Application of a Randomized Algorithm for Extracting a Shallow Low-Rank Structure in Low-Frequency Reverberation

Abstract

The detection performance of active sonar is often hindered by the presence of seabed reverberation in shallow water. Separating the reverberations from the target echo and noise in the received signal is a crucial challenge in the field of underwater acoustic signal processing. To address this issue, an improved Go-SOR decomposition method is proposed based on the subspace-orbit-randomized singular value decomposition (SOR-SVD). This method successfully extracts the low-rank structure with a certain striation pattern. The results demonstrate that the proposed algorithm outperforms both the original Go algorithm and the current state-of-the-art (SOTA) algorithm in terms of the definition index of the low-rank structure and computational efficiency. Based on the monostatic reverberation theory of the normal mode, it is established that the low-rank structure is consistent with the low-frequency reverberation interference striation. This study examines the interference characteristics of the low-rank structure in the experimental sea area and suggests that the interferences of the fifth and seventh modes mainly control the low-rank structure. The findings of this study can be applied to seafloor exploration, reverberation waveguide invariant (RWI) extraction, and data-driven reverberation suppression methods.

1. Introduction

The study of shallow-water active detection requires careful consideration of seabed reverberation, which poses significant challenges and is of utmost importance [1]. In shallow water, seafloor reverberation presents a significant obstacle to active sonar performance, as weak target echoes can easily be overwhelmed by reverberation interference. Reverberation clutter, which is caused by the inhomogeneity of the water column, discrete scatters on the seafloor, non-discrete seabed structures [2], and fish schools [3], exacerbates the undesired constant false alarm rate (CFAR) of active sonar.

To address these challenges, researchers have developed reverberation suppression and signal separation technologies. However, as seabed reverberation and transmitted signals are inseparable in the Fourier domain, processing methods in the wavelet domain [4] and fractional Fourier domain [5,6,7] have become increasingly significant. To enhance data quality, singular value decomposition (SVD) [8] and mode decomposition (MD) [9,10,11,12,13,14] techniques can be implemented at different stages of processing. Some studies suggest that the reverberation waveguide invariant (RWI) $\beta$, calculated by reverberation interference striation (RIS) that represents the regular structure in the time–frequency domain, can be used to develop RWI-based methods for CFAR and estimating the reverberation level (RL) [15,16]. However, these methods focus on a single measurement, ignoring the internal relationship between multiple measurements and resulting in difficulties with data-driven requirements.

In the processing of reverberation data using optimization and matrix computation-based theories, the structural properties of multi-ping data matrices are of significant importance due to their inherent low-rank features. In [17], a unitary transformation is applied to the block Toeplitz matrix of received signals, and a fast approximate algorithm is devised for the approximation of the reverberation low-rank subspace. For the localization problem of receiver failure in reverberant environments, a matrix completion algorithm that relies on the low-rank characteristic is developed for the Hankel matrices of received signals [18]. The non-negative matrix factorization (NMF) algorithm is utilized to calculate the low-rank structure based on the non-negative multiplier criterion and has been demonstrated to be effective for underwater blind-separation cases of linear frequency modulation (LFM) signals [19,20]. Background separation algorithms, such as robust principal component analysis (RPCA) [21] and Go decomposition [22], have been introduced to separate reverberation, target echo, and noise in multi-ping received signals [23,24] when reverberation signals are viewed as interference backgrounds. Furthermore, a combination method of NMF and Go decomposition is presented to detect underwater targets [25].

However, the current multi-ping methods of reverberation separation [23,24,25] face three significant challenges. Firstly, the method’s application is restricted to high-frequency, short-range reverberation, while cases of low-frequency distant reverberation have not been validated. In contrast to high-frequency short-range reverberation, low-frequency distant reverberation, which can be regarded as a two-way long-range propagation process, is more susceptible to the environmental factors of water velocity disturbance [26], geoacoustics [27], waveguide effects [28,29], and environmental uncertainty [30], resulting in more intricate reverberation pings. Consequently, the low-rank structure of reverberation obtained from multi-ping reverberation will exhibit distinct characteristics due to the impact of frequency and propagation distance. Secondly, although the low-rank structure between different reverberation pings is naturally introduced, a strictly physical explanation in specific scenarios has not been observed, resulting in the inability to guarantee the accuracy and effectiveness of the algorithm’s applicability and separation results. Thirdly, the existing methods mainly achieve signal separation through simulation and only a few have been validated at the level of experimental data, which may lead to uncertainties in the practical application of the algorithms.

To tackle the aforementioned issues, an efficient Go-SOR decomposition algorithm is proposed for separating reverberation signals. The algorithm is demonstrated to be effective for low-frequency and distant cases in underwater reverberant environments. Additionally, the results indicate that the proposed algorithm outperforms both the original Go algorithm and the state-of-the-art (SOTA) algorithm in terms of the definition index of the low-rank structure and computational efficiency. Furthermore, the physical significance of the low-rank structure obtained by the algorithm for low-frequency reverberation in shallow water is explicated as follows: under the monostatic reverberation model of the normal mode, the low-rank structure reverberation aligns with the theoretically predicted RIS. The experimental findings in the South China Sea reveal that the low-rank structure reverberation is primarily formed by the interference of the fifth and seventh modes, thereby ensuring the existence of a low-rank structure in shallow seabed reverberation.

1.1. Notations

