Sonar Model Based Matched Field Signal Processing

1. Introduction

Passive underwater sonar array systems were subject to development in the past decades for accomplishing tasks in naval surveillance and geo-acoustical investigation of the ocean bottom. Modern sonar passive array systems are using adaptive matched field signal processing techniques. Passive sonar array signal processing based on ocean physical propagation model is known as “Matched Field Processing” (MFP). It is a specific application of the more general case of space-time adaptive signal processing in the area of sonar array signals processing (Klemm, 2002). It is a cross-correlation technique developed for matching the values of the computed with a propagation model sound pressure with the measured values at the output of a passive sonar vertical array. The maximum in the cross-correlation or ambiguity surface gives an estimation of the position of the underwater source. This technique is first proposed by Clay in 1966 and later on theoretically and experimentally developed by Hinich and Bucker in 1972 and 1975 (Sullivan & Middleton, 1993). In 1980 and 1981 Klemm introduced in the matched field signal processing algorithms the Capon’s minimum variance and Burg’s maximum entropy estimators for better resolution in shallow water environments (Sullivan & Middleton, 1993). The acoustic field in an ocean horizontally stratified waveguide may be presented as a sum of “normal waves” or “normal modes”

standing waves

, each satisfying the boundary conditions and propagating down the layer with its own speed, wave constant and attenuation. The coherent modal sum is used in matched field signal processing for localization of an underwater source in range and depth. It can be used for localization of the source also with a horizontal array because the angular modal spread is proportional to the direction of arrival relative to the axis of the array (Klemm, 2002). The matched field processing for source position localization is presented in most of the investigations as a parameter estimation problem. It is assumed that the detection has been carried out and the frequency of the source is known. Some authors consider joint detection and localization problem [Booth et al., 2000]. In the beginning MFP of narrowband signals was used which was later modified with the more realistic for practical applications broadband case where MFP results for several frequencies are averaged. The aim of this chapter is to introduce the reader to the principles of matched field signal processing and present narrowband MFP experimental results. An available on the web normal mode propagation code and vertical sonar array signals data from a controlled source localization experiment are used. A comparison is made between MFP results of three known matched field processors applied to vertical array sonar experimental data.

2. Model based matched field processing of vertical sonar array signals

In conventional antenna array systems signal processing the usual methods of beamforming include phase control of the signals in receive or transmit mode in order to form and steer antenna beam pattern in given direction (see Van Trees, 2002). In MFP the steering vector is determined by the predicted with a propagation model sound pressure for a range of source coordinates in a waveguide with known depth sound speed profile (SSP) and bottom acoustical parameters. A flow diagram of MFP is given on fig. 1.

Figure 1.

Matched Field Processing Flow Diagram for Source Position Estimation

From the time domain array sampled signals the cross-spectral density matrix (CSDM) is estimated. As a steering vector for the matched filed processor the calculated sound pressure level with normal mode propagation model is used. An ambiguity surface is formed through matching or cross-correlating of the steering vector with the received sound pressure field. The maximum of the ambiguity surface gives the most likely position of the sound source in range and depth. There are several general signal processing algorithms applicable for MFP.

2.1. Cross-spectral density matrix estimation

Cross spectral density matrix

The name covariance matrix is also used

(CSDM) is formed through estimation of the cross-spectral functions of the signal at the output of the sonar array. The cross-spectral function may be defined with a Furier transform of the cross-correlation function or directly via finite Furier transform of the input time domain signal (Bendat, 2010). If sample records are available from a stationary random processes in the time domain

For example signals at the output of two hydrophones from the array

$$\left\{x_k\left(t\right)\right\}$$,$$\left\{y_k\left(t\right)\right\}$$ for a finite time interval $$0\le t\le T$$the cross-spectal function of the two processes is defined:

$$S_{xy}\left(\omega ,T,k\right)=\frac{1}{T}{X}_k^{*}\left(\omega ,T\right)Y_k\left(\omega ,T\right)$$E1

where

$$X_k\left(\omega ,T\right)={\displaystyle \underset{0}{\overset{T}{\int }}{x}_k\left(t\right)e^{-j\omega t}dt}$$E2

