Harmonic Source Depth Estimation by a Single Hydrophone under Unknown Seabed Geoacoustic Property

Abstract

The passive estimation of harmonic sound source depth is of great significance for underwater target localization and identification. Passive source depth estimation using a single hydrophone with an unknown seabed geoacoustic property is a crucial challenge. To address this issue, a harmonic sound source depth estimation algorithm, seabed independent depth estimation (SIDE) algorithm, is proposed. This algorithm combines the estimated mode depth functions, modal amplitudes, and the sign of each modal to estimate the sound source depth. The performance of the SIDE algorithm is analyzed by simulations. Results show that the SIDE is insensitive to the initial range of the sound source, the source depth, the hydrophone depth, the source velocity, and the type of the seabed. Finally, the effectiveness of the SIDE algorithm is verified by the SWellEX-96 data.

1. Introduction

Acoustic methods for localizing a source, particularly a submerged source, in an underwater acoustic waveguide is a topic of great interest. As an important part of source localization, the depth estimation has particular difficulties but is of great significance to the identification of underwater targets.

Many methods have been proposed for source depth estimation. One of the most popular methods is the matched field processing (MFP) [1,2], which matches the measured field and the replica field for source depth estimation. The MFP can estimate source depth accurately without considering any environmental mismatches. However, when the sound propagation is not exactly modeled, the performance of depth estimation by the MFP degrades greatly, which is referred to as the environmental mismatch problem (EMP). In the MFP, the depth estimation depends on the range estimation, which indicates that an erroneous estimation of the source range would have serious consequences for the depth estimation. Compared with the MFP, matched mode processing (MMP) [3,4,5] can estimate the source range and depth separately and is less sensitive to environmental mismatches. Noting that source depth determines the amplitudes of excited modes, the MMP can determine the source depth by matching the mode amplitudes extracted from real data with the simulated excitation modal amplitudes of the sound source at different depths. However, the MMP is still bothered by the EMP in applications. Moreover, the resolution in the depth direction is poor when only mode amplitudes are used for depth estimation. To extract the mode amplitudes accurately, a full-depth vertical linear array (VLA) [6,7] or a horizontal linear array (HLA) with a large aperture [8,9] is required. Nevertheless, it is difficult to meet the above conditions in many application scenarios. In order to alleviate the effect of the EMP, machine learning (ML) methods have been introduced for underwater source location in recent years [10,11,12,13,14], and they show better performance than the MFP or the MMP in both ranging and depth estimation. Although the ML methods perform well in terms of positioning accuracy and adaptability to the environment in underwater acoustics, they often require a large number of training samples. It is difficult to meet that in many underwater applications. In addition, the interpretability issue of ML has not been resolved, which also limits its promotion in the field of underwater acoustics. An alternative method for overcoming the EMP is MMP based on data-derived modal parameters (mode depth functions and wavenumbers) [15,16]. Nevertheless, a larger aperture is also required to guarantee accuracy. To solve the problem of insufficient aperture, some research works have tried to use the idea of synthetic aperture [17,18] to obtain long-distance spatial sampling via the source motion, but the speed of a moving source is generally unknown. Although the smooth-averaged cross-correlation function of signals recorded by two horizontally separated hydrophones can be used for depth estimation of a source when the speed of the source in unknown [19], this method requires a broadband signal of the sound source and is not suitable for a harmonic source.

Recently, Hongxu Cui et al. proposed a method called the mode depth function constrained modal Doppler velocity estimation (MDFMD-v) method to estimate the velocity of a harmonic source in shallow water [20]. In this method, the harmonic source radiation signal only needs to contain two different frequencies, and no seabed geoacoustic information is required. Combining the MDFMD-v method, this paper introduces a algorithm for harmonic source depth estimation using a single hydrophone. In our algorithm, only the sound speed profile (SSP) of shallow water needs to be known. To address the problem of low depth direction resolution caused by the traditional MMP using only modal amplitudes for depth estimation, this paper proposes a method called seabed independent depth estimation (SIDE). In the SIDE algorithm, a one-dimensional mode sign search algorithm called depth-sign search (DSS) is introduced to obtain the sign information of each mode. Combining the searched modal signs and the extracted modal amplitudes, the SIDE algorithm can improve the depth estimation accuracy effectively. The performance of the SIDE algorithm is analyzed by numerical simulation in this paper, and the simulation results show that the performance of the SIDE algorithm is insensitive to source depth, hydrophone depth, range between source and hydrophone, source velocity, and type of the seabed. When the time domain signal-to-noise ratio (TSNR) is 0 dB, the average relative error of the source depth estimation by the SIDE algorithm is less than 7%, which is far less than that of the conventional MMP with an average relative error of more than 20% under the same condition. Moreover, the effect of TSNR on the SIDE algorithm is analyzed, and the results show that the average relative error of the depth estimation results is less than 12% when the TSNR is greater than −20 dB. Finally, the effectiveness of the SIDE algorithm is validated by the SwellEx-96 data [21].

The innovations of this paper are as follows:

• The seabed independent depth estimation (SIDE) algorithm is proposed, which can realize the source depth estimation without the information of the seabed parameters and the source velocity.
• The depth-sign search (DSS) method is proposed, which can transform the multi-dimensional mode-sign search problem into a one-dimensional search problem, improving the search efficiency significantly.
• The performance of the SIDE algorithm is analyzed by simulation, and the effectiveness of the SIDE algorithm is verified by the Swell-Ex96 data.

The structure of this paper is as follows. Section 2 is the preliminary, where the necessary preparatory knowledge for the SIDE algorithm will be introduced. The theory of the SIDE algorithm is introduced in Section 3, and the DSS is also introduced in this section. In Section 4, the performance of the SIDE algorithm is analyzed by simulations. The effectiveness of the SIDE algorithm is validated by the date in the Swell-Ex96 in Section 5. Finally, conclusions are given in Section 6.

2. Preliminary

In order to ensure the theoretical integrity of the article, before introducing the SIDE algorithm, we first introduce the modal Doppler effect and briefly describe the MDFMD-v, which utilizes the modal Doppler effect to estimate the velocity and frequency of a moving source.

Consider a range-independent waveguide as shown in Figure 1 and assume that a point source is approaching the hydrophone with radial velocity $v_0\left(v_0\ll c\right)$. The source transmits signal $s\left(t,\omega \right)=e^{-i\omega t}$, with $\omega =2\pi f$ being the intrinsic angular frequency of the source. The range between the source and the hydrophone is unknown, and the source depth is assumed to be a constant. According to the normal mode theory, the sound field received by the hydrophone can be expressed as [22]

$$\begin{array}{cc}\hfill p\left(t,\omega \right)& \approx i\sqrt{2\pi }{e}^{-i\frac{\pi }{4}}\sum _{m=1}^M\left(\frac{{\psi }_m\left(\omega ,z_s\right){\psi }_m\left(\omega ,z_r\right)}{\sqrt{k_{rm}\left(\omega \right)R\left(t\right)}}{e}^{\left[i\omega t-i{k}_{rm}\left(\omega \right)R\left(t\right)\right]}\right)\hfill \\ \hfill & \approx i\sqrt{2\pi }{e}^{-i\frac{\pi }{4}}\sum _{m=1}^M\left(\frac{{\psi }_m\left(\omega ,z_s\right){\psi }_m\left(\omega ,z_r\right)}{\sqrt{k_{rm}\left(\omega \right)R\left(t\right)}}{e}^{\left[i{\omega }_{m}t-i{k}_{rm}\left(\omega \right)R_0\right]}\right),\hfill \end{array}$$

(1)

where M denotes the number of propagating modes; $R\left(t\right)=R_0-v_{0}t$ is the range between the source and the hydrophone at time t, with $R_0$ being the initial range; $k_{rm}\left(\omega \right)$ and ${\psi }_m\left(\omega ,z\right)$ are the horizontal wavenumber and the mode depth function of mode m, respectively; and $z_s$ and $z_r$ represent the depth of source and hydrophone, respectively. In Equation (1),

$$\begin{array}{c}\hfill {\omega }_m=\omega +k_{rm}\left(\omega \right)v_0,\hspace{0.17em}m=1,2,\dots ,M,\end{array}$$

(2)

which means that the motion of the sound source causes a frequency shift in the received signal. And the frequency shift of different modes is different, which is the modal Doppler effect.

The MDFMD-v method can estimate the harmonic sound velocity when seabed parameters are unknown. This method combines the modal Doppler and the homomodal similarity (HS) proposed in [20] to realize the estimation of the source velocity. The HS means that the mode depth functions usually vary slowly with frequency, so the mode depth functions excited by two adjacent frequencies are almost the same with high correlation. From Equation (1), one can find that $k_{rm}$ and $v_0$ are coupled. In the MDFMD-v, the HS is utilized to decouple $k_{rm}$ and $v_0$, in which the intrinsic relationship between the horizontal wavenumbers and the mode depth functions are needed. Taking the $m\mathrm{th}$ normal mode as an example, the horizontal wavenumber and the corresponding mode depth function are related by the following differential equation:

$$\begin{array}{c}\hfill \frac{{\partial }^2{\psi }_m\left(\omega ,z\right)}{\partial z^2}+k_{zm}^2\left(\omega ,z\right){\psi }_m\left(\omega ,z\right)=0,\end{array}$$

(3)

where

$$\begin{array}{c}{k}_{zm}\left(\omega ,z\right)=\sqrt{{k_0}^2\left(\omega ,z\right)-k_{rm}^2\left(\omega \right)},\\ k_0\left(\omega ,z\right)=\frac{\omega }{c\left(z\right)},\end{array}$$

and $c\left(z\right)$ is the SSP of the water column. According to reference [7], the above equation can be discretized to the following difference equation:

$$\begin{array}{c}\hfill {\psi }_m\left(\omega ,z+\Delta z\right)+{\psi }_m\left(\omega ,z-\Delta z\right)-2cos\left(\Delta z\sqrt{{\left(\frac{\omega }{c\left(z\right)}\right)}^2-k_{rm}^2\left(\omega \right)}\right){\psi }_m\left(\omega ,z\right)=0,\end{array}$$

(4)

where $\Delta z$ is the difference step. Therefore, ${\psi }_m(\omega ,z)$ can be calculated when $k_{rm}\left(\omega \right)$ and $c\left(z\right)$ are known by the following equation:

$${\psi }_m(\omega ,n\Delta z)=2cos\left(\Delta z\sqrt{{\left(\frac{\omega }{c\left(z\right)}\right)}^2-k_{rm}^2\left(\omega \right)}\right){\psi }_m\left(\omega ,(n-1)\Delta z\right)-{\psi }_m\left(\omega ,(n-2)\Delta z\right),$$

(5)

where $m=1,2,\dots ,M,n=2,3,\dots ,N$, ${\psi }_m(\omega ,0)=0$, ${\psi }_m(\omega ,\Delta z)=1$. Define

$$\begin{array}{c}\hfill {\mathbf{\Psi }}_m\left(\omega ,\mathbf{z}\right)={\left[{\psi }_m\left(\omega ,z_0\right),{\psi }_m\left(\omega ,z_1\right),{\psi }_m\left(\omega ,z_2\right),\cdots ,{\psi }_m\left(\omega ,z_N\right)\right]}^{\mathrm{T}},\end{array}$$

(6)

where $z_n=n\Delta z$ with $n=0,1,\dots ,N$, $N\Delta z\approx H$, H is the water depth, $\mathbf{z}={\left[z_0,z_1,\cdots ,z_N\right]}^{\mathrm{T}}$, and ‘T’ represents the matrix transpose. Then, the HS can be judged by calculating the correlation coefficient of the same mode depth function at two different frequencies, which can be approximated by

$$\begin{array}{c}\hfill C{C}_m\left({\omega }_1,{\omega }_2\right)=\frac{{\mathbf{\Psi }}_m^{\mathrm{H}}\left({\omega }_1,\mathbf{z}\right){\mathbf{\Psi }}_m\left({\omega }_2,\mathbf{z}\right)}{{∥{\mathbf{\Psi }}_m\left({\omega }_1,\mathbf{z}\right)∥}_2{∥{\mathbf{\Psi }}_m\left({\omega }_2,\mathbf{z}\right)∥}_2},\end{array}$$

(7)

where ‘H’ represents the conjugate transpose.

Assume that a moving source emits signals at two adjacent frequencies ${\omega }_{10}$ and ${\omega }_{20}$. In passive detection, both the velocity and the frequencies of the moving source are unknown. Due to the modal Doppler effect, the frequencies of the received signal consist of ${\mathit{\omega }}_1={\left[{\omega }_{11},{\omega }_{12},\cdots ,{\omega }_{1M_1}\right]}^{\mathrm{T}}$ and ${\mathit{\omega }}_2={\left[{\omega }_{21},{\omega }_{22},\cdots ,{\omega }_{2M_2}\right]}^{\mathrm{T}}$, where ${\omega }_{1m_1}\approx {\omega }_{10}+k_{r{m}_1}\left({\omega }_{10}\right)v_0$ and ${\omega }_{2m_2}\approx {\omega }_{20}+k_{r{m}_2}\left({\omega }_{20}\right){\nu }_0$ with $m_1=1,2,\dots M_1$ and $m_2=1,2,\dots ,M_2$. Although it is difficult to accurately obtain the frequencies of the source, an approximate range of the two frequencies can be determined by spectral analysis. From the time–frequency analysis, one can determine the source frequencies search intervals ${\omega }_1\in \left[{\omega }_{\mathrm{lmin}},{\omega }_{\mathrm{lmax}}\right],{\omega }_2\in \left[{\omega }_{2\min },{\omega }_{2\max }\right].$ According to experience, one can determine the velocity search interval $\nu \in \left[{\nu }_{\min },{\nu }_{\max }\right].$ Then, the mode depth functions corresponding to each ${\omega }_1,{\omega }_2$ and $\nu$ are calculated according to Equations (1) and (5). Let

$${\mathbf{\Psi }}_{m_1}\left({\omega }_1,\mathbf{z}\right)={[{\psi }_{m_1}\left({\omega }_1,z_0\right),{\psi }_{m_1}\left({\omega }_1,z_1\right),\cdots ,{\psi }_{m_1}\left({\omega }_1,z_N\right)]}^{\mathrm{T}}$$

and

$${\mathbf{\Psi }}_{m_2}\left({\omega }_2,\mathbf{z}\right)={[{\psi }_{m_2}\left({\omega }_2,z_0\right),{\psi }_{m_2}\left({\omega }_2,z_1\right),\cdots ,{\psi }_{m_2}\left({\omega }_2,z_N\right)]}^{\mathrm{T}}$$

represent the extracted mode depth functions at ${\omega }_1$ and ${\omega }_2$ respectively, where $m_1,m_2\in {\mathbb{N}}^{+}$. Based on the fact that not all modes satisfy the HS when $\left\{{\omega }_1,{\omega }_2,v\right\}\ne \left\{{\omega }_{10},{\omega }_{20},v_0\right\}$, the loss function for the moving source velocity estimation can be constructed as:

$$\begin{array}{c}\hfill Loss\left({\omega }_1,{\omega }_2,v\right)=-\sum _{m_1=m_2}\frac{{\mathbf{\Psi }}_{m_1}^{\mathrm{H}}\left({\omega }_1,\mathbf{z}\right){\mathbf{\Psi }}_{m_2}\left({\omega }_2,\mathbf{z}\right)}{{∥{\mathbf{\Psi }}_{m_1}\left({\omega }_1,\mathbf{z}\right)∥}_2{∥{\mathbf{\Psi }}_{m_2}\left({\omega }_2,\mathbf{z}\right)∥}_2}.\end{array}$$

(8)

When the loss function is minimal, the correlation between ${\mathbf{\Psi }}_{m_1}({\omega }_1,\mathbf{z})$ and ${\mathbf{\Psi }}_{m_2}({\omega }_2,\mathbf{z})$ is the highest:

$$\begin{array}{c}\hfill \left\{{\widehat{\omega }}_1,{\widehat{\omega }}_2,\widehat{v}\right\}=\underset{{\omega }_1,{\omega }_2,v}{\arg \min } Loss\left({\omega }_1,{\omega }_2,v\right).\end{array}$$

(9)

3. Theory of the SIDE Algorithm

In Section 2, we introduce the MDFMD-v method, which offers an innovative approach to estimating the velocity of a moving harmonic source without seabed information. Combining Equations (2) and (5), we can also calculate the horizontal wavenumbers and the corresponding mode depth functions.

After obtaining the horizontal wavenumbers and the mode depth functions, the MMP can be used to estimate the source range and the source depth at the same time, theoretically. However, considering the fact that target ranging errors will lead to target depth estimation errors, it is not wise to use the MMP method to estimate the target distance at the same time while only the target depth information is focused. Therefore, T. C. Yang proposed a method to estimate the depth independently [4]. This method utilizes only the normal mode amplitudes, and defines the source depth ambiguity function as

$$\begin{array}{c}\hfill D_a(z,\omega )=\left|\sum _{m=1}^M\left|{\mathbf{\Psi }}_m\left(\omega ,z\right)\right|\left|\frac{A_m}{{\psi }_m\left(\omega ,z_r\right)}\right|\right|,\end{array}$$

(10)

where $\left|A_m\right|$ represents the amplitude of the $m\mathrm{th}$ normal mode,

$$\begin{array}{c}\hfill A_m\propto \frac{1}{\sqrt{k_{rm}\left(\omega \right)}}{\psi }_m\left(\omega ,z_s\right){\psi }_m\left(\omega ,z_r\right).\end{array}$$

(11)

And the depth corresponding to the maximum value of $D_a\left(z,\omega \right)$ is the estimation of the source depth. When the velocity of a moving source is known, $\left|A_m\right|$ can be estimated by the Hankel transform (HT) [23], and the source depth can be estimated by Equation (10). Due to only the mode amplitudes being used in Equation (10), the depth estimation method proposed by T. C. Yang is sensitive to the estimation results of the amplitude of each mode, which shows a high sidelobe. In addition, when the estimation results of the amplitude of each mode are biased due to seabed attenuation and other factors, the sidelobe peak value may be higher than that of the main lobe. This will lead to the wrong depth estimation.

Different from the method proposed by T. C. Yang, not only the mode amplitudes but also the mode signs are used for source depth estimation in the SIDE algorithm. Actually, the source depth information is modulated simultaneously on the amplitude and sign of each normal mode. Moreover, the sign of each mode is coupled with the propagation phase of each mode, which is difficult to obtain when the distance of the sound source is unknown. A naive method for obtaining modal signs is to try all possible combinations of modal signs by ergodic search. However, the number of all possible combinations of mode signs increases exponentially with the number of modes M. And the efficiency of the ergodic search method is low when the number of modes in the ocean waveguide is large. In addition, in order to obtain the real sign of each mode, it is also necessary to set the sign discrimination standard. Based on the current research, neither the efficient search method of modal signs nor the correctness judgment method of modal signs have been studied. The DSS method, a fast search method for modal sign, which is a component of the SIDE algorithm, is introduced first. And then we give the correctness judgment criteria of the searches’ modal signs.

From Equation (11), it is easy to find that the signs of $A_m$ are only influenced by ${\psi }_m(\omega ,z_s)$ and ${\psi }_m(\omega ,z_r)$, where $m=1,2,\dots ,M$. Generally, the hydrophone depth $z_r$ is a known parameter. When the mode depth functions have been extracted, ${\psi }_m(\omega ,z_r)$ are known, and the signs of $A_m$ are only determined by $z_s$, which means that an ergodic search for all meaningful combinations of modal signs can be achieved by searching only for a single parameter $z_s$. So, this modal sign search method is named depth-sign search (DSS).

To judge the correctness of the modal signs, we combine the DSS with the source depth estimation process. Firstly, we construct the sign scalar

$$\begin{array}{c}\hfill {\delta }_m\left(\omega ,q\right)=\mathrm{sign}\left[{\psi }_m\left(\omega ,z_r\right){\psi }_m\left(\omega ,z_q\right)\right],m=1,\cdots ,M,q=1,\cdots ,N,\end{array}$$

(12)

where $z_q=q\Delta z$, and 