• c, $c_b$: minimum sound velocity in the water column, the sound velocity of seabed sediment.
• $\mathrm{d}r$, $\mathrm{d}S(t,\theta ,\tau )$, $\mathrm{d}\theta$: width of the elliptical ring, area increment of bottom scattering, increment of polar angle.
• e: base of the natural logarithm.
• $E\left(\beta \right)$: distribution of RWI.
• $\mathbb{E}(·)$: expectation of $(·)$.
• f, $f_{m{m}^{\prime }}$: frequency point, interference frequency of mode m and $m^{\prime }$.
• $g_{mn}(f,t,\theta )$, $g_{mn}(f,t)$: scattering amplitude of mode m and n at polar angle $\theta$, scattering amplitude of mode m and n.
• i: square root of $-1$.
• $I(f,t)$, $\mathit{I}$, ${\mathit{I}}_{\mathrm{mean}}$: reverberation intensity, reverberation intensity matrix, mean value of reverberation intensity matrix.
• $\stackrel{\mathrm{q}}{{\overline{\mathit{I}}}_T}$, $\stackrel{\mathrm{q}}{\overline{\mathit{L}}}$, $\stackrel{\mathrm{q}}{\overline{\mathit{S}}}$, $\stackrel{\mathrm{q}}{\overline{\mathit{G}}}$: q-th ping of reverberation intensity, low-rank structure, sparse signal, and residual signal.
• j: iteration index in the algorithms.
• $k_{rm}$, $k_{rn}$, $k_{r{m}^{\prime }}$, $k_{r{n}^{\prime }}$: horizontal wavenumbers (eigenvalues) of different modes.
• $2l$: focal length of an ellipse.
• $\overline{\mathit{L}}$: low-rank matrix without the ping number (each ping is the same).
• m, $m^{\prime }$, n, $n^{\prime }$: mode numbers.
• $max(·)$, $min(·)$: maximum value of $(·)$, minimum value of $(·)$.
• N: total number of modes.
• p: oversampling the dimension parameter.
• $P_{\theta }(f,t)$, $P(f,t)$: reverberation pressure at polar angle $\theta$, total reverberation pressure.
• $\mathbb{P}{}_{\mathsf{\Omega }}^k(·)$, $\mathbb{P}{}_{\omega }^{\sigma }(·)$: retain the first k-largest absolute entries of a matrix $(·)$ and set other entries to zero (k-largest non-negative projection on matrix domain), hard threshold operator with threshold $\sigma$ on the entries of a matrix $(·)$.
• q, $\mathrm{q}$: index of the power scheme in the paper, index of reverberation pings.
• Q: total number of reverberation pings.
• ${\mathit{Q}}_{\mathit{1}}$,${\mathit{Q}}_{\mathit{2}}$: unitary\orthogonal matrix of QR decomposition.
• r, $r_0$, $r_1$, $r_2$: rank of a matrix, the horizontal distance of monostatic reverberation, the horizontal distance from scatters to the receiver of bistatic reverberation, the horizontal distance from the source to the scatters of bistatic reverberation.
• $\mathbb{R}$: real field.
• ${\mathit{R}}_{\mathit{1}}$, ${\mathit{R}}_{\mathit{2}}$: upper triangular matrix of QR decomposition.
• $R_{mn}(f,t)$: element of random amplitude in the scattering matrix.
• $\mathrm{rank}(·)$, $\mathrm{card}(·)$: rank of a matrix $(·)$, number of non-zero elements in a matrix $(·)$.
• $S_{\eta }(·)$: soft threshold operator for $(·)$ with threshold $\eta$.
• t, $t_{\max }$: reverberation time, maximum iterations.
• $T_{mn}(f,t)$: element of the random scattering matrix without area term.
• $\mathit{U}$, ${\mathit{U}}_r$, ${\mathit{U}}_2$, ${\mathit{V}}_2$, ${\mathit{V}}_r$, ${\mathit{V}}_2$: unitary\orthogonal matrix of SVD.
• $\mathrm{vec}(·)$, $\mathrm{mat}(·)$: vectorization of a matrix $(·)$, matricization of a vector $(·)$.
• $\mathit{X}$, $\mathit{L}$, $\mathit{S}$, $\mathit{G}$: input matrix, low-rank matrix, sparse matrix, residual matrix.
• $z_s$, $z_r$, $z_b$: depth of the source, depth of the receiver, depth of the bottom.
• $\varphi \left(t\right)$: random variable with uniform distribution in $[0,2\pi ]$ at reverberation time t.
• $\mu$, ${\mu }_1$, ${\mu }_2$: constants of Lambert amplitude.
• $\alpha$: sound absorption coefficient.
• $\beta$: value of RWI.
• $\rho \left(\theta \right)$, ${\rho }_b$: circular similarity of the ellipse at polar angle $\theta$, bottom density.
• $\theta$, ${\theta }_m$, ${\theta }_n$: polar angle, grazing angle of mode m, grazing angle of mode n.
• $\mathsf{\sigma }$: hard threshold parameter.
• $\mathsf{\Sigma }$, ${\mathsf{\Sigma }}_r$, ${\mathsf{\Sigma }}_2$: diagonal matrix of SVD.
• $\tau$: width of the transmitted pulse.
• $ϵ$: error margin.
• ${\phi }_m\left(z_s\right)$, ${\phi }_m\left(z_b\right)$, ${\phi }_n\left(z_r\right)$, ${\phi }_n\left(z_b\right)$: eigenfunctions of different depths and mode numbers.
• ${\parallel (·)\parallel }_F$: Frobenius norm of a matrix $(·)$.
• ${(·)}^{*}$, ${(·)}^H$, ${(·)}^{-1}$, ${(·)}^{†}$: complex conjugation of $(·)$, transpose\Hermitian transpose of a real\complex matrix $(·)$, inverse of $(·)$, Moore–Penrose inverse of a matrix $(·)$.
• $〈·〉$: ensemble operation.

2. Theoretical Basis of Seabed Distant Reverberation

2.1. Reverberation Model

A depth-projection schematic of the bistatic bottom reverberation model in shallow water is shown in Figure 1a. The black solid and cyan circles represent the source and receiver, respectively. The geometric prior yields an elliptical ring of width $\mathrm{d}r$ for the received reverberation at time t, with a fixed focal length $2l$. An area element on the elliptical ring is shown as a red square, and the solid blue line with an arrow indicates the projected path from the source to the receiver. Then, it is shown

$$\mathrm{d}r\simeq \frac{c\tau }{2},\mathrm{d}S(t,\theta ,\tau )=r_1\mathrm{d}r\mathrm{d}\theta \simeq \frac{c\tau }{2}{r}_1\mathrm{d}\theta .$$

(1)

Based on the normal mode reverberation theory [31,32], the seabed reverberation field $P_{\theta }(f,t)$ corresponding to the unit scattering area $\mathrm{d}S(t,\theta ,\tau )$ can be expressed as

$$\begin{array}{c}\hfill P_{\theta }(f,t)=\frac{i}{8\pi }\sum _{m=1}^N\sum _{n=1}^N\frac{{\phi }_m\left(z_s\right){\phi }_m\left(z_b\right){\phi }_n\left(z_r\right){\phi }_n\left(z_b\right)g_{mn}(f,t,\theta )}{\sqrt{k_{rm}{k}_{rn}}}\frac{e^{i(k_{rm}{r}_1+k_{rn}{r}_2)}}{\sqrt{r_1{r}_2}}.\end{array}$$

(2)

For simplicity, we omit the index of the polar angle $\theta$ in the scattering amplitude $g_{mn}(f,t,\theta )$:

$$\begin{array}{c}\hfill g_{mn}(f,t)=g_{mn}(f,t,\theta )=T_{mn}(f,t)\sqrt{\mathrm{d}S(t,\theta ,\tau )},\\ \hfill T_{mn}(f,t)=R_{mn}(f,t)e^{i\varphi \left(t\right)}\sqrt{sin{\theta }_{m}sin{\theta }_n},\end{array}$$

(3)

where the bottom scattering amplitude $R_{mn}(f,t)$ can be a empirical constant $\sqrt{\mu }$ in the Lambert scattering case [33], or follows a Rayleigh or K distribution in the stochastic scattering case [34]; $\varphi \left(t\right)$ is a random variable with uniform distribution in $[0,2\pi ]$ and ${\theta }_m$ is the grazing angle of the m-th mode. More information can be found in the Notations section (Section 1.1). Integrating along angle $\theta$ in Equation (2), the reverberation field can be expressed as follows:

$$\begin{array}{cc}\hfill & P(f,t)=\frac{i}{8\pi }\int \sum _{m=1}^N\sum _{n=1}^N\frac{{\phi }_m\left(z_s\right){\phi }_m\left(z_b\right){\phi }_n\left(z_r\right){\phi }_n\left(z_b\right)T_{mn}(f,t)e^{i(k_{rm}{r}_1+k_{rn}{r}_2)}\sqrt{\mathrm{d}S(t,\theta ,\tau )}}{\sqrt{k_{rm}{k}_{rn}{r}_1{r}_2}}\hfill \\ \hfill & =\frac{i\sqrt{c\tau }}{8\sqrt{2}\pi }\sum _{m=1}^N\sum _{n=1}^N\frac{{\phi }_m\left(z_s\right){\phi }_m\left(z_b\right){\phi }_n\left(z_r\right){\phi }_n\left(z_b\right)T_{mn}(f,t)}{\sqrt{k_{rm}{k}_{rn}}}{\int }_0^{2\pi }\frac{e^{i(k_{rm}{r}_1+k_{rn}{r}_2)}}{\sqrt{r_2}}\sqrt{\mathrm{d}\theta }.\hfill \end{array}$$

(4)

According to the geometric properties of an ellipse, $r_1$ and $r_2$ can be written as follows:

$$\begin{array}{c}\hfill r_1=\frac{l^2-{(ct/2)}^2}{lcos\theta -ct/2},r_2=\frac{lctcos\theta -l^2-{(ct/2)}^2}{lcos\theta -ct/2}.\end{array}$$

(5)

The bistatic seabed reverberation in shallow water can be deduced from Equations (4) and (5). We define the circular similarity $\rho \left(\theta \right)$ as

$$\begin{array}{c}\hfill \rho \left(\theta \right)=\frac{ct/2-l}{ct/2-lcos\theta }=\frac{1-2l/\left(ct\right)}{1-2lcos\theta /\left(ct\right)}.\end{array}$$

(6)

When the ratio of the focal length $2l$ to the length of the scattering path $ct$ approaches zero, $\rho \left(\theta \right)\simeq \rho =1$, and Equation (5) can be approximated as

$$\begin{array}{c}\hfill r_1\simeq r_2\simeq \frac{ct}{2}=r_0,\end{array}$$

(7)

which means that the scattering elliptic ring of bistatic reverberation can be regarded as the scattering circle ring of the monostatic case. Substituting Equation (7) into (4) yields the monostatic reverberation field

$$\begin{array}{c}\hfill P(f,t)=\frac{i\sqrt{c\tau }}{8\sqrt{\pi r_0}}\sum _{m=1}^N\sum _{n=1}^N\frac{{\phi }_m\left(z_s\right){\phi }_m\left(z_b\right){\phi }_n\left(z_r\right){\phi }_n\left(z_b\right)T_{mn}(f,t)e^{i(k_{rm}+k_{rn})r_0}}{\sqrt{k_{rm}{k}_{rn}}}.\end{array}$$

(8)

The corresponding schematic of monostatic reverberation is shown in Figure 1b, and the legends remain the same as in Figure 1a.

Figure 1. Schematic of the shallow seabed reverberation. (a) Projection of the in-depth direction of the bistatic reverberation: elliptic ring. (b) Monostatic reverberation: circle ring.

Figure 1. Schematic of the shallow seabed reverberation. (a) Projection of the in-depth direction of the bistatic reverberation: elliptic ring. (b) Monostatic reverberation: circle ring.

2.2. Simulations of Low-Frequency Bottom Reverberation

A sound speed profile (SSP) of the South China Sea is shown in Figure 2a. The SSP represents a negative gradient structure with a water depth of $z_b=88.75$ m. The bottom velocity $c_b$, density ${\rho }_b$, and absorption coefficients are 1664 m/s, 1.9 g/cm${}^3$, and 0.2 dB/$\lambda$, respectively. The depths $z_s$ of the source and $z_r$ of the receiver are 50 m and 65 m, and the reverberation time and frequency bands are 2–6 s and 300–555 Hz, respectively. For the bistatic reverberation case, the focal length $2l$ is set to 200 m.

Considering the reverberation field $P(f,t)$, the reverberation intensity can be expressed as follows:

$$\begin{array}{c}\hfill I(f,t)=P(f,t)P^{*}(f,t).\end{array}$$

(9)

