Reevaluation of Algorithmic Basics for ZUPT-Based Pedestrian Navigation

During the last 10 to 15 years, pedestrian navigation based on zero velocity updates (ZUPT) has become very popular. One of the main reasons for this is the increasing availability of small, low-cost inertial measurement units. However, the processing of the data from these units for pedestrian navigation is almost exclusively based on algorithmic features that originate from classical inertial navigation with high-grade sensors. In addition, the historical background of the ZUPT approach presupposes also sensors with high accuracy. This leads to the problem that neither the ZUPT approach nor the algorithmic features mentioned are consistent with the accuracy level of the inertial measurement units used normally for pedestrian navigation. Therefore, the question arises of whether the usual algorithmic basics and numerical procedures employed by ZUPT-based pedestrian navigation are adequate. Supported by a literature review showing the state of the art, this study investigates the effect of basics such as the system states and the usage of Runge-Kutta method on the navigation accuracy from pedestrian inertial data aided by ZUPT. To this end, a comparative processing of a well-known, published dataset of a short walk and of own data from the authors using different data fusion algorithms was employed. The important results of this study are that the often-used omission of inertial sensor biases cannot be recommended, that a continuous-discrete Kalman filter in combination with a Runge-Kutta algorithm performs better than traditional filter configurations, and that total navigation states are an interesting alternative to classical navigation error-states.


I. INTRODUCTION
For more than two decades, the availability of inertial measurement units (IMUs) being designed as microelectronic mechanical systems (MEMS) has strongly increased. These devices offer a favorable price-performance relationship combined with low weight, volume, and power consumption. Furthermore, their use benefits from cheap and small microprocessors of growing performance for signal processing. This fact is reflected in the remarkable emergent popularity of IMUs for biomechanical applications [1]. One such usage is pedestrian inertial navigation.
The associate editor coordinating the review of this manuscript and approving it for publication was Halil Ersin Soken .
On the other hand, there has been great success in designing increasingly small satellite navigation (GNSS) receivers for about three decades. As an essential application, these units were early adapted to the guidance of pedestrians but are of course restricted to areas of good satellite visibility. IMUs are particularly suitable for overcoming the latter limitation and have led to considerable research efforts to expand pedestrian navigation to urban and indoor areas for many years. These activities effected not only to a strong miniaturization of pedestrian inertial navigation systemsfrom the initial size of a handcart [2] to the integration into a shoe heel [3] -but also to a great variety of methods and technical approaches used. Nevertheless, this development exhibits some algorithmic basics that have remained virtually VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ unchanged for many years. These basics are the focus of the paper at hand. To first introduce and classify the algorithmic basics, the paper continues with a typing of pedestrian inertial navigation systems, with explaining the background of the so-called ZUPT method, and with an overview of relevant research topics. The introduction concludes with an overview of the study presented in the main body of the paper.

A. PEDESTRIAN INERTIAL NAVIGATION SYSTEMS
During the years three types of pedestrian inertial navigation systems developed. The first one is based on wearable sensor networks consisting of many MEMS IMUs and magnetometers. These sensor units are distributed over the human body, which is simultaneously modeled as a multibody system like a chain [4]. Such networks serve as motion capture systems (not only for pedestrian motion) and exist in numerous variants [5]. Nevertheless, their inertial units are mostly designed as local, single Attitude and Heading Reference Systems (AHRS) [6]. Using the outputs of one IMU and of one magnetometer, each AHRS determines the attitude (i.e. the pitch, roll, and yaw angle) of that body part (foot, lower leg, thigh, etc.) it is attached to. Based on this, the motion of the multibody system is then composed from the partial movements of all single body parts. (Remark: Alternatively, the signals of all IMUs could be fused as a whole and together with the kinematical structure of the entire multibody system [7], but a realization of this integrated approach is not yet state of the art.) The second type is a strongly simplified version of the first one. It is based on a single AHRS, a kinematical model of a two-dimensional pedestrian motion on a horizontal surface, a step detector, and an algorithm for step length estimation. The pedestrian motion is reconstructed from the length, the direction, and the summation of all single steps [8]. This approach is called Pedestrian Dead Reckoning (PDR), and it aims, as an example, to smartphone applications. Because of its underlying simplifications, its accuracy potential is limited. Although PDR as well as wearable sensor networks are not subject of the paper at hand, some findings of the following sections may nevertheless advance also these approaches.
The third type, which is especially popular, pursues also a relatively simple inertial sensor arrangement: Generally, only one foot of a pedestrian is equipped with a single MEMS IMU being strapped down (abbreviated: ''strapdown'') or rather fixed to the shoe or the ankle. This measurement approach intends to reconstruct the pedestrian movement from recording only the motion of the foot. In this case, however, the IMU is not the signal source for an AHRS but for a complete inertial navigation system (INS), which determines additionally the three-dimensional position and the velocity of the IMU. That technical form or type of pedestrian inertial navigation is the subject of the paper at hand; and it includes the so-called Zero Velocity Update (ZUPT) aiding.

B. THE BACKGROUND OF ZUPT
It is well known that the mathematical procedures inside an INS are subject to an inherent numerical long-term instability. The type of pedestrian inertial navigation being outlined last forms no exception to this effect, which leads to growing errors of position, velocity, and attitude. For limiting these inaccuracies, a so-called aiding of the INS is required. Traditionally, radio navigation devices like GNSS receivers or doppler radar provide the necessary aiding data. Another aiding approach, which can be utilized especially for pedestrian navigation, takes advantage of the mentioned ZUPT technique. This method is based on the idea that during the resting phases (i.e. stance phases) of the shoe (or, more precisely, of the IMU), artificial velocity measurements of zero values can be generated to aid the INS.
Originally, ZUPT was developed in geodesy for surveying with precise standalone inertial navigation systems, which stopped at regular, well-defined points for a few minutes to update the INS. Applications date back to the 1970ies [9]. The tradition of ZUPT-based pedestrian navigation is much shorter, only about 15 years [10], and it has been especially promoted by a few papers, such as that of Foxlin [11] and Jiménez et al. [12], as well as the open-source software OpenShoe of Nilsson et al. [3]. In contrast to geodesy, the detection of the -short -resting phases of the inertial sensors is more difficult, as foot-mounted IMUs always show a remaining, inevitable motion during stance phases. Although a certain number of heuristic procedures have been developed to identify stance phases, their reliability is limited [10], [13]. Consequently, aiding by ZUPT in pedestrian navigation has significantly lower precision than that of geodesy.
Consequently, the accuracy level of ZUPT-based pedestrian navigation differs from classical inertial surveying, not only in the IMU precision grade but also in the ZUPT quality. Nonetheless, the mechanization of ZUPT-based pedestrian navigation systems, that is the mathematical description for processing the IMU output and aiding measurements, is normally taken from classical strapdown systems with much higher accuracy levels [6]. The mechanization is simplified only in a subordinate step for considering the lower accuracy of consumer-grade IMUs. This modification, as used in [3], [11], and [12], is explained in detail in a new textbook by Wang and Shkel, who additionally review, systematize, and analyze the state of the art of ZUPT-based pedestrian navigation [14]. Nevertheless, the simplifications do not affect the general system features of classical, high-grade strapdown systems -features that are also very familiar in integrated navigation, such as combining an INS with a GNSS receiver.
Against this background, one can ask to what extent the usual mechanization of ZUPT-based pedestrian navigation systems is at all compatible with the low accuracy level of the ZUPT aiding and the MEMS IMU precision. This issue is the main question of the following sections.