$$\mathrm{sign}\left[x\right]=\left\{\begin{array}{cc}1\hfill & x\ge 0\hfill \\ -1\hfill & x<0\hfill \end{array}\right..$$

${\delta }_m(\omega ,q)$ can be considered an estimation of the mth modal sign. Based on Equation (10), we define a new source depth ambiguity function by considering ${\delta }_m\left(\omega ,q\right)$:

$$\begin{array}{c}\hfill D(z,q,\omega )=\left|\sum _{m=1}^M{\mathbf{\Psi }}_m\left(\omega ,z\right)\frac{\left|A_m\right|{\delta }_m(\omega ,q)}{{\psi }_m\left(\omega ,z_r\right)}\right|,\end{array}$$

(13)

According to reference [24], $D(z,q,\omega )$ closely approximates the Dirichlet kernel function when ${\delta }_m\left(\omega ,q\right)$ is the correct modal signs. Leveraging this property, a loss function can be devised to measure the discrepancy between $D(z,q,\omega )$ and that of the Dirichlet kernel. The KL divergence is chosen here to measure the difference between $D(z,q,\omega )$ and the Dirichlet kernel:

$$\begin{array}{c}\hfill KL(q,\omega )={\int }_0^{H}D(z,q,\omega )ln\frac{D(z,q,\omega )}{\left|Ds(z,q)\right|}dz,q=1,\cdots ,N,\end{array}$$

(14)

where the integral operation can be approximated by discrete summation,

$$Ds(z,q)=\frac{sin\left[\left(M+1\right)\pi \left(z-z_q\right)/H\right]}{sin\left[\pi \left(z-z_q\right)/\left(2H\right)\right]}+\epsilon$$

is the Dirichlet kernel. $\epsilon$ is a small quantity to avoid $Ds(z,q)$ being zero. In this paper, we take $\epsilon ={10}^{-6}$. Then, the correct modal signs can be determined by

$$\begin{array}{c}\hfill {\delta }_m\left(\omega ,q_0\right)=\mathrm{sign}\left[{\psi }_m\left(\omega ,z_r\right){\psi }_m\left(\omega ,z_{q_0}\right)\right],m=1,\cdots ,M,\end{array}$$

(15)

where $q_0=\underset{q}{\arg \min } KL(q,\omega ),q=1,\cdots ,Q$. When both the modal amplitudes and modal signs are obtained, a joint source depth estimation function can be defined as Equation (16) by considering at least two different frequencies which are needed in the MDFMD-v:

$$\begin{array}{c}\hfill {\widehat{z}}_{joint}=\underset{z}{\arg \max }\left\{\frac{D\left(\mathbf{z},q_1,{\omega }_1\right)}{{∥D\left(\mathbf{z},q_1,{\omega }_1\right)∥}_2^2}+\frac{D\left(\mathbf{z},q_2,{\omega }_2\right)}{{∥D\left(\mathbf{z},q_2,{\omega }_2\right)∥}_2^2}\right\},\end{array}$$

(16)

where $q_1=\underset{q}{\arg \min } KL\left(q,{\omega }_1\right),q_2=\underset{q}{\arg \min } KL\left(q,{\omega }_2\right)$. Since the aforementioned process of sound source depth estimation does not incorporate any seabed information, we refer to the method of sound source depth estimation as the seabed independent depth estimation method (SIDE) algorithm. The pseudocode implementation of the SIDE algorithm is presented in Algorithm 1. In the following sections, the performance of the SIDE algorithm will be analyzed through simulations in Section 4 and subsequently validated with experimental data presented in Section 5.

Algorithm 1 SIDE

Input:  time domain pressure $p\left(t\right)$, source frequencies search intervals ${\omega }_1\in \left[{\omega }_{1\min },{\omega }_{1\max }\right]$, ${\omega }_2\in \left[{\omega }_{2\min },{\omega }_{2\max }\right]$, velocity search interval $v\in \left[v_{\min },v_{\max }\right]$, hydrophone depth $z_r$ Output:  source depth $z_r$  1: compute modal Doppler frequencies ${\omega }_1$, ${\omega }_2$ and $A_m$ using the Fourier transform of $p\left(t\right)$: ${\omega }_1={\left[{\omega }_{11},{\omega }_{12},\cdots ,{\omega }_{1M_1}\right]}^{\mathrm{T}}{\omega }_2={\left[{\omega }_{21},{\omega }_{22},\cdots ,{\omega }_{2M_2}\right]}^{\mathrm{T}}$  2: compute mode depth functions ${\mathbf{\Psi }}_{m2}\left({\omega }_1,\mathbf{z}\right)$ and ${\mathbf{\Psi }}_{m2}\left({\omega }_2,\mathbf{z}\right)$ of each $\left\{{\omega }_1,{\omega }_2,v\right\}$ using Equation (5)  3: compute source frequencies and velocity: $\left\{{\widehat{\omega }}_1,{\widehat{\omega }}_2,\widehat{v}\right\}=\underset{{\omega }_1,{\omega }_2,v}{\arg \min }\left\{-{\displaystyle \sum _{m_1=m_2}}\frac{{\mathbf{\Psi }}_{m_1}^{\mathrm{H}}\left({\omega }_1,\mathbf{z}\right){\mathbf{\Psi }}_{m_2}\left({\omega }_2,\mathbf{z}\right)}{{∥{\mathbf{\Psi }}_{m_1}\left({\omega }_1,\mathbf{z}\right)∥}_2{∥{\mathbf{\Psi }}_{m_2}\left({\omega }_2,\mathbf{z}\right)∥}_2}\right\}$  4: compute horizontal wavenumber: $k_{r{m}_1}\left({\omega }_1\right)=\frac{{\omega }_{1m_1}-{\widehat{\omega }}_{10}}{\widehat{v}}$, $k_{r{m}_2}\left({\omega }_2\right)=\frac{{\omega }_{2m_2}-{\widehat{\omega }}_{20}}{\widehat{v}}$  5: construct the sign scalar: ${\delta }_m\left(\omega ,q\right)=\mathrm{sign}\left[{\psi }_m\left(\omega ,z_r\right){\psi }_m\left(\omega ,z_q\right)\right]$  6: compute the minimum value of the loss function: $q_0=\underset{q}{\arg \min }{\int }_0^{H}D(z,q,\omega )ln\frac{D(z,q,\omega )}{\left|Ds(z,q)\right|}dz$  7: compute source depth: ${\widehat{z}}_s=\underset{z}{\arg \max }\left\{\left|{\displaystyle \sum _{m=1}^M}{\mathbf{\Psi }}_m\left(\omega ,z\right)\frac{\left|A_m\right|{\delta }_m(\omega ,q_0)}{{\psi }_m\left(\omega ,z_r\right)}\right|,\right\}$

4. Simulation

In this section, the performance of the SIDE algorithm will be analyzed by simulations. In our analysis, we focus on the effects of hydrophone depth, source depth, source velocity, initial distance between source and hydrophone, signal-to-noise ratio, and seabed type on the SIDE algorithm.

Firstly, consider the waveguide as depicted in Figure 1. The SSP is selected from the SWellEX-96 experiment, illustrated in Figure 2b, with specific seabed information denoted as $c_b=1750\mathrm{m}/\mathrm{s}$, ${\rho }_b=1.9\mathrm{g}/{\mathrm{cm}}^3$, ${\alpha }_b=0.15\mathrm{dB}/\lambda$, and the depth of water $H=200\mathrm{m}$. The source, whose depth is $z_s=54\mathrm{m}$, approaches the hydrophone from position $R_0=10\mathrm{km}$ at a constant velocity $v_0=2.5\mathrm{m}/\mathrm{s}$, and emits a signal $s\left(t\right)=cos\left(2\pi f_{1}t\right)+cos\left(2\pi f_{2}t\right)$, where $f_1=170\mathrm{Hz}$, $f_2=180\mathrm{Hz}$. The time domain waveform (within a 1 s observation window) is shown in Figure 3a. And Figure 3b is the normalized spectrum of the emitted signal. The hydrophone is positioned at a depth of $z_r=120\mathrm{m}$, and the initial distance between the source and the hydrophone is 10 km. The received signal with a 1 s observation window is shown in Figure 3c. In order to simulate the influence of environmental noise, white Gaussian noise has been added to the received signal, and the TSNR is 0 dB. Figure 3d shows the normalized spectrum of the received signal. In order to observe the modal Doppler, the length of the time window of the Fourier transform (FT) is set as 2000 s to ensure adequate frequency resolution. From Figure 3d, one can clearly see the line spectrum splitting phenomenon caused by the mode Doppler effect. In addition, the amplitude of each modal Doppler frequency is approximated as the amplitude of each mode. It is assumed that the source is a non-cooperative source with unknown velocity and frequency during the data processing.

The velocity and frequencies of the source are estimated using the MDFMD-v method, and the results are shown in Figure 4. Note that the result of Equation (8) is a three-dimensional matrix corresponding to $f_1$, $f_2$ and v, which will make it difficult to visualize the estimation results. To facilitate the presentation of the results, Figure 4 shows the results for the $f_1-v$ and the $f_2-v$ planes at the location of the minimum loss. $f_1$ searches from 160 Hz to 180 Hz at 0.1 Hz intervals, $f_2$ searches from 170 Hz to 190 Hz at 0.1 Hz intervals, and v searches from 0 m/s to 8 m/s at 0.1 m/s intervals. In Figure 4, the white circles represent the locations of the real source velocity and frequencies. And the frequencies and velocity estimation results are $f_1=170\mathrm{Hz}$, $f_2=180\mathrm{Hz}$ Hz and $v_0=2.5\mathrm{m}/\mathrm{s}$, which are well matched with the real.

Upon determining the source velocity and frequencies, the horizontal wavenumbers and the corresponding mode depth functions can be computed by Equations (2) and (5), and the amplitude of each mode can be estimated using the FT. Figure 5 illustrates the estimated horizontal wavenumbers and mode depth functions for the frequency $f_1=170\mathrm{Hz}$, and all subsequent results pertain to the frequency $f_1$ unless explicitly specified. In Figure 5a, the red circle and blue cross symbols represent the estimated and the real horizontal wavenumbers, respectively. Similarly, in Figure 5b, the blue line and the red dashed lines represent the estimated and the real mode depth functions, respectively. Almost all of the correlation coefficients $C{C}_m$ are larger than 0.99, except for the fourth mode with a correlation coefficient of 0.66. The poor estimation result for the lowest-order mode is mainly due to the fact that the source travel distance is smaller, resulting in inaccurate horizontal wavenumber estimation of the lowest normal mode. Inaccurate estimation of the horizontal wavenumber will lead to inaccurate estimation of the mode depth function.

Figure 6a illustrates the normalized amplitudes of modal Doppler frequencies. Here, the normalized amplitudes mean that the largest amplitude is one. As a reference, normalized modal amplitudes calculated by Equation (11) for different modes are also given in Figure 6a. Figure 6b shows the loss function defined by Equation (14) evolution during the optimization process utilizing DSS. It can be found that the value of the loss function is the smallest near the true depth of the sound source.

Figure 6c exhibits the depth ambiguity function computed by the SIDE algorithm under various scenarios: individual estimation for $f_1=170\mathrm{Hz}$, $f_2=180\mathrm{Hz}$, and the joint estimation. For comparison, the result estimated by traditional MMP is also given in Figure 6c. In the cases of individual estimation for $f_1$, $f_2$, and joint estimation, the sign combinations are chosen based on the minimization of the loss function. In Figure 6c, we can find that the results estimated by the SIDE algorithm, $z_s=54\mathrm{m}$, have very low sidelobes for both the single line spectrum and joint estimation. In contrast, the conventional MMP method, despite the presence of a peak at the depth of the sound source, has high sidelobes, and the highest peak occurs at $z_s=198\mathrm{m}$, with a large estimation error. Notice that the loss function minimum in Figure 6b is located at the source depth, but this is not necessarily the case in practice. As in Figure 6d, the signs of the mode depth functions in the red boxes are the same.

To evaluate the performance of the SIDE algorithm, we analyze the source depth estimation performance of the SIDE algorithm under different initial ranges of the source, source depth $z_s$, hydrophone depth $z_r$, and source velocity $v_0$. The initial range of the source is varied from 6 km to 15 km in increments of 0.2 km, the source velocity varied from 1 m/s to 8 m/s in increments of 0.1 m/s, and both source depth and hydrophone depth are varied from 1 m to 200 m with increments of 1 m. When we analyze the effect of one of the parameters on the SIDE algorithm, the values of the other parameters are consistent with the initial simulation. The simulation results presented in Figure 7 are the root mean square error (RMSE) obtained from 400 independent Monte Carlo experiments, and the RMSE is defined as

$$\begin{array}{c}\hfill \mathrm{RMSE}=\sqrt{\frac{1}{N_{mon}}\sum _{i=1}^{N_{mon}}{\left({\widehat{z}}_i-z_s\right)}^2}\end{array}$$

(17)

where $N_{mon}$ is the number of independent Monte Carlo experiments, ${\widehat{z}}_i$ is the depth estimation result for the ith independent Monte Carlo experiment, and $z_s$ is the real source depth.

In Figure 7, the blue line represents the RMSE of frequency $f_1$, the black line represents the RMSE of frequency $f_2$, and the red line represents the RMSE of the joint estimate. In addition, the blue, black, and red regions are the range of standard deviations of the results for $f_1$, $f_2$, and joint estimation, respectively. From Figure 7, it is evident that the SIDE algorithm achieves low RMSE in depth estimation when $\mathrm{SNR}=0\mathrm{dB}$.

The results in Figure 7 suggest that the initial range of the source, the source velocity, the hydrophone depth, and the source depth have little impact on the performance of the SIDE algorithm, and the average relative error $Ea$ in the source depth estimation is less than 7% when $\mathrm{SNR}=0\mathrm{dB}$. The average relative error is defined by Equation (18). For comparison, the traditional MMP was conducted to estimate the source depth under the same conditions. Results are shown in Figure 8, where the standard deviation of the estimation results is not given as in Figure 7 because the traditional MMP estimation results are very unstable. Compared with Figure 7, it is easy to find that the accuracy of the sound source depth determination by the traditional MMP is not high, and the average relative error $Ea$ is more than 20%. This indicates that the traditional MMP is sensitive to the source velocity, the hydrophone depth, and the source depth:

$$\begin{array}{c}\hfill Ea=\frac{1}{z_s{N}_{mon}}\sum _{i=1}^{N_{mon}}\left|{\widehat{z}}_i-z_s\right|\times 100\%.\end{array}$$

(18)

For a more comprehensive analysis of the performance of the SIDE algorithm, it is essential to examine the impact of SNR on the SIDE algorithm. Figure 9 illustrates the source depth estimation results across different SNR and compares them with the traditional MMP, with the time-domain SNR varying from −40 dB to 0 dB in increments of 2.5 dB.

In Figure 9, the blue, black, and red lines show the results estimated using the SIDE algorithm, indicating the RMSE for frequencies $f_1$, $f_2$, and the joint estimation, respectively. Similarly, the cyan, pink, and green lines are the RMSEs for the conventional MMP. Figure 9 shows that the SIDE algorithm outperforms the traditional method. The SIDE algorithm has a RMSE of less than 10 m at $\mathrm{SNR}>-20\mathrm{dB}$ for both the single line spectrum and the joint estimation. Additionally, joint estimation exhibits lower errors compared to single line spectrum estimation.

To demonstrate the impact of the seafbed on the SIDE algorithm, we tested the source depth estimation performance of the SIDE algorithm under different seabed parameter conditions. Consider a multi-layered seabed with sedimentary layer, the environment in this simulation are shown in Figure 10. The SSP is kept the same as the previous one. And the seabed consists of two layers and an elastic half-space. Layer 1 can be considered the sedimentary layer. And the parameters of Layer 1 have six cases as shown in

Table 1. The estimation results of source depth with different Layer 1.

Type	${\mathit{\rho }}_1$ ($\mathbf{g}/{\mathbf{cm}}^3$)	${\mathit{c}}_{\mathit{p}1}$ ($\mathbf{m}/\mathbf{s}$)	${\mathit{\alpha }}_{\mathit{p}1}$ (dB/$\mathit{\lambda }$)	${\mathit{c}}_{\mathit{s}1}$ ($\mathbf{m}/\mathbf{s}$)	${\mathit{\alpha }}_{\mathit{s}1}$ (dB/$\mathit{\lambda }$)	${\mathit{z}}_{\mathit{s}}$ ($\mathbf{m}$)	${\widehat{\mathit{z}}}_{\mathit{s}}$ ($\mathbf{m}$)	$\mathit{\sigma }$ ($\mathbf{m}$)
sand	1.98	1742	0.4	382	12.1	54	54.304	0.098
slit	1.83	1677	0.5	467	12.2	54	54.512	0.050
sandy silt	1.56	1552	0.5	379	14.1	54	53.417	0.166
sand-mud-clay	1.58	1578	0.3	409	14.0	54	53.707	0.042
clayey silt	1.43	1535	0.2	312	14.5	54	53.763	0.229
silty clay	1.42	1519	0.2	287	15.0	54	53.192	0.404

. The seabed type can be changed by changing the parameters of Layer 1. For Layer 2, which has fixed parameters, its densities ${\rho }_2=2.06\mathrm{g}/{\mathrm{cm}}^3$, compression wave velocity $c_{p2}=1700\mathrm{m}/\mathrm{s}$, attenuation of compression wave ${\alpha }_{p2}=0.06\mathrm{dB}/\lambda$, shear wave velocity $c_{s2}=460\mathrm{m}/\mathrm{s}$, shear wave attenuation ${\alpha }_{s2}=10\mathrm{dB}/\lambda$, and thickness is 40 m. The elastic half-space density ${\rho }_h=2.54\mathrm{g}/{\mathrm{cm}}^3$, compression wave velocity $c_{ph}=3500\mathrm{m}/\mathrm{s}$, attenuation of compression wave ${\alpha }_{ph}=0.03\mathrm{dB}/\lambda$, shear wave velocity $c_{sh}=1750\mathrm{m}/\mathrm{s}$, and shear wave attenuation ${\alpha }_{sh}=3\mathrm{dB}/\lambda$. During the test, SNR is fixed at 0 dB. In addition, the results in Table 1 are obtained by averaging the results of 400 independent Monte Carlo experiments, and $z_s$ is true source depth, and ${\widehat{z}}_s$ and $\sigma$ are the mean and standard deviation of the 400 Monte Carlo experiments, respectively.

Table 1 indicates that the seabed type has little impact on the SIDE algorithm. To further discuss the effect of the ocean waveguide on the SIDE algorithm, we vary the SSP and water depth of the waveguide. Figure 11 illustrates a typical negative gradient (NG) SSP waveguide (Figure 11a) and a Pekris (Figure 11b) waveguide with multi-layered sedbed. Layer 1 is slit, and its parameters are shown in Table 1. The sound speed of the Pekris waveguide is 1500 m/s. In addition, the water depths H are set to be 50 m, 100 m, and 200 m, respectively. Similarly, 0 dB Gaussian noise is added, and 400 Monte Carlo experiments are performed. The source depth, receiver depth, and results are shown in

Table 2. The estimation results of source depth with different SSP and water depth.

	$\mathbf{H}=200\mathbf{m}$				$\mathbf{H}=100\mathbf{m}$				$\mathbf{H}=50\mathbf{m}$
SSP	${\mathit{z}}_{\mathit{s}}$ ($\mathbf{m}$)	${\mathit{z}}_{\mathit{r}}$ ($\mathbf{m}$)	${\widehat{\mathit{z}}}_{\mathit{s}}$ ($\mathbf{m}$)	$\mathit{\sigma }$ ($\mathbf{m}$)	${\mathit{z}}_{\mathit{s}}$ ($\mathbf{m}$)	${\mathit{z}}_{\mathit{r}}$ ($\mathbf{m}$)	${\widehat{\mathit{z}}}_{\mathit{s}}$ ($\mathbf{m}$)	$\mathit{\sigma }$ ($\mathbf{m}$)	${\mathit{z}}_{\mathit{s}}$ ($\mathbf{m}$)	${\mathit{z}}_{\mathit{r}}$ ($\mathbf{m}$)	${\widehat{\mathit{z}}}_{\mathit{s}}$ ($\mathbf{m}$)	$\mathit{\sigma }$ ($\mathbf{m}$)
NG	54	120	54.200	0.201	54	80	54.395	0.206	9	40	8.978	0.987
Pekris	54	120	56.486	0.439	54	80	54.380	0.173	9	40	8.701	0.929

. The rest of the simulation conditions are the same as before.

From Table 2, we can find that the variations in the SSP and depth of the waveguide have little impact on the SIDE algorithm, and the relative errors in the depth estimation are all less than 7%.

5. Experiment

In this section, the efficacy of the SIDE algorithm is validated using data from the SWellEX-96 experiment, which includes the S5 and S59 events. Due to S5 events having more stable Doppler phenomena than S59 events, this section will use S5 events for method validation. The movement path of the two sources for the S5 event and the locations of the hydrophone arrays are illustrated in Figure 2a. The deep source is 54 m underwater, and the shallow source is 9 m underwater. The SSP is depicted in Figure 2b. The sources started their track south of all of the arrays and proceeded northward at a speed of 5 knots, and transmitted numerous tonals of various source levels between 49 Hz and 400 Hz. The VLA recorded the full 75 min event.

This section will perform method validation using 30 min of data, 5 min after the start of event S5, recorded by one hydrophone, which is the one closest to the seabed in the VLA (at a depth of 212.25 m). Along the trajectory, water depths predominantly range between 200 m and 220 m, with minimal variation in water depth, allowing it to be treated as a homogeneous waveguide. In Figure 2c, the range between the source and the vertical array is depicted, with the blue dotted line representing the data segment being analyzed. The shallow source frequencies of 127 Hz and 145 Hz, along with the deep source frequencies of 130 Hz and 148 Hz, are selected for examination. Throughout the data analysis process, we assume that the source velocity, frequency and range are unknown, and only the seawater SSP information is available.

The velocity and frequency of the two sources are estimated using the MDFMD-v method as shown in Figure 12. The $f_1$ and $f_2$ search interval is 0.1 Hz, and v search interval is 0.1 m/s. The frequencies and velocity estimation results for the shallow source are $f_1=127\mathrm{Hz}$, $f_2=145\mathrm{Hz}$ and $v_0=2.6\mathrm{m}/\mathrm{s}$. The parameters estimation for the deep source are $f_1=130\mathrm{Hz}$, $f_2=148\mathrm{Hz}$ and $v_0=2.6\mathrm{m}/\mathrm{s}$. The estimation results of the parameters are well matched with the experiment.

Based on the velocity estimation result, the horizontal wavenumbers and the corresponding mode depth functions are computed by Equations (2) and (5), respectively. Then, the DSS is used to determine the optimal sign of each mode. The source depth can be estimated by combining the estimated mode signs, the mode depth functions and the modal amplitudes, and the estimated results are shown in Figure 13. Figure 13a shows the loss function defined by th evolution of Equation (14) during the optimization process, utilizing DSS. Figure 13b exhibits the depth ambiguity function computed under various scenarios: individual estimation for $f_1=127\mathrm{Hz}$, $f_2=145\mathrm{Hz}$, joint estimation, and traditional MMP. Similarly, Figure 13c,d show the estimation results for the deep source. Similar to the simulation, the traditional MMP has high sidelobes and large estimation error. For the SIDE algorithm, despite the poor estimation result for the deep source at $f_2$, the joint estimation results are accurate. The joint estimation estimates depths of shallow source $z_s=8\mathrm{m}$ and deep source $z_s=58\mathrm{m}$, and the relative error of the estimation results is less than 11%. In contrast, the results estimated by the traditional MMP method are 134 m and 45 m, respectively. And the relative errors are about 1388% and 16%, respectively, which are significantly larger than those of the SIDE algorithm.

6. Conclusions

The SIDE algorithm is introduced to address the problem of depth estimation for a non-cooperative moving harmonic source when the seabed information is unknown. The SIDE algorithm enables passive depth estimation for a source that emits signals at two adjacent frequencies when only the SSP is known in a range-independent waveguide. In the SIDE algorithm, a one-dimensional sign search method called DSS is proposed to efficiently determine the modal amplitudes signs. DSS takes advantage of the property that the depth ambiguity function distribution closely resembles the Dirichlet kernel function when the signs align with those at the source depth. Simulations and SWellEx-96 experimental data validate the effectiveness of the SIDE algorithm for source depth estimation. In the simulations, we can find unexplained fluctuations in the depth estimation results; although these fluctuations have small peaks, they still need to be further explored in the future. Furthermore, with this method, it takes a long time to observe the modal Doppler effect, and it will be necessary to consider reducing this time in the next works. Finally, the current method focuses solely on scenarios where the source moves at a consistent speed in a horizontally uniform waveguide with minimal variations in water depth. Further investigations are warranted for sources moving at uniform or non-uniform velocities in range-dependent waveguides.