Superdirective beamforming applied to SWellEx96 horizontal arrays data for source localization

The case where array aperture is short, on the order of several wavelengths, is of interest in many practical applications. For example, for hearing aids, where strong constraints are present on the maximum aperture of the array. A short microphone array, which can be built into an eyeglass frame, is attractive if it can improve speech intelligibility in noise.¹ “Advanced deployable arrays” use short horizontal arrays, say of 16 elements, for target localization and detection, as compared with traditional design of tens to hundreds of elements. Small arrays are limited in their performance using conventional beamforming (CBF) and require superdirectivity.

Superdirectivity is not precisely defined, but it generally refers to the fact that a small size array, using superdirective beamforming, can perform as well as a large size array using CBF. Superdirective array processing has received a lot of interest recently as it has the potential of providing not only high angular discrimination, but also high directivity index (DI) or directivity factor (DF), defined as the array gain (AG) in isotropic noise. In previous work, superdirective beamforming exploits the spatial dependence of the acoustic noise field to improve the angular resolution and DF.² Theoretically, DF can approach M² for a spherical array using M number of modes,^3,4 and DF for a linear array can also reach N² along the endfire direction in spherically isotropic noise for an array of N closely spaced elements.⁵ This phenomenon has been called supergain and is most noticed at frequencies lower than the designed frequency, defined as the frequency whose wavelength is twice the element spacing for a uniform array. However, the nice features derived theoretically for superdirective arrays are difficult to achieve in practice because superdirective beamforming is highly sensitive to signal mismatch due to imbalance in sensors characteristics, errors in array configuration, signal distortion due to scattering or imprecision of the signal model, etc. How to minimize the sensitivities to these errors is not a trivial task and many methods have been proposed in this regard. These methods impose constraints to reduce the sensitivity to the signal mismatch at the expense of the array gain.² 

The approach used here to achieve superdirective beamforming is by deconvolution.^6,7 Given that the beam power output of conventional (delay and sum) beamforming can be expressed as the convolution of the source distribution and the response of the beamformer to a point source, the latter referred to as the beam pattern or the point spread function (PSF), the original source distribution can be recovered by deconvolution, yielding not only improved accuracy in the direction of arrival (DOA) estimation and suppressed side lobes, but also high array gain compared with CBF (the deconvolution gain) as demonstrated by Yang for a uniformly spaced horizontal line array (HLA)⁶ and circular arrays.⁷ Deconvolution can be an ill-posed problem depending on the number of receivers, the nature of the signals, including the signal-to-noise ratio (SNR), etc. But, deconvolving the CBF beam power [deconvolved CBF (dCv) in short] has been shown to be well-behaved and produce a unique solution since the beam power (including the beam pattern) is positive.⁶ As a result, unique pencil-like beams were obtained for point sources. Another advantage of deconvolved CBF is that it retains the advantages of CBF in terms of the sensitivity and snapshots requirement. The method is insensitive to signal mismatch compared with other high resolution methods such as minimum variance distortionless response (MVDR) or superdirective mode decomposition methods reported previously. The reader is referred to Refs. 6, 7 for details.

This paper applies deconvolved CBF to the SWellEx96 horizontal array data and demonstrates that short horizontal arrays (aperture ∼ 4–8 wavelength) can achieve performance (superdirectivity) equivalent to that of an array of much (order of magnitude) larger aperture using CBF. It is a first underwater experimental demonstration of superdirectivity to the best of the author's knowledge.

The SWellEx96 horizontal array data have been analyzed by Tollesfen et al. for source localization using matched field processing.⁸ Source localization in latitude and longitude (or x-y coordinate) is coupled to and hence influenced by the estimated source depth. Multiple (5–9) frequencies, or multi-tonal signals, are needed to yield robust range and depth estimates. In this paper, only a single frequency tonal signal is used to localize the source in x-y.

It is well known that the CBF beam power is the convolution of the beam pattern with the source (including noise) distribution,⁶ as is the case for a linear system,

