Off-Grid Underwater Acoustic Source Direction-of-Arrival Estimation Method Based on Iterative Empirical Mode Decomposition Interval Threshold

To solve the problem that the hydrophone arrays are disturbed by ocean noise when collecting signals in shallow seas, resulting in reduced accuracy and resolution of target orientation estimation, a direction-of-arrival (DOA) estimation algorithm based on iterative EMD interval thresholding (EMD-IIT) and off-grid sparse Bayesian learning is proposed. Firstly, the noisy signal acquired by the hydrophone array is denoised by the EMD-IIT algorithm. Secondly, the singular value decomposition is performed on the denoised signal, and then an off-grid sparse reconstruction model is established. Finally, the maximum a posteriori probability of the target signal is obtained by the Bayesian learning algorithm, and the DOA estimate of the target is derived to achieve the orientation estimation of the target. Simulation analysis and sea trial data results show that the algorithm achieves a resolution probability of 100% at an azimuthal separation of 8° between adjacent signal sources. At a low signal-to-noise ratio of −9 dB, the resolution probability reaches 100%. Compared with the conventional MUSIC-like and OGSBI-SVD algorithms, this algorithm can effectively eliminate noise interference and provides better performance in terms of localization accuracy, algorithm runtime, and algorithm robustness.


Introduction
In recent years, with the increasing demand for marine resource development, environmental protection, and national defense, hydrophones have been widely used in various fields such as oil and gas exploration, subsea geological surveys, marine ecological monitoring, noise pollution control, and submarine detection [1].However, despite significant advancements in hydrophone technology, many challenges remain.The noise in marine environments is complex and variable, comprising both natural sources (such as waves, marine life, and rain) and artificial sources (such as ships and military sonar) [2].The non-stationary and complex characteristics of this noise greatly increase the difficulty of signal processing.Furthermore, marine noise often consists of a mixture of multiple noise sources, forming a complex noise field that further complicates the extraction and identification of target signals.Against this backdrop, it is crucial to develop effective denoising algorithms and methods to improve the positioning accuracy of hydrophone signals and address the challenges posed by complex underwater environments and noise interference [3].
In the field of underwater acoustic DOA estimation, hydrophone arrays have been widely employed for source localization [4].Traditional spatial spectrum estimation algorithms include Conventional Beamforming (CBF) [5,6], Minimum Variance Distortionless Response (MVDR), and subspace algorithms with super-resolution capabilities such as Sensors 2024, 24, 5835 2 of 26 MUSIC and ESPRIT [7,8].However, in unfavorable environments such as low signal-tonoise ratios, the performance of these algorithms in estimating the target orientation can be significantly reduced.In order to enhance the precision of location determination, a number of sparse signal processing algorithms, including Bayesian inference (SBI), have been put forth as potential solutions [9].Nevertheless, these methods rely on predefined spatial sampling grids.Despite the high grid density, the fact that actual signal sources do not precisely align with these grid points limits the angle resolution.To address this issue, gridless algorithms and algorithms that incorporate off-grid errors have been widely studied.The authors of [10] proposed a gridless convex optimization problem to address this issue by recovering the covariance matrix of a virtual array.The authors of [11] introduced a perturbed Sparse Signal Representation (SSR) model, incorporating bias parameters into the DOA estimation framework to address these challenges.The authors of [12] tackled the angle resolution limitation caused by sources not precisely aligning with grid points by constructing and optimizing a co-array tensor embedded with displacement information.The authors of [13] introduced the off-grid sparse Bayesian inference (OGSBI) algorithm, which incorporates off-grid errors into the array's manifold matrix, thereby avoiding the substantial computational burden associated with mesh refinement and overcoming the angle resolution limitations.
The non-stationarity of underwater acoustic signals coupled with the presence of interfering noises renders traditional stationary signal processing methods inadequate.To address this, algorithms such as empirical mode decomposition (EMD) and wavelet transform are utilized due to their proficient handling of non-stationary signals.In 1998, Hilbert Huang et al. proposed an adaptive signal processing method, empirical mode decomposition (EMD), for the analysis of non-linear non-stationary time series, which is based on the characteristic time scale of the signal itself.The method does not require the setting of a basis function.The intrinsic modal function is generated adaptively based on the analyzed signal [14].In the studies by [15,16], ensemble empirical mode decomposition (EEMD) and modified ensemble EMD (MEEMD) algorithms have been proposed to reduce the aliasing problem in empirical modal decomposition by adding noise to the original data and also to extend it to the field of high-dimensional processing.In the study by [17], an iterative EMD interval thresholding denoising algorithm is proposed, which is based on the idea of EMD with translation-invariant wavelet thresholding, and builds different noisy versions of the original signal by changing the position of the first intrinsic mode function (IMF) sample through multiple iterations, and then obtains multiple denoised versions of the signal by the IMF thresholding method and averaging them to enhance the denoising ability.
Specifically, the contributions of our work can be summarized as follows: a.
The algorithm leverages EMD for the adaptive decomposition of non-stationary signals into their respective IMFs.To overcome EMD's limitations in handling complex noise, IIT is applied to iteratively threshold the IMFs, effectively filtering out low-energy noise.This combined approach significantly enhances noise suppression, ensuring accurate DOA estimation even under low signal-to-noise ratio (SNR) conditions.b.
The integration with OGSBI further optimizes signal reconstruction, enhancing target resolution and addressing grid mismatch issues, which is particularly effective when dealing with closely spaced signal sources.c.
The use of singular value decomposition (SVD) reduces computational complexity, making this method not only more accurate than traditional algorithms but also more efficient in real-time applications.
The paper is structured as follows: Section 2 introduces the background of EMD (empirical mode decomposition) and discusses the traditional DOA (direction of arrival) estimation method in detail.It also establishes the signal reception model and elaborates on the principles and processes of the EMDIIT-OGSBI algorithm.Section 3 compares the EMDIIT-OGSBI algorithm with the MUSIC algorithm and the OGSBI-SVD algorithm through simulation analysis and analyzes the superiority of the algorithm in this paper.Section 4 presents the source of the experimental data from the sea trial and uses it to validate the algorithm.The conclusions are given in Section 5.

