Inﬁnite Weighted p -Norm Sparse Iterative DOA Estimation via Acoustic Vector Sensor Array under Impulsive Noise

: Recently, many direction-of-arrival (DOA) estimation techniques based on sparse representation have been proposed. However, these techniques often suffer from performance degradation issues in the presence of impulsive noise. This paper aims to overcome this challenge in conventional sparse-based techniques on an acoustic vector sensor array (AVSA). Firstly, to remove high outliers from the array output data, the output information of the AVSA is weighted by using the inﬁnite norm. To further suppress outliers, a p -order cost function is formulated by extending the Frobenius norm to lower order, and then the expression of the signal power is quantiﬁed. Lastly, the DOA is approximated on the signal power by a spectral peak search mechanism. DOA estimation results based on Monte Carlo simulations validate the accuracy and robustness of the proposed techniques herein compared to the current, available methods in the context of intense impulsive noise, low generalized signal–to–noise ratio (GSNR), and a smaller number of snapshots.


Introduction
The estimation of direction of arrival (DOA) using an acoustic vector sensor (AVS) has been crucial in the fields of underwater environment detection, underwater target tracking, underwater target feature extraction, underwater acoustic communication, ocean energy exploration, and seabed geoacoustic parameter estimation [1][2][3][4][5][6][7][8][9].Unlike the classic acoustic scalar sensor (ASS), that can simply measure acoustic pressure information, a single AVS has the ability to assess acoustic element velocity and acoustic pressure information in two or three quadratic orientations at a single point in the field [10].Nehorai and Paldi demonstrated in [11] that a single AVS can attain DOA without left-right ambiguity.Compared to the acoustic scalar sensor array (ASSA), the AVS array (AVSA) can utilize more available acoustic data to improve DOA estimation performance.Consequently, AVS and AVSA are extensively used for DOA estimation.Some well-known examples of subspace-based algorithms comprise Multiple Signal Classification (MUSIC) [12] and Estimated Signal Parameter via Rotational Invariance Techniques (ESPRIT) [13].In addition, numerous improved algorithms have been proposed in [14][15][16][17].Despite the fact that these subspace-based methods can achieve accurate DOA estimation by manipulating the orthogonality among the signal subspace and the noise subspace, they have a number of obvious drawbacks.For instance, they need to be aware of the number of sources.Unfortunately, it is usually a challenging task to obtain the number of sources in practice.Furthermore, these subspace-based methods are ineffective in the existence of small snapshots or associated signals.Recently, the DOA estimation technology based on sparse signal reconstruction has attracted considerable attention due to its high resolution and robustness towards insufficient snapshots or correlated sources.The L1SVD method proposed in [18] is one of the most well-known sparse signal reconstruction methods.However, it has to be aware of the number of transmitters beforehand, and its DOA computation precision is only good in conditions with high signal-to-noise ratios (SNR).In complex underwater environments, selecting a suitable normalizing limit is a difficult task to degrade the DOA estimation performance of the L1SVD method.As stated in [19], the acoustic pressure and the acoustic particle velocity of the noise measured from AVS are uncorrelated when the ambient noise is isotropic.Based on this fact, the augmented cross-covariance vectorization DOA estimation techniques based on sparse representation has been proposed in [20,21].The above sparse representation-based methods have the advantages of higher angular resolution and estimation accuracy compared to traditional high-resolution methods.However, they are subject to the convex optimization problem, which has a significant computational burden.A user parameter-free and nonparametric method, called the iterative adaptive approach (IAA), was introduced in [22] to perform DOA estimation in order to figure out this issue, and high-precision DOA estimation performance was demonstrated even with a small sample size of snapshots.However, its DOA estimation performance suffers significantly when analyzing signals with a small angular interval.Furthermore, some sparse iterative DOA estimation based on AVSA have been developed.Ref. [23] proposed a weighted power compensation sparse iterative algorithm via AVSA, which used the current solution to compensate for the next iteration, and the accuracy was effectively improved.In addition, by compensating for the signal with l 2 -norm to improve the sparsity of signals in the spatial domain, the weighted l 2 -norm based IAA was proposed in [24] to provide accurate DOA.However, these sparse iterative algorithms are based on the assumption of Gaussian noise.In complex underwater acoustic environments, their DOA estimation performance will be severely degraded when dealing with data-containing outliers.
It is important to note that the ambient noise may be of different types and reveal non-Gaussian attributes in real-world scenarios such as wireless radio systems and underwater acoustic communication.Impulsive noise and burst noise are two common types of non-Gaussian noise [25,26].The probability density function (PDF) of impulsive noise has larger tails compared to the Gaussian distribution [27].When compared to the Gaussian distribution, heavy-tailed distributions have a greater probability of offering results with multiple standard deviations.These significant values can be regarded as outliers since they are highly improbable to occur within the standard Gaussian noise model.Several frameworks have been developed to model the impulse property, the most common of which are the Gaussian mixture model (GMM) [28], generalized Gaussian distribution (GGD) [29], and α-stable distribution [30].The GMM consists of two or more Gaussian distributions and has a waveform that exhibits randomness and Gaussian characteristics.It is widely used to deal with large-scale data sets and uncertainty problems.GGD is a generalization of Cauchy distribution that can describe probability distributions with heavier tails.It has the advantages of controllable noise shape and non-negative value and broad applicability in digital image processing, communication, signal processing, and other fields.α-stable distribution, which is famous for its heavy-tailed distribution and non-Gaussian model, is especially suitable for dealing with data involving extreme outliers.Thus, standard symmetric α-stable (SαS) distribution may be the most appropriate to describe impulsive noise [31,32].
When there are outliers in the data output of the ASS or AVS array, it may cause problems such as signal distortion, noise enhancement, and signal detection difficulties.To overcome outliers, various subspace-based DOA approximation techniques following fractional lower-order statistics (FLOS) exist in the literature.For example, manipulation of impulsive noise as a complex symmetric α-stable, ref. [33] has been used to formulate a robust covariation model, which is identical to the covariance matrix of Gaussian distribution.Then, the robust covariation matrix was used to improve DOA estimation performance under impulsive noise.Applying the fractional lower-order moment (FLOM) to the problems where the signal is circular and the additive noise is SαS distribution, a FLOM-MUSIC algorithm was formulated in [34].However, the DOA approximation performance of FLOM-MUSIC can be improved in the existence of impulsive noise when α falls within the range of 1-2.When α is less than 1, its performance will suffer.To solve this issue, a new subspace algorithm involving phased fractional lower-order moment (PFLOM) was proposed in [35], showing a higher resolution capability within the range of 0-2.The weighted array data-based MUSIC (WARD-MUSIC) algorithm is proposed in [36] by weighting the output data snapshot by snapshot.However, in order to achieve a satisfactory performance, the FLOS subspace-based methods necessitate a large sample size.To address this issue, the FLOS vector sparse representation method was proposed in [37] by combining sparse signal representation technology and the FLOS theory.Formed on the hypothesis that both parts of noise measurement are independent and uniformly dispersed random variables, then the original complex-valued measurement model is transformed and represented as an augmented real-valued measurement.An effective method involving sparse mechanism is proposed in [38] to handle the DOA approximation challenge under impulsive noise.Since the outliers were sparse, the noise matrix was separated into outliers and complex Gaussian noise matrix terms, and a series of sparse Bayesian learning methods were proposed in [39][40][41][42].It is well known that these Bayesian algorithms face the challenge of high computational complexity.Recently, deep learning was applied to DOA estimation under impulsive noise in [43], revealing that the DOA prediction ability was reasonably enhanced by utilizing the powerful learning ability of neural networks.However, deep learning algorithms heavily depend on excessive training datasets, and models trained for specific conditions may perform poorly when applied to data models in other new environments.It cannot be ignored that the above-mentioned DOA estimation methods can be applied to the AVSA to deal with the performance degradation issues.However, as indicated above, they either face the problem of time consumption and off-grid.To solve the problems of off-grid and impulsive noise in sparse signal processing, by using the PFLOM to suppress impulsive noise, a three-step alternating iterative orthogonal matching pursuit algorithm is proposed in [44] to achieve the off-grid DOA estimation for an AVSA.
Inspired by de-outlier techniques and lower-order processing methods, an infinite weighted p-norm sparse iterative algorithm (IWPN-SIA) is designed to predict the DOA via the AVSA while considering impulsive noise in this paper.Noting that the participation of high outliers in the sparse iteration will produce large errors in the results, the infinite norm is used to weigh each snapshot of the sensors' output to shrink the high outliers.Although the output data of AVSA no longer has ultra-high outliers after the weighting process, the outliers are still significant values in the output dataset since all data in the output dataset is scaled at the same level.To further reduce the influence of impulsive noise, the Frobenius norm is extended to a lower order and then a p-order (1 < p < 2) cost function is constructed to obtain the expression of sparse signal power.The simulation results verify that the suggested algorithm can attain more precise and robust DOA approximation under strong impulsive noise.
In the following, (•) −1 , (•) H , and (•) * denote the matrix inverse, conjugate transpose and conjugate operation, respectively.⊗ indicates the Kronecker product.D{•} is a transverse matrix that yields the row vector of the key diagonal variables of the matrix in parentheses.tr{•} is the trace variable, [•] ηζ denotes the η-th row, and ζ-th column elements of the matrix.abs{•} means taking the absolute value of each element in the matrix.

Data Model of AVSA
Assuming a uniform linear array consisting of M AVSs with the unit spacing between sensors d as depicted in Figure 1, K far-field equal power narrowband signals from different directions θ k (k = 1, 2, • • • , K) are collected at the AVSA.The k-th signal object is incident on the m-th AVS, and its array manifold can be written as [11]  Then, the array diverse matrix of AVSA can be represented as is the acoustic pressure response parameter of the AVSA, and λ is representing signal period.Then at time t, the output vector by the AVSA can be written as and T is the impulsive noise data gained by the AVSA, and L is representing the snapshot.
For ease of calculation, the concept of the sparse grid is introduced into the model.The angle space −180 • ∼ 180 • is divided into N identical discrete nets, i.e., θ = { θ1 , θ2 , • • • , θN } and N K, assuming that the DOAs of signal sources are placed at the discrete grids, i.e., θ k ∈ θ.At time t the sparse signal vector of the AVSA can be calculated as and likewise is the array manifold matrix of 3M × N, and ] T is the waveform vector of the signal with dimensions N × 1.

α-Stable Distribution Model
The α-stable distribution is utilized as the noise model since it is famous for its heavy tail, stability and non-Gaussian characteristics, which enable it to perfectly simulate extreme anomalies in the signal.Unfortunately, there is no PDF closed-form expression for α-stable spreading.Thus, its fitness criterion is often utilized to express noise as [45] with and and α (0 < α ≤ 2) is the characteristic exponent, that manages the thickness of the tail of the PDF, γ (γ > 0) is the diffusion variable, which is identical to the variance of Gaussian distribution, µ denotes the position parameter, and β (−1 ≤ β ≤ 1) is the proportionality parameter.Noting that the α-stable spreading is proportional about µ, thus, it is named symmetric α-stable (SαS) spreading when β = 0.The SαS spreading represents a different PDF after different parameters are selected.Therefore, it is equivalent to the Gaussian distribution when α = 2 and β = 0, and it denotes the Cauchy distribution when α = 1 and β = 0.

The Proposed IWPN-SIA Algorithm
In this section, the IWPN-SIA algorithm is projected to realize the efficient estimation of DOA via the AVSA in the presence of impulsive noise.Firstly, by utilizing an infinite norm weighting technique, the impact of high outliers is mitigated through extensive numerical scaling.Subsequently, the Frobenius norm is extended to lower orders, and then a new cost function based on p-order is formulated to obtain the DOA estimation.Simultaneously, eigenvalue decomposition is incorporated to further attenuate noise.

Weighted Preconditioning with the Infinite Norm
To mitigate the impact of outliers and obtain a bounded covariance matrix, the weighted term is defined as where x m (t) ∞ denotes infinite norm of the signal vector x(t) ∈ 3M × 1.The weighted signal vector of the AVSA at time t is rewritten as so T is the weighted waveform vector, and T is the noise weighted vector.
It can be seen from ( 7) and ( 8) that the weighted processing technique is blind since it does not necessitate any previous information of impulsive noise data.Hence, it is also adaptable for use with various impulsive noise designs.In addition, the higher the outlier in the observed data, the greater the inhibition of these outliers in the weighted processing, demonstrating the dominance of this approach in processing high outliers.
Built on (8), the weighted correlation information of the AVSA is obtained as as with Y = [y(1), denoting the power of the n-th weighted signal, and η 2 v is the weighted noise strength of the respective observer.As mentioned in [36,46], the weighted array noises from separate observers exhibit no correlation, the weighted array signals from different sensors also show no correlation, and there is no correlation between the weighted array noises and signals.Consequently, the eigenvalue division of matrix R is subsequently performed by where • • • , a 3M are the corresponding eigenvectors, Λ s denotes a diagonal matrix consisting of K large eigenvalues, U s signifies the signal segment formed from eigenvectors corresponding to the large eigenvalues, Λ v represents the diagonal matrix of 3M − K small eigenvalues, U v expresses are called the noise subspace, which is also formed from the eigenvectors corresponding to the small eigenvalues, signal correlation matrix is denoted with R s , and R v stands for the noise correlation matrix.In practical applications, the estimated weighted array covariance matrix R usually needs to be swapped by the prototype correlation matrix

A Sparse Iterative Algorithm Based on p-Norm
Based on the concept of FLOS [33], the low-order covariance matrix is defined as where is the low-order weighted matrix.It is obvious that the diagonal element of the low-order covariance R c represents the signal power, and its value is closely related to p. Thus, R c can be rewritten as where |s wn (t)| p .Based on (8), the signal correlation matrix from the n-th grid is given as where the n-th column input of A( θ) is represented with A n .Next, the correlation matrix of interference combined with noise on the n-th grid can be represented as Based on [22,47], we can obtain the inverse of the matrix Θ n To further enhance the suppression ability of outlier, the Frobenius norm-based method in [22] is extended to p-norm.The cost function is constructed as follows Based on the statement in [48] and the operational relationship between the matrix norm and the trace in [47], for any matrix B ∈ M × M, it holds that According to ( 17), ( 16) can be further transformed into where Ψ is a low-order weighted matrix, which is defined as with , and Γ ζη denoting the entity in η-th column and ζ-th row of Γ.To simplify the calculation, its diagonal elements are extracted to create a new weighted matrix Ξ, which is defined as Then, (18) can be further simplified as The initial and secondary derivatives of g( qn ) in terms of qn are and Based on (23), we can have g ( qn ) > 0, showing that g ( qn ) has a unique solution.Setting g ( qn ) = 0, we obtain By substituting ( 15) into (24), the following can be obtained, where It is obvious that estimating qn by (25) requires the prior information of Ξ, and then solving Ξ also requires the prior information of qn .Therefore, an iterative algorithm must be used to estimate qn .Assuming q(j+1) n is the estimation of qn in the (j + 1)-th time.Next, q(j+1) where and Ultimately, the presented method is consolidated as Algorithm 1.
Algorithm 1: Infinite Weighted p-norm Sparse Iterative Algorithm 1: Calculate the weighting term (t) and the weighted signal vector y(t) with ( 7) and ( 8), respectively.2: Obtain the low-order sample covariance matrix R c with (11).3: Perform the eigenvalue decomposition to obtain the signal covariance matrix R s with (10).j) with ( 29).

Simulation Results and Discussions
Here, we compare the designed IWPN-SIA paradigm with classical algorithms MU-SIC [12], IAA [22], and advanced algorithms for impulsive noise including FLOM-MUSIC [34] and WARD-MUSIC [36].Unless otherwise stated, a ULA with M = 5 AVSs are used in the following simulations, the scanning range −180 • ∼ 180 • is divided into 361 grids with a 1 • interval, the error tolerance υ is set to 10 −3 , and the maximum iteration times T max is 25.The impulsive noise is modeled as SαS distribution, i.e., β = 0, and the diffusion variable γ is fixed to 1. Due to the features of α-stable distribution, the commonly used SNR becomes meaningless.To quantify the signal power and noise diffusion ratio, the generalized signal-to-noise ratio (GSNR) is computed with [27] as The two incoming waves of DOAs are θ 1 = 60 • , θ 2 = 134 • , respectively.For each simulation, 2000 Monte Carlo trials (N m = 2000) are executed.The resolution likelihood and root mean square error (RMSE) are examined for these algorithms.Among them, it is believed that the DOA prediction of two objects can be successfully distinguished whether the DOA estimation errors of these objects are lower than τ = 2 • .Suppose that the number of successful resolutions is N c in N m Monte Carlo trials, then the resolution probability is shown as N c /N m .The RMSE is defined as in the i-th Monte Carlo trial, and θk,i is the k-th signal estimated DOA.

The Influence of p on Resolution Probability and RMSE
Here, we provide the relationship between p and resolution probability and RMSE to determine the value of p.
The simulation results are depicted in Figures 2 and 3, respectively, with GSNR = −2 dB and L = 300 where we consider three cases: α = 1, α = 1.3, and α = 2.It can be seen from Figure 2 that the resolution probability of the IWPN-SIA method can reach 100% when α = 2 and p > 1.4.Although the resolution probability of the IWPN-SIA method cannot reach 100% within the range of 1 to 2 for p when α = 1.3, it has the highest resolution probability when p is approximately 1.7.It is obvious that compared to α = 1.3 and α = 2, the resolution probability of the IWPN-SIA method is the worst when α = 1, but it can be seen that its resolution probability is relatively high when p is in the range of 1.4 to 1.7.In addition, it can be seen in Figure 3 that the RMSE of the IWPN-SIA method tends to stabilize when α = 2 and p > 1.3.Obviously, its RMSE is the smallest when α = 1.3 and p is around 1.8.Moreover, it can be seen that the RMSE is relatively large throughout the entire range of p when α = 1, but its RMSE is the smallest when p = 1.7.Therefore, in the following simulation experiments, p is fixed to 1.7.

Comparison of the Resolution Probability
In this subsection, we compared the resolution probability of these algorithms from three aspects of GSNR, snapshots, and characteristic exponent α.
Resolution probability over GSNR with α = 1.3 and L = 300 can be seen in Figure 4.It is noticeable that the resolution probability of MUSIC and IAA methods cannot reach 100% with the gradual rise of GSNR, mainly due to the fact that both MUSIC and IAA methods are the DOA prediction techniques built on the Gaussian noise model.On the contrary, FLOM-MUSIC can effectively enhance the resolution probability by obtaining the covariance matrix through low-order processing.Furthermore, the performance of WARD-MUSIC is significantly improved by successively weighting snapshots to remove outliers and can achieve a resolution probability of 100% when the GSNR is not less than 6 dB.In general, the performance of IWPN-SIA, which combines outlier suppression and low-order processing, is significantly better than that of FLOM-MUSIC and WARD-MUSIC, and it is obvious that the proposed IWPN-SIA paradigm is able to attain a resolution likelihood of 100% when the GSNR is not less than 0 dB.This reveals that the proposed IWPN-SIA method has a higher resolution probability at low GSNR regions compared to other conventional methods.In Figure 4, it can be seen that the proposed IWPN-SIA method has a high resolution probability when the DOA values of the signal sources are placed onto the discrete values of the angle grid, considering the fact that the aggregate signal will no longer be sparse in the grid domain if the DOAs are not placed onto the discrete values of the grid.To illustrate the impact of the DOAs not being on the grid on the performance, the result is given in Figure 5, where the DOAs of the incident waves are set to 60.5 • and 134.1 • , and the other parameter settings remaining consistent.It can be seen from Figure 5 that compared to the DOAs of the signal on the grid, the resolution probability of the IWPN-SIA and WARD-MUSIC methods can still reach 100% when GSNR ≥ 0 dB and GSNR ≥ 6 dB, respectively, but their resolution probability is reduced when GSNR < 0 dB and GSNR < 6 dB, respectively.The resolution probability of other compared methods is decreased throughout the entire GSNR range when the DOAs of the signal are not placed onto the discrete values of the grid.To validate the potency of the designed IWPN-SIA method with a reduced quantity of snapshots, we compare the resolution probability of the above mentioned techniques across varying numbers of snapshots, ranging from 10 to 500, with the GSNR static at 6 dB, and the remaining coefficients being kept constant with Figure 4.The resolution probability results of these algorithms are presented in Figure 6, where it can be noticed that the resolution probability of the MUSIC and IAA algorithms shows a decreasing trend as the amount of snapshots extends.The reason is that the quantity of outliers increases with the number of snapshots received in the impulsive noise atmosphere, significantly impacting the resolution capability of these algorithms based on Gaussian assumptions.Although the resolution probability of the FLOM-MUSIC, WARD-MUSIC, and IWPN-SIA methods has been increased with the rise in snapshot quantity, the suggested IWPN-SIA method has the highest resolution probability in the whole snapshot region, and its resolution performance is particularly significant compared to other methods in small snapshots.The characteristic exponent α is an essential parameter of α-stable dispersion, which has the ability to measure the strength of impulsive noise.Figure 7 depicts the resolution probability of these methods versus the characteristic exponent α where α ranges from 0.4 to 2, L = 300, and the other parameters remain the same as Figure 6.It is noticed that the resolution accuracy of these techniques improves as α rises.The MUSIC, IAA, and FLOM-MUSIC methods fail when α < 1.In contrast to the MUSIC, IAA, and FLOM-MUSIC algorithms, while the resolution probability of the WARD-MUSIC and IWPN-SIA methods is significantly improved when 0.4 < α < 1, the proposed IWPN-SIA method still has a superior resolution probability compared to WARD-MUSIC.

Comparison of RMSE
Here, we simulate the RMSE of the designed method and the comparison methods mentioned above from three aspects, including GSNR, snapshots, and the characteristic exponent α.
The RMSEs of DOA estimates with α = 1.3 and L = 300 are shown in Figure 8.It is obvious that the RMSE efficiency of IAA and MUSIC algorithms is poor, especially at the points when GSNR is low.This reveals that the impulsive noise has a serious impact on their DOA estimation performance.Although the RMSE of FLOM-MUSIC, WARD-MUSIC, and IWPN-SIA has improved in the case of low GSNR, the approximation precision of the designed IWPN-SIA algorithm is the highest in the given GSNR region.We similarly consider the impact of the off-grid case on the RMSE.The DOAs of the incident wave are set to 60.5 • and 134.1 • , and the other parameter settings are consistent with Figure 8; the result is shown in Figure 9.It can be seen that compared to Figure 8, the RMSE of all compared methods has increased, especially in areas with low GSNR, but the performance of the proposed IWPN-SIA method is the best.For the purpose of assessing the effectiveness of the designed IWPN-SIA technique when dealing with small snapshots, Figure 10 shows the effect of the value of snapshots, where the GSNR is fixed to 6 dB and other parameters remain consistent with Figure 8, from which it can be observed that as the volume of snapshots rises, the RMSE of the MUSIC and IAA techniques no longer declines.This also indicates that the impulsive noise has a meaningful impact on their estimation behavior.Although the RMSE of the FLOM-MUSIC, WARD-MUSIC, and IWPN-SIA methods increases with the number of snapshots, their RMSEs are never lower than that of IWPN-SIA.Even under small snapshot conditions, the proposed IWPN-SIA algorithm continues to exhibit superior estimation performance.The comparison of different parameters indicates that the effectiveness of the proposed IWPN-SIA techniques is far better than the existing methodologies.When examining the RMSE of all techniques regarding characteristic parameter α, α ranges from 0.4 to 2, L = 300, and the other parameters remain the same as Figure 10; this is shown in Figure 11.It is worth noting that the α-stable dispersion gradually tends to the Gaussian spread with the increase of α.As a result, the estimation accuracy of each compared method is significantly enhanced.Clearly, the MUSIC, IAA, and FLOM-MUSIC methods have high RMSE when α < 1.In contrast, the RMSE of the WARD-MUSIC and IWPN-SIA methods is significantly lower than MUSIC and IAA methods when 0.4 < α < 1.Although the DOA prediction accuracy of MUSIC, IAA, and WARD-MUSIC methods is higher than the designed IWPN-SIA algorithm when α = 2, the proposed IWPN-SIA method has the lowest RMSE when α < 2.