C. RESEARCH TOPICS ON ZUPT-BASED PEDESTRIAN NAVIGATION
To assess the importance of that question, the study being presented below included an extensive literature research. Section III presents details about the procedure of this paper review and some special statistical results. To classify the algorithmic basics mentioned above, however, a general overview of current research topics related to the advancement of ZUPT-based pedestrian navigation is already given in the following paragraphs. Due to an extensive, broad literature situation, the papers cited for this overview can, however, only be a selection of preferably recent articles. Furthermore, a subdivision into six thematic groups is made for this overview: Stance phase detection: To ensure a good aiding by ZUPT, a clear detection of stance phases of natural walk is an obvious topic for research. Besides the classical techniques like those of Skog et al. [13], the available methods for this task occupy meanwhile a wide space in the literature of pedestrian navigation and range from rather general approaches like artificial intelligence [15], machine learning [16], and neural networks [17] to more targeted techniques like motion classification [18], detection of special environments like escalators [19], decomposition of the single steps [20], or usage of insole sensors [21].
Complementing ZUPT aiding: Another obvious research area is the usage of additional aiding techniques. Such methods are particularly numerous in the literature and aim at compensating shortcomings of ZUPT but also at enhancing the reliability and accuracy of indoor navigation in general [6]. At this, an important methodological subgroup is characterized by waiving additional sensors. Using motion typing or motion constraints, the techniques of this subgroup try instead to identify e.g. the walk on horizontal surfaces, along corridors, or between marked waypoints in buildings. A constant height or specific heading is then used as a derived, also artificial measurement [22], [23], [24]. Furthermore, the aiding technique Zero Angular Rate Update (ZARU) is a central member of this subgroup. This approach is addressed in more detail in Subsection II.E below. A second subgroup covers traditional aiding devices like radio receivers for GNSS (at indoor-outdoor transition areas), UWB, WLAN, Bluetooth etc. [25], [26], [27] as well as magnetometers [28] and cameras for visual odometry and ambient sensing [29]. Also, map matching approaches for positioning, collision avoidance, and identifying preferred walking directions pertain to this subgroup [29], [30].
IMU: The next group of research topics relates to the quality of the IMU data. To increase the IMU accuracy, calibration [31] or modeling systematic sensor errors [32] form a natural approach. Arrays of redundant gyroscopes and accelerometers with non-orthogonal input axis are also known for increasing the navigation precision [33]. Other research is concerned with optimal IMU positioning on the shoe [34] and optimal anti-aliasing filtering of the IMU data [35]. A classical topic of these research approaches is the modeling of inertial sensor biases. This problem is addressed in the Subsection II.F below.
Second inertial unit: The fourth group of research activities borrows elements from wearable sensor networks. Typically, both feet are equipped with an IMU each, and the IMU outputs feed the data processing of a dual INS. The INS aiding comprises not only the ZUPT technique for each foot separately but also the utilization of the distance between the two feet. The latter aiding is implemented by a pure mathematical evaluation of distance constraints [36] or by using additional measurements of sensors such as a pair of cameras [37], radar [38], or ultrasound units [39]. An alternative variation of this dual IMU configuration is to instead wear one unit on the chest to form a local AHRS. The attitude of the AHRS is regularly transmitted down to the ZUPT-based INS [40].
Applications outside the pedestrian area: A smaller group of ZUPT-based research focuses on transferring this type of inertial navigation to mechatronic applications such as production engineering [41], robotics [42], or specialized vehicles [43].
System aspects: Finally, also a rather small group of research topics addresses issues such as system accuracy [44] or fusion techniques for the IMU and aiding data [45], [46].
To these results of the literature review must also be added the already indicated, important observation that -despite the wide thematic spectrum of research topics -almost all papers share some common methodological characteristics. These features relate to the numerics and the system design of ZUPT-based pedestrian navigation and are the algorithmic basics mentioned above. They fit thematically into the last group of research topics and are explained in Subsections II.A through II.D. However, their current use appears like a standard, like a norm that is not questioned. This circumstance may have its origin in the fact that the ZUPT-based pedestrian navigation can be attributed to a few studies such as [3], [11], and [12], while tutorials such as [47] consolidated this effect.
Against this background, Section II presents the algorithmic basics and possible alternatives in detail. In addition, it explains the ZARU method and a modeling of the inertial sensor biases, as these two topics will be investigated in terms of numerical improvements related to the algorithmic basics. Section III follows with a statistical analysis for the literature review to show how established the algorithmic basics are in the state of the art. Section IV focuses on processing pedestrian inertial navigation data from experiments. This is done to evaluate numerical improvements with respect to alternatives to the algorithmic basics. Finally, Section V discusses the results of the data processing and derives recommendations for pedestrian inertial navigation.

II. ALGORITHMIC BASICS
The algorithmic basics mentioned repeatedly concern the system architecture, system states, Kalman filter type, and numerical integration of ordinary differential equations required for determining the motion of a strapdown IMU. The next subsections each describe one of these four topics and present alternatives to the current state of the art. Two subsections on ZARU and modeling of inertial sensor biases complete these explanations.

