A Scalable Framework for Map Matching Based Cooperative Localization

Localization based on scalar field map matching (e.g., using gravity anomaly, magnetic anomaly, topographics, or olfaction maps) is a potential solution for navigating in Global Navigation Satellite System (GNSS)-denied environments. In this paper, a scalable framework is presented for cooperatively localizing a group of agents based on map matching given a prior map modeling the scalar field. In order to satisfy the communication constraints, each agent in the group is assigned to different subgroups. A locally centralized cooperative localization method is performed in each subgroup to estimate the poses and covariances of all agents inside the subgroup. Each agent in the group, at the same time, could belong to multiple subgroups, which means multiple pose and covariance estimates from different subgroups exist for each agent. The improved pose estimate for each agent at each time step is then solved through an information fusion algorithm. The proposed algorithm is evaluated with two different types of scalar field based simulations. The simulation results show that the proposed algorithm is able to deal with large group sizes (e.g., 128 agents), achieve 10-m level localization performance with 180 km traveling distance, while under restrictive communication constraints.


Introduction
Cooperative multi-agent systems have become increasingly popular due to the wide range of applications that they support, such as surveillance [1], search and rescue [2], and exploration [3]. In these applications, high-quality localization, the ability for agents to reliably and accurately estimate their poses (i.e., positions and orientations) with respect to the surrounding environment or to a geographic coordinate system, is crucial. The GNSS, as a traditional solution, is not always available or reliable, due to reasons such as signal blockages, multipath reflection, and jamming [4]. One potential solution for localization in GNSS-denied environments is to utilize map matching techniques, given a prior map represented as a scalar field. Scalar fields associate a scalar value with every point in space, and applications include gravity anomaly [5], magnetic anomaly [6,7], topographic [8], and olfaction [9], to name a few. Methods utilizing scalar fields for localization regulate agents' dead-reckoning error growth through matching the information measured by on-board sensors with the prior given scalar field maps, such as terrain-aid navigation [10] and magnetic anomaly-based navigation [7]. However, these methods are sensitive to the characteristic information available in the local area near the agent, sensor noises, and the resolution and accuracy of the given maps. For example, for a single agent localization, the on-board sensor measurements could match to multiple positions on a scalar field map, creating ambiguity. This can be alleviated through matching the past sensor measurements along the agent's trajectory to the map [11] but is still often not robust in real-world applications. Fortunately, compared to a single agent, a group of collaborative agents may provide several navigational benefits, such as tolerance against individual sensor failures. This can be achieved through sharing observations across a large spatial area on the scalar proposed algorithm is introduced in Section 5, with results discussed in Section 6. The paper concludes in Section 7.

Fusion of Multiple Estimates Full Group of Agents
Examples of Local Subgroups The geometry of the subgroups (i.e., the relative positions inside the subgroups) are estimated, using range-only measurements, then the geometry is used to extract measurements of the scalar field to estimate the pose and associated uncertainty of each member in the subgroup. (Lower right) An agent, as an example, receives multiple copies of its pose estimate through its membership in several subgroups and fuses them to reduce pose error.

