A Joint Graph-Based Approach for Simultaneous Underwater Localization and Mapping for AUV Navigation Fusing Bathymetric and Magnetic-Beacon-Observation Data

: Accurate positioning is the necessary basis for autonomous underwater vehicles (AUV) to perform safe navigation in underwater tasks, such as port environment monitoring, target search, and seabed exploration. The position estimates of underwater navigation systems usually suffer from an error accumulation problem, which makes the AUVs difficult use to perform long-term and accurate underwater tasks. Underwater simultaneous localization and mapping (SLAM) approaches based on multibeam-bathymetric data have attracted much attention for being able to obtain error-bounded position estimates. Two problems limit the use of multibeam bathymetric SLAM in many scenarios. The first is that the loop closures only occur in the AUV path intersection areas. The second is that the data association is prone to failure in areas with gentle topographic changes. To overcome these problems, a joint graph-based underwater SLAM approach that fuses bathymetric and magnetic-beacon measurements is proposed in this paper. In the front-end, a robust dual-stage bathymetric data-association method is used to first detect loop closures on the multibeam bathymetric data. Then, a magnetic-beacon-detection method using Euler-deconvolution and optimization algorithms is designed to localize the magnetic beacons using a magnetic measurement sequence on the path. The loop closures obtained from both bathymetric and magnetic-beacon observations are fused to build a joint-factor graph. In the back-end, a diagnosis method is introduced to identify the potential false factors in the graph, thus improving the robustness of the joint SLAM system to outliers in the measurement data. Experiments based on field bathymetric datasets are performed to test the performance of the proposed approach. Compared with classic bathymetric SLAM algorithms, the proposed algorithm can improve the data-association accuracy by 50%, and the average positioning error after optimization converges to less than 10 m.


Introduction
The underwater navigation system is critical in order for underwater vehicles to perform activities that require long-duration and accurate positioning [1,2].Given the absence of Global Navigation Satellite System (GNSS) signals underwater, the principal navigation systems predominantly rely on the Inertial Navigation System (INS), complemented by Doppler Velocity Log (DVL), underwater acoustic transponders, and geophysical field data [3,4].The INS provides autonomous position estimation but is susceptible to error accumulation due to its integral solution method.DVL can help to reduce the INS accumulated error greatly, but the navigation output still contains approximate linear-growth errors.
Acoustic positioning systems, such as long baseline, short baseline, and ultra-short baseline systems, offer precise positions within their working ranges.Nonetheless, their deployment necessitates intricate advance preparation and is quite costly and time-consuming.
Geophysical field-aided navigation methods have garnered significant interest for enhancing the reliability of underwater navigation systems [5][6][7].In areas where prior geophysical field maps are available, the positioning errors can be bounded through datamatching approaches.Conversely, in unknown environments, a Simultaneous Localization and Mapping (SLAM) approach is usually taken to correct the accumulated errors through data association based on the measured data in real time.Several studies involving the utilization of bathymetric data for underwater navigation have been conducted to evaluate practical performance.The Norwegian Defense Research Establishment has carried out a series of simulations and experimental works on terrain-matching algorithms using the HUGIN AUV [8][9][10].Stanford University and the Monterey Bay Aquarium Research Institute collaborated on the development of a long-range and low-cost AUV terrain navigation system [11,12].By using low-precision dead-reckoning (DR) and filtering algorithms, a positioning accuracy of 4-10 m is achieved.Roman C and Singh H proposed an improved SLAM approach for multibeam seafloor mapping [13,14].To create more accurate topographic maps, the multibeam-sonar data is divided into submaps that are used for registration.Meanwhile, the delayed state extended Kalman filter (EKF) is used to re-navigate vehicle trajectories.
The research into geomagnetic field-aided navigation methods for marine vehicles has garnered significant attention.Huang et al. proposed a method for long-range navigation of underwater autonomous underwater vehicles using magnetic anomalies generated by magnetic dipoles [15].In their study, the magnetic-gradient tensor and the draft depth are used to achieve continuous underwater positioning.The relative positioning error is less than 4% on average in simulation experiments.Alimi et al. demonstrated the viability of using passive magnets for underwater navigation, estimating the position of the sensor with a relative error of 1% to 3% through real-world systematic studies [16].Chang et al. proposed a magnetic beacon-based simultaneous localization and mapping (SLAM) approach [17,18].This approach based on particle filter is designed to correct the vehicle-position errors by identifying magnetic beacons in unknown background magnetic fields.Du et al. devised a program to track a mobile-magnetic-dipole target and estimated its magnetic moments using measurements from multiple scalar magnetometers.This program enables rapid detection of the direction, velocity, elevation angle, and magnetic moment of the target [19].Wu and Yao proposed an adaptive SLAM-based method for magnetic-gradient-aided navigation [20].The magnetic-gradient values of landmarks are measured by magnetometers that are installed on an underwater vehicle.Subsequently, an adaptive-unscented-Kalman filter-SLAM method is applied to derive optimal estimations of the vehicle poses.
In underwater SLAM systems, the AUVs detect the seabed using sonar or cameras.The images or bathymetric data obtained are then used to make data associations for positioning-error correction [21,22].The main steps include the following: (1) accurate environmental-data observation; (2) effective data association; and (3) robust state estimation.Magnetic SLAM is mainly based on the magnetic-dipole model, and in order to establish the observation equation, the measurements used are usually magnetic-field intensity or magnetic-field gradient.Lee et al. used the gradient-based magnetic tensor to identify and locate ferromagnetic objects in the environment.Two different types of dense objects are used to verify the effectiveness and accuracy of a proposed method [23].In bathymetric SLAM, the multibeam echo sonar (MBES) method is used to collect bathymetric data in pings, and the iterative closest point (ICP) method is commonly employed for data association.Palomer et al. used MBES to collect three-dimensional bathymetric point clouds and used heuristic ICP to eliminate the influence of point cloud uncertainty on registration accuracy [24,25].Jung et al. proposed a regularized ICP method by re-scaling bathymetric values to improve the data-association accuracy in areas with poor topographic features [26].Bichucher et al. proposed an xy-alignment generalized ICP (GICP) algorithm to create registration between bathymetric submaps [27].This proposed algorithm aims to enforce consistent bathymetric structures and constrain long-term navigation drift.
The SLAM solutions can be broadly categorized into filter-based and graph-based methods.Filtering methods primarily comprise extended Kalman filter (EKF) and particle filter (PF) methods [28,29].For instance, Callmer et al. presented a silent underwater SLAM scheme employing triaxial magnetometers as sensors to locate a vessel with a known magnetic dipole, with trajectory and sensor positions estimated simultaneously using an extended Kalman filter (EKF) [30].Torroba et al. proposed a Rao-Blackwellized particle filter (RBPF) SLAM solution for online large-scale SLAM in an embedded platform, incorporating implicit data association for terrain-loop closure detection and utilizing the particle filter to fuse terrain constraints and estimate AUV poses [31].Additionally, He et al. integrated the particle swarm optimization (PSO) method into Fast SLAM to enhance AUV pose estimation accuracy, adapting PSO to optimize particles and prevent particle degeneracy and impoverishment [32].The filtering methods often have a high computational complexity when the state vector is of high dimension, and their accuracy can be easily affected by unstable measurement noises.The graph-based SLAM approach utilizes factor graphs to represent the relationships between nodes in the SLAM system, where vehicle poses and observed objects are represented as nodes and odometry and loop closures are represented as edges.The optimization of trajectories and maps is achieved using a cost function composed of the edges in the factor graph [33].Graph optimization methods have been applied in BSLAM and magnetic-beacon SLAM studies.Ling et al. used a sparse pseudo-input Gauss process (SPGP) to establish weak data association between adjacent moments and used a terrain-matching method to obtain closed loops among survey moments, and the underwater vehicle positions were estimated using a graph optimization method [34].In order to solve the effect of invalid loop closures on SLAM results, some robust algorithms have been proposed.Sunderhauf and Protzel added switch variables to the objective function to obtain a robust solution for pose graph SLAM.The switch variables can enable or disable the factors associated with them [35].Latif et al. proposed a realizing, reversing, and recovering (RRR) algorithm based on the consensus that correct loop closures tend to agree among themselves and agree with the sequential constraints from odometry, to create robust loop closures in SLAM [36].
This paper proposes a joint graph-based underwater SLAM approach that fuses bathymetric and magnetic beacon measurements.In the front-end, magnetic beacon observations and bathymetric loop closures are generated based on the Euler-deconvolution method and the dual-stage-data-association method, respectively; in the back-end, the generated constraints in the front-end are used to construct the factor graph, and a robust optimizer is used to estimate the AUV path and construct a bathymetric map.The main contributions are as follows: 1.
A novel dual-stage bathymetric data-association method is introduced.The curvature feature vectors on different scales are used to make robust transformations in the first stage, and the GICP algorithm is taken to make fine registrations in the second stage.