A. SYSTEM ARCHITECTURE
The following discussion is based on Fig. 1 and 2, which have been published in a similar form by Wagner and Wieneke [48], who investigated system simplifications for integrating a strapdown INS and a GNSS receiver for airborne applications. These simplifications are now tailored to the pedestrian inertial navigation.
First, two coordinate systems, as illustrated in Fig. 3, have to be introduced: For navigation, the gravity-fixed system (or frame) n shall be used. It represents a Cartesian north-eastdown system with coordinate axes n x, n y, and n z. Furthermore, it is assumed that the origin of this system is set as required by the respective navigation task and that the local surface of the earth is flat. This simplification is justified by the fact that customary MEMS IMUs for pedestrian navigation cannot sense low angular velocities like Earth's rotation rate. Therefore, system n can be considered as an inertial frame. The other coordinate system b is fixed to the respective moving body or object, which is typically a shoe or an ankle in pedestrian navigation. In addition, an IMU is fixed to this body, and the IMU measurement center is the origin of system b. The coordinate axes of the system are b x, b y, and b z.
Next, several vectors and subvectors like position r (Fig. 3) have to be introduced (vectors are indicated by bold letters). All the vector components are functions of time t. In addition, the components of all vectors with a physical meaning have to be represented using the axes of one coordinate system. The used system is indicated again by the prefixed index b or n.
aiding vector y = velocity vector n v .
In Fig. 1, the block x represents the moving body, which is subject to an acceleration and angular rate vector (equation (1a)). From this input u, a motion results, which can be characterized by a state x consisting primarily of the position, velocity, and angular attitude of the object (equation (1b)). If x must be known but cannot be measured directly, it must be reconstructed. For this purpose, an IMU measures u, and a strapdown algorithm calculates an estimate (notation^) of x in block z from the IMU signals [6], [14].
The IMU consists usually of three accelerometers and three gyroscopes. (In the following, the term gyro is used as an abbreviation for gyroscope.) Both sensor sets have pairwise orthogonal input axes. To modeling the IMU error characteristics, x is traditionally enhanced by a bias for each of the six sensors, as shown in (1b) [14], [48].
Historically, the block z was not a strapdown system but a classical stabilized mechanical platform (small enough for a vehicle but much too large for a shoe). To increase its accuracy numerically, the only possibility was emulating the error behavior of the INS in a separate navigation computer and calculating an estimate of the error δ of thex. This is performed in block | (and is used to update the INS, dotted line). Today, this separation is still standard, but no longer necessary, as all strapdown and error model computations can be performed using one navigation computer with an IMU as the data source. (Remark: For δx, the quaternion q of x should be replaced by the vector α of the three errors of the Euler angles φ, θ, and ψ for roll, pitch, and yaw. Roll and yaw are illustrated in Fig. 3, pitch is the turn around an intermediate axis between the n y and the b y axis.) Furthermore, former navigation computers had a very limited accuracy as they were based on processors of 16 bits or less. The same applies to the output of classical INSs mentioned. From that viewpointx can be understood as the more significant digits and δx as the less significant digits of the common estimatex+δx of x. Such a separation, however, is also not necessary when using the numerics of a modern microprocessor with 64 bits.
Based on these considerations, blocks z and | can be merged for strapdown systems if the strapdown calculations are performed with comparatively high accuracy. The merger is represented by block z in Fig. 2.
The same line of thought is possible for the right part of Fig. 1. To limit the aforementioned inevitable numerical instability of an INS or a strapdown calculation, an aiding of x is required: First, a suitable measuring device y with output y is modeled by an appropriate function of x. By means of this function, an estimateŷ can be calculated fromx in block {. The difference y −ŷ is a measure for x −x. Clearly, y −ŷ contains only the less significant digits of the aiding part. The error δŷ of this difference can be estimated from δx using a suitable error model (block }). The comparison of δŷ with y −ŷ can be used for a suitable feedback device (block~) to correct δx (which, in turn, correctsx via the dotted line). In addition, the two comparison points on the right side of Fig. 1 lead to the expression y−(ŷ+δŷ) showing that the more significant digits and the less significant digits of the aiding estimate can be merged for microprocessors with significant accuracy. The merger is represented by block { in Fig. 2.
The structure of Fig. 2 and the lower part of Fig. 1 (blocks z to~) represent the principle of an observer; accordingly, the theory of Kalman filters is normally used for designing the feedback unit~. Fig. 2 not only has the advantages of a simpler system structure but also of less computing effort. It also leads to a noteworthy reduction in the estimation error covariance of the Kalman filter [48]. This indicates a gain in the filter stability.
With respect to ZUPT-based pedestrian navigation, two additional aspects are important.
First, the error models in blocks | and } result from a linearization of the equations used inside blocks z and {. Therefore, they are bound by the premise that δx and δŷ are small. For the high error level of ZUPT-based pedestrian navigation, it is not certain whether this assumption is still reasonable. Fig. 2 shows the direct usage of nonlinear observers, such as the extended Kalman filter (EKF).
Second, the error models in blocks | and } havex as a linearization reference (see the dashed lines in Fig. 1), which is kept intermittently constant during the operation of the scheme shown in Fig. 1. This is acceptable for vehicles such as ships and aircraft with moderate dynamics, but is not acceptable for abrupt motion phases such as the impact of a foot on the floor. Such quick changes in x can be handled easily using Fig. 2 in combination with the theory presented in Subsection II.D.
To assess the impact of the system architecture described above on ZUPT-based pedestrian navigation, Section IV compares of the processing of real data with the system shown in Fig. 1 and the system in Fig. 2.

B. STATE AUGMENTATION BY IMU SENSOR BIASES
To further develop the discussion of Figs. 1 and 2, the mathematical content of the blocks z to } must be considered. The block z contains a set of ordinary, coupled, and nonlinear differential equations:ẋ which must be numerically solved. The block { consists of the following set of algebraic equations: which simplifies toŷ =v together with y = 0 in the case of pure ZUPT. The blocks | and } are formed with Jacobian matrices F and H of f and h with respect tox: In the last two equations,x and u have the function of linearization reference.
To achieve robust performance of the systems in Fig. 1 and 2, it is necessary that y contains sufficient information about every element of x. Then, the aiding acts on all elementsx and δx, respectively. In this case, x is completely observable by y. However, assuring full observability is not trivial and requires that the matrix has a full rank, where n is the number of elements of x. This condition must not be fulfilled permanently, but at least during regular intervals to provide satisfactory system performance. An example of aircraft navigation is the aiding of an INS with a single-antenna GNSS receiver. In this case, turning flights are required occasionally [49]. It is a well-accepted fact that aiding an INS using a pure ZUPT technique generally cannot reach full observability. As the authors could not find a simple proof for this circumstance, the appendix contains a compact mathematical derivation, which provides a more detailed insight and shows that the first three lines of B remain empty. This implies that position r is not observable. Similar shortcomings also exist for the yaw angle, gyro bias of the vertical axis, and accelerometer biases of the horizontal axes if one axis of the IMU is (as mostly usual) oriented vertically.
Incomplete observability means that the estimated position r diverges, and other state variables also show symptoms of instability. Against this background, Skog et al. [13] proposed omitting the six inertial sensor biases of b a and b ω in x and VOLUME 10, 2022 δx, which reduces the degree of missing observability by deleting three unobservable elements of x, but prevents the compensation of the biases.
Skog et al. argued that the biases are a minor source of error. However, the literature lacks a comparison of the extent to which omitting sensor biases affects system performance. Therefore, such an evaluation is made in Section IV.

C. DISCRETE-DISCRETE AND CONTINUOUS-DISCRETE KALMAN FILTER
The next point of discussing Fig. 1 and 2 is the assumption made so far that the systems depicted operate continuously over time. For blocks y, {, }, and~however, digital data processing technology is common today. Therefore, these blocks operate in a sampling mode, which is based on discrete time points t i . By contrast, the block x represents a process that is principally time-continuous. For such circumstances, the theory of the extended Kalman filter offers a continuousdiscrete filter variant, which can be outlined as a recursive algorithm (with x(t i ) =: x i , y(t i ) =: y i etc.): 1) Predictor step from t i to t i+1 : • numerical solution of (2) in block z with the initial vector x i and result x i+1 , as well as for the estimation error covariance matrix P(t) in the block~, where P i is the initial value, P i+1 is the result, Q is the covariance matrix of the system noise w(t) of the block x, and G is the matrix describing the effect of w on x [50]. Note: The model of block x which corresponds to (2) isẋ 2) Corrector step at t i+1 (in block~): • calculating the Kalman gain matrix K, • then correcting of δx andx, respectively, . 1) , • and updating of P (more accurate Joseph form), where R is the covariance matrix of the measurement noise v(t) of the block y, and I is the unit matrix of appropriate dimension [50]. Remark: Matrices R and Q are typically constant over time. They help for tuning the Kalman filter to adjust it to the sensor quality of u and y. In addition, ZUPT-based pedestrian navigation requires careful tuning to adapt the filter to motion characteristics that vary with the walking task and personal walking style [51].
Although the continuous-discrete extended Kalman filter is well known [52], [53], in ZUPT-based pedestrian navigation, the discrete-discrete filter variant, in combination with Fig. 1, is common [14]. Using the step size t = t i −t i−1 in time, the discrete-discrete filter equations approximate the continuous model of (4) and (7) as These equations alter the time-continuous block | to a seemingly discrete-time unit. However, it is only an artifice. Thus, comparisons of the different performances of the continuous-discrete and discrete-discrete Kalman filters are required to assess the effects of approximations (12) and (13). Such studies exist for some mathematical cases [52], but investigations of the special conditions of ZUPT-based pedestrian navigation were not found by the authors. Therefore, Section IV presents such a comparison.
Remark: In (13), the term G i Q G T i t is often replaced with Q i or Q t 2 [14], [47].

D. RUNGE-KUTTA NUMERICAL INTEGRATION
The conditions of the continuous-discrete extended Kalman filter and the classification of (12) require additional remarks about the numerical solution of ordinary differential equations (2) and (7). Owing to the updates of (10) and (11),x and P are discontinuous functions of time at all t i . To numerically solve the initial value problems of (2) and (7), one-step methods are therefore preferred to multistep methods [54]. Many such approaches exist that differ from each other with respect to accuracy and computational cost. Fig. 4 outlines three possibilities of the one-step method for an arbitrary function s(t). The first is the (explicit) Euler method. It is especially important as (12) is a realization of this approach. The method assumes that the time-derivativeṡ of s is a function of s: s(t) = f (t, s). Starting at time t i , the Euler method calculates an approximation of s with error e (see Fig. 4) at time t i+1 = t i + t: Error e is inevitable, but it can be reduced by iteration. For this, (14) is supposed to be the first prediction of s(t i+1 ), denoted by the upper-right index p,1. From s p,1 i+1 , an estimate of the time derivative of s at t i+1 can be calculated using f . The mean value of this derivative and of the derivative of (14) is then used as the approximated, averaged derivative for the entire time interval. This leads to a better accuracy of s(t i+1 ) but can also be used as the beginning of k iterations for further improvement: If k = 1 is chosen, then the Heun method results. This is the 2 nd possibility in Fig. 4 and represents a Runge-Kutta approach of 2 nd order.
Instead of using higher values of k, more efficient Runge-Kutta algorithms of a higher order can be used. The best-known variant, a Runge-Kutta form of 4 th order, is also outlined in Fig. 4. It starts similar to the Heun method; however, it uses only half of the step size and erects an additional grid point in the middle of the interval at t i+1/2 = t i + t/2. Next, an iteration is used to calculate the derivative of s at this extra point before computing the derivative at the right end of the interval. Finally, all calculated derivatives of s are combined to a weighted average of the entire time interval, and this mean value is used to improve the calculation of s(t i+1 ) [50].
As the Runge-Kutta theory is well developed, a large family of algorithms, such as that of (16), exist [54]. Fenta and Derese [55] presented, for example, an algorithm of 6 th order, which is utilized together with (16) and (14) in Section IV. Owing to its size, it is not reproduced here, and it roughly duplicates the computing effort compared with (16).
In addition to analyses for subsections II.A and II.B, Section IV presents also a comparison between the Euler and the Runge-Kutta method for ZUPT-based pedestrian navigation. This is combined with the comparison of a continuousdiscrete and discrete-discrete Kalman filter by assigning the Euler method via (12) to the discrete-discrete variant, and the Runge-Kutta method to the continuous-discrete variant.
Remark 1: The above equations are also valid if f and s are vector functions.
Remark 2: The above discussion assumes that u is continuous in time, but is sampled by the IMU only at each point t i . Therefore, u must be interpolated if it is used at the additional grid points of the Runge-Kutta method while numerically integrating (2), (4), and (7). This requires that the IMU outputs be instantaneous measurements of the components of a and ω. This is customary for modern MEMS units and differs from classical inertial sensors, which typically generate velocity and angular increments as outputs. In the latter case, blocks z and | should be handled like time discrete units with, however, more accurate sculling and coning algorithms instead of (12) [56], [57].
Remark 3: It is a matter of course that the step size t has to be sufficiently small for all numerical methods explained above. For flight tests, a minimum sampling frequency 1/ t of 100 Hz is an appropriate practical value [58]. However, the high dynamics of a foot or shoe hitting the floor require a higher rate. Based on a thorough frequency analysis, Munoz Diaz et al. [59] recommended a minimum rate of approximately 250 Hz.

E. ZERO ANGULAR RATE UPDATE
As introduced above, the ZUPT aiding can be complemented by the ZARU method [12]. This simple, well-known technique is explained in the next subsection, since the alternative algorithmic approaches of subsections II.A to II.D indicate numerical improvements -and the condition of better accuracy suggests a re-evaluation of ZARU.
The aiding by ZARU is also based on the idea that the IMU is at rest during stance phases. Therefore, the IMU generates gyro measurements that are nominally zero during these periods. This means that during such phases, the real gyro signals consist only of the respective possible sensor bias, which, in practice, adds itself to the gyro output. Hence, the IMU measures directly b ω during stance phases. For using ZARU aiding, (1c) and (3) must be accordingly expanded as follows: According to the Appendix, ZARU provides the missing observability of the gyro bias of the vertical axis. This stabilizes the Kalman filter, to a certain extent. The algorithmic implementation of ZARU is also very simple, Therefore, ZUPT-based pedestrian navigation is often combined with ZARU. However, the benefits of ZARU are uncertain. Some papers show accuracy enhancements, e.g. [12], while authors of other papers express reservations [51], [60]. Hence, as mentioned, Section IV considers the possible accuracy improvements of ZARU in conjunction with the other algorithmic features described above. VOLUME 10, 2022

F. GAUSS-MARKOV MODEL FOR IMU SENSOR BIASES
Owing to the partially poor observability of IMU sensor biases, there are modeling approaches to improve the estimates of b a and b ω . Also, this research topic was already mentioned and should be re-evaluated under the precondition of numerical improvements.
Traditionally, all six IMU biases are individually described in (2) with the same simple type of differential equation: where w is a single white noise signal, and there is no correlation between the six white noise signals from all these equations. This simple model results in an angular random walk for each gyro, and a velocity random walk for each accelerometer. (Both random walk types are important quality indicators of inertial sensors.) Another well-known and easy-to-implement option is to model the bias by a Gauss-Markov process of 1 st order with time constant T :ḃ Quinchia et al. [61] showed that (19) leads to a better accuracy than (18) when integrating a GNSS receiver and an INS with a MEMS IMU. As the authors of the paper at hand could not find an equivalent assessment for ZUPTbased pedestrian navigation, they added such a comparison to Section IV.

