Journal of Ocean University of China  2024, Vol. 23 Issue (4): 871-882  DOI: 10.1007/s11802-024-5742-6

Citation  

GAO Siyu, LI Weilu, ZHANG Yinquan, et al. Extraction of Acoustic Normal Mode Depth Functions Using Range-Difference Method with Vertical Linear Array Data[J]. Journal of Ocean University of China, 2024, 23(4): 871-882.

Corresponding author

LI Weilu, E-mail: liweilu@nmdis.org.cn; ZHANG Yinquan, E-mail: zhyq_nmdis@126.com.

History

Received May 29, 2023
revised September 27, 2023
accepted November 8, 2023
Extraction of Acoustic Normal Mode Depth Functions Using Range-Difference Method with Vertical Linear Array Data
GAO Siyu1),2) , LI Weilu2) , ZHANG Yinquan2) , LI Xiaolei1) , and WANG Ning1)     
1) College of Marine Technology, Ocean University of China, Qingdao 266100, China;
2) National Marine Data and Information Service, Tianjin 300171, China
Abstract: Data-derived normal mode extraction is an effective method for extracting normal mode depth functions in the absence of marine environmental data. However, when the corresponding singular vectors become nonunique when two or more singular values obtained from the cross-spectral density matrix diagonalization are nearly equal, this results in unsatisfactory extraction outcomes for the normal mode depth functions. To address this issue, we introduced in this paper a range-difference singular value decomposition method for the extraction of normal mode depth functions. We performed the mode extraction by conducting singular value decomposition on the individual frequency components of the signal’s cross-spectral density matrix. This was achieved by using pressure and its range-difference matrices constructed from vertical line array data. The proposed method was validated using simulated data. In addition, modes were successfully extracted from ambient noise.
Key words: range difference    depth function extraction    normal mode    
1 Introduction

Acoustic propagation in shallow water can be conveniently described by normal mode theory, which posits that the acoustic field can be understood as the sum of several modal components, with each mode propagating in a dispersive manner. Modal information is highly valuable in underwater acoustics because it contains environmental data that can be used by other algorithms to localize transient sound sources (Yang, 1987, 2014; Nicolas et al., 2006) and infer environmental information along the propagation path (Shang and Wang, 1994; Wang and Song, 2020). The goal of modal filtering is to extract the information of each mode from the received acoustic field. The dispersion characteristics of the acoustic field can be simplified by applying modal filtering (Touze et al., 2009), resulting in the increased temporospatial correlation of the acoustic field. Consequently, modal filtering has widespread practical applications because it can significantly enhance the performance of underwater acoustic signal processing.

Thus far, researchers have proposed various modal filtering methods depending on the scenario. For example, the traditional approach involves measuring the acoustic field on a dense and synchronized water-column-spanning vertical line array (VLA) and using the orthogonality over the depth of the modes at each frequency (Brown et al., 1996; Udovydchenkov and Brown, 2006). This method requires a priori waveguide information about the sound speed profile (SSP) in the water column or the geo-acoustic profile of the bottom. However, a critical issue in using this method is that waveguide information and dense VLAs are often unavailable. Hence, several methods have been proposed to address this problem, including data-derived (Ferris, 1972; Wolf, 1987; Hursky et al., 2001; Roux et al., 2004; Niu et al., 2020), time-warping (Touze, 2009; Bonnel et al., 2010, 2012, 2013; Bonnel and Chapman, 2011; Brown, 2020), and compressive sensing (Yang, 2017; Park et al., 2019; Liang et al., 2020) methods.

The data-derived method (Hursky et al., 2001) is achieved by performing a singular value decomposition (SVD) on the individual frequency components of the temporally averaged, spatial cross-spectral density matrix of the signal (Neilsen et al., 1997; Neilsen and Westwood, 2002). Here, the SVD generates a matrix containing a mutually orthogonal set of basis functions, which are proportional to the depth-dependent normal modes, and a diagonal matrix containing the singular values, which are proportional to the modal source excitations and mode eigenvalues. Time-warping, which is a specific case of unitary similarity transformation, utilizes a unique, nonuniform temporal sampling of the measured signal (Baraniuk and Jones, 1995; Bonnel et al., 2017; Michael, 2020). In this method, the contributions from individual mode numbers are isolated in the frequency spectrum of the time-warped signal, thus allowing the isolation of broadband contributions corresponding to fixed mode numbers using only measurements made at a single receiver. Meanwhile, compressive sensing is a signal processing technique for efficiently acquiring and reconstructing a signal by finding solutions to under-determined linear systems as long as the original signal and the measurement matrix meet certain mathematical conditions. Here, normal modes can be separated on the plane of frequency and azimuth by leveraging compressive sensing, which offers high resolution in azimuth estimation. In turn, the waveform can be recovered through inverse Fourier transformation.

However, a significant issue with the abovementioned data-derived method is that if two or more singular values are nearly equal, the corresponding singular vectors (depth-dependent normal modes) are not uniquely determined. To address this problem, the present paper proposes a new approach to extract depth-dependent normal modes based on VLA data. Here, we introduce a finite range-difference pressure matrix aside from the pressure matrix measured on a VLA. The singular vectors of the two matrices theoretically become the same, except for the order of arrangement, through the use of the SVDs of these two matrices. Furthermore, the singular values of the two matrices are distinct from each other in terms of the difference factors that depend on the normal modes. In this way, we can distinguish the nearly equal singular values of the pressure matrix so that the corresponding singular vectors would be uniquely determined.

