Underwater Acoustic Matched Field Imaging Based on Compressed Sensing

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.


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 Bandexcluded 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.

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.

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 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 0 0 r r , z z = = , denoted as ( 0 r , 0 z ), with unit strength and an angle frequency ω will generate acoustic pressure field at a far-field position ( r, z ) as where ⋅     is the ceiling operator and z k , ς are vectors of length N representing the eigenvalues of wave modes, It can be shown that ς and z k 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, z0) shown in Figure 1a as the point in the left. The reader can be referred to [24] for a full presentation of propagation theory. The pressure field intensity generated by a source.

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' A uniformly linear sensor array with P elements is set in the far field at positions S S ( , ) r z to collect data denoted as y . The sensors can receive all the existing sound pressure generated by these targets. If we set ( ) S S r,z ( , ) = r z and 0 0 T T

( r ,z ) ( , ) = r z in
Equation (2), then ( ) p r,z can be also represented as Hence, y is a linear superposition of the entire sound pressure intensities generated from all the existing candidate pixels according to Equation (2) as where A represents the discrete version of Equation (2) and e is an observation error. With slight abuse of notations we have where i = 1,2, ···, P, j = 1,2,···,Q 2 , vec(·) is a function piling a matrix to a vector column by column.
Clearly, Equation (6) can be viewed as the weighted sum of the atoms in A with the weights representing the strengths of the targets. 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 Ax 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., 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 as μ .
Here we use i, j instead of n to distinguish two subscriptions. According to the definition, the Gram 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 0 z but a ∆ r distance in r direction will generate two atoms with a correlation coefficient as 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 Now let us focus on the second part of Equation (11). Similar to the above discussion, we have for two matrix i A and j A ( i j ≠ ), two targets with the same r and a ∆ z distance in the z direction, For two targets with the same z and a ∆ r distance in the r direction, in Equation (12). Similarly, ( )=0 z i, j μ and ( )=0 r i, j μ 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 in Equation (12) According to [25], if the target number is less than then in Equation (9) can be recovered uniquely using greedy algorithms. This Theorem determines the performance of CS-MFP model.

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 And the denominator is

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 rz space, frequency dimension is added as the third dimension, whose dimensionality is denoted as R . In this way, we obtain a more intricate matrix  to replace A in Equation (9) with 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.

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 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 Figures 1b and 3b). We then modify the BLOOMP method to solve this problem. Assuming 0 1 < η < , we first define the η-coherence index set of atom k to be the set and the η-coherence index set of an index set  is hence defined as Here we assume the targets are not distributed on positions with correlation coefficient larger than η, or for any 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.  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. 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. 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. 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. 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.

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.