2.1. Underlying principle of MAP

The basic idea underlying MAP is use of measurement data to estimate the system status. It is similar to Maximum Likelihood Estimation (MLE) as proposed by Fisher (Geisser, Reference Geisser, Kotz and Johnson1992), and can be seen as regularised MLE (Ding and Xiao, Reference Ding and Xiao2014).

If the posterior Probability Density Function (PDF) p(x|Z) of the estimated amount x is given, the estimated value ${\bi{\hat x}}$ of MAP is the value of x which maximises the posterior PDF. From the Bayesian formula, if p(Z|x) is known, the posterior PDF of the system status can be calculated using Equation (1):

(1)$$p\lpar {\bi Z}\vert {\bi x}\rpar =\displaystyle{{p\lpar {\bi Z}\vert {\bi x}\rpar p\lpar {\bi x}\rpar }\over{p\lpar {\bi Z}\rpar }} $$

where p(x) is the prior PDF of x, p(Z) is the PDF of measurement data Z and p(Z|x) can be obtained via calculation or experience. Let ${\bi{\hat x}}_{MAP}$ be the estimated value of MAP, taking the natural logarithms of both sides of Equation (1) results in Equation (2):

(2)$$\ln p\lpar {\bi x}\vert {\bi Z}\rpar =\ln p\lpar {\bi Z}\vert {\bi x}\rpar +\ln p\lpar {\bi x}\rpar -\ln p\lpar {\bi Z}\rpar $$

The physical significance of MAP is as follows: if measurement data Z are determinate, the probability that x is in the neighbourhood of ${\bi{\hat x}}_{MAP}$ is higher than that of other areas. Obviously, MAP should satisfy the following equation:

(3)$$\displaystyle{{\partial}\over{\partial {\bi x}}} \ln p\lpar {\bi x}\vert{\bi Z}\rpar \vert_{{\bi x}={\bi{\hat x}}_{MAP}}=0 $$

Equation (3) is called the posterior equation, which we can solve to obtain ${\bi{\hat x}}_{MAP}$. As p(Z|x) and x are independent, combining Equations (2) and (3), the posterior equation can be rewritten as follows:

(4)$$\left.\left[\displaystyle{{\partial}\over{\partial {\bi x}}}\ln p\lpar {\bi Z}\vert{\bi x}\rpar +\displaystyle{{\partial}\over{\partial {\bi x}}}\ln p\lpar {\bi x}\rpar \right]\right\vert_{{\bi x}={\bi{\hat x}}_{MAP}}=0 $$

If there is no prior information of x – that is, the probability value of x is equivalent – the prior PDF p(x) can be considered as a normal distribution with variance tending to infinity, and p(x) is as follows:

(5)$$p\lpar {\bi x}\rpar =\displaystyle{{1}\over{\lpar 2\pi\rpar ^{n/2}\vert {\bi P}_{{\bi x}}\vert ^{1/2}}} \exp\left\{-\displaystyle{{1}\over{2}}\lpar {\bf x}-\mu\rpar^T{\bi P}_{{\bi x}}^{-1}\lpar {\bi x}-\mu\rpar \right\}$$

where n is the number of measured data, P_x is the covariance matrix of the measured error and μ is the expected value. As the variance tends to infinity P_x → ∞, ${\bi P}_{{\bi x}}^{-1} \to 0$. Taking the natural logarithms of both sides of Equation (5):

(6)$$\eqalign{\displaystyle{{\partial}\over{\partial {\bi x}}}\ln p\lpar {{\bi x}} \rpar & =-\displaystyle{{1}\over{2}}\ln \vert {{\bi P}_{x}} \vert -\ln \lpar {{\bi P}_{{\bi x}}^{-1}}\rpar -{\bi P}_{{\bi x}}^{-1} \lpar {{\bi x}-\mu } \rpar \cr & =-{\bi P}_{{\bi x}}^{-1} \lpar {{\bi x}-\mu } \rpar \to 0} $$

At this point, ${\bi{\hat x}}_{MAP}$ satisfies the following condition:

(7)$$\displaystyle{{\partial}\over{\partial{\bi x}}} \ln p\lpar {\bi z}\vert{\bi x}\rpar \vert_{{\bi x}={{{\bi{\hat x}}}_{MAP}}}=0 $$

According to the basics of MLE, Equation (7) is the judging condition of MLE. Therefore, if there is no prior information of x, MAP will degenerate into MLE. In other words, MLE is a special MAP. However, as MAP considers the a priori statistical information of the estimated amount, MAP is better than MLE.

2.2. UTPM model

The vertical coordinates of an AUV can be measured using an acoustic altimeter and a water pressure sensor in real time, resulting in no cumulative error. Thus, the vertical position need not be considered in terrain matching. Based on the above analysis, in this paper, the positioning accuracy considered is horizontal accuracy.

Without considering the coordinate changes of the AUV in the vertical plane, and assuming the navigation system performs a data update at each sampling time, the UTPM model is stated as follows:

(8)$${\bi X}_{t+1} = {\bi X}_t+{\bf U}_t+{\bi v}_t $$

(9)$${\bi Y}_{t} = {\bi H}_{t} \lpar {\bi x}_{t}\rpar +{\bi E}_{t} $$

where X_t is the INS position at time t, U_t is the position change between adjacent positioning times, ν_t is the error of INS, Y_t is the real-time sounding data at time t, H_t is the local DTM data at x_t and E_t is the sounding error, which is defined as Gaussian white noise. To simplify the analysis, Equation (9) can be rewritten as follows:

(10)$${\bi y}_{t} ={\bi h}_{t} \lpar {\bi x}_{t} \rpar +{\bi e}_{t} $$

Assuming that the measurement errors are independent of the water depth, and the errors are independent and identically distributed, the possibility of measuring sequence y_t at x_t is p(y_t|x_t) = p(e_t), and the prior PDF can be expressed as Equation (11):

(11)$$p\lpar {{\bi y}_{t} \vert {\bi x}_{t} } \rpar =p\lpar {{\bi y}_{t} -{\bi h}_{t} \lpar {{\bi x}_{t} } \rpar } \rpar =\displaystyle{{1}\over{\left({\sqrt {2\pi} \sigma_{e} } \right)^{n}}}\exp \left[{-\displaystyle{{1}\over{2\sigma_{e}^{2}}}\lpar {{\bi y}_{t} -{\bi h}\lpar {{\bi x}_{t} } \rpar } \rpar^{2}} \right]$$

Assuming the positioning error of INS obeys v_t ~ N(0, Σ), p(x_t) can be expressed as follows:

(12)$$p\lpar {\bi x}_t\rpar =\displaystyle{{1}\over{2\pi\sqrt{\det\, {\bi \Sigma}}}}\exp\left\{\displaystyle{{1}\over{2}}\lpar {\bi x}_t-{\bi X}\rpar {\bi \Sigma}^{-1}\lpar {\bi x}_t-{\bi X}\rpar^T\right\}$$

In Equation (12), X = (X, Y) is the indicating position of INS. Substituting Equations (11) and (12) into Equation (4), we obtain the posterior equation used for terrain matching positioning. The positioning point which satisfies the a posteriori equation is the MAP estimated position of terrain matching.

2.3. Pseudo-anchor points phenomenon

In principle, the positioning point which satisfies the a posteriori equation can always be calculated, and the true position can be fixed. Due to the complex underwater environment and the interference of measurement errors, the position is not always unique and the posterior PDF p(x_t|y_t) has multiple peaks. The peak which represents the actual position of the AUV is called the true positioning. The peaks outside the true positioning are called pseudo-peaks, and the corresponding anchor points are pseudo-anchor points. Pseudo-anchor points are a common phenomenon in underwater terrain matching positioning and negatively influence the matching result. This influence is particularly evident in flat areas and causes matching accuracy to decrease or even fail.

Let us analyse the existence characteristics of a pseudo-peak. According to Equation (9), we obtain the following:

(13)$${\bi y}_{t} =\lsqb {\bi h}\lpar x_{t} \rpar +{\bi E}_{t} \rsqb =\left\{{{\matrix{ {h_{t\comma 1} =h_{1} \lpar x_{t} \rpar +e_{1} } \cr \vdots \cr {h_{t\comma k} =h_{k} \lpar x_{t} \rpar +e_{k} } \cr \vdots \cr {h_{t\comma n} =h_{n} \lpar x_{t} \rpar +e_{n} } \cr }}} \right\}$$

where e _k is the measurement error of the k-th real-time sounding data and h(x _k) is the extracted depth data from DTM at x _k. Suppose that x ₀ is the real matching position, x ₁ is the pseudo-anchor point, and y₀ and y₁ are the sounding sequence at x ₀ and x ₁, respectively. Let Δ y = [Δ y ₁, Δ y ₂, …, Δ y _n] ^T be the difference sequence between y₀ and y₁; thus, y₁ = y₀ − Δ y. Further, assuming that the measurement errors are independent and identically distributed, we obtain Equation (14):

(14)$${\bi y}_{1} ={\bi y}_{0} -\Delta {\bi y}=\lsqb {\bi h}\lpar x_{0} \rpar -\Delta {\bi h}+{\bi E}_{k} \rsqb =\left\{{{\matrix{ {h_{t\comma 1} =h_{1} \lpar x_{t} \rpar -\Delta h_{1} +e_{1} } \cr \vdots \cr {h_{t\comma k} =h_{k} \lpar x_{t} \rpar -\Delta h_{k} +e_{k} } \cr \vdots \cr {h_{t\comma n} =h_{n} \lpar x_{t} \rpar -\Delta h_{n} +e_{n} } \cr } }} \right\}$$

$p\lpar {\bi x}\vert {\bi y}\rpar =\displaystyle{{p\lpar {\bi x}\vert {\bi y}\rpar p\lpar {\bi x}\rpar }\over{p\lpar {\bi y}\rpar }}\propto p\lpar {\bi y}\vert {\bi x}\rpar $. Let J(x_t) represent the power term of the exponential term part in Equation (11), then:

(15)$${\bi J}\lpar {\bi x}_t\rpar =\lpar {\bi y}_t-{\bi h}\lpar {\bi x}_t\rpar \rpar ^2=\sum_{i=1}^{n}\lpar y_{t\comma i}-h_i\lpar x_t\rpar \rpar ^2 $$

Thus:

(16)$$\eqalign{J\lpar x_{1} \rpar -J\lpar x_{0} \rpar &=\lpar {\bi y}_{1} -{\bi h}\lpar x_{1} \rpar \rpar ^{{\bf T}}\lpar {\bi y}_{1} -{\bi h}\lpar x_{1} \rpar \rpar -\lpar {\bi y}_{0} -{\bi h}\lpar x_{0} \rpar \rpar ^{{\bf T}}\lpar {\bi y}_{0} -{\bi h}\lpar x_{0} \rpar \rpar \cr &=\sum\limits_{k=1}^n {\lpar e_{k} -\Delta h_{k} \rpar ^{2}} -\sum\limits_{k=1}^n {e_{k}^{2}} \cr &=\sum\limits_{k=1}^n {\lpar \Delta h_{k}^{2}-2} e_{k} \Delta h_{k} \rpar \cr } $$

As can be seen from Equation (16), when $\sum_{k=1}^{1}\lpar \Delta h_{k}^{2}-2e_{k}\Delta h_{k}\rpar \gg 0$, there is no pseudo-peak on the posterior PDF and x ₀ is the best matching position. In the presence of pseudo-peaks, with the real-time sounding data increasing, the value of $\sum_{k=1}^{n} \lpar \Delta h_{k}^{2}-2 e_{k} \Delta h_{k}\rpar $ will gradually increase. Thus, theoretically, when the sounding data are sufficient, there will be no pseudo-positioning. However, owing to the error of the measurement terrain and terrain similarity, the pseudo-anchor points will not be adequately removed by increasing the amount of sounding data. Further, too much sounding data results in heavy calculation, which affects the real-time performance of the terrain matching algorithm.

The pseudo-anchor points are derived from terrain similarity and measurement error. If the errors are independent and identically distributed, the pseudo-anchor points are decided by the sum of all measurement data errors. The terrain matching algorithm is based on the hypothesis that the effect of each set of measurement data for terrain matching is the same. In this hypothesis, the overall error is the average of individual errors, and differences in individual errors are not considered. When terrain features are flat, the differences in the overall errors are not obvious. However, with increasing measurement data, the effect of ‘error averaging’ also increases, which means that the result of terrain matching positioning cannot necessarily improve with increasing quantity of measurement data. Thus, if data with large error can be removed effectively, the effect of ‘error averaging’ can be effectively suppressed, and the terrain matching performance can be improved. To achieve this, a pseudo-anchor points discriminant method is proposed in this paper. In the proposed method, pseudo-anchor points are effectively removed by screening the data that have smaller errors in the measurement data, then discriminating all matching points.

3.1. PCNN model

By studying the cerebral cortex of mammals, Eckhorn et al. (Reference Eckhorn, Reitboeck, Arndt and Dicke1990) established a neuronal conduction characteristics model for the visual area, which eventually developed into PCNN (Lindblad and Kinser, Reference Lindblad and Kinser2013; Mohammed et al., Reference Mohammed, Badr and Abdelhalim2015). PCNN is a single-layer artificial neural network model that does not need training on sample data as the network implementation is dominated by an iterative algorithm. In contrast to traditional multi-layer artificial neural networks, PCNN can conduct self-learning and self-supervision.

As shown in Figure 1, the PCNN model includes three mutually coupled functional units: linking component, feeding component and step function. As shown in Equation (17), the input to the linking component is determined by the output from the previous iteration:

(17)$${\bi L}_{ij} \lpar n \rpar =e^{-\alpha_{L} }{\bi L}_{ij} \lpar {n-1} \rpar +V_{L} \times {\bi Y}_{ij} \lpar {n-1}\rpar $$

where, L_ij is the link status related to neuron (i, j), V _L and α _L are an amplification factor and a decay time constant, Y_ij is the output of neuron (i, j) and n is the number of iterations.

Figure 1. The PCNN model.

The input to the feeding component is determined by the output from the previous iteration, which is shown as Equation (18):

(18)$${\bi F}_{ij} \lpar n\rpar =e^{-\alpha_{F} }{\bi F}_{ij} \lpar {n-1} \rpar +V_{F} \times {\bi Y}_{ij} \lpar {n-1} \rpar +A_{ij} $$

where F_ij(n) is the feed status related to neuron (i, j), V _F and α _F are the amplification factor and decay time constant and A _ij is the normalised value of sounding data (i, j). In contrast to the original PCNN model, it considers a single weight value β, which is a constant value. The combined output is shown as Equation (19):

(19)$${\bf U}_{ij} \lpar n\rpar ={\bi F}_{ij} \lpar n\rpar +\lpar {1+\beta {\bi L}_{ij} \lpar n\rpar } \rpar $$

where U_ij is the status of neuron (i, j). As shown in Figure 1, the neuron output is defined by the step function shown in Equation (20):

(20)$${\bi Y}_{ij}\lpar n\rpar \left\{\matrix{1\comma \; & {\bf U}_{ij}\lpar n\rpar \gt {\bi T}_{ij}\lpar n\rpar \cr 0\comma \; & otherwise}\right. $$

where T_ij(n) is the threshold value, which can be dynamically calculated and updated through Equation (21):

(21)$${\bi T}_{ij} \lpar n\rpar =e^{-\alpha_{T} }{\bi T}_{ij} \lpar {n-1} \rpar +V_{T} {\bi Y}_{ij} \lpar n\rpar $$

where V _T and α _T are an amplification factor and a decay time constant, respectively. As noted from Equations (17)–(21), the output of the PCNN is affected by its parameters: V _L, α _L, V _F, α _F, β, V _T and α _T. Therefore, to improve the performance of PCNN, the above parameters must be adjusted to the optimum.

The key features of PCNN are a nonlinear modulation coupling mechanism, the core of the PCNN model and threshold exponential decay mechanism. In the traditional PCNN model, the threshold changes repeatedly. Specifically, after a period of decline, owing to the activation of neurons, the threshold rises sharply, falls sharply and then rises sharply again. Thus, in PCNN, a significant amount of feature information is stored after processing in the activation cycle and activation phase, and the information in the direct output is poor. In order to overcome the above defect and reduce the amount of calculation, we improve the threshold function as follows:

(22)$${\bi F}_{ij} \lpar n\rpar ={\bi A}_{ij} $$

(23)$${\bi L}_{ij} \lpar n\rpar =V_{L} \times {\bi Y}_{ij} \lpar {n-1} \rpar $$

(24)$${\bi T}_{ij} \lpar n\rpar =e^{-\alpha_{T} }{\bi T}_{ij} \lpar {n-1} \rpar $$

The improved model reduces the operating parameters of PCNN, with the meaning of the parameters remaining the same as in the standard PCNN model. The number of parameters is reduced from nine to five, and the complexity of their setting is also reduced. The input F_ij is the value of the sounding data.

When the PCNN model is used to distinguish pseudo-anchor points, let the expressions for U_ij(n) and Y_ij(n) remain unchanged and Equations (22), (23) and (24) respectively replace Equtions (18), (17) and (21). Further, we define the following three concepts:

1. (1) The external input of the PCNN neurons is the water depth value associated with the coordinates, which is F_ij(n) = A_ij;

2. (2) All PCNN neurons have the same structure and parameter;

3. (3) Each PCNN neuron can only be activated once.

3.2. Discriminant method of the pseudo-anchor points

The discriminant process is divided into two steps: (1) matching unit data selection using the PCNN model and (2) matching judgement by the selected data - the position which has the highest probability is the real matching position.

3.2.1. Data selection method based on PCNN

The input to the traditional PCNN is a two-dimensional matrix, whereas the input parameters of terrain matching are at least three two-dimensional matrices. To simplify the description process, suppose that the real-time sounding terrain is also a regular grid terrain. Then, the terrain matching process can be as shown in Figure 2. If the number of matching position and pseudo-anchor points is m, the number of terrain matrices processed by PCNN will be m + 1. The PCNN model can only process one data matrix and the output is Boolean data. Thus, when designing the PCNN structure, the output of the PCNN can be used to determine the similarity of the real-time sounding terrain and the DTM extracted data: ‘0’ means ‘not similar’, ‘1’ means ‘similar’. The inputs to the PCNN are the similarity measure of real-time sounding terrain and the DTM extracted data. Normalising the data to [0, 1], we obtain the similarity measure sequence $\lcub {\lambda_{i} } \rcub _{j=1}^{n} $, where λ _i is defined as follows:

(25)$$\lambda_i=\exp\lsqb -\lpar {\bi y}_{t\comma i}-{\bi h}_{t\comma i}\ \lpar {\bi x}_{t\comma i}\rpar \rpar ^{2}/\sigma_e^2\rsqb $$

Figure 2. Sketch of the terrain searching and matching process.

In the output matrix of PCNN, a ‘1’ indicates that the sounding error of the corresponding data is small, the data should be retained in the second discrimination. Otherwise, if the data has output ‘0’, the data should be removed in the second discrimination. In other words, the ignition process is the similar terrain selection process. Figure 3 is a schematic diagram of the data selection process.

Figure 3. The data selection process based on PCNN.

The PCNN ignition calculation is performed for each node in the matching area, and the number of ignitions is counted. We assume that there are N sounding data points in real-time terrain, and there are m positioning points satisfying the MAP condition. The flowchart of the data selection process is presented in Figure 4.

Figure 4. Flowchart of the data selection process based on PCNN.

3.2.2. The matching judgement method

After the PCNN ignition calculation, the data with the smaller differences are extracted, then further judgement calculations are performed. Before judgement calculations, the data in (y_t − h_t(x_t)) ² have to undergo weight processing using the ignition matrix $\lcub {\mathop{\lambda}\limits^{\frown}}_j\rcub _{j=1}^{N}$, which is defined as follows:

(26)$$W=\sum_{j=1}^{N}{{\mathop{\lambda}\limits^{\frown}}_j}\lpar y_{t\comma j}-h_{t\comma j}\lpar {\bi x}_t\rpar \rpar ^2 $$

Let W replace (y_t − h_t(x_t)) ² in Equation (11), the weighted PDF can be defined as Equation (27):

(27)$$p\lpar {\bi y}_t\vert {\bi x}_t\rpar =\lpar \sqrt{2\pi}\sigma_e\rpar ^{-K}\exp \left[-\displaystyle{{1}\over{2\sigma_e^2}}W\right]$$

where K is the number of nodes with value ‘1’ after PCNN ignition calculation. Calculate the prior PDF of m positioning points and the positioning point which has the maximum value is the best position estimate.

4.1 Search area

For a terrain matching navigation system, the larger the search area, the higher the terrain matching reliability, and the greater the matching time. Therefore, a reasonable search area is necessary.

The search area is related to the performance of the INS. Considering the terrain resolution of DTM, the search area can be defined as a rectangular area determined by a confidence ellipse (Nygren, Reference Nygren2005), with the centre of the rectangular area (p _c,x, p _c,y) indicating the position from the INS. The index numbers of the local terrain can be calculated using Equation (28):

(28)$$R=\left\{{\lpar {x\comma \; y} \rpar \left\vert {\displaystyle{{\vert {x-p_{c\comma x} } \vert }\over{R_{x}}}\le U\wedge \displaystyle{{\vert {y-p_{c\comma y} } \vert}\over{R_{y}}}\le V} \right.} \right\}$$

where R _x is the grid resolution along the x direction, R _y is the grid resolution along the y direction, and, usually, R _x = R _y. U and V are the coefficients associated with l _x and l _y, which depend on the error accumulation characteristics of INS; usually, l _x = 2UR _x, l _y = 2VR _y. The search area is shown in Figure 5, where the w ₁ − −w ₂ coordinate system is the confidence ellipse coordinate system defined by Nygren (Reference Nygren2005).

Figure 5. Confidence ellipsoid and rectangular search area.

4.2. The UTPM process

The flowchart of the UTPM algorithm is shown in Figure 6.

Figure 6. UTPM algorithm flowchart.

The instructions for the calculation process are as follows:

1. (1) Multi-beam data must be processed in real time, which includes removing outlier data and extracting terrain profile information.

2. (2) Multi-beam data include multiple pings. Considering the influence of the amount of computation and the cumulative errors of INS, two to ten pings are usually used.

3. (3) According to the requirements of underwater terrain matching navigation and basic knowledge of probability theory (Xie, Reference Xie2005; Chang et al., Reference Chang, Yang, Kou and Zhang2005), the confidence of the confidence ellipse must not be less than 0·95.

4. (4) The UTPM result can be used to correct the cumulative errors of INS diametrically or fused with other navigational data.

5.1. Underwater semi-physical simulation system

The simulation tests were executed in an underwater semi-physical simulation system. The system centre was a PC104 embedded computer. The structure of the simulation system is shown in Figure 7.

Figure 7. Structure of the underwater semi-physical simulation system.

As shown in Figure 7, the simulation system comprised three components: monitoring computer, environment simulation computer and PC104 embedded computer. The monitoring computer was the actual AUV monitoring computer, which is responsible for issuing task-level parameters and startup instructions. At the same time, the ‘AUV’ status information was displayed. The environment simulation computer simulated the AUV movement in the marine environment and the sensors carried by the AUV. A PC104 embedded computer ran the AUV control system, and it had the same parameters as the master AUV computer, with the only difference being that the sensor information and the data packet were sent by the environment simulation computer, and the actuator information was also sent to the environment simulation computer through the Ethernet. The specifications of the PC104 computer are listed in Table 1.

Table 1. Specifications of the PC104 computer

Based on the simulation system, the simulation test was divided into two parts: Test 1, algorithm validation based on an electronic chart, which verified the feasibility of the proposed algorithm; Test 2, playback simulation test using real multi-beam sounding data, which verified the applicability of the proposed algorithm in a real marine environment.

5.2. Algorithm validation

In Test 1, a visual simulation toolkit issued by MultiGen-Paradigm Inc., Vega software, was used for used to simulate the motion of AUV, Figure 8 shows the initial interface of the Vega-based simulation module. The intersecting rays (The green line in Figure 8) of Vega (were used to simulate a multi-beam sounding system (Li et al., Reference Li, Chen and Dong2011). Adding a Vega-based simulation module to the environment simulation computer, ‘AUV’ motion and data acquisition were conducted in this module.

Figure 8. Initial interface of the Vega-based simulation module.

The DTM for the simulation test was obtained by interpolating the electronic chart. The size of the DTM was 5 km  ×  5 km, and the grid spacing was 10 m. The DTM is shown in Figure 9(a). A simulation model of DTM was then built and input to the Vega-based simulation module.

Figure 9. DTM and matching area selection (a) DTM; (b) matching area selection.

As shown in Figure 9(b), five different terrain feature areas in DTM were selected. During Test 1, the selected local underwater map was interpolated with 2 m grid spacing. The systematic errors of attitude and heading would influence measurement data. However, these impacts are always negligible in terrain matching navigation because a Fibre Optic Gyrocompass (FOG) is very accurate (Hagen and Ånonsen, 2014; Zhao et al., Reference Zhao, Gao, Huang, Wang and Zhou2015). Therefore, the measurement error of the attitude sensor can be ignored in the simulation test. In order to quantify the terrain feature richness, the terrain entropy of the selected local underwater map was calculated (Wang et al., Reference Wang, Yan, Qian and Zhu2007). The interval distance between two adjacent ping samplings was 10 m. Table 2 shows the simulation parameters.

Table 2. Simulation parameters

After setting the simulation parameters, the simulation tests were conducted based on different terrain profile combinations. For comparison, simulation tests using the MLE method (Peng et al., Reference Peng, Zhou, Li and Zhang2016) were conducted with the same simulation conditions. The simulation results obtained after 500 simulation tests are shown in Table 3.

Table 3. Algorithm simulation results based on electronic chart data

In Table 3, the average positioning error is the average error in each matching area, the number of pseudo-anchor points is the maximum value in the simulation, 60 × 2 and 60 × 5 are terrain combined modes. For example, 60 × 2 is a terrain with two terrain profiles comprising 60 units of depth data in each profile.

In the simulation tests, the terrain tended to be gentle from Area 1 to Area 5. According to the definition of terrain entropy, the smaller the terrain entropy value, the richer are the terrain features and the higher the terrain adaptability. Comparing the MAP and MLE positioning results, in the area with smaller terrain entropy, the performance of the two algorithms is similar. As the terrain entropy increased, the matching positioning accuracy decreased and the pseudo-anchor points appeared. As can be seen in Table 3, under the same test conditions, the number of pseudo-anchor points of MAP is less than that of MLE. Without removal of the pseudo-anchor points, the matching accuracy of MAP is higher than that of MLE. Figure 10 shows the terrain matching results of MAP and MLE in Area 3. As can be seen, under the same simulation conditions, the pseudo-anchor points of MAP are less than those of MLE, and the probability distribution is more reasonable; hence, MAP is superior to MLE.

Figure 10. Terrain matching results (a) MAP (b) MLE.

Comparison of the matching positioning results of different terrain profile combinations shows that the matching positioning accuracy of 60 × 5 is higher than that of 60 × 2. As can be seen in Table 3, in Areas 3 and 4, the number of pseudo-anchor points is significantly reduced. For example, there are five pseudo-anchor points using the MAP method in Area 4 by the 60 × 2 terrain, and only one pseudo-anchor point by the 60 × 5 terrain. However, in areas where the terrain is flatter (Area 5), the number of pseudo-anchor points is not significantly reduced (from 11 to ten using MAP and from 30 to 28 using MLE), and the calculation time is significantly increased. As can be seen, in the flat terrain area, because there are more pseudo-anchor points, the real-time performance of MLE is worse.

Figure 11 shows the posterior PDF distribution with different real-time terrain combinations in Area 5. As can be seen, in the flat terrain area, the posterior PDF distribution of 60 × 5 is not changed significantly compared with 60 × 2. In other words, the matching performance cannot be significantly improved by increasing real-time terrain data. After the discrimination of pseudo-anchor points by PCNN, the true positioning point can be found, and the positioning accuracy increases significantly. Therefore, the proposed pseudo-anchor points discriminant method is effective. However, it should also be noted that the calculation time of PCNN is proportional to the number of pseudo-anchor points and the amount of real-time terrain data. As the number of pseudo-anchor points and the amount of real-time terrain data increases, the calculation time also increases, which can affect the real-time performance of terrain matching.

Figure 11. Probability distribution with different real-time terrains in Area 5 (a) 60 × 2 (b) 60 × 5.

In theory, by screening and determining each search point using PCNN, then calculating the PDF, the true matching position can be obtained. Figure 12 shows the probability distribution after PCNN screening with 60 × 5 real-time terrain data. After screening all search points, the probability distribution is a single peak distribution, which means that the effect of ‘error averaging’ is effectively suppressed. However, in the simulation, the calculation time is over 100 seconds, which is unacceptable for real-time performance. Therefore, combining MAP and PCNN is effective for terrain matching positioning. With advancements in computer performance, PCNN terrain matching positioning has greater advantage.

Figure 12. Matching positioning based solely on PCNN.

5.3. Playback simulation

In Test 2, the data source of DTM was multi-beam sounding data measured in a sea trial. The multi-beam sounding data were acquired using GeoSwath Plus (GS+), a phase interferometry multi-beam sounding side scanner produced by the GeoAcoustics Company. The size of the measured area was 1,000 m  ×  900 m, and the depth was from 5 m to 40 m. After filtering and gridding processing, the grid spacing of DTM was 1 m  ×  1 m (Figure 13).

Figure 13. The DTM and independent survey line.

The navigational data in the GS+ original file were obtained using the Real-time Kinematic (RTK) technique (GeoAcoustics Limited, 2007). Thus, the navigation accuracy was up to several centimetres (Yao et al., Reference Yao, Hu and Xu2016). Therefore, RTK data can be used for the real position.

The path indicated by the arrow in Figure 13 is an independent multi-beam survey line, which is used to simulate real-time sounding data. The survey line is perpendicular to the lines marking the DTM; thus, data independence is ensured.

The playback simulation means that the data in the simulation are from sea trials and the data are read according to the state when the data were acquired in simulation tests. As the test data are from real multi-beam sounding, the playback simulation could verify the applicability of the proposed method in the real marine environment. The initial error of INS is established as σ _x = σ _y = 15 m, and the error characteristic satisfies σ _x = σ _y = N(0, (0 · 05t) ²). Calculation was completed of the terrain features of the independent line passing through the area, and five matching areas can be selected sequentially along the survey line, as shown in Figure 12. Similarly, the terrain entropy in the matching area is calculated separately and is shown in Table 4.

Table 4. Terrain entropy in the matching area.

Multiple playback simulation tests were conducted, with the results in Table 5 being obtained. Matching length is defined as the distance between the two furthest real-time terrain profiles. For example, there are five terrain profiles in the real-time terrain data, the length between adjacent terrain profiles is 10 m and the matching length is 40 m. The accuracy requirement is that positioning errors be less than 5 m.

Table 5. Playback simulation results

As can be seen from Table 5, in Area 1, which has the richest terrain features, the true positioning point can be found using only 80 × 2 real-time terrain data and increasing the data does not improve the positioning accuracy. In Areas 3 and 4, which have relatively rich terrain features, there are pseudo-anchor points using 80 × 2 real-time terrain data, but increasing the data to 80 × 5 causes the pseudo-anchor points to disappear; the positioning accuracy of MAP using 80 × 5 is close to the PCNN result for 80 × 2. In Areas 2 and 5, which have the poorest terrain features, the pseudo-anchor points still exist with 80 × 5 real-time terrain data; pseudo-anchor points discrimination must be conducted. In Area 5, the positioning accuracy between 80 × 5 and 80 × 2 is approximately the same, whereas in Area 2, the positioning accuracy meets the requirements only with 80 × 5 real-time terrain data. This is because the terrain is very flat in Area 2, which results in the terrain data using 80 × 2 being too small to describe the terrain features after PCNN ignition calculation. Therefore, before terrain matching positioning, the matching suitability of the local terrain must first be assessed, which is used to determine the real-time terrain combination mode in the terrain matching calculation.

As can be seen from the above analysis, the true positioning point can be found using only a small amount of real-time terrain data by the proposed terrain matching positioning method in most areas, and the pseudo-anchor points can be removed effectively. However, in excessively flat terrain areas, an overly small amount of terrain data cannot describe the terrain features, and so the amount of real-time terrain data has to be increased. In the playback simulation tests, because of the high resolution of DTM, the positioning error was less than 5 m in each matching area, and the accuracy was unrelated to the sailing distance and sailing time. Therefore, the proposed method is applicable for AUV underwater navigation.