Active Sonar Detection in Reverberation via Signal Subspace Extraction Algorithm

This paper presents a new algorithm called Signal Subspace Extraction (SSE) for detecting and estimating target echoes in reverberation. The new algorithm can be taken as an extension of the Principal Component Inverse (PCI) and maintains the benefit of PCI algorithm and moreover shows better performance due to a more reasonable reverberation model. In the SSE approach, a best low-rank estimate of a target echo is extracted by decomposing the returns into short duration subintervals and by invoking the Eckart-Young theorem twice. It was assumed that CW is less efficiency in lower Doppler than broadband waveforms in spectrum methods; however, the subspace methods show good performance in detection whatever the respective Doppler is. Hence, the signal emitted by active sonar is CW in the new algorithm which performs well in detection and estimation even when low Doppler is low. Further, a block forward matrix is proposed to extend the algorithm to the sensor array problem. The comparison among the block forwardmatrix, the conventional matrix, and the three-mode array is discussed. Echo separation is also provided by the new algorithm. Examples are presented using both real, active-sonar data and simulated data.


Introduction
A major problem in moving platform active sonar systems is the detection of targets in reverberation. Reverberation is caused mainly by the multiple reflections, diffusions, or diffractions of the transmitted signal by the surface and bottom interfaces. When the target is close to one interface, the target echo is hidden in the reverberation resulting in a low signal-to-reverberation ratio (SRR). Moreover, since the reverberation is strongly correlated with the signal, classical detection methods like matched filtering (MF) are inefficient. In order to improve detection, we can use a model of reverberation, but a correct model is difficult to find because reverberation contains both diffuse components (which look like noise) and also more discrete components (which look like signal). If a global statistical description of reverberation is available, like a reverberation scattering function, the structure of the theoretical optimal receiver is known [1]. Generally, this is not the case and a simplified model is used.
One often used statistical model considers reverberation as nonstationary, colored noise. This approach is used in [2] for monochromatic transmitted signals and in [3] for wideband signals. In [4], they show that algorithms based on this approach have some problems when the Doppler shifts of reverberation and target echo are similar. Further, this model does not take advantage of the connection between reverberation and the transmitted signal.
In [4], they propose to use a simplified model which is deterministic: reverberation is considered as a sum of undesirable echoes. The method for detection consists in estimating these echoes and deleting them before applying the classical MF. It is important to choose a metric to distinguish reverberation echoes from target echoes, and since the target echo power is often lower than reverberation power, they choose echo power as a metric. Next, they need to find an algorithm which is able to separate echoes with different power. They used the Principal Component Inverse (PCI) algorithm which was introduced in [4][5][6].
This algorithm originally assumes that noise is completely different from the searched signal. But [4] shows that PCI algorithm can separate several similar echoes (which means echoes with slight time shift or/and Doppler shift) differing in powers.
PCI is applied to detection in presence of reverberation by taking reverberation as a sum of echoes with higher power than target echoes. PCI algorithm separates the received data into two parts: reverberation and target echoes. By this means we detect targets. However, even when reverberation power is high, there are still some reverberation echoes with lower power. The sum of these lower reverberation echoes sometimes makes a strong confusion with targets.
In this paper, we present a new algorithm named Signal Subspace Extraction (SSE) based on a more real-life model. The SSE algorithm divides the reverberation with target echoes into three parts: higher reverberation echoes, target echoes, and lower reverberation echoes. It makes use of a low rank characteristic of the target echo subspace and separates the signal subspace via the singular value decomposition (SVD) method. PCI separates the reverberation and the target echo by invoking Eckart-Young theorem [7] while SSE extracts the signal by invoking Eckart-Young theorem twice.
Broadband waveforms are generally preferred to continuous wave (CW) in low Doppler [8] when using spectrum methods. However, the subspace algorithms are efficient in whatever the respective Doppler is [4,9]. Hence, in this paper, CW is brought back into use and shows good performance in signal extraction and estimation by experiments with real and simulated reverberation.
In [10], a three-mode array is brought for PCI for sensor array problem, and the detection is improved by the threemode array. However, the three-mode array is a kind of a three dimension matrix. To make it work, the Eckart-Young theorem and SVD have to be extended to a three-dimension problem, too. In this paper, we provide a block forward matrix which is a two-dimension matrix, but this matrix still extends SSE into sensor array problem. SSE with the block forward matrix can be regarded as a new algorithm of Space Time Adaptive Process [11] which jointly processes received data in angle and Doppler to improve the separation of target echo and reverberation. The comparison among block forward matrix, traditional matrix, and three-mode array is also presented. Section 2 quickly reviews the classical detection/estimation hypothesis, the Block Normalized Matched Filter (BNMF). Section 3 presents the SSE algorithm and gives results with adapted real temporal data. Section 4 presents the block forward matrix and extends the algorithm to sensor array problem and discusses the property of the new matrix in comparison with the conventional matrix in [4] for PCI. In Section 4.3, experiments are taken by simulation of space time reverberation in comparison between PCI and SSE. We also give the examples of comparison results among block forward matrix, traditional matrix, and three-mode array.

