Range Localization of a Moving Source Based on Synthetic Aperture Beamforming Using a Single Hydrophone in Shallow Water

: To localize a moving source in shallow water with a single hydrophone, a passive range localization method based on synthetic aperture beamforming is proposed. First, the horizontal wavenumber spectrum excited by the source is obtained by synthetic aperture beamforming. Then, according to the theoretical derivation (when the integration time is short, the maximum value of the horizontal wavenumber spectrum is related to the average horizontal wavenumber and the radial velocity of the source), the radial velocity can be obtained after obtaining the average horizontal wavenumber. Finally, in the case where there is a closest point of approach (CPA), the range can be recovered from estimation of the range and time of CPA, and from the constant source speed alone the linear track by ﬁtting the source velocity with the model of radial velocity. The only a priori information required is the sound velocity in water. The processing results using simulated data and SWellEx-96 experimental data show that the proposed method can e ﬀ ectively estimate the range of a moving source in shallow sea.


Introduction
Acoustic methods for localizing a source in an underwater acoustic waveguide, including direction finding and range and depth estimation, are a topic of great interest in the field of underwater acoustic signal processing. Popular methods for estimating the source range and depth are matched field processing (MFP) [1,2] and matched mode processing (MMP) [3,4], which require multiple hydrophones to form a vertical or horizontal array to receive the source radiated noise. Knowledge of the environment, such as the sound speed profile (SSP) and the geoacoustic parameters, is also needed to calculate the copy field vector for MFP or the modal coefficient vector for MMP using a propagation model to estimate source range and depth. The MFP and MMP methods are often plagued by environmental mismatch problems: if the assumed environment differs from the true environment, there may be large error in the final position estimate [5,6]. A number of robust localization algorithms have been developed to solve the mismatch problem [7][8][9][10], and will not be described here.
In actual use, the hydrophone array is expensive, and it is cumbersome to deploy at sea. Therefore, the idea of using a single hydrophone to perform source range localization emerged. A single hydrophone can be deployed more quickly and flexibly, at a lower cost, and monitoring of a large area

The Range Localization Method
We consider a range-independent shallow ocean environment and a moving source with a depth z s and constant speed v 0 during the observation time T. The radial velocity is given by v s (t) = δr δt (t) where r 0 and t 0 are the range at the closest point of approach (CPA) and the time of CPA, respectively, and δt is the variation of time associated to δr, the radial range variation for each acquisition time t. The normal mode theory expression of the pressure field received at a receiver at a depth z r and time t is where r t is the source range at time t, φ m (z) is the mth mode depth function, k m is the mth horizontal wavenumber, and α m is the mth mode attenuation coefficient. M is the number of modes that can be propagated in the current environment.
The pressure field at different times is beamformed as follows [17,18]: where t 1 and t 2 represent the start and end moments of the observation time T, respectively, and e ik r t is the steering vector and k r is a variable.
In [17] the observation time T should set to be long enough to separate different horizontal wavenumbers of modes. Here, we set the observation time in the opposite way. When the observation time T is short, the radial velocity can be considered to be constant, which is represented by v s . Within the observation time T, the source range r t can be expressed as where r t 1 represents the range at time t 1 . After substituting the normal mode expression into Equation (3), one obtains where and the wavenumbers can be expressed as k m = k + δk m , where k is the average horizontal wavenumber of the modes since δk m k (typically, δk m /k ∼ O(10 −2 ) for discrete mode wavenumbers, which are between k min = ω/c b and k max = ω/min(c w ) where c b and c w are the sound speeds at the bottom surface and in the water, respectively [14]. For some lower frequencies' cases, δk m /k ∼ 10 −1 ). Then g(k r ) can be written as It can be seen from Equation (6) that the horizontal wavenumber spectrum g(k r ) will have maxima for k r = kv s , noting that the mode attenuation coefficient is generally orders of magnitude smaller than the horizontal wavenumber. The average phase speed v p can be approximated by the sound speed in water: v p ≈ c w , and we can easily obtain the average horizontal wavenumber since we have that k( f ) = 2π f /v p for a particular frequency (narrowband case). Furthermore, the average radial velocity within the observation time T can be obtained as In the calculation process, the term √ r t in Equation (3) can be considered to be constant within a short observation period T, and thus can be set to 1 without affecting the subsequent results. Further, Equation (3) can be conveniently calculated using a method such as fast Fourier transform (FFT).
After obtaining the radial velocity, we can easily obtain the source speed v 0 , the range and time of CPA, r 0 and t 0 , according to Equation (1) using the least squares criterion and then calculate the source range as follows For a broadband case, the frequency processing gain can be obtained in a manner of sub-band division and non-coherent accumulation. Assume that the number of sub-bands divided is L, and the center frequency corresponding to each sub-band is f l , where l = 1, 2, · · · , L. Each sub-band is processed according to Equation (3) to obtain the sub-band horizontal wavenumber spectrum g(k r ) f l , where l = 1, 2, · · · , L. Since the average horizontal wavenumbers corresponding to different frequencies are different, the peak positions of the sub-band horizontal wavenumber spectrums are also different. Therefore, the sub-band horizontal wavenumber spectrum cannot be directly added together. Here, we choose any one of the sub-band center frequencies as the reference frequency f ref , and convert the horizontal wavenumber spectrum of each sub-band into a form relative to the reference frequency by Appl. Sci. 2020, 10, 1005 4 of 12 stretching/compression transformation ( g(k r ) f l ⇒ g(k r ( f l / f ref )) f l ). Finally, we can obtain the following broadband horizontal wavenumber spectrum where k r,ref is the variable corresponding to the reference frequency f ref . The broadband horizontal wavenumber spectrum is similar to the narrowband horizontal wavenumber spectrum at a frequency f ref , so that the average radial velocity can be estimated as shown in Equation (10), similarly to the process of obtaining Equation (7) v s = k r,ref where the average horizontal wavenumber It can be seen from Equation (7) that the approximation of the average horizontal wavenumber by the wavenumber in water causes a bias on the estimation of the average radial velocity, which will in turn affect the accuracy of the range estimation. Assuming that the error of the average horizontal wavenumber is δk, then the estimation error of the average radial velocity is given here [14] ∆v and the maximum error of the horizontal wavenumber is δk = k max − k min . According to the previous description, δk/k is small since δk k, so the radial velocity estimation error is also small, ensuring the accuracy of the range localization.
So far, we have completed the theoretical derivation of the source range localization. The effectiveness of the method will be verified using simulation and at-sea data in the next two sections.

Simulation and Results
We consider an 88 m depth waveguide with a negative gradient SSP as shown in Figure 1. It consists of a 10 m isothermal mixed layer with a sound speed of 1533 m/s. The sound speed reduces linearly to 1478 m/s at a depth of 40 m and then remains unchanged until the bottom at 88 m. The bottom sound speed is 1650 m/s, the bottom density is 1.76 g/cm 3 , and the bottom attenuation is 0.8 dB/λ. We assume that the source speed is v 0 = 2.0 m/s, the range of CPA is r 0 = 1000 m, the time of CPA is t 0 = 1000 s, the source depth is 50 m, a narrowband signal of 350 Hz is emitted, and the depth of the receiver is 70 m. The received sound pressure is calculated using the Kraken normal mode model [19]. The relationship between the real part of the sound pressure and time is shown in Figure 2a. After zooming in to reveal the details, it can be seen from Figure 2b that the real part of the sound pressure exhibits obvious periodicity. According to the above theoretical derivation, this period is related to the average horizontal wavenumber and the radial velocity. The received sound pressure is calculated using the Kraken normal mode model [19]. The relationship between the real part of the sound pressure and time is shown in Figure 2a. After zooming in to reveal the details, it can be seen from Figure 2b that the real part of the sound pressure exhibits obvious periodicity. According to the above theoretical derivation, this period is related to the average horizontal wavenumber and the radial velocity. The received sound pressure is calculated using the Kraken normal mode model [19]. The relationship between the real part of the sound pressure and time is shown in Figure 2a. After zooming in to reveal the details, it can be seen from Figure 2b that the real part of the sound pressure exhibits obvious periodicity. According to the above theoretical derivation, this period is related to the average horizontal wavenumber and the radial velocity. In the simulation, the observation period is set to be 0.5s T = , and Equation (3)  can be obtained according to Equation (7). Following this method, the radial velocity in the observation time of 6000 seconds is obtained as shown in Figure 4, from which we can see that the estimated radial velocity remains consistent with its theoretical value (the black solid line in the figure). After calculating the source speed 0 v , and the range and time of CPA, 0 r and 0 t , using Equation (1), we can obtain the source range estimation using Equation (8).
The result is shown in Figure 5. In the simulation, the observation period is set to be T = 0.5 s, and Equation (3) is calculated using the 128-point FFT. The horizontal wavenumber spectrum of the first-time step is obtained and shown in Figure 3. It can be seen that the maximum value of the spectrum appears at k r = 2.651. We set the average phase velocity to be v p = 1478 m/s and obtain k = 1.4879 m −1 . Then an estimated value of the radial velocity v s = 1.7817 m/s can be obtained according to Equation (7). Following this method, the radial velocity in the observation time of 6000 seconds is obtained as shown in Figure 4, from which we can see that the estimated radial velocity remains consistent with its theoretical value (the black solid line in the figure). After calculating the source speed v 0 , and the range and time of CPA, r 0 and t 0 , using Equation (1), we can obtain the source range estimation using Equation (8). The result is shown in Figure 5.          It can be seen from the above analysis that the range localization result obtained by the simulation data is good. The processing result of the SWellEx-96 experimental data is given below.

Experiment Description
Data from the SWellEx-96 experiment were used. The SWellEx-96 experiment was conducted by the Marine Physical Lab in the littoral waters just outside the port of San Diego in May 1996. The environmental information and data are available online [20]. Therefore, they are widely used to verify the performance of underwater target detection [21], localization [22], tracking [23], and underwater acoustic inversion [24] methods.
The receiving arrays used in the experiment include a vertical line array (VLA), a tilt line array (TLA), and two horizontal line arrays (HLA). The experiment is divided into two events: S5 and S59. Sound sources in both events travel along the contour line and are suitable for signal processing in a range-independent sound field environment, and event S59 has one more interference source than event S5. Here, the data of the vertical array from event S5 are used to verify the performance of the proposed range localization method. Figure 6 shows the ship trajectory for event S5, with two sources being towed at depths of 54 m and 9 m, respectively. The deeper source emits comb-like single-frequency signals at various signal-to-noise ratios (SNR) and frequencies from 49 to 400 Hz, and the shallower source emits 9 single-frequency signals with a frequency range of 109 to 385 Hz. The sampling rate is 1500 Hz. The SSP was measured using a CTD several times during the 75 min duration of the experiment. underwater acoustic inversion [24] methods.
The receiving arrays used in the experiment include a vertical line array (VLA), a tilt line array (TLA), and two horizontal line arrays (HLA). The experiment is divided into two events: S5 and S59. Sound sources in both events travel along the contour line and are suitable for signal processing in a range-independent sound field environment, and event S59 has one more interference source than event S5. Here, the data of the vertical array from event S5 are used to verify the performance of the proposed range localization method. Figure 6 shows the ship trajectory for event S5, with two sources being towed at depths of 54 m and 9 m, respectively. The deeper source emits comb-like single-frequency signals at various signal-to-noise ratios (SNR) and frequencies from 49 to 400 Hz, and the shallower source emits 9 single-frequency signals with a frequency range of 109 to 385 Hz. The sampling rate is 1500 Hz. The SSP was measured using a CTD several times during the 75 min duration of the experiment.

Data Processing Results and Analysis
In the following, the source range is estimated from the analysis of the signal received at the deepest array sensor. The results given below are similar when processing the signal measured at

Data Processing Results and Analysis
In the following, the source range is estimated from the analysis of the signal received at the deepest array sensor. The results given below are similar when processing the signal measured at the other array elements. As shown in Figure 8, the time-frequency diagram of the measured signal shows that there is a suspected strong interference in the frequency band below 100 Hz around the 10th minute, which will have an impact on the radial velocity estimation based on single-frequency signal processing. Around the 60th minute, there are strong interference patterns in the~500 to 700 Hz frequency band, which is caused by the noise radiated by the source ship.

Data Processing Results and Analysis
In the following, the source range is estimated from the analysis of the signal received at the deepest array sensor. The results given below are similar when processing the signal measured at the other array elements. As shown in Figure 8, the time-frequency diagram of the measured signal shows that there is a suspected strong interference in the frequency band below 100 Hz around the 10th minute, which will have an impact on the radial velocity estimation based on single-frequency signal processing. Around the 60th minute, there are strong interference patterns in the ~500 to 700 Hz frequency band, which is caused by the noise radiated by the source ship.

Narrowband Source Range Localization
For the narrowband case, ( ) pt in Equation (2) is calculated and the 127 Hz line is extracted from the spectrogram of the time series that was acquired continuously during the 75 min of the experiment, which is shown in Figure 9. The size of the window for the FFT is 1 s and the number of FFT points is 1500.

Narrowband Source Range Localization
For the narrowband case, p(t) in Equation (2) is calculated and the 127 Hz line is extracted from the spectrogram of the time series that was acquired continuously during the 75 min of the experiment, which is shown in Figure 9. The size of the window for the FFT is 1 s and the number of FFT points is 1500. can be obtained according to Equation (7). Following this method, the radial velocity in the 75 min duration of the experiment is obtained as shown in Figure 11, wherein the black solid line is the maximum point. It can be seen that the estimation of the radial velocity is good for most of The time step is set to be T = 1 s, and Equation (3) is calculated using the 128-point FFT. The horizontal wavenumber spectrum of the first-time step is obtained and shown in Figure 10. The maximum value of the spectrum appears at k r = 1.325. We set the average phase velocity to be v p = 1490 m/s and obtain k = 0.5355 m −1 . Then an estimated value of the radial velocity v s = 2.4741 m/s can be obtained according to Equation (7). Following this method, the radial velocity in the 75 min duration of the experiment is obtained as shown in Figure 11, wherein the black solid line is the maximum point. It can be seen that the estimation of the radial velocity is good for most of the time. Around the 10th minute, some points with large deviations appear, caused by the strong interference mentioned above.
The time step is set to be 1s T = , and Equation (3) is calculated using the 128-point FFT. The horizontal wavenumber spectrum of the first-time step is obtained and shown in Figure 10. The maximum value of the spectrum appears at can be obtained according to Equation (7). Following this method, the radial velocity in the 75 min duration of the experiment is obtained as shown in Figure 11, wherein the black solid line is the maximum point. It can be seen that the estimation of the radial velocity is good for most of the time. Around the 10th minute, some points with large deviations appear, caused by the strong interference mentioned above.    The source range is estimated and shown in Figure 13. It can be seen that the estimated source ranges are in good agreement with the ranges from GPS data, and exhibit a bias of 6.23%. The source range is estimated and shown in Figure 13. It can be seen that the estimated source ranges are in good agreement with the ranges from GPS data, and exhibit a bias of 6.23%.
Appl. Sci. 2020, 10, 1005 11 of 13 Figure 13. Comparison of moving source range from GPS data and from estimated radial velocity.

The Broadband Case
In the experiment, the sound source emits signals of multiple frequencies, so the broadband range localization method can also be used to estimate the source range. Three frequencies of 109, 127, and 130 Hz are selected to estimate the source range, and the remaining parameter settings are the same as the narrowband case. According to the calculation of Equation (9), the radial velocity estimation result is obtained and shown in Figure 14. It can be seen from the figure that the interference in the 10th minute has no effect on the radial velocity estimation result, showing the advantage of broadband processing. After obtaining the radial velocity using the broadband method, the source range estimation can be accomplished following the same steps as the narrowband case, which will not be described herein again.

The Broadband Case
In the experiment, the sound source emits signals of multiple frequencies, so the broadband range localization method can also be used to estimate the source range. Three frequencies of 109, 127, and 130 Hz are selected to estimate the source range, and the remaining parameter settings are the same as the narrowband case. According to the calculation of Equation (9), the radial velocity estimation result is obtained and shown in Figure 14. It can be seen from the figure that the interference in the 10th minute has no effect on the radial velocity estimation result, showing the advantage of broadband processing. After obtaining the radial velocity using the broadband method, the source range estimation can be accomplished following the same steps as the narrowband case, which will not be described herein again. the same as the narrowband case. According to the calculation of Equation (9), the radial velocity estimation result is obtained and shown in Figure 14. It can be seen from the figure that the interference in the 10th minute has no effect on the radial velocity estimation result, showing the advantage of broadband processing. After obtaining the radial velocity using the broadband method, the source range estimation can be accomplished following the same steps as the narrowband case, which will not be described herein again.

Conclusions
In this paper, the estimation of the source radial velocity is realized based on the synthetic aperture beamforming method. Then, the least squares criterion is used to estimate the range and time of CPA and the source velocity. Finally, source range localization is accomplished. The proposed method is validated to be effective by both simulation data and SWellEx-96 experimental data. This paper presents a method for estimating the radial velocity in both narrowband and broadband forms. The broadband form is more stable in the presence of interference in the environment or low SNR case. Because the FFT algorithm is used in the calculation of the horizontal

Conclusions
In this paper, the estimation of the source radial velocity is realized based on the synthetic aperture beamforming method. Then, the least squares criterion is used to estimate the range and time of CPA and the source velocity. Finally, source range localization is accomplished. The proposed method is validated to be effective by both simulation data and SWellEx-96 experimental data. This paper presents a method for estimating the radial velocity in both narrowband and broadband forms. The broadband form is more stable in the presence of interference in the environment or low SNR case. Because the FFT algorithm is used in the calculation of the horizontal wavenumber spectrum, the running speed of the method can be increased and a shorter running time can be obtained. Since only the sound speed in water is needed as a priori information, the proposed method is robust to the environmental uncertainty. It should be noted that for the case of no CPA present, after the radial velocity is obtained, the estimation of the source range can still be realized based on the waveguide invariant method. In this case, some broadband signal information is needed, and the specific process can be found in [14].