where the beam pattern $Bp(θ|ϑ)$ is the beam power (response) at angle $θ$ to a far-field point source at direction $ϑ$⁠. Deconvolution is a method to solve for $S(ϑ)$ given the data $BCBF(θ)$⁠. The solution will be denoted as $BdCv(θ)$⁠. Equation (1) can be written in a general form as $r(y)=∫−∞∞P(y|x)s(x) dx$⁠, where $P(y|x)$ is the PSF. The Richardson-Lucy (R-L) algorithm^9,10 is a Bayesian-based iterative method to solve for a normalized source term $s(x)$⁠, satisfying $∫−∞∞s(x) dx=1$⁠, given the measured data $r(y)$⁠, and the PSF $P(y|x)$⁠, which satisfies $∫−∞∞P(y|x) dy=1$⁠. It follows consequently that $∫−∞∞r(y) dy=1$⁠. Interpreting $P(y|x)dy$ as the probability that y lies in the interval $(y,y+dy)$⁠, conditioned on a given x, the R-L solution is given by $s(n+1)(x)=∫−∞∞r(y) Q(n)(x|y) dy$⁠, where $Q(n)(x|y)$ is the probability that $r(y)$ lies in $(y,y+dy)$ and is given, based on the Bayes theorem, by $Q(n)(x|y)=P(y|x)s(n)(x)/r(n)(y)$⁠, where $r(n)(y)=∫−∞∞P(y|x)s(n)(x) dx$⁠. The above algorithm was applied to a uniformly distributed line array and circular array, for which the beam pattern is shift invariant, i.e., $P(x|y)=P(x−y)$⁠. The readers are referred to Refs. 6 and 7 for details, and corresponding references. For an array of arbitrary configuration, this condition is not satisfied. For a shift variant PSF, where $P(x|y)≠P(x−y)$⁠, the integral $M(y)=∫−∞∞P(y|x) dx$ is a function of y and not a constant, thus $M(y)≠1$⁠. For this case, let us define $P¯(y|x)=P(y|x)/M(y)$ and $r¯(y)=r(y)/M(y)$⁠, one finds $r¯(y)=∫−∞∞P¯(y|x)s(x) dx$⁠, which satisfies the condition set by the R-L algorithm, namely, $∫−∞∞s(x) dx=1$ and $∫−∞∞P¯(y|x) dx=1$⁠. The solution is then given, after some manipulation, by

The HLA's in the SWellEx96 experiments are curved line arrays. The extended R-L algorithm as described above will be used for deconvolution.

To simulate superdirectivity and sugergain (assuming high input SNR, e.g., >30 dB), we shall consider a sub-array of the HLA North array deployed during the SWellEx96 experiment, consisting of 14 elements with an aperture of 107 m as shown in the upper left corner in Fig. 1(a). (See Ref. 11 for details.) Figure 1(b) shows the −3 dB beam width for this array using CBF and dCv as a function of frequency covering 40 to 160 Hz (array aperture L increasing from 2.9 λ to 11.4 λ, where λ is the wavelength) for a source along the track (shown by the solid line). One notes that the beam width using CBF increases rapidly as frequency decreases; the beam width is 4.5° at $L=11.4 λ$⁠, and reaches 16° at $L=2.9 λ$⁠. On the other hand, the beam width using dCv remains relatively small, 1°–2°, despites the decreasing frequency or decreasing L/λ. This is superdirectivity. Likewise, for a L = 3λ, the maximum array gain using CBF for isotropic noise is 10 log(7) ≈ 8.5 dB (for seven receivers spaced at half-wavelength). Using dCv, one obtains an additional deconvolution gain (DG) of ∼15 dB as shown in Fig. 1(c) (on top of the CBF AG). DG is the ratio of the directivity index between dCv and CBF, and is given by $DG=DIdCv/DICBF=∫−ππBCBF(θ|ϑ)dϑ/∫−ππBdCv(θ|ϑ)dϑ$⁠. DG has been studied for a HLA and circular array below design frequencies. It is the “supergain” obtained by deconvolution.

In this section, we apply CBF and dCv to the SWellex96 HLA data¹¹ to localize the moving source. We select a period when the source travels in between the HLA South array and vertical line array (VLA) similar to Ref. 8, where one notes that the source has a high bearing rate with respect to HLA North (∼1.2°/s), for which MVDR works poorly due to the snapshot deficient problem. [Signals arriving close to the endfire of a line array may appear in different beams (the split beams problem) due to multipaths arriving at different elevation angles and will be treated in a separate publication.] For each array, we first calculate the normalized CBF beam power given by $BCBF(θ)=|aH(θ)p|2/(|a(θ)|2|p|2)$⁠, where $p$ is a vector of the pressure field received on the array, and $a(θ)=exp(−ikr⋅u(θ))$ is the steering vector, where k is the wavenumber, $r$ is a vector of the two-dimensional coordinates of the array elements (or a displacement vector of element positions relative to that of a reference element), and $u$ is a unit vector pointing to the angle $θ$⁠. Next, we define an ambiguity function $Amb(x,y)=BCBF(θ)$⁠, where the beam power at position (x, y) is determined by the beam power at angle $θ$⁠, where $tan θ=(y−y0)/(x−x0)$⁠, $(x0,y0)$ denoting the position of the first element of the array. Finally, the source position (x,y) ambiguity function is obtained by summing the ambiguity functions from (each of) the north and south (sub-) arrays, i.e., $AmbT(x,y)=AmbN(x,y)+AmbS(x,y)$⁠. The peak of the ambiguity function can be used to determine the source position. The result is shown in Fig. 2(a) for CBF for a snapshot of 1.2 s 127 Hz signal. The ambiguity function in Fig. 2(a) is not symmetric since the source is close to HLA North and located on the broadside direction of the array, where the beam width is narrower; note that the cross-range resolution is proportional to range times the beam width. To compare with CBF, we apply deconvolution to the CBF beam power and create a similar normalized spatial distribution for each array (based on the dCv beam power), and sum the ambiguity functions from both arrays to create the dCv source position ambiguity function as shown in Fig. 2(b). To compare the resolution in source localization, we zoom in on the peak of the ambiguity function as shown by the insets in Figs. 2(a) and 2(b). The insert in Fig. 2(a) shows an area of size 500 m × 500 m and the insert in Fig. 2(b) shows an area of size 100 m × 100 m. One finds the resolution (the area at −3 dB contour) for dCv is order of magnitude smaller than that of CBF. The resolution varies from ping to ping. The time variations of the source localization results (the −3 dB contour of the source ambiguity functions) are shown in Fig. 2(c) using CBF and Fig. 2(d) using dCv (d) for the data analyzed, at a time increment of 12.5 s for the 127 Hz signal. Statistically speaking, the localization uncertainty as measured by the area of the −3 dB contours is order of magnitude (30–100 times) smaller for dCv than CBF, which is consistent with the beam resolution shown in Fig. 2(b).