Detection/Estimation Problem in the Presence of Reverberation
The detection problem is written as follows: where x(t) is the observed or received signal, r(t) is the reverberation noise generated by the transmitted signal, and n(t) represents white noise. The signal emitted by the active sonar is assumed to be a CW. s(t) is the signal to be detected. We assume here that it is linked to the emitted signal e(t) in a simple way: s(t) is differed from the emitted signal by a time delay τ, a Doppler shift f d , and an amplitude attenuation A in the block where signal presents and let x n be the sampled vector of x(t) and let s n be the sampled signal, where T = 1/ f s is the the sampling interval and we sampled at time t = nT: or where f d / f s is named as normalized Doppler Frequency. All signals are complex valued and represent the sonar output after complex demodulation. We work with time-sampled signals.

One-Dimensional Signal Subspace Extraction
We model reverberation as a sum of echoes issued from the transmitted signal which implies that reverberation and the target echoes share almost the same properties.
By cutting x into X i , the forward matrix Y i is generated by where N is the block length and p is chosen close to N/2. It is well known that a vector, which is a linear combination of k complex exponentials, can be made into the forward matrix above, and the matrix will have rank k, if min(p, N − p + 1) ≥ k [5]. However, since for the reverberation echoes k min(p, N − p + 1), reverberation echoes span the full space of matrix Y i . In the context of the real data of interest, if one assumes that the reverberation or target echoes are well approximated by a series of CW suitably scaled, the rank of target echoes subspace is low because for target echoes k ≤ min(p, N − p + 1) in the matrix in most detection cases.
The SSE algorithm consists of decomposing Y i into three matrices Y r1 i , Y s i , and Y r2 i : where Y r i = Y r1 i + Y r2 i is the reverberation plus white noise subspace and Y s i is the received signal dominant subspace. As reverberation power is stronger than received signal in most cases, according to Eckart-Young theorem and [4], Y r1 i is the best r-rank approximation of Y i if r is the rank of dominant reverberation subspace. After we delete the dominant reverberation Y r1 i , the residual matrix contains the target echoes, residual reverberation, and noise, and target echoes become the principal component in the residual matrix. Then we use Eckart-Young theorem for the second time. Y s i is the best s-rank approximation of Y o i if s is the rank of target echo subspace. The result is obtained via the Singular Value Decomposition(SVD) of Y i : where U is the left singular-vector matrix of Y i , V is the right singular-vector matrix of Y i , and Σ is a diagonal matrix which contains the decreasing singular values of Y i , {σ i } (σ 1 > σ 2 > · · · ). Vector X s i is then collected from Y s i . The subspace signal estimation is obtained and then X s i contains only the signal. The detection processing is done on the vector X s i .
The rank used to partition the matrix is not known and must be estimated. In the SSE approximation it is determined by following the method suggested in [13]. This procedure uses the partial sums of squared singular values from the SVD of the data matrix as its test statistic. We start testing from the smallest sum and work our way upwards till, for some I, the partial sum exceeds a specified threshold. The singular values are in descending order, {σ i } (σ 1 > σ 2 > · · · ). We seek the smallest I, I min for which where R Y is the rank of Y i . Following this method, we also seek the largest J, J max for which where Q and P are the SSE threshold values. Q is related to the higher power of reverberation. The sum of Q and P is related to the whole power of reverberation [4]. The rank is then chosen as r = J max and s = R Y − J max − I min . From real cases studied, usually Q is simplified to the higher power of reverberation, since Q P. Hence, the first step of SSE is the same as the PCI procedure including the threshold. And s approximates to the number of target echoes since for target echoes k ≤ min(p, N − p +1) in the matrix, if the transmitted signal is CW, and the target echo is present in the block.
The SSE thresholds used here are based upon the background reverberation power and may be set using prior knowledge or derived from the data. If I min + J max ≥ R Y , SSE does not treat this block. And only when I min + J max ≤ R Y , the SVD is required to determine the signal subspace.
A hypothesis is necessary for a correct running of SSE: the rank r1 of Y r1 i must be small. This hypothesis is the same to PCI and indicates that SSE will fail when SRR is extremely low.