Signal Processing Model for Hydrophone Arrays
Multiple hydrophone array elements arranged in space according to a certain geometric position constitute a base array, which has better directivity than a single array element.It is assumed that there are M uniformly deployed hydrophone arrays for hydroacoustic signal localization [18], with neighboring arrays spaced d = λ/2 and λ being the signal wavelength.Under the plane wave assumption [19], the hydrophone base array with an array streamlines of a(θ) ∈ C M×1 receives K far-field uncorrelated narrowband signals.As shown in Figure 1, the incident azimuth is , and [•] T denotes transposition, where θ K denotes the orientation of the Kth signal.The orientation of the hydrophone array is measured relative to the zero angle, which corresponds to the broadside angle of the array.At the instantaneous moment t, the base array output vector x(t) ∈ C M×1 is expressed as [20] x where

Iterative EMD Interval Thresholding Methodology
The EMD-IIT denoising algorithm combines EMD with a translation-invariant wave let thresholding method [21], and the section below describes the specific principles in volved in this denoising algorithm.If each source signal has T snapshots, the model for multiple snapshots can be written as where

Iterative EMD Interval Thresholding Methodology
The EMD-IIT denoising algorithm combines EMD with a translation-invariant wavelet thresholding method [21], and the section below describes the specific principles involved in this denoising algorithm.

Empirical Mode Decomposition
The basic idea of EMD [22] is that for all signals, the time series can be decomposed into a finite number of IMFs with different characteristic scales, each representing a characteristic oscillation in a time scale and with local orthogonality and adaptive properties.In the analysis of non-stationary signals, the interference of the same frequency components between different components over time can be effectively avoided.
In the signal model established above for array reception, where the array element received signal is X ∈ C M×T , we perform EMD of each array element signal separately to obtain IMFs and residuals for each row vector, as follows: . . .
where r mL (t) is the residual, which is a slowly varying function of the non-zero mean with only a few extreme values.
The EMD algorithm proceeds as follows, taking as an example the first array received signal x 1 (t) ∈ C 1×T of the array received signal as the input signal: (1) Find the local maxima and minima of a, obtain the sequence of maxima and minima, and interpolate the maxima and minima, respectively, to obtain the upper and lower envelopes of x 1 (t).The two envelopes are asymmetric to each other.(2) Average the maxima and minima envelopes, and obtain the envelope mean m 1 (t).
Subtract the envelope mean from the original signal to obtain the first component ) Use h 1 (t) as the input and repeat steps 1 and 2 to obtain h 2 (t), iterating continuously.
Use the standard deviation of the two adjacent decomposition components as the criterion for stopping the iteration, which is generally taken as 0.2 ≤ S ≤ 0.
The iterative flow chart of the EMD algorithm is shown in Figure 2.

Wavelet Threshold Denoising
Wavelet threshold denoising has the feature of suppressing the useless part of th signal and enhancing the useful part [23,24].The following describes the principle o wavelet threshold denoising: The following noise signal ( ) y t is given by ( ) ( ) ( ) y t x t n t σ = + (5

Wavelet Threshold Denoising
Wavelet threshold denoising has the feature of suppressing the useless part of the signal and enhancing the useful part [23,24].The following describes the principle of wavelet threshold denoising: The following noise signal y(t) is given by where x(t) is the original signal, n(t) is the white noise signal following a standard Gaussian distribution, and σ is the noise variance.The wavelet threshold denoising method starts with the selection of a suitable wavelet basis W, and the discrete wavelet transform (DWT) is obtained as c = Wy (6) where is the wavelet decomposition coefficient.Then, a threshold is set for quantization, i.e., λ = σC, where C is a constant.The basic principle of wavelet thresholding is to set all wavelet coefficients below that threshold to zero, and to keep those above that threshold directly or to process them accordingly.Currently, the most common thresholding functions are the hard and soft threshold values proposed by Donoho [25], which are defined as and Sensors 2024, 24, 5835 6 of 26 For the selection of the threshold value, the common threshold formula is λ = σ √ 2 ln T, which guarantees a high probability that all noisy components will have a low amplitude.Here, T is the length of the sampled signal, and σ is the standard deviation of the noisy signal, which is estimated as [26] After processing with hard or soft thresholding, the denoised signal is obtained by inverting the processed wavelet coefficients, as follows:

Implementation of Iterative EMD Interval Thresholding
The IMF threshold denoising method in this paper is related to the wavelet threshold denoising method.Wavelet denoising is the thresholding of wavelet components, whereas in the EMD case, IMF thresholding is performed on samples from each IMF.However, the fact is that the noise contained in each IMF is colored, i.e., each mode has a different energy in it.In this sense, EMD denoising is most closely related to wavelet denoising of signals corrupted by colored noise, where the threshold must be scale-dependent, and by adapting the threshold function to the nature of each IMF, threshold denoising of the IMF obtained by EMD can locally exclude the low-energy IMF parts disturbed by high noise.Based on the idea of wavelet thresholding, the IMF threshold is given by for hard thresholding and given by for soft thresholding.In both thresholding cases, we can see that the thresholds are also different for the different IMFs, where h I MF i (t) represents the ith thresholded IMF.The thresholds we study become multiples of the independent thresholds for each IMF, i.e., , where E k represents the energy of the kth IMF, which is given by where E 2 1 is the energy of the first IMF.According to [27], the values of parameters β and ρ are 0.719 and 2.01.
The reconstruction of the denoised signal is given by In particular, the addition of the M 1 and M 2 parameters gives us more flexibility in excluding noisy low-order IMFs and in selecting thresholds for higher-order IMFs.Regarding the parameter M 1 , it is calculated by conventional EMD denoising [28][29][30].The denoised signal must be reconstructed from IMFs of order J and higher; in other words, the energy of the noise is lower than the energy of the signal in IMFs of order J and higher.The best choice for the M 1 parameter is And the best choice for the parameter is For an independent IMF sample amplitude, it is not possible to infer whether it corresponds to a noisy or useful signal.However, it is possible to guess whether the signal within this interval is noise-dominated or signal-dominated, based on the single extreme value h I MF i (r (j) If a strong signal is present within this interval, the absolute value of the extreme value will exceed the threshold; conversely, if the signal is small, the absolute value of the extreme value will be below the threshold.Thus, the modified hard and soft thresholds, referred to as the EMD interval threshold (EMD-IT), are and where represents the sample values from interval t in the ith IMF.Based on the idea of translation-invariant wavelet thresholding [28], multiple denoised versions of the signal are obtained by iteration, and their denoising power is enhanced by averaging them.In the EMD case, different denoised versions of the noisy signal can only be obtained by thresholding different versions of the IMF.We know that under Gaussian white noise conditions, the first IMF is mainly noise; in other words, it contains more noise than the others.By changing the position of the first IMF sample and then adding the newly generated noise signal to the sum of the remaining IMFs, we can obtain a different noisy version of the original signal.In fact, when the first IMF contains only noise, the total noise variance of the newly generated noise signal is the same as that of the original noise signal.Figure 3 shows the flowchart of the algorithm referred to as iterative EMD interval thresholding (EMD-IIT), in the following steps: (1) EMD expansion of the initial noisy signal x(t).
(5) Performing EMD processing on the newly obtained noisy signal x (a) (t).(6) The denoised version of the original signal, x 1 (t) of x, is obtained by denoising the IMFs of the newly obtained noisy signal x (a) (t) by Formula ( 16) or Formula ( 17).(7) Iterating Q − 1 more times in steps 3-6 (typically, Q is set to 20) to obtain q denoised versions of x, i.e.,

The Algorithm for Iterative EMD Interval Thresholding and Off-Grid Sparse B Learning
The DOA estimation algorithm in this paper consists of the following p part is to perform EMD-IIT denoising on the noisy signal received by each a to obtain the signal with noise interference removed; the second part is to sp nal space by building an off-grid sparse grid model and to decompose the de into singular values to further reduce the sensitivity to noise; and the third p by Bayesian continuous iteration, updating the hyperparameters, and final state of convergence to obtain the DOA estimates.The overall framework of is shown in Figure 4.

The Algorithm for Iterative EMD Interval Thresholding and Off-Grid Sparse Bayesian Learning
The DOA estimation algorithm in this paper consists of the following parts: the first part is to perform EMD-IIT denoising on the noisy signal received by each array element to obtain the signal with noise interference removed; the second part is to specify the signal space by building an off-grid sparse grid model and to decompose the denoised signal into singular values to further reduce the sensitivity to noise; and the third part is to learn by Bayesian continuous iteration, updating the hyperparameters, and finally reaching a state of convergence to obtain the DOA estimates.The overall framework of the algorithm is shown in Figure 4.   for ease of calculation.We write the noise-bearing signal in the f Then, the EMD-IIT algorithm is used to denoise each row of the noi rately to obtain the denoised signal ˆˆˆ[ ] , , , ⋅⋅⋅ 1 2 m X = X X X .Then, the ne vector after removing the noise is where X is the vector signal of an M T × dimension.

Off-Grid Sparse Model and Singular Value Decomposition
The spatial angle range [ 2,2] π π with each point representing a potential incident direction, such as

EMD-IIT Denoising
According to Formula (3), we obtain the noise-bearing signal received by the hydrophone arrays, i.e., X = A(Θ s )S+N, where X ∈ C M×T is the signal of T snapshots.Here, with the EMD-IIT algorithm above, we denoise the noisy signal.Since each array is subjected to EMD-IIT, the signal matrix received by the mth array element is defined as X m ∈ C 1×T for ease of calculation.We write the noise-bearing signal in the following form: Then, the EMD-IIT algorithm is used to denoise each row of the noisy signal separately to obtain the denoised signal X=[ X1 , X2 , • • • , Xm ].Then, the new array signal vector after removing the noise is where X is the vector signal of an M × T dimension.

Off-Grid Sparse Model and Singular Value Decomposition
The spatial angle range [−π/2, π/2] is uniformly divided into N grid points, with each point representing a potential incident direction, such as Θ = θ 1 , θ 2 , • • • , θ N , and It is evident that the target orientation can be reinterpreted as an overcomplete sparse representation across the N divided grid points.In the off-grid sparse model, there is an issue where the target incidence direction does not align perfectly with the grid points, resulting in a mismatch, represented as generally a denser grid point can reduce the error, but this does not contain all possible incidence directions and increases the computational effort.To address this issue, an off-grid error is incorporated Sensors 2024, 24, 5835 10 of 26 into the array's prevalence matrix by performing a first-order Taylor expansion between two neighboring grid points, allowing the steering vector to be approximated as where θ n k denotes the nearest grid point to , the grid error can be expressed as According to Formula (20), the overcomplete array popularity matrix can be expressed as Considering the case where the background noise has been filtered out above, we have X = Φ Ŝ, where K ≤ T and Rank( X ≤ Rank( Ŝ ≤ K. where V 1 and V 2 are matrices consisting of the first K columns of V and the remaining T − K columns, respectively.By singular value decomposition, we obtain matrices XV = [ XSV XV 2 ] con- taining all signal information, where XSV = XV 1 ∈ C M×K , and where the first part XSV retains most of the signal information and is used in the signal reconstruction process later, and the second part is discarded.Let ŜSV = ŜV 1 ; then, we can represent the signal after singular value decomposition as XSV = Φ ŜSV (23) where the signal ŜSV still has joint sparsity [29].

Sparse Bayesian Inference
The Bayesian inference method is used for the estimation of hydroacoustic target orientation, and the optimal estimate can be obtained.Assuming that the noise signal obeys the complex Gaussian distribution, the following likelihood function of XSV is given by p XSV ŜSV ; α 0 = CN XSV Φ ŜSV , α −1 0 (24) where CN represents the complex Gaussian distribution, α 0 = σ −2 , where σ 2 is the noise variance and α 0 is usually unknown, which is assumed to obey the Gamma prior distribution.
Let the prior probability of the sparse signal ŜSV be where γ is a set of hyperparameters, γ = [γ 1 , γ 2 , • • • , γ N ] T , which represents the source signal power incident to the array in each direction and also affects the sparsity of the sparse signal.Here, Γ = diag(γ) represents the covariance matrix of the sparse signal.
Combining the prior information and the likelihood function, the joint probability density function is obtained as The posterior probability density function of ŜSV can be obtained from Bayesian inference as p ŜSV XSV , α 0 , γ, β = p( XSV| ŜSV ,α 0 ,β)p( ŜSV| γ) where the posterior mean µ and posterior covariance matrix Σ of the sparse signal are given separately by µ = α 0 ΣΦ H XSV (28) and where µ and Σ are functions with respect to the three hyperparameters α 0 , γ, and β.We use a maximum posterior probability criterion to maximize the probability p α 0 , γ, β XSV , which leads to the derivation of two parameters α 0 and γ.These two parameters are given by and where Ξ t ≜ µ(t)(µ(t)) H , c, d → 0 , and ρ is a positive constraint taking small values.
The grid error β determines the accuracy of the target orientation estimation.We can use the expectation maximization criterion to find the grid error so that the expectation E ln p XSV ŜSV , α 0 , β p(β) is maximized, which is equivalent to minimizing where C is a constant depending on β.P is a semi-positive definite matrix whose expression is given by The angle correction vector is obtained by the above derivation as [30] Next, we obtain the expression for β by taking the derivative of β by Formula (32) and setting it to zero as where β −n is β without the nth entry for a vector β.By means of the constraint we have In the above Bayesian inference, we maximize the posterior probability to find out the updated formula of the two hyperparameters α 0 and γ.Then, through these two hyperparameters, we use the expectation maximization criterion to find out the update formula of the angle correction vector hyperparameter β and obtain the final grid error β n .Finally, we initialize α 0 , γ, and β and keep iterating Formulas (30), (31), and (37) until convergence, so that we can calculate the estimated value of the K DOAs as The specific flow of the algorithm for iterative EMD interval thresholding and off-grid sparse Bayesian learning is as follows: (1) Pass the received signal from the array through EMD-IIT to obtain the denoised signal X. (2) Construct the off-grid sparse model to obtain the overcomplete sparse dictionary Φ and perform the singular value decomposition of X to obtain XSV .(3) Initialize the hyperparameters α 0 and γ for the noise and signal, respectively, and initialize the angle correction vector β, the mean µ, and the variance Σ to zero.(4) Use Formulas ( 28) and ( 29) to solve for the mean and variance, respectively.

Simulation Analysis
In order to verify the feasibility of the algorithms in this paper, simulation analysis is performed in this section and the estimated performance of the three algorithms is compared.If not specified, the following parameters and initial values were used in the DOA estimation of the objectives: In OGSBI-SVD and EMDIIT-OGSBI, ρ = 0.01 and , and β=0 are initialized.The uniform line array element M = 8 is initialized with a grid distance of r = 1, τ = 10 −3 is set, and the maximum number of iterations is 800.The original signal is a frequency band signal with a center frequency of f m = 1000 Hz, λ = v/ f m , the sampling frequency is 10 kHz, v = 1500 m/s is set as the underwater sound velocity, and the adjacent array element spacing is d = λ/2 = 0.75 m.The signal-to-noise ratio is calculated as SNR = 10lg(P s /P n ), where P s is the signal power, P n is the noise power, and the spatial angle is divided into [−π/2, π/2], and the noise is complex Gaussian white noise.The root mean square error (RMSE) is defined as RMSE = where S denotes the number of Monte Carlo experiments, K is the number of source signals, and θk s denotes the orientation estimate of the kth signal source in the first experiment.

EMD-IIT Denoising Analysis
To visualize the denoising performance of the EMD-IIT algorithm, the noisy signals received by each array element are compared and analyzed with the original signals in this section.The time-frequency spectrum of the original signal is shown in Figure 5.Let the incident directions of the two target signals be [−3.6 • 11 • ]; the number of snapshots is T = 1024.The time-frequency spectrum of the received signals of each array element at a signal-to-noise ratio of 5 dB is given in Figure 6. Figure 7 gives the time-frequency spectrum of the received signals of each array element after the EMD-IIT denoising algorithm.We can see that after the EMD-IIT denoising algorithm, a portion of the frequency interference has been effectively removed from the time-frequency spectrum of the array element.This is due to the enhanced noise reduction capability of the EMD-IIT algorithm through multiple noises averaging by the translation-invariant wavelet threshold.The signal-to-noise ratio after denoising and the signal-to-noise ratio before denoising are given in Table 1 for each array element, which further illustrates the good denoising performance of EMD-IIT.

Spatial Power Spectrum Estimation Analysis
To validate the proposed method and highlight its superiority, we compare it under identical conditions with several existing algorithms, including the MUSIC algorithm, the off-grid sparse Bayesian inference algorithm (OGSBI), and the MUSIC-like algorithm.This comparative analysis demonstrates the advantages of our algorithm in terms of performance and accuracy.The orientation of the two target signals is [−3.6 • 11 • ], SNR = 0 dB, and the number of snapshots is T = 1024.Figure 8 illustrates that the EMDIIT-OGSBI algorithm, as proposed in this paper, is capable of matching the target incident orientation, exhibiting higher spatial-spectral gain and a narrower main flap width.The OGSB-SVD algorithm is able to distinguish the target orientation; however, the main flap width is wider and comprises four pseudo-peaks.The MUSIC algorithm is capable of approximating the target orientation; however, the main flap width is wider, and the peak is not discernible.The MUSIC-like algorithm is able to align with the target orientation; however, its main lobe width is wider, and the spatial gain is inferior in comparison to the algorithm proposed in this paper.

Spatial Power Spectrum Estimation Analysis
To validate the proposed method and highlight its superiority, we compare it under identical conditions with several existing algorithms, including the MUSIC algorithm, the off-grid sparse Bayesian inference algorithm (OGSBI), and the MUSIC-like algorithm.This comparative analysis demonstrates the advantages of our algorithm in terms of performance and accuracy.The orientation of the two target signals is [−3.6°11°], SNR 0 dB = , and the number of snapshots is 1024 T = . Figure 8 illustrates that the EMDIIT-OGSBI algorithm, as proposed in this paper, is capable of matching the target incident orientation, exhibiting higher spatial-spectral gain and a narrower main flap width.The OGSB-SVD algorithm is able to distinguish the target orientation; however, the main flap width is wider and comprises four pseudo-peaks.The MUSIC algorithm is capable of approximating the target orientation; however, the main flap width is wider, and the peak is not discernible.The MUSIC-like algorithm is able to align with the target orientation; however, its main lobe width is wider, and the spatial gain is inferior in comparison to the algorithm proposed in this paper.

RMSE of the Algorithm at Different Numbers of Monte Carlo Trials
Figure 9 shows the variation in RMSE with the number of Monte Carlo trials.The two target orientations in this experiment are [−3.6°11°], the signal-to-noise ratio is 0 dB, the number of snapshots is 1024, and other parameters are unchanged.It can be seen from the graph that when the number of Monte Carlo trials is less than 100, the RMSE values of the four algorithms converge less, and the DOA estimates of the algorithms have a chance at this time.To avoid this problem, the number of Monte Carlo trials in the following root mean square error analysis is 200.

RMSE of the Algorithm at Different Signal-to-Noise Ratios
Figure 10 shows the variation in root mean square error with signal-to-noise ratios for the four algorithms.The number of snapshots is 1024, SNR = −10 : 1 : 10 dB, and the two target signal orientations are [−14.3• 6 • ], respectively.As shown in Figure 9, the MUSIC algorithm has the relatively highest RMSE value at low SNRs, and its orientation estimation performance is limited.The OGSBI-SVD algorithm and the MUSIC-like algorithm have higher RMSE values; however, their performance improves as the signal-to-noise ratio increases, leading to better orientation estimation under the influence of weak noise.In contrast, the EMDIIT-OGSBI algorithm in this paper has a relatively lower RMSE value and a relatively more outstanding DOA estimation accuracy and still has a high DOA estimation accuracy at low signal-to-noise ratios.This is due to the suppression of Gaussian noise by the EMDIIT denoising algorithm, and the singular value decomposition of the off-grid sparse observation matrix also reduces the effect of noise interference.], the number of snapshots is 1024, the signal-to-noise ratio range is SNR = −7 : 1 : 10 dB, and other conditions remain unchanged.As can be seen from Figure 12, the RMSE values of the four different grid distances gradually decrease as the signal-to-noise ratio increases.The finer the grid distance is divided, the smaller the RMSE value is.We can also see that the difference in RMSE between coarse and fine grid distance is insignificant, which shows that the proposed algorithm still has high accuracy in estimating the target orientation under coarse grid distance.It is worth noting that the algorithm in this paper can still maintain a high estimation accuracy even with coarse grid spacing., and other conditions remain unchanged.As can be seen from Figure 12, the RMSE values of the four different grid distances gradually decrease as the signalto-noise ratio increases.The finer the grid distance is divided, the smaller the RMSE value is.We can also see that the difference in RMSE between coarse and fine grid distance is insignificant, which shows that the proposed algorithm still has high accuracy in estimating the target orientation under coarse grid distance.It is worth noting that the algorithm in this paper can still maintain a high estimation accuracy even with coarse grid spacing.

Analysis of the Discriminative Probability of Compact Sound Sources
This section focuses on the analysis of the spatial resolution probability of compact sound sources.Under a uniform line array, the spatial resolution capability of the algorithm in this paper for compact sound sources is analyzed.We set the incident directions of the two target signals as [2.4 8.5 ] °°, the signal-to-noise ratio as 5 dB, the number of snapshots as 1024, and the grid distance as 1°, with other conditions being constant.Figure

Analysis of the Discriminative Probability of Compact Sound Sources
This section focuses on the analysis of the spatial resolution probability of compact sound sources.Under a uniform line array, the spatial resolution capability of the algorithm in this paper for compact sound sources is analyzed.We set the incident directions of the two target signals as [2.4 • 8.5 • ], the signal-to-noise ratio as 5 dB, the number of snapshots as 1024, and the grid distance as 1 • , with other conditions being constant.Figure 13 gives the spatial-spectral estimation plots of the four algorithms for compact sound sources.From Figure 13, it can be seen that compared with the MUSIC, MUSIC-like, and OGSBI-SVD algorithms, the algorithms in this paper have a better spatial resolution for the compact sound sources and have narrower main flap widths and sharper peaks, indicating that the orientation estimation of the compact sound sources also has higher accuracy.To further investigate the discrimination ability of the algorithm proposed in this paper for spatially tight signal DOA estimation, the discrimination probabilities of the three algorithms are analyzed at different DOA intervals with a signal-to-noise ratio of 0 dB with other conditions being constant.The DOAs of the two target signals are defined as satisfied, the discrimination of the compact sound sources is successful.As can be seen from Figure 14, the EMDIIT-OGSBI algorithm in this paper achieves a 100% resolution probability at source intervals ∆θ ≥ 8 • .The OGSBI-SVD algorithm achieves 100% resolution probability at ∆θ = 22 • ; the MUSIC algorithm achieves 100% resolution probability at ∆θ = 26 • ; and the MUSIC-like algorithm achieves 100% resolution probability at ∆θ = 12 • .It can be seen that the algorithm has a good spatial resolution of spatially immediate neighboring signals.

Analysis of the Discriminative Probability at Different Signal-to-Noise Ratios
To investigate the algorithms' ability to discriminate between targets at low signalto-noise ratios, the following simulations are made.The incident directions of the two target signals are randomly generated in the range of [−90 • 90 • ], the separation between them is 14 • , and the signal-to-noise ratio range is SNR = −10 : 1 : 10 dB, with all other conditions being equal.Figure 15 shows the resolution probabilities of the four algorithms as the signal-to-noise ratio varies.From this figure, it can be seen that the MUSIC and OGSBI-SVD algorithms have a poor resolution of the target at a low S/N ratio, and their resolution ability is improving as the S/N ratio increases.Compared to the MUSIC and OGSBI-SVD algorithms, the MUSIC-like algorithm has better target resolution capability at low signal-to-noise ratios, with its resolution probability reaching 100% at a low SNR of −4 dB.In contrast to the other three algorithms, the proposed EMDIIT-OGSBI algorithm still has a very strong resolution of the target at low SNRs, and the resolution probability reaches 100% only at a low SNR of −9 dB, which effectively shows that the algorithm still has an excellent resolution of the target signal at low SNRs.
at low signal-to-noise ratios, with its resolution probability reaching 100% at a low SNR of −4 dB.In contrast to the other three algorithms, the proposed EMDIIT-OGSBI algorithm still has a very strong resolution of the target at low SNRs, and the resolution probability reaches 100% only at a low SNR of −9 dB, which effectively shows that the algorithm stil has an excellent resolution of the target signal at low SNRs.

Algorithm Runtime Analysis
Figure 16 shows the comparison of the algorithm running time of the four algorithms with different grid spacing.SNR = 5 dB was set, the number of snaps is 1024, and the two target incidence directions are [ 13.4 5.6 ] − °°, respectively.We can see that the running time of the four algorithms decreases as the grid spacing increases.The MUSIC algorithm has the shortest runtime because it does not involve iterative operations; however, this also results in its poor orientation estimation performance, as discussed earlier.The MU-SIC-like algorithm requires the calculation of fourth-order cumulants, which increases its runtime.The EMDIIT-OGSBI algorithm includes iterative operations not only in the noise reduction process but also in the azimuth angle estimation, leading to a longer runtime than the MUSIC algorithm.Nevertheless, its runtime is still shorter than that of the OGSBI-SVD algorithm.This is because, in the EMDIIT-OGSBI algorithm, the EMDIIT method effectively reduces the impact of background noise, significantly decreasing the number of iterations required by OGSBI.Additionally, singular value decomposition is employed to reduce the signal's dimensionality and discard unnecessary signal matrices further reducing the runtime.

Algorithm Runtime Analysis
Figure 16 shows the comparison of the algorithm running time of the four algorithms with different grid spacing.SNR = 5 dB was set, the number of snaps is 1024, and the two target incidence directions are [−13.4• 5.6 • ], respectively.We can see that the running time of the four algorithms decreases as the grid spacing increases.The MUSIC algorithm has the shortest runtime because it does not involve iterative operations; however, this also results in its poor orientation estimation performance, as discussed earlier.The MUSIC-like algorithm requires the calculation of fourth-order cumulants, which increases its runtime.The EMDIIT-OGSBI algorithm includes iterative operations not only in the noise reduction process but also in the azimuth angle estimation, leading to a longer runtime than the MUSIC algorithm.Nevertheless, its runtime is still shorter than that of the OGSBI-SVD algorithm.This is because, in the EMDIIT-OGSBI algorithm, the EMDIIT method effectively reduces the impact of background noise, significantly decreasing the number of iterations required by OGSBI.Additionally, singular value decomposition is employed to reduce the signal's dimensionality and discard unnecessary signal matrices, further reducing the runtime.

Validating Algorithms with Sea Trial Data
The data obtained from a sea trial experiment conducted in a sea area are used t verify the performance of the algorithm.The experimental sea area on that day was calm in terms of ocean noise, with no other passing vessels on the surface and low wind speed The sound velocity profiles of the sea surface on that day are shown in Figure 17, whic were measured at 13:52 pm, 14:57 pm, and 16:41 pm.In the sea trial, the sound source tran mitting equipment UW350 was lifted at a depth of 5 m and the transmitting signal was

Validating Algorithms with Sea Trial Data
The data obtained from a sea trial experiment conducted in a sea area are used to verify the performance of the algorithm.The experimental sea area on that day was calm in terms of ocean noise, with no other passing vessels on the surface and low wind speeds.The sound velocity profiles of the sea surface on that day are shown in Figure 17, which were measured at 13:52 p.m., 14:57 p.m., and 16:41 p.m.In the sea trial, the sound source transmitting equipment UW350 was lifted at a depth of 5 m and the transmitting signal was a broadband long pulse signal of 200-600 Hz, with a signal length of 15 s and a sampling frequency of 10 kHz.The transmitted signal is shown in Figure 18.The time-frequency spectrum of the sound source is shown in Figure 19.The number of hydrophone array elements is eight, the array element spacing is half a wavelength, and the experimental sea depth is 25.5 m.The above equipment deployment depth and seawater depth are measured by the depth sensor.The schematic layout of the sea trial is shown in Figure 20.

Validating Algorithms with Sea Trial Data
The data obtained from a sea trial experiment conducted in a sea area are used verify the performance of the algorithm.The experimental sea area on that day was c in terms of ocean noise, with no other passing vessels on the surface and low wind spee The sound velocity profiles of the sea surface on that day are shown in Figure 17, wh were measured at 13:52 pm, 14:57 pm, and 16:41 pm.In the sea trial, the sound source tra mitting equipment UW350 was lifted at a depth of 5 m and the transmitting signal w broadband long pulse signal of 200-600 Hz, with a signal length of 15 s and a samp frequency of 10 kHz.The transmitted signal is shown in Figure 18.The time-frequency sp trum of the sound source is shown in Figure 19.The number of hydrophone array eleme is eight, the array element spacing is half a wavelength, and the experimental sea dept 25.5 m.The above equipment deployment depth and seawater depth are measured by depth sensor.The schematic layout of the sea trial is shown in Figure 20.It can be seen that the noise of the day interferes with the original signal.The results of DOA spatial spectrum estimation at different target orientations are given in Figures 23 and 24, respectively.In the orientation estimation of the sea trial data, the number of snapshots chosen is 512 and 1024.From these figures, it can be seen that the MUSIC algorithm can roughly distinguish the target orientation, has some deviation in the target orientation estimation, has a wide main flap width and poor spatial power spectrum gain, and the OGSBI-SVD algorithm can approximate the target orientation, with a relatively narrow width of the main flap and a more pronounced peak.However, there are pseudo-peaks in the spatial power spectrum, which will seriously affect the discrimination of the true orientation; the MUSIC-like algorithm is relatively stable and is capable of distinguishing target direction angles with a narrower main lobe width.However, it exhibits poor spatial gain, whereas the EMDIIT-OGSBI algorithm in this paper has a narrower main flap width and higher spatial-spectral gain, which can estimate the target orientation more accurately.In summary, the DOA estimates of the four algorithms are approximately equal to the actual deployed orientation.In comparison, both the EMDIIT-OGSBI algorithm and the MUSIC-like algorithm demonstrate high accuracy at low snapshot counts, are less affected by the number of snapshots, and exhibit better robustness.
the number of snapshots chosen is 512 and 1024.From these figures, it can be seen that the MUSIC algorithm can roughly distinguish the target orientation, has some deviation in the target orientation estimation, has a wide main flap width and poor spatial power spectrum gain, and the OGSBI-SVD algorithm can approximate the target orientation, with a relatively narrow width of the main flap and a more pronounced peak.However, there are pseudo-peaks in the spatial power spectrum, which will seriously affect the discrimination of the true orientation; the MUSIC-like algorithm is relatively stable and is capable of distinguishing target direction angles with a narrower main lobe width.However, it exhibits poor spatial gain, whereas the EMDIIT-OGSBI algorithm in this paper has a narrower main flap width and higher spatial-spectral gain, which can estimate the target orientation more accurately.In summary, the DOA estimates of the four algorithms are approximately equal to the actual deployed orientation.In comparison, both the EMDIIT-OGSBI algorithm and the MUSIC-like algorithm demonstrate high accuracy at low snapshot counts, are less affected by the number of snapshots, and exhibit better robustness.To better compare the performance of the four algorithms for the orientation estimation of the sea trial data, 200 Monte Carlo experiments are conducted on the sea trial experimental data of the three positions, and the number of snapshots is chosen as 1024, and the mean value and root mean square error of the DOA estimation results are recorded in Table 2.The estimation results of the algorithm in this paper are closer to the target bearing To better compare the performance of the four algorithms for the orientation estimation of the sea trial data, 200 Monte Carlo experiments are conducted on the sea trial experimental data of the three positions, and the number of snapshots is chosen as 1024, and the mean value and root mean square error of the DOA estimation results are recorded in Table 2.The estimation results of the algorithm in this paper are closer to the target bearing and have a stronger suppression ability to the background noise interference of the ocean, which is generally consistent with the numerical simulation analysis.

Conclusions
This paper investigates the application of EMD-IIT-based denoising with off-grid sparse Bayesian learning algorithms for hydroacoustic direction-of-arrival (DOA) estimation.The proposed algorithm effectively reduces noise, improves target resolution, and enhances estimation accuracy.
The EMDIIT-OGSBI algorithm demonstrates strong robustness, shorter runtime, and reduced sensitivity to grid spacing and the number of snapshots compared to conventional algorithms.It maintains high resolution for target signals and neighboring signals in the presence of strong noise interference.Overall, the algorithm addresses the challenges of low precision and poor resolution in target DOA estimation under low signal-to-noise ratios, making it a relevant and beneficial approach in hydroacoustic signal processing.

2 Figure 1 .
Figure 1.Model diagram of array received signal.

Figure 1 .
Figure 1.Model diagram of array received signal.

( 3 )
Randomly changing the sample position of the first IMF, IMF

Figure 4 .
Figure 4. Flowchart of algorithm for iterative EMD interval thresholding and off-g ian learning.

Figure 4 .
Figure 4. Flowchart of algorithm for iterative EMD interval thresholding and off-grid sparse Bayesian learning.

Figure 6 .Figure 7 .
Figure 6.The time-frequency spectrum of each array when the noisy signal is received.(a) The timefrequency spectrum of the first to the fourth array when the noisy signal is received.(b) The timefrequency spectrum of the fifth to the eighth array when the noisy signal is received.

Figure 7 .
Figure 7.The time-frequency spectrum of each array after EMD-IIT denoising.(a) The time-frequency spectrum of the first to the fourth array after EMD-IIT denoising.(b) The time-frequency spectrum of the fifth to the eighth array after EMD-IIT denoising.

Figure 8 .
Figure 8.The spatial power spectrum of three algorithms.

Figure 8 .
Figure 8.The spatial power spectrum of three algorithms.

3. 3 . 28 Figure 9 .
Figure 9 shows the variation in RMSE with the number of Monte Carlo trials.The two target orientations in this experiment are [−3.6 • 11 • ], the signal-to-noise ratio is 0 dB, the number of snapshots is 1024, and other parameters are unchanged.It can be seen from the graph that when the number of Monte Carlo trials is less than 100, the RMSE values of the four algorithms converge less, and the DOA estimates of the algorithms have a chance at this time.To avoid this problem, the number of Monte Carlo trials in the following root mean square error analysis is 200.Sensors 2024, 24, x FOR PEER REVIEW 18 of 28

Figure 10 . 28 Figure 11 .
Figure 10.RMSE vs. signal-to-noise ratios.3.3.3.RMSE of Algorithms at Different Snap Counts Figure 11 shows a comparison of the RMSEs of the four algorithms for different numbers of snapshots.The range of the number of snapshots is T = [32, 64, 128, 256, 512, 1024], the signal-to-noise ratio is −3 dB, the two target orientations are [−17.2• 2.4 • ], and other conditions remain unchanged.As can be seen from Figure 11, the RMSE values of the four algorithms gradually decrease as the number of snapshots increases, while the RMSE values of the proposed algorithm are relatively lower, and the target orientation can be accurately estimated under the low snapshot condition.Sensors 2024, 24, x FOR PEER REVIEW 19 of 28

Figure 12 .
Figure 12.Variation plots of RMSE with signal-to-noise ratios under different grid distances.

Figure 12 .
Figure 12.Variation plots of RMSE with signal-to-noise ratios under different grid distances.

θ 1 =
u • and θ 2 = (u + ∆θ) • , and their corresponding spatial power spectrum values are ζ 1 and ζ 2 ; the DOA interval ∆θ was varied from 2 • to 30 • in steps of 2 • , and 200 Monte Carlo experiments are performed at each DOA interval.The intermediate values of the two target DOAs are set to θ µ = (θ 1 + θ 2 )/2, and θ µ corresponds to a spatial power spectrum value of ζ

Figure 13 .
Figure 13.Spatial power spectrum of compact sound source.

Figure 17 .
Figure 17.The profile of sound speed.

Figure 19 .
Figure 19.The time-frequency spectrum of the sound source.

Figure 19 .
Figure 19.The time-frequency spectrum of the sound source.

Figure 19 .
Figure 19.The time-frequency spectrum of the sound source.

Figure 20 .
Figure 20.Sea trial deployment diagram.DOA verification of the target was carried out by selecting data from two different locations from the sea trials.The target orientations were [−8.2 • 16.7• ], which were placed at 14:57 p.m. and 16:41 p.m., respectively.Figures21 and 22show the signal's timefrequency spectrum of the eight array elements collected at two time periods, respectively.It can be seen that the noise of the day interferes with the original signal.The results of DOA spatial spectrum estimation at different target orientations are given in Figures23 and 24, respectively.In the orientation estimation of the sea trial data, the number of snapshots chosen is 512 and 1024.From these figures, it can be seen that the MUSIC algorithm can roughly distinguish the target orientation, has some deviation in the target orientation estimation, has a wide main flap width and poor spatial power spectrum gain, and the OGSBI-SVD algorithm can approximate the target orientation, with a relatively narrow width of the main flap and a more pronounced peak.However, there are pseudo-peaks in the spatial power spectrum, which will seriously affect the discrimination of the true orientation; the MUSIC-like algorithm is relatively stable and is capable of distinguishing

Figure 21 .
Figure 21.The time-frequency spectrum of each array at 14:57 pm.(a) The time-frequency spectrum of the first to the fourth array.(b) The time-frequency spectrum of the fifth to the eighth array.

Figure 21 .Figure 21 .
Figure 21.The time-frequency spectrum of each array at 14:57 p.m. (a) The time-frequency spectrum of the first to the fourth array.(b) The time-frequency spectrum of the fifth to the eighth array.

Figure 22 .Figure 23 .
Figure 22.The time-frequency spectrum of each array at 16:41 pm.(a) The time-frequency spectrum of the first to the fourth array.(b) The time-frequency spectrum of the fifth to the eighth array.

Figure 22 .Figure 22 .Figure 23 .
Figure 22.The time-frequency spectrum of each array at 16:41 p.m. (a) The time-frequency spectrum of the to the fourth array.(b) The time-frequency spectrum of the fifth to the eighth array.

Figure 23 .
Figure 23.Estimation of the spatial power spectrum at the second position.(a) The snapshot count is 512.(b) The snapshot count is 1024.

Figure 23 .Figure 24 .
Figure 23.Estimation of the spatial power spectrum at the second position.(a) The snapshot count is 512.(b) The snapshot count is 1024.

Figure 24 .
Figure 24.Estimation of the spatial power spectrum at the third position.(a) The snapshot count is 512.(b) The snapshot count is 1024.

Author Contributions:
Conceptualization, C.X. and G.T.; methodology, C.X.; software, S.D.; validation, S.D.; formal analysis, C.X.; resources, S.D.; data curation, S.D.; writing-original draft preparation, G.T.; writing-review and editing, G.T.; visualization, S.D.; supervision, C.X.; project administration, C.X.; funding acquisition, C.X.All authors have read and agreed to the published version of the manuscript.Funding: The authors would like to acknowledge the National Natural Science Foundation of China (Grant No. 61761048) and the Basic Research Special General project of Yunnan Province, China (Grant No. 202101AT070132).Institutional Review Board Statement: Not applicable.Informed Consent Statement: Not applicable.
(t) is the first IMF, noted as I MF 11 (t), and subtract the input signal from I MF 11 (t) to obtain the residual quantity r 11 (t) = x 1 (t) − I MF 11 (t).(5)User11(t) as the new input signal and re-execute steps 1 to 5 to obtain a new residual r 12 (t) and a second I MF 12 (t), and so on, to obtain I MF 13 (t), I MF 14 (t), • • • , I MF 1L (t).At this point, the residual r 1L (t) becomes a monotonic function and cannot be decomposed into IMF again.The final result is x 1 (t) =The iterative flow chart of the EMD algorithm is shown in Figure2.

Table 1 .
SNR of each array before and after denoising.

Table 1 .
SNR of each array before and after denoising.

Table 2 .
Mean and root mean square errors of DOA estimate.