Cooperative localization by dual foot-mounted inertial sensors and inter-agent ranging

The implementation challenges of cooperative localization by dual foot-mounted inertial sensors and inter-agent ranging are discussed, and work on the subject is reviewed. System architecture and sensor fusion are identified as key challenges. A partially decentralized system architecture based on step-wise inertial navigation and step-wise dead reckoning is presented. This architecture is argued to reduce the computational cost and required communication bandwidth by around two orders of magnitude while only giving negligible information loss in comparison with a naive centralized implementation. This makes a joint global state estimation feasible for up to a platoon-sized group of agents. Furthermore, robust and low-cost sensor fusion for the considered setup, based on state space transformation and marginalization, is presented. The transformation and marginalization are used to give the necessary flexibility for presented sampling-based updates for the inter-agent ranging and ranging free fusion of the two feet of an individual agent. Finally, the characteristics of the suggested implementation are demonstrated with simulations and a real-time system implementation.


I. INTRODUCTION
High accuracy, robust, and infrastructure-free pedestrian localization is a highly desired ability for, among others, military, security personnel, and first responders.Localization together with communication are key capabilities to achieve situational awareness and to support, manage, and automatize individual's or agent group actions and interactions.See [1][2][3][4][5][6][7][8] for reviews on the subject.The fundamental information sources for the localization are proprioception, exteroception, and motion models.Without infrastructure, the exteroception must be dependent on prior or acquired knowledge about the environment [9].Unfortunately, in general, little or no prior knowledge of the environment is available and exploiting acquired knowledge without revisiting locations is difficult.Therefore, preferably the localization should primarily rely on proprioception and motion models.Proprioception can take place on the agent level, providing the possibility to perform dead reckoning; or on inter-agent level, providing the means to perform cooperative localization.Pedestrian dead reckoning can be implemented in a number of different ways [10].However, foot-mounted inertial navigation, with motion models providing zero-velocity updates, constitute a unique, robust, and high accuracy pedestrian dead reckoning capability [11][12][13][14].With open-source implementations [15][16][17] and several Fig. 1: Illustration of the considered localization setup.A group of agents are cooperatively localizing themselves without any infrastructure.For this purpose, each agent is equipped with dual foot-mounted inertial sensors, a ranging device, and a communication (com.) and processing (proc.)device.
products appearing on the market [18][19][20][21], dead reckoning by foot-mounted inertial sensors is a readily available technology.In turn, the most straight-forward and well studied interagent measurement, and mean of cooperative localization, is ranging [22][23][24][25].Also here, there are multiple (radio) ranging implementations available in the research literature [26][27][28][29][30] and as products on the market [31][32][33].Finally, suitable infrastructure-free communication equipment for inter-agent communication is available off-the-shelf, e.g.[34][35][36][37], and processing platforms are available in abundance.Together this suggests that the setup with foot-mounted inertial sensors and inter-agent ranging as illustrated in Fig. 1 is suitably used as a base setup for any infrastructure-free pedestrian localization system.However, despite the mature components and the in principle straight-forward combination, cooperative localization with this sensor setup remains challenging, and only a few similar systems can be found in the literature [38][39][40][41][42] and no off-the-shelf products are available.
The challenges with the localization setup lie in the system architecture and the sensor fusion.The inter-agent ranging and lack of anchor nodes mean that some global state estimation is required with a potentially prohibitive large computational cost.The global state estimation, the distributed measurements, and the (required) high sampling rates of the inertial sensors mean that a potentially substantial communication is needed and that the system may be sensitive to incomplete or varying connectivity.The feet are poor placements for inter-agent ranging devices and preferably inertial sensors are used on both feet meaning that the sensors of an individual agent will not be collocated.This gives a high system state dimensionality and means that relating the sensory data from different sensors to each other is difficult and that local communication links on each agent are needed.Further, interagent ranging errors as well as sensor separations, often have far from white, stationary, and Gaussian characteristics.Together, this makes fusing ranging and dead reckoning in a high integrity and recursive Bayesian manner at a reasonable computational cost difficult.
Unfortunately, the mentioned challenges are inherent to the system setup.Therefore, they have to be addressed for any practical implementation.However, to our knowledge, the implementation issues have only been sparsely covered in isolation in the literature and no complete satisfactory solution has been presented.Therefore, in this article we present solutions to key challenges to the system setup and a complete localization system implementation.More specifically, the considered overall problem is tracking, i.e recursively estimating, the positions of a group of agents with the equipment setup of Fig. 1.The available measurements for the tracking are inertial measurements from the dual foot-mounted inertial sensors and inter-agent range measurements.The position tracking is illustrated in Fig. 2. The measurements will be treated as localized to the respective sensors and the necessary communication will be handled as an integral part of the overall problem.However, we will not consider specific communication technologies, but only communication constraints that naturally arise in the current scenario (low bandwidth and varying connectivity).See [43][44][45][46][47] and references therein for treatment of related networking and communication technologies.Also, for brevity, the issues of initialization and time synchronization will not be considered.See [48,49] for the solutions used in the system implementation.
To arrive at the key challenges and the solutions, initially, in Section II, the implementation challenges are discussed in more detail and related work is reviewed.Following this, we address the key challenges and present a cooperative localization system implementation based on dual foot-mounted inertial sensors and inter-agent ranging.The implementation is based on a partially decentralized system architecture and statistical marginalization and sampling based measurement updates.In Section III, the architecture is presented and argued to reduce the computational cost and required communication by around two orders of magnitude and to make the system robust to varying connectivity, while only giving negligible information loss.Thereafter, in Section IV, the sampling based measurement updates with required state space transformation and marginalization are presented and shown to give a robust and low computational cost sensor fusion.Subsequently, in Section V, the characteristic of the suggested implementation is illustrated via simulations and a real-time system implementation.The cooperative localization is found to give a bounded relative position mean-square-error (MSE) and an absolute position MSE inversely proportional to the number of agents, in the worst case scenario; and a bounded position MSE, in the best case scenario.Finally, Section VI concludes the article.

II. IMPLEMENTATION CHALLENGES
The lack of anchor nodes, the distributed nature of the system, the error characteristics of the different sensors, and the non-collocated sensors of individual agents poses a number of implementation challenges for the cooperative localization.Broadly speaking, these challenges can be divided into those related to: designing an overall system architecture to minimize the required communication and computational cost, while making it robust to varying connectivity and retaining sufficient information about the coupling of non-collocated parts of the system; and fusing the information of different parts of the system given the constraints imposed by the system architecture and a finite computational power, while retaining performance and a high system integrity.In the following two subsections, these two overall challenges are discussed in more detail and related previous work is reviewed.

A. System architecture and state estimation
The system architecture has a strong connection to the position/state estimation and the required communication.The range of potential system architectures and estimation solutions goes from the completely decentralized, in which each agent only estimates its own states, to the completely centralized, in which all states of all agents are estimated jointly.
A completely decentralized architecture is often used in combination with some inherently decentralized beliefpropagation estimation techniques [38,50,51].The advantage of this is that it makes the localization scalable and robust to varying and incomplete connectivity between the agents.Unfortunately, the belief-propagation discards information about the coupling between agents, leading to reduced performance [51][52][53][54].See [52] for an explicit treatment of the subject.Unfortunately, as will be shown in Section V, in a system with dead reckoning, inter-agent ranging, and no anchor nodes, the errors in the position estimates of different agents may become almost perfectly correlated.Consequently, discarding these couplings/correlations between agents can significantly deteriorate the localization performance and integrity.
In contrast, with a centralized architecture and estimation, all correlations can be considered, but instead the state dimensionality of all agents will add up.Unfortunately, due to the lack of collocation of the sensors of individual agents, the state dimensionality of individual agents will be high.Together this means computationally expensive filter updates.Further, the distributed nature of the system means that information needs to be gathered to perform the sensor fusion.Therefore, communication links are needed, both locally on each agent as well as on a group level.Inter-agent communication links are naturally wireless.However, the foot-mounting of the inertial sensors makes cabled connections impractical opting for battery powering and local wireless links for the sensors as well [55,56].Unfortunately, the expensive filter updates, the wireless communication links, and the battery powering combines poorly with the required high sampling rates of the inertial sensors.With increasing number of agents, the computational cost and the required communication bandwidth will eventually become a problem.Moreover, an agent which loses contact with the fusion center cannot, unless state statistics are continually provided, easily carry on the estimation of its own states by itself.Also, to recover from an outage when the contact is restored, a significant amount of data would have to be stored, transferred and processed.
Obviously, neither of the extreme cases, the completely decentralized nor the completely centralized architectures, are acceptable.The related problems suggest that some degree of decentralization of the estimation is required to cope with the state dimensionality and communication problems.However, some global book keeping is also required to handle the information coupling.Multiple approximative and exact distributed implementations of global state estimation have been demonstrated, see [54,[57][58][59] and references therein.However, these methods suffer from either a high computational cost or required guaranteed and high bandwidth communication, and are not adapted to the considered sensor setup with high update rates, local communication links, and lack of sensor collocation.Therefore, in Section III we suggest and motivate a system architecture with partially decentralized estimation based on a division of the foot-mounted inertial navigation into a step-wise inertial navigation and dead reckoning.This architecture does not completely solve the computational cost issue, but makes it manageable for up to a platoon-sized group of agents.For larger groups, some cellular structure is needed [39,58].However, the architecture is largely independent of how the global state estimation is implemented and a distributed implementation is conceivable.
The idea of dividing the filtering is not completely new.A similar division is presented in an application specific context in [60] and used to fuse data from foot-mounted inertial sensors with maps, or to build the maps themselves, in [61][62][63].However, the described division is heuristically motivated and the statistical relation between the different parts is not clear.Also, no physical processing decentralization is exploited to give reduced communication requirements.

B. Robust and low computational cost sensor fusion
The sensor fusion firstly poses the problem of how to model the relation between the tracked inertial sensors and the range measurements.Secondly, it poses the problem of how to condition the state statistic estimates on provided information while retaining reasonable computational costs.
The easiest solution to the non-collocated sensors of individual agents is to make the assumption that they are collocated (or have a fixed relation) [38,[64][65][66].While simple, this method can clearly introduce modeling errors resulting in suboptimal performance and questionable integrity.Instead, explicitly modeling the relative motion of the feet has been suggested in [67].However, making an accurate and general model of the human motion is difficult, to say the least.As an alternative, multiple publications suggest explicitly measuring the relation between the sensors [14,[68][69][70].The added information can improve the localization performance but unfortunately introduces the need for additional hardware and measurement models.Also, it works best for situations with line-of-sight between measurement points, and therefore, it is probably only a viable solution for foot-to-foot ranging on clear, not too rough, and vegetation/obstacle free ground [71].Instead of modeling or measuring the relation between navigation points of an individual agent, the constraint that the spatial separation between them has an upper limit may be used.This side information obviously has an almost perfect integrity, and results in [72] indicate that the performance loss in comparison to ranging is transitory.For inertial navigation, it has been demonstrated that a range constraint can be used to fuse the information from two foot-mounted systems, while only propagating the mean and the covariance [73,74].Unfortunately, the suggested methods depend on numerical solvers and only apply the constraint on the mean, giving questionable statistical properties.Therefore, in Section IV, based on the work in [72], we suggest a simpler and numerically more attractive solution to using range constraints to perform the sensor fusion, based on marginalization and sampling.
The naive solution to the sensor fusion of the footmounted inertial navigation and the inter-agent ranging is simply using traditional Kalman filter measurement updates for the ranging [38].However, the radio ranging errors are often far from Gaussian, often with heavy tails and nonstationary and spatially correlated errors [75][76][77][78][79][80].This can cause unexpected behavior of many localization algorithms, and therefore, statistically more robust methods are desirable [79][80][81].See [82] and references therein for a general treatment of the statistical robustness concept.The heavy tails and spatially correlated errors could potentially be solved by a flat likelihood function as suggested in [75,83].However, while giving a high integrity, this also ignores a substantial amount of information and requires multi-hypothesis filtering (a particle filter) with unacceptable high computational cost.Using a more informative likelihood function is not hard to imagine.Unfortunately, only a small set of likelihood functions can easily be used without resorting to multi-hypothesis filtering methods.Some low cost fusion techniques for special classes of heavy-tailed distributions and H ∞ criteria have been suggested in the literature [84][85][86][87][88].However, ideally we would like more flexibility to handle measurement errors and noncollocated sensors.Therefore, in Section IV, we propose a marginalization and sample based measurement update for the inter-agent ranging, providing the necessary flexibility to handle an arbitrary likelihood function.A suitable likelihood function is proposed, taking lack of collocation, statistical robustness, and correlated errors into account, and shown to providing a robust and low computational cost sensor fusion.

III. DECENTRALIZED ESTIMATION ARCHITECTURE
To get around the problems of the centralized architecture, the state estimation needs somehow to be partially decentralized.However, as previously argued, some global state estimation is necessary.Consequently, the challenge is to do the decentralization in a way that does not lead to unacceptable loss in information coupling, leading to poor performance and integrity, while still solving the issues with computational cost, communication bandwidth, and robustness to varying connectivity.In the following subsections it is shown how this can be achieved by dividing the filtering associated with footmounted inertial sensors into a step-wise inertial navigation and step-wise dead reckoning, as illustrated in Fig. 3. Pseudocode for the related processing is found in Alg. 1 and the complete system architecture is illustrated in Fig. 5.

A. Zero-velocity-update-aided inertial navigation
To track the position of an agent equipped with footmounted inertial sensors, the sensors are used to implement an inertial navigation system aided by so called zero-velocity updates (ZUPTs).The inertial navigation essentially consists of the inertial sensors combined with mechanization equations.In the simplest form, the mechanization equations are where k is a time index, dt is the time difference between measurement instances, p k is the position, v k is the velocity, f k is the specific force, g = [0, 0, g] is the gravity, and ω k is the angular rate (all in R 3 ).Further, q k is the quaternion describing the orientation of the system, the triple product q k−1 f k q k−1 denotes the rotation of f k by q k , and Ω(•) is the quaternion update matrix.For a detailed treatment of inertial navigation see [89,90].For analytical convenience we will interchangeably represent the orientation q k with the equivalent Euler angles (roll,pitch,yaw) . .] is used to denote a column vector.The mechanization equations (1) together with measurements of the specific force fk and the angular rates ωk , provided by the inertial sensors, are used to propagate position pk , velocity vk , and orientation qk state estimates.Unfortunately, due to its integrative nature, small measurement errors in fk and ωk accumulate, giving rapidly growing estimation errors.Fortunately, these errors can be modeled and estimated with ZUPTs.A first-order error model of ( 1) is given by  where δ(•) k are the error states, I and 0 are 3 × 3 identity and zero matrices, respectively, and [•] × is the crossproduct matrix.As argued in [91], one should be cautious about estimating systematic sensor errors in the current setup.Indeed, remarkable dead-reckoning performance has been demonstrated, exploiting dual foot-mounted sensors without any sensor error state estimation [92].Therefore, in contrast to many publications, no additional sensor bias states are used.
Together with statistical models for the errors in fk and ωk , ( 2) is used to propagate statistics of the error states.To estimate the error states, stationary time instances are detected based on the condition Z({ fκ , ωκ } W k ) < γ Z , where Z(•) is some zero-velocity test statistic, { fκ , ωκ } W k is the inertial measurements over some time window W k , and γ Z is a zerovelocity detection threshold.See [93,94] for further details.The implied zero-velocities are used as pseudo-measurements which are modeled in terms of the error states as where H = [0 I 0] is the measurement matrix and n k is a measurement noise, i.e. ỹk = δv k + n k .A similar detector is also used to lock the system when completely stationary.See [95] for further details.Given the error model ( 2) and the measurements model ( 4), the measurements (3) can be used to estimate the error states with a Kalman type of filter.See [11,94,96,97] for further details and variations.See [98] for a general treatment of aided navigation.Since there is no reason to propagate errors, as soon as there are any non-zero estimates δ pk , δv k , or δ θk , they are fed back correcting the navigational states, pk vk := pk vk + δ pk δv k and qk := Ω(δ θk )q k (5) and consequently the error state estimates are set to zero, i.e. δ pk := 0 3×1 , δv k := 0 3×1 , and δ θk := 0 3×1 , where := indicates assignment.Unfortunately, all (error) states are not observable based on the ZUPTs.During periods of consecutive ZUPTs, the system (2) becomes essentially linear and time invariant.Zero-velocity for consecutive time instances means no acceleration and ideally f k = q k gq k .This gives the system and observability matrices Obviously, the position (error) is not observable, while the velocity is.Since the roll and pitch are observable while the heading (yaw) of the system is not.Ignoring the process noise, this implies that the covariances of the observable states decay as one over the number of consecutive ZUPTs.Note that there is no difference between the covariances of the error states and the states themselves.Consequently, during stand-still, after a reasonable number of ZUPTs, the state estimate covariance becomes where P x,y = cov(x,y), P x = cov(x) = cov(x,x), (•) denotes the transpose, and 0 n×m denotes a zero matrix of size n × m.

B. Step-wise dead reckoning
The covariance matrix (6) tells us that the errors of pk and ψk are uncorrelated with those of vk and [ φk , θk ].Together with the Markovian assumption of the state space models and the translational and in-plan rotation invariance of ( 1)-( 4), this means that future errors of vk and [ φk , θk ] are uncorrelated with those of the current pk and ψk .Consequently, future ZUPTs cannot be used to deduce information about the current position and heading errors.In turn, this means that, considering only the ZUPTs, it makes no difference if we reset the system and add the new relative position and heading to those before the reset.However, for other information sources, we must keep track of the global (total) error covariance of the position and heading estimates.
Resetting the system means setting position pk and heading ψk and corresponding covariances to zero.Denote the position and heading estimates at a reset by dp and dψ .These values can be used to drive the step-wise dead reckoning where x and χ are the global position in 3 dimensions and heading in the horizontal plan of the inertial navigation system relative to the navigation frame, is the rotation matrix from the local coordinate frame of the last reset to the navigation frame, and w is a (by assumption) white noise with covariance, The noise w in (7) represents the accumulated uncertainty in position and heading since the last reset, i.e. the essentially non-zero elements in (6) transformed to the navigation frame.The dead reckoning (7) can trivially be used to estimate x and χ and their error covariances from dp and dψ and related covariances.The relation between the step-wise inertial navigation and dead reckoning is illustrated in Fig. 3.
Step-wise dead reckoning Step-wise inertial navigation Fig. 3: Illustration of the step-wise inertial navigation and the step-wise dead reckoning.The displacement and heading change over a step given by the inertial navigation is used to perform the step-wise dead-reckoning.
To get [dp , dψ ] from the inertial navigation, reset instances need to be determined, i.e. the decoupled situation (6) needs to be detected.However, detecting it is not enough.If it holds for one time instance k, it is likely to hold for the next time instance.Resetting at nearby time instances is not desirable.Instead we want to reset once at every step or at some regular intervals if the system is stationary for a longer period of time.The latter requirement is necessary to distinguish between extended stationary periods and extended dynamic periods.Further, to allow for real-time processing, the detection needs to be done in a recursive manner.The longer the stationary period, the smaller the cross-coupling terms in (6).This means that the system should be reset as late as possible in a stationary period.However, if the stationary period is too short, we may not want to reset at all, since then the cross terms in (6) may not have converged.In summary, necessary conditions for a reset are low enough cross-coupling and minimum elapsed time since the last reset.If this holds, there is a pending reset.In principle, the cross-coupling terms in (6) should be used to determine the first requirement.However, in practice, all elements fall off together and a threshold γ p on e.g. the first velocity component can be used.To assess the second requirement, a counter c p which is incremented at each time instance is needed, giving the pending reset condition where c min is the minimum number of samples between resets.A pending reset is to be performed if the stationary period comes to an end or a maximum time with a pending reset has elapsed.To assess the latter condition, a counter c d is needed which is incremented if (9) holds.Then a reset is performed if where c max is the maximum number of samples of a pending reset.Together ( 9) and ( 10) make up the sufficient conditions for a reset.When the reset is performed, the counters are reset, c p := 0 and c d := 0. This gives a recursive step segmentation.Pseudo-code for the inertial navigation with recursive step segmentation (i.e.step-wise inertial navigation) and the stepwise dead reckoning is found in Alg. 1.Not to lose performance in comparison with a sensor fusion approach based on centralized estimation, the step-wise inertial navigation combined with the step-wise dead reckoning needs Algorithm 1 Pseudo-code of the combined step-wise inertial navigation and step-wise dead reckoning.The ZUPT-aided inertial navigation and the step-wise dead reckoning refers to the effect of (1)-( 5) and (7), respectively, combined with Kalman type of filtering.For notational compactness, below P k = P [p k ,v k ,q k ] and P = P ZUPT-aided inertial navigation Step-wise dead reckoning end if 22: end loop to reproduce the same state statistics (mean and covariance) as those of the indefinite (no resets) ZUPT-aided inertial navigation.If the models (1), (2), and (7) had been linear with Gaussian noise and the cross-coupling terms of (6) were perfectly zero, then the divided filtering would reproduce the full filter behavior perfectly.Unfortunately, they are not.However, as shown in the example trajectory in Fig. 4, in practice the differences are marginal and the mean and covariance estimates of the position and heading can be reproduced by only [dp , dψ ] and the corresponding covariances.Due to linearization and modeling errors of the ZUPTs, the step-wise dead reckoning can even be expected to improve performance since it will eliminate these effects to single steps [91,99].Indeed, resetting appropriate covariance elements (which has similar effects as of performing the step-wise dead reckoning) has empirically been found to improve performance [100].

C. Physical decentralization of state estimation
The step-wise inertial navigation and dead reckoning as described in Alg. 1 can be used to implement a decentralized architecture and state estimation.The ranging, as well as most additional information, is only dependent on position and not on the full state vector [p k , v k , θ k ].Further, as argued in the previous subsection, the errors of vk and [ φk , θk ] are weakly correlated with those of pk and ψk .Therefore, only the states [x , χ ] (for all feet) have to be estimated jointly and only line 19 need to be executed centrally.The step-wise inertial navigation, i.e.Alg. 1 apart from line 19, can be implemented locally in the foot-mounted units, and thereby, only [dp , dψ ] and related covariances need to be transmitted from the feet.This way, the required communication will be significantly lower compared to the case in which all inertial data would have to be transmitted.Also, since the computational cost of propagating ( 7) is marginal, this can be done both locally on the processing device of each agent and in a global state estimation.This way, if an agent loses contact with whomever who performs the global state estimation, it can still perform the dead reckoning, and thereby, keep an estimate of where it is.Since the amount of data in the displacement and heading changes is small, if contact is reestablished, all data can easily be transferred and its states in the global state estimation updated.The other way around, if corrections to the estimates of [x , χ ] are made in the central state estimation, these corrections can be transferred down to the agent.Since the recursion in ( 7) is pure dead reckoning (no statistical conditioning), these corrections can directly be used to correct the local estimates of [x , χ ].This way, the local and the global estimates can be kept consistent.The straight-forward way of implementing the global state estimation is by one (or multiple) central fusion center to which all dead reckoning data are transmitted (potentially by broadcasting).The fusion center may be carried by an agents, or reside in a vehicle or something similar.Range measurements relative to other agents only have a meaning if  the position estimate and its statistics are known.Therefore, all ranging information is transferred to the central fusion center.This processing architecture with its three layers of estimation (foot, processing device of agent, and common fusion center) is illustrated in Fig. 5.However, the division in step-wise inertial navigation and dead reckoning is independent of the structure with a fusion center and some decentralized global state estimation could potentially be used.

D. Computational cost and required communication
The step-wise dead reckoning is primarily motivated and justified by the reduction in computational cost and required communication bandwidth.With a completely centralized sensor fusion [ fk , ωk ], 6 measurement values in total, need to be transferred to the central fusion center at a sampling rate f IMU in the order of hundreds of [Hz] with each measurement value typically consisting of some 12-16 bits.With the step-wise dead reckoning, [dp , dψ ], P p , P p ,ψ , and P ψ ,ψ , in total 14 values, need to be transferred to the central fusion center at a rate of f sw ≈ 1 [Hz] (normal gait frequency [101]).In practice, the 14 values can be reduced to 8 values since crosscovariances may be ignored and the numerical ranges are such that they can reasonably be fitted in some 12-16 bits each.The other way around, some 4 correction values need to be transferred back to the agent.Together this gives the ratio of the required communication of (6 • f IMU )/(12 • f sw ) ≈ 10 2 , a two orders magnitude reduction.
In turn, the computational cost scales linearly with the update rates f IMU and f sw .In addition, the computational cost has a cubic scaling (for covariance based filters) with the state dimensionality.Therefore, the reduction in the computational cost at the central fusion center is at the most f IMU /f sw • ( 9 /4) 3 ≈ 10 3 .However, at higher update rates, updates may be bundle together.Consequently, a somewhat lower reduction may be expected in practice giving a reduction of again around two orders of magnitude.

IV. ROBUST AND LOW-COST SENSOR FUSION
The step-wise dead reckoning provides a low dimensional and low update rate interface to the foot-mounted inertial navigation.With this interface, the global state of the localization system (the system state as conceived by the global state estimation) becomes where x j and χ j are the positions and headings of the agents' feet with dropped time indices.Other auxiliary states may also be put in the state vector.Our desire is to fuse the provided dead-reckoning with that of the other foot and that of other agents via inter-agent ranging.This fusion is primarily challenging because: 1) The high dimensionality of the global system.2) The non-collocated sensors of individual agents.
3) The potentially malign error characteristic of the ranging.The high dimensionality is tackled by only propagating mean and covariance estimates and by marginalization of the state space.The lack of collocation is handled by imposing range constraints between sensors.Finally, the error characteristic of the ranging is handled by sampling based updates.In the following subsections, these approaches are described.Pseudocode for the sensor fusion is found in Algs.2-3 in the final subsection.

A. Marginalization
New information (e.g.range measurements) introduced in the systems is only dependent on a small subset of the states.Assume that the state vector can be decomposed as z = [z 1 , z 2 ], such that some introduced information π is only dependent on z 1 .Then with a Gaussian prior with mean and covariance ẑ = ẑ1 ẑ2 and P z = P z1 P z1z2 the conditional (with respect to π) mean of z 2 and the conditional covariance can be expressed as [72] ẑ2|π = V + U ẑ1|π where and C z1|π is the conditional second order moment of z 1 .Note that this will hold for any information π only dependent on z 1 , not just range constraints as studied in [72].Consequently, the relations (12) give a desired marginalization.To impose the information π on (11), only the first and second conditional moments, ẑ1|π and C z1|π , need to be calculated.If π is linearly dependent on z 1 and with Gaussian errors, this will be equivalent with a Kalman filter measurement update.This may trivially be used to introduce information about individual agents.However, as we will show in the following subsections, this can also be used to couple multiple navigation points of individual agents without any further measurement and to introduce non-Gaussian ranging between agents.

B. Fusing dead reckoning from dual foot-mounted units
The position of the feet x a and x b of an agent (in general two navigation points of an agent) has a bounded spatial separation.This can be used to fuse the related dead reckoning without any further measurements.In practice, the constraint will often have different extents in the vertical and the horizontal direction.This can be expressed as a range constraint where D γ is a diagonal scaling matrix with γ = [1, 1, γ xy /γ z ] on the diagonal, and γ xy and γ z are the constraints in the horizontal and vertical direction.Unfortunately, there is no standard way of imposing such a constraint in a Kalman like framework [102].Also, the position states being in arbitrary locations in the global state vector, i.e. x = [. . ., x a , . . ., x b , . . .], means that the state vector is not on the form of z.Further, since the constraint ( 13 In contrast to (13), this constraint has a bounded support.Therefore, as suggested in [72], the conditional means can be approximated by sampling and projecting sigma points.where ≈ denotes approximate equality and Here s (i) and w (i) are sigma points and weights where l i is the ith column of the Cholesky decomposition P z1 = LL and the scalar η reflects the portion of the prior to which the constraint is applied.See [72] for further details.The application of the constraint for the two-dimensional case is illustrated in Fig. 6.

C. Inter-agent range measurement updates
Similar to the range constraint between the feet of an individual agent, the geometry of the inter-agent ranging gives the constraints where r is the (true) range between agents' ranging devices and γ a and γ b are the maximum spatial separation of respective ranging device and x a and x b ; where in this case x a and x b are the positions of a foot of each agent.The range only being dependent on x a − x b means that the state transformation z = T 1 x where 1 = [1, 1, 1] and z 1 = x a − x b , and the corresponding mean and covariance transformations as explained in the previous subsection, can be used to let us exploit the marginalization (12).The inter-agent ranging gives measurements r of the range r.As reviewed in Section II-B, the malign attributes of r which we have to deal with are potentially heavy-tailed error distributions and non-stationary spatially correlated errors due to diffraction, multi-path, or similar.This can be done by using the model r = r + v + v where v is a heavy-tailed error component and v is a uniformly distributed component intended to cover the assumed bounded correlated errors in a manner similar to that of [75].Combining the model with (16) and the state transformation T 1 gives the measurement model where γ r is chosen to cover the bounds in (16), the asynchrony between xa and xb , and the correlated errors v .In practice γ r will be a system parameter trading integrity for information.
To update the global state estimate with the the range measurement r, the state ẑ1 and covariance estimates P z1 must be conditioned on r via (17).Due to the stochastic term v, we cannot use hard constraints as with the feet of a single agent.However, by assigning a uniform prior to the constraint in (17), the likelihood function of r given ẑ1 becomes where U(−γ r , γ r ) is a uniform distribution over the interval is the distribution of v with mean ẑ1 − r and some scale σ r , and * denotes convolution.Then, with the assumed Gaussian prior z 1 ∼ N (ẑ 1 , P z1 ), the conditional distribution of z 1 given r, ẑ1 , and P z1 is Since z 1 is low-dimensional, the conditional moments ẑ1|r and C z1|r can be evaluated by sampling.With the marginalization (12) and the inverse transformation T −1 1 , this will give the conditional mean and covariance of x.
Since the likelihood function ( 18) is typically heavy-tailed, it cannot easily be described by a set of samples.However, since the prior is (assumed) Gaussian, the sampling of it can efficiently be implemented with the eigenvalue decomposition.With sample points u (i) of the standard Gaussian distribution, the corresponding sample points of the prior is given by where P z1 = QΛQ is the eigenvalue decomposition of P z1 .With the sample points s (i) , the associated prior weights only become dependent on u (i) (apart from normalization) since and can therefore, be precalculated.Reweighting with the likelihood function, ) and normalizing the weights w , with suitable chosen u (i) the conditional moments can be approximated by Consequently, as long as the likelihood function can be efficiently evaluated, any likelihood function may be used.
For analytical convenience, we have typically let V(•, •) be Cauchy-distributed giving the heavy-tailed likelihood function  7: Illustration of the suggested range update in two dimensions.From left to right, the plots show: the Gaussian prior given by the mean ẑ1 (blue dot) and the covariance Pz 1 (blue ellipse) and with indicated samples (red dots) and eigenvectors/-values, the (one dimensional) likelihood function given by the range measurement r and used to reweight the samples, and the resulting posterior with conditional mean ẑ1|r and covariance P z 1 |r calculated from the reweighted samples.The sampling based range update with this likelihood function and u (i) from a square sampling lattice is illustrated in Fig. 7. Potential more elaborate techniques for choosing the sample points can be found in [103].
The presented ranging update gives a robustness to outliers in the measurement data.In Fig. 8, the influence functions for the sample based update and the traditional Kalman measurement update are shown for the ranging likelihood function (20) with γ r = 2 [m] and σ r = 0.5 [m] and position covariance values of P z1 = I [m 2 ] and P z1 = 0.3 I [m 2 ].By comparing the blue solid and the red dashed-dotted lines, it is seen that when the position and ranging error covariances are of the same size, the suggested ranging update behaves like the Kalman update up to around three standard deviations, where it gracefully starts to neglect the range measurement.In addition, by comparing the blue dashed and the red dotted lines, it is seen that for smaller position error covariances, in contrast to the Kalman update, the suggested range update neglects ranging measurements with small errors (flat spot in the middle of the influence function).This has the effect that multiple ranging updates will not make the position error covariance collapse, which captures the fact that due to correlated errors, during standstill, multiple range measurements will contain a diminishing amount of information; and during motion the range measurements should only "herd" the dead reckoning.
With slight modifications, the ranging updates can be used to incorporate information from many other information sources.Ranging to anchor nodes whose positions are not Algorithm 2 Pseudo-code for imposing the range constraint (13), between navigation points xa and x b , on the global state estimate x.
1: ẑ := T γ x and P z := T γ P x T γ 2: L := chol(P z1 ) 3: Sample and project sigma points/weights (z 5: Calculate conditional mean and covariance by marginalization 1: ẑ := T 1 x and P z := T 1 P x T 1 2: (Q, Λ) := eig(P z1 ) 3: s (i) := ẑ1 + QΛ 1/2 u (i) ∀i 4: po s (i) (s (i) ) 7: Calculate conditional mean and covariance by marginalization (ẑ |r , P z|r ) ← (ẑ 1|r , ẑ2 , C z1|r , P z ) 8: x|r := T −1 1 ẑ|r and P x|r := T −1 1 P z|r T − 1 kept in x or position updates (from a GNSS receiver or similar) may trivially be implemented as range updates (zero range in the case of the position update) with z 1 = x a − x b replaced with z 1 = x a − x c where x c is the position of the anchor node or the position measurement.Fusion of pressure measurements may be implemented as range updates in the vertical direction, either relative to other agents or relative to a reference pressure.

D. Summary of sensor fusion
The central sensor fusion, as described in Section III-C, keeps the position and heading of all feet in the global state vector x.From all agents, it receives dead reckoning updates, [dp , dψ ], P p , P p ,ψ , and P ψ ,ψ , and inter-agent range measurements r.The dead reckoning updates are used to propagate the corresponding states and covariances according to (7).At each dead reckoning update, the range constraint is imposed on the state space as described in subsection IV-B, and corrections are sent back to the agent.The inter-agent range measurements are used to condition the state space as described in in subsection IV-C.Pseudo-code for conditioning the state mean and covariance estimates on the range constraint and range measurements is shown in Algs.2-3.

V. EXPERIMENTAL RESULTS
To demonstrate the characteristics of the sensor fusion presented in the previous section, in the following subsection we first show numerical simulations giving a quantitative description of the fusion.Subsequently, to demonstrate the practical feasibility of the suggested architecture and sensor fusion, a real-time localization system implementation is briefly presented.

A. Simulations
The cooperative localization by foot-mounted inertial sensors and inter-agent ranging is nonlinear and the behavior of the system will be highly dependent on the trajectories.Therefore, we cannot give an analytical expression for the performance.Instead, to demonstrate the system characteristics, two extreme scenarios are simulated.For both scenarios the agents move with 1 [m] steps at 1 [Hz].Gaussian errors with standard deviation 0.01 [m] and 0.2[ • ] were added to the step displacements and the heading changes, respectively, and heavy-tailed Cauchy distributed errors of scale 1 [m] where added to the range measurements.The ranging is done timemultiplexed in a Round-robin fashion at a total rate of 1 [Hz].
1) Straight-line march: N agents are marching beside each other in straight lines with agent separation of 10 [m].The straight line is the worst case scenario for the dead reckoning and the position errors will be dominated by errors induced by the heading errors.In Fig. 9, an example of the estimated trajectories of the right (blue) and left (green) feet are shown from 3 agents without any further information, with range constraints between the feet, and with range constraints and inter-agent ranging.The absolute and relative root-meansquare-error (RMSE) as a function of the walked distance, and for different number of agents, are shown in Fig. 10.The relative errors are naturally bounded by the inter agent ranging.However, the heading RMSE grows linearly with time/distance, and therefore, the absolute position error is seen to grow with distance.Similar behavior can be observed in the experimental data in [64,91].Naturally, the heading error, and therefore, also the absolute position RMSE drops as 1/ √ N where N is the number of agents.This is shown in Fig. 11.We may also note that the position errors of different agents become strongly correlated.The correlation coefficients for two agents as a function of distance are shown in Fig. 12.  √ N (dashed red line), the position error is seen to be decaying as the square-root of the number of agents.
2) Three static agents: Three non-collinear agents are standing still.This will be perceived by the foot-mounted inertial navigation, and therefore, they essentially become anchor nodes.This is obviously the best-case scenario.A fourth agent walk around them in a circle.An example of an estimated trajectory is shown in Fig. 13 and the RMSE as a function of time is shown in Fig. 14.Since anchor nodes are essentially present in the system, the errors are bounded.See [104] for further discussions.The non-zero RMSE reflects the range constraints in the system.
From the two scenarios, we can conclude that the relative position errors are kept bounded by the inter-agent ranging while the absolute position errors (relative starting location) are bounded in the best-case (stationary agents); and that the error growth is reduced by a factor of 1/ √ N in the worst case.

B. Real-time implementation
The decentralized system architecture has been realized with OpenShoe units [15] and Android smart-phones and tablets (Samsung Galaxy S III and Tab 2 10.1) in the inhouse developed tactical locator system TOR.The communication is done over an IEEE 802.11WLAN.Synthetic inter-agent ranging is generated from position measurements from a Ubisense system (Ubisense Research & Development Package), installed in the KTH R1 reactor hall [105].The intension is to replace the Ubisense system with in-house developed UWB radios [26].The equipment for a single agent  is shown in Fig. 15.The multi-agent setup with additional equipment for sensor mounting is shown in Fig. 16.
The step-wise inertial navigation and the associated transfer of displacements and heading changes has been implemented in the OpenShoe units.The agent filtering has been implemented as Android applications together with graphical user interfaces.A screen-shot of the graphical user interface with trajectories from a ∼10 [min] search in the reactor hall and adjacent rooms (built up walls not displayed) by 3 smoked divers is shown in Fig. 17.The central sensor fusion has been implemented as a separate Android application running on one agent's Android platform.Recently, voice radio communication and 3D audio has been integrated into the localization system [106].VI.CONCLUSIONS Key implementation challenges of cooperative localization by foot-mounted inertial sensors and inter-agent ranging are: designing an overall system architecture to minimize the required communication and computational cost, while retaining the performance and making it robust to varying connectivity; and fusing the information from the system under the constraint of the system architecture while retaining a high integrity and accuracy.A solution to the former problem has been presented in the partially decentralized system architecture based on the division and physical separation of the stepwise inertial navigation and the step-wise dead reckoning.A solution to the latter problem has been presented in the marginalization and sample based spatial separation constraint and ranging updates.By simulations, it has been shown that in the worst case scenario, the absolute localization RMSE improves as the square-root of the number of agents and the relative errors are bounded.In the best case scenario, both the relative and the absolute errors are bounded.Finally, the feasibility of the suggested architecture and sensor fusion has been demonstrated with simulations and a real-time system implementation featuring multiple agents and a meter level accuracy for operation times of tenth of minutes in a harsh industrial environment.

