Subspace-based compressive sensing algorithm for raypath separation in a shallow-water waveguide

Multipath acoustic waves propagate in shallow water due to the reflection and (or) refraction of acoustic eigenrays. In ocean acoustic tomography(OAT)Munk1995 , we use the variation in the arrival time of rays to inversely derive fluctuations of sound speed profile in the water column. Recent works have emphasized the use of angle variations to image sound speed perturbations using the sensitivity kernel approach AulanierShallow2013 ; AulanierTime2013 . Multiple paths cover the different parts of the ocean, so they can provide spatially-extended information for the inversion process. However, they also produce interferences between them. Thus, successful separation of these raypaths based on direction-of-arrival (DOA) estimation is the first step of OAT. Researchers have paid great attention to solving this problem, and many algorithms have been proposedschimidt1986 ; jiang2017raypath ; jiang2017Active . Conventional beamforming (CBF) and Capon BF are two such conventional algorithms, but they fail in separating raypaths that arrive close to each other at the receiver array. Moreover, Iturbe and co-authors proposed a double-beamforming algorithm in a double-array configuration to improve the performance of the conventional beamforming algorithm Iturbe2008 . On the other hand, high-resolution algorithms are also introduced to improve the angular resolution of beamforming. The MUltiple SIgnal Classification (MUSIC) algorithm schimidt1986 is a representative of the high-resolution algorithms, which is based on the orthogonality of the signal and noise subspace. Following the strategy of high-resolution algorithms, exploring the spectral information of the emitted signal, Jiang et al proposed a smoothing-MUSICAL algorithm jiang2017raypath to separate the correlated or coherent raypaths and a higher-order algorithmjiang2017Active in the frequency domain to enable the separation the raypaths in an ocean environment in the presence of colored noise. Experimental results illustrate that both algorithms achieve resolution improvement compared to beamforming-like algorithms. However, the noise and signal subspaces are not always totally orthogonal in actual operationthomas1995 . Thus, the DOA estimation based on subspace decomposition still faces great challenges,especially when multiple raypaths propagate in a strongly noisy ocean environment. Recently, the DOA estimation based on CS has been very attractive to many researchers in this area. CS theory as a new information technology makes use of the fact that received signals with inherent sparsity can be reconstructed perfectly from very few measurements through the solution of a convex minimization problem. In underwater acoustics, Edelmann and Gaumond Edelmann2011Beamforming first applied a CS algorithm to estimate DOAs using ship-towed horizontal arrays. Compared with conventional beamforming, the algorithm achieved a resolution improvement and interference suppression. Xenaki and Gerstoft Xenaki2014Compressive developed a compressive beamforming algorithm and further discussed its performance, especially with coherent arrivals, single-snapshot data and random array geometries. In addition, the authors indicated that the limitations of CS-based algorithms are related to the beampattern. Gerstoft et al. gerstoft2015multiple extended the algorithm to the case of multiple snapshots combining a maximum a posteriori method with the least absolute shrinkage and selection operator (LASSO) method. The experiment shows that this method also outperforms conventional high-resolution methods even under challenging scenarios such as coherent arrivals and single-snapshot data. Recently, a compressive sensing beamforming method (CSB-II) zhong2013compressive is developed based on sampling covariance matrix in the application context of aeroacoustic experiments containing strong background noise at broadband frequencies. The experiment results illustrate that the CSB-II method achieves a higher resolution and is robust to sensing noise for presumably spatially sparse and incoherent signals. However, in the context of ocean acoustic tomography, the multiple raypaths are produced by the emitted signal through the refraction or the reflection of the ocean surface and the continental shelf. Thus, they are always fully correlated or coherent. The scattering at different points on the sea surface or the eigenpaths crossing different areas of the water column may attenuate these correlations; however, there is still the possibility that the case of fully correlated or coherent raypaths exists. This is a very difficult separation case. In addition, even in the case of partly correlated raypaths, the correlations between the raypaths impede the successful separation or decrease the separation resolution of raypaths. However, the existing CS-based algorithms do not focus on solving this problem.

