# Terrain Matching Localization for Underwater Vehicle Based on Gradient Fitting

Authors: Jiaqi Gao; Dongdong Peng; Tian Zhou; Tianhao Wang; Chao Xu
Journal: Journal of Sensors
 (2018)

---

## 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.


---

## Body

## 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 Section2, 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 byz=zx,y. Under a certain assumption of the map precision, the DEM can be converted into the M×N matrix with interpolation function H.
(1)Ai,j=Hzx,y,1≤i≤M,1≤j≤N,where M is the number of rows and N is the number of columns of matrix A.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 positionx0,y0 is the center and terrain data z′=z′x,y is the grouping of m consecutive profiles. The measured terrain elevation matrix can be written as
(2)Zx0y0′i,j=Hz′x,y,1≤i≤m,1≤j≤n,where m and n are, respectively, determined by the numbers of profiles and beams (m≪M,n≪N). A reference matrix from the estimated position of the vehicle is determined in matrix A, as shown in Figure 1 and the following process.Figure 1
Schematic of TERCOM.ConsiderA’s top left m×n submatrix and its central x1,y1
(3)Zx1y1i,j=Ai,j,1≤i≤m,1≤j≤n.There is a relationship between the arbitrary positionx,y and the proposed position x1,y1
(4)x=x1+sd′y=y1+td′,0≤s≤M−m,0≤t≤N−n,s,t∈N∗,where d′ is the matching interval. The reference matrix can be written as
(5)Zxyi,j=Zx1+sd1,y1+td2i,j=Ai+td′d,j+sd′d=Ai+tk,j+sk,where d is the width of the unit cell in the map and d′=kdk∈N∗.Thus, the MAD used for correlating the measured terrain elevation matrix with each size corresponding reference matrix is defined as follow.(6)Dxy=1mn∑i=1m∑j=1nZxyi,j−Zx0y0′i,j=1mn∑i=1m∑j=1nAi+tk,j+sk−Zx0y0′i,j.Finally, find indexs0,t0 of the minimum Dxy
(7)Dmins0,t0=min0≤s≤M−m,0≤t≤N−n.Dxyand determine x^0,y^0 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(8)Zx0y0′i,j=Zx0y0i,j+Ni,j,1≤i≤m,1≤j≤n,where Ni,j is the bias at point i,j. According to (6), the MAD at position x0,y0 can be written as
(9)Dx0y0=1mn∑i=1m∑j=1nNi,j.Moreover, the relationship of the reference matrixes between positionx0,y0 and adjacent position x0+α,y0+β can be depicted as
(10)limα⟶0,β⟶0Zx0+α,y0+βi,j=Zx0y0i,j.In the flat terrain or repetitive terrain, even ifα and β are considerably large, there is still
(11)Zx0+α,y0+βi,j≈Zx0y0i,j,and with the influence of Ni,j, the value of the difference Dx0y0 is not small enough to be distinguished from Dx0+α,y0+β, which ultimately leads to a false position fix.