Fig. 2 :
Fig. 2: Illustration of the localization estimation problem and the desired output from the localization system.The problem is to track (i.e.recursively estimate) the positions of multiple agents in three dimensions by inter-agent range measurements and inertial measurements from the foot-mounted inertial sensors.

3 Fig. 4 :
Fig.4: The plots show the trajectory (upper), the position error covariances (middle), and the covariance between the position and heading errors (lower) as estimated by an extended Kalman filter based indefinite ZUPT-aided inertial navigation (solid lines) and a step-wise inertial navigation and dead reckoning (crossed lines).The agreement between the systems are far below the accuracy and integrity of the former system.
d χ d p , d ψ , r d x , d χ Inter-agent ranging device Foot-mounted IMU with step-wise inertial navigation Local processing device with step-wise dead reckoning.

Fig. 5 :
Fig. 5: Illustration of the decentralized system architecture.Step-wise inertial navigation is done locally in the foot-mounted units.Displacement and heading changes are transferred to a local processing device, where stepwise dead reckoning is performed, and relayed together with ranging data to a central fusion center.The fusion center may be carried by an agents, reside in a vehicle or something similar, or be distributed among agents.
) has unbounded support, the conditional means xa|γ , xb|γ and covariances cov [x a|γ , x b|γ ] cannot easily be evaluated.Moreover, since the states are updated asynchronously (steps occurring at different time instances), the state estimates xa and xb may not refer to the same time instance.The latter problem can be handled by adjusting the constraint γ xy by the time difference of the states.In principle, this means that an upper limit on the speed by which an agent moves is imposed.The former problems can be solved with the state transformation z = T γ x where T γ = D γ −D γ D γ D γ ⊕ I m−6 Π (14) where I m−6 is the identity matrix of size m − 6, m is the dimension of x, ⊕ denotes the direct sum of matrices, Π is a permutation matrix fulfilling Π[. . ., x a , . . ., x b , . . .] = [x a , x b , . . .] and z 1 = D γ (x a − x b ).With the state transformation T γ , the mean and covariance of z become ẑ = T γ x and P z = T γ P x T γ .Inversely, since T γ is invertible, the conditional mean and covariance of x become x|γ = T −1 γ ẑ|γ and P x|γ = T −1 γ P z|γ T − γ .Therefore, if ẑ1|γ and C z1|γ are evaluated, (12) gives ẑ|γ and P z|γ and thereby also x|γ and P x|γ .Fortunately, with z 1 = D γ (x a − x b ), the constraint (13) becomes z 1 ≤ γ xy .