In order to perform the separation between eigenrays, the present paper proposes a subspace-based compressive sensing algorithm that formulates the statistic information of the signal subspace in a CS framework. The signal subspace is first obtained through an eigenvalues decomposition to the spectral matrix of the received signal. Then, to suppress the interferences produced by the correlations between the raypaths, the sparse signal reconstruction is achieved by designing a subspace-based $L_{1}$ regularization optimization.

The paper is organized as follows. In Section II, the signal model is established. In Section III, the principle of the comparative method–reweighed compressive sensing is described. In Section IV, we introduce our proposed algorithm. In Section V, we illustrate the performance of the proposed algorithm based on simulations, a tank experiment and an ocean experiment.

In a shallow-water waveguide, we assume that $P$ raypaths emitted by a far-field source $e(t)$ are received by each element of a vertical array composed of $M$ sensors. $s_{p}(t)=a_{p}e(t)$ denotes the $p^{th}$ raypath, where $a_{p}$ is the random amplitude of the $p^{th}$ path. The received signal on the $m^{th}~{}(m=1,2,...,M)$ sensor, denoted by $y_{m}(t)$, can be modeled in the time domain as follows.

where $s_{p}$ denotes the amplitude of the $p^{th}$ source on the $m^{th}$ sensor; and $n_{m}(t)$ is the white Gaussian noise received at the $m^{th}$ sensor. $\tau_{mp}$ represents the propagation delay of the $p^{th}$ raypath between the $m^{th}$ sensor and the reference sensor.

In the frequency domain, Eq. 1 can be rewritten as

with $G_{mp}=exp(-j\nu(m-1)2\pi\frac{dsin\theta_{p}}{c})$, $d$ is the distance between two adjacent sensors and $c$ is the propagation velocity of the acoustic signal, which is assumed to be constant in this paper. The term $S_{p}(\nu)=a_{p}e_{p}(\nu)$, where $a_{p}$ is the random amplitude of the $p^{th}$ path and $e_{p}(\nu)$ is the deterministic Fourier Transform of the emitted signal $e(t)$. ${N_{m}}(\nu)$ is the received white Gaussian noise at frequency $\nu$. and $\theta_{p}$ is the direction of arrival of the $p^{th}$ raypath at the reference sensor.

Considering the received signals of all the $M$ sensors and constructing the sparse signal, Eq. 2 can be rewritten in matrix form as follows.

where $\mathbf{Y}(\nu)=[Y_{1}(\nu),Y_{2}(\nu),\ldots,Y_{M}(\nu)]^{T}$ and $\mathbf{Y}(\nu)$ represents the received signal on all the $M$ sensors, whose size is $M\times 1$. $\mathbf{N}(\nu)=[n_{1}(\nu),n_{2}(\nu),\ldots,n_{M}(\nu)]^{T}$ is a vector of dimension $M\times 1$, which is obtained by the concatenation of the noise at each sensor. $\mathbf{G}=[\mathbf{G}_{1},\mathbf{G}_{2},\ldots,\mathbf{G}_{M}]^{T}$, and $\mathbf{G}_{m}=[e^{-2i\pi\nu\tau_{m1}},\ldots,e^{-2i\pi\nu\tau_{mQ}}]$, which is an over-complete dictionary. $\mathbf{G}$ contains the terms $e^{-2i\pi\nu\tau_{mp}},p=1,2,\ldots,P$, which describe the transfer functions between the source and the sensors. $\mathbf{S}=[S_{1},S_{2},\cdots,S_{P},\cdots,S_{Q}]^{T}$ denotes the sparse-signal vector containing the angle and the amplitude information of the raypaths, which are produced from the emitted signal. $Q$, (${P<M},P<<Q$ and $M<Q$) is the dimension of the sparse-signal vector, and $T$ means “transposed”. The so-called sparsity means that the number of nonzero entries in $\bf S$ is pretty small, i.e., ${\parallel\bf S\parallel}_{0}<Q$.

Based on the signal model in frequency domain described above, we first review the compressive sensing beamforming algorithm. The signal $\bf S$ is sparse, as it has only $P$ nonzero components, with $P<<Q$. Thus, the seemingly underdetermined problem is considered. The sparse vector estimate $\mathbf{S}$ contains the angle and amplitude information of the raypaths, and it can be estimated by the following $l_{1}$ minimization problem:

where the noise is assumed to be spatially white and $\varepsilon$ is an upper bound for the noise norm. , such that $\|\mathbf{N}\|_{2}\leq\epsilon$.

The reweighed CS Xenaki2014Compressive improves the performance of the conventional CS algorithm through the use of a diagonal weight matrix ${\bf W}_{nn}\in{\mathbb{R}}$ to enhance the sparsity of the solution:

The weight matrix, $\mathbf{W}=\textstyle{diag}[w_{1},w_{2},\cdots,w_{Q}]$, is initialized with the identity matrix $\mathbf{I}_{Q}$, and the elements $w_{n}^{k+1}$ at the $(k+1)^{th}$ iteration are updated as

where $\widehat{\mathbf{S}}_{n}^{k}=\textstyle{diag}[\hat{s}_{n}^{1},\hat{s}_{n}^{2},\cdots,\hat{s}_{n}^{k}]$ is the solution in the previous iteration. Setting parameter $\xi>0$ ensures that a null element in the current estimate, $\widehat{\mathbf{S}}_{n}^{k}$, prevents instability in the iteration. The algorithm will proceed until solutions in adjacent iterations finally reach approximately equal values; i.e., $\widehat{\mathbf{S}}_{n}^{k+1}=\widehat{\mathbf{S}}_{n}^{k}$, where

which results in its convergence.

In addition, for stationary sources, which do not move in space, and for the whole $L$ snapshots, which is a single observation vector of an array data, in the time domain, due to the consistency of the sparsity the problem is solved by minimizing the $l_{1}$-norm of the product of the weight $\bf W$ and the vector ${\bf S}_{l_{1}}$ combined with the $l_{1}$-norm of the row vectors in $\bf S$,

where the sparse signal $\bf S$ is a $Q\times L$ matrices and the received data $\bf Y$ has dimensions $M\times L$. The equation is transformed into the equation (9),(10).

where $(s_{q}^{k})_{l_{1}}$ is the $q^{th}$ element of the vectors ${\bf S}_{l_{1}}$ at the $k^{th}$ iteration.

In this paper, the broadband acoustic pulse emission combined to a wavelength that is small compared to water depth allows us to describe the waveguide propagation with the ray model. In the case of tank data, the propagation is dominated by reflections on the top and bottom waveguide boundaries while both refraction and reflections in shallow waters are observed with the ocean data used in the present study (Fig. 3c). This means that the $P$ largest eigenvalues of the spectral covariance matrix can be associated to the $P$ main eigenrays, which may not be the case in oceanic waveguides where modal dispersion is dominating. We first compute the spectral matrix of the received signals as follows.

where $\hat{\mathbf{R}}\in\mathbb{R}^{M\times M}$. $E\{\cdot\}$ denotes the expectation. The expectation is computed by a frequency-smoothing operation jiang2017raypath . $\cdot^{\textbf{H}}$ denotes the conjugate transpose. In addition, the multiple raypaths are produced by the same source, and there exist correlations between the raypaths, which inevitably lead to rank deficiency, so we use frequency smoothing for preprocessing. For the wideband signal, two-sided correlation (TSC) matrices are designed at each frequency. Each transforms the spectral matrix at a certain frequency into the focusing spectral matrix at central frequency Valaee1995Wideband . Then, we apply an eigenvalue decomposition to the spectral matrix $\hat{\mathbf{R}}$. That is,

where $\lambda_{1},\lambda_{2},\cdots,\lambda_{M}$ are the eigenvalues of $\hat{\mathbf{R}}$. The eigenvalues are sequentially arranged as $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{M}$, and $\mathbf{\mu}_{1},\cdots,\mathbf{\mu}_{M}$ denote the corresponding eigenvectors; Therefore, the signal subspace is projected by the largest $P$ eigenvalues (We assume that the number of raypaths $P$ is known in this paperJiang2014Automatic .) and the corresponding eigenvectors.

$\hat{\mathbf{R}_{s}}$ is rearranged into a long vector $\hat{\mathbf{r}_{V}}\in\mathbb{R}^{M^{2}\times 1}$ by orderly stacking each row of $\hat{\mathbf{R}_{s}}$ behind the previous ones. That is, $\hat{\mathbf{r}}_{V}=[R_{s,1,1},R_{s,1,M},R_{s,M,1}\cdots,R_{s,M,M}]^{T}$, where $R_{s,i,j}$ denotes the $i^{th}$ row $j^{th}$ column element of $\hat{\mathbf{R}_{s}}$.

According to Eq. 3, we have

Assuming that the signals and the noise are uncorrelated, we obtain

and the signal subspace is denoted as follows:

In the context of OAT, the raypaths are produced by the same source. They may be de-correlated when they scatter at different points on the sea surface or cross different areas of the water column; however, there is still the possibility that the case of correlated or coherent raypaths exists. To further improve separation resolution, it is necessary to consider and suppress the interferences produced by the correlations between the raypaths. Thus, the above equation can be rewritten as follows:

where

and

Stacking each row of $\mathbf{D}_{1}$ below the previous ones, we rearrange $\mathbf{D}_{1}$ into a long vector $\mathbf{d_{1V}}$ and further reshape it as follows:

where ${\mathbf{G}}^{\prime}=\begin{bmatrix}G_{11}G_{11}^{H}&G_{12}G_{12}^{H}&\cdots&G_{1N}G_{1N}^{H}\\ G_{11}G_{21}^{H}&G_{12}G_{22}^{H}&\cdots&G_{1N}G_{2N}^{H}\\ &&\vdots&\\ G_{M1}G_{M1}^{H}&G_{M2}G_{M2}^{H}&\cdots&G_{MN}G_{MN}^{H}\end{bmatrix}\in\mathbb{C}^{M^{2}\times N}$ and
$\mathbf{S^{{}^{\prime}}}=\begin{bmatrix}S_{1}S_{1}^{H}&S_{2}S_{2}^{H}&\cdots&S_{N}S_{N}^{H}\end{bmatrix}^{T}\in\mathbb{R}^{N\times 1}$. $\cdot^{T}$ denotes the matrix transpose.

Similarly, we rearrange $\mathbf{D}_{2}$ into a long vector $\mathbf{{d}_{2V}}$ by stacking each row of $\mathbf{D}_{2}$ below the previous ones, and $\mathbf{{d}_{2V}}$ will be set as a variable in the calculation with the convex optimization tools.

Finally, we obtain

Based on CS theory, an accurate reconstruction is assured under two conditions: sparsity of the underlying signal and sufficient incoherence of the process that maps the underlying signals to the observations. In our case, $Q>>P$ assures sparsity, as there are only $P$ nonzero values. The incoherence of the process is assured, since the first column to the $M^{th}$ column of the sensing matrix ${\mathbf{G}}^{\prime}$ constitute a Vandermonde matrix (The element $a_{ij}$ at $i^{th}$ row and the $j^{th}$ column of a Vandermonde matrix can be expressed as $a_{ij}=(a_{i})^{j-1}$.). If any two $P$ columns extracted from the sensing matrix are orthogonal, then the Restricted Isometry Property (RIP) condition is satisfied. At this time, the rank of the sub-matrix formed by the extracted arbitrary two $P$ columns is $2P$. As a Vandermonde matrix satisfies the orthogonality of any two $P$ columns mentioned above under certain $P$Jafarpour2011Deterministic ; Cohen2009Compressed , Vandermonde matrixusually satisfies the RIP condition of compressed sensing.

Thus, the sparse signal can then be reconstructed by solving the following $L_{1}$ regularization optimization.

We solved the above programming using available convex optimization tools, such as CVX Chen2001Atomic . The peak positions of $\hat{\mathbf{S^{{}^{\prime}}}}$ represent the DOA estimations of the raypaths.

In this section, we illustrate the performance of the subspace-based CS algorithm with simulation data in a multipath environment, a set of real data obtained in an ultrasonic waveguideRoux2011Travel and, finally, a set of ocean shallow water dataRoux2008The . These tests are all performed in the same configuration, which is composed of a point source and a vertical array of receivers. We used $150$, $300$ and $60$ snapshots in the simulations, tank experiments and at-sea experiments, respectively. The sound speed is uniform in the experiments, $1500$ m/s for the simulation, $1473$ m/s for tank data, $1509$ m/s for the real data. We denote the wavelength as $\lambda$, the height of the water column as $H$, the interval between two adjacent sensors as $d$ and the distance between the point source and the reference sensor as $D$. Their specific values corresponding to each test are shown in Table 1. In order to test the robustness of the subspace-based CS algorithm to the variation of sensor number and the interval between the two adjacent sensors, we chose vastly different values of $M$, and $d/\lambda$ for the three cases.

