We assume that $K$ uncorrelated far field narrowband sources impinge on a linear array of $M$ equally spaced sensors with sensor spacing $d$, which is less than or equal to half the wavelength. The $M \times 1$ received signal vector can be expressed as

where ${\text{s}}(t) = {[\begin{array}{*{20}{c}} {{s_1}(t)}& ... &{{s_k}(t)} \end{array}]^{\rm{T}} }$, s_k(t) is the incident signal from the $k {\rm{th}}$ source, ${(\cdot)^{\rm{T}}}$ denotes the transpose; A= [a(θ₁) ... $ \begin{array}{*{20}{c}}{{\text{a}}({\theta _k})} \end{array}]$, ${\text{a}}({\theta _k})$ is the array manifold vector, ${\theta _k}$ denotes the signal arriving from the direction of the $k {\rm{th}}$ source; and ${\text{n}}(t)$ denotes the noise vector. For plane wave propagation, we have

where $\lambda $ is the wavelength. If we suppose that the noise is zero-mean white Gaussian with variance $\sigma _n^2$ and uncorrelated with the source signals, then the covariance of the received signal vector is

where ${{\text{I}}_M}$ denotes an $M \times M$ identity matrix, ${(\cdot)^{\rm{H}}}$ denotes the complex conjugate transpose; and ${{\text{R}}_{{\text{ss}}}}$ is the source covariance matrix of rank $K$. Estimating the covariance requires sufficient statistical samples of both the signals and noise, which is particularly demanding when the number of sensors is large. In addition, the computational complexity is related to the number of sensors.

First, we divide the entire array into $N$ identical subarrays, each with $P$ sensors that are equally separated by $D$, which is an integer multiple of the sensor spacing. We denote the set of $L$ beam pointing angles as ${\bar {\text{θ}} _l}, \; l = 1, ..., L$. The ESPRIT method with multiple-angle subarray beamforming can then be applied.

The commonly used total least squares (TLS) version of ESPRIT with multiple-angle subarray beamforming is denoted as TLS-ESPRIT-SB and its implementation can be summarized as follows (Jiang, 2011).

(1) Divide the entire array into $N$ identical subarrays, each with $P$ sensors that are equally separated by $D$, and obtain the signal received ${{\text{x}}_i}(t), \; i = 1, ...,\, N$ by each subarray. We denote a set of $L$ beam pointing angles as ${\bar {\text{θ}} _l}, \; l = 1, ...,\, L$.

(2) Obtain the beamforming outputs ${{\text{y}}_i}(t), \; i = 1, ...,\, N$ from each subarray, and compute ${\text{z}}(t) = {[\begin{array}{*{20}{c}} {{\text{y}}_1^{\rm{T}} (t)}& ..., &{{\text{y}}_N^{\rm{T}} (t)} \end{array}]^{\rm{T}} }$.

(3) Obtain the sample covariance ${{\hat{\text {R}}}_{{\text{zz}}}}$ as an estimate of ${{\text{R}}_{{\text{zz}}}}$.

(4) Compute the eigendecomposition of ${{\hat{\text {R}}}_{{\text{zz}}}}$ as

(5) If necessary, estimate the number of sources $K$.

(6) Partition ${\hat{\text {E}}}{\rm{ = [}}{{\hat{\text {E}}}_{\text{s}}}|{{\hat{\text {E}}}_{\text{n}}}{\rm{]}}$, where ${{\hat{\text {E}}}_{\text{s}}}$ are the principal eigenvectors corresponding to the $K$ largest eigenvalues.

(7) Compute the eigendecomposition of the $2K \times 2K$ matrix as

where ${{\hat{\text {E}}}_{{{\rm{z}}_{\rm{1}}}}}$ and ${{\hat{\text {E}}}_{{{\rm{z}}_{\rm{2}}}}}$ pick the first and last $L(N - 1)$ rows of ${{\hat{\text {E}}}_{\text{s}}}$, respectively.

(8) Partition ${\text{E}}$ into $K \times K$ submatrices, as follows:

(9) Calculate the eigenvalues ${\gamma _k}$ of $ - {{\text{E}}_{12}}{\text{E}}_{22}^{ - 1}$, then

where $\arg (\cdot)$ retrieves the angle information.

Note that a simple choice of $L = 2$ or 3 is often sufficient to improve the beamspace DOA estimation performance. The dimension in the beamspace can be far less than that in the element space, thereby reducing the computational complexity and decreasing the number of snapshots needed to obtain a statistically robust covariance estimate. Hence, for $L \ll P$, the ESPRIT-SB exhibits significant computational advantages compared with the original ESPRIT method.

In practical multibeam system applications, the above described method can be applied to the selected beams following conventional beamforming to generate beams at a set of pre-selected angles and for subsequent source signal detection at each beam.

The implementation of HRBD in the proposed MBES system for the beamforming outputs and sensor outputs from every time slice can be summarized as follows.

(1) Estimate the angle range of the receiving signals and determine the corresponding subarray. Choose an angle sector ${Q_h}$ and determine the structure of the virtual ULA and the corresponding transformation matrices $\left\{ {{{\text{B}}_h}} \right\}$, which are pre-computed off-line.

(2) Obtain the sensor output vector ${{\text{y}}_{\rm{v}}}$ for the virtual ULA with ${{\text{y}}_{\rm{v}}} = {{\text{B}}_h}{\text{y}}$, where ${\text{y}}$ is the sensor output vector from the subarray. If necessary, estimate the number of sources $K$.

(3) Choose $L$ angles ${\bar {\text{θ}} _l}, \; l = 1, ..., \, L$ in the angle sector ${Q_h}$ as the beam pointing angles and separate the virtual ULA into ${N_{\rm{v}}}$ subarrays.

(4) Operate TLS-ESPRIT-SB as before and then obtain the estimated DOA.

(5) Apply motion compensation to the estimated DOA and calculate the estimated TOA.

To assess the performance of the ESPRIT-SB for a U-shaped array (which we denote as via-ESPRIT-SB), we considered an example of DOA estimation using the TLS versions of both via-ESPRIT (we denote the ESPRIT with the virtual array transformation as via-ESPRIT) and via-ESPRIT-SB. The number of sources was assumed to be known where the individual source processes were uncorrelated and corrupted by white Gaussian noise. The simulation was evaluated based on 1 000 Monte-Carlo runs.

Two equal-power incident signals came from ${\theta _1} = {34.6^ \circ }$ and ${\theta _2} = {36.8^ \circ }$. For the via-ESPRIT-SB, we divided the virtual ULA into $N({\rm{ = 4)}}$ subarrays with $P({\rm{ = 49)}}$ sensors, where each was equally separated by 5d and three given-pointing beams were formed. For the via-ESPRIT, the virtual ULA was divided into two subarrays, where each had 63 sensors and they were equally separated by d. Thirty snapshots were used for both algorithms.

Average root mean squared errors (RMSE) for all of the sources are plotted as a function of the signal-to-noise ratio (SNR) in Fig. 4. When the SNR approached –15 dB, the RMSE approached the variance of a uniform random variable (Van Trees, 2002). As expected, the RMSE exhibited a threshold behavior, where it increased dramatically below some SNR due to ambiguity sidelobes (Van Trees, 2002). A SNR gain in the order of $10\log P$ was clearly obtained, where the via-ESPRIT entered the threshold region earlier and the via-ESPRIT-SB performed better in the high SNR region.

Clearly, the virtual transformation performed efficiently when operating high resolution methods on a U-shaped array. The performance of the via-ESPRIT-SB was much better with a few snapshots compared with that of the via-ESPRIT for the following two reasons. First, the subarray configurations changed the performance where the via-ESPRIT-SB could increase the separation among the subarrays without any prior information to improve the performance. Second, the beamforming processing in the via-ESPRIT-SB could obtain a SNR gain and lower the threshold. Therefore, the via-ESPRIT-SB is more suitable for use in practical applications with a nonlinear array.

An experiment was performed in the Qiandao Lake of Hangzhou, China. The bottom of this lake has a bowl shape and the depth decreases gradually from the edge to the middle of the lake, where it changes from about 10 to 60 m. The MBES system was installed underwater at a depth of 3 m in the well of a ship and the ship was anchored far off the shore. The system could be rotated 180° in the horizontal plane. The transmitting frequency was 165 kHz and the transmitting pulse width was 0.5 ms. The sample frequency for receiving was 4/3 times the transmitting frequency (Xu et al., 2009). Each record file comprised the receiving data according to one emission that included six pings. The MBES system formed 290 uniform receiving beams spanning from –85° to 85°. We assumed that the acoustic speed was a simple, fixed value such as 1 500 m/s.

The receiving signals were demodulated and a sampling rate for the baseband signal was 20 kHz. The results obtained by conventional beamforming processing using the baseband signal for one ping are shown in Fig. 5, where the abscissa is the beam angle and the ordinate is the number of snapshots. The value in the color bar denotes the quantization value of the energy after beamforming. The darker colored spots show the rough contours of the bottom. The BDI and the HRBD were applied to the beamforming outputs. A threshold of three was selected, so a signal with an energy three times higher than the average energy in this beam (the SNR was about 5 dB) was detected as the bottom signal. The estimates of the bottom obtained by the two methods are plotted in Figs 6a and b. In Fig. 6, the abscissa is the number of snapshots, where it corresponds to the TOA from the bottom, and the ordinate is the beam angle. The major differences between Figs 6a and 6b are in terms of the normal and oblique incident angles. In brief, at [–20° 20°], the data obtained by the whole array of sensors were used by the BDI, but only data from half the array of sensors were used for the virtual array transformation by the HRBD. The aperture of the array determines the resolution of beam forming and the number of array elements is related to the array gain. The results show that the HRBD obtained fewer hits than the BDI (Fig. 6). Thus, the performance of the HRBD was poorer in the case of normal incidence. At [–85° –70°] and [70° 85°], the hits in Fig. 6b formed a better trace than those in Fig. 6a, and the difference was even more significant in Fig. 7. Figure 7 shows the depth estimated based on the hits from three pings after bottom detection, where it was assumed that the acoustic speed was fixed, and the red and the black indicate the results obtained by the BDI and the HRBD, respectively. The estimated depths were the original results without the elimination of outliers and curve fitting. The left of Fig. 7 shows that when the range was more than 50 m, the beam angle was larger than ${70^ \circ }$. The BDI missed the lake bottom whereas the HRBD continued following. The right of Fig. 7 shows that when the range was more than 150 m, the beam angle was larger than ${75^ \circ }$. In addition, the depth in another direction is shown in Fig. 8 where the MBES was rotated nearly 100°. According to Fig. 8, there were two differences between the results obtained by the two algorithms: the results obtained by the HRBD were on the either side of the graph and the variance of the results obtained by the BDI was larger than that of the HRBD results, where both were caused by the higher resolution of the HRBD. When the beam angle was larger, the footprint area was wider and the resolution of the BDI was lower, and thus the variance of the hits was greater with the BDI. As shown in Fig. 5, the footprint spanned nearly 3° with 1 000 snapshots and the depth results obtained by the BDI ranged over more than 5 m. The HRBD followed the bottom better and the range was more than 50 m on both sides. In the edge region, the HRBD had a greater capacity for terrain following.

In summary, the performance of the HRBD in the normal incident area was poorer than that of the BDI due to the employment of aperture segmentation. However, on the edges of the beams, the HRBD could estimate the DOAs better and provide more details of the bottom depth. The HRBD method is more sensitive to the SNR, thereby leading to more false hits, and more accurate discrimination of outliers is required. However, the false hits and outliers can be eliminated easily after depth tracking according to the depth information coherence.