Matched-field localization using a virtual time-reversal processing method in shallow water

Time-reversal processing (TRP) is an implementation of matched-field processing (MFP) where the ocean itself is used to construct the replica field. This paper introduces virtual time-reversal processing (VTRP) that is implemented electronically at a receiver array and simulates the kind of processing that would be done by an actual TRP during the reciprocal propagation stage. MFP is a forward propagation process, while VTRP is a back-propagation process, which exploits the properties of reciprocity and superposition and is realized by weighting the replica surface with the complex conjugate of the data received on the corresponding element, followed by summation of the processed received data. The number of parabolic equation computational grids of VTRP is much smaller than that of MFP in a range-dependent waveguide. As a result, the localization surface of VTRP can be formed faster than its MFP counterpart in a range-dependent waveguide. The performance of VTRP for source localization is validated through numerical simulations and data from the Mediterranean Sea.

Time-reversal processing (TRP) [1][2][3][4][5][6][7] is a process of retransmitting a received signal in a time-reversed fashion by a source-receiver array. In the frequency domain, time reversal is equivalent to phase conjugation. When TRP is applied to source localization, known as virtual time-reversal processing (VTRP), it is unnecessary to send a signal back and forth between the source and receiver. Instead, assuming that the acoustic channel is sufficiently stable in time, the retransmission of the temporal dispersed signals in a time reversed fashion will be done by a computer. A passive array is used in VTRP instead of a source-receiver array.
When TRP is applied to source localization, there is one drawback to this back-propagation approach. Because of sound attenuation, the acoustic intensity close to the "source-receiver" array is much stronger than the focusing peak at the real source location. This means that the focusing peak is only a local maximum around the real source *Corresponding author (email: walternwpu@gmail.com) location rather than a global one. To suppress the stronger peaks close to the passive array and make the focusing peak the global maximum, an appropriate normalization factor should be introduced to the back-propagation field. Therefore, the localization surface formed by VTRP contains the focusing peak, which provides the greatest likelihood of the source location.
Matched-field processing (MFP) and VTRP are conceptually similar, but differ in their physical meaning and physical implementation. MFP [8][9][10][11][12][13] is a forward propagation process that consists of systematically placing a test point source at each point of a search grid, computing the replica vectors on the array and then correlating these replicas with the data from the real source. A search is performed in a region of possible target positions. VTRP is a back-propagation process which exploits the properties of medium reciprocity and superposition. VTRP can be realized by weighting the replica surface with the complex conjugate of the data received on the corresponding element, followed by summation of the processed received data.
The normal mode (NM) method [14] provides an accurate and computationally efficient propagation model to solve range-independent ocean acoustic problems. Although the computational orders of MFP and VTRP are different, their total numbers of NM computational grids are the same in a range-independent waveguide. However, many realistic environments cannot be adequately described by rangeindependent models. The parabolic equation (PE) method [14] is effective for solving non-adiabatic range-dependent ocean acoustic problems. To compute the replica vector at the array corresponding to a test point source at a range search grid, the acoustic field from the source and the array must first be computed. Thus, the number of PE computational grids of VTRP is much smaller than that of MFP in a range-dependent waveguide. As a result, the localization surface of VTRP can be formed much faster than its counterpart of MFP in a range-dependent waveguide.
1 Theoretical model

Matched field processing
MFP deals with target localization by matching the data of the target radiated acoustic field, acquired by an array to a model-based replica vector at a test target position. The replica vector on the array for each candidate position (r, z) is normalized to the unit norm and is denoted as w MFP (r, z) here. For a Bartlett matched field processor, the ambiguity function (or surface) is given by [12]  where d(r s , z s ) is a data vector observed on the array for the real source at (r s , z s ), and R is the data covariance matrix. Superscript ( ) H denotes the Hermitian or conjugate transpose of the matrix.

