Fast raypath separation based on low-rank matrix approximation in a shallow-water waveguide

In shallow water, acoustic waves propagate in multi-paths due to refraction and reflection at the ocean surface and the continental shelf. Multi-path propagation has been widely used in many applications, such as active sonar, ocean acoustic tomography,¹ and source localization, because more raypaths provide more information. However, the successful use of the information is based on the premise of first eliminating the raypaths' interferences with each other produced by the multi-path propagation. Among the existing algorithms to separate raypaths, subspace-based algorithms are noticeable for their large resolution improvement. In particular, the multiple signal classification (MUSIC) algorithm,² as a representative subspace-based method, was applied to separate raypaths under the assumption of an uncorrelated narrow-band signal. However, MUSIC fails to separate raypaths when they are fully correlated or coherent. To overcome this drawback, Jiang et al.^3,4 developed a high-resolution method called the smoothing-multiple signal classification active large band (MUSICAL) algorithm, which greatly improves the separation performance, especially when two fully correlated raypaths arrive close together. Although the smoothing-MUSICAL algorithm was developed under the assumption of white Gaussian noise, colored noise exists in real ocean environments. Recently, Jiang et al.⁵ presented a higher-order wideband active algorithm to separate raypaths interrupted by colored noise. Naturally, it takes advantage of fourth-order cumulants containing non-Gaussian information in non-Gaussian processes. The experimental results obtained from tank data at small scale⁵ illustrate that the higher-order algorithm can obtain more accurate results and suppress noise; however, it has large computational cost (both time and memory space) for computing the singular value decomposition of a huge cumulants matrix. This cost is prohibitive from a practical perspective. To overcome this drawback, we propose in this study a fast higher-order separation algorithm, which has the advantages of low computational cost and high resolution. In addition, its performance will be tested using synthetic data, experimental data obtained in a tank and ocean data. The most important contribution is that the proposed algorithm provides a strategy to find the signal subspace efficiently, which can be extended to all the subspace-based algorithms. Based on the strategy, their computational cost can be largely reduced especially when the cumulants matrix has large dimensions.

This paper is organized as follows. Our fast raypath separation algorithm is introduced in Sec. 2. In Secs. 3 and 4 we discuss the performance and the error of the proposed method. In Sec. 5, its performance test using synthetic data, experimental data obtained in a tank, and ocean data is conducted, respectively. We draw conclusions in Sec. 6.

We consider the case where P raypaths arrive on a vertical array composed of M sensors. The received signal of the mth $(m=1,2,…,M)$ sensor in the frequency domain is modeled by

where the term $e(ν)$ is the deterministic amplitude of the emitted signal at frequency ν. The amplitude of each raypath, a_p, is considered random, where T_p represents the arrival time of the $pth$ raypath on the reference sensor and $tm(θp)$ is the delay between the reference sensor and the $mth$ sensor. Additionally, $tm(θp)=d sin(θp)/c$⁠, where θ_p is the arrival direction of raypath p on the array, d is the distance between two adjacent sensors, c is the propagation velocity of the acoustic signal, and $nm(ν)$ is the additive noise received at the $mth$ sensor.

We use F frequency bins of the wideband signal. Hence, Eq. (1) is written in matrix form as follows:

where $xg=[x+(ν1),x+(ν2),…,x+(νF)]+$ with $x(νi)=[x1(νi),x2(νi),…,xM(νi)]+$$(i=1,…,F)$⁠; that is, $xg$ is a long-vector of size M × F, which is arranged according to all the frequency bins on all the sensors. $ng$ is similarly the concatenation of the observed noise vectors at each frequency. Specifically, $ng=[n+(ν1),n+(ν2),…,n+(νF)]+$ with $n(νi)=[n1(νi),n2(νi),…,nM(νi)]+$⁠. $A=[a1,a2,…,aP]+$ is a vector of dimension P. $H=[h1,h2,…,hP]+$ with $hp=[e(νi)e−2jπν1τ1p,…,e(νF)e−2jπνFτMp]+$⁠, where $e−2jπνiτmp(i=1,…,F;m=1,…,M)$ is the transfer function between the source and the mth sensor, and + indicates transposition.

We first apply spatial-frequency smoothing^4,5 as a preprocessing step to eliminate the interferences among raypaths. That is, for a spatial smoothing factor K_s, an antenna composed of M sensors is divided into $2Ks+1$ sub-antennas. The $ks$th sub-antenna is composed of the sensors in the range of $[ks,ks+1,…,M−2Ks]$⁠. At the same time, a similar frequency-smoothing process is also applied, where K_f denotes the frequency-smoothing factor. Finally, using the recurrences produced by the spatial-frequency smoothing, we can estimate the wideband trispectral matrix with the following formula:

