Multi-resolution estimation of the interference spectrum per pair of modes in the frequency domain Open Access

In a shallow water waveguide, a small number of normal modes can be excited for low frequency and long-range acoustic propagation. Due to the mutual interference between different modes, each pair of normal modes will produce a specific interference spectrum fluctuating with frequency in the way of quasi-cosine function. In this letter, “interference spectrum” means the sound intensity fluctuation caused by a pair of normal modes interfering with each other. Accurate knowledge of the interference spectrum is often required in many practical applications. For example, the fluctuation quasi-period of each interference spectrum is determined by both the wavenumber difference of the corresponding normal mode pair and the source-receiver range. When a single hydrophone and a guide source are available, even if there is only one pair of modes in the received signals, the range of unknown source can also be data-derived estimated according to its fluctuation period.¹ However, in most cases, there are more than two significant normal modes contributing to the sound intensity in a shallow water waveguide, the observed fluctuation in the frequency domain is a superposition of quasi-cosine signals with different periods, and these quasi-cosine signals correspond to the interference spectra for different mode pairs. In this case, the superposition effect will make identification of the quasi-periods difficult, and then the traditional source ranging method mentioned in Ref. 1 will be limited or will fail. Hence it is critically important to separate the interference spectrums from the sound intensity fluctuation for further applications, such as source ranging.

The interference spectrum can be expressed as a quasi-cosine function based on normal mode theory, which is the foundation of interference spectrum separation. However, both the amplitude and phase of the interference spectrum strongly depend on the numbers of considerable mode-pairs and environment parameters. In addition, the interference spectrum will show non-stationary oscillation sometimes. In addition, only one-dimensional (1-D) data with a finite bandwidth can be used for interference spectrum separation. Actually, interference spectrum separation is similar to extraction of horizontal wavenumbers in ocean waveguides, which has been thoroughly studied.^2–8 The intensity fluctuation along the frequency is the superposition of different interference spectra just like the sound field fluctuation along the range is the superposition of different normal modes. A classical method for horizontal wavenumber extraction is Fourier or Hankel transform,^2,3 which has low resolution and is not suitable for interference spectrum separation. In contrast to Hankel transform, the resolution of the Prony method is very high.⁴ However, the Prony method requires the signal to be in the form of an exponential signal, and interference spectra generally cannot meet this condition. Recently, compressed sensing (CS) and Bayesian learning (BL) methods have been used to extract horizontal wavenumbers.^5,6 However, both CS and BL cannot be used in interference spectrum separation directly for the follow reasons. Because of the nonstationarity of the interference spectrum, it is difficult to express the intensity fluctuation with a set of fixed basis functions sparsely. In addition, it is difficult to give the prior probability of an interference spectrum.

In this letter, a multi-resolution estimation method (MREM) of the interference spectrum per pairs of modes is presented and verified by simulations for the low-frequency narrow-band signal. The idea of MREM is as follows: (i) a $2×L$ Hankel matrix is constructed using the received sound intensity data of length L + 1 in the frequency domain; (ii) then the recursion reconstruction and singular value decomposition (SVD) is combined to process this kind of data matrix and separate the interference spectra with different fluctuation periods. The theoretical derivation of MREM can be found in the supplementary material.⁹ The main advantage and quality of MREM include the unbiased estimation of interference spectra, only requiring a single receiver and without any prior environmental information.

When considering low-frequency acoustic propagation at sufficiently long range in shallow water, the received acoustic intensity at depth z_r created by an omni-directional point source with circular frequency ω at depth z_s and range r can be expressed as

where

$P(ω)$ is the power spectral density of the source, $N(ω)$ is the number of propagation modes, $kn(ω)$ and $ϕn(z)$ are the horizontal wavenumber and the mode depth function for the nth mode, and $Δkmn(ω)=km(ω)−kn(ω)$ with m < n. The first term in Eq. (1) is the non-interference intensity, which varies slowly with range and frequency. $AnAm* cos[Δkmn(ω)r]$ corresponds to the interference spectrum of $m−n$ mode pair, and the second term in Eq. (1) is a summation of many interference spectra for different mode pairs, which gives rise to the stripe-pattern in the $r−ω$ plane.^10,11 Combining Eqs. (1) and (2), it is easy to find that the amplitude of interference spectrum is a weakly depend on ω and the phase is depend on horizontal wavenumber difference $Δkmn$ and range r. If more than two propagation modes exist in the ocean waveguide, the observed broadband acoustic intensity by a single hydrophone is usually the superposition of non-interference intensity and these interference spectra of different mode pairs. It is difficult to observe interference spectrum of a mode pair directly from the intensity fluctuation along frequency. And the main goal of this letter is to separate interference spectra of different mode pairs, which will be useful in source ranging.¹ Because only the summation term that contains interference spectra will be considered in this letter, Eq. (1) can be rewritten as