III. REVIEW OF LITERATURE
Section I already discussed that the practical implementation of ZUPT-based pedestrian navigation contains key elements that are handled like a standard. After introducing these elements in Section II, an estimation lends itself how prevalent they are in state-of-the-art. Also to assess the significance of the numerical results from Section IV, the authors statistically analyzed therefore the literature review of Subsection I.C. This investigation was based on the following sequential steps.
• Several publisher websites were selected to cover important journals for pedestrian navigation and inertial sensors such as GPS • From these papers, all articles were removed that were erroneously selected, which contained reviews instead of describing a specific research project, which were not detailed enough for an analysis, were based essentially on a handheld device, were pure AHRSs, followed an observer principle that was not comparable with a Kalman filter, or referred to unpublished papers.
• After sorting out these papers, 186 remained. Their research topics were already summarized in Subsection I.C. The content and the number of papers were deemed sufficient to allow a statistical assessment of the state of the art.
• All these papers were examined with respect to six questions. In case one article was not detailed enough or was too ambiguous regarding a single question, this question was skipped for that paper. These questions were guided by the subsections of Section II. They are now listed together with the possible answers and the number of papers which could be assigned to the answers. Assignment of the answers was easy for most papers; occasionally, a guess seemed reasonable. 1) Were the sensor biases b a and b ω used for x? -Yes: 119 -No: 60 -Not specified: 7 2) Did the authors use mainly the pseudo-discrete Euler method of (12) or a numerical procedure of higher order such as Runge-Kutta for time-continuous problems? -Euler method: 160 -Higher-order method: 0 -Not specified: 26 3) Which sample rate 1/ t was used for u (following Remark 3 in Section II.D)? -1/ t < 100 Hz: 16 -100 Hz ≤ 1/ t < 250 Hz: 86 -1/ t ≥ 250 Hz: 27 -Not specified: 57 4) Did the authors use an error-state filter (Fig. 1) or total-state filter (Fig. 2) Following the number of the questions, the result of this review can be interpreted as follows: 1) The situation regarding the omission of the sensor biases is heterogeneous. Only a part of the papers complies with the recommendation of Skog et al. [13]. This circumstance lightens up with the results of Section IV: Omitting the sensor biases is indeed a simple means of limiting the estimation error. However, using the system structure of Fig. 2 together with a Runge-Kutta algorithm of sufficiently high order achieves at least the same effect and additionally allows the retention of bias compensation, as is traditional in inertial navigation. 2) Numerical procedures of higher order are seemingly unusual or unknown, and a significant number of the authors even considered a specification of the numerical method unnecessary. This is a critical finding in light of the results in Section IV, and the finding is also surprising compared with the reputation of sculling and coning algorithms as well as with the reference to Runge-Kutta algorithms in well-known textbooks such as [6] or [50].
3) The critical outcome of the answers to Question 2 intensified when the IMU sample rate was considered. A predominant part of all papers (85%) does not fulfill the minimum frequency requirement of [59] or does not provide any frequency information at all. There is a clear preference for a sample rate of about 100 Hz, which may result from the preset value of the customary IMUs [12] or from the practical values of other INS applications. 4) The usage of the total-state filter in Fig. 2 is not unknown, but is much less established than the errorstate filter. Furthermore, it is interesting that the papers with total states do not elaborate the reasons for choosing this filter variant instead of the error-state filter. Also, comparisons between both filter variants are missing. This effect could be caused by the fact that the authors of these articles seem to be ''lateral entrants'' to inertial technology as they come from areas such as computer science or forestry [62], [63]. 5) The papers selected for the literature review show a heterogeneous situation in the use of ZARU. This reflects the uncertain benefit of ZARU for pedestrian navigation.
6) Traditional random walk for modeling sensor biases still dominates the state of the art. A significant advantage of the Gauss-Markov model is not obvious.

IV. EXPERIMENTAL RESULTS
After explaining the algorithmic basics in Section II and assessing their prevalence in the state of the art in Section III, a numerical investigation of these key elements is an essential next step. For this task, real pedestrian IMU data were processed using several software variants. ZUPT and ZARU were the only aiding techniques used in block y to exclude the influence from other sensors. The IMU data were obtained from experiments conducted by the authors and from the OpenShoe group (see below).
No. 7: The Runge-Kutta algorithm of (16) is replaced by the higher-order algorithm of [55].

B. ORIGINAL OPENSHOE DATA
The first trial to check the performance of all seven software variants was performed using the real test data included in the software package from www.openshoe.org. These data were obtained using a MicroStrain 3DM-GX2 IMU mounted in the sole of the right-side shoe of a test person [13]. The IMU measurement ranges were adjusted to ±18 g and ±1,200 deg/s. The sampling interval was t = 0.004 s. This value fulfils the minimum requirement of approximately 250 Hz for the sampling frequency of pedestrian navigation [59]. Due to the relative smooth character of the IMU signals, it seems that an additional anti-aliasing filtering had been applied to the data. VOLUME 10, 2022 However, details about this signal preprocessing, could not be found in the documentation of the OpenShoe package.
The OpenShoe file package contains several sets of similar real IMU data. To test the seven software variants described above, one was selected, which is located in the folder ''Measurement_100521_44''. The main results of this data evaluation are presented below. Beforehand some statements regarding Kalman filter tuning should be made: First, the filter matrices R and Q were mainly the same for all seven variants to obtain a definite basis for comparing results. Both matrices have a diagonal shape, and the values of their elements are included in the downloaded MATLAB R code. (There is only one set of values in the file package.) The only change regarding Q was made for the transition from (13) to (7) in variant no. 3 when Q had to be replaced by Q/ t to ensure the equivalence of both equations.
Second, ZUPT-based pedestrian navigation requires careful filter tuning, as previously mentioned. This was the main reason for using the given values of R and Q from the OpenShoe MATLAB R code. During the processing of the IMU data, however, it turned out that OpenShoe tuning aims mainly at good results of the position only, but not of the height. This can be seen below by examining Fig. 5-7.
Third, the transition from variant nos. 4 to 5 entailed an extension of R for ZARU. A diagonal shape with three identical values for the variance of ZARU noise was assigned to this matrix part. Systematic sequential attempts were made to find a good choice for the numeric values of these diagonal elements. The selection criterion was a reduction in the position error at the endpoint (according to the determination of Q and the remaining part of R).
Fourth, the transition from variant no. 5 to no. 6 entailed the introduction of six time constants T -one for each IMU sensor model, following (19). The numerical values were different for each sensor owing to the different characteristics of the mechanical loads of the sensors (e.g., a much higher effect of gravity on the vertical accelerometer axis than on the horizontal axes). Again, systematic sequential attempts were made to find good choices for the value of T for each sensor. The selection criterion was once more a reduction in position error at the endpoint.
Finally, the software of the ZUPT unit for detecting the stance phases in the blocks y of Figs. 1 and 2 remained unchanged for all variants. This includes the parameters for the detector tuning. The detector type was the stance hypothesis optimal detector (SHOE) [13].
The discussion of the results of the software variant test is limited to the position, height, and yaw angle, as these motion components are generally not observable using ZUPT and ZARU (see Appendix). Therefore, their error rates should be particularly characteristic. Fig. 5 shows the ground track of the motion recorded using the OpenShoe data. (For r, the north-east-down coordinate system n with n r = [ n r x n r yn r z ] T is used, but the prefixed index n is omitted below for a more compact view.) It reflects a walk in a closed-loop trajectory with an approximate shape of 8.
The starting and endpoint were identical. This was realized during the test using a shoe imprint [3] and included the same shoe attitude. The entire walk lasted for approximately 42 s. The results of no. 3 and no. 4 cannot be seen, as they are covered by the curve of no. 1, but no. 2 obviously led to a significant decrease in long-term accuracy. The results of nos. 5 to 7 are not included in Figs. 5, 6, and 7 in order not to overload the diagrams by curves with no essential new information. Fig. 6 shows a close-up view of the estimated endpoint. The result of no. 2 is outside the range of the diagram, but the results of all other variants are close together and at similar distances (see the circles) to the true endpoint. Considering the height, Fig. 7 reveals additional significant details. First, the accuracy of no. 1 is one order of magnitude worse than that of the position. In addition, the accuracy decreased again for no. 2. Variants no. 3 and no. 4, which are difficult to distinguish, not only compensate for this loss of precision but also significantly increase accuracy.
For all seven software variants, Fig. 8 contains in the upper two diagrams the distance error and the height error at the end point (i.e. the difference to the starting point). In both cases, it is visible that variant no. 3 is better than variant no. 1, and the accuracy can additionally be slightly enhanced by variant no. 4. ZARU and Gauss-Markov modeling in no. 5 and no. 6 also improve the position accuracy but impair the height estimate. This last effect is probably a consequence of the filter tuning strategy. The usage of a Runge-Kutta algorithm of the 6 th order apparently yields no consistent improvement, too.
The bottom diagram in Fig. 8 shows the yaw difference between the start and end values. An interpretation of this result is difficult because the values of all variants vary obviously around an offset. It seems as if there is a systematic yaw deviation of the shoe between the start and end. This difference may result from a certain clearance of the shoe imprint.
As an interim conclusion, it can be stated that omission of the inertial sensor biases is only a limitedly effective means against filter instability owing to the incomplete observability of δx. The use of a better numerical integration instead of the Euler method is obviously a more appropriate tool and challenges the argument of Skog et al. [13] that the biases b a and b ω are of little importance. This finding also illuminates the heterogeneous result regarding b a and b ω in Section III: part of the authors of the selected papers followed the suggestion of Skog et al. to omit the inertial sensor biases, and the remaining part of the authors apparently tried to limit the filter instability by other methods.
Another interim conclusion is that the use of total states instead of error states can lead to at least small accuracy improvements, whereas the question of whether ZARU or Gauss-Markov modeling has a benefit remains open. In addition, the effort of the Runge-Kutta algorithm of 6 th order instead of 4 th order does not lead to consistent improvements.

