Passive Array-Invariant-Based Localization for a Small Horizontal Array Using Two-Dimensional Deconvolution

Abstract

Recently, the array-invariant method was proposed to passively localize sources of opportunity in shallow water. It exploits multiple arrivals which are different in terms of beam angle and travel time. Conventional plane-wave beamforming in the existing array-invariant method is used to obtain beam-time migration. The resolution capability of conventional plane-wave beamforming is determined by array aperture, which, however, limits the localization accuracy of the existing array-invariant method. To improve the localization accuracy, this study proposes the use of two-dimensional (2D) deconvolution to obtain a better beam-time migration than in conventional plane-wave beamforming. Our simulation with a small horizontal array showed that the range estimation error of the proposed array-invariant method based on 2D deconvolution was only one-third of that of the existing method. The experiment also demonstrated the validity of our proposed method.

1. Introduction

Underwater source localization is a hot topic in underwater acoustic signal processing [1,2,3,4,5,6,7,8]. The array-invariant method, based on the dispersion characteristics of wideband signals in waveguides, was proposed by Lee and Makris in 2006 for shallow sea source localization [2]. With the development of blind deconvolution [3,4,5], the array-invariant method has been successfully applied to passively localize surface ships or submerged sources, using a vertical [6,7,8,9] or horizontal [10,11] array without knowledge of the ocean environment. The blind deconvolution method extracts the Green’s function, namely the channel impulse response, between a source of opportunity (e.g., ship) and an array. Subsequently, conventional plane-wave beamforming is employed to obtain the beam-time migration of the Green’s function. Finally, the array-invariant method estimates the source range from the beam-time migration using the least-squares approach [2,11].

Following the steps of the existing array-invariant method, it can be found that beam-time migration plays a key role in the source–range estimation. To date, only conventional plane-wave beamforming has been utilized to obtain the beam-time migration in localization applications [7,8,9,11]. However, the resolution capability of conventional plane-wave beamforming is determined by the array aperture, which influences the localization accuracy of the existing array-invariant method, especially when using a small horizontal array. The relatively wide beam-time migration should be one of the main reasons causing the source–range estimation error of the existing array-invariant method. This study attempts to decrease the source–range estimation error by obtaining better beam-time migration by using two-dimensional (2D) deconvolution.

In this paper, with some mathematical derivations, we employ frequency-domain beamforming to process the extracted Green’s function to obtain the beam-time migration, which is expressed as a 2D convolution. Based on the deduced 2D convolution expression, we propose a 2D deconvolution method to deconvolve the beam-time migration obtained with conventional plane-wave beamforming to obtain the high-resolution beam-time migration. The more accurate source range can thus be estimated from the high-resolution beam-time migration. The existing array-invariant method [11] involving blind deconvolution and conventional plane-wave beamforming is referred to as the reference method in this study. The simulation with a small horizontal array in shallow water shows that the source–range estimation error of our proposed method based on the 2D deconvolution is only one third of that of the reference method. The experiment with a small horizontal array in the Yellow Sea of China also proves the validity of our proposed method. Additionally, the results in this study also demonstrate that the wide beam-time migration obtained with conventional plane-wave beamforming in the array-invariant method is one of the main reasons behind source–range estimation errors.

This brief report is organized as follows: Section 2 describes the Green’s function and blind deconvolution for the horizontal array. In Section 3.1, conventional plane-wave beamforming is reviewed. In Section 3.2, the 2D deconvolution method is proposed to deconvolve the beam-time migration obtained with conventional plane-wave beamforming. The 2D deconvolution method can obtain high-resolution beam-time migration. The simulation and the experiment are described in Section 4 and Section 5, respectively. The simulation and the experiment demonstrate the validity of our proposed method. Finally, our conclusion is provided in Section 6.

2. Extracting Green’s Function

For the array-invariant method, extracting the Green’s function with blind deconvolution [8,9,11] from the received signals is necessary. Assume that an unknown source signal is $s(t)$. The frequency spectrum of $s(t)$ is $S(\omega )=\left|S(\omega )\right|e^{j{\mathsf{\Phi }}_s(\omega )}$, where $\omega$ is the radian frequency and ${\mathsf{\Phi }}_s(\omega )$, is the corresponding phase. In this study, we focus on using a horizontal array in shallow water to realize passive localization. Consider a horizontal array with N elements, as shown in Figure 1. The received signal of the nth element can be expressed as

$$P_n(\omega )=G(r_n,r_s,\omega )\left|S(\omega )\right|e^{j{\mathsf{\Phi }}_s(\omega )},$$

(1)

where $r_n$ and $r_s$ are the position vectors of the nth element of the array and the unknown source, respectively, and $G(r_n,r_s,\omega )$ is the Green’s function between the source at $r_s$ and the nth element at $r_n$. In a ray path approximation [11], the Green’s function $G(r_n,r_s,\omega )$ can be written as

$$G(r_n,r_s,\omega )={\displaystyle \sum _{k=1}^K{A}_k{e}^{j\omega \{{\tau }_n({\theta }_k,\varphi )-T({\theta }_k,\varphi )\}}},$$

(2)

where $\varphi$ is the horizontal azimuth shown in Figure 1, ${\theta }_k$ is the grazing angle of the kth ray, ${\tau }_n({\theta }_k,\varphi )$ is the local time delay of the kth ray at the nth element, $T({\theta }_k,\varphi )$ is travel times of the kth ray, and $A_k$ is the amplitude of the kth ray.

For the horizontal array in Figure 1, the local time delays depend mainly on the horizontal azimuth due to a small grazing angle, thus ${\tau }_n({\theta }_k,\varphi )\approx {\tau }_n(\varphi )$ [11]. Subsequently, by steering the beam to $\varphi$, the output of the conventional plane-wave beamforming is identified by

$$\begin{array}{ll}F(\omega ;{\theta }_k,\varphi )& ={\displaystyle \sum _{n=1}^N{e}^{-j\omega {\tau }_n({\theta }_k,\varphi )}{P}_n(\omega )}\\ & =\left|S(\omega )\right|e^{j{\mathsf{\Phi }}_s(\omega )}{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{k_q=1}^K{A}_{k_q}{e}^{j\omega \{{\tau }_n({\theta }_{k_q,}\varphi )-T_{k_q}\}}}}{e}^{-j\omega {\tau }_n({\theta }_{k,}\varphi )}\\ & \approx \left|F(\omega ;{\theta }_k,\varphi )\right|e^{\{j{\mathsf{\Phi }}_s(\omega )-j\omega T_k\}}\end{array}$$

(3)

Blind convolution uses the beamforming output phase, $\psi (\omega ,\varphi )={\mathsf{\Phi }}_s(\omega )-\omega T({\theta }_k,\varphi )$, and removes the phase component ${\mathsf{\Phi }}_s(\omega )$ of the unknown source signal from $P_n(\omega )$ by a phase rotation, expressed as:

$$\widehat{G}(r_n,r_s,\omega )=\frac{P_n(\omega )}{\sqrt{{{\displaystyle \sum _{n=1}^N\left|P_n(\omega )\right|}}^2}}{e}^{-j\psi (\omega ,\varphi )}=\frac{G(r_n,r_s,\omega )}{\sqrt{{{\displaystyle \sum _{n=1}^N\left|G(r_n,r_s,\omega )\right|}}^2}}{e}^{j\omega T({\theta }_k,\varphi )}$$

(4)

The beam-time migration in the existing array-invariant method [11] can be obtained by processing the extracted $\widehat{G}(r_n,r_s,\omega )$ with conventional plane-wave beamforming.

3. 2D Deconvolution for Obtaining Better Beam-Time Migration

To ensure our proposed method based on 2D deconvolution is easily understandable, frequency-domain conventional plane-wave beamforming is briefly introduced first in this section. Then, our proposed method is presented.

3.1. Conventional Plane-Wave Beamforming

Considering a far-field source, the angle of incidence $\tilde{\theta }$ is shown in Figure 1. When using uniform horizontal array, as shown in Section 2, with interelement spacing denoted by d, the received signal of the array can be expressed as

$$X_S(\omega )=p_{a}S(\omega )+n=\left[\begin{array}{c}{p}_{a,1}\\ p_{a,2}\\ ⋮\\ p_{a,N}\end{array}\right]S(\omega )+n$$

(5)

where $p_a$ is the array response vector given as $p_a={[p_{a,1},p_{a,2},\cdots p_{a,N}]}^T$, $p_{a,n}=e^{-j\omega nd\sin \tilde{\theta }/c}$, and $n$ is the noise vector.

Conventional beamforming [12,13,14] multiplies the array signals $X_S(\omega )$ by a steering vector $\mu =[{\mu }_1,{\mu }_2,\cdots {\mu }_N{]}^T$ for the direction $\phi$, where ${\mu }_n=(1/N)e^{-j\omega nd\sin \phi /c}$. The beam power as a function of $\sin \phi$, which is obtained with conventional beamforming, is written as

$$X_p(\sin \phi )={\left|X_S(\omega ){\mu }^T\right|}^2={\left|\frac{S(\omega )}{N}{\displaystyle \sum _{n=1}^N{e}^{j2\pi n\frac{d}{\lambda }(\sin \phi -\sin \tilde{\theta })}}\right|}^2={\left|S(\omega )\right|}^2{\left|\frac{1}{N}{\displaystyle \sum _{n=1}^N{e}^{j2\pi n\frac{d}{\lambda }(\sin \phi -\sin \tilde{\theta })}}\right|}^2$$

(6)

3.2. 2D Deconvolution

The local time delays ${\tau }_n({\theta }_k,\varphi )$ in Equation (2) can be written as ${\tau }_n(\sin {\tilde{\theta }}_k)$, which is the function of the beam angle measured $\sin {\tilde{\theta }}_k$, and $T({\theta }_k,\varphi )$ can be abbreviated as $T_k$. Then, the expression of the Green’s function in (2) can be rewritten as

$$G(r_n,r_s,\omega )={\displaystyle \sum _{k=1}^K{A}_k{e}^{j\omega \{{\tau }_n(\sin {\tilde{\theta }}_k)-T_k\}}}$$

(7)

Let $G(r_s,\omega )=[G(r_1,r_s,\omega ),G(r_2,r_s,\omega ),\cdots G(r_N,r_s,\omega )]$. We define the 2D rotation matrix as

$$\Psi =[\mu (t_0),\mu (t_1),\dots \mu (t_M)],$$

(8)

where $\mu (t_m)=[{\mu }_1(t_m),{\mu }_2(t_m),\cdots {\mu }_N(t_m){]}^T$, and ${\mu }_n(t_m)=(1/N)e^{j\omega ({\tau }_n(\sin \phi )-t_m)}$.

Assuming that we process M samples, $t_0$ and $t_M$ are the start and end time of the sampling segment, respectively. Considering (7), after the time-beam rotation, the output can be written as

$$\begin{array}{c}Y(\sin \phi ,t)={\displaystyle \sum _{{\omega }_q={\omega }_{\min }}^{{\omega }_{\max }}G(r_s,{\omega }_q)}\Psi ={\displaystyle \sum _{{\omega }_q={\omega }_{\min }}^{{\omega }_{\max }}{\displaystyle \sum _{k=1}^K\frac{A_k}{N}}{\displaystyle \sum _{n=0}^N{e}^{-j{\omega }_q(({\tau }_n(\sin \phi )-{\tau }_n(\sin {\tilde{\theta }}_k))+(T_k-t))}}}\\ ={\displaystyle \sum _{k=1}^K\frac{A_k}{N}{\displaystyle \sum _{{\omega }_q={\omega }_{\min }}^{{\omega }_{\max }}{\displaystyle \sum _{n=1}^N{e}^{-j{\omega }_q((\frac{nd(\sin \phi -\sin {\tilde{\theta }}_k)}{c}+(T_k-t))}}}}\end{array}$$

(9)

where ${\omega }_{\min }$ and ${\omega }_{\max }$ are the highest and lowest frequencies of the Green’s function $G(r_s,{\omega }_q)$ in the discrete frequency domain, respectively. The power of $Y(\sin \phi ,t)$ is

$$B(\sin \phi ,t)={\left|Y(\sin \phi ,t)\right|}^2={{\displaystyle \sum _{k=1}^K\frac{{\left|A_k\right|}^2}{N^2}\left|{\displaystyle \sum _{{\omega }_q={\omega }_{\min }}^{{\omega }_{\max }}{\displaystyle \sum _{n=1}^N{e}^{-j{\omega }_q((\frac{nd(\sin \phi -\sin {\tilde{\theta }}_k)}{c}+(T_k-t))}}}\right|}}^2.$$

(10)

After a simple derivation, Equation (10) can be written as a 2D convolution of the beam-time pattern, denoted as PSF, with the source power distribution $S$. $PSF$ depends only on the difference in bearing and time, so $PSF$ is shift-invariant in the steering angle and time, and can be written as

$$B(\sin \phi ,t)={\displaystyle {\int }_{t_0}^{t_M}{\displaystyle {\int }_{-1}^{1}PSF((\sin \phi -\sin \vartheta ),(t-\widehat{t}))}\cdot S(\sin \vartheta ,\widehat{t})d\sin \vartheta d\widehat{t}},$$

(11)

where

$$PSF((\sin \phi -\sin \vartheta ),(t-\widehat{t}))={\left|\frac{1}{N}{\displaystyle \sum _{{\omega }_q={\omega }_{\min }}^{{\omega }_{\max }}{\displaystyle \sum _{n=0}^N{e}^{-j{\omega }_q((\frac{nd(\sin \phi -\sin \vartheta )}{c}+(t-\widehat{t}))}}}\right|}^2,$$

(12)

and

$$S(\sin \vartheta ,\widehat{t})={\displaystyle \sum _{k=1}^K{\left|A_k\right|}^2\delta ((\sin \vartheta -\sin {\tilde{\theta }}_k)+(\widehat{t}-T_k))}.$$

(13)

It should be noted that $S(\sin \vartheta ,\widehat{t})$ separates the multiple arrivals, in terms of both the bearing and time domains, much better than $B(\sin \phi ,t)$ does. If we estimate $S(\sin \vartheta ,\widehat{t})$ from $B(\sin \phi ,t)$ using 2D deconvolution, high-resolution beam-time migration can be obtained, which helps improve the localization accuracy.

Some well-known deconvolution algorithms, such as the Richardson–Lucy (R–L) [12,13,14,15,16] algorithm, have been widely applied in various fields. We propose the use of the R–L algorithm to deconvolve $B(\sin \phi ,t)$ both the bearing and time domains to improve the localization accuracy of the array-invariant method.

After obtaining the high-resolution beam-time migration, the range is estimated with the following formula:

$$r_0=\frac{-c\sin \tilde{\theta }}{\tilde{\chi }}$$

(14)

where $\tilde{\chi }=d(\sin \tilde{\theta })/dt$, which is extracted from the beam-time migration with the least-squares approach [2,11]. The source direction $\tilde{\theta }$ is estimated by searching for the maximum value in the beam-time migration.

4. Simulation

A simulation in a typical shallow-water waveguide was conducted to demonstrate our proposed method based on 2D deconvolution. The KRAKEN toolbox was used. The environmental parameters of the employed Pekeris [2,17] waveguide are summarized in Figure 2a. The water depth was 100 m and the speed of sound was assumed to be 1500 m and the same at different depths. The seabed was typically sand bottom.

A horizontal array with 16 elements was employed. The array length was 2.8 m. The bandwidth of the noise radiated by the source of opportunity was 2–4 kHz. The sampling frequency was 20 kHz. The existing array-invariant method [11] was used as the reference method. The source and the horizontal array were placed at 3 m and 97 m under water, respectively. The azimuth angle of the sound source was 65°. Additive white noise was used. The signal-to-noise ratio at each array element was set at 10 dB. The length of simulation data is 2 s. The Green’s function of the 16th element at the range of 1000 m was extracted with blind deconvolution and is shown in Figure 2b.

The beam-time migrations obtained using the reference method and proposed methods are shown in Figure 2c,d, respectively, which have the same scale. The peak position of beam-time migrations is marked by the red square. In Figure 2c, it can be seen that the relative arrival time of the two paths was 13.95 ms and 40.1 ms, respectively, and the sine of the beam angle of the two paths was 0.879 ms and 0.839 ms, respectively. In Figure 2d, it can be seen that the relative arrival time of the two paths was 13.2 ms and 40.1 ms, respectively, and the sine of the beam angle of the two paths was 0.883 ms and 0.849 ms, respectively.

It is clear that our proposed method based on 2D deconvolution, shown in Figure 2d, had a much higher resolution in both the bearing and time domains, compared to the reference method. Based on the high-resolution beam-time migration, it is clear that the source–range estimation error can be decreased.

For passive shallow-water localization with a horizontal array, the relative error of range estimation determines the localization accuracy. The source range was estimated according to Equation (14). In the case of the true source range of 1000 m, the range estimated with the reference method, shown in Figure 2c, was 850.2 m. The range estimated with the proposed method, shown in Figure 2d, was 1038.7 m. The relative estimation errors of the proposed and reference methods were 3.9% and 15.0%, respectively. The estimation performances of the two methods at different ranges are summarized in

Table 1. The range estimates and relative errors of the proposed and reference methods at different ranges in the simulation.

True Range	Proposed (Estimated Range, Relative Error)	Reference (Estimated Range, Relative Error)
1000 m	1038.7 m, 3.9%	850.2 m, 15.0%
1500 m	1615.1 m, 7.7%	1216.3 m, 18.9%
2000 m	2146.8 m, 7.3%	2454.3 m, 22.7%
3000 m	3232.3 m, 7.7%	2251.0 m, 24.9%

. On average, the estimation error of our proposed method was only one third of that of the reference method. The results also prove that the wideband beam-time migration obtained using conventional plane-wave beamforming in the reference array-invariant method [11] is one of the main factors influencing the accuracy of the source range estimation.

5. Experiment

The data from an experiment in the Yellow Sea of China in 2019 were processed for testing the proposed and reference methods. A 2.8 m movable horizontal array was used. The horizontal array had 16 elements. The sampling frequency was 20 kHz. The source was fixed as shown in Figure 3a. The source transmitted wideband noise within a frequency band of 2–4 kHz. Both the source and the horizontal array were located at a depth of 5 m. The trajectory of the horizontal array is shown in Figure 3a.

The extracted Green’s function at Position #1 (988.9 m) in Figure 3a is shown in Figure 3b. The beam-time migrations obtained using the reference and proposed methods are shown in Figure 3c,d, respectively, which have the same scale. The peak position of beam-time migrations is marked by the red square. Based on the beam-time migrations in Figure 3c,d, the range estimated with the reference and proposed methods were 1244.3 m and 1130.2 m, respectively, and the relative estimation errors of the reference and proposed methods were 25.9% and 14.4%, respectively.

The received signals at Positions #2 (718.6) and #3 (580.1) in Figure 3a were also processed with the proposed and reference methods.

Table 2. The range estimates and relative errors of the proposed and reference methods at different ranges in the experiment.

True Range	Proposed (Estimated Range, Relative Error)	Reference (Estimated Range, Relative Error)
#1 (988.9 m)	1130.2 m, 14.4%	1244.3 m, 25.9%
#2 (718.6 m)	604.2 m, 15.9%	535.3 m, 25.3%
#3 (580.1 m)	515.2 m, 11.2%	425.5 m, 26.6%

shows the range estimates and relative errors of the proposed and reference methods at different positions. It can be seen that the relative error of the proposed method is about half of that of the reference method.

6. Conclusions

This brief report proposed an array-invariant method based on 2D deconvolution for robust passive source localization in shallow water. Our proposed method achieved much better beam-time migration than the existing array-invariant method, which improved the range estimation accuracy significantly. Both the simulation and the experiment show that our proposed method based on 2D deconvolution achieves a much higher accuracy than the existing array-invariant method. This brief report also reveals that the relatively wideband beam-time migration obtained using conventional plane-wave beamforming in the existing array-invariant method is one of the main factors that influences the localization accuracy.