where $⊗$ represents the Kronecker product, $*$ indicates conjugation, and H indicates the conjugate transpose. Based on the preprocessing, the rank of $Ĉ$ is equal to $K=(2Ks+1)(2Kf+1)$ (it is necessary to select K greater than P).

The traditional subspace-based algorithms use methods such as the singular value decomposition (SVD) to span a signal subspace, however, they require memory and time which are super-linear in the size of the trispectral matrix (or spectral matrix). It is prohibitive for a practical application, especially when the sizes are very large. Inspired by Frieze et al.,^6–8 we introduce a randomized strategy to span the signal subspace at low computation cost.

The strategy is based on the following three principles.

• Principle 1: We can approximate the product of two matrices by an algorithm called the basic matrix multiplication.⁷ That is, when we input two matrices, $Ĉ∈ℝ(MF)2×(MF)2$ and $ĈT$⁠, a probability distribution ${qi}i=1(MF)2$ and a number $c≤(MF)2$⁠, two output matrices $Ĉ1$ and $Ĉ1T$ are returned, such that $Ĉ1Ĉ1T≈ĈĈT$⁠, where $Ĉ1∈ℝ(MF)2×c$ is a matrix whose columns are c randomly chosen columns of $Ĉ$ with suitably rescale.

• Principle 2: An appropriate choice of the probability distribution is crucial steps of the proposed algorithm. Although a uniform distribution could be always used, the optimal probabilities ${qi}i=1MF2$ are proved and defined as the following equation:⁷ Suppose $Ĉ∈ℝ(MF)2×(MF)2, c∈ℤ+$⁠, such that $1≤c≤n$⁠, and ${qi}i=1(MF)2$ are such that $qi≥0,∑i=1nqi=1$⁠, and such that for some positive constant $β≤1$⁠: $qi≥β|Ĉi|2/‖Ĉ‖F2$⁠.

• Principle 3: According to linear algebra, for symmetric matrices the left and right singular vector are the same. The singular values of $Ĉ$ are the nonnegative square roots of the eigenvalues of $ĈTĈ$ and of $ĈĈT$⁠; in addition, the left singular vectors of $Ĉ$ are eigenvectors of $ĈĈT$ and the right singular vectors are eigenvectors of $ĈTĈ$⁠.

The general process of the proposed algorithm is first described: (1) based on principle 3, finding the eigenvectors of the matrix $Ĉ$ is equivalent to find the eigenvector of the matrix $ĈĈT$⁠; (2) based on principles 1 and 2, the matrix $ĈĈT$ can be approximated by $C1̂C1̂T$ with the optimal probabilities, where $Ĉ1$ consists of those c columns of $Ĉ$ (after appropriate rescaling); it has been proven that by choosing appropriate sampling probabilities and enough columns the approximation error can be made arbitrarily small; (3) based on principle 3, finding the eigenvectors of the matrix $C1̂C1̂T$ can be transformed to find the eigenvectors of the matrix $C1̂TC1̂$⁠, then from $Ĉ1T Ĉ1$ approximations to the top singular values and corresponding singular vectors of $Ĉ1 Ĉ1T$ may be computed, which are the approximations of the eigenvectors of the matrix $ĈĈT$⁠. Thus, the fast algorithm are finally operated in four specific steps.

1. To sample c columns of $Ĉ$ (⁠$Ĉ∈ℝ(MF)2×(MF)2$⁠) randomly according to a judicious probability q, we first compute $qi=|Ĉ(i)|2/||Ĉ||F2,i=1,2,…,(MF)2$ for all $(MF)2$ columns (where $Ĉ(i)$ is the ith column of matrix $Ĉ$ and $||·||F$ indicates Frobenius norm). Subsequently, we get a random number r, where $0<r<1$⁠. If $∑j=1j=i−1qj<r<∑j=1j=iqj$⁠, we retain the ith column. The sampling step will be repeated another c – 1 times and finally we will get c columns of $Ĉ$⁠.