2.
Magnetic beacons are used to aid bathymetric data to give better SLAM performance.Also, a new approach that combines Euler-deconvolution and clustering methods is presented that can be used to locate the magnetic beacons accurately in real time.

3.
A false loop-closure-factor-diagnosis mechanism that considers positioning uncertainty, heading, and bathymetric consistency is designed and introduced to the novel back-end optimizer.The robustness of the SLAM system to bad data-association results is subsequently improved.
The remaining sections of the paper are organized as below.The robust dual-stage multibeam-bathymetric-point-cloud-registration method is introduced in Section 2, followed by the magnetic-beacon-detection approach introduced in Section 3. The back-end graph optimizer with a false loop-closure-diagnosis-mechanism design is introduced in Section 4.Last of all, experiments and analysis are given in Section 5.

Dual-Stage Bathymetric Registration Using Curvature-Based Method and GICP
There are usually some outliers in the multibeam bathymetric dataset.To perform bathymetric SLAM in real time, a denoising method is introduced first in this section.At the same time, to improve the data-association accuracy using bathymetric point clouds, the Gaussian-process-regression method is used to model and resample the topographic data.On the basis of this, a dual-stage bathymetric registration approach is presented to make data associations between bathymetric point clouds.A statistical method is used to remove the outliers from the bathymetric point cloud by detecting the bathymetric difference among neighbor points.For each bathymetric point in the cloud, the mean u and standard deviation σ of the K-nearest neighbors is calculated.If a bathymetric point differs from u by greater than 3σ, then it will be taken as an outlier.Specific steps are as follows: (1) Read the bathymetric point cloud S = {b i , i = 1, 2, . . ., n}; (2) Use the k-dimensional tree to index the bathymetric points in the cloud to speed up the searching efficiency for the closest points.
(3) Search the K-nearest neighbors for each bathymetric point b i ∈ S and put them in a set Nb(b i ).
(4) Calculate the bathymetric difference D mid between b i and the mean value of Nb(b i ) as in Equation (1).
If D mid (b i ) > 3σ, then b i will be taken as an outlier and removed.Otherwise, it will be retained.