Related Works
Scalar field based localization system designs have been researched in applications such as gravity-aid navigation [11,18,19], magnetic anomaly-based navigation [7,20], and terrain-based navigation [21]. In order to perform robust localization in featureless areas or with low-quality sensors, cooperative multi-agent localization systems are proposed to achieve accurate estimations [22][23][24]. Distributed multi-agent localization methods were first formulated based on Kalman filters [25,26]. Even though these methods allow to perform an observation update and data exchange when agents are within the communication range, each agent in the group is required to estimate the poses of all agents, which does not scale well to large groups. Meanwhile, the Kalman filter-based estimation methods assume that the pose estimate can be presented by a unimodal Gaussian distribution. However, the scalar field-based estimation error distributions are usually multi-modal and difficult to be approximated by the Gaussian distribution [7,12].
Canciani et al. formulated the magnetic anomaly-based cooperative navigation problem as a particle filter [14]. The method does not scale for large groups, due to the use of a centralized particle filter. Our previous works in [6,15] broke the localization process into two steps: the relative poses between agents are estimated using inter-agent ranging measurements through an extended Kalman filter (EKF), and then each agent estimates its pose using all magnetic anomaly measurements and relative poses of the group through a particle filter. Although the particle filter in [6,15] only contains four states, the EKF formulation, which includes all agents' poses, does not scale to large group sizes.
A decentralized cooperative bathymetry based localization method was proposed in [17]. In [17], each agent is able to estimate its pose through matching altimeter measurements with a bathymetric map, using a marginalized particle filter. Then, the Gaussian belief, estimated based on the inter-agent ranging measurement and the other agent's position estimate, is applied to update the particles in the filter. Although this method is able to achieve scalable cooperative localization, it ignores the correlation of the information, which may lead to over convergence. Rui et al. then presented an extended information filter to address the issues about the correlation of the information [27]. However, the method described in [27] is reliable to GPS measurement or highly accurate bathymetric information-based estimations for the prediction update, which leads this method to be non-feasible in an underwater environment.
Wiktor et al. presented a decentralized CI based collaborative multi-agent localization algorithm applied in natural terrain-aid navigation [12]. Similar as [17], each agent is assumed to perform terrain-aid navigation to estimate its own pose and related covariances. The pose estimates are updated using inter-agent ranging measurements and agents' pose and covariances through a CI filter, which can fuse estimates with unknown correlation. One potential limitation of this method is the robustness within feature-poor regions due to the single measurement used in map matching for each agent at each time step. Compared with those methods presented in [12,17,27], which only use the information from the immediate neighbors, in our proposed algorithm, a subgroup strategy, instead of utilizing single measurements, is applied to improve the robustness of the localization system.
Active multi-agent navigation algorithms [28,29], which combine localization and active path planning algorithm, are interesting research directions to improve the robustness of the pose estimation. However, to the best knowledge of the authors, there is so far no active multi-agent navigation algorithms focused on localization based on map matching using scalar field information.

Problem Statement and Notations
In this study, the case is considered where a large group of N agents (e.g., space, aerial, ground, surface, or underwater vehicles) is entering an environment without GNSS. At the moment entering the GNSS-denied environment, the initial pose in the global frame for each agent is assumed to be known with a small uncertainty (e.g., provided by GNSS). The main objective is to achieve accurate and robust global localization for all agents in a large group mainly based on local communication and inter-agent ranging measurements along with scalar field information. In this paper, a one-based numbering system is used for indexing. Let the group of N agents be denoted by Φ = {φ 1 , φ 2 , · · · , φ N } where φ i denotes the ith agent in the group. The pose of each agent is partially observable; thus, at each time step, the pose is represented as a Gaussian belief and φ i [P k ] are the mean vector and covariance matrix of the state of the ith agent at the kth time step.
Each agent is assumed to be equipped with radios that enable communication with nearby agents to exchange information, as well as to perform undirected inter-agent ranging. Due to communication limitations (e.g., range or number of channels), each agent can only communicate and perform ranging with a limited number of agents inside a certain range. Therefore, the agents are divided into subgroups based on the communication constraints as follows: where ζ(φ i , φ j ) is a user defined function indicating if agent i is capable of communicating with agent j and is a user defined threshold. For example, ζ(φ i , φ j ) may be based on signal strength, and may define a threshold determining whether the signal strength is adequate for agent i to communicate with agent j. Based on this definition, the group has at least N subgroups, each with M members, where subgroup Φ i corresponds to the agents capable of communicating with agent φ i . The agent i is named the fusion center agent of the subgroup i. Furthermore, the subgroups are not disjointed (i.e., Φ = N i=1 Φ i ), as each agent is included in multiple subgroups. In other words, the situations of isolated agents and isolated subgroups are not considered in this study. The beliefs of subgroup i at time step k are denoted by the following: are similarly the set of means and covariance matrices of subgroup i at time step k. Each agent in a subgroup measures the distance between other agents in the subgroup. Let the distance between agent i and agent j be denoted by d[φ i , φ j ], and let the set of distance measurements between agents in a subgroup be denoted by Each agent is assumed to be loaded with a prior scalar field map covering the operation area. Since the types of scalar fields evaluated in this study change very little in a short time frame [7,8], the prior loaded scalar field map is assumed to be deterministic. Each agent also performs real-time measurements of the scalar field at the current location with sensor noises, such as gravity anomaly, magnetic anomaly, or altimeter measurements, at each time step, which can be exchanged through the communication links. In general, scalar fields may vary with altitude depending on the applications (e.g., magnetic anomaly) and the available data at different altitudes, which may not be available or dense enough; however, this study only deals with the navigation problem in 2D (i.e., assuming the agents are moving at a constant altitude). The scalar field measurements of the ith agent are denoted by φ i [I], and the sets of these measurements in a subgroup are denoted by Φ i [I].
Meanwhile, each agent is assumed to be able to measure or estimate its velocity at each time step, which could be achieved through utilizing the Doppler velocity log [30], wheel odometry [31], or vehicle dynamic model [32], along with yaw rate measurements, using gyroscopes. Similar to before, the velocity and yaw rate measurements of the ith agent are denoted by φ i [v] and φ i [ω], and the sets of measurements in a subgroup are denoted by