where

and it has been assumed that $⟨AnAm* cos[Δkmn(ω)r]⟩≈0$⁠.

To facilitate the computer processing, we discretize the one-dimensional broadband acoustic intensity recorded by a hydrophone in the frequency domain with frequency interval of $Δω$⁠. Then $I=[I0,I1,…,IL]$⁠, where $In=I(r,zr,zs,ω0+nΔω)$ with r, z_r, and z_s are fixed. In order to realize interference spectra separation, a $2×L$ Hankel matrix H is reconstructed by I, i.e.,

According to SVD,

where $u1,2∈C2×1$ represent the left singular vectors, $v1,2∈CL×1$ represent the right singular vectors, “C” represents complex number, $σa≥σd≥0$ are singular values, $Ha=σau1v1H$⁠, which describes the outline of I, $Hd=σdu2v2H$⁠, which describes the details of I and “H” represents conjugate transpose. I can be expressed as the summation of an outline vector A and a detail vector D,

where

according to Eqs. (4) and (5), $Ha,d(m,n)$ is the m-th row and n-th column element of $Ha,d$⁠. Replacing I by A, and repeating the process from Eq. (4) to Eq. (7)J times, one can get

where $AJ$ is the outline vector after iteration J times. Assuming that $Δkmn(ω)≈amnω+b$ when $ω∈[ω0,ω0+LΔω]$⁠, the interference spectrum for m-n mode pair can be approximated as a cosine signal with an angular frequency of $amnr$⁠. In the supplementary material, we have shown that compared with I, the amplitude of the signal with the angular frequency of $amnr$ in $AJ$ is approximately reduced by $cos2J(amnrΔω/2)$ times,⁹ which means that only the interference spectrum with $amnr$ being the minimum is significant in $AJ$ after a certain number J iterations. In addition, the phase of each interference spectrum is unbiased in each iteration according to the supplementary material.⁹ How to choose J will be discussed in Sec. 3. Assume that different interference spectra are orthogonal to each other, i.e.,

where $amm′nn′$ is a constant and δ_mn is the Kronecker δ function. Then the interference spectrum with $amnr$ being the minimum can be extracted as follows:

where $ÃJ=AJ/||AJ||2$ and “T” means transpose. Generally, I is the superposition of cosine type interference spectra with different periods and the interference spectrum of the biggest period in I can be separated by the above process, i.e., Eqs. (4)–(10). Replacing I by $I−I1$ and repeating the process from Eq. (4) to Eq. (10), one can separate the interference spectrum of the lowest frequency in $I−I1$⁠, which is also the interference spectrum with the second lowest frequency in I. By repeating this process, the interference spectrum of different frequencies can be separated one by one. The details of the decomposition of I is shown in Table 1, where $H(1,α)$ is the 1st row and α-th column element of H and $H(2,α)$ is similar.

A specific numerical simulation serves to illustrate the performances of MREM for interference spectrum separation presented in Sec. 2. Figure 1 shows the simulated environmental model with a thermocline from depth 15–30 m overlying a homogeneous basement. The source and receiver are located at 40 and 50 m below the sea surface, respectively. The horizontal range between the source and the receiver is 21 km. The synthetic acoustic data is computed by kraken¹² with