Experiments.
The experiments are performed by comparison. They are based on the real data taken from a sea trial in South China Sea. The transmitted signal is CW. Data is received by active sonar. A moving target presents in this trial. The sampling frequency f s is 5 kHz. The normalized Doppler frequency f d / f s due to the moving target is 0.049. Here we use the normalized Doppler frequency to plot the detection/estimate results. Reverberation is mainly caused by bottom echoes. The BNMF algorithm is applied after PCI and SSE to see the detection improvement. The time series is cut into 0.1 s in each block.
Three experiments are performed with different SRRs. We add weighted adjacent block without target echoes into the block in which target echo is present to obtain data with different SRRs. This is reasonable based on the temporally local stationarity [14]. The results are shown in Figures 1-3. Here we only plot the result of the block in which the target echo is present. We observed that without PCI or SSE, the target echo could hardly be detected. The detection and estimation are improved by PCI and SSE.  Figure 1 shows the results of BNMF outputs without processing and after PCI and SSE with SRR −12 dB. For PCI and SSE, they both detect the target and give right estimation of normalized Doppler frequency f d / f s of the target echo which is 0.049. However, many false alarms appear with PCI while no false alarm is present with SSE in this experiment. Figure 2 shows the results of BNMF outputs without processing and after PCI and SSE with SRR −17 dB. It is observed that the target would not be detected with the PCI if only the largest peak is chosen. The false alarms are the consequence of a lack of lower reverberation removal. Figure 3 shows the results of BNMF outputs without processing and after PCI and SSE with SRR −22 dB. When SRR is lower, the detection becomes worse: both PCI and SSE give false alarms. But the false alarms are still less and lower with SSE than PCI. The target would not be detected with the PCI if only the largest peak is chosen.

Two-Dimensional Signal Subspace Extraction
where x n,m is the received signal on the mth sensor at time sample n. Let us consider a linear array of M sensors with equally interelement spacing d. The signal emitted by the active sonar is CW. Then each element of s n,m is written as where β is the direction of target and λ is the wavelength. e(t) is the emitted signal. τ is the time delay. f d is the Doppler shift. f s is the sampling frequency. Hence f d / f s is the normalized Doppler frequency. The detection and estimation are also performed block by block. The temporal length of each block is N, and the spacial length is M which is equal to the number of sensors. We use the classical generalized likelihood ratio test (GLRT) to build the algorithms. GLRT does not only choose H 0 and H 1 but also estimates the Doppler frequency f d and azimuth β. The estimation of the delay τ is quantified by the shift between two blocks. For the ith block, the statistic test of the BNMF for space time GLRT is This detector requires the estimation of new parameter β. H 1 is chosen if max fd,β L i ( f d , β) > η. η is a given threshold.

Extension of SSE to Senor Array Data.
For the sensor array problem, the SSE algorithm is the same. Only the arrangement of matrix Y i changes. Here we propose a block forward matrix for SSE. The block forward matrix is similar to the block Hankel matrix [15,16]. The block forward matrix of x n,m is defined as where for m = 1, 2, . . . , M, where p is chosen close to M/2 and q is chosen close to N/2. If x n,m is composed of one complex exponential, the block forward matrix has rank one, because each row can be expressed as a complex scale factor times the first row. Matrices (14) and (6) share similar structure: for matrix (6), the shift between two rows or two columns is equal to one sample, and so it is for matrix (14). Hence the rank analysis is the same. The additional degree of freedom given by the spatial dimension leads to easier separation of different echoes.
We are now interested in the separation of two echoes issued from the transmitted signal. As we have known, two different target echoes are represented by different time delays, different directions, or different Doppler frequencies.
The consideration of different time delays of echoes turns to the signal present or not after we divide the received data into short time duration. If time delays of two different echoes are the same, which means that two target echoes appear in the same block, it is obvious that two echoes have to be described by different singular values in order to separate them. We have shown that rank of signal subspace is related to number of target echoes differed by frequencies and directions since for CW the target echo vector is composed of different complex exponentials. Then it is easy to separate different target echoes. Here we also present the matrix derived in [4]. It is built from the data received on all sensors and has the general form: where N is the block length, and M is the number of sensors. Every column corresponds to the output of one sensor. The algorithm described in this section is applied to this matrix. We present the arrangement of the first block to illustrate the structure of this matrix. The matrix is then written as follows: In matrix (17), the shift between two rows is the same and so is between two columns.
The comparison between these two matrices will be presented by experiments in the next section. We analyze them theoretically in this section. First, the dimension of  the block forward matrix is higher than the second matrix. Hence for the full-rank matrix, echoes could be represented by more singular values in block forward matrix. The separation of echoes will be easier. Then, the lengths of the column and the row are nearly the same in block forward matrix but can be quite different in the second matrix in which they completely depend on the number of the sensors and the block length. The number of singular values depends on the shorter one, which means that the data cannot be used in the most efficient way in (16), if the lengths of column and row are different.

