Open access peer-reviewed chapter

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

Written By

Salah Bourennane and Caroline Fossati

Submitted: 13 June 2017 Reviewed: 26 September 2017 Published: 22 November 2017

DOI: 10.5772/intechopen.71272

From the Edited Volume

Advances in Underwater Acoustics

Edited by Andrzej Zak

Chapter metrics overview

1,394 Chapter Downloads

View Full Metrics


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.


  • array processing
  • source localization
  • wide-band
  • fast algorithm

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 [24] 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 P objects. 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 N sensors receives K signals 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]:


where A(f) = [a(f, θ1), a(f, θ2), …, 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) = [s1(f), s2(f), …, sK(f)]T is the vector of signals, b(f) = [b1(f), b2(f), …, bK(f)]T is the vector of Gaussian white noise. We define the matrix interspersals by:


This matrix is estimated by Γ̂f=1Lrl=1Lrrlfrl+f where Lr represents the number of realizations. Thus the spectral matrix Γf is formed:


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


where Λ(f) = diag {λ1(f), …, λP(f)} and V(f) = [v1(f), …, vP(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) ≅ σ2(f).

  • the eigenvectors corresponding to the minimal eigenvalues are orthogonal to the columns of A(f), namely, Vb(f) is equal to by the definition of {VP + 1(f)…VN(f)} orthogonal to {a(f, θ1)…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:


where afθ=1e2jπfdsinθce2jπfN1dsinθc is the steering vector of plane wave model, Vb is the eigenvectors of the noise subspace, c is the sound speed, d is the interspacing of the sensors and j is the complex operator. In the presence of P objects, the 1-D MUSIC(f, θ) algorithm cannot solve all the P angles 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 θinc as 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 h and the nature of the sediment are known or can be determined. Also, the speed of wave propagation in the sediment c2 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 θ1 and ρ1, but in terms of five unknown parameters θ11, ρ11, θ12, ρ12 and yc (the depth of buried object). So we will express θ11, ρ11, θ12, ρ12 and yc based on θ1 and ρ1 (see Figure 2). We use the law of Snell-Descartes and generalize the Pythagorean theorem to obtain the expressions:


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 Pcyl(f, θk1, ρk1) contains three acoustic pressure components:


where Pinwater − cyl = ejk1(−(ρk1 sin(θk1)) sin(θinc) + hcos(θinc)) is the pressure incident in the water, Prefwater − cyl = R(θinc)ejk1((ρ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, Pdiffcyl=m=+ξTcIDc1Ψcylt is the diffused acoustic pressure wave transmitted in the water, where I is the identity matrix, Dc is a linear operator, Tc is the transition diagonal matrix, Ψcylt is the vector of transmitted wave and ξ = [ξ1, ξ2, …, ξm] is defined by ξm = Twater − sed(θinc)ejk2yc cos(θk11)jmejm(π − θk11), where Twater − sed(θinc) is the transmission coefficient.

3.2. Spherical shell

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


where Pinwatersph=m=+jmPmcosθince2jcosθincyc is the incident wave generates in the water, Prefwatersph=m=+RθincPmcosθince2jcosθincyc 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 Psphsediment=Ts1Cm1Pinwatersph, where Ts is the transition matrix and C is the matrix containing the conversion coefficients. The diffused acoustic pressure wave transmitted in the water is given by: Pdiffsph=m=0εmYcosmθk1θ12, where Y = [PsphsedimentPm(cos(θk1))hm(k2ρk1) + Xm], εm = 2 for m > 0, ε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:


Eq. (11) or (12) give the first component of the vector. The other Pcyl(f, θki, ρki) for i = 1, 2, …, N associated 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:


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 Ls overlapping 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 1 Proposed Algorithm for Bearing and Range Estimation of Buried Objects

  1. use the beamformer method to find an initial estimate of θ0k̂, where k = 1, …, P0, with P0P.

  2. compute the initial values of ρk=Xcosθk for k = 1, …, P0, where X = h + yc represents the distance between the receiver and the bottom of the tank (seabed),

  3. fill the transfer matrix, Âf=afθ1ρ1afθ2ρ2afθKρK, where each source steering vector is filled using Eq. (11),

  4. estimate the spectral matrix Γf=Erfr+f=1Lrl=1Lrrlfrl+f, where Ll is 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 σ2f=1NP0i=P0+1Nλi, where λi is the ith eigenvalue of Γ(f). Then we calculate Γb(f) = σ2I, where I is the identify matrix,

  6. calculate the spectral matrix of the signals reflected on the objects by Γsf=Â+fÂf1Â+f×ΓfΓnfÂfÂ+fÂf1,

  7. compute the average of the spectral matrices Γ¯sf=1LsS=1LSΓsf, where Ls represents the number of subarrays, then calculate V¯sf by SVD,

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


where 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) = σ2(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 qo a complex vector of N elements, ∥qo = 1;

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

  3. Calculate xl : Γxl = ql − 1

    μl =  ∥ xl


It is shown that: limlμl=1λN, where λN is 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. Vb(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 N is 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 Vb(f) exactly. We can calculate signal subspace Vs(f) = [v1(f), v2(f), …, vP(f)] whose columns are the P orthonormal basis vectors. The projector onto the noise subspace spanned by the (N − P) eigenvectors associated with the (N − P) smallest eigenvalues is VbfVb+f, given by:


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 P orthonormal 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 pth basis vector. The algorithm for pth basis vector will converge when the new value vp+ and old value vp are such that vp+vp=1. It is usually economical to use a finite tolerance error to satisfy the convergence criterion vp+vp1<η where η is a prior fixed threshold. The proposed algorithm is given as follows:

Algorithm 2 Fixed Point Algorithm

  1. Choose P, the number of principal axes or eigenvectors required to estimate. Consider covariance matrix Γ and set p ← 1.

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

  3. while vpHvp1<η

    1. Update vp as vpΓvp;

    2. Do the Gram-Schmidt orthogonalization process vpvpj=1j=p1vpTvjvj;

    3. Normalize vp by dividing it by its norm: vpvpvp.

  4. Increment counter p ← p + 1 and go to step 2 until p equals 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)).


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 n samples. The estimation is given as follows:


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:


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


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 3 PAST 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)


    Pt=1β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)


The operator Tri indicates 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


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 4 PASTD Algorithm

  1. Initialization:

    Wi(0) and di(0)

  2. for n = 1,2…

    r1(t) = r(t)

    for i = 1,2,…,P yit=Wi+t1rit

    di(t) = βdi(t − 1) + |yi(t)|2

    ei(t) = ri(t) − Wi(t − 1)yi(t)


    ri + 1(t) = ri(t) − Wi(t)yi(t)



where estimates are made of P eigenvectors with the largest eigenvalues. In practice since P <  < N, this indicates an important optimization compared to the eigendecomposition or SVD. The eigenvector projection estimates Wi are initialized to the columns of some nonzero orthogonal matrix. di(t) is initialized to arbitrary nonzero constants. When PASTD has converged, the Wi(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 di(t) by 1ββ.

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:


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


where I is the identity matrix, W+(t − 1)p(t) = 0 and r=defptqt. Thus




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


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

Algorithm 5 OPAST Algorithm

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


W(t) = τ(t)W(t − 1)q(t) + (1 + τ(t) ∥ q(t)∥2)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 N3 computations. The computational complexity of fixed-point algorithm, PAST, PASTD and OPAST is (NP2 + N2P), 3NP + O(P2), 4NP + O(P) and 4NP + O(P2) respectively. If the number of sensors N is larger compared to the number of objects P, the computational complexity can be estimated to be around N2P for fixed-point algorithm, 3NP for PAST, 3NP for PASTD and 4NP for 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).

Time SVD (second)0.951.
Time fixed point (second)
Time PAST (second)
Time PASTD (second)
Time OPAST (second)
Ratio SVD/fixed point1.
Ratio SVD/PAST2.

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 ∅a the outer radius (the inner radius ∅b = ∅a − 0.001 m):

  1. the 1st couple (O1, O2): spherical shells, ∅a = 0.3 m, δ = 0.33 m, full of air,

  2. the 2nd couple (O3, O4): cylindrical shells, ∅a = 0.01 m, δ = 0.13 m, full of air,

  3. the 3rd couple (O5, O6): cylindrical shells, ∅a = 0.018 m, δ = 0.16 m, full of water,

  4. the 4th couple (O7, O8): 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 D2 = 1800 kg/m3, the speed of the wave in the water c1 is 1500 m/s and in the sediment c2 is 1700 m/s, the longitudinal and transverse-elastic wave velocities inside the shell medium are cl = 6300 m/s and ct = 3200 m/s, respectively. The speed of the wave in the water c1 is 1500 m/s and in the sediment c2 is 1700 m/s, The external fluid is water with density D1 = 1000 kg/m3 and the internal fluid is water or air with density D3air = 1.2·10−6 kg/m3 or D3water = 1000 kg/m3.

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 Xmin = 50 mm and position Xmax = 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 Ei(Oii, Oii + 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 E1(O1, O2), …, E4(O7, O8) associated to the 1st, 2nd, 3rd and the 4th couple, respectively. Then we performed the other four experiments E5(O1, O2), …, E8(O7, O8) 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 [fmin = 150, fmax = 250] kHz and the center frequency f0 is f0 = 200 kHz. The sampling rate is 2 MHz. The duration of the received signal was 700 μs. The variance of Gaussian white noise σ2 is 100 and the angle of incidence θinc is 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:


where Xiexp (respectively Xiest) 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 RE calculated between the expected and the estimated bearings and ranges errors is weep.


  1. 1. Bourennane S, Frikel M. Localization of wideband sources with estimation of an antenna shape. In: Proc. IEEE Workshop on Statistical Sig. Array Proc.; 1996. pp. 97-100
  2. 2. Bermani E, Boni A, Caorsi S, Massa A. An innovative real-time technique for buried object detection. IEEE Transactions on Geoscience and Remote Sensing. 2003;41:927-931
  3. 3. Caorsi S, Anguita D, Bermani E, Boni A, Donelli M, Massa A. A comparative study of NN and SVM-based electromagnetic inverse scattering approaches to on-line detection of buried objects. Journal of the Applied Computational Electromagnetics Society. Jul. 2003;18:1-11
  4. 4. Massa A, Boni A, Donelli M. A classification approach based on SVM for electromagnetic subsurface sensing. IEEE Transactions on Geoscience and Remote Sensing. Sept. 2005;43(9):2084-2093
  5. 5. Zhang B, O’Neill K, Kong JA, Grzegorczyk TM. Support vector machine and neural network classification of metallic objects using coefficients of the spheroidal MQS response modes. IEEE Transactions on Geoscience and Remote Sensing. Jan. 2008;46(1):159-171
  6. 6. Ho KC, Carin L, Gader PD, Wilson JN. An investigation of using the spectral characteristics from ground penetrating radar for landmine/clutter discrimination. IEEE Transactions on Geoscience and Remote Sensing. Apr. 2008;46(4):1177-1191
  7. 7. Guillermin R, Lasaygues P, Sessarego JP, Wirgin A. Characterization of buried objects by a discretized domain integral equation inversion method using born approximation. In: Proc. 5th European Conference of Underwater Acoustics. Vol. 2; Jul. 10-13, 2000. pp. 863-868
  8. 8. Hetet A, Amate M, Zerr B. SAS processing results for the detection of buried objects with a ship-mounted sonar. In: the 7th European Conference on Underwater Acoustics (ECUA ’04); Delft, The Netherlands; July. 2004. pp. 1127-1132
  9. 9. Kumareasan R, Tufts DW. Estimating the angles of arrival of multiple source plane waves. IEEE Transactions on Aerospace and Electronic Systems. 1983;19:134-139
  10. 10. Baggeroer AB, Kuperman WA, Mikhalevsky PN. An overview of matched field methods in ocean acoustics. IEEE Journal of Oceanic Engineering. 1993;18(4):401-424
  11. 11. Sahin A, Miller EL. Object-based localization of buried objects using high resolution array processing techniques. Proc. SPIE Aerosense. May 1996;2765
  12. 12. Fawcett JA, Lim R. Evaluation of the integrals of target/seabed scattering using the method of complex images. The Journal of the Acoustical Society of America. Sep. 2003;114(3):1406-1415
  13. 13. Fawcett JA, Fox WLJ, Maguer A. Modeling of scattering by objects on the seabed. The Journal of the Acoustical Society of America. Dec. 1998;104(6):3296-3304
  14. 14. Bourennane S, Bendjama A. Locating wide band acoustic sources using higher-order statistics. Applied Acoustics. 2002;63:235-251
  15. 15. Saidi Z, Bourennane S. Cumulant-based coherent signal subspace method for bearing and range estimation. EURASIP Journal on Advances in Signal Processing. 2007. p. 9. Article ID 84576. Doi:10.1155/2007/84576
  16. 16. Doolittle RD, Uberall H. Sound scattering by elastic cylindrical shells. The Journal of the Acoustical Society of America. 1966;39(2)
  17. 17. Ye Z. Recent developments in underwater acoustics: Acoustic scattering from single and multiple bodies. Proceedings of the National Science Council, Republic of China, Part A: Physical Science and Engineering. 2001;25(3):137-150
  18. 18. Lim R. Acoustic scattering by a partially buried three dimensional elastic obstacle. The Journal of the Acoustical Society of America. 1998;104(2):769-782
  19. 19. Tesei A, Maguer A, Fox WLJ, Schmidt H. Measurements and modeling of acoustic scattering from partially and completely buried spherical shells. The Journal of the Acoustical Society of America. Nov. 2002;112(5):1817-1830
  20. 20. Bourennane S, Han D, Fossati C. Fast coherent signal subspace-based method for bearing and range of buried objects estimation in the presence of colored noise. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2009;2:329-338
  21. 21. Tayem N, Kwon HM. Azimuth and elevation angle estimation with no failure and no eigen decomposition. EURASIP Signal Processing. Jan. 2006;86(1):8-16
  22. 22. Bourennane S, Fossati C, Marot J. About noneigenvector source localization methods. EURASIP Journal on Advances in Signal Processing. 2008 Article ID 480835
  23. 23. Hyvärinen A, Oja E. A fast fixed-point algorithm for independent component analysis. Neural Computation. Oct. 1997;9(7):1483-1492
  24. 24. Yang B. Projection approximation subspace tracking. IEEE Transactions on Signal Processing. Jan. 1995;44(1):95-107
  25. 25. Yang B. Asymptotic convergence analysis of the projection approximation subspace tracking algorithms. Signal Processing. Apr. 1996;50:123-136
  26. 26. Yang B. An extension of the PASTd algorithm to both rank and subspace tracking. IEEE Signal Processing Letters. Mar. 2000;7(3):60-62
  27. 27. Golub GH, Van Loan CF. Matrix Computations. third ed. Baltimore and London: The Johns Hopkins University Press; 1996
  28. 28. Hua Y, Xiang Y, Chen T, Abed-Meraim K, Miao Y. A new look at the power method for fast subspace tracking. Digital Signal Processing. Oct. 1999:297-314
  29. 29. Abed-Meraim K, Chkeif A, Hua Y. Fast orthogonal PAST algorithm. IEEE Signal Processing Letters. 1995;2:179-182
  30. 30. Badeau R, Abed-Meraim K, Richard G, David B. Sliding window orthonormal PAST algorithm. In: IEEE Int. Conf. on Acoustic, Speech and Signal Processing; 2003
  31. 31. Gonen E, Mendel JM. Applications of cumulants to array processing - Part III: Blind beamforming for coherent signals. IEEE Transactions on Signal Processing. Sept. 1997;45(9):2252-2264
  32. 32. Mendel JM. Tutorial on higher-order statistics (spectra) in signal processing and system theory: Theoretical results and some applications. IEEE-Proceedings. Mar. 1991;79(3):278-305
  33. 33. Villemin G, Fossati C, Bourennane S. Efficient time of arrival estimation in the presence of multipath propagation. The Journal of the Acoustical Society of America (JASA). Oct. 2013;134(4):301-306
  34. 34. Villemin G, Fossati C, Bourennane S. Spatio-temporal-based joint range and angle estimation for wideband signals. EURASIP Journal on Advances in Signal Processing. 2013. Article ID 131

Written By

Salah Bourennane and Caroline Fossati

Submitted: 13 June 2017 Reviewed: 26 September 2017 Published: 22 November 2017