Bathymetric Down-Sampling Using Gaussian Process Regression
The purpose of bathymetric down-sampling is to reduce the density of the points without reducing the matching suitability.Gaussian process regression (GPR) is a nonparametric model that uses Gaussian process priors to perform regression analysis [37].Assuming that the bathymetric points obey the multidimensional Gaussian distribution, the GPR can achieve high-precision and uniform down-sampling on the bathymetric point cloud.Let z i indicate the bathymetric value of the point x i = b i,x , b i,y .The fitted function between x i and z i can be expressed as Equation (2).
where f(x) = GP(m(x), k(x, x′)) is the Gaussian process with known mean and variance, and σ 2 n is the variance of measurement noise.k(x, x′) is the kernel function that computes the covariance between two bathymetric points x and x′ as, where the first term on the right side represents the square exponential covariance kernel function, and the second term represents the Kronecker delta function.The l, σ f , and σ n values are the hyperparameters of the kernel function.
Using the disordered bathymetric point cloud as the training set to find the optimal hyperparameters, the coordinate x * of an interpolation point is input into the model, and the depths z * can be estimated.
The joint distribution of the training data z and estimated value z * is expressed as where C is the covariance matrix of training data x, C * is the covariance matrix between the train data x and the interpolation points x * , and C * * is the covariance matrix of x * .The z * can then be estimated as shown in Equation (5), In which m(z An example of the denoising and down-sampling of a multibeam bathymetric dataset is shown in Figure 1.
where the first term on the right side represents the square exponential covariance kernel function, and the second term represents the Kronecker delta function.The ,  , and  values are the hyperparameters of the kernel function.
Using the disordered bathymetric point cloud as the training set to find the optimal hyperparameters, the coordinate  * of an interpolation point is input into the model, and the depths z * can be estimated. The

Dual-Stage Bathymetric Data Association
The initial position of two bathymetric point clouds strongly affects the iterative closest point (ICP) results [38], especially for bathymetric maps with low overlap ratio and low matching suitability.To achieve accurate data association, a rough registration method based on feature vectors is proposed in the first stage.Then, in the second stage, the generalized ICP (GICP) is used to obtain accurate loop closure.The detailed process is shown in Algorithm 1.

Dual-Stage Bathymetric Data Association
The initial position of two bathymetric point clouds strongly affects the iterative closest point (ICP) results [38], especially for bathymetric maps with low overlap ratio and low matching suitability.To achieve accurate data association, a rough registration method based on feature vectors is proposed in the first stage.Then, in the second stage, the generalized ICP (GICP) is used to obtain accurate loop closure.The detailed process is shown in Algorithm 1.

Coarse Registration Based on Multi-Scale Feature Vectors
A curvature is selected to represent the topographic change between adjacent bathymetric points.First, the multi-scale curvature feature vectors are established for the points.For each point, the k-nearest neighbor points form a surface S : z = z(x, y).With a definition → r ≜ {x, y, z(x, y)}, the first-order and second-order partial derivative of → r with respect to x and y can be calculated as The unit vector of the principal normal to the surface is and N = → r yy • → n, the average curvature H and Gaussian curvature K of surface S can be calculated by Equations ( 9) and (10), respectively.
where k 1 and k 2 are the maximum and minimum curvatures of the surface S, also called the principal curvature.They can be calculated using Equations ( 11) and ( 12), respectively.
Then, the Gaussian curvature, mean curvature, and principal curvature form a feature vector c as With each bathymetric point as a center, three different-scale square windows are created that contain n 1 , n 2 , and n 3 bathymetric points, respectively, as shown in Figure 2. The feature vector F is formed as Equation (14).
where  and  are the maximum and minimum curvatures of the surface S, also called the principal curvature.They can be calculated using Equations ( 11) and ( 12), respectively.
Then, the Gaussian curvature, mean curvature, and principal curvature form a feature vector  as With each bathymetric point as a center, three different-scale square windows are created that contain  ,  , and  bathymetric points, respectively, as shown in Figure 2. The feature vector  is formed as Equation ( 14). =  ,  ,  (14) where  can be computed using Equation (13). ,  , and  are set to 9, 16, and 25 in this research, respectively.The cosine similarity is used to evaluate the similarity of two bathymetric point clouds, as shown in Equation ( 15) P and Q are used to indicate the source and target bathymetric point clouds, respectively. is the number of bathymetric points in a cloud.S  ,  is the cosine value of the angle between two n-dimensional vectors.The Random Sample Consensus (RANSAC) algorithm [39] is used to select the best transformation between the two clouds.Specific steps are given as follows: (1) Randomly select a translation  between two matched points, between which the curvature feature vectors have the highest similarity.
(2) Use  to transform the source point clouds P. P and Q are used to indicate the source and target bathymetric point clouds, respectively.n is the number of bathymetric points in a cloud.S F P , F Q is the cosine value of the angle between two n-dimensional vectors.The Random Sample Consensus (RANSAC) algorithm [39] is used to select the best transformation between the two clouds.Specific steps are given as follows: (1) Randomly select a translation t i between two matched points, between which the curvature feature vectors have the highest similarity.(2) Use t i to transform the source point clouds P.
(3) Find the closest point pairs in between Q and P, and count the number of point pairs with bathymetric differences less than the consistency threshold τ. (4) Repeat steps (1)-( 3) until enough iterations are performed.Then, the translation t * with the largest n is selected as the final transformation.
This method is illustrated in Figure 3. (3) Find the closest point pairs in between Q and P, and count the number of point pairs with bathymetric differences less than the consistency threshold .(4) Repeat steps (1)-( 3) until enough iterations are performed.Then, the translation  * with the largest  is selected as the final transformation.
This method is illustrated in Figure 3.

Fine Registration based on GICP
Generalized ICP (GICP) is an extension of the traditional ICP algorithm.GICP improves the distance metric of the ICP algorithm with a probability model and is more suitable for dealing with noisy point cloud registration problems [40].
Base on the coarse registration result, the closest point pairs between the source point cloud P = { } ,…, and target point cloud Q =  ,…, can be determined as an initialization for the fine registration.
The transformation between P and Q can be established as in Equation (16).
The bathymetric values of the source points and target points are supposed to obey Gaussian distribution, as  ∼   ,  ,  ∼   ,  .Denoting the true transformation by  =  * •  , the transformation error of  can be calculated as in Equation (17), and the probabilistic distribution can be represented as in Equation (18).
The optimal transformation is obtained by minimizing the transformation error.The  can be calculated using maximum likelihood estimation The complete dual-stage bathymetric data-association approach is summarized in Algorithm 1 below.

Fine Registration Based on GICP
Generalized ICP (GICP) is an extension of the traditional ICP algorithm.GICP improves the distance metric of the ICP algorithm with a probability model and is more suitable for dealing with noisy point cloud registration problems [40].
Base on the coarse registration result, the closest point pairs between the source point The transformation between P and Q can be established as in Equation (16).
The bathymetric values of the source points and target points are supposed to obey Gaussian distribution, as b the transformation error of T can be calculated as in Equation ( 17), and the probabilistic distribution can be represented as in Equation (18).
The optimal transformation is obtained by minimizing the transformation error.The T can be calculated using maximum likelihood estimation The complete dual-stage bathymetric data-association approach is summarized in Algorithm 1 below.
Step3: Computing the similarity use cosine distance, and using KNN algorithm to select the best match index = KNNsearch F P , F Q , ′cosine′ Step4: Selecting the best transform with RANSAC

Detection of Magnetic Beacons
To use the magnetic beacons as landmarks to obtain more loop closures in SLAM, an accurate detection method is essential in order to obtain relative localizations of the beacons.A new magnetic-beacon-detection method is proposed in this paper.First, the Euler-deconvolution method is utilized to obtain the general location of potential beacons near the AUV track.Second, the Euler-deconvolution results are taken as initializations for the nonlinear optimization approach, which is used to obtain the accurate positions of the beacons simultaneously.The detailed process is given in Algorithm 2.
Locating the relative positions of the magnetic beacons and clustering the locations of magnetic beacons for j = 1 : k As the distances between the vehicle and magnetic beacons deployed on seabed are much larger than the beacon size, the beacons can be treated as magnetic dipoles [41].Placing the magnetic beacon on the coordinate origin, the magnetic field it generates at point P(x, y, z) can be simplified as Equation ( 20) where B x , B y and B z are magnetic-field component of the beacons in the x, y and z directions, respectively.µ is the permeability of the medium, and µ 0 = 4π × 10 −7 V•s/(A•m) in vacuum.r is the vector from the beacon to point P, and ∥r∥ = x 2 + y 2 + z 2 .The magnetic-gradient tensor G describes the changes of B x , B y , and B z in the x, y and z directions, as shown in Equation (21).
When the AUV navigates in the magnetic-beacon deployment area, the magnetic strength and gradient it measures are the sum of the magnetic fields of the Earth the and beacons, as shown in Equations ( 22) and (23).
where The tensor Euler-deconvolution (TED) method [42] can be used by AUVs to estimate the position of the magnetic sources.The magnetic anomaly field satisfies the homogeneous Euler equations where n is the structure index of the source.For magnetic-dipole sources, n = 3.Then, we obtain Furthermore, to improve the robustness of TED against measurement noise in source localization, the Eigen analysis of the magnetic-gradient tensor (EGT) is used together with TED to improve the location accuracy.As G is a real symmetrical matrix, the eigenvector set {b 1 , b 2 , b 3 } of G can be obtained.The b 3 corresponds to the small- est eigenvalue and is perpendicular to the plane m − r, in which m is the magnetic moment vector of the beacon.This means that the r to be solved is orthogonal to b 3 .
T and r = r x r y r z T , then Then, we obtain Equation ( 27) by combining Equations ( 25) and (26).Compared with Equation ( 25), the solution space of r is significantly reduced and can be solved with the least square method.
The measurement noise also affects the TED+EGT performance greatly.With an increase in distance between the vehicle and magnetic beacons, the influence of measurement noise on inversion results increases.For each beacon, the inversion results obtained from the vehicle positions close to it show obvious convergence, whereas those obtained from far vehicle positions contain large errors and differ from each other greatly.We use the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) [43] method to cluster the estimated locations of the beacons.
Supposing there are n beacon locations R= {r 1 ,r 2 ,. . .,r n }, the distance between r i and r j can be computed as For a location r i , the points in R that are shorter than ε from the location form a set N (r i ), as in Equation ( 29) If the number of elements in N (r i ) is larger than a given threshold, then r i will be taken as a core point and a new cluster N (r i ) is formed.Then, each element in N (r i ) will be taken as a core point for cluster expansion.All of the locations in R should be processed in this way and those do not fall into any clusters will be neglected.The centers of the clusters will be taken as initial locations of potential magnetic beacons.