Virtual time-reversal processing
The VTRP is an implementation of the back-propagation process which is done not by the ocean itself as a TRP but by a computer. We define a vector of search depths as z=[z 1 ,z 2 ,···z ND ] T and a vector of search ranges as r=[r 1 ,r 2 ,···r NR ] T , where ND and NR are the numbers of depth grids and range grids. Superscript ( ) T denotes the transpose of a matrix or vector. The localization surface S VTRP (r,z) can propagate outward from the array toward the potential target location(s) in a search region (r, z): where d(r i ,d i ;r s ,z s ) represents the received acoustic pressure of the ith element at (r i , z i ) propagated from the source at (r s , z s ). w VTRP (r,z;r i ,z i ) is the normalized replica surface generated by the virtual source corresponding to the ith element at (r i , z i ). N is the element number. Superscript ( ) * denotes the complex conjugate. The superposition of these N surfaces yields the final localization surface whose peak value provides the greatest likelihood that the target is present. The amplitude and phase of the signal fluctuate with time even when stringent efforts are made to fix all controllable parameters. The effects of such variability can be reduced by considering the mean behavior of the fields rather than individual samples of the fields themselves. When using VTRP to process the real data, the acoustic pressure d(r i ,z i ;r s ,z s ) used for back-propagation can be obtained by following the procedure described in [5]. The signal matrix X is constructed by gathering multiple signal vectors (or snapshots). Using singular value decomposition [15,16], we can obtain where U=[u 1 ,···,u k ,···,u K ] is a N×K matrix whose columns u k are left singular vectors, Σ=diag(σ 1 ,···σ k ,···σ K ) is a K×K matrix whose diagonal elements are singular values, and V=[v 1 ,···v k ,···v K ] is a N×K matrix whose columns v k are right singular vectors. The data vector d(r s ,z s ) can be obtained from the combination of the left singular vectors with the corresponding singular values Now, the acoustic pressure d(r j ,z j ;r s ,z s ) used in eq. (2) corresponds to the ith element of the data vector d(r s ,z s ).

Comparison between MFP and VTRP
(i) Range-independent waveguide. In most cases, the application of MFP has been restricted to range-independent waveguides. The NM method is convenient because it is accurate and computationally efficient for low frequency and far field in range-independent waveguides. According to NM theory, the acoustic pressure field can be expressed in terms of an NM expansion. Once the geo-acoustic parameters are defined and the source frequency is given, the mode shape functions and horizontal wavenumbers can be calculated numerically by a single evaluation of the propagation model. The total acoustic pressure field is then the weighted sum of the contributions from each mode.
The physical implementations of MFP and VTRP are different as shown in Figure 1. MFP is a forward propagation process and a test point source is systematically placed at each point on a search grid. Then the N-dimensional replica vector RI MFP P at the array is computed as shown in Figure 1(a). The corresponding number of NM computational grids is N. For a search region of ND×NR grids, we will run ND×NR times to calculate all these replica vectors. Thus, the total number of NM computational grids is ND×NR×N.
VTRP is a back-propagation process. It is assumed that there are N virtual sources located at each of the N element positions on the array. Each virtual source generates a replica surface RI VTRP P which corresponds to a search region of ND×NR grids as shown in Figure 1(b). The corresponding number of NM computational grids is ND×NR. Therefore, the total number of NM computational grids is N×ND×NR.
Although the computational orders of MFP and VTRP are different, their total numbers of NM computational grids are the same in a range-independent waveguide.
(ii) Range-dependent waveguide. For real shallow water waveguides, the supposition that the depth, the sound speed profile and the bottom are constants is usually too coarse. In a majority of cases, the variations in the waveguide parameters are very significant. Variations in depth in the coastal regions are particularly marked. The PE method is very effective for solving non-adiabatic range-dependent ocean acoustic problems.
The PE computation range step and depth step are dr and dz. The maximum computational depth is Z max . A search region over the source range and depth is performed from Z 1 to Z 2 in depth and from R 1 to R 2 in range (m). The range-depth search grid is ΔR×ΔZ. Then, the number of range search grids is NR=(R 2 -R 1 )/ΔR+1, and the number of depth search grids is NR=(Z 2 -Z 1 )/ΔZ+1. Let the largest search range R 2 be the 1st range search grid, and the nearest search range R 1 be the NRth range search grid. Figure 2(a) shows the implementation procedure of MFP. To compute the replica vector RD MFP P at the array corresponding to a test point source at the kth range search grid, the acoustic field before the kth range search grid must first be computed. The number of PE computational grids will be [ ] There are ND test point sources at each range search grid. Thus, the number of PE computational grids corresponding to the kth range search grid is ND×C k . As a result, the total number of PE computational grids for all the NR range search grids is region of ND×NR grids. The total number of PE computational grids corresponding to the N replica surfaces is The ratio of the number of PE computational grids between MFP and VTRP is Because the largest search range R 2 is often hundreds of times larger than ΔR, the number of PE computational grids of VTRP is much smaller than that of MFP. For example, if R 1 =1000 m, R 2 =10000 m, ΔR=30 m, Z 1 =5 m, Z 2 =115 m, ΔZ=2 m, and the element number N=60, then C MFP /C VTRP =155. As a result, the localization surface of VTRP can be formed faster than its counterpart of MFP in a range-dependent waveguide. Figure 3 illustrates the range-dependent shallow-water channel used for the simulations. A range-dependent water layer overlies a constant thickness sediment layer and sub-bottom half-space layer. The values of bathymetry at the receiver and the source are 140 and 130 m. The speed of sound in the water-column is modeled with a summer profile having an almost isovelocity layer extending to 60 m and a strong thermocline spanning 60-80 m depth. The ocean surface is treated as a flat pressure-release surface.

Numerical simulation 2.1 Simulation conditions
It is assumed that the whole environment is time invariant. Thus the sound channel is reciprocal. For simplicity, the ambient noise is ignored. A search region over a source range and depth is performed from 5 to 115 m in depth and from 1000 to 10000 m in range. The range-depth search grid is 30 m×2 m. The localization surface is normalized so that the peak value is 0 dB. The location of the maximum (0 dB) on the localization surface is used as the source location parameter estimate. The acoustic field used in this study is calculated by RAM [17][18][19].
To quantify the localization performance, the output signal to interference noise ratio (SINR) [20] is used. In the following simulations, once the maximum of the ambiguity surface is within the neighborhood (±500 m in range and ±10 m in depth) of the real source location, then the localization is correct, and vice versa. The SINR is defined as the maximum output signal (normalized to 0 dB) minus the 75th percentile of the normalized localization surface powers after sorting in ascending order (100th percentile means the peak of the localization surface). Once the localization failed, the SINR is meaningless, and is left at 0 dB.

Simulation results
The point source is narrow-band at a frequency of 170 Hz and locates at 5000 m in range and 75 m in depth. The VLA consists of 60 elements spanning the water column from 20 to 138 m with 2 m inter element spacing. Figure 4(a) shows the localization surface formed by MFP. The focused peak illustrates that the source localization is correct, and its output SINR is 15.3 dB. Figure 4(b) shows the case for VTRP. The source is correctly localized, and its output SINR is also 15.3 dB. Note that the localization performances of MFP and VTRP are identical.
The average computer processing times of MFP and VTRP are about 1365 and 8 min, respectively. These times are reported for a personal computer (model number: p6515cn; Intel(R) Core(TM) i3 540@3.07 GHz). Using this computer, VTRP proceeds about 170 times faster than MFP. The reason is that the number of PE computational grids of VTRP is much smaller than that of MFP in a range-dependent waveguide.

Experiment description
The approach proposed in this study is also evaluated using vertical array data collected during the October 1993 Mediterranean Sea Trials [21,22]. The environment model and experimental configurations are similar to Figure 3  spacing. The source signal was pseudo-random noise in a 20 Hz band around 170 Hz, whose −3 dB bandwidth was approximately 12 Hz, and the sampling frequency was 1 kHz. Because the narrow-band problem was of interest here, only the frequency bin near 169.9 Hz was used to obtain a narrowband snapshot of the sensor outputs. A detailed description of the experimental dataset may be found in [22]. The data are also available at http://spib.rice.edu/spib/saclant.html.

Stationary source localization
The source level of stationary source was approximately 163 dB re: 1 μPa Hz . Based on the known uncertainties of the GPS position for the vertical array and the source buoy, the source range with respect to the vertical array was predicted to be (5600±200) m. The accuracy of the knowledge about the source depth leads to a prediction of (80±2) m. The input signal to noise ratio was high, at about 10 dB. The data at each hydrophone was transformed into the frequency domain using fast Fourier transform. Each snapshot had 1024 data points. Both the data covariance matrix R for MFP and the signal matrix X for VTRP are averaged over 60 snapshots.
Ambiguity surfaces were computed using the range-dependent model parameters from [22]. Figure 5 shows the localization surfaces of the MFP and VTRP using the first minute of data. The estimated source range and depth are (r, z) MFP =(5530 m, 76 m), and (r, z) VTRP =(5530 m, 76 m). This is consistent with the estimated value (5560 m, 77 m) [22]. All of them are able to estimate the source position correctly. The output SINRs are 10.6 dB and 10.9 dB for the MFP and VTRP. The output SINR of VTRP is a little higher than that of MFP. The reason is that the data used for VTRP is obtained through singular value decomposition of the signal matrix. The average CPU times of MFP and VTRP are 1305 and 6 min. As the number of PE computational grids for VTRP is much smaller than that for MFP, VTRP proceeds about 217 times faster than MFP.

Conclusions
In the forward propagation MFP, the replica vectors corresponding to the search grids should be calculated one at a time. In this paper, the concept of VTRP for source localization is introduced. It is a back-propagation process which exploits the properties of medium reciprocity and superposition. Though the computational orders of MFP and VTRP are different, their total numbers of NM computational grids are the same in a range-independent waveguide. However, many realistic environments cannot be adequately described by range-independent models. The number of PE computational grids of VTRP is much smaller (hundreds of times) than that of MFP in a range-dependent waveguide.
Compared with MFP, VRTP can achieve the same localization performance while using much less CPU time in a range-dependent waveguide. Furthermore, VTRP is also verified by previously published experimental data, and similar results are obtained.