Research Advances and Prospects of Underwater Terrain-Aided Navigation

: Underwater terrain-aided navigation (TAN) can obtain high-precision positioning independently and autonomously under the conditions of a communication rejection space, which is an important breakthrough for the autonomous and refined operation of deep-sea autonomous under-water vehicles near the seabed. Although TAN originated in the aviation field, the particularity of the underwater physical environment has led to the formation of a different theoretical and technical system. In this article, the application background, operating principles, and most important technical aspects of underwater TAN are introduced. Then, the relevant algorithms involved in the two main modules (the terrain-aided positioning module and the iterative filtering estimation module) of the underwater TAN are reviewed. Finally, other cutting-edge issues in the field of underwater TAN are summarized. The purpose of this article is to provide researchers with a comprehensive understanding of the current research status and possible future developments in the TAN field.


Introduction
With the advancement of ocean exploration and development into deep-sea spaces, underwater vehicles face increasing challenges in terms of communication-limited and communication-denied operational environments [1,2].This often results in the inability of human intervention to provide timely and accurate information to guide their mission progress.Therefore, underwater vehicles tasked with deep-sea operations must possess autonomy and intelligent operational capabilities to meet the long-term and precision requirements of future deep-sea operations.High-precision navigation technology serves as the foundation for the safe navigation and autonomous intelligent operation of underwater vehicles, and its theoretical and technical levels are directly related to the level of underwater vehicle intelligence.
Owing to the inability of satellite signals to penetrate seawater and reach the deep sea, the scope of acoustic navigation systems for underwater vehicles is limited.Inertial navigation error accumulates infinitely over time.Traditional navigation methods can no longer support underwater vehicles in performing long-term and fine operation tasks in deep-sea spaces [3,4].Therefore, the research and development of high-precision, long-term autonomous navigation technology that is adaptable to deep-sea environments has become an urgent need in deep-sea exploration and development.
Deep-sea terrain-aided navigation (TAN) utilizes underwater terrain features as positioning references.It achieves the real-time positioning of underwater vehicles in a priori terrain maps by actively detecting and tracking characteristic terrains with positioning errors bounded in the time and spatial domains [5].Therefore, the TAN is a significant breakthrough in achieving the long-term, high-precision positioning of underwater vehicles in the deep sea.As the demand for deep-sea and polar seabed exploration and vehicles.In the 1980s, the Norwegian Defence Research Establishment (FFI, Forsvarets forskningsinstitut) introduced TAN technology into the field of underwater navigation.It successfully developed the underwater TAN system "TerrNav" and the experimental platform "TerrLab" [46,47].During the past four decades, the development of underwater TAN technology has led to the formation of a relatively complete theoretical and technical system (Figure 1).The theoretical and technical system of underwater TAN mainly includes three parts: terrain measurement and restructuring [48], TAP [48,49], and the path planning of underwater TAN [50].Among them, terrain measurement and reconstruction include measuring and reconstructing the global a priori DEM and the real-time measured terrain map, on which the relevant research is already very comprehensive.TAN/TAP and TAN path planning are the current research hotspots of TAN.TAP research mainly involves single-point TAP estimation [49], continuous recursive positioning estimation [51,52], and positioning in terrain-unknown environments [53].TAN path planning mainly involves Terrain Adaptability Analysis [54], path search, and optimization methods [55].

Composition of Two TAN Systems
The underwater TAN system mainly consists of three components: the main navigation system, the terrain-aided positioning (TAP) module (including the terrain measurement module), and the localization information fusion module.The main navigation system provides the basic position reference information for the underwater TAN system.

Composition of Two TAN Systems
The underwater TAN system mainly consists of three components: the main navigation system, the terrain-aided positioning (TAP) module (including the terrain measurement module), and the localization information fusion module.The main navigation system provides the basic position reference information for the underwater TAN system.The primary navigation localization module is generally based on the dead reckoning (DR) or inertial navigation systems (INSs), in which the navigation positions have recursive relationships and strong correlations.However, the positioning errors of the main navigation system accumulate over time without bounds.The TAP module is used to calculate the position probability distribution of the real-time measured map (RTM) matched to the a priori terrain map in the spatial domain.Because the RTM corresponds to the autonomous underwater vehicle (AUV) position, the position probability distribution of the RTM can also describe the positioning probability distribution of the AUV in the a priori digital elevation model (DEM).The TAP module includes the terrain measurement module, which is primarily used to obtain real-time terrain information beneath the AUV.Various sensors can be used for this purpose, such as single-beam echo sounders [46,56,57], multibeam echo sounders [3,4], the Doppler velocity log (DVL) [35,58], and structured light [59][60][61].The TAN filter estimation module serves as the data fusion center for the entire system.It is responsible for receiving position information from both the main navigation module and the TAP module and then obtaining new positioning information through fusion estimation and iterative computation.The main modules and information flow of the TAN system are shown in Figure 2. The primary navigation localization module is generally based on the dead reckoning (DR) or inertial navigation systems (INSs), in which the navigation positions have recursive relationships and strong correlations.However, the positioning errors of the main navigation system accumulate over time without bounds.The TAP module is used to calculate the position probability distribution of the real-time measured map (RTM) matched to the a priori terrain map in the spatial domain.Because the RTM corresponds to the autonomous underwater vehicle (AUV) position, the position probability distribution of the RTM can also describe the positioning probability distribution of the AUV in the a priori digital elevation model (DEM).The TAP module includes the terrain measurement module, which is primarily used to obtain real-time terrain information beneath the AUV.Various sensors can be used for this purpose, such as single-beam echo sounders [46,56,57], multibeam echo sounders [3,4], the Doppler velocity log (DVL) [35,58], and structured light [59][60][61].The TAN filter estimation module serves as the data fusion center for the entire system.It is responsible for receiving position information from both the main navigation module and the TAP module and then obtaining new positioning information through fusion estimation and iterative computation.The main modules and information flow of the TAN system are shown in Figure 2.

Mathematical Model of the TAN System
The operation of an underwater TAN system can be described using a set of state transition equations.However, the specific forms of the state transition equations vary depending on the structure of the TAN system.The main structural modes of the TAN system include tightly and loosely coupled filtering modes [35].The loosely coupled filtering mode, as shown in Figure 3, is the most commonly used.In this structure, the main navigation system outputs the positioning information only to the TAN system.The filter receives the positioning information from the main navigation system and the TAP information, performs fusion estimation for positioning integration, and outputs the fused positioning information.The main navigation system no longer provides independent positioning results.The state transition equations for this structure are as follows.Equation (1) describes the state transition of the main navigation system, and the second equation represents the state observation equation based on terrain matching.

Mathematical Model of the TAN System
The operation of an underwater TAN system can be described using a set of state transition equations.However, the specific forms of the state transition equations vary depending on the structure of the TAN system.The main structural modes of the TAN system include tightly and loosely coupled filtering modes [35].The loosely coupled filtering mode, as shown in Figure 3, is the most commonly used.In this structure, the main navigation system outputs the positioning information only to the TAN system.The filter receives the positioning information from the main navigation system and the TAP information, performs fusion estimation for positioning integration, and outputs the fused positioning information.The main navigation system no longer provides independent positioning results.The state transition equations for this structure are as follows.Equation (1) describes the state transition of the main navigation system, and the second equation represents the state observation equation based on terrain matching.
Here, x t−1 represents the positioning output of the TAN system at time t, u t−1 represents the input of the TAN system at time t, x t represents the positioning output of the main navigation system at time t after receiving the input u t−1 , v t−1 represents the system input error, y t represents the sequence of real-time sonar measurements of the terrain height obtained at time t, h(x t ) represents the sequence of the interpolated terrain height at the positioning point x t , and e t represents the terrain measurement error.

Mathematical Model of the TAN System
The operation of an underwater TAN system can be described using a set of state transition equations.However, the specific forms of the state transition equations vary depending on the structure of the TAN system.The main structural modes of the TAN system include tightly and loosely coupled filtering modes [35].The loosely coupled filtering mode, as shown in Figure 3, is the most commonly used.In this structure, the main navigation system outputs the positioning information only to the TAN system.The filter receives the positioning information from the main navigation system and the TAP information, performs fusion estimation for positioning integration, and outputs the fused positioning information.The main navigation system no longer provides independent positioning results.The state transition equations for this structure are as follows.Equation (1) describes the state transition of the main navigation system, and the second equation represents the state observation equation based on terrain matching.Here, x t−1 represents the positioning output of the TAN system at time t, u t−1 represents the input of the TAN system at time t, x t represents the positioning output of the main navigation system at time t after receiving the input u t−1 , v t−1 represents the system input error, y t represents the sequence of real-time sonar measurements of the terrain height obtained at time t, h(x t ) represents the sequence of the interpolated terrain height at the positioning point x t , and e t represents the terrain measurement error.
Another structural mode of the TAN system is the tightly coupled mode, as shown in Figure 4.In this structure, the main navigation system outputs attitude information, velocity information, and other parameters to the TAN system.The TAN filter utilizes velocity information, attitude information, and terrain measurement information to perform positioning estimation.Furthermore, the main navigation system independently performs positioning estimation and outputs the positioning information separately, without interfering with the TAN positioning information.The state transition equations for this structure are shown in Equation (2), where the first equation describes the state transition of the main navigation system, and the second equation represents the state observation equation based on terrain matching.
where x k and x k−1 represent the system state vectors obtained at times k and k − 1, respectively.In addition, , where p k = [x N , x E , δ Z ] represents the coordinates of the AUV main navigation system, and q k = [φ, θ, ψ] represents the AUV attitude angles.Moreover, p t k = x t N , x t E , z t represents the sequence of the coordinates of the intersection points where the terrain is measured, x N represents the north coordinate in the map coordinate system, x E represents the east coordinate in the map coordinate system, δ Z represents the elevation difference between the measured terrain and the prior terrain map, ω the terrain height obtained at time , ℎ  represents the sequence of the interpolated terrain height at the positioning point  , and  represents the terrain measurement error.
Another structural mode of the TAN system is the tightly coupled mode, as shown in Figure 4.In this structure, the main navigation system outputs attitude information, velocity information, and other parameters to the TAN system.The TAN filter utilizes velocity information, attitude information, and terrain measurement information to perform positioning estimation.Furthermore, the main navigation system independently performs positioning estimation and outputs the positioning information separately, without interfering with the TAN positioning information.The state transition equations for this structure are shown in Equation (2), where the first equation describes the state transition of the main navigation system, and the second equation represents the state observation equation based on terrain matching.
where  and  represent the system state vectors obtained at times  and  − 1, respectively.In addition,  =  ,  ,  , ,  , , where  =  ,  ,  represents the coordinates of the AUV main navigation system, and  = , ,  represents the AUV attitude angles.Moreover,  =  ,  ,  represents the sequence of the coordinates of the intersection points where the terrain is measured,  represents the north coordinate in the map coordinate system,  represents the east coordinate in the map coordinate system,  represents the elevation difference between the measured terrain and the prior terrain map,  , and  , denote the angular velocities of the AUV around the yaxis and z-axis, respectively,   represents the rotation matrix formed by the AUV The tightly coupled navigation filtering mode is generally used for TAN estimation with low-cost sensor systems.For AUVs equipped with high-precision main navigation systems, the loosely coupled navigation filtering mode is typically employed.

TAP in Underwater TAN Systems
Underwater TAN systems utilize terrain feature matching and tracking for position estimation.Prior to providing accurate positioning information, the navigation coordinate system must be aligned with the spatial coordinate system of the terrain map.Subsequently, the TAN system performs continuous position estimation and outputs stable results by matching and tracking the underwater terrain features.In a previous study [5], the operation of underwater TAN systems was divided into two stages: initial positioning and tracking navigation.The primary objective of the initial positioning stage is to correct the accumulated time error of the main navigation system using high-precision TAP methods.This ensures the rapid convergence and stable output of the system within a short period and distance.However, the tracking positioning stage involves fusing the positioning information from the main navigation system with the TAP information, enabling continuous recursive estimation of the position and the provision of reliable positioning information.With the increasing diving depth and range of AUVs, larger initial positioning deviations can lead to challenges, such as slow convergence or even the divergence of the filtering process during the initial stage [3,4,48].Therefore, the issue of convergence during the initial stage of underwater TAN systems has attracted significant attention.Comparative studies [3,4] have investigated the effects of different initial positioning and filtering initialization methods on navigation accuracy.These studies have shown that convergence errors and speeds during the initial stage affect the final positioning accuracy.Moreover, filter convergence serves as the foundation for reliable system operation and the generation of valid positioning information.Thus, considering system reliability, the operational stages of underwater TAN systems should encompass the initial positioning and filter initialization, the navigation filter convergence stage, and the tracking navigation filter stage, as depicted in Figure 5.The Primary algorithms utilized in the three operational phases of TAN is shown in Table 2.
Moreover, filter convergence serves as the foundation for reliable system operation and the generation of valid positioning information.Thus, considering system reliability, the operational stages of underwater TAN systems should encompass the initial positioning and filter initialization, the navigation filter convergence stage, and the tracking navigation filter stage, as depicted in Figure 5.The Primary algorithms utilized in the three operational phases of TAN is shown in Table 2.

Main Algorithms Involved in TAP
Underwater TAN systems exhibit strong nonlinearity and non-Gaussian input error characteristics.The state equation of the system contains insolvable terms, making it challenging to apply existing state estimation methods to underwater TAN systems.Over the years, researchers have attempted to develop high-precision and reliable technologies for underwater TAN systems, focusing on such areas as initial positioning and system tracking filtering.The research achievements can be categorized into underwater TAP techniques, underwater TAN filtering techniques, and other cutting-edge issues in underwater TAN systems.In this section, the development of underwater TAN technologies is reviewed and explored based on these aspects.

TAP Estimation
Underwater TAP estimation involves approximating the probability distribution of an AUV at a search point based on the likelihood estimation of the RTM terrain and the prior terrain.This estimation can be considered the position points of the AUV and their probability distribution in the terrain space, serving as a primary source of positioning information for underwater TAN systems.During the initial stage of system operation, TAP estimation can be used to obtain the initial position points and probability distribution estimation of the system.In the tracking filtering stage, the filter combines information from TAP with the main navigation system to output TAN positioning information.Commonly used underwater TAP methods include point-cloud-based-and image-based matching methods.

Terrain Point-Cloud-Based Matching
TERCOM is the most commonly used batch-processing TAP method based on terrain elevation (Figure 6).This method quantifies the likelihood measure using the elevation difference between two terrain profiles and normalizes it into a probability distribution function using an exponential function [5].
where Z represents the elevation sequence of the measured terrain, Σ represents the error sequence of the measured terrain, H() represents the interpolation function for the terrain surface elevation, which can provide the elevation information at a given point based on its horizontal coordinates, and X s ij represents the coordinates of the search point within the search interval, where ij represents its index sequence number.The TERCOM algorithm can calculate the likelihood  of each search point corresponding to the prior terrain map using RTM data [49].By computing the maximum value of Equation ( 4), the position of the RTM in the DEM can be estimated.
Here, m and n represent the row and column indices of the measurement footprint sequence of the RTM, respectively, M and N represent the numbers of rows and columns in the measurement footprint sequence of the RTM, respectively,  represents the measured elevation of the footprint with index ,  , ℎ represents the terrain surface function, which describes the relationship between the horizontal coordinates and the elevation at each grid point on the terrain surface,  represents the horizontal coordinates corresponding to the index point ,  within the search interval  , and  represents the measurement error of the terrain elevation.
The mapping relationship between the RTM measurement points and the AUV position is unique.As shown in Figure 7, the RTM measurement point  obtained in the bathymetric sonar coordinate system O − x y z can be transformed to the vehicle coordinate system G − xyz through geometric relationships and then transformed to the AUV navigation coordinate system E − ξζη through the navigation system, which is represented by Equation (5).Therefore,  describes the probability distribution estimation of the AUV positioning points in the prior terrain map space.

𝑿
=      (5) where  represents the coordinate sequence of the RTM measurement footprints in the navigation coordinate system,  represents the coordinate vector of the RTM The TERCOM algorithm can calculate the likelihood L ij of each search point corresponding to the prior terrain map using RTM data [49].By computing the maximum value of Equation ( 4), the position of the RTM in the DEM can be estimated.
Here, m and n represent the row and column indices of the measurement footprint sequence of the RTM, respectively, M and N represent the numbers of rows and columns in the measurement footprint sequence of the RTM, respectively, z mn represents the measured elevation of the footprint with index (m, n), h( ) represents the terrain surface function, which describes the relationship between the horizontal coordinates and the elevation at each grid point on the terrain surface, X s ij represents the horizontal coordinates corresponding to the index point (i, j) within the search interval X s , and σ represents the measurement error of the terrain elevation.
The mapping relationship between the RTM measurement points and the AUV position is unique.As shown in Figure 7, the RTM measurement point c i obtained in the bathymetric sonar coordinate system O − x ′ y ′ z ′ can be transformed to the vehicle coordinate system G − xyz through geometric relationships and then transformed to the AUV navigation coordinate system E − ξζη through the navigation system, which is represented by Equation ( 5).Therefore, L ij describes the probability distribution estimation of the AUV positioning points in the prior terrain map space.
where X sonar mn represents the coordinate sequence of the RTM measurement footprints in the navigation coordinate system, X b mn represents the coordinate vector of the RTM measurement footprints in the sonar coordinate system, R b s represents the rotation matrix for transforming the RTM footprints from the sonar coordinate system to the vehicle coordinate system, X 0 represents the coordinate vector of the sonar installation position in the vehicle coordinate system, R n b represents the attitude rotation matrix for transforming from the vehicle coordinate system to the navigation coordinate system, and X n t represents the navigation output positioning at the measurement time.In one study [5], the characteristics of the TERCOM algorithm were investigated extensively.It was found that the likelihood function obtained from terrain-matching positioning using the TERCOM algorithm exhibits asymptotic Gaussian distribution properties.However, the TERCOM algorithm has been found to have issues with poor stability and reliability, which has led many researchers to investigate improvements to the algorithm.Li et al. addressed these issues by introducing geometric features as descriptors to construct a multivariate feature vector for terrain representation.This approach increased the amount of terrain-matching information, leading to improved positioning accuracy and stability [49].In addition to the TERCOM algorithm, the iterative closest point (ICP) algorithm and its variants, such as the improved convergence-constraint point (ICCP) algorithm, are commonly used for TAP based on point cloud data [62,63,85].Because terrain elevation data can be treated as point cloud data, many researchers have applied the ICP algorithm to estimate the TAP.The principle of the ICP algorithm is illustrated in Equation ( 6 where  represents the rotation matrix,  represents the translation matrix,  represents the error vector, and  represents the sum of the squared distance residuals between point sets  and  after rotation and translation.The optimal rotation and translation In one study [5], the characteristics of the TERCOM algorithm were investigated extensively.It was found that the likelihood function obtained from terrain-matching positioning using the TERCOM algorithm exhibits asymptotic Gaussian distribution properties.However, the TERCOM algorithm has been found to have issues with poor stability and reliability, which has led many researchers to investigate improvements to the algorithm.Li et al. addressed these issues by introducing geometric features as descriptors to construct a multivariate feature vector for terrain representation.This approach increased the amount of terrain-matching information, leading to improved positioning accuracy and stability [49].In addition to the TERCOM algorithm, the iterative closest point (ICP) algorithm and its variants, such as the improved convergence-constraint point (ICCP) algorithm, are commonly used for TAP based on point cloud data [62,63,85].Because terrain elevation data can be treated as point cloud data, many researchers have applied the ICP algorithm to estimate the TAP.The principle of the ICP algorithm is illustrated in Equation ( 6), where the point set X = {x 1 , . . . ,x N } represents the reference point set, and P = {p 1 , . . . ,p N } represents the point set to be aligned.Here, N represents the number of elements in the point set [86].Point set P can be aligned with point set X through rotation matrix R and translation matrix T.
   where R represents the rotation matrix, T represents the translation matrix, e i represents the error vector, and J represents the sum of the squared distance residuals between point sets P and X after rotation and translation.The optimal rotation and translation matrices are obtained when J is minimized.The ICP algorithm is widely used in outdoor mobile robot navigation for point cloud data matching and localization based on point cloud information [87][88][89].Therefore, the ICP algorithm and its improved versions have been widely applied in underwater TAP.A previously proposed ICCP algorithm [64] provides an improved method for matching and positioning accuracy in scenarios with significant initial errors.This method models the terrain distribution feature as a two-dimensional Gaussian distribution to reduce the interference of the heading error, and a terrain feature descriptor with rotation invariance is established.Other researchers [90,91] addressed the problem of filter divergence that inertial navigation systems face when encountering large track deviations.They improved the trajectory accuracy in ICCP using the difference between the matching result and the output of the inertial navigation system as the measurement value of the Kalman filter (KF).In another work [62], a multipath parallel ICCP algorithm was proposed, which considers the error differences in the measurement frames of the multibeam sonar, thus enhancing the stability of the ICCP algorithm.However, the ICP/ICCP algorithm is sensitive to initial deviations, may converge slowly, and is prone to local minima when dealing with large initial deviations, leading to mismatches.
TERCOM localization remains the most commonly used matching and localization method for terrain-matching navigation systems.The TERCOM algorithm is less sensitive than the ICP and ICCP methods to errors and exhibits higher robustness.In addition, TER-COM describes the similarity between terrains using Euclidean distance and an exponential function, and the likelihood function obtained by TERCOM exhibits asymptotic Gaussian distribution characteristics [5].These advantages provide significant convenience for the fusion of terrain-matching localization data and primary navigation localization data.

Image Matching Positioning
In addition to such methods as TERCOM, ICP, and ICCP, which utilize point cloud information for matching and positioning, underwater terrain matching and positioning algorithms based on image-matching techniques have been studied, including side-scan sonar images and optical camera images.In previous studies [92,93], matching and positioning algorithms based on side-scan images were investigated.They utilized log-polar transformation to solve the problem of estimating the rotation angle between real-time maps and reference images, thereby improving the matching accuracy and stability when there is an angular deviation in the measurement images.
In another work [94], a detailed comparison was made between the speeded-up robust features (SURF), oriented FAST and rotated BRIEF (ORB), binary robust invariant scalable keypoints (BRISK), and SURF-Harris algorithms, which were originally developed for optical image analysis.It was concluded that the further development of current feature description methodologies might be required for underwater acoustic image analysis.Other researchers [95] used an image attribute transfer algorithm and an advanced local feature descriptor to solve the problem of underwater acousto-optic image matching.The technique made it possible to preprocess multimodal images to obtain accurate matching results.
In one study [96], an underwater acoustic image simulation system was designed specifically for underwater acoustic image navigation.Building upon this system, research has been conducted on navigation systems based on forward-looking sonar images [97].The system utilizes deep-learning methods for image processing and registration, as well as for the estimation of the attitude of the AUV [98], demonstrating the advantages of deep-learning methods in underwater acoustic image matching, positioning, and navigation.In other research [99], a novel 3D preprocessing filter, additions to line and line segment detection methods, a system description using concepts from the theory of mul-tilevel hierarchical systems, and extensive evaluation tests using synthetic and real data were proposed.
Although the statistical moment information of the images enhances the noise resistance of the system to a certain extent, it reduces the dimensionality of the terrain elevation information, thereby weakening the amount of matching information.In the context of matching and positioning based on optical images, researchers [100][101][102] have investigated techniques for the simultaneous matching, stitching, and localization of images.These techniques make it possible to expand and map terrains in unknown underwater environments.The research conducted in these studies is important for achieving fine-grained measurements of unknown underwater topography and the high-precision positioning of AUVs.These technologies are commonly employed for precise measurements in such fields as underwater archaeology, underwater geological surveys, and underwater ecosystem assessments, as well as for the high-precision positioning of underwater robots [103][104][105].In addition, image feature-based terrain matching and positioning techniques are widely applied in planetary rover landing positioning and navigation [106][107][108][109].

Iterative Filtering Methods for TAN
The filter serves as the data fusion center in TAN systems, where it combines information from the TAP and primary navigation positioning at the time step of the navigation solution.The TAP error is bounded in the spatial and temporal domains, whereas the propagated navigation error accumulates indefinitely.However, there is a strong correlation between the navigation points.Underwater TAN filtering algorithms integrate the boundedness of the TAP error and the strong correlation of the propagated navigation positioning, thereby harnessing the advantages of both positioning information types and obtaining more stable and accurate positioning information.Filtering methods for TAN systems have become an important research direction for underwater TAN technology.With the development of underwater mapping technology, computer hardware technology, and nonlinear estimation theory, filtering algorithms for TAN systems have gone through three major stages: Kalman, nonparametric, and robust filtering.Particle filtering has become the primary filtering method for underwater terrain-matching navigation.

Kalman Filter and Its Improvements
The initial iterative filtering estimation method is Kalman filtering [110].However, the KF is a linearized filter that requires linearization of the terrain during the filtering processspecifically, linearizing the H( ) term in the state observation equation of the underwater TAN system and estimating the terrain slope [111].Owing to the strong nonlinearity of the terrain, significant linearization errors are introduced, leading to filter divergence.Some researchers [73,112] have studied a TAN method based on the innovation transformation Kalman filtering technique, which mitigated the adverse effects of linearization errors to some extent.Comparative experiments with a particle filter (PF) and point mass filter showed that the PMF achieved the highest positioning accuracy when the measurement information was limited.As the measurement information increased, the PF approached the accuracy of the PMF, but both outperformed the KF and unscented KF in terms of positioning accuracy.Owing to the strong nonlinearity of the terrain and the non-Gaussian nature of the measurement errors, the Kalman filtering algorithm struggles to adapt to TAN systems with strong nonlinearity and non-Gaussian characteristics.PF and PMF methods, based on posterior Bayesian estimation theory and their direct numerical solution techniques, offer higher robustness and adaptability, making them particularly suitable for addressing the positioning estimation problem in TAN systems.Consequently, they have gradually replaced Kalman filtering methods as the most commonly used methods for state estimation in TAN systems.

Posterior Bayesian Estimation and Its Numerical Solution
The posterior Bayesian estimation equation for the underwater TAN states, as shown in Equation (7), consists of the state transition equation of the primary navigation system as the first term and the posterior Bayesian estimation equation of the system state as the second term.Solving the equation p( x t |Y t ) requires an evaluation of the integral term R 2 p( y t |x t )p( x t |Y t−1 )dx t .However, p( y t |x t ) does not have an analytical form, mainly because the mapping relationship between the AUV positioning and the terrain surface cannot be represented analytically.Although under the assumption of linearization and Gaussian error distribution, the Kalman filtering equations can be derived through derivation [5], the instability of the system caused by terrain linearization errors and non-Gaussian errors has become an insurmountable bottleneck.PF and PMF methods based on numerical solution techniques have been proposed to solve the state estimation problem of nonlinear dynamic systems [68,69,74,113].The PMF and PF approximate the probability distribution function p( y t |x t ) of the TAN positioning using finite grids and particle sets, making pos- sible a direct numerical solution of the posterior Bayesian estimation of the TAN system positioning probability distribution in Equation ( 2) and significantly reducing the influence of terrain nonlinearity.Numerical solution methods based on Bayesian theory, such as PMF and PF, significantly improve filtering accuracy and stability [114], thereby advancing the development and application of underwater TAN technology.
where p( x t |Y t−1 ) represents the predicted positioning probability distribution of the under- water TAN system at time t, p( x t |Y t ) represents the posterior corrected probability of TAP in the underwater system, α t = R 2 p( y t |x t )p( x t |Y t−1 )dx t is the normalization parameter related to the terrain measurement error, p(x t − u t−1 − x t−1 ) represents the probability distribution of the propagated navigation system estimated position (where u t−1 is the system state input information from time t − 1 to t), p( x t−1 |Y t−1 ) is the system state at time t − 1, R 2 dx t−1 represents the integration range for solving the probability distribution of underwater TAN positioning, and p( y t |x t ) represents the conditional distribution probabil- ity of observation value y t obtained by x t at time t.Equations ( 7) and ( 8) make it possible to estimate the TAP points.
Similarly, the estimation of the underwater TAP error can be obtained using Equation ( 9): (1) Particle Filter The FFI was one of the earliest research institutions to apply TAN technology to underwater navigation [46,47,115,116].In the 1980s, FFI researchers introduced nonparametric filtering methods for state iteration estimation in TAN systems.The FFI also developed the TERLAB simulation platform specifically for the research on and application of underwater TAN technology [116].Subsequently, an increasing number of researchers have studied underwater TAN technology based on PF methods.PFs approximate probability distributions with randomly distributed particles that convert an integral operation into a summation operation (Figure 8).Therefore, it can be very convenient to solve non-Gaussian problems and integral operations without analytic representations.
The operation of a PF typically involves three steps: prediction, updating, and resampling.During the prediction step, the particle set is moved to the next navigation point based on the primary navigation information.In the updating step, the particle weights are adjusted based on the probability distribution function of TAP.The resampling step involves adjusting the distribution and weights of particles based on their weights to prevent particle degeneration.Strong nonlinearity and non-Gaussian input errors in TAN systems often lead to significant discrepancies between the true and proposed particle distributions.Furthermore, the problem of pseudopeaks in the probability distribution functions caused by terrain self-similarity further affects the stability and accuracy of the PF.Several studies [75,117,118] have addressed the issues of non-Gaussian and multimodal distributions in particle filtering owing to terrain self-similarity by proposing PF algorithms based on regularized representations and clustering analysis.To address the problem of intractable state prediction in TAN systems caused by observation equations, a Gaussian process PF method based on historical observation data and Gaussian models was proposed [119] for terrain prediction in underwater vehicle motion.A data-driven particle-filtering framework [75] was also proposed to avoid the explicit modeling of system disturbances.Another approach [76] addressed the non-Gaussian nature of the particle distribution by approximating the Bayesian function using multiple Gaussian components.To mitigate the computational burden caused by an increasing number of particles, an edge PF based on a more computationally efficient marginalized particle filter (MPF) was applied to state estimation in TAN systems [70].A combined MPF that integrates a KF and a PF was proposed [81] to reduce the computational load of the PF process.An improved PF method based on the Huber function for modifying the particle weights was presented [77].Furthermore, to address the problem of an explosive increase in particle quantity during the initial stage, an adaptive PF TAN filtering method utilizing the Kullback-Leibler distance measure was introduced [78].A three-dimensional marginalized PF was proposed to enhance the filtering robustness and reduce the computational load in underwater TAN under uncertain conditions of tidal depth deviation [66,77].An initialization sampling method for a PF based on constraints from TAP intervals was proposed [5] to reduce the particle spread range and quantity and improve the computational efficiency and convergence speed.A PF-based TAN method [120] that incorporates a nonrigid terrain transformation model that can correct nonrigid terrain transformation errors caused by the absence of a sound velocity profile and eliminate its impact on TAN was proposed.Because of its advantages in handling nonlinear and non-Gaussian noise in dynamic system state estimation problems, the PF algorithm has become the primary estimation method in TAN.The operation of a PF typically involves three steps: prediction, updating, and resampling.During the prediction step, the particle set is moved to the next navigation point based on the primary navigation information.In the updating step, the particle weights are adjusted based on the probability distribution function of TAP.The resampling step involves adjusting the distribution and weights of particles based on their weights to prevent particle degeneration.Strong nonlinearity and non-Gaussian input errors in TAN systems often lead to significant discrepancies between the true and proposed particle distributions.Furthermore, the problem of pseudopeaks in the probability distribution functions caused by terrain self-similarity further affects the stability and accuracy (2) Marginalized Particle Filter Similar to the PF algorithm, the MPF is a numerical method for solving posterior Bayesian estimation.It employs a finite grid approximation of the probability distribution function and is well suited for nonlinear and non-Gaussian state estimation in dynamic systems.In recent years, it has been widely applied to TAN.The MPF algorithm was initially proposed by Bucy and first applied to underwater TAN by Bergman [114,121].As shown in Figure 9, the MPF algorithm divides the integral continuous probability function surface into volumes by using rectangular prisms with a base length of δ and a height equal to the probability value at that point.This transforms the integration of the continuous probability function surface into a volume summation problem.A comparative study of the MPF and PF methods in TAN state estimation was conducted [51,122], demonstrating that the MPF has a higher computational accuracy, but the computational complexity of the algorithm increases rapidly with the dimensionality of the state.However, although the PF has a lower accuracy than the MPF, it exhibits a higher computational efficiency when dealing with high-dimensional problems.
nonrigid terrain transformation model that can correct nonrigid terrain transformation errors caused by the absence of a sound velocity profile and eliminate its impact on TAN was proposed.Because of its advantages in handling nonlinear and non-Gaussian noise in dynamic system state estimation problems, the PF algorithm has become the primary estimation method in TAN.
(2) Marginalized Particle Filter Similar to the PF algorithm, the MPF is a numerical method for solving posterior Bayesian estimation.It employs a finite grid approximation of the probability distribution function and is well suited for nonlinear and non-Gaussian state estimation in dynamic systems.In recent years, it has been widely applied to TAN.The MPF algorithm was initially proposed by Bucy and first applied to underwater TAN by Bergman [114,121].As shown in Figure 9, the MPF algorithm divides the integral continuous probability function surface into volumes by using rectangular prisms with a base length of δ and a height equal to the probability value at that point.This transforms the integration of the continuous probability function surface into a volume summation problem.A comparative study of the MPF and PF methods in TAN state estimation was conducted [51,122], demonstrating that the MPF has a higher computational accuracy, but the computational complexity of the algorithm increases rapidly with the dimensionality of the state.However, although the PF has a lower accuracy than the MPF, it exhibits a higher computational efficiency when dealing with high-dimensional problems.Because a standard grid is used in the implementation of the PMF algorithm, the grid quality can affect the estimation accuracy.The adaptive grid PMF algorithm was proposed [123,124] to achieve higher accuracy with a lower computational load.In another study [66], the application performance of the PMF algorithm under complex conditions was investigated, and it was found that the algorithm is significantly affected by measurement depth errors and initial position deviations.In other work [66], a marginalized PMF in three dimensions was developed to concurrently estimate and compensate for tidal depth bias.In this method, the tidal depth bias is extended as a state variable and estimated using a KF, whereas the horizontal position state is estimated using the original two-dimensional point mass filter.Using multibeam sonar, simulation experiments in a real Because a standard grid is used in the implementation of the PMF algorithm, the grid quality can affect the estimation accuracy.The adaptive grid PMF algorithm was proposed [123,124] to achieve higher accuracy with a lower computational load.In another study [66], the application performance of the PMF algorithm under complex conditions was investigated, and it was found that the algorithm is significantly affected by measurement depth errors and initial position deviations.In other work [66], a marginalized PMF in three dimensions was developed to concurrently estimate and compensate for tidal depth bias.In this method, the tidal depth bias is extended as a state variable and estimated using a KF, whereas the horizontal position state is estimated using the original two-dimensional point mass filter.Using multibeam sonar, simulation experiments in a real DEM demonstrated that the method can accurately estimate the tidal depth bias and obtain a robust navigation solution in suitable terrain.
Although the PF-and PF-based filters have significant advantages in the state estimation of TAN systems, strong nonlinearity and non-Gaussian input errors still lead to filter instability.Currently, the techniques for measuring and reconstructing underwater terrain maps primarily include gravity field inversion, acoustic measurement, and optical measurement.

Other Relevant Technological Advancements
The practical application of underwater TAN technology requires high-precision prior terrain maps, which makes the practical application of underwater TAN difficult.However, the greatest bottleneck affecting the practical application of underwater TAN technology is the insufficient reliability of underwater TAN systems.The main application prospects of the underwater TAN technology lie in long-term, long-distance detection and development in deep-sea and remote areas.Because of the high investment, difficulty, and risk associated with deep-sea exploration and development, the autonomous equipment used for deepsea exploration must possess high reliability to ensure the safety and smoothness of the detection equipment and processes.This requires underwater navigation systems with sufficiently high reliability.However, underwater terrain-aided navigation systems face several challenges and unresolved difficulties.

Prior Digital Elevation Model
Obtaining high-precision prior DEMs is crucial for effective and reliable positioning information of underwater TAN systems.However, the acquisition of underwater terrain, particularly large-scale high-precision bathymetric maps in deep and remote ocean areas, remains a global challenge.Currently, the techniques used for bathymetric map measurement and reconstruction primarily include gravity field inversion, acoustic ranging, and optical imaging.
Researchers [71,125,126] have proposed methods for obtaining bathymetric maps based on gravity field inversion, which can achieve global bathymetric inversion.However, the accuracy and resolution of the inversion results based on the gravity field are low, making them unsuitable for underwater TAN systems.High-precision bathymetric detection mainly relies on acoustic sounding techniques, such as single-beam or multibeam echo sounding.Based on the carrier of the acoustic sensor, these techniques can be categorized as shipborne, towed, or underwater-robot sounding systems.Currently, most global bathymetric data are derived from shipborne sounding systems.However, less than 20% of the seafloor has been adequately surveyed.
In 2017, the Japan Society for the Promotion of Science and the General Bathymetric Chart of the Oceans jointly discussed how marine mapping can support the United Nations' sustainable development goals and proposed the Seabed 2030 project, which aims to achieve 100% global seafloor coverage by 2030.As of 2023, 24.9% of the global bathymetric mapping was completed [127].However, high-precision and high-efficiency deep-sea bathymetric mapping remains a bottleneck that must be overcome.
With the rapid development of underwater robotics technology, an increasing number of research institutions have focused on developing deep-sea terrain mapping techniques based on underwater robotic platforms.These techniques show promise for improving the accuracy and efficiency of deep-sea terrain mapping.In one study [128], a technique for underwater terrain mapping based on an AUV platform was devised.This technique utilizes a multibeam sonar system mounted on an AUV to measure seabed profiles.The obtained profile information is then integrated with the AUV navigation system to compute the terrain profiles in the geodetic coordinate system, making autonomous underwater terrain mapping possible by using an AUV.However, because of the accumulation of errors in the navigation coordinate system of an AUV over time, it is challenging to achieve long-term and large-scale terrain mapping.In another study [129], a methodology was proposed for processing 3D multibeam sonar big data based on the stepwise processing of a dataset with 3D models and isoline map generation that makes fast processing possible, and the obtained DEMs were of good quality.
To address the issue of cumulative errors in AUV navigation systems over time, researchers worldwide have focused on underwater simultaneous localization and mapping (SLAM) technology.This technology aims to enable AUVs to construct global terrain maps autonomously and perform navigation and localization in unknown underwater environments without external assistance [130][131][132].This includes bathymetric SLAM (BSLAM) based on multibeam bathymetry information [12,16,17,[133][134][135], sensor-based acoustic and optical image SLAM utilizing acoustic and optical image data [136][137][138][139][140], and SLAM techniques based on structured light point cloud matching [133,141].In BSLAM, multibeam bathymetry systems are commonly used to acquire three-dimensional point cloud information of the underwater terrain for matching and mapping.Image-based SLAM techniques utilize the image information obtained from underwater optical or acoustic devices for matching and mapping.Point cloud SLAM techniques employ stereo cameras or line-structured light sensors for data acquisition.Figure 10 shows the seabed terrain data obtained using multibeam sonar, an optical camera, and structured light.
to achieve long-term and large-scale terrain mapping.In another study [129], a methodology was proposed for processing 3D multibeam sonar big data based on the stepwise processing of a dataset with 3D models and isoline map generation that makes fast processing possible, and the obtained DEMs were of good quality.
To address the issue of cumulative errors in AUV navigation systems over time, researchers worldwide have focused on underwater simultaneous localization and mapping (SLAM) technology.This technology aims to enable AUVs to construct global terrain maps autonomously and perform navigation and localization in unknown underwater environments without external assistance [130][131][132].This includes bathymetric SLAM (BSLAM) based on multibeam bathymetry information [12,16,17,[133][134][135], sensor-based acoustic and optical image SLAM utilizing acoustic and optical image data [136][137][138][139][140], and SLAM techniques based on structured light point cloud matching [133,141].In BSLAM, multibeam bathymetry systems are commonly used to acquire three-dimensional point cloud information of the underwater terrain for matching and mapping.Image-based SLAM techniques utilize the image information obtained from underwater optical or acoustic devices for matching and mapping.Point cloud SLAM techniques employ stereo cameras or line-structured light sensors for data acquisition.Figure 10 shows the seabed terrain data obtained using multibeam sonar, an optical camera, and structured light.The comparison in Figure 10 reveals that only the camera and multiline structured light [143] can provide information about a planar region in a single measurement, whereas the multibeam sonar and single-line structured light can only capture a single-line measurement.Because multibeam bathymetry can only obtain a strip-shaped terrain in a single measurement, it cannot construct continuous matching and correlation information.Therefore, submap partitioning is commonly used for map segmentation and feature matching [12,131].Figure 11 illustrates the typical workflow for multibeam bathymetric SLAM.In contrast to single-frame multibeam bathymetry, a single-frame optical image can cover a limited planar region, making it possible to construct matching correlations between adjacent measurement frames.Therefore, there are significant differences in the algorithm design between bathymetric SLAM and image-based SLAM.However, regardless of the mapping approach, SLAM techniques rely on matching the correlations of environmental features to establish equations for AUV pose estimation.These equations serve as constraints related to AUV pose and feature matching, and they form the optimization equations for global mapping and AUV pose estimation.
In the equations, f (X i−1 , u i ) represents the state transition equation with an input covariance matrix Q.In addition, h(X k , X l ) represents the measurement equation with a measurement error of Γ k,l , where (k, l) represents the indices of two different measure- ment frames that form a loop closure.Moreover, LO represents the set of indices of all measurement frames that form loop closures, z k,l represents the actual terrain elevation measurement sequence, X i−1 represents the position information of the AUV, X * represents the optimized trajectory of the AUV, and M represents the length of the AUV state sequence.lations between adjacent measurement frames.Therefore, there are significant differences in the algorithm design between bathymetric SLAM and image-based SLAM.However, regardless of the mapping approach, SLAM techniques rely on matching the correlations of environmental features to establish equations for AUV pose estimation.These equations serve as constraints related to AUV pose and feature matching, and they form the optimization equations for global mapping and AUV pose estimation.
In the equations,   ,  represents the state transition equation with an input covariance matrix .In addition, ℎ  ,  represents the measurement equation with a measurement error of  , , where ,  represents the indices of two different Due to the large measurement range and relatively low measurement accuracy of acoustic sensors, the SLAM technology is mainly used in the initial stage of seabed surveys, for large-scale seabed topography and geomorphology surveys, to obtain information on the large-scale changes in the seabed elevation and seabed substrate.Although SLAM technology solves the problem of unknown underwater terrain detection, the process of measuring the large-scale underwater terrain still faces the challenges of a low efficiency and high cost, owing to the low speed and limited range of AUVs.To address the efficiency issue, underwater terrain measurement techniques based on multi-robot collaboration have been rapidly developed [143][144][145][146].In one study [147], large-scale seabed mapping using multiple AUVs was developed based on a graph-based cooperative bathymetric SLAM system that can compress many bathymetric measurements into small-scale acoustic packets and yield accurate navigation results with a 10% loss in acoustic packets caused by unreliable acoustic communication.In another work [148], the problem of multi-AUV collaborative online target detection was studied, and the simulation results showed that the algorithm can improve the detection efficiency by at least 40% compared with a single AUV.In another study [149], a path-planning method for a multi-AUV TAN was investigated by considering constraints, such as collision avoidance, arrival times, energy minimization, and acoustic communication between AUVs.A multi-AUV collaborative TAN technology [150] was studied, and it was verified through simulations that the proposed collaborative localization technique could extend the AUV operating time.Underwater SLAM technology was investigated [151] based on multi-robot collaboration, improving the terrain measurement efficiency and mapping accuracy.In other research [152], a multi-robot collaborative measurement method suitable for larger-scale underwater terrain measurements was developed.As shown in Figure 12, Robots A and B are both equipped with acoustic positioning devices.During the measurement process, Robot A is first submerged in a predetermined area and remains stationary as an acoustic positioning reference.Robot B, with Robot A as the positioning reference, performs terrain measurements in a certain seabed area.Once the measurement area is covered, Robot B sub-merges to the predetermined area and remains stationary as the new positioning reference, and Robot A, with Robot B as the positioning reference, performs terrain measurement in another seabed area.This cycle is repeated, theoretically allowing for the measurement of the entire underwater terrain.Compared to SLAM technology based on acoustic measurement sensors, SLAM technology based on optical information is more mature and has been widely applied in autonomous vehicles, drones, and outdoor mobile robots [153].While structured light and stereo vision measurements have a limited measurement range, they can perform closerange observation with high accuracy, making them suitable for small-scale and high-precision measurement tasks.They have been widely adopted in underwater archaeology for observation [154,155], the detection of underwater sediments or structures [156,157], and other underwater tasks.Reference [156] uses 3D visual reconstruction methods to study the sea floor hydrothermal systems and identify the habitats of Shinkaia crosnieri squat lobsters and their population density.Reference [36] presents an extensive review of the sensors and the methodologies used in archaeological underwater 3D recording and mapping, together with relevant highlights of well-renowned projects in 3D recording underwater.Furthermore, the literature [158][159][160] has investigated the common issue of light refraction correction in underwater optical measurement and SLAM.In reference [101,161], visual SLAM with a focus on a stereo camera system is presented to estimate the motion of autonomous underwater vehicles (AUVs) and build the feature map of the surrounding environment in real time.Reference [105] details the operations, discusses the results of the ancient shipwreck survey, and identifies the specific challenges of adapting AUV technology for deep water archaeology.

Terrain Adaptability Analysis and Path Planning
Underwater TAN technology can provide AUVs with positional information that has bounded errors.In theory, as long as the terrain adaptability meets the accuracy requirements of matching navigation, it can provide sufficiently accurate position information for AUVs over a long period of time.However, flat regions are likely to exist in underwater terrain maps, which have low adaptability and result in poor matching localization accuracy.This is particularly true for long-range AUV TAN, where large flat areas may exist in the extensive prior terrain map.Once the AUV enters such an area, there is a risk of filter divergence.Furthermore, AUVs have limited onboard energy; therefore, it is necessary to select the optimal path considering such factors as the mission, energy, and navigation accuracy, while ensuring that the energy reserves are sufficient to meet the consumption during the planned path tracking.Therefore, navigation path planning is a crucial research area in underwater TAN systems.Research on path planning for underwater TAN has mainly focused on terrain adaptability analysis, adaptability quantification, and Compared to SLAM technology based on acoustic measurement sensors, SLAM technology based on optical information is more mature and has been widely applied in autonomous vehicles, drones, and outdoor mobile robots [153].While structured light and stereo vision measurements have a limited measurement range, they can perform close-range observation with high accuracy, making them suitable for small-scale and highprecision measurement tasks.They have been widely adopted in underwater archaeology for observation [154,155], the detection of underwater sediments or structures [156,157], and other underwater tasks.Reference [156] uses 3D visual reconstruction methods to study the sea floor hydrothermal systems and identify the habitats of Shinkaia crosnieri squat lobsters and their population density.Reference [36] presents an extensive review of the sensors and the methodologies used in archaeological underwater 3D recording and mapping, together with relevant highlights of well-renowned projects in 3D recording underwater.Furthermore, the literature [158][159][160] has investigated the common issue of light refraction correction in underwater optical measurement and SLAM.In reference [101,161], visual SLAM with a focus on a stereo camera system is presented to estimate the motion of autonomous underwater vehicles (AUVs) and build the feature map of the surrounding environment in real time.Reference [105] details the operations, discusses the results of the ancient shipwreck survey, and identifies the specific challenges of adapting AUV technology for deep water archaeology.

Terrain Adaptability Analysis and Path Planning
Underwater TAN technology can provide AUVs with positional information that has bounded errors.In theory, as long as the terrain adaptability meets the accuracy requirements of matching navigation, it can provide sufficiently accurate position information for AUVs over a long period of time.However, flat regions are likely to exist in underwater terrain maps, which have low adaptability and result in poor matching localization accuracy.This is particularly true for long-range AUV TAN, where large flat areas may exist in the extensive prior terrain map.Once the AUV enters such an area, there is a risk of filter divergence.Furthermore, AUVs have limited onboard energy; therefore, it is necessary to select the optimal path considering such factors as the mission, energy, and navigation accuracy, while ensuring that the energy reserves are sufficient to meet the consumption during the planned path tracking.Therefore, navigation path planning is a crucial research area in underwater TAN systems.Research on path planning for underwater TAN has mainly focused on terrain adaptability analysis, adaptability quantification, and path-planning algorithms.Adaptability analysis involves analyzing and parameterizing terrain features while identifying the most important terrain feature parameters that affect the matching accuracy.In one study [54], directional feature factors that affect the matching accuracy based on error analysis were derived, showing that the positioning accuracy of a point in the terrain map varies in different directions.Other research [162] addressed the selection of terrain feature parameters in adaptability analysis by constructing a binary logistic regression classifier based on factor analysis results, confirming the rationality of the factor interpretation and the effectiveness of the analysis results.
Adaptability quantification aims to evaluate the adaptability of the terrain quantitatively based on the terrain feature parameters and to divide the terrain map optimally into highly and poorly adaptable regions.In a previous study [54], an optimal partitioning method for terrain adaptability based on a terrain spatial standard grid was investigated.This method divides the terrain space using standard grids and obtains the optimal gridbased terrain map by maximizing the allocation of highly and poorly adaptable regions to different grids [54].Figure 13 shows the grid-based results of the prior terrain map obtained using the method proposed [54].The grid is determined based on the principle of optimal segmentation according to adaptability, in which highly and poorly adaptable regions in the terrain map are optimally allocated to different grids.Figure 14 shows the quantified values of terrain adaptability in each grid under the optimal grid partitioning.path-planning algorithms.Adaptability analysis involves analyzing and parameterizing terrain features while identifying the most important terrain feature parameters that affect the matching accuracy.In one study [54], directional feature factors that affect the matching accuracy based on error analysis were derived, showing that the positioning accuracy of a point in the terrain map varies in different directions.Other research [162] addressed the selection of terrain feature parameters in adaptability analysis by constructing a binary logistic regression classifier based on factor analysis results, confirming the rationality of the factor interpretation and the effectiveness of the analysis results.
Adaptability quantification aims to evaluate the adaptability of the terrain quantitatively based on the terrain feature parameters and to divide the terrain map optimally into highly and poorly adaptable regions.In a previous study [54], an optimal partitioning method for terrain adaptability based on a terrain spatial standard grid was investigated.This method divides the terrain space using standard grids and obtains the optimal gridbased terrain map by maximizing the allocation of highly and poorly adaptable regions to different grids [54].Figure 13 shows the grid-based results of the prior terrain map obtained using the method proposed [54].The grid is determined based on the principle of optimal segmentation according to adaptability, in which highly and poorly adaptable regions in the terrain map are optimally allocated to different grids.Figure 14 shows the quantified values of terrain adaptability in each grid under the optimal grid partitioning.Path-planning algorithms primarily focus on finding optimal navigation paths considering multiple constraints, such as terrain adaptability, distance, and energy consumption.In one study [9], the navigation path-planning problem of guiding AUVs through the Arctic ice cover using underwater terrain-matching navigation was investigated.Building upon this, researchers [28] further considered such constraints as Arctic Ocean current path-planning algorithms.Adaptability analysis involves analyzing and parameterizing terrain features while identifying the most important terrain feature parameters that affect the matching accuracy.In one study [54], directional feature factors that affect the matching accuracy based on error analysis were derived, showing that the positioning accuracy of a point in the terrain map varies in different directions.Other research [162] addressed the selection of terrain feature parameters in adaptability analysis by constructing a binary logistic regression classifier based on factor analysis results, confirming the rationality of the factor interpretation and the effectiveness of the analysis results.
Adaptability quantification aims to evaluate the adaptability of the terrain quantitatively based on the terrain feature parameters and to divide the terrain map optimally into highly and poorly adaptable regions.In a previous study [54], an optimal partitioning method for terrain adaptability based on a terrain spatial standard grid was investigated.This method divides the terrain space using standard grids and obtains the optimal gridbased terrain map by maximizing the allocation of highly and poorly adaptable regions to different grids [54].Figure 13 shows the grid-based results of the prior terrain map obtained using the method proposed [54].The grid is determined based on the principle of optimal segmentation according to adaptability, in which highly and poorly adaptable regions in the terrain map are optimally allocated to different grids.Figure 14 shows the quantified values of terrain adaptability in each grid under the optimal grid partitioning.Path-planning algorithms primarily focus on finding optimal navigation paths considering multiple constraints, such as terrain adaptability, distance, and energy consumption.In one study [9], the navigation path-planning problem of guiding AUVs through the Arctic ice cover using underwater terrain-matching navigation was investigated.Building upon this, researchers [28] further considered such constraints as Arctic Ocean current Path-planning algorithms primarily focus on finding optimal navigation paths considering multiple constraints, such as terrain adaptability, distance, and energy consumption.In one study [9], the navigation path-planning problem of guiding AUVs through the Arctic ice cover using underwater terrain-matching navigation was investigated.Building upon this, researchers [28] further considered such constraints as Arctic Ocean current models, high-latitude navigation estimation errors, and terrain depth measurement errors.They employed the optimized rapidly exploring random tree algorithm to obtain the optimal navigation path and conducted simulations using the Autosub Long-Range 1500 (ALR1500) AUV platform.
Currently, underwater TAN path planning only considers constraints in the horizontal plane and does not consider the vertical motion of the AUV and constraints related to vertical terrain variations.In practice, it is necessary to consider the vertical height of the AUV during underwater TAN path planning.The vertical height of the AUV from the seabed affects the range and accuracy of terrain detection, which directly impacts the matching localization accuracy.Moreover, the AUV must adjust its navigation depth based on terrain variations, affects energy consumption.Moreover, the AUV itself has constraints on the navigation depth.Therefore, underwater TAN path planning must consider not only adaptability constraints in the horizontal plane but also depth constraints for path points and even differential constraints for the path.

Initial Positioning and Accuracy Evaluation
Previous studies have mostly focused on cases with small initial localization errors, where the use of TERCOM for initial localization can provide high-precision initial positions and PF initialization results.However, as AUV missions extend to deep-sea and longrange submerged navigation, navigation systems experience significant error accumulation because deep diving and extended travel distances are necessary.In such cases, the pseudopeaks caused by the initial localization errors can lead to severe filter initialization errors and affect the stability and convergence of the filter.
To address the high-precision initial localization requirement in underwater TAN systems, a method [65] was proposed called the "terrain point diversity consistency test" (TPDCT) based on multiple TERCOM localization points to determine the effectiveness of TAP points.In other work [163], this method was applied to geomagnetic matching navigation, and practical experiments were conducted.However, underwater terrain acquisition is challenging, and the region suitable for matching positioning is limited.As a passive filtering method, TPDCT requires a considerable number of TAP points.In addition, the errors in the TAP points are related to the terrain features, and the probability distribution of different TAP points varies, resulting in low data utilization and poor applicability in underwater environments.Therefore, a multiple transmission-reception point fusion positioning (MTFP) positioning method [3] was studied based on data fusion estimation and trajectory validity determination.The MTFP and TPDCT methods use multiple TERCOM points for positioning estimation.Unlike the TPDCT method, the MTFP method actively uses data to obtain positioning estimates for trajectory segments.Based on this, the estimation cost of each trajectory segment is calculated, and the optimal trajectory is selected based on the lowest cost condition, thereby transforming the positioning estimation problem into a local trajectory segment estimation problem [48].Furthermore, to address the problem of excessive pseudo-peaks in TAP causing a high computational load, researchers [39,164] proposed a PF initialization method based on confidence interval constraints.By estimating the confidence interval of the initial localization probability distribution function and using it as the distribution boundary for initializing the particle set, the number and range of initial particles are reduced, improving the convergence speed.
Regarding the accuracy evaluation of underwater TAP, most studies have used the Fisher information matrix as the lower bound for the error distribution of TAP [5,35,58].However, the error distribution of the TAP estimates differs fundamentally from that of the recursive positioning estimates.The confidence interval is not elliptical.In one study [49], it was demonstrated that TAP points have distinct directional distribution characteristics, particularly in areas with significant terrain directionality.A stochastic jump model was proposed for the confidence interval estimation of TAP points (Figure 15), and a new equation was derived for the confidence interval estimation of TAP points, significantly improving the accuracy of the confidence interval estimation for TAP points: Here,   represents the sum of squared residuals of the TAP positions,  denotes the variance in terrain measurement errors, and 1 −  represents the confidence level, typically in the range of 0.95,0.97 .The TAP points are randomly fluctuating points associated with terrain features in the spatial domain, and their probability distribution does not have a unique central peak but multiple pseudo-peaks, as shown in Figure 16.Therefore, the confidence intervals for TAP points cannot be estimated using the traditional elliptical confidence interval theory.Other position confidence interval estimation methods of TAN in different application scenarios list in Table 3.Here, S(X p ) represents the sum of squared residuals of the TAP positions, (σ p ) 2 denotes the variance in terrain measurement errors, and 1 − α represents the confidence level, typically in the range of [0.95, 0.97].The TAP points are randomly fluctuating points associated with terrain features in the spatial domain, and their probability distribution does not have a unique central peak but multiple pseudo-peaks, as shown in Figure 16.Therefore, the confidence intervals for TAP points cannot be estimated using the traditional elliptical confidence interval theory.Other position confidence interval estimation methods of TAN in different application scenarios list in Table 3.
Here,   represents the sum of squared residuals of the TAP positions,  denotes the variance in terrain measurement errors, and 1 −  represents the confidence level, typically in the range of 0.95,0.97 .The TAP points are randomly fluctuating points associated with terrain features in the spatial domain, and their probability distribution does not have a unique central peak but multiple pseudo-peaks, as shown in Figure 16.Therefore, the confidence intervals for TAP points cannot be estimated using the traditional elliptical confidence interval theory.Other position confidence interval estimation methods of TAN in different application scenarios list in Table 3.   Table 3.Primary method for representing confidence intervals in statistics.

Confidence Interval Estimation Equation Estimation
Theory Remarks References where σ 2 e represents the variance in terrain measurement errors, and where L 1−α represents the lower bound of the likelihood function value of the TAP points at a confidence level of 1 − α, σ p represents the standard deviation of terrain measurement errors, and S(X)

Simulation of an Underwater TAN System
Underwater TAN technology is primarily used for long-duration, long-distance navigation and positioning in the deep sea and remote underwater spaces.The deep-sea environment is a communication-denied space, where it is difficult for human intervention to occur quickly and effectively during operational processes.Therefore, underwater TAN systems must be sufficiently stable and reliable.Field trials and simulation experiments are important means of testing system performance.However, conducting field trials in deep-sea environments is challenging and risky.Therefore, shallow-water tests are often used as substitutes.This leads to significant differences between the experimental environment and the actual operating conditions, and the terrain features may also be limited, resulting in limited practical reference values for actual experiments.Simulation experiments offer such advantages as a low risk, high efficiency, and safety control, and they have been widely adopted in research on underwater autonomous systems.In recent years, they have been extensively used in underwater TAN research.
Obtaining system input information with high fidelity is essential to simulating underwater TAN systems, particularly when simulating input errors, which must have probabilistic completeness.This means that the input errors must encompass all possible error input patterns with a certain level of confidence.In previous research [166], a surface vessel was employed as a simulation platform to conduct sea trials using underwater TAN and SLAM techniques.FFI developed the NavLab integrated navigation simulation system and the TerrLab terrain-matching navigation simulation system [47,167].These systems can be used for the algorithm verification and data postprocessing of underwater TAN systems, making them the earliest underwater TAN simulation platforms.The University of Southampton, UK, conducted research on the feasibility of using TAN technology to support AUVs in performing ultralong-range missions using the 6000AUV as the carrier.They also conducted a simulation verification using underwater TAN technology to assist AUVs in crossing the Arctic ice cover [9,28,168].Currently, almost all the research institutions engaged in underwater TAN technology have developed their own simulation systems.However, most simulation studies have been limited to algorithm verification platforms based on data playback.In general, an underwater TAN simulation system should include a basic terrain measurement simulation, an AUV dynamics simulation, a TAN algorithm simulation, and other components (Figure 17).In particular, it should possess probabilistic completeness simulation capabilities for different terrain conditions and measurement systems, making it possible to simulate almost all possible operating conditions in a limited number of simulation experiments.These tasks cannot be accomplished through actual sea trials and are the core content of underwater TAN simulation systems.Some progress has been made in related research.For example, in terms of bathymetric system error analysis, researchers [169][170][171][172] have studied the error composition of multibeam bathymetric measurement systems with multiple sensors.They analyzed the distribution characteristics of depth measurement errors with respect to the wave speed and incident angle in a single measurement frame.Figure 18 shows the estimated measurement errors of multibeam bathymetric measurements at different beam angle positions as output by this simulation system.These results provide a theoretical basis for simulating bathymetric measurement errors using multibeam systems.
Currently, accurately simulating input errors for TAP systems is difficult.This is primarily because of the highly nonlinear characteristics of terrain measurement systems, where measurement errors are the result of coupling among various sensors in the measurement system.Figure 19 shows the residual sequence of the matching points in TAP.The signal in Figure 19b exhibits high autocorrelation and does not follow a Gaussian distribution, mainly because the prior terrain map was obtained through low-resolution map interpolation.Although it has been proved [5] that the residuals of terrain-matching positioning systems exhibit asymptotic Gaussian characteristics, the residual sequence of individual measurement beams still demonstrates a strong autocorrelation.Furthermore, matching residuals are closely related to terrain features, prior terrain maps, underwater acoustic environments, seabed geological conditions, and other factors.Therefore, the primary challenges are simulating non-Gaussian signals that match the measurement system, terrain features, acoustic environment characteristics, etc. and ensuring that the limited number of simulated residual signals satisfies the statistical completeness requirements of error signals.This is essential for conducting effective and reliable simulations of terrain-matching positioning and is a problem that must be addressed urgently.
multibeam bathymetric measurement systems with multiple sensors.They analyzed the distribution characteristics of depth measurement errors with respect to the wave speed and incident angle in a single measurement frame.Figure 18 shows the estimated measurement errors of multibeam bathymetric measurements at different beam angle positions as output by this simulation system.These results provide a theoretical basis for simulating bathymetric measurement errors using multibeam systems.Currently, accurately simulating input errors for TAP systems is difficult.This is primarily because of the highly nonlinear characteristics of terrain measurement systems, where measurement errors are the result of coupling among various sensors in the measurement system.Figure 19 shows the residual sequence of the matching points in TAP.The signal in Figure 19b exhibits high autocorrelation and does not follow a Gaussian distribution, mainly because the prior terrain map was obtained through low-resolution map interpolation.Although it has been proved [5] that the residuals of terrain-matching positioning systems exhibit asymptotic Gaussian characteristics, the residual sequence of individual measurement beams still demonstrates a strong autocorrelation.Furthermore, matching residuals are closely related to terrain features, prior terrain maps, underwater Depth error (m) mary challenges are simulating non-Gaussian signals that match the measurement system, terrain features, acoustic environment characteristics, etc. and ensuring that the limited number of simulated residual signals satisfies the statistical completeness requirements of error signals.This is essential for conducting effective and reliable simulations of terrain-matching positioning and is a problem that must be addressed urgently.

Conclusions and Outlook
Underwater TAN performs tracking and positioning navigation with non-Gaussian inputs, strong nonlinearity, and analytically intractable system state transition equations.Compared with other navigation methods, such as dead reckoning and acoustic positioning, TAN not only shares common issues related to position information fusion and estimation but also faces specific system characteristic problems owing to the uniqueness of the system input information (terrain information, terrain measurement errors, etc.).These issues result in non-Gaussian positioning errors, weak correlations between TAP points, initial instability during system operation, and the poor generalization performance of the algorithms.In addition, underwater TAN faces technical challenges, such as insufficient prior terrain maps and their accuracy, inadequate adaptability to terrain features, the accuracy and effectiveness evaluation of positioning, and high-precision system simulation.This article summarized the status of these issues, and the technical challenges that are currently faced were discussed.
Although underwater TAN technology has undergone more than 40 years of development, it has not yet been widely applied, and several important technical challenges must be overcome.Specifically, the following aspects should be addressed.
(1) The stability and reliability of filters remain significant challenges, particularly for state estimation in non-Gaussian and nonlinear dynamic systems.Therefore, there is an urgent need to develop robust and reliable filtering algorithms and filters.In addition, filters should possess self-awareness and self-correction capabilities to enhance the generalization performance of the algorithms.For example, the work [97,98] has introduced machine learning into image matching, the work [173] has incorporated intelligent algorithms into Particle Filters (PFs), and the work [75] has proposed data-driven methods that can learn approximate proposal distributions from previous data.These are all good attempts.(2) With the diversification and refinement of AUV sensing information, there are differences in the granularity, data structures, and physical characteristics of the matching information.This involves matching and assimilating different resolutions, granularities, and data structure information.(3) The integration of TAN technology with underwater robot planning and control is

Conclusions and Outlook
Underwater TAN performs tracking and positioning navigation with non-Gaussian inputs, strong nonlinearity, and analytically intractable system state transition equations.Compared with other navigation methods, such as dead reckoning and acoustic positioning, TAN not only shares common issues related to position information fusion and estimation but also faces specific system characteristic problems owing to the uniqueness of the system input information (terrain information, terrain measurement errors, etc.).These issues result in non-Gaussian positioning errors, weak correlations between TAP points, initial instability during system operation, and the poor generalization performance of the algorithms.In addition, underwater TAN faces technical challenges, such as insufficient prior terrain maps and their accuracy, inadequate adaptability to terrain features, the accuracy and effectiveness evaluation of positioning, and high-precision system simulation.This article summarized the status of these issues, and the technical challenges that are currently faced were discussed.
Although underwater TAN technology has undergone more than 40 years of development, it has not yet been widely applied, and several important technical challenges must be overcome.Specifically, the following aspects should be addressed.
(1) The stability and reliability of filters remain significant challenges, particularly for state estimation in non-Gaussian and nonlinear dynamic systems.Therefore, there is an urgent need to develop robust and reliable filtering algorithms and filters.In addition, filters should possess self-awareness and self-correction capabilities to enhance the generalization performance of the algorithms.For example, the work [97,98] has introduced machine learning into image matching, the work [173] has incorporated intelligent algorithms into Particle Filters (PFs), and the work [75] has proposed datadriven methods that can learn approximate proposal distributions from previous data.These are all good attempts.
(2) With the diversification and refinement of AUV sensing information, there are differences in the granularity, data structures, and physical characteristics of the matching information.This involves matching and assimilating different resolutions, granularities, and data structure information.(3) The integration of TAN technology with underwater robot planning and control is crucial.The incorporation of TAN information into intelligent decision-making and control systems for AUVs should be explored.(4) The availability of large-scale, high-precision prior terrain maps remains a major bottleneck in the development of underwater TAN technology.Promising breakthroughs can be achieved through multi-AUV cooperative underwater positioning and terrain mapping techniques, which have the potential to improve the measurement accuracy and efficiency for unknown seafloor terrain.For instance, the work [174] and the work [152] proposed two different strategies for collaborative SLAM techniques involving multiple Autonomous Underwater Vehicles (AUVs).

35 Figure 1 .
Figure 1.Underwater TAN theory and technical system.

Figure 1 .
Figure 1.Underwater TAN theory and technical system.

Figure 2 .
Figure 2. Main functional modules of TAN.

Figure 2 .
Figure 2. Main functional modules of TAN.
b y,k and ω b z,k denote the angular velocities of the AUV around the y-axis and z-axis, respectively, R qk−1 represents the rotation matrix formed by the AUV attitude matrix qk−1 at time k − 1, v k−1 represents the velocity matrix in the body coordinate system, ∆∅ k−1 and ∆θ k−1 represent the changes in the roll and pitch angles measured by the attitude sensors, respectively, ω y,k−1 and ω z,k−1 represent the angular velocities, a represents the unit vector in the sonar beam direction, r k represents the distance measured by the beam at time k, and r •,k−1 represents the measurement error of the sensor.

Figure 5 .
Figure 5. Operation stages of the TAN system.

Figure 5 .
Figure 5. Operation stages of the TAN system.

Figure 6 .
Figure 6.Search and matching process of TERCOM.

Figure 6 .
Figure 6.Search and matching process of TERCOM.
), where the point set  =  , … … ,  represents the reference point set, and  =  , … … ,  represents the point set to be aligned.Here, N represents the number of elements in the point set [86].Point set  can be aligned with point set  through rotation matrix  and translation matrix . =       = ‖ −   ‖

Figure 8 .
Figure 8. PF solution for nonanalytical and non-Gaussian posterior Bayesian estimation.

Figure 17 .
Figure 17.Diagram of a TAP simulation system.Figure 17.Diagram of a TAP simulation system.

Figure 17 .
Figure 17.Diagram of a TAP simulation system.Figure 17.Diagram of a TAP simulation system.

Figure 19 .
Figure 19.TRP residual signals of a multibeam measurement system.(a) Overlapping of multiple terrain-matching residual signals.(b) Individual terrain-matching residual signal.

Figure 19 .
Figure 19.TRP residual signals of a multibeam measurement system.(a) Overlapping of multiple terrain-matching residual signals.(b) Individual terrain-matching residual signal.

Table 2 .
Primary algorithms utilized in the three operational phases of a TAN.

Table 2 .
Primary algorithms utilized in the three operational phases of a TAN. ) ∂x 1 and ∂h i ∂x 2 represent the partial derivatives of the terrain surface at the i-th footpoint in the x 1 and x 2 directions, 1−α represents the sum of squared residuals of the matching residuals at a confidence level of 1 − α.