Fine Detection with Nonlinear Optimization
The magnetic beacon location and magnetic moments of the beacons are contained in a variable vector m = [m 1 m 2 . . .m n . . .m N ], in which m n = x n , y n , z n , M x,n , M y,n , M z,n , and M x,n , M y,n , M z,n is the magnetic moment vector of the beacon.The m is to be solved using optimization methods.
Given a section of AUV path x = x 1 x 2 . . .x j . . .x J , the magnetic strength and gradient data are used as measurements for optimization, as shown in Equations ( 30) and ( 31) where B j = B x,j B y,j B z,j and G j = G xx,j G xy,j G xz,j G yy,j G yz,j represent the magnetic strength vector and magnetic-gradient tensor measured by the AUV at the j-th position.
Then, the objective function can be built for optimization using the measurements above where b(m) is the function to calculate the magnetic strength vector at the measurement point as Equation ( 20), and g(m) represents the function to calculate the magnetic-gradient vector at the measurement point as Equation (21).Ω B and Ω G are noise covariance matrixes.B b is the background magnetic-field strength obtained by consulting the World Magnetic Model (WMM).

Joint-Factor-Graph Construction and Robust Optimization 4.1. Joint-Factor-Graph Construction
There are two main kinds of nodes in the joint graph, which are the AUV pose nodes x = {x i , i = 1, 2, . . ., I} and the magnetic-beacon nodes m = m j , j = 1, 2, . . ., J .There are three kinds of factors among the nodes, which are the dead-reckoning (DR) factors u r between temporally successive AUV pose nodes, the magnetic-beacon-observation factors u m between AUV pose nodes and magnetic-beacon nodes, and the bathymetric loopclosure factors u c obtained from multibeam-bathymetric-data association, as illustrated in Figure 4. where b() is the function to calculate the magnetic strength vector at the measurement point as Equation ( 20), and g() represents the function to calculate the magnetic-gradient vector at the measurement point as Equation (21).Ω and Ω are noise covariance matrixes. is the background magnetic-field strength obtained by consulting the World Magnetic Model (WMM).

Joint-Factor-Graph Construction
There are two main kinds of nodes in the joint graph, which are the AUV pose nodes  = { ,  = 1,2, … , I} and the magnetic-beacon nodes  = { ,  = 1,2, … , J} .There are three kinds of factors among the nodes, which are the dead-reckoning (DR) factors  between temporally successive AUV pose nodes, the magnetic-beacon-observation factors  between AUV pose nodes and magnetic-beacon nodes, and the bathymetric loop-closure factors  obtained from multibeam-bathymetric-data association, as illustrated in Figure 4.The bathymetric loop-closure factors are obtained using the dual-stage bathymetric data-association approach described in Section 2 and the objective function is shown in Equation (36).
where  is the loop-closure measurement,  is the covariance matrix of the loop-closure-observation model's noise.The computation of a bathymetric loop-closure factor is illustrated in Figure 5 and formulated in Equation (37).The DR factors are obtained from the AUV autonomous navigation systems and the objective function is where J r represents the DR factor, u r i = x DR i+1 − x DR i is the DR increment between time i and i + 1, ∥•∥ 2 represents the 2-norm of the factor error vector, and Σ r i is the covariance matrix of the motion-model noise.
The bathymetric loop-closure factors are obtained using the dual-stage bathymetric data-association approach described in Section 2 and the objective function is shown in Equation (36).
where u c ij is the loop-closure measurement, Σ c i is the covariance matrix of the loop-closureobservation model's noise.The computation of a bathymetric loop-closure factor is illustrated in Figure 5 and formulated in Equation (37).
where  is the loop-closure measurement,  is the covariance matrix of the loop-closure-observation model's noise.The computation of a bathymetric loop-closure factor is illustrated in Figure 5 and formulated in Equation (37).
The x DR i and x DR j are the AUV positions at moments i and j, respectively.The c P and c Q are the centers of the source and target submaps.The x i ′ and c P ′ are the AUV position and center of the source submap after registration, respectively.
The magnetic beacon observation factors can be obtained by the magnetic-beaconlocation approach described in Section 3, and the corresponding objective function is given in Equation (38). where is the relative position of the magnetic beacons to the vehicle pose nodes.Σ m i is the covariance matrix of the magnetic beacon observation noise.m k Indicates the position of the k-th magnetic beacon that will be optimized together with the AUV poses in the back-end.
By combining the three kinds of factors together, the general objective function for back-end optimization can be achieved as

Robust Back-End Optimizer
The bathymetric-loop closures and magnetic-beacon observations generated in the front-end form the observation factor set F = {f 1 , f 2 , . . . ,f m }.To identify and exclude potential false factors in F , a diagnosis mechanism is designed and introduced into the back-end optimizer to make it able to detect false factors during optimization.First, some symbols should be defined for clear explanation.
t n : the moment that the n-th observation factor (bathymetric or magnetic-beacon observation) is obtained.
L n : the whole optimized path obtained at t n using the first n observation factors.L − n : the whole optimized path obtained at t n using only the first n − 1 observation factors.L − n−1,n : the dead-reckoning path obtained between t n−1 and t n .L n−1,n : the optimized path section between t n−1 and t n .P n : the vehicle position at t n on L n .P − n : the vehicle position at t n on L − n .P On the basis of this, a robust back-end optimizer for the joint SLAM is introduced.Once a new observation factor is added into the factor graph and the optimization is updated, a diagnosis mechanism with three conditions will be performed to determine whether the new factor is true and should be reserved.The conditions are given below and illustrated in The first condition is to guarantee that the new factor optimizes the DR path in the uncertainty range.If the optimization result is completely outside the error ellipse of the DR estimate, then it is probably a false factor that results in extra errors.The second condition is based on the fact that the underwater navigation systems of AUVs usually have good accuracy.The heading angles of the optimized path should be consistent with the DR path.Otherwise, the path must have been seriously distorted by the newly added factor.The third condition follows the truth that the bathymetric data observed on different paths in a specific area should have strong consistency.False factors lead to wrong positioning estimates and make bad bathymetric consistency at path-cross areas.In summary, only when all of the three conditions are satisfied can the new observation factor be determined as true, and the optimization will be saved for later SLAM solution.

Experiment and Analysis
Simulation experiments based on field-multibeam-bathymetric data are conducted to test the joint SLAM approach above.The bathymetric data were surveyed by the authors' group using a ship equipped with a SONIC 2024 multibeam echo sounder (R2 SONIC, Austin, TX, USA) and a differential GPS receiver (Hi-Target, Guangzhou, China).The GPS path is about 20 km long and include 2868 poses, as shown in Figure 7.In this experiment, the AUV bathymetric survey process and dead-reckoning estimates are simulated using the bathymetric dataset and GPS path.
An ideal AUV DR path is generated first to follow the GPS path, as shown in Equation (40).Then, velocity noise (0, 0.2 ) m/s and heading angle noise (0, 0.02 ) rad/s are added to simulate an DR path with accumulated errors.
where   ,   and  represent the position, velocity, and heading angle of the AUV at the -th moment, respectively.∆ is the sampling interval and is set to 1s.
The multibeam-bathymetric data used in simulation is surveyed by us using a The first condition is to guarantee that the new factor optimizes the DR path in the uncertainty range.If the optimization result is completely outside the error ellipse of the DR estimate, then it is probably a false factor that results in extra errors.The second condition is based on the fact that the underwater navigation systems of AUVs usually have good accuracy.The heading angles of the optimized path should be consistent with the DR path.Otherwise, the path must have been seriously distorted by the newly added factor.The third condition follows the truth that the bathymetric data observed on different paths in a specific area should have strong consistency.False factors lead to wrong positioning estimates and make bad bathymetric consistency at path-cross areas.In summary, only when all of the three conditions are satisfied can the new observation factor be determined as true, and the optimization will be saved for later SLAM solution.

Experiment and Analysis
Simulation experiments based on field-multibeam-bathymetric data are conducted to test the joint SLAM approach above.The bathymetric data were surveyed by the authors' group using a ship equipped with a SONIC 2024 multibeam echo sounder (R2 SONIC, Austin, TX, USA) and a differential GPS receiver (Hi-Target, Guangzhou, China).The GPS path is about 20 km long and include 2868 poses, as shown in Figure 7.In this experiment, the AUV bathymetric survey process and dead-reckoning estimates are simulated using the bathymetric dataset and GPS path.
experiment area in simulation, as shown in Figure 7.The magnetic moments of the beacons are  = 10 A • m ,  = 2 * 10 A • m , and  = 10 A • m .The magnetic strength and gradient data generated by the beacons in the area can be computed using Equations (20) and (21).The magnetometer and gradiometer make measurement at 1 Hz and the measurement noise is set to (0, 0.1 )nT.

Joint-Factor-Graph Construction
During AUV navigation, the magnetic beacons are detected in real time.The magnetic-beacon-location results and the magnetic-beacon-observation factors are shown in Figure 8a,b, respectively.The nine sub-frames in Figure 8a are used to show the magneticbeacon-detection results clearly.The clustering method based on Euler deconvolution makes good estimations and the optimization approach can further improve the location accuracy.In total, 158 magnetic-beacon observations are reserved in the factor graph and will be taken for back-end optimization.An ideal AUV DR path is generated first to follow the GPS path, as shown in Equation ( 40).Then, velocity noise N 0, 0.2 2 m/s and heading angle noise N 0, 0.02 2 rad/s are added to simulate an DR path with accumulated errors.
where x k , v k and φ k represent the position, velocity, and heading angle of the AUV at the k-th moment, respectively.∆t is the sampling interval and is set to 1 s.The multibeam-bathymetric data used in simulation is surveyed by us using a SONIC 2024 multibeam echo sounder in a 2 km × 2 km area in the Huanghai Sea, China.The water-depth range of the data set is 5.2~69.8m and the bathymetric point resolution is about 5 m.Multibeam sonar collects a ping every second and the bathymetric survey process is simulated correspondingly.Nine magnetic beacons are placed in the experiment area in simulation, as shown in Figure 7.The magnetic moments of the beacons are M x = 10 6 A•m 2 , M y = 2 * 10 5 A•m 2 , and M z = 10 5 A•m 2 .The magnetic strength and gradient data generated by the beacons in the area can be computed using Equations (20) and (21).The magnetometer and gradiometer make measurement at 1 Hz and the measurement noise is set to N 0, 0.1 2 nT.

Joint-Factor-Graph Construction
During AUV navigation, the magnetic beacons are detected in real time.The magneticbeacon-location results and the magnetic-beacon-observation factors are shown in Figure 8a,b, respectively.The nine sub-frames in Figure 8a are used to show the magnetic-beacon-detection results clearly.The clustering method based on Euler deconvolution makes good estimations and the optimization approach can further improve the location accuracy.In total, 158 magnetic-beacon observations are reserved in the factor graph and will be taken for back-end optimization.
netic-beacon-location results and the magnetic-beacon-observation factors are sho Figure 8a,b, respectively.The nine sub-frames in Figure 8a are used to show the mag beacon-detection results clearly.The clustering method based on Euler deconvo makes good estimations and the optimization approach can further improve the lo accuracy.In total, 158 magnetic-beacon observations are reserved in the factor will be taken for back-end optimization.The multibeam bathymetric strips are divided into 125 submaps, as shown in Figure 9.The minimum submap length L min and the maximum submap length L max are set to 100 m and 300 m, respectively.The variance of curvature at all bathymetric points in a submap is used to measure the matching suitability α.During the multibeam survey, once the submap length satisfies L min ≤ L ≤ L max and the matching suitability α is larger than a threshold α thd the submap construction is completed.The α thd is set to 0.01 in this research.
submap length satisfies  ≤  ≤  and the matching suitability  is larger than a threshold  the submap construction is completed.The  is set to 0.01 in this research.Two bathymetric loop closures are taken as examples to show the effectiveness of the data-association method.The pure GICP method is taken for comparison here.The loop closures labeled with 1  ⃝, 4 ⃝, 7 ⃝ and 9 ⃝ in Figure 9 are taken to show the dataassociation results in Figure 10.The green and blue point clouds indicate the true positions of target and source bathymetric submaps, respectively.The red point cloud represents the registration results obtained using corresponding approaches.Obviously, the closer the red point cloud is to the green point cloud, the higher the bathymetric registration accuracy.For the four loop-closure areas, the dual-stage data-association errors are 4.2 m, 4.8 m, 2.6 m, and 1.5 m, respectively, which are much smaller than those for the GICP results.Then, 100 Monte Carlo data-association experiments are made on the ten submap pairs.The average matching errors of GICP, coarse registration, and dual-stage data-association methods are 14.2 m, 9.1 m, and 4.4 m, respectively.Compared with GICP, the proposed dual-stage method improves the data-association accuracy by 51.7%.

