Survey of Terrain-Aided Navigation Methods for Underwater Vehicles

The prerequisite for underwater vehicles to accomplish various underwater tasks is to obtain real-time position information about themselves in the underwater environment using internal or external sensors, which is known as underwater positioning and navigation technology. This paper presents a general review of underwater vehicle positioning and navigation technology and illustrates the importance of underwater terrain-aided navigation technology. It also summarizes its general structure and describes the principles and functions of each component. The paper mainly analyzes various mainstream methods applied to underwater terrain-aided navigation, focuses on terrain suitability analysis and terrain matching algorithms, which play a central role in the analysis. Additionally, the research progress of terrain-aided navigation technology is described, and the research methods of terrain suitability analysis, terrain matching navigation, and the latest technology applications are compared. Some shortcomings of the current underwater terrain-aided navigation technology are presented, and future research directions are foreseen to provide a basis for subsequent research.


I. INTRODUCTION
From a flat boat to a giant ship, and from the development of water routes to underwater diving, human beings have been exploring the ocean. With the continuous development and breakthrough of sensor technology, control science, and artificial intelligence, the Autonomous Underwater Vehicle (AUV) has been used to explore the ocean. In recent years, countries worldwide have increased their investment in research. The navigation and positioning system is a crucial part of the AUV, as it confirms the vehicle's position and reaches the target point. Unlike on land and in the air, electromagnetic waves have serious attenuation problems in water. Therefore, the Global Navigation Satellite System (GNSS), which has high positioning accuracy and is frequently used on land and air, cannot be extended for use underwater [1].
The associate editor coordinating the review of this manuscript and approving it for publication was Rosalia Maglietta .
Finding a navigation and positioning method that is passive and autonomous, with high positioning accuracy and good concealment, is the focus of AUV research and a significant challenge. The terrain-aided navigation system has gradually become the main research direction for AUV positioning and navigation due to its advantages of passive autonomy, wide navigation range, high positioning accuracy, and strong antiinterference. The main contributions of this paper are the following three points: 1) Previous review-type papers on terrain-aided navigation have either been less comprehensive or simpler [2], [3] for technical application methods, and rarely provide comprehensive inductive discussions. This paper provides a general overview of terrain-aided navigation systems, summarizes the methods for each component, and provides a detailed discussion that provides the basis for subsequent research.
2) This paper provides a comprehensive summary and generalization of the principles and methods involved in terrain suitability analysis. It divides the methods for characterizing the inherent feature information of regional physical fields in a dimensionality reduction manner into feature extraction methods for invisible characterization and feature selection methods for explicit characterization. The suitability region is classified into direct classification methods and classification methods for pattern recognition fields, highlighting the latest research directions represented by data-driven approaches. 3) This paper defines terrain matching algorithm as an optimization problem to solve the similarity between matching data, comprehensively summarizes and compares various terrain matching algorithms, and provides an in-depth analysis of terrain matching algorithms based on statistical theory from the direction of algorithm optimization, which provides a reliable reference for future research. The general framework of the paper is described below. Section II describes several major navigation and positioning methods, analyzes and compares the characteristics of each navigation and positioning method, and points out the importance of terrain-aided navigation technology. Section III describes the components of the terrain-aided navigation system, presenting the functions of each unit individually and analyzing the main methods and their characteristics. Considering that the suitability analysis of the terrain matching area and the terrain matching algorithm are the core of the terrain-aided navigation system, these two parts are separately described in Sections IV and V, respectively. Finally, Section VI summarizes the paper.

II. UNDERWATER NAVIGATION METHODS
Underwater navigation and positioning mainly involve several navigation methods such as Inertial Navigation System (INS), Acoustic Navigation System (ANS), Geophysical Aided Navigation System (GANS), Simultaneous Localization and Mapping (SLAM) and Multi-AUVs Collaborative Navigation (MACN). These methods are briefly described below.

A. INERTIAL NAVIGATION SYSTEM
The basic principle of INS is to obtain the position direction, rotational angular velocity, and acceleration of the underwater vehicle using the 3-axis gyroscope, 3-axis accelerometer, and Doppler Velocity Log (DVL). It then integrates time to deduce the current position [4]. For short-range AUVs, INS can achieve navigation and positioning more accurately by itself. However, over time, the accumulated error of INS will gradually increase, leading to deviation from the target position. Therefore, when using INS as the primary navigation unit, other navigation and positioning systems must be used as auxiliary units to correct INS's positioning error in time.

B. ACOUSTIC NAVIGATION SYSTEM
The principle of ANS is to use the transmission time of acoustic signals and their phase difference to determine the position of the underwater submersible, which is more accurate than the inertial navigation system. Based on their different baseline lengths, ANS can be divided into three types: long baseline (LBL), short baseline (SBL), and ultra-short baseline (USBL) as shown in Figure 1. Sture et al. [5] used the ultra-short baseline acoustic positioning system to help one or more underwater vehicles achieve more accurate navigation.
All three navigation types require the deployment and calibration of the base array, which is time-consuming and laborintensive. Furthermore, their navigation ranges are influenced by the acoustic base array, which is not conducive to navigation across the entire sea.

C. GEOPHYSICAL AIDED NAVIGATION SYSTEM
GANS uses the relative positions of observable geophysical features to obtain the AUV's own positioning [6]. The positioning method is defined by the inertial system, and then the data matching algorithm is used to match the geophysical measurements acquired by the sensors with the priori geophysical database to obtain the real-time position of the underwater submersible, update the positioning error of the INS, and thus estimate the position of the AUV itself. Zhang et al. [7] proposed a navigation path that relies on the surrounding geomagnetic environment, uses the heading angle predicted by the geomagnetic gradient to constrain the navigation path, and improves the evaluation function based on the principle of simultaneous convergence of multiple parameters, thus improving the reliability and accuracy of navigation. VOLUME 11, 2023 47511 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.

