Robust Direction Finding via Acoustic Vector Sensor Array with Axial Deviation under Non-Uniform Noise

: To minimize the major decline in direction of arrival (DOA) estimation performance for an acoustic vector sensor array (AVSA) with the coexistence of axial deviation and non-uniform noise, a two-step iterative minimization (TSIM) method is proposed in this paper. Initially, the axial deviation measurement model of an AVSA is formulated by incorporating the disturbance parameter into the signal model, and then a novel AVSA manifold matrix is deﬁned to estimate the sparse signal power and noise power mutually. After that, to mitigate a joint optimization problem to achieve the sparse signal power, the noise power and the axial deviation matrix, two auxiliary cost functions, are presented based on the covariance matrix ﬁtting (CMF) criterion and the weighted least squares (WLS), respectively. Furthermore, their analytical expressions are also derived. In addition, to further enhance their prediction accuracy, the estimated axial deviation matrix is modiﬁed based on its speciﬁc structural properties. The simulation results demonstrate the superiority and robustness of the proposed technique over several conventional algorithms.


Introduction
Accurate direction of arrival (DOA) estimation is a fundamental problem in a passive acoustic system. Acoustic scalar sensor (ASS), which measures the acoustic pressure information generated by underwater sound waves, has been intensively applied in various fields. However, it is impossible to analyze DOA estimation for a single ASS, which may limit its application on small-scale platforms. Acoustic vector sensor (AVS), which is another promising type of sensor used in the passive acoustic system, comprises an omni-directional ASS and two or three acoustic particle velocity (APV) components [1][2][3]. Compared with ASS, which is only capable of measuring the acoustic pressure information, AVS can also measure additional APV information in two or three orthogonal directions [4,5]. Since the direction information of acoustic sources is contained in the APV components, the acoustic vector sensor array (AVSA) is capable of providing more direction information on acoustic sources relative to their equivalent-aperture acoustic scalar sensor array (ASSA) [6,7]. By using the additional APV component information, the AVSA is able to estimate both elevation and azimuth without left-right ambiguity. Consequently, it is the core of many emerging applications, such as underwater acoustic communication [8], underwater target tracking [9], speech signal processing [10], sonar [11], and geoacoustic inversion problems [12].
The DOA estimation technologies based on the ASSA, such as CBF [13], MVDR [14], MUSIC [15], ESPRIT [16], and sparse representation [17] methods, have been extended to the AVSA. Among them, the typical CBF method is robust against mismatches. However, it usually encounters the problem of large side-lobe fluctuations and main-lobe spectral bandwidth. The MVDR method has low side lobes and a narrow beam width, but its performance may be degraded dramatically under considerable ambient noise. The MUSIC and ESPRIT methods are subspace-based algorithms that achieve super-resolution ability. However, their performance may deteriorate seriously under small snapshots, low signal-tonoise ratio (SNR), or correlated signals since they depend on the eigenvalue decomposition of the covariance matrix. To tackle this issue, the spatial smoothing technology [18] is proposed to decorrelate the correlated sources. However, the aperture of the array is usually lost, which may limit the performance of DOA estimation. Moreover, the particle-velocityfield (PVF) smoothing [19] and PVF difference smoothing [20] methods are also formulated to overcome this problem. Unlike the spatial smoothing technology, the aperture of the array is not lost for the PVF smoothing and PVF difference smoothing methods. However, the data dimension output of these methods is reduced due to the differential operation between additional APV information, which may cause the loss of the degree of freedom of the array.
In addition, with the vigorous development of compressed sensing technology, several sparse signal reconstruction (SSR) methods, such as 1 -norm based singular value decomposition ( 1 -SVD) [21], sparse Bayesian learning (SBL) [22], iterative adaptive approach (IAA) [23], and their variants [17,24,25], have been developed to deal with the DOA estimation problem. By using compressed sensing technology, the DOA estimation performance of these sparse methods is usually better than subspace-based methods under lower SNR, limited snapshots, and correlation sources. Among them, the 1 -SVD is a favorable SSR algorithm because of its guaranteed recovery accuracy. However, it usually requires prior knowledge of the number of sources and the resolution performance if it is limited under the small angle spacing of the sources. To mitigate this problem, two augmented cross-covariance matrices are reconstructed to estimate the DOA of acoustic sources [26,27]. However, in a complex underwater environment, the selection of the appropriate regularization parameter is still a challenging problem, which might decrease the DOA estimation accuracy of these algorithms. The SBL method is another famous following SSR technology. However, the disadvantage of SBL is that it guarantees less signal recovery accuracy than the 1 optimization because it is a probabilistic method. Moreover, these methods encounter the issue of high computational complexity. The IAA method, which is proposed in [23], is an amplitude estimation method without user-defined parameters. Since it can accurately estimate the DOA and amplitude information of acoustic sources, it has been widely applied in various fields [28,29].
Note that the superior performance of these above methods is based on the fact that the noise power of all components in the AVSA is the same, called the uniform noise. However, in a practical passive acoustic system, the main noise received by underwater AVSA comes from ambient noise rather than the internal noise of sensors. Consequently, the noise powers of the APV components are usually unequal to the acoustic pressure components in the condition of ambient noise. It is because APV components filter out some noise due to the directivity of AVS [30]. In this case, their DOA estimation performance seriously deteriorates because the difference in noise power received between the acoustic pressure and APV components is not taken into account. For an ASSA, since the channel response and the mutual coupling are not perfect, the noise received by different sensors may have different powers, called non-uniform noise. To solve this issue, many methods have been proposed in the last decade. In [31,32], deterministic and stochastic non-uniform maximum-likelihood (ML) estimators are proposed to estimate the DOAs and noise nuisance parameters by iteration. However, their practical application may be significantly restricted due to the nonlinear optimization problem involved in each iteration. To reduce the computational complexity, a novel DOA estimator, based on a power domain method, is presented in [33]. Although this method is free of estimating the nuisance parameters, the nonlinear optimization problem is still involved in the power domain. In [34], a linear transformation technology is exploited to deal with the inequality of noise powers, while this method suffers from the problem of being time-consuming since it is efficiently solved by second-order cone programming. In [35], two iterative methods called the iterative ML subspace estimation (IMLSE) and iterative least-squares subspace estimation (ILSSE) methods, are developed to estimate the DOA by combining traditional estimators. Moreover, a linear prediction DOA estimation approach is proposed in [36]. However, the DOA estimation performance of it is inferior to the iterative subspace methods when the SNR is high. Considering the advantages of the SSR technique, a new DOA estimation method is presented in [37] by using the tail minimization technique. These methods can be applied to the AVSA to improve the DOA estimation performance under the inequality of the acoustic pressure and APV components. However, they usually face the problem of high computational complexity. In [30], an augmented subspace MUSIC (ASMUSIC) is proposed to improve the DOA estimation performance of AVSA in the case of non-uniform noise by extending the number of virtual sources to the signal subspace. However, the different noise power of acoustic pressure and APV components in a single AVS is only considered. The noise power of each channel of AVS array may be different when any two sensors encounter imperfect responses in the AVSA array, in which case, the DOA estimation performance of the ASMUSIC method is limited. In addition, considering the characteristics of non-Gaussian noise in shallow water, a robust variational Bayesian inference method is proposed in [38] to address the problem of DOA estimation performance degradation under impulse noise. Then, by using the common sparsity of signal and the element-wise sparsity of the impulsive noise for different snapshots, the posterior estimation of signal and impulsive noise components is obtained based on the variational Bayesian inference.
The above-mentioned DOA estimation techniques can be extended to the ideal AVSA, in which the axes of any two sensors point in the same direction. In practical applications, the axes of any two sensors can not be guaranteed to point in the same direction because of errors in installation and the directivity of AVS. These problems may deteriorate the performance of existing DOA estimation techniques dramatically. To handle this problem, by utilizing the specific structural properties of the axial deviation matrix, a modification method is presented in [39] to estimate the DOA accurately by combining existing technologies. In summary, almost all the above-mentioned DOA estimation techniques consider either the axial deviation or the non-uniform noise for an AVSA. However, the DOA estimation technology in the presence of axial deviation and non-uniform noise is rare. Therefore, we are committed to overcoming this situation under the condition of axial deviation and non-uniform noise via an AVSA.
In this paper, a robust two-step iterative minimization (TSIM) method is proposed to estimate the DOA via an AVSA with the coexistence of axial deviation and non-uniform noise. The main contributions of this paper are summarized as follows:

•
We present the axial deviation measurement model of an AVSA by incorporating the disturbance parameter and formulating a novel AVSA manifold matrix to estimate the noise vector as a potential signal. • Based on the covariance matrix fitting (CMF) criterion and weighted least squares (WLS), the cost functions with respect to sparse signal power, noise power, and axial deviation parameters are constructed, and the properties of concave functions are used to transform nonlinear optimization problems into linear optimization problems. Then, their analytical expressions are derived.

•
We formulate the modified methods of the covariance matrix and axial deviation matrix in non-uniform noise and axial deviation, respectively, to improve the estimation accuracy of sparse signal power in each iteration.

•
Extensive numerical examples illustrate the superiority and robustness of our proposed method over the existing ones in the condition of axial deviation and nonuniform noise via an AVSA.
The remaining of this paper is organized as follows. The classical AVSA model and the AVSA model with axial deviation are briefly described in Section 2. A TSIM approach is formulated to estimate sparse signal power, noise power, and axial deviation matrix in Section 3. Numerical simulation examples are given in Section 4. Finally, the conclusion of the whole paper is summarized in Section 5.
Notations: A T , A H , A −1 , and A † denote transpose, conjugate transpose, inverse, and pseudo-inverse operations, respectively. [A] τη is an element of τ-th row and η-th column in A. (A) takes the real part of every element in A. Tr(·) denotes the trace operation of the matrix. max(A) means to take the maximum value of A. ||A|| F stands for the Frobenius norm of A. ⊗ indicates the Kronecker product operation. diag(·) and blkdiag(·) are the diagonalization and block diagonalization operations, respectively. D(·) represents the differential operation.

The Classical Measurement Model of an AVSA
As given in Figure 1a, it is assumed that a uniform linear array (ULA) is composed of M-element AVS, which receives K far-field narrowband planer wave acoustic sources with unknown DOAs where each AVS consists of three components, including an omni-directional ASS as well as two orthogonal APV components named u x and u y , respectively, and d indicates the spacing between adjacent AVS. It is assumed that the number of sources has been correctly chosen by exploiting the proposed detection method in [40]. According to [4], the 3 × 1 array manifold vector of the m-th AVS for the k-th source under the classical AVSA model can be expressed as Then, the 3M × K manifold matrix of the AVSA can be formulated as where λ is the signal wavelength. The output of the AVSA at time l can be written as where being the outputs of the m-th AVS pressure channel and the two APV channels, s(l) = [s 1 (l), s 2 (l), · · · , s K (l)] T denotes the sources vector, and e(l) = [e 1 (l), e 2 (l), · · · , e 3M (l)] T is the Gaussian noise vector. With L snapshots collected, the array output can be given as

The Axial Deviation Measurement Model of an AVSA
It can be seen from Figure 1a that the u x (or u y ) axes of adjacent AVS or non-adjacent AVS are arranged in the same direction. However, in practical applications, it is difficult to arrange the u x axes of adjacent AVS or non-adjacent AVS in the same direction due to errors in installation and the directivity of AVS. This case, which is called the axial deviation model of an AVSA in this paper, as given in Figure 1b with µ m denoting the axial angle deviation parameter of the m-th AVS relative to the axis of the first AVS at the coordinate origin, degrades the DOA estimation accuracy abruptly because of the involvement of the parameter µ m . In the following, we consider the axial deviation model of an AVSA.
Let µ = [µ 1 , µ 2 , · · · , µ M ] T be the axial deviation vector of the AVSA. In this case, the array manifold vector of the m-th AVS for the k-th source can be rewritten as with which is the axial deviation matrix of the m-th AVS. Then, the steering matrix of the AVSA can be rewritten as A(θ, µ) = [a(θ 1 , µ), a(θ 2 , µ), · · · , a(θ K , µ)], where a(θ k , µ) = CA p (θ k )H(θ k ) with C denotes the axial deviation matrix of the AVSA, Therefore, (7) can be rewritten as Considering the sparse characteristics of signals in the spatial domain, the spatial domain of a range (−180 • , 180 • ] can be uniformly divided into N potential DOAs θ = {θ 1 ,θ 2 , . . . ,θ N } with θ k ∈θ, where N is much greater than the number of acoustic sources and the number of AVS. Then, the manifold dictionary of the classical AVSA can be reformulated as A(θ) = [a(θ 1 ), a(θ 2 ), · · · , a(θ N )]. Thus, (11) can be reformulated as withS = [s 1 ,s 2 , · · · ,s N ] T being the extension of S from θ toθ, where the non-zero entries ofS indicates the true sources locations. Assuming that the signal and the noise are uncorrelated, then the expected covariance matrix of the sparse signal can be written as withP and Q being the sparse signal covariance matrix and the noise covariance matrix, respectively,P and Ls ns H n , U m represents the m-th column data of I 3M , σ 2 mp , σ 2 mx , and σ 2 my are the noise powers of the acoustic pressure and APV components of the m-th AVS, respectively, and q m = [σ 2 mp , σ 2 mx , σ 2 my ] is the noise power vector of the m-th AVS. In the actual environment, the noise collected by AVS is usually ambient noise rather than internal sensor noise. Furthermore, due to the directivity of APV components, as well as the non-ideal channel response and mutual coupling from pressure and APV channels, there often exists (15) is reduced to the uniform noise model with

The Proposed Method
Noting thatθ is the row-sparse vector due to the number of grids being much greater than the number of acoustic sources. Thus, the DOAs of acoustic sources can be addressed by estimating the row parameters ofS. However, the involvement of the unknown axial deviation matrix C in (12) may cause the estimated DOA to deviate from the true DOA. In addition, the non-uniform noise in (12) may deteriorate the DOA estimation performance of AVSA. To obtain the DOAs of acoustic sources in the presence of the axial deviation and the non-uniform noise, the sparse signal power, the axial deviation matrix, and the noise power need to be estimated simultaneously. In this section, a TSIM algorithm is proposed to jointly estimate them based on the CMF criterion and WLS.
Let us define a new manifold matrix of AVSA to jointly estimate the sparse signal power and the noise power, where A n represents the n-th column data of A(θ). Then, (12) can be further rewritten as It can be seen from (18) that B is regarded as a new array manifold matrix composed of a real array manifold matrix [CA 1 , · · · , CA N ] and the identity matrix I 3M , in which the former is used to estimate the sparse signal and the latter is used to estimate the noise vector. In addition, it is noted that the noise vector E in (12) is estimated as a potential sparse signal vector in (18), which is formulated in this paper. Thus, it can be understood that the signal is divided into N + 3M discrete grids in the spatial angle space, where the first N sparse signal vectors correspond to the angle space of the actual sparse signal, and the following 3M sparse signal vectors are the noise vector. IfS can be recovered from (18), the sparse signal power and the noise power can be estimated based onp r = 1 Ls rs H r with r = 1, 2, · · · , N, N + 1, · · · , N + 3M, where {p r } N r=1 = {p n } N n=1 and {p r } N+3M r=N+1 are the sparse signal power and the noise power, respectively. Then, based on (18), (13) can be further transformed into with ∆ being (21).
From (20) and (21), it is clear that the data from the (N + 1)-th column to the (N + 3M)-th column of B in (17) is related to the noise vector if {p r } N+3M r=N+1 is treated as the potential signal powers when estimating the true signal powers {p r } N r=1 in ∆. The signal covariance matrix corresponding to the r-th column of B can be formulated as R r =p r B r B H r , r = 1, 2, · · · , N, N + 1, · · · , N + 3M.
Then, the interference plus noise covariance matrix corresponding to the r-th column of B can be given by Obviously, if the sparse signal power, the noise power, and the axial deviation matrix can be accurately estimated from the interference plus noise covariance matrix Γ r , the DOAs can be determined by searching the spectral peak of the signal powers {p r } N n=1 . As stated in [41], since the second-order statistics model has a higher array output SNR, the accuracy of DOA estimation can be improved by utilizing the covariance matrix of the array output. It can be seen from (23) that it is difficult to directly estimate the sparse signal power and the noise power due to the covariance matrix R being unknown.
Therefore, the covariance matrix R in (23) is replaced bŷ which is the sample covariance matrix of AVSA, then the interference plus noise covariance matrix corresponding to the r-th column of B can be further rewritten asR − R r . Clearly, in the noiseless environment, the product ofR − R r and R −1 r should be the identity matrix. Therefore, in the first step, the regularized weighted CMF method is developed to jointly estimate the sparse signal power and the noise power while keeping the axial deviation matrix fixed,p (j+1) r = arg min Fp(p r ) with where the first part of (26) indicates the CMF term, the second part represents a penalty term further enforcing the sparsity ofs r , λ s is the regularization parameter balancing the relationship between the sparse signal and fitting error, and 0 < q ≤ 1 denotes a user parameter constraining the sparsity of signals. The second step optimizes over the axial deviation matrix based on the regularized WLS while fixing the sparse signal power and the noise power, with where the first part Γ is a WLS term, the solution of which is consistent with IAA [23] when C is the identity matrix, the second part λ c ||C|| q F is the penalty term with λ c > 0 being the regularization parameter controlling the tradeoff between the WLS term and the Frobenius norm term, and q constraining the efficiency of the axial deviation matrix.
It is proven that both (26) and (28)  and C from (26) and (28). In the next subsection, how to estimate {p r } N+3M r=1 and C will be presented.

Estimating the Sparse Signal Power and the Noise Power
As stated in [42], for any concave function F y (y) of y > 0, the following inequality is constant y (y) ≤ y (y 0 ) + y (y 0 )(y − y 0 ), (29) where y 0 > 0 is a given value. Suppose y (y) = y q 4 , then it is observed that y (y) is a concave function. Let y = s r Noting that the first-order derivative of y (y) = y q 4 with respect to y can be expressed as Inserting As a consequence, substituting y = s r Multiplying λ s and then adding Γ on both sides of (33), we get Fp(p r ) ≤ Dp(p r ) (34) with It is interesting to note thats (j) r , which has been obtained in the j-th iteration, is known in (35). Thus, (25) can be further equivalent tō with Then, the estimate ofp r can be updated at the (j + 1)-th iteration as (see Appendix A) , r = 1, 2, · · · , N, N + 1, · · · , N + 3M Furthermore, according to (17) and (21), the sparse signal power and the noise power can be obtained as, respectively, with and In view of (21) and (40), the noise power of the acoustic pressure and APV channels of the m-th AVS can be written asσ Obviously, according to the estimated signal power in (39) and noise power in (40), the covariance matrix in (A12) is reconstructed to approximate the expected covariance matrix. The DOA estimation performance of AVSA may be further improved; however, it cannot be ignored that the performance of AVSA may be severely deteriorated due to the involvement of the axial deviation matrix C in (41) and (42). To overcome this problem, how to obtain the desired axial deviation matrix is further analyzed in the following subsection.

Estimating and Modifying the Axial Deviation Matrix
As mentioned above, (28) is also a nonlinear function with respect to C due to the involvement of q, so it is difficult to solve C by optimizing (27). To deal with the optimization problem in (27), a coarse estimate for C is first implemented. Then, based on the specific structural properties of the axial deviation matrix presented previously in (9) and (10), the coarsely estimated axial deviation matrix is further modified.

Coarsely Estimating the Axial Deviation Matrix
Similar to the derivation of the sparse signal power in Section 3.1, it is assumed that z (z) = z q 2 . Noting that z (z) is a concave function of z > 0. Then, we can get where z (z 0 ) represents the function value of the derivative of function z (z) when z = z 0 . Assuming that z = C 2 F and z 0 = C (j) 2 F . Then, we hold Inserting Multiplying λ c and then adding Γ on both sides of (53), we can obtain with In (55), since C (j) has been obtained in the j-th iteration, the minimization of (27) with respect to C can be further rewritten as Significantly, the axial deviation matrix coarsely estimated by minimizing (56) is related to the corresponding potential sparse signals {s r } N r=1 . Moreover, it can be seen that the inverse operation of matrix Γ r needs to be performed N times in each iterative, resulting in a high computational burden. To overcome this problem, Γ r can be replaced by R based on the derivation of (A4) to (38). In addition, to improve the estimation accuracy of the axial deviation matrix, we are only concerned with the sparse signals corresponding to K estimated DOAs. Then, (57) can be further formulated as where A indicates the DOAs corresponding to the maximum K Assuming thatC is the coarse estimation of C. In Appendix B, it is proven that by taking the partial derivative of G c (C) with respect to C, we haveC where and