Back-End Optimization Results
There are 2868 pose nodes, 10 bathymetric loop-closure factors, 2867 DR factors, and 158 magnetic-beacon-observation factors in the factor graph.To verify the significance of magnetic beacons, the factor graph is optimized in two situations.First, all of the factors and nodes in the joint graph are put into the optimizer designed in Section 4.2.Second, only the bathymetric loop-closure and DR factors are used for back-end optimization of pose nodes.The optimization results and errors computed are shown in Figure 11a,b.It can be seen that both situations can give good optimization results that are much closer to the truth than is DR result.The magnetic beacons can further improve the optimization accuracy compared with the case of using only bathymetric loop closures.
The RRR method [36] and Fast SLAM algorithm [44] are compared at last to evaluate the accuracy and efficiency of the proposed algorithm.Fast SLAM uses 100 particles for filtering estimation.Both the Fast SLAM and RRR methods are performed using the same factor graph with the proposed algorithm.The optimized paths and positioning errors of the three methods are shown in Figure 12.The Fast SLAM method can reduce the DR error to a great extent, but the overall error tends to diverge.The RRR method is able to detect and exclude bad loop closures and has a better performance than Fast SLAM but is inferior to the proposed algorithm.The path estimated by proposed method is the closest to the truth and the average error is less than 10 m.Statistics of the running time of three methods are also made on the computer that has 16G RAM and Intel i7-11700 CPU with 2.5 GHz main frequency.The running time of the proposed SLAM, RRR, and Fast SLAM methods were 521.4 s, 665.6 s, and 752.0 s, respectively.

Coarse registration results
Fine registration results GICP results

Back-End Optimization Results
There are 2868 pose nodes, 10 bathymetric loop-closure factors, 2867 DR factors, 158 magnetic-beacon-observation factors in the factor graph.To verify the significanc magnetic beacons, the factor graph is optimized in two situations.First, all of the fac and nodes in the joint graph are put into the optimizer designed in Section 4.2.Seco only the bathymetric loop-closure and DR factors are used for back-end optimizatio pose nodes.The optimization results and errors computed are shown in Figure 11a, can be seen that both situations can give good optimization results that are much cl to the truth than is the DR result.The magnetic beacons can further improve the opt zation accuracy compared with the case of using only bathymetric loop closures.The RRR method [36] and Fast SLAM algorithm [44] are compared at last to ev the accuracy and efficiency of the proposed algorithm.Fast SLAM uses 100 partic filtering estimation.Both the Fast SLAM and RRR methods are performed using the factor graph with the proposed algorithm.The optimized paths and positioning er the three methods are shown in Figure 12.The Fast SLAM method can reduce t error to a great extent, but the overall error tends to diverge.The RRR method is a  With the true bathymetric map of the research area shown in Figure 13a as a ence, the mapping errors for pure DR, RRR, Fast SLAM, and the proposed algorith computed and shown in Figure 13b-e, respectively.It is obvious that the pro method has the highest mapping accuracy.With the true bathymetric map of the research area shown in Figure 13a as a reference, the mapping errors for pure DR, RRR, Fast SLAM, and the proposed algorithm are computed and shown in Figure 13b-e, respectively.It is obvious that the proposed method has the highest mapping accuracy.

Conclusions
To overcome the shortcomings of using pure seabed topographic data in underwater navigation in unknown environments, a joint graph-based underwater SLAM approach that fuses bathymetric and magnetic-beacon observations is proposed in this paper.In the dual-stage bathymetric data association method, the coarse registration method can provide good initialization for the fine registration method using the GICP algorithm.More accurate bathymetric loop closures have been obtained compared with using pure GICP.The combination of Euler-deconvolution method and clustering algorithms make the vehicle able to identify and localize the magnetic beacons deployed on the seafloor accurately.A robust back-end optimizer is also designed to identify the false loop closures in the factor graph.Experiments verified that much better SLAM estimates can be achieved using the proposed joint-factor graph than using classic bathymetric methods.However,

Conclusions
To overcome the shortcomings of using pure seabed topographic data in underwater navigation in unknown environments, a joint graph-based underwater SLAM approach that fuses bathymetric and magnetic-beacon observations is proposed in this paper.In the dual-stage bathymetric data association method, the coarse registration method can provide good initialization for the fine registration method using the GICP algorithm.More accurate bathymetric loop closures have been obtained compared with using pure GICP.The combination of Euler-deconvolution method and clustering algorithms make the vehicle able to identify and localize the magnetic beacons deployed on the seafloor accurately.A robust back-end optimizer is also designed to identify the false loop closures in the factor graph.Experiments verified that much better SLAM estimates can be achieved using the proposed joint-factor graph than using classic bathymetric methods.However, there are still some limitations in the paper that need further exploration.The first is that the magnetic-beacon detection is sensitive to environmental noise.Shielding the environmental magnetic-field noise generated by the AUV electronic devices is a challenge.The second is that as the magnetic-field strength decays rapidly with the increase of vehicle-to-beacon distance, it will be difficult to use this method in very large sea areas.We will focus on two works in the near future.The first is combining the geomagnetic anomaly field and magnetic beacons to further enrich the factor graph.The second is improving the robustness of the magnetic-beacon-detection algorithm against environment noises.

Figure 3 .
Figure 3. Procedure of coarse bathymetric point cloud registration.