D. SIMULTANEOUS LOCALIZATION AND MAPPING
SLAM is also a hot topic of current research. In an unknown environment without the assistance of an acoustic marker base array and without a priori geophysical database, researchers have attempted to allow vehicles to autonomously create maps using sensors such as sonar and cameras to collect information about the surrounding environment and to implement localization and navigation based on this [8], [9]. This approach is still essentially a physical environmentbased navigation method [10], a method that uses geophysical information from the surrounding environment to achieve autonomous navigation and localization. Li et al. [11] successfully implemented real-time position correction for SLAM systems using forward-looking sonar (FLS) images by solving the problems of data association, key frame selection, and outlier rejection.

E. MULTI-AUVS COLLABORATIVE NAVIGATION
The MACN approach solves this problem better in the case where the underwater vehicles are loaded with low sensor accuracy and a single robot cannot accomplish the navigation task. MACN can effectively reduce the accuracy error of the inertial guidance system and improve the intelligence and efficiency of multiple AUVs to accomplish tasks that are difficult to be accomplished by a single submersible [12]. There are 2 types of MACN, namely parallel and primary-secondary, depending on their structure. The primary-secondary approach is cheaper and more accurate than the parallel approach, which makes it the mainstream approach for MACN. Allotta et al. [13] designed a primary-secondary MACN system in which multiple AUVs are equipped with low-cost inertial guidance units and depth sensors, only one of which has a high-precision DVL, and the DVL-equipped primary AUV is used to determine the positions of other secondary AUVs. position and communicate with the mother ship through GPS to obtain the specific position of each AUV. Chen et al. [14] proposed two new fault-tolerant control techniques for leaderless multi-AUV systems and lead-following multi-AUV systems, which make each AUV exchange information with its neighboring AUVs, i.e., a distributed consensus control strategy to tolerate information transmission failures of multiple AUV systems.
Several underwater positioning and navigation methods are analyzed and compared in Table 1.
The INS cannot meet the long-endurance and highprecision navigation and positioning requirements due to the inherent defects of the system. ANS has higher positioning accuracy, but it cannot be applied to a wide range of navigation and positioning tasks due to the high complexity of base array deployment and small deployment range. GANS has high positioning accuracy and wide navigation range, but it needs INS guidance to matchable area before it can be used. Therefore, the future navigation system for underwater intelligent robots incorporates a combination of technologies for navigation [15]. INS has the characteristics of passive, fully autonomous, stealthy, and all-sea navigation. And GANS relies on external sensors that do not cause accumulated errors over time, which is the best choice to avoid INS drift problems and can correct INS errors for a long time [16]. A geophysically assisted navigation system that combines the features of INS passive, fully autonomous, good stealth, and all-sea navigation with the advantages that the external sensors of GANS do not lead to accumulated errors over time is the future trend of underwater positioning and navigation [17].

III. TERRAIN AIDED NAVIGATION SYSTEM
Geophysical Aided Navigation is a data-driven aided navigation system that uses information from geophysical data of the surrounding environment. Geophysical data is mainly obtained by measuring the physical environment surrounding an underwater vehicle using specific sensors. The types of geophysical data include visual imaging data, sonar imaging data, geomagnetic data, gravity data, terrain elevation data, and others. Each of these data types has its own characteristics and is suitable for positioning and navigation in different underwater scenarios, Table 2 provides an analysis and comparison of these data types. Visual imaging data are more accurate [18], [19] and are often used for end-target localization at close range, with a smaller range of application and are only suitable for use in clear water environments, and their effectiveness is poor in turbid waters. Sonar imaging data is an alternative to visual imaging data [20], [21], but it often produces blurred images and is not suitable for accurate matching localization because most of the seabed images are not distinctly characterized. The geomagnetic data have a relatively large range of matching areas and low accuracy. The measurement instruments for gravity data are relatively expensive and large, making them unsuitable for widespread installation on small underwater vehicles. With the improvement and development of sonar rangefinder over the years, the accuracy and range of terrain elevation data have been greatly improved, and its small size and moderate cost make 47512 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.  it suitable for ''high-precision, long-range'' underwater positioning and navigation. This paper will mainly focus on Terrain Aided Navigation System (TANS). TANS was first applied as a navigation method for land vehicles and air vehicles. This method can provide navigation and positioning functions very well in the absence of GNSS. Hollowell [22] provided a TAN algorithm for helicopter navigation systems using radar altimeter ground clearance measurements combined with conventional navigation systems and stored digital terrain elevation maps. A set of single-state Kalman filters employing a multi-model adaptive estimation (MMAE) technique ensures reliable position estimates even in the face of large initial position errors. The algorithm reduces the corrected center radial error to less than 50m. Xian et al. [23] analyzed the isometric flight scheme and the inertial system-based absolute altitude solution when the barometric altimeter does not work properly. This allows for the realization of the absolute altitude solution based on strapdown inertial navigation system (SINS), and the simulation results show that this solution basically meets the accuracy requirements of the terrain matching-aided navigation system. This method has a height error of no more than 10.52m per unit time with gyroscopic error and accelerometer error.
Bergem [24] was one of the most early proponents, who used the contours from multibeam sonar to match with a pre-stored reference map to obtain an estimate of absolute position. Massa et al. [25] proposed a Kalman filterbased terrain-aided navigation method suitable for real-time position estimation, but could not resolve the uniqueness of the matching results. Strauss et al. [26] proposed a terrain-referenced localization method based on matching rough sonar data with digital elevation maps. They compared two matching algorithms based on statistical and fuzzy methods, respectively, and used the Kalman filter to predict the position of underwater divers more accurately in simulations. In the simulation experiment, this method reduces the average error in both directions in the plane right-angle coordinate system from 80m to 65m.Ingemar Nygren [27] analyzed correlation-based terrain navigation using 3D sonar. One of the most important results is the theoretical demonstration that the maximum likelihood estimator is optimal under certain conditions, and that the likelihood measurement function converges asymptotically to a Gaussian probability distribution as the number of measurement beams increases.
Terrain-assisted navigation systems have been developed and studied by researchers in many countries as an applied VOLUME 11, 2023 technology in engineering. In 2016, the Institute of Industrial Science of the University of Tokyo built an AUV called the Highly Agile Terrain Tracker for Ocean Research and Investigation (HATTORI) [28], equipped with a scanning sonar for seafloor exploration with a beam frequency of 675 kHz, a maximum range of 100 m, and a beam width of 3 × 30 • . Maki et al. [28] successfully tracked 90m of rocky terrain in 130 seconds during the 2016 experiment and 40m of coral reef terrain in 85 seconds during the 2017 experiment in two at-sea trials using this AUV. FFI developed a HUGIN 1000 HUS AUV [29], and conducted two sea trials in Oslofjord, where the error of terrain-matching positioning with hydroacoustic positioning was 4m and 5m with GPS in a 50km trip, showing the high positioning accuracy of terrainmatching positioning. Rock et al. [30] improved MBAIR and conducted a sea trial in 2014 at Portuguese Ledge in Monterey Bay to verify the possibility of successful AUV recovery using terrain-assisted navigation. The results showed that the AUV successfully returned to the recovery point after a four-stage test path. Zhang et al. [31] conducted an experimental validation of the proposed improved interated closest contour point (ICCP) algorithm at the Xin'an River Reservoir in Hangzhou, China, equipped with an MBES iBeam 8140 sonar sensor and a fiber optic gyroscope for INS heading projection, using GPS for creating topographic maps and reference paths. The experimental study shows that the improved ICCP algorithm achieves higher positioning accuracy than the conventional ICCP in the case of small initial errors, and relatively better positioning accuracy in the case of larger initial errors where the conventional ICCP cannot match the positioning. Ding et al. [32] experimentally validated the proposed terrain-assisted navigation algorithm using navigation-grade SINS, bidirectional antenna real-time kinematic global navigation satellite system (RTK-GNSS), and multibeam bathymetry, and the experimental results showed that the matching accuracy of the algorithm improved by 23.97% over the ICCP algorithm.
TANS consists of three components, as shown in Figure 2, which are the inertial navigation unit, the terrain matching unit, and the path planning unit.

A. INERTIAL NAVIGATION UNIT
The inertial navigation unit is the basic navigation unit of the system and provides the determination of the initial local matching region. The range of the matching region can be determined in two methods: the matching region determined by the error ellipse and the matching region determined by the INS drift.

1) THE MATCHING REGION DETERMINED BY THE ERROR ELLIPSE
The first method requires calculating the maximum error of the INS and defining an ellipse with the position provided by the INS as the center of the ellipse. The search region is then defined as a rectangular region, which is the minimum enclosing rectangle of the ellipse region, as shown in Figure 3. If σ x , σ y are the standard deviation of the INS positioning system in the eastward and northward directions; σ 2 x , σ 2 y , σ xy , σ yx are the variance and covariance of the INS positioning system, respectively, the error ellipse can be determined by the following equation: where a and b are the long and short semi-axes of the ellipse, respectively; θ is the angle between the long semi-axis of the ellipse and the positive north direction;σ 0 is the expansion factor.
Since the position error of the AUV is almost proportional to the circular error probability (CEP) at that position [33],the distribution of the positioning error of the INS is nearly a standard normal distribution, according to the ''3σ '' principle, ifσ 0 = 3.03, then the probability that the measured point lies within the positioning error ellipse is 95%.
The error elliptic exterior tangent rectangle is calculated as: where x m and y m are the lengths of the long and wide sides of the error elliptic exterior tangent rectangle.

2) THE MATCHING REGION DETERMINED BY THE INS DRIFT
Another method [34], [35] assumes that the drift produced by INS since the last update is σ . This method centers on the inertial guidance indication position with the ±3σ error of the system as the magnitude, extract the terrain elevation parallel to the navigation indication path in a certain area around the navigation position from the terrain reference map area, i.e., the matching search area, as shown in Figure 4, and traverse the search on the base map. Based on this, Liu et al. [36] proposed a method to construct an adaptive grid model to improve the average positioning accuracy to 23.6m for the situation where the regular grid model uses equidistant dissection of the terrain. This method results in a less detailed description of complex terrain areas and an overly cumbersome description of simple terrain areas.

B. TERRAIN MATCHING UNIT
The terrain matching unit first performs a suitability analysis of the matching area based on the prior terrain elevation database. Then, it uses the terrain data acquired by the sonar sensors and the prior terrain database as inputs for the terrain matching algorithm to estimate the true position of the vehicle.

1) DIGITAL ELEVATION MODEL
Due to the unique characteristics of the underwater environment, obtaining high-resolution underwater topographic images is very challenging. Therefore, sonar rangefinders and Digital Elevation Model (DEM) are commonly used to acquire topographic elevation data for topographic matching, and the position of the AUV is estimated based on the matching results. The DEM is a special grid data model that can represent continuous spatial undulation changes. The position of each grid is represented by a 2D vector expressed by latitude and longitude coordinates, with the elevation value being the attribute value of each grid. The grid size is called spatial resolution [37]. The DEM is particularly suitable for constructing underwater terrain elevation models and for terrain matching navigation because of its ease of storage, updating, propagation, and automatic computer processing, as well as its multi-scale characteristics, which make it suitable for various quantitative analyses and 3D modeling [38].

2) TERRAIN MATCHING ALGORITHM
Terrain matching algorithms can be classified into direct image matching-based approach, image feature point-based approach, and sensor-measured terrain elevation data-based approach, depending on the data used. Currently, most terrain matching algorithms are implemented based on DEM. According to different matching methods, they can be further classified into correlation-based approach, approach based on extended Kalman filtering, approach based on direct probability criterion, and approach based on image matching [3]. Terrain suitability analysis and terrain matching algorithms are the core components of the terrain matching unit, and these aspects will be discussed separately in Sections IV and V of this paper.

3) TERRAIN DEPTH SOUNDER
The seafloor topography measurement mainly consists of two parts, as in Figure 5, if h 1 denotes the depth from sea level measured by pressure transducer, h 2 denotes the measured depth from sonar transducer, then the absolute water depth value H can be expressed as: Sonar sensors are mainly used to derive seafloor elevation values by transmitting and receiving sound waves returned from the seafloor and calculating the phase angle and delay time between them [39]. Additionally, the underwater sound velocity is related to various factors such as temperature and salinity of seawater. Therefore, only considering these factors together can result in more accurate measurement of seafloor topographic elevation values. Kamolov et al. [40] investigated the application of the fuzzy c-means (FCM) clustering algorithm on crowdsourced bathymetry and obtained more accurate results, with an average absolute error of 1.67 m and an average value of 2.09 m between its predicted depth and the real data. Sonar sensors can be classified as single-beam sonar sensors, DVL sonar sensors, side-scan sonar sensors, and multi-beam sonar sensors. DVL sonar sensors are sensors that use the Doppler effect to measure Vehicle's velocity VOLUME 11, 2023 and are generally not used directly for topographic elevation measurements. Side-scan sonar sensors do not provide bathymetric data, but rather provide a relatively high-resolution acoustic image of the seafloor, which can be interpreted as a monochromatic representation of the physical properties of the seafloor [41].
Both single-beam and multibeam sonar sensors are used to measure topographic elevation values directly, and single-beam sonar achieves a ''point-line'' measurement range in the case of continuous measurements. Compared with single-beam sonar, multibeam sonar is a qualitative leap, multibeam sonar sensors form hundreds of irradiated footprints in the strip area perpendicular to the ship's direction of travel [42], thus further forming a ''line-surface'' measurement range, as shown in Figure 6. Multi-beam sonar has a high resolution compared to single-beam [43], and therefore, multi-beam sonar sensors have been routinely used to measure large areas of the seafloor [44]. Multi-beam bathymetric systems are integrated systems composed of several subsystems. They generally consist of four parts: the acoustic subsystem, data acquisition subsystem, data processing subsystem, and peripheral auxiliary equipment. These are responsible for tasks such as acoustic wave transmission and reception, data correction, data conversion, and comprehensive data processing, respectively [45]. Coherent research has mechanically solved the problem of limited resolution in multibeam bathymetry, caused by the number of beams [46], greatly improving the resolution of multibeam bathymetry data.

C. PATH PLANNING UNIT
AUV path planning can be implemented based on terrain suitability analysis and sonar detection sensors in five steps: environmental information acquisition, underwater environment modeling, global path planning, underwater situational awareness, and local path planning [47]. The first three steps rely on the suitability analysis of a priori terrain maps, while the last two steps depend on the detection of sonar sensors to achieve path planning with accurate localization, shortest distance, and real-time obstacle avoidance.
The current path planning algorithms mainly include heuristic optimization algorithms and machine learning algorithms. The ideas behind heuristic optimization algorithms are derived from long-term observation and practice of physical and biological phenomena, as well as the summary and generalization of natural phenomenon [48], including Genetic Algorithm (GA), Ant Colony Optimization (ACO), Particle Swarm Optimization (PSO), Firefly Algorithm (FA), A * algorithm, RRT * algorithm, etc. Due to the limited space, the principles of some of these algorithms are given as examples. The particle swarm algorithm originated from the study of bird feeding behavior, in which birds search for food randomly, not knowing where their current position is from the food, the simplest and most effective search strategy is to search the area around the bird that is currently closest to the food. The simplest and most effective search strategy is to search the area around the bird that is currently closest to the food and constantly adjust the state of the particles to find the food and reach the global optimum [49]. The basic principle of the ant colony algorithm is that ants release pheromone in the process of searching for food. Combined with the property of its volatilization, the pheromone concentration of the shortest path is higher than other paths as time increases, Yan et al. [50] proposed a good information matrix that differentially distributes the initial pheromone concentration, avoiding blind search, reducing the search range, and shortening the search time. Compared with the traditional ant colony algorithm, this algorithm reduces the planned path from 40.28km to 38.63km, and the number of iterations to convergence is reduced from 155 to 68. Cong et al. [51] compared four different path planning methods: genetic algorithm, A * algorithm, RRT * algorithm and ant colony algorithm, taking into account the terrain complexity and the motion characteristics of the autonomous underwater vehicle, they verified that all four algorithms could provide a reliable path for the TANS.
Machine learning algorithms view the path planning problem as a Markovian decision process with penalties or rewards in constant trial to optimally update the path. Sun et al. [52] proposed a 3D path planning for AUVs using hierarchical deep Q networks (HDQN) combined with prioritized experience replay, which divides the path planning task into three layers to achieve dimensionality reduction of the state space and solve the dimensional disaster problem. Patle et al. [53] applied the fuzzy algorithm as a fuzzy function to a mobile robot navigation system. They set up two parallel rules, the first to classify the metric and the second to set up an optimal decision. The fuzzy function provides the robot the ability to decide the choice of paths and rank them. The paths consist of roots, and a new path can be set up for the target using the roots and their affiliation functions. The algorithm controls the error rate to within 5%. Kahn et al. [54] proposed a self-supervised learning-based mobile robot navigation system that can be trained with automatically labeled non-strategic data collected in a real environment without any simulation or human supervision. This system can be improved as more data is collected to continuously optimize the path. The algorithm achieved a successful arrival rate in experiments that is nearly 40% higher than conventional algorithms. Li et al. [55] constructed an Occupancy Grid Map (OGM) for path planning. It is a grid mapping in which each grid is assigned a unique value using a specific potential function. In each step, the unsearched neighboring grid with the smallest value is selected. When the vehicle detects an obstacle, it resets all grid values that are not 1 (the grid occupied by the obstacle will be assigned a value of 1). Using this resettable OGM, the vehicle will search for an alternative path to approach the target point while avoiding collisions with the obstacle.
There are many kinds of path planning methods, and they are widely used, mature, and selectable. For a terrain-aided navigation system, it is enough to choose the appropriate method according to the data and algorithm. Therefore, this paper will not discuss this in-depth.

