Open AccessArticle

Underwater Acoustic Matched Field Imaging Based on Compressed Sensing

Huichen Yan

¹,

Jia Xu

^2,*

Teng Long

² and

Xudong Zhang

Department of Electronic Engineering, Tsinghua University, Beijing 100084, China

School of Information and Electronics, Beijing Institute of Technology, Beijing 100081, China

Author to whom correspondence should be addressed.

Sensors 2015, 15(10), 25577-25591; https://doi.org/10.3390/s151025577

Submission received: 4 August 2015 / Revised: 6 September 2015 / Accepted: 18 September 2015 / Published: 7 October 2015

(This article belongs to the Section Physical Sensors)

Download

Browse Figures

Figure 1
(a) The sketch of underwater circumstances; (b) The pressure field intensity generated by a source. "> Figure 2
(a) Conventional MFP model Conventional MFP model; (b) CS-MFP model. "> Figure 3
(a) A real part of the matrix A example; (b) A’s average correlation coefficients. "> Figure 4
Percentage of high coherent atoms (coefficient ≥ 0.5) in the dictionary. "> Figure 5
(a) Relative imaging errors w.r.t. numbers of targets under different depth resolutions with ρr = 2 m fixed; (b) Relative imaging errors w.r.t. numbers of targets under different range resolutions with ρz = 2 m fixed. "> Figure 6
(a) Recovery errors w.r.t. numbers of targets and sensors with CS-MFP model; (b) Recovery errors w.r.t. numbers of targets and sensors with traditional MFP model. "> Figure 7
(a) A Gram matrix of <math display="inline"> <semantics> <mi mathvariant="script">A</mi> </semantics> </math> in a multi-frequency case; (b) The targets of interest; (c) The recovery result. "> Figure 8
The successful recovery rate in terms of SNR. ">

Versions Notes

Abstract

Matched field processing (MFP) is an effective method for underwater target imaging and localizing, but its performance is not guaranteed due to the nonuniqueness and instability problems caused by the underdetermined essence of MFP. By exploiting the sparsity of the targets in an imaging area, this paper proposes a compressive sensing MFP (CS-MFP) model from wave propagation theory by using randomly deployed sensors. In addition, the model’s recovery performance is investigated by exploring the lower bounds of the coherence parameter of the CS dictionary. Furthermore, this paper analyzes the robustness of CS-MFP with respect to the displacement of the sensors. Subsequently, a coherence-excluding coherence optimized orthogonal matching pursuit (CCOOMP) algorithm is proposed to overcome the high coherent dictionary problem in special cases. Finally, some numerical experiments are provided to demonstrate the effectiveness of the proposed CS-MFP method.

Keywords:

matched field processing (MFP); compressed sensing (CS); wave propagation; coherence parameter; coherence-excluding coherence optimized orthogonal matching pursuit (CCOOMP)

1. Introduction

Matched field processing (MFP) [1,2,3] is an attractive technique for imaging and localizing underwater targets, and its effectiveness has been widely verified in many real scenarios. Without emitting signals, it is known that the underwater targets can be recovered via MFP in the range-depth plane by using a sensor array passively. MFP collects sound pressure field data from targets and background and compares the data in the Bartlett processor with a series of pre-calculated data of all possible locations from models with greedy algorithms to finally determine the existence and positions of targets [4,5]. However, due to the limited number of sensors, the small array size and the complicated background, the recovery performance of existing MFP is not guaranteed in real application. Normally, the false MFP recoveries may occur on both number and positions for the targets.

These problems above are primarily owing to the fact that the MFP solves an underdetermined problem in practice [6]. To deal with this problem, the time-reversal technique [7] has been proposed as an alternative, but it also needs to be improved. As most MFP recovers a very small number of targets, this reminds us of the compressed sensing (CS) technique. To solve the underdetermined problem, the “small number”, or “sparsity”, property of the underwater targets in the scene should be effectively exploited. During the past decade, CS theories and techniques have been progressing at an overwhelming pace. For an underdetermined linear problem y = Ax, where y is a given observation and A is the so-called dictionary with far more columns than rows and each of its columns called atoms, it is known that the CS theory can provide a unique and stable solution as long as x is sparse and A satisfies the restricted isometrics property (RIP) condition. The RIP condition can be characterized by its coherence parameter [8,9,10,11], and the coherence parameter of A can be viewed as the largest normalized cross-correlation coefficient of the atom pairs. As computational imaging problems are mostly underdetermined, CS methods have been widely discussed and applied [6,12,13,14,15] in cases where imaging scene is sparse. CS imaging has been also proposed for MFP in [16], where it focused more on the computational complexity reduction rather than the imaging model and recovery performance.

We give a performance analysis of MFP from the view of CS. Instead of accelerating computation speed in Bartlett processing, in this paper, CS is hinged in the physical propagation model as [6] to improve the MFP performance. We focus on the shallow sea case and exploit the normal-mode acoustic propagation model to analyze the seabed and surface effects. In this case, the water is viewed as isotropic medium and the seabed and surface act as sound propagation boundaries. Boundaries affect the wave propagation, resulting in the problems of multi-path effects, interference effects, “blind spot effect” and so on [17]. According to wave propagation theory, the partial differential equation (PDE) is used for describing the propagation of sound in water from the wave equation, and the boundaries determine the constraints of the PDE directly. Based on the normal-mode acoustic propagation model, the shallow water is then regarded as a waveguide and a monochromatic wave is described as the sum of waves with different wave modes. Therefore, the MFP can be related to CS theory in essence and the CS-MFP model is proposed. In addition, this model can be easily transformed from the monochromatic one into a multi-frequency case.

