Open AccessArticle

An Enhanced Data-Driven Array Shape Estimation Method Using Passive Underwater Acoustic Data

Qisong Wu

Hao Zhang

Zhichao Lai

Youhai Xu

Shuai Yao

and

Jun Tao

Key Laboratory of Underwater Acoustic Signal Processing of Ministry of Education, Southeast University, Nanjing 210096, China

Author to whom correspondence should be addressed.

Remote Sens. 2021, 13(9), 1773; https://doi.org/10.3390/rs13091773

Submission received: 19 March 2021 / Revised: 18 April 2021 / Accepted: 26 April 2021 / Published: 2 May 2021

Download

Browse Figures

Graphical abstract
"> Figure 1
An example of a distorted hydrophone array. "> Figure 2
Spectrum example of real radiated noise. "> Figure 3
Time-delay differences varying with narrow-band frequencies in real data case. "> Figure 4
Performance comparisons. (a) Beamforming result based on the hypothetical array. (b) Estimation of array shape. (c) Comparison of the estimates of time-delay differences. (d) Time-delay differences in these detected narrow-band signals. (e) Spectrum of the improved output signal. "> Figure 5
Estimated time-delay difference variances. "> Figure 6
Amplitude error versus bearings of the radiated noise source. "> Figure 7
Average amplitude errors versus SINR in narrow-band signals. "> Figure 8
Performance comparisons. (a) Time-delay differences in 4 detected narrow-band signals. (b) Comparisons of the estimates of time-delay differences. "> Figure 9
Performance comparisons of Short Time Fourier Transform using 20 s hamming window with a window shift of 10 s. (a) Proposed method. (b) Average method. (c) Subspace method. (d) CBF method. "> Figure 10
Gains of the detected narrow-band frequencies. ">

Versions Notes

Abstract

Beamforming-based signal enhancement technologies in passive sonar array processing often have poor performance due to 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 formulation is proposed using passive underwater acoustic data. Beamforming based on a hypothetically ideal array is firstly employed to perform the detection of narrow-band components from sources of opportunity, and the corresponding phases of these detected narrow-band components are subsequently extracted to acquire time-delay differences. Then, a weighted outlier-robust Kalman smoother is proposed to acquire enhanced estimates of the time-delay differences, since the underlying properties of slowly changing time-delay differences in the hydrophone array and diverse signal to interference and noise ratios in multiple narrow-band components are explored; and its Cramer–Rao Lower Bound is also provided. Finally, the hydrophone array shape is estimated based on the estimated time delay differences. The proposed formulation fully exploits directional radiated noise signals from distant underwater acoustic targets as sources of opportunity for real-time array shape estimation, and thus it requires neither the number nor direction of sources to be known in advance. The effectiveness of the proposed method is validated in simulations and real experimental data.

Keywords:

array shape estimation; passive sonar system; radiated noise; time delay estimation; beamforming

Graphical Abstract

1. Introduction

Hydrophone arrays are widely used in underwater acoustics to perform direction-of-arrival (DOA) estimation and beamforming-based signal enhancement in the passive sonar system. Most array signal processing techniques are based on the assumption that the shape of the array remains unchanged. Unfortunately, the shape of a hydrophone array, such as a towed array, is often perturbed due to fluctuations in ship tactical maneuvering, oceanic currents, nonneutral buoyancy, and so on, particularly in the demand of the long-towed arrays for the acquisition of a low-frequency signal source at long range. These would lead to uncertainty in array shape and deviation from the hypothetical one [1,2,3,4]. However, beamforming as a typical signal enhancement technique requires knowledge of accurate array positions to properly process the measurements, which means an accurate estimate of the array shape is required to properly beamform the acoustic arrivals. The uncertainty in array shape inevitably degrades the beamforming performance [2]. A theoretical derivation of the performance loss in a distorted array shape due to the tactical maneuvering is provided [5]. Therefore, it is requisite to implement an array shape estimation to improve the acoustic system performance.

1.1. Related Works

A huge amount of effort has been devoted to array shape estimation. In general, these methods can be classified into source-independent methods and source-dependent methods. For source-independent methods, one effective method is to place depth-measurement sensors and compasses at several points within the array. These additional sensors provide localized respectively vertical and horizontal information on the transverse displacements of the array, whereas it is less economical due to the use of additional sensors. Further work is to model the array shape by a low order polynomial and to use the available sensors’ information to determine the coefficients of this polynomial [6]. A Kalman filter technique with a small angle approximation is used to estimate the array shape evolution over time [7], and an extended Kalman filter approach is also considered with the predicted state update following a higher-order nonlinear model [8].

As for source-dependent methods, acoustic hydrophone data would be used to estimate the array shape when there is no sensors’ information within the array or there are too few working sensors, and the shape can be estimated by measuring the relative time or phase delays among sensors due to one or more sources with known or unknown directions [9,10,11,12,13,14,15,16,17,18,19]. Historically, the array shape model employed for acoustic-based estimation is often considered to be a horizontal axis with perturbations at regular intervals on the vertical axis. Under this formulation, the Cramer–Rao Lower Bound (CRLB) for the narrow-band plane-wave sources is provided for source localization with a single known source direction or a known angle between the first receiver and the second one [9]. The maximum-likelihood (ML) approach is used to jointly estimate both the array and source parameters for multiple narrow-band or broadband sources with the known number of sources through a two-stage iteration procedure [10]. Subspace-based eigenvector techniques form another way to estimate array shapes when a single source at a known frequency and direction is known [11]. The eigenvector method for estimating the positions of the array elements is based on the eigen-decomposition of the array covariance matrix to acquire the phase delays between adjacent array elements. The statistic of phase delays is analyzed in [12]. The CRLB of the subspace-based array shape estimation method is provided in [3]. The eigenvector method is similar to the method [13], and it can be also regarded as a special case of the maximum-likelihood method [14]. Reference [15] proposes an improved array shape estimation method by jointly exploiting the acoustic hydrophone data and sensor data.

1.2. Contributions

In this paper, an enhanced data-driven array shape estimation formulation is proposed using passive underwater acoustic data. Instead of using the traditional method with a deterministic signal model and a known number of discrete sources, a novel data-driven approach is proposed by exploiting these directionals radiated noise sources, such as shipping and cargo as sources of opportunity for real-time array shape estimation. The conventional beamforming based on the hypothetical uniform linear array is firstly employed to acquire an improved waveform due to spatial gain from the hydrophone array. Since the radiated noise sources involve the amount of narrow-band components that originated from the vibration of mechanical components like diesel generator and air conditioning equipment, we perform the narrow-band frequency detection and extract their corresponding phases for the acquisition of noisy time-delay differences. A novel weighted outlier-robust Kalman smoother (WORKS) method is then proposed to acquire enhanced time-delay estimates by exploring the properties of slowly changing time-delay difference in hydrophone array as well as the diverse signal to interference and noise ratios (SINR) in multiple narrow-band components, and its CRLB is also analyzed and provided. Finally, the hydrophone array shape is obtained using the estimated time delay differences, and thus beamforming-based improved radiated noise signal is acquired. The proposed formulation fully exploits the directional radiated noise source signals as sources of opportunity for real-time array shape estimation, and thus it requires neither the number nor direction of sources to be known in advance. Furthermore, the proposed WORKS approach has the capability of acquiring the enhanced time-delay difference estimation by treating the weights associated with each observed noisy time-delay difference data in a probabilistic approach.

The contributions of the paper lie in three aspects: (1) An enhanced data-driven array shape estimation formulation using passive underwater acoustic data is provided. The proposed formulation fully exploits the directional radiated noise source signals as sources of opportunity for real-time array shape estimation, and thus it requires neither the number nor direction of sources to be known in advance. (2) A novel weighted outlier-robust Kalman smoother method is proposed to acquire enhanced time-delay estimates by the exploitation of slowly changing time-delay differences in hydrophone array and has the capability of alleviating the effects of the outliers resulted from the diverse noise levels in the multiple frequencies. (3) The CRLB of the proposed method is analyzed.

1.3. Organization and Notation

The paper is organized as follows: The underwater acoustic radiated noise signal model is introduced in Section 2. A novel array shape estimate technique is described in Section 3, and the weighted robust-outlier Kalman smoother method is also proposed in this section. Performance bound is also provided in Section 3. Simulations and experiments demonstrate the robustness of the proposed algorithms compared to conventional methods in Section 4. Concluding remarks are given in Section 5.

Notation: We use lower-case (upper-case) bold characters to denote vectors (matrices).

p (\cdot)

denotes the probability density function (pdf), and

N (x | a, b)

denotes that random variable x follows a Gaussian distribution with mean a and variance b. In addition,

{(\cdot)}^{T}

and

{(\cdot)}^{H}

denote transpose and conjugate transpose, respectively,

I_{N}

denotes the

N \times N

identity matrix.

\nabla_{x} = [\frac{\partial}{\partial x_{1}}, \dots, \frac{\partial}{\partial x_{N}}]

denotes the gradient operator respect to the variable vector

x = {[x_{1}, \dots, x_{N}]}^{T} \in R^{N \times 1}

, and

E (x)

represents the expectation operator of the random variable x.

2. Underwater Acoustic Radiated Noise Model in the Distorted Array

In this section, the underwater acoustic radiated noise model in a distorted hydrophone array is introduced. We assume that a thin flexible cable with M omnidirectional hydrophone sensors mounted at fixed inter-element spacing d is employed. Due to the transverse motion of ship tactical maneuvering, oceanic currents, and hydrodynamic effects, it would result in transverse displacements of the array from a straight line, as shown in Figure 1. The positions of the first two sensors define a baseline for the coordinate system and consequently are always assumed to reside on the x axis. It is observed that a deformed array occurs when any one of the inner sensor y positions are no longer zero.

The deformation is caused by bending at the sensor junctions due to the soft and elastic external shell of hydrophone array, which implies that the inter-sensor distances remain constant, and thus we have,

\begin{matrix} d = \sqrt{{(x_{m + 1} - x_{m})}^{2} + {(y_{m + 1} - y_{m})}^{2}}, m = 0, \dots, M - 1, \end{matrix}

(1)

where the location of the mth sensor is denoted by

(x_{m}, y_{m})

in the

x - y

Cartesian plane.

With respect to K far-field sources of opportunity impinging on the hydrophone array and the source with the bearing or direction of

θ_{k}

, the range difference between the first two sensors becomes, as shown in Figure 1,

\begin{matrix} Δ r_{k, 0} = r_{k, 1} - r_{k, 0} \approx d cos θ_{k} = (x_{1} - x_{0}) cos θ_{k}, \end{matrix}

(2)

where

r_{k, 1}

and

r_{k, 0}

denote the distances from the kth source to the second and reference sensors, respectively, and the result above is based on the far-field hypothesis. However, the locations of sensors deviate from the nominal ones, when the hydrophone array shape is distorted or deformed. Assuming the deviation angle formed by the line between the mth sensor and

(m + 1)

th sensor along x axis is denoted as

ϕ_{m}

, the range difference becomes,

\begin{matrix} Δ r_{k, m} \approx d cos (θ_{k} - ϕ_{m}) = (x_{m + 1} - x_{m}) cos θ_{k} + (y_{m + 1} - y_{m}) sin θ_{k} . \end{matrix}

(3)

For the first two sensors as a baseline,

ϕ_{0} = 0

. The range difference between the mth sensor and the reference can be written by,

\begin{matrix} r_{k, m} - r_{k, 0} = \sum_{i = 0}^{m - 1} Δ r_{k, i} = \sum_{i = 0}^{m - 1} d cos (θ_{k} - ϕ_{i}) = \sum_{i = 0}^{m - 1} (x_{i + 1} - x_{i}) cos θ_{k} + (y_{i + 1} - y_{i}) sin θ_{k} . \end{matrix}

(4)

The corresponding time-delay difference is

\begin{matrix} τ_{k, m} = \frac{\sum_{i = 0}^{m - 1} Δ r_{k, i}}{v} = \sum_{i = 0}^{m - 1} Δ τ_{k, i}, \end{matrix}

(5)

where v is the speed of acoustic wave in water, and

Δ τ_{k, i}

denotes the time-delay difference between the

(i + 1)

th and ith sensors with respect to the source of opportunity. It should be noted that the accuracies of the array shape estimates will be confined to the accuracies of each estimated time-delay difference for adjacent sensors, giving rise to error accumulations of the time-delay difference according to Equation (5).

Assuming the kth radiated noise source signal as

s_{k} (t)

with the bearing of

θ_{k}

, the received signal

y_{m} (t)

in the mth sensor can be written as

\begin{matrix} y_{m} (t) = \sum_{k = 1}^{K} ζ_{k, m} s_{k} (t - τ_{k, m}) + e_{m} (t), \end{matrix}

(6)

where

ζ_{k, m}

and

τ_{k, m}

, represent respectively the amplitude and time delay from the source to the mth sensor, and

e_{m} (t)

is an additive noise.

According to the statistical model of the underwater acoustic radiated noise [20,21], the corresponding time-domain

s_{k} (t)

is expressed by,

\begin{matrix} s_{k} (t) = p_{k l} (t) + (1 + p_{k m} (t)) p_{k c} (t), \end{matrix}

(7)

where

p_{k l} (t)

is a narrow-band signal arised from the vibration of machinery,

p_{k m} (t)

represents a modulation function due to the propeller movement, and

p_{k c} (t)

denotes a stationary hydrodynamic noise.

The narrow-band signal

p_{k l} (t)

, which is often generated from the vibration of mechanical components like diesel generator and air conditioning equipment, can be denoted by the superposition of multiple sinusoid signals as,

\begin{matrix} p_{k l} (t) = \sum_{l = 1}^{L} α_{k l} sin (2 π f_{k l} t), \end{matrix}

(8)

where

α_{k l}

is an amplitude of the lth single-tone frequency from the source, L is the number of the sinusoid signals, and

f_{k l}

is the corresponding frequency.

In the real sonar system, the narrow-band signals in the underwater acoustic radiation noise often have high source levels owing to inevitable mechanical vibration, and thus they are easier to be detected and recognized, compared to the broadband hydrodynamic noise signals [20,21]. Figure 2 shows a typical spectrum example of radiated noise from a ship. It is observed that there exist several narrow-band frequencies that have higher powers and are prominent in their nearby frequency ranges, their phases involve the information of time delays between the sources of opportunity and sensors in the deformed array. It possibly provides a way to estimate the distorted array shape.

3. Array Shape Estimation

In this section, a novel data-driven array shape estimation formulation will be provided, and the weighted robust-outlier Kalman smoother method is proposed to acquire an enhanced time-delay difference estimation.

3.1. Phase Extraction Based on Detected Narrow-Band Frequencies

Assume a passive sonar system with the sample rate of

F_{s}

and the time interval of T, and thus the total number of sample is

N = F_{s} T

. In accordance with Equation (6), the observed signal

y_{m} (t)

at the mth sensor is decomposed into N narrow-band components using Discrete Fourier Transform (DFT). The resulting narrow-band components corresponding to any one of the frequency points may be expressed as,

\begin{matrix} y_{m} (f_{n}) = \sum_{k = 1}^{K} ζ_{k, m} (f_{n}) s_{k} (f_{n}) exp (- j 2 π f_{n} τ_{k, m}) + e_{m} (f_{n}), n \in (1, \dots, N), \end{matrix}

(9)

where

f_{n} = n F_{s} / N

is the nth frequency grid,

ζ_{k, m} (f_{n})

is the coefficient at the

f_{n}

frequency with a high value, when a strong sinusoidal noise signal at the frequency of

f_{n}

exists, and K is the number of the far-field sources of opportunity.

According to Equation (9), the wideband beamforming is turned to be multiple narrow-band beamforming through the DFT procedure. The weight vector with the pre-defined

ϑ_{q}

in the nth frequency point is written by,

\begin{matrix} w (f_{n}, ϑ_{q}) = {[1, exp (\frac{- j 2 π f_{n} d cos ϑ_{q}}{v}) \dots, exp (\frac{- j 2 π f_{n} (M - 1) d cos ϑ_{q}}{v})]}^{T}, \end{matrix}

(10)

where

ϑ_{q} \in (0, 180)

q \in 1, \dots, Q

, and Q is the number of the pre-defined bearing grid. According to weights in Equation (10), beamforming is achieved by

\begin{matrix} ϱ (ϑ_{k}) = \sum_{n = 1}^{N} {|w {(f_{n}, ϑ_{k})}^{H} y (f_{n})|}^{2}, k \in 1, \dots, Q, \end{matrix}

(11)

where

y (f_{n}) = {[y_{1} (f_{n}), \dots, y_{M} (f_{n})]}^{T} \in C^{M}

is a received signal vector in the nth frequency point. According to the result in Equation (11), the

{\hat{ϑ}}_{k}

is estimated based on the maximum-value detection criterion. The output signal based on the estimated

{\hat{ϑ}}_{k}

is,

\begin{matrix} {\tilde{y}}_{k} (n) = \sum_{m = 1}^{M} y_{m} [n - \frac{(m - 1) d cos {\hat{ϑ}}_{k}}{v}], n \in (1, \dots, N), \end{matrix}

(12)

where M is the number of the hydrophone sensors.

The operations above in Equations (9)–(12) 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

{\tilde{y}}_{k} = {[{\tilde{y}}_{k} (1), \dots, {\tilde{y}}_{k} (N)]}^{T}

can be improved by M times, when the array is nominally uniform linear in the Gaussian noise case. 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 wrong. Therefore, the performance of beamforming would worsen dramatically, when the hydrophone array greatly deviates from the nominal one. Taking a uniform linear towed array of

M = 128

hydrophone sensors with a spatial interval of

1.6

m for example, the shape of the distorted array approximately becomes a parabola and deviation angle

ϕ_{m}

sin (π / (2 M))

, the average loss of the narrow-band signals would reach up to about 15 dB.

The output signal

{\tilde{y}}_{k}

based on the nominal 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 ones. Considering the fact that an important component in radiated noise originates from the machinery noise and can be regarded as the superposition of multiple narrowband signals, we would utilize the phases of these narrow-band signals to estimate time delays in the distorted array instead of the time delays in the nominal array. Instead of directly performing the element-space narrow-band detection, we carry out the detection on the spectrum of

P_{y k} (f)

, which is acquired by performing the DFT operation on the output signal

{\tilde{y}}_{k}

. By exploiting the background calibration and threshold-based detection techniques, the detected narrow-band frequencies radiated by the kth source of opportunity can be denoted as

{\hat{f}}_{k l}

with

l \in 1, \dots, L_{0}

, where

L_{0}

is the number of detected narrow-band frequencies. The corresponding unwrapped phase of the lth narrow-band frequency from the kth source to the mth sensor can be acquired by,

\begin{matrix} ψ_{m} ({\hat{f}}_{k l}) = 2 π {\hat{f}}_{k l} τ_{k, m} . \end{matrix}

(13)

In a similar way, we obtain the phase of

ψ_{m + 1} ({\hat{f}}_{k l})

in the

(m + 1)

th hydrophone, and thus estimate the adjacent phase difference by

\begin{matrix} Δ {\hat{ψ}}_{k, l, m} = ψ_{m + 1} ({\hat{f}}_{k l}) - ψ_{m} ({\hat{f}}_{k l}), l \in 1, \dots, L_{0} . \end{matrix}

(14)

\begin{matrix} Δ {\hat{τ}}_{k, l, m} = \frac{ψ_{m + 1} ({\hat{f}}_{k l}) - ψ_{m} ({\hat{f}}_{k l})}{2 π {\hat{f}}_{k l}}, l \in 1, \dots, L_{0} . \end{matrix}

(15)

3.2. Proposed Weighted Outlier-Robust Kalman Smoother

Ideally, the estimated time-delay difference

Δ {\hat{τ}}_{k, l, m}

between the

(m + 1)

th hydrophone sensor the mth hydrophone sensor expressed by Equation (14) should be independent of narrow-band frequencies, and be identical among

L_{0}

narrow-band frequencies. The reason is that its phase difference is only determined by the spatial interval d between adjacent hydrophones and the bearing of impinging source. However, the time-delay differences

Δ {\hat{τ}}_{k, l, m}

with

l \in 1, \dots, L_{0}

are often diverse in the real case. Figure 3 shows an example in a public SWellEx-96 HLA dataset [22]. It is observed that the time-delay differences in the detected narrow-band frequencies are completely different across sensors. The possible reason is that there are diverse SINRs in these detected narrow-band frequencies, and these time-delay differences will deviate from the true values at these frequencies.

One way to alleviate the fluctuations of

Δ {\hat{τ}}_{k, l, m}

is average these estimated time-delay differences and acquire [23]

\begin{matrix} Δ {\hat{τ}}_{k, m} = \frac{1}{L_{0}} \sum_{l = 1}^{L_{0}} Δ {\hat{τ}}_{k, l, m}, m \in 0, \dots, M - 2 . \end{matrix}

(16)

However, the average operator being a linear estimator is sensitive to outliers and easily deviates from the true value when there exists a small number of outliers. The

Δ {\hat{τ}}_{k, l, m}

acquired by the adjacent phase differences easily deteriorates in the low SINR cases, and thus become abnormal. In this case, these outliers would lead to a large deviation in the average operator. Therefore, the average operator is not an optimal way to acquire the time-delay estimation in the distorted array. To address this issue, we propose to use a weighted robust-outlier Kalman smoother.

Instead of averaging all the

Δ {\hat{ψ}}_{k, l, m}

in multiple narrow-band frequencies, we stack

Δ {\hat{ψ}}_{k, l, m}

with

l = 1, \dots, L_{0}

into a vector as

z_{k, m} = {[Δ {\hat{ψ}}_{k, 1, m}, \dots, Δ {\hat{ψ}}_{k, L_{0}, m}]}^{T}

, and regard

z_{k, m}

as the observation data sample. In conformity with the analysis above, the observation equation is expressed by,

\begin{matrix} z_{k, m} = α (f) \times Δ τ_{k, m} + ϵ_{m}, \end{matrix}

(17)

where

α (f) = {[2 π f_{1}, \dots, 2 π f_{L_{0}}]}^{T} \in R^{L_{0} \times 1}

is a vector,

Δ τ_{k, m}

is the true time-delay difference between the

(m + 1)

th and the mth sensors, and

ϵ_{m} \in R^{L_{0} \times 1}

is an additive noise.

According to Equation (17), the observed time-delay differences in multiple narrow-band frequencies can be considered as the true time-delay difference plus noise. To alleviate the effects of outliers and diverse noise levels in the multiple frequencies, we introduce a novel Bayesian algorithm that treats the weights associated with each observation data samples probabilistically as proposed by [24,25]. However, unlike [24,25] which use a scalar to weight the observed data sample and all elements in observed data sample are equally treated, we introduce a vector weight

γ_{m}

for each observed data sample

z_{k, m}

such that the variances of

z_{k, m}

are weighted with

γ_{m}

, owing to diverse SINRs in the corresponding narrow-band frequencies, and have,

\begin{matrix} z_{k, m} | γ_{m}, r \sim N (z_{k, m} | α (f) Δ τ_{k, m}, Γ_{m}^{- 1} R), \end{matrix}

(18)

where

R = diag (r)

, which is a covariance matrix for observation noise, is a diagonal matrix with

r = {[r_{1}, \dots, r_{L_{0}}]}^{T} \in R^{L_{0} \times 1}

on its diagonal.

Γ_{m} = diag (γ_{m})

is a diagonal matrix with

γ_{m} = {[γ_{1 m}, \dots, γ_{L_{0} m}]}^{T} \in R^{L_{0} \times 1}

on its diagonal. It should be noted that each elements in

r

is weighted by the corresponding

γ_{m}

. A Gamma prior distribution for the weights in order to ensure they remain positive is imposed on the weight vector

γ_{m}

as,

\begin{matrix} γ_{l m} \sim Gamma (γ_{l m} | a_{0}, b_{0}), l = 1, \dots, L_{0}, m = 0, \dots, M - 2, \end{matrix}

(19)

where Gamma(·) represents a Gamma distribution, and

a_{0}

and

b_{0}

are hyperparameters respectively, and often set to

a_{0} = b_{0}

Considering the fact that external shell of hydrophone array is soft and elastic, it is reasonable to assume the time-delay differences

Δ τ_{k, m}

are slowly changing across sensors (These

{Δ τ}_{k, m}

with

m = 0, \dots, M - 2

should be identical, when the array is ideally straight), and thus we have the state equation as

\begin{matrix} Δ τ_{k, m + 1} = Δ τ_{k, m} + ε_{m}, m = 0, \dots, M - 2, \end{matrix}

(20)

where

Δ τ_{k, m + 1}

is the time-delay difference between the

(m + 2)

th sensor and the

(m + 1)

th sensor, and

ε_{m}

is an additive noise. For the

m = 0

, the initialization of

Δ τ_{k, 0}

is assumed to follow the Gaussian distribution of

Δ τ_{k, 0} \sim N (Δ τ_{k, 0} | 0, ν_{0})

. Without loss of generality, a Gaussian prior is considered to model this noise, and can be expressed by,

\begin{matrix} ε_{m} | ν \sim N (ε_{m} | 0, ν), \end{matrix}

(21)

where

ν

is the variance of state noise.

In accordance with these time-delay differences in the M sensors, we have observed data samples

Δ τ_{k, 0 : M - 2}

at one time. To estimate the posterior distributions of the random variables and parameter values, we should consider the log evidence of the data samples observed, i.e.,

log p (z_{k, 0 : M - 2}, Δ τ_{k, 0 : M - 2}, γ_{0 : M - 2})

. It is a typical Kalman smoother.

We can treat this entire problem as an Expectation Minimization (EM) learning problem [26]. The expectation of the complete data likelihood should be taken with respect to the true posterior distribution of all hidden variables

q (Δ τ_{k, 0 : M - 2}, γ_{0 : M - 2})

. However, since this is an analytically intractable expression, we make a factorial approximation of the true posterior as follows [27],

\begin{matrix} q (Δ τ_{k, 0 : M - 2}, γ_{0 : M - 2}) = \prod_{m = 0}^{M - 2} q (γ_{m}) \prod_{m = 0}^{M - 3} q (Δ τ_{k, m + 1} | Δ τ_{k, m}) q (Δ τ_{k, 0}) . \end{matrix}

(22)

This factorization of

Δ τ_{k, 0 : M - 2}

considers the influence of each

Δ τ_{k, m}

within its Markov chain. According to this factor approximation, all resulting posterior distributions over hidden variables become analytically tractable. We can derive the final EM update equations from standard manipulations of Gaussian and Gamma distributions and acquire the following,

E - step :

Forward recursion,

\begin{matrix} N (Δ τ_{k, m} | μ_{m}, v_{m}) \propto N (z_{k, m} | α (f) Δ τ_{k, m}, Γ^{- 1} R) N (Δ τ_{k, m} | μ_{m - 1}, 〈ν〉 + v_{m - 1}), \end{matrix}

(23)

\begin{matrix} \to μ_{m} = v_{m} (\sum_{l = 1}^{L_{0}} \frac{〈γ_{l m}〉}{〈r_{l}〉} z_{k, m, l} + \frac{2 π f_{l} μ_{m - 1}}{〈ν〉 + v_{m - 1}}), \end{matrix}

(24)

\begin{matrix} \to v_{m} = {(\sum_{l = 1}^{L_{0}} \frac{〈γ_{l m}〉}{〈r_{l}〉} + \frac{1}{〈ν〉 + v_{m - 1}})}^{- 1}, \end{matrix}

(25)

\begin{matrix} \to 〈γ_{l m}〉 = \frac{a_{0} + 1 / 2}{\frac{1}{2 〈r_{l}〉} {(z_{k, m, l} - 2 π f_{l} μ_{m})}^{2} + b_{0}} l \in {1, \dots, L_{0}} . \end{matrix}

(26)

Backward recursion,

\begin{matrix} \hat{γ} (Δ τ_{k, m}) = N (Δ τ_{k, m} | {\hat{μ}}_{m}, {\hat{v}}_{m}) \propto N ({\hat{μ}}_{m + 1} | Δ τ_{k, m}, 〈ν〉) N (Δ τ_{k, m} | μ_{m}, v_{m}), \end{matrix}

(27)

\begin{matrix} \to {\hat{μ}}_{m} = {\hat{v}}_{m} (\frac{μ_{m}}{v_{m}} + \frac{{\hat{μ}}_{m + 1}}{〈ν〉}) = μ_{m} + \frac{v_{m}}{v_{m} + 〈ν〉} ({\hat{μ}}_{m + 1} - μ_{m}), \end{matrix}

(28)

\begin{matrix} \to {\hat{v}}_{m} = v_{m} + \frac{v_{m}^{2}}{{(v_{m} + 〈ν〉)}^{2}} [{\hat{v}}_{m + 1} - (v_{m} + 〈ν〉)] . \end{matrix}

(29)

M - step :

\begin{matrix} 〈r_{l}〉 & = \sum_{i = 1}^{m} \frac{〈γ_{l i}〉}{m} {(z_{k, i, l} - 〈Δ τ_{k, i}〉)}^{2}, \end{matrix}

(30)

\begin{matrix} 〈ν〉 & = \sum_{i = 1}^{m} \frac{1}{m} {(〈Δ τ_{k, i}〉 - 〈Δ τ_{k, i - 1}〉)}^{2} . \end{matrix}

(31)

Once these time-delay differences

{μ_{m}}_{m = 0}^{M - 1}

have been obtained, we can estimate the array shape. Since the first two sensors define the the baseline, the estimated bearing

{\hat{ϑ}}_{k}

of the kth noise source can be calculated based on the estimated time-delay

{\hat{μ}}_{0}

according to Equation (3),

\begin{matrix} {\hat{ϑ}}_{k} = arccos (\frac{v {\hat{μ}}_{1}}{d}) . \end{matrix}

(32)

Therefore, the estimated deviation angle

{\hat{ϕ}}_{m}

can be calculated as according to Equation (3),

\begin{matrix} {\hat{ϕ}}_{m} = arccos (\frac{v {\hat{μ}}_{1}}{d}) - arccos (\frac{v {\hat{μ}}_{m}}{d}) \end{matrix}

(33)

The location

({\hat{x}}_{m + 1}, {\hat{y}}_{m + 1})

(m + 1)

th sensor can be estimated by,

\begin{matrix} {\hat{x}}_{m + 1} & = {\hat{x}}_{m} + v {\hat{μ}}_{m} cos ({\hat{ϕ}}_{m}) + \sqrt{d^{2} - v^{2} {\hat{μ}}_{m}^{2}} sin ({\hat{ϕ}}_{m}) \end{matrix}

(34)

\begin{matrix} {\hat{y}}_{m + 1} & = {\hat{y}}_{m} + v {\hat{μ}}_{m} sin ({\hat{ϕ}}_{m}) + \sqrt{d^{2} - v^{2} {\hat{μ}}_{m}^{2}} cos ({\hat{ϕ}}_{m}), \end{matrix}

(35)

where

m = 2, \dots, M - 1

Finally, we perform the beamforming and acquire the improved output signal based on the estimated array shape as

\begin{matrix} {\bar{y}}_{k} (n) = y_{1} + \sum_{m = 1}^{M - 1} y_{m} [n - \sum_{i = 1}^{m} {\hat{μ}}_{i - 1}], \end{matrix}

(36)

where

n = [1, \dots, N]

3.3. Summary of Proposed Distorted Array Shape Estimation Scheme

The proposed distorted array shape estimation scheme is summarized in Table 1.

3.4. Performance Bound Analysis

To quantify the estimation performance of the proposed method, we perform the CRLB analysis in this section. Considering the problem of estimating the time-delay difference vector

Δ τ_{k} = {[Δ τ_{k, 0}, \dots, Δ τ_{k, M - 2}]}^{T}

from a set of observed data vectors of

z_{k} = {[z_{k, 0}^{T}, \dots, z_{k, M - 2}^{T}]}^{T}

with

z_{k, m} = {[Δ {\hat{ψ}}_{k, 1, m}, \dots, Δ {\hat{ψ}}_{k, L_{0}, m}]}^{T}

, we can acquire the joint probability distribution function (pdf)

p (Δ τ_{k}, z_{k})

from the kth source as based on Equations (18) and (20),

\begin{matrix} p (Δ τ_{k}, z_{k}) = \int_{γ_{m}} & \prod_{m = 0}^{M - 2} N (z_{k, m} | α (f) Δ τ_{k, m}, Γ_{m}^{- 1} R) N (Δ τ_{k, m + 1} | Δ τ_{k, m}, ν) \\ \times N (Δ τ_{k, 0} | 0, ν_{0}) Gamma (γ_{m} | a_{0}, b_{0}) d γ_{m} . \end{matrix}

(37)

The Fisher information matrix (FIM)

J (Δ τ_{k})

for the parameter vector

Δ τ_{k}

is defined as follows,

\begin{matrix} J (Δ τ_{k}) = - E \{\nabla_{Δ τ_{k}} {[\nabla_{Δ τ_{k}} ln p ({Δ τ}_{k}, z_{k})]}^{T}\} . \end{matrix}

(38)

When the derivatives and expectation exist, we acquire an alternative formula for the FIM as,

\begin{matrix} J (Δ τ_{k}) = E \{{[\nabla_{Δ τ_{k}} ln p ({Δ τ}_{k}, z_{k})]}^{T} \nabla_{Δ τ_{k}} ln p (Δ τ_{k}, z_{k})\} . \end{matrix}

(39)

Let

Δ {\hat{τ}}_{k} = Δ {\hat{τ}}_{k} (z_{k})

be an estimated of the parameter vector

Δ τ_{k}

. The mean-square error matrix (MSEM) is defined as

\begin{matrix} Ξ_{k} (Δ {\hat{τ}}_{k}) = E [(Δ τ_{k} - Δ {\hat{τ}}_{k}) {(Δ τ_{k} - Δ {\hat{τ}}_{k})}^{T}] . \end{matrix}

(40)

In fact, the MSEM is bounded by an inverse of the FIM, and we acquire [28],

\begin{matrix} Ξ_{k} (Δ {\hat{τ}}_{k}) \geq J^{- 1} (Δ τ_{k}) . \end{matrix}

(41)

The inverse of

J^{- 1} (Δ τ_{k})

is just the well-known posterior or Bayesian Cramer–Rao bound [29], and it will be denoted as

C (Δ τ_{k})

. According to the proposed generative model in Section 3.2, the logarithm of joint pdf of the state and measurement histories

p (z_{k}, Δ τ_{k})

can be expressed as,

\begin{matrix} ln p (z_{k}, Δ τ_{k}) = \sum_{m = 0}^{M - 2} ln p (z_{k, m} | Δ τ_{k, m}) + \sum_{m = 1}^{M - 2} ln p (Δ τ_{k, m} | Δ τ_{k, m - 1}) + ln p (Δ τ_{k, 0}) . \end{matrix}

(42)

Motivated by [28], we follow a similar recursive approach for the smoothing CRLBs, and express CRLBs for the filtering as in the forward recursion,

\begin{matrix} C_{m | m}^{- 1} & = C_{m | m - 1}^{- 1} + S_{m}^{m}, \end{matrix}

(43)

\begin{matrix} C_{m + 1 | m}^{- 1} & = G_{m + 1}^{m + 1} - G_{m + 1}^{m + 1, m} {(G_{m + 1}^{m} + C_{m | m}^{- 1})}^{- 1} G_{m + 1}^{m + 1, m}, \end{matrix}

(44)

\begin{matrix} C_{0 | - 1}^{- 1} & = G_{0}^{0}, \end{matrix}

(45)

where

m = 0, \dots, M - 2

. We acquire the CRLBs of the smoother in the backward recursion [28],

\begin{matrix} C_{m | M - 2}^{- 1} = C_{m | m}^{- 1} + G_{m + 1}^{m} - G_{m + 1}^{m + 1, m} {(G_{m + 1}^{m + 1} + C_{m + 1 | M - 2}^{- 1} + C_{m + 1 | m}^{- 1})}^{- 1} G_{m + 1}^{m + 1, m} . \end{matrix}

(46)

where,

\begin{matrix} G_{m + 1}^{m} & = E \{- \nabla_{Δ τ_{k, m}} [\nabla_{Δ τ_{k, m}} ln p (Δ τ_{k, m + 1} | Δ τ_{k, m})]\}, \end{matrix}

(47)

\begin{matrix} G_{m + 1}^{m, m + 1} & = E \{- \nabla_{Δ τ_{k, m + 1}} [\nabla_{Δ τ_{k, m}} ln p (Δ τ_{k, m + 1} | Δ τ_{k, m})]\}, \end{matrix}

(48)

\begin{matrix} G_{m + 1}^{m + 1} & = E \{- \nabla_{Δ τ_{k, m + 1}} [\nabla_{Δ τ_{k, m + 1}} ln p (Δ τ_{k, m + 1} | Δ τ_{k, m})]\}, \end{matrix}

(49)

\begin{matrix} S_{m}^{m} & = E \{- \nabla_{Δ τ_{k, m}} [\nabla_{Δ τ_{k, m}} ln p (z_{k, m} | Δ τ_{k, m})]\} . \end{matrix}

(50)

According to Equations (18) and (19) in the generative model,

p (z_{k, m})

would follow a multivariate Student distribution by integrating the weight vector

γ_{m}

, and can be expressed by,

\begin{matrix} z_{k, m} | Δ τ_{k, m}, \sim st (α (f) Δ τ_{k, m}, R, υ), \end{matrix}

(51)

where

st (\cdot)

denotes the student-t distribution, and the freedom

υ = 2 a_{0}

According to [28], we derive the recursive formula for the CRLB of the proposed WORKS method as in the forward recursion,

\begin{matrix} C_{m | m}^{- 1} & = C_{m | m - 1}^{- 1} + \frac{υ + L}{υ + L + 2} α {(f)}^{H} R^{- 1} α (f), \end{matrix}

(52)

\begin{matrix} C_{m + 1 | m}^{- 1} & = ν^{- 1} - ν^{- 1} {(ν^{- 1} + C_{m | m}^{- 1})}^{- 1} ν^{- 1}, \end{matrix}

(53)

\begin{matrix} C_{0 | - 1}^{- 1} & = ν_{k, 0} . \end{matrix}

(54)

and we acquire the corresponding recursive formula for the CRLB as in the backward recursion,

\begin{matrix} C_{m | M - 2}^{- 1} = C_{m | m}^{- 1} + ν^{- 1} - ν^{- 1} {(ν^{- 1} + C_{m + 1 | M - 2}^{- 1} + C_{m + 1 | m}^{- 1})}^{- 1} ν^{- 1}, \end{matrix}

(55)

where

m = 0, \dots, M - 2

4. Simulations and Experiments

In this section, both simulations and experiments will verify the effectiveness and correctness of the proposed method. Several typical source-dependent approaches, such as subspace-based eigenvector method [11] and average operation [23], are used for comparisons in simulations and experiments.

4.1. Simulations

Suppose that a uniform linear towing array has

M = 128

hydrophone sensors with a spatial interval of

1.6

m, and the bearing of a radiated noise source is

60^{\circ}

. The deviation angle

ϕ_{m}

is set to be

sin (π / (2 M))

, and thus the shape of the distorted array approximately becomes a parabola, as the black line shown in Figure 4a. The time duration is

T = 5

s and the sampling rate

F_{s} = 4

kHz. According to the characteristics of radiated noise [20,21], a three parameters model with

f_{c} = 500

Hz,

f_{m} = 200

Hz and

κ = 0

is used to model the continuous spectrum component. In addition, seven narrow-band (sinusoid) signals with frequencies of 59 Hz, 97 Hz, 125 Hz, 163 Hz, 198 Hz, 232 Hz and 280 Hz are considered. Without loss of generality, additive Gaussian noise and additional interferences in three narrow-band frequencies with 59 Hz, 97 Hz, and 163 Hz are considered with the SINR of

- 10

dB.

Figure 4b shows the beamforming result based on the hypothetical uniform linear array. The maximum value in the beam energy diagram is

55^{\circ}

, which is deviated from the true bearing of

60^{\circ}

, due to the time-delay mismatches in the distorted array. Only five narrow-band frequencies with 59 Hz, 97 Hz, 198 Hz, 163 Hz, and 232 Hz are acquired after the narrow-band detection based on the beamforming with the hypothetical array shape, and the narrow-band frequencies with 125 Hz and 280 Hz are missing. The corresponding time-delay differences of these detected narrow-band frequencies with 59 Hz, 97 Hz, 163 Hz, 198 Hz, and 232 Hz are shown in Figure 4d. It is observed that the time-delay differences are generally different, and show significant fluctuations particularly in the frequencies of 59 Hz, 97 Hz, and 163 Hz, due to the diverse SINRs on these narrow-band frequencies. Several outliers appear in the estimated time-delay differences. The time-delay differences should be independent of the considered narrow-band frequencies and be identical for the five narrow-band frequencies.

The ground truth of time-delay differences is shown in Figure 4c with the black line, while the average of these phase differences is shown in Figure 4c with the blue line. The estimation of time-delay differences in the subspace method is shown in Figure 4c with the green line. It is found that the estimates of time-delay differences in both the subspace-based method and the average operation method exhibit relatively large deviations from the ground truth and lead to quite similar estimates of the time-delay differences. The possible reason is that both the subspace-based method and the average operation approach are based on the narrow-band phase differences between adjacent elements, and both of them utilize the average operators for the estimates of time-delay differences in essence. The estimated time-delay differences in the proposed WORKS method, which significantly inhibits outliers, are quite closer to the true values, as shown in Figure 4c with the red line. According to the detected narrow-band frequencies, the array shape estimation is acquired by Equations (34) and (35), as shown in Figure 4a. The estimated array shape in the proposed method is closer to the real one, compared to those of the other two methods, even if gaps exist between them, especially in the latter part of the array. It is due to the error accumulations of time-delay differences in the array shape estimation, according to Equation (5).

The output spectrum based on the estimated array shape is provided in Figure 4e. It is reasonable that the output spectrum has higher gain, and is closer to the spectrum in the known array shape, compared to these in the subspace method and average operation method. It is also observed that the missing narrow-band frequency of 125 Hz becomes quite clear and can easily be detected in this enhanced spectrum. To quantify the superior performance, we use the amplitude error as the performance index. The amplitude error is defined as the absolute error between the amplitude of the detected narrow-band frequency and the amplitude in the corresponding frequencies in the spectrum based on the known array shape. The amplitude errors of the three methods are presented in Table 2. It is obvious that the proposed method is superior over the average operation and the subspace-based methods, owing to the ability to mitigate the effects of the outliers. It is also observed that the proposed method generally has the highest gain improvements among these three methods, and its amplitude errors in these detected narrow-band frequencies are generally less than

0.3

dB.

4.1.1. Variance Estimates of Time-Delay Difference

In this subsection, the variance estimates of time-delay differences are analyzed. All the parameters are the same as the simulations above. Figure 5 shows the variance results of the estimated

Δ τ

. It is obvious that the proposed approach has the lowest estimated variance across 128 sensors, and its average variance is 0.023, which is quite closer to the CRBL of 0.006, compared to 0.842 and 0.821 in the average operation and the subspace-based methods, respectively. The main reason is that the proposed method is quite robust to these outliers in the estimated time-delay differences possibly caused by the low-SNR phases, whereas both average operation and the subspace-based method are sensitive to these outliers, and would lead to relatively poor performance.

4.1.2. Performance Comparison versus Bearings of the Radiated Noise Source

The effectiveness of the proposed method has been verified thanks to the aforementioned simulations. In the following simulation, the performance versus bearings of the noise source will be analyzed. The source bearings change from

5^{\circ}

175^{\circ}

10^{\circ}

steps. Other parameters are the same as in the previous simulations. It is observed that all three approaches seem robust to the bearings of targets and their average gain errors of narrow-band frequencies in the end-fire bearing are comparable to those in the upright bearing. However, the average gain errors of results in the proposed method are much smaller than in the average operation and the subspace-based methods, as shown in Figure 6.

4.1.3. Performance Comparison versus SINR

In this simulation, we focus on the performance comparison versus diverse SINRs. The narrow-band frequency interferences of 59 Hz, 97 Hz, and 163 Hz with different amplitudes, which would cause diverse outliers in time-delay differences, are considered to form diverse SINRs from

- 3

dB to

- 14

dB with the step of

- 2

dB. Figure 7 shows average amplitude errors versus diverse SINRs. The proposed method has superior performance over the average operation and the subspace-based methods. It is observed that the average amplitude errors of the subspace-based method and the average operation approach gradually increase as the SINR decreases, whereas the proposed method is robust to different SINRs, which means the proposed method could mitigate the effects of different SINRs.

4.2. Experiments in Lake Trial Data

In this section, a set of lake trial data based on the distorted array is used to verify the effectiveness of the proposed method. In the experiment, the power amplifier as the acoustic source with 15 m in depth actively transmits the 83 Hz, 104 Hz, and 109 Hz low-frequency narrow-band signals and a chirp signal whose initial frequency is 90 Hz and bandwidth 8 Hz. The other ship equipped with a towing array of 64 hydrophone sensors with a uniform interval of

1.5

m, 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 second, and the sample rate is 4 kHz.

The three narrow-band frequencies of 83 Hz, 104 Hz, and 109 Hz are acquired by using the narrow-band frequency detection. The corresponding time-delay differences are shown in Figure 8a. It is observed that some perturbations of time-delay differences exist among these detected frequencies, and they would result in large derivations from each other and lead to the appearances of outliers. The possible reason is that the SINRs in each detected narrow-band frequencies are quite different. The time-delay differences estimation in the average method, the subspace-based method, and the proposed WORKS method are shown in Figure 8b, respectively. The time-delay differences in the proposed method are quite smooth and slowly change, compared to thees in the average method and the subspace method. The output time-frequency spectrum based on the enhanced signals during 300 s are given in Figure 9a–d. It is found that the spectrum of enhanced signals based on the proposed method has higher amplitudes in both narrow-band signals and chirp signals, compared to these in the other two methods. The proposed method is superior to the average operation and the subspace-based methods. The quantitative performance concerning the average gains of three narrow-band frequencies are provided as shown in Figure 10. Note that the gains are calculated by the narrow-band frequency amplitude in the output signal over that in the conventional beamforming assuming a hypothetical uniform linear array due to the unknown distorted array shape. It is found the gains in the our proposed method are generally above 15 dB, even up to 30 dB, compared to about 5 dB gains in both the average operation and the subspace-based methods.

5. Conclusions

A robust and data-driven signal enhancement scheme is proposed to take into account the array distortion in this paper. The proposed scheme firstly extracts the phases of narrow-band signals from the conventional beamforming result using narrow-band signal detection. Then, we proposed a novel weighted outlier-robust Kalman smoother to estimate time-delay differences to mitigate outliers in the time-delay differences. Finally, an enhanced output signal is achieved by beamforming based on the estimated time-delay differences. The proposed scheme is validated by both simulations and experiments. The proposed method is a data-driven approach that fully exploits the directional radiated noise signal due to distant underwater acoustic targets as sources of opportunity for real-time array shape estimation. It requires neither the number nor direction of sources to be known in advance.

Author Contributions

Conceptualization, Q.W. and H.Z.; methodology, Q.W., H.Z. and Z.L.; software, H.Z. and Y.X.; validation, H.Z., Z.L. and S.Y.; formal analysis, Q.W. and H.Z.; investigation, Q.W., Y.X. and S.Y. writing—original draft preparation, Q.W. and H.Z.; writing—Q.W. and J.T.; supervision, Q.W., Y.X., S.Y. and J.T.; funding acquisition, Q.W. All authors have read and agreed to the published version of the manuscript.

Funding

The work was supported in part National Natural Science Foundation under Grants No. 61701109 and 91938203, in part by Science and Technology on Sonar Laboratory under Grant No. 6142109180202, and in part by National Defense Basis Scientific Research program of China under Grant No. JCKY2019110C143, in part supported by National Natural Science Foundation of Jiangsu Province under Grant No BK20160701.

Data Availability Statement

The data presented in this paper are available after contacting the corresponding author.

Acknowledgments

We thank anonymous reviewers for their comments towards improving this manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

Schurman, I. Reverberation rejection with a dual-line towed array. IEEE J. Ocean. Eng. 1996, 21, 193–204. [Google Scholar] [CrossRef]
Felisberto, P.; Jesus, S.M. Towed-array beamforming during ship’s maneuvering. Proc. Inst. Elect. Eng. Radar Sonar Navig. 1996, 143, 210–215. [Google Scholar] [CrossRef]
Park, H.Y.; Lee, C.; Kang, H.G.; Youn, D.H. Generalization of the subspace-based array shape estimations. IEEE J. Ocean. Eng. 2004, 29, 847–856. [Google Scholar] [CrossRef]
Hodgkiss, W. The effects of array shape perturbation on beamforming and passive ranging. IEEE J. Ocean. Eng. 1983, 8, 120–130. [Google Scholar] [CrossRef]
Groen, J.; Beerens, S.P.; Doisy, Y. Shape and Doppler correction for a towed sonar array. In Proceedings of the Oceans’04 MTS/IEEE Techno-Ocean ’04, Kobe, Japan, 9–12 November 2004; Volume 2, pp. 613–620. [Google Scholar]
Owsely, N.L. Shape estimation for a flexible underwater cable. IEEE EASCON 1981, 14, 16–19. [Google Scholar]
Gray, D.; Anderson, B.; Bitmead, R. Towed array shape estimation using Kalman filters-theoretical models. IEEE J. Ocean. Eng. 1993, 18, 543–556. [Google Scholar] [CrossRef]
Newhall, B.; Jenkins, J.; Dietz, J. Improved estimation of the shape of towed sonar arrays. In Proceedings of the IEEE Instrumentation and Measurement Technology Conference, Como, Italy, 18–20 May 2004; pp. 873–876. [Google Scholar]
Rockah, Y.; Schultheiss, P.M. Array shape calibration using sources in unknown locations—Part I: Far-field sources. IEEE Trans. Acoust. Speech Signal Process. 1987, 35, 286–299. [Google Scholar] [CrossRef] [Green Version]
Weiss, A.J.; Friedlander, B. Array shape calibration using sources in unknown locations-a maximum likelihood approach. IEEE Trans. Acoust. Speech Signal Process. 1989, 37, 1958–1966. [Google Scholar] [CrossRef]
Gray, D.A.; Wolf, W.O.; Riley, J.L. An eigenvector method for estimating the positions of the elements of an array of receivers. In Proceedings of the Conference Signal Processing Applic, Adelaide, Australia, 17–19 April 1989; pp. 391–393. [Google Scholar]
Smith, J.J.; Leung, Y.H.; Cantoni, A. Sptatistics of the phase delays between array receivers estimated from the eigendecomposition of the signal correlation matrix. In Proceedings of the Conference IEEE International Conference on Acoustics, Speech and Signal Processing, Adelaide, SA, Australia, 19–22 April 1994; pp. 325–329. [Google Scholar]
Bouvet, M. Beamforming of a distorted line array in the presence of uncertainties on the sensor positions. J. Acoust. Soc. Am. 1987, 81, 1833–1840. [Google Scholar] [CrossRef]
Owsley, N.L.; Swope, G.R. Array shape determination using time delay estimation procedures. IEEE EASCON 1980, 13, 158–164. [Google Scholar]
Odom, J.L.; Krolik, J.L. Passive towed array shape estimation using heading and acoustic data. IEEE J. Ocean. Eng. 2015, 40, 465–474. [Google Scholar] [CrossRef]
Ferguson, B.G. Sharpness applied to the adaptive beamforming of acoustic data from a towed array of unknown shape. J. Acoust. Soc. Am. 1990, 88, 2695–2701. [Google Scholar] [CrossRef]
Quinn, B.G.; Barrett, R.F.; Kootsookos, P.J.; Searle, S.J. The estimation of the shape of an array using a hidden Markov model. IEEE J. Ocean. Eng. 1993, 18, 557–564. [Google Scholar] [CrossRef] [Green Version]
Wahl, D.E. Towed array shape estimation using frequency wavenumber data. IEEE J. Ocean. Eng. 1993, 18, 582–590. [Google Scholar] [CrossRef]
Viberg, M.; Swindlehurst, A.L. A Bayesian approach to autocalibration for parametric array signal processing. IEEE Trans. Signal Process. 1994, 42, 3495–3507. [Google Scholar] [CrossRef] [Green Version]
Ross, D.; Kuperman, W.A. Mechanics of underwater noise. J. Acoust. Soc. Am. 1989, 86, 1626. [Google Scholar] [CrossRef]
Wales, S.C.; Heitmeyer, R.M. An ensemble source spectra model for merchant ship-radiated noise. J. Acoust. Soc. Am. 2002, 111, 1211–1231. [Google Scholar] [CrossRef] [PubMed]
Information on the Swellex96 Experiment Data. Available online: http://swellex96.ucsd.edu/ (accessed on 15 July 2020).
Wu, Q.; Xu, P.; Li, T.; Fang, S. Feature enhancement technique with distorted towed array in the underwater radiated noise. In Proceedings of the INTER-NOISE and NOISE-CON Congress and Conference Proceedings. Institute of Noise Control Engineering, Hong Kong, China, 27–30 August 2017; Volume 255, pp. 3824–3830. [Google Scholar]
Ting, J.; Evangelos, T.; Stefan, S. A Kalman filter for robust outlier detection. In Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, San Diego, CA, USA, 29 October–2 November 2007; pp. 1514–1519. [Google Scholar]
Ting, J.; Evangelos, T.; Stefan, S. Learning an outlier-robust Kalman filter. In Proceedings of the European Conference on Machine Learning, Warsaw, Poland, 17–21 September 2007; pp. 748–756. [Google Scholar]
Dempster, A.; Laird, N.; Rubin, D. Maximum likelihood from incomplete data via the EM algorithm. J. R. Stat. Soc. Ser. B 1977, 39, 1–38. [Google Scholar]
Bishop, C.M. Pattern Recognition and Machine Learning; Springer: New York, NY, USA, 2006; pp. 423–517. [Google Scholar]
Šimandl, M.; Královec, J.; Tichavsky, P. Filtering, predictive, and smoothing Cramer–Rao bounds for discrete-time nonlinear dynamic systems. Automatica 2001, 37, 1703–1716. [Google Scholar] [CrossRef]
Tichavsky, P.; Muravchik, C.H.; Nehorai, A. Posterior Cramer–Rao bounds for discrete-time nonlinear filtering. IEEE Trans. Signal Process. 1998, 46, 1386–1396. [Google Scholar] [CrossRef] [Green Version]

Figure 1. An example of a distorted hydrophone array.

Figure 2. Spectrum example of real radiated noise.

Figure 3. Time-delay differences varying with narrow-band frequencies in real data case.

Figure 4. Performance comparisons. (a) Beamforming result based on the hypothetical array. (b) Estimation of array shape. (c) Comparison of the estimates of time-delay differences. (d) Time-delay differences in these detected narrow-band signals. (e) Spectrum of the improved output signal.

Figure 5. Estimated time-delay difference variances.

Figure 6. Amplitude error versus bearings of the radiated noise source.

Figure 7. Average amplitude errors versus SINR in narrow-band signals.

Figure 8. Performance comparisons. (a) Time-delay differences in 4 detected narrow-band signals. (b) Comparisons of the estimates of time-delay differences.

Figure 9. Performance comparisons of Short Time Fourier Transform using 20 s hamming window with a window shift of 10 s. (a) Proposed method. (b) Average method. (c) Subspace method. (d) CBF method.

Figure 10. Gains of the detected narrow-band frequencies.

Table 1. Summary of Proposed Distorted Array Shape Estimation Scheme.

Given the observed array data

{y_{m} (t)}_{m = 1}^{M}

(1) Calculate

{\tilde{y}}_{k}

with respect to the kth source based on the hypothetical uniform linear array.

(2) Acquire the detected narrow-band frequencies

{\hat{f}}_{k l}

with

l \in 1, \dots, L_{0}

based on

y_{m} (t)

(3) Calculate time-delay differences vector

z_{k, m} = {[Δ {\hat{ψ}}_{k, 1, m}, \dots, Δ {\hat{ψ}}_{k, L_{0}, m}]}^{T}, m = 0, \dots, M - 2

(4) Estimate the time-delay differences using proposed Weighted Robust-Outlier Kalman Smoother:

a. Initialize hyperparameters

a_{0}

and

b_{0}

b. Initialize covariance for observation noise

R

, variance of state noise

ν

c. Estimate time-delay differences

{\hat{μ}}_{m}

using EM algorithm Equations (23)–(30).

(5) Perform array shape estimates

({\hat{x}}_{m}, {\hat{y}}_{m})

with

m = 2, \dots, M - 1

using Equations (34) and (35).

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

{\bar{y}}_{k} (t)

Table 2. Amplitude errors of detected narrow-band frequencies.

Frequency/Hz	59	97	125	163	198	232	280
Amplitude error in average method/dB	25.3	15.0	43.2	24.6	37.9	11.4	25.2
Amplitude error in subspace method/dB	24.9	15.7	49.7	24.8	32.8	11.3	22.6
Amplitude error in proposed method/dB	0.1	0.1	0.1	0.1	0.2	0.3	0.3

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Wu, Q.; Zhang, H.; Lai, Z.; Xu, Y.; Yao, S.; Tao, J. An Enhanced Data-Driven Array Shape Estimation Method Using Passive Underwater Acoustic Data. Remote Sens. 2021, 13, 1773. https://doi.org/10.3390/rs13091773

AMA Style

Wu Q, Zhang H, Lai Z, Xu Y, Yao S, Tao J. An Enhanced Data-Driven Array Shape Estimation Method Using Passive Underwater Acoustic Data. Remote Sensing. 2021; 13(9):1773. https://doi.org/10.3390/rs13091773

Chicago/Turabian Style

Wu, Qisong, Hao Zhang, Zhichao Lai, Youhai Xu, Shuai Yao, and Jun Tao. 2021. "An Enhanced Data-Driven Array Shape Estimation Method Using Passive Underwater Acoustic Data" Remote Sensing 13, no. 9: 1773. https://doi.org/10.3390/rs13091773

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

An Enhanced Data-Driven Array Shape Estimation Method Using Passive Underwater Acoustic Data

Abstract

1. Introduction

1.1. Related Works

1.2. Contributions

1.3. Organization and Notation

2. Underwater Acoustic Radiated Noise Model in the Distorted Array

3. Array Shape Estimation

3.1. Phase Extraction Based on Detected Narrow-Band Frequencies

3.2. Proposed Weighted Outlier-Robust Kalman Smoother

3.3. Summary of Proposed Distorted Array Shape Estimation Scheme

3.4. Performance Bound Analysis

4. Simulations and Experiments

4.1. Simulations

4.1.1. Variance Estimates of Time-Delay Difference

4.1.2. Performance Comparison versus Bearings of the Radiated Noise Source

4.1.3. Performance Comparison versus SINR

4.2. Experiments in Lake Trial Data

5. Conclusions

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI