A Nonlinear Data-Driven Towed Array Shape Estimation Method Using Passive Underwater Acoustic Data

Abstract

Large-aperture towed linear hydrophone array has been widely used for beamforming-based signal enhancement in passive sonar systems; however, its performance can drastically decrease due to the array distortion caused by rapid tactical maneuvers of the towed platform, oceanic currents, hydrodynamic effects, etc. In this paper, an enhanced data-driven shape array estimation scheme is provided in the passive underwater acoustic data, and a novel nonlinear outlier-robust particle filter (ORPF) method is proposed to acquire enhanced estimates of time delays in the presence of distorted hydrophone array. A conventional beamforming technique based on a hypothetical array is first used, and the detection of the narrow-band components is sequentially carried out so that the corresponding amplitudes and phases at these narrow-band components can be acquired. We convert the towed array estimation problem into a nonlinear discrete-time filtering problem with the joint estimates of amplitudes and time-delay differences, and then propose the ORPF method to acquire enhanced estimates of the time delays by exploiting the underlying properties of slowly changing time-delay differences across sensors. The proposed scheme fully exploits directional radiated noise targets as sources of opportunity for online array shape estimation, and thus it requires neither the number nor direction of sources to be known in advance. Both simulations and real experimental data show the effectiveness of the proposed method.

1. Introduction

Beamforming-based signal enhancement technique in a passive sonar system is a hot topic in underwater acoustic array signal processing, and has been paid much attention by researchers. Large-aperture towed linear hydrophone array, which is typically formed by trailing a hydrophone array behind a towing platform in a nominally straight line, is often used to perform high-resolution direction-of-arrival (DOA) estimation and signal enhancement with high spatial gain [1,2]. However, the towed linear array would be often perturbed and deviated from the nominally straight line, due to inevitable oceanic currents, hydrodynamics, and tactical maneuvers of the towing platform [3,4,5]. It is a significant challenge for beamforming-based signal enhancement based on the known array shape in the underwater acoustic array signal processing [3,4,5,6,7].

A huge amount of effort has been devoted to array shape estimation. These methods can be generally classified into source-independent methods and source-dependent methods. For source-independent methods, one naive method is to install depth sensors and compasses on the array for the acquisition of vertical and horizontal information of the transverse displacements of the hydrophone sensors. The translation invariance of the array curve is used to establish and solve the underdetermined equation of the segment curve on the array state function with nonacoustic auxiliary sensors [8]. However, such solutions are not economical due to the burden of extra sensors. In addition, the limited precision of these auxiliary sensors makes it difficult to accurately estimate the array shape in real-time application [9,10]. The source-dependent methods, in essence, measure the relative time or phase delays among sensors by using one or more sources with known or unknown directions based on underwater acoustic data [11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28]. The single and double near-field source-dependent method is used for towed array shape estimation [11]. Error sources in the localization systems are analyzed to reduce the impact on final estimation of array shape [12]. The Kalman filters theoretical model is used to address the towed array shape estimation [13]. The approximately optimal distribution of a number of depth sensors over the discretized hydrophone sensors yields the approximately best achievable performance in terms of the minimum space average mean square error [14]. A hidden Markov model is used to estimate the slowly changing array shape, and the Viterbi algorithm enables a maximum likelihood estimate of the array shape to be obtained efficiently [15]. The Cramer–Rao lower bound (CRLB) for the narrow-band plane-wave source is used for the array shape estimation with a single known source direction [16]. Subspace-based eigenvector techniques are introduced to estimate array shapes with a known single source at a known frequency and direction [17]. An iterative maximum likelihood method for simultaneously estimating directions of arrival (DOA) and sensor locations is developed [18]. Both the acoustic hydrophone data and sensor data are used for the distorted array shape estimation [25]. The phase of the narrow-band frequency in the passive acoustic data, which involves the information of time delays in the presence of array distortion, can be used to perform the array shape estimation. The phases in the multiple narrow-band frequencies from directional radiated noise source signals as sources of opportunity were used to acquire time-delay estimation [26]; an improved time-delay estimation method based on the phases of multiple narrow-band frequencies is developed by exploiting the underlying properties of slowly changing time-delay difference in the hydrophone array [27]; furthermore, the time-delay estimation method based on the phases of multiple narrow-band frequencies was developed for the improvement of hydrophone array estimation by exploiting the underlying properties of slowly-changing time-delay difference over time, and the Chinese remainder theorem technique is used to effectively address the phase difference ambiguity issue [28]. However, the proposed method directly uses the complex values of the detected narrow-band frequencies for the time-delay difference estimation, and effectively avoids the issue of the phase-difference ambiguity issue, which often occurs in the phase extraction operator, especially in the low signal-to-noise-ratio (SNR) case.

In this paper, an enhanced data-driven array shape estimation scheme is proposed using passive underwater acoustic data, and a novel nonlinear outlier-robust particle filter (ORPF) method is proposed to acquire enhanced estimates of time delays in the presence of distorted hydrophone array. Considering the fact that the radiated noise sources have an amount of narrow-band components that originated from the vibration of mechanical equipment, we carry out the narrow-band detection for the acquisition of the information of time delays after beamforming-based signal enhancement in the hypothetical uniform linear array shape. Unlike these methods in [26,27,28], where only phase information of the multiple narrow-band frequencies is exploited to perform the time-delay estimation, a nonlinear ORPF method is proposed to acquire improved estimates of time delays by jointly exploiting the amplitude and phases of detected multiple narrow-band frequencies, and the posterior probability distribution inference is provided based on the expectation–maximization technique. The proposed scheme fully exploits the directional radiated noise source signals as sources of opportunity for online array shape estimation, and thus it requires neither the number nor direction of sources to be known in advance.

The contributions of the paper lie in three aspects: (1) An enhanced data-driven array shape estimation scheme using passive underwater acoustic data is provided; it exploits the directional radiated noise source signals as sources of opportunity for online array shape estimation, and thus it requires neither the number nor direction of sources to be known in advance. (2) An enhanced nonlinear ORPF method is proposed by jointly exploiting the amplitude and phases of detected multiple narrow-band frequencies. (3) The posterior probability distribution inference of the proposed method is provided based on the expectation–maximization technique.

The paper is organized as follows: the passive underwater acoustic signal model in the distorted towed array is first introduced in Section 2. The distorted towed hydrophone array shape estimation method is provided in Section 3. In this section, the time delays in the narrow-band components are analyzed; the outlier-robust particle filter model is given, and the posterior probability distribution inference is provided. Simulations and experiments demonstrate the robustness of the proposed algorithms compared to conventional methods in Section 4. The summarization of this paper is presented in Section 5.