2. According to the rule provided by the first step, we sample c columns where $1<P<c<(MF)2$⁠, rescale each column by the factor $1/cqit,t=1,2,…,c$⁠, and arrange all c columns to construct the approximation matrix of $Ĉ$⁠, which is denoted as $C1̂∈ℝ(MF)2×c$⁠. We would like to approximate the eigenvectors of square matrix $Ĉ$⁠. Normally, we can obtain these vectors directly through performing a SVD on matrix $C1̂C1̂T∈ℝ(MF)2×(MF)2$⁠, but its size is still the same as that of $Ĉ$⁠. Thus, we first perform a SVD on $C1̂TC1̂∈ℝc×c$⁠, whose size is quite smaller, to find its eigenvectors, which are also the right singular vectors of $C1̂$⁠; then, they can be transformed into the left singular vectors of $C1̂$⁠, which are also the eigenvectors of the matrix $C1̂C1̂T$⁠. Therefore in the next steps first we carry out the SVD on matrix $C1̂TC1̂$ to obtain right singular vectors of $C1̂$⁠; Second we compute P² left singular vectors of $C1̂$⁠, which successively approximate P² eigenvectors of square matrix $Ĉ.$

3. SVD is performed on matrix $Ĉ1T Ĉ1$⁠, and the P² eigenvector $yt(t=1,2,…,P2)$⁠, corresponding to the P² largest eigenvalues $σt(t=1,2,…,P2)$⁠, is identified.

4. Eigenvectors $ht(t=1,2,…,P2)$ are calculated using the equation $ht=Ĉ1yt/σt(Ĉ1)$⁠.

Based on the eigenvectors computed previously, the signal subspace $Ĉs$ is projected as follows:

where $HP2$ is composed of $ht(t=1,2,…,P2)$⁠.

Finally, the direction of arrival (DOA) and the arrival time of each raypath are obtained by maximizing the following estimator:

where the wideband steering vector $a(θ,T)$ is the Kronecker product of $b(θ,T)$⁠,

where $b(θ,T)$ is the concatenation of the classical steering vectors $d(νi,θ)$⁠, that is,

The error analysis of the low-rank matrix approximation is performed in two steps: (1) define the bound of the approximation error and (2) choose c and define the sampling probabilities.

The low-rank approximation incurs two error components: first, the error between $Ĉ$ and $Ĉ1P2$ (assuming $Ĉ1P2$ is the most optimal P²-rank approximation to matrix $Ĉ$⁠); second, $HP2HP2TĈ$ incurs an error depending on $‖ĈĈT−C1̂C1̂T‖F$ or, equivalently, $‖Ĉ‖F2$⁠. It has been proven⁷ that the error appears entirely within the following two bounds with respect to the Frobenius norm and spectral norm, respectively,

where $‖‖F$ indicates the Frobenius norm and $‖‖2$ indicates the spectral norm. The optimal P²-rank approximation matrix is $C1̂P2=UP2UP2TC1̂$⁠, where $UP2$ is composed of the top P² left singular vectors of $C1̂$⁠.

To control the approximation error, the method used to select c and to define the sampling probabilities is important. According to the theory of Drineas et al.,⁷ we can determine the sampling probabilities as follows. We assume $HP2$ is constructed by sampling c columns of $Ĉ$ with probabilities $qi (i=1,…,(MF)2)$ to make $qi≥β|Ĉi|2/‖Ĉ‖F2$ for some positive $β≤1$⁠. Set $η=1+(8/β) log(1/δ)$⁠. ϵ represents an additional error coefficient and δ denotes the failure probability.

• If a Frobenius norm bound is desired: set $qi=|Ĉ(i)|2/‖Ĉ‖F2, c≥4P2η2/ϵ2$⁠; then, with probability of at least $1−δ$⁠,

  $‖Ĉ−HP2HP2TĈ‖F2≤‖Ĉ−Ĉ1P2|F2+ϵ‖Ĉ‖F2.$

  (10)
• If a spectral norm bound is desired: set $qi=|Ĉ(i)|2/‖Ĉ‖F2, c≥4η2/ϵ2$⁠; then, with probability of at least $1−δ$⁠,

  $‖Ĉ−HP2HP2TĈ‖22≤‖Ĉ−Ĉ1P2|22+ϵ‖Ĉ‖F2.$

  (11)

It is necessary note that the approximation error can be made arbitrarily small by choosing enough columns.

In the proposed algorithm, we can sample one column in one pass over the matrix without previously knowing its size. With the algorithm running in parallel, we can obtain more samples in a single pass in O(c) additional space and time. Given the index to be sampled, we need additional $O((MF)2c)$ space and time to construct matrix $C1̂$⁠. Computing $C1̂TC1̂$ requires $O((MF)2c)$ additional space and $O(((MF)2c)2)$ additional time, and we consider the SVD of $C1̂TC1̂$ to require $O(c2)$ additional space and $O(c3)$ additional time. Then, computing $HP2$ requires P² matrix-vector multiplications for $O((MFP)2c)$ additional space and time. Overall, the proposed fast algorithm requires $O((MF)2c+c2)$ additional space and $O((MF)2c2+c3)$ additional time, which have linear relationships with the number of rows in the high-order cumulants square matrix.

In this section, we test the performance of the proposed algorithm with many sets of simulation data, a set of real data obtained in a tank⁹ (the experiment is based on the principle that if we amplify the frequency of the signals by a factor while reducing the spatial distance by the same factor, the actual physical phenomena occurring in nature can be reproduced on a smaller scale) and finally a set of ocean data.¹⁰ (The experiment was performed in July 2005 north of Elba Island, Italy.) The conventional higher-order algorithm in Ref. 5 is provided for comparison. These tests are performed on a server with 6 Intel Xeon CPUs(X7542) with 48 cores @2.67 GHz and 256.0 GB RAM on the Redhat Linux 5.6 operating system. All the tests are performed between a point source and a vertical array, which are composed of M equally spaced sensors. We show the parameters including the number of sensors (M), the source location (z_s), the range of the array covering the ocean from the position at $zr1$ m to the position at $zrM$⁠, the interval between two adjacent sensors d and the distance between the source and the reference sensor D in each configuration in Table 1. We give the central frequency (ν_c), the bandwidth of the emitted signals (ν_w), the number of the frequencies (⁠$νn:$ from ν_b to ν_e) and the taps of the broadband time-domain signal. (N) used in the tests in Table 2.

The proposed fast algorithm is applied to the simulation data, a set of real data obtained in a tank and finally a set of ocean data and the separation results are shown in Fig. 1 and Fig. 2, respectively. The expected positions are marked with black crosses. Figure 1(A) illustrates the separation result with a set of simulation data. We can see that the five raypaths are all successfully separated with few errors compared with their expected values [(0.5, 0), (0.9, −0.09), (1.2, 0.15), (1.9, −0.254), and (2.5, 0.235)]. Moreover, the efficiency of the proposed algorithm is tested with many sets of simulation data for various numbers of snapshots. Ten realizations are used for each number of snapshots. Figure 1(B) shows that the time reduction ratio increases sharply with increasing number of snapshots. Figure 1(C) shows the separation result with a set of real data obtained in a tank. The proposed algorithm successfully separates the seven raypaths with few errors. Their expected values are (0.6789, 0), (0.6799, 3.1), (0.68, −3.2), (0.683, 6.3), (0.683, −6.3), (0.688, 9.3), and (0.6882, −9.4). Figures 2(A) and 2(B) show the separation results for the at-sea experimental data with the beamforming algorithm and our proposed algorithm, respectively. The proposed fast algorithm are more efficient for noise suppression and obtain the smaller spots representing the positions of the raypaths compared to the beamforming algorithm. Figure 2(C) displays the three raypaths in Figs. 2(A) and 2(B), which propagate between the centers of source and receive arrays. The centers of source and receive arrays both located at 96 m in the water. Due to the sound speed variation described in the Ref. 10, the three raypaths arrived at close time. There are a surface-reflected ray and two refracted rays. Each refracted ray has a turning point near the depth of maximum sound-speed variability. Figures 2(a) and 2(b) illustrate that the proposed fast algorithm is robust to the fluctuation of environmental factors, for instance the water column sound speed variation we have described. The performance of the 4-smoothing-MUSICAL algorithm and the proposed fast algorithm are almost the same, which are not included in the paper due to space limitations. Table 3 compares the time and space consumed by the 4-smoothing-MUSICAL algorithm and the proposed fast algorithm, illustrating that with the proposed algorithm, the cost of the time in decomposition step has been reduced by approximately 99.62% (simulation), 99.7% (small-scale experiment), and 99.9% (ocean data), respectively. The peak memory has been reduced by approximately 95.50% (simulation), 99.2% (small-scale experiment), and 99.75% (ocean data), respectively.

In this paper we proposed a fast raypath separation method that achieves large reductions in time and space consumption compared with the conventional high-order cumulant-based algorithm. In the future, the method will be extended to improve the efficiency of the multi-dimensional subspace-based algorithms.

This paper was supported by the National Natural Science Foundation of China (Grant Nos. 61401085 and 31400842), the State Key Laboratory of Acoustics, Chinese Academy of Sciences (Grant No. SKLA201604), and the Scientific Research Foundation for the Returned Overseas Chinese Scholars. The Focused Acoustic Forecasting experiment (FAF05) experiment was performed in collaborative experiments with the NATO Underwater Research Centre (NURC), La Spezia, Italy, with Mark Stevenson as Chief Scientist. Scientists who contributed to these experiments include Tuncay Akal, W. A. Kuperman, W. H. Hodgkiss, H. C. Song, B. D. Cornuelle, Piero Boni, Piero Guerrini, other NURC staff, and the officers and crew of the RV Alliance.