Research on Passive Ranging Technology of Moving Ship Based on Vertical Hydrophone Array

: This paper introduces a method for the range localization of a moving ship based on the vertical hydrophone line array. The main implementation steps of the proposed algorithm are as follows. First, the stable low-frequency line spectrum component of the broadband radiated noise from the moving ship was extracted through the Detection of Envelope Modulation on Noise (DEMON) spectrum analysis method. Second, the pressure di ﬀ erence between the two di ﬀ erent ranges was derived, and the corresponding interference fringes were observed in the plane of time and time interval. Then, the radial velocity of the moving ship could be obtained based on the period of the pattern oscillations of the interference fringes. Further, we estimated the time and range information of the Closest Point of Approach (CPA) and computed the ship range versus time. Finally, each element of the vertical hydrophone line array was processed by the method proposed above, and data fusion technology was adopted to reduce the impact of ine ﬀ ective elements and improve the range estimation accuracy. The results of the simulation and experiment of a 16-element vertical array performed in the South China Sea veriﬁed the e ﬀ ectiveness of the algorithm.


Introduction
Range localization is a significant component of the Sound Navigation and Ranging (SONAR) system in the field of underwater target detection. Currently, the most popular methods for ranging are three subarrays, the matched field processing (MFP) and the waveguide invariant method. The three subarray passive ranging method mainly uses the principle of the curvature of the spherical wavefront. The target range is estimated by measuring the relative time delay of the received signal on each of the three subarrays [1]. However, this method is easily affected by the time delay estimation accuracy, whose error cannot be ignored, especially in the complex ocean environment [2]. The matching field passive positioning technology is an environment-dependent method based on ocean channel modelling [3]. The copy sound field is first calculated through the ocean channel model, and the matching processing is applied in the copy field and signals received by the hydrophone array. We can obtain the range and depth information of the sound source after the parameter search processing [4,5]. However, the mismatch problem of the ocean channel environment cannot be solved effectively in practice, and thus the matching field positioning accuracy becomes poor. Comparatively, the waveguide invariant technology is more tolerant of the environmental demands, which is an important over the time interval. Finally, the weighted least square (WLS) fusion method is applied to fuse the multi-hydrophone ranging results. The flow diagram is as shown in Figure 1.

DEMON Spectrum Analysis for A Single Hydrophone
As mentioned above, the DEMON spectral analysis is widely used in the feature extraction of ship-radiated noise. Generally speaking, the high-frequency band of the vessel-radiated noise has a periodic modulation, and the modulation frequency is composed of low-frequency characteristic signals such as the shaft frequency and blade frequency. The shaft frequency and blade frequency contain the information on the propeller speed and blade number. Normally, the ship-radiated noise is a periodic and locally stable process, whose expression in the time domain can be expressed as: where   ft represents the modulation signal and is a low-frequency, slowly varying periodic function.

 
gt is the carrier frequency and is a stable White Gaussian random process with a limited bandwidth.   xt is a periodic stationary stochastic Gaussian process for ship-radiated noise. The modulation signal   ft can be obtained by demodulating the ship noise   xt . This paper adopts the absolute value demodulation technology of the incoherent demodulation method. The specific flow diagram has been presented in Figure 2. The process includes four steps: band-pass filtering, the operation of taking the absolute value, the low-pass filter, the spectrum analysis.
We take a simple situation as an example to explain the operation principle. In this case, the amplitude of the high-frequency carrier is modulated, and the carrier signal is a single frequency signal. The amplitude modulation (AM) signal is expressed as:

DEMON Spectrum Analysis for A Single Hydrophone
As mentioned above, the DEMON spectral analysis is widely used in the feature extraction of ship-radiated noise. Generally speaking, the high-frequency band of the vessel-radiated noise has a periodic modulation, and the modulation frequency is composed of low-frequency characteristic signals such as the shaft frequency and blade frequency. The shaft frequency and blade frequency contain the information on the propeller speed and blade number. Normally, the ship-radiated noise is a periodic and locally stable process, whose expression in the time domain can be expressed as: where f (t) represents the modulation signal and is a low-frequency, slowly varying periodic function. g(t) is the carrier frequency and is a stable White Gaussian random process with a limited bandwidth.
x(t) is a periodic stationary stochastic Gaussian process for ship-radiated noise. The modulation signal f (t) can be obtained by demodulating the ship noise x(t). This paper adopts the absolute value demodulation technology of the incoherent demodulation method. The specific flow diagram has been presented in Figure 2.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 16 over the time interval. Finally, the weighted least square (WLS) fusion method is applied to fuse the multi-hydrophone ranging results. The flow diagram is as shown in Figure 1.

