Moving harmonic source depth estimation by MDFMD-v in a weakly range-depend waveguide

To estimate the depth of a non-cooperative moving sound source in a range-depend waveguide, a passive depth estimation method combining depth-polarity search for a harmonic moving source is proposed in this paper. Firstly, the method obtains each mode parameters at the source and horizontal linear array (HLA) locations by utilizing the modal depth function constrained modal Doppler velocity estimation (MDFMD-v) and matrix pencil (MP) methods, respectively. Subsequently, mode amplitudes are calculated by compressive sensing. To determine the polarity of each mode amplitude, a depth-polarity search method is given in this paper. Finally, a depth ambiguity function with mode amplitude polarity is used to obtain the depth estimation results. Simulation results verify the effectiveness of the method.


Introduction
Depth estimation is important in many underwater applications.One classical method for sound source depth estimation is the matched field processing (MFP) [1]- [3], which requires marine environmental information to compute the replica fields, and can utilize the joint search of ambiguity function to estimate source depth and range.However, the traditional MFP method is restricted in practical applications due to the environmental mismatch problem.To alleviate the impact of environmental mismatch, one idea is to simplify the source depth estimation problem.Premus and Conan et al. simplify the source depth estimation to a binary classification problem [4,5], i.e., to discriminate whether the source is near the surface or submerged.These approaches address depth classification problem by extracting depth-depend acoustic features and combining them with algorithms such as the matchedsubspace method [4] or the trapped energy ratio [5].However, these methods are usually applicable to low frequencies and cannot estimate the source depth accurately.Another idea is to extract the normal modal parameters from the sound field data received by an array.A representative one is the data-based MMP method proposed by Yang [6,7], which estimates the mode depth functions through a vertical line array (VLA), then estimates the mode amplitudes using a virtual horizontal array synthesized by a moving source, and finally constructs a depth ambiguity function to achieve source depth estimation.The limitation of the method proposed by Yang is that the velocity and the frequency of the moving source is required, which is ineffective for non-cooperative targets.Besides, the most of the above methods are studied in range-independent waveguides.There are few studies on source depth estimation in range-depend waveguides with uncertain environmental information.
This paper proposes a method for estimating the depth of a non-cooperative target in a range-depend waveguide with unknown seabed parameters based on a horizontal linear array (HLA).And only waveguides satisfying the adiabatic approximation is considered for simplicity.In this case, to estimate the source depth, it is necessary to know the mode depth functions of source and HLA locations to decouple the amplitude of each mode excited by the source.Under the assumption that sound speed profile (SSP) is known, the mode depth functions at the HLA and the source locations can be extracted based on matrix pencil (MP) [8] and modal depth function constrained modal Doppler velocity estimation (MDFMD-v) [9] respectively by combing the difference equation introduced in Ref. [10].And then the amplitude of each mode excited by the source can be extracted by using compressive sensing (CS) [11].Considering the fact that the polarity of each mode amplitude has an impact on the depth estimation, a depth-polarity search method is designed to estimate the polarity of each mode amplitude.Finally, the depth of the sound source can be obtained by the combination of depth ambiguity function and optimal amplitude polarity.

Mode depth function constrained modal Doppler velocity estimation
A range-depend waveguide is shown in Figure 1, with a maximum value of water column depth 1 h and a minimum depth of 2 h .A Cartesian coordinate system is established as shown in Figure 1, the sound speed at coordinates ( ) , rz is ( ) , c r z , and the parameters of the seabed is unknown.In a shallow water waveguide, the sound field can be calculated by the adiabatic approximation theory when the bottom changes slowly [12].Consider a moving harmonic source at the bottom of a slope departing from away from the HLA at a speed 0 v with 0 vc , the sound field at frequency  received by the sensor located at coordinates ( ) , rr rz in the shallow water waveguide can be represented as: 0  .
Equation ( 1) has the same form as the equation ( 4) in Ref. [10].If the SSP at the source location is known and the source contains the adjacent frequencies 1 f and 2 f , the MDFMD-v method in Ref. [10] can be used to jointly estimate the velocity 0 v , the horizontal wavenumber ( , ) rm s kr  and the source frequencies 12 , ff .Then the difference equation proposed in Ref. [8], see equation ( 2), can be utilized to estimate the mode depth function.
where z  is the difference step.

Depth estimation
As shown in Figure 1, the HLA is arranged at the top of the slope with L elements; the distance between two adjacent elements is d ；the first element is located at the (0, ) r z .Therefore, the sound field emitted by the source located at ( , ) ss rz and received by the th n element in the HLA can be expressed as where ( ) mr kr  at this location can be estimated by the MP method [9].When ( , )   mr kr  are estimated, the corresponding mode depth functions


, respectively.After estimating the horizontal wavenumbers and mode depth functions at the source and HLA locations, one still needs to estimate amplitude of the normal modes.CS is used to estimate the complex mode amplitudes in this paper:  , , , ] ,  4) can be efficiently solved by convex optimization, such as the CVX toolkit [13].However, due to the lack of information such as the initial range of the source, it is difficult to obtain an accurate polarity of each mode amplitude, which has a significant effect on the estimation of the source depth.A depth-polarity search method is given in this paper to solve the polarity estimation of each mode amplitude.Before introducing the depthpolarity search method, we can restrict the search range of polarity combinations.Firstly, the preestimation of the source depth is given by , then we can calculate the polarity combinations of the mode amplitudes where the peaks located and the polarity matrix S is modeled as ( , , ) ( , ) , (0, , , , The response of ( , , ) D z n  is concentrated at the depth of the source when the polarity combination is searched to match the actual mode amplitudes polarity combination.To confirm the actual polarity combination of the mode amplitudes, a cost function is defined as The polarity combination where the cost function taken to its minimum value can be considered as the actual polarity of the mode amplitudes.Finally, the source depth can be calculated as Furthermore, a joint estimation function that uses the estimation results of both frequencies can be defined as arg max , , , , , where .The effectiveness of the proposed method will be verified by simulation in Section 3.

Simulation
The range-depend waveguide is shown in Figure 1, and the aperture of the HLA is 1200 m with spacing 5m.The sound source is moving away from the HLA at a velocity , and emits a combination of harmonic signal of frequencies 1 150 Hz f = and 2 160 Hz f = .The signal with a duration of 1000s is selected for depth estimation.Noting that the source velocity and source frequencies are unknown parameters in the simulation process.The sound field is calculated by equation (1), in which the mode depth functions and the horizontal wavenumbers can be calculated by Kraken, and in order to be more realistic, Gaussian white noise with signal-to-noise ratio (SNR) of 0 dB is added to the sound field.
With the source parameters estimation method introduced in Section 2, one can get the estimation of the velocity and two frequencies of the source as shown in the Once the source velocity and frequencies have been estimated, it is possible to compute the horizontal wavenumbers ( , ) m s k rf that correspond to the frequencies and then determine the mode depth functions ( , ) ms rz  at the source location.Figure 2(a) shows the comparison of the estimated mode depth functions with the Kraken-calculated mode depth functions at a frequency of 160Hz (Figure 2 illustrates the estimation results for modes 8, 10, 11 and 12, and the other modes estimated perform similarly.Besides, the results are similar for 150 Hz) when the source depth is 50m.From figure 2, it can be seen that the estimated mode depth functions tally well with the Kraken calculation.The estimation and comparison results of 160Hz are shown in figure 2(b).As mentioned in Section 2, the complex mode amplitudes are computed using the CS method, and the polarities are estimated using the depth-polarity search method.Figure 3 demonstrates the comparison of depth estimation results between the real mode amplitudes (with attenuation and without attenuation) and the mode amplitudes solved by the CS method.In the CS method, the regularization parameter  is set as 0.01 to ensure a better fit to the data when the horizontal wavenumbers are determined.As depicted in figure 3, the influence of low-number mode amplitudes attenuation on depth estimation can be negligible within a certain range.Moreover, if the polarity combination of mode amplitudes is known, estimated depth using the mode amplitudes calculated by CVX program is reliable.To further validate the performance of the method, several simulations with different source depths were performed, and Figure 5 shows the results obtained by using the two frequencies alone and joint estimation.The average relative errors of the source depths estimated at 150Hz, 160Hz and the joint estimate are 26.20%,21.77% and 19.27%.And the depths estimation in the intervals of 0-10m and 150-200m result in the main estimation error (the average relative errors of the estimated source depths are 8.07%, and 3.48% for the 150 Hz, 160 Hz and joint estimation in the intervals of 10-150m, respectively).Therefore, compared with the estimation results using individual frequencies, the errors of the joint estimation are smaller.The depth estimation errors for shallower sources are mainly due to the less highnumber modes are used (the maximum mode number used during the simulation is 15, while the total number of modes is about 30).The depth estimation errors of the deeper sources are mainly due to the estimation of these depths depends on the low-number modes.And the ability of the MDFMD-v method to estimate the wavenumbers of the low-number modes at the source location is affected by the size of the synthetic aperture, which is 2500m in this paper and it is not able to separate the low-number modes at the source location (the first two modes interference span is about 7,000m).

Conclusion
In this paper, the combination of the MDFMD-v and depth-polarity search method achieves the depth estimation of a non-cooperative target sound source in a weakly range-depend waveguide using a HLA.Only two adjacent frequencies emitted by the source and a priori of the sound speed profile in the water column is required in this method.The simulation results show that the method has a good performance in source depth estimation.However, only the moving source at the bottom of the slope in weakly rangedepend waveguide is simulated.In future work, we will extend this method to more complex waveguides and perform experimental data processing and validation.

R
is the horizontal coordinate of the sound source at 0 t = .s z is depth of the source,  are the mode depth function and horizontal wavenumber of mode m at the location ( ), rz , and M is the total number of propagation modes.The dependence of the mode depth function on frequency is neglected in(1).

Figure 1 .
Figure 1.In a weakly range-depend waveguide, the source moves away from the horizontal linear array (HLA) at the speed 0 v.
a fixed and large source range s r for different HLA elements.The frequency domain field received by the HLA is approximated as a linear superposition of exponential functions, and the horizontal wavenumber( , )

.
The estimation of the horizontal wavenumber and mode depth function of the th m normal mode is denoted by ( , ) rm r kr  and ( , ) mr rz controls the relative importance between the data fit and the 1 l norm of the solution.Equation ( and assuming the polarity combination of one peak matches that at the actual source depth.If the peaks of ( , ) pre Dz  occur at' , 1, 2, , represents the sign-taking operation.So, a depth estimation function can be defined as , 1

Figure 2 .
Figure 2. Comparison between the estimated mode depth function with that calculated by Kraken at a frequency of 160Hz.(a) is the result of the mode depth function at source location.(b) is the result of the mode depth function at HLA location.The blue solid line is the estimated mode depth function and the red dashed line is the mode depth function calculated by Kraken.

Figure 3 .
Figure 3.Comparison of estimation results in different cases at a source depth 50 m s z = and a frequency of 160 Hz.The red solid line shows the depth estimation results without considering attenuation, the blue dotted line shows the depth estimation results considering attenuation and the black dashed line shows the depth estimation results using the mode amplitudes estimated by the CVX.

Figure 4 Figure 4 .
Figure 4 illustrates that a peak occurs in the spectrogram of

Figure 5 .
Figure 5. Estimation results at different depths.(a) Depth estimation results estimated by the frequencies 150Hz (black scatter) and 160Hz (red scatter) separately when the source is located at different depths.(b) Depth estimation results estimated by joint estimation when the source is located at different depths.The solid black line is the reference line.

Table 1 .
Results of velocity and frequencies estimated by the MDFMD-v method.