where l is a nature number and $Δω=0.5π rad/s$⁠. Gaussian white noise is added to the acoustic data to simulate the noise in the real data and the signal to noise ratio (SNR) is 20 dB. Figure 2(a) shows the simulated intensity, which fluctuates with frequency. Because the fluctuation of intensity with frequency is the result of superposition of multiple interference spectra, there is no obvious periodic oscillation in Fig. 2(a). One can use the program in Table 1 to separate different interference spectra once the parameters N and J are chosen. Generally, N is consistent with the estimated number of interference spectra. In addition, it is easy to verify that $σa/σd$ increases with the number of SVD iterations J, and when only one interference spectrum dominates, $σa/σd$ is approximate constant.⁹ So, $σa/σd$ can be selected as an indicator which is used to determine the value of J. Taking the separation of the interference spectrum corresponding to the mode 1–2 pair as an example, Fig. 2(b) shows $σa/σd$ as a function of J. One can find that $σa/σd$ is almost a constant when $J≥180$ and J = 180 (marked by an arrow) in the separation of the interference spectrum corresponding to the mode 1–2 pair. Figures 2(c) and 2(d) show the estimated interference spectra for the mode 1–2 pair and 1–3 pair, respectively. For comparison, the corresponding truth interference spectra are also shown in Figs. 2(c) and 2(d), which are consistent with the estimated results. Note that the interference spectrum for the mode 1–3 pair is a non-stationary signal where the amplitude of the interference spectrum changes with ω slowly, and the SVD based separation method still works. The other interference spectra are submerged by noise and cannot be separated. Assuming that the interference spectrum is a stable cosine signal, i.e., the amplitude does not change with ω, the oscillation frequency t_m of the mth interference spectrum can be estimated by Fourier transform (FT) or multiple signal classification (MUSIC), where m is a positive integer. When t_m is known, the following equation can be used to estimate the mth interference spectrum:

where I_pm represents the mth estimated interference spectrum, ω₁ and ω₂ are the lower and upper bounds of acoustic signal frequency, respectively,

and

Figures 2(c) and 2(d) also show the estimated interference spectra based on FT and MUSIC. One finds that the oscillation frequency and initial phase of the interference spectra estimated by a and B are different than the real situation from Figs. 2(c) and 2(d). The main reason for this phenomenon is that the amplitude of the interference spectra is not constant but changes slowly with ω.

In order to further verify the effectiveness of interference spectrum separation method, a second simulation will be presented in the letter. In the simulation below, all parameters are the same with the simulation above except that

where l is a nature number and $Δω=0.5π rad/s$⁠. Figure 3(a) shows the simulated intensity, which fluctuates with frequency again. Figures 3(b)–3(d) show the estimated interference spectra for the mode 1–2 pair, 1–3 pair, and 1–6 pair, respectively. For comparison, the corresponding truth interference spectra are also shown in Figs. 3(b) to 3(d), which are consistent with the estimated results again. The above two simulations verify the effectiveness of MREM for interference spectrum separation. For comparison, the estimated interference spectra based on FT and MUSIC are also shown in Figs. 3(b) to 3(d). Compared with the real situation, the oscillation frequency and initial phase of the interference spectra estimated by FT and MUSIC are quite different, which is consistent with the results of the first simulation.

Finally, it should be stressed that the performance of MREM degrades with the decrease of SNR. When SNR decreases to 5 dB, the interference spectra for the mode 1–3 pair in Fig. 2(d) and for the mode 1–6 pair in Fig. 3(d) cannot be separated anymore. When SNR decreases to 0 dB, the interference spectra for the mode 1–2 pair in Fig. 2(c) and for the mode 1–3 pair in Fig. 3(c) cannot be separated anymore. When SNR decreases to −3 dB, the interference spectrum for the mode 1–2 pair in Fig. 3(a) cannot be separated anymore.

In a broadband acoustic intensity, different interference spectra for each pair of modes are generally superimposed over one another. Thus, the detailed features of the interference spectra are often unable to be observed due to this superimposing effect. In this letter, we investigated the corresponding inverse problem, and proposed the MREM for interference spectrum separation. First, a $2×L$ Hankel matrix is constructed by a broadband intensity data of length L + 1. Second, the interference spectrum per pairs of modes is estimated by the recursion reconstruction and SVD of this data matrix. Finally, simulation examples are carried out to illustrate the validity of the proposed method. The limitation of the proposed method is the moderate requirement of signal-to-noise ratio. Future research, which include experimental tests, a more effective method for complex cases with higher frequency or deeper sea-depth, and a method to determine the mode numbers of estimated interference spectrum, will be studied in the near future.

This work was supported by the Taishan Scholars under Grant No. tsqn201909053 and Fundamental Research Funds for the Central Universities under Grant Nos. 202065005 and 862001013102. The authors thank Ning Wang for useful suggestions and discussions. W.G. and X.L. contributed equally to this work.