In accordance with the RIP condition of CS theory, it is required that the atoms should be approximately orthogonal, or the coherence parameter of different atoms should be as small as possible. However, the dictionary in CS-MFP model cannot be designed arbitrarily. Furthermore, it is found that the dictionaries in CS-MFP are highly coherent by verifying the lower bound coherence parameter. This problem is related to the well-known off-grid problem and spectral CS (SCS) problem in CS recovery. In these problems, the nonzero elements of targets may not be exactly on the predefined grids, which will inevitably cause recovery errors. Nevertheless, refining the grids will generate the high coherent dictionaries. To solve this above problems, [18] and [19] discussed the off-grid CS recovery, and [20,21] proposed Spectral Iterative Hard Thresholding (SIHT) and Band-excluded Locally Optimized Orthogonal Matching Pursuit (BLOOMP) algorithms, respectively. BLOOMP can use fewer samples yet performed far better than SIHT. On the other hand, more complicated methods have been proposed recently for Fourier dictionaries [22,23]. In this paper, as the dictionary in CS-MFP is more complicated, we modify the BLOOMP method and provide a coherence-excluding coherence optimized orthogonal matching pursuit (CCOOMP) algorithm to solve the high coherent dictionary problem.

The contribution of this paper is as follows,

(1): The model of CS-MFP is established from wave propagation angle considering boundaries;
(2): The recovery performance of the model is discussed from a CS angle by examining the lower bound of the CS coherent parameter of the above model;
(3): The effective solution is given when targets are sparsely distributed.

The remainder of this paper is arranged as follows. In Section 2, the CS-MFP model is established for monochromatic wave cases. In Section 3, we propose a CCOOMP algorithm to solve the CS-MFP problem with high coherent dictionaries. In Section 4, we present some numerical experiments to verify the model and algorithm and extend the model to multi-frequency case. In Section 5, some conclusions are drawn.

2. Compressed Sensing MFP Model Formulation

This section describes the formulation of underwater acoustic CS-MFP imaging model in a shallow sea case from the angle of wave propagation. To articulate the model clearly, we consider simple boundaries in the formulation.

2.1. Wave Propagation Theories

In propagation theory, a monochromatic acoustic wave generated from a point source in a boundless field will propagate in a form of spherical wave in the range-depth plane (r–z plane) according to the wave equation. As Figure 1a indicates, in a shallow sea case, assuming that sound propagation velocity c is a constant, sea surface (z = 0) is a soft boundary and seabed (z = H) is a hard boundary, the propagation of the acoustic wave can be modeled as a waveguide model in r–z plane. The propagation follows a wave equation and the two constraints for the wave equation are

Z (n) |_{z = 0} = 0, {\frac{d Z (n)}{d z} |}_{z = H} = 0, 1 \leq n \leq N,

(1)

where

Z (n)

is the eigenfunction and

N

is the number of normal modes of the wave propagating in the water. Constrained by Equation (1), an acoustic source at

r = r_{0}, z = z_{0}

, denoted as (

r_{0}

z_{0}

), with unit strength and an angle frequency

ω

will generate acoustic pressure field at a far-field position ( r, z ) as

p (r, z) = - \frac{2 j}{H} \sum_{n = 1}^{N} {\sqrt{\frac{2 π}{ς (n) (r - r_{0})}} s i n (k_{z} (n) z) s i n (k_{z} (n) z_{0}) \cdot e x p [- j (ς (n) (r - r_{0}) - \frac{π}{4})]}

(2)

in which,

j = \sqrt{- 1}

N = ⌊ \frac{H ω}{π c} + \frac{1}{2} ⌋

(3)

where

⌊ \cdot ⌋

is the ceiling operator and

k_{z}

ς

are vectors of length N representing the eigenvalues of wave modes,

k_{z} (n) = (n - \frac{1}{2}) \frac{π}{H}

(4)

and

ς (n) = \sqrt{{(\frac{ω}{c})}^{2} - k_{z} {(n)}^{2}}

(5)

It can be shown that

ς

and

k_{z}

are the horizontal and vertical components of the wave number ω/c. Figure 1b presents an example of the generated sound pressure field intensity from a source at position (0, z₀) shown in Figure 1a as the point in the left. The reader can be referred to [24] for a full presentation of propagation theory.

Figure 1. (a) The sketch of underwater circumstances; (b) The pressure field intensity generated by a source.

2.2. CS-MFP Model Formulation

As Figure 2a shows, in the conventional MFP, a two-dimensional observation area of interest with hundreds of “candidate” pixels is determined in advance. For simplicity, we assume the targets only exist on a square lattice of dimension Q × Q with grid spacing

Δ r

and

Δ z

and denote the grid points’ coordinates as

r_{T}

z_{T}

). Hence,

r_{T} (i) = r_{T} (1) + (i - 1) Δ r,

z_{T} (i) = z_{T} (1) + (i - 1) Δ z

i = 1, 2, \cdot \cdot \cdot, Q

. Matrix

G

denotes the strengths of the grid points to be determined, i.e.,

G (i, j) (i, j = 1, 2, \cdot \cdot \cdot, Q)

represents the strength at coordinate

(r_{T} (i), z_{T} (j))

. A uniformly linear sensor array with

P

elements is set in the far field at positions

(r_{S}, z_{S})

to collect data denoted as

y

. The sensors can receive all the existing sound pressure generated by these targets. If we set

(r, z) = (r_{S}, z_{S})

and

(r_{0}, z_{0}) = (r_{T}, z_{T})

in Equation (2), then

p (r, z)

can be also represented as

p (r_{S}, z_{S} | r_{T}, z_{T})

. Hence,

y

is a linear superposition of the entire sound pressure intensities generated from all the existing candidate pixels according to Equation (2) as

y = A \cdot vec (G) + e

(6)

where

A

represents the discrete version of Equation (2) and

e

is an observation error. With slight abuse of notations we have

A (i, j) = p (r_{S} (i), z_{S} (i) | r_{T} (⌈ j / Q ⌉), z_{T} (j - (⌈ j / Q ⌉ - 1) Q))

(7)

where

i = 1, 2, \cdot \cdot \cdot, P, j = 1, 2, \cdot \cdot \cdot, Q^{2}

, vec(·) is a function piling a matrix to a vector column by column. Denoting

{(\cdot)}^{T}

as matrix transpose, we have

vec (G) = {[G {(:, 1)}^{T} G {(:, 2)}^{T} \dots G {(:, Q)}^{T}]}^{T}

(8)

Clearly, Equation (6) can be viewed as the weighted sum of the atoms in

A

with the weights representing the strengths of the targets.

Figure 2. (a) Conventional MFP model Conventional MFP model; (b) CS-MFP model.

In the conventional MFP framework, Bartlett processing is adopted therefrom to determine the most possible target positions among the candidate pixel positions. As the pixel number is greater compared to the number of sensors, it is obvious that the problem is underdetermined. Consequently, we obtain a linear model with the form of CS as

y = A x + e

(9)

In CS theory, A can be Gaussian, Bernoulli and partial Fourier random matrices, as they satisfy the RIP condition [10]. The randomness in these matrices plays a crucial role. In this new CS-MFP model, we propose to use a random distributed sensor array in a 2D plane instead of the traditional uniformly distributed linear array as Figure 2b plots. We will see the advantage of this setup of matrix A from the angle of compressed sensing in the following analysis, especially for resolution in

r

direction. In practice, the random array is also applicable. In fact, we can randomly deploy the sensors and keep track of their positions separately just as a linear array. In addition, we don’t have to stick to the uniform distance requirement. However, as Equation (7) indicates, it is not trivial to analyze the RIP or even for the coherence parameter of

A

. To comprehend the characteristics of

A

in Equation (7) to certain extent, let us separately observe each of the

N

components in Equation (2), i.e.,

A = \sum_{n = 1}^{N} A_{n}

(10)

and analyze the lower bound of the coherence parameter. Note in reality the separation is also possible by employing sound pressure sensor subarrays to form spatial bandpass filters for certain bandwidth at particular observation locations to form space filters. Denote the coherence parameter of

A

μ

. Here we use

i, j

instead of

n

to distinguish two subscriptions. According to the definition, the Gram matrix of A is

Gram (A) = A^{H} A = \sum_{i = 1}^{N} A_{i}^{H} \sum_{j = 1}^{N} A_{j} = \sum_{i = 1}^{N} A_{i}^{H} A_{i} + \sum_{i, j = 1, i \neq j}^{N} A_{i}^{H} A_{j}

(11)

The coherence parameter of A is the maximum of the off-diagonal entry of Gram (A). Let us firstly analyze the first part of the right hand side of Equation (11). For a fixed

i

, two targets with the same

z_{0}

but a

Δ r

distance in

r

direction will generate two atoms with a correlation coefficient as

\begin{array}{l} μ_{r} (i, i) \\ = \frac{| {[\frac{1}{\sqrt{r}} s i n (k_{z} (i) z) s i n (k_{z} (i) z_{0}) e x p (- j ς (i) r)]}^{H} [\frac{1}{\sqrt{r + 1 Δ r}} s i n (k_{z} (i) z) s i n (k_{z} (i) z_{0}) e x p (- j ς (i) (r + 1 Δ r))] |}{| \frac{1}{\sqrt{r}} s i n (k_{z} (i) z) s i n (k_{z} (i) z_{0}) e x p (- j ς (i) r) | | \frac{1}{\sqrt{r + 1 Δ r}} s i n (k_{z} (i) z) s i n (k_{z} (i) z_{0}) e x p (- j ς (i) r + 1 Δ r) |} \\ = \frac{\sum_{p = 1}^{P} \frac{1}{\sqrt{r_{p} (r_{p} + Δ r)}} {s i n}^{2} (k_{z} (i) z_{p})}{\sqrt{\sum_{p = 1}^{P} \frac{1}{r_{p}} {s i n}^{2} (k_{z} (i) z_{p})} \sqrt{\sum_{p = 1}^{P} \frac{1}{r_{p} + Δ r} {s i n}^{2} (k_{z} (i) z_{p})}} \end{array}

(12)

On the other hand, two targets with the same

r

but a

Δ z

distance in

z

direction will generate two atoms with a correlation coefficient

