Maximum Likelihood Deconvolution of Beamformed Images with Signal-Dependent Speckle Fluctuations from Gaussian Random Fields: With Application to Ocean Acoustic Waveguide Remote Sensing (OAWRS)

Abstract

Wide area acoustic remote sensing often involves the use of coherent receiver arrays to determine the spatial distribution of sources and scatterers at any instant. The resulting acoustic intensity images are typically corrupted by signal-dependent noise from Gaussian random field fluctuations arising from the central limit theorem and have a spatial resolution that depends on the incident direction, sensing array aperture and wavelength. Here, we use the maximum likelihood method to deconvolve the intensity distribution measured on a coherent line array assuming a discrete angular distribution of incident plane waves. Instantaneous wide area population density images of fish aggregations measured with Ocean Acoustic Waveguide Remote Sensing (OAWRS) are deconvolved to illustrate the effectiveness of this approach in improving angular resolution over conventional planewave beamforming.

1. Introduction

Coherent receiver arrays are often used to generate acoustic images of the undersea environment to chart the spatial distribution of acoustic scatterers or sources, such as fish, marine mammals, plankton, submerged objects and bathymetric features [1,2,3]. In a typical wide-area application, such as Ocean Acoustic Waveguide Remote Sensing (OAWRS), involving a horizontal, coherent line array, an acoustic image is formed by charting the acoustic intensity received from scatterers or sources in range by travel time analysis and in azimuth by plane-wave beam forming [4,5,6,7]. Conventional time-harmonic plane wave beamforming [2,3] on a line array by Fourier transform from the spatial to the wavenumber domain can non-linearly blur the angular distribution of incident plane waves for array geometries that lack spherical symmetry due to the nonlinear relationship between the incident direction and wavenumber components along the array caused by foreshortening.

Here, we develop and apply a maximum likelihood (ML) approach for the spatial deconvolution of acoustic intensity images generated by a coherent acoustic array. We take advantage of the relatively pervasive circular complex Gaussian random (CCGR) statistics of acoustic fields in the ocean that follow from the central limit theorem for many typical scenarios of propagation, radiation or scattering [4,5,6,8,9,10,11,12,13,14]. Instantaneous intensity then follows a negative exponential distribution with a standard deviation proportional to the mean [9]. This leads to intensity images that have inherent signal-dependent noise known as speckle [8,9,11,15], as well as possible independent additive noise in the CCGR fields. Given these statistics, the ML estimator is then derived to resolve the angular distribution of incident plane waves from coherent array measurements. We focus on OAWRS applications [4,5,6,11,13,14] where side lobes are sufficiently low such that distinct beams are statistically independent of each other and far field signal and noise sources dominate. Many past statistical models for beamforming in the ocean [16,17,18,19,20,21] are less well suited to OAWRS applications, because they assume deterministic signals that are inconsistent with observations and additive noise with spatial correlations that are inconsistent with nonuniform far field origins. While the approach presented here is general, a line array geometry is used in examples both for simplicity and because line arrays are typically used in OAWRS applications. The examples show that the ML method leads to significant increases in angular resolution over conventional space-wavenumber Fourier transform-based planewave beamforming.

2. Methods: Theoretical Formulation

2.1. Beamformed Intensity Measurement on a Receiver Array

Let $P_i(\overline{r},t)=P({\overline{k}}_i,\omega )e^{j({\overline{k}}_i·\overline{r}-\omega t)}$ denote the pressure field at point $\overline{r}$ and time t due to a time-harmonic plane wave propagating in the direction of the acoustic wavenumber vector ${\overline{k}}_i=k_x^i{\widehat{i}}_x+k_y^i{\widehat{i}}_y+k_z^i{\widehat{i}}_z$, where $P({\overline{k}}_i,\omega )$ is the complex plane wave amplitude, $k=|{\overline{k}}_i|=\omega /c$ is the acoustic wavenumber, ω is the radial frequency, c is the speed of sound and i is an index for $i=1,2,3,...,M$ directions. A plane wave from any of these directions is incident on a discrete receiver array of length L containing N elements. The pressure field from any direction received on the n-th element of the array located at ${\overline{r}}_n=(0,y_n,0)$ (Figure 1) is $P_i({\overline{r}}_n,t)=P({\overline{k}}_i,\omega )e^{j(k_y^i{y}_n-\omega t)}$.

The y-component of the wavenumber vector can be expressed as a function of the vertical inclination angle ψ and horizontal azimuthal angle θ of the incident plane wave as $k_y^i=ksin{\psi }_{i}sin{\theta }_i$. Then, in the presence of M incident plane waves, assuming the inclination angle ${\psi }_i=\pi /2$, the total pressure field received on the n-th element of the receiver array is:

$$P({\overline{r}}_n,t)=\sum _{i=1}^{M}P({\overline{k}}_i,\omega )e^{j(k{y}_{n}sin{\theta }_i-\omega t)}$$

(1)

The total beamformed pressure field at scan angle θ in terms of plane wave incident angle ${\theta }_i$ is:

$$\begin{array}{ccc}\hfill P_B(\theta ,t,\omega )& =& \sum _{i=1}^M\sum _{n=-\frac{N-1}{2}}^{\frac{N-1}{2}}P({\overline{k}}_i,\omega )e^{j(k{y}_n[sin{\theta }_i-sin\theta ]-\omega t)}\hfill \\ \mathrm{(2)}& & =& \sum _{i=1}^{M}P({\overline{k}}_i,\omega )B(sin\theta -sin{\theta }_i,\omega )e^{-j\omega t}\hfill \end{array}$$

where $B(sin\theta ,\omega )$ is the discrete array beam pattern given [1,2,3] by $B(sin\theta ,\omega )=\frac{1}{N}\frac{sin\left(\frac{k}{2}N\Delta ysin\theta \right)}{sin\left(\frac{k}{2}\Delta ysin\theta \right)}$ for θ spanning $[-\pi ,2,\pi /2]$ and array element spacing $\Delta y=\frac{L}{N-1}$.

The beamformed pressure field in Equation (2) is the convolution of the complex pressure amplitude of incident plane waves from different directions and the receiver array beam pattern. In the case of a single incident plane wave from direction ${\theta }_0$, Equation (2) simplifies to $P_B(\theta ,t,\omega )=P({\overline{k}}_0,\omega )B(sin\theta -sin{\theta }_0,\omega )e^{-j\omega t}$, which is proportional to the beam pattern steered in the direction of the incident plane wave, as illustrated in Figure A1 for a single incident plane wave at 30${}^{\circ }$.

The beamformed plane wave intensity measured on the array is then:

$$W(\theta ,\omega )={\left|P_B(\theta ,t,\omega )\right|}^2=|\sum _{i=1}^{M}P({\overline{k}}_i,\omega )B(sin\theta -sin{\theta }_i,\omega ){|}^2$$

(3)

By the central limit theorem, the incident plane wave complex amplitude $P({\overline{k}}_i,\omega )$ follows circular complex Gaussian random (CCGR) statistics with zero mean due to random fluctuations arising from waveguide propagation, scattering, source mechanisms or a combination of these [4,5,6,11,12,13,14], such that $\langle P({\overline{k}}_i,\omega )\rangle =0$. Moreover, incident plane wave fields from different directions are assumed to be independent of each other as found in typical OAWRS applications [4,5,6,13,14], such that $\langle P({\overline{k}}_i,\omega )P^{*}({\overline{k}}_j,\omega )\rangle =\langle P({\overline{k}}_i,\omega )\rangle \langle P^{*}({\overline{k}}_j,\omega )\rangle$ for $i\ne j$. The expected beamformed intensity is then given by:

$$\sigma (\theta ,\omega )=\langle W(\theta ,\omega )\rangle =\sum _{i=1}^{M}S({\overline{k}}_i,\omega ){\left|B(sin\theta -sin{\theta }_i,\omega )\right|}^2$$

(4)

where $S({\overline{k}}_i,\omega )=\langle |P({\overline{k}}_i,\omega ){|}^2\rangle$ is the expected intensity of the incident plane wave propagating in the direction ${\overline{k}}_i$.

2.2. Deconvolution of Beamformed Intensity from a Line Array: Maximum Likelihood Estimate of Expected Source Intensity

Beamformed intensity measurements (Equation (3)) at $j=1,2,3,...,J$ discrete scan angles can be expressed as:

$$W_j=W({\theta }_j,\omega )=|\sum _{i=1}^{M}P({\overline{k}}_i,\omega )B(sin{\theta }_j-sin{\theta }_i,\omega ){|}^2$$

(5)

which, in matrix form, can be written as:

$$W_j=|\sum _{i=1}^M{P}_i{B}_{ij}{|}^2$$

(6)

where $P_i=P({\overline{k}}_i,\omega )$ are the elements of a vector $\mathbf{P}$ and $B_{ij}=B(sin{\theta }_j-sin{\theta }_i,\omega )$ are the elements of the scanning matrix $\mathit{B}$. The expected beamformed intensity from Equation (4) then becomes:

$${\mathit{\sigma }}^T={\mathbf{S}}^T\mathbb{B}$$

(7)

where the vector $\mathit{\sigma }$ contains expected beamformed intensities $\langle W_j\rangle$, the vector $\mathbf{S}$ contains the expected plane wave intensities $S_i=S({\overline{k}}_i,\omega )=\langle P({\overline{k}}_i,\omega )P^{*}({\overline{k}}_i,\omega )\rangle$ from all horizontal azimuths and the elements of matrix $\mathbb{B}$ are ${\mathbb{B}}_{ij}={\left|B_{ij}\right|}^2$.

The expected values of measurements $\langle W_j\rangle ={\sigma }_j(\mathbf{S})$ depend on the expected incident plane wave intensity vector $\mathbf{S}$, which in general varies in space, and is to be estimated from $\mathbf{W}$. In typical operational scenarios for OAWRS imaging, the combination of a sufficiently large, densely-sampled and spatially-tapered aperture [4,5,6,7] leads to low enough side lobes that non-overlapping beams are effectively independent for far field signal and noise sources. The conditional probability distribution for all measurements in vector $\mathbf{W}$ given $\mathbf{S}$ is then the product of gamma distributions [9,15]:

$$P(\mathbf{W}|\mathbf{S})=\prod _{j=1}^J\frac{{\left[\frac{{\mu }_j}{{\sigma }_j(\mathbf{S})}\right]}^{{\mu }_j}{W}_j^{{\mu }_j-1}exp\left[-{\mu }_j\frac{W_j}{{\sigma }_j(\mathbf{S})}\right]}{\Gamma ({\mu }_j)}$$

(8)

where ${\mu }_j$ is the number of coherence cells in the intensity average derived from a sum of independent measurements, which is also equal to the signal-to-noise ratio ${\langle W_j\rangle }^2/(\langle W_j^2\rangle -{\langle W_j\rangle }^2)$ [9,15].

The log-transformed vector $\mathbf{L}$ defined by $L_j=ln(W_j/I_{ref})$, where $I_{ref}$ is the reference intensity, obeys the conditional distribution [15]:

$$P(\mathbf{L}|\mathbf{S})=\prod _{j=1}^J\frac{{\left[\frac{{\mu }_j}{{\sigma }_j^{\prime }(\mathbf{S})}\right]}^{{\mu }_j}exp\left[-{\mu }_j\frac{exp(L_j)}{{\sigma }_j^{\prime }(\mathbf{S})}+{\mu }_j{L}_j\right]}{\Gamma ({\mu }_j)}$$

(9)