Experiments.
Simulations are performed to check the proposed algorithms in array data in this section. We first consider the reverberation containing one target echo with block forward matrix. Then, reverberation with two target echoes will be used to perform the separation. Finally, the comparison of different matrices will be presented.

Space Time Reverberation Model.
Consider a narrowband, M element linear sonar array with a constant intersensor spacing d towed along the x-direction with a velocity v. The complex envelope of the Doppler-shifted reverberation data received at the mth sensor at (x m , y m ) = ((m − 1)d, 0), (x 0 , y 0 ) = (0, 0) at time t n = τ 0 + nT, can be written as [17] r mn = θiφl α θ i , φ l e j(2π/λ) cos φl(sin θi(m−1)d+2v sin θinT) , where T = 1/ f s is the the sampling interval, azimuth −π ≤ θ i < π, 1 ≤ i ≤ M θ , and elevation angle |φ l | ≤ φ max , where φ max is the multipath elevation angle spread defined by the critical angle of the ocean acoustic channel. α(θ i , φ l ) is the complex scatter amplitude from a reverberation patch at range cτ 0 /2, where c is the propagation speed of sound in water. The total number of reverberation patches M θ N φ MN.

Signal Extraction.
In this experiment, suppose that there is one target echo with an SRR of −18 dB, normalized Doppler frequency f d / f s = 0.05, and azimuth β = π/4. The number of sensors is M = 16 and number of time samples is N = 64 in each block. The results of the detector in which block the target is detected are shown in Figure 4. The target is detected after both PCI and SSE. Both show right estimates of the true parameters. However, a lot more false alarms appear with PCI. The superiority of SSE is easily shown.

ROC.
The superiority of the proposed detection scheme is demonstrated from the experiments above. However, to make this claim more precise, we evaluate the experimental performance of the detectors by receiver operating characteristic (ROC) curve where the detection rate is plotted versus the false alarm probability in Figure 6. Monte Carlo simulations were performed comprising 1500 realizations with one target echo present with an SRR of −12 dB and equally many with reverberation and noise only. The number of sensors is M = 16 and number of time samples is N = 16 in each block. The ROC curves are shown in Figure 6 for the detectors using PCI and for the proposed detectors using SSE. Comparing the two curves, we see that SSE has a higher probability of detection when probability of false alarm is low.

Separation.
In the following experiments, the separation performance of SSE will be demonstrated.
In the first experiment of separation, suppose that there are two target echoes in one block. They are one target echo with normalized Doppler frequency 0.05 and azimuth π/4, and the other with normalized Doppler frequency −0.05 and azimuth 3π/8. The amplitude of two target echoes is slightly different with 1 : 0.97 ratio. The SRR is −16 dB. The result of the detector after SSE is shown in Figure 5 In the second experiment of separation, the ratio of amplitude of two target echoes is changed into 1 : 0.5. The result of the detector after SSE is shown in Figure 7(a). When we apply SSE, a few false alarms appear. The echo with higher power is easy to detect, but the less powerful echo is no stronger than the false alarms. However, after the separation performance in Figures 7(b) and 7(c), the two target echoes are well detected and estimated. This step of performance requires the preknowledge of the power level of each target echo.

Matrix Comparison.
We use the conventional matrix in (16) with the same data as in the first experiment for Section 4.3.4 above. Results are shown in Figure 8. False alarms appear with PCI in Figure 8 We also show the result of PCI with the three-mode array [10] in Figure 9. The experiment is arranged in the same condition with Section 4.3.2. Comparing with Figure 4, the detection/estimation capability is equally the same with the block forward matrix using PCI. And since the three-mode array is a three-dimension matrix, SSE is too complicated to be applied to it and so is the echo separation which [10] is not mentioned either. Hence block forward matrix still performs the best among the three matrices in efficiency and detection/estimation capability.

Conclusions
In this paper, we present a new algorithm Signal Subspace Extraction to extract the signal subspace from reverberation. SSE is tested by adapted real signal-channel data and shows good results. Then we derive a block forward matrix and extend the method to the sensor array problem. Experiments by simulations show the block forward matrix works well with the new algorithm not only in detection of target echoes but also in separation of target echoes.

A. Singular Value Decomposition
Given a matrix A m×n whose rank is r and m × n, there exist two orthogonal matrixes U m×m = (u 1 , u 2 , . . . , u n ) and V n×n = (v 1 , v 2 , . . . , v n ):

B. Eckart-Young Theorem
Let the SVD of A be given by (A.1) with r = rank(A) ≤ p = min{m, n} and the singular values are in descending order, {σ i } (σ 1 > σ 2 > · · · ), and define where k < p; then A k is the optimal approximation of A in the view of