C. ADDITIONAL EXPERIMENTS BY THE AUTHORS
The last statements are only the results of a single experiment, and a better assessment requires further tests. Therefore, the authors executed additional experiments. To promote the generality of the findings from these tests, they deliberately introduced several differences compared to the OpenShoe trials above: • Use of newer IMUs with different accuracy levels and different measuring ranges of ±8 g and ±500 deg/s in common: -Bosch BMI085 [64] and -Analog Devices ADIS16465-2 [65].
• Higher sample rate of 2000 Hz. • No anti-aliasing filtering except for the inherent filtering of the IMUs.
• Different shoe type ( Fig. 9) with authors as test persons. • Test track for a different trajectory ( Fig. 9 and 10).
• Position, height, and yaw checks not only at the endpoint but also at several intermediate positions (Fig. 10). Fig. 9 shows the laboratory room with the test track. This is a rectangular field with a special cradle in each corner. These four devices provided a position and attitude reference for a shoe with the outer shape of a cuboid. The shoe is also shown in Fig. 9 and carries in the forefoot area the two IMUs being plugged into breakout boards of the manufacturers. Fig. 10 shows the horizontal part of the estimated trajectory of one test walk. (Only the results of two software variants are shown to avoid overloading the diagram.) This test walk was selected because of its smooth test sequence. Its route started in the lower left corner, moved to the upper left corner, then to the upper right corner and subsequently to the lower right corner. Then, two diagonals followed, with a repetition of the segment from the upper left to the upper right corner in between. The walk ended with the second diagonal to the lower left corner. The covered distance was approximately 25 m, the duration was approximately 120 s.
Before discussing the results of the selected walk, three general remarks about the parameter tuning of the ZUPT detector and Kalman filter must be made.
First, the ZUPT detector type was again SHOE as employed before; however, the usage of significantly other footwear and different IMU types required to adapt the detector parameters in order to achieve a usable step detection. The variances for acceleration and angular rate noise were adapted to the sensors used. This implied similar values as for the OpenShoe IMU. However, the window size for identifying the resting phases was partly increased, and the accompanying threshold level had to be strongly aggravated by lowering the value for the admitted residual motion.
The parameter modification of the ZUPT detector reflects a stronger narrowing of the stance phases. This approach turned out to be appropriate because of the omitted antialiasing filtering; and it was possible because of the large, stiff footprint of the test shoe (leading to temporarily steady stance phases). However, the stiffness of the test shoe made the shoe motion more abrupt with additional impacts, while inserting the shoe into the cradles. This led to spikes along  the trajectory, which are visible in Fig. 10 near the cradles (''targets''), particularly for software variant no. 3.
Second, to explore the influence of parameter tuning and IMUs, four different cases of tuning were applied for the test walk shown in Fig. 10: AD1: Tuning for the Analog Devices IMU data to achieve a good horizontal position accuracy with software variant no. 1 (exactly the same tuning approach as for the OpenShoe data). B1: Equivalent tuning for the Bosch IMU data to achieve a good horizontal position accuracy with variant no. 1. AD4: Tuning for the Analog Devices IMU data to achieve a good horizontal position accuracy, now with variant no. 4. B4: Tuning for the Bosch IMU data to achieve a good horizontal position accuracy, also with variant no. 4.
Third, the tuned initial values of the covariance matrix P for the Kalman filter were, in most cases, slightly higher than those of OpenShoe. However, the effect of these parameters on the estimation results was rather weak. The OpenShoe values for the three diagonal elements of R for the ZUPT aiding were varied in a range of 100 times higher to 100 times lower to obtain the horizontal position results as accurately as possible. The influence of the R-element values on the estimation results was also weak (without a clear regularity). The strongest influence on the estimation results came from the diagonal elements of the system noise matrix Q (or more precisely, from the ratio between the elements of Q and R). These were mainly lower than the OpenShoe values and possibly reflected the sharper narrowing of the ZUPT phases.
With each of these four sets of parameters, all seven software variants were used again to estimate the IMU trajectory, that is, the parameters were kept constant from variant to variant, as was performed with the OpenShoe data. The tuning of the additional parameters of variants no. 5 and no. 6 was also performed in all four cases in the same way as for the OpenShoe data.
In each of the software runs (i.e., four tuning approaches with seven software variants), the position error was checked not only at the end of the walk but, as aforementioned, every time the cradle was reached (see Fig. 9). Six of these events (with index j) took place during the walk, and each of them led to a distance error in the estimated position relative to the true position, as indicated by the circles in Fig. 10 (compare also Fig. 6). The position error was then computed as the RMS value of the six distance errors: position error = 1 6 6 j=1 ((r x,j − r x,cradle,j ) 2 +(r y,j − y,cradle,j ) 2 ). (20) For the height and yaw errors, the RMS values of all six events were analogously calculated from the estimates at the cradles. The formula for the height is: (r z,j − r z,cradle,j ) 2 with r z,cradle,j = 0.  Fig. 11 shows the height trajectories of the first four variants.
To continue by comparing the different software variants, diagrams as shown in Fig. 8 are presented in Fig. 12-14 for the cases AD1, B1, and B4. (To limit the number of figures, AD4 was omitted, as it did not reveal additional findings.) The most remarkable difference of Fig. 12-14 to the Open-Shoe results is the strong increase in the yaw error (bottom diagrams of Fig. 12-14) for software variant no. 2. (For display reasons, the scales of the ordinate axes of Fig. 12-14 are different from the diagrams in Fig. 8.) This effect manifests itself by an increasing turn of the horizontal motion pattern with accordingly impaired accuracy (top diagrams of Fig. 12-14). As described by Jiménez et al. [12], the turn of the horizontal motion pattern is well known for ZUPT-based pedestrian navigation for data processing according to no. 2 (compare no. 2 in Fig. 5). Without extra aiding, this effect goes completely back to the error level of no. 1 when using only the Runge-Kutta method in combination with the totalstates approach. Partially, ZARU can provide additional small improvements.
In contrast to the horizontal accuracy, which behaves similarly to the yaw error, the height shows only moderate error growth from variant no. 1 to no. 2, if any. The transition to no. 3 does not always lead to improvements, but after the transition to no. 4, there is consistently a strong error reduction with an accuracy that is much better than that of no. 1.
Furthermore, the Bosch IMU (a low-cost device) led to the same horizontal accuracy level as the OpenShoe data, whereas the vertical accuracy in Fig. 12-14 is even one order of magnitude better -a statement that is, however, only valid for variants no. 4 to no. 7. A significant difference between cases B1 and B4 for these variants was not observed. Therefore, it seems of minor importance if software variants no. 1 or no. 4 are used for filter tuning.
Unexpectedly, the more precise IMU from Analog Devices did not lead to an increase in accuracy compared with the Bosch IMU. The vertical accuracy in Fig. 12-14 is one order of magnitude worse. On the other hand, the IMU reacts more sensitively to the harsh walking style owing to the higher bandwidth of its vertical accelerometer (600 Hz [65] compared with 353 Hz of the Bosch sensor [64]) and higher internal sample rates. In addition, Kronenwett [66] observed in a similar way a rather low influence of IMU quality on the results of ZUPT-based pedestrian navigation. He compared two IMUs with comparable accuracy to the Bosch and Analog Devices units, but employed additional aiding techniques.
The uncertain benefit of aiding by ZARU, as reflected in the literature review, remains uncertain given the results for variant no. 5 above. For the Bosch IMU, ZARU turned out to be helpful; for the data of the Analog Devices IMU and the OpenShoe data, ZARU did not show an advantage. Apparently, the statement of Wang and Shkel [14] (p. 68) about ZARU is also valid against an improved numerical background: ZARU requires phases in which the IMU particularly ''stands still''.
Variant no. 6 for the Gauss-Markov bias model and no. 7 for the Runge-Kutta algorithm of the 6 th order delivered neither significant accuracy improvements nor clear VOLUME 10, 2022 deterioration. Therefore, these results firstly reflect the outcome of the literature review that a significant advantage of the Gauss-Markov bias model is not obvious. Second, the results show that the high computational cost of the alternative Runge-Kutta algorithm does not pay off.