$$Y_k\left(\omega ,T\right)={\displaystyle \underset{0}{\overset{T}{\int }}{y}_k\left(t\right)e^{-j\omega t}dt}$$E3

The functions$$X_k\left(\omega ,T\right)$$, $$Y_k\left(\omega ,T\right)$$represent finite Furier transforms of the time domain samples $$x_k\left(t\right)$$$$y_k\left(t\right)$$, (*) is the complex conjugate operator, k is sample index. An estimate of the cross spectral density function is

$$S_{xy}\left(\omega ,k\right)=\underset{T\to \infty }{\lim }E\left[S_{XY}\left(\omega ,T,k\right)\right]$$E4

In practise the record length is always finite and the expectation operation $$E\left[\right]$$is taken from a finite number of ensemple elements. If the signals at the outputs of the M elements antenna array are stationary the data can be presented after the finite Furier transforms in the frequency domain with the following matrix:

$$\begin{array}{ccccc}{X}_1\left({\omega }_1\right)& X_1\left({\omega }_2\right)& X_1\left({\omega }_3\right)& \cdots & X_1\left({\omega }_N\right)\\ X_2\left({\omega }_1\right)& X_2\left({\omega }_2\right)& X_2\left({\omega }_3\right)& \cdots & X_2\left({\omega }_N\right)\\ X_3\left({\omega }_1\right)& X_3\left({\omega }_2\right)& X_3\left({\omega }_3\right)& \cdots & X_3\left({\omega }_N\right)\\ ⋮& ⋮& ⋮& \cdots & ⋮\\ X_M\left({\omega }_1\right)& X_M\left({\omega }_2\right)& X_M\left({\omega }_3\right)& \cdots & X_M\left({\omega }_N\right)\end{array}$$E5

The row in (5) represents the frequency spectrum of the signal from a given hydrophone (channel).

The array cross-spectral matrix for the n-th frequency bin is:

$$\begin{array}{ccccc}{X}_1\left({\omega }_n\right)X_1^{\ast }\left({\omega }_n\right)& X_1\left({\omega }_n\right)X_2^{\ast }\left({\omega }_n\right)& X_1\left({\omega }_n\right)X_3^{\ast }\left({\omega }_n\right)& \cdots & X_1\left({\omega }_n\right)X_M^{\ast }\left({\omega }_n\right)\\ X_2\left({\omega }_n\right)X_1^{\ast }\left({\omega }_n\right)& X_2\left({\omega }_n\right)X_2^{\ast }\left({\omega }_n\right)& X_2\left({\omega }_n\right)X_3^{\ast }\left({\omega }_n\right)& \cdots & X_2\left({\omega }_n\right)X_M^{\ast }\left({\omega }_n\right)\\ X_3\left({\omega }_n\right)X_1^{\ast }\left({\omega }_n\right)& X_3\left({\omega }_n\right)X_2^{\ast }\left({\omega }_n\right)& X_3\left({\omega }_n\right)X_3^{\ast }\left({\omega }_n\right)& \cdots & X_3\left({\omega }_n\right)X_M^{\ast }\left({\omega }_n\right)\\ ⋮& ⋮& ⋮& \cdots & ⋮\\ X_M\left({\omega }_n\right)X_1^{\ast }\left({\omega }_n\right)& X_M\left({\omega }_n\right)X_2^{\ast }\left({\omega }_n\right)& X_M\left({\omega }_n\right)X_3^{\ast }\left({\omega }_n\right)& \cdots & X_M\left({\omega }_n\right)X_M^{\ast }\left({\omega }_n\right)\end{array}$$E6

Where (*) is the complex conjugate operator. The cross-spectral matrix for an M element array is square MxM matrix for a given frequency$${\omega }_n$$. The estimation is defined with the following expectation operator for the frequency bin$${\omega }_n$$:

$${\widehat{S}}_X\left({\omega }_n\right)={\rm E}\left\{X\left({\omega }_n\right)X^{*}\left({\omega }_n\right)\right\}$$E7

An estimate of the CSDM $${\stackrel{}{\stackrel{⌢}{S}}}_X\left({\omega }_n\right)$$is formed through averaging ensemble cross-spectrum data. On figure 2 a data cube of CSDM estimates is presented for three frequency bins.