MUSIC and reweighed CS beamforming Xenaki2014Compressive are provided as comparative algorithms. Fig. 1 shows the received signals in simulations, tank experiments and at-sea experiments. Figs. 2, 3 and 4 illustrate the separation contrast results with simulation data in a multipath environment, a set of real data obtained in an ultrasonic waveguide and, finally, a set of ocean shallow water data, respectively. The SNR is defined as the ratio of signal power to the noise power in the frequency band of the signal for the array. It is denoted by the equation $SNR=10*lg(P_{s}/P_{n})$, where $P_{s}$ and $P_{n}$ represent the power of signal and noise in the frequency band of the signal, respectively.

Using different line types, Fig. 2(a) shows the separation results with simulation data when the SNR is equal to -5 dB, and the black crosses denote the theoretical positions. The theoretical positions indicate the correct arrival angles (or arrival times), which are computed using ray theory Cornuelle1987Compressive ; Skarsoulis2004Travel . It is clear that the subspace-based CS algorithm can successfully separate all the five raypaths, while the MUSIC algorithm totally fails in separating them. In addition, the reweighed CS algorithm detects four peaks, but it is difficult to identify them by their expected positions. When the SNR is equal to 0 dB, five raypaths are successfully separated with the subspace-based CS algorithm, and the corresponding results are shown in Fig. 2(b). In addition, the MUSIC algorithm does not resolve any raypath. The reweighed CS algorithm provides four peaks, and each of them is located between two theoretical values. Thus, it is also difficult to know the corresponding theoretical values of these peaks and whether these peaks are pseudo-peaks. Moreover, we test the performance of the subspace-based CS algorithm at different SNRs and compare its separation root-mean-square errors (RMSE) with those of other algorithms. We define the average RMSE for the direction of arrival ${\hat{\theta}_{p}}$ for the raypath $p$ by the quantity:

Where $K_{i}$ is the number of the trials, which is equal to 20 for the simulation. That is, the simulations are performed 20 times for different SNRs. Generally, the subspace-based CS algorithm is more robust to noise, as shown in Figs. 2(a) and (b). In Fig. 2(c), even for these raypaths, which can be separated by all three algorithms, the subspace-based CS algorithm achieves the minimum RMSEs.

Fig. 4 shows the separation contrast with the data obtained in an ultrasonic waveguide immersed in a water tank. If the frequency of the signals is multiplied by a factor and the spatial distances, including both the one between the source and the receivers and the one between the adjacent receivers, are divided by the same factor, then the physical phenomena occurring in the environment remain the same. Thus, the small-scale experiment can reproduce the actual physical phenomena occurring in nature at a smaller scale inside the laboratory, which achieves a reduced cost and a totally controlled experiment. In the present work, a steel bar acts as the bottom, for which the boundary conditions are nearly perfect at the water-bottom interface, and a  1.10 m long, 5.4 cm deep acoustic waveguide is constructed Roux2014Inverting . Two coplanar 64-element vertical line arrays (VLAs) are placed, and A broadband 1-$\mu$s-long ultrasonic pulse is emitted by each source successively at the 3 MHz central frequency of the transducers Roux2011Travel . The transducer dimensions $0.75~{}mm\times 12~{}mm$ are used to make the linear arrays omnidirectional in the plane defined by the source-receiver arrays, and the beams are collimated in the plane perpendicular to the waveguide axis. We chose five point-to-array configurations (one emitted source and 31 receivers ) and we described their specific positions for the benchmark configuration in Table I. Five eigenrays are expected between the source-receiver arrays with this experimental configuration. To further test the performance of the subspace-based algorithm, we study statistically the outputs of the beamforming algorithms using a set of sources around the actual source [the two above and the two below in the source array (SA)] and sets of subarrays among the vertical received array that are close to the center of the receiver array (RA). The parameters of the benchmark configuration are listed in the second line of Table I. In addition, the noise (SNR=-3 dB) has been added to the received signal for the ultrasound scenario. The same processing using the subspace-based CS algorithm, the CS beamforming algorithm and the MUSIC algorithm was performed to produce the beamformer outputs. In Fig. 3(a), the average beamformer outputs are shown. The purple line corresponds to the average output of the CS beamforming algorithm; the blue line with circles corresponds to that of the subspace-based CS algorithm; and the yellow line with very short vertical lines indicates that of the MUSIC algorithm. Clearly, the subspace-based CS algorithm correctly separates all the five raypaths with minimal bias compared to the theoretical position. At the same time, the MUSIC algorithm fails in separating three of the five raypaths, and the CS beamforming algorithm gives additional pseudo-peaks. In Fig. 3(b), only the peaks of these spectra are retained. The purple triangles denote the peaks of the CS beamforming algorithm; the blue circles represent those of the subspace-based CS algorithm; and the yellow squares correspond to those of the MUSIC algorithm. It can be seen that the subspace-based algorithm shows the smallest dispersion to the theoretical values of the DOAs, which means that it is more robust to slight changes in the experimental configuration than the other two algorithms.

Fig. 3 shows the contrasts of separation results with five data sets obtained in an at-sea experiment Roux2013Analyzing . The FAF05 (Focused Acoustic Field 2005) experiment was conducted in July 2005, to the north of Elba Island, Italy with repetitive data collection over more than eight consecutive hours between two source-receiver vertical arrays that were separated by a distance 4.071 km in a 123 m deep waveguide. The FAF05 at-sea experimental setup is similar to that of the small-scale ultrasonic experiment, although at a much larger scale. The configuration is composed of two equally spaced VLAs. The source array (SA) has 29 transducers spanning 78 m in 120 m water, and the receiver array (RA) has 32 hydrophones covering 62 m. The central frequency of the transducers is 3.2 kHz, with 1 kHz bandwidth. The source signals were 200-ms linear frequency modulated chirps that were compressed after reception to their pulse equivalent by cross-correlation. This provided 40 dB of signal-to-noise ratio for reception with power-limited transmission. Similar to the case of the tests with tank data, we use five different subarrays that will provide different experimental configurations with arrival angles that are close to each other. Fig. 4(c) shows the five eigenrays (including one surface-reflected ray and four refracted rays) propagating between two elements of the SA and the RARoux2013Analyzing . The eigenrays are either refracted at thermoclines or reflected at the air-water interface. All of the eigenrays that interact with the bottom and the surface-reflected eigenray have lower amplitudes than the refracted eigenrays because of reflection loss. The FAF05 data set has been used to develop a coordinated source-receiver array processing procedure and illustrate its effectiveness in an example of tracking raylike arrivals in a fluctuating ocean environment Roux2008The . In addition, using these data, Roux et el. Roux2013Analyzing analyzed the time-dependent sound speed fluctuations in shallow water from group-velocity versus phase-velocity data representation. The subspace-based CS algorithm, the CS beamforming algorithm and the MUSIC algorithm are applied to the five sets of data. We calculate the expected DOAs based on the propagation of the raypaths. The average spectra for the five sets of data, shown in Fig. 4(a), are used to statistically study the performance of each algorithm. The subspace-based CS algorithm provides five peaks around the expected theoretical DOAs, while the CS beamforming algorithm gives multiple peaks with relatively large deviation. In addition, we only keep the peaks for each algorithm and plot them in Fig. 4(b). It can be seen that there is no dispersion in the results of the subspace-based CS algorithms and that there is much more dispersion with those of the CS beamforming algorithm. Additionally, although the results of the MUSIC algorithm gather together around the expected values, it fails in detecting some of the raypaths for some subarray configurations. Based on these results, we can conclude that the subspace-based CS algorithm provides more reliable results with higher angular resolution.

In this paper, a subspace-based CS algorithm is proposed for separating raypaths, which introduces the subspace-based algorithm into a CS framework. Compared to MUSIC and CS beamforming, the proposed algorithm achieved a better performance in terms of resolution and robustness to ocean fluctuations, especially when in a strongly noisy environment.