1   Introduction

In recent years, the increasing demand for applications such as radar, sonar, audio signal processing, and underwater communications has drawn considerable research attention to time delay estimation, a critical step in source localization (Li et al., 2019; Nicos et al., 2024; Shi et al., 2019). The matched filter (MF) is a traditional time delay estimation algorithm known for its low computational complexity, making it widely adopted in engineering practices (Hama and Ochiai, 2019; Challinor and Cegla, 2024). However, its time resolution is susceptible to signal bandwidth limitations, and its performance is vulnerable to reverberation interference in active sonar systems.

To address the above problems, Jiang extended the inverse folding technique to the time domain by convolving the MF output and proposed a high-resolution two-dimensional inverse folding method to enhance temporal resolution (Jiang et al., 2020).

In many practical applications, source signals are often non-Gaussian, while the additive noise is typically Gaussian or even correlated. In this case, Wu (2000) proposed a delay estimation method based on the third-order statistics of the signal. Tang et al. (2004) introduced a delay estimation method using the four-order statistics of the signal.

The multiple signal classification (MUSIC), originally proposed by Schmidt, is a super-resolution delay estimation method based on feature decomposition theory. This method compensates for the limited accuracy caused by the bandwidth-constrained noise of transmitting signals under special conditions (Drake et al., 2023; Zhang et al., 2020). However, the performance of the traditional MUSIC algorithm is limited by the dimension of the autocorrelation matrix M. A large M value may lead to substantial deviations in the estimation value, while a M value can degrade the temporal resolution of the algorithm due to a reduction in M. In recent years, domestic and foreign scholars have introduced various improvements to the MUSIC algorithm. The Root-MUSIC algorithm enhances performance when processing small-sample received data (Yan et al., 2020). The TST-MUSIC algorithm improves performance under low signal-to-noise ratio (SNR) conditions (Wang et al., 2001). Meanwhile, the PM-MUSIC algorithm eliminates the need for complex matrix calculations such as eigenvalue decomposition and singular value decomposition of the received data’s autocorrelation matrix (Qasaymeh et al., 2009), thereby substantially reducing the complexity of the algorithm.

Xia et al. (2016) transformed the geometric acoustic scattering with different time delays in the target echo into multiple single-frequency components and applied the MUSIC algorithm to enhance spectral resolution when the frequency intervals were smaller than the resolution limit of the Fourier transform. Compared to traditional time-delay estimation methods, this approach offers improved performance, characterized by a sharper main flap and reduced side flap interference.

Multipath signals arrive at the receiver with varying angles and path delays. Considering this phenomenon, Wang et al. (2023) proposed a low-complexity framework based on sparse Bayesian learning for joint estimation of the angle and delay of received signals. This approach offers considerable advantages over traditional methods, such as higher accuracy and resolution. While the dual-sensor differential denoising approach is effective for suppressing pump noise, accurately identifying the time delay remains challenging. Yan et al. (2023) introduced an improved variable step-size adaptive time delay estimation method, TVSS-LMSTDE. Additionally, Ding et al. (2023) proposed a novel time-delay estimation method within the framework of variational Bayesian inference. Unlike conventional grid-less compressive time-delay estimation methods, this method does not require parameter adjustment and can automatically estimate the number of time delays.

The frequency sliced wavelet transform (FSWT) algorithm offers superior resolution in the time-frequency images of time delay signals, making it widely applicable in fields such as medical signal processing and power system analysis (Zhao et al., 2019). Sheng et al. (2021) leveraged the unrestricted filtering properties of FSWT to reconstruct and extract composite spectral boundaries. Additionally, Hu et al. (2021) successfully extracted subcomponents containing fault information from bearing vibration signals by integrating FSWT with an improved particle swarm optimization algorithm and empirical wavelet transform.

In this paper, based on the principles of linear geometric acoustics, a bright spot model for typical underwater target echoes is constructed. The scattering characteristics of echo signals at various incidence angles are analyzed. On this basis, echo signals in a reverberation background are simulated and processed using the fractional order Fourier transform–frequency sliced wavelet transform (FrFT–FSWT) algorithm, which indirectly enhances the accuracy of time delay estimation.

The rest of the paper is organized as follows:

1) Section 2 introduces the characteristics of geometric acoustic scattering by integrating methods from physical acoustics. This section derives mathematical expressions for the relative time delays of individual bright spots in the experimental target model, based on the geometric parameters of typical bright spot shapes. Furthermore, the effects of different incidence angles on the echoes from geometric bright spots are simulated and analyzed.

2) Section 3 introduces the theoretical framework of the FSWT. On this basis, the fractional-order Fourier transform (FrFT) is introduced, and the optimal parameters in the fractional-order domain are determined. The performance improvements of the FrFT–FSWT method, particularly in terms of interference suppression and time–frequency aggregation, are also analyzed.

3) Section 4 constructs the time–frequency observation matrix of the bright spot model using the FrFT–FSWT and provides a detailed description of the proposed delay estimation algorithm.

4) Section 5 verifies the superiority of the proposed delay estimation algorithm over the conventional MF through simulation.

5) Section 6 further validates the reliability and effectiveness of the proposed delay estimation algorithm through the analysis of experimental data.

6) Section 7 summarizes the key contributions and findings of this study.

2   Fundamentals of the target highlight model

Active sonar systems emit acoustic signals to stimulate underwater targets, generating echo responses. These echoes result from a physical interaction that depends not only on the intrinsic vibrations or fluctuating characteristics of the target but also on the characteristics of the incident acoustic wave. In most cases, the incident wave acts as a small-amplitude excitation source, and the resulting echo follows linear acoustic principles and is independent of the initial transmission time. For engineering applications, the detailed generation mechanism of the echo is often simplified by modeling the system as linear and time-invariant, wherein the echo represents the response of the targets to the incident acoustic wave. In the high-frequency domain, the main echo components typically include mirror reflection waves, angular scattering waves, and various elastic scattering waves. These components can be approximated as echoes from discrete scattering centers, commonly referred to as bright spots. Accordingly, the entire target can be modeled as a spatial distribution of bright spots, a representation known as the bright spot model. The wavelet generated by each bright spot coherently superimposes to form the complex overall echo of the target. According to their formation mechanism, bright spots are generally categorized into two types: elastic bright spots and geometric bright spots, which correspond to elastic scattering waves and geometric scattering waves, respectively (Zhou et al., 2024).

An elastic bright spot is an equivalent, rather than a physical, bright spot. Elastic echoes can arise from various mechanisms; for example, when elastic scattered waves or surface traveling waves are present, elastic bright spots may appear at specific incident angles and positions, often accompanied by dispersion phenomena. The propagation characteristics and energy radiation modes of different elastic highlights are inconsistent, and their analysis typically relies on wave theory or geometric diffraction theory.

Geometric highlights generally fall into two categories. The first type is the mirror reflection highlights, generated by reflections on convex smooth surfaces, providing the most contribution when the curvature radius is large. The second type includes scattering bright spots generated at the edges of the target (Shen et al., 2022). The sound center of geometric highlights aligns with the geometric center of the target, and their characteristics are mainly determined by the target’s geometric scale. Among the various types of bright spots in target echoes, geometric bright spots are relatively stable, have simple causes, and do not exhibit dispersion. Therefore, this study focuses on the analysis of geometric bright spots, disregarding elastic bright spots.

The target comprises a set of bright spots, each associated with a transfer function that reflects the properties of the target. The transfer function of a single bright spot is defined in Eq. (1).

where r represents the straight-line distance between the target and the bright spot, and A(r, ω) represents the amplitude reflection factor, which is typically a function of the frequency. For narrowband signals, the amplitude reflection factor can be approximated by its value at the center frequency. The parameter τ represents the delay factor, determined by the equivalent scattering center and the acoustic range between a reference point ξ, defined as τ= 2ξ/c. The acoustic range ξ is a function of the acoustic wave incidence angle θ. The parameter φ denotes the phase factor, representing the phase shift that occurs during the formation of the target echo. This phase shift depends on the shape of the target and the characteristics of the bright spot.

The echo from any complex target can be viewed as the coherent superposition of multiple wavelets. By combining the transfer functions of individual bright spots, the transfer overall function of H(r, ω) of the target is defined in Eq. (2).

where M represents the total number of bright spots in the target echo. A_m, τ_m, and φ_m are the values of the amplitude reflection, time delay, and phase factors, respectively, for different bright spots.

When the transmitted signal is a linear frequency modulated (LFM) signal, the target echo can be defined based on the transfer function of the highlight model, as shown in Eq. (3).

where f₀ is the center frequency of the LFM pulse, k is the frequency modulation slope, and T is the signal pulse width. The experimental target model is shown in Figure 1.

Figure  1  Experimental target model

Download: Full-Size Img

The experimental target model comprises seven geometric scattering centers, with the labels and specific positions of each scattering center shown in Figure 1. The target body is a cylindrical shell with a hollow interior. Scattering centers #1, #2, #3, and #4 are associated with angular reflection, while scattering center #5 is associated with mirror reflection. Scattering center #6 corresponds to hemispherical scattering and is located at the intersection point between the hemisphere and the line passing through the center of the hemisphere. Scattering center #7 is a planar scatterer, located at the center of the bottom of the model. The incident sound field is a plane wave; thus, the incident angle of the sound wave for all scattering centers is θ. Based on the properties of geometric bright spots mentioned earlier, the main observable scattering centers in the target echo follow the transformation rules with respect to the incident angle θ as outlined below.

1) When θ=0°, the main scattering center #6 can be observed;

2) When 0° < θ < 90°, the main scattering centers #1, #4, and #6 can be observed;

3) When θ=90°, the main scattering center #5 can be observed;

4) When 90° < θ < 180°, the main scattering centers #1, #2, and #4 can be observed;

5) When θ=180°, the main scattering center #7 can be observed.

As the incident angle θ changes, the time delay of each scattering center on the target also varies, contributing to the dynamic characteristics of the geometric bright spot structure in the target echoes. Taking the center O of the experimental target model as the reference point, the time delay expression for each acoustic scattering center relative to O is calculated and is defined in Table 1.

Based on the time delay expressions provided in Table 1, and assuming the experimental target l= 210 cm and a bottom diameter a= 53.3 cm, the variation in relative time delay for each acoustic scattering center with respect to the reference point O, as the acoustic wave incidence angle changes from 0° to 180°, is shown in Figure 2.

Figure  2  Relationship between time delay of acoustic scattering center and incidence angle of acoustic wave

Download: Full-Size Img

As shown in Figure 2, the time delay curve as a function of the incident angle θ is asymmetric at 90°. Assuming the transmitted signal is an LFM signal, the effect of different incident angles on the geometric bright spot echoes is shown in Figure 3.

Figure  3  Time-domain signals of bright spot models with different angles of incidence

Download: Full-Size Img

According to Figures 3(a) to (d), as the incident angle of the sound wave increases sequentially, the delay difference between highlights #1 and #2 gradually decreases, while that between the two and highlight #4 increases. Therefore, in the early part of the waveform (Figures 3(a) to (d)), the duration of the LFM signal shortens, whereas in the later part of the waveform, the duration becomes longer. Notably, in Figure 3(d), the echo from highlight #4 is completely separated from those of #1 and #2. Based on the incident angle of the sound wave in Figure 1, when the incident sound wave approaches the transverse plane of the target, the amplitude of the geometric highlight echoes in Figure 3 decreases, and the echo energy gradually weakens.

3   Signal preprocessing

From the energy viewpoint, target echo signals exhibit strong energy aggregation characteristics in the time-frequency domain, whereas reverberation signals display a random distribution with poor energy aggregation. This distinction enables effective separation of the two in the time-frequency domain. Aiming to suppress the influence of reverberation and noise on the echoes, this section proposes a joint FrFT–FSWT algorithm designed to improve the time-frequency clustering of the target signals.

3.1   Frequency slice wavelet transform

FSWT is an advanced time–frequency signal analysis method that combines the advantages of the short-time Fourier transform (STFT) and the wavelet transform (WT), while also making improvements to both. FSWT introduces the frequency slice function (FSF), which overcomes the fixed time and frequency resolution limitations of the STFT and addresses the challenge of selecting an appropriate wavelet basis in WT. Additionally, FSWT avoids the cross-term problem encountered in the Wigner–Ville distribution (WVD) when analyzing multiple frequency components (Yan et al., 2011). Therefore, FSWT is used in this paper to reconstruct and separate the desired signal components.

For the signal x(t), if it satisfies ∫|x(t)|²dt < ∞, then the FSWT is defined in Eq. (4).

where the scale factor σ and energy factor λ are constants or functions of ω, u, and t; ˆ represents the frequency spectrum after the corresponding Fourier transform; * represents conjugate; ω and t are the observation frequency and observation time, respectively; u represents the evaluation frequency; and p(t) represents the FSF.

FSWT is reversible, allowing for signal decomposition in the time and frequency domains, followed by filtering or segmenting the signals for processing. According to Yan et al. (2011), if $\hat{p}(\omega) $ satisfied $\hat{p}(\omega) $=1, then the frequency slice inverse wavelet transform (IFSWT) can be simplified as shown in Eq. (5).

As shown in Eq. (5), signal reconstruction based on IFSWT is independent of FSF selection. Compared to WT, IFSWT avoids the issue of selecting an appropriate wavelet basis during the signal reconstruction process.

From Eqs. (4) and (5), the FSWT–IFSWT realizes the function of a bandpass filter in the time and frequency domains. Therefore, the SNR of the reconstructed signal, obtained after applying FSWT–IFSWT, is improved.

The scale factor σ plays a more substantial role than the energy factor λ in the FSWT. Typically, λ=1, the scale factor is not constant. According to the concept of Morlet transformation, σ∝ω.

According to Yan et al. (2011), the range of κ is shown in Eq. (7).

where µ and v are related to FSF. The FSF typically takes $\hat{p}(\omega)=\mathrm{e}^{-(1 / 2) \omega^2}, $; thus, $\mu=1 / 2, \quad v=\mathrm{e}^{-(1 / 2) \mu} \approx 0.7788 $, and κ≈0.707/η. The frequency resolution η is shown in Eq. (8).

where Δω denotes the frequency interval between two adjacent points in the frequency domain in the time–frequency diagram and ω_x is the center frequency of the signal x(t).

For single-component CW signals, the parameter κ can be solved accurately. However, if the signal x(t) to be measured contains multiple frequency components, such as an LFM signal whose frequency changes over time, then κ cannot be calculated accurately. In such cases, κ is typically taken as an empirical value, and the value of κ is determined based on the expected response ratio η, which is also empirical. This approach has limitations in underwater acoustic signal processing. As shown in Figure 4, when using FSWT to process an LFM signal, the energy of the signal is not constant across the bandwidth range; it increases from weak to strong as time progresses. Additionally, the time–frequency aggregation of the signal deteriorates.

Figure  4  FSWT processing results of LFM signals under different conditions

Download: Full-Size Img

3.2   Definition of fractional fourier transform

This subsection introduces the FrFT, which has a basis decomposition property for LFM signals (Omair and Ahmet, 2024). When the appropriate fractional order is selected, the LFM signal behaves as an impulse function in the fractional order domain (Sun et al., 2015; Guo and Yang, 2022). This condition indicates that the LFM signal is transformed into a single-frequency signal, allowing for the precise calculation of the value κ_u of κ in the u-domain.

From a linear integration perspective, let the signal to be measured be denoted as x(t), and its FrFT is defined in Eq. (9).

where p is the order of FrFT, which can take any real value, and its relationship with the rotation angle α is shown in Eq. (10).

where F^P is the operator symbol of FrFT, and K_p(t, u) is the transform kernel. The specific expression is shown in Eq. (11).

In Eq. (11), when p= 4n, the rotation angle is 2nπ, $K_p(t, u)=\delta(u-t) $; when p= 4n± 2, then α=(2n+1)π, $K_p(t, u)=\delta(u+t) $. Considering the above derivation, the order in FrFT is based on a period of 4. Eq. (11) is replaced with a variable substitution. When $u=u / \sqrt{2 \pi} $, t=$t / \sqrt{2 \pi} $, the expansion of Eq. (11) is shown in Eq. (12).

Based on Eqs. (12) and (9), FrFT is a linear transformation. Its linearity indicates that the FrFt of the echo is equal to the sum of the FrFts of each component, simplifying the analysis of the target echo and reverberation. Additionally, p is based on a period of 4; thus, only the case − 2≤p≤2 needs to be considered in practical applications. When p=1, FrFT corresponds to the traditional Fourier transform.

The inverse transformation of FrFT is given in Eq. (13).

where X_p(u) is the function space of x(t) in the inverse transformation kernel K_-p(t, u), which is an orthogonal LFM signal group in the u-domain.

3.3   Parameter solving for fractional domain

From the perspective of the time–frequency plane, FrFT denotes the result of rotating the time–frequency plane by an angle of α with the origin as the center. This rotation maps the signal into a new domain, known as the u-domain. As shown in Figure 5, the LFM signal is distributed along a straight line in the time–frequency plane, where its projection on the time axis t corresponds to the pulse width T of the signal, and its projection on the frequency axis ω represents the bandwidth B of the signal. The expression of the LFM signal in the time domain is given by the following: x(t)=exp(i2πf₀t+iπkt²), t∈[0, T], where f₀ is the signal center frequency and k is the frequency modulation slope. x(t) is substituted into Eq. (9) to obtain the p-order FrFT mathematical expression of the LFM pulse, as shown in Eq. (14).

Figure  5  Fractional Fourier transform under the rotation angle of α

Download: Full-Size Img

Eq. (14) provides the analytical solution for the energy accumulation position of the LFM pulse in the u-domain. By selecting the optimal order p, the peak point u₀x(t) in the u-domain can be determined. That is, the energy of the LFM signal, which is dispersed across multiple frequency points, will concentrate at the position of u₀ in the u-domain. Therefore, the LFM pulse similarly behaves like a narrowband signal in the u-domain. Substituting Eq. (14) into Eq. (13) yields the narrowband signal x_u(v).

Two methods are used to calculate the angle of rotation: the first method involves calculating the theoretical value, as shown in Eq. (15).

where f_s is the sampling frequency, N is the number of data sampling points, and k is the frequency modulation slope. The second method aims to determine the location of the peak point on the u-v plane. However, the first method is used in this paper.

For the LFM signal shown in Figure 4, Gaussian white noise with different SNRs and reverberation with different signal-to-reverberation ratios (SRRs) are added to generate the signal x(t). FSWT and IFSWT are then successively applied to the signal x(t), and the SNR and SRR of the processed signals are calculated to evaluate the algorithm’s performance. Simultaneously, using FrFT, the signal x(t) is rotated into the u-domain to obtain a single-frequency signal x_u(t). The simulation result of the LFM signal after fractional Fourier transform is shown in Figure 6.

Figure  6  Change of signal before and after FrFT

Download: Full-Size Img

As shown in Figure 6, due to differences in time–frequency concentration between noise, reverberation, and the LFM signal in the u-domain, the total energy of the signal remains unchanged before and after the FrFT. However, the distribution of this energy is altered. As described in Section 3.1, applying the FSWT–IFSWT to a signal is equivalent to its passage through a bandpass filter within a specific frequency band. Meanwhile, FrFT concentrates the signal’s energy at a specific frequency point. Based on the conditional distribution characteristics of reverberation and noise, combining FrFT with FSWT–IFSWT enables highly effective suppression of reverberation and noise. For x_u(t), FSWT and IFSWT are applied sequentially, and the SNR and SRR of the processed signal are calculated. The results are shown in Figure 7.

Figure  7  Change curves of SNR and SRR before and after signal processing

Download: Full-Size Img

As shown in Figure 7, the joint FrFT–FSWT time–frequency analysis algorithm proposed in this paper demonstrates superior interference suppression capabilities compared to the standalone FSWT algorithm.

4   Algorithm of time delay estimation

Aiming to suppress noise in target echo signals, the FSWT, known for its noise immunity, is chosen. However, in active sonar systems, the main source of interference is reverberation. Therefore, a time–frequency observation matrix based on geometric acoustic scattering characteristics can be constructed. By leveraging the differences in the time–frequency distribution between reverberation and target echoes, reverberation can be effectively suppressed. The target echo signals, after interference suppression, can then be reconstructed using IFSWT. Finally, the arrival time of the geometric echo signals at the active sonar receiver can be estimated (Zhou et al., 2022; Zhao et al., 2022).

The schematic of the underwater target and sonar system is shown in Figure 8. In this figure, H denotes the water depth, and θ is the angle between the incident acoustic wave and the horizontal plane. The distance between the target and the transceiver (combined transducer) is cτ/2, where c denotes the speed of sound and τ is the absolute time delay measured by the transducer. The relative time delay τ0τ0 between different target prisms depends only on the angle of incidence θ and the target’s physical dimensions, l and a.

Figure  8  Schematic of the relative position of the target and the transceiver-combined sonar

Download: Full-Size Img

The incident acoustic wave is assumed to be a far-field plane wave. Therefore, the arrival time at each prism is denoted as τ. The target is placed on the seafloor, which is assumed to be flat, allowing the depth of the target to be known. The target is located in a shallow sea area. Thus, the effects of the sound speed gradient and reflections from the bottom of the sea surface on the acoustic wave are neglected.

In this paper, the experimental target model is not left–right symmetric. If the target and the transceiver-combined sonar’s position relationship is shown in Figure 8(a), then the formula represents the first half of Eq. (16). If the target and the transceiver-combined sonar’s position relationship is shown in Figure 8(b), then the formula represents the second half of Eq. (16). The specific steps for the time-delay estimation algorithm are as follows:

1) The target echo signal is first subjected to the FrFT to obtain the domain signal. The FrFT preserves the total energy of the signal but alters its energy distribution. Thus, the FrFT process helps suppress some of the reverberation and noise.

2) The FSF is selected, and the parameter κ_u is computed to obtain the FSWT’s time–frequency image W_xu(v, u) of x_u(v).

3) The range of τ, τ_i∈[t₁, t₂], and τ_i is determined in the steps of τ_s=0.01N/f_s, where N is the number of sampling points of the echo signal and f_s is the sampling frequency. The acoustic wave incidence angle θ_i for the two cases corresponding to each τ_i is determined in accordance with Eq. (16). Meanwhile, in the case of known target size information, Table 1 is used to obtain the relative time delay information between each prism, and different bright spot models $S_{h_i}\left(t, \tau_i\right) $ are identified, as defined in Eq. (17). The calculations of τ₁, τ₂, τ₄, and τ₆ in Eq. (17) are combined with Table 1, where the two expressions in Eq. (17) use crown #6 and prong #1 as reference points. The transmitting signal is set to be an LFM signal.

4) S_{h_i}(t, τ_i) is rotated using FrFT to obtain a bright spot model in the domain S_{h_i_u}(v, τ_i). The FSWT corresponding to S_{h_i_u}(v, τ_i) is calculated to obtain its time–frequency observation matrix W_{Sh_i_u}(v, u, τ_i).

5) The interval in the time-frequency image W_xu(v, u) for signal reconstruction is determined based on the time–frequency observation matrix W_{Sh_i_u}(v, u, τ_i), as shown in Eq. (18).

6) The filtered time–frequency image Y_i(v, u, τ_i) is calculated and then substituted into Eq. (5) to obtain the reconstructed signal y_i. The definition of y_i is given in Eq. (19).

7) Two evaluation parameters are introduced for estimating the absolute delay τ of the echo reaching the receiver:

(ⅰ) The first term is the mean square error, defined in Eq. (20) as follows:

(ⅱ) The second term is the waveform similarity coefficient, as shown in Eq. (21).

8) Steps 3) to 7) are repeated in accordance with the range and step size in Step 3).

(ⅰ) The mean square error curve is plotted based on the results of e_{MSE_i} at each time, and the lowest point of the curve is identified, as shown in Eq. (22).

(ⅱ) The waveform similarity coefficient curve is plotted based on the results of ρ_{WSC_i} at each time, and the highest point of the curve is then determined, as shown in Eq. (23).

τ_ρ is the optimal delay estimation result based on the waveform similarity coefficient. If τ_e and τ_ρ are equal to τ_i, then τ_i is the arrival time of the target echo at the receiving end of the active sonar. When the SNR or SRR is notably low, a certain difference may exist between τ_ρ and τ_e in the ranges of 0°≤θ < 90° and 0°≤θ < 90°. During this time, the angle range corresponding to the small difference between τ_ρ and τ_e is selected, and τ_i is recorded as the average of τ_ρ and τ_e. Once τ_i is determined, the distance between the target and the transducer, as well as the acoustic incidence angle θ, can be determined.

In shallow water environments, the influence of sound velocity gradients and seabed roughness on time delay estimation is negligible. When the target tilts slightly in the vertical direction relative to the seabed, the true θ value can be recorded as θ_true=θ+Δθ (Δθ≪θ). According to Eq. (16), the difference between the actual delay value and the measured delay value can be expressed as Eq. (24).

In shallow sea environments, H is relatively small. Thus, when Δθ is substantially smaller than θ, the value of ΔτΔτ is negligible.

5   Simulation analysis in the context of strong reverberation

The unit scattering model is used, with the number of scatterers set to 100. The transmitted signal is an LFM signal with the following parameters: a center frequency of 37.5 kHz, a bandwidth of 25 kHz, a pulse width of 2 ms, and a reverberation duration ranging from 0.02 s to 0.08 s. The target echo model at the receiving end follows the bright spot model S_h(t), with an acoustic incidence angle of 120°. The target model and dimensions are consistent with those described in Section 2. The SRR is set to -3 dB, and S_h(t) is located in the reverberation at 0.034 6 s, at which time S_h(t) is submerged within the reverberation. As shown in Figure 8(b), the water depth H is set to 300 cm. When the bottom surface of the target is facing the receiving end, the acoustic wave incidence angle is 60°, which corresponds to an angle of θ= 120° between the incident wave and the target surface. The straight-line distance between the target and the receiving end is 346 cm, and the speed of sound c= 1 500 m/s. Under these conditions, the target echo arrives at the receiving end at time τ= 0.004 6 s. Aiming to simulate the effect of noise, Gaussian noise with an SNR of -5 dB is added to the signal in the interval between 0.03 s and 0.04 s, resulting in the signal to be processed x(t).

First, x(t) is subjected to the FrFT to obtain a single-frequency signal x_u(v) in the u-domain. However, the time–frequency intervals of the signal x_u(v) cannot be observed.

The water depth H is 300 cm, and the sound wave is incident on the bottom of the water. The time it takes for the echo to reflect and reach the receiving end is 0.004 s. Therefore, the absolute time delay τ of the echo arriving at the receiving end should satisfy τ_i∈[0.004 1 s, 0.01 s]. However, in this range, the bright spot model corresponds to the longest time width of 0.004 2 s0.004 2 s. Consequently, the range of τ is adjusted to τ_i∈[0.004 1 s, 0.005 8 s], and the step size is set to 0.000 1 s. As shown in Figure 7, once the distance cτ/2 is determined, the target may be in two possible states: one where the target’s hemispherical head is facing the sonar, and the other where the target’s bottom surface is facing the sonar. Therefore, constructing the bright spot model S_h(t) in the two cases is necessary. When the acoustic wave incidence angle is between 0° and 90°, the scattering centers that contribute to the bright spot model are #1, #4, and #6. The bright spot model is S_{h_i_0_90}(t, τ_i). Similarly, when the acoustic wave incidence angle is between 90° and 180°, the scattering centers contributing to the bright spot model are #1, #2, and #4. In this case, the bright spot model is S_{h_i_90_180}(t, τ_i). Correspondingly, two time–frequency observation matrices, W_{i_0_90}(v, u, τ_i) and W_{i_90_180}(v, u, τ_i), need to be constructed. Therefore, two mean square error curves and two waveform similarity coefficient curves will be obtained in Step 8) of Section 4, as shown in Figure 9.

Figure  9  Mean square error curves and waveform similarity coefficient curves for different scenarios

Download: Full-Size Img

When the angle of acoustic incidence is within 0° < θ < 90°, the observation results of W_{i_0_90}(v, u, τ_i) are shown in Figures 9(a) and (c). The lowest point of the mean square error curve is located at τ_e=0.004 6 s, with a corresponding error value of e_{MSE_0.004 6}≈0.09. Meanwhile, the highest point of the waveform similarity coefficient is located at τ_ρ=0.005 s, with a corresponding waveform similarity coefficient ρ_{WSC_0.005}≈0.09. τ_ρ and τ_e are not the same; therefore, the angle of acoustic incidence is not within 0° < θ < 90°. When the angle of acoustic incidence is within 90° < θ < 180°, the observation results of W_{i_90_180}(v, u, τ_i) at this time, as shown in Figures 9(b) and (d). The lowest point of the mean square error is located at τ_e=0.004 6 s, corresponding to e_{MSE_0.004 6}≈0.09, while the highest point of the waveform similarity coefficient is located at τ_ρ=0.004 6 s, corresponding to ρ_{WSC_0.046}≈0.85. τ_e and τ_ρ are consistent. Furthermore, eMSEeMSE is smaller than that at 0° < θ < 90°, and ρ_WSE is larger than that at 0° < θ < 90°, indicating that the acoustic wave incidence angularity is within 90° < θ < 180°. Therefore, the time for the target echo to arrive at the receiving end of the active sonar is approximately 0.004 6 s0.004 6 s.

Aiming to verify the effectiveness of the time-delay estimation algorithm proposed in this paper, this algorithm is compared with the MF method. The relative errors of the time delay estimation are discussed under different interference sources, including reverberation and noise. The calculation of the relative error ε is shown in Eq. (25).

where τ_estimate is the estimated delay, and τrτr is the true delay, Monte Carlo 500 times and observe the results.

As shown in Figure 10, the "blue squares" indicate the delay estimation performance of the FrFT–FSWT joint time–frequency analysis method using mean square error as the evaluation parameter. The "red circles" indicate the performance of the same method evaluated by the waveform similarity coefficient. The "yellow asterisks" denote the delay estimation performance of the MF method. The "red circle" represents the delay estimation performance of the joint FrFT–FSWT time–frequency analysis method using the waveform similarity coefficient as the evaluation parameter, while the "yellow asterisk" denotes the performance of the MF method. As shown in Figure 10(a), the performance of MF deteriorates substantially when the SRR falls below -3 dB. Under an SRR = -10 dB, the relative error of the joint FrFT–FSWT algorithm is reduced by 66.7%–88.9% compared to MF. Similarly, Figure 10(b) shows that MF performance deteriorates when the SNR drops below -19 dB. At an SNR = -25 dB, the FrFT–FSWT joint algorithm achieves a 65.9%–91.5% lower relative error than MF. These results indicate that the delay estimation performance of MF deteriorates substantially under lower SRR/SNR conditions, whereas the FrFT–FSWT algorithm, particularly when evaluated using the mean square error, exhibits more stable and robust performance.

Figure  10  Proposed method performance under different SNR and SRR values

Download: Full-Size Img

6   Analysis of experimental data

Experiments conducted in a lake are performed to verify the effectiveness of the proposed time-delay estimation algorithm in a strong reverberation background. The transmitted signal parameters and target size information used in these tests are the same as those used in Section 4. The lake has a depth of 3 m. The target was towed to a designated location by a small sampan and buried in the lakebed. The sensor was suspended by a cable from a survey vessel, which was anchored to the shore with a mooring line. The transmitting array was connected to the receiving array to form a transceiver-combined sonar system. A schematic of the test instrumentation is shown in Figure 11, which includes a signal generator, power amplifier, data collector, transceiver-combined transducer, and personal computer. The array is highly directional, with a horizontal directivity opening angle of 6° and a vertical directivity opening angle of 30°. The main flap of the transmitting array is oriented toward the target to minimize reverberation intensity, given the strong reverberation in the echo signals when using active sonar to detect a submerged target.

Figure  11  Schematic of test apparatus

Download: Full-Size Img

As shown in Figure 12, the reverberation in the experimental data is strong, and the time frame in which the signal appears cannot be clearly observed. As shown in Figure 13, the delay value estimated by the MF is approximately 0.002 3 s. The delay estimation result using the FrFT–FSWT is shown in Figure 14, where the mean square error curve indicates that τ_e occurs at 0.055 s, with the acoustic incidence angle ranging between 0° and 90°. As shown in Figure 15, the highest point of the curve of the waveform τ_ρ similarity coefficient also occurs at 0.055 s. Compared with the real delay value of 0.056 6 s, the relative error is 2.83%. As shown in Figure 14, the mean square error remains low until approximately 0.04 s, and in Figure 15, the waveform similarity coefficient remains high until 0.04 s. This finding indicates that the reverberation is strong before 0.04 s. Therefore, under strong reverberation conditions, the MF estimate corresponds to the time when the reverberation first appears. Based on the time-delay estimation result τ_estimate=0.055 s from the FrFT–FSWT, the acoustic wave incidence angle is calculated to be approximately equal to 4.2°. The hemispherical head side of the underwater target is facing the sonar, and the straight-line distance from the sonar to the target is approximately 41.25 m. The relative position is shown in Figure 8(a).

Figure  12  Received signal in the time and time–frequency domains

Download: Full-Size Img

Figure  13  Matched filter delay estimation results

Download: Full-Size Img

Figure  14  Mean square error curve

Download: Full-Size Img

Figure  15  Waveform similarity coefficient curve

Download: Full-Size Img

7   Conclusions

In this paper, based on geometric acoustics theory and the parameters of a typical bright spot model, a bright spot parameter model for the experimental target is deduced.

Compared with the WT, STFT, and WVD algorithms, the FSWT algorithm is selected for its superior performance. FSWT is used to reconstruct multiple geometrical target echo signals in the time–frequency domain.

The FSWT algorithm alone cannot accurately determine the κ-value for broadband signals. Thus, this paper introduces the FrFT algorithm to address this issue. By rotating the LFM pulse into the u-domain using a certain angle, the broadband signal is transformed into a single-frequency signal, allowing for the derivation of an expression for the corresponding κ-value. Simulation results further verify that the FrFT–FSWT algorithm exhibits better anti-interference performance compared to the standalone FSWT algorithm.

In this paper, a time–frequency observation matrix for different incident angles is constructed based on the known geometric acoustic scattering characteristics of the target. This matrix is combined with the FrFT–FSWT algorithm to suppress reverberation and noise in the received signal. The signal is reconstructed using IFSWT from the de-temporalized time–frequency observation matrix. The reconstructed signals are compared with the theoretical highlight model under different delay conditions using the mean square error and waveform similarity coefficients. The optimal time delay estimation is determined based on the results.

Simulation analysis and experimental data show that, compared to the MF method, the algorithm proposed in this paper exhibits superior performance in time delay estimation. This improvement offers a novel approach for estimating the azimuth of shallow-water submerged targets.