According to the theory of normal-mode, for a point source excitation with circular frequency ω and depth z_s in shallow water, the received SIS I_T at a depth of z_r after propagation over a long distance d can be expressed approximately as in Equation 1 (Grachev, 1993; Gao et al., 2021):

where B_n is shown in Equation 2

P(ω) is the power spectral density of the source, and B_n (ω, z_r, z_s ) and B_m (ω, z_r, z_s ) represent the amplitudes of the nth and mth normal-modes, respectively. M is the number of the propagating normal-modes. Δκ_nm (ω) is the horizontal wavenumber difference between the nth and mth normal-modes, and Δκ_nm (ω) = κ_n (ω) − κ_m (ω). ψ_n (z_r ) and ψ_n (z_s ) are the mode depth functions for the nth normal-mode receiver and source, respectively. ‘∗’ represents conjugation. N_s (ω) is the additive noise. ∑ n = 1 M | B n ( ω ,   z r ,   z s ) | 2 represents the non-interference components of SIS, which vary slowly with the frequency. I n m ( d ,   z r ,   z s ,   ω )   = B n ( ω ,   z r ,   z s ) B m ∗ ( ω ,   z r ,   z s ) exp   ( i Δ κ n m   ( ω )   d )   is called the NMIS for the n–m pair of normal modes, which oscillates with frequency. I ( d ,   z r ,   z s ,   ω )   =   ∑ n ≠ m B n ( ω ,   z r ,   z s ) B m ∗ ( ω ,   z r ,   z s ) exp   ( i Δ κ n m   ( ω )   d ) corresponds to the sum of the different NMIS and represents the interference components of the SIS. When multiple propagation modes exist in the ocean waveguide, the SIS received by a single hydrophone is typically a superposition of the NMIS, non-interference components, and environmental noise. In general, it is difficult to observe accurate fluctuation periods of the NMIS directly owing to the ocean ambient noise and the limitations of the received equipment. It should be noted that we are mainly concerned with NMIS, and the non-interference components can be considered as approximately constant and removed to obtain I_c according to the Equation 3 :

where 〈 f ( x ) 〉 = ∑ i = 0 Z e i x i ( Z = 2 ) represents the polynomial fit term obtained using the least-squares method (Björck, 1990), where e is the fitting coefficient. N ^′ _s (ω) represents the noise after the removal of non-interference components. In subsequent work, our goal is to extract different NMIS from I_c .

The SSA is a classical signal decomposition/extraction method. Assuming that the input sequence is y = [y ₁ ,y ₂ ,.,y_N ], the steps of the SSA algorithm are as follows:

1. The first step is embedding. Embedding creates a Hankel trajectory matrix. It can be regarded as a mapping that transfers a one-dimensional vector y to Hankel matrix H, the H is shown in Equation 4 :

where the numbers of the row and column vectors are L and K, respectively. L + K = N + 1.

2. The second step is to decompose the matrix H using SVD. That is, matrix H is decomposed into the Equation 5 :

where U ∈ ℝ L × R , its column vectors are orthogonal and normalized. V † ∈ ℝ R × K , its row vectors are orthogonal and normalized. Σ ∈ ℝ R × R is a singular value matrix consisting of zero off-diagonal entries and obvious non-zero singular values on the main diagonal entries (Lin and Wu, 2022). The symbol ‘†’ represents the conjugate transpose. Assuming the rank of a matrix H is R, then Σ = diag(σ ₁ ,σ ₂ ,…,σ_R ) (σ ₁ > σ ₂ >… > σ_R ), the matrix H can also be expressed as the Equation 6 :

3. The third step is grouping. The grouping step splits the Hankel matrix H into several groups and sums up the matrices within each group. For I_c , if there are r NMIS, matrix H can be divided into two group, as shown in Equation 7 :

where A = H 1 + H 2 + … + H r denotes the interference components of SIS (I) (It is also possible to obtain Σ r = diag ( σ 1 , σ 2 , … , σ r , 0 , … , 0 ) by setting all the elements in Σ except the first r singular values to 0, and then A = U Σ r V † , A ∈ ℝ L × K ) , B = H r + 1 + H r + 2 … + H R denotes the noisy spectrum ( N s ' ( ω ) ) .

4. The fourth step was diagonal averaging. Taking A as an example, the mapped vector y ^ is obtained using the following steps,as shown in Equation 8 :

The SSA has various grouping methods (Unnikrishnan and Jothiprakash, 2022). The grouping process above divides the matrix into signal and noise subspaces, which completes the denoising operation of I_c . The SSA (Gao, 2016) method divides the matrix H into r + 1 groups (i.e., H ₁, H ₂, H ₃ to H _r , and B) according to the singular value, and then uses diagonal averaging to obtain the corresponding NMIS. This study adopts the first grouping method, and the SSA is unrolled and designed as a neural network model called SSANet. The specific design ideas and basic structure of SSANet are described in the following sections.

It should be noted that the selection of the effective rank r in the above process is crucial. In this study, r was determined based on the discontinuity of singular values. When white noise is present, the distribution of singular values (σ_i ) is characterized by sudden changes. That is, the singular values corresponding to the interference components of SIS(I) are relatively large, whereas those corresponding to the noisy spectrum ( N s ' ( ω ) ) are much smaller. Therefore, the sequence of relative differences between adjacent singular values is denoted as shown in Equation 9:

is referred to as the difference spectrum of the singular value (Wax and Kailath, 1985; Fishler et al., 2002). The effective rank r is then determined using the peak of the difference spectrum, which is the estimated value of the number of NMIS.

Based on the steps of SSA algorithm, we unroll it and design it as SSANet for the extraction of NMIS. The unrolling process of the SSA algorithm and the corresponding relationship with the SSANet are illustrated in Figure 1 . Specifically:

Based on the unrolling of the SSA algorithm mentioned above, we designed a six-layer network structure called SSANet. Assuming the input is the received and normalized I ¯ c ∈ ℝ 1 × N , which undergoes SVD in the SSA algorithm to obtain the unitary matrices U ∈ ℝ L × R and V † ∈ ℝ R × K , and the singular value matrix Σ ∈ ℝ R × R . The structure of SSANet is illustrated in Figure 2 (the parameters of the SSA algorithm are shown in Figure 2 ). The dimensions of the network input, weight, and output parameters for each layer are presented in Table 1 . The specific steps are as follows:

The implementation of SSANet was carried out using Pytorch 1.6.0, with an NVIDIA Quadro GV100 GPU. The weight parameters W ₁ and W ₂ were initialized with Kaming initialization (He et al., 2015). The mean absolute error was selected as the loss function and optimized using the Adam optimizer (Kingma and Ba, 2014), with the learning rate (Smith, 2017) set to 1e−4. The batch size was set as 128. To prevent network overfitting, an early stop strategy was adopted in training; the network training was stopped when the loss on the validation set did not drop within 10 epochs (Liang et al., 2019).

This study utilizes the typical sound field simulation software Kraken (Porter, 1992), which is used to simulate typical winter (isovelocity) and summer (thermocline) waveguide sound speed profiles. The sound speed and medium parameters are shown in Figures 3A, B . In these two waveguide environments, the emission frequencies of the sound source were f ∈ [300,360] Hz and f ∈ [180,220] Hz, respectively. The frequency resolution are 0.3 Hz and 0.2 Hz, respectively. It is assumed that the depths of the known hydrophones are 15 m and 30 m, respectively, whereas the depths and distances of the sound sources are unknown.

For the training samples in the isovelocity waveguide, the depth of the sound source was set as z_s = {1,2,…,30} m. At each source depth, 500 data points were randomly generated within the distance r ∈ [10,15] km. A total of 15,000 samples are generated. Similarly, in the thermocline waveguide, the depth of the sound source was set as z_s = {1,2,…,50} m. At each source depth, 500 data points were randomly generated within r ∈ [20,25] km. In total, 25,000 samples were collected. The two training samples are randomly added with noise at SNR = [−10,10] dB and are divided into training and validation sets in an 8:2 ratio for network training. It should be noted that the SSANet model in this study can extract multiple pairs of NMIS by setting the value of num. In this study, we use the extraction of two pairs of NMIS as an example to introduce the SSANet method. Therefore, for the training set, the true labels are the first two pairs of NMIS with larger interference amplitudes, i.e., num = 2 in the SSANet.

According to the method based on Gao et al. (2021), in an isovelocity environment, when z_r = 17 m, z_s = 7 m, and r = 13 km, the singular value distribution of I is shown in Figure 4A . In a thermocline environment, when z_r = 30 m, z_s = 29 m, and r = 23 km, the singular value distribution of I is shown in Figure 4B . As can be seen in Figure 4 , there are two large singular values, corresponding to the 1st–3rd NMIS(I ₁₃) and the 1st–2nd NMIS (I ₁₂). Based on the method for determining the effective rank r in this study, it is known that only two sets of NMIS dominate in the above environment. Therefore, a test set was generated at the aforementioned sound source depth. The isovelocity waveguide set z_r = 17 m, z_s = 7 m, and r = {10,10.025,10.5,…,15} km, and the thermocline waveguide set z_r = 30 m, z_s = 29 m, and r = {20,20.025,20.5,…,25} km. Each of them generates 201 data, which are added to the noise of SNR = [−10,10] dB with an interval of 2 dB, serving as the test set. The datasets used are presented in Table 2 . In the data preprocessing stage, all the above data samples were normalized using maximum value normalization.

The root-mean-square error (RMSE) and Mean Absolute Error (MAE) were used to evaluate the performance of the methods. A smaller value for both metrics indicates a minor error between the true and predicted values, indicating better extraction of the NMIS. RMSE and MAE are defined in Equation 11 and Equation 12, as follows:

where I ^ n m and I ¯ n m represent the predicted and true NMIS, respectively.

This study compares the SSANet method with the traditional methods: FT, MUSIC, and SSA methods (Gao, 2016). For the SSANet method, we assumed the construction of a Hankel matrix as a square matrix; thus, the network parameters are denoted as L = K = R = N + 1 2 =101. The MAE and RMSE obtained by applying the four mentioned methods for extracting I ₁₃ and I ₁₂ in the isovelocity waveguide under varying SNR are shown in Figures 5A–D . Similarly, the MAE and RMSE obtained by applying the four methods for extracting I ₁₃ and I ₁₂ in a thermocline waveguide under varying SNR are shown in Figures 6A–D . In addition, a random sample was selected from the test dataset in both environments. Figure 7 shows the input I_c for the isovelocity and thermocline waveguides. Figures 8A–F , respectively, represent the results of extracting I ₁₃ and I ₁₂ in the isovelocity waveguide using SSANet compared to other methods. Figures 9A–F , respectively, represent the results of extracting I ₁₃ and I ₁₂ in the thermocline waveguide using SSANet compared to other methods.

From Figures 5 , 6 , 8 , 9 , it can be observed that SSANet provides the overall best extraction results for the phase, amplitude, and oscillation period of the NMIS. However, the FT, MUSIC, and SSA methods exhibit a significant decrease in performance when the amplitude and phase of the NMIS exhibit nonlinear variations with frequency under low SNR conditions. This is because the traditional extraction algorithms are designed by analyzing the physical processes and through handcrafting, while SSANet attempts to automatically discover model information and incorporate NMIS information by optimizing network parameters that are obtained from the training samples (Monga et al., 2021). On the one hand, when there is noise, traditional extraction algorithms do not consider the prior information of the noise, whereas SSANet learns the prior information of the noise during network training, and thus SSANet has stronger noise robustness. On the other hand, when the amplitude and phase of the NMIS exhibit nonlinear changes, the prior assumption of traditional extraction algorithms makes it difficult to extract nonlinear information. The SSANet can learn nonlinear information through training. Overall, it can be said that the trained SSANet is a parameter-optimized version of the SSA algorithm and therefore outperforms traditional extraction algorithms. These results confirm the effectiveness of the SSANet proposed in this study.