System Overview
The proposed algorithm is developed for estimating each agent's global pose in a scalable group with communication constraints. The pipeline of the proposed multisubgroup based algorithm is shown in Figure 2. Each agent in the group is assigned to different subgroups based on the communication range or the number of channels, and these subgroups overlap, as discussed earlier in Section 3. Note that in this study, in order to conveniently evaluate the performance of the proposed algorithm for each agent in the group, the subgroups are pre-set at the beginning, which means that the memberships of each subgroup are assumed to be constant during an operation. Without pre-setting the subgroups, each agent in the group cannot be guaranteed to belong to the same subgroups during the operation, which could lead the performance of each agent to be unequal and difficult to evaluate. Each subgroup takes each agent's velocity, yaw rate, scalar field measurements, as well as the ranging measurements between the agents in the subgroup as inputs. The outputs of each subgroup are pose estimates of the agents belonging to this subgroup. As mentioned earlier, each agent could belong to multiple subgroups. Thus, multiple pose and covariance estimates from different subgroups exist for each agent in the group at each time step. For example, each subgroup contains M agents and each agent is a member of L subgroups. Therefore, L estimates of the pose and covariance can be communicated to that agent. For each agent, the information fusion algorithm, explained in Section 4.3, is developed for estimating the agent's pose based on this redundant information.

Cooperative
Ranging Localization

Subgroup Cooperative Localization
The cooperative localization algorithm that was developed in our previous work [6] is applied to obtain the estimated poses and covariances of all agents in a subgroup at each time step. For completeness, the main approach and key equations are explained in this subsection. The pipeline of the subgroup process is described in Figure 3. Specifically, two sequential algorithms, named cooperative ranging localization and cooperative scalar field localization, are included in this process and performed in the fusion center agent in each subgroup. The cooperative ranging localization step estimates the relative poses between each pair of agents in the subgroup. The cooperative scalar field localization step then estimates the pose and covariances of the fusion center agent in the subgroup through matching the multiple scalar field measurements with relative pose constraint to the prior scalar field map. Since the processes in the rest of this subsection are identical for each subgroup, in order to present equations clearly, all notations listed in the rest of this subsection do not include the index of the subgroup.