Fig. 6 :
Fig. 6: Illustration of the separation constraint update for two feet in the horizontal plane.The plots show (upper-left) the prior position and covariance estimates, (upper-right) the transformed system with the sigma points (blue crosses) and the constraint indicated with the dashed circle, (lower-left) the projected sigma points and the conditional mean and covariance, and finally (lower-right) the conditioned result in the original domain with the prior covariances indicated with thinner lines.
Fig.7: Illustration of the suggested range update in two dimensions.From left to right, the plots show: the Gaussian prior given by the mean ẑ1 (blue dot) and the covariance Pz 1 (blue ellipse) and with indicated samples (red dots) and eigenvectors/-values, the (one dimensional) likelihood function given by the range measurement r and used to reweight the samples, and the resulting posterior with conditional mean ẑ1|r and covariance P z 1 |r calculated from the reweighted samples.

Fig. 8 :
Fig. 8: Influence function of the ranging update for ẑ1 = 10 [m] performed with the suggested method (blue solid/dashed lines) and a traditional Kalman measurement update (red dashed-dotted/dotted lines) with Pz 1 = I [m 2 ] and Pz 1 = 0.3I [m 2 ], respectively.For the suggested update γr = 2 [m] and v 2 was Cauchy distributed with σr = 0.5 [m] and for the Kalman measurement update the measurement error variance was 1 [m 2 ].