which is also the likelihood function, where ${\sigma }_j^{\prime }(\mathbf{S})={\sigma }_j(\mathbf{S})/I_{ref}$. The log-likelihood function is then:

$$lnP(\mathbf{L}|\mathbf{S})=\sum _{j=1}^J{\mu }_{j}ln\left[\frac{{\mu }_j}{{\sigma }_j^{\prime }(\mathbf{S})}\right]+\left[\frac{-{\mu }_{j}exp(L_j)}{{\sigma }_j^{\prime }(\mathbf{S})}+{\mu }_j{L}_j\right]-ln\Gamma ({\mu }_j)$$

(10)

The value of $\mathbf{S}$ for which the log-likelihood function attains its global maxima is the maximum likelihood estimate (MLE) of $\mathbf{S}$,

$$\widehat{\mathbf{S}}=arg\underset{\mathbf{S}}{max}\left[lnP(\mathbf{L}|\mathbf{S})\right]$$

(11)

which is also the value of $\widehat{\mathbf{S}}$ for which the derivative of $lnP(\mathbf{L}|\mathbf{S})$ with respect to $\mathbf{S}$ vanishes:

$$\frac{dlnP(\mathbf{L}|\mathbf{S})}{d\mathbf{S}}=\sum _{j=1}^J-\frac{{\mu }_j}{{\sigma }_j^{\prime }(\mathbf{S})}\frac{d{\sigma }_j^{\prime }(\mathbf{S})}{d\mathbf{S}}+\frac{{\mu }_j}{{\sigma }_j^{\prime }{(\mathbf{S})}^2}\frac{W_j}{I_{ref}}\frac{d{\sigma }_j^{\prime }(\mathbf{S})}{d\mathbf{S}}=0{|}_{\mathbf{S}=\widehat{\mathbf{S}}}$$

(12)

The MLE of $\mathbf{S}$ is then:

$$\widehat{\mathbf{S}}={\left[\mathbb{B}{\widehat{\mathit{d}}}_{\sigma }{\mathbb{B}}^T\right]}^{-1}\mathbb{B}{\widehat{\mathit{d}}}_{\sigma }\mathbf{W}$$

(13)

where the elements of diagonal matrix ${\widehat{\mathit{d}}}_{\sigma }$ are ${\left[{\widehat{d}}_{\sigma }\right]}_{ji}={\delta }_{ji}{\mu }_j/{W_j^{\prime }}^2$ and $W_j^{\prime }=W_j/I_{ref}$, as shown in Appendix A.2.

When the number of scan angles J is equal to the number of incident plane wave directions M, then matrices $\mathbb{B}$ and ${\widehat{\mathit{d}}}_{\sigma }$ are square matrices, and Equation (13) simplifies to:

$$\widehat{\mathbf{S}}={{\mathbb{B}}^T}^{-1}\mathbf{W}$$

(14)

and the solution of $\widehat{\mathbf{S}}$ is linear in $\mathbf{W}$.

2.3. Statistics of the Maximum Likelihood Estimator

A criterion to determine the performance of an estimator $\widehat{\mathbf{S}}$ obtained from measurement vectors $\mathbf{W}$ or $\mathbf{L}$ is the Cramer–Rao Lower Bound (CRLB) on the parameter’s variance [22,23] given by:

$$\text{Var}({\widehat{S}}_i)\ge {\left[{\mathbb{I}}^{-1}((\mathbf{S}))\right]}_{ii}$$

(15)

where $S_i$ is the true parameter value and $\mathbb{I}(\mathbf{S})$ is the Fisher information matrix with elements:

$${\left[\mathbb{I}\right]}_{in}(\mathbf{S})=-〈\frac{{\partial }^2}{\partial S_i\partial S_n}lnP(\mathbf{L}|\mathbf{S})〉$$

(16)

The Fisher information matrices for $\mathbf{W}$ and $\mathbf{L}$ are identical and equal to [11,15]:

$${\left[\mathbb{I}\right]}_{in}(\mathbf{S})=\sum _{j=1}^J\frac{{\mu }_j}{{\sigma }_j^{\prime }{(\mathbf{S})}^2}\frac{\partial {\sigma }_j^{\prime }(\mathbf{S})}{\partial S_i}\frac{\partial {\sigma }_j^{\prime }(\mathbf{S})}{\partial S_n}$$

(17)

which can be further simplified using Equation (7) to:

$$\begin{array}{ccc}\hfill {\left[\mathbb{I}\right]}_{in}(\mathbf{S})& =& \sum _{j=1}^N\frac{{\mu }_j}{{\sigma }_j{(\mathbf{S})}^2}{\mathbb{B}}_{ij}{\mathbb{B}}_{nj}\hfill \\ \hfill \mathbb{I}& =& \mathbb{B}\Xi {\mathbb{B}}^T\hfill \end{array}$$

(18)

where Ξ is a diagonal matrix with elements ${\Xi }_{jj}=\frac{{\mu }_j}{{\sigma }_j{(\mathbf{S})}^2}$.

In the case of $M=J$, $\text{Var}({\widehat{S}}_i)$ can be obtained from Equation (14) as:

$$\begin{array}{ccc}\hfill \text{Var}({\widehat{S}}_i)& =& \text{Var}({{\mathbb{B}}^T}^{-1}\mathbf{W})\hfill \\ & =& \sum _{j=1}^M{\mathfrak{B}}_{ij}^2\frac{{\sigma }_j{(\mathbf{S})}^2}{{\mu }_j}\hfill \\ \mathrm{(19)}& & =& \text{Diagonal}(\mathfrak{B}{\Xi }^{-1}{\mathfrak{B}}^T)\hfill \end{array}$$

where the matrix $\mathfrak{B}={{\mathbb{B}}^T}^{-1}$ and $\text{Var}(W_j)=\frac{{\sigma }_j{(\mathbf{S})}^2}{{\mu }_j}$ from the Gamma distribution of $W_j$ in Equation (8). The variance of the MLE when $J=M$ then attains the CRLB, since it equals the diagonal elements of the inverse of $\mathbb{I}$ from Equation (18):

$$\begin{array}{ccc}\hfill {\mathbb{I}}^{-1}& =& {{\mathbb{B}}^T}^{-1}{\Xi }^{-1}{\mathbb{B}}^{-1}\hfill \\ \mathrm{(20)}& & =& \mathfrak{B}{\Xi }^{-1}{\mathfrak{B}}^T\hfill \end{array}$$

For sufficiently large μ, the ML estimator’s probability distribution $\widehat{\mathbf{S}}$ is asymptotically Gaussian given by [22,23]:

$$\widehat{\mathbf{S}}{\sim }^a\mathcal{N}(\mathbf{S},J^{-1}(\mathbf{S}))$$

(21)

where the variance is asymptotically equal to the CRLB given in Equation (15), and the mean converges to the true parameter value making the ML estimator asymptotically optimal since it becomes unbiased and attains minimum variance [22,23].

3. Results: Illustrative Examples

Here, we use the maximum likelihood method to estimate the angular intensity distribution of incident plane waves from coherent line array measurements with analytic, synthetic and actual wide-area OAWRS imagery data [4,5,6]. In all of the examples shown in this section, $N=64$ array elements spaced $0.75$ m apart are used at a frequency of 950 Hz, as is consistent with some typical OAWRS imaging scenarios [4,5]. The array element spacing corresponds to half the acoustic wavelength for a cutoff frequency of 1000 Hz and speed of sound in water c = 1500 m/s.

3.1. Analytic Solutions

For a single incident plane wave with no signal dependent noise, the ML deconvolution of the conventional beamformed measurement results in the ground truth incident plane wave by invariance of the ML estimator [22,23], as illustrated for both broadside and end-fire incidence examples in Figure 2. This can be seen by noting that in the absence of randomness, the beamformed measurement is its expected value $\mathbf{W}=\mathit{\sigma }$. The ML estimate can be then found analytically either via Equation (14) or via Equation (7), where resulting ML estimate $\widehat{\mathbf{S}}$ becomes $\mathbf{S}$, which is unity in the direction of the plane wave and zero elsewhere. Reversing the process, the convolution of this pressure field with the receiver array beam pattern yields the conventional beamformed data, which is the beam pattern itself shifted to the direction of the incident pressure field.

3.2. Synthetic Data

We now illustrate ML deconvolution of angular or spatial features sensed by nonuniform distributions of incident CCGR fields with synthetic data and numerical simulation (Figure 3, Figure 4 and Figure 5).

For prominent and relatively narrow angular box-like features, a roughly 50% improvement in angular resolution is found (Figure 3 and Figure 4). For a box-like feature of roughly 1${}^{\circ }$ width at broadside, the 3-dB beamwidth of the deconvolved plane wave intensity distribution is found to be roughly 1.2${}^{\circ }$ (Figure 3). This is roughly two times smaller than that from conventional beamforming of roughly 3${}^{\circ }$. Similarly, near array end-fire for a box-like feature of roughly an 8${}^{\circ }$ width, the 3-dB beamwidth of the deconvolved plane wave intensity distribution is roughly 9${}^{\circ }$ (Figure 4), which is roughly two-times smaller than the 20${}^{\circ }$ width found in conventional beamforming. Only at much lower values at least 5–7 dB down from the peak intensity does the match between ground truth and ML deconvolved estimate deteriorate. This deterioration, however, is barely perceivable in a linear intensity scale. The expected intensity of the ML deconvolved estimate matches the ground truth well, within roughly ±1 dB for the broadside case and ±0.5 dB for the end-fire case. This comprises roughly an order of magnitude improvement over conventional beamforming, as shown in Figure 3 and Figure 4. In Figure 3 and Figure 4, a spatial Hanning window taper is applied; beamformed intensity is averaged over five independent and instantaneous realizations, such that $\mu =5$ for each $W({\theta }_j,\omega )$, as is consistent with typical OAWRS imaging [4,5]; simulations are performed as described in Appendix A.4; and the beamformed data $W({\theta }_j,\omega )$ is deconvolved using the algorithm in Appendix A.5.

For a prominent and relatively wide two-dimensional box-like feature that extends both in azimuth and range from the receiver array center, a roughly 40% improvement in angular resolution over conventional beamforming is found (Figure 5). Such features are representative of actual fish distributions in the ocean [4,5,6]. The improvement in angular resolution is measured in terms of the spatial span or the area within the 3-dB down contour from the peak intensity of the feature after conventional beamforming and subsequent ML deconvolution with respect to the true areal coverage of the feature. For a feature centered at 50${}^{\circ }$ from an array broadside of roughly a 9${}^{\circ }$ width and a 1-km extent in range, the error in the spatial span of the ML deconvolved feature is roughly 8% of the true feature area. In comparison, the error in the spatial span of the feature from conventional beamforming is several times higher at roughly 50% of the true feature area. In Figure 5, a spatial Hanning window taper is applied; beamformed intensity is averaged over 10 independent and instantaneous realizations, such that $\mu =5$ for each $W({\theta }_j,\omega )$, as is consistent with typical OAWRS imaging [4,5]; simulations are performed as described in Appendix A.4; and the beamformed data $W({\theta }_j,\omega )$ are deconvolved using the algorithm in Appendix A.5.

3.3. Ocean Acoustic Waveguide Remote Sensing Images of Fish Shoals

In the fall of 2006, vast shoals of spawning Atlantic herring were imaged using OAWRS in the Gulf of Maine (Figure 6) [5,6]. During evening hours, small clusters or schools of fish aggregated together to form large shoals spanning several tens of kilometers in length.

Maximum likelihood deconvolution of OAWRS fish shoal imagery leads to a roughly 10%–50% improvement in angular resolution over conventional beamforming (Figure 7, Figure 8, Figure 9 and Figure 10). For features located 20–50${}^{\circ }$ from the array end-fire with population density >0.5 fish/m${}^2$ (Figure 7, Figure 8, Figure 9 and Figure 10), a 40%–50% improvement in resolution is found. Similarly, for features located close to the array broadside, a 10%–25% improvement in resolution is found. The results indicate that ML improvements in angular resolution over conventional beamforming depend strongly on array beamwidth and the actual angular extent of the features. These results are consistent with synthetic data deconvolution results shown in Figure 3, Figure 4 and Figure 5.

In the OAWRS images shown in this section, $N=64$ receiver array elements are spaced 0.75 m apart; the sensing frequency is 950 Hz; and a spatial Hanning window taper is applied to the receiver array elements [4,5]. The array element spacing corresponds to half the acoustic wavelength for a cutoff frequency of 1000 Hz and c = 1500 m/s. The OAWRS beamformed measurements are deconvolved using the algorithm in Appendix A.5.

4. Discussion

Following experimental observations from OAWRS measurements, we assume the demodulated field of incident plane waves from a discrete set of directions follows circular complex Gaussian random statistics (CCGR) arising from random wave propagation, scattering or radiation via the central limit theorem [4,5,6,11,12,13,14]. Intensity measurements then have signal-dependent noise. By time-harmonic plane-wave beamforming [2,3,22], the received intensity at the array center from different scan directions is determined and obeys a Gamma distribution [15,24]. Such standard plane-wave beamforming by Fourier transform from the spatial to the wavenumber domain, however, can non-linearly blur the angular distribution of incident plane waves for array geometries that lack spherical symmetry due to the nonlinear relationship between the incident direction and wavenumber components along the array caused by foreshortening.

By employing a sufficiently large aperture and applying a spatial taper function across the array as in OAWRS applications, distinct directional beams are effectively statistically independent [4,5,6,7] for the far field signals and noise sources typically encountered in OAWRS. Maximizing the likelihood function for these independent beams then yields the maximum likelihood estimate of intensity incident from a given parameterized angular sector, which acts as a deconvolution process to extract the angular distribution of the incident planewave field from blurring caused by the array’s diffraction limited beam pattern.

The far field signal and noise sources encountered in OAWRS typically have correlation along the array [4,5,6,13,14] that is not known a priori. Past statistical models for acoustic beamforming, including past maximum likelihood approaches, in the ocean are not well suited to OAWRS applications, because they have overwhelmingly assumed a deterministic signal model in additive Gaussian noise with a priori knowledge of the noise’s spatial correlation, which is most often assumed to be spatially uncorrelated for horizontal apertures’ element spacing exceeding half the wavelength [16,17,18,19,20,21]. Neither of these assumptions is consistent with those observed in OAWRS sensing, where both signal and noise fields are random and highly correlated in space across the array due to their far field origins. Additionally, it has been noted [25] that some persistent confusion and misidentification has occurred between maximum likelihood [22,23] and other methods [19,20,21] in beamforming applications.

5. Conclusions

A maximum likelihood method is developed to deconvolve underwater acoustic intensity images corrupted with signal-dependent noise. The images are assumed to be generated by planewave beamforming of data received from a coherent spatial receiving array. The method is shown to enable significant improvements in resolving the directions of incident plane waves over conventional planewave beamforming by theoretical analysis, numerical simulation and the application to experimental data. Deconvolution of experimentally-measured Ocean Acoustic Waveguide Remote Sensing (OAWRS) images of fish areal population density show an angular resolution improvement of 40%–50% for fish clusters close to the array end-fire and 10%–25% for fish clusters close to the array broadside. The approach developed is generally applicable to planewave beamforming on multidimensional arrays of arbitrary geometry. Line array implementations are presented here for simplicity and to provide examples relevant to typical OAWRS scenarios.

multiple