Notation: We use lower-case (upper-case) bold characters to denote vectors (matrices). $p(·)$ denotes the probability density function (pdf), and $\mathcal{N}\left(x\right|a,b)$ denotes that random variable x follows a Gaussian distribution with mean a and variance b. $\mathrm{Gamma}\left(x\right|\alpha ,\beta )$ denotes that the random variable x follows the gamma function with the parameter $\alpha$ and $\beta$. ${(·)}^T$ and ${(·)}^H$ denote transpose and conjugate transpose, respectively, ${\mathbf{I}}_N$ denotes the $N\times N$ identity matrix. ${}^{*}$ donates a conjugate operator, and $\langle x\rangle$ represents the expectation operator of the random variable x. $lnx$ denotes the natural logarithm of x.

2. Passive Underwater Acoustic Signal Model in the Distorted Towed Array

In this section, the passive underwater acoustic signal model in a distorted towed array is introduced [26]. Supposing a towed linear array consisting of M hydrophone sensors, a uniform linear array (in red) with the interelement spacing d is shown in Figure 1. Due to the transverse motion of ship tactical maneuvering, oceanic currents, and hydrodynamic effects, it would result in transverse displacements of the array away from a straight line, and an example of distorted array shape is shown in Figure 1 (in blue). In the paper, the directional line formed by the first two sensors is defined to be the baseline in the $x-y$ plane, and x-axis is defined along the directional line, as shown in Figure 1. However, the distortion of the linear array would make the position of the hydrophone sensor locate away from the x-axis.

Due to the flexible material of the hydrophone array, the interelement spacing between adjacent sensors is constant even in the distorted array, and can be written by

$$\begin{array}{c}\hfill d=\sqrt{{(x_{m+1}-x_m)}^2+{(y_{m+1}-y_m)}^2},m=0,\cdots ,M-1,\end{array}$$

(1)

where the location of the mth sensor is denoted by $(x_m,y_m)$ in the $x-y$ plane.

Considering a far-field target as a source of opportunity impinging on the towed linear array with the bearing of $\varphi$, we have the range offset between interelement sensors in the hypothetical linear array

$$\begin{array}{c}\hfill \mathsf{\Delta }{R}_0=R_1-R_0\approx dcos\varphi =\left(x_1-x_0\right)cos\varphi ,\end{array}$$

(2)

where $R_1$ denotes the distance from the acoustic source to the second sensor, and $R_0$ denotes the distance from the acoustic source to the first sensor, which is often regarded as the reference. Assuming the deviation angle formed by the line between the mth sensor and $(m+1)$th sensor along the x-axis is denoted as ${\phi }_m$ in the distorted array, the range offset $\mathsf{\Delta }{R}_m$ becomes, in the presence of array distortion,

$$\begin{array}{c}\hfill \mathsf{\Delta }{R}_m\approx dcos(\varphi -{\phi }_m)=\left(x_{m+1}-x_m\right)cos(\varphi -{\phi }_m)+\left(y_{m+1}-y_m\right)sin(\varphi -{\phi }_m).\end{array}$$

(3)

Given the first sensor as the reference, the range offset between the mth sensor and the reference can be written by

$$\begin{array}{c}\hfill \mathsf{\Delta }{R}_m=R_m-R_0=\sum _{i=0}^{m-1}dcos(\varphi -{\phi }_i)=\sum _{i=0}^{m-1}\left(x_{i+1}-x_i\right)cos(\varphi -{\phi }_m)+\left(y_{i+1}-y_i\right)sin(\varphi -{\phi }_m)\end{array}$$

(4)

The corresponding time-delay offset is

$$\begin{array}{c}\hfill {\tau }_m=\frac{{\sum }_{i=0}^{m-1}\mathsf{\Delta }{R}_i}{c}=\sum _{i=0}^{m-1}\mathsf{\Delta }{\tau }_i\end{array}$$

(5)

where c is the speed of an underwater acoustic wave, and $\mathsf{\Delta }{\tau }_i$ denotes the acoustic source time-delay difference between the $(i+1)$th and ith sensors. It is observed that the accuracy of the array shape estimation will depend on the accuracy of the estimated delay difference $\mathsf{\Delta }{\tau }_i$ by the adjacent sensors, and the adjacent errors would result in the accumulation in the distorted array shape estimation.

Assuming the acoustic source signal as $s\left(t\right)$, the received signal ${\mathrm{sr}}_m\left(t\right)$ in the mth sensor can be written by

$$\begin{array}{c}\hfill {\mathrm{sr}}_m\left(t\right)={\alpha }_{m}s\left(t-{\tau }_m\right)+e_m\left(t\right)\end{array}$$

(6)

where ${\alpha }_m$ and ${\tau }_m$ represent, respectively, the amplitude and time delay from the source to the mth sensor, and $e_m\left(t\right)$ is additive Gaussian noise.

According to the statistic model of the underwater acoustic radiated noise, the time-domain $s\left(t\right)$ can be expressed by [24,29]

$$\begin{array}{c}\hfill s\left(t\right)=\left(1+p_{\mathrm{mod}}\left(t\right)\right)p_{\mathrm{cont}}\left(t\right)+p_{\mathrm{nar}}\left(t\right),\end{array}$$

(7)

where $p_{\mathrm{cont}}\left(t\right)$ denotes the hydrodynamic noise yielding a stationary continuous spectrum, $p_{\mathrm{mod}}\left(t\right)$ represents a modulation function caused by the propeller noise, and $p_{\mathrm{nar}}\left(t\right)$ denotes the narrow-band noise. The narrow-band signal $p_{\mathrm{nar}}\left(t\right)$, which is quite an important component in the radiated noise, is often generated from the vibration of mechanical components, such as diesel generators and air conditioning equipment, and can be denoted by the superposition of multiple sinusoid signals as

$$\begin{array}{c}\hfill p_{\mathrm{nar}}\left(t\right)=\sum _{l=1}^L{\gamma }_{l}sin\left(2\pi f_{l}t\right),\end{array}$$

(8)

where ${\gamma }_l$ is the amplitude of the lth single-tone frequency from the acoustic source of oppotunity, L is the total number of the narrow-band signals radiated from the acoustic source, and $f_l$ is the narrow-band frequency.

In the real passive sonar system, the narrow-band components generated from the inevitable vibration of mechanical equipment, such as the diesel generator and air conditioning system, are quite important and useful elements in the underwater ship-radiated noise [30]. Generally, the power of the narrow-band components are at least 10 dB higher than that of their nearby continuous spectrum such that the narrow-band components are ready to be detected [31]. Figure 2 shows an example of the real radiated noise spectrum, and it is observed that the radiated noise spectrum involves multiple clear narrow-band components due to their prominent noise levels. The peaks with the red stars could be the narrow-band components from the vibration of the mechanical equipment of the ship. It also demonstrates the claim that the narrow-band components are often higher than that of their local continuous spectrum and are easily detected. Therefore, it provides a feasible way to estimate the distorted hydrophone array shape by exploiting these narrow-band components.

3. Distorted Towed Hydrophone Array Shape Estimation Method

In this section, an enhanced data-driven array shape estimation scheme is proposed to acquire improved array shape estimation by exploiting the narrow-band components, and a novel nonlinear ORPF method is proposed to acquire an enhanced time-delay difference estimation.

3.1. Time-Delay Model in the Narrow-Band Components

Assuming the received signal in the mth hydrophone sensor is ${\mathrm{sr}}_m\left(t\right)$, the sampling frequency and sampling interval of data collection are $F_s$ and T, respectively, and the total number of sampling points is $N=F_{s}T$. According to Equation (6), the spectrum of the received signal ${\mathrm{sr}}_m\left(t\right)$ at the mth sensor by using discrete Fourier transform (DFT) can be written by

$$\begin{array}{c}\hfill {\mathrm{sr}}_m\left(f_n\right)={\gamma }_m\left(f_n\right)s\left(f_n\right)exp(-j2\pi f_n{\tau }_m)+e_m\left(f_n\right),n\in (1,\cdots ,N),\end{array}$$

(9)

where $f_n=n{F}_s/N$ is the nth frequency grid, ${\gamma }_m\left(f_n\right)$ is the corresponding amplitude at the frequency of $f_n$, and $e_m\left(f_n\right)$ denotes a complex value at the frequency of $f_n$ after the DFT operator of the noise $e_m\left(t\right)$.

According to Equation (9), the wideband beamforming becomes multiple narrow-band beamforming through the DFT procedure. The weight vector with the predefined ${\varphi }_q$ in the nth frequency point is written by

$$\begin{array}{c}\hfill \mathbf{w}(f_n,{\varphi }_q)={\left[1,exp\left(\frac{-j2\pi f_{n}dcos{\varphi }_q}{c}\right)\cdots ,exp\left(\frac{-j2\pi f_n(M-1)dcos{\varphi }_q}{c}\right)\right]}^T,\end{array}$$

(10)

where ${\varphi }_q\in (0,180)$, $q\in 1,\cdots ,Q$, and Q is the number of the predefined bearing grid. Based on weights in Equation (10), the beamforming method is achieved by

$$\begin{array}{c}\hfill \rho \left({\varphi }_q\right)=\sum _{n=1}^N{\left|\mathbf{w}{(f_n,{\varphi }_q)}^H\mathrm{Sr}\left(f_n\right)\right|}^2,q\in 1,\cdots ,Q,\end{array}$$

(11)

where $\mathrm{Sr}\left(f_n\right)={[{\mathrm{sr}}_1\left(f_n\right),\cdots ,{\mathrm{sr}}_M\left(f_n\right)]}^T\in {\mathcal{C}}^M$ is a received signal vector in the nth frequency point. According to the result in Equation (11), the $\widehat{\varphi }$ is estimated based on the maximum-value detection criterion, and can be expressed by

$$\begin{array}{c}\hfill \widehat{\varphi }=\underset{\varphi }{argmax}\rho \left({\varphi }_q\right),q\in 1,\cdots ,Q.\end{array}$$

(12)

The output signal based on the estimated $\widehat{\varphi }$ is

$$\begin{array}{c}\hfill \widetilde{\mathrm{sr}}\left(n\right)=\sum _{m=1}^M{\mathrm{sr}}_m\left[n-\frac{(m-1)dcos\widehat{\varphi }}{c}\right]n\in (1,\cdots ,N).\end{array}$$

(13)

The operations above in Equations (9)–(13) are involved in the basic wideband frequency beamforming method. The weight vector in Equation (10) is optimal, and the signal-to-noise ratio (SNR) of the output signal $\widetilde{\mathrm{sr}}\left(n\right)={[{\mathrm{sr}}_1\left(n\right),\cdots ,{\mathrm{sr}}_N\left(n\right)]}^T$ can be improved by M times when the towed array is ideally uniform linear. However, the array should be deformed severely and not remain a straight line: due to the transverse motion of the towing vessel, oceanic currents, and hydrodynamic effects, weights in Equation (10) would be not optimal. Therefore, the performance of beamforming would worsen dramatically when the hydrophone array greatly deviates from the nominal one.

The output signal $\widetilde{\mathrm{sr}}\left(n\right)$ based on hypothetical array shape would not be the optimal one in the distorted array cases, because the time delays in the weight vector Equation (10) do not match the true one. Instead of directly performing the element-space narrow-band detection, we carry out the detection on the spectrum, which is acquired by performing the DFT operation on the output signal $\widetilde{\mathrm{sr}}\left(n\right)$ in Equation (13). By exploiting the background calibration and threshold-based detection techniques, the detected narrow-band frequencies radiated by the source of opportunity can be denoted as ${\widehat{f}}_l$ with $l\in 1,\cdots ,L_0$, where $L_0$ is the number of detected narrow-band frequencies. The corresponding complex value $g_{l,m}$ of the lth narrow-band frequency from the acoustic source to the mth sensor can be acquired by

$$\begin{array}{c}\hfill g_{l,m}\left({\widehat{f}}_l\right)={\gamma }_m\left({\widehat{f}}_l\right)exp\left(-j2\pi {\widehat{f}}_l{\tau }_m\right).\end{array}$$

(14)

In a similar manner, we have the complex value of ${\widehat{f}}_l$ in the $(m+1)$th hydrophone, and thus estimate the adjacent time-delay difference by

$$\begin{array}{c}\hfill \mathsf{\Delta }{g}_{l,m}\left({\widehat{f}}_l\right)=g_{l,m}\left({\widehat{f}}_l\right)g_{l,m+1}^{*}\left({\widehat{f}}_l\right)={\tilde{\gamma }}_{m}exp\left(j2\pi {\widehat{f}}_l\mathsf{\Delta }{\tau }_m\right),\end{array}$$

(15)

where ${\tilde{\gamma }}_m={\gamma }_m{\gamma }_{m+1}$ denotes the amplitude, and $\mathsf{\Delta }{\tau }_m={\tau }_{m+1}-{\tau }_m$ is the time-delay difference between the $(m+1)$th sensor and the mth sensor.

3.2. Proposed Outlier-Robust Particle Filter Method

Unlike the methods in [26,27,28], where only phase information of the detected narrow-band components is used to estimate time-delay differences, the complex narrow-band components including both amplitudes and phases would be exploited to acquire enhanced time-delay estimates in the paper. The benefits of the complex values lie in two aspects: on the one hand, although the information of time delays is involved only in the phases, the corresponding amplitudes can represent the credibility and stability of these phases, therefore, the joint exploitation of both the phase and amplitude information would provide a better way for the improvement of the time-delay estimation compared to those methods with only the use of the phase information. On the other hand, the time-delay estimate based on the complex values of narrow-band components would alleviate the phase ambiguity issue. In the following section, a novel outlier-robust particle filter algorithm is proposed to acquire the enhanced time-delay difference estimates by jointly exploiting the amplitude and phases of detected multiple narrow-band frequencies.

Instead of directly extracting their phases of multiple detected narrow-band components, we stack the complex value $\mathsf{\Delta }{g}_{1,m}$ into a vector as ${\mathbf{z}}_m={\left[\mathsf{\Delta }{g}_{1,m},\cdots ,\mathsf{\Delta }{g}_{L_0,m}\right]}^T$, and regard ${\mathbf{z}}_m$ as the observation data sample. In conformity with the analysis above, the observation equation is expressed by

$$\begin{array}{c}\hfill {\mathbf{z}}_m={\widehat{\gamma }}_{m}exp\left(-j2\pi \widehat{\mathbf{f}}\mathsf{\Delta }{\tau }_m\right)+{\mathbf{n}}_m,\end{array}$$

(16)

where ${\widehat{\gamma }}_m={\left[{\tilde{\gamma }}_{1,m},\cdots ,{\tilde{\gamma }}_{L_0,m}\right]}^T\in {\mathcal{R}}^{L_0\times 1}$ denotes an amplitude vector, $\widehat{\mathbf{f}}={\left[{\widehat{f}}_1,\cdots ,{\widehat{f}}_{L_0}\right]}^T\in {\mathcal{R}}^{L_0\times 1}$ represents a detected narrow-band component vector, and ${\mathbf{n}}_m\in {\mathcal{R}}^{L_0\times 1}$ is an additive Gaussian noise.

Considering the fact that the external shell of a hydrophone array is soft and elastic, it is reasonable to assume the time delays in the adjacent hydrophone sensors are slowly changing (${\left\{\mathsf{\Delta }{\tau }_m\right\}}_{m=0}^{M-1}$ should be identical, when the array is ideally straight), and thus we model the state equation by introducing a joint state vector ${\mathbf{\theta }}_m=\left[{\widehat{\gamma }}_m,\mathsf{\Delta }{\tau }_m\right]\in {\mathcal{R}}^{L_0+1}$ at the mth hydrophone array:

$$\begin{array}{c}\hfill {\mathbf{\theta }}_m={\mathbf{\theta }}_{m-1}+{\mathbf{u}}_m,m=1,\cdots N.\end{array}$$

(17)

Without generality of loss, a Gaussian process noise ${\mathbf{u}}_m\in {\mathcal{R}}^{L_0+1}$ is considered in the model.

According to Equations (16) and (17), the time-delay difference estimate problem can be regarded as a typical nonlinear discrete-time filtering problem, and can be expressed as

$$\begin{array}{cc}\hfill & {\theta }_{\mathbf{m}}=f\left({\theta }_{m-1}\right)+{\mathbf{u}}_m,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & {\mathbf{z}}_m=h\left({\theta }_m\right)+{\mathbf{n}}_m\hfill \end{array}$$

(19)

where ${\theta }_{\mathit{m}}\in {\mathit{R}}^{L_0+1}$ and ${\mathit{z}}_{\mathit{m}}\in {\mathit{R}}^{L_0}$ represent the state vector and observation vector at the mth hydrophone sensor, and ${\mathbf{u}}_m$ and ${\mathbf{n}}_m$ are zero-mean Gaussian white noise, satisfying ${\mathbf{u}}_m\sim \mathcal{N}\left({\mathbf{u}}_m|\mathbf{0},\mathit{Q}\right)$ and ${\mathbf{n}}_m\sim \mathcal{N}\left(\mathbf{0},\mathit{R}\right)$, respectively. $\mathit{R}$ is a diagonal matrix with diagonal elements $\mathbf{r}={[r_1,...,r_{L_0}]}^T$. $\mathit{Q}$ and $\mathit{R}$ are, respectively, the covariance matrix of process noise and observation noise. In our model, the function $f(·)$ is a linear function in Equation (17), and the function $h(·)$ is the nonlinear exponential function in Equation (16).

Following the outlier-robust Kalman filter method in [27,32,33,34], the weight vector is also introduced to alleviate the outlier. The probability distribution of observation noise is modified to

$$\begin{array}{c}\hfill {\mathbf{n}}_m\sim \mathcal{CN}\left({\mathbf{n}}_m|\mathbf{0},{\mathbf{W}}_m^{-1}\mathit{R}\right),\end{array}$$

(20)

where ${\mathbf{W}}_m\in {\mathit{R}}^{L_0\times L_0}$ is a diagonal array, and the weight vector ${\mathbf{w}}_m={[w_{1m},...,w_{L_{0}m}]}^T\in {\mathit{R}}^M$ is diagonal at the mth sensor. ${\mathbf{W}}_m^{-1}$ is the inverse of ${\mathbf{W}}_m$. Each element in the weight vector ${\mathbf{w}}_m$ obeys the Gamma distribution to ensure that each element is a positive value. Therefore, the prior probability of the model is

$$\begin{array}{c}\hfill {\mathbf{z}}_m|{\mathbf{\theta }}_m,{\mathbf{W}}_m\sim \mathcal{CN}(h\left({\mathbf{\theta }}_m\right),{\mathbf{W}}_m^{-1}\mathbf{R})\end{array}$$

(21)

$$\begin{array}{c}\hfill {\mathbf{\theta }}_m|{\mathbf{\theta }}_{m-1}\sim \mathcal{N}(f\left({\mathbf{\theta }}_{m-1}\right),\mathit{Q})\end{array}$$

(22)

$$\begin{array}{c}\hfill w_{lm}\sim \mathrm{Gamma}\left(w_{lm}|v_l/2,v_l/2\right),m=1,2,\dots ,M,l=1,\cdots ,L_0\end{array}$$

(23)

where $v_l$ is the hyperparameter of the Gamma distribution. By introducing the weight vector ${\mathbf{W}}_m$, the model has the ability to address the effect of the outliers [27]. Considering that the local linearization scheme, such as the extended Kalman filter, will bring about the loss of filtering accuracy, the expectation–maximum (EM) framework is used to perform the posterior inference. According to the expectation–maximum (EM) technique [35,36,37], it is necessary to carry out the posterior distribution of hidden variables in E-step, and to estimate the parameters in the system by maximizing the likelihood function in M-step. E-Step and M-Step are performed iteratively until the algorithm converges.

Considering that latent variables ${\mathbf{\theta }}_{1:M}$ and ${\mathit{w}}_{1:M}$ are independent of each other, the complete data has the following posterior distribution:

$$\begin{array}{c}\hfill p({\mathbf{\theta }}_{1:M},{\mathit{w}}_{1:M}|{\mathit{z}}_{1:M})=\frac{p\left({\mathit{z}}_{1:M}\right|{\mathbf{\theta }}_{1:M},{\mathit{w}}_{1:M})p\left({\mathbf{\theta }}_{1:M}\right)p\left({\mathit{w}}_{1:M}\right)}{p\left({\mathit{z}}_{1:M}\right)}\end{array}$$

(24)

Using the Markov property of variables, the prior distribution of states can be written as

$$\begin{array}{c}\hfill p({\mathbf{\theta }}_{1:M},{\mathit{w}}_{1:M}|{\mathit{z}}_{1:M})=\frac{p\left({\mathit{z}}_{1:M}\right|{\mathbf{\theta }}_{1:M},{\mathit{w}}_{1:M})p\left({\mathbf{\theta }}_{1:k}\right)p\left({\mathit{w}}_{1:M}\right)}{p\left({\mathit{z}}_{1:M}\right)}\end{array}$$

(25)

Using the conditional independence property, the joint distribution can be written as

$$\begin{array}{c}\hfill p({\mathbf{\theta }}_{1:M},{\mathit{w}}_{1:M}|{\mathit{z}}_{1:M})\propto p\left({\mathbf{\theta }}_1\right)\prod _{i=1}^{M}p\left({\mathit{z}}_i\right|{\mathbf{\theta }}_i,{\mathit{w}}_i)\prod _{i=1}^{M}p\left({\mathbf{\theta }}_i\right|{\mathbf{\theta }}_{i-1})\prod _{i=1}^{M}p\left({\mathit{w}}_i\right).\end{array}$$

(26)

Assuming that ${\mathbf{\theta }}_1\sim \mathcal{N}\left({\mathbf{\theta }}_1\right|{\mathbf{\mu }}_0,{\mathbf{V}}_0)$ with the initial mean and covariance of ${\mathbf{\mu }}_0$ and ${\mathbf{V}}_0$, the posterior distribution of hidden variables can be analytically acquired using the conjugate distribution property.

3.2.1. Posterior Probability Distribution Inference

The general form of the posterior estimation of the system state function $g\left({\mathbf{\theta }}_{\mathbf{m}}\right)$ can be expressed by

$$\begin{array}{cc}\hfill E\left[g\left({\mathbf{\theta }}_m\right)\right|{\mathit{z}}_{1:m},{\mathit{w}}_{1:m}]& =\int g\left({\mathbf{\theta }}_m\right)p\left(\mathbf{\theta }\right|{\mathit{z}}_{1:m},{\mathit{w}}_{1:m})d{\mathbf{\theta }}_m\hfill \end{array}$$

(27)

$$\begin{array}{c}=\int g\left({\mathbf{\theta }}_m\right)\frac{p\left({\mathbf{\theta }}_m\right|{\mathit{z}}_{1:m},{\mathit{w}}_{1:m})}{{\pi }_{{\mathbf{\theta }}_m}\left({\mathbf{\theta }}_m\right|z_{1:m})}{\pi }_{{\mathbf{\theta }}_m}\left({\mathbf{\theta }}_m\right|{\mathit{z}}_{1:m})d{\mathbf{\theta }}_m\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & =\frac{\int g\left({\mathbf{\theta }}_m\right)\frac{p({\mathbf{\theta }}_m,{\mathit{z}}_{1:m}|{\mathit{w}}_{1:m})}{{\pi }_{{\mathbf{\theta }}_m}\left({\mathbf{\theta }}_m\right|{\mathit{z}}_{1:m})}{\pi }_{{\mathbf{\theta }}_m}\left({\mathbf{\theta }}_m\right|\mathit{z}{}_{1:m})d{\mathbf{\theta }}_m}{\int \frac{p({\mathbf{\theta }}_m,{\mathit{z}}_{1:m}|{\mathit{w}}_{1:m})}{{\pi }_{{\mathbf{\theta }}_m}\left({\mathbf{\theta }}_m\right|{\mathit{z}}_{1:m})}{\pi }_{{\mathbf{\theta }}_m}\left({\mathbf{\theta }}_m\right|{\mathit{z}}_{1:m})d{\mathbf{\theta }}_m},\hfill \end{array}$$

(29)

where ${\pi }_{{\mathbf{\theta }}_m}\left({\mathbf{\theta }}_m\right|{\mathit{z}}_{1:m})$ represents the importance sampling of the system state ${\mathbf{\theta }}_m$. In our case, $g\left({\mathbf{\theta }}_m\right)={\mathbf{\theta }}_m$. According to the Monte Carlo sampling method, assuming that the particles draw from a nominal distribution ${\mathbf{\theta }}_m^{\left(s\right)}\sim {\pi }_{{\mathbf{\theta }}_m}\left({\mathbf{\theta }}_m\right|{\mathbf{z}}_{1:m}),s=1,2,...,S$, the weight of the particles can be written as

$$\begin{array}{c}\hfill w\left(\mathbf{\theta }\right)=\frac{p({\mathbf{\theta }}_m,{\mathbf{z}}_{1:m}|{\mathbf{w}}_{1:m})}{{\pi }_{{\mathbf{\theta }}_m}\left({\mathbf{\theta }}_m\right|{\mathbf{z}}_{1:m})}.\end{array}$$

(30)

The posterior expectation of estimated ${\mathbf{\theta }}_m$ can be written as

$$\begin{array}{c}\hfill {\widehat{\mathbf{\theta }}}_m=\frac{{\displaystyle \sum _{s=1}^S}{\mathbf{\theta }}_m^{\left(s\right)}w\left({\mathbf{\theta }}_m^{\left(s\right)}\right)}{{\displaystyle \sum _{s=1}^S}w\left({\mathbf{\theta }}_m^{\left(s\right)}\right)},\end{array}$$

(31)

where ${\widehat{\mathbf{\theta }}}_m$ is the estimated state vector, and the last element ${\widehat{\theta }}_{L_0+1,m}$ in the vector corresponds to the estimated time-delay difference $\mathsf{\Delta }{\widehat{\tau }}_m$.

According to the Bayes’ theorem, the recursive form of the particle weight $w\left({\mathbf{\theta }}_m^{\left(s\right)}\right)$ can be expressed by

$$\begin{array}{c}\hfill w\left({\mathbf{\theta }}_m^{\left(s\right)}\right)\propto w\left({\mathbf{\theta }}_{m-1}^{\left(s\right)}\right)p\left({\mathit{z}}_m|{\mathbf{\theta }}_m^{\left(s\right)}\right).\end{array}$$

(32)

Using the abovementioned particle update method and the recursive calculation method of particle weights, the posterior expectation calculation of the function of the system state can be realized.

For the learning of the weight vector ${\mathit{w}}_m$, the posterior obeys the Gamma distribution, therefore the posterior distribution of ${\mathit{w}}_m$ can be analytically obtained due to the conjugate distribution between Gaussian and Gamma distributions, and can be written by

$$\begin{array}{cc}\hfill p\left(w_m|z_{1:m},{\mathbf{\theta }}_{1:m},\right)& \propto p\left(z_m\right|{\mathbf{\theta }}_m,w_m)p\left(w_m\right)\hfill \\ \hfill & \propto {\left|{\mathit{W}}_m\right|}^{1/2}·exp\left\{-\frac{1}{2}{[{\mathit{z}}_m-h\left({\mathbf{\theta }}_m\right)]}^H{\mathit{R}}^{-1}{\mathit{W}}_m[z_m-h\left({\mathbf{\theta }}_m\right)]\right\}·\prod _{l=1}^{L_0}{w}_{lm}^{v_l/2-1}{e}^{-\frac{v_l{w}_{lm}}{2}}.\hfill \end{array}$$

(33)

Using ${\mathbf{\delta }}_m=z_m-h\left({\mathbf{\theta }}_m\right)$, the upper expression can be written as

$$\begin{array}{c}\hfill p\left({\mathit{w}}_m|{\mathit{z}}_{1:m},{\mathbf{\theta }}_{1:m}\right)\propto \prod _{l=1}^{L_0}{w}_{lm}^{(v_l/2+1/2)-1}exp\left[\left(-\frac{v_l}{2}-\frac{{\delta }_{lm}^2}{r_l}\right)w_{lm}\right],\end{array}$$

(34)

where ${\mathbf{\delta }}_m={[{\delta }_{1m},{\delta }_{2m},\dots ,{\delta }_{L_{0}m}]}^T$. It is observed that the posterior distribution of $w_{lm}$ still obeys Gamma distribution, so we can deduce the posterior distribution of $w_{lm}$, and have the following form:

$$\begin{array}{c}\hfill w_{lm}|z_{1:m},{\mathbf{\theta }}_{1:m}\sim \mathrm{Gamma}\left(w_{lm}\left|\frac{v_l}{2}+\frac{1}{2},\frac{v_l}{2}+\frac{{〈{\delta }_{lm}〉}^2}{2r_l}\right.\right).\end{array}$$

(35)

According to the conjugate distribution property, the posterior expectation of $w_{lm}$ can be provided by

$$\begin{array}{c}\hfill 〈w_{lm}〉=\frac{r_l{v}_l+r_l}{r_l{v}_l+{〈{\delta }_{lm}〉}^2},\end{array}$$

(36)

where $〈w_{lm}〉$ represents the expectation of a posterior distribution $w_{lm}$, and $〈{\delta }_{lm}〉$ represents the expectation of the lth element in the vector $[z_m-h\left({\mathbf{\theta }}_m\right)]$.

The above operation in the E-step stage is used to estimate the hidden variables of the system. Considering that the model parameter vector $\beta =\{{\mathbf{\mu }}_0,{\mathit{V}}_0,\mathit{Q},\mathit{R},\mathit{v}\}$ is often unknown in practical application, this learning process would be achieved in the M-step stage of the EM algorithm.

The learning of parameters in the EM algorithm is obtained by the maximum objective:

$$\begin{array}{c}\hfill {\beta }^{new}=\underset{\beta }{argmax}Q\left(\beta ,{\beta }^{old}\right),\end{array}$$

(37)

where $Q(\beta ,{\beta }^{old})$ represents the cost objective, and ${\beta }^{old}$ represents the value of the parameter set before learning, and ${\beta }^{new}$ denotes the updated parameter vector in the M-step stage. The cost objective function can be written by

$$\begin{array}{c}\hfill Q\left(\beta ,{\beta }^{old}\right)=\int p\left({\mathbf{\theta }}_{1:m},{\mathit{w}}_{1:m}|{\mathit{z}}_{1:m},{\beta }^{old}\right)lnp\left({\mathbf{\theta }}_{1:m},{\mathit{w}}_{1:m},{\mathit{z}}_{1:m}|\beta \right)d{\mathbf{\theta }}_{1:m}d{\mathit{w}}_{1:m}.\end{array}$$

(38)

To learn the parameter ${\mathbf{\mu }}_0$, we preserve these terms related to ${\mathbf{\mu }}_0$ in $Q(\beta ,{\beta }^{old})$, and express them as

$$\begin{array}{c}\hfill Q\left(\beta ,{\beta }^{old}\right)=〈-\frac{1}{2}{({\mathbf{\theta }}_1-{\mathbf{\mu }}_0)}^T{{\mathit{V}}_0}^{-1}({\mathbf{\theta }}_1-{\mathbf{\mu }}_0)〉+\mathrm{const},\end{array}$$

(39)

where the term of const represents all formulas unrelated to the estimated vector ${\mathbf{\mu }}_0$. The partial derivative of ${\mathbf{\mu }}_0$ is calculated by the above equation, and the partial derivative is equal to zero, and we obtain

$$\begin{array}{c}\hfill {\mu }_0^{new}=〈{\mathbf{\theta }}_1〉.\end{array}$$

(40)

To learn the parameter ${\mathit{V}}_0$, we re-express the $Q\left(\beta ,{\beta }^{old}\right)$, only keeping the terms related to ${\mathbf{\mu }}_0$, and denote

$$\begin{array}{c}\hfill Q\left(\beta ,{\beta }^{old}\right)=-\frac{1}{2}ln\left|{\mathit{V}}_0\right|-\frac{1}{2}〈{({\mathbf{\theta }}_1-{\mathbf{\mu }}_0)}^T{{\mathit{V}}_0}^{-1}({\mathbf{\theta }}_1-{\mathbf{\mu }}_0)〉+\mathrm{const}.\end{array}$$

(41)

The upper formula asks for the partial derivative of ${\mathit{V}}_0$ and makes the partial derivative equal to zero, and by taking the partial derivative of ${\mathit{V}}_0$ in the $Q\left(\beta ,{\beta }^{old}\right)$, we obtain

$$\begin{array}{c}\hfill \frac{\partial Q\left(\beta ,{\beta }^{old}\right)}{\partial {\mathit{V}}_0}=0\Rightarrow {\mathit{V}}_0^{new}=〈{\mathbf{\theta }}_1{{\mathbf{\theta }}_1}^T〉-〈{\mathbf{\theta }}_1〉{〈{\mathbf{\theta }}_1〉}^T.\end{array}$$

(42)

In a similar manner, we obtain the learning $\mathit{Q}$ by taking the the partial derivative of $\mathit{Q}$ in the cost objective function $Q\left(\beta ,{\beta }^{old}\right)$:

$$\begin{array}{c}\hfill {\mathit{Q}}^{new}=\frac{1}{m-1}\sum _{i=2}^m\left[〈{\mathbf{\theta }}_i{\mathbf{\theta }}_i^T〉-〈{\mathbf{\theta }}_i{\mathbf{\theta }}_{i-1}^T〉{\mathit{F}}_{i-1}^T-{\mathit{F}}_{i-1}〈{\mathbf{\theta }}_{i-1}{\mathbf{\theta }}_i^T〉+{\mathit{F}}_{i-1}〈{\mathbf{\theta }}_{i-1}{\mathbf{\theta }}_{i-1}^T〉{\mathit{F}}_{i-1}^T\right],\end{array}$$

(43)

where ${\mathit{F}}_{\mathbf{i}}$ represents the Jacobian matrix of the state function $f(·)$ at ${\mathbf{\theta }}_i$. Considering that the real state of the system represented by ${\mathbf{\theta }}_i$ is unknown, it can be replaced by $f\left(〈{\mathbf{\theta }}_{k-1}〉\right)$, which is expressed as

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathbf{i}}=\frac{\partial f\left(\mathbf{\theta }\right)}{\partial \mathbf{\theta }}{|}_{\mathbf{\theta }={\mathbf{\theta }}_i}.\end{array}$$

(44)

The upper derivation of $r_l^{-1}$ is obtained, the derivation is equal to zero, and we obtain

$$\begin{array}{c}\hfill r_l^{new}=\frac{1}{k}\sum _{i=1}^m〈w_{li}〉〈{\delta }_{li}^2〉.\end{array}$$

(45)

In order to learn the parameter vector $\mathit{v}$, we re-express the items related to $\mathit{v}$ in $Q\left(\beta ,{\beta }^{old}\right)$, and denote this as

$$\begin{array}{c}\hfill Q\left(\beta ,{\beta }^{old}\right)=\sum _{i=1}^m〈\frac{v_l}{2}ln\frac{v_l}{2}+\frac{v_l}{2}ln{w}_{li}-\frac{v_l}{2}{w}_{li}-ln\mathsf{\Gamma }\left(\frac{v_l}{2}\right)〉+\mathrm{const}.\end{array}$$

(46)

Finding the upper derivation of $v_l$, and letting the derivation be zero,

$$\begin{array}{c}\hfill ln\frac{v_l}{2}+2-\psi \left(\frac{v_l}{2}\right)=\frac{1}{m}\sum _{i=1}^m[〈w_{li}〉-\frac{v_l}{2}〈ln{w}_{li}〉].\end{array}$$

(47)

By solving the above equation, the estimated value $v_l^{new}$ can be obtained, and the parameter vector ${\mathbf{v}}^{new}={[v_1^{new},v_2^{new},\dots ,v_{L_0}^{new}]}^T$ can be obtained.

Once the state vector of ${\widehat{\mathbf{\theta }}}_m$ in Equation (31) is acquired, the estimated time-delay difference $\mathsf{\Delta }{\widehat{\tau }}_m$ simply corresponds to the last element in the vector. According to ${\left\{\mathsf{\Delta }{\widehat{\tau }}_m\right\}}_0^{M-1}$, we can calculate the locations of the distorted array shape. As the first two sensors define the the baseline, the estimated bearing $\widehat{\varphi }$ of the noise source of the opportunity can be calculated based on the estimated time-delay $\mathsf{\Delta }{\widehat{\tau }}_0$, according to Equation (3):

$$\begin{array}{c}\hfill \widehat{\varphi }=arccos\left(\frac{c\mathsf{\Delta }{\widehat{\tau }}_0}{d}\right).\end{array}$$

(48)

Therefore, the estimated deviation angle ${\widehat{\varphi }}_m$ can be calculated, as according to Equation (3),

$$\begin{array}{c}\hfill {\widehat{\phi }}_m=arccos\left(\frac{c\mathsf{\Delta }{\widehat{\tau }}_0}{d}\right)-arccos\left(\frac{c\mathsf{\Delta }{\widehat{\tau }}_m}{d}\right)\end{array}$$

(49)

The location $({\widehat{x}}_{m+1},{\widehat{y}}_{m+1})$ of $(m+1)$th sensor can be estimated by

$$\begin{array}{cc}\hfill {\widehat{x}}_{m+1}& ={\widehat{x}}_m+c\mathsf{\Delta }{\widehat{\tau }}_{m}cos\left({\widehat{\phi }}_m\right)+\sqrt{d^2-c^2\mathsf{\Delta }{\widehat{\tau }}_m^2}sin\left({\widehat{\phi }}_m\right)\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill {\widehat{y}}_{m+1}& ={\widehat{y}}_m+c\mathsf{\Delta }{\widehat{\tau }}_{m}sin\left({\widehat{\phi }}_m\right)+\sqrt{d^2-c^2\mathsf{\Delta }{\widehat{\tau }}_m^2}cos\left({\widehat{\phi }}_m\right),\hfill \end{array}$$

(51)

where $m=1,\cdots ,M-1$. Until now, we acquire the accurate estimates of the distorted linear hydrophone array.

3.2.2. Summary of the Proposed ORPF Method

According to the EM technique, the parameter estimates of the hidden variables ${\widehat{\mathbf{\theta }}}_m$ and ${\mathit{w}}_m$ are acquired in the E step, and other parameter learning is carried out by maximizing the objective function in the M step. The two steps are iteratively carried out to achieve the learning of model parameters while completing the system state estimation until the algorithm converges. The proposed distorted array shape estimation scheme using the ORPF algorithm is summarized in

Table 1. Proposed ORPF method.

Given the observed array data ${\left\{{\mathrm{sr}}_m\left(t\right)\right\}}_{m=1}^M$.

(1) Calculate $\widetilde{\mathrm{sr}}$ based on the hypothetical uniform linear array.

(2) Acquire the detected narrow-band frequencies ${\widehat{f}}_l$ with $l\in 1,\cdots ,L_0$ based on ${\mathrm{sr}}_m\left(t\right)$.

(3) Stacking observed measurement vector consisting of the complex values of each narrow-band component.

(4) Estimate the time-delay differences using proposed outlier-robust particle filter algorithm:

a. Initialize initial parameter set $\beta$.

b. Initialize covariance for observation noise $\mathbf{n}$, variance of state noise $\mathbf{u}$.

c. Estimate state vector ${\widehat{\theta }}_m$ and ${\widehat{\mathbf{w}}}_m$.

d. Estimate parameter set $\beta$.

(5) Perform array shape estimates $({\widehat{x}}_m,{\widehat{y}}_m)$.

(6) Perform the beamforming and acquire the improved output signal.

.

4. Simulations and Experiments

In this section, the proposed ORPF algorithm is verified and analyzed in the simulation and the real data. Three underwater acoustic signal enhancement methods, including subspace feature vector method [17], average delay difference estimation method [26], and outlier-robust Kalman filter (ORKF) method [27], are used for performance comparisons.

4.1. Simulations

In the simulations, the number of hydrophone sensors $M=128$ is used, and the interelement spacing between adjacent sensors d is $1.5$ m. The far-field underwater acoustic source of opportunity is located at ${60}^{\circ }$. The distorted towed hydrophone array is shown in the black curve in Figure 3d. The sampling frequency is $F_s=4000$ Hz. The sampling time of the simulated signal is 10 s. A three-parameter model with $f_c=500$ Hz, $f_m=200$ Hz and $\kappa =0$ is used to simulate the continuous spectrum. The five sinusoid waveforms with the frequencies of 32 Hz, 62 Hz, 96 Hz, 190 Hz, and 310 Hz are added in the radiated noise to model the narrow-band components, an additive Gaussian noise with the SNR of $-10$ dB is taken into consideration in the simulation, and additional noise is added in the narrow-band frequencies so that outliers occur.

Figure 3a shows the DOA estimation results based on the hypothetical uniform linear array. The DOA estimation result shows that the maximum bearing in the beam energy is ${54}^{\circ }$, which deviates from the bearing of the true one due to the distortion of the towed array. The narrow-band frequency detection is carried out, and five simulated narrow-band frequencies are detected. The time-delay differences in the detected narrow-band components are shown in Figure 3d. It is found that the time-delay differences are generally different, and show drastic fluctuations particularly in the frequencies of 32 Hz and 62 Hz, owing to the diverse SNRs on these narrow-band components. A number of outliers of time-delay differences appear, shown in Figure 3b. However, the time-delay differences should be identical among five narrow-band frequencies, because the true value of the time-delay difference is only determined by the bearing of the target source and the interelement spacing d according to Equation (2), and is independent of these detected narrow-band frequencies. The estimated time-delay differences are provided by Figure 3c. The ground truth of time-delay differences are shown in Figure 3c with the black line, and the estimated time-delay difference curves based on the average method and subspace-based method are, respectively, indicated by the blue and green lines. It is observed that the estimated time-delay differences drastically fluctuate across hydrophone sensors; the possible reason may be that both the subspace-based method and the average operation approach are sensitive to outliers. However, both ORKF and ORPF methods, which are robust to outliers by introducing the weight vector for the observation covariance matrix, acquire enhanced time-delay difference estimates. In addition, the proposed ORPF method has the lowest estimation error owing to the joint use of complex values, and the estimated time-delay difference curve is closest to the ground truth. Based on these estimated time-delay differences, the estimated distorted array shapes are shown in Figure 3d. The estimated array shape in the proposed method is closer to the real one, compared to those in the other methods.

The output spectra based on the estimated array shape are shown in Figure 3e. It is found that the output spectrum is closer to the spectrum in the known array shape, compared to those in other methods. To quantify the performance comparison, we introduce the gain loss, which is defined as the difference between the gains of the narrow-band components in the output signal with the estimated array shape and the gains with a known array shape, as the performance index, which is presented in

Table 2. Gain loss of detected narrow-band frequencies.

Frequency/Hz	32	62	96	190	310
Gain loss in the average method/dB	−20.3	−15.0	−13.2	−14.6	−27.9
Gain loss in the subspace-based method/dB	−20.9	−15.7	−13.7	−14.8	–27.8
Gain loss in the ORKF method/dB	−1.0	−1.7	−1.7	−1.8	−1.8
Gain loss in the proposed method/dB	−0.2	−0.2	−0.2	−0.3	−0.4

. It is clear that the proposed method is superior to other methods, the proposed method generally has the lowest gain losses among these four methods, and its average gain loss in these detected narrow-band frequencies is generally less than $-0.3$ dB, compared to the average losses of $-1.6$ dB in the ORKF method, $-18.6$ dB in the subspace-based method, and $-18.2$ dB in the average method.

The performance versus bearings of the noise source will be analyzed. The bearing of the underwater acoustic source of opportunity varies from ${10}^{\circ }$ to ${180}^{\circ }$ with the interval step of ${10}^{\circ }$, and the corresponding results are shown in Figure 4. It can be observed that the gain losses seem robust to the bearings of sources. However, the average gain losses of $-0.8$ dB in the proposed method are lowest, compared to $-1.3$ loss in the ORKF method, $-12.3$ dB loss in the subspace method, and $-13.1$ dB loss in the average method.

The performance comparisons versus SNR are also analyzed for the evaluation of the proposed method. The SNRs from $-15$ dB to $-3$ dB with a step of 2 dB are considered in this simulation. The results are shown in Figure 5. It is found that the gain losses generally decrease with the increase of SNR, and the ORKF and ORPF methods have lower gain losses than those in the average and subspace methods, due to the consideration of outlier robustness in the model. However, the proposed method has the lowest gain loss across the SNR, and the average gain loss is about $-0.5$ dB, compared to $-1.5$ dB in the ORKF method, $-4.8$ dB in the subspace-based method, and $-5.3$ dB in the average method.

4.2. Experiments

In this section, real data is used to verify the proposed algorithm. Figure 6 shows the real experiment scenarios on Qingjiang Lake, China. The target ship, as an acoustic source with a power amplifier, actively transmits low-frequency sinusoid waveforms with the frequencies of 83 Hz, 104 Hz, and 109 Hz, and a chirp signal whose initial frequency is 90 Hz and bandwidth 8 Hz are 15 m underwater. The other, ship-mounted by a towing array of 64 hydrophone sensors with a uniform interval of $1.5$ m, which is located 1 km away, is used to acquire the radiated noise data. There is a slight bow in the array. The data duration of T is 300 s, and the sample rate is 4 kHz.

The three narrow-band frequencies of 83 Hz, 104 Hz, and 109 Hz are obtained by using the narrow-band frequency detection technique, and the corresponding phases and amplitudes corresponding to these narrow-band components are shown in Figure 7a,b, respectively. It is observed that there indeed exist a number of outliers in both phases and amplitudes. Figure 7c shows the comparison of the time-delay difference estimation results of the proposed method with those of the average method, the subspace method, and the ORKF method. It can be found that the time-delay difference in the proposed ORPF method is relatively smooth, whereas they drastically fluctuate in both the average method and the subspace-based method. Although the ORKF algorithm can suppress the influence of some outliers, the fluctuation of the time-delay differences is still higher than those in the ORPF algorithm. The comparisons of the time-frequency spectrum obtained by the short-time Fourier transform of the 300-second output signal are provided in Figure 8a–e. The proposed ORPF method has significant enhancement in the time-frequency domain compared to the other three methods. To qualify the performance of the proposed method, we introduce the gain as an extra performance index, which is defined as the narrow-band frequency amplitude in the output signal over that in the conventional beamforming based on the hypothetical uniform linear array due to the unknown distorted array shape. The gains in our proposed method are generally more than 16 dB, even up to 25 dB, compared to the average gain of $12.1$ dB in the ORKF method, $5.8$ dB in the subspace-based method, and $5.0$ dB in the average method, as presented in

Table 3. Gains of detected narrow−band frequencies in real data.

Frequency/Hz	83	104	109
Gains in the average method/dB	5.3	5.0	4.8
Gains in the subspace-based method/dB	6.5	5.7	5.2
Gains in the ORKF method/dB	10.0	10.7	15.7
Gains in the proposed method/dB	16.1	17.2	25.2

.

5. Conclusions

Beamforming-based signal enhancement technologies in passive sonar array processing often have poor performances in the presence of array shape distortion. In this paper, an enhanced data-driven shape array estimation scheme is provided using passive underwater acoustic data, and a novel nonlinear outlier-robust particle filter method is proposed to acquire enhanced estimates of time delays in the presence of distorted hydrophone arrays. The proposed method takes full advantage of both amplitude and phase information of these narrow-band components for the improvement of the time-delay estimates; in addition, it also effectively alleviates the effect of phase ambiguities in the phase extraction in the only phase-based methods. The proposed scheme fully exploits directional radiated noise targets as sources of opportunity for online array shape estimation, and thus it requires neither the number nor direction of sources to be known in advance. Both simulations and experimental data show the superiority of the proposed method over other methods.