A Time-Frequency Joint Time-Delay Difference Estimation Method for Signal Enhancement in the Distorted towed Hydrophone Array

Abstract

To acquire the enhanced underwater ship-radiated noise signal in the presence of array shape distortion in a passive sonar system, the phase difference of the line-spectrum component in ship-radiated noise is often exploited to estimate the time-delay difference for the beamforming-based signal enhancement. However, the time-delay difference estimation performance drastically degrades with decreases of the signal-to-noise ratio (SNR) of the line-spectrum component. Meanwhile, although the time-delay difference estimation performance of the high-frequency line-spectrum components is generally superior to that of the low-frequency one, the phase difference measurements of the high-frequency line-spectrum component often easily encounter the issue of modulus 2$\pi$ ambiguity. To address the above issues, a novel time-frequency joint time-delay difference estimation method is proposed in this paper. The proposed method establishes a data-driven hidden Markov model with robustness to phase difference ambiguity by fully exploiting the underlying property of slowly changing the time-delay difference over time. Thus, the phase difference measurements available for time-delay difference estimation are extended from that of low-frequency line-spectrum components in a single frame to that of all detected line-spectrum components in multiple frames. By jointly taking advantage of the phase difference measurements in both time and frequency dimensions, the proposed method can acquire enhanced time-delay difference estimates even in a low SNR case. Both simulation and at-sea experimental results have demonstrated the effectiveness of the proposed method.

1. Introduction

Beamforming-based signal enhancement is a crucial problem in array signal processing, and plays a significant role in feature extraction and target recognition [1,2,3,4]. In a passive sonar system, an important topic is to acquire the enhanced underwater ship-radiated noise signal from the received data of the hydrophone array [5,6,7]. A large aperture is often required to achieve accurate localization and high array gain [8,9]. The large aperture is typically formed by trailing a hydrophone array behind a towing platform in a nominally straight line [10,11]. However, the array is often deformed or distorted due to inevitable oceanic currents, hydrodynamics, and tactical maneuvers of the towing platform, resulting in time-delay mismatch in beamforming-based signal enhancement, which seriously degrades the signal enhancement performance [12,13,14,15].

Over the past several decades, a tremendous amount of effort has been devoted to time-delay difference estimation in the distorted towed hydrophone array. One intuitive method is to install the compasses and depth sensors at several points within the towed array, providing localized horizontal and vertical information on the transverse displacements of the array, respectively [11,16]. Although this type of method can acquire the array shape straightforwardly, the limited accuracy and information update rate of these auxiliary sensors make it difficult to accurately estimate the array shape in real-time [17,18]. The generalized cross-correlation (GCC) estimator, consisting of a pair of prefilters and a cross-correlator, determines the time-delay difference by locating the peak of the cross-correlator output [19,20]. Note that the towed array with a large aperture mainly focuses on the weak targets. However, the correlation between the wideband components of ship-radiated noise signals received by different hydrophones decreases significantly in the low signal-to-noise ratio (SNR) scenario, drastically degrading the time-delay difference estimation performance of the GCC method [21,22].

Line-spectrum components generated from the inevitable vibration of mechanical equipment, such as the diesel generator and air conditioning system, are a quite important and useful element in the underwater ship-radiated noise signal [23,24]. Generally, the power of the line-spectrum components is at least 10 dB higher than that of their nearby continuous spectrum such that the line-spectrum components are ready to be detected and recognized [25]. The phases of these relatively strong line-spectrum components involve the information of time-delays from the target to hydrophones, and thus can be exploited to estimate the time-delay difference of the radiated noise signal received by different hydrophones. The methods developed in [26,27] detect the line-spectrum components by locating the peaks of the power spectrum of the received acoustic signal. An adaptive discrete-time nonlinear frequency locked-loop filter is properly designed to extract the line-spectrum components from the acquired acoustic signal [28,29]. The method in [6] first modeled the line-spectrum components as the superposition of sinusoidal signals and narrow-band noise, then estimated the time-delay difference of each line-spectrum component utilizing the phase difference measurement, and finally averaged the results to acquire the final time-delay difference estimates. The minimum time-delay difference estimation variance using a sinusoidal signal with unknown frequency has been derived in [30]. The method in [31] regarded the time-delay difference as the slope of the regression line of the phase difference and applied regression analysis to acquire the time-delay differences estimates. Note that first, the accuracy of time-delay difference estimation exploiting phase difference of line-spectrum component is sensitive to noise, and drastically degrades with decreases of the SNR of the line-spectrum component [3]. Furthermore, regarding the line-spectrum components with the same SNR, the time-delay difference estimation accuracy is proportional to the frequency of the line-spectrum component, i.e., the estimation accuracy of the high-frequency line-spectrum component outperforms that of the low-frequency one [30]. However, for the phase difference measurements of high-frequency line-spectrum component whose half-wavelength is less than the inter-hydrophone spacing, there often exists the issue of modulus 2$\pi$ ambiguity. If one ignores the phase difference ambiguity and still estimates the time-delay exploiting the wrapped phase difference measurements, the resulting time-delay estimates deviate from the truth values significantly. The algorithms mentioned above estimate the time-delay difference utilizing the unambiguous low-frequency line-spectrum phase difference measurement from a single frame, and thus require a relatively high SNR to achieve satisfactory performance. Unfortunately, this requirement cannot be met for the underwater ship-radiated noise signal due to the development of the vibration control and noise suppression techniques [5] resulting in a much weaker radiated noise level. Moreover, the phase difference measurements are acquired in each hydrophone element without any array gain [3].

Several improved algorithms have been proposed to acquire the enhanced time-delay difference estimates. For example, an effective outlier-robust Kalman smoother (WORKS) method has been developed in [3] to acquire enhanced time-delay estimates by exploiting the properties of slowly changing time-delay difference in the hydrophone array. This method can alleviate the effects of the outliers resulting from the diverse noise levels in the multiple line-spectrum components. Chinese remainder theorem (CRT) has been reformulated to improve the time-delay difference estimation accuracy by unwrapping the phase difference measurements of line-spectrum components [32,33,34,35,36,37]. This type of method was used to determine an integer (multiple of the time-delay difference) from its remainders (wrapped and noisy phase difference measurements) of several moduli (wavelengths of the line-spectrum components) [32]. For the conventional CRT, a small error in the remainder may cause a large error in the solution of an integer [33]. A robust CRT and its fast implementation have been proposed to reduce the reconstruction error caused by the remainder error [34,35]. In addition, a maximum likelihood estimation-based robust CRT has been developed to address the scenario when the remainder noises do not have the same variance [36]. However, the CRT-based methods have to choose the frequency of the line-spectrum component properly to ensure that there is a positive integer such that the products of this positive integer and the wavelengths of different line-spectrum components are coprime [37]. This is difficult for a passive sonar system since the frequency of line-spectrum components in ship-radiated noise is random, rather than designed artificially. The least-squares phase unwrapping estimator (LSPUE) converted the phase unwrapping problem into the problem of finding the nearest points in a lattice and solved the problem utilizing the lattice point theory in the least-squares sense [38,39]. Unfortunately, a relatively high SNR is required to achieve a satisfying phase unwrapping success rate. Overall, time-delay difference estimation exploiting phase difference measurements of line-spectrum components of the underwater ship-radiated noise signal is still an open problem for beamforming-based signal enhancement in the presence of array shape distortion, particularly in the low SNR case.

In this paper, a novel time-frequency joint time-delay difference estimation method is proposed for signal enhancement in the distorted towed hydrophone array. Unlike these conventional methods, which estimate the time-delay difference utilizing only the unambiguous low-frequency line-spectrum phase difference measurements in a single frame, the proposed method is able to acquire enhanced time-delay difference estimates by fully exploiting phase difference measurements of all detected line-spectrum components in the multiple frames. First, the line-spectrum detection is performed on the pre-enhanced signal based on the hypothetical linear array. Then, a technology combining weighted least-square (WLS) and hidden Markov model (HMM) is developed to obtain the coarse time-delay difference estimates utilizing the unambiguous phase difference measurements of low-frequency line-spectrum components. Next, a data-driven HMM robust to phase difference ambiguity is established to acquire the refined time-delay difference estimates exploiting the phase difference measurements of all detected line-spectrum components. Finally, beamforming based on the refined time-delay difference estimates is employed to achieve the signal enhancement. The performance of the proposed method is verified by both numerical simulations and at-sea experiments. The contributions of the proposed method are summarized as follows.

(1)

A data-driven HMM with robustness to phase difference ambiguity is established to acquire enhanced time-delay difference estimates by taking advantage of the underlying property of slowly changing time-delay difference over time;

(2)

The phase difference measurements available for time-delay difference estimation are extended from that of low-frequency line-spectrum components in a single frame to that of all detected line-spectrum components in multiple frames;

(3)

The signal enhancement performance of the proposed method with a distorted array is close to that of the existing approach with a known array shape, even if the SNR of all the line-spectrum components is as low as 4 dB.

The remainder of this paper is organized as follows. Section 2 introduces the signal model. Section 3 analyses the method for time-delay difference estimation exploiting the phase difference of line-spectrum components. The proposed time-frequency joint time-delay difference estimation method is presented in Section 4. Simulation and experiment results are provided in Section 5. Section 6 concludes the paper.

2. Signal Model

In this section, we first introduce the signal model for underwater ship-radiated noise received by a distorted towed hydrophone array, then present the problem of signal enhancement in the presence of array shape distortion.

2.1. Signal Model of Underwater Ship-Radiated Noise in Distorted Towed Hydrophone Array

Consider a flexible cable with M omnidirectional hydrophones mounted at a fixed spacing d. Because of oceanic currents, hydrodynamics, and tactical maneuvers of the towing platform, the hydrophone array towed behind a maneuvering platform cannot be kept as a straight line [11], as shown in Figure 1. Suppose that the array shape distortion only occurs on the horizontal plane [40]. Take the position of the hydrophone nearest from the platform as the origin of the coordinate system. The hydrophone in the distorted array is denoted by red circles in Figure 1 as opposed to that in the original linear array denoted by blue squares.

Suppose that an underwater acoustic noise signal radiated by a far-field source impinges on the towed array with an incident angle $\vartheta$, which is defined as the angle between the incident direction and the negative x-axis, as shown in Figure 1. The data received by the $m\mathrm{th}$ hydrophone can be expressed as

$$r_m\left(t\right)=s\left(t-{\tau }_m\right)+n_m\left(t\right),m=1,2,\cdots ,M,$$

(1)

where $s\left(t\right)$ represents the waveform of the radiated noise signal received by the reference hydrophone, ${\tau }_m$ denotes the time-delay of the signal propagating from the reference hydrophone to the $m\mathrm{th}$ hydrophone, and $n_m\left(t\right)$ is the additive noise uncorrelated with the signal, respectively. When the first hydrophone is taken as the reference, ${\tau }_m$ is given by

$${\tau }_m=\frac{x_{m}cos\vartheta -y_{m}sin\vartheta }{c},m=1,2,\cdots ,M,$$

(2)

where $\left(x_m,y_m\right)$ represents the coordinates of the $m\mathrm{th}$ hydrophone, and c denotes the sound speed in water.

The spectrum of the received hydrophone data can be written as

$$R_m\left(f\right)=S\left(f\right)e^{-j2\pi f{\tau }_m}+N_m\left(f\right),m=1,2,\cdots ,M,$$

(3)

where $R_m\left(f\right)$, $S\left(f\right)$, and $N_m\left(f\right)$ represent the spectrum of $r_m\left(t\right)$, $s\left(t\right)$, and $n_m\left(t\right)$, respectively. According to the statistical model of the underwater ship-radiated noise [24], the spectrum $S\left(f\right)$ can be decomposed into three components, i.e.,

$$S\left(f\right)=S_{ch}\left(f\right)+S_{mp}\left(f\right)+S_{lm}\left(f\right),$$

(4)

where $S_{ch}\left(f\right)$ represents the stationary continuous spectrum originated from the hydrodynamic noise, $S_{mp}\left(f\right)$ denotes the modulation spectrum caused by the propeller noise, and $S_{lm}\left(f\right)$ is the line-spectrum generated by the machinery noise, respectively. The machinery noise $s_{lm}\left(t\right)$, generated from the inevitable vibration of mechanical equipment such as diesel generator and air conditioning system, is generally denoted by the collection of multiple sinusoidal signals, i.e., [3],

$$s_{lm}\left(t\right)=\sum _{k=1}^K{A}_{k}cos(2\pi f_{k}t+{\varphi }_k),$$

(5)

where $A_k$, $f_k$, and ${\varphi }_k$ represent the amplitude, frequency, and initial phase of the $k\mathrm{th}$ sinusoidal signal, and K denotes the number of the sinusoidal signals. Thus, the machinery noise is also termed as the line-spectrum components of the ship-radiated noise.

2.2. Beamforming-Based Signal Enhancement in the Presence of Array Shape Distortion

Beamforming-based signal enhancement has been widely used in passive sonar systems to acquire the enhanced underwater ship-radiated noise signal from the received data of the hydrophone array [13,25]. It can be written as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y\left(t\right)& =\frac{1}{M}\sum _{m=1}^M{r}_m(t+{\widehat{\tau }}_m)=\frac{1}{M}\sum _{m=1}^{M}s\left(t-{\tau }_m+{\widehat{\tau }}_m\right)+\frac{1}{M}\sum _{m=1}^M{n}_m\left(t+{\widehat{\tau }}_m\right),\hfill \end{array}\end{array}$$

(6)

where ${\widehat{\tau }}_m$ represents the estimate of the $m\mathrm{th}$ time-delay. As indicated by (6), the effect of signal enhancement depends on the estimation accuracy of the time-delay. Specifically, when the estimated time-delay ${\widehat{\tau }}_m$ is coinciding with the true time-delay ${\tau }_m$, the source signal received by various hydrophones of the towed array are coherently summed. On the contrary, noise is incoherently summed when the noise in different hydrophones are uncorrelated. Consequently, the SNR of the source signal in the enhanced signal is M times that of the one in the received hydrophone data [41]. According to (2), the estimate of the $m\mathrm{th}$ time-delay can be given by

$${\widehat{\tau }}_m=\frac{x_{m}cos\widehat{\vartheta }-y_{m}sin\widehat{\vartheta }}{c},$$

(7)

where $\widehat{\vartheta }$ denotes the estimate of the target direction. As indicated by (7), the exact coordinates of each hydrophone in the array are required to estimate the time-delay accurately. However, it is difficult to obtain the exact coordinates for a distorted towed array. If one ignores the array shape distortion and estimates the time-delay based on the hypothetical uniform linear array instead, i.e.,

$${\widehat{\tau }}_m=\frac{\left(m-1\right)dcos\widehat{\vartheta }}{c},$$

(8)

the resulting time-delay estimates deviate from the truth values significantly. Thus, the performance of the beamforming-based signal enhancement degrades severely.

Generally, the power of the line-spectrum components in ship-radiated noise is at least 10 dB higher than that of their nearby continuous spectrum. Thus, the line-spectrum components are readily to be detected and recognized [25], as shown in Figure 2. According to (3) and (5), the phase difference of the $k\mathrm{th}$ line-spectrum component between the $m\mathrm{th}$ and $\left(m-1\right)\mathrm{th}$ hydrophones is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \Delta {\phi }_{m,k}& ={\phi }_{m,k}-{\phi }_{m-1,k}=-2\pi f_k\Delta {\tau }_m,\hfill \end{array}\end{array}$$

(9)

where ${\phi }_{mk}$ represents the phase of the $k\mathrm{th}$ line-spectrum component of the radiated noise signal received by the $m\mathrm{th}$ hydrophone, $\Delta {\tau }_m={\tau }_m-{\tau }_{m-1}$ is the time-delay difference of the radiated noise signal between the $m\mathrm{th}$ and $\left(m-1\right)\mathrm{th}$ hydrophones for $m=2,3,\cdots ,M$, and $\Delta {\tau }_1=0$. The relation between time-delay ${\tau }_m$ and time-delay difference $\Delta {\tau }_m$ is given by

$${\tau }_m=\sum _{i=1}^m\Delta {\tau }_i.$$

(10)

It is noted from (9) and (10) that the time-delay required for beamforming-based signal enhancement can be acquired from the phase differences of line-spectrum components in the underwater ship-radiated noise signal.

3. Time-Delay Difference Estimation Exploiting Phase Difference of Line-Spectrum Components of Underwater Ship-Radiated Noise

In this section, we first introduce the time-delay difference estimation utilizing phase difference of line-spectrum components in underwater ship-radiated noise. Then, we analyze the performance degradation caused by modulus 2$\pi$ ambiguity in phase difference measurement of the high-frequency line-spectrum component.

To acquire the spatial gain of the hydrophone array, we pre-enhance the ship-radiated noise signal exploiting conventional beamforming (CBF) based on the hypothetical uniform linear array, i.e.,

$$y\left[n\right]=\frac{1}{M}\sum _{m=1}^M{r}_m\left[n+\frac{\left(m-1\right)dcos\widehat{\vartheta }}{c}{f}_s\right],n=0,1,\cdots ,N-1,$$

(11)

where $r_m\left[n\right]$ represents the sample of $r_m\left(t\right)$ with a sampling rate $f_s$, and N denotes the number of samples in a single frame of observation. Suppose that K line-spectrum components are detected from the power spectrum of the pre-enhanced signal, and the estimated frequency of the $k\mathrm{th}$ line-spectrum component is denoted as ${\widehat{f}}_k$ for $k=1,2,\cdots ,K$.

The received hydrophone data $r_m\left[n\right]$ with $n=0,1,\cdots ,N-1$ can be decomposed into N narrow-band frequency bins using discrete Fourier transform. The frequency bin containing the line-spectrum component can be approximately modeled as the superposition of the sinusoidal signal and narrow-band noise [3]. Exploiting maximum likelihood, the phase estimation for the $k\mathrm{th}$ line-spectrum component received by the $m\mathrm{th}$ hydrophone can be expressed as [42]

$${\widehat{\phi }}_{m,k}=arctan2\frac{-{\displaystyle \sum _{n=0}^{N-1}}{r}_m\left[n\right]sin\left(2\pi {\widehat{f}}_{k}n/\phantom{2\pi {\widehat{f}}_{k}n{f}_s}{f}_s\right)}{{\displaystyle \sum _{n=0}^{N-1}}{r}_m\left[n\right]cos\left(2\pi {\widehat{f}}_{k}n/\phantom{2\pi {\widehat{f}}_{k}n{f}_s}{f}_s\right)},-\pi \le {\widehat{\phi }}_{m,k}<\pi .$$

(12)

The corresponding probability density function can be expressed as [43]

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & p\left({\widehat{\phi }}_{m,k}\left|{\phi }_{m,k}\right.\right)=\frac{1}{2\pi }{e}^{-{\eta }_{m,k}}\hfill \\ \hfill & +\sqrt{\frac{{\eta }_{m,k}}{\pi }}cos\left({\widehat{\phi }}_{m,k}-{\phi }_{m,k}\right)e^{-{\eta }_{m,k}{sin}^2\left({\widehat{\phi }}_{m,k}-{\phi }_{m,k}\right)}Q\left[-\sqrt{2{\eta }_{m,k}}cos\left({\widehat{\phi }}_{m,k}-{\phi }_{m,k}\right)\right],\hfill \end{array}\end{array}$$

(13)

where $Q\left[·\right]$ represents the right-tail integral of the standard Gaussian distribution, and ${\eta }_{mk}$ denotes the SNR of the $k\mathrm{th}$ line-spectrum component received by the $m\mathrm{th}$ hydrophone, which is defined as the power ratio of the sinusoidal signal to the narrow-band noise with a bandwidth $f_s/\phantom{f_{s}N}N$, i.e.,

$${\eta }_{m,k}=\frac{N{A}_k^2}{2N_m\left(f_k\right)f_s},$$

(14)

where $N_m\left(f_k\right)$ is the power spectral density at frequency $f_k$ of the additional noise received by the $m\mathrm{th}$ hydrophone.

According to (9), the phase difference measurement of the $k\mathrm{th}$ line-spectrum component between the $m\mathrm{th}$ and $\left(m-1\right)\mathrm{th}$ hydrophones is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \Delta {\widehat{\phi }}_{m,k}& ={\widehat{\phi }}_{m,k}-{\widehat{\phi }}_{m-1,k}\hfill \\ \hfill & =\left(\Delta {\phi }_{m,k}+n_{\Delta {\phi }_{m,k}}\right)\left(mod2\pi \right)\hfill \\ \hfill & =\Delta {\phi }_{m,k}+n_{\Delta {\phi }_{m,k}}+2q_{m,k}\pi ,-\pi \le \Delta {\widehat{\phi }}_{m,k}<\pi ,\hfill \end{array}\end{array}$$

(15)

where $\Delta {\phi }_{m,k}=-2\pi f_k\Delta {\tau }_m$ represents the actual phase difference, $n_{\Delta {\phi }_{m,k}}\in \left[-\pi ,\pi \right)$ denotes the noise in phase difference measurement, $\left(\Delta {\phi }_{m,k}+n_{\Delta {\phi }_{m,k}}\right)$ is termed noisy phase difference, $q_{m,k}$ is an integer to ensure that $-\pi \le \Delta {\widehat{\phi }}_{m,k}<\pi$, respectively. It is noted from (15) that the obtained phase difference measurement is the result of the noisy phase difference modulo 2$\pi$. When the obtained phase difference measurement is unambiguous, i.e., $q_{m,k}=0$ or $\Delta {\widehat{\phi }}_{m,k}=\Delta {\phi }_{m,k}+n_{\Delta {\phi }_{m,k}}$, time-delay difference estimation utilizing phase difference measurement of the $k\mathrm{th}$ line-spectrum component is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \Delta {\widehat{\tau }}_{m,k}& =\frac{\Delta {\widehat{\phi }}_{m,k}}{-2\pi {\widehat{f}}_k}=\frac{\Delta {\phi }_{m,k}\left(1+e_{\Delta {\phi }_{m,k}}/\phantom{e_{\Delta {\phi }_{m,k}}\Delta {\phi }_{m,k}}\Delta {\phi }_{m,k}\right)}{-2\pi f_k\left(1+e_{f_k}/\phantom{e_{f_k}{f}_k}{f}_k\right)}\hfill \\ \hfill & \approx \Delta {\tau }_m-\frac{e_{\Delta {\phi }_{m,k}}/\phantom{e_{\Delta {\phi }_{m,k}}2\pi }2\pi }{f_k}-\frac{e_{f_k}\Delta {\tau }_m}{f_k},k=1,2,\cdots ,K,\hfill \end{array}\end{array}$$

(16)

where $e_{\Delta {\phi }_{m,k}}$ and $e_{f_k}$ denote the estimation error of $\Delta {\phi }_{m,k}$ and $f_k$, respectively. It can be noted in (16) that, for a given $\Delta {\tau }_m$, the time-delay difference estimation error is proportional to the estimation errors of both phase difference and frequency of the line-spectrum component, and is inversely proportional to the frequency of the line-spectrum component. The variance of $\Delta {\widehat{\tau }}_{m,k}$ is approximately given by [30]

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{var}\left(\Delta {\widehat{\tau }}_{m,k}\right)& \simeq \frac{{\widehat{\eta }}_k\left({\widehat{\eta }}_{m,k}+{\widehat{\eta }}_{m-1,k}\right)\left(N^2-1\right)+12{\widehat{\eta }}_{m,k}{\widehat{\eta }}_{m-1,k}{\left(f_s\Delta {\tau }_m\right)}^2}{{\left(2\pi f_k\right)}^2\left(N^2-1\right){\widehat{\eta }}_k{\widehat{\eta }}_{m,k}{\widehat{\eta }}_{m-1,k}},\hfill \end{array}\end{array}$$

(17)

where ${\widehat{\eta }}_k$ represents the estimate for SNR of the $k\mathrm{th}$ line-spectrum component in the pre-enhanced signal.

Note that the detected line-spectrum components generally have diverse SNRs and frequencies. To accurately estimate the time-delay difference exploiting the phase difference measurements of all detected line-spectrum components, we can apply WLS on the phase difference measurements, i.e.,

$$\Delta {\widehat{\tau }}_m^{\mathrm{w}}=\frac{-{\displaystyle \sum _{k=1}^K}2\pi {\widehat{f}}_k\Delta {\widehat{\phi }}_{m,k}{\widehat{\eta }}_{m,k}{\widehat{\eta }}_{m-1,k}/\phantom{2\pi {\widehat{f}}_k\Delta {\widehat{\phi }}_{m,k}{\widehat{\eta }}_{m,k}{\widehat{\eta }}_{m-1,k}\left({\widehat{\eta }}_{m,k}+{\widehat{\eta }}_{m-1,k}\right)}\left({\widehat{\eta }}_{m,k}+{\widehat{\eta }}_{m-1,k}\right)}{{\displaystyle \sum _{k=1}^K}{\left(2\pi {\widehat{f}}_k\right)}^2{\widehat{\eta }}_{m,k}{\widehat{\eta }}_{m-1,k}/\phantom{{\left(2\pi {\widehat{f}}_k\right)}^2{\widehat{\eta }}_{m,k}{\widehat{\eta }}_{m-1,k}\left({\widehat{\eta }}_{m,k}+{\widehat{\eta }}_{m-1,k}\right)}\left({\widehat{\eta }}_{m,k}+{\widehat{\eta }}_{m-1,k}\right)},$$

(18)

where the superscript w denotes WLS. Substituting (17) into (18) yields the approximated variance of the WLS time-delay difference estimation, i.e.,

$$\mathrm{var}\left(\Delta {\widehat{\tau }}_m^{\mathrm{w}}\right)\approx \frac{1}{{\displaystyle \sum _{k=1}^K}{\left(2\pi f_k\right)}^2{\eta }_{m,k}{\eta }_{m-1,k}/\phantom{{\sum }_{k=1}^K{\left(2\pi f_k\right)}^2{\eta }_{m,k}{\eta }_{m-1,k}\left({\eta }_{m,k}+{\eta }_{m-1,k}\right)}\left({\eta }_{m,k}+{\eta }_{m-1,k}\right)}.$$

(19)

The estimator in (18) is optimal in terms of minimizing the sum of squared errors, when there is no modulus 2$\pi$ ambiguity in phase difference measurements of all detected line-spectrum components, i.e., $q_{m,k}=0$ or $\Delta {\widehat{\phi }}_{m,k}=\Delta {\phi }_{m,k}+n_{\Delta {\phi }_{m,k}}$ for $k=1,2,\cdots ,K$.

As indicated by (16), for the line-spectrum components with the same SNR, the time-delay difference estimation accuracy is proportional to the frequency of the line-spectrum component. However, for the line-spectrum component with frequency larger than $c/\phantom{c\left(2d\right)}\left(2d\right)$, the absolute value of the actual phase difference may larger than $\pi$. Thus, the obtained phase difference measurement is ambiguous and wrapped by 2$\pi$ radians, i.e., $q_{m,k}\ne 0$ and $\Delta {\widehat{\phi }}_{m,k}\ne \Delta {\phi }_{m,k}+n_{\Delta {\phi }_{m,k}}$. In the remaining parts of this paper, unless otherwise stated, the line-spectrum component with a frequency larger than $c/\phantom{c\left(2d\right)}\left(2d\right)$ is termed high-frequency line-spectrum component, otherwise termed low-frequency line-spectrum component. If one ignores the phase difference ambiguity and still estimates the time-delay difference according to (18) exploiting the wrapped phase difference measurements, the resulted time-delay difference estimation accuracy degrades significantly. Figure 3 shows the time-delay difference estimation results utilizing phase difference measurements from line-spectrum components with frequencies of 30 Hz, 120 Hz, and 480 Hz. It can be noted that, although the time-delay difference estimates obtained from the line-spectrum component with a frequency of 480 Hz have a smaller fluctuation compared to those obtained from the line-spectrum components with frequencies of 30 Hz and 120 Hz, they exhibit significant deviations from the actual values.

In addition, note that the phase difference measurement of the line-spectrum component is sensitive to noise. Moreover, the time-delay difference estimation performance degrades significantly in the low SNR case. Thus, a relatively high SNR is required to achieve satisfactory time-delay difference estimation accuracy, which is difficult in practice, as discussed in Section 1. Therefore, time-delay difference estimation exploiting phase difference measurements of line-spectrum components of the underwater ship-radiated noise signal is still an open problem for beamforming-based signal enhancement in the presence of array shape distortion, particularly in the low SNR case.

4. Proposed Time-Frequency Joint Time-Delay Difference Estimation Method for Signal Enhancement in the Distorted Towed Hydrophone Array

In this section, we propose a time-frequency joint time-delay difference estimation method to obtain the enhanced time-delay difference estimates in the low SNR case. First, we reformulate the HMM for time-delay difference estimation in Section 4.1. Then, we show the proposed method in detail, which is divided into two stages, i.e., the coarse estimation stage in Section 4.2 and the fine estimation stage in Section 4.3. Next, we summarized our scheme for signal enhancement in the distorted towed hydrophone array in Section 4.4. Eventually, we analyze the calculation complexity of the proposed method in Section 4.5.

4.1. HMM for Time-Dealy Difference Estimation

Note that the time-delay difference for ship-radiated noise signal received by adjacent hydrophones of a towed array always changes slowly and continuously. Thus, it is reasonable to model the change of time-delay difference as a first-order hidden Markov process. The HMM is characterized by $\lambda =\left(\mathbf{A},\mathbf{B},\mathbf{\Pi }\right)$, where $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{\Pi }$ represent the state transition probability matrix, observation probability matrix, and initial state probability vector, respectively [44,45]. Let $\mathbf{u}=\left\{u_1,u_2,\cdots ,u_L\right\}$ denote the set of L hidden states (time-delay differences). $u_l$ is uniformly distributed over $\left[\Delta {\tau }_m^{\mathrm{low}},\Delta {\tau }_m^{\mathrm{up}}\right]$ with an interval ${\delta }_m$ for $l=1,2,\cdots ,L$, where $\Delta {\tau }_m^{\mathrm{low}}$ and $\Delta {\tau }_m^{\mathrm{up}}$ represent the lower and upper bounds of the hidden states, respectively. The dimension for the set of hidden states is given by

$$L=⌈\frac{\Delta {\tau }_m^{\mathrm{up}}-\Delta {\tau }_m^{\mathrm{low}}}{{\delta }_m}⌉+1,$$

(20)

where $⌈·⌉$ rounds up to an integer. The state sequence with length T is denoted as ${\mathbf{I}}_m=\left[i_{m,1},i_{m,2},\cdots ,i_{m,T}\right]$, and $i_{m,t}\in \mathbf{u}$ is the time-delay difference state at frame t for $t=1,2,\cdots ,T$. The observation sequence is represented by ${\mathbf{Z}}_m=\left[{\mathbf{z}}_{m,1},{\mathbf{z}}_{m,2},\cdots ,{\mathbf{z}}_{m,T}\right]$, and ${\mathbf{z}}_{m,t}$ is the observation obtained at frame t.

The state transition probability matrix is denoted as $\mathbf{A}={\left[a_{i,j}\right]}_{L\times L}$, where $a_{i,j}$ represents the probability of the state transitioning to $u_j$ at frame t when the state is $u_i$ at frame $t-1$, i.e.,

$$a_{i,j}=p\left(i_{m,t}=u_j\left|i_{m,t-1}=u_i\right.\right).$$

(21)

The state equation of $i_{m,t}$ can be expressed as

$$i_{m,t}=i_{m,t-1}+{\beta }_{m,t}+\epsilon ,t=2,3,\cdots ,T,$$

(22)

where ${\beta }_{m,t}$ represents the change in time-delay difference between adjacent observations, which is an unknown non-random variable determined by the change rate of the target direction, $\epsilon$ denotes the state noise caused by the random change for the positions of target and hydrophones, and is assumed to be normally distributed with mean zero and variance ${\sigma }_{\epsilon }^2$. The $\left(i,j\right)\mathrm{th}$ element of $\mathbf{A}$ is consequently given by

$$a_{i,j}=\frac{{\xi }_i}{\sqrt{2\pi {\sigma }_{\epsilon }^2}}exp\left[\frac{{\left(i_{m,t}-i_{m,t-1}-{\beta }_{m,t}\right)}^2}{2{\sigma }_{\epsilon }^2}\right],$$

(23)

where ${\xi }_i$ is a scaling factor to ensure that $a_{i,j}$ satisfies ${\sum }_{j=1}^L{a}_{i,j}=1$ for $i=1,2,\cdots ,L$. The observation probability matrix is denoted as $\mathbf{B}=\left[b_1\left({\mathbf{z}}_{m,t}\right),b_2\left({\mathbf{z}}_{m,t}\right),\cdots ,b_L\left({\mathbf{z}}_{m,t}\right)\right]$, where $b_j\left({\mathbf{z}}_{m,t}\right)$ represents the probability of observing ${\mathbf{z}}_{m,t}$ given that the state is $u_j$ at frame t, i.e.,

$$b_j\left({\mathbf{z}}_{m,t}\right)=p\left({\mathbf{z}}_{m,t}\left|i_{m,t}=u_j\right.\right).$$

(24)

The initial state probability vector is denoted as $\mathbf{\Pi }=\left\{{\pi }_1,{\pi }_2,\cdots ,{\pi }_L\right\}$, where ${\pi }_i$ denotes the probability of the starting state being $u_i$, which is defined as

$${\pi }_i=p\left(i_{m,1}=u_i\right),i=1,2,\cdots ,L.$$

(25)

Generally, there is no prior information about the time-delay difference at the starting stage of signal enhancement in the presence of array shape distortion. Hence, $i_{m,1}$ is assumed to be uniformly distributed over $\mathbf{u}$, i.e., ${\pi }_i=1/\phantom{1L}L$ for $\forall i$.

To obtain the time-delay difference estimation using HMM, we need to determine the state sequence that maximizes the conditional probability function $p\left({\mathbf{I}}_m\left|{\mathbf{Z}}_m\right.\right)$ given the model $\lambda =\left(\mathbf{A},\mathbf{B},\mathbf{\Pi }\right)$. Therefore, the estimates of time-delay difference (optimal state sequence) is obtained by

$${\widehat{\mathbf{I}}}_m=arg\underset{{\mathbf{I}}_m}{max}p\left({\mathbf{I}}_m\left|{\mathbf{Z}}_m,\lambda \right.\right),$$

(26)

which can be efficiently computed utilizing a dynamic programming method called the Viterbi algorithm [46,47] as presented in Appendix A. The obtained time-delay difference estimates are denoted as $\Delta {\widehat{\tau }}_m=\left[i_{m,1}^{*},i_{m,2}^{*},\cdots ,i_{m,T}^{*}\right]$.

4.2. Coarse Time-Delay Difference Estimation

The proposed method contains two stages, i.e., the coarse and fine estimation stages, present in this and the next subsection, respectively. In the coarse estimation stage, the observation is the WLS time-delay difference estimation exploiting the unambiguous phase difference measurements from low-frequency line-spectrum components, i.e., ${\mathbf{z}}_{m,t}=\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}$. According to (18), $\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}$ is given by

$$\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}=\frac{-{\displaystyle \sum _{k=1}^{K_1}}2\pi {\widehat{f}}_{k,t}\Delta {\widehat{\phi }}_{m,k,t}{\widehat{\eta }}_{m,k,t}{\widehat{\eta }}_{m-1,k,t}/\phantom{{\widehat{\eta }}_{m,k,t}{\widehat{\eta }}_{m-1,k,t}\left({\widehat{\eta }}_{m,k,t}+{\widehat{\eta }}_{m-1,k,t}\right)}\left({\widehat{\eta }}_{m,k,t}+{\widehat{\eta }}_{m-1,k,t}\right)}{{\displaystyle \sum _{k=1}^{K_1}}{\left(2\pi {\widehat{f}}_{k,t}\right)}^2{\widehat{\eta }}_{m,k,t}{\widehat{\eta }}_{m-1,k,t}/\phantom{{\widehat{\eta }}_{m,k,t}{\widehat{\eta }}_{m-1,k,t}\left({\widehat{\eta }}_{m,k}+{\widehat{\eta }}_{m-1,k,t}\right)}\left({\widehat{\eta }}_{m,k}+{\widehat{\eta }}_{m-1,k,t}\right)},t=1,2,\cdots ,T,$$

(27)

where $\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}$ represents the WLS estimation of time-delay difference at frame t, ${\widehat{f}}_{k,t}$, ${\widehat{\eta }}_{m,k,t}$, and $\Delta {\widehat{\phi }}_{m,k,t}$ denote the estimate of frequency, SNR, and phase difference for the $k\mathrm{th}$ line-spectrum component at frame t, T stands for the number of frames (observation sequence length), and $K_1$ is the number of the detected low-frequency line-spectrum components, respectively. According to (19), the estimated variance of $\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}$ is approximately given by

$${\sigma }_{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}}^2\approx \frac{1}{{\displaystyle \sum _{k=1}^{K_1}}{\left(2\pi {\widehat{f}}_{k,t}\right)}^2{\widehat{\eta }}_{m,k,t}{\widehat{\eta }}_{m-1,k,t}/\phantom{{\left(2\pi {\widehat{f}}_{k,t}\right)}^2{\widehat{\eta }}_{m,k,t}{\widehat{\eta }}_{m-1,k,t}\left({\widehat{\eta }}_{m,k,t}+{\widehat{\eta }}_{m-1,k,t}\right)}\left({\widehat{\eta }}_{m,k,t}+{\widehat{\eta }}_{m-1,k,t}\right)}.$$

(28)

Then, the lower and upper bounds of the hidden states can be set as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \Delta {\tau }_m^{\mathrm{low}}=\max \left[\underset{1\le t\le T}{min}\left(\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}-3{\sigma }_{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}}\right),-\frac{d}{c}\right],\hfill \\ \hfill & \Delta {\tau }_m^{\mathrm{up}}=\min \left[\underset{1\le t\le T}{max}\left(\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}+3{\sigma }_{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}}\right),\frac{d}{c}\right].\hfill \end{array}\end{array}$$

(29)

The interval of the hidden states can be set as

$${\delta }_m=\frac{{\sigma }_{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}}}{3\sqrt{T}},$$

(30)

which is small enough to reduce the impact of time-delay difference discretization on estimation accuracy to a negligible level. It can be seen from (20) and (28)–(30) that the dimension for the set of the hidden states of HMM is a function of the observation sequence length and the number, frequency, and SNR of the line-spectrum components.

Although the value of time-delay difference change ${\beta }_{m,t}$ varies from time to time, it can be approximately regarded as a constant value within a short time. Therefore, when the window length of data fitting is set to be $\left(2T_0+1\right)$, the estimate of ${\beta }_{m,t}$ utilizing least-square linear fitting is given by

$${\widehat{\beta }}_{m,t}=\frac{3{\displaystyle \sum _{i=t-T_0}^{t+T_0}}\left(i-t\right)\left(\Delta {\widehat{\tau }}_{m,i}^{\mathrm{w}}-\Delta {\overline{\widehat{\tau }}}_{m,t}^{\mathrm{w}}\right)}{T_0\left(T_0+1\right)\left(2T_0+1\right)},$$

(31)

where $\Delta {\overline{\widehat{\tau }}}_{m,t}^{\mathrm{w}}={\sum }_{i=t-T_0}^{t+T_0}\Delta {\widehat{\tau }}_{m,i}^{\mathrm{w}}/\phantom{{\sum }_{i=t-T_0}^{t-T_0}\Delta {\widehat{\tau }}_{m,i}^{\mathrm{w}}(2T_0+1)}(2T_0+1)$ denotes the mean of WLS time-delay difference estimation results within the sliding window. Then, the state transition probability matrix $\mathbf{A}$ can be acquired according to (23). As for the observation probability matrix, it is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill b_j\left({\mathbf{z}}_{m,t}\right)& =p\left(\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}\left|i_{m,t}=u_j\right.\right)=\frac{1}{\sqrt{2\pi {\sigma }_{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}}^2}}exp\left[\frac{{\left(\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}-u_j\right)}^2}{2{\sigma }_{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}}^2}\right],\hfill \end{array}\end{array}$$

(32)

when $\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}$ is assumed to follow a Gaussian distribution.

After acquiring the HMM parameters for the coarse estimation stage as the methods described above, we can perform the Viterbi algorithm to estimate the time-delay difference. The resulted estimates of time-delay difference are denoted as $\Delta {\widehat{\tau }}_m^{\mathrm{c}}=\left[\Delta {\widehat{\tau }}_{m,1}^{\mathrm{c}},\Delta {\widehat{\tau }}_{m,2}^{\mathrm{c}},\cdots ,\Delta {\widehat{\tau }}_{m,T}^{\mathrm{c}}\right]$ (the superscript c denotes coarse).

4.3. Fine Time-Dealy Difference Estimation

Once the coarse time-delay difference estimates $\Delta {\widehat{\tau }}_m^{\mathrm{c}}$ are obtained, the HMM parameters for the fine estimation stage of the proposed method can be acquired as follows. The lower and upper bounds of the hidden state can be set as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \Delta {\tau }_m^{\mathrm{low}}=\max \left[\underset{1\le t\le T}{min}\left(\Delta {\widehat{\tau }}_{m,t}^{\mathrm{c}}-\sqrt{3}{\sigma }_{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}}\right),-\frac{d}{c}\right],\hfill \\ \hfill & \Delta {\tau }_m^{\mathrm{up}}=\min \left[\underset{1\le t\le T}{max}\left(\Delta {\widehat{\tau }}_{m,t}^{\mathrm{c}}+\sqrt{3}{\sigma }_{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}}\right),\frac{d}{c}\right].\hfill \end{array}\end{array}$$

(33)

The interval of the hidden state is related to all detected line-spectrum components and can be set as

$${\delta }_m=\frac{1}{3\sqrt{T{\displaystyle \sum _{k=1}^K}{\left(2\pi {\widehat{f}}_{k,t}\right)}^2{\widehat{\eta }}_{m,k,t}{\widehat{\eta }}_{m-1,k,t}/\phantom{{\left(2\pi {\widehat{f}}_{k,t}\right)}^2{\widehat{\eta }}_{m,k,t}{\widehat{\eta }}_{m-1,k,t}\left({\widehat{\eta }}_{m,k,t}+{\widehat{\eta }}_{m-1,k,t}\right)}\left({\widehat{\eta }}_{m,k,t}+{\widehat{\eta }}_{m-1,k,t}\right)}}.$$

(34)

After updating the estimate of ${\beta }_{m,t}$ according to (31) utilizing the coarse time-delay difference estimates ${\left\{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{c}}\right\}}_{t=1}^T$, the state transition probability matrix $\mathbf{A}$ can be obtained according to (23).

Unlike the coarse estimation stage, where the observation sequence is the WLS time-delay difference estimates from the unambiguous phase difference measurements of low-frequency line-spectrum components, the observation sequence in the fine estimation stage is the phase difference measurements of all detected line-spectrum components, i.e., ${\mathbf{z}}_{m,t}={\left[\Delta {\widehat{\phi }}_{m,1,t},{\widehat{\phi }}_{m,2,t},\cdots ,{\widehat{\phi }}_{m,K,t}\right]}^T$. Thus, the $j\mathrm{th}$ element of the observation probability matrix is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill b_j\left({\mathbf{z}}_{m,t}\right)& =p\left({\mathbf{z}}_{m,t}\left|i_{m,t}=u_j\right.\right)=\prod _{k=1}^{K}p\left(\Delta {\widehat{\phi }}_{m,k,t}\left|i_{m,t}=u_j\right.\right).\hfill \end{array}\end{array}$$

(35)

In (35), it is assumed that the phase difference measurements for different line-spectrum components are independent of each other. According to (9) and (12), the phase difference $\Delta {\phi }_{m,k,t}$ is a function of ${\phi }_{m,k,t}$ and ${\phi }_{m-1,k,t}$, and the estimation of ${\phi }_{m,k,t}$ and ${\phi }_{m-1,k,t}$ are mutually independent. Thus, the conditional probability of the phase difference measurement can be obtained by calculating the marginal probability density function with respect to ${\phi }_{m-1,k,t}$, i.e.,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & p\left(\Delta {\widehat{\phi }}_{m,k,t}\left|i_{m,t}=u_j\right.\right)\hfill \\ \hfill & =\underset{-\pi }{\overset{\pi }{\int }}p\left({\widehat{\phi }}_{m-1,k,t}\left|{\phi }_{m-1,k,t}\right.\right)p\left(\left({\widehat{\phi }}_{m-1,k,t}+\Delta {\widehat{\phi }}_{m,k,t}-2q\pi \right)\left|\left({\phi }_{m-1,k,t}-2\pi f_{k,t}{u}_j\right)\right.\right)d{\phi }_{m-1,k,t}\hfill \\ \hfill & =p\left(\left(\Delta {\widehat{\phi }}_{m,k,t}-2q\pi \right)\left|-2\pi f_{k,t}{u}_j\right.\right).\hfill \end{array}\end{array}$$

(36)

More specifically, it is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & p\left(\left(\Delta {\widehat{\phi }}_{m,k,t}-2q\pi \right)\left|-2\pi f_{k,t}{u}_j\right.\right)\hfill \\ \hfill & =\frac{1}{2\pi }exp\left(-{\widehat{\eta }}_{m-1,k,t}-{\widehat{\eta }}_{m,k,t}\right)\left\{1+\underset{-\pi }{\overset{\pi }{\int }}\sqrt{{\widehat{\eta }}_{m-1,k,t}/\phantom{{\widehat{\eta }}_{m-1,k,t}\pi }\pi }cos\theta exp\left({\widehat{\eta }}_{m-1,k,t}{cos}^2\theta \right)Q\left(-\sqrt{2{\widehat{\eta }}_{m-1,k,t}}cos\theta \right)d\theta \right.\hfill \\ \hfill & +\underset{-\pi }{\overset{\pi }{\int }}\sqrt{{\widehat{\eta }}_{m,k,t}/\phantom{{\widehat{\eta }}_{m,k,t}\pi }\pi }cos\theta exp\left({\widehat{\eta }}_{m,k,t}{cos}^2\theta \right)Q\left(-\sqrt{2{\widehat{\eta }}_{m,k,t}}cos\theta \right)d\theta \hfill \\ \hfill & \left.+\underset{-\pi }{\overset{\pi }{\int }}2\sqrt{{\widehat{\eta }}_{m-1,k,t}{\widehat{\eta }}_{m,k,t}}Dcos\theta exp\left({\widehat{\eta }}_{m-1,k,t}{cos}^2\theta +{\widehat{\eta }}_{m,k,t}{D}^2\right)Q\left(-\sqrt{2{\widehat{\eta }}_{m-1,k,t}}cos\theta \right)Q\left(-\sqrt{2{\widehat{\eta }}_{m,k,t}}D\right)d\theta \right\},\hfill \end{array}\end{array}$$

(37)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D& =cos\left(\Delta {\widehat{\phi }}_{m,k,t}-2q\pi +2\pi f_{k,t}{u}_j\right)cos\theta -sin\left(\Delta {\widehat{\phi }}_{m,k,t}-2q\pi +2\pi f_{k,t}{u}_j\right)sin\theta \hfill \\ \hfill & =cos\left(\Delta {\widehat{\phi }}_{m,k,t}+2\pi f_{k,t}{u}_j\right)cos\theta -sin\left(\Delta {\widehat{\phi }}_{m,k,t}+2\pi f_{k,t}{u}_j\right)sin\theta .\hfill \end{array}\end{array}$$

(38)

Note that the conditional probability obtained in (37) associates the wrapped phase difference measurement with the actual time-delay difference. Furthermore, it is independent of the modulus 2$\pi$ ambiguity of the phase difference measurements due to the cosine operation shown in (38). Therefore, the data-driven HMM established in the fine estimation stage is robust to phase difference ambiguity, and can be adopted to estimate the time-delay difference from the phase difference measurements of all detected line-spectrum components.

After obtaining the HMM parameters for the fine estimation stage, we can perform the Viterbi algorithm to estimate the time-delay difference. The resulted refined estimates of time-delay difference are denoted as $\Delta {\widehat{\tau }}_m^{\mathrm{r}}=\left[\Delta {\widehat{\tau }}_{m,1}^{\mathrm{r}},\Delta {\widehat{\tau }}_{m,2}^{\mathrm{r}},\cdots ,\Delta {\widehat{\tau }}_{m,T}^{\mathrm{r}}\right]$ (the superscript r denotes refined).

4.4. Summary of the Signal Enhancement in the Distorted Towed Hydrophone Array

Once the refined time-delay difference estimates are obtained, the beamforming-based signal enhancement is conducted as follows

$$y_t\left[n\right]=\frac{1}{M}\sum _{m=1}^M{r}_{m,t}\left[n+f_s\sum _{i=1}^m\Delta {\widehat{\tau }}_{i,t}^{\mathrm{r}}\right],t=1,2,\cdots ,T,n=0,1,\cdots ,N-1,$$

(39)

where $\Delta {\widehat{\tau }}_{1,t}^{\mathrm{r}}=0$ for $t=1,2,\cdots ,T$. The proposed time-frequency joint time-delay difference estimation method for signal enhancement in the distorted towed hydrophone array is summarized in Algorithm 1.

Algorithm 1: Summary of the Signal Enhancement in the Distorted Towed Hydrophone Array

Given the received hydrophone data ${\left\{r_{m,t}\left[n\right]\right\}}_{n=0}^{N-1}$ for $m=1,2,\cdots ,M$ and $t=1,2,\cdots ,T$. (1) Pre-enhance the signal based on the hypothetical uniform linear array according to (11). (2) Detect the line-spectrum components from the pre-enhanced signal, i.e., ${\widehat{f}}_{k,t}$, $k=1,2,\cdots ,K$. (3) Calculate the phase difference of the detected line-spectrum components according to (15), i.e., $\Delta {\widehat{\phi }}_{m,k,t}$, $m,=2,3,\cdots ,M$. (4) Estimate the time-delay difference utilizing the proposed method: Acquire the WLS time-delay difference estimates ${\left\{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}\right\}}_{t=1}^T$ exploiting ${\left\{\Delta {\widehat{\phi }}_{m,k,t}\right\}}_{k=1}^{K_1}$ according to (27). Obtain the coarse time-delay difference estimates ${\left\{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{c}}\right\}}_{t=1}^T$ utilizing ${\left\{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{w}}\right\}}_{t=1}^T$ according to (28)–(32) and the Viterbi algorithm. Acquire the refined time-delay difference estimates ${\left\{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{r}}\right\}}_{t=1}^T$ using ${\left\{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{c}}\right\}}_{t=1}^T$ and ${\left\{\Delta {\widehat{\phi }}_{m,k,t}\right\}}_{k=1}^K$ according to (33)–(38) and the Viterbi algorithm. (5) Perform the signal enhancement exploiting ${\left\{\Delta {\widehat{\tau }}_{m,t}^{\mathrm{r}}\right\}}_{t=1}^T$ according to (39).

4.5. Calculation Complexity

In this subsection, we analyze the calculation complexity of the proposed method and compare it with those of the existing methods, including the CBF, average method [6], and WORKS [3]. In complexity analysis, we neglect all operations not involving M, N, $L^c$, or $L^r$, where $L^c$ and $L^r$ denote the numbers of the hidden states in coarse and fine estimation stages, respectively.

The calculation complexity of the proposed method mainly lies in the direction-finding and signal pre-enhancement, line-spectrum component detection, and phase difference estimation, and coarse and fine time-delay difference estimation. The beam power maximization-based direction-finding with $3M$ candidate beam directions requires the calculation of $\left({log}_{2}N+3M\right)MN$ complex multiplications and additions [40]. The signal pre-enhancement, in which the precise time-delay is achieved by an $N_{fo}$ order fractional time-delay filter, requires the calculation of $\left(N_{fo}+1\right)MN$ real multiplications and additions. The detection and phase-difference estimation of K line-spectrum components requires the calculation of $4\left[\left(K+1\right){log}_{2}N+KM\right]N$ real multiplications and $2\left[2\left(K+1\right){log}_{2}N+KM\right]N$ real additions [48]. For T frames of observation, the coarse and fine time-delay difference estimation requires the calculation of $\left(M-1\right)\left[\left(7T{L}^c+3T+1\right)L^c+\left(7T{L}^r+30KT+1\right)L^r\right]$ real multiplications and $\left(M-1\right)$$\left[\left(3L^c+1\right)T{L}^c+\left(3L^r+8K\right)T{L}^r\right]$ real additions. The complex multiplication can be obtained by four real multiplications and two real additions; the complex addition requires two real additions [48]. Therefore, the overall calculation of the proposed method with T frames of observation is approximately equal to

$$\begin{array}{c}\hfill \begin{array}{c}\hfill T\left\{\left[\left(4{log}_{2}N+12M+4K+N_{fo}\right)M+4K{log}_{2}N\right]N+\left[\left(7L^c+3\right)L^c+\left(7L^r+30K\right)L^r\right]M\right\}\end{array}\end{array}$$

(40)

real multiplications and

$$\begin{array}{c}\hfill \begin{array}{c}\hfill T\left\{\left[\left(4{log}_{2}N+12M+2K+N_{fo}\right)M+4K{log}_{2}N\right]N+\left[\left(3L^c+1\right)L^c+\left(3L^r+8K\right)L^r\right]M\right\}\end{array}\end{array}$$

(41)

real additions.

For the CBF method, the calculation complexity mainly lies in the beam power maximization-based direction-finding. Thus, its overall calculation complexity is approximately equal to

$$\begin{array}{c}\hfill \begin{array}{c}\hfill 4T\left({log}_{2}N+3M\right)MN\end{array}\end{array}$$

(42)

real multiplications and

$$\begin{array}{c}\hfill \begin{array}{c}\hfill 4T\left({log}_{2}N+3M\right)MN\end{array}\end{array}$$

(43)

real additions.

For the average method, the calculation complexity mainly lies in the direction-finding and signal pre-enhancement, line-spectrum component detection, and phase-difference estimation. Therefore, its overall calculation complexity is approximately equal to

$$\begin{array}{c}\hfill \begin{array}{c}\hfill T\left[\left(4{log}_{2}N+12M+4K+N_{fo}\right)M+4K{log}_{2}N\right]N\end{array}\end{array}$$

(44)

real multiplications and

$$\begin{array}{c}\hfill \begin{array}{c}\hfill T\left[\left(4{log}_{2}N+12M+2K+N_{fo}\right)M+4K{log}_{2}N\right]N\end{array}\end{array}$$

(45)

real additions.

For the WORKS method, the calculation complexity mainly lies in the direction-finding and signal pre-enhancement, line-spectrum component detection and phase-difference estimation, and weighted robust-outlier Kalman smoother. The weighted robust-outlier Kalman smoother requires the calculation of $2(M-1)\left[5M(M-1)+17K+11\right]$ real multiplications and $2(M-1)\left[M(M-1)+7K+7\right]$ real additions [3]. Thus, its overall calculation complexity is approximately equal to

$$\begin{array}{c}\hfill \begin{array}{c}\hfill T\left\{\left[\left(4{log}_{2}N+12M+4K+N_{fo}\right)M+4K{log}_{2}N\right]N+2M\left(5M^2+17K\right)\right\}\end{array}\end{array}$$

(46)

real multiplications and

$$\begin{array}{c}\hfill \begin{array}{c}\hfill T\left\{\left[\left(4{log}_{2}N+12M+2K+N_{fo}\right)M+4K{log}_{2}N\right]N+2M\left(M^2+7K\right)\right\}\end{array}\end{array}$$

(47)

real additions.

It is difficult to compare the calculation complexity intuitively from (40)–(47). A specific example is presented to compare the calculation complexity visually. In this example, we take $M=60$, $d=5.0$ m, $\vartheta ={45}^{\circ }$, $N=20000$, $f_s=20$ kHz, $N_{fo}=10$ and $T=50$. In addition, $K=6$ line-spectrum components with frequencies of 53 Hz, 86 Hz, 139 Hz, 382 Hz, 847 Hz, and 1451 Hz are considered, and all the SNRs of the line-spectrum components are 4 dB. Thus, the numbers of the hidden states in coarse and fine estimation stages are $L^c=136$ and $L^r=745$, respectively. The calculation complexities for different methods are listed in

Table 1. The comparison of calculation complexity.

Method	Real Multiplications	Real Additions
CBF	$4.68\times {10}^{10}$	$4.68\times {10}^{10}$
Average	$4.92\times {10}^{10}$	$4.85\times {10}^{10}$
WORKS	$4.94\times {10}^{10}$	$4.86\times {10}^{10}$
Proposed	$6.17\times {10}^{10}$	$6.09\times {10}^{10}$

. It is noted that the calculation complexity of the proposed method is larger than those of the other three methods. However, since the largest calculation lies in the beam power maximization-based direction-finding required by all methods, the calculation complexities of all methods are at the same order of magnitude. Furthermore, the proposed method can realize the enhanced time-delay difference estimation in real-time on an AMD Ryzen 7 4800H CPU, 16G RAM laptop. Under the above conditions, it takes 16.83 s to process 50 frames of observation (the observation time span is 50 s).

5. Numerical Simulations and At-Sea Experiments

In this section, we evaluate the performance of the proposed method utilizing both simulated data and real sea trial data, and compare it with existing methods, including the CBF method, average method [6], and WORKS method [3]. In addition to the results corresponding to the overall algorithm output (fine estimation), we also present the results corresponding to the coarse estimation (marked with WLS-HMM) for comparison.

5.1. Simulation Results

Suppose that the data are collected from a distorted towed array comprising of $M=60$ hydrophones with inter-hydrophone spacing $d=5.0$ m. The shape of the distorted towed array is modeled as a sinusoid, which is the most common distortion form and can be expressed as [49]

$$y=Asin\left(\gamma x\right),$$

(48)

where A and $\gamma$ are a pair of parameters to describe the shape distortion. The x-coordinate for each hydrophone of the distorted towed array can be recursively acquired as [40]

$$x_m=\left\{\begin{array}{ccc}\hfill & 0,\hfill & \hfill m=1,\\ \hfill & x_{m-1}+\frac{d}{\sqrt{1+A^2{\gamma }^2{cos}^2\left(\gamma x_{m-1}\right)}},\hfill & \hfill m=2,\cdots ,M.\end{array}\right.$$

(49)

The corresponding y-coordinate can then be acquired according to (49). In the following simulations, the parameters for array shape distortion are set to be $A=15$ and $\gamma =\pi /\phantom{\pi 240}240$, respectively. The hydrophone coordinates of the actual distorted array are shown in Figure 4 and compared with that of the hypothetical uniform linear array.

Suppose that the underwater acoustic signal radiated by a far-field target impinges on the towed array with an initial angle $\vartheta ={45}^{\circ }$, and the change rate of the target angle is set to be $\mu =0.025{}^{\circ }/\phantom{{}^{\circ }\mathrm{s}}\mathrm{s}$. According to the characteristics of ship-radiated noise, a three-parameter model with $f_m=200$ Hz, $f_c=500$ Hz, and $\kappa =0$ is employed to model the continuous spectrum component [50]. Besides, six sinusoidal signals with frequencies of 53 Hz, 86 Hz, 139 Hz, 382 Hz, 847 Hz, and 1451 Hz are used to model the line-spectrum components. The sound velocity in water is assumed to be 1500 m/s. The observation sequence length is set to be $T=50$, and the sampling number for each observation is $N=20000$ with a sampling rate $f_s=20$ kHz.

5.1.1. Effectiveness Evaluation of the Proposed Method

In this simulation, the effectiveness of the proposed method is investigated in detail. The SNR of all the line-spectrum components is set to be 4 dB. Figure 5a shows the acquired beam power pattern based on the hypothetical uniform linear array. Due to the time-delay mismatch in the distorted towed array, the main beam is split into two peaks, and the maximal peak is located at $\widehat{\vartheta }=37.4^{\circ }$, which deviates from the actual target direction $\vartheta ={45}^{\circ }$. Six line-spectrum components are detected from the power spectrum of the obtained pre-enhanced signal. Figure 5b presents the estimates of the time-delay difference between the first 2 hydrophones for 50 frames of observation. Figure 5c shows the estimates for the inter-hydrophone time-delay difference of the distorted towed array at the 17th observation. As seen from Figure 5b,c, the time-delay difference estimates are fluctuated around the truths for all methods except the CBF method. Moreover, the estimates in the proposed method exhibit a smaller fluctuation and are closer to the actual values, compared to those in other three methods. This is because the proposed method establishes a data-driven HMM with robustness to phase difference ambiguity by fully exploiting the underlying property of slowly changing time-delay difference over time. Thus, it is able to estimate the time-delay difference utilizing the phase difference measurements of all detected line-spectrum components from multiple frames of observation.

Figure 5d shows the average for the power spectrum of the 50 frames of enhanced signal. It can be seen that the power spectrum of the proposed method outperforms those of its counterparts in terms of fidelity, especially for the high-frequency parts of the power spectrum. Moreover, the power spectrum of the proposed method agrees well with that acquired when the accurate hydrophone coordinates are known. To quantify the signal enhancement performance, the amplitude errors between the line-spectrum components in the power spectrum for the considered five methods and those in the power spectrum for a known array shape are listed in

Table 2. Amplitude errors of line-spectrum components.

Frequency (Hz)	53	86	139	382	847	1451
CBF (dB)	5.20	4.75	5.64	13.80	16.15	15.30
Average (dB)	0.20	0.44	1.25	7.90	15.68	14.56
WLS-HMM (dB)	0.05	0.33	0.52	1.98	7.69	14.30
WORKS (dB)	0.02	0.05	0.18	1.40	4.92	9.33
Proposed (dB)	0	0	0.01	0.01	0.04	0.09

. It is observed that the improvement by using the proposed method is evident, and the amplitude errors of line-spectrum components in the proposed method are generally less than 0.1 dB.

5.1.2. Performance Comparison versus the SNR of the Line-Spectrum Component

The effectiveness of the proposed method has been verified in the previous simulation. In this simulation, we analyze the impact of SNR of the line-spectrum component on the time-delay difference estimation accuracy of the proposed method. The estimation accuracy of time-delay difference is evaluated in terms of the root mean square error (RMSE) defined as

$$RMS{E}_{\Delta \tau }=\sqrt{\frac{1}{N_{r}TM}\sum _{i=1}^{N_r}\sum _{t=1}^T\sum _{m=1}^M{\left(\Delta {\widehat{\tau }}_{m,t,i}-\Delta {\tau }_{m,t}\right)}^2},$$

(50)

where $\Delta {\widehat{\tau }}_{m,t,i}$ denotes the estimate for time-delay difference between the $m\mathrm{th}$ and $\left(m-1\right)\mathrm{th}$ hydrophones at frame t of the $i\mathrm{th}$ run, $\Delta {\tau }_{m,t}$ represents the actual time-delay difference between the $m\mathrm{th}$ and $\left(m-1\right)\mathrm{th}$ hydrophones at frame t, $N_r$ stands for the number of Monte Carlo runs and is set to be 200 in this simulation, respectively.

In this simulation, the SNR of the line-spectrum component varies over the interval [−15, 25] dB. Other parameters are the same as those in the preceding simulation. Figure 6 shows the RMSEs of time-delay difference estimates versus the SNR of line-spectrum component for different methods. It is noted that the RMSEs of the proposed method, WORKS method, WLS-HMM method, and average method decrease as the SNR increases for SNR $\ge -5$ dB. However, the RMSE of the CBF method is almost independent of SNR and fixed at a constant. The incorrect assumption of hydrophone coordinates is the main reason resulting in time-delay difference estimation errors in this scenario. Additionally, note that the RMSEs of the proposed method are much smaller than those of its counterparts throughout the SNR interval. This is to be expected, since the proposed method extends the phase difference measurements available for time-delay difference estimation from that of low-frequency line-spectrum components in a single frame of observation to that of all detected line-spectrum components in multiple frames of observation.

5.2. Experimental Results of Real Sea Trial Data

In this subsection, the real data collected at sea are employed to verify the effectiveness of the proposed method. Figure 7 shows the diagram of the sea trial conducted in the South China Sea. As indicated in Figure 7, the depths of the acoustic source, the receiving hydrophone array and the seafloor are 25 m, 50 m, and 105 m, respectively. A signal is transmitted from the power amplifier (acoustic source) below the anchored transmitting ship. The power spectrum of the transmitted signal consists of continuous spectrum component and eight line-spectrum components with frequencies of 36 Hz, 97 Hz, 157 Hz, 310 Hz, 847 Hz, 1500 Hz, 3000 Hz, and 5000 Hz. The spectrum level of the continuous spectrum component is 130 dB at 1 kHz, and the powers of the line-spectrum components are 15 dB higher than their nearby continuous spectrum. The sound speed profile acquired in the sea trial area shows a weak negative gradient. The towing ship runs according to the pre-designed route with a speed of 6 knots. The distance between the acoustic source and the receiving array is about 12 km. The data are collected from a towed array comprising 60 hydrophones uniformly spaced at 6.0 m. Other parameters are the same as those in Figure 5.

Figure 8a presents the direction time record during a turn of the towing ship. It is observed that the towing ship starts the turn at the time of 300 s and ends the turn at the time of 1250 s. Without loss of generality, we select 120 s of data starting from the time of 1000 s as an example to verify the effectiveness of the proposed method. Five line-spectrum components with frequencies of 36 Hz, 97 Hz, 157 Hz, 310 Hz, and 3000 Hz are detected from the pre-enhanced signal based on the hypothetical uniform linear array. Figure 8b shows the estimates of the time-delay difference between the first two hydrophones for 120 frames of observation. Figure 8c presents the estimates for the inter-hydrophone time-delay difference of the towed array at the 67th observation. As seen from Figure 8b,c, whether in time dimension or space dimension, the time-delay difference estimates in the proposed method exhibit a smaller fluctuation compared to those in other three methods.

Figure 8d presents the average for the power spectrum of the enhanced signal in 120 frames of observation. To quantify the signal enhancement performance, the amplitude gains between the line-spectrum components in the power spectrum for the considered four methods and those in the power spectrum based on the hypothetical uniform linear array are listed in

Table 3. Amplitude gains of line-spectrum components.

Frequency (Hz)	36	97	157	310	847	1500	3000	5000
Average (dB)	1.09	3.27	4.0	5.33	1.99	0.23	0.69	0.07
WLS-HMM (dB)	1.16	3.81	4.72	6.91	5.53	0.54	0.86	0.41
WORKS (dB)	1.36	3.84	4.87	7.50	8.53	2.59	2.98	0.57
Proposed (dB)	1.62	4.20	6.95	10.80	12.76	7.64	11.51	6.27

. Note that the improvement by using the proposed method is evident, particularly for the line-spectrum components with frequencies larger than 100 Hz. The time-frequency spectrum of the enhanced signal for different methods are presented in Figure 9a–e. It is observed that the line-spectrum components in the time-frequency spectrum of the proposed method have higher amplitudes than those in the time-frequency spectrum of other four methods. It is consistent with the results shown in Figure 8d and Table 3.

6. Conclusions

In this paper, we proposed a time-frequency joint time-delay difference estimation method for signal enhancement in the distorted towed hydrophone array. The proposed method fully exploits the underlying property of slowly changing time-delay difference over time by modeling the change of time-delay difference as a first-order hidden Markov process. Furthermore, a data-driven HMM with robustness to phase difference ambiguity is established in the proposed method to estimate the time-delay difference. Thus, the phase difference measurements available for time-delay difference estimation are extended from that of low-frequency line-spectrum components in a single frame to that of all detected line-spectrum components in multiple frames. Since the phase difference measurements in both time and frequency dimensions are exploited jointly, the proposed method has the capability of acquiring enhanced time-delay difference estimates even in the low SNR case. Simulation results show that the signal enhancement performance of the proposed method with a distorted array is close to that of the existing approach with a known array shape, even if the SNR of all the line-spectrum components is as low as 4 dB. At-sea experimental results prove that, even though the signal of interest is contaminated by towed platform noise and real ocean ambient noise, the proposed method still achieves a better signal enhancement performance compared to these existing state-of-the-art approaches.