Abstract
Terrain matching positioning is a promising method to overcome the problem that the inertial navigation error of the underwater vehicle accumulates over time. In the conventional terrain matching method, all measurement points are commonly used for matching and positioning. However, this method fails to be taken into a balanced consideration on both the computation complexity and the positioning accuracy. To reduce the computation and ensure the accuracy at the same time, an improved terrain matching method based on the gradient fitting is proposed in this paper. In the method, the gradient distributions of multiple terrain regions are firstly analyzed. Then, normal distribution is used to fit them, and according to the distribution, points with larger gradient values are selected as matching points. Finally, minimum absolute difference matching is chosen to match for positioning. Simulation results using multibeam sonar show that the improved terrain matching localization method not only reduces the computational complexity but also improves the accuracy of positioning.
1. Introduction
Underwater vehicles are becoming more and more popular in marine fields. At present, most underwater navigation systems are inertial navigation systems (INS). As the error of INS inevitably accumulates over time [1], many auxiliary measures need to be adopted, including Global Positioning System- (GPS-) aided methods [2, 3], acoustic navigation methods [4, 5], and Doppler velocity log- (DVL-) aided methods [6, 7]. However, the limitation of these techniques is that the navigation system must rely on either electromagnetic wave signal (e.g., GPS signal), artificial beacons (e.g., long baseline), complex mathematical models (e.g., flow prediction model), or the accessibility of the sea bottom (e.g., DVL bottom tracking). It would be difficult to satisfy the demands of efficiency, mobility, and invisibility at the same time.
As an autonomous underwater vehicle (AUV) is generally equipped with pressure measurement sensor and multibeam sonar, topographic information becomes an available and reliable resource in aided navigation. The terrain-aided navigation (TAN) system does not require any other external equipment, with high localization precision, favorable invisibility, complete autonomy, and no cumulative time error [8, 9]. TAN was first applied to cruise missiles and aircraft, such as the famous terrain contour matching (TERCOM) [10, 11] and Sandia inertial terrain-aided navigation (SITAN) systems [12]. In recent years, TAN’s far-reaching application potential in marine fields has aroused widespread concern of many scholars. Algorithms for underwater topographic matching positioning and navigation have been proposed one after another, such as iterated closest contour point- (ICCP-) based [13], cross-correlation-based [14], and maximum likelihood estimation-based [15] terrain navigation. The ICCP-based terrain navigation algorithm works by single-beam sounding sonar, whereas the terrain navigation method which is based on correlation and maximum likelihood estimation is applicable to not only single-beam sonar but also multibeam sonar. Compared with the single-beam sonar, multibeam sonar can measure a survey line perpendicular to the track at a time, obtain more detailed topographical information, and apply much more to terrain matching navigation [16].
The nonlinear Bayesian method is regarded as popular and state-of-the-art in TAN’s algorithms. Literature [17, 18] demonstrates the ability of the TAN system using this method to sustain a high navigation accuracy, without the need of any external position measurements. The TAN system is built based on the motion of an AUV and the measurement of the bottom depth. Two common recursive Bayesian filter equations are used recursively to update the distribution of the position offsets for each time step, and the particle filter may be used as a sampled-based approximation to solve the optimal Bayesian filter equations. However, since each multibeam measurement is performed on hundreds of points and calculated thousands of times, the conventional matching algorithm that traverses all positions and calculates with all measured points is computational expensive and unconducive to real-time positioning. Meanwhile, using uniform sampling points for matching may not guarantee the integrity of topographic information and it is inevitably affected by multibeam sonar measurement noise. In order to improve the performance of matching, terrain matching localization methods based on image processing are proposed in [19, 20]. However, the multibeam measurement noise may result in the distortion of natural terrain image, and it is uncertain which features of terrain would maximize the precise topographical information. Moreover, the feature extraction process is very time-consuming and greatly increases the extra computation.
To overcome the problem that the conventional terrain matching method fails to be taken into a balanced consideration in the computation complexity and the positioning accuracy, an effective sampling method of the matching points based on gradient distribution fitting is proposed in this paper, which is carried out from the strategy of selecting matching points. As the probability density function (PDF) of some common distributions (e.g., normal, gamma, and Weibull) is extensively used to characterize the normalized acquired data [21–23], the distribution which yields the lowest mean absolute error (MAE) would be chosen as the representation of the gradient data. The selection criteria of the matching points are proposed based on the quantile property of the fitting distribution, making the selection of the matching points more efficient and convenient. In this study, the fitting universality and computational efficiency of the normal distribution is proved, the criterion with optimum threshold is given, and according to the fitting distribution and best criterion, the measured depths with larger gradient values are selected as matching points for positioning. Different from the matching with uniform sampling points or all points, the matching based on gradient fitting can not only reduce the computation but also retain the main features of the terrain. Most importantly, the improved method reduces the influence of the adjacent or repetitive terrain and finally decreases the positioning error.
The rest of the paper is organized as follows. In Section 2, the conventional minimum absolute difference- (MAD-) based TERCOM algorithm and the essence of its matching error are analyzed. Section 3 provides the comparison between three different distributions and gives the selection criteria. An improved algorithm is proposed under the criteria. Then, the simulation results are given in Section 4. The optimal threshold of the selection criteria is given, and the difference between the conventional algorithms and the improved algorithm in both error performance and computation is analyzed. Finally, some conclusions are given in Section 5.
2. Conventional TERCOM Algorithm
Most modern underwater terrain localization algorithms, based on mean absolute difference (MAD) or mean square difference (MSD), process the data from the multibeam sonar in batches to derive an accurate estimate of the AUV’s position. This section will review the basic matching approach of the TERCOM method and analyze the essence of its limitation of localization performance.
2.1. Matching Approach
The depth values are measured as terrain data using multibeam sonar and then stored into batch files, called a measured profile, during a specific update period and are processed in batches using a correlation algorithm.
The digital elevation map (DEM) provided by bathymetric mapping sonar is a depth value function of discrete positions, represented by . Under a certain assumption of the map precision, the DEM can be converted into the matrix with interpolation function . where is the number of rows and is the number of columns of matrix .
In real-time navigation, the measured profile intervals provided by the multibeam sonar in the matching region are integrated and calculated. In the matching region, assume that position is the center and terrain data is the grouping of consecutive profiles. The measured terrain elevation matrix can be written as where and are, respectively, determined by the numbers of profiles and beams (). A reference matrix from the estimated position of the vehicle is determined in matrix , as shown in Figure 1 and the following process.
Consider ’s top left submatrix and its central
There is a relationship between the arbitrary position and the proposed position where is the matching interval. The reference matrix can be written as where is the width of the unit cell in the map and .
Thus, the MAD used for correlating the measured terrain elevation matrix with each size corresponding reference matrix is defined as follow.
Finally, find index of the minimum and determine as the position fix according to (4).
2.2. Applicability and Limitation
As the most widely known TAN algorithm, TERCOM is a correlation-oriented method using the collection of measured profiles rather than individual measurement recursive processing, which makes it insensitive to the initial position error and more stable to provide an accurate position fix to compensate for the navigation error of the INS. Meanwhile, the MAD correlation matching algorithm works as an efficient and robust way to determine the most relevant location in the DEM. With the increasing measurement range and accuracy of multibeam sonar, the position estimation of the AUV becomes more reliable owing to a large number of accurate calculation points. Since the DEM provides underwater terrain depth values in a matrix form, the batch processing works on the assumption that the AUV moves along a rectilinear trajectory in the matching region. This MAD-based TERCOM algorithm is beneficial due to high accuracy and simplicity of usage compared to the recursive algorithms like SITAN.
However, a false position fix of the AUV resulting from the wrong MAD estimation may occur while moving over a flat terrain or repetitive terrain which has no terrain uniqueness. Although the expansion of the matching region leads to comparable localization performance, the bias caused by the attitude change of the AUV would inevitably be enlarged and so is the computation time. In order to obtain more suitable solutions for these limitations, the essence of the correlated estimation errors is expected to be revealed as follows.
Most of the bias in the measured terrain elevation matrix is caused by ambient noise, in addition to the sensor accuracy, interpolation function, tidal fluctuations, and other factors. The measured terrain elevation matrix can be represented by where is the bias at point . According to (6), the MAD at position can be written as
Moreover, the relationship of the reference matrixes between position and adjacent position can be depicted as
In the flat terrain or repetitive terrain, even if and are considerably large, there is still and with the influence of , the value of the difference is not small enough to be distinguished from , which ultimately leads to a false position fix.
3. Improved Algorithm
In the following sections, we formulate a solution both for the improvement of localization accuracy and for the reduction of computing time utilizing calculating point selection based on depth gradient fitting. The formulation and process of the improved algorithm are described in Section 3.1 and 3.2, respectively. Section 3.3 lists the preconditions of the proposed algorithm. Section 3.4 realizes the depth gradient fitting and puts forward a feasible distribution, and Section 3.5 provides selection criteria of matching points based on the property of the proposed distribution.
3.1. Formulation
Consider the real-time position and its adjacent positions and . Due to the large overlap of the depth matrices in the adjacent positions, the reference matrixes are related in most of and , which can be depicted by
According to (8), since is relatively small compared to the depth, the absolute difference in the adjacent positions can be approximated as where and represent the absolute difference of the positions adjacent to in horizontal and vertical matching directions, respectively. Their average values determine the value of defined in (6) and ultimately influence the localization accuracy. Suppose if some specific , which have relatively large values of or , can be efficiently sampled out and calculated, the position bias would be reduced together with the interference from the adjacent position.
3.2. Process Flow
Therefore, the following processes are constructed for the improved algorithm.
Step 1. Construct the gradient matrixes. The -direction gradient matrix is defined as Similarly, the -direction gradient matrix is defined as
Step 2. Sample out specific points in the two matrixes. In , find out the points with the large gradient value (let us say for ), and according to (14a), add and to the point set . Similarly, in , according to (14b), add and to the point set . Finally, take the element union of two sets. where is the final set of points for matching.
Step 3. Match. The MAD is rewritten as
where represents the number of elements in the set .
In this way, the improved TERCOM algorithm is constructed by replacing (6) with steps 1–3.
3.3. Preconditions for Application
From the comparison between (13) and (14a) and (14b), or is proportional to the value of their own gradient matrix, only when , as shown in Figure 2. Due to that or calculated with the points in is relatively large, the side lobe interference in the real-time position would be comparably smaller, which ultimately improves the reliability of the position estimation.
Therefore, the algorithm works under the preconditions that (1)The matching direction of the algorithm must be consistent with the direction of gradient calculation; that is, the horizontal matching direction corresponds to the direction, and the vertical matching direction corresponds to the direction.(2)The intervals must be the same on matching and gradient calculation. It can be seen from (14a) and (14b) that . Certainly, the interval of the gradient calculation in (14a) and (14b) can be changed to satisfy the engineering precision demand.(3)The distribution of gradient values have been learnt, and according to the property of the fitting distribution, the corresponding point-sampling criteria have been set up, which makes the determination of the matching points more feasible and efficient.
3.4. Gradient Fitting Simulation
In order to determine the best distribution which fits the gradients of depth, the gradient data were compared with various standard distributions, i.e., normal, gamma, and Weibull. In the map, seven topographic regions are selected, numbered from 1 to 7, as shown in Figure 3.
In fact, the mean of the gradient values is approximately zero, and the gradient values in the matrixes defined in (14a) and (14b) are symmetrically distributed at zero. However, considering the properties of the gamma and Weibull distribution, a positive value is added to all the gradient values to ensure that they are distributed in the positive section.
For the sake of clarity, taking region no. 1, no. 2, and no. 4 for example, their landforms are shown in Figures 4 and 5 which depict their frequencies of the gradient data in both and direction as well as the PDF of normal, gamma, and Weibull distributions, which are given in (17), (18), and (19), respectively. where is the value of the positive-transformed gradient, is the mean, is the standard deviation, and are the scale parameters, and are shape parameters, and is the gamma function. These parameters are calculated for each region and tabulated in Tables 1 and 2.
(a) No. 1
(b) No. 2
(c) No. 4
(a) No. 1 direction fitting
(b) No. 1 direction fitting
(c) No. 2 direction fitting
(d) No. 2 direction fitting
(e) No. 4 direction fitting
(f) No. 4 direction fitting
Then, to test how well the distribution resembles the gradient values, MAE between the frequencies and PDF values is calculated in (20). where is the total number of gradient histograms, is the gradient level of each histogram, is the frequency, and is the probability density. The most fitting distribution has the lowest , and the results of the testing distributions are summarized in Tables 1 and 2.
3.5. Selection Criteria
As shown in the plots in Figure 5, three distribution curves are similar in shape and have a certain degree of fitting to the gradient values. As Figure 5(e) shows, when the number of negative gradients is more than the positive ones, the histogram is in a shape narrower on the left and wider on the right, and the gradients tend to obey the gamma distribution. In Figure 5(c), when the number of positive gradients is more than the negative ones, the shape of the histogram tends to be narrower on the right and wider on the left, and the gradients are inclined to obey the Weibull distribution. In Figure 5(b), when the number of positive and negative gradients is equal, the histogram is left-right symmetric, and the gradients fit the normal distribution best.
Although the result in Tables 1 and 2 reveals that Weibull distribution is generally smaller in , normal distribution appears to be more robust and universal for both positive and negative gradients. In addition, as is shown in Figures 6 and 7, 400 different regions are tested. The statistics in Table 3 shows that the of normal distribution has the smallest mean and variance.
Normal distribution is common in usage and simple in calculation as well. In conclusion, it is feasible to use normal distribution as the fitting curve of gradients and criteria of selecting matching points can be developed as follows.
First, estimate the normal distribution parameters.
According to the 3-sigma rule, it can be approximately assumed that all gradients are distributed in interval . The criterion can be given by
Under this criterion, in both and directions, 31.74% of the points with larger gradient values can be screened out as matching points through decision threshold . Finally, the set of matching points is determined in (15).
In many cases, a specific amount of matching points is required, and the criterion can be generalized as
Therefore, the threshold can be changed to approximately determine the amount of matching points.
4. Simulation
In this section, in order to evaluate the proposed algorithm, the estimation bias of three MAD-based TERCOM algorithms, respectively, matching with all points, uniform sampling points, and gradient-fitting sampling points, is firstly compared in different terrains, both flat and rough. Then, the influence of different thresholds in the generalized selection criterion on performance is studied in three different terrains. Finally, the optimal threshold of the selection criteria is given and the difference between the conventional algorithms and the improved algorithm in both error performance and computation is analyzed.
4.1. Localization Error Comparison
Assuming that 2025 of depth points are provided by the measurement in the matching region, the measured terrain elevation matrix in (2) with 45 rows and 45 columns can be built for matching. The purpose of MAD matching is to find the maximal from the matrix in (1) and determine the estimated position. The matching simulations of three sampling forms in a rough and flat terrain (shown in Figure 8) are shown in Figures 9 and 10, respectively.
(a) Rough terrain
(b) Flat terrain
(a) All 2025 points
(b) 1444 uniform sampling points
(c) 1432 gradient-fitting sampling points
(d) 841 uniform sampling points
(e) 795 gradient-fitting sampling points
(a) All 2025 points
(b) 1444 uniform sampling points
(c) 1397 gradient-fitting sampling points
(d) 841 uniform sampling points
(e) 804 gradient-fitting sampling points
In the rough terrain, Figure 9(a) shows the matching result of all 2025 points. The highest peak represents the estimated position as well as the true position. An interference position stands at the second-high peak, caused by the repetitive terrain. Figure 9(b) shows the matching result with 1444 points, uniformly sampled from measured data in the matching region. With the reduction in the amount of matching points, the peaks spread and the second-high peak becomes higher. While in Figure 9(c), matching with 1432 points sampled under the criterion defined in (22), the second-high peak collapses to a relatively low level and the uncertainty of the estimated position is narrowed. Figure 9(d) represents the matching result with 841 uniform sampling points. As shown in the figure, a new disturbance occurs as the price of the point reduction. Moreover, Figure 9(e) represents the matching result with 795 points sampled under the generalized criterion (23) with . Even if the number of matching points is further reduced, the certainty of the true position is not lowered compared to Figure 9(a). Therefore, in the rough terrain, the matching method under the gradient fitting criteria can effectively reduce the influence of repetitive terrain.
In the flat terrain, however, the risk of a false position fix increases as there are more interferential peaks, which is mainly caused by the similarity of nearby regions. Figure 10(a) shows the matching result of all 2025 points. Compared with Figure 10(a), the highest peak also represents the estimated position as well as true position, whereas the number of interference positions increases. Figure 10(b) represents the matching result of 1444 uniform sampling points, and there is a serious expansion of the highest peak, which ultimately results in the bias of the position estimation, while Figure 10(c), representing the result of 1397 points sampled under the criterion defined in (22), shows the same certainty of the position estimation as Figure 10(a). As the number of uniform sampling points further reduces to 841, Figure 10(d) shows a great deviation between estimated and true positions, while Figure 10(e) represents the matching result with 804 points sampled under the generalized criterion (23) with . Even though its number of matching points is comparable to Figure 10(d) and far less than Figure 10(a), its certainty and accuracy of the estimated position are higher than Figure 10(d) and comparable to Figure 10(a). Therefore, in the flat terrain, the matching method under the gradient fitting criteria can effectively reduce the influence of the adjacent terrain.
4.2. Optimum Threshold and Computation Complexity
The robustness of TERCOM based on gradient fitting is already confirmed in Section 4.1, and the performance of this method varies with Δ in (23), which determines the number of matching points and so is the computation. It is noteworthy that represents the matching with all points, and with the increase in Δ, the matching points decrease together with its compute. In this section, the localization performance under different is studied by Monte-Carlo simulation, which can be written as where represents the simulation times, is the true position, is the estimated position for the experiment, and the result is evaluated as matching error, indicating the mean deviation of the position estimation. The curves in Figure 11 reveal the movement of with increasing in three different terrains ().
(a) Region no. 3
(b) Region no. 5
(c) Region no. 6
All of the curves in Figure 11 present the same tendency of the matching error, monotonically decreasing in , reaching the minimum at , which is called min level, and monotonically increasing in . Moreover, when , the error approximately reaches the level at , which is called all-point level. Although the number of matching points decreases as increases, the matching error would become huge when or larger. Therefore, endlessly increasing to reduce computational complexity is not feasible.
However, due to the low fluctuation of terrain no. 3, the matching errors are generally large, and there would be a vast horizontal disparity between all-point level and min level, as shown in Figure 11(a). In Figure 11(b), terrain no. 5 has a higher degree of fluctuation, and the matching error is generally small, so the difference between the two levels is smaller. The performance in the moderate undulating terrain no. 6 is sort of in between, as shown in Figure 11(c).
In summary, the limitation of error reduction varies by terrain, but the choices of remain to be consistent. For the balance between matching error and computation complexity, the optimum choice of should be within (1,1.2), where the amount of matching points is relatively small, and the matching error is also close to the theoretical minimum. In this way, it can be approximated that (22) is the best criterion with optimum threshold. Under this criterion, the matching performance of the proposed gradient-fitting method is compared with the two conventional methods mentioned above. The results are listed in Table 4.
In Table 4, and represent the standard deviations of the gradient in the and directions, respectively, so as to determine the degree of terrain fluctuation. The result in the table shows that the higher fluctuation of the terrain is that, the lower the matching error becomes, the smaller the difference is between the three. However, the proposed matching algorithm has less matching error than the two conventional algorithms, regardless of the degree of fluctuation. Furthermore, compared with the all-point matching method, the proposed algorithm filters out an average of 609 points and retains an average of 1416 points for matching to achieve a larger error reduction. It explains that the filtered points have a certain interference to the position estimation and the remaining points are enough to characterize the topographic information, which radically reduces the estimation error as well as the compute. This also verifies the previous theoretical analysis and superiority of the proposed algorithm.
In terms of the compute, the complexity of the proposed algorithm is mainly determined by the amount of matching points and the number of reference positions . The amount of computation can be approximated as , that is, proportional to with a proportionality coefficient . Therefore, the change in computational complexity can be evaluated by the change of points with , which is shown in Figure 12.
The best is within (1,1.2), and the matching error is closest to minimum when . As is shown in Figure 9, when , the compute is reduced by , and when , the computation is reduced by . More importantly, in actual matching, would be larger, thus the reduction of the compute would be more significant.
5. Conclusion and Future Work
In this paper, an improved matching algorithm is proposed to overcome the shortcomings of the conventional matching algorithm. The simulation results show that the improved algorithm not only has better matching performance but also greatly reduces the computational cost. Meanwhile, the applicability of the normal distribution in gradient fitting is verified, corresponding criteria for sampling are given, and the algorithm gains the smallest localization error under the criterion with . Under this sampling criterion, the improved algorithm can better account for “computation” and “accuracy” to improve the matching positioning performance.
In fact, the proposed method can be a step in the current trend of navigation like using Bayesian methods to process dense terrain sensor measurements. Real-time positioning is hardly questionable with today’s computational power as well as the reduction in matching points. Also, compared to the conventional matching method, the proposed method gains higher positioning accuracy due to its unique sampling criterion. However, in our current study, there is no single evidence yet to prove the applicability and validity of the algorithm in mainstream navigation methods. Therefore, we will improve this part in our future work.
Data Availability
The terrain depth data used to support the findings of this study are available from the corresponding author upon request. All other data are shown above in the form of tables or charts.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was sponsored by the National Natural Science Foundation of China under Grants U1709203 and 41376103 and Research Foundation of Qingpu District 2018-14 (corresponding author: Tian Zhou) The authors are with the Underwater Acoustic Science and Technology Laboratory, College of Underwater Acoustic Engineering, Harbin Engineering University, Harbin, China (e-mail: [email protected]).