\begin{array}{l} μ_{z} (i, i) \\ = \frac{| {[\frac{1}{\sqrt{r}} s i n (k_{z} (i) z) s i n (k_{z} (i) (z_{0} + Δ z)) e x p (- j ς (i) r)]}^{H} [\frac{1}{\sqrt{r}} s i n (k_{z} (i) z) s i n (k_{z} (i) z_{0}) e x p (- j ς (i) r)] |}{| [\frac{1}{\sqrt{r}} s i n (k_{z} (i) z) s i n (k_{z} (i) (z_{0} + Δ z)) e x p (- j ς (i) r)] | | [\frac{1}{\sqrt{r}} s i n (k_{z} (i) z) s i n (k_{z} (i) z_{0}) e x p (- j ς (i) r)] |} \\ = \frac{| \sum_{p = 1}^{P} \frac{1}{r_{p}} {s i n}^{2} (k_{z} (i) z_{p}) s i n (k_{z} (i) z_{0}) s i n (k_{z} (i) (z_{0} + Δ z)) |}{\sqrt{\sum_{p = 1}^{P} \frac{1}{r_{p}} {s i n}^{2} (k_{z} (i) z_{p}) {s i n}^{2} (k_{z} (i) (z_{0} + Δ z))} \sqrt{\sum_{p = 1}^{P} \frac{1}{r_{p}} {s i n}^{2} (k_{z} (i) z_{p}) {s i n}^{2} (k_{z} (i) z_{0})}} = 1. \end{array}

(13)

Now let us focus on the second part of Equation (11). Similar to the above discussion, we have for two matrix

A_{i}

and

A_{j}

(

i \neq j

), two targets with the same

r

and a

Δ z

distance in the

z

direction,

\begin{array}{l} μ_{z} (i, j) \\ = \frac{| \sum_{p = 1}^{P} \frac{1}{\sqrt{ς (i) ς (j)} r_{p}} s i n (k_{z} (i) z_{p}) s i n (k_{z} (j) z_{p}) s i n (k_{z} (i) z_{0}) s i n (k_{z} (j) (z_{0} + Δ z)) e x p (- j [ς (i) - ς (j)] r_{p}) |}{\sqrt{\sum_{p = 1}^{P} \frac{1}{ς (i) r_{p}} {s i n}^{2} (k_{z} (j) z_{p}) {s i n}^{2} (k_{z} (j) z_{0})} \sqrt{\sum_{p = 1}^{P} \frac{1}{ς (j) r_{p}} {s i n}^{2} (k_{z} (i) z_{p}) {s i n}^{2} (k_{z} (i) (z_{0} + Δ z))}} \\ = \frac{| \sum_{p = 1}^{P} \frac{1}{r_{p}} s i n (k_{z} (i) z_{p}) s i n (k_{z} (j) z_{p}) e x p (- j [ς (i) - ς (j)] r_{p}) |}{\sqrt{\sum_{p = 1}^{P} \frac{1}{r_{p}} {s i n}^{2} (k_{z} (i) z_{p})} \sqrt{\sum_{p = 1}^{P} \frac{1}{r_{p}} {s i n}^{2} (k_{z} (j) z_{p})}} \end{array}

(14)

For two targets with the same

z

and a

Δ r

distance in the

r

direction,

\begin{array}{l} μ_{r} (i, j) \\ = \frac{| \sum_{p = 1}^{P} \frac{1}{\sqrt{ς (i) ς (j) r_{p} (r_{p} + Δ r)}} s i n (k_{z} (i) z_{p}) s i n (k_{z} (j) z_{p}) s i n (k_{z} (i) z_{0}) s i n (k_{z} (j) z_{0}) e x p (- j [ς (i) - ς (j)] r_{p}) |}{\sqrt{\sum_{p = 1}^{P} \frac{1}{ς (i) r_{p}} {s i n}^{2} (k_{z} (i) z_{p}) {s i n}^{2} (k_{z} (i) z_{0})} \sqrt{\sum_{p = 1}^{P} \frac{1}{ς (j) (r_{p} + Δ r)} {s i n}^{2} (k_{z} (j) z_{p}) {s i n}^{2} (k_{z} (i) z_{0})}} \\ = \frac{| \sum_{p = 1}^{P} \frac{1}{\sqrt{r_{p} (r_{p} + Δ r)}} s i n (k_{z} (i) z_{p}) s i n (k_{z} (j) z_{p}) e x p (- j [ς (i) - ς (j)] r_{p}) |}{\sqrt{\sum_{p = 1}^{P} \frac{1}{r_{p}} {s i n}^{2} (k_{z} (i) z_{p})} \sqrt{\sum_{p = 1}^{P} \frac{1}{(r_{p} + Δ r)} {s i n}^{2} (k_{z} (j) z_{p})}} \end{array}

(15)

Let us consider the case where z_p is uniform distribution on [0, H] and r_p is uniform distribution on [r_min, r_max]. Then we have for large P,

\begin{array}{l} \sum_{p = 1}^{P} \frac{1}{\sqrt{r_{p} (r_{p} + Δ r)}} {s i n}^{2} (k_{z} (i) z_{p}) = \sum_{p = 1}^{P} \frac{1}{2 \sqrt{r_{p} (r_{p} + Δ r)}} [1 - c o s (2 k_{z} (i) z_{p})] \\ = \int_{r = r_{m i n}}^{r_{m a x}} \int_{z = 0}^{H} \frac{1}{2 \sqrt{r (r + Δ r)}} [1 - c o s (π \frac{2 i - 1}{H} z)] d z d r \approx \frac{H}{2} \ln (\frac{r_{m a x} + \frac{1}{2} Δ r}{r_{m i n} + \frac{1}{2} Δ r}) \end{array}

(16)

\begin{array}{l} \sqrt{\sum_{p = 1}^{P} \frac{1}{r_{p}} {s i n}^{2} (k_{z} (i) z_{p})} = \sqrt{\sum_{p = 1}^{P} \frac{1}{2 r_{p}} [1 - c o s (2 k_{z} (i) z_{p})]} \\ = \sqrt{\int_{r = r_{m i n}}^{r_{m a x}} \int_{z = 0}^{H} \frac{1}{2 r} [1 - c o s (π \frac{2 i - 1}{H} z)] d z d r} = \sqrt{\frac{H}{2} \ln (\frac{r_{m a x}}{r_{m i n}})} \end{array}