DEMON Spectrum Analysis for A Single Hydrophone
As mentioned above, the DEMON spectral analysis is widely used in the feature extraction of ship-radiated noise. Generally speaking, the high-frequency band of the vessel-radiated noise has a periodic modulation, and the modulation frequency is composed of low-frequency characteristic signals such as the shaft frequency and blade frequency. The shaft frequency and blade frequency contain the information on the propeller speed and blade number. Normally, the ship-radiated noise is a periodic and locally stable process, whose expression in the time domain can be expressed as: where   ft represents the modulation signal and is a low-frequency, slowly varying periodic function.

 
gt is the carrier frequency and is a stable White Gaussian random process with a limited bandwidth.   xt is a periodic stationary stochastic Gaussian process for ship-radiated noise. The modulation signal   ft can be obtained by demodulating the ship noise   xt . This paper adopts the absolute value demodulation technology of the incoherent demodulation method. The specific flow diagram has been presented in Figure 2. The process includes four steps: band-pass filtering, the operation of taking the absolute value, the low-pass filter, the spectrum analysis.
We take a simple situation as an example to explain the operation principle. In this case, the amplitude of the high-frequency carrier is modulated, and the carrier signal is a single frequency signal. The amplitude modulation (AM) signal is expressed as: where m A and c A are the amplitude of the carrier signal and modulation signal, respectively.
Moreover, we assume that m A and c A are both greater than one. m f and c f are the frequency of the carrier signal and modulation signal, respectively. Taking the absolute value of Equation (2), the AM signal can be written as: Then, the content of the absolute value is expanded into a Fourier series, that is: The process includes four steps: band-pass filtering, the operation of taking the absolute value, the low-pass filter, the spectrum analysis.
We take a simple situation as an example to explain the operation principle. In this case, the amplitude of the high-frequency carrier is modulated, and the carrier signal is a single frequency signal. The amplitude modulation (AM) signal is expressed as: where A m and A c are the amplitude of the carrier signal and modulation signal, respectively. Moreover, we assume that A m and A c are both greater than one. f m and f c are the frequency of the carrier signal and modulation signal, respectively. Taking the absolute value of Equation (2), the AM signal can be written as: Then, the content of the absolute value is expanded into a Fourier series, that is: where A n and B n are the constant coefficients. The multiplication of Equations (4) and (5) will produce various harmonic components. After passing through the designed low-pass filter, only the modulated signal remains. Then, the spectrum analysis of the low-frequency modulation signal is carried out after a Fourier transform. According to the spectrum distribution, we extract the line spectrum components representing the characteristics of the ship for the next step of the ranging analysis.

Range Estimation by Sound Intensity of The Pressure Difference
In the acoustic field, the intensity of the pressure difference received at different ranges contains information about the range interval, from which the radial velocity of the moving target can be obtained. Then, the range over time can be estimated.

Radial Velocity Derivation
According to the normal wave model theory, the sound pressure at a certain receiving depth z and range r in the sound field can be expressed as [20]: where z s is the source depth, ρ is the density of sea water and m stands for the pattern number. Ψ m (z) and Ψ m (z s ) are the eigenfunction of the m-th pattern for the source depth and receiving depth, respectively. k rm represents the eigenvalue. Similarly, the sound pressure at a close range r + ∆r can be expressed as: where ∆r is a small amount. We define the middle range as R = r + ∆r/2 and assume ∆r/2 R. Therefore, the pressure difference between the range r and r + ∆r in the same receiving depth is [7]: ∆p(R, ∆r) = p(r + ∆r) − p(r) B m e ik rm R sin k rm is a constant, and B m = is unrelated to the range interval ∆r. The sound intensity of the pressure difference can be indicated as: It is noted that the period of the cos (k rm − k rn ) ∆r 2 is much longer than that of the term cos (k rm + k rn ) ∆r 2 . Thus, the sound intensity of the pressure difference is mainly dominated by the latter term. The simplified expression is as follows: Further, we define k ri = δk i + k r , i = m, n, where δk i is a minimal quantity and k r is the average wave number of the modes. Its value is between k min = ω c b and k max = ω min(c w ) , where c b and c w are the sound speed at the bottom surface and in the water layer, respectively. Thus, Equation (10) can be further expressed as: (11) where the range interval is also expressed by ∆r(t) = v s (t)∆t, and v s (t) stands for the radial velocity of the moving target varying from time t. From Equation (10), we find that the interference pattern is mainly determined by the term cos k r v s (t)∆t . When meeting the following relationship: I ∆ (R, ∆r) will achieve the maximum. The sound intensity of the pressure difference can be calculated in another form. From Equation (9), we find that the intensity can be described as: where I r and I r+∆r stand for the intensity of p(r) and p(r + ∆r), respectively. The real part of the cross term p(r + ∆r) * · p(r) has the same interference pattern as Equation (11) [7]. Thus, we can conclude the following relationship: The average wave number k r is related to the frequency f and satisfies the expression: where v p represents the average phase speed. In this paper, we only consider the single-frequency condition. Substituting Equation (15) into Equation (12), we obtain the radial velocity of the moving target as: where d t is the peak interval in the (∆t, t) plane and can be obtained from the interference fringes. Generally speaking, v p approximates the speed of the sound field.

Range Estimation
We consider one normal situation where the source moves from far to near with a constant moving speed v 0 . When the source reaches the point nearest to the receiving point, the range is denoted as r 0 and the corresponding time is denoted as t 0 . At this time, the point is named as Closest Point of Approach (CPA). The trajectory of the moving source is shown in Figure 3.
When r 0 and t 0 are known, the range r of any point can be expressed as follows according to the geometric relationship: where t is the moving time.
The time at any point on the trajectory is recorded as t, and the time at a very close range is denoted as t − δt. The range between this point and the receiving point is denoted as r(t − δt). Additionally, the radial velocity v s (t) can be denoted as: Appl. Sci. 2020, 10, 7396 6 of 16 First, the appropriate range of search values are set for v 0 , t 0 and r 0 . Then, a set of values are substituted into Equation (18) to calculate v s . This value is compared with the radial velocity obtained from Equation (16), and the mean square error (MSE) is calculated. v 0 , t 0 and r 0 corresponding to the minimum MSE are the estimated moving velocity, time and range of the CPA point. Last, depending on the values of the CPA point, the moving range over time is obtained by Equation (17).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 16 where t is the moving time.
The time at any point on the trajectory is recorded as t , and the time at a very close range is denoted as tt   . The range between this point and the receiving point is denoted as   Additionally, the radial velocity   s vt can be denoted as:

Data Fusion Technology Based on Weighted Least Squares
Due to the complex ocean environment, the range localization result for a single hydrophone is not accurate. Therefore, we adopt data fusion technology based on the weighted least squares to improve the localization accuracy.
Here, we assume that a vertical hydrophone array is composed of N elements. Each element collects sound pressure information on the ship noise separately. The signal received by each hydrophone is processed by the DEMON spectrum, and then the estimated range is obtained by the sound intensity of the pressure difference.
We regard ˆj  as the estimated range of the j -th hydrophone, and the estimation range vector of the vertical hydrophone line array can be written as: The estimated ranges of different sensors are fused to obtain a more precise target location. According to the linear fusion model, the estimated range from the source to the sensors can be

Data Fusion Technology Based on Weighted Least Squares
Due to the complex ocean environment, the range localization result for a single hydrophone is not accurate. Therefore, we adopt data fusion technology based on the weighted least squares to improve the localization accuracy.
Here, we assume that a vertical hydrophone array is composed of N elements. Each element collects sound pressure information on the ship noise separately. The signal received by each hydrophone is processed by the DEMON spectrum, and then the estimated range is obtained by the sound intensity of the pressure difference.
We regardθ j as the estimated range of the j-th hydrophone, and the estimation range vector of the vertical hydrophone line array can be written as: The estimated ranges of different sensors are fused to obtain a more precise target location. According to the linear fusion model, the estimated range from the source to the sensors can be denoted as [21]: where x denotes the state parameters to be measured. e = [e 1 , e 2 , · · · e N ] T represents the N × 1 dimensional error and satisfies the Gaussian distribution, namely E e j = 0 , j = 1, 2 · · · N. H is the N × 1 dimensional mapping matrix between r and x. Generally, H = [1, 1, · · · 1] T because of the same type of each hydrophone. However, the estimation range of each hydrophone has a different accuracy due to complex environments in the deep sea. Thus, we adopt the weighted least squares (WLS) method to eliminate the uncertainty. That is, for the high-accuracy element the weight is set to be larger, but for the low-accuracy element the weight is smaller. After using the WLS method, the estimation range becomes more accurate.
We denote the weighted coefficient of N hydrophones as w = [w 1 , w 2 , · · · , w N ], which satisfies the relation N j=1 w j = 1, 0 ≤ w j ≤ 1. In addition, the weight needs to satisfy the relationship E w T r − x = 0. Moreover, the estimation error is lower than the variance of estimated ranges of any single element. The sum of the square of weighted errors can be expressed as: where W is a weight vector and satisfies the expression W = diag(w 1 , w 2 , · · · , w N ). According to the WLS rule, the optimal weight W can be acquired when J w (x) reaches the minimum. Take the partial derivative of Equation (21) and set it to zero, namely: The value of the estimated range after optimal WLS fusion can be written as: Further, the variance after WLS data fusion is expressed as: where σ j 2 is the estimated variance of the j-th sensor. Based on the extreme value theory, the weight corresponding to the minimum total mean square error (MSE) is: Substituting Equation (25) into Equation (23), the ultimate range result after fusion is derived. At this time, the corresponding fusion result has the highest accuracy, and its root mean square error (RMSE) is: Equation (26) indicates that the variance after fusion is lower than that of any single hydrophone. The WLS data fusion method can effectively get rid of the elements with inaccurate range estimation results and improve the estimated accuracy significantly.

Simulation
In this section, a simulation is carried out to verify the performance of the above-proposed method. The concrete simulation setting is as follows. The source radiates a broadband signal with a frequency band from 200 Hz to 1000 Hz and moves with a constant speed of 3.5 m/s. The source depth is 50 m. The range and time of the CPA point are 500 m and 500 s, respectively. The moving time of the ship is from 0 s to 1800 s, and the corresponding range between the sound source and the receiving hydrophone is from 500 m to 4568 m. Figure 4 shows the curve of the target range versus time.
The concrete simulation setting is as follows. The source radiates a broadband signal with a frequency band from 200 Hz to 1000 Hz and moves with a constant speed of 3.5 / ms . The source depth is 50 m. The range and time of the CPA point are 500 m and 500 s, respectively. The moving time of the ship is from 0 s to 1800 s, and the corresponding range between the sound source and the receiving hydrophone is from 500 m to 4568 m. Figure 4 shows the curve of the target range versus time. A vertical hydrophone array is adopted to collect the pressure information of the ship-radiated noise. It consists of 16 elements. The hydrophones at different depths receive different signals due to a variety of scattering and reflection in the sound field. The range estimation results obtained by each element are also different. In order to reflect the differences in the received signal, we designed different signal-to-noise ratios (SNR) for the hydrophones at different depths. The algorithm proposed in this paper makes no requirements on the sound field model and is not sensitive to environmental parameters. In order to prove that the method is adaptable to different sound field environments, we used a sound speed profile (SSP), which was measured at a sea depth of 2200 m. The concrete parameter settings such as the water density and sea floor attenuation coefficient are based on the experience value and presented in Figure 5.  A vertical hydrophone array is adopted to collect the pressure information of the ship-radiated noise. It consists of 16 elements. The hydrophones at different depths receive different signals due to a variety of scattering and reflection in the sound field. The range estimation results obtained by each element are also different. In order to reflect the differences in the received signal, we designed different signal-to-noise ratios (SNR) for the hydrophones at different depths. The algorithm proposed in this paper makes no requirements on the sound field model and is not sensitive to environmental parameters. In order to prove that the method is adaptable to different sound field environments, we used a sound speed profile (SSP), which was measured at a sea depth of 2200 m. The concrete parameter settings such as the water density and sea floor attenuation coefficient are based on the experience value and presented in Figure 5.
frequency band from 200 Hz to 1000 Hz and moves with a constant speed of 3.5 / ms . The source depth is 50 m. The range and time of the CPA point are 500 m and 500 s, respectively. The moving time of the ship is from 0 s to 1800 s, and the corresponding range between the sound source and the receiving hydrophone is from 500 m to 4568 m. Figure 4 shows the curve of the target range versus time. A vertical hydrophone array is adopted to collect the pressure information of the ship-radiated noise. It consists of 16 elements. The hydrophones at different depths receive different signals due to a variety of scattering and reflection in the sound field. The range estimation results obtained by each element are also different. In order to reflect the differences in the received signal, we designed different signal-to-noise ratios (SNR) for the hydrophones at different depths. The algorithm proposed in this paper makes no requirements on the sound field model and is not sensitive to environmental parameters. In order to prove that the method is adaptable to different sound field environments, we used a sound speed profile (SSP), which was measured at a sea depth of 2200 m. The concrete parameter settings such as the water density and sea floor attenuation coefficient are based on the experience value and presented in Figure 5.  As described above, the high-frequency signal is modulated by a low-frequency signal. However, the low-frequency line spectrum in the wideband spectrum is unstable, so a DEMON spectrum analysis is used to extract the low-frequency line spectrum components. In this case, we set the modulation frequency to f = 14 Hz and employ the normal mode model to simulate the pressure signal. Figure 6 presents the wave of the broadband noise in the time domain and the DEMON spectrum result demodulated for a single element. It is clear that the wideband signal has a low-frequency envelope, and the demodulated line spectrum is stable at 14 Hz.
analysis is used to extract the low-frequency line spectrum components. In this case, we set the modulation frequency to 14 f Hz  and employ the normal mode model to simulate the pressure signal. Figure 6 presents the wave of the broadband noise in the time domain and the DEMON spectrum result demodulated for a single element. It is clear that the wideband signal has a lowfrequency envelope, and the demodulated line spectrum is stable at 14 Hz. Then, the estimated radial velocity versus time is computed by Equation (16). As presented in Figure  7b, the yellow line represents the radial velocity estimation result, and the blue line stands for the actual radial velocity. For the CPA point, the estimated radial velocity, the range and the time are 3.7 m/s, 555 m and 496 s. Note that there is a widening of the estimated radial velocity near the CPA point because the radial component of the moving velocity in the region near the CPA point is quite small and even close to zero. Figure 7c compares the estimated range and actual range, which are represented by the solid red line and the blue dotted line, respectively. The comparison result shows that the two curves basically coincide. The two curves are relatively close, especially at the CPA point. However, the deviation increases over time. In order to describe the difference between the estimation value and true value at each time, we introduce a parameter, namely the relative error of the range estimation. It is defined as:  The interference fringes with the changing time can be observed when the time interval ∆t varies from 1 s to 200 s, as shown in Figure 7a. The period of the interference oscillations d t can be obtained through calculating the Fourier transform of the sound intensity of the pressure difference. Then, the estimated radial velocity versus time is computed by Equation (16). As presented in Figure 7b, the yellow line represents the radial velocity estimation result, and the blue line stands for the actual radial velocity. For the CPA point, the estimated radial velocity, the range and the time are 3.7 m/s, 555 m and 496 s. Note that there is a widening of the estimated radial velocity near the CPA point because the radial component of the moving velocity in the region near the CPA point is quite small and even close to zero. Figure 7c compares the estimated range and actual range, which are represented by the solid red line and the blue dotted line, respectively. The comparison result shows that the two curves basically coincide. The two curves are relatively close, especially at the CPA point. However, the deviation increases over time. In order to describe the difference between the estimation value and true value at each time, we introduce a parameter, namely the relative error of the range estimation. It is defined as: where r is the range estimation, and R true is the true range versus time. The corresponding calculation result is shown in Figure 7d. One can see that the relative error of the estimated result is stable at about 0.1 and changes little with time.
In the application, the range estimation result using only a single element may be inaccurate because the ocean environment is very complicated. Hence, the vertical line array composed of multiple hydrophones at different depths is used in the simulation, and further the data fusion method is adopted. The comparison results of the estimated ranges for each hydrophone and the fusion result are shown in Figure 8.
The solid black lines represent the ranging results of each element. The solid red line is the 16-element fusion result, and the dotted blue line is the actual range with time. This figure indicates that the estimation results of individual elements have a large bias. In particular, the estimated CPA point is very different from the true point, which seriously affects the accuracy of the range estimation. The data fusion method assigns a small weight to inaccurate results to improve the estimation accuracy. After data fusion, the fusion curve is in better agreement with the actual curve. In the application, the range estimation result using only a single element may be inaccurate because the ocean environment is very complicated. Hence, the vertical line array composed of multiple hydrophones at different depths is used in the simulation, and further the data fusion method is adopted. The comparison results of the estimated ranges for each hydrophone and the fusion result are shown in Figure 8.   In the application, the range estimation result using only a single element may be inaccurate because the ocean environment is very complicated. Hence, the vertical line array composed of multiple hydrophones at different depths is used in the simulation, and further the data fusion method is adopted. The comparison results of the estimated ranges for each hydrophone and the fusion result are shown in Figure 8.

Experiment Situation
In 2016, our research group conducted a sea trial in the South China Sea. The experimental configuration diagram is shown in Figure 9. This experiment employed a 16-element vertical hydrophone line array. A moving ship was chosen as the sound source, and the hydrophone array received the broadband-radiated noise of the vessel. The concrete experiment parameters were set as follows: • The vertical hydrophone line array was composed of 16 elements with an interval of 5 m. The depth of the top hydrophone was 370 m, and that of the bottom sensor was about 445 m. The depth of water over the experimental area was about 1800 m.

•
The position where the heavy block attached below the array finally entered the water was regarded as the real array position recorded by the Global Positioning System (GPS) during the experiment. The moving vessel had GPS modules to record its real-time spatial positions. The corresponding trajectory map has been shown in Figure 10a. The real range from the ship to the array was calculated by GPS data and was then shown in Figure 10b. From this figure, we can see that the range and time of the CPA point were 344 m and 130 s, respectively. The moving velocity of the ship was about 5 m/s.

•
The array may have a small tilt due to the disturbance of ocean currents. However, this small angular offset has little effect (only within 10%) on the final ranging result.

•
The sound velocity profiler (SVP) was used to measure the sound speed profile (SSP) at the beginning of the experiments. The measured SSP is shown in Figure 11.

•
The sampling rate was 16 kHz. The sea condition was level 3, and no other vessels passed during the measurement time.
is very different from the true point, which seriously affects the accuracy of the range estimation. The data fusion method assigns a small weight to inaccurate results to improve the estimation accuracy. After data fusion, the fusion curve is in better agreement with the actual curve.

Experiment Situation
In 2016, our research group conducted a sea trial in the South China Sea. The experimental configuration diagram is shown in Figure 9. This experiment employed a 16-element vertical hydrophone line array. A moving ship was chosen as the sound source, and the hydrophone array received the broadband-radiated noise of the vessel. The concrete experiment parameters were set as follows:  The vertical hydrophone line array was composed of 16 elements with an interval of 5 m. The depth of the top hydrophone was 370 m, and that of the bottom sensor was about 445 m. The depth of water over the experimental area was about 1800 m.  The position where the heavy block attached below the array finally entered the water was regarded as the real array position recorded by the Global Positioning System (GPS) during the experiment. The moving vessel had GPS modules to record its real-time spatial positions. The corresponding trajectory map has been shown in Figure 10a. The real range from the ship to the array was calculated by GPS data and was then shown in Figure 10b. From this figure, we can see that the range and time of the CPA point were 344 m and 130 s, respectively. The moving velocity of the ship was about 5 m/s.  The array may have a small tilt due to the disturbance of ocean currents. However, this small angular offset has little effect (only within 10%) on the final ranging result.  The sound velocity profiler (SVP) was used to measure the sound speed profile (SSP) at the beginning of the experiments. The measured SSP is shown in Figure 11.  The sampling rate was 16 kHz. The sea condition was level 3, and no other vessels passed during the measurement time.

Results Analysis
The top hydrophone is chosen as the reference element for the below discussion, and the other 15 elements are processed in the same way. Then, the range estimation results of 16 hydrophones are fused to improve the estimation accuracy. We choose a piece of data with a length of 900 s for processing. Figure 12 shows the time domain diagram and time-frequency domain diagram of the received signal for the reference element.

Results Analysis
The top hydrophone is chosen as the reference element for the below discussion, and the other 15 elements are processed in the same way. Then, the range estimation results of 16 hydrophones are fused to improve the estimation accuracy. We choose a piece of data with a length of 900 s for processing. Figure 12 shows the time domain diagram and time-frequency domain diagram of the received signal for the reference element. From Figure 12a, we can see that the received signal has a low-frequency envelope. Figure 12b validates that the ship-radiated noise contains multiple frequency components. In addition, adjacent spectrum components are prone to aliasing, and the low-frequency line spectrum components directly extracted from the spectrum are unstable. Therefore, the DEMON spectrum analysis method is employed to extract the stable low-frequency line spectrum. As shown in Figure 13a, the signal frequency after the DEMON analysis is 14 Hz. From Figure 12a, we can see that the received signal has a low-frequency envelope. Figure 12b validates that the ship-radiated noise contains multiple frequency components. In addition, adjacent spectrum components are prone to aliasing, and the low-frequency line spectrum components directly extracted from the spectrum are unstable. Therefore, the DEMON spectrum analysis method is employed to extract the stable low-frequency line spectrum. As shown in Figure 13a, the signal frequency after the DEMON analysis is 14 Hz.
The time interval between two points ranges from 1 s to 200 s, and the corresponding interference fringes of the sound intensity of the pressure difference is shown in Figure 13b. The period of the interference oscillations is calculated through Fourier transform and is then applied to calculate the radial velocity. The estimated radial velocity changing over time is shown in Figure 13c. From this figure, we can find that the estimated radial velocity is about 4.6 m/s. Then, when the MSE between the estimated value and true value reaches the minimum, the corresponding range and time of the CPA point can be obtained. Through the above processing, the range and time of the CPA point are 190 m and 115 s, respectively. Based on the above analysis, the range versus time is estimated by Equation (17). Then, the estimated range and the range recorded by the GPS data are compared in Figure 13d. The two curves have similar trends, but there is a slight difference in the range r at the same time t. One can see that the range estimation result of a single hydrophone is inaccurate. Our experiment adopts a 16-element vertical hydrophone line array, and all of the elements are processed separately by the proposed method. Figure 14 shows the range estimation results of the 16 elements and the data fusion result. From Figure 12a, we can see that the received signal has a low-frequency envelope. Figure 12b validates that the ship-radiated noise contains multiple frequency components. In addition, adjacent spectrum components are prone to aliasing, and the low-frequency line spectrum components directly extracted from the spectrum are unstable. Therefore, the DEMON spectrum analysis method is employed to extract the stable low-frequency line spectrum. As shown in Figure 13a, the signal frequency after the DEMON analysis is 14 Hz. The time interval between two points ranges from 1 s to 200 s, and the corresponding interference fringes of the sound intensity of the pressure difference is shown in Figure 13b. The period of the interference oscillations is calculated through Fourier transform and is then applied to calculate the radial velocity. The estimated radial velocity changing over time is shown in Figure 13c. From this figure, we can find that the estimated radial velocity is about 4.6 m/s. Then, when the MSE between the estimated value and true value reaches the minimum, the corresponding range and time of the CPA point can be obtained. Through the above processing, the range and time of the CPA point are 190 m and 115 s, respectively. Based on the above analysis, the range versus time is estimated by Equation (17). Then, the estimated range and the range recorded by the GPS data are compared in Figure 13d. The two curves have similar trends, but there is a slight difference in the range r at the same time t . One can see that the range estimation result of a single hydrophone is inaccurate. Our experiment adopts a 16-element vertical hydrophone line array, and all of the elements are processed separately by the proposed method. Figure 14 shows the range estimation results of the 16 elements and the data fusion result. From this figure, we can see that the estimation results of most elements are close to the GPS true value and that there are three obvious curves with large deviations. When performing WLS fusion, we assign smaller weights to the estimation results with larger range deviations, and larger weights to the elements with smaller deviations. Therefore, the fusion result is in good agreement with the actual curve. From this figure, we can see that the estimation results of most elements are close to the GPS true value and that there are three obvious curves with large deviations. When performing WLS fusion, we assign smaller weights to the estimation results with larger range deviations, and larger weights to the elements with smaller deviations. Therefore, the fusion result is in good agreement with the actual curve.
To further measure the performance of this range localization method, we calculate the parameter of the relative mean square error (RMSE), defining RMSE as: where r i t j represents the estimated range at the j-th second by the i-th element. R t j is the real range at the j-th second. σ i 2 is the RMSE of the i-th element. M and N are the length of time and the total number of elements. In this paper, we are more concerned about the accuracy of the estimation results near the CPA point, because the range, time and radial velocity of the CPA point determine the total estimated range. Therefore, we choose the estimation results of the first 300 s of Figure 14 for statistical calculations. Table 1 gives the RMSE of each hydrophone and the RMSE after the 16-element fusion. From this table, we can see that the RMSEs of the 16 elements are slightly different. The RMSEs of the 7th, 10th and 13th element are significantly higher than those of the other elements, which is consistent with the curves in Figure 14. Moreover, the RMSE of the fusion results is 4.9 × 10 −4 , which is lower than the RMSE of most elements. Hydrophones at different depths receive signals with different SNRs. Inaccurate range estimation results can be obtained by a low SNR signal, as shown in the 7th, 10th and 13th element, because the sound signals at different depths go through different paths, undergoing different reflections and scattering. The signals reaching the receivers are attenuated differently, and the SNRs are also various.
The performance evaluation of the algorithm proposed in our paper is based on a comparison with real GPS data and the range estimation results obtained by a single hydrophone element. The comparison of these three results infer that this method can improve the estimation accuracy effectively and obtain a more accurate range estimation result than a single hydrophone.
It should be noted that the target and array should avoid being located in the shadow zone. In the shadow zone, the sound field is mainly composed of bottom-reflected rays rather than direct rays. The sound rays reflected by the sea bottom are greatly attenuated, and the SNRs of the signal received by the hydrophones are low. Thus, the range localization performance may be inaccurate.

Conclusions
In this paper, we proposed an algorithm to estimate the range of a moving ship. First, the DEMON spectrum analysis method is adopted to extract the low-frequency line spectrum components. Based on the sound intensity of the pressure difference, the extracted single-frequency signal at the given time interval interferes and forms the interference fringes in the plane of time and time interval. Then, the radial velocity of the moving source is derived from the period of interference oscillations. By fitting the time, velocity and range of the CPA point, we estimate the range from the source to the receiving hydrophone over time. Finally, range estimation results of multiple elements of the vertical line array are fused to obtain a more accurate estimation result. The simulation and deep-sea experiment analysis all verify that this method can estimate the range effectively. Moreover, the estimation accuracy is significantly improved since it fuses all the estimation results of each element and eliminates the impact of inaccurate results. The relative mean square error of the fusion result is lower than that of most hydrophones. In addition, this method does not require prior knowledge of marine environmental parameters and sound field modelling and thus demonstrates good universality. This method provides a new idea for the passive ranging of moving ships.

Conflicts of Interest:
The authors declare no conflict of interest.