The remainder of this paper is organized as follows. In Section 2, we introduce the basic theory discussed in this paper, including the normal mode theory, the data-derived method, and the proposed RD-SVD method based on VLA data. In Section 3, the effectiveness of the proposed method is verified through numerical simulations. Next, we address the degeneracy problem appearing in the traditional data-derived method by extracting the modal depth functions for each order of the normal mode – a task that can be achieved by processing the acoustic pressure and its range-difference matrices. Then, we compare the results obtained by the two methods with those calculated by the Kraken code. The effect of the signal-to-noise ratio on mode extraction is also discussed. Finally, we provide conclusions and suggest future research directions.

2 Theory 2.1 Normal Mode Theory

The propagation of low-frequency acoustic waves in shallow water is described by the normal mode theory (Wang and Shang, 2013). Considering the two-dimensional Helmholtz equation, where sound speed and density depend only on depth z, this is expressed as:

$ \frac{1}{r}\frac{\partial }{{\partial r}}\left({r\frac{{\partial P}}{{\partial z}}} \right) + \rho (z)\frac{\partial }{{\partial z}}\left({\frac{1}{{\rho (z)}}\frac{{\partial P}}{{\partial z}}} \right) + \frac{{{\omega ^2}}}{{{c^2}(z)}}P = - \frac{{\delta (r)\delta (z - {z_s})}}{{2{\rm{ \mathsf{ π} }}r}} . $ (1)

Using the technique of separation of variables, a solution of the unforced equation is given in the form of P(r, z) = Ψ(r)Φ(z). After substituting into the above equation and dividing through byΨ(r)Φ(z), we arrive at

$ \frac{1}{\mathit{\Psi } }\left[ {\frac{1}{r}\frac{{\text{d}}}{{{\text{d}}r}}\left({r\frac{{{\text{d}}\mathit{\Psi } }}{{{\text{d}}r}}} \right)} \right] + \frac{1}{\mathit{\Phi } }\left[ {\rho (z)\frac{{\text{d}}}{{{\text{d}}z}}\left({\frac{1}{{\rho (z)}}\frac{{{\text{d}}\mathit{\Phi } }}{{{\text{d}}z}}} \right) + \frac{{{\omega ^2}}}{{{c^2}(z)}}\mathit{\Phi } } \right] = 0 . $ (2)

The modal depth function Φ depends on both frequency and depth. For the sake of simplicity, modal depth functions are referred to as ‘modes’ in the following text. The modes are classically introduced as a solution to the Sturm-Liouville eigenvalue problem. As such, they are mutually orthogonal and are often used for mode estimation over a VLA of NZ hydrophones. We measured the pressure time series at each of the NZ receiver depths and at each of the NR source-receiver ranges. Then, we obtained the frequency components of the time series by performing a fast Fourier transform (FFT). For each frequency component, the FFT yields an NZ × NR complex matrix P (Neilsen and Westwood, 2002), which is written as follows:

$ \boldsymbol{P} = \left[ {\begin{array}{*{20}{c}} {p({z_1}, {r_1})}&{p({z_1}, {r_2})}& \cdots &{p({z_1}, {r_{{N_R}}})} \\ {p({z_2}, {r_1})}&{p({z_2}, {r_2})}& \cdots &{p({z_2}, {r_{{N_R}}})} \\ \vdots & \vdots & \ddots & \vdots \\ {p({z_{{N_Z}}}, {r_1})}&{p({z_{{N_Z}}}, {r_2})}& \cdots &{p({z_{{N_Z}}}, {r_{{N_R}}})} \end{array}} \right] . $ (3)

The elements of P may be expressed according to the normal mode theory as follows:

$ p({z_i}, {r_j}) = \frac{{\sqrt {2{\rm{ \mathsf{ π} }}} {{\text{e}}^{i{\rm{ \mathsf{ π} }}/4}}}}{{\rho ({z_s})}}\sum\limits_{m = 1}^{{N_M}} {{\phi _m}({z_s})} {\phi _m}({z_i})\frac{{{{\text{e}}^{i{k_m}{r_j}}}}}{{\sqrt {{k_m}{r_j}} }}, $ (4)

where NM is the number of propagating modes; ϕm(z) are the orthonormal, depth-dependent mode functions; km is the horizontal wavenumber of the m-th mode; and the depths z1, z2, ···, zi represent the locations of the corresponding hydrophones of the VLA wherein the source is assumed to remain at a fixed depth of ZS. For the trapped modes, the following equation is a good approximation:

$ \int_0^H {{\phi _n}(z)\phi _m^*(z)} \approx {\delta _{nm}}, $ (5)

where H indicates the depth of the water column, and * represents the conjugate operation. Several excellent papers, such as some on the extraction of modes in different experimental setups (Hursky et al., 2001; Neilsen and Westwood, 2002), have been published in recent years. The SVD of the matrix plays a key role in these discussions.

2.2 Singular Value Decomposition and Its Connection with Normal Modes

Recognizing the similarity between the orthonormality conditions of the modes and the singular vectors, the pressure matrix P can be defined by its elements in Eq. (4) as follows:

$ P = {{\text{e}}^{i{\rm{ \mathsf{ π} }}/4}}\boldsymbol{\varPhi } \boldsymbol{\varLambda } \boldsymbol{F}, $ (6)

where Φ, Λ, $and F are respectively given by the following equations:

$ \boldsymbol{\varPhi } = \frac{1}{{\sqrt {\rho ({z_s})} }}\left[ {\begin{array}{*{20}{c}} {{\phi _1}({z_1})}&{{\phi _2}({z_1})}& \cdots &{{\phi _{{N_M}}}({z_1})} \\ {{\phi _1}({z_2})}&{{\phi _2}({z_2})}& \cdots &{{\phi _{{N_M}}}({z_2})} \\ \vdots & \vdots & \ddots & \vdots \\ {{\phi _1}({z_{{N_Z}}})}&{{\phi _2}({z_{{N_Z}}})}& \cdots &{{\phi _{{N_M}}}({z_{{N_Z}}})} \end{array}} \right], $ (7)
$ \boldsymbol{\varLambda } = \frac{{\sqrt {2{\rm{ \mathsf{ π} }}{N_R}} }}{{\sqrt {\rho ({z_s})} }}\left[ {\begin{array}{*{20}{c}} {\frac{1}{{\sqrt {{k_1}} }}{\phi _1}({z_s})}&0& \cdots &0 \\ 0&{\frac{1}{{\sqrt {{k_2}} }}{\phi _2}({z_s})}& \cdots &0 \\ \vdots & \vdots & \ddots & \vdots \\ 0&0& \cdots &{\frac{1}{{\sqrt {{k_{{N_M}}}} }}{\phi _{{N_M}}}({z_s})} \end{array}} \right], $ (8)
$ \boldsymbol{F} = \frac{1}{{\sqrt {{N_R}} }}\left[ {\begin{array}{*{20}{c}} {\frac{{{{\text{e}}^{i{k_1}{r_1}}}}}{{\sqrt {{r_1}} }}}&{\frac{{{{\text{e}}^{i{k_1}{r_2}}}}}{{\sqrt {{r_2}} }}}& \cdots &{\frac{{{{\text{e}}^{i{k_1}{r_{{N_R}}}}}}}{{\sqrt {{r_{{N_R}}}} }}} \\ {\frac{{{{\text{e}}^{i{k_2}{r_1}}}}}{{\sqrt {{r_1}} }}}&{\frac{{{{\text{e}}^{i{k_2}{r_2}}}}}{{\sqrt {{r_2}} }}}& \cdots &{\frac{{{{\text{e}}^{i{k_2}{r_{{N_R}}}}}}}{{\sqrt {{r_{{N_R}}}} }}} \\ \vdots & \vdots & \ddots & \vdots \\ {\frac{{{{\text{e}}^{i{k_{{N_M}}}{r_1}}}}}{{\sqrt {{r_1}} }}}&{\frac{{{{\text{e}}^{i{k_{{N_M}}}{r_2}}}}}{{\sqrt {{r_2}} }}}& \cdots &{\frac{{{{\text{e}}^{i{k_{{N_M}}}{r_{{N_R}}}}}}}{{\sqrt {{r_{{N_R}}}} }}} \end{array}} \right] . $ (9)

The cross-spectral density matrix of the acoustic pressure field is given by

$ \boldsymbol{C} = {\boldsymbol{PP}^{\text{H}}} = \boldsymbol{\varPhi } \boldsymbol{\varLambda } {\boldsymbol{FF}^{\text{H}}}{\boldsymbol{\varLambda } ^{\text{H}}}{\boldsymbol{\Phi } ^{\text{H}}}, $ (10)

where ‘H’ indicates the conjugate transpose.

If Φ and F are unitary matrices satisfying FFHId, where Id is a unit matrix, then the cross-spectral density matrix can be approximated as follows:

$ \boldsymbol{C} = {\boldsymbol{PP}^{\text{H}}} = \boldsymbol{\varPhi } \boldsymbol{\varLambda } {\boldsymbol{FF}^{\text{H}}}{\boldsymbol{\varLambda } ^{\text{H}}}{\boldsymbol{\varPhi } ^{\text{H}}} \approx \boldsymbol{\varPhi } \boldsymbol{\varLambda } {\boldsymbol{\varLambda } ^{\text{H}}}{\boldsymbol{\varPhi } ^{\text{H}}} . $ (11)

Here, ΛΛH is a diagonal matrix whose elements on the main diagonal are the eigenvalues of the cross-spectral density matrix. When FFHId is satisfied, the cross-spectral density matrix C is similar to the diagonal matrix ΛΛH. Therefore, the SVD of the cross-spectral density matrix C yields the matrix Φ. In turn, this can be used to extract the modes.

The corresponding eigenvectors of these eigenvalues of the matrix Φ are not unique, especially when the cross-spectral density matrix C of the acoustic pressure field has two or more eigenvalues that are close to each other. Thus, in this case, the normal mode depth functions corresponding to these eigenvalues are a linear combination of the eigenvectors of the matrix Φ, which is called ‘normal mode degeneracy’.

Physically, the acoustic source excites several orders of normal modes. In this case, the phenomenon of normal mode degeneracy is likely to occur during SVD if two or more orders of normal modes have similar eigenvalues. At this point, the resulting poor extraction effect of the normal mode depth functions with modal degeneracy significantly limits the practical application of the data-derived normal mode filtering method.

2.3 RD-SVD Based on the VLA Data

The analysis in the previous section shows that the primary reason for normal mode degeneracy is the existence of close eigenvalues. Thus, the key to solving the degeneracy problem is to determine close eigenvalues. Based on this idea, this paper proposes the RD-SVD method to solve the degeneracy problem, and the principle of this method is described below.

In terms of the normal mode expansion, the pressure matrix can be written in the following matrix-product form:

$ {P_{i, j}} = p({z_i}, {r_j}) = {{\text{e}}^{i{\rm{ \mathsf{ π} }}/4}}\sum\limits_{m = 1}^{{N_M}} {\sum\limits_{m' = 1}^{{N_M}} {{\Phi _{i, m}}{\Lambda _{m, m'}}{F_{m', j}}} }, $ (12)

where the elements of the matrices, Φ, Λ, and F are respectively defined as follows:

$ {[\boldsymbol{\varPhi } ]_{i, m}} = \frac{1}{{\sqrt {\rho ({z_s})} }}{\phi _m}({z_i}), $ (13)
$ {[\boldsymbol{\varLambda } ]_{m, m'}} = \frac{{\sqrt {2{\rm{ \mathsf{ π} }}{N_R}} }}{{\sqrt {\rho ({z_s})} }}\frac{1}{{\sqrt {{k_m}} }}{\phi _m}({z_s}){\delta _{m, m'}}, $ (14)
$ {[\boldsymbol{F}]_{m', j}} = \frac{1}{{\sqrt {{N_R}} }}\frac{{{{\text{e}}^{i{k_{m'}}{r_j}}}}}{{\sqrt {{k_{m'}}{r_j}} }} . $ (15)

We introduce the following range-difference acoustic pressure matrices PD, which can be expressed as the products of the matrices:

$ P_{i, j}^D = \frac{{P_{i, j}^ + + P_{i, j}^ - }}{2}, $ (16)
$ P_{i, j}^ + = {P_{i, j + 1}} = \left[ {\begin{array}{*{20}{c}} {p({z_1}, {r_1})}&{p({z_1}, {r_2})}& \cdots &{p({z_1}, {r_{{N_R} - 2}})} \\ {p({z_2}, {r_1})}&{p({z_2}, {r_2})}& \cdots &{p({z_2}, {r_{{N_R} - 2}})} \\ \vdots & \vdots & \ddots & \vdots \\ {p({z_{{N_Z}}}, {r_1})}&{p({z_{{N_Z}}}, {r_2})}& \cdots &{p({z_{{N_Z}}}, {r_{{N_R} - 2}})} \end{array}} \right], $ (17)
$ P_{i, j}^ - = {P_{i, j - 1}} = \left[ {\begin{array}{*{20}{c}} {p({z_1}, {r_3})}&{p({z_1}, {r_4})}& \cdots &{p({z_1}, {r_{{N_R}}})} \\ {p({z_2}, {r_3})}&{p({z_2}, {r_4})}& \cdots &{p({z_2}, {r_{{N_R}}})} \\ \vdots & \vdots & \ddots & \vdots \\ {p({z_{{N_Z}}}, {r_3})}&{p({z_{{N_Z}}}, {r_4})}& \cdots &{p({z_{{N_Z}}}, {r_{{N_R}}})} \end{array}} \right] . $ (18)

Substituting Eqs. (12), (17), and (18) into Eq. (16) yields the following:

$ P_{i, k}^D = {{\text{e}}^{i{\rm{ \mathsf{ π} }}/4}}\sum\limits_{m = 1}^{{N_M}} {\sum\limits_{m' = 1}^{{N_M}} {{\mathit{\Phi } _{i, m}}{\mathit{\Lambda } _{m, m'}}\frac{{{F_{m', j + 1}} + {F_{m', j - 1}}}}{2}} }, $ (19)

where

$\begin{array}{l} {F_{m', j + 1}} + {F_{m', j - 1}} = \\ \frac{1}{{\sqrt {{N_R}} }}\frac{{{{\text{e}}^{i{k_{m'}}{r_j}}}}}{{\sqrt {{k_{m'}}{r_{j + 1}}} }}{{\text{e}}^{i{k_{m'}}\Delta r}} + \frac{1}{{\sqrt {{N_R}} }}\frac{{{{\text{e}}^{i{k_{m'}}{r_j}}}}}{{\sqrt {{k_{m'}}{r_{j - 1}}} }}{{\text{e}}^{ - i{k_{m'}}\Delta r}} . \end{array}$ (20)

In Eq. (20) above, ∆r denotes the step size of the source moving distance. Given that the distance step size is smaller than the source-receiver spacing, Eq. (20) is thus approximated as follows:

$ {F_{m', j + 1}} + {F_{m', j - 1}} = \frac{1}{{\sqrt {{N_R}} }}\frac{{{{\text{e}}^{i{k_{m'}}{r_j}}}}}{{\sqrt {{k_{m'}}{r_j}} }}\left[ {{{\text{e}}^{i{k_{m'}}\Delta r}} + {{\text{e}}^{ - i{k_{m'}}\Delta r}}} \right] . $ (21)

According to Euler’s formula, substituting Eq. (21) into Eq. (19) yields the following:

$ P_{i, k}^D = {{\text{e}}^{i{\rm{ \mathsf{ π} }}/4}}\sum\limits_{m = 1}^{{N_M}} {\sum\limits_{m' = 1}^{{N_M}} {{\mathit{\Phi } _{i, m}}{\mathit{\Lambda } _{m, m'}}\frac{1}{{\sqrt {{N_R}} }}\frac{{{{\text{e}}^{i{k_{m'}}{r_k}}}}}{{\sqrt {{k_{m'}}{r_k}} }}\cos ({k_{m'}}\Delta r)} } . $ (22)

Eq. (22) above shows that PD can be written in the matrix form as follows:

$ \boldsymbol{P}^{\boldsymbol{D}} = {{\text{e}}^{i{\rm{ \mathsf{ π} }}/4}}\boldsymbol{\varPhi } \boldsymbol{K}\boldsymbol{\varLambda } \boldsymbol{F}, $ (23)

where

$ {K_{m, m'}} = \cos ({k_m}\Delta r){\delta _{m, m'}} . $ (24)

The cross-spectral density matrix of the finite-distance differential acoustic field is expressed as

$ \boldsymbol{P}^{\boldsymbol{D}}{[\boldsymbol{P}^{\boldsymbol{D}}]^{\text{H}}} = \boldsymbol{\varPhi } {\boldsymbol{KK}^{\text{H}}}\boldsymbol{\varLambda } {\boldsymbol{\varLambda } ^{\text{H}}}{\boldsymbol{\varPhi } ^{\text{H}}} . $ (25)

Comparing Eqs. (10) and (25), we find that the SVD of the cross-spectral density matrix of the original acoustic field PPH yields the eigenmatrix Φ and the diagonal matrix ΛΛH. Here, the diagonal matrix is KKHΛΛH while the eigenmatrix obtained by performing SVD on the cross-spectral density matrix of the range-difference acoustic pressure field PD[PD]H is also Φ.

The analysis in this section reveals that the SVD and RD-SVD methods can be used to extract the depth function Φ, with the shape of the modes being the same. Furthermore, compared with the SVD method, the RD-SVD method does not increase the complexity of the experimental settings and does not require a priori information of ∆r. In particular, the RD-SVD method introduces an eigenvalue weighting factor cos2(kmr), which can adjust the magnitudes of the eigenvalues, thus affecting the sequence of the eigenvectors of Φ.

Notably, the degeneracy problem can be solved once the weight factor cos2(kmr) generates differences between the originally close eigenvalues obtained by SVD. This is the core idea of using the RD-SVD method to address the normal mode degeneracy problem in this paper. The degeneracy phenomenon can be improved by adjusting ∆r, once the weight factors of the normal modes cos2(kmr) are nearly equal. The effectiveness of this method will be verified by the simulations presented below.

3 Numerical Simulations

The basic principles of the data-derived and RD-SVD methods have been introduced in the previous section. In the current section, we verify the proposed method by conducting simulations in three waveguide environments. First, we obtained the depth functions of normal modes by calculating the acoustic pressure matrices and range-difference pressure matrices of moving acoustic sources at different depths, as well as by diagonalizing the corresponding cross-spectral density matrices. Then, to verify the feasibility of the RD-SVD method, we compared these modes with the normal mode eigenfunctions calculated by the Kraken code. In addition, we discuss in the following sections the influence of noise on mode extraction upon using the RD-SVD method.

3.1 Environment Settings

This section primarily verifies the effectiveness of the proposed approach and its dependence on the parameters. We simulated the acoustic field in a 40-m-deep, range-independent ocean using the Kraken normal mode model. The field was generated by a 350 Hz source. The source was towed on an aperture L from an initial range of R0 = 4 km, and the geometric parameter of the moving source interval and the number of moving sources were respectively represented by ∆r and NR. The bottom sound speed, density, and acoustic attenuation were set at 1753 m s−1, 1900 km m−3, and 0.5 dB λ−1, respectively, where λ is the wavelength. The acoustic field received from the source was sampled by NZ hydrophones at equidistant depths. For example, when NZ = 40 and NR = 64, the 1st to (NR−2)-th sources were used for $ P_{i, j}^ - $ and the 3rd to NR-th sources were used for $ P_{i, j}^ + $. With this setting, there were 16 propagating modes, and only the first eight modes significantly contributed to the received pressure field.

To perform the simulation validation, three different SSPs of the water column were adopted in this study. First, we introduced a typical waveguide environment, which we call ‘Environment A’, with a negative sound speed gradient in summer (Fig.1). This SSP consists of three segments: a constant sound speed of 1520 m s−1 from 0 to 10 m; a linear thermocline with a negative gradient between 10 and 20 m, and the remaining part with a constant sound speed of 1500 m s−1 from 20 to 40 m. The number of moving sources was set at 64, and there were 40 elements in the VLA. The distance between the elements in the VLA was set at 1 m. The range apertures in the simulations were taken as L = 2λmax = 2070.6 m (where λmax represented the maximum modal cycle for the first eight modes).

Fig. 1 Environment A: summer sound speed profile in the water over a half-space.

The second simulation environment is the Pekeris waveguide, and we call this ‘Environment B’. The sound speed is constant at 1510 m s−1 from 0 to 40 m. Everything, except the water sound speed, is the same between Environments A and B.

Fig. 2 Environment B: isovelocity in the water over the same half-space.

The third simulation environment is a waveguide with a springtime SSP measured in the South China Sea (SCS). The waveguide model consisted of a water column with a measured SSP obtained from expendable bathythermograph (XBT) in the SCS, overlaying a two-layer bottom. This two-layer bottom model contains eight parameters: boundary velocities of 1600 m s−1 and 1700 m s−1, a density of 1900 km m−3, acoustic attenuation of 0.2 dB λ−1, and a depth of 10 m for the upper layer. The corresponding values for the bottom layer are 1753 m s−1, 1900 km m−3, and 0.1 dB λ−1 for the bottom layer, and we call this waveguide environment ‘Environment C’. The schematic diagram of this model is shown in Fig.3.

Fig. 3 Environment C: spring sound speed profile in the water with a sediment layer over the half-space.
3.2 Mode Extraction Using the Traditional SVD and RD-SVD Methods Without the Degeneracy Phenomenon

Numerical simulations were first performed for Environment A (Fig.1). The field was generated by a 350 Hz source at a depth of 31 m. The modes extracted by the SVD method and the calculation results of the Kraken code are compared in Fig.4a, where the modes of Nos. 1, 2, 5, 6, 4, 8, 7, and 3 are plotted from left to right, respectively. As shown in the figure, the modes extracted by the SVD method agree well with the calculation results of the Kraken code. The eigenvalues obtained by the SVD method are shown in Fig.4b.

Fig. 4 (a) Depth-dependent normal modes for Environment A, with ZS = 31 m computed using the SVD on the simulated data (red dots) compared with those calculated with the Kraken model (black lines). (b) Singular values from the SVD.

Next, we compared the modes extracted by the RD-SVD method and the calculation results of the Kraken code, where the modes of Nos. 1, 2, 5, 4, 6, 8, 3, and 7 are plotted from left to right (Fig.5a). The modes extracted by the RD-SVD method agree well with the calculation results of the Kraken code. In addition, the eigenvalues obtained by the RD-SVD method are shown in Fig.5b.

Fig. 5 (a) Depth-dependent normal modes for Environment A, with ZS = 31 m computed using the RD-SVD on the simulated data (magenta dots) compared with those calculated with the Kraken model (black lines). (b) Singular values from the RD-SVD.

The comparison between the SVD and RD-SVD methods reveals that both methods can effectively extract the modes and that the shapes of the extracted modes are the same. The only difference between them lies in the arrangement order of the modes, which can be attributed to the varying eigenvalues.

The following numerical simulations were performed in the Pekeris waveguide environment (Environment B), where the source depth was set at 16 m. Comparing the modes extracted by the SVD method and the results of the Kraken code calculation, where the modes of Nos. 1, 4, 2, 6, 7, 3, 5, and 8 are plotted from left to right (Fig.6a). The figure shows that the modes extracted by the SVD method and the calculation results of the Kraken code are in good agreement. The eigenvalues obtained by the SVD method are shown in Fig.6b. Next, we compared the modes extracted by the RD-SVD method and the results of the Kraken code calculation, where the modes of Nos. 1, 2, 4, 7, 3, 6, 5, and 8 are plotted from left to right (Fig.7a). The modes extracted by the RD-SVD method also coincide with the results of the Kraken code calculation. The specific eigenvalues obtained by the RD-SVD method are shown in Fig.7b.

Fig. 6 (a) Depth-dependent normal modes for Environment B, with ZS = 16 m computed using the SVD on the simulated data (red dots) compared with those calculated with the Kraken model (black lines). (b) Singular values from the SVD.
Fig. 7 (a) Depth-dependent normal modes for Environment B, with ZS = 16 m computed using the RD-SVD on the simulated data (magenta dots) compared with those calculated with the Kraken model (black lines). (b) Singular values from the RD-SVD.

We conducted the following numerical simulations under Environment C (Fig.3), in which the source depth was set at 11 m. Then, we compared the modes extracted by the SVD method and the calculation results of the Kraken code, where the modes of Nos. 4, 5, 3, 6, 8, 7, 2, and 1 are plotted from left to right (Fig.8a). As shown in the figure, the modes extracted by the SVD method are in good agreement with the calculation results of the Kraken code. The eigenvalues obtained by the SVD method are presented in Fig.8b.

Fig. 8 (a) Depth-dependent normal modes for Environment C, with ZS = 11 m computed using the SVD on the simulated data (red dots) compared with those calculated with the Kraken model (black lines). (b) Singular values from the SVD.

In addition, the images shown in Fig.9a compare the modes extracted by the RD-SVD method and the calculation results of the Kraken code, where the modes of Nos. 4, 5, 3, 6, 8, 2, 1, and 7 are plotted from left to right. This figure shows that the modes extracted by the RD-SVD method also agree well with the calculation results of the Kraken code, thereby validating the effectiveness of the RD-SVD method in extracting the normal mode depth functions. The specific eigenvalues obtained by the RD-SVD method are shown in Fig.9b.

Fig. 9 (a) Depth-dependent normal modes for Environment C, with ZS = 11 m computed using the RD-SVD on the simulated data (magenta dots) compared with those calculated with the Kraken model (black lines). (b) Singular values from the RD-SVD.
3.3 Mode Extraction Using the RD-SVD Method to Solve the Degeneracy Problem Occurring in the Traditional SVD

Numerical simulations were also performed under Environment A (Fig.1), with the acoustic source depth adjusted to 6 m while all other simulation conditions remained unchanged. The images shown in Fig.10a compare the modes extracted by the SVD method and the theoretical values. The modes of Nos. 4, 3, 5, 6, 2, 7, 8, and 1 are plotted from left to right, respectively. As shown in the figure, the modes of Nos. 3 and 5 that we extracted using the SVD method exhibit certain deviations from the theoretical values. The main reason is that the eigenvalues k2 and k3 obtained by diagonalizing the density matrix are too close to each other (Fig.10b). In turn, this leads to a poor extraction effect of the modes of Nos. 3 and 5, thus indicating the occurrence of modal degeneracy. The modes extracted using the RD-SVD method and the theoretical values are compared in Fig.11a, where the modes of Nos. 4, 3, 5, 6, 2, 1, 8, and 7 are plotted from left to right. The RD-SVD method shows great performance in extracting the normal mode depth functions. The comparison between Figs.11b and 10b indicate that the eigenvalues k2 and k3 have been numerically separated, thereby addressing the degeneracy problem.

Fig. 10 (a) Depth-dependent normal modes for Environment A, with ZS = 6 m computed using the SVD on the simulated data (red dots) compared with those calculated with the Kraken model (black lines). (b) Singular values from the SVD.
Fig. 11 (a) Depth-dependent normal modes for Environment A, with ZS = 6 m computed using the RD-SVD on the simulated data (magenta dots) compared with those calculated with the Kraken model (black lines). (b) Singular values from the RD-SVD.

The following numerical simulations were performed under Environment B (Fig.2), with the acoustic source depth adjusted to 11 m and all other simulation conditions remaining the same. The images shown in Fig.12a compare the normal mode depth functions extracted by the SVD method and the theoretical values, where the modes of Nos. 2, 1, 3, 5, 6, 7, 4, and 8 are plotted from left to right. As can be seen, the modes of Nos. 5 and 6 extracted using the SVD method significantly deviate from the theoretical values. The main reason for this is that the eigenvalues k3, k4, and k5 are close to each other (Fig.12b), resulting in the poor extraction effects of the modes of Nos. 3, 5, and 6, thus demonstrating the occurrence of modal degeneracy. The normal mode depth functions extracted using the RD-SVD method and the theoretical values are compared in Fig.13a, where the modes of Nos. 2, 1, 3, 6, 7, 5, 4, and 8 are plotted from left to right. The figure shows the excellent performance of the RD-SVD method in mode extraction. The main reason for this is that the eigenvalues k3, k4, and k5 corresponding to the modes of Nos. 3, 5, and 6 have been numerically separated (Fig.13b).

Fig. 12 (a) Depth-dependent normal modes for Environment B, with ZS = 11 m computed using the SVD on the simulated data (red dots) compared with those calculated with the Kraken model (black lines). (b) Singular values from the SVD.
Fig. 13 (a) Depth-dependent normal modes for Environment B, with ZS = 11 m computed using the RD-SVD on the simulated data (magenta dots) compared with those calculated with the Kraken model (black lines). (b) Singular values from the RD-SVD.

The following numerical simulations were performed under Environment C (Fig.3). First, we adjusted the acoustic source depth to 10 m, and all other simulation conditions were kept the same. The images shown in Fig.14a compare the modes extracted by the SVD method and the theoretical values. The modes of Nos. 4, 5, 3, 6, 7, 2, 8, and 1 are plotted from left to right, respectively. As can be seen, the modes of Nos. 3 and 6 extracted using the SVD method exhibit deviations from the theoretical values, mainly because the eigenvalues k3 and k4 are relatively close to each other (Fig.14b). The modes extracted using the RD-SVD method and the theoretical values are presented in Fig.15a. In this figure, the modes of Nos. 4, 5, 3, 6, 2, 8, 7, and 1 are plotted from left to right, and as can be seen, the RD-SVD method shows good performance in mode extraction. This is because the eigenvalues k3 and k4 corresponding to the modes of Nos. 3 and 6, which have led to modal degeneracy, have been numerically separated with the RD-SVD method (Fig.15b), thus addressing the degeneracy problem.

Fig. 14 (a) Depth-dependent normal modes for Environment C, with ZS = 10 m computed using the SVD on the simulated data (red dots) compared with those calculated with the Kraken model (black lines). (b) Singular values from the SVD.
Fig. 15 (a) Depth-dependent normal modes for Environment C, with ZS = 10 m computed using the RD-SVD on the simulated data (magenta dots) compared with those calculated with the Kraken model (black lines). (b) Singular values from the RD-SVD.
3.4 Mode Extraction in Dependence on the SNR

We tested the proposed approach under different SNR conditions. Here, SNR is defined as

$ SNR = 10{\log _{10}}\left({\frac{{\sigma _s^2}}{{\sigma _b^2}}} \right), $ (26)

where $ \sigma _s^2 $ is the signal power over the source bandwidth, and $ \sigma _b^2 $ is the noise power for the same bandwidth. Multiple parameters can affect the extraction of modes, including the waveguide environment, acoustic source depth, and SNR conditions. Thus, the current paper uses a simulation example to briefly illustrate the performance of the RD-SVD method in extracting modes under different SNR conditions. The simulation used the Pekeris profile (Environment B), in which the acoustic source depth was set to 11 m and other simulation conditions were unchanged. The images in Figs.16a and 16b show the results of the SVD and RD-SVD methods, respectively, in extracting modes compared with theoretical values under the condition without noise. The images in Figs.16c and 16d show the results of the SVD and RD-SVD methods, respectively, in extracting modes compared with theoretical values under SNR = 35 dB. The images in Figs.16e and 16f show the results of the SVD and RD-SVD methods, respectively, in extracting modes compared with theoretical values under SNR = 20 dB. As expected, the modal estimation errors increase as the SNR decreases for both methods. However, under SNR = 20 dB, the modal estimation errors remain small with the RD-SVD method, especially when modal degeneracy occurs and when the extraction errors of the SVD method are relatively large.

Fig. 16 Depth-dependent normal modes for Environment B, with ZS = 11 m computed using the SVD (red dots) and the RD-SVD (blue dots) compared with those calculated with the Kraken model (black lines). (a) and (b), results with no noise; (c) and (d), results with SNR = 35 dB; and (e) and (f), results with SNR = 20 dB.
4 Conclusions

With the goal of addressing the problem of modal degeneracy that arises due to close eigenvalues, this paper proposed a new approach for extracting acoustic normal mode depth functions based on VLA data in a range-independent waveguide. In this work, we conducted simulations under three waveguide environments, and the effectiveness of the RD-SVD method in extracting normal mode depth functions was verified. We also explained the phenomenon of modal degeneracy and discussed the reasons behind it through an analysis of the modal extraction effect and eigenvalue distribution pattern of the traditional SVD method. The proposed method has been proven to have better separating ability via comparisons between processing results obtained by using the RD-SVD and the traditional SVD methods. We also applied the RD-SVD method to solve the degeneracy problem, which does not increase the complexity of experimental settings and does not require a priori information of ∆r. We believe that by adjusting ∆r when the weight factors of the two degeneracy eigenvalues are nearly equal, the degeneracy phenomenon can be further improved. Finally, the paper illustrates the effectiveness of the RD-SVD method in extracting normal mode depth functions under different signal-to-noise ratios. In the future, we aim to further validate the RD-SVD method using more marine acoustic experimental data.

Acknowledgements

This research was supported in part by the Young Scientists Fund of National Natural Science Foundation of China (No. 42206226), and the National Key Research and Development Program of China (No. 2021YFC3101603).

References
Baraniuk, R. G., and Jones, D. L., 1995. Unitary equivalence: A new twist on signal processing. IEEE Transactions on Signal Processing, 43(10): 2269-2282. DOI:10.1109/78.469861 (0)
Bonnel, J., and Chapman, N. R., 2011. Geoacoustic inversion in a dispersive waveguide using warping operators. Journal of the Acoustical Society of America, 130(2): EL101-EL107. DOI:10.1121/1.3611395 (0)
Bonnel, J., Caporale, S., and Thode, A., 2017. Waveguide mode amplitude estimation using warping and phase compensation. Journal of the Acoustical Society of America, 141(3): 2243-2255. DOI:10.1121/1.4979057 (0)
Bonnel, J., Dosso, S. E., and Chapman, N. R., 2013. Bayesian geoacoustic inversion of single hydrophone light bulb data using warping dispersion analysis. Journal of the Acoustical Society of America, 134(1): 120-130. DOI:10.1121/1.4809678 (0)
Bonnel, J., Gervaise, C., Nicolas, B., and Mars, J. I., 2012. Single-receiver geoacoustic inversion using modal reversal. Journal of the Acoustical Society of America, 131(1): 119-128. DOI:10.1121/1.3664083 (0)
Bonnel, J., Nicolas, B., Mars, J. I., and Walker, S. C., 2010. Estimation of modal group velocities with a single receiver for geoacoustic inversion in shallow water. Journal of the Acoustical Society of America, 128(2): 719-727. DOI:10.1121/1.3459855 (0)
Brown, M. G., 2020. Time-warping in underwater acoustic waveguides. Journal of the Acoustical Society of America, 147(2): 898-910. DOI:10.1121/10.0000693 (0)
Brown, M. G., Viechnicki, J., and Tappert, F. D., 1996. On the measurement of modal group time delays in the deep ocean. Journal of the Acoustical Society of America, 100(4): 2093-2102. DOI:10.1121/1.417919 (0)
Ferris, R. H., 1972. Comparison of measured and calculated normal-mode amplitude functions for acoustic waves in shallow water. Journal of the Acoustical Society of America, 52(3B): 981-989. DOI:10.1121/1.1913204 (0)
Hursky, P., Hodgkiss, W. S., and Kuperman, W. A., 2001. Matched field processing with data-derived modes. Journal of the Acoustical Society of America, 109(4): 1355-1366. DOI:10.1121/1.1353592 (0)
Liang, Y. Q., Zhou, S. H., and Gong, Z. X., 2020. Normal mode separation based on compressive sensing with a horizontal array. Acta Acustica, 45(5): 609-624. (0)
Michael, G. B., 2020. Time-warping in underwater acoustic waveguides. Journal of the Acoustical Society of America, 147(2): 898-910. DOI:10.1121/10.0000693 (0)
Neilsen, T. B., and Westwood, E. K., 2002. Extraction of acoustic normal mode depth functions using vertical line array data. Journal of the Acoustical Society of America, 111(2): 748-756. DOI:10.1121/1.1432982 (0)
Neilsen, T. B., Westwood, E. K., and Udagawa, T., 1997. Mode function extraction from a VLA using singular value decomposition. Journal of the Acoustical Society of America, 101(5): 3025. (0)
Nicolas, B., Mars, J. I., and Lacoume, J. L., 2006. Source depth estimation using a horizontal array by matched-mode processing in the frequency-wavenumber domain. EURASIP Journal on Applied Signal Processing, 2006: 1-16. (0)
Niu, H. Q., Gerstoft, P., Ozanich, E., Li, Z. L., Zhang, R. H., and Gong, Z. X., 2020. Block sparse Bayesian learning for broadband mode extraction in shallow water from a vertical array. Journal of the Acoustical Society of America, 147(6): 3729-3739. DOI:10.1121/10.0001322 (0)
Park, Y., Gerstoft, P., and Seong, W., 2019. Grid-free compressive mode extraction. Journal of the Acoustical Society of America, 145(3): 1427-1442. DOI:10.1121/1.5094345 (0)
Roux, P., Cassereau, D., and Roux, A., 2004. A high-resolution algorithm for wave number estimation using holographic array processing. Journal of the Acoustical Society of America, 115(3): 1059-1067. DOI:10.1121/1.1648321 (0)
Shang, E. C., and Wang, Y. Y., 1994. Tomographic inversion of the El Niño profile by using a matched-mode processing (MMP) method. IEEE Journal of Oceanic Engineering, 19(2): 208-213. DOI:10.1109/48.286643 (0)
Touze, G. L., Nicolas, B., Mars, J. I., and Lacoume, J. L., 2009. Matched representations and filters for guided waves. IEEE Transactions on Signal Processing, 57(5): 1783-1795. DOI:10.1109/TSP.2009.2013907 (0)
Udovydchenkov, I. A., and Brown, M. G., 2006. Modal group time spreads in weakly range-dependent deep ocean. Journal of the Acoustical Society of America, 120(5S): 3062-3062. (0)
Wang, D. Z., and Shang, E. C.,, 2013. Underwater Acoustics. 2nd edition. Science Press, Beijing, 105-106. (0)
Wang, P. Y., and Song, W. H., 2020. Matched-field geoacoustic inversion using propagation invariant in a range-dependent waveguide. Journal of the Acoustical Society of America, 147(6): 491-497. DOI:10.1121/10.0000966 (0)
Wolf, S. N., 1987. Experimental determination of modal depth functions from covariance matrix eigenfunction analysis. Journal of the Acoustical Society of America, 81(S1): S64-S64. (0)
Yang, T. C., 1987. A method of range and depth estimation by modal decomposition. Journal of the Acoustical Society of America, 82(5): 1736-1745. DOI:10.1121/1.395825 (0)
Yang, T. C., 2014. Data-based matched-mode source localization for a moving source. Journal of the Acoustical Society of America, 135: 1218-1230. DOI:10.1121/1.4863270 (0)
Yang, T. C., 2017. Estimating the mode wavenumbers, depth functions, and amplitudes from moving source data using compressive sensing. OCEANS 2017Anchorage, IEEE. (0)