Simulated by Equations (8) and (9), the Lambert RIS for the monostatic case is shown in Figure 3a. By Equations (4) and (9), the Lambert RIS for the bistatic case is shown in Figure 3b. It can be seen that the slope information of the bistatic RIS is nearly the same as that of the monostatic RIS, but the intensity information differs. To measure the RISs, a two-dimensional Fourier transform method [35] (see Appendix A) is applied to calculate the distribution $E\left(\beta \right)$ of RWI. The results are presented in Figure 3c, where the blue solid line and the red dashed line represent $E\left(\beta \right)$ of the monostatic Lambert RIS in Figure 3a and bistatic Lambert RIS in Figure 3b, respectively, and are in good agreement with the peaks and trends. The peak distribution $E\left(\beta \right)$ is conventionally defined as the real value of RWI, i.e., $\beta$. Furthermore, when the simulating condition satisfies the reverberation time $t>2$ s, the circular similarity $\rho \left(\theta \right)$ in Equation (8) satisfies $\rho \left(\theta \right)\approx 1$. Thus, the bistatic bottom reverberation can be well-depicted by a much simpler monostatic bottom reverberation model in both theoretical and simulated bases under the given conditions.

3. Randomized Algorithms and Signal Separation

3.1. Randomized Form of Low-Rank Approximation

The Eckart–Young theorem [36] states that the optimal solution for low-rank approximation can be achieved through SVD under the constraint of the minimum Frobenius norm or spectral norm. However, the efficiency of SVD is limited by the size of the data matrix, and it is computationally expensive for large-dimensional matrices. To overcome this limitation, randomized algorithms based on the Johnson–Lindenstrauss theorem [37] have been introduced into SVD or low-rank approximation methods. Randomized SVD (RSVD) was first proposed by Halko et al. [38], then some researchers developed the subspace-orbit RSVD (SOR-SVD) [39,40] and theoretically provided the approximate error. Additionally, the low-rank approximation derived directly by the bilateral random projection (BRP) approach was first described in Equation (8) in [41], and optimized to the tight rank-r form by Zhou et al. [42], with a theoretical analysis of the approximate error under the constraint of the spectral norm.

The pseudocode implementations of the two randomized low-rank approximation methods, BRP and SOR-SVD, are presented in Algorithms 1 and 2, respectively. BRP and SOR-SVD are demonstrated to be strictly equivalent when $q=0$ and approximately equivalent when $q>0$ (refer to Appendix B). However, the efficiency of SOR-SVD is proved much better than that of BRP in the following content, especially when the power scheme parameter q increases.

Algorithm 1 BRP with the power scheme

Input: $\mathit{X}\in {\mathbb{R}}^{m\times n}\left(m\ge n\right)$, standard Gaussian matrix ${\mathit{A}}_1\in {\mathbb{R}}^{n\times r}$, r, q Output: low-rank matrix ${\mathit{L}}_r$ of rank r 1: $\tilde{\mathit{X}}={\left(\mathit{X}{\mathit{X}}^H\right)}^q\mathit{X}$ 2: ${\mathit{A}}_2=\tilde{\mathit{X}}{\mathit{A}}_1,{\mathit{Y}}_2={\tilde{\mathit{X}}}^H{\mathit{A}}_2,{\mathit{Y}}_1=\tilde{\mathit{X}}{\mathit{Y}}_2$ 3: compute QR decompositions: ${\mathit{Y}}_2={\mathit{Q}}_2{\mathit{R}}_2,{\mathit{Y}}_1={\mathit{Q}}_1{\mathit{R}}_1$ 4: compute low-rank approximation: ${\mathit{L}}_r={\left(\tilde{\mathit{L}}\right)}^{\frac{1}{2q+1}}={\mathit{Q}}_1{\left[{\mathit{R}}_1{\left({\mathit{A}}_2^H{\mathit{Y}}_1\right)}^{-1}{\mathit{R}}_2^H\right]}^{\frac{1}{2q+1}}{\mathit{Q}}_2^H$

Algorithm 2 SOR-SVD with power scheme

Input: $\mathit{X}\in {\mathbb{R}}^{m\times n}\left(m\ge n\right)$, standard Gaussian matrix ${\mathit{A}}_2\in {\mathbb{R}}^{n\times r}$, p, r, q Output: low-rank matrix ${\mathit{L}}_p$ of rank p, rank-r approximate SVD of ${\mathit{L}}_r$ 1: for$j=1:q+1$do 2:    $\mathit{A}{}_1=\mathit{X}{\mathit{A}}_2$ 3:    ${\mathit{A}}_2={\mathit{X}}^H{\mathit{A}}_1$ 4: end for 5: compute QR decompositions: ${\mathit{A}}_1={\mathit{Q}}_1{\mathit{R}}_1,{\mathit{A}}_2={\mathit{Q}}_2{\mathit{R}}_2$ 6: compute the small-dimensional intermediate quantity: $\mathit{M}={\mathit{Q}}_1^H\mathit{X}{\mathit{Q}}_2$ 7: compute the rank-r truncated SVD of ${\mathit{M}}_r={\mathit{U}}_r{\mathsf{\Sigma }}_r{\mathit{V}}_r^H$ 8: compute the rank-p approximation: ${\mathit{L}}_p={\mathit{Q}}_1\mathit{M}{\mathit{Q}}_2^H$ 9: compute the rank-r approximation SVD form: ${\mathit{L}}_r=\left({\mathit{Q}}_1{\mathit{U}}_r\right){\mathsf{\Sigma }}_r{\left({\mathit{Q}}_2{\mathit{V}}_r\right)}^H$

3.2. Algorithms of Signal Separation

As shown in Algorithm 3, the original Go decomposition is a signal separation algorithm based on alternating projection and the BRP–low-rank approximation [22]. For the following optimization problem:

$$\begin{array}{cc}\underset{L,S}{min}& {\parallel \mathit{X}-\mathit{L}-\mathit{S}\parallel }_F^2,\\ \mathrm{s}.\mathrm{t}.& rank\left(\mathit{L}\right)\le r,\\ & card\left(\mathit{S}\right)\le k,\end{array}$$

(10)

the solution of the Go decomposition yields the low-rank background $\mathit{L}$, sparse target signal $\mathit{S}$, and residual noise $\mathit{G}$ corresponding to the input matrix $\mathit{X}$.

As discussed in Section 3.1, the efficiency of BRP in the Go decomposition is lower than to that of SOR-SVD. Therefore, replacing BRP with SOR-SVD can reduce the computing time without sacrificing approximate precision. An improved algorithm, named Go-SOR, is presented in Algorithm 4. The SOTA randomized-separation algorithm, ALM-SOR-SVD RPCA [40], is also introduced in Algorithm 5 for comparison.

Algorithm 3 Original Go decomposition

Input: $\mathit{X}\in {\mathbb{R}}^{m\times n}\left(m\ge n\right)$, standard Gaussian matrix ${\mathit{A}}_1\in {\mathbb{R}}^{n\times r}$, r, q, k, $ϵ$, $t_{\max }$ Output: low-rank matrix $\mathit{L}$, sparse matrix $\mathit{S}$, residual matrix $\mathit{G}$ 1: while true do 2:    $t=t+1$ 3:    ${\mathit{L}}_t=\mathrm{BRP}\left(\mathit{X}-{\mathit{S}}_{t-1},{\mathit{A}}_1,r,q\right)$    (Algorithm 1) 4:    ${\mathit{S}}_t=\mathbb{P}{}_{\mathsf{\Omega }}^k(\mathit{X}-{\mathit{L}}_t)$    (Equation (11)) 5:    if $\parallel \mathit{X}-{\mathit{S}}_t-{\mathit{L}}_t{\parallel }_F^2/{\parallel \mathit{X}\parallel }_F^2\le ϵ$ or $t\ge t_{\max }$ then 6:      $\mathbf{\text{break}}$ 7:    end if 8: end while 9: compute $\mathit{G}=\mathit{X}-\mathit{S}-\mathit{L}$

Algorithm 4 Go-SOR decomposition

Input: $\mathit{X}\in {\mathbb{R}}^{m\times n}\left(m\ge n\right)$, standard Gaussian matrix ${\mathit{A}}_1\in {\mathbb{R}}^{n\times r}$, r, q, k, $ϵ$, $t_{\max }$ Output: low-rank matrix $\mathit{L}$, sparse matrix $\mathit{S}$, residual matrix $\mathit{G}$ 1: while true do 2:    $t=t+1$ 3:    ${\mathit{L}}_t=\mathrm{SOR-SVD}\left(\mathit{X}-{\mathit{S}}_{t-1},{\mathit{A}}_1,p,r,q\right)$, where $p=r$   (Algorithm 2) 4:    ${\mathit{S}}_t=\mathbb{P}{}_{\omega }^{\sigma }(\mathit{X}-{\mathit{L}}_t)$    (Equation (12)) 5:    if $\parallel \mathit{X}-{\mathit{S}}_t-{\mathit{L}}_t{\parallel }_F^2/{\parallel \mathit{X}\parallel }_F^2\le ϵ$ or $t\ge t_{\max }$ then 6:      $\mathbf{\text{break}}$ 7:    end if 8: end while 9: compute $\mathit{G}=\mathit{X}-\mathit{S}-\mathit{L}$

Algorithm 5 ALM-SOR-SVD RPCA

Input: $\mathit{X}\in {\mathbb{R}}^{m\times n}\left(m\ge n\right)$, standard Gaussian matrix ${\mathit{A}}_1\in {\mathbb{R}}^{n\times r}$, p, r, q, $ϵ$, $t_{\max }$ Output: low-rank matrix $\mathit{L}$, sparse-plus-noise matrix $\mathit{E}$ 1: while true do 2:    $t=t+1$ 3:    $\left(\mathit{U},\mathsf{\Sigma },\mathit{V}\right)=\mathrm{SOR-SVD}\left(\mathit{X}-{\mathit{E}}_{t-1}+{\mu }_t^{-1}{\mathit{Y}}_{t-1},{\mathit{A}}_1,p,r,q\right)$   (Algorithm 2) 4:    ${\mathit{L}}_t=\mathit{U}{S}_{{\mu }_{t-1}^{-1}}\left(\mathsf{\Sigma }\right){\mathit{V}}^H$ 5:    ${\mathit{E}}_t=S_{\lambda {\mu }_{t-1}^{-1}}(\mathit{X}-{\mathit{L}}_t+{\mu }_{t-1}^{-1}{\mathit{Y}}_{t-1})$ 6:    ${\mathit{Y}}_t={\mathit{Y}}_{t-1}+{\mu }_{t-1}(\mathit{X}-{\mathit{L}}_t-{\mathit{E}}_t)$ 7:    ${\mu }_t=max(\rho {\mu }_{t-1},\overline{\mu })$ 8:    if $\parallel \mathit{X}-{\mathit{E}}_t-{\mathit{L}}_t{\parallel }_F^2/{\parallel \mathit{X}\parallel }_F^2\le ϵ$ or $t\ge t_{\max }$ then 9:      $\mathbf{\text{break}}$ 10:    end if 11: end while

Note that the output matrix $\mathit{E}$ of the ALM-SOR-SVD RPCA differs from that of the Go-SOR decomposition. Alternatively, the sparse matrix $\mathit{S}$ and residual matrix $\mathit{G}$ can be separated from the sparse-plus-noise matrix $\mathit{E}$ by applying a k-largest non-negative projection $\mathbb{P}{}_{\mathsf{\Omega }}^k(·)$ as follows:

$$\mathbb{P}{}_{\mathsf{\Omega }}^k:{\mathit{S}}_{ij}=\left\{\begin{array}{c}{\mathit{E}}_{ij}\left|{\mathit{E}}_{ij}=max\left(\left|\mathit{E}\right|,l\right)\right.,1\le l\le k,\\ 0,\end{array}\right.\mathit{G}=\mathit{X}-\mathit{S}-\mathit{L}.$$

(11)

where the subscript “$ij$” and the operation “max” with a parameter “l” indicates the index of elements and the l-th maximum in the matrix, respectively. Similarly, the element-wise hard threshold operator $\mathbb{P}{}_{\omega }^{\sigma }(·)$ with a threshold $\sigma$ on $\mathit{X}-\mathit{L}$ is applied in Go-SOR and is defined as follows:

$$\mathbb{P}{}_{\omega }^{\sigma }:{\mathit{S}}_{ij}=\left\{\begin{array}{c}{(\mathit{X}-\mathit{L})}_{ij},\left|{(\mathit{X}-\mathit{L})}_{\mathit{ij}}\right|\ge \sigma ,\\ 0,\left|{(\mathit{X}-\mathit{L})}_{\mathit{ij}}\right|<\sigma .\end{array}\right.$$

(12)

4. Results

4.1. Low-Rank Structure of the Shallow Reverberation Experiment

Experimental data of shallow seabed reverberation were obtained from the 2015 reverberation experiment in the South China Sea. The SSP of the sea area is shown in Figure 2a, and the geoacoustic parameters are the same as those used in the simulation. The receiver was fixed at a depth of 65 m, and 31 explosive sources were released at a depth of 50 m around the receiver. The relative positions of the sources and the receiver are shown in Figure 2b, where the symbols of the blue circle and red asterisks indicate the receiver and sources, respectively. The distances between the sources and receiver are approximately 200 m. According to the results in Section 2, it is reasonable to simplify the experimental bistatic reverberation by monostatic reverberation when the reverberation time is $t>2$ s and the focal length is $2l\sim 200$ m.

The moving signal was not predetermined during the sea trial. Low-frequency analysis and recording (LOFAR) spectra of the experimental reverberation data were obtained by short-time Fourier transform (STFT) and labeled LOFARgram. LOFARgrams $\stackrel{1}{{\overline{\mathit{I}}}_T}$, $\stackrel{11}{{\overline{\mathit{I}}}_T}$, $\stackrel{21}{{\overline{\mathit{I}}}_T}$, and $\stackrel{31}{{\overline{\mathit{I}}}_T}$ of ping numbers 1, 11, 21, and 31 are displayed in Figure 4a, Figure 4b, Figure 4c, and Figure 4d, respectively. It can be observed that these spectrograms vary with different pings, and there is no obvious regularity for the environmental uncertainty and difference in positions. The structural similarity index measure (SSIM) values of every two pings are shown in

Table 1. SSIM of every two different reverberation pings.

Ping Number	1	11	21	31
1	-	0.0845	0.1396	0.1030
11	0.0845	-	0.1002	0.0860
21	0.1396	0.1002	-	0.1080
31	0.1030	0.0860	0.1080	-

, and the results are concentrated around 0.1, which indicates the non-correlation of different reverberation pings.

The input matrix $\mathit{X}$ is obtained by the vectorization of different LOFARgrams

$$\mathit{X}=\left[\mathrm{v}\mathrm{e}\mathrm{c}\left(\stackrel{1}{{\overline{\mathit{I}}}_T}\right),\mathrm{v}\mathrm{e}\mathrm{c}\left(\stackrel{2}{{\overline{\mathit{I}}}_T}\right),\dots ,\mathrm{v}\mathrm{e}\mathrm{c}\left(\stackrel{31}{{\overline{\mathit{I}}}_T}\right)\right].$$

(13)

The initialization of Go-SOR (found in Algorithm 4), Go (found in Algorithm 3), and ALM-SOR-SVD RPCA (found in Algorithm 5) are $r=p=1$, $q=1$, and $t_{max}=100$. Because there is no preset target for the experiment, the values of k and $\sigma$ are meaningless. The low-rank matrix $\stackrel{\mathrm{q}}{\overline{\mathit{L}}}$, sparse matrix $\stackrel{\mathrm{q}}{\overline{\mathit{S}}}$, and residual matrix $\stackrel{\mathrm{q}}{\overline{\mathit{G}}}$ of q-th ping LOFARgram $\stackrel{\mathrm{q}}{{\overline{\mathit{I}}}_T}$ are obtained by the matricization of the output

$$\stackrel{\mathrm{q}}{\overline{\mathit{L}}}=mat\left[\mathit{L}\left(:,\mathrm{q}\right)\right],\stackrel{\mathrm{q}}{\overline{\mathit{S}}}=mat\left[\mathit{S}\left(:,\mathrm{q}\right)\right],\stackrel{\mathrm{q}}{\overline{\mathit{G}}}=mat\left[\mathit{G}\left(:,\mathrm{q}\right)\right].$$

(14)

The computing times of the three algorithms are listed in

Table 2. Reverberation data processing times for $r=p=1$, $q=1$, and $t_{max}=100$.

Algorithms	Go-SOR	Original Go	ALM-SOR-SVD RPCA
Computing time (s)	0.44	55.20	0.76

. These algorithms ran on a personal desktop computer equipped with an Intel i7 10700 CPU and 64 GB RAM, using Matlab 2022b software for implementations. It can be seen that the proposed Go-SOR algorithm (computing time 0.44 s) has a certain improvement compared to the SOTA performance of ALM-SOR-SVD RPCA (computing time 0.76 s). However, the original Go algorithm has the poorest performance (computing time 55.20 s) and takes two orders of magnitude longer than the other two algorithms. The improvement in efficiency of the proposed Go-SOR algorithm can be attributed to two factors: (a) the use of SOR-SVD instead of BRP in the low-rank approximation, and (b) the use of the hard threshold operator $\mathbb{P}{}_{\omega }^{\sigma }(·)$ instead of the non-negative projection operator $\mathbb{P}{}_{\mathsf{\Omega }}^k(·)$ in each iteration [43]. The first factor was explained in detail earlier, and the second factor will be briefly discussed below. The non-negative projection operator $\mathbb{P}{}_{\mathsf{\Omega }}^k(·)$ involves time-consuming element-sorting, while the hard threshold operator $\mathbb{P}{}_{\omega }^{\sigma }(·)$ is a more efficient operator that involves element subtraction. Additionally, the parameters of ALM-SOR-SVD RPCA are difficult to determine, whereas those of Go-SOR are relatively simple, which is a significant advantage in practical applications.

The low-rank structures $\stackrel{\mathrm{q}}{\overline{\mathit{L}}}$ for the three algorithms are shown in Figure 5. Figure 5a–c represent the results of Go-SOR, original Go, and ALM-SOR-SVD RPCA, respectively. The graph in Figure 5d symbolizes the mean LOFARgram of 31 reverberation pings. Although all four graphs show similar interference striation patterns around 2 s to 4.5 s, the low-rank structure of Go-SOR is the most distinct. The definition index, which is commonly used in optical image quality assessment, is calculated using four different matrix gradient methods: Energy of gradient (EOG), Tenengrad, Brenner, and Laplacian [44,45]. The results in

Table 3. Normalized definition indices of the low-rank structure reverberation via different algorithms.

	Go-SOR	Original Go	ALM-SOR-SVD RPCA
EOG	1	0.8352	0.9100
Tenengrad	1	0.8398	0.9304
Brenner	1	0.7472	0.6825
Laplacian	1	0.8406	0.8974

demonstrate that the proposed Go-SOR algorithm outperforms the other two algorithms in all four gradient methods, with the definition indices normalized to 1 for ease of understanding and comparison. Therefore, the proposed Go-SOR algorithm achieves the best performance in terms of both the image quality of the generated low-rank structure and computing efficiency.

The corresponding distributions $E\left(\beta \right)$ of RWIs in Figure 5 are shown in Figure 6a. The blue solid, red dashed, green dashed, and black dotted lines represent $E\left(\beta \right)$ for the low-rank structure in Figure 5a (via Go-SOR), Figure 5b (via Go), Figure 5c (via ALM-SOR-SVD RPCA), and mean LOFARgram in Figure 5d, respectively. The curves show similar distributions of RWIs, indicating that the striation features of the four subfigures in Figure 5 are almost consistent. The coincidence of the peaks demonstrates the close $\beta$ values, making the processing of low-rank structure reverberations reliable.

4.2. Uniformity of RIS, Mean LOFARgram, and Low-Rank Structure

The contents in Figure 5 and Figure 6a illustrate the relationship between the RIS, mean LOFARgram, and low-rank structure, highlighting their uniformity.

For a simplified monostatic reverberation case in the experimental environment, the stochastic scattering amplitude $R_{mn}(f,t)$ in Equation (3) is assumed to be stationary, and the first and second moments exist [46]

$$\mathbb{E}\left[R_{mn}(f,t)\right]=〈R_{mn}(f,t)〉=\sqrt{{\mu }_1},\mathbb{E}\left\{{\left[R_{mn}(f,t)\right]}^2\right\}=〈{\left[R_{mn}(f,t)\right]}^2〉={\mu }_2.$$

(15)

Using Equation (8), the ensemble reverberation intensity element $〈I(f,t)〉$ of intensity matrix $〈\mathit{I}〉$ with frequency index f and time index t can be expressed as follows:

$$\begin{array}{cc}\hfill & 〈I(f,t)〉=〈P(f,t)P^{*}(f,t)〉=\frac{c\tau }{64\pi r_0}\sum _{m=1}^N\sum _{n=1}^N\frac{{\phi }_m\left(z_s\right){\phi }_m\left(z_b\right){\phi }_n\left(z_r\right){\phi }_n\left(z_b\right)}{\sqrt{k_{rm}{k}_{rn}}}\times \hfill \\ & \sum _{m^{\prime }=1}^N\sum _{n^{\prime }=1}^N\frac{{\phi }_{m^{\prime }}\left(z_s\right){\phi }_{m^{\prime }}\left(z_b\right){\phi }_{n^{\prime }}\left(z_r\right){\phi }_{n^{\prime }}\left(z_b\right)〈T_{mn}(f,t)T_{m^{\prime }{n}^{\prime }}^{*}(f,t)〉e^{i(k_{rm}+k_{rn}-k_{r{m}^{\prime }}-k_{r{n}^{\prime }})r_0}}{\sqrt{k_{r{m}^{\prime }}{k}_{r{n}^{\prime }}}}.\hfill \end{array}$$

(16)

where

$$\begin{array}{cc}\hfill 〈T_{mn}(f,t)T_{m^{\prime }{n}^{\prime }}^{*}(f,t)〉& =〈R_{mn}(f,t)R_{m^{\prime }{n}^{\prime }}^{*}(f,t)e^{i\left[\varphi \right(t)-\varphi (t\left)\right]}〉\sqrt{sin{\theta }_{m}sin{\theta }_{n}sin{\theta }_{m^{\prime }}sin{\theta }_{n^{\prime }}}\hfill \\ \hfill & =\left\{\begin{array}{c}{\mu }_1\sqrt{sin{\theta }_{m}sin{\theta }_{n}sin{\theta }_{m^{\prime }}sin{\theta }_{n^{\prime }}},ifm\ne m^{\prime }orn\ne n^{\prime },\\ {\mu }_2\sqrt{sin{\theta }_{m}sin{\theta }_{n}sin{\theta }_{m^{\prime }}sin{\theta }_{n^{\prime }}},ifm=m^{\prime }andn=n^{\prime }.\end{array}\right.\hfill \end{array}$$

(17)

Equations (16) and (17) show that the ensemble expression $〈\mathit{I}〉$ of stochastic reverberation $\mathit{I}$ loses randomness and is theoretically consistent with a certain Lambert reverberation ${\mathit{I}}_{\mathrm{Lambert}}$ with parameters ${\mu }_1$ and ${\mu }_2$. In the experimental data treatment, the ensemble reverberation intensity $〈\mathit{I}〉$ is approximated by the mean LOFARgram ${\mathit{I}}_{\mathrm{mean}}$ of the Q ($Q=31$) pings. Thus, we naturally obtain the following equation:

$${\mathit{I}}_{\mathrm{Lambert}}=〈\mathit{I}〉\approx {\mathit{I}}_{\mathrm{mean}}=\frac{1}{Q}\sum _{\mathrm{q}=1}^Q\stackrel{\mathrm{q}}{{\overline{\mathit{I}}}_T}.$$

(18)

Therefore, the uniformity of the RIS (symbolizing the Lambert reverberation ${\mathit{I}}_{\mathrm{Lambert}}$) and mean LOFARgram ${\mathit{I}}_{\mathrm{mean}}$ is demonstrated. When considering a low-rank structure, the initialized condition $r=1$ in Go-SOR leads to the following relation:

$$\stackrel{{\mathrm{q}}_1}{\overline{\mathit{L}}}=\stackrel{{\mathrm{q}}_2}{\overline{\mathit{L}}}=\overline{\mathit{L}},$$

(19)

which means that the low-rank structures of different ping numbers ${\mathrm{q}}_1$ and ${\mathrm{q}}_2$ are consistent. Another condition of the non-existence of a moving target implies

$$\mathit{0}=〈\stackrel{\mathrm{q}}{\overline{\mathit{S}}}+\stackrel{\mathrm{q}}{\overline{\mathit{G}}}〉,$$

(20)

where $\mathit{0}$ is the zero matrix. Meanwhile, each ping of the reverberation intensity matrix can be written as

$${\stackrel{\mathrm{q}}{\overline{\mathit{I}}}}_T=〈\mathit{I}〉+\left[{\stackrel{\mathrm{q}}{\overline{\mathit{I}}}}_T-〈\mathit{I}〉\right].$$

(21)

It can be observed that factorized terms $〈\mathit{I}〉$ and $\left[{\stackrel{\mathrm{q}}{\overline{\mathit{I}}}}_T-〈\mathit{I}〉\right]$ on the right-hand side (RHS) of Equation (21), which satisfy Equations (19) and (20), respectively, are unique. It is demonstrated that the low-rank structure obtained by the proposed algorithm is exactly the unique solution of the ensemble form of experimental RIS (mean LOFARgram ${\mathit{I}}_{\mathrm{mean}}$) under specific prior conditions. By combining these equations, the consistency between the low-rank structure and Lambert RIS can be theoretically established. Therefore, the RIS ${\mathit{I}}_{\mathrm{Lambert}}$, mean LOFARgram ${\mathit{I}}_{\mathrm{mean}}$, and low-rank structure $\overline{\mathit{L}}$ are related by the following expressions:

$$\overline{\mathit{L}}={\mathit{I}}_{\mathrm{Lambert}}\approx {\mathit{I}}_{\mathrm{mean}}.$$

(22)

4.3. Interference Properties of the Low-Rank Structure

Due to the influence of the time-varying environment, the interference properties of single-ping reverberation change sharply. As shown in Figure 5, obvious periodic interference phenomena can be observed in the low-rank matrices for multi-ping reverberation from 2 s to 4.5 s. By the uniformity of RIS and the low-rank structure, the interference frequency f of experimental data can be theoretically analyzed through Equation (16), which is the simplified monostatic model for shallow seabed reverberation.

The concept of interference frequency $f_{m{m}^{\prime }}$ between modes m-th and $m^{\prime }$-th is introduced to analyze the one-way process of bottom reverberation interference patterns:

$$f_{m{m}^{\prime }}=\frac{\left|k_{rm}-k_{r{m}^{\prime }}\right|}{2\pi }\frac{c}{2}.$$

(23)

The dominant interference frequencies for the reverberation experimental data are determined through spectral analysis of the low-rank structure shown in Figure 5. For accurate frequency estimation, the reverberation intensity curves at different frequencies are detrended during pretreatment. We use the multiple signal classification (MUSIC) method to improve the resolution of estimated results [47]. The estimation of interference frequencies for the low-rank structure is shown in Figure 7a, where the blue solid line, black dashed line, and red dotted line represent the results of 300, 400, and 500 Hz, respectively. The frequencies corresponding to the highest peaks match the direct calculations of the fifth and seventh modes by Equation (23) under all three different frequencies.

The simulated fifth-and-seventh-only RIS, which is similar to the low-rank structure shown in Figure 5a between 2 s to 4.5 s, is described in Figure 7b. For comparison, distributions $E\left(\beta \right)$ of the fifth-and-seventh-only RIS, monostatic Lambert RIS in Figure 3a, and low-rank structure in Figure 5a are described in Figure 6b, with the black dotted line, red dashed line, and blue solid line, respectively. We find that $E\left(\beta \right)$ of the simulated fifth-and-seventh-only RIS matches well with the experimental low-rank structure. Because of the difference in the mode number in the simulation, the monostatic Lambert RIS, which is calculated under all modes, presents a more complex pattern compared to the fifth-and-seventh-only RIS, and the difference of $E\left(\beta \right)$ in Figure 6b is apparent. The differences in the RWI corresponding to the peaks in Figure 6b are small enough, despite the deviation of the monostatic $E\left(\beta \right)$ curve from the others. This demonstrates the conclusive control of the fifth and seventh modes in the experiment.

5. Discussion

In this study, we employed randomized separation algorithms to extract the low-rank structure of shallow seabed reverberation. We propose a modified Go-SOR decomposition method based on the equivalence of SOR-SVD and BRP, which exhibits superior comprehensive performance compared to the SOTA algorithm when processing low-frequency reverberation experimental data. We demonstrated the uniformity between the simulated Lambert RIS, low-rank structure reverberation, and ensemble intensity structure based on the simplified monostatic reverberation model. The spectral analysis of the low-rank structure reverberation and theoretical calculation of interference frequencies indicate that the interference structure is primarily controlled by the fifth and seventh modes in the experimental environment of the South China Sea.

The original contributions of this work are:

• A novel signal separation technique, named Go-SOR, was proposed and evaluated for processing reverberation experimental data. The results show that Go-SOR outperforms the SOTA algorithm ALM-SRO-SVD RPCA in terms of the computing time and definition index. Furthermore, we established the conditional equivalence of low-rank approximation between the SOR-SVD of the proposed algorithm and the BRP of the original Go algorithm.
• A bistatic low-frequency reverberation simulation model based on the normal mode theory was developed. When the transition indicator circular similarity is $\rho \left(\theta \right)\to 1$, the model can be transformed into a monostatic distant seabed reverberation model.
• The equivalence between the low-rank structure by the proposed algorithm and RIS was described, which provides the interpretability of the algorithm’s processing results. Our findings suggest that the study of RIS can be transformed into the study of the low-rank structure reverberation obtained from the data. This provides the possibility for the study of the data-driven RIS and other methods based on the data-driven RWI, and facilitates the acquisition of the dominant modes of the experimental sea area from the data.

The presented method has potential applications in seafloor exploration, RWI extraction, and active sonar dereverberation. Additionally, the proposed algorithm may facilitate the detection of strong moving targets by extracting sparse components and it may be a promising approach for detecting weak signals in the future.