Fig. 9 :
Fig.9: Illustration of the gain of dual foot-mounted sensors and inter-agent ranging.The upper plot shows the step-wise dead reckoning of the individual feet (in blue and green) without any further information.The middle plot shows the step-wise dead reckoning with range constraints between the feet of individual agents.The lower plot shows the complete cooperative localization with step-wise dead reckoning, range constraints, and inter-agent ranging.

Fig. 10 :
Fig. 10: Absolute position RMSE (blue lines) and relative position RMSE (red lines) as a function of distance for 100 Monte-Carlo runs.The different blue lines correspond, in ascending order, to increasing number of agents.Clearly, the relative error is bounded by the inter-agent ranging while the absolute error grows slower the larger the number of agents.The final position RMSEs as a function of the number of agents are shown in Fig. 11.

Fig. 11 :
Fig. 11: Final position RMSE as a function of the number of agents (blue crossed line) for 100 Monte-Carlo runs.From the fit to 1/

Fig. 12 :Fig. 13 :
Fig.12: Position error correlation coefficient for two agents in the straightline march scenario in the x direction (solid blue), y direction (dashed green), and z direction (dotted/dashed red).Clearly, the positions errors of different agents become strongly correlated with increasing distance traveled.

Fig. 14 :
Fig.14:The plot shows the RMSE of the mobile unit for the three static agents scenario over 100 Monter-Carlo runs.Being static, the three stationary agents essentially become anchor nodes, and therefore, the RMSE is bounded.

Fig. 16 :
Fig. 16: Four agent with equipment displayed in Fig. 15.The OpenShoe units are integrated in the soles of the shoes and the radio tags are attached to the helmets.The cables and the USB hubs are not displayed.

Fig. 17 :
Fig.17: Screen-shot of the cooperative localization system user interface overlayed on a photography of the Android phone.The interface shows the recursively estimated positions (trajectories) of 3 smoke divers, a smoke diving pair (magenta) and a smoke diving leader (blue), during a ∼10 minutes search of the R1 reactor hall and adjacent rooms.The built up rooms/walls are not displayed but can clearly be seen in the search pattern.At the time of the screen-shot, the smoke diving pair has advanced to the second floor as can be seen by the hight estimates displayed on the left side.