D. FURTHER ASPECTS OF ACCURACY
To complete the discussion of Figs. 8, 12, 13, and 14, the horizontal accuracy, which is of particular interest for pedestrian navigation, is looked at even more closely.
As a first point, variants no. 4 to no. 6 match or exceed the performance of the classical OpenShoe approach, which still reflects the state of the art in pure ZUPT-based pedestrian navigation. Therefore, a gain in accuracy could be observed for the alternative algorithmic basics described in Subsection II.A to II.D.
A comparison with results of other research projects that also employ pure ZUPT-aiding is useful as well. To this, the variants no. 4 to no 6 of the cases B1 and B4 are used, which reached an absolute horizontal accuracy of the order of 0.08 m corresponding to 0.3% relative to the distance covered. Bagchi et al. [41] achieved an absolute accuracy of 0.05 m to 0.1 m for the conditions of production engineering with shorter distances covered, whereas Wang et al. [34] report a relative accuracy of about 0.4% for a distance travelled being ten times longer than that of Fig. 10.
The last two papers confirm, as a second point, an adequate accuracy level reached by the authors. Comparisons beyond this with results of research projects using more than ZUPT and ZARU aiding are less certain because added navigation system elements typically alter accuracy. Nevertheless, an attempt will be made to draw a more general picture. The comparisons are based on rather randomly selected recent publications, which however (like the last references [34]  and [41]) largely satisfy the minimum IMU sampling rate of about 250 Hz being recommended by Munoz Diaz et al. [59].
Moderate changes to the navigation system design, as a third point, do not seem to significantly change the accuracy level as long as ZUPT aiding is of central importance. Ma et al. [67] used an advanced detection of ZUPT intervals as well as an additional height aiding by a pressure sensor and obtained a relative accuracy of also about 0.3%. Zhang et al. [68] employed an additional position aiding by an UWB receiver and reported an absolute accuracy around 0.35 m. The fourth point considered is the transition to other types of pedestrian navigation systems that still include inertial sensors. Interestingly, this change also shows no strong change in the accuracy level: Qiu et al. [4] achieved a relative accuracy of 0.3% as well with a wearable sensor network. Experiments of Kronenwett [66] with a hand worn integrated system consisting of an IMU, a sophisticated aiding by a camera, and a pressure sensor yielded absolute accuracies between 0.05 m and 0.09 m. The last point is the accuracy level of pedestrian navigation techniques that do not employ inertial sensors. Pure Map matching is such an approach. Using a combination of indoor floor plans and WLAN fingerprints, Yu et al. [69] reported an absolute accuracy of 1.24 m. Employing instead magnetic field patterns as a map, Ashraf et al. [70] yielded absolute accuracies of about 0.5 m to 1.0 m. Another exemplary technique without inertial sensors is the traditional use of radio signals. Josephy et al. [71] were able to reach absolute VOLUME 10, 2022 accuracies of 13 m to 54 m by combining GNSS with ADS-B messages from aircraft.
On this last point, El-Sheimy and Li [8] note that the techniques and accuracies of non-inertial indoor navigation are diverse. Therefore, additional data are required for a meaningful comparison with pure ZUPT-based pedestrian navigation. However, the other four points show that the accuracy achieved in the authors' experiments is state of the art even when compared to inertial pedestrian navigation systems with additional components.
Two other aspects of the relevant literature on the accuracy of ZUPT-based pedestrian navigation should also be mentioned. The first is the presumption of Munoz Diaz et al. [72], who attribute similar accuracy to the total-state approach for pedestrian navigation as to the error-state approach. However, the above experiments show that the total-state Kalman filter performs better than the error-state filter.
Second, a continuous-discrete Kalman filter, as used above, together with an accurate numerical integration, can show the same performance as an unscented Kalman filter [73], which is repeatedly used in integrated navigation as a (supposedly) more accurate alternative (e.g. [74]). Based on variant no. 4, a comparison of an extended Kalman filter with an unscented one could therefore be an interesting study.

V. CONCLUSION AND RECOMMENDATIONS
The accuracy comparisons at the end of the last section showed that the results of the paper at hand are in line with the usual accuracy level of pedestrian inertial navigation with MEMS IMUs. The direct comparison with the results of the OpenShoe software proved furthermore that the alternatives to the classical algorithmic basics as described in Subsection II.A through II.D offer a potential of significant numerical improvements. Against this background, the following conclusions can be drawn.
The introduction of the classical algorithmic basics by Foxlin [11] and their promotion by projects like OpenShoe [3] created the conditions for an extensive and broad research work on pedestrian inertial navigation with ZUPT aiding. This exceptional, important success, however, neglected the circumstance that Foxlin's formulation of ZUPT-based pedestrian navigation is only one possible realization of strapdown inertial technology -an implementation being characterized by simplicity but not by exploiting the numerical accuracy.
The second point of the conclusions concerns the effect of omitting the inertial sensor biases from the motion state x (Subsection II.B). This proposal by Skog et al. [13] has also significantly supported the research on ZUPT-based pedestrian navigation as it reduced the system instability caused by the incomplete observability of x. Using total states instead of error states for the Kalman filter in combination with a Runge-Kutta method of sufficient order (Subsections II.A and II.D) realizes, however, at least the same stabilizing effect but keeps additionally the compensation of the inertial sensor biases. This shows that the instability being caused by the missing observability is mainly ''fueled'' by avoidable numerical inaccuracies and not by the sensor biases.
Third, the better performance of the Runge-Kutta method compared with the usual Euler method is natural and needs no additional explanation. Instead, the current prevalence of the Euler method is peculiar and little appropriate. (The high popularity of the discrete-discrete Kalman filter (Subsection II.C) may be the reason for this effect.) Based on the results in Section IV, the statement seems therefore necessary that the low accuracy of MEMS IMUs (so called ''consumer grade'' and ''industrial grade'' IMUs [14]) requires nevertheless a signal processing with numerical procedures of high quality, which should be limited only by the computational effort.
The fourth point relates to the better performance of the total state Kalman filter compared with the error state Kalman filter (Subsection II.A). This effect indicates that the linear error models (4), (5) of the filter are not adequate for the high error level of the MEMS IMUs. There is a long tradition of using linear error models for the more accurate ''navigation grade'' and ''tactical grade'' IMUs [14], but MEMS IMUs require a nonlinear modeling in combination with adequate computer numerics.
Fifth, a careful parameter tuning of the ZUPT detector and the data fusion filter is still necessary to individually consider the walking style, footwear, and IMU characteristics. The differences between the software variants no. 1 and no. 4 do not seem to have a primary influence on the tuning. The improved navigation system accuracy provided by the alternative algorithmic basics require, as a sixth point of the conclusions, to reconsider also the performance of other system design features. This has been done above for the aiding by ZARU and the Gauss-Markov modelling of the inertial sensor biases. In both cases, the results did not change the assessment of the current literature: ZARU is applicable only in pronounced resting phases of the IMU, and a significant advantage of the Gauss-Markov modelling is not visible. Further research should address the IMU sampling rate and IMU calibration (as already done by the authors [75]) the usage of other data fusion techniques like an unscented Kalman filter or like smoothing, advanced detection of ZUPT phases, dual IMU configurations etc.
As a last point, the study presented above has proven again that ZUPT is an effective aiding technique. Nevertheless, the application of ZUPT for pedestrian navigation is not very reliable because of the difficult realization and detection of real IMU resting phases. Accordingly, Wahlström and Skog [10] suggest to use ZUPT aiding only to the minimum extent. Besides additional aiding techniques for complementing ZUPT, the improved numerics described above supports this recommendation.
Based on the conclusions above, other recommendations are: 1) The inertial sensor biases should be included in the system error-state. However, this step should not be combined with the Euler method according to (12) and (13). 2) Equations (12) and (13) should be replaced with (4) and (7). Solving these differential equations requires a Runge-Kutta algorithm of sufficiently high order. 3) For further numerical improvements, it is suggested to use a total-state Kalman filter according to Fig. 2 instead of an error-state Kalman filter. 4) Depending on the walking style, footwear, IMU, and aiding methods, aiding by ZARU can be useful, but requires careful tuning and testing. 5) Recommendations 1 through 3 and an IMU sampling rate greater than 250 Hz should form the usual basis for future research on ZUPT-based pedestrian navigation. The authors hope that these recommendations will advance the state of the art in pedestrian and indoor inertial navigation.

APPENDIX: SYSTEM OBSERVABILITY
In the following, the system observability is approximately checked during a single stance phase (outside the stance phases, there is no system observability at all, because h does not exist for pure ZUPT and ZARU aiding).
For the analysis, the coordinate systems n and b (Fig. 3) are used again. Furthermore, it is approximately assumed again that system n represents an inertial frame. In this case, the matrix F of (4) with x of (1b) and the replacement of δq by the error angle vector α is [14]: The additional symbols are the gravity g, transformation matrix nb T from the system b to to the system n, unit matrix I of dimension 3 × 3, and zero matrix 0. Furthermore, for pedestrian navigation, the pitch and roll angles can be assumed to be small during the stance phase. In addition, the yaw angle can be assumed also to be small, as in the navigation frame there is no distinct horizontal direction, so North can be taken as an arbitrary orientation and can be set as required during the stance phase. Under these conditions, nb T becomes approximately equal to I.
Next, matrix H for the combined aiding of ZUPT and ZARU follows from (17) Using these matrices, we get from (6) Considering the different lines of B, the following conclusions can be drawn. Prof. Wagner is a member of the Institute of Navigation ION, the German Institute of Navigation DGON, the German Society for Aeronautics and Astronautics DGLR, and the Association of Applied Mathematics and Mechanics GAMM.
MICHAEL KOHL received the B.Sc. and M.Sc. degrees in aerospace engineering from the University of Stuttgart, Germany, in 2020, with focus on mathematical and physical as well as drive and energy systems.
In 2019, he has completed the semester abroad in Samara and Russia from the Samara State University of Aerospace Engineering. Since 2021, he has been working as a Research Assistant at the Chair of Flight Measuring Technology, University of Stuttgart. His research interests include structural and system dynamics and inertial and integrated measurement systems.
BENEDIKT GYÖRFI received the B.Eng. degree in mechanical engineering from the University of Applied Science Pforzheim, Germany, in 2015, and the M.S. degree in product development from the Faculty of Mechanical Engineering, University of Applied Science Pforzheim.
In 2015, he was a Research Assistant with the research project BikeSafe at the University of Applied Science Pforzheim. He joined the Chair of Flight Measuring Technology, University of Stuttgart, Germany, as a Test Engineer, in 2018. VOLUME 10, 2022