IV. SUITABILITY ANALYSIS
Suitability, also known as navigability, is a physical concept used to reflect the richness of a region's inherent geophysical fields (e.g., topography, geomagnetism, gravity) containing navigational information features. It is a property inherent to any finite region that reflects the ability of a region to provide planar position information along a particular directional profile [56]. The terrain suitability analysis is based on the terrain feature parameters derived from the terrain elevation values, and the terrain suitability analysis is a prerequisite for TANS, which is based on the pre-acquired DEM. The purpose of terrain suitability analysis is to find the suitable adaptation area for terrain matching algorithm, which is roughly divided into two parts: the extraction and selection of features and the selection of the adaptation area using the processed terrain features.

A. EXTRACTION AND SELECTION OF TERRAIN FEATURES
Terrain features are metrics used to characterize the geometric relationship between terrain elevation data and its similarity. The more obvious terrain features that can be found in a region, the better the adaptation performance and the higher the accuracy of terrain matching. The basic terrain adaptation features include elevation mean, elevation standard deviation, slope standard deviation, correlation length, correlation coefficient, kurtosis coefficient, skewness coefficient, cumulative gradient mean, Feche information, terrain roughness, abundance coefficient, terrain entropy, etc. [57]. Because all these parameters characterize the terrain features of the terrain elevation map from different perspectives, it is difficult to describe the adaptation performance of this region in a more comprehensive way using a single feature parameter for terrain feature characterization. If all these parameters are used to characterize the terrain features, any point of the feature parameters will be a very high-dimensional feature vector, which is not conducive to the subsequent calculation, and there are redundant relationships among these feature parameters, which will also cause a waste of computational resources. Therefore, it is necessary to reduce the dimensionality of these feature parameters.
Dimensionality reduction is a method to solve dimensional catastrophes and other irrelevant properties in high-dimensional spaces [58], which is divided into linear and nonlinear dimensionality reduction. Principal Component Analysis (PCA) is a linear dimensionality reduction method, which reduces dimensionality by retaining low-order principal components and ignoring high-order principal components through linear transformation. In contrast, Roweis et al. [59] proposed a method to recover the global nonlinear structure from a local linear fit, achieving a better nonlinear dimensionality reduction. Feature dimension reduction is a process of selecting a low-dimensional feature set from an initial high-dimensional feature set to optimally reduce the feature space according to certain evaluation criteria [60], which is divided into two ways: feature extraction and feature selection.

1) FEATURE SELECTION
Feature selection (FS) is the process of selecting some features from the original data feature set based on the relevance and redundancy of the features so that they satisfy specific metrics, and is an explicit expression of the integrated features. The selected features are a subset of the original features, hence it is also called feature subset selection [61]. This process of selecting can also be seen as a gradual merit-seeking process, which contains elements such as merit-seeking strategies, evaluation metrics, and stopping conditions. The optimization strategy is mainly based on heuristic algorithms, such as [62], which uses genetic algorithms to select high-performance synthetic features from the synthetic feature space by performing selecting, crossover and variation operations on multiple randomly generated synthetic feature trees individuals using the effectiveness evaluation value as a metric. The evaluation metrics are the basic indicators for judging the performance of regional fitness. The matching probability is the most commonly used metric [63], and its formula is defined as follows: where CMA is the candidate matching area;CMM is the matching degree of a point p in the CMA, and CMM(p) = 1 if the match is successful, and CMM(p) = 0 if the match fails; PCM is the matching probability; d denotes the square matching area of a given size in the neighborhood of point p; NCMA is the total number of matches in the CMA area. Numerous studies have been conducted by scholars to investigate the matching performance of terrain features in regions to be matched. Wang et al. [64] fused a hierarchical statistical significance detection method and a frequency domain decorrelation significance detection method to create a fused terrain feature map, and then matched the map using normalized cross correlation (NCC) to quickly select landmark points. This method achieved a correct matching rate of 73.9%-88.3% in actual terrain matching experiments. Zhang et al. [65] constructed a BP neural network to establish the mapping relationship between input terrain feature covariates and matching performance, thereby achieving automatic classification of the matching area with a classification rate of around 90%, more than 10% higher than traditional classification methods. Teng et al. [66] used terrain entropy and terrain variance entropy as terrain features to analyze the matching performance of a region, this method uses the joint criterion and fuzzy criterion to calculate the weight of entropy (criterion) in the cost function, and the path planning time in the experiment was 16.606 seconds, with a matching error in each matching region of less than 10 meters. Li et al. [67] combined the shortest arc law in spherical geometry with attitude control theory in space and marine environments to propose a new geodesic-based measurement method that reduces the radius of the search matching region. The search matching time was reduced from 9.84 seconds to 1.29 seconds (about 7.6 times) with high matching efficiency.

2) FEATURE EXTRACTION
Feature extraction (FE) is a process of reducing the high dimensionality of the feature vector by creating new features that depend on the original input feature set. This is done to provide an implicit representation of the integrated features. The new features created maintain the original relative distance between the features without losing a large amount of information during the feature transformation [68]. Among the latest research directions on FE, data-driven terrain navigation is increasingly becoming a hot topic of interest. This approach takes data as the input of the model, learns from the data, and optimizes the model in continuous iterations, which makes the FE capability of the model increase and gradually converge the model to the realistic state [69] . The process schematic is shown in Figure 7. Most of the data-driven FE methods are based on deep learning algorithms that are currently very widely used in the field of artificial intelligence. This method avoids the limitations of manually designed features, and its FE performance is closely related to the training data set and model structure.
In recent studies, more and more researchers are focusing on this approach. References [70], [71] used deep learning networks to extract sonar image features, trained using images generated by FLS simulators, and tested the model using images from real detections to verify the effectiveness of the method to achieve underwater navigation. Zhang et al. [72] proposed a new data-driven framework for template matching of underwater topographic images, which introduces contrast learning algorithms to terrain matching to achieve data enhancement by simulating grayscale and texture differences between images, enabling selfsupervised end-to-end learning without additional data annotation. Although the matching time of the model is slightly longer compared to the traditional SSD and NCC, the matching accuracy reaches up to 98%, which cannot be achieved by other traditional methods. Contante et al. [73] proposed a new deep neural network (DNN) architecture that predicts robot position and self-motion by processing image sequences to compute relative pose predictions while providing a measure of uncertainty about these estimates. Cohen et al. [74] proposed an end-to-end deep learning approach to estimate the velocity of an AUV instead of a traditional filter, using only INS and DVL, which improved the effectiveness by 66-87% over the traditional model-based approach. Melo et al. [75] proposed a data-driven particle filter, where sonar sensors are limited and cannot be modeled. The method learns from previous data to obtain an approximate estimated distribution, and it has relatively high accuracy and efficiency with a 40% lower data resampling frequency and 17% less running time than the traditional method.
Both FE and FS are methods for processing terrain data, and Table 3 summarizes and compares these two methods.
With the significant improvement in computer computing power, the rapid development of artificial intelligence technology, and the continuous enhancement of sensor measurement accuracy, there is an increasing trend towards using a data-driven approach in FE methods to handle terrain data and its matching algorithm. As a result, this method has gradually become a research hotspot.

B. TERRAIN REGION CLASSIFICATION
The purpose of terrain area classification algorithm is to classify the area to be matched in the DEM into suitable and non-suitable areas, which is a prerequisite for path planning and helps to improve the accuracy of terrain matching. It can be roughly divided into two methods. The first is the direct classification method (DCM), which involves classifying using a salient feature or combination of features, setting a suitable threshold value, and dividing the fitness region based on it. Wang et al. [76] achieved the division of terrain suitability areas and non-suitability areas by quantifying the terrain suitability features. Xu et al. [77] proposed navigation coefficients based on gray fuzzy theory for multi-feature parameter fusion and used it as the basis for terrain adaptation domain selection.
Another approach is implemented based on classification methods in the field of pattern recognition (PRCM), which can be divided into supervised and unsupervised classification according to the type of data and into linear and nonlinear classification models based on the type of model. Since terrain features are nonlinear and manual labeling is too costly, the current mainstream approach is unsupervised classification using nonlinear classification models. Cao et al. [78] used Gabor wavelet features of images and fused the image fitness parameters, then combined with the support vector machine (SVM) classification method to transform the matching probability estimation problem of any region within the image into the classification problem of pixels in the region to be estimated, and measured the matching performance of each region in the navigation baseline map from a quantitative perspective. The prediction probability of this method is within 5% of the actual probability error in more than 80% of the total, and the prediction estimation is good. Liu et al. [79] transferred the geomagnetic matching problem to statistical pattern recognition and analyzed it from the perspective of pattern recognition. They proposed an SVM-based matching algorithm, which greatly improved the matching accuracy. This method improved the matching rate by 74.43% and 52.86% in comparison experiments with Mean Square Deviation (MSD) method in two regions, respectively.
Since terrain features are nonlinear and manual labeling is too costly, the current mainstream approach is unsupervised classification using nonlinear classification models. It will be very interesting in future research if combined with data-driven feature extraction-based methods for terrain area classification.

V. TERRAIN MATCHING ALGORITHM
The terrain matching algorithm is the core of TANS, and its performance directly determines the accuracy of TANS, which is key to successful navigation. The terrain matching algorithm is based on the terrain elevation value or the terrain features derived from it, upon which the real-time terrain elevation data measured by the sonar sensor is matched with the underwater digital terrain database to estimate the current position of the vehicle. The workflow of the algorithm is shown in Figure 8.
The terrain matching algorithm is essentially an optimization problem that aims to solve the similarity between matching data, with the ultimate goal of finding the matching data with optimal similarity. Therefore, terrain matching algorithms are mainly divided into two categories: methods based on sequential data correlation and methods based on statistical theory, according to the different research methods.

A. METHOD BASED ON SEQUENTIAL DATA CORRELATION
The method based on coherent data correlation (CDCM) is used to find the matching point that has the highest correlation with the real-time elevation data. This is achieved by calculating the similarity between the real-time terrain elevation data and the area to be matched through a continuous search or iterative process to realize the positioning correction of INS. The methods mainly include Terrain Contour Matching (TERCOM) and ICCP algorithms [80].
TERCOM involves searching and matching real-time terrain elevation data with a terrain elevation database and calculating the correlation by a certain similarity metric. The resulting position of the correlated extreme point is the current position of the vehicle [81]. There are many ways to measure similarity, including Euclidean Distance, Manhattan Distance, Cross Correlation (COR), Mean Absolute Difference (MAD), Mean Square Error (MSE), Cosine, Terrain Entropy, Hausdorff Distance, and so on. Xu et al. [82] proposed the Hausdorff distance between the real-time bathymetry sequence and the search sequence on the topographic map of the seafloor as the metric of spatial curve matching, and gave the approximate real position of the vehicle through simulation. The algorithm rapidly reduced the 50% Circular Error Probability (CEP) from 2500m to 250m and kept the navigation stable in the simulation experiment. According to the basic principle of TERCOM, Rao et al. [83] proposed the terrain bathymetric slope sequence as a terrain feature parameter, and theoretically proved that the larger its mean value is, the smaller the matching error is and the better the matching navigation performance is. References [84], [85] introduced the particle swarm optimization algorithm into the matching area search, used the function that is related to normalized product as the particle fitness function, and compared the similarity between the baseline submap and the real-time map profile by the maximum metric value of the fitness. The simulation results showed that the algorithm was better than the traditional TERCOM algorithm. In [86], a moving least squares (MLS) method is used to reflect the curvature of the MLS surface at each water depth point, and a local curvature binary identification (LCBD) method descriptor is proposed to describe the shape of the local MLS surface. The points satisfying the LCBD case are divided into the set of sampling points, after which a GA is introduced to replace the traditional traversal method (TM) from the underwater a priori graph to search for the optimal estimated position of the vehicle. This method improves the matching localization efficiency over the conventional method by more than 80% at 7 localization points.
The basic idea of ICCP is to measure the topographic elevation data Hi of the point in real time based on a certain position i of the indicated place of the INS, find the nearest point on the isobath with depth Hi near it for each point, rotate and translate the measurement path to minimize the sum of squares of the distances between the set of points on the path after translation and the set of nearest points found, and repeat this process iteratively until the distance reaches a set threshold or the maximum number of iterations, then the calculation is finished [87]. Chen et al. [88] used Hausdorff distance to calculate the position error and particle swarm optimization algorithm to optimize the reference route. They proposed an improved ICCP algorithm to enhance the positioning accuracy and robustness of the algorithm, which reduced the mean value of the error from 83.2m to 14.4m for the trajectory position indicated by INS.
Wang et al. [89] proposed an improved multi-path parallel ICCP algorithm based on multibeam bathymetry data collected by multibeam echosounder. The algorithm selects data points located at the center and both sides of the edge of the strip sounding data, a total of three paths, to effectively solve the matching dispersion problem caused by the large initial error. Simulation results verify the effectiveness of the multi-path parallel ICCP algorithm, where the matching error between the matching trajectory and the real path of the algorithm is less than 100 m, even with an initial error of INS up to 1200 m.
Compared with TERCOM, ICCP has the process of rotation and translation, resulting in higher positioning accuracy. However, it is extremely sensitive to the error of the initial position, while TERCOM is not. Therefore, [35] divided the matching process into two stages, first using TERCOM for coarse matching and later using ICCP for fine matching. Simulation results show that after using this algorithm, when the initial position error of INS is 4.46km, the average matching error of each point still converges to 206.94m, achieving a high matching accuracy. Therefore, the direction of future research should be to combine the advantages of both TER-COM and ICCP methods, taking into account the tolerance for initial position errors and the requirements for positioning accuracy.

B. METHOD BASED ON STATISTICAL THEORY
The method based on statistical theory (BSTM) assumes that terrain features have some definite form of probability distribution and achieves the estimation of the probability distribution parameters by computing a specific criterion on the sample data. Statistics offers two approaches to this problem, namely, the maximum likelihood estimation method advocated by the frequentist school and the Bayesian estimation method advocated by the Bayesian school.
In the following, the principles of these two methods and their applications in terrain matching are described.

1) RECURSIVE BAYESIAN ESTIMATION
Bayesian theory was first proposed by the British mathematician Thomas Bayes in 1763. After more than 200 years of development, it has become a crucial theory in the field of statistics and plays an important role in artificial intelligence, information technology, finance, and other fields. The main idea of Bayesian estimation is to assume that the sample parameters follow a prior distribution, and then calculate the posterior distribution of the parameters based on the observed data [90].
Recursive Bayesian estimation (RBE) is a generalization of Bayesian estimation and an optimal estimation method for data. The application of RBE in filter estimation is extremely widespread. For linear systems, the Kalman filter is a special instance of RBE, which is an optimal estimation algorithm for linear systems. Many scholars at home and abroad have studied its application in terrain-matching navigation. Reference [91] proposed a navigation filter to update the state estimation by range or position measurement using the positioning system on the support ship. The algorithm reduced the maximum position estimation error from 1317m to 76m compared to DR. Xu et al. [92] proposed a method to add conditional constraints and confidence evaluation to the Extended Kalman Filter (EKF) to filter the localization values of the USBL, making the filtering results more robust and smooth. Li et al. [93] used a combined topographic matching algorithm based on ICCP and Kalman filtering, using underwater depth measurements, digitally stored topographic maps, and SINS output to generate position estimates for AUV. They used the output of TAN as the observation of the Kalman filter and used the filtering results to correct the output of SINS with feedback to obtain a more accurate indication of the ICCP algorithm trajectory. In the simulation experiments, the algorithm controls the matching position error always below 50m.
Since terrain features are nonlinear, terrain matching navigation can be considered as a nonlinear state-optimal estimation problem. Therefore, a nonlinear improvement of the Kalman filter is required to deal with nonlinear systems, such as terrain features. Currently, the most commonly used algorithm in terrain matching is the SITAN algorithm proposed by Sandia Laboratories in the 1970s [94]. It mainly consists of two stages, search mode, and tracking mode, to achieve coarse matching and fine matching of TANS, respectively. Essentially, it is a route correction method that uses the extended Kalman filter. However, it requires linearization of the terrain height matching model and a linear minimum variance estimation criterion. Since the TERCOM algorithm matches terrain elevation values only after measuring a distance, which is poor in real time, while SITAN is extremely sensitive to initial position errors although it is a real-time position estimation, [95] used the TERCOM algorithm for coarse matching in the search phase and the SITAN algorithm in the tracking phase to obtain a higher positioning accuracy. The search time of this algorithm is 0.281s, which is much smaller than the 1.023s of the inertial terrain aided navigation algorithm (BITAN) II. It consistently keeps the positioning error within 100m in the offline phase. Yang et al. [96] used the quadratic interpolation (QI) method and the VB method to derive the predicted error covariance matrix and the measurement noise matrix. The Kalman filter was improved by nonlinearization to estimate the state vector and the observation vector more accurately. Reference [97] proposed a mapping matching method based on the sliding window iterative nearest contour point (ICCP) algorithm, combining ICCP with SITAN to improve the matching navigation accuracy, which reduces the longitude error by nearly 60% and improves the overall matching accuracy by 50.3%.
However, these evolutionary filters may exhibit poor performance, such as filter scatter, when the topographic random linearization error or initial uncertainty is large. Particle filters (PF) or prime point filters (PMF) addresses this problem much better by adopting a more sophisticated approach to deal with highly nonlinear estimation problems [98]. Taking PF as an example, the basic principle is briefly explained [99]: the reference position and weights contained in each particle are initialized using the initial position of the vehicle and the error of the INS. The position of each particle at the next moment is calculated using the equation of motion, and the corresponding reference depth is obtained with the help of the DEM. The weights are corrected according to the difference between the measured depth and the reference depth. The position of the navigator is estimated by the particle swarm information. Repeating the above position recursion and weight update process, the estimation of full trajectory can be completed. Nordlund et al. [100] proposed the marginalized particle filter MPF. They used the PF to estimate the nonlinear part and the Kalman filter for its linear part. This greatly reduced the number of samples, making the applied MPF computationally easy to handle. Long et al. [101] proposed an algorithmic framework consisting of acquisition mode, tracking mode, and mode control logic. It first uses TERCOM for coarse localization to reduce the initial position error range, followed by fusion filters for precise localization. The fusion filters include particle filter and Kalman filter, which estimate the nonlinear and linear parts of the measurement data respectively, improving the real-time positioning accuracy of the system with an initial positioning error of 100m, which is smaller than the 150m of BITAN II. Claus et al. [102] used a particle filter to localize the AUV and calculated the particle weights based on the terrain complexity. In [103], the PF uses the computational solution of the EKF to propagate the particles during range measurement updates, and fuses the position updates back into the EKF solution. In the EKF and PF, state estimation and range measurements can be used as inputs to a deviation estimator for calculating speed and time synchronization deviations, and in experiments, the algorithm reduces the position root-mean-square error from 51m to 27m.

2) MAXIMUM LIKELIHOOD ESTIMATION
RBE is a continuous filter with large energy dissipation for multibeam sonar sensors, and the prior probability of the algorithm is difficult to solve with large initial uncertainty, VOLUME 11, 2023 47521 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.   which makes it unsuitable for long-duration diving missions of underwater vehicles. Therefore, the RBE degenerates to the maximum likelihood estimation (MLE). The maximum likelihood estimation utilizes discrete measurements of multiple sonar pulses to accomplish localization rather than long continuous measurements from a multibeam sonar sensor, which is simpler and more effective in practice [98].
Based on the above advantages, Chen [104] proposed an underwater terrain matching localization method based on maximum likelihood estimation, introduced the Fisher criterion to constrain the pseudo-peaks appearing in the likelihood function of flat terrain regions. He also proposed a Bayesian one-Fisher estimation method with the maximum directed feature parameter as the basis for quadratic discrimination, effectively eliminating the pseudo-peak points and enhancing 47522 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
the recognition rate of terrain features. Reference [98], which used multibeam echo sounder to obtain bathymetry data and used maximum likelihood estimation algorithm to match the real-time bathymetry with the analysis of the reference digital map, and the measurement beam of multibeam sonar in the simulation experiment was 10 × 21, and the Circular error was minimum 0.2m. In [32], the maximum likelihood estimation was used to achieve coarse matching by searching the surface terrain averaging factor to improve the matching efficiency under the assumption that the matching residuals obey Gaussian distribution. In the fine matching stage, inertia constraints were established to help determine the optimal set of matching reference points and match weight constraints were used to reduce the contribution of remote matching points in the overall cost function. The algorithm improved the matching accuracy in the root mean square errors (RMSE) by 66.54% over ICCP in the experiment on the lake.
The maximum likelihood estimation method always finds the maximum value of the likelihood function in a given topographic dataset, and uses it as the best estimate. When the dataset is large enough with certain regularization conditions, the maximum likelihood estimation has a Gaussian probability density function (PDF) and is asymptotically optimal. However, compared to Bayesian estimation, the maximum likelihood estimation lacks the regularization term and is often not the optimal estimate. Additionally, the integral calculation in RBE is complicated, making it difficult to compute posterior probabilities. A more appropriate approach is to use an approximate method, such as the maximum a posteriori probability (MAP), to find the posterior probability by approximating the computation of the RBE. This method, while less accurate than the traditional recursive Bayesian algorithm, greatly improves computational efficiency and is more compatible with practical engineering applications. Figure 9 shows the structural relationship between them.
Thus, in future research, we can introduce a regularization term equivalent to the prior probability in the maximum likelihood estimation to obtain the maximum a posteriori probability estimate. This approach has been successfully applied in current deep learning algorithms, which are very popular in artificial intelligence. In these algorithms, the loss function is the target of model optimization, and the optimized posterior estimate is obtained by introducing the regularization term in the learning process. Table 5 summarizes and compares all current mainstream terrain matching algorithms, providing a comprehensive reference for future research.

VI. CONCLUSION
This paper provides an overview of the framework of underwater terrain-aided navigation techniques, discussing the working principles of component units and comparing the characteristics of various research methods, including their advantages and disadvantages. Innovative discussions focus on the methods involved in terrain suitability analysis and terrain matching algorithms, using the ideas of dimensionality reduction to analyze terrain features and statistics to predict terrain matching results. These methods are summarized in a way that is useful for scholars who are new to this field, as well as those who have already studied it.
Although many TANS methods have been validated by numerous researchers, the following methods still require improvement: 1) Most of the current optimization for feature selection and matching algorithm parameters focuses on heuristic optimization algorithms. However, there are many other excellent optimization algorithms in the field of artificial intelligence, and it is interesting to try using these algorithms for feature selection and parameter tuning. 2) In terms of feature reduction for fitness analysis, feature selection methods are mostly used to form a composite feature with fixed parameters and then classify the region to be matched using a classifier. However, feature extraction methods are less applied. Future work should focus on feature dimensionality reduction using a data-driven approach based on the model for autonomous learning and variable parameter feature extraction methods. Combining this with nonlinear models for unsupervised classification will be a hot spot and difficult area for future research. 3) Topographic matching of terrain features using RBE requires integration of the posterior probability density, a process that is usually difficult to achieve. The maximum likelihood estimation lacks a regularization term compared to Bayesian estimation and is often not the most optimal estimate. Therefore, future research should try to introduce a regularization term equivalent to the prior probability in the maximum likelihood estimation to obtain the maximum a posteriori probability estimate.