Figure 2.

Data cube of a real CSDM estimates for different frequency bins.

2.2. Steering vector calculation

The frequency range and the parameters of the ocean waveguide determine the sound propagation model for the matched field processing steering vector. In the application a low frequency source positioned in a shallow water ocean waveguide is considered which specifies the normal mode propagation model. For sound propagation numerical calculations in a range-independant environment the normal mode representation of the acoustical pressure actuated by a harmonic source in a horizontally stratified medium is given by the following infinite radical sum (Ferla et al., 1993):

$$p(r,z)=\frac{i}{4\rho (z_s)\sqrt{8\pi r}}{\displaystyle \sum _{m=1}^{\infty }{\Psi }_m(z_s){\Psi }_m(z)H_0^1\left(k_{m}r\right)}$$E8

Where r is the range, z the depth, ρ the density, z_s the source depth, Ψ_m the mode amplitude, $$H_0^1$$the zero-order Hankel function of the first kind, k_m is the eigenvalue. In the solution a large argument asimptotic approximation of the Hankel function is used. When the number of modes is limited to M the final expression for the sound pressure is:

$$p(r,z)=\frac{j}{\rho (z_s)\sqrt{8\pi r}}{e}^{-j\frac{\pi }{4}}{\displaystyle \sum _{m=1}^M{\Psi }_m(z_s){\Psi }_m(z)\frac{e^{j{k}_{m}r}}{\sqrt{k_m}}}$$E9

,

where dependence of the pressure – p from distance – r and depth z is obvious.

Each of the modes has unique wave number, attenuation, phase speed and group speed. The sum of the modes for different source positions and distances from the source forms unique vertical distribution of the pressure which is subject to matching with the experimentally measured pressure (see fig.7). The calculated sound pressure is used as a steering vector in the matched field processor.

Figure 3.

Graphics of nine normal modes representing the sound field in a shallow water waveguide considered in the examples.

Figure 4.

Data cube of the calculated sound pressure level

The steering vector is the sound pressure level for a given possible position of the known source in depth. This forms a sound pressure level data cube presented on fig. 4. For every possible position of the source in depth $$S_1,S_2,S_3,\dots ,S_{end}$$ the sound pressure level is calculated in a 2D region. For every distance from the source $$R_1,R_2,R_3,\dots ,R_{\max }$$a depth sound pressure vector $$P_1,P_2,P_3,\dots ,P_m$$is formed for the known depth hydrophone receiver positions$$D_1,D_2,D_3,\dots ,D_{\max }$$. This vector is later used as a steering vector for correlation with the estimated CSM and forming the ambiguity surface. A maximum in the ambiguity surface gives an estimation of the source position in range and depth. On figure 5 a data cube with real pressure fields respectively for sources: S1, S10 and S68 considered in chapter 3.

Figure 5.

Data cube of the calculated sound pressure fields

2.3. Matched field signal processors

Two matched field processors are subject to investigation in the application. The first one is the classical or Bartlett matched field processor and the second one is the so called minimum variance distortion less processor (MVDR).

2.3.1. Bartlett (conventional) matched field signal processor

In conventional beamforming the sonar array pattern is steered through phase control of the signals at the output of the antenna array (see fig. 1) through steering vector. In matched field processing the steering vector represents the computed pressure in the points of the array elements. The processor output or as in our case a matched field filter output is:

$$Y\left(\omega \right)=v^H\left(\omega ,r,z\right)X\left(\omega \right)$$E10

Where $$v^H\left(\omega ,r,z\right)$$is the pressure vector, $$X\left(\omega \right)$$is the signal at the input of the procerssor and $$Y\left(\omega \right)$$ is the output signal. We want to estimate the power of the signal at the output of the filter:

$$P\left(\omega ,r,z\right)=E\left[v^H\left(\omega ,r,z\right)X{X}^{*}v\left(\omega ,r,z\right)\right]=v^H\left(\omega ,r,z\right)S_X\left(\omega \right)v\left(\omega ,r,z\right)$$E11

where $$S_X\left(\omega \right)$$ corresponds to CSDM given with (7). This conventional processor is known as a Bartlett processor. Bartlett matched field processor forms ambiguity surface as a direct dot product of model replica vectors $$\vartheta \left(\omega ,r,z\right)$$and the sample CSDM -$$S_X\left(\omega \right)$$. The maximum of the power estimate in the ambiguity surface gives the probable position of the source.

$$P_{BART}=v^H\left(\omega ,r,z\right)S_X\left(\omega \right)v\left(\omega ,r,z\right)$$E12

This processor is a benchmark for the other more sophisticated processors because it has much power and low probability of no detection.

2.3.2. MVDR matched field signal processor

The MVDR processor is derived with quadratic constraints (Van Trees, 2002). An optimal filter is given with the block diagram on fig. 6. $$X\left(\omega \right)$$and$$Y\left(\omega \right)$$represent the signals at the input and at the output of the filter in the frequency domain. $$W^H\left(\omega \right)$$is the filter frequency response.

Figure 6.

General optimal filter (processor) structure.

The signal vector at the input can be written as:

$$X_S\left(\omega \right)=F\left(\omega \right)v\left(\omega ,r,z\right)$$E13

,

where $$F\left(\omega \right)$$ is the frequency-domain snapshot of the source signal and $$v\left(\omega ,r,z\right)$$ is the ocean waveguide manifold vector. It is required that in the absence of noise$$Y\left(\omega \right)=F\left(\omega \right)$$. The constraint of no distortion

Or no gain

implies $$W^H\left(\omega \right)v\left(\omega ,r,z\right)=1$$ in order to receive the input signal at the output of the processor without distortion (with unity gain). Minimization of the mean square of the output noise leads to MVDR beamformer first derived by Capon. Optimal weight of the processor is given (Van Trees, 2002):

$$W_{MVDR}^H\left(\omega \right)=\frac{v^H\left(\omega ,r,z\right)S_N^{-1}\left(\omega \right)}{v^H\left(\omega ,r,z\right)S_N^{-1}\left(\omega \right)v\left(\omega ,r,z\right)}$$E14

where $$S_N^{-1}\left(\omega \right)$$ is the inverse of noise CSDM for a given angular frequency, $$v^H\left(\omega ,k_S\right)$$is the hermitian transpose operator. When instead of $$S_N^{}\left(\omega \right)$$ the signal CSDM $$S_X^{}\left(\omega \right)$$is used in (14) the estimator is known as minimum power distortionsless response (MPDR) (Van Trees, 2002). Usually it is more difficult to estimate $$S_n^{}\left(\omega \right)$$ to form the weight, because it is colored with multipath ocean waveguide interference. It is not possible however also to estimate a pure signal matrix $$S_X^{}\left(\omega \right)$$ beacause it is mixed with noise. The optimum weight of MPDR

In the majority of the literature it is also referred as MVDR (Van Trees, 2002)

has the form:

$$W_{MPDR}^H\left(\omega \right)=\frac{v^H\left(\omega ,r,z\right)S_X^{-1}\left(\omega \right)}{v^H\left(\omega ,r,z\right)S_X^{-1}\left(\omega \right)v\left(\omega ,r,z\right)}$$E15

Later in the text we are using the name MVDR for weights (14), (15) and specify if the CSDM of the signal or noise is used. This estimator (beamformer) in the literature is known also as Minimum Variance Estimator (MVE) and Maximum Likehood Estimator (Klemm, 2002). The power of MFP-MVDR is given with (Van Trees, 2002):

$$P_{MVDR}=W_{}^H\left(\omega \right)S_X\left(\omega \right)W_{}\left(\omega \right)$$E16

where $$S_X\left(\omega \right)$$ is the signal CSDM and $$W_{}\left(\omega \right)$$ is the optimal weight defined with (14) or (15). When we have acceptable for array processing signal to noise ratio (SNR) it is relatively easy to estimate$$S_X\left(\omega \right)$$ and use (15). Subject to optimization in the described processors is the match between the function of the computed sound pressure over the array elements and the measured output signal power for a given frequency, given with CSDM. On fig. 7 the computed normalized values of the pressure and the measured signal power are superimposed at distance and depth of good match.

Figure 7.

Predicted and measured signal power (diagonal values of CSDM) for a frequency bin of interrest over the array at distance and depth of good match (red – measured signal power, blue – computed squared pressure).

One way to build robust to some degree of mismatch weight vector and improve the numerical stability in matrix inversion is through diagonal loading of the CSDM. This is the so called method of “diagonal loading”

Van Trees shows that imposition of quadratic constrain leads to diagonal loading

(Van Trees, 2002). It is also referred in the literature as a “white-noise-constraint” (WNC). Adding noise to diagonal elements expands the noise space and eliminates small eigenvalues in the inverse CSDM. Usually the CSDM is not well conditioned matrix

The analysis of the experimental averaged CSDM showed that the ratio between largest singular value to the smallest is more than 1 000 000

and the diagonal loading is a way toward more stable matrix inversion procedure. The MVDR weight with diagonal loading of the CSDM is (Van Trees, 2002):

$$W_{MVDRDL}^{}\left(\omega \right)=\frac{v^H\left(\omega ,r,z\right){\left[S_X\left(\omega \right)+{\sigma }^{2}I\right]}^{-1}\left(\omega \right)}{v^H\left(\omega ,r,z\right){\left[S_X\left(\omega \right)+{\sigma }^{2}I\right]}^{-1}v\left(\omega ,r,z\right)}$$E17

Another way to overcome problems with CSDM inversion is to apply singular value decomposition (Van Trees, 2002). With singular value decomposition the square CSDM is presented in the following form:

$$S_X=U\Lambda V^{*}$$E18

where Λ is diagonal matrix of the eigen values:

$$\Lambda =\left|\begin{array}{ccccc}{\lambda }_1& & & & \\ & {\lambda }_2& & & \\ & & .& & \\ & & & .& \\ & & & & {\lambda }_M\end{array}\right|$$E19

where $${\lambda }_1>{\lambda }_2>\mathrm{...}>{\lambda }_M$$ and U and V* are unitary with size MxM. Choosing the rank “r” of the decomposed CSDM directly imposes constraints on the spatial signal subspace and is powerful instrument for multipath interference rejection (r is smaller than the number of sensors M which in our case is 48). The inverse matrix is formed with the following matrix multiplication:

$$S^{-1}{}_X=\left|\begin{array}{cccc}{U}_{11}& .& .& U_{1r}\\ .& .& & \\ .& & .& \\ U_{M1}& & & U_{Mr}\end{array}\right|\left|\begin{array}{ccc}1/{\lambda }_1& & 0\\ & .& \\ 0& & 1/{\lambda }_r\end{array}\right|\left|\begin{array}{cccc}{V}_{11}& .& .& V_{1M}\\ .& .& & \\ .& & .& \\ V_{r1}& & & V_{rM}\end{array}\right|$$E20

The optimal weight for MVDR is formed:

$$W_{MVDRPINV}\left(\omega \right)=\frac{v^H\left(\omega ,r,z\right)\left(U{\Lambda }^{-1}{V}^{*}+{\epsilon }^{2}I\right)}{v^H\left(\omega ,r,z\right)\left(U{\Lambda }^{-1}{V}^{*}+{\epsilon }^{2}I\right)v\left(\omega ,r,z\right)}$$E21

where the inverse matrix

Also known in the literature as pseudo-inverse

is formed with rank reduction and diagonal loading is applied (adding noise to the inverse matrix through diagonal loading), controlled by the parameter ε. A similar form of this processor weight was used in (Kolev & Georgiev, 2007).

The rank in the application is determined if the corresponding singular value is under 1/100 from the maximum value. On fig. 8 the singular values of a CSDM are plotted and the rank is determined to be r=5.

Figure 8.

Typical eigenvalues of CSDM determining the matrix rang – r=5.

3. Vertical sonar array MFP experimental results

Passive array sonar data from SACLANTC 1993 North Elba experiment available on Internet (Saclant sonar data, 1996) is used for processing. Additional information about the experiments and signal processing examples is given in (Gingras & Gerstoft,1995, Gerstoft & Gingrass, 1996). The array data was collected in shallow-water off the Italian west coast by the NATO SACLANT Center in La Spezia, Italy. SACLANT Center has made this data available to the public for the purposes of fostering signal processing research. The original SACLANT time series has been converted to a series of MATLAB files each of which contains a matrix data file from the 48 sensors. Each file represents about 1 minute of data The vertical array consists of 48 hydrophones with spacing 2 m between elements at total aperture length 94 m (18.7 m to 112.7 m in depth). On 27-th of October a source was towed in the surveillance area at a speed approximately 3.5 kn. The source emitted PRN signal with a center frequency of 170 Hz. For this source 10 data files are available. The data for the moving source was subject to MFP. The available on the site sound speed profiles were used for replica computations with KRAKEN normal mode program (Porter, 2007). Technical manual of the software was used from (Porter, 1992). On fig. 9 a datagram of sampled time domain array signals and their spectrum is given.

The CSDM was formed from averaged FFT data. The FFT length was 4096 resulting in 0.2441 Hz FFT binwidth for 1000 Hz sampling data rate. A Hamming window with 50 % overlapping was used for weighting. It was observed from fig. 11 the relatively low level of the power of CSDM in the frequency window of interest with center frequency 170 Hz compared to the low frequency noise (50 – 100 Hz). For estimation of signal CSDM the frequency bin for 169.9 Hz is used.

Figure 9.

A datagram of time domain array sampled multichannel signals and their power spectrum.

Figure 10.

Measured Sound Speed Profile of SCALANTC 1993 North Elba experiment.

Figure 11.

CSDM for 1 minute moving source array signal.

On fig. 12 the ambiguity surfaces of three matched field processors: Bartlet, MVDR with diagonal loading and MVDR with reduced rank SVD and pseudoiverse matrix are plotted. The position of the source is determined with the maximum in the ambiguity surface. All data was equally treated with rank of the reduced matrix r = 5, and diagonal loading constant ε² = σ² = 10. Comparison of the three processors side lobes for an equal depth cut of the ambiguity surfaces is given on fig. 13. The level of the side lobes is lower when introducing spatial constraints with rank reduction of the inverse matrix. In the source localization results on fig. 12 the Bartlett processor is used as a benchmark. Even it has more side lobes in the ambiguity surface there is global maximum corresponding to the source position. During the processing of the available array data it was established that the accuracy of passive localization with the narrow band MFP is very sensitive to environment parameters change. Little inaccuracy of SSP or bottom geo-acoustical parameters leads to substantial errors in localization. The images on fig. 12 are for the moving source data recorded on the 27-th of October. The available environment file on the web site is used in the computations (Saclant Sonar Data, 1996). If in the environment file the SSP is changed with the data from another file - “sspprofiles.txt” for the 27-th of October more than one and halve thousand meters offset in the source range occurs.

Figure 12.

Comparison of ambiguity surfaces localization results of the three processors

Figure 13.

Comparison of the normalized power of the three processors MFP-Bartlet - blue, MFP-MVDRDL - green and MFPPINV - red.

4. Conclusion

The basics of matched field signal processing were presented. Three matched field processors were applied on real experimental vertical array sonar data for source localization purposes. The MFP-MVDR processor has a lower side lobe levels especially applied with rank reduction and pseudo-inverse matrix diagonal loading (see fig. 13). For the investigated ranges and shallow water waveguides however the sound structure in range is reverberant and repetitive resulting in multiple ambiguity surface local maximums or false alarms. The current state of the technology leads to development of underwater sensor networks with application of bistatic and multistatic sonar concepts. Results from measurements with a bistatic sonar experimental setup are presented in (Kolev et al. 2009), (Kolev et al. 2010). The interest for developing adaptive beamforming methods is as with uniform but also with non-uniform and random sparse (ad-hoc) arrays with distance between nodes bigger than a half-wavelength (Hodgkiss, 1981, Gerstoft, P & Hodgkiss, 2011). MFP methods are subject to development for application in these underwater networks.

Acknowledgments

The authors are grateful to SACLANT centre scientists who have made experimental sonar array data and normal mode propagation codes available on the web.