Localization of Buried Objects Using Reflected Wide-Band Underwater Acoustic Signals

Abstract

This chapter deals with the localization of wide-band underwater acoustic sources. A combination of high resolution methods with scattering acoustic model are proposed. The bearing and the range sources at each sensor are expressed as a function to those at the first sensor. We present the noneigendecomposition methods fixed-point algorithm, projection approximation subspace tracking (PAST) algorithm, PAST with deflation (PASTD) algorithm and orthogonal PAST (OPAST) algorithm to track the signal subspace to compute leading eigenvectors. The proposed algorithms are faster than singular value decomposition (SVD) for MUSIC. The spatial smoothing operator is used to decorrelate the received signals and to estimate the coherent signal subspace. The performance of the different methods are evaluated by both computer simulations and experimental and data recorded during underwater acoustic experiments.

1. Introduction

Non-invasive detection and localization of sources is an important application area in many application domains, such as radar, sonar, seismology and communications. Thus there has been a growing interest in developing techniques for the estimation wavefronts of the direction-of-arrival (DOA) in order to detect and localize the emitting sources [1]. Support vector machine (SVM) based on electromagnetic approach [2–4] and conventional neural networks (NN) based on inverse scattering technique [5] are proposed for buried object detection. Ground penetrating radar (GPR) is used to improve the detection of weak-scattering plastic mines [6]. But electromagnetic filed inversion require more computational effort. The inversion of measured scattered acoustical waves is used to image buried objects, but it needs high frequencies and the application in a real environment is difficult [7]. Therefore, the acoustic imagery technique is not suitable because the high frequencies are strongly attenuated inside the sediment. Using a low frequency, synthetic aperture sonar (SAS) has been recently applied on partially and shallowly buried cylinders in a sandy seabed [8]. The bearing and the range estimation using correlated signals scattered from nearfield and farfield objects, in a noise environment, still a challenging problem. The MUSIC algorithm is one of the most thoroughly studied and best understood subspace based high resolution methods. It divides the observation space into two signal-subspaces: the signal subspace and the noise subspace [9]. MUSIC uses the orthogonality property between the two areas to locate sources. Different approaches exist to detect and localize buried objects but acoustic techniques will be considered in our study. Match field processing (MFP) [10] has been successfully used for localization sources in ocean acoustic. We discuss the proposed approach based on MUSIC associated with acoustic scattering model referred to MFP [10]. We take into account the water-sediment interface [11]. This means that we attempt to combine both the reflection and refraction of wave in the model [12]. From the exact solution of the acoustic scattered field [13], we have derived a new source steering vector including both the ranges and the bearings of the objects. This source steering vector is employed in objective function instead of the classical plane wave model [14, 15] which have extended the 1-D MUSIC to 2-D MUSIC. The acoustic scatter field model has been addressed in many published researches with different configurations. For example, the configurations can be single [16] or multiple objects [17], buried or partially buried objects [18] with cylindrical [16] or spherical shape [19]. All those scattering models can be used with the proposed source steering vector. In this chapter a spatial smoothing operator is proposed to estimate the coherent signal subspace [20]. Inverse power method, which allows to find an approximate eigenvector when an approximation to corresponding eigenvalue is already known, is proposed to estimate the required noise variance. In high resolution method, we use singular value decomposition (SVD) in music for obtaining the eigenvectors noise subspace. However, the main drawback is the inherent complexity and computational time load [21]. So a large number of approaches have been introduced for fast subspace tracking in order to overcome this difficulty. We propose to replace SVD by Fixed Point for computing leading eigenvectors from the spectral matrix [22, 23]. We propose another methods to accelerate MUSIC, such as projection approximation subspace tracking (PAST) [24, 25], which makes the expectation of square difference between the input vector and the projected vector minimum. With proper projection approximation, the PAST derives a recursive least squares (RLS) algorithm for tracking the signal subspace. The PAST algorithm computes an asymptotically orthogonal basis of the signal subspace. PAST with deflation (PASTD) is derived from PAST by applying the deflation technique in order to get the signal eigenvectors and eigenvalues [24, 26]. It has been shown that these subspace trackers are closely linked to the classical power iterations method, but does not guarantee the orthonormality at each iteration [27, 28]. Orthogonal PAST (OPAST) algorithm is another fast implementation of the power method which outperforms both PAST and PASTD to reduce computation time [29, 30]. The performance of the proposed algorithms are evaluated by several numerical simulations and the data has been recorded using an experimental water tank.

The remainder of the chapter is as follows: Section 2 introduces the problem formulation. Section 3 presents the scattering acoustic model of generating the received signals. Then proposed algorithm for fast localization of underwater acoustic in presence of correlated noise is presented in Section 4. Section 5 proposes the new versions of “MUSIC” without eigendecomposition. Some numerical results and Experimental tests are shown in Sections 6 and 7, respectively. Finally, Section 8 summarizes the main conclusions of this chapter.

Throughout the chapter, we use to denote: transpose operation “T,” complex conjugate transpose “+, ” complex conjugate “∗,” expectation operator E[.], cumulant Cum(.), Kronecker product ⊗, determinant ∣·∣ and Frobenius norm ∥·∥_F.

2. Problem formulation

Consider a transmitter that generates a plane wave with an angle θ_inc. The incident plane wave will propagate and be reflected by Pobjects. For example, when it is located in the bottom of a tank filled with sand and water. We name the objects, which reflect the signals, is the sources. An array composed of Nsensors receives Ksignals emitted by the sources (P < N). The received signals are grouped into a vector r(f), which is the Fourier transform of the array output vector at frequency f, is written as [31, 32]:

$\mathbf{r}\left(f\right)=\mathbf{A}\left(f\right)\mathbf{s}\left(f\right)+\mathbf{b}\left(f\right)$E1

where A(f) = [a(f, θ₁), a(f, θ₂), …, a(f, θ_K)], matrix of dimensions (N × P) is the transfer matrix of the source-sensor array systems with respect to some chosen reference point, s(f) = [s₁(f), s₂(f), …, s_K(f)]^Tis the vector of signals, b(f) = [b₁(f), b₂(f), …, b_K(f)]^Tis the vector of Gaussian white noise. We define the matrix interspersals by:

$\mathbf{\Gamma }\left(f\right)=E\left[\mathbf{r}\left(f\right){\mathbf{r}}^{+}\left(f\right)\right]$E2

This matrix is estimated by $\widehat{\mathbf{\Gamma }}\left(f\right)=\frac{1}{L_r}{\sum }_{l=1}^{L_r}{\mathbf{r}}_l\left(f\right){\mathbf{r}}_l^{+}\left(f\right)$where L_rrepresents the number of realizations. Thus the spectral matrix $\mathbf{\Gamma }\left(f\right)$is formed:

$\mathbf{\Gamma }\left(f\right)=\mathbf{A}\left(f\right){\mathbf{\Gamma }}_s\left(f\right){\mathbf{A}}^{+}\left(f\right)+{\mathbf{\Gamma }}_b\left(f\right),$E3

where Γ_b(f) = E[b(f)b⁺(f)] is the spectral matrix of noise vector, the spectral matrix of signal vector is given as:

${\mathbf{\Gamma }}_s\left(f\right)=E\left[\mathbf{s}\left(f\right){\mathbf{s}}^{+}\left(f\right)\right]=\mathbf{V}\left(f\right)\mathbf{\Lambda }\left(f\right){\mathbf{V}}^{+}\left(f\right)+{\mathbf{\Gamma }}_b\left(f\right)$E4

where Λ(f) = diag {λ₁(f), …, λ_P(f)} and V(f) = [v₁(f), …, v_P(f)]. Assuming that the columns of A(f) are linearly independent, in other words, A(f) is full rank, it follows that for nonsingular Γ_s(f), the rank of A(f)Γ_s(f)A⁺(f) is P. This rank property implies that:

• the (N − P) multiplicity of its smallest eigenvalues: λ_P + 1(f) = … = λ_N(f) ≅ σ²(f).

• the eigenvectors corresponding to the minimal eigenvalues are orthogonal to the columns of A(f), namely, V_b(f) is equal to by the definition of {V_P + 1(f)…V_N(f)} orthogonal to {a(f, θ₁)…a(f, θ_P)}.

The eigenstructure-based techniques are based on the exploitation of these properties. When the objects are far away from the array, the wavefront is assumed to be plane. Then DOA of the sources are obtained, at the frequency f, by the peak positions in a so-called spectrum (MUSIC) defined as:

$\mathbf{Z}\left(f,\theta \right)=\frac{1}{{\mathbf{a}}^{+}\left(f,\theta \right){\mathbf{V}}_b\left(f\right){\mathbf{V}}_b^{+}\left(f\right)\mathbf{a}\left(f,\theta \right)}$E5

where $\mathbf{a}\left(f,\theta \right)=\left[1,e^{-2\mathit{j\pi f}\frac{dsin\left(\theta \right)}{c}},\dots ,e^{-2\mathit{j\pi f}\left(N-1\right)\frac{dsin\left(\theta \right)}{c}}\right]$is the steering vector of plane wave model, V_bis the eigenvectors of the noise subspace, cis the sound speed, dis the interspacing of the sensors and jis the complex operator. In the presence of Pobjects, the 1-D MUSIC(f, θ) algorithm cannot solve all the Pangles because the signals are correlated. In the following sections, we use simultaneously all the information contained in the signals to estimate the coherent signal subspace which extends the conventional 1-D MUSIC(f, θ) algorithm to 2-D MUSIC(f, θ, ρ) for joint range ρand DOA θestimation when the objects are buried in the sand with small depth.

3. Scattering acoustic model: to generate the received signals

In this section, we will present how to fill the vector of the scattering model. We consider a sedimentary covered with water and the interface is treated as a plane. An object of cylindrical or spherical shell is buried in the sediment. An incident plane wave propagating in the water reaches the interface with an angle of incidence θ_incas show in Figure 1. The incident plane wave generates a wave reflecting plane in the water and refracted plane wave propagating in the sediment. So the array located in the water receives three components [18]:

• the incident plane wave,

• the reflecting plane wave,

• the transmitter plane wave diffused by the object.

Figure 1.

Geometry configuration of buried object.

The array-interface height hand the nature of the sediment are known or can be determined. Also, the speed of wave propagation in the sediment c₂ is assumed to be known. Because the object is buried, the pressure in the water and sediment will not be expressed directly in terms of θ₁ and ρ₁, but in terms of five unknown parameters θ₁₁, ρ₁₁, θ₁₂, ρ₁₂ and y_c(the depth of buried object). So we will express θ₁₁, ρ₁₁, θ₁₂, ρ₁₂ and y_cbased on θ₁ and ρ₁ (see Figure 2). We use the law of Snell-Descartes and generalize the Pythagorean theorem to obtain the expressions:

$y_c={\rho }_{1}cos\left({\theta }_1\right)-h$E6

${\theta }_{12}=arcsin\left(\frac{c_2}{c_1}sin\left({\theta }_{\mathit{inc}}\right)\right)$E7

${\rho }_{12}=\frac{{\rho }_{1}cos\left({\theta }_1\right)-h}{cos\left[arcsin\left(\frac{c_2}{c_1}sin\left({\theta }_{\mathit{inc}}\right)\right)\right]}$E8

${\theta }_{11}=arctan\left(\frac{{\rho }_{1}cos\left({\theta }_1\right)-{\rho }_{12}cos\left({\theta }_{12}\right)}{{\rho }_{1}sin\left({\theta }_1\right)-{\rho }_{12}sin\left({\theta }_{12}\right)}\right)$E9

${\rho }_{11}=\frac{h}{cos\left[arctan\left(\frac{{\rho }_{1}cos\left({\theta }_1\right)-{\rho }_{12}cos\left({\theta }_{12}\right)}{{\rho }_{1}sin\left({\theta }_1\right)-{\rho }_{12}sin\left({\theta }_{12}\right)}\right)\right]}$E10

Figure 2.

Configuration of the buried object-1st sensor.

3.1. Cylindrical shell

Assume a cylindrical shell long enough which is buried in the sediment with axis parallel to the interface plane. Thus the acoustic pressure wave received by the first sensor of the array P_cyl(f, θ_k1, ρ_k1) contains three acoustic pressure components:

$P_{\mathit{cyl}}\left(f,{\theta }_{k1},{\rho }_{k1}\right)=P_{\mathit{\text{inwater}}-\mathit{cyl}}+P_{\mathit{\text{refwater}}-\mathit{cyl}}+P_{\mathit{\text{diff}}-\mathit{cyl}}$E11

where P_{inwater − cyl} = e^{jk₁(−(ρ_k1 sin(θ_k1)) sin(θ_inc) + hcos(θ_inc))} is the pressure incident in the water, P_{refwater − cyl}= R(θ_inc)e^{jk₁((ρ_k1 sin(θ_k1)) sin(θ_inc) − hcos(θ_inc))} is the pressure reflected by the sediment-water interface, where R(θ_inc) is the reflection coefficient of the interface, $P_{\mathit{\text{diff}}-\mathit{cyl}}={\sum }_{m=-\infty }^{+\infty }\xi {\mathbf{T}}_c{\left(\mathbf{I}-{\mathbf{D}}_c\right)}^{-1}{\Psi }_{\mathit{cyl}}^t$is the diffused acoustic pressure wave transmitted in the water, where Iis the identity matrix, D_cis a linear operator, T_cis the transition diagonal matrix, ${\Psi }_{\mathit{cyl}}^t$is the vector of transmitted wave and ξ = [ξ₁, ξ₂, …, ξ_m] is defined by ξ_m = T_{water − sed}(θ_inc)e^{jk₂y_c cos(θ_k11)}j^me^{−jm(π − θ_k11)}, where T_{water − sed}(θ_inc) is the transmission coefficient.

3.2. Spherical shell

In this section P_sph(f, θ_k1, ρ_k1) is the acoustic pressure wave received by the first sensor and is expressed as follows [12, 18]:

$P_{\mathit{sph}}\left(f,{\theta }_{k1},{\rho }_{k1}\right)=P_{\mathit{\text{inwater}}-\mathit{sph}}+P_{\mathit{\text{refwater}}-\mathit{sph}}+P_{\mathit{\text{diff}}-\mathit{sph}}$E12

where $P_{\mathit{\text{inwater}}-\mathit{sph}}={\sum }_{m=-\infty }^{+\infty }{j}^m\left[P_m\left(cos\left({\theta }_{\mathit{inc}}\right)\right)e^{2jcos\left({\theta }_{\mathit{inc}}\right)y_c}\right]$is the incident wave generates in the water, $P_{\mathit{\text{refwater}}-\mathit{sph}}={\sum }_{m=-\infty }^{+\infty }R\left({\theta }_{\mathit{inc}}\right)P_m\left(-cos\left({\theta }_{\mathit{inc}}\right)\right)e^{2jcos\left({\theta }_{\mathit{inc}}\right)y_c}$is the reflected wave, where R(θ_inc) is the reflection coefficient of the interface. We define the acoustic wave of the diffused by spherical shell by $P_{\mathit{\text{sphsediment}}}={\left({\mathbf{T}}_s^{-1}-{\mathbf{C}}^m\right)}^{-1}{P}_{\mathit{\text{inwater}}-\mathit{sph}}$, where T_sis the transition matrix and Cis the matrix containing the conversion coefficients. The diffused acoustic pressure wave transmitted in the water is given by: $P_{\mathit{\text{diff}}-\mathit{sph}}={\sum }_{m=0}^{\infty }{\epsilon }_{m}Ycos\left(m\left({\theta }_{k1}-{\theta }_{12}\right)\right)$, where Y = [P_sphsedimentP_m(cos(θ_k1))h_m(k₂ρ_k1) + X_m], ε_m = 2 for m > 0, ε₀ = 1.

The vector a(f, θ_K, ρ_K) is filled with cylindrical or spherical scattering model considering the sources shape. For example, when the sources are cylindrical shells, the vector is given by:

$\mathbf{a}\left(f,{\theta }_K,{\rho }_K\right)=\left[P_{\mathit{cyl}}\left(f,{\theta }_{k1},\mathit{\rho k}1\right),\dots ,P_{\mathit{cyl}}\left(f,{\theta }_{\mathit{kN}},{\rho }_{\mathit{kN}}\right)\right]$E13

Eq. (11) or (12) give the first component of the vector. The other P_cyl(f, θ_ki, ρki) for i = 1, 2, …, Nassociated with the ith sensors can be formed by a geometric recursive relationship. The relationship allows to express (θ_ki, ρ_ki) according to (θ_ki − 1, ρ_ki − 1) (see Figure 1). This recursive calculate is done as follows:

${\rho }_{\mathit{ki}}=\sqrt{{\rho }_{\mathit{ki}-1}^2+d^2-2{\rho }_{\mathit{ki}-1}dcos\left(\frac{\pi }{2}+{\theta }_{\mathit{ki}-1}\right)}$E14

${\theta }_{\mathit{ki}}=-\frac{\pi }{2}+{cos}^{-1}\left(\frac{d^2+{\rho }_{\mathit{ki}}^2-{\rho }_{\mathit{ki}-1}^2}{2{\rho }_{\mathit{ki}-1}d}\right),i=2,\dots ,N$E15

These equations are employed in Eq. (5) to estimate simultaneously range and bearing of the objects.

In the following section, we summarize the proposed algorithm for fast localization of underwater acoustic sources using a wide-band transmitter to receive the signals at different frequencies, then the coherent signal subspace can be applied to decorrelate the source signals.

4. Proposed algorithm for fast localization of underwater acoustic sources

We use spatial smoothing to deal with narrow band correlated signal, we divide the array into L_soverlapping subarrays. The spatially smoothed covariance matrix is the average of the subarray covariances [33]. The step-by-step proposed algorithm for fast localization of underwater acoustic sources is given as following:

Algorithm 1Proposed Algorithm for Bearing and Range Estimation of Buried Objects

1. use the beamformer method to find an initial estimate of $\widehat{{\theta }_{0k}}$, where k = 1, …, P₀, with P₀ ≤ P.

2. compute the initial values of ${\rho }_k=\frac{X}{cos\left({\theta }_k\right)}$for k = 1, …, P₀, where X = h + y_crepresents the distance between the receiver and the bottom of the tank (seabed),

3. fill the transfer matrix, $\widehat{\mathbf{A}}\left(f\right)=\left[\mathbf{a}\left(f,{\theta }_1,{\rho }_1\right),\mathbf{a}\left(f,{\theta }_2,{\rho }_2\right),\dots ,\mathbf{a}\left(f,{\theta }_K,{\rho }_K\right)\right]$, where each source steering vector is filled using Eq. (11),

4. estimate the spectral matrix $\mathbf{\Gamma }\left(f\right)=E\left[\mathbf{r}\left(f\right){\mathbf{r}}^{+}\left(f\right)\right]=\frac{1}{L_r}{\sum }_{l=1}^{L_r}{\mathbf{r}}_l\left(f\right){\mathbf{r}}_l^{+}\left(f\right)$, where L_lis the realization number,

5. estimate noise covariance matrix Γ_b(f) = E[b(f)b⁺(f)]. When the noise is white noise, that is, estimate noise variance ${\sigma }^2\left(f\right)=\frac{1}{N-P_0}{\sum }_{i=P_0+1}^N{\lambda }_i$, where λ_iis the i^theigenvalue of Γ(f). Then we calculate Γ_b(f) = σ²I, where Iis the identify matrix,

6. calculate the spectral matrix of the signals reflected on the objects by ${\mathbf{\Gamma }}_s\left(f\right)={\left({\widehat{\mathbf{A}}}^{+}\left(f\right)\widehat{\mathbf{A}}\left(f\right)\right)}^{-1}{\widehat{\mathbf{A}}}^{+}\left(f\right)\times \left(\mathbf{\Gamma }\left(f\right)-{\mathbf{\Gamma }}_n\left(f\right)\right)\widehat{\mathbf{A}}\left(f\right){\left({\widehat{\mathbf{A}}}^{+}\left(f\right)\widehat{\mathbf{A}}\left(f\right)\right)}^{-1}$,

7. compute the average of the spectral matrices ${\overline{\mathbf{\Gamma }}}_s\left(f\right)=\frac{1}{L_s}\sum _{S=1}^{L_S}{\mathbf{\Gamma }}_s\left(f\right)$, where L_srepresents the number of subarrays, then calculate ${\overline{\mathbf{V}}}_s\left(f\right)$by SVD,

8. calculate the spatial spectrum of the MUSICmethod for bearing and range object estimation:

$\mathit{\text{MUSIC}}=\frac{1}{\mid \mid \mathbf{a}{\left(f_0,{\theta }_k,{\rho }_k\right)}^{+}{\overline{\mathbf{V}}}_b\left(f_0\right){\overline{\mathbf{V}}}_b^{+}\left(f_0\right)\mathbf{a}\left(f_0,{\theta }_k,{\rho }_k\right)\mid \mid },$E16

where ${\overline{\mathbf{V}}}_b$is the eigenvector matrix of noise subspace associated with the (N − P) smallest eigenvalues.

Inverse Power method can be used to estimate the noise power λ_N(f) = σ²(f). This variance can be used at step 5 in Alg. 1 for estimating the noise variance in case of white noise. The principle of this method is recalled in the following, using the maximum norm.

1. Let q_oa complex vector of Nelements, ∥q_o∥_∞ = 1;

2. For l = 1, 2, 3, …;

3. Calculate x_l : Γx_l = q_l − 1

   μ^l =  ∥ x_l∥_∞

   ${\mathbf{q}}_l=\frac{{\mathbf{x}}_l}{{\mu }^l}$;

It is shown that: ${lim}_{l\to \infty }{\mu }^l=\frac{1}{{\lambda }_N}$, where λ_Nis the smallest eigenvalue of Γ.

In the high resolution noise subspace based methods, the DOA’s are given by the local maximum points of a cost function, for example Eq. (16) of MUSIC. V_b(f) is the orthogonal projector onto the noise subspace given by the eigenvectors associated with the smallest eigenvalues of the covariance matrix of the received data. This requires an enormous computational load which limits its use for tracking by SVD [34]. In many applications, only a few eigenvectors are required. Since the number of sensor Nis often larger than the number of sources P. It means that the vector dimension of noise subspace is larger than signal subspace. It is more efficient to work with the lower dimensional signal subspace than with the noise subspace. That is to say, it is not necessary to obtain V_b(f) exactly. We can calculate signal subspace V_s(f) = [v₁(f), v₂(f), …, v_P(f)] whose columns are the Porthonormal basis vectors. The projector onto the noise subspace spanned by the (N − P) eigenvectors associated with the (N − P) smallest eigenvalues is ${\mathbf{V}}_b\left(f\right){\mathbf{V}}_b^{+}\left(f\right)$, given by:

${\mathbf{V}}_b\left(f\right){\mathbf{V}}_b^{+}\left(f\right)=\mathbf{I}-{\mathbf{V}}_s\left(f\right){\mathbf{V}}_s^{+}\left(f\right)$E17

On the other hand, the additive noise is assumed to be white. But in practice, the noise is not always spatially white noise. In generally, the noise is correlated or unknown.

So in the next two sections, we will introduce the algorithms to replace SVD in MUSIC for reduce computation times and propose a new algorithm for estimating the spectral matrix of an unknown limited length spatially correlated noise.

5. MUSIC without eigendecomposition

In this section, we propose the noneigenvector versions of “MUSIC” to replace SVD to accelerate computation times.

5.1. Fixed point algorithm

One way to compute the Porthonormal basis vectors is to use Gram-Schmidt method. The eigenvector with dominant eigenvalue will be measured first. Similarly, all the remaining P − 1 basis vectors will be measured one by one in a reducing order of dominance. The previously measured (p − 1)^th basis vectors will be utilized to find the p^th basis vector. The algorithm for p^th basis vector will converge when the new value ${\mathbf{v}}_p^{+}$and old value v_pare such that ${\mathbf{v}}_p^{+}{\mathbf{v}}_p=1$. It is usually economical to use a finite tolerance error to satisfy the convergence criterion $\left|\left|{\mathbf{v}}_p^{+}{\mathbf{v}}_p-1\right|\right|<\eta$where ηis a prior fixed threshold. The proposed algorithm is given as follows:

Algorithm 2Fixed Point Algorithm

1. Choose P, the number of principal axes or eigenvectors required to estimate. Consider covariance matrix $\mathbf{\Gamma }$and set p ← 1.

2. Initialize eigenvector v_pof size d × 1, e.g. randomly;

3. while $\left|\left|{\mathbf{v}}_p^H{\mathbf{v}}_p-1\right|\right|<\eta$

   1. Update v_pas ${\mathbf{v}}_p\leftarrow \mathbf{\Gamma }{\mathbf{v}}_p$;

   2. Do the Gram-Schmidt orthogonalization process ${\mathbf{v}}_p\leftarrow {\mathbf{v}}_p-{\sum }_{j=1}^{j=p-1}\left({\mathbf{v}}_p^T{\mathbf{v}}_j\right){\mathbf{v}}_j$;

   3. Normalize v_pby dividing it by its norm: ${\mathbf{v}}_p\leftarrow \frac{{\mathbf{v}}_p}{\left|\left|{\mathbf{v}}_p\right|\right|}$.

4. Increment counter p ← p + 1 and go to step 2 until pequals K.

5.2. Projection approximation subspace tracking (PAST) algorithm

Suppose that we have an estimation of the signal subspace W(t) where each column is an eigenvector. The Linear Principal Analysis Criterion gives the definition of the scalar cost function J(W(t)).

$J\left(\mathbf{W}\left(t\right)\right)=E\left\{{||\mathbf{r}\left(t\right)-\mathbf{W}\left(t\right){\mathbf{W}}^{+}\left(t\right)\mathbf{r}\left(t\right)||}^2\right\}$E18

where W(t)W⁺(t)r(t) is the projection of r(t) into the subspace W(t). The error surface of the function has several local minimal and one global minimum. When W(t) is equal to a basis for the signal subspace, J(W(t)) has a global minimum which can estimate the signal subspace by Eq. (18). Note that W(t) is not equal to the signal subspace itself, but merely provides a possible basis. If W(t) is a signal column vector, it does indeed become equal to the Principal Component (dominant eigenvector) under minimization.

The cost function J(W(t)) can be minimized by the application of a gradient-descent technique or recursive least squares variant. We can replace the expectation operator in Eq. (18) by an exponentially weighted sum over nsamples. The estimation is given as follows:

$\widehat{J}\left(J\left(\mathbf{W}\left(t\right)\right)\right)=\sum _{t=1}^n{\beta }^{n-t}{||\mathbf{r}\left(t\right)-\mathbf{W}\left(t\right)\mathbf{W}{\left(t\right)}^{+}\mathbf{r}\left(t\right)||}^2$E19

where βis the forgetting factor (0 < β < 1). The forgetting factor allows the subspace estimation to track geostationary signal over time.

We get another cost function by approximating W(t)⁺r(t) with W(t − 1)⁺r(t), using the previous value of W(t)⁺ in the iteration, giving:

${\widehat{J}}^{\text{'}}\left(\mathbf{W}\left(t\right)\right)=\sum _{t=1}^n{\beta }^{n-t}{||\mathbf{r}\left(t\right)-\mathbf{W}\left(t\right)\mathbf{W}{\left(t-1\right)}^{+}\mathbf{r}\left(t\right)||}^2$E20

This function resembles the cost function used to define a recursive least squares (RLS) filter:

$C\left(\mathbf{W}\left(t\right)\right)=\sum _{t=1}^n{\beta }^{n-t}{\left|\mathbf{e}\left(t\right)\right|}^2$E21

where e(t) is the error signal. The error signal is the difference between the “desired” signal r(t) and its projection into the subspace W(t)W(t − 1)⁺r(t). Consequently, the PAST algorithm may be summarized by the following equations:

Algorithm 3PAST Algorithm

1. Initialization:

   W(0) and P(0)

2. for t = 1,2…

   y(t) = W⁺(t − 1)r(t)

   h(t) = P(t − 1)^∗y(t)

   $\mathbf{g}\left(t\right)=\frac{\mathbf{h}\left(t\right)}{\left(\beta +\mathbf{y}\left(t\right)\right)\mathbf{h}\left(t\right)}$

   $\mathbf{P}\left(t\right)=\frac{1}{\beta }$Tri{P(t − 1) − g(t)r⁺(t)}

   e(t) = r(t) − W(t − 1)y(t)

   W(t) = W(t − 1) + e(t)g⁺(t)

   end

The operator Triindicates that only the upper (or lower) triangular part of the matrix is calculated and its Hermitian transposed version is copied to the another lower (or upper) triangular part.

5.3. Projection approximation subspace tracking with deflation (PASTD) algorithm

The PAST algorithm provides a method to estimate only a basis for the dominant subspace. The exact eigenvectors (singular vector) are not calculated unless W(t) is a column vector in which case only the dominant eigenvector (principal component) is estimated. We present a second subspace tracking algorithm - PAST with deflation (PASTD) which is derived from the PAST approach. The PASTD algorithm is based on the deflation technique which is the sequential estimation of the principal components. According to the Karhunen-Loève expansion

$\mathbf{r}\left(t\right)=\sum _{i=1}^N{\mathbf{W}}_i^{+}\left(t\right)\mathbf{r}\left(t\right){\mathbf{W}}_i\left(t\right)$E22

tells that any r(t) may be expressed as a linear combination of the eigenvector of the correlation matrix.

The first step of PASTD is to update the most dominant eigenvector by applying PAST with i = 1, then the contribution of the dominant eigenvector in Eq. (22) is removed by subtraction. So the second dominant eigenvector becomes the most dominant and can be extracted in the same way. Then we repeat the procedure until all desired eigencomponents are estimated. This iterative process is called deflation. So the algorithm may be summarized as follows:

Algorithm 4PASTD Algorithm

1. Initialization:

   W_i(0) and d_i(0)

2. for n = 1,2…

   r₁(t) = r(t)

   for i = 1,2,…,P ${\mathbf{y}}_i\left(t\right)={\mathbf{W}}_i^{+}\left(t-1\right){\mathbf{r}}_i\left(t\right)$

   d_i(t) = βd_i(t − 1) + |y_i(t)|²

   e_i(t) = r_i(t) − W_i(t − 1)y_i(t)

   ${\mathbf{W}}_i\left(t\right)={\mathbf{W}}_i\left(t-1\right)+{\mathbf{e}}_i\left(t\right){\mathbf{y}}_i^{\ast }\left(t\right)/d_i\left(t\right)$

   r_i + 1(t) = r_i(t) − W_i(t)y_i(t)

   end

   end

where estimates are made of Peigenvectors with the largest eigenvalues. In practice since P <  < N, this indicates an important optimization compared to the eigendecomposition or SVD. The eigenvector projection estimates W_iare initialized to the columns of some nonzero orthogonal matrix. d_i(t) is initialized to arbitrary nonzero constants. When PASTD has converged, the W_i(t) will contain estimates of the eigenvector of the correlation matrix of the data in r(t). The corresponding eigenvalues may be calculated by multiplying the d_i(t) by $\frac{1-\beta }{\beta }$.

5.4. Orthogonal projection approximation subspace tracking (OPAST) algorithm

The OPAST algorithm is the modification of PAST. The weight matrix W(t) is forced to be orthonormal to each iteration. So we can get:

$\mathbf{W}\left(t\right)=\mathbf{W}\left(t\right){\left({\mathbf{W}}^{+}\left(t\right)\mathbf{W}\left(t\right)\right)}^{-1/2}$E23

where (W⁺(t)W(t))^−1/2 denotes an inverse square root of (W⁺(t)W(t)). (W⁺(t)W(t))^−1/2 can be calculated by using the updating equation of W(t). Note that W(t − 1) is now an orthonormal matrix, we have

${\mathbf{W}}^{+}\left(t\right)\mathbf{W}\left(t\right)=\mathbf{I}+\parallel \mathbf{p}\left(t\right){\parallel }^2\mathbf{q}\left(t\right){\mathbf{q}}^{+}\left(t\right)$E24

where Iis the identity matrix, W⁺(t − 1)p(t) = 0 and $\mathbf{r}\stackrel{\mathit{def}}{=}\parallel \mathbf{p}\left(t\right)\parallel \mathbf{q}\left(t\right)$. Thus

$\begin{array}{ll}{\left({\mathbf{W}}^{+}\left(t\right)\mathbf{W}\left(t\right)\right)}^{-1/2}\hfill & =\mathbf{I}+\frac{1}{\parallel \mathbf{r}{\parallel }^2}\left(\frac{1}{\sqrt{1+\parallel \mathbf{r}{\parallel }^2}}-1\right){\mathbf{xx}}^H\\ & =\mathbf{I}+\tau \left(t\right)\mathbf{q}\left(t\right){\mathbf{q}}^{+}\left(t\right)\hfill \end{array}$E25

where

$\tau {\left(t\right)}^{def}=\frac{1}{\parallel q\left(t\right){\parallel }^2}\left(\frac{1}{\sqrt{1+\parallel p\left(t\right){\parallel }^2\parallel q\left(t\right){\parallel }^2}}-1\right)$E26

Using Eqs. (23) and (26), and the updating equation of W(t), we obtain

$\begin{array}{ll}\mathbf{W}\left(t\right)\hfill & =\left(\mathbf{W}\left(t-1\right)+\mathbf{p}\left(t\right){\mathbf{q}}^{+}\left(t\right)\right)\left(\mathbf{I}+\tau \left(t\right)\mathbf{q}\left(t\right){\mathbf{q}}^{+}\left(t\right)\right)\\ & =\mathbf{W}\left(t-1\right)+{\mathbf{p}}^{\mathbf{\prime }}\left(t\right){\mathbf{q}}^{+}\left(t\right)\hfill \end{array}$E27

where p^'(t) = τ(t)W(t − 1)q(t) + (1 + τ(t) ∥ q(t)∥²)p(t). Thus, the OPAST algorithm can be written as the PAST (see Algorithm 3):

Algorithm 5OPAST Algorithm

W(t) = W(t − 1) + p^'(t)q⁺(t)

$\tau \left(t\right)=\frac{1}{\parallel \mathbf{q}\left(t\right){\parallel }^2}\left(\frac{1}{\sqrt{1+\parallel \mathbf{p}\left(t\right){\parallel }^2\parallel \mathbf{q}\left(t\right){\parallel }^2}}-1\right)$

W(t) = τ(t)W(t − 1)q(t) + (1 + τ(t) ∥ q(t)∥²)p(t)

6. Simulations results

6.1. Complexity

The traditional MUSIC method estimate the noise subspace eigenvectors by SVD. From the computational point of view, the well-known SVD method is the cyclic Jacobi’s method which requires around N³ computations. The computational complexity of fixed-point algorithm, PAST, PASTD and OPAST is (NP² + N²P), 3NP + O(P²), 4NP + O(P) and 4NP + O(P²) respectively. If the number of sensors Nis larger compared to the number of objects P, the computational complexity can be estimated to be around N²Pfor fixed-point algorithm, 3NPfor PAST, 3NPfor PASTD and 4NPfor OPAST. Several experiments are carried out with various numbers of sensors, to study the computational load of the proposed algorithm with SNR = 0 dB and P = 8 sources. DOA values are: 5^∘, 10^∘, 20^∘, 25^∘, 35^∘, 40^∘, 50^∘ and 55^∘.

The number of realizations is 1000, and the number of observations is 1000. Choosing a number of snapshots equal to 100, such as in [14, 21, 22, 29, 30], does not change the results. The mean computational load is then up to 2.5 times less with fixed point algorithm than with SVD (see Figure 3 and Table 1, N = 10 up to 30). Both versions of MUSIC provide the same results (see Figure 4, take fixed-point algorithm for example).

Figure 3.

Computational times, SVD (red), Fixed–point (green), PAST (Black), OPAST(blue) and PASTD (pink).

	10	15	20	25	30
Time SVD (second)	0.95	1.3	2.4	4.4	7.1
Time fixed point (second)	0.5	0.6	1.1	1.8	2.8
Time PAST (second)	0.4	0.5	0.8	1.4	2.2
Time PASTD (second)	0.3	0.4	0.7	1.3	2.0
Time OPAST (second)	0.5	0.7	1.3	2.1	2.9
Ratio SVD/fixed point	1.9	2.2	2.2	2.4	2.5
Ratio SVD/PAST	2.4	2.6	3.0	3.1	3.2
Ratio SVD/PASTD	3.2	3.3	3.4	3.4	3.6
Ratio SVD/OPAST	1.9	1.9	1.8	2.1	2.4

Table 1.

Computational time needed to run MUSIC for various numbers of sensors.

Figure 4.

(a) Pseudospectrum of MUSIC obtained using fixed point, (b) pseudospectrum of MUSIC obtained using SVD.

7. Experimental setup

The studied signals are recorded during an underwater acoustic experiment in order to estimate the developed method performance. The experiment is carried out in an acoustic tank under the conditions similar to those in a marine environment. The bottom of the tank is filled with sand. The experimental device is presented in Figure 5. The tank is topped by two mobile carriages. The first carriage supports a transducer issuer and the second supports a transducer receiver pilot by the computer.

Figure 5.

Experimental setup: (a) Data acquisition system, (b) Experimental tank.

Four couples of spherical and cylindrical shells (see Figure 6) are buried between 0 and 0.05 m under the sand. The considered objects have the following characteristics, where δrepresents the distance between the two objects of the same couple and ∅_athe outer radius (the inner radius ∅_b = ∅_a − 0.001 m):

1. the 1st couple (O₁, O₂): spherical shells, ∅_a = 0.3 m, δ = 0.33 m, full of air,

2. the 2nd couple (O₃, O₄): cylindrical shells, ∅_a = 0.01 m, δ = 0.13 m, full of air,

3. the 3rd couple (O₅, O₆): cylindrical shells, ∅_a = 0.018 m, δ = 0.16 m, full of water,

4. the 4th couple (O₇, O₈): cylindrical shells, ∅_a = 0.02 m, δ = 0.06 m, full of air,

Figure 6.

Experimental objects.

The considered objects are made of dural aluminum with density D₂ = 1800 kg/m³, the speed of the wave in the water c₁ is 1500 m/s and in the sediment c₂ is 1700 m/s, the longitudinal and transverse-elastic wave velocities inside the shell medium are c_l = 6300 m/s and c_t = 3200 m/s, respectively. The speed of the wave in the water c₁ is 1500 m/s and in the sediment c₂ is 1700 m/s, The external fluid is water with density D₁ = 1000 kg/m³ and the internal fluid is water or air with density D_3air = 1.2·10⁻⁶ kg/m³ or D_3water = 1000 kg/m³.

In addition to estimate the performance of the propose method, the signal source a spatially correlated noise is emitted with K = 10. The objective is to estimate the directions of arrival of the signals during the experiment. The signals are received on one uniform linear array. The observed signals come from various reflections on the objects being in the tank. Generally the aims of acousticians is the detection, localization and identification of these objects. In this experiment we have recorded the reflected signals by a single receiver. This receiver is moved along a straight line between position X_min = 50 mm and position X_max = 150 mm with a step of α= 1 mm in order to create a uniform linear array. The experimental setup is shown in Figure 7. We have measured eight times E_i(O_ii, O_ii + 1) with i = 1, …, 8 and ii = 1, 3, 5, 7. At first, the receiver horizontal axis XX^' is fixed at 0.2 m, we performed the experiments E₁(O₁, O₂), …, E₄(O₇, O₈) associated to the 1st, 2nd, 3rd and the 4th couple, respectively. Then we performed the other four experiments E₅(O₁, O₂), …, E₈(O₇, O₈) with XX^' fixed at 0.4 m. RR^' is the vertical axis which goes through the center of the first object of each couple. For each experiment, the transmitted signal had the following properties: pulse duration is 15 μs, the frequency band is, the frequency of the band is [f_min = 150, f_max = 250] kHz and the center frequency f₀ is f₀ = 200 kHz. The sampling rate is 2 MHz. The duration of the received signal was 700 μs. The variance of Gaussian white noise σ² is 100 and the angle of incidence θ_incis 60^∘.

Figure 7.

Experimental setup.

At each sensor, time-domain data corresponding only to target echoes are collected with signal to noise ratio equal to 20 dB. The typical sensor output signals recorded during one experiment is shown in Figure 8.

Figure 8.

Example of observed signals during experiment Exp. 1.

The proposed algorithms were applied on each experimental data set. Forming the directional vector by the model of acoustic diffusion appropriate to locate the objects and the spectral matrix of the simulated data. We use the focusing operator on the signals by dividing the frequency band [150, 250] kHz in 11 frequencies and Alg. 2, Alg. 3, Alg. 4 and 5 to calculate the noise subspace, respectively. Finally we apply Eq. (16) to estimate DOA of objects and object-1st sensor distance.

As is shown in Figure 9, X axis is the object-1st sensor distance ρ, Y axis is the DOA of object-1st sensor θ. The white points and the peak positions, which present the maximum values, correspond to the coordinate of 2 objects (29^∘, 0.31 m) and (33^∘, 0.34 m). The bearings and the ranges of buried objects are (28.1^∘, 0.298 m) and (33.9^∘, 0.361 m) if we use conventional SVD (see Figure 9(a) and (b)). The results of the proposed algorithms are (29.5^∘, 0.301 m) and (33.3^∘, 0.351 m) for fixed-point algorithm (see Figure 9(c) and (d)), (29.2^∘, 0.312 m) and (32.8^∘, 0.343 m) for PAST algorithm (see Figure 9(e) and (f)), (29.8^∘, 0.313 m) and (33.3^∘, 0.355 m) for PASTD algorithm (see Figure 9(g) and (h)) and (29.3^∘, 0.306 m) and (32.5^∘, 0.334 m) for OPAST algorithm (see Figure 9(i) and (j)).

Figure 9.

Example of object localization with different methods: (a)-(b) SVD, (c)-(d) Fixed-point, (e)-(f) PAST, (g)-(h) PASTD and (i)-(j) OPAST.

We have done statistical study in order to a posteriori verify the quality of estimation of the proposed method. Standard Deviation (Std) is defined as follows:

$\mathit{std}=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(X_{\mathit{\text{iexp}}}-X_{\mathit{\text{iest}}}\right)}^2}$E28

where X_iexp(respectively X_iest) represents the ith expected (respectively the ith estimated) value of θor ρ. Standard deviation of the bearing and the range estimation at different signal-to-noise ratio (SNRs) (from −10 to 20 dB) are given in Figure 10.

Figure 10.

Standard deviation versus SNR of the bearing and the range estimation.

8. Conclusion

The main target of array processing is the estimation of the parameters: DOA of objects and the objects-sensors distance. In this chapter, we have proposed a new fast localization algorithm to estimate both the ranges and the bearings of buried sources underwater acoustic in presence of correlated noise. This algorithm takes into account both the reflection and refraction of water-sediment interface. We develop fixed point algorithm in MUSIC instead of SVD to keep the small computational time load. A new focusing operator is proposed to estimate the coherent signal subspace. Some simulations have been done to test our method. We compare the computation time of MUSIC with SVD and fixed point, it shows that fixed point is faster than SVD. The proposed method performance was investigated through scaled tank tests associated with some cylindrical and spherical shells buried in an homogenous fine sand. The obtained results are promising and the REcalculated between the expected and the estimated bearings and ranges errors is weep.