## 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 byz=zx,y. Under a certain assumption of the map precision, the DEM can be converted into the M×N matrix with interpolation function H.
(1)Ai,j=Hzx,y,1≤i≤M,1≤j≤N,where M is the number of rows and N is the number of columns of matrix A.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 positionx0,y0 is the center and terrain data z′=z′x,y is the grouping of m consecutive profiles. The measured terrain elevation matrix can be written as
(2)Zx0y0′i,j=Hz′x,y,1≤i≤m,1≤j≤n,where m and n are, respectively, determined by the numbers of profiles and beams (m≪M,n≪N). A reference matrix from the estimated position of the vehicle is determined in matrix A, as shown in Figure 1 and the following process.Figure 1
Schematic of TERCOM.ConsiderA’s top left m×n submatrix and its central x1,y1
(3)Zx1y1i,j=Ai,j,1≤i≤m,1≤j≤n.There is a relationship between the arbitrary positionx,y and the proposed position x1,y1
(4)x=x1+sd′y=y1+td′,0≤s≤M−m,0≤t≤N−n,s,t∈N∗,where d′ is the matching interval. The reference matrix can be written as
(5)Zxyi,j=Zx1+sd1,y1+td2i,j=Ai+td′d,j+sd′d=Ai+tk,j+sk,where d is the width of the unit cell in the map and d′=kdk∈N∗.Thus, the MAD used for correlating the measured terrain elevation matrix with each size corresponding reference matrix is defined as follow.(6)Dxy=1mn∑i=1m∑j=1nZxyi,j−Zx0y0′i,j=1mn∑i=1m∑j=1nAi+tk,j+sk−Zx0y0′i,j.Finally, find indexs0,t0 of the minimum Dxy
(7)Dmins0,t0=min0≤s≤M−m,0≤t≤N−n.Dxyand determine x^0,y^0 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(8)Zx0y0′i,j=Zx0y0i,j+Ni,j,1≤i≤m,1≤j≤n,where Ni,j is the bias at point i,j. According to (6), the MAD at position x0,y0 can be written as
(9)Dx0y0=1mn∑i=1m∑j=1nNi,j.Moreover, the relationship of the reference matrixes between positionx0,y0 and adjacent position x0+α,y0+β can be depicted as
(10)limα⟶0,β⟶0Zx0+α,y0+βi,j=Zx0y0i,j.In the flat terrain or repetitive terrain, even ifα and β are considerably large, there is still
(11)Zx0+α,y0+βi,j≈Zx0y0i,j,and with the influence of Ni,j, the value of the difference Dx0y0 is not small enough to be distinguished from Dx0+α,y0+β, 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 Section3.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 positionx0,y0 and its adjacent positions x0±d′,y0 and x0,y0±d′. Due to the large overlap of the depth matrices in the adjacent positions, the reference matrixes are related in most of i and j, which can be depicted by
(12)Zx0±d′,y0i,j=Zx0y0i,j±k,Zx0,y0±d′i,j=Zx0y0i±k,j.According to (8), since Ni,j is relatively small compared to the depth, the absolute difference in the adjacent positions can be approximated as
(13)D1i,j=Zx0±d′,y0i,j−Zx0y0′i,j≈Zx0y0i,j±k−Zx0y0i,j,D2i,j=Zx0,y0±d′i,j−Zx0y0′i,j≈Zx0y0i±k,j−Zx0y0i,j,where D1 and D2 represent the absolute difference of the positions adjacent to x0,y0 in horizontal and vertical matching directions, respectively. Their average values determine the value of Dxy defined in (6) and ultimately influence the localization accuracy. Suppose if some specific i,j, which have relatively large values of D1 or D2, 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. TheX-direction gradient matrix is defined as
(14a)Gradxi,j=Zx0y0′i,2−Zx0y0′i,1d,j=1,Zx0y0′i,n−Zx0y0′i,n−1d,j=n,Zx0y0′i,j+1−Zx0y0′i,j−12d,else.
Similarly, theY-direction gradient matrix is defined as
(14b)Gradyi,j=Zx0y0′2,j−Zx0y0′1,jd,i=1,Zx0y0′m,j−Zx0y0′m−1,jd,i=m,Zx0y0′i+1,j−Zx0y0′i−1,j2d,else.Step 2.
Sample out specific points in the two matrixes. InGradx, find out the points with the large gradient value (let us say for i,j), and according to (14a), add i,j+1 and i,j−1 to the point set U1. Similarly, in Grady, according to (14b), add i+1,j and i+1,j to the point set U2. Finally, take the element union of two sets.
(15)U=U1⋃U2,where U is the final set of points for matching.Step 3.
Match. The MAD is rewritten as(16)Dxy=1sizeU∑i,j∈UZxyi,j−Zx0y0′i,j,where sizeU represents the number of elements in the set U.
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), D1 or D2 is proportional to the value of their own gradient matrix, only when d′=2d, as shown in Figure 2. Due to that D1 or D2 calculated with the points in U 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.Figure 2
Illustration of the preconditions.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 theX direction, and the vertical matching direction corresponds to the Y direction.
(2)
The intervals must be the same on matching and gradient calculation. It can be seen from (14a) and (14b) that d′=2d. 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 Figure3.Figure 3
DEM and seven topographic regions.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 Figures4 and 5 which depict their frequencies of the gradient data in both X and Y direction as well as the PDF of normal, gamma, and Weibull distributions, which are given in (17), (18), and (19), respectively.
(17)pXx;μ,σ=12πσexp−x−μ22σ2,(18)pXx;a,b=baΓaxa−1exp−bx,x>0,0,x≤0,(19)pXx;λ,k=kλxλk−1exp−xλk,x≥0,0,x<0,where x is the value of the positive-transformed gradient, μ is the mean, σ is the standard deviation, b and λ are the scale parameters, a and k are shape parameters, and Γa is the gamma function. These parameters are calculated for each region and tabulated in Tables 1 and 2.Figure 4
Landforms of three regions.
(a)
No. 1
(b)
No. 2
(c)
No. 4Figure 5
Results of gradient fitting.
(a)
No. 1X direction fitting
(b)
No. 1Y direction fitting
(c)
No. 2X direction fitting
(d)
No. 2Y direction fitting
(e)
No. 4X direction fitting
(f)
No. 4Y direction fittingTable 1
Parameters and MAE ofX direction gradient fitting.
Region
1
2
3
4
5
6
7
Parameters
μ
1.0742
0.9434
0.9828
1.0229
0.8960
1.0427
0.8773
σ
0.0874
0.1671
0.1010
0.1324
0.1886
0.1289
0.2354
a
69.602
31.886
41.458
78.289
32.975
94.024
13.608
b
64.796
33.800
42.185
76.535
36.802
90.178
15.512
λ
1.1293
1.0034
1.0402
1.0752
0.9563
1.0901
0.9528
k
9.4194
7.8748
8.0858
8.6665
7.1263
10.1016
5.4273
dk
Normal
0.2025
0.3192
0.3061
0.2108
0.3784
0.1589
0.3657
Gamma
0.2253
0.3953
0.3702
0.1852
0.4458
0.1760
0.4407
Weibull
0.2180
0.2214
0.2635
0.3691
0.3196
0.2533
0.3002Table 2
Parameters and MAE ofY direction gradient fitting.
Region
1
2
3
4
5
6
7
Parameters
μ
1.0119
0.9583
0.9630
1.0353
0.9164
1.0185
0.9420
σ
0.0773
0.1339
0.0950
0.1160
0.1457
0.1312
0.1643
a
79.846
52.816
43.202
78.720
59.301
81.237
16.181
b
78.904
55.114
44.861
76.036
64.714
79.759
17.177
λ
1.0612
1.0105
1.0208
1.0895
0.9659
1.0690
1.0268
k
9.5668
8.9917
7.6033
8.1490
8.9526
8.8259
5.0442
dk
Normal
0.1402
0.2812
0.2700
0.2261
0.3005
0.1550
0.2183
Gamma
0.1529
0.3425
0.3243
0.1982
0.3566
0.1405
0.2740
Weibull
0.2423
0.2199
0.2594
0.4166
0.2831
0.3053
0.1812Then, to test how well the distribution resembles the gradient values, MAE between the frequencies and PDF values is calculated in (20).
(20)dk=1M∑t=1Mht−pt,where M is the total number of gradient histograms, t is the gradient level of each histogram, ht is the frequency, and pt is the probability density. The most fitting distribution has the lowest dk, and the results of the testing distributions are summarized in Tables 1 and 2.


### 3.5. Selection Criteria

As shown in the plots in Figure5, 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 Tables1 and 2 reveals that Weibull distribution is generally smaller in dk, 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 dk of normal distribution has the smallest mean and variance.Figure 6
X-direction gradient fitting in 400 regions.Figure 7
Y-direction gradient fitting in 400 regions.Table 3
Statistics ofdk in 400 regions.
Normal
Gamma
Weibull
X direction
Mean
0.2758
0.2916
0.3202
Variance
0.1085
0.1185
0.1248
Y direction
Mean
0.2712
0.5069
0.3313
Variance
0.1099
0.3152
0.1443Normal 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.(21)μ=1mn∑i=1m∑j=1nGradi,j,σ2=1mn∑i=1m∑j=1nGradi,j−μ2.According to the 3-sigma rule, it can be approximately assumed that all gradients are distributed in intervalμ−3σ,μ+3σ. The criterion can be given by
(22)PGradi,j−μ>σ=31.74%.Under this criterion, in bothX and Y 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(23)PGradi,j−μ>Δσ=PΔ.Therefore, the thresholdΔσ can be changed to approximately determine the amount of matching points.


## 3.1. Formulation

Consider the real-time positionx0,y0 and its adjacent positions x0±d′,y0 and x0,y0±d′. Due to the large overlap of the depth matrices in the adjacent positions, the reference matrixes are related in most of i and j, which can be depicted by
(12)Zx0±d′,y0i,j=Zx0y0i,j±k,Zx0,y0±d′i,j=Zx0y0i±k,j.According to (8), since Ni,j is relatively small compared to the depth, the absolute difference in the adjacent positions can be approximated as
(13)D1i,j=Zx0±d′,y0i,j−Zx0y0′i,j≈Zx0y0i,j±k−Zx0y0i,j,D2i,j=Zx0,y0±d′i,j−Zx0y0′i,j≈Zx0y0i±k,j−Zx0y0i,j,where D1 and D2 represent the absolute difference of the positions adjacent to x0,y0 in horizontal and vertical matching directions, respectively. Their average values determine the value of Dxy defined in (6) and ultimately influence the localization accuracy. Suppose if some specific i,j, which have relatively large values of D1 or D2, 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. TheX-direction gradient matrix is defined as
(14a)Gradxi,j=Zx0y0′i,2−Zx0y0′i,1d,j=1,Zx0y0′i,n−Zx0y0′i,n−1d,j=n,Zx0y0′i,j+1−Zx0y0′i,j−12d,else.
Similarly, theY-direction gradient matrix is defined as
(14b)Gradyi,j=Zx0y0′2,j−Zx0y0′1,jd,i=1,Zx0y0′m,j−Zx0y0′m−1,jd,i=m,Zx0y0′i+1,j−Zx0y0′i−1,j2d,else.Step 2.
Sample out specific points in the two matrixes. InGradx, find out the points with the large gradient value (let us say for i,j), and according to (14a), add i,j+1 and i,j−1 to the point set U1. Similarly, in Grady, according to (14b), add i+1,j and i+1,j to the point set U2. Finally, take the element union of two sets.
(15)U=U1⋃U2,where U is the final set of points for matching.Step 3.
Match. The MAD is rewritten as(16)Dxy=1sizeU∑i,j∈UZxyi,j−Zx0y0′i,j,where sizeU represents the number of elements in the set U.
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), D1 or D2 is proportional to the value of their own gradient matrix, only when d′=2d, as shown in Figure 2. Due to that D1 or D2 calculated with the points in U 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.Figure 2
Illustration of the preconditions.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 theX direction, and the vertical matching direction corresponds to the Y direction.
(2)
The intervals must be the same on matching and gradient calculation. It can be seen from (14a) and (14b) that d′=2d. 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 Figure3.Figure 3
DEM and seven topographic regions.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 Figures4 and 5 which depict their frequencies of the gradient data in both X and Y direction as well as the PDF of normal, gamma, and Weibull distributions, which are given in (17), (18), and (19), respectively.
(17)pXx;μ,σ=12πσexp−x−μ22σ2,(18)pXx;a,b=baΓaxa−1exp−bx,x>0,0,x≤0,(19)pXx;λ,k=kλxλk−1exp−xλk,x≥0,0,x<0,where x is the value of the positive-transformed gradient, μ is the mean, σ is the standard deviation, b and λ are the scale parameters, a and k are shape parameters, and Γa is the gamma function. These parameters are calculated for each region and tabulated in Tables 1 and 2.Figure 4
Landforms of three regions.
(a)
No. 1
(b)
No. 2
(c)
No. 4Figure 5
Results of gradient fitting.
(a)
No. 1X direction fitting
(b)
No. 1Y direction fitting
(c)
No. 2X direction fitting
(d)
No. 2Y direction fitting
(e)
No. 4X direction fitting
(f)
No. 4Y direction fittingTable 1
Parameters and MAE ofX direction gradient fitting.
Region
1
2
3
4
5
6
7
Parameters
μ
1.0742
0.9434
0.9828
1.0229
0.8960
1.0427
0.8773
σ
0.0874
0.1671
0.1010
0.1324
0.1886
0.1289
0.2354
a
69.602
31.886
41.458
78.289
32.975
94.024
13.608
b
64.796
33.800
42.185
76.535
36.802
90.178
15.512
λ
1.1293
1.0034
1.0402
1.0752
0.9563
1.0901
0.9528
k
9.4194
7.8748
8.0858
8.6665
7.1263
10.1016
5.4273
dk
Normal
0.2025
0.3192
0.3061
0.2108
0.3784
0.1589
0.3657
Gamma
0.2253
0.3953
0.3702
0.1852
0.4458
0.1760
0.4407
Weibull
0.2180
0.2214
0.2635
0.3691
0.3196
0.2533
0.3002Table 2
Parameters and MAE ofY direction gradient fitting.
Region
1
2
3
4
5
6
7
Parameters
μ
1.0119
0.9583
0.9630
1.0353
0.9164
1.0185
0.9420
σ
0.0773
0.1339
0.0950
0.1160
0.1457
0.1312
0.1643
a
79.846
52.816
43.202
78.720
59.301
81.237
16.181
b
78.904
55.114
44.861
76.036
64.714
79.759
17.177
λ
1.0612
1.0105
1.0208
1.0895
0.9659
1.0690
1.0268
k
9.5668
8.9917
7.6033
8.1490
8.9526
8.8259
5.0442
dk
Normal
0.1402
0.2812
0.2700
0.2261
0.3005
0.1550
0.2183
Gamma
0.1529
0.3425
0.3243
0.1982
0.3566
0.1405
0.2740
Weibull
0.2423
0.2199
0.2594
0.4166
0.2831
0.3053
0.1812Then, to test how well the distribution resembles the gradient values, MAE between the frequencies and PDF values is calculated in (20).
(20)dk=1M∑t=1Mht−pt,where M is the total number of gradient histograms, t is the gradient level of each histogram, ht is the frequency, and pt is the probability density. The most fitting distribution has the lowest dk, and the results of the testing distributions are summarized in Tables 1 and 2.


## 3.5. Selection Criteria

As shown in the plots in Figure5, 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 Tables1 and 2 reveals that Weibull distribution is generally smaller in dk, 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 dk of normal distribution has the smallest mean and variance.Figure 6
X-direction gradient fitting in 400 regions.Figure 7
Y-direction gradient fitting in 400 regions.Table 3
Statistics ofdk in 400 regions.
Normal
Gamma
Weibull
X direction
Mean
0.2758
0.2916
0.3202
Variance
0.1085
0.1185
0.1248
Y direction
Mean
0.2712
0.5069
0.3313
Variance
0.1099
0.3152
0.1443Normal 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.(21)μ=1mn∑i=1m∑j=1nGradi,j,σ2=1mn∑i=1m∑j=1nGradi,j−μ2.According to the 3-sigma rule, it can be approximately assumed that all gradients are distributed in intervalμ−3σ,μ+3σ. The criterion can be given by
(22)PGradi,j−μ>σ=31.74%.Under this criterion, in bothX and Y 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(23)PGradi,j−μ>Δσ=PΔ.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 1/D 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.Figure 8
Two types of terrain for testing.
(a)
Rough terrain
(b)
Flat terrainFigure 9
Matching performance with different matching points in rough 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 pointsFigure 10
Matching performance with different matching points in a flat terrain.
(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 pointsIn the rough terrain, Figure9(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 Δ=1.5. 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. Figure10(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 Δ=1.5. 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 Section4.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 Δ=0 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
(24)e=1L∑l=1Lx^0l−x02+y^0l−y02,where L represents the simulation times, x0,y0 is the true position, x^0l,y^0l is the estimated position for the lth experiment, and the result e is evaluated as matching error, indicating the mean deviation of the position estimation. The curves in Figure 11 reveal the movement of e with increasing Δ in three different terrains (L=5000).Figure 11
Error curves of three terrains.
(a)
Region no. 3
(b)
Region no. 5
(c)
Region no. 6All of the curves in Figure11 present the same tendency of the matching error, monotonically decreasing in 0,1, reaching the minimum at Δ=1, which is called min level, and monotonically increasing in 1,3. Moreover, when Δ∈1.2,1.4, the error approximately reaches the level at Δ=0, which is called all-point level. Although the number of matching points decreases as Δ increases, the matching error would become huge when Δ=2 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 Figure11(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.Table 4
Matching error comparison of three methods.
Region
1
2
3
4
5
6
7
Fluctuation
σx
0.0874
0.1671
0.1010
0.1324
0.1886
0.1289
0.2354
σy
0.0773
0.1339
0.0950
0.1160
0.1457
0.1312
0.1643
Number of matching points
All
2025
2025
2025
2025
2025
2025
2025
Uniform
1444
1444
1521
1444
1444
1521
1444
Grad-fit
1385
1435
1450
1376
1389
1452
1428
e (m)
All
2.8293
0.6474
1.2387
1.2078
0.1680
2.1053
0.2680
Uniform
6.4865
2.0072
3.5304
2.2676
0.3186
6.8306
0.5573
Grad-fit
2.6418
0.5448
1.1272
1.0042
0.1106
1.9543
0.2153In Table4, σx and σy represent the standard deviations of the gradient in the X and Y 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 pointsJ and the number of reference positions K. The amount of computation can be approximated as J×K, that is, proportional to J with a proportionality coefficient K. Therefore, the change in computational complexity can be evaluated by the change of points with Δ, which is shown in Figure 12.Figure 12
The regulation of point reduction.The bestΔ is within (1,1.2), and the matching error is closest to minimum when Δ=1. As is shown in Figure 9, when Δ=1, the compute is reduced by 2025−1416×K=609K, and when Δ=1.2, the computation is reduced by 2025−1150×K=875K. More importantly, in actual matching, K would be larger, thus the reduction of the compute would be more significant.


## 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 1/D 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.Figure 8
Two types of terrain for testing.
(a)
Rough terrain
(b)
Flat terrainFigure 9
Matching performance with different matching points in rough 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 pointsFigure 10
Matching performance with different matching points in a flat terrain.
(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 pointsIn the rough terrain, Figure9(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 Δ=1.5. 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. Figure10(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 Δ=1.5. 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 Section4.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 Δ=0 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
(24)e=1L∑l=1Lx^0l−x02+y^0l−y02,where L represents the simulation times, x0,y0 is the true position, x^0l,y^0l is the estimated position for the lth experiment, and the result e is evaluated as matching error, indicating the mean deviation of the position estimation. The curves in Figure 11 reveal the movement of e with increasing Δ in three different terrains (L=5000).Figure 11
Error curves of three terrains.
(a)
Region no. 3
(b)
Region no. 5
(c)
Region no. 6All of the curves in Figure11 present the same tendency of the matching error, monotonically decreasing in 0,1, reaching the minimum at Δ=1, which is called min level, and monotonically increasing in 1,3. Moreover, when Δ∈1.2,1.4, the error approximately reaches the level at Δ=0, which is called all-point level. Although the number of matching points decreases as Δ increases, the matching error would become huge when Δ=2 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 Figure11(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.Table 4
Matching error comparison of three methods.
Region
1
2
3
4
5
6
7
Fluctuation
σx
0.0874
0.1671
0.1010
0.1324
0.1886
0.1289
0.2354
σy
0.0773
0.1339
0.0950
0.1160
0.1457
0.1312
0.1643
Number of matching points
All
2025
2025
2025
2025
2025
2025
2025
Uniform
1444
1444
1521
1444
1444
1521
1444
Grad-fit
1385
1435
1450
1376
1389
1452
1428
e (m)
All
2.8293
0.6474
1.2387
1.2078
0.1680
2.1053
0.2680
Uniform
6.4865
2.0072
3.5304
2.2676
0.3186
6.8306
0.5573
Grad-fit
2.6418
0.5448
1.1272
1.0042
0.1106
1.9543
0.2153In Table4, σx and σy represent the standard deviations of the gradient in the X and Y 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 pointsJ and the number of reference positions K. The amount of computation can be approximated as J×K, that is, proportional to J with a proportionality coefficient K. Therefore, the change in computational complexity can be evaluated by the change of points with Δ, which is shown in Figure 12.Figure 12
The regulation of point reduction.The bestΔ is within (1,1.2), and the matching error is closest to minimum when Δ=1. As is shown in Figure 9, when Δ=1, the compute is reduced by 2025−1416×K=609K, and when Δ=1.2, the computation is reduced by 2025−1150×K=875K. More importantly, in actual matching, K 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Δ=1. 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.



---

*Source: 3717430-2018-12-30.xml*