(17)

\begin{array}{l} \sqrt{\sum_{p = 1}^{P} \frac{1}{r_{p} + Δ r} {s i n}^{2} (k_{z} (i) z_{p})} = \sqrt{\sum_{p = 1}^{P} \frac{1}{2 (r_{p} + Δ r)} [1 - c o s (2 k_{z} (i) z_{p})]} \\ = \sqrt{\int_{r = r_{m i n}}^{r_{m a x}} \int_{z = 0}^{H} \frac{1}{2 (r + Δ r)} [1 - c o s (π \frac{2 i - 1}{H} z)] d z d r} = \sqrt{\frac{H}{2} \ln (\frac{r_{m a x} + Δ r}{r_{m i n} + Δ r})} \end{array}

(18)

According to Equations (16)–(18),

μ_{r} (i, i) = \frac{\ln (r_{m a x} + \frac{1}{2} Δ r) - \ln (r_{m i n} + \frac{1}{2} Δ r)}{\sqrt{\ln (\frac{r_{m a x}}{r_{m i n}}) \ln (\frac{r_{m a x} + Δ r}{r_{m i n} + Δ r})}}

(19)

in Equation (12). Similarly,

μ_{z} (i, j)=0

and

μ_{r} (i, j)=0

in this case.

From the above discussion, we have four conclusions.

(1): the distance $Δ z$ has no impact on the coherence parameter.
(2): the random deployment of sensors helps lower the coherence parameter in Equations (12) and (15), according to Cauchy-Schwarz inequality. Traditionally, if we use uniform array, we will have ${\tilde{μ}}_{r} (i, i) = 1$ in Equation (12) and greater ${\tilde{μ}}_{r} (i, j)$ as

${\tilde{μ}}_{r} (i, j) = \sum_{p = 1}^{P} s i n (k_{z} (i) z_{p}) s i n (k_{z} (j) z_{p}) / \sqrt{\sum_{p = 1}^{P} {s i n}^{2} (k_{z} (i) z_{p})} \sqrt{\sum_{p = 1}^{P} {s i n}^{2} (k_{z} (j) z_{p})}$
(3): in the special case, the high coherent part is mainly contributed by signals of the same modes.
(4): if all the $A_{n}$ ( $n = 1, 2, \dots, N$ ) are compensated to have equal energy, we have a lower bound

$\begin{array}{l} μ \geq m a x (μ_{r}, μ_{z}) \\ = \frac{1}{N} m a x (μ_{r} (i, i), μ_{z} (i, i)) = \frac{1}{N} \end{array}$

(20)

According to [25], if the target number is less than

(μ^{- 1} + 1) / 2,

(21)

then

x

in Equation (9) can be recovered uniquely using greedy algorithms. This Theorem determines the performance of CS-MFP model.

2.3. Robustness of Coherence

In this subsection, we discuss the robustness of the coherence parameter. We analyze the correlation coefficient change w.r.t. the displacement of the sensors. As the sensors are deposed in the flowing water, this discussion becomes necessary. When there is a

Δ z

displacement in the

z

direction, for the

n

th component in Equation (2), the correlation coefficient between the new measurement and the theoretical measurement can be expressed as

μ_{z} = \frac{| {[sin (k_{z} (n) z)]}^{T} [sin (k_{z} (n) (z + 1 Δ z))] |}{{‖ sin (k_{z} (n) z) ‖}_{2} {‖ sin (k_{z} (n) (z + 1 Δ z)) ‖}_{2}}

(22)

When vector

z

is set randomly, the sensors are uniformly distributed in the

z

direction, we have

μ_{z}^{2} = \frac{{[\sum_{p = 1}^{P} cos [2 k_{z} (n) z_{p} + k_{z} (n) Δ z] - P cos (k_{z} (n) Δ z)]}^{2}}{[P - \sum_{p = 1}^{P} cos (2 k_{z} (n) z_{p})] [P - \sum_{p = 1}^{P} cos (2 k_{z} (n) (z_{p} + Δ z))]}

(23)

The summation is approximated as integral for large

P

. Hence, the numerator can be simplified as

\begin{array}{l} P^{2} {cos}^{2} (k_{z} (n) Δ z) - \frac{P cos (k_{z} (n) Δ z)}{k_{z} (n)} {s i n [2 k_{z} (n) H + k_{z} (n) Δ z] - s i n [k_{z} (n) Δ z]} \\ + \frac{2 c o s [2 k_{z} (n) H + 2 k_{z} (n) Δ z] - c o s [4 k_{z} (n) H + 2 k_{z} (n) Δ z] - c o s [2 k_{z} (n) Δ z] - 2 c o s [2 k_{z} (n) H] + 2}{8 k_{z} {(n)}^{2}} \end{array}

(24)

And the denominator is

\begin{array}{l} P^{2} - \frac{P}{2 k_{z} (n)} s i n (2 k_{z} (n) H) - \frac{P}{2 k_{z} (n)} s i n (2 k_{z} (n) H) \\ + \frac{1}{4 k_{z} {(n)}^{2}} s i n (2 k_{z} (n) H) [s i n (2 k_{z} (n) (H + Δ z)) - s i n (2 k_{z} (n) Δ z)] \end{array}

(25)

Hence, we have

μ_{z} = | c o s (k_{z} (n) Δ z) |

(26)

This result indicates the robustness deteriorates with the increase of

n

. On the other hand, the robustness w.r.t.

r

has the same form of Equation (12).

2.4. Multi-Frequency CS-MFP Case

In model Equation (6), the angle frequency of the targets is fixed as

ω

. It is a reasonable assumption as in real applications the bandwidths of the transducers are relatively narrow. However, in some wide band or multi-frequency cases, the model can be extended. Besides a two-dimensional search in

r

z

space, frequency dimension is added as the third dimension, whose dimensionality is denoted as

R

. In this way, we obtain a more intricate matrix

A

to replace

A

in Equation (9) with

A = [A_{ω_{1}} A_{ω_{2}} \dots A_{ω_{R}}]

Similarly, we can add even more dimensions to the model, for instance, a time series search. However these dimensions deteriorate the RIP of the dictionary, making it more difficult to accurately recover the targets with the same number of sensors.

3. CS-MFP Recovery with High Coherent Dictionaries via CCOOMP

However, as it can be seen from Equations (12)–(15) and the simulations in Section 4, in general, it is very difficult to give an explicit expression of

μ

. It also may be high. This is a common problem faced in CS imaging applications, as in conventional CS theories matrix

A

can be meticulously designed, whereas in imaging applications,

A

is linked to physical quantities without much design freedom as Equation (6) shows. This problem resembles an off-grid problem in CS and spectral CS problems. Off-grid problems occur universally in CS imaging applications [13,14,15,21,26], as in reality, targets are not distributed exactly as the lattice is divided, or a discretized dictionary cannot fully express an analog signal

\hat{x}

from Equation (9). BLOOMP is a successful tool to solve this problem. In these problems, high coherent atoms form a neighborhood or the so-called “band”, whereas in this problem, the high coherent atoms are randomly distributed in the dictionary (see Figure 1b and Figure 3b). We then modify the BLOOMP method to solve this problem.

Figure 3. (a) A real part of the matrix A example; (b) A’s average correlation coefficients.

Assuming

0 < η < 1

, we first define the η-coherence index set of atom

k

to be the set

C_{η} (k) = {i | μ (i, k) > η}

(27)

and the η-coherence index set of an index set

S

is hence defined as

C_{η} (S) = \cup_{k \in S} C_{η} (k)

(28)

and

C_{η}^{2} (S) = C (C_{η} (S))

(29)

Here we assume the targets are not distributed on positions with correlation coefficient larger than η, or for any

i \in supp {x}

i \cap C_{η}^{2} (supp (x) \ i) = \emptyset

(30)

Then, we obtain the Coherence-excluding Coherence Optimized Orthogonal Matching Pursuit (CCOOMP) method as Algorithm 1 shows. CCOOMP is similar to BLOOMP [21] except in the implementing, CCOOMP virtually excludes the atoms according to the correlation coefficients rather than the physical distances. Coherence Optimization (CO) (Algorithm 2) is embedded in CCOOMP.

Algorithm 1 Coherence-Excluding Coherence Optimized Orthogonal Matching Pursuit (CCOOMP)

1: Input:  $A, y, η > 0, s, t, e p$ 
2: Initialization:  $x_{0} = 0, r = y$  and  $S^{0} = Ø$ 
3: Iteration: For  $n = 1, \dots, s$ 
4:  $i_{m a x} = \arg \underset{i}{m a x} | 〈 r^{n - 1}, A (:, i) 〉 |, i \notin C_{η}^{2} (S^{n - 1})$ 
5:  $S^{n} = CO (S^{n - 1} \cup {i_{m a x}})$ 
6:  $x^{n} = \arg \underset{z}{m i n} {‖ A z - y ‖}_{2} s . t . supp (z) \in S^{n}$ 
7:  $r^{n} = y - A x^{n}$ 
8: return:  $x^{s}$

Algorithm 2 Coherence Optimization (CO)

1: Input:  $A, y, η > 0, S^{0} = {i_{1}, i_{2}, \dots, i_{k}}$ 
2: Iteration: For  $n = 1, \dots, k$ 
3:  $x^{n} = \arg \underset{i}{m a x} {‖ Ax - b ‖}_{2}, supp (x) = (S^{n - 1} \ {i_{n}} \cup {j_{n}}), j_{n} \in C_{η} ({i_{n}})$ 
4:  $S^{n} = supp (x^{n})$ 
5: return:  $S^{k}$

4. Numerical Experiments and Performance Analysis

In this section, numerical experiments are conducted to demonstrate the proposed model and algorithm. We still consider a shallow sea case. In this case H = 200 m. We are going to observe an area at a distance r

\in

[100 m, 132 m] and depth z

\in

[2 m, 17 m] with 16 × 16 pixels. Range resolution

ρ_{r}

is set 2 m, depth resolution ρz = 1 m. We uniformly randomly place 80 transducers along r

\in

[0 m, 10 m] and z

\in [0, H]

, and the central frequency of the transducers is 6 kHz. Note according to Equation (3), N = 25. Then matrix

A

is determined. Figure 3a illuminates an example of real part of matrix

A

. In Figure 3b we calculate the normalized correlation coefficients of each atom with all the other atoms in the dictionary and average it. In this case, high correlation coefficients occur randomly rather than forming a neighbor band around 0 as in SCS. This is a result of the scattering of the boundaries.

To obtain a better understanding of the dictionary, we analyze the coherence property of the matrix w.r.t., the resolution in r–z dimensions. In this particular experiment, the average percentages of high coherent atoms are drawn in Figure 4 after 10 iterations. As can be expected, the pairs of high coherent atoms increase with the resolution of imaging. Now, from the above experiment, we can see the high coherence of the dictionaries.

Figure 4. Percentage of high coherent atoms (coefficient ≥ 0.5) in the dictionary.

In the next experiment, the CCOOMP algorithm is applied to the imaging process. The sparse targets are set randomly in the observing area with coherence no greater than 0.5. For comparison, Figure 5 shows the imaging errors with CCOOMP in terms of number of targets under different resolutions. For each number of targets, 100 simulations are performed and the average error is plotted. It is seen that by properly setting the coherent parameter, the recovery result can be satisfactory.

Figure 5. (a) Relative imaging errors w.r.t. numbers of targets under different depth resolutions with ρ_r = 2 m fixed; (b) Relative imaging errors w.r.t. numbers of targets under different range resolutions with ρ_z = 2 m fixed.

In the next experiment, the numbers of sensors and targets change simultaneously, and the targets are set with random positions and phases regardless of the coherence. These targets are recovered with the CS-MFP sensor model and the traditional MFP model. The average relative square errors of 100 trials are shown in Figure 6. With greater number of targets, the CS-MFP model outperform traditional model overwhelmingly.

Figure 6. (a) Recovery errors w.r.t. numbers of targets and sensors with CS-MFP model; (b) Recovery errors w.r.t. numbers of targets and sensors with traditional MFP model.

In the next experiment, the multi-frequency case is simulated. In the experiment, three frequencies 4 kHz, 8 kHz, 16 kHz are chosen as the frequency grids. The numbers of normal modes are 16, 32 and 67 respectively. Figure 7a presents a picture of the Gram matrix of matrix A. Three blocks are clearly seen along the diagonal. As can be seen, the high coherent atoms are distributed mainly in atoms with the same frequency. Hence, expanding dictionary will not obviously worsen the coherence parameter in this case. CS is applicable nevertheless. Figure 7b,c presents a result. Three targets with different frequencies are generated in the area of interest, and they are recovered with CCOOMP.

Figure 7. (a) A Gram matrix of

A

in a multi-frequency case; (b) The targets of interest; (c) The recovery result.

Figure 7. (a) A Gram matrix of

A

in a multi-frequency case; (b) The targets of interest; (c) The recovery result.

In the last experiment, we examine the performance of the algorithm w.r.t. noise. In this experiment, two targets are placed randomly with random amplitudes and phases. The coherence of the targets is less than 0.5. The measurement noise is assumed as White Gaussian noise. We use bottleneck distance as the measurement, i.e., the recovery is counted as a success when the support of x is recovered. Figure 8 shows the successful recovery rate in terms of signal to noise ratio (SNR) in dB with 100 trials in each case. From Figure 8, with cases of SNR > 30 dB, good performance can be obtained. This constraint applies to the sea with low Beaufort scales and strong underwater targets.

Figure 8. The successful recovery rate in terms of SNR.

5. Conclusions

Based on acoustic propagation in medium with boundaries, this paper presented a CS-MFP model in the application of underwater sparse target imaging and localizing. By exploiting the sparsity of targets in underwater environment, the feasibility of CS method in MFP application is discussed. Both the coherence of the constructed matrix and the robustness of the coherence w.r.t. the displacement of the sensors were analyzed. However, the direct application of CS may result in a high coherence parameter problem. The proposed CCOOMP algorithm was effective in solving high coherent dictionary recovery problems, especially when an atom's high coherent atoms do not form a neighborhood around it. Finally, numerical experiments were provided to demonstrate the effectiveness of the proposed CS-based MFP method.

Acknowledgments

This work was partly supported by 111 Project of China under Grant B14010, partly supported by the Program for New Century Excellent Talents in University of Ministry of Education of China under Grant NCET-13-0049), partly supported by supported by the Chang Jiang Scholars Program under Grant T2012122, partly supported in part by China National Science Foundation under Grant 61271391 and China Ministry Research Foundation under Grant 9140A21050114HT05338. The authors thank Albert Fannjiang with the discussion of the BLOOMP algorithm.

Author Contributions

All of the authors contributed extensively to the preparation of this manuscript. Huichen Yan conceived the work, designed the algorithms and wrote the manuscript. Jia Xu designed the experiments, wrote the manuscript and commented on the manuscript. Teng Long provided the analysis in terms of the robustness and made suggestions. Xudong Zhang conceived the work and commented on the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

Baggeroer, A.B.; Kuperman, W.A.; Schmidt, H. Matched field processing: Source localization in correlated noise as an optimum parameter estimation problem. J. Acoust. Soc. Am. 1998, 83, 571–587. [Google Scholar] [CrossRef]
Tolstoy, A.; Diachok, O.; Frazer, L.N. Acoustic tomography via matched field processing. J. Acoust. Soc. Am. 1991, 89, 1119–1127. [Google Scholar] [CrossRef]
Thode, A.M.; D’Spain, G.L.; Kuperman, W.A. Matched field processing, geoacoustic inversion, and source signature recovery of blue whale vocalizations. J. Acoust. Soc. Am. 2000, 107, 1286–1300. [Google Scholar] [CrossRef] [PubMed]
Tolstoy, A. Matched Field Processing in Underwater Acoustics; World Scientific: Singapore, 1993; pp. 93–153. [Google Scholar]
Baggeroer, A.B.; Kuperman, W.A.; Mikhalevsky, P.N. An overview of matched field methods in ocean acoustics. IEEE J. Ocean Eng. 1993, 18, 401–424. [Google Scholar] [CrossRef]
Fannjiang, A.C.; Strohmer, T.; Yan, P. Compressed remote sensing of sparse objects. SIAM J. Imaging Sci. 2010, 3, 595–618. [Google Scholar]
Fink, M.; Prada, C. Acoustic time-reversal mirror. Inverse Probl. 2001, 17. [Google Scholar] [CrossRef]
Donoho, D.L. Compressed sensing. IEEE Trans. Inform. Theory 2006, 52, 1289–1306. [Google Scholar] [CrossRef]
Baraniuk, R.G. Compressive sensing. IEEE Signal Process. Mag. 2007, 24, 118–124. [Google Scholar] [CrossRef]
Candes, E.J.; Romberg, J.K.; Tao, T. Stable signal recovery from incomplete and inaccurate measurements. Commun. Pure Appl. Math. 2006, 59, 1207–1223. [Google Scholar]
Tropp, J.A.; Gilbert, A. Signal recovery from partial information via orthogonal matching pursuit. IEEE Trans. Inf. Theory 2007, 53, 4655–4666. [Google Scholar] [CrossRef]
Lustig, M.; Donoho, D.L.; Pauly, J.M. Sparse MRI: The application of compressed sensing for rapid MR imaging. Magn. Reson. Med. 2007, 58, 1182–1195. [Google Scholar] [CrossRef] [PubMed]
Herman, M.A.; Strohmer, T. High-resolution radar via compressed sensing. IEEE Trans. Signal Process. 2009, 57, 2275–2284. [Google Scholar] [CrossRef]
Yan, H.; Xu, J.; Peng, S.; Zhang, X. A compressed sensing method for a wider swath in synthetic aperture imaging. In Proceedings of the IET Interantional Radar Conference, Xi’an, China, 14–16 April 2013; pp. 1–6.
Yan, H.; Xu, J.; Xia, X.G.; Liu, F.; Peng, S.; Zhang, X.; Long, T. Wideband underwater sonar imaging via compressed sensing with scaling effect compensation. Sci. China Inf. Sci. 2015, 58, 1–11. [Google Scholar]
Mantzel, W.; Romberg, J.; Sabra, K. Compressive matched-field processing. J. Acoust. Soc. Am. 2012, 132, 90–102. [Google Scholar] [CrossRef] [PubMed]
Katsnelson, B.G.; Petnikov, V.G. Shallow Water Acoustics; Springer Science & Business Media: Berlin, Germany, 2002; pp. 165–180. [Google Scholar]
Herman, M.A.; Strohmer, T. General deviants: An analysis of perturbations in compressed sensing. IEEE J. Sel. Top. Signal Process. 2010, 4, 342–349. [Google Scholar] [CrossRef]
Chi, Y.; Scharf, L.L.; Pezeshki, A. Sensitivity to basis mismatch in compressed sensing. IEEE Trans. Signal Process. 2011, 59, 2182–2195. [Google Scholar]
Duarte, M.F.; Baraniuk, R.G. Spectral compressive sensing. Appl. Comput. Harmonic Anal. 2013, 35, 111–129. [Google Scholar] [CrossRef]
Fannjiang, A.C.; Liao, W. Coherence pattern-guided compressive sensing with unresolved grids. SIAM J. Imaging Sci. 2012, 5, 179–202. [Google Scholar] [CrossRef]
Tang, G.; Bhaskar, B.N.; Shah, P.; Recht, B. Compressive sensing off the grid. IEEE Trans. Inf. Theory 2012, 59, 7465–7490. [Google Scholar] [CrossRef]
Chen, Y.; Chi, Y. Robust spectral compressed sensing via structured matrix completion. IEEE Trans. Inf. Theory 2014, 60, 6576–6601. [Google Scholar]
Stephens, R.W.B. Underwater Acoustics; John Wiley & Sons: New York, NY, USA, 1970. [Google Scholar]
Tropp, J.A. Greed is good: Algorithmic results for sparse approximation. IEEE Trans. Inf. Theory 2004, 50, 2231–2242. [Google Scholar] [CrossRef]
Yu, Y.; Petropulu, A.P.; Poor, H.V. MIMO radar using compressive sampling. IEEE J. Sel. Top. Signal Process. 2010, 4, 146–163. [Google Scholar]

© 2015 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Yan, H.; Xu, J.; Long, T.; Zhang, X. Underwater Acoustic Matched Field Imaging Based on Compressed Sensing. Sensors 2015, 15, 25577-25591. https://doi.org/10.3390/s151025577

AMA Style

Yan H, Xu J, Long T, Zhang X. Underwater Acoustic Matched Field Imaging Based on Compressed Sensing. Sensors. 2015; 15(10):25577-25591. https://doi.org/10.3390/s151025577

Chicago/Turabian Style

Yan, Huichen, Jia Xu, Teng Long, and Xudong Zhang. 2015. "Underwater Acoustic Matched Field Imaging Based on Compressed Sensing" Sensors 15, no. 10: 25577-25591. https://doi.org/10.3390/s151025577

Article Menu

Underwater Acoustic Matched Field Imaging Based on Compressed Sensing

Abstract

1. Introduction

2. Compressed Sensing MFP Model Formulation

2.1. Wave Propagation Theories

2.2. CS-MFP Model Formulation

2.3. Robustness of Coherence

2.4. Multi-Frequency CS-MFP Case

3. CS-MFP Recovery with High Coherent Dictionaries via CCOOMP

4. Numerical Experiments and Performance Analysis

5. Conclusions

Acknowledgments

Author Contributions

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI