2.1. Formation of virtual LBL matrix window

The proposed method in this study is based on a virtual LBL matrix window. In this system, a single sound source is placed at the bottom of the sea and sends periodic signals and a single hydrophone is installed on the AUV. In the AUV navigation trajectory, four selected recent positions of the AUV are regarded as four virtual hydrophones of the LBL matrix, which constitute a virtual LBL matrix window, as shown in Figure 1.

Figure 1. Formation of virtual LBL matrix.

2.2. Positioning algorithm based on a moving virtual LBL matrix window

In this method, the single sound source is fixed on the seabed in which its absolute position is measured as (x, y, z) in advance in Earth-centred Earth-fixed coordinates. This sound source sends periodic pulse signals at regular intervals, with the time period defined as Ts. A virtual LBL matrix window is created when the AUV enters the area where the source signal can be received. The current position of the AUV is set as the initial position P ₀, where the period of receiving acoustic signals is denoted as x₀. Then, the AUV continues to move forward, and its position is denoted as P ₁ at the beginning of the eleventh period. A period of receiving acoustic signals is denoted as x₁. In this way, the AUV position is denoted as P _i(i = 2, 3,···), and the receiving acoustic signal is denoted as x_i(i = 2, 3,···) for every ten periods. Each of the four consecutive positions of the AUV forms a virtual LBL matrix window. For example, the first virtual LBL matrix window is composed of the initial four positions (P ₀, P ₁, P ₂, P ₃), which are marked by the blue rectangle in Figure 2. The second virtual LBL window is formed by removing position P ₀ and adding position P ₄ when the AUV arrives at position P ₄; this window is shown as the red rectangle in Figure 2. In this way, the moving virtual LBL windows are established by removing the previous position and adding the more recent position to iteratively calculate the current position of the AUV.

Figure 2. Time window model based on periodic movement.

3.1. Iterative updating of position based on moving virtual LBL windows

Considering that the same algorithm is used to calculate the position of the AUV based on each window, we describe the calculation method for the AUV by taking the first window as an example.

The initial position of the AUV is denoted as P ₀(x ₀, y ₀, z ₀) in the inertial coordinate system. From the above description, each new position of the AUV is obtained every 10Ts, for example, P ₁(x ₁, y ₁, z ₁), P ₂(x ₂, y ₂, z ₂) and P ₃(x ₃, y ₃, z ₃). The four positions form the first virtual LBL matrix window. The distances between position P _i(i = 0, 1, 2) and P ₃ in the three-axis direction (x, y, z) are defined as Δ x _3i, Δ y _3i and Δz _3i respectively; their values can be measured by dead reckoning using INS. P _i can be described by recent position P ₃ as follows; its position is considered unknown and should be calculated:

(1)$$\left\{\matrix{\lpar x_{0}{,y}_{0}{,z}_{0} \rpar = \lpar \lpar x_{3}{-\Delta x}_{30} \rpar ,\lpar y_{3}{-\Delta y}_{30} \rpar ,\lpar z_{3}{-\Delta z}_{30} \rpar \rpar \hfill \cr \lpar x_{1}{,y}_{1}{,z}_{1} \rpar =\lpar \lpar x_{3}{-\Delta x}_{31} \rpar , \lpar y_{3}{-\Delta y}_{31} \rpar , \lpar z_{3}{-\Delta z}_{31} \rpar \rpar \hfill \cr \lpar x_{2}{,y}_{2}{,z}_{2} \rpar = \lpar \lpar x_{3}{-\Delta x}_{32} \rpar , \lpar y_{3}{-\Delta y}_{32} \rpar , \lpar z_{3}{-\Delta z}_{32} \rpar \rpar \hfill}\right.$$

Time difference τ_3i of the received signal at positions P₃ and P _i can be obtained through a cross-correlation operation. Two factors lead to the time difference. The first factor is the difference in geographic location (Δ t _3i), and the second factor is the time interval of the sent signal, as indicated in Equation (2):

(2)$${\Delta t}_{3i} = \tau_{3i} - \lpar 30-10i \rpar t, \quad \lpar {\hbox{i} = 0,1,2} \rpar $$

The sound propagation velocity underwater is assumed to be a constant value denoted as c, and three equations are established as follows:

(3)$$\left\{\matrix{{\Delta t}_{30}c = \sqrt {{(x_{3}-x)}^{2}+{(y_{3}-y)}^{2}+{(z_{3}-z)}^{2}}\hfill \cr -\,\sqrt {{(x_{0}-x)}^{2}+{(y_{0}-y)}^{2}+{(z_{0}-z)}^{2}} \hfill \cr {\Delta t}_{31}c = \sqrt {{(x_{3}-x)}^{2}+{(y_{3}-y)}^{2}+{(z_{3}-z)}^{2}} \hfill \cr -\,\sqrt {{(x_{1}-x)}^{2}+{(y_{1}-y)}^{2}+{(z_{1}-z)}^{2}} \hfill \cr {\Delta t}_{32}c=\sqrt {{(x_{3}-x)}^{2}+{(y_{3}-y)}^{2}+{(z_{3}-z)}^{2}} \hfill \cr -\,\sqrt {{(x_{2}-x)}^{2}+{(y_{2}-y)}^{2}+{(z_{2}-z)}^{2}}\hfill}\right.$$

(x ₀, y ₀, z ₀),(x ₁, y ₁, z ₁) and (x ₂, y ₂, z ₂) can be expressed as expressions that contain (x ₃, y ₃, z ₃) through Equation (1), and recent position P₃(x ₃, y ₃, z ₃) can be calculated by Equation (3). Subsequently, the recent geodetic coordinates can be obtained through the coordinate transformation of latitude and longitude.

The AUV coordinates in the geodetic coordinate and Earth rectangular coordinate systems are set as (λ, L, H) and (x, y, z), respectively. The formula for converting the geodetic coordinates into the Earth frame's rectangular coordinates is expressed as follows:

$$\left\{\matrix{x = \lpar N + H \rpar cosLcos\lambda \hfill \cr y = (N + H) cosLsin \lambda \hfill \cr z = [N\lpar \hbox{1-}e^{2} \rpar +H] sinL} \right.$$

where e is the first eccentricity and N is the radius of the unitary circle at latitude L.

$$N = a / \sqrt{\lpar \hbox{1-}e^{{\rm 2}}{sin}^{{\rm 2}}L \rpar }$$

where a is the long semi-axis of the Earth.

Therefore, the formula for converting the Earth's rectangular coordinates to geodetic coordinates is:

$$\left\{\matrix{\lambda = arctan (y / x) \hfill \cr L = \hbox{arctan} \left(\displaystyle{z+be_{{2}}^{{2}}{sin}^{{\rm 3}}U \over \sqrt{x^{{\rm 2}}+y^{{\rm 2}}} -ae^{{2}}{cos}^{{\rm 3}}U}\right)\hfill \cr H = \sqrt {x^{{\rm 2}}+y^{{\rm 2}}} / {cosL-N} \hfill} \right.$$

Among the coordinates, b is the short semi-axis of the Earth and e ₂ is the second eccentricity, where $U=arctan\left(z / {\sqrt{x^{{\rm 2}}+y^{{\rm 2}}} \sqrt {\hbox{1-}e^{{\rm 2}}}} \right)$.

The first and second virtual LBL windows are used as examples to explain the iterative updating of the AUV position. The AUV continues to move and obtains the current position P ₄(x ₄, y ₄, z ₄) after 10Ts when the recent position P₃(x₃, y ₃, z ₃) is obtained. P ₀(x ₀, y ₀, z ₀) is removed from the window and P₄ is added in the new window with positions P₁, P₂ and P₃. The AUV updates the current position by using this method, which corrects the errors of the SINS.

3.2. Calculation of Time Difference of Arrival (TDOA)

A TDOA positioning method provides the sound source position by measuring the TDOA at different targets from the sound source. Given the difficulty of synchronising underwater signals, the time difference is easy to obtain through the cross-correlation processing of the received signals.

Assuming that the received signals at positions P_i and P_j are:

(4)$$\matrix{\hbox{x}_{i}=\hbox{a}_{i}x\lpar t-{\rm \tau }_{i} \rpar +\hbox{n}_{i}{\rm (}t{\rm )} \cr \hbox{x}_{j}=\hbox{a}_{j}x\lpar t-{\rm \tau }_{j} \rpar +\hbox{n}_{j}{\rm (}t{\rm )}\hfill}$$

where a_i and a_j are the attenuation coefficients of the acoustic signals that propagate underwater, n_i(t) and n_j(t) are the uncorrelated noise signals and τ_i, τ_j are the propagation times.

The cross-correlation function of x _i(t) and x_j(t) is:

(5)$$\hbox{R}_{{\rm x}_{i}{\rm x}_{j}} ({\rm \tau}) {\rm =E} [\hbox{x}_{i} (\hbox{t}) \hbox{x}_{j}^{\ast} (\hbox{t}-\tau)] = \displaystyle{1 \over \hbox{T}-\tau} \int_{\rm \tau}^{\rm T} \hbox{x}_{i} (\hbox{t}) \hbox{x}_{j} (\hbox{t}-\tau) \hbox{dt}$$

where τ = τ _j − τ _i indicates the TDOA signals and T denotes the observation time. In theory, the peak of $\hbox{R}_{{\rm x}_{i}{\rm x}_{j}}$ corresponds to the time difference of x _i(t) and x_j(t), which is τ.

3.3. Smooth Coherent Transformation (SCOT)-weighted generalised cross-correlation

As a result of complicated underwater environments, refractions and reflections are common, and they are called multipath effects. Figure 3 shows a simplified multipath channel model of underwater acoustic propagation, in which only the direct propagation path (P _iD,(i=1, 2, … )), reflection propagation path of the sea surface (P _iS,(i=1, 2, …)) and reflection propagation path of the seabed (P _iB,(i=1, 2, …)) are considered. In addition, noise exists in the underwater environment, and includes the radiated noise of the AUV, oceanic reverberation and other factors that interfere with acoustic signals. This condition results in a number of similar peaks that correspond to different time differences. However, the positioning error is amplified after the time differences are multiplied by the speed of the acoustic propagation (Deng et al., Reference Deng, Bao and Zhao2009).

Figure 3. Simplified channel mode for underwater acoustic multipath propagation.

Mahdinejad and Seghaleh (Reference Mahdinejad and Seghaleh2013) used a different weighted generalised cross-correlation method to calculate time delays and compared the experimental results in a reverberation room under an acoustic environment. They determined that the SCOT-weighted generalised cross-correlation algorithm obtains good results in acoustic signal processing and time delay estimation, and this method has better accuracy and precision than other methods. In this study, the SCOT-weighted generalised cross-correlation is used to make the relevant peaks prominent to effectively enhance the spectral component of the source signal in the received signal. Therefore, the signal-to-noise ratio is improved, and the cross-correlation function is enhanced. These improvements increase the accuracy of the estimation of the time difference.

The received signals of the hydrophone at positions P_i and P_j are transformed from the time domain into the frequency domain by using Fourier transformation:

(6)$$\matrix{F_{x_{i}}\lpar \omega \rpar =\int_{-\infty }^{+\infty } {x_{i}\lpar t \rpar } e^{-j\omega t}dt\cr F_{x_{j}}\lpar \omega \rpar =\int_{-\infty }^{+\infty } {x_{j}\lpar t \rpar } e^{-j\omega t}dt}$$

The function of cross-power spectral density is obtained by the following equation:

(7)$$G_{x_{i}x_{j}}\lpar \omega \rpar = F_{x_{i}}\lpar \omega \rpar F_{x_{j}}^{\ast }\lpar \omega \rpar $$

where $F_{x_{j}}^{\ast} \lpar \omega \rpar $ is the conjugate of $F_{x_{j}}\lpar \omega \rpar $.

According to the theorem of Wiener-Khinchin (Cohen, Reference Cohen1998), the relationship between the function of cross-power spectral density and cross correlation is:

(8)$$R_{x_{i}x_{{\rm j}}}\lpar \tau \rpar =\int_0^\pi {G_{x_{i}x_{j}}\lpar \omega \rpar } e^{-j\omega \tau }d\omega$$

To reduce the multipath interference and noise in the selection of correct correlation peaks, the received signals are processed by the above mentioned SCOT-weighted method. Signals x_i and x_j are passed through a pre-filter, and the SCOT-weighted general cross-correlation operation is performed on the pre-filtered signals. The cross-correlation function based on weighted function $\varphi_{x_{i}x_{j}}\lpar \omega \rpar $ in the frequency domain is:

(9)$$R_{x_{i}x_{{\rm j}}} \lpar \tau \rpar = \int_{\rm 0}^{\pi} \varphi_{x_{i} x_{j}} \lpar \omega \rpar G_{x_{i} x_{j}} \lpar \omega \rpar e^{-j\omega \tau} d \omega$$

The flow chart of the SCOT-weighted generalised cross correlation is shown in Figure 4.

Figure 4. Flow chart of the SCOT frequency domain-weighted generalised cross-correlation algorithm.

The SCOT-weighted generalised cross-correlation function improves the proportion of effective spectra in the received signal, which suppresses the noise interference and improves the accuracy of the estimated time difference. The abscissa value that corresponds to the maximum peak in the SCOT-weighted generalised cross-correlation function is the time difference τ.

4.1. Simulation to verifying the validity of the SCOT weighted function

The Bellhop model (Michael, Reference Michael2011) is widely used in underwater acoustic channel simulation. This model is equally effective in channel analysis and modelling. It simulates underwater environment noise, multipath effects, and Doppler frequency shifts by introducing sound field data. It thus forms a real channel model and obtains various practical data, such as eigen line, impulse response, propagation loss and arrival time sequence of sound.

To verify the validity of the SCOT-weighted cross-correlation function in the frequency domain, we performed simulations and compare the method with traditional cross-correlation methods by establishing a propagation model of underwater acoustic signals using the Bellhop software. The sound source signal was selected as a type of amplitude modulation signal with a bandwidth of 40 kHz and a centre frequency of 22 kHz. The sound source was placed at a depth of 50 m, and two hydrophones were placed at a depth of 25 m underwater at different locations. The results of the two methods of cross-correlation functions are shown in Figures 5 and 6. The simulations were conducted three times, and the results are subsequently compared (Table 1).

Figure 5. Basic cross-correlation function in simulation.

Figure 6. SCOT-weighted cross-correlation function in simulation.

Table 1. Comparisons of time difference errors in simulations.

Figures 5 and 6 and Table 1 show that the SCOT-weighted cross-correlation function in the frequency domain reduces the amplitude of pseudo-peaks near the ideal peak that corresponds to the time difference by eliminating the inference of multiple correlation peaks caused by the multipath propagation effect. The result reveals that the time difference estimated by this method is more accurate than that obtained by the traditional method.

4.2. Semi-physical experiment to verify the validity of the SCOT-weighted function

To evaluate the performance of the above method, we conducted a semi-physical experiment involving an ultrasonic sound source and five receivers in an open space with low noise. The four receivers formed a square shape with a side length of 5 m, and the remaining receiver was placed in the centre of the square. The equipment layout is shown in Figure 7, and the experimental environment is shown in Figure 8.

Figure 7. Equipment layout.

Figure 8. Experimental environment.

The sound source signal is a type of pulse signal with a frequency of 40 kHz and a pulse width of 20 ms. The traditional cross-correlation and SCOT-weighted cross-correlation methods in the frequency domain were applied for the signals received by two receivers. The results are shown in Figures 9 and 10 and Table 2.

Figure 9. Semi-physical basic cross-correlation function graph.

Figure 10. Semi-physical SCOT-weighted cross-correlation function graph.

Table 2. Comparison of time difference errors in the semi-physical experiment.

Table 3. Rocking parameters of the AUV.

The result of the traditional cross-correlation method reveals several pseudo-peaks with similar peaks. By contrast, the SCOT-weighted cross-correlation method in the frequency domain effectively reduces all the amplitudes of the pseudo-peaks and attains the maximum peak. The results of the semi-physical experiment show that the SCOT-weighted cross-correlation method significantly enhances the real correlation peak and effectively estimates the time difference.

4.3. Simulation of dynamic positioning

A simulation was conducted to verify the dynamic positioning effect of the time window model based on periodic movement. A single sound source was placed at a longitude of 120 · 01°, latitude of 40° and depth of 50 m. The sound source transmitted the pulse signal with a frequency of 0·025 Hz and a pulse width of 100 ms. The AUV moved along the northeast direction from its initial position of 119·995° longitude and 32·995° latitude at the speed of 1 m/s and depth of 25 m. The AUV swung in accordance with the three-axis sine model. The random and constant drifts of the gyroscopes were 0·05° h, and the random and constant biases of the accelerometers were 50 μg. The initial misalignment angles at the three axes were 1·5° and the period of the time window was set as 400 s. The simulation time was 3,600 s. In the simulation, the AUV position was estimated only on the basis of the SINS for the first 1,200 s. The algorithm based on a virtual LBL matrix window was introduced to calculate the position for the remaining 2,400 s.

Three algorithms, namely the algorithms based on SINS with no LBL aiding, virtual LBL matrix window using traditional cross-correlation function and SCOT-weighted cross-correlation function in the frequency domain, are compared in Figure 11. Table 4 shows the positioning errors of the three algorithms. As shown in Figure 11 and Table 4, the position obtained by SINS significantly deviates from the real trajectory. This result indicates that the accumulated errors of inertial sensors have a considerable influence on the positioning accuracy of the AUV. In the initial stage of the algorithm based on the virtual LBL matrix window, the accuracy of the position is rapidly obtained to within 3 m. The positioning accuracy was improved as the number of iterations increased for the virtual LBL matrix window algorithm. The positioning result became accurate and reliable when the number of iterations reached four and five. Then, the positioning error gradually increases because the distance from the AUV to the sound source becomes large. Compared with the traditional cross-correlation function, the SCOT-weighted cross-correlation function in the frequency domain performs better and meets the high accuracy requirements of underwater vehicle navigation positioning.

Figure 11. Comparison of two tracks with the ideal track.

Table 4. Comparison of pure inertia and acoustic cycle time window distance error.