Figure 3 .
Figure 3. Procedure of coarse bathymetric point cloud registration.
cloud P = b P i i=1,...,N and target point cloud Q = b Q i i=1,...,N can be determined as an initialization for the fine registration.

Algorithm 2 : 1 .
the set of magnetic beacons and logging magnetic-beacon observations for m * k in m * opt Optimization of magnetic-beacon positions if m * k not in M m = m, m * Fast Identification Using Euler-Deconvolution and Clustering Methods the superscripts b and e indicate the beacon and Earth, respectively.The B e ≜ from global or local magnetic models accurately.Due to the weak spatial change of the Earth's magnetic field, G e is much smaller than G b and thus can be neglected.

Figure 4 .
Figure 4. Factor graph of the joint SLAM.The arrows indicate the AUV motion direction.The DR factors are obtained from the AUV autonomous navigation systems and the objective function is J = ‖ +  −  ‖

Figure 4 .
Figure 4. Factor graph of the joint SLAM.The arrows indicate the AUV motion direction.
n : the j-th position of the vehicle on L − n−1,n .P j n−1,n : the j-th position of the vehicle on L n−1,n .E(P): the error ellipse of P. b 1 i,z , b 2 i,z : the bathymetric values of the i-th closest point pair in the overlapped area of a bathymetric loop-closure.

21 Condition 3 :Figure 6 .
Figure 6.Schematic diagram of false factor detection in back-end.

Figure 6 .
Figure 6.Schematic diagram of false factor detection in back-end.

Figure 8 .Figure 7 .
Figure 8. Data association by magnetic-beacon observations.The nine sub-frames in (a) are used to show the relative positions of magnetic-beacon-detection results in enlarged views.(a) Magneticbeacon-location result; and (b) the magnetic-beacon-observation factors.The multibeam bathymetric strips are divided into 125 submaps, as shown in Figure 9.The minimum submap length  and the maximum submap length  are set to 100 m and 300 m, respectively.The variance of curvature at all bathymetric points in a submap is used to measure the matching suitability .During the multibeam survey, once the

Figure 8 .Figure 8 .
Figure 8. Data association by magnetic-beacon observations.The nine sub-frames in (a) are show the relative positions of magnetic-beacon-detection results in enlarged views.(a) Ma beacon-location result; and (b) the magnetic-beacon-observation factors.The multibeam bathymetric strips are divided into 125 submaps, as shown in 9.The minimum submap length  and the maximum submap length  are set m and 300 m, respectively.The variance of curvature at all bathymetric points in a su is used to measure the matching suitability .During the multibeam survey, on

Figure 9 .
Figure 9. Multibeam submaps and data association.The symbols ① to ⑩ indicate the numbers of loop closures.Two bathymetric loop closures are taken as examples to show the effectiveness of the dual-stage data-association method.The pure GICP method is taken for comparison here.The loop closures labeled with ①, ④, ⑦ and ⑨ in Figure9are taken to show the dataassociation results in Figure10.The green and blue point clouds indicate the true positions of target and source bathymetric submaps, respectively.The red point cloud represents the registration results obtained using corresponding approaches.Obviously, the closer the red point cloud is to the green point cloud, the higher the bathymetric registration accuracy.For the four loop-closure areas, the dual-stage data-association errors are 4.2 m, 4.8 m, 2.6 m, and 1.5 m, respectively, which are much smaller than those for the GICP results.Then, 100 Monte Carlo data-association experiments are made on the ten submap pairs.The average matching errors of GICP, coarse registration, and dual-stage data-association methods are 14.2 m, 9.1 m, and 4.4 m, respectively.Compared with GICP, the proposed dual-stage method improves the data-association accuracy by 51.7%.

Figure 9 .
Figure 9. Multibeam submaps and data association.The symbols 1 ⃝ to 10 ⃝ indicate the numbers of loop closures.

Figure 10 .
Figure 10.Bathymetric data-association example.The subfigures (a,d,g,j) are the coarse registra results of the loop closures ①, ④, ⑦ and ⑨ in Figure 9. Correspondingly, the subfigures (b,e are the fine registration results, and (c,f,i,l) are the registration results using the GICP method.

Figure 10 .Figure 11 .
Figure 10.Bathymetric data-association example.The subfigures (a,d,g,j) are the coarse registration results of the loop closures 1⃝, 4 ⃝, 7 ⃝ and 9 ⃝ in Figure9.Correspondingly, the subfigures (b,e,h,k) are the fine registration results, and (c,f,I,l) are the registration results using the GICP method.J. Mar.Sci.Eng.2024, 12, x FOR PEER REVIEW

Figure 11 .
Figure 11.Comparison of back-end results with and without magnetic beacons.(a) Comparison of optimized path; and (b) comparison of path errors.

Figure 12 .
Figure 12.Comparison of AUV paths using different SLAM methods.(a) Comparison of opt paths; (b) Comparison of path errors.

Figure 12 .
Figure 12.Comparison of AUV paths using different SLAM methods.(a) Comparison of optimized paths; (b) Comparison of path errors.

Figure 13 .
Figure 13.Bathymetric map errors from different methods.(a) True bathymetric map; (b) mapping error of DR; (c) mapping error of Fast SLAM; (d) mapping error of RRR; and (e) mapping error of proposed SLAM.

Figure 13 .
Figure 13.Bathymetric map errors from different methods.(a) True bathymetric map; (b) mapping error of DR; (c) mapping error of Fast SLAM; (d) mapping error of RRR; and (e) mapping error of proposed SLAM.
joint distribution of the training data z and estimated value z * is expressed as  is the covariance matrix of training data ,  * is the covariance matrix between the train data  and the interpolation points  * , and  * * is the covariance matrix of  * .The z * can then be estimated as shown in Equation (5),