The localization results at 64 Hz are shown in Fig. 2(e) using CBF and Fig. 2(f) using dCv (f). One observes a good resolution at 64 Hz for an array of aperture L = 4.6λ. To achieve the same beam width using CBF would require an array of L ≥ 46λ according to a figure for beam width for 64 Hz (not shown here) similar to that shown in Fig. 2(b). [Note that the beam width result, Fig. 2(b), scales as a function of aperture over wavelength.] It demonstrates “superdirecivity.”

Figure 3 shows an example of the signal beam power distribution measured on the HLA south sub-array, for a snapshot of ∼5 s at ∼ 1000 s into the recording of HLA data, using CBF and dCv for the 127 Hz signal [Fig. 3(a)] and 64 Hz signal [Fig. 3(c)]. The dCv shows a much narrower beam width and smaller sidelobe levels. Figures 3(b) and 3(d) show the beam spectra using CBF and dCv along the source look direction at 127 and 64 Hz, respectively. Also shown is the phone spectrum; all spectra are normalized by the peak power. One observes from Fig. 3(b) that the signal at 127 Hz has a SNR at ∼20 dB for individual hydrophones (based on an eyeball average of the noise level). The CBF beam output shows a SNR at ∼30 dB, yielding a CBF array gain at about 10 dB; a total of 30 dB AG is consistent with the Peak-to-sidelobe-level ratio (PSR) seen in Fig. 3(a), assuming the sidelobe level is noise. The beam output after deconvolution shows a much higher SNR than CBF. The exact number is dependent on where the noise level is measured. For the noise in the neighborhood of the signal, the SNR is 60 dB, which is consistent with the PSR of the dCv beam output shown in Fig. 3(a). If the noise level is measured by averaging the non-signals in Fig. 3(b), one obtains a SNR > 40 dB and hence a DG > 10 dB. Exact measurement of the dCv output SNR from the beam output spectrum (and DG) requires a statistical average of many snapshots of data (from a stationary source), since the noise background is fluctuating (rapidly) from snapshots to snapshots. The moving source data are non-stationary; the measurement of DG requires additional analysis and is beyond the scope of this paper. One further notes that the 127 and 130 Hz signals show equal strengths at the phone level and likewise in the CBF output power. The 127 Hz is observed in the dCv output, but the 130 Hz signal is not in Fig. 3(b). The reason is that the 127 Hz shows up at a bearing of 85.5° whereas the 130 Hz shows up at a bearing of 88.5° based on the dCv output. They fall in the same beam using CBF but fall in two separate beams using dCv. [Figure 3(b) is focused at 85.5° bearing.] Note that the 127 Hz signal is transmitted from a shallow source at 9 m depth, and 130 Hz signal is transmitted from a deep source at 54 m depth. They arrive at the HLA with different vertical angles which influence the bearing estimations.

In Fig. 3(d), one finds the 64 Hz signal has a SNR of ∼18–20 dB at the phone level, ∼25–30 dB at the CBF output, and a ∼35 dB SNR at the dCv output, leading to an estimated DG of ∼10 dB. Again, this is a rough estimate and more quantitative analysis is needed to estimate the DG. The at-sea data, Fig. 3, demonstrates the potential of “supergain” using dCv.

Classical delay-and-sum beamforming provides a wide beam and small amount of array gain at decreasing frequencies given a limited aperture, since the array performance is mainly determined by the ratio of aperture over wavelength, which decreases with decreasing frequencies. While adaptive array processing, such as MVDR, has been shown to yield narrow beam width, it usually does not work well at below design frequencies, where it was not intended for, and also not for a fast bearing-changing target due to the snapshot deficient problem. Superdirective arrays offer a heretofore overlooked solution in which optimal performance can be obtained for a stationary random noise field. It provides a directivity (beam width) with less aperture than that required by conventional arrays with the same directivity index (array gain). In other words, it reduces the size of the receiving array without sacrificing array performance. The above features are studied and illustrated using the SWellEx96 HLA data using dCv. It is noted that many proposed superdirective methods are sensitive to signal mismatch, such as that induced by sensor position errors, limiting their use in practice.² Deconvolution is relatively insensitive to such problems as previously studied in Refs. 6 and 7. Also, the beam width varies slowly with frequency and provides a frequency-invariant beam pattern for broadband signal processing. These are two useful features associated (uniquely) with deconvolution.

This work is supported by a grant from the National Science Foundation of China No. 61531017.