1. INTRODUCTION
An Autonomous Underwater Vehicle (AUV) is an indispensable piece of engineering equipment in ocean exploration. Because of the complexity of the underwater environment, the research and development of AUVs has encountered many difficulties, and precise navigation is one of them. At present, underwater navigation is mainly based on the Inertial Navigation System (INS), but the cumulative error of INS is an unsolvable problem. Regularly breaking the surface to “zero out” the position error greatly increases energy consumption, and for long-endurance AUVs, this is particularly undesirable. Terrain positioning is therefore an important problem. To solve this problem, researchers have proposed terrain-aided navigation methods that were first used in the field of aviation. Initial studies focused on improvements based on land terrain matching navigation systems. There are already a number of countries and institutes that have developed their own underwater terrain matching navigation system.
For instance, the Norwegian Defense Research Establishment developed the TerrLab simulation system (Hagen, Reference Hagen2006) and have conducted several sea trials based on the HUGIN AUV. In 2010, they conducted a sea trial for five hours in Oslo, and the navigation error was only 5 m (Anonsen and Hagen, Reference Anonsen and Hagen2010; Hagen et al., Reference Hagen, Ånonsen and Mandt2010). Stanford University and the Monterey Bay Aquarium Research Institute developed a low-cost terrain navigation method that is suitable for long-range AUVs. They took an MBARI Dorado AUV out into Monterey Bay in 2008 to verify the feasibility of using low-precision measuring equipment for terrain-aided navigation and used low-accuracy stern projections and a terrain navigation filter, achieving an accuracy of 4–10 m (Meduna et al., Reference Meduna, Rock and McEwen2008; Reference Meduna, Rock and McEwen2010). The Swedish Royal Institute of Technology (Nygren, Reference Nygren2008), University of Southampton (Carreno et al., Reference Carreno, Wilson, Ridao and Petillot2010), and University of Tokyo with the Japan Institute of Marine Engineering (Nakatani et al., Reference Nakatani, Ura, Sakamaki and Kojima2009) have also developed their own underwater terrain matching navigation systems. In addition, Peking University (Zhang et al., Reference Zhang, Si, Yan and Ge2003), Wuhan University (Zhang et al., Reference Zhang, Zhao and Shao2011; Reference Zhang, Zhao and Chen2014), and Harbin Engineering University (Chen, Reference Chen2013; Tian, Reference Tian2007) have also studied underwater terrain matching navigation systems. However, most studies are based on improving the aviation terrain matching algorithm or single beam terrain matching system. Related studies include the improved TERCOM algorithm (Liu, Reference Liu2003) and Image Matching Technology (Zhang et al., Reference Zhang, Zhao and Shao2011).
In contrast to land-based environments, the underwater environment makes it difficult for current bathymetric technology to obtain a terrain map that is as precise as a land-based map. In addition, the underwater terrain is flatter and AUVs move slowly, which decrease the signal-to-noise ratio and reduce positioning accuracy, reliability, and stability. This problem is particularly prominent in flat regions (Chen, Reference Chen2013; Xie, Reference Xie2005). In this paper, the Node Multi-Information Fusion (NMIF) algorithm is proposed, which effectively improves the accuracy of flat terrain positioning and measurement stability. This method was verified using a real underwater terrain map.
2. CONSTRUCTION OF NODE DATA PACKETS
In terrain matching, positioning methods are based on the Digital Elevation Map (DEM), which stores the height of each grid node sampled from the ocean floor. In the Measured Terrain Map (MTM), all heights form a sequence of height values. Matching algorithms determine the position of a point in the DEM by matching an MTM height sequence with the height sequence from the DEM. This point is called a positioning point, but it is never the actual position of the AUV. In a flat region, false position peaks caused by low numbers of measured terrain nodes and limited feature information are inevitable (Chen, Reference Chen2013; Xie, Reference Xie2005; Nygren, Reference Nygren2005), so the focus of this study is to extract more information from limited topographic terrain features. In fact, multi-beam methods increase the number of nodes measured per measurement period, and this increases the amount of position information, which does improve the accuracy and efficiency of terrain matching. However, there is no improvement in the utilisation of node information. There is a difference between a terrain profile formed from a single beam and a topographic profile formed from multiple beams. The surface includes significantly more topographic features than a terrain profile does, and it does not significantly improve accuracy and efficiency to only use height information. A single node includes only height information, but if it is observed using a multi-beam MTM, it is possible to add surrounding surface information into a node data packet, causing it to become more distinct. The construction of such a data packet, the Single Node Data Packet (SNDP), is shown in Figure 1.
In order to make full use of the node information, this paper proposes a new terrain matching method based on the SNDPs. Suppose a ij denotes the data packet of node (i, j) in the DEM, expressed as an ξ-dimensional vector, a ij = [σ 1, σ 2 · · · σ ξ ].
In Figure 2 an SNDP contains the information surrounding one node and consists of four information units. How much information is included in these units is decided by many factors, but the reliability of the packet information must be ensured. In contrast to a single node, which only contains height information, more information is included in the SNDP. Hence, one can think of the SNDP as an expansion of the node information, or a fusion of multiple types of information regarding a certain single node.
3. TERRAIN MATCHING BASED ON PACKETS
Similarity is the key to positioning, and a packet contains a variety of information with which to calculate it. We assume that every item of information is independent of each other. A sub-element matching area is a part of the search area, and its location and size is provided by inertial or reckoning navigation and the corresponding navigation error. The search area, sub-element matching area, node, and surrounding node information are shown in Figure 3.
3.1. Matching SNDPs
The MTM contains SNDPs to be matched with the corresponding nodes in the DEM, and the matching algorithm decides which SNDP in the DEM is most similar to the SNDP from the MTM. For example, as shown in Figure 4, suppose s ij is a node in the DEM and there are k × l corresponding nodes from the MTM that must be compared with it. Further suppose that the matching algorithm considers node a ij from sub-element areas A 11 and A ij to be similar to s ij . In this case, A 11 and A ij have SDNP s ij added to them. Figure 4 shows the process that determines the SNDPs to be matched in the DEM and MTM and the corresponding comparison of the SNDPs. If an SNDP in the DEM qualifies as similar, as defined in the algorithm, then the node in the DEM has that SNDP added to it. The centre of the matched sub-element search area in the DEM records the number of added SNDPs, and a higher number indicates greater similarity. Of course, if a node in a certain sub-element in the matched area meets the similarity condition, it can always add node s ij . In Figure 4, s 11 is similar to SNDP α 11 in sub-element areas A 11, A 12, and A ij , so A 11, A 12, and A ij will have one SNDP added each. As shown in Figure 4, A 11 adds one SNDP, A 12 adds one SNDP, A ij adds three SNDPs, and A kl adds no SNDPs. In this case, sub-element area A 11 is the most similar to the MTM, and the searching point represented by A ij is marked by a red dot.
3.2. Algorithm implementation
The number of SNDPs added at a certain sub-element area in the DEM is a measure of the level of similarity between the MTM and the sub-element search area. In Figure 5, it can be seen that only the similarity between the SNDPs from the sub-element search area and the SNDPs from the MTM exceed the threshold defined in the algorithm. The node in the sub-element search area adds the packet, and the value of A ij is incremented. At this point, A ij represents both the added data packets and the centre of the sub-element search area. Details of the procedure are as follows.
3.2.1. Node similarity measurement -Algorithm
(1) For example, if there is an m × n node DEM (m and n represent the number of nodes in the x and y directions, respectively), the feature vector of each node must be calculated, so the algorithm must calculate m × n feature vectors (σ 1, σ 2 … σ ξ ), where ξ represents the number of features of each node. The algorithm marks the centre of the sub-element search area with A ij (i = 1, 2, …m, j = 1, 2, …n). Assuming that every sub-element search area contains k × l nodes and each feature vector has ξ dimensions, the algorithm marks the sub-element search area with a ij i = 1, 2, …k, j = 1, 2 …l, and a ij = (σ 1, σ 2, …σ ξ ). Different nodes will have different values for σ 1 σ 2 σ 3…σ ξ . The relation between A ij , a ij , and (σ 1, σ 2 …σ ξ ) is shown below.
The MTM is processed in the same way to obtain its feature vector. It is marked with (σ 1, σ 2, σ 3 · · · σ ξ ) because a MTM has only one unit S. Hence, the relation between (σ 1, σ 2, σ 3 · · · σ ξ ) and S is as shown below.
(2) Using a similarity decision function, the algorithm calculates the similarity between two feature vectors. For example, given a node s(i, j) from MTM S and node a(i, j) from sub-element search area A ij i = 1, 2, …m, j = 1, 2, …n, which is the node corresponding to s(i, j), the algorithm sets λ ij i = 1, 2, …m, j = 1, 2, …n, to be the similarity measurement between a(i, j) and s(i, j), where λ threshold represents the similarity threshold. That is, if the similarity λ ij > λ threshold , s(i, j) and a(i, j) are considered to represent the same node. At the same time, sub-element search area A ij adds an SNDP and its value is incremented by one.
(3) Step 2 is repeated until all nodes have been compared with the corresponding node. Finally, the algorithm calculates the sum of all the SNDPs added to sub-element search area A ij as follows.
The value of δ ij is defined by the following formula:
The algorithm next determines the maximum A max from all A ij (i = 1, 2, …m, j = 1, 2, …n). Here, A ij represents the centre of the sub-element search area and simultaneously the sum of the SNDPs that sub-element search area A ij has added. Hence, A ij contains two elements, one is a coordinate that represents the centre of A ij and another is a number that represents the sum of the SNDPs that A ij has added. As mentioned, higher values of A ij indicate a higher similarity between A ij and S. Hence, A max represents the sub-element search area that is the most similar to the real-time DEM S, and the coordinate of A max is the final position.
3.2.2. Similarity measurement function
The similarity function is a measurement function that represents similarity (Zhang et al., Reference Zhang, Liu and Ji2009; Shi, Reference Shi2003). We define a function that transforms the deviation between two vectors into a measurement scale λ. We use the Gaussian function to measure similarity. Suppose Θ is an n-dimensional vector, Θ i i = 1, 2, …m is a sample space that contains m samples, and Θ i also has n dimensions. We formulate Θ as follows: $\Theta = \left( {\matrix{ {\sigma _1} \cr {\sigma _2} \cr \vdots \cr {\sigma _n} \cr}} \right)$ , $\Theta _i = \left( {\matrix{ {\sigma _{i1}} \cr {\sigma _{i2}} \cr \vdots \cr {\sigma _{in}} \cr}} \right)$ , n = 1, 2, … . Suppose Δ i is the deviation between Θ and Θ i . Using the Euclidean distance to measure the deviation between Θ and Θ i , Δ i can be represented by the following equation.
Our goal is now to determine the Θ i that are similar to Θ using the value Δ i . The problem is how to represent the level of similarity and decide whether Θ i is similar. The solution is to use a function f(Δ i ) to transform Δ i into a value λ i , where λ i = f(Δ i ) is the data range [0, 1], then define a threshold λ threshold . If λ ij > λ threshold , then Θ i is considered similar to Θ. From this analysis, λ should have the following properties:
-
1. If Δ i = 0 (i = 1, 2, 3 …), i.e., the deviation between Θ and Θ i is zero, Θ and Θ i are exactly similar, then λ i = 1. If Δ i → ∞, Θ and Θ i are exactly dissimilar, at this point λ i = 0.
-
2. As the similarity value increases, λ i increases monotonically.
-
3. The value of λ i is within the range [0, 1].
Given the above analysis, we constructed a similarity measurement function based on the Gaussian function (Zhang and Yang, Reference Zhang and Yang2004). The properties of a Gaussian function mean that it can transform deviation value Δ i into a number between [0, 1], achieving our goals. The normal form of a Gaussian function is:
Suppose x is parameter Δ i . At point b = 0, λ = 1 must be true, so a = 1. After reformulating the new Gaussian function by replacing f(Δ i ) with λ, we have:
Parameter c is an important element that controls the shape of function λ. Hence it clearly influences the value of λ. For various c, the shape of function λ becomes slimmer as the value increases. Hence, λ is high in a small area surrounding “0”. The situation will be the reverse if parameter c is small and we decrease its value. Figure 6 shows the shape of function λ for different c, indicating the value of λ at point Δ = 1.
This result shows that when c is small, the ability to measure the similarity between two vectors is significantly weakened. Hence, in practice, we should choose larger values. We discuss how to choose value c below. First, we establish the principle of similarity measurement, as shown in Figure 7.
Suppose X and Y are two measurement results at the same point. The measurement noise of the two measurements is denoted by ω x and ω y , respectively. Because measurements X and Y come from the same point, given our definition of similarity, we examine the value of λ, which must exceed threshold λ threshold . Suppose λ threshold = 0.9. The similarity can be calculated as follows:
If X and Y come from the same node without measurement error, X − Y = 0, and clearly λ = 1. However, we must take the error into account (Chen, Reference Chen2013). The similarity function is reformulated as follows:
Usually, we assume that the measurement noise is zero-mean Gaussian white noise. As it is preferable to use a larger value, |ω x − ω y | should result in a larger value, that is, σ x + σ y . Here, σ x + σ y represents the sum of the standard deviations of the measurement noise for both measurements. We further suppose that each measurement noise is independent. If X and Y are vectors with n dimensions, then σ x + σ y will be a diagonal matrix with n dimensions (Xie, Reference Xie2005). The above explanation describes the value of parameter c when λ threshold = 0.9. If the threshold is changed, the procedure is unchanged. In order to improve stability and efficiency, choosing a higher threshold is better. Hence, we always use λ threshold ≥ 0.9.
4. EXPERIMENTAL ANALYSIS
Using theoretical analysis, we obtained the structure and necessary parameters of the algorithm. We next examine the stability and robustness of the algorithm by conducting a series of experiments using an actual ocean floor DEM located in the South China Sea (Figure 8). The data was obtained using a GS+ bathymetric system and bilinear interpolation to create a DEM with a grid size of 4 × 4 m.
4.1. Best feature selection
The DEM has an area of 800 × 700 m, with a grid size of 4 × 4 m. The grid size in the MTM is 2 × 2 m; hence, we used a linear interpolation algorithm to transform the DEM into a 2 × 2 m grid. A simulation MTM was obtained from the DEM and covers an area of 200 × 200 m. The centre was located at (500, 500) in the local coordinate system, the same as below. The position of nodes in the global coordinate system can be determined by the following formula:
where x 0, y 0 represent the coordinates of the origin of the local coordinates in the global coordinate system, d is the length of a grid side, and i, j represents the index of the nodes.
At the same time, we take the measurement of the circular error probability and height error into consideration. We set the circular error probability to be σ r = 0.5 m and the height error probability to σ = 1 m. They are both modelled as zero mean Gaussian white noise. Considering the deviation in MNS, we chose an area of 200 × 200 m centred at point (500, 500). The search area and MTM are shown in Figures 9 and 10.
4.1.1. Node packets containing height only
The similarity measurement function is expressed as follows:
In this equation, σ 1(s) represents the height of a node in the DEM, and σ 1(a) represents the height of a node in the MTM. We set the threshold to λ threshold = 0.9 and interpolated the DEM with grids respectively equal to 2, 3, 4, 5, and 6 m into a 2 m grid. The final simulation results are shown in Figure 11 for the 2, 4, and 6 m grids. Figures 11(a)–(c) show the distribution of the sum of the SNDPs in the search area using height only. Figure 12 shows the number of SNDPs at the positioning point in the three tests. From the results, it can be seen that as the grid length increases, the number of added SNDPs surrounding the located point increases rapidly. The false peaks appear when the grid length is set to 6 m.
4.1.2. Node packets containing surface features only
Consider the case where there is more than one node in the DEM that is similar to the node in the MTM. In such a case, we must judge whether this node comes from the positioning area. The algorithm cannot decide by height alone because it is too similar to other positions; hence, the surface information becomes useful for distinguishing many matching nodes in a series of similar nodes. The surface information can be represented in many ways, such as by entropy information, slope, gradient, and so on. Briefly, we must choose a method that can express the features of the terrain surface effectively and accurately. In this paper, we use the sum of the height of a node and its height standard deviation in the local area that is centred on that node. The area should be extensive enough so that it can distinguish the node from all other nodes at a similar height. The feature equation is as follows:
where σ 2 represents the feature parameter of a node that contains height and surface information, and σ 0 represents the standard deviation of height in the local area that is centred on the matching node. In this paper, the area of the local region is 10 × 10 m2. Furthermore, σ 1 represents the height of the matching node. The similarity is measured as follows:
In Equation (16), σ 2(s) represents the blended surface feature of the node in the MTM and σ 2(a) represents the blended surface feature surrounding the node in the DEM. We set the threshold to λ threshold = 0.9. We also chose the DEM with grids respectively equal to 2, 3, 4, 5, and 6 m interpolated into a 2 m grid. The final simulation results are shown in Figures 13(a)–(c). These figures show the distribution of the sum of SNDPs in the sub-element search area centre for the 2, 4, and 6 m grids using the surface feature only. Figure 14 shows the number of SNDPs at this position point in the three tests.
The test results show that the surface feature in the local area can identify similar nodes more easily, as the final point had more SNDPs added. On the other hand, the ability to inhibit false peaks was weakened.
4.2. Node Multi-Information Fusion (NMIF)
From these results, it can be seen that both height and surface information have their own advantages. Both types of information converge near matches, and using surface features makes it easier to identify similar nodes but does not inhibit false point matches. Hence, we next fuse the two types of position information in two ways and compare them.
4.2.1. Fusion based on logic
For logical “and” fusion, the result is true only if all of the parameters meet the criteria. In this method, both the height and the surface height characteristics must exceed the threshold before the node in the MTM receives the SNDP.
In this equation, $\lambda _{ij}^{\sigma _1} $ and $\lambda _{ij}^{\sigma _2} $ are at a similar scale for the height and surface features.
We next test the algorithm on the DEM 6 m grid, comparing it with the results using height only. The results are shown in Figure 15.
Figure 15(a) shows the distribution of the sum of SNDPs using logic fusion, Figure 15(b) is the result using height value only, and Figure 15(c) is the result using surface features only. This comparison shows that for the same grid precision and measurement noise, the results based on the fusion of logic “and” inhibits false peaks around the located point and the number of added packets is more focused. The advantage of NMIF is clear.
4.2.2. Fusion based on vector similarity
Vector similarity fusion, that is, thinking of the SNDP as a vector, compares two node similarities as two vectors (Sun, Reference Sun2002; Zhang et al., Reference Zhang, Liu and Ji2009). Measuring the similarity of two vectors is usually equal to comparing their angle or length. First, we transform the two characteristic parameters of the node σ 1 and σ 2 into a two dimensional vector Γ = (σ 1, σ 2). We use the Euclidean distance as the similarity measurement of two vectors, expressed as follows.
where Γ s and Γ a are two-dimensional vectors composed of the node characteristic parameters, and c is a two-dimensional diagonal matrix representing the measurement noise of σ 1 and σ 2, the two elements of Γ s and Γ a . Therefore, the similarity is measured as follows:
The experimental results for the DEM on a 6 m grid are shown in Figure 16. Figure 16(a) shows the distribution of the sum of SNDPs using vector similarity fusion, Figure 16(b) is the result using height value only, and Figure 16(c) is the result using surface features only.
From the results, we can see that for the same grid precision and measurement noise, the vector similarity method inhibited the false peaks around the located point better than the result using height value or surface feature information alone. The sum of the added SNDPs at the position point is more focused.
5. ACCURACY ANALYSIS
As can be seen from the above results, NMIF can make full use of the information of a single node and its surrounding surface. The nearby pseudo-peaks surrounding the positioning node have been inhibited effectively. These results show that for the same terrain, NMIF obtains a more concentrated distribution of the number of packets, which is advantageous for positioning point identification, and inhibits pseudo-peaks. We now focus on analysing the precision and stability before and after using data fusion for positioning. In order to compare the precision and stability of the proposed method and the three comparison methods (no fusion, fusion based on vector similarity, and fusion based on logic) in a DEM with low precision and high measurement noise, we chose a DEM with a grid of 6 m and interpolated it into a 2 m grid. In addition, we added Gaussian measurement noise with zero mean, standard variance σ = 1, and circular error probability σ r = 0.5. We computed the deviation for every experiment using $\Delta d = 2 \times \sqrt {\Delta i^2 + \Delta j^2} $ . The results of the simulation output are as follows.
Figures 17(a)–(c) are the results for “height value,” “logic fusion,” and “vector similarity” for a grid of 6 m interpolated into a 2 m grid. Figure 18 compares the positioning error of the three positioning methods, and Figure 19 shows the mean and standard deviation of the positioning error for the three positioning method tests.
From the simulation positioning results, we can see that the positioning point distribution based on the NMIF method, both for vector similarity fusion and logic fusion, are more concentrated around the true location, and the mean and standard deviation of the positioning deviation are far smaller than the positioning result using height value only. It is clear that the NMIF method has better accuracy and stability than the height-based matching position method. NMIF extracts more than one terrain feature and the additional positioning information from the terrain surface increases the accuracy of the matching position stability.
5. CONCLUSIONS
The essence of terrain matching positioning is a process that positions terrain features in the DEM. Previous methods consider the measured data as a discrete point sequence and hence cannot take full advantage of the measured information. As a result, they obtain a poor position result because of the lower amount of measured data or terrain features. NMIF is a good way to extract the terrain features, and the test results have shown the great advantages of NMIF at terrain matching positioning. We hence obtain the following conclusions. First, the MTM and DEM are not only point-wise discrete sequences, but form three-dimensional surfaces that include a large amount of terrain feature information. There is a determinate position relationship between the nodes, so we can extract more positioning information. Second, the experiments show that NMIF can extract diverse characteristics of the topography, increasing the utilisation ratio of the measured terrain to improve the accuracy and stability of the terrain matching position. Finally, as the measurement noise and topography grid side lengths increase, the position deviation will increase, but NMIF performs better than the position result using height or surface features alone. That is to say, the NMIF is more robust to measurement error.
ACKNOWLEDGEMENT
This project was supported by the National Natural Science Foundation of China (51279221).