Broadband Sound Intensity Interference Frequency Periodicity and Pulse Source Localization

Abstract

In order to analyze the frequency periodicity characteristics of acoustic field interference and realize acoustic source ranging (ASR), the normal mode model is used to analyze the interference characteristics of the broadband acoustic field under the condition of horizontally layered medium; the broadband received signal field when the broadband pulse signal passes through the acoustic field is also simulated. The variation of interference patterns with frequency is analyzed, and their spatial interference characteristics and mechanisms are analyzed. Based on the interference theory, the relation between the acoustic source range and the frequency periodicity of the broadband acoustic intensity interference is derived. Simulation and experimental results show that this relation can accurately estimate the far-field acoustic source range, and the estimation accuracy and real-time performance are greatly improved compared with previous methods. Besides, simulation shows that the method combined with multiple-receiver ranging obtains high-precision direction of arrival (DOA) estimation as well as ASR. The relation between acoustic source position and broadband acoustic field interference frequency periodicity can be used to improve far-field ASR and DOA estimation, which is of great value for oceanography, marine engineering, and marine military. In addition, this relation can also be extended to that between the modal interference frequency periodicity and other related parameters in other physical fields for parameter inversion.

1. Introduction

Matched field processing (MFP) is a classic underwater acoustic source localization (ASL) method, and a series of improved methods [1,2,3,4,5,6] based on this method are often used for underwater ASL. However, this method requires consistency between marine environment and the model, but marine environment is time variant and spatial variant, leading to a degeneration in the ASL accuracy of MFP, or even complete failure. A series of methods have been developed to solve the problem of model mismatch. Richardson [7] used the maximum a posterior estimation and realized narrow-band ASL in the deep-sea uncertain sound velocity field, the result of which was more stable than that of MFP. Krolik [8] used the minimum variance beamforming to realize ASL in a random ocean channel, which had high stability to the random change of the sound speed profile (SSP) between the acoustic source and the receiver. Mourad [9] used reference sound waves to correct marine environment variation, realized ASL, and eliminated distance integral effect of marine environment. Fialkowski [10] used Bartlett processing to eliminate noise, and used focusing to reduce the impact of marine environment uncertainty, realizing ASL in noisy and uncertain ocean. In the presence of internal waves, Yoo [11] adopted uncoupled normal modes to achieve ASL, while Book [12] used the hybrid Cramer-Rao lower bound and significantly improved the performance of low-frequency broadband coherent ASL. Song [13] adopted maximum likelihood estimation to realize ASL in a random medium, with an ASL relative error of 20%. Yoo [14] adopted broadband matched beam processing and realized shallow water broadband ASL in the presence of internal waves. Collison [15] used the normalized matched mode processing for ASL and used the replica field as a priori estimation to improve adaptation to environmental mismatch. Soares [16] found that MFP for ASL would become robust to temporal and spatial changes of environment, if the parameters related to mismatch were clearly identified, appropriately modeled and estimated by focusing methods. Finette [17] used horizontal array beamforming to deal with anisotropic internal wave fields and weakened the horizontal refraction effect by internal waves.

Zhang [18] adopted a seabed horizontal line array (HLA) MFP to realize shallow water ASL, the performance of which was better than that of HLA at other depths. Based on a reliable acoustic path time-reverse ASL method, Wang [19] used an HLA to realize deep-sea ASL, and enhanced the focusing effect by increasing the number and spacing of the array elements. For broadband ASL in shallow water, based on the waveguide invariant, Qi [20] proposed an operator according to the autocorrelation function of the received signal, effectively realizing passive ASL. By extracting the peak points of the autocorrelation function of the experimentally received signal, Thomassin [21] obtained the corresponding time delay of multipath arrivals, realizing ASL. Based on autocorrelation function of the broadband signal received by a single hydrophone, Zhang [22] realized shallow water ASL adapting to marine environment mismatch.

In addition, the waveguide invariant method is also a commonly used broadband ASL method. Based on the relation between the slope of broadband sound field interference striation and the acoustic source range, this method realizes underwater acoustic source ranging (ASR). Yu [23] achieved passive ranging of moving targets based on waveguide invariants. Based on vector sound field interference structure, Lin [24] used waveguide invariant to carry out passive ASR. Wang [25] used HLA to estimate the target velocity, using shallow ocean waveguide invariant for ASR. Combined with the virtual source method, Liu [26] proposed a waveguide invariant ASR method, improving the source ranging accuracy. Under the condition of linearly changed ocean depth with distance, based on waveguide invariant and warping transform, Shao [27] used an iterative algorithm to realize ASR, which reduced the ranging error.

However, in the broadband sound field interference pattern, as an important piece of information, the quasi-periodic variation of the sound intensity interference with frequency is often ignored. In our previous research [28], the quasi-periodicity in broadband sound intensity fluctuation with frequency was discovered, and the quantitative relation between the quasi-periodicity and the position of the soliton internal wave was derived, realizing soliton internal wave localization. This paper focuses on the phenomenon that the broadband sound intensity fluctuates quasi-periodically with frequency. The quantitative relation between the quasi-periodicity in broadband sound intensity fluctuation with frequency and the acoustic source range is found, and the underwater acoustic source is localized. The accuracy of the results obtained is significantly better than that of MFP. Moreover, this work is remarkable in that the results coming from the study of interference frequency periodicity suggest that the frequency periodicity can be used for high-accuracy underwater ASR even with a single receiving range. Further, this new scientific thinking model will be applied to other fields involving frequency periodicity of physical field modal interference, bringing more practical application results. For instance, underwater broadband vector acoustic field ranging and acoustic source ranging in the air can be realized in the same way, which will be studied next.

Based on the principle of frequency synthesis, we adopted the normal mode method and obtained the relevant results of the broadband pulse signal propagating in shallow water. The theoretical derivation leads us to the relation between the acoustic source range and the interference period of the two-normal-mode interference curves in frequency. This relation can accurately estimate the acoustic source range when the receiver is far from the acoustic source (further than 10 km), and is simulated and verified by experimental data processing. Moreover, compared with the traditional “waveguide invariant-interference striation slope” ranging method, the frequency periodicity method has better real-time performance and lower computational complexity for the acoustic source range monitoring. The frequency periodicity method can not only accurately estimate the acoustic source range, but also can be used as a scientific idea to be extended to the application of other physical fields for normal mode broadband interference parameter measurement.

2. Materials and Methods

Simulation and Processing of Broadband Received Signals

The selected SSP is shown in Figure 1a. The sea depth is 60 m, the liquid seabed density is 1.7 g/cm³, the seabed longitudinal wave sound velocity is 1650 m/s, and the longitudinal wave attenuation coefficient is 0.8 dB/wavelength. The acoustic source depth is z_s = 50 m, and the receiver depth z_r = 51 m (the purpose is to ensure that the received signal intensity of the first two normal modes is relatively large), and the reference sound velocity C₀ = 1500 m/s.

The broadband pulse source signal shown in Figure 2a is $s_{in}\left(t\right)=0.5\sin \left({\omega }_{c}t\right)\left[1-\cos \left(0.2{\omega }_{c}t\right)\right]$, in which ${\omega }_c=2\pi f_c$,$f_c=350\mathrm{Hz}$ is the central frequency for broadband source signal, and $t\in \left[0,5/f_c\right]$. Performing Fourier transform on the signal, we obtained the broadband source spectrum $S_{in}\left(\omega \right)$ as shown in Figure 1b.

According to the principle of frequency synthesis, the normal mode method is used to simulate the received pulse sequence of the broadband source signal. According to the normal mode method, the broadband acoustic pressure field is as follows:

$$P\left(\omega ,R\right)={\displaystyle \sum _{m=1}^M{A}_m\left(\omega ,R\right)\frac{e^{i{k}_m\left(\omega \right)R}}{\sqrt{k_m\left(\omega \right)R}}}$$

(1)

where $A_m\left(\omega ,R\right)$ and $k_m\left(\omega \right)$ are the amplitude and eigenvalue of the m-th normal mode, respectively, and $\omega$ is the angular frequency of the sound wave. For bottom with attenuation, $k_m\left(\omega \right)$ can be written as $k_m\left(\omega \right)=K_m\left(\omega \right)+i{\alpha }_m\left(\omega \right)$, in which ${\alpha }_m\left(\omega \right)$ is the attenuation coefficient. Then, broadband acoustic intensity is as follows:

$$\begin{array}{ll}I\left(\omega ,R\right)& ={\displaystyle \sum _{m=1}^M{\left|A_m\left(\omega ,R\right)\frac{e^{\left[i{K}_m\left(\omega \right)-{\alpha }_m\left(\omega \right)\right]R}}{\sqrt{k_m\left(\omega \right)R}}\right|}^2}\\ & +{\displaystyle \sum _{\begin{array}{l}m=1\\ m\ne n\end{array}}^M{\displaystyle \sum _{n=1}^M\frac{A_m\left(\omega ,R\right)A_n^{\ast }\left(\omega ,R\right)e^{-\left[{\alpha }_m\left(\omega \right)+{\alpha }_n\left(\omega \right)\right]R}}{\sqrt{k_m\left(\omega \right)k_n^{\ast }\left(\omega \right)}R}\exp \left[i{K}_{mn}\left(\omega \right)R\right]}}\end{array}$$

(2)

where $K_{mn}\left(\omega \right)=K_m\left(\omega \right)-K_n\left(\omega \right)$.

Equation (2) is the broadband sound intensity, where the interference between normal modes can be seen from $\exp \left[i{K}_{mn}\left(\omega \right)R\right]$. In the waveguide invariant method, by the partial derivative of the broadband sound intensity with respect to distance and frequency [29], the relation between the slope of broadband sound intensity interference striation and the acoustic source position is obtained:

$$\frac{\partial R}{\partial \omega }=\frac{\partial I\left(\omega ,R\right)/\partial \omega }{\partial I\left(\omega ,R\right)/\partial R}=-\frac{r}{\omega }\frac{\mathrm{d}{S}_g}{\mathrm{d}{S}_p}$$

(3)

Thus, by measuring the slope of broadband sound intensity interference striation, the broadband underwater ASR is realized.

The waveguide invariant method can obtain accurate results for broadband ASR. However, the method relies on the slope of interference striation, which requires a change in R. For a fixed source, it is inconvenient to obtain a broadband sound intensity interference pattern through moving receivers and realize ASR.

However, it is noticed that the broadband interference in Equation (2) is the interference between normal modes, and the interference term is $\exp \left[i{K}_{mn}\left(\omega \right)R\right]$. For such an interference term, we make the following analysis:

When R is fixed, $K_{mn}\left(\omega \right)$ changes with $\omega$. Generally, for low-order normal modes, $K_{mn}\left(\omega \right)$ becomes smaller as $\omega$ increases. Additionally, when $K_{mn}\left(\omega \right)R=2n\pi$ ($n\in {\mathrm{N}}^{+}$), $\exp \left[i{K}_{mn}\left(\omega \right)R\right]=1$, and the interference between the m-th and n-th normal modes takes a peak. We also assume that the frequencies corresponding to two adjacent interference peaks are ${\omega }_1$ and ${\omega }_2$. Then,

$$\begin{array}{ll}\{\begin{array}{l}{K}_{mn}\left({\omega }_1\right)R=2n\pi \\ K_{mn}\left({\omega }_2\right)R=2\left(n\pm 1\right)\pi \end{array}& \\ \Rightarrow \left|K_{mn}\left({\omega }_1\right)-K_{mn}\left({\omega }_2\right)\right|R=2\pi & \\ \Rightarrow R=2\pi /\left|K_{mn}\left({\omega }_1\right)-K_{mn}\left({\omega }_2\right)\right|& ,(n\in {\mathrm{N}}^{+})\end{array}$$

(4)

Thus, R can be obtained through Equation (4). Compared with the waveguide invariant method, the advantage of the proposed method is that R can be obtained by the broadband sound intensity at a single distance. To use the Equation (4) for underwater broadband ASR, we need to filter the received sound signal and only retain the m-th and n-th modes. Generally, the amplitudes of the received first two modes are the strongest, so we keep the first two modes to be conducive to use Equation (4) for ASR.

Next, we take the broadband pulse acoustic source as an example, and consider the actual received signal field, to test the effectiveness of the above ASL algorithm, Equation (4). Based on the received signal spectrum

$$S_{out}\left(\omega ,R\right)=S_{in}\left(\omega \right)P\left(\omega ,R\right)$$

(5)

and the corresponding received broadband interference spectrum

$$\begin{array}{ll}{S}_{out,I}\left(\omega ,R\right)& ={\left|S_{out}\left(\omega ,R\right)\right|}^2\\ & ={\left|S_{in}\left(\omega \right)\right|}^2\times \{{\displaystyle \sum _{m=1}^M{\left|A_m\left(\omega ,R\right)\frac{e^{\left[i{K}_m\left(\omega \right)-{\alpha }_m\left(\omega \right)\right]R}}{\sqrt{k_m\left(\omega \right)R}}\right|}^2}\\ & +{\displaystyle \sum _{\begin{array}{l}m=1\\ m\ne n\end{array}}^M{\displaystyle \sum _{n=1}^M\frac{A_m\left(\omega ,R\right)A_n^{\ast }\left(\omega ,R\right)e^{-\left[{\alpha }_m\left(\omega \right)+{\alpha }_n\left(\omega \right)\right]R}}{\sqrt{k_m\left(\omega \right)k_n^{\ast }\left(\omega \right)}R}\exp \left[i{K}_{mn}\left(\omega \right)R\right]}}\}\end{array}$$

(6)

as well as inverse Fourier transform, the received pulse sequence is obtained:

$$s_{out}\left(t,R\right)={\displaystyle {\int }_{{\omega }_1}^{{\omega }_2}{S}_{out}\left(\omega ,R\right)e^{-i\omega t}d\omega }$$

(7)

Figure 2b–e show the received pulse sequence at different distances. As shown in Figure 2a–b, the pulse with an initial length of 0.02 s diverges strongly, and the final length becomes 0.2 s. Different normal mode separation methods, such as mode filtering [30] and warping transform [31], can handle these received pulses sequence. Next, we intercepted the first two largest wave packets in Figure 2c–e, starting from 0 s and intercepted to 0.62, 0.42 and 0.22 s, respectively. Figure 3 shows the intercepted received temporal signals.

3. Broadband Interference Simulation and Theoretical Analysis

After intercepting the first two largest wave packets, we obtained a new received temporal signal, written as $s_{out,m+n}\left(t,R\right)$ (for the first two modes, m = 1, n = 2). Next, $s_{out,m+n}\left(t,R\right)$ is processed with Fourier transform, to obtain frequency domain received signal as follows:

$$S_{out,m+n}\left(\omega ,R\right)={\displaystyle {\int }_{-\infty }^{\infty }{s}_{out,m+n}\left(t,R\right)e^{i\omega t}dt}$$

(8)

and the corresponding received broadband interference spectrum is

$$\begin{array}{ll}{S}_{out,m+n,I}\left(\omega ,R\right)& ={\left|S_{out,m+n}\left(\omega ,R\right)\right|}^2\\ & ={\left|S_{in}\left(\omega \right)\right|}^2\times \{{\left|A_m\left(\omega ,R\right)\frac{e^{\left[i{K}_m\left(\omega \right)-{\alpha }_m\left(\omega \right)\right]R}}{\sqrt{k_m\left(\omega \right)R}}\right|}^2+{\left|A_n\left(\omega ,R\right)\frac{e^{\left[i{K}_n\left(\omega \right)-{\alpha }_n\left(\omega \right)\right]R}}{\sqrt{k_n\left(\omega \right)R}}\right|}^2\\ & +\frac{A_m\left(\omega ,R\right)A_n^{\ast }\left(\omega ,R\right)e^{-\left({\alpha }_m\left(\omega \right)+{\alpha }_n\left(\omega \right)\right)R}}{\sqrt{k_m\left(\omega \right)k_n^{\ast }\left(\omega \right)}R}\exp \left[i{K}_{mn}\left(\omega \right)R\right]\}\end{array}$$

(9)

Combined with the broadband acoustic source frequency spectrum $S_{in}\left(\omega \right)$ (Figure 1b), only the 250–450 Hz band with a large amplitude is reserved, to obtain the frequency response of the broadband sound field

$$P_{m+n}\left(\omega ,R\right)=S_{out,m+n}\left(\omega ,R\right)/S_{in}\left(\omega \right)$$

(10)

as well as the broadband sound intensity containing only two normal modes:

$$\begin{array}{ll}{I}_{m+n}\left(\omega ,R\right)& ={\left|P_{m+n}\left(\omega ,R\right)\right|}^2\\ & =\{{\left|A_m\left(\omega ,R\right)\frac{e^{\left[i{K}_m\left(\omega \right)-{\alpha }_m\left(\omega \right)\right]R}}{\sqrt{k_m\left(\omega \right)R}}\right|}^2+{\left|A_n\left(\omega ,R\right)\frac{e^{\left[i{K}_n\left(\omega \right)-{\alpha }_n\left(\omega \right)\right]R}}{\sqrt{k_n\left(\omega \right)R}}\right|}^2\\ & +\frac{A_m\left(\omega ,R\right)A_n^{\ast }\left(\omega ,R\right)e^{-\left[{\alpha }_m\left(\omega \right)+{\alpha }_n\left(\omega \right)\right]R}}{\sqrt{k_m\left(\omega \right)k_n^{\ast }\left(\omega \right)}R}\exp \left[i{K}_{mn}\left(\omega \right)R\right]\}\end{array}$$

(11)

Frequency domain received signals for various R are shown in Figure 4a–c, and the interference pattern $I\left(\omega ,R\right)$ of the broadband acoustic pressure field $P\left(\omega ,R\right)$ containing all normal modes is shown in Figure 4e. The 250–450 Hz band is chosen because the source amplitude of the band is much larger and more flat than that of other bands, resulting in high ASL accuracy. This is a limitation of the proposed algorithm: the received signal spectrum is required to be relatively flat throughout the band. If not, you need to choose a relatively flat band from the broadband source spectrum, such as the 250:1:450 Hz band in Figure 1b.

As shown in Figure 4e, the interference pattern of the sound field obtained by superposing all modes tends to be stable as R increases, and the sound field interference at a far field is dominated by the interference between the first two modes (interference spatial period d_mn is close to $d_{12}\left(\omega \right)=2\pi /K_{12}\left(\omega \right)$). As shown in Figure 4a–c, the frequency response of the broadband sound field fluctuates with frequency quasi-periodically, and since only the first two modes are retained, the interference in Figure 4a–c is dominated by the spatial interference of the two modes. In other words, the interference pattern $I\left(\omega ,R\right)$ is actually dominated by the interference pattern $I_{1+2}\left(\omega ,R\right)$ of the first two modes, and, according to Equation (11),

$$\begin{array}{ll}{I}_{1+2}\left(\omega ,R\right)& =\{{\left|A_1\left(\omega ,R\right)\frac{e^{\left[i{K}_1\left(\omega \right)-{\alpha }_1\left(\omega \right)\right]R}}{\sqrt{k_1\left(\omega \right)R}}\right|}^2+{\left|A_2\left(\omega ,R\right)\frac{e^{\left[i{K}_2\left(\omega \right)-{\alpha }_2\left(\omega \right)\right]R}}{\sqrt{k_2\left(\omega \right)R}}\right|}^2\\ & +\frac{A_1\left(\omega ,R\right)A_2^{\ast }\left(\omega ,R\right)e^{-\left({\alpha }_1\left(\omega \right)+{\alpha }_2\left(\omega \right)\right)R}}{\sqrt{k_1\left(\omega \right)k_2^{\ast }\left(\omega \right)}R}\exp \left[i{K}_{12}\left(\omega \right)R\right]\}\end{array}$$

(12)

Taking source spectrum into consideration, the received interference spectrum

$$\begin{array}{ll}{S}_{out,1+2,I}\left(\omega ,R\right)={\left|S_{in}\left(\omega \right)\right|}^2& \times \{{\left|A_1\left(\omega ,R\right)\frac{e^{\left[i{K}_1\left(\omega \right)-{\alpha }_1\left(\omega \right)\right]R}}{\sqrt{k_1\left(\omega \right)R}}\right|}^2+{\left|A_2\left(\omega ,R\right)\frac{e^{\left[i{K}_2\left(\omega \right)-{\alpha }_2\left(\omega \right)\right]R}}{\sqrt{k_2\left(\omega \right)R}}\right|}^2\\ & +\frac{A_1\left(\omega ,R\right)A_2^{\ast }\left(\omega ,R\right)e^{-\left[{\alpha }_1\left(\omega \right)+{\alpha }_2\left(\omega \right)\right]R}}{\sqrt{k_1\left(\omega \right)k_2^{\ast }\left(\omega \right)}R}\exp \left[i{K}_{12}\left(\omega \right)R\right]\}\end{array}$$

(13)

As shown by Equation (12), when R is constant, $\exp \left[i{K}_{12}\left(\omega \right)R\right]$ term causes the fluctuation of sound intensity in frequency. Combining the fluctuation property with Equation (4), if the frequencies of two adjacent peak points for the sound intensity fluctuation are ${\omega }_A$ and ${\omega }_B$, the relation between R and ${\omega }_A$ and ${\omega }_B$ is as follows:

$$R=2\pi /\left|K_{12}\left({\omega }_A\right)-K_{12}\left({\omega }_B\right)\right|$$

(14)

According to Equation (14), the broadband pulse ASR can be done. If the adjacent valley values ${\omega }_C$ and ${\omega }_D$ of the sound intensity interference in frequency are taken, $R=2\pi /\left|K_{12}\left({\omega }_C\right)-K_{12}\left({\omega }_D\right)\right|$ still holds, and the derivation is not repeated here. So, $K_{12}\left(\omega \right)$ changes with frequency, resulting in a quasi-periodical change in frequency for $S_{out,1+2,I}\left(\omega ,R\right)$, which can be used to obtain R.

In addition, the farther R is, the smaller $\left|K_{12}\left({\omega }_A\right)-K_{12}\left({\omega }_B\right)\right|$ required for frequency interference to complete a full cycle is, and the smaller corresponding $\left|{\omega }_A-{\omega }_B\right|$ is, which is easier to satisfy in a narrower frequency band, so as to achieve ASR more easily. On the contrary, when R is too small, it is not easy to meet the above condition, causing ranging failure. Therefore, the performance of the proposed method on far-field ranging is theoretically better than that of near-field. The poor accuracy of near-field ranging is the second drawback of the proposed method.

However, in practice, for other ASL methods, near-field ASL is often relatively simple, while far-field precise ASL is more difficult, so this algorithm can solve the far-field precise ASL problem that is difficult for other methods, and near-field ASL can be solved by other methods.

In addition, to reduce the ranging error, we used the method below. For example, ${\omega }_{A1}$ and ${\omega }_{B1}$ are frequency points corresponding to the two peaks of the broadband sound field interferogram. The acoustic intensity fluctuates with frequency, and the band between ${\omega }_{A1}$ and ${\omega }_{B1}$ contains n interference cycles ($n\in {\mathrm{N}}^{+}$). Then,

$$\begin{array}{l}\left|K_{12}\left({\omega }_{A1}\right)-K_{12}\left({\omega }_{B1}\right)\right|R=2n\pi \\ R=2n\pi /\left|K_{12}\left({\omega }_{A1}\right)-K_{12}\left({\omega }_{B1}\right)\right|\end{array}$$

(15)

Therefore, R can be accurately obtained. Taking Figure 4a–c, for example, according to multiple interference periods in the broadband interferogram of the sound field at three distances, the estimated R is shown in

Table 1. Relative error for acoustic source ranging (ASR) at three ranges.

Real Range R/km	Frequency Periods n	Estimated Range $\widehat{\mathit{R}}$/km	Relative Error/%
20	5	19.814	–1.18
30	7	30.004	0.01
40	10	40.036	0.09

.:

As shown in Table 1, ranging errors at the three distances are all less than 2%, and the ASR errors at 30 and 40 km are even less than 0.1%. Therefore, the method realizes high accuracy in ASR for the far field. Next, considering the limitation of the proposed method in the near field, we used the method to analyze broadband sound intensity interference $I\left(\omega ,R\right)$ and draw the estimated source range when the actual source range is R = 30:0.1:100 km, as shown in Figure 5. Considering that high-order normal mode attenuation is large in the far field, far-field sound intensity interference is dominated by low-order normal modes, so there is no need to truncate the received temporal signal and keep the first two wave packets only. To obtain the broadband receiving spectrum, we conducted Fourier transform on the entire received temporal signal; then, combining broadband receiving spectrum with Equation (15), we obtained the ASR result shown in Figure 5. In addition, as R increases, the interference period in frequency will become shorter, and the wrong peak points will be found due to insufficient sampling frequency points, bringing ranging errors. Therefore, we avoided the interference peaks at lower frequencies, and, starting from the third interference peak point counting from low to high frequency, we counted the number of interference peaks and the corresponding number of interference cycles n and calculated R.

As shown in Figure 5a, at all R, the relative ranging error is far less than 10%. To further show the high precision of the ASR algorithm, we reduced the error limitation to 1%, and the relative ranging error with R is shown in Figure 5b, where the error is less than 1% at most distances and between 1 and 2% at a small part of distances, and higher than 2% and lower than 10% for only one measured point. Moreover, when R is further than 40 km, the error of all distances is within 1%, illustrating the high accuracy of the proposed method on far-field ASR.

It is often impossible for us to predict the frequency spectrum of the acoustic source, so we can only perform ASR based on the received signal. Under such a premise, we still used the above method, select the 250:1:450 Hz band, and use the fluctuation periodicity of the received intensity spectrum $S_{out,1+2,I}\left(\omega ,R\right)$ in frequency to perform ASR, and the ranging results are shown in Figure 6.

As shown in Figure 6b, when $S_{out}\left(\omega ,R\right)$ is used for ASR, the relative ranging error is within 2% at most distances, and 2–5% for a small number of distance points, which also achieves accurate far-field ASR.

Next, we determined the acoustic source azimuth based on the multi-receiver ranging. Theoretically, only three non-collinear receivers are required to locate an acoustic source, and ASR relative to the three receivers can be calculated based on the geometric relation. To test the source azimuth measurement performance of the proposed method, we kept the marine environment parameters consistent with those in Figure 1a, and choose the layout shown in Figure 7a, where A, B, and C are three receivers and D is the acoustic source, and their relative positions are shown in Figure 7a. The distances from D to A, B, and C are, respectively,

$$R_{AD}\approx 78.102\mathrm{km},R_{BD}\approx 50.990\mathrm{km},R_{CD}\approx 102.956\mathrm{km},$$

Combining $S_{out}\left(\omega ,R\right)$, the estimated range $\widehat{r}$ and relative error $E_r$ obtained by the frequency periodicity method are as follows:

$$\{\begin{array}{l}\widehat{R_{AD}}=79.447\mathrm{km},\widehat{R_{BD}}=51.752\mathrm{km},\widehat{R_{CD}}=104.620\mathrm{km}\\ E_{r,AD}=1.72\%,E_{r,BD}=1.49\%,E_{r,CD}=1.62\%\end{array}$$

As demonstrated by the the above results, the ASR errors at three receivers are all within 2%, achieving high-precision ranging. Next, according to the geometric relation, we drew three circles with the points A, B, and C as the center and $\widehat{R_{AD}}$,$\widehat{R_{BD}}$,$\widehat{R_{CD}}$ as the radius, respectively, to determine the acoustic source position, as shown in Figure 7b.

As shown in Figure 7b, the three circles centered on A, B, and C have a common intersection point close to the acoustic source, D. The details of the intersection point in Figure 7b are enlarged and shown in Figure 7c, where there are actually three intersections for the three circles near D, namely (50.29,51.74) km, (50.42,51.67) km, and (50.44,51.56) km. If the average of the horizontal and vertical co-ordinates of the three points, (50.383,51.657) km, is used as the positioning result, selecting A(–10,0) km as the reference point and positive x axis as the 0° direction, the actual acoustic source position by polar co-ordinate is (78.102 km, 39.806°), the estimated acoustic source position is (79.464 km, 40.547°), the relative ranging error is 1.74%, and the error of direction estimation is 0.741°. Thus, the distance and direction of the acoustic source have been determined accurately. Here, we still emphasize that, due to theoretical limitation, the proposed method is more accurate for far-field ASL but has insufficient performance for near-field ASL. Moreover, the received signal spectrum is required to be relatively flat.

4. Experimental Verification

In order to test the underwater ASL performance of the proposed method in the actual marine environment, the data from the Qingdao sea test in the summer of 2007 by the Institute of Acoustics of the Chinese Academy of Sciences is processed using the above method, and the results are as follows.

The SSP measured in the experiment is shown in Figure 8a. The sea depth is H = 36 m, and the 32-element vertical array was used to receive data. Among the elements, three were faulty and 29 were effective. The effective receiving depth z_r = 2.5:1:30.5 m, the acoustic source depth z_s = 16.5 m, and 38 explosion acoustic source tests were performed in the experiment. The 38 distances from acoustic sources to the receiver are shown in Figure 8b. The bottom geometry parameters were not measured in the experiment. According to the results in Reference [32], supposing the bottom is half infinite space and liquid seabed, the bottom sound speed c_b = 1588 m/s, bottom density ${\rho }_b=1.85\mathrm{g}/{\mathrm{cm}}^3$, and bottom longitudinal wave attenuation coefficient ${\alpha }_b$ is fitted base on the results in Reference [31]: ${\alpha }_b=5.1525\times {10}^{-7}{f}^2+3.2252\times {10}^{-4}f+0.04347$, where the unit of frequency f is Hz, the unit of ${\alpha }_b$ is $\mathrm{dB}/\lambda$, and $\lambda$ is the wavelength of the sound wave. Combined with the normal mode method, the corresponding eigenvalue and eigenfunction are calculated, and the needed sound wave frequency band is determined according to the actual calculation requirements. The calculated $K_{12}\left(\omega \right)$ combined with the test environment is shown in Figure 8c.

In the experiment, an amount of 0.1 kg TNT was used as the explosion acoustic source, and the Fourier transform was performed on the received temporal sound signal to obtain the received signal spectrum:

$$S_{out}\left(z_r;\omega \right)={\displaystyle {\int }_{-\infty }^{\infty }{s}_{out}\left(z_r;t\right)e^{i\omega t}dt}$$

(16)

The normal mode method is used to obtain the eigenvalue and eigenfunction, and then, combined with the sound field interference frequency periodicity method, the acoustic source range is estimated. According to the KRAKEN program [33] and the environment in Figure 8a, only two propagating normal modes are obtained when the sound wave frequency f = 105:0.1:175 Hz. The sound field interference structure is simple, which is convenient for subsequent processing. Therefore, this frequency band is selected for ASR. The received signal spectrum $S_{out}\left(z_r;\omega \right)$ is used to realize ASR. To obtain high-precision results, mode filtering [30] is performed on $S_{out}\left(z_r;\omega \right)$ to obtain received spectrum $S_{out,1}\left(z_r;\omega \right)$ and $S_{out,2}\left(z_r;\omega \right)$ of the first two normal modes, and the filtered received signal spectrum is as follows:

$$S_{out,1+2}\left(z_r;\omega \right)=S_{out,1}\left(z_r;\omega \right)+S_{out,2}\left(z_r;\omega \right)$$

(17)

Taking the 38th bombing location as an example, the actual distance r₃₈ = 31.348 km. Selecting the received acoustic signal at the 27th receiving depth (z_r = 28.5 m, at which the peak of the received spectrum is sharper that that at other depths), and the corresponding $S_{out,1+2}\left(z_r;\omega \right)$ is shown in Figure 9.

As shown in Figure 9, $S_{out,1+2}\left(z_r;\omega \right)$ fluctuates with frequency. To improve the accuracy of ASR, the 128.4–160.3 Hz band is selected, which has a long and stable received signal fluctuation period, and there are eight interference cycles in this band, so the estimated source range $r_{38,est}$ and relative error $E_{r,38}$ are, respectively,

$$\{\begin{array}{l}{r}_{38,est}=8\times 2\pi /\left|\mathrm{Re}\left(k_{12}\left(128.4\mathrm{Hz}\right)-k_{12}\left(160.3\mathrm{Hz}\right)\right)\right|=31.098\mathrm{km}\\ E_{r,38}=\left(r_{38,est}-r_{38}\right)/r_{38}\times 100\%=-0.79\%\end{array}$$

Thus, the frequency periodicity method obtains an accurate ASR result. Next, the method is applied to other bombing locations. To reduce the error, after the mode filtering is used to obtain $S_{out,1+2}\left(z_r;\omega \right)$, the $S_{out,1+2}\left(z_r;\omega \right)$ at each receiving depth is added to obtain the sum of received signal spectrum:

$$S_{out,sum}\left(\omega \right)={\displaystyle \sum _{z_r}{S}_{out,1+2}\left(z_r;\omega \right)}$$

(18)

The ASR is performed according to the interference period of $S_{out,sum}\left(\omega \right)$ in frequency, and the ASR result is shown in Figure 10b; while the ASR result by the MFP algorithm (Equation (1b) in Reference [14]) is shown in Figure 10a.

Comparing Figure 10a,b, it can be seen that the accuracy of the ASR result by the frequency periodicity method is significantly higher than the result by MFP. In the ranging results of MFP, the ranging error of a considerable part of the locations is far more than 10%, while the ranging results of the frequency periodicity method are all within 10%, in which a considerable part are far less than 10%, showing the superiority of the proposed method in ASL accuracy. In addition, compared with MFP, the frequency periodicity method does not need to calculate the entire replica field that changes with distance, but only processes the actual received data, which greatly reduces the amount of calculation for ASR. Moreover, the method only needs the signal at a single receiving distance, and does not rely on the slope of the interference striation of the sound field interference pattern composed of multiple distances to realize ASR. Compared with the waveguide invariant ranging method, the proposed method greatly improves the real-time nature of ranging.

In addition, there are only 36 estimated distance points in Figure 10b, and the estimated distances of explosion locations No.1 and No.2 do not appear. The reason for this is that the actual acoustic source ranges of the two explosion locations are 1.496 and 2.071 km, respectively, which are too close for the proposed method, causing ASR to fail, and they are not marked in the Figure 10b.

5. Conclusions

According to the principle of frequency synthesis, we adopted the normal mode method to simulate and analyze the interferogram of the received signal, when the broadband pulse signal propagates in the sound field. The results show that when there are only two normal modes in the sound field, the sound field interferogram is dominated by the two normal modes, and the periodicity of the sound field interference in frequency is also dominated by the two normal modes, thus obtaining a simple and clear relation between source range and the frequency periodicity of normal mode interference. According to the relation, a method for underwater ASR with high accuracy, good real-time performance and small calculation amount has been obtained. Theoretical simulation and experimental data have verified the above advantages of the proposed method.

The frequency periodicity method not only helps us realize ASR, but can also be extended to the field of frequency periodicity analysis caused by other physical field mode interference, so that real-time monitoring can be done for other types of targets and parameters, such as ASL in air acoustics, and all other fields involving mode interference and ASL analysis.

However, the proposed ranging method is currently only suitable for broadband pulse acoustic sources due to theoretical limitation. The spectrum energy of the received signal should change relatively smoothly with frequency; the near-field ASR accuracy becomes worse, or even completely invalid; the algorithm is only verified with the 2007 experiment in Qingdao, where the medium is nearly horizontally layered, and with a shallow ocean depth of no more than 40 m. In the future, the algorithm will be improved to be suitable for near-field ASR and other types of ASR, including complicated situations, such as the effect of a range-dependent environment and deep ocean.