Cooperative Ranging Localization
The cooperative ranging localization algorithm is formulated as an EKF with the state vector given by the following: (1) , π (2) , . . . , π (M) ] T where π (i) = [x (i) , y (i) , θ (i) ] T is the pose of the ith agent in the subgroup in the global frame where x and y are the Cartesian coordinates for the position, θ is the heading, and M is the number of agents in the subgroup. The states, shown in Equation (3), at the next time step are predicted using each agent's velocity and yaw rate measurements and current states. The inter-agent ranging measurements are applied to perform observation updates.
The observation function is given by the following: where (x (i) , y (i) ) and (x (j) , y (j) ) are Cartesian coordinates of an arbitrary pair of agents in the subgroup. At each time step, the relative pose of each agent with respect to other agents in the Cartesian coordinates can be calculated as the following: where [x (i/j) ,ŷ (i/j) ,θ (i/j) ] T are the relative pose for agent i with respect to agent j in the subgroup, [x (i) ,ŷ (i) ,θ (i) ] T and [x (j) ,ŷ (j) ,θ (j) ] T are the pose estimates of agent i and agent j in the subgroup from EKF, respectively.

Cooperative Scalar Field Localization
With the estimated relative poses in the subgroup, the scalar field measurements from each agent can be treated as points on a fix geometric shape to match to the map. The cooperative scalar field localization method is formulated as a particle filter involving only four states as follows: where x and y are the Cartesian coordinates for the global position of the fusion center agent in the subgroup, θ is the heading of the agent, and γ is the rotation error of the subgroup, which is explained in detail in [15]. In brief, the rotation error is defined as the rotation angle in the global frame of the estimated geometric structure of the subgroup from the cooperative ranging localization. This is necessary since the EKF in the cooperative ranging localization is only capable of preserving the pairwise distances, and as a result, only the geometric structure of the subgroups are maintained through the EKF. The relative position estimates of each agent with respect to the other agents in the subgroup can be corrected based on γ by the following: where [x (i/j) ,ỳ (i/j) ] T are the relative position coordinates for agent i with respect to agent j in the subgroup after applying the update utilizing γ, and [x (i/j) ,ŷ (i/j) ] T are the relative position coordinates for agent i with respect to agent j in the subgroup from Equation (5).
Since the particle filter is performed in the fusion center agent, for each particle, the states, shown in Equation (6), at the next time step are predicted, using the fusion center agent's velocity and yaw rate measurements and the state estimate at the current time step. Note that since the rotation error changes slightly and is difficult to predict, the state γ is propagated by a random walk. The observation model for each particle is given by the following: where y is the observation vector with length m and η is the measurement noise vector. The relative position r is calculated from Equation (5). Therefore, in the subgroup, the other agents' predicted global positions can be extracted by adding the updated relative positions from Equation (7) to the predicted states from the particle filter. The vector-valued observation function h M is used to extract the predicted scalar field measurements from the given map based on each agent's predicted position.

Information Fusion
The goal of the information fusion step is to estimate an improved pose for each agent. Since the correlation among the pose and covariance estimates is unknown due to information reuse between subgroups, the CI algorithm [33] is applied to fuse the estimates of each agent from multiple subgroups. In general, the CI algorithm fuses information based on a convex combination of the information matrices (i.e., the information states and the inverse of corresponding covariance matrices). Therefore, to apply CI, each agent must have knowledge of the expectations and associated covariance matrices estimated in the subgroups it belongs to. The procedure for obtaining these values is described in detailed in this section.

Covariance Estimates From Particles
The expectation of the state at time k for the fusion center agent in each subgroup, denoted asp k , can be obtained from the particle filter by computing the weighted average of the particles given the particle weights. Let l p(i) denote ith state of the lth particle, and letp(i) denote the expected value of the ith state. Now, the covariance matrix of the state denoted by C is calculated from the particles as follows: where iw is the weight of the ith particle, n is the total number of particles, and l and p represent the indices of the elements in the matrix. Then, the covariance matrix of the pose estimate of each agent at time step k is given by the following: where φ i [P k ](l, p) denotes the element in the lth row and pth column of the covariance matrix of ith agent. The pose estimates of other agents in the subgroup are able to be derived by utilizing the fusion center agent's pose estimate, the relative pose from Equation (5), and the subgroup rotation error γ. According to [34], the covariance matrix of the pose estimate for each agent in the subgroup at one time step is the same as the covariance of the fusion center agent.

Covariance Intersection
In order to present the details of CI applied in this problem clearly, denotes the jth estimate of pose and covariance of ith agent at time step k. The implementation of CI in this situation is presented as follows: where α = [α 1 , α 2 , · · · ] are the weighting coefficients, which need to be solved, and L is the total number of solutions of the ith agent (i.e., the total number of subgroups, including agent i). Based on [33], α is assumed to be computed by solving a nonlinear optimization problem, which is to find α to minimize the determinant or trace of φ i [P k ]. However, to solve the nonlinear optimization problem is computationally expensive, especially, when the ith agent has a large number of solutions. To reduce the computation cost, the fast CI method presented in [35], which is designed to solve the weighting coefficients quickly, is applied as follows: By solving for α from Equation (13), the improved pose of the agent at time step k (i.e., φ i [x k ]) can be computed through Equations (11) and (12).

Simulation
In order to verify that the proposed algorithm is suitable for different types of scalar fields, a magnetic anomaly map and a bathymetric map are utilized in the simulation study.
For the ease of simulating a large number of agents, the agents are simulated with a bicycle model [36] in the 2D plane. A feedback controller steers an agent to follow a reference trajectory and reference velocity based on the pose estimates from the proposed algorithm. All reference trajectories are set to be parallel to latitude lines in a geographic coordinate system. New parallel reference trajectories are added with a certain distance when the group size increases. Agents are assumed to be moving from west to east with the same initial longitude and small uncertainties (a Gaussian white noise with zero mean and 1 m standard deviation) on the initial poses as discussed in Section 3. In order to simulate situations when all agents are moving at different speeds, the reference velocity for each agent varies slowly, following a sine function with a random offset throughout the duration of the mission.
In this simulation, the communication among agents is restricted by the limited number of communication channels available for each agent. The fusion center agent and agents it can reach are assigned as a subgroup at the beginning of the simulation. Meanwhile, different scalar fields are related to different experimental environments and platforms. For example, magnetic anomaly maps can be used by the localization system of aerial vehicles on Earth [7], or space probes exploring other planets, and bathymetric maps are used for underwater vehicles' localization [12]. Therefore, different parameter settings are introduced corresponding to two map types in the following.

Magnetic Anomaly Map
The Earth's magnetic anomaly information that presents the variations of the Earth's magnetic field is stable and distinctive at different locations in a certain range of altitude [7]. The magnetic anomaly map utilized in this simulator, shown in Figure 4, is obtained from the United States Geological Survey [37] and contains the magnetic anomaly information at 305 m altitude from the area around Columbus, Ohio, United States. The noise of the velocity measurement and yaw rate measurement are drawn from a Gaussian distribution with zero mean and standard deviation σ v = 0.3 m/s for velocity and σ g = 0.005 deg/s for the yaw rate. A turn-on bias b v ∼ N(0, 0.1σ v ) for velocity and b g ∼ N(0, 0.1σ g ) for yaw rate is added separately. According to the references [7,38], the noise of the ranging measurements and the magnetic anomaly measurements are drawn from a Gaussian distribution with zero mean and standard deviation σ r = 1 m and σ m = 10 nT, respectively. The initial distance between each pair of neighbor agents in latitude is set to 1000 m. The reference velocity for each agent varies from 40 to 60 m/s. The mission duration is 1 h, and the trajectory length of each agent is about 180 km.

Bathymetric Map
The bathymetry, also known as submarine topography, presents the depths of the underwater terrain. The bathymetric map applied in this simulator, shown in Figure 5, is acquired from the United States Geological Survey [39]. According to the descriptions of parameter setting in underwater experiments in [12,13], the standard deviation of the velocity measurement noise and yaw rate measurement noise are selected as σ v = 0.1 m/s and σ g = 0.1 deg/s. A turn-on bias b v ∼ N(0, 0.1σ v ) for velocity and b g ∼ N(0, 0.1σ g ) for the yaw rate is also added separately. The noise of the altimeter measurements is drawn from a Gaussian distribution with zero mean and standard deviation σ a = 1 m. The noise of the inter-agent ranging measurements are drawn from a Gaussian distribution with zero mean and standard deviation σ r = 1 m. The initial distance between each pair of neighboring agents in latitude is set to be 200 m. The reference velocity for each agent varies from 0.5 to 1.5 m/s. The mission duration is 1 h and the trajectory length of each agent is about 3600 m.

Results and Discussions
The proposed algorithm is designed to remove the constraints of the group size in our previous works in performing cooperative localization [15]. Therefore, the algorithms presented in [15] are utilized as the comparison to evaluate the proposed algorithm. The performance of the proposed algorithm is compared with simulations under the full communication (FC) assumption, where all agents can communicate with each other agent at each time step. The performances and sensitivity analysis of FC could be found in [15]. The average position root mean squared errors (RMSE) in each simulation from both FC and the proposed algorithm are evaluated. In each case, multiple Monte Carlo simulations are performed, consisting of 160 trials each.
The proposed algorithm and FC are first evaluated, using amagnetic anomaly mapbased simulation environment. Figure 6 shows the cumulative distribution function (CDF) of each simulation's average position error for all agents in the group at each time step from the Monte Carlo simulations.  Table 1 shows the statistical data of the Figure 6 along with performance of the simulation with dead-reckoning (DR) without using the ranging and magnetic anomaly information. From Figure 6 and Table 1, it is clear that the performance of both the proposed algorithm and FC improves as the group size increases, and both are much better than doing DR alone. Note that since the number of particles used in different scenarios is the same at the cooperative scalar field localization step, the performance comparison is suitable to be made between FC and the proposed algorithm, when the group size of FC is equal to the subgroup size of the proposed algorithm. In this case, with the same number of particles, the particle filter in both situations deals with the same number of scalar field measurements at each time step. Meanwhile, the computation of each agent applying the proposed algorithm is similar to the computation of each agent working in FC because of the fast CI method [35] applied. The performance of the proposal algorithm shows improvement over the FC in Figure 6 and Table 1.
Similar results are acquired using a bathymetric map-based simulation environment. Figure 7 shows the CDF of each simulation's average position error for all agents in the group at each time step. Due to the smaller bathymetric map, a maximum of up to 16 in group sizes were simulated.  Table 2 shows the statistical data of the Figure 7 along with performance of the DR. The performance of the proposed algorithm applied with different group sizes and the same subgroup size using the magnetic anomaly map is shown in Table 3. It shows that the proposed algorithm works with large group sizes (e.g., N = 128). The computation of each agent in the group, which is only spent on locally centralized cooperative localization and information fusion, does not increase when the group size grows. However, since each agent uses the same amount of information in different sizes of groups (i.e., the sizes of the subgroups are the same), the performance of these simulations are similar.

Conclusions and Future Work
In this paper, a scalable cooperative localization framework based on scalar field prior maps and real-time measurements is presented. In order to satisfy the communication constraints, a large agent group is separated into several subgroups, where each agent is treated as the fusion center in each subgroup. A locally centralized cooperative localization is performed to estimate the agents' poses in each subgroup through matching multiple scalar field measurements constrained by relative positions to the given map. In order to avoid over-convergence due to using correlated information, a fast CI algorithm is applied to estimate an improved pose for each agent based on its multiple pose and covariance estimates from its membership in multiple subgroups.
The simulation results show that the proposed algorithm is able to deal with large groups (e.g., N = 128), and to achieve higher performance, even under more restrictive communication constraints, compared to our previous works [6,15]. Additionally, this approach is suitable for missions with different types of scalar fields.
There are several limitations in this work that need to be addressed in the future, such as extending the algorithm to agents distributed in 3D spaces and better integration of the uncertainty in the group geometry estimates into the map-matching process. A great challenge is to find decentralized localization solutions that can utilize all agents' measurements in the group in an efficient and robust manner, under communication constraints. Currently, the performance of the proposed algorithm is limited by the subgroup size, instead of the full group size. Our future work will focus on allowing information to flow beyond the immediate neighbors while maintaining the stability of the pose estimation algorithm in the agent network.