Modifying the Axial Deviation Matrix
It is interesting to note that the estimatedC (j+1) does not satisfy the expected axis deviation matrix, mainly because of the following reasons: (a) the magnitude of elements in C (j+1) may be greater than one; (2) the element values inC (j+1) may be a complex number due to the disturbance of the calculation and the involvement of the complex signal. To deal with this problem, the following operation onC (j+1) is performed by Let η = 3m − 1. Based on (9) and (10), it is obvious that both Φ where and g (j+1) Since the coarsely estimatedC (j+1) is carried out in a noisy environment, the obtained G (j+1) m may contain the disturbance information. As a result, the modified axial deviation matrix is not accurate enough. To solve this problem, the estimated axial deviation matrix is modified again. It is assumed that the v x and v y directions at the coordinate origin, as described in Figure 1b, are used as the reference directions of other APV components. To eliminate the axial angle disturbance information, the relative axial deviation matrix from the reference AVS can be performed by Obviously, the axial deviation parameter of the m-th AVS is determined bŷ with ϕ denoting the DOAs of K signals determined by , and υ indicating a user tolerance parameter. Then, according to (9) and (10), the desired axial deviation matrix of the AVSA can be reconstructed bŷ As a consequence, when the estimated sparse signal power has a negligible change compared with the previous estimation, the final DOA can be realized by searching the spectral peak ofp (j+1) .

Computational Complexity of the Proposed Algorithm
The computational complexity of our proposed TSIM method with MUSIC, ASMUSIC, IMLSE, and AIAA methods are compared in this subsection. MUSIC has a complexity of O{M 2 L + M 2 N}. The computation complexity of ASMUSIC is O{M 2 L + M 2 N + 3M 3 }. By considering T iter (T iter < T max ) iterations, for IMLSE, we have O{M 2 L + M 2 N + T iter M 3 }, AIAA is O{(13T iter + 2)M 2 N + (8T iter + 7)MN + (3T iter + 2L)MN + T iter KN + 42M 3 }. For our proposed TSIM method, it is about O{31T iter M 2 N + (12T iter + 2)MN + (5T iter + 4L)N + T iter KN + 42M 3 }. It is noted that the computational complexity of our proposed TSIM method is the highest. However, the DOA estimation performance of it is best in the compared methods, which is depicted in the following simulation results.

Numerical Results
In this section, the DOA estimation performance of our proposed method is evaluated by simulations, and MUSIC [15], IAA [23], ASMUSIC [30], IMLSE [35], and alternating iterative adaptive approach (AIAA) [39] are taken into comparison. Further, the Cramer Rao Lower Bound (CRLB) [4] is plotted to compare the performance of these algorithms. In our simulation experiments, a four-element AVS ULA is utilized, and the adjacent AVS spacing is fixed to one wavelength [43]. To simplify, the axis orientation of APV at the coordinate origin position is set to µ 1 = 0 • and the axial deviation parameter of others {µ m } 4 m=2 is generated with a uniform distribution in each simulation experiment, where the mean value is µ, and the variance is σ 2 . The grid interval is fixed to 2 • . For our proposed method, T max is fixed to 20, τ is fixed to 10 −3 , and both λ s and λ c are fixed to 0.625 [44]. According to the analysis in [39,45], p is set to 0.5 and υ is fixed to 45. In addition, it is assumed that the non-uniform noise covariance matrix in (13) is Q = diag{3, 1.5, 2, 5, 2, 1.5, 4, 2, 5, 2, 30, 2, 1.5}.
As stated in [35,36], the worst noise power ratio (WNPR) and the SNR are defined as, respectively, and where σ 2 max and σ 2 min represent the maximum and minimum values in Q, respectively, and σ 2 s denotes the signal power. The average root mean square error (RMSE) is introduced as whereθ wj indicates the estimated result of θ k by the w Monte Carlo simulation experiment, and W = 3000 is the total number of Monte Carlo simulation experiments. First, it is assumed that two far-field equal-power uncorrelated sources impinge on the AVSA from θ 1 = −14 • and θ 2 = −14 • + ∆, where ∆, denoting the DOA separation of the two sources, varies from 25 to 55. The RMSE of these compared methods is given in Figure 2 with SNR = 5dB and L = 300, where we consider four cases: (a) uniform noise with µ = 0 and σ = 0; (b) uniform noise with µ = 15 and σ = 4; (c) non-uniform noise with µ = 0 and σ = 0; (d) non-uniform noise with µ = 15 and σ = 4. It can be seen from Figure 3 that there is a clear gap between the IMLSE and CRLB curves in the range of small DOA intervals in four cases. As is depicted in Figure 3a, when ∆ ≥ 30, the DOA estimation performance of all the compared methods except the IMLSE method is within an acceptable range if the background noise environment is uniform and there is no axial deviation in the AVSA. Compared with the simulation result in Figure 3a, it can be seen from Figure 3b that the performance of MUSIC and ASMUSIC methods deteriorates when ∆ < 35, revealing the high importance of axial deviation in the MUSIC and ASMUSIC methods. In addition, from Figure 3c, it is noted that the performance of the AIAA method further deteriorates when ∆ < 30 compared to the results shown in Figure 3b. The reason lies in the fact that the AIAA method only considers the influence of axial deviation on DOA estimation without considering the non-uniform noise. On the other hand, the results in Figure 3b,c indicate that the axial deviation has a more significant effect on the performance of the ASMUSIC method compared with the non-uniform noise. Further, it is witnessed from Figure 3d that the TSIM method can still offer a satisfactory performance with ∆ < 45, while the performance of other compared methods further deteriorates under axial deviation and non-uniform noise, showing the robustness and effectiveness of the TSIM method under the scenario of the axial deviation and the non-uniform noise. Then, we compare the RMSE of these methods with the SNR varying from −15 to 10 dB, where all simulation parameters maintain unchanged, as shown in Figure 3d, except that ∆ is fixed to 52. The simulation result is given in Figure 2. It can be seen that the performance of all the compared methods is improved as the SNR grows except MUSIC. The RMSE of MUSIC is limited as the SNR continues to increase. It is clear that ASMUSIC and AIAA show worse performance in the low SNR region. In fact, this is mainly because ASMUSIC and AIAA take into account either non-uniform noise or axial deviation. Although the RMSE of the IMLSE method and the TSIM method tends to be flat in the whole SNR region, the TSIM method can still provide a superior DOA estimation performance.
We also test the advantage of our proposed TSIM method for different µ with the previous settings except for the SNR being fixed to 5. The simulation result is shown in Figure 4. It is noted that the performance of MUSIC, ASMUSIC, and IMLSE methods deteriorate to some extent when the µ is large. The main reason is that the axial deviation parameter is ignored in the MUSIC, ASMUSIC, and IMLSE methods. On the contrary, the performance of the AIAA method and the TSIM method is further improved by modifying the observation data through the estimated axial deviation matrix. Obviously, the TSIM method is capable of maintaining superior performance compared to the AIAA method.  Finally, we verify the ability of the TSIM method for different WNPR, where the noise power of the pressure channel of the last AVS is varied from 10 to 50 while the noise power of other channels remains unchanged. According to (71), we can see that the WNPR varies from 20 to 100. Following the previous setting.s except for µ = 15, the RMSE of DOA estimation for different WNPRs is calculated. The result obtained by these compared methods is illustrated in Figure 5. It is noted that the performance of MUSIC, ASMUSIC, and AIAA severely deteriorates when the WNPR becomes large. Although IMLSE and TSIM methods can maintain DOA estimation performance in the whole range of WNPR, the DOA estimation accuracy of the proposed TSIM method is significantly higher than that of the IMLSE method. The reason is that the IMLSE method only considers the non-uniform noise. However, the TSIM method considers the coexistence of the axial deviation and the non-uniform noise via an AVSA.

Conclusions
In this paper, we proposed a robust TSIM method for DOA estimation via an AVSA with the coexistence of axial deviation and non-uniform noise. In our proposed TSIM method, based on the estimated results of the sparse signal power and the noise power in the previous iteration, the covariance matrix is reconstructed to approximate the expected covariance matrix. Then, by utilizing the structural characteristics of the axial bias matrix, the estimated axial deviation matrix is further modified to estimate the sparse signal power and the noise power more accurately in the next iteration. Extensive numerical simulations have verified that our proposed method can provide robust DOA estimation for an AVSA with the coexistence of axial deviation and non-uniform noise.

Appendix B
Based on the Frobenius norm properties in [47], we have Tr CC H (A13) Using the matrix derivative properties [47], the first-order differential of the last three terms of (A13) are respectively given as D Tr R −0.5(j+1) CA The partial derivatives of G c (C) can be expressed as C. (A17)