The research of direction of arrival (DOA) estimation plays a crucial role in the field of array signal processing [1-3], which is widely used in military and economic fields, such as radar [4,5], communication [6,7], and sonar [8]. In the past few decades, subspace decomposition methods have developed rapidly and achieved great success. This method’s core process is obtaining the signal subspace from the sampling covariance matrix of the incident signals. A representative method of this type is the multiple signal classification (MUSIC) [9] method. In recent years, signal sparsity theory has developed rapidly and has been well applied in the research of array signal processing, such as DOA estimation [10-12]. Compared with subspace decomposition techniques, sparse DOA estimation methods can provide better accuracy in more harsh environments, such as low signal-to-noise ratio (SNR) or small snapshot situations [13,14]. However, the above methods did not consider non-ideal factors, such as the mutual coupling between array elements. Weiss et al. [15] demonstrated that mutual coupling significantly impacts the accuracy of the incident signal’s DOA when the array manifold matrix is affected by it.

Recently, several papers have emerged on DOA estimation under array mutual coupling. The solution methods for this problem mainly consist of three categories. The first category is preprocessing methods. This category of methods construct a preprocessing operator based on the number of mutual coupling coefficients, which can effectively reduce the negative impact of array mutual coupling. After preprocessing, the non-ideal array manifold matrix affect by mutual coupling can be restored to the ideal Vandermonde matrix [16]. However, the construction of the preprocessing operator comes at the cost of sacrificing the array aperture [17,18]. The second category involves the unique deformation of the production of the mutual coupling matrix (MCM) and the array manifold matrix. Since the MCM is a Toeplitz matrix and the array manifold matrix is a Vandermonde matrix, their product can be represented as a new transformation. Additionally, based on the orthogonality between the noise subspaces and the steering vector of the incident signal [14], Liao et al. [19] and Ge et al. [20] proposed a rank deficiency method, respectively. Since the noise subspace is indispensable in the methods above, these methods are not applicable in scenarios with coherent signals. The third category of methods is the block sparse method. For example, according to the preprocessing method proposed by [17], Dai et al. [21] restored the array manifold matrix affected by mutual coupling to the ideal Vandermonde matrix and established a ${l_{2,0}}$ mixed norm to solve the model, where ${l_{2,0}}$ is the mixed norms of ${l_2}$ and ${l_0}$; Wang et al. [22] used the unique deformation of the product of the MCM and the array manifold matrix to establish a new block sparse model about mutual coupling. Unlike [21], Wang et al. [22] solved this problem using the ${l_{F,0}}$ norm, which provided a novel perspective for us, where ${l_{F,0}}$ is the mixed norms of ${l_F}$ and ${l_0}$. Since singular value decomposition (SVD) is used for preprocessing, the above methods still apply to coherent signals. However, the DOA estimation accuracy of the above two methods is not ideal; in order to further improve the DOA estimation performance of the method proposed by [22], Wang et al. [23] used a set of weights that can enhance the sparsity of the ${l_{F,1}}$ norm, which improved the DOA estimation accuracy, where ${l_{F,1}}$ is the mixed norms of ${l_F}$ and ${l_1}$. Meng et al. [24,25] weighted the signal subspace affected by mutual coupling appropriately, which also improves the accuracy of DOA estimation. However, the above block sparse methods have high computational complexity, so they cannot achieve DOA estimation quickly.

In this paper, a low computational complexity Newton-like method is designed to solve the ${l_{F,0}}$ norm, achieving DOA estimation under unknown mutual coupling. Firstly, a block sparse model under array mutual coupling is established via signal subspace. Then, we use exponential to represent the ${l_{F,0}}$ norm and establish the Lagrangian function ${L_\sigma }(\bar {\boldsymbol{B}})$. However, the Lagrangian function’s Hessian matrix cannot always guarantee positive definiteness, and the dimension of the Hessian matrix of ${L_\sigma }(\bar {\boldsymbol{B}})$ is very large. The dimensionality of the Hessian matrix is too large, which means that the calculation of its inverse matrix will be a difficult task. To maximize computational efficiency as much as possible, we design a Newton-like method based on mapping $\zeta (\bar {\boldsymbol{B}})$ to achieve DOA estimation and reduce its computational complexity by matrix inversion lemma. Finally, simulations demonstrate the effectiveness of the proposed method in DOA estimation.

Assuming the number of signal sources is Q, and the incident DOA of the signals are ${\boldsymbol{\theta }} = {[}{\theta _1},{\theta _2},\cdots ,{\theta _Q}{]}$. The array consists of M elements with a uniform spacing of half signal wavelength, i.e., $d = \bar \lambda /2$, $\bar \lambda $ is the signal’s wavelength. As shown in Fig. 1, the received signal of the array at time t under array mutual coupling is represented aswhere ${\boldsymbol{X}}(t)$ is the received signal from far-field at time t under mutual coupling, ${\boldsymbol{a}}({\theta _q}) = {[1,\beta ({\theta _q}),\cdots ,\beta {({\theta _q})^{M - 1}}]^{\mathrm{T}}}$ is the steering vector unaffected by mutual coupling of the qth far-field signal, $\beta ({\theta _q}) = {\text{exp}}\left( - {\mathrm{j}}{{2\text{π} d\sin {\theta _q}}}/{{\bar \lambda }}\right)$, and ${\boldsymbol{A}} = [{\boldsymbol{a}}({\theta _1}),{\boldsymbol{a}}({\theta _2}),\cdots ,{\boldsymbol{a}}({\theta _Q})]$ is the array manifold matrix unaffected by mutual coupling, ${\boldsymbol{s}}(t) = [{{\boldsymbol{s}}_1}(t), {{\boldsymbol{s}}_2}(t), \cdots , {{\boldsymbol{s}}_Q}(t)]^\text{T}$ is the source signal, ${\boldsymbol{n}}(t)$ represents Gaussian additive white noise $ \sim {{\mathrm{CN}}}(0,\sigma _n^2{{\boldsymbol{I}}_M})$, $\sigma _n^2$ is the noise power, ${\boldsymbol{C}} = {\text{Toeplitz}}({\boldsymbol{c}},{\boldsymbol{0}}_{M - P}^\text{T})$ is the MCM, ${\boldsymbol{C}} \in {{{\bf{C}}}^{M \times M}}$ is a Toeplitz matrix, where c is the row vector composed of mutual coupling coefficients, ${\boldsymbol{c}} = [1,{{{c}}_1}, \cdots ,{{{c}}_{P - 1}}]$ and satisfies $\begin{array}{*{20}{c}} {1 \gt \mid {{{c}}_1}\mid \gt \mid {{{c}}_2}\mid \gt \cdots \gt \mid {{{c}}_{P - 1}}\mid } \end{array}$, P is the length of ${\boldsymbol{c}}$.

Therefore, the received signal under mutual coupling can be represented aswhere $\tilde {\boldsymbol{A}}$ is the array manifold matrix affected by mutual coupling, $\tilde {\boldsymbol{A}} = [{\boldsymbol{Ca}}({\theta _1}),{\boldsymbol{Ca}}({\theta _2}), \cdots ,{\boldsymbol{Ca}}({\theta _Q})]$. The ideal covariance matrix can be represented aswhere ${\text{E}}\{ \cdot \} $ is the mathematical expectation.

When the number of snapshots is limited, the covariance matrix can be expressed aswhere T represents the number of snapshots, $ {\text{(}\cdot \text{)}}^{\text{H}} $ is the conjugate transpose, ${{\boldsymbol{\varSigma }}_s}$ is a diagonal matrix composed of the first Q large eigenvalues, ${{\boldsymbol{\varSigma }}_s} = {\mathrm{diag}}\{ {\lambda _1},{\lambda _2},\cdots ,{\lambda _Q}\} $, and ${{\boldsymbol{\varSigma }}_n}$ is a diagonal matrix composed of the last Q+1 to M small eigenvalues, ${{\boldsymbol{\varSigma }}_n} = {\mathrm{diag}}\{ {\lambda _{Q + 1}},{\lambda _{Q + 2}},\cdots ,{\lambda _M}\} $; ${{\boldsymbol{U}}_s}$ represents the signal subspace, ${{\boldsymbol{U}}_s}$ is composed of eigenvectors corresponding to the first Q large eigenvalues; ${{\boldsymbol{U}}_n}$ represents the noise subspace, ${{\boldsymbol{U}}_n}$ is composed of eigenvectors corresponding to the last Q+1 to M small eigenvalues.

When the signal is incoherent, there exists a nonsingular matrix ${\boldsymbol{B}} \in {{{\bf{C}}}^{Q \times Q}}$, and the signal subspace can be represented as

According to [26], the product of MCM and the incident signal’s steering vector can be represented aswhere${\boldsymbol{\varLambda}} ({\theta _q}) \in {{{\bf{C}}}^{M \times (2P - 1)}},\;{\boldsymbol{\varLambda}} ({\theta _q})$ is only related to the DOA of the incident signal and is independent of the mutual coupling coefficients, andwhere $\varOmega ({\theta _q})$ is a constant. In this paper, it is assumed that $\varOmega ({\theta _q}) \ne 0$, and the expression of ${\boldsymbol{\varGamma}} ({\theta _q})$ iswhere ${\boldsymbol{\varGamma }}({\theta _q}) \in {{{\bf{C}}}^{2P - 1}}$ andwhere $l = 1,2, \cdots ,P - 1$.

By substituting the above results into (5), we can obtain thatwhere $\hat {\boldsymbol{A}} = [{\boldsymbol{\varLambda}} ({\theta _1}),{\boldsymbol{\varLambda}} ({\theta _2}),\cdots ,{\boldsymbol{\varLambda}} ({\theta _Q})] \in {{{\bf{C}}}^{M \times Q(2P - 1)}}$ can be regarded as a block array manifold matrix, ${\boldsymbol{\varLambda}} ({\theta _q}) \in {{{\bf{C}}}^{M \times (2P - 1)}}$ can be regarded as a block steering matrix that is only related to the DOA and independent of the mutual coupling coefficients, and assuming $\tilde Q = 2P - 1$, ${\boldsymbol{\varLambda}} ({\theta _q}) \in {{{\bf{C}}}^{M \times \tilde Q}}$, andwhere ${\boldsymbol{\varDelta }} \in {{{\bf{C}}}^{Q\tilde Q \times Q}}$. Since ${\boldsymbol{\varDelta }}$ is a block diagonal matrix, thereforewhere ${{\hat {\boldsymbol{B}}}} = {{{\boldsymbol{\varDelta}} {\boldsymbol{B}}}} = {{\text{[}}{{\hat {\boldsymbol{B}}}}_1^\text{T},{{\hat {\boldsymbol{B}}}}_2^\text{T}, \cdots ,{{\hat {\boldsymbol{B}}}}_Q^\text{T}{\text{]}}^\text{T}} \in {{{\bf{C}}}^{Q\tilde Q \times Q}}$ is the block signal after deformation, ${\hat {\boldsymbol{B}}_q} = {{\varOmega }}({\theta _q}){\boldsymbol{\varGamma}} ({\theta _q}) \otimes {{\boldsymbol{B}}_{(q,:)}} \in {{{\bf{C}}}^{\tilde Q \times Q}}$, $q = 1,2, \cdots ,Q$, ${{{\hat {\boldsymbol{B}}}}_q}$ represents the qth block matrix of ${{\hat {\boldsymbol{B}}}}$, ${{\boldsymbol{B}}_{(q,:)}}$ represents the qth row’s elements of the matrix ${\boldsymbol{B}}$, $ \otimes $ is the Kronecker product. Since ${{\hat {\boldsymbol{A}}}}$ is only related to the angle, and the DOA of the incident signal in the spatial domain is sparse, a block sparse model about mutual coupling can be constructed according to (14).

For the purpose of DOA estimation, the spatial region is evenly partitioned into N grids, i.e., the spatial domain is ${{{\boldsymbol{\varTheta}} }} = [{\theta _1},{\theta _2}, \cdots ,{\theta _N}]$, the block sparse model about the mutual coupling of ${{\boldsymbol{U}}_s}$ can be represented aswhere $\bar {\boldsymbol{A}} = [{\boldsymbol{\varLambda}} ({\theta _1}),{\boldsymbol{\varLambda}} ({\theta _2}),\cdots ,{\boldsymbol{\varLambda}} ({\theta _N})] \in {{{\bf{C}}}^{M \times N(2P - 1)}}$ is the overcomplete dictionary for the block sparse model, and ${\boldsymbol{\varDelta }}$ is represented by ${\boldsymbol{\bar \varDelta }}$ in the sparse model, ${\boldsymbol{\bar \varDelta }} \in {{{\bf{C}}}^{N\tilde Q \times N}}$, and${{{\bar {\boldsymbol{B}}}}_s}$ is the representation of B in the sparse model, and whether the elements contained in the nth row of ${{{\bar {\boldsymbol{B}}}}_s}$ are all 0, which is contingent on the presence of a signal incident from the DOA corresponding to the grid point, $ n = 1,2,\cdots,{\text{ }}N $. In addition, ${{\bar {\boldsymbol{B}}}} = {{\bar{\boldsymbol{ \varDelta}} }}{{{\bar {\boldsymbol{B}}}}_s} = [{{\bar {\boldsymbol{B}}}}_1^\text{T},{{\bar {\boldsymbol{B}}}}_2^\text{T}, \cdots , {{\bar {\boldsymbol{B}}}}_N^\text{T}]^\text{T} \in {{{\bf{C}}}^{N\tilde Q \times Q}}$ is the block sparse model, ${\bar {\boldsymbol{B}}_n} = {{\varOmega }}({\theta _n}){\boldsymbol{\varGamma}} ({\theta _n}) \otimes {\bar {\boldsymbol{B}}_{s(n,:)}} \in {{{\bf{C}}}^{\tilde Q \times Q}}$ is the nth block matrix of ${{\bar {\boldsymbol{B}}}}$, $ n = 1,2,\cdots,N $. Therefore, whether the elements of ${\bar {\boldsymbol{B}}_n} = \bar {\boldsymbol{B}}((n - 1)\tilde Q + 1, n\tilde Q,:)$ are all 0 depends on whether the element of ${\bar {\boldsymbol{B}}_{s(n,:)}}$ are all 0, where $\bar {\boldsymbol{B}}((n - 1)\tilde Q + 1,n\tilde Q,:)$ are the elements from row $(n - 1)\tilde Q + 1$ to row $n\tilde Q$ of ${{\bar {\boldsymbol{B}}}}$, ${\bar {\boldsymbol{B}}_{s(n,:)}}$ represents the elements of the nth row of the matrix ${\bar {\boldsymbol{B}}_s}$.

According to (15), ${{\bar {\boldsymbol{B}}}}$ is the block sparse signal, and the index n of ${{\bar {\boldsymbol{B}}}}$’s block matrix ${{{\bar {\boldsymbol{B}}}}_n}$ corresponds to ${\theta _n}$. Therefore, after recovering the block sparse signal ${{\bar {\boldsymbol{B}}}}$, solve for all ${\left\| {{{{{\bar {\boldsymbol{B}}}}}_n}} \right\|_F}$ and find the first Q maximum values to obtain the true DOA estimation. From [10], we know that the most direct sparse recovery method for ${{\bar {\boldsymbol{B}}}}$ is to use ${l_0}$ ${\text{norm}}$ penalty, i.e.where ${\bar {\boldsymbol{B}}^{{l_F}}} = {\left[{\left\| {{{\bar {\boldsymbol{B}}}_1}} \right\|_F},{\left\| {{{\bar {\boldsymbol{B}}}_2}} \right\|_F},\cdots ,{\left\| {{{\bar {\boldsymbol{B}}}_N}} \right\|_F}\right]^\text{T}}$, $\xi $ is the threshold parameter, which determines the upper bound of the fitting error. However, the ${l_0}$ ${\text{norm}}$ question is an NP-hard problem, which is a challenging and intractable optimization problem. Therefore, the following work focuses on solving the problem with greater efficiency. $\parallel {{{\bar {\boldsymbol{B}}}}^{{l_F}}}{\parallel _0}\;$ can be represented aswhere the function expression for $\mathcal{I}(\alpha ) $ isconsidering the non-smooth function of $\parallel {\bar {\boldsymbol{B}}^{{l_F}}}{\parallel _0}$, this paper employs an approximation technique for $\parallel {\bar {\boldsymbol{B}}^{{l_F}}}{\parallel _0}$ by utilizing a series of exponential functions. The expression of the exponential function is as follows:so we can get $ \lim_{\sigma \to 0} {f_\sigma }(\alpha ) = 1 - \mathcal{I}(\alpha )$, when $\alpha \ne 0$, and

So

When $\sigma \to 0$, then there is

According to (21), (17) can be represented as

Our upcoming research will focus on solving (24). By solving this Lagrangian function, we can obtainwhere ${\bar {\boldsymbol{B}}_ * }(\sigma )$ is the optimal solution that satisfies $\bar{\boldsymbol{B}}_*(\sigma)=\arg \min _{\bar{B}} L_\sigma(\bar{\boldsymbol{B}})$. The parameter $\lambda $ regulates the balance between the signal’s sparsity and the residual energy [27], which is a constant. The choice of the value of $\sigma $ will have different effects: the smaller $\sigma $, the better approximation of $N - {F_\sigma }(\bar {\boldsymbol{B}})$ to $\parallel {\bar {\boldsymbol{B}}^{{l_F}}}{\parallel _0}$, but there is negative effect of many local minima; in addition, the larger $\sigma $, the smoother $N - {F_\sigma }(\bar {\boldsymbol{B}})$. To regulate this process, a set of gradually decreasing $\sigma $ is selected until the iteration termination condition is satisfied [28], i.e., the value of $\sigma $ gradually decreases from large to small. Additionally, Mohimani et al. [28] demonstrated thatwhere $\bar {\boldsymbol{B}}_ * ^{(0)}(\sigma )$ is the optimal initial value in the iterative solution (24).

It should be emphasized that (24) cannot use the Newton method because the Hessian matrix of $L_\sigma(\bar{\boldsymbol{B}})$ cannot always guarantee to be positive definite [29]. In addition, the calculation of $L_\sigma(\bar{\boldsymbol{B}})$’s Hessian matrix is very complex when the block sparse signal $\bar {\boldsymbol{B}}$ is no longer a real-valued column vector. For example, when $\bar {\boldsymbol{B}} \in {{{\bf{C}}}^{N\tilde Q \times Q}}$, the Hessian matrix dimension of the real-valued scalar ${L_\sigma }(\bar {\boldsymbol{B}})$ is $2N\tilde QQ \times 2N\tilde QQ$, the computational complexity of the inverse matrix of this Hessian matrix is extremely large. Therefore, this paper develops a Newton-like method based on mapping ${\xi }(\bar {\boldsymbol{B}})$ to solve (24), which can minimize computational complexity as much as possible.

Lemma 1　Define the mapping $\zeta: {{\bf{C}}}^{N \tilde{Q} \times Q} \rightarrow {{\bf{C}}}^{N \tilde{Q} \times Q}$, i.e.,

In reality, $\zeta (\bar {\boldsymbol{B}})$ is the expression for $\bar {\boldsymbol{B}}$ when ${\boldsymbol{G}}(\bar {\boldsymbol{B}})=0$, where ${\boldsymbol{G}}(\bar {\boldsymbol{B}})$ is the gradient of $ L_\sigma(\bar{\boldsymbol{B}}) $ at $\bar {\boldsymbol{B}}$:the specific derivation process of ${\boldsymbol{G}}(\bar {\boldsymbol{B}})$ can be found in Proof 1. ${\boldsymbol{W}}(\bar {\boldsymbol{B}})$ is a diagonal matrix:and ${{\boldsymbol{W}}_n}(\bar {\boldsymbol{B}}) = {f_\sigma }(\parallel {\bar {\boldsymbol{B}}_n}{\parallel _F}){{\boldsymbol{I}}_{\tilde Q}},n = 1,2,\cdots ,N$.

In addition, for any $\bar {\boldsymbol{B}}$, there exists a scalar real number $\kappa \geqslant 0$ that satisfies

In reality, (30) reveals that the value of ${L_\sigma }(\bar {\boldsymbol{B}})$ along the $\zeta (\bar {\boldsymbol{B}}) - \bar {\boldsymbol{B}}$ direction gradually decreases until the iteration termination condition is satisfied. The rationality of direction $\zeta (\bar {\boldsymbol{B}}) - \bar {\boldsymbol{B}}$ is given in as follows.

Proof 1　According to (24), the Lagrangian function can be obtainedand set ${L_{1\sigma }}(\bar{\boldsymbol{ B}}) = {F_\sigma }(\bar {\boldsymbol{B}})$, ${L_{2\sigma }}(\bar {\boldsymbol{B}}) = \parallel {{\boldsymbol{U}}_s} - \bar {\boldsymbol{A}}\bar {\boldsymbol{B}}\parallel _F^2$. Taking the derivative of $\bar {\boldsymbol{B}}$:where $\bar {\boldsymbol{B}}[i,j]$ is the element of the jth column and ith row of $\bar {\boldsymbol{B}}$, where $i \in (n - 1)\tilde Q + 1,\cdots ,n\tilde Q$, sothe specific expression of ${\boldsymbol{W}}(\bar {\boldsymbol{B}})$ can be found in (29).

Based on the above discussion,where ${\boldsymbol{G}}(\bar {\boldsymbol{B}})$ is the gradient of $L(\bar {\boldsymbol{B}})$ at $\bar {\boldsymbol{B}}$, considering ${\boldsymbol{G}}(\bar {\boldsymbol{B}}) = 0$, so we can get$\zeta (\bar {\boldsymbol{B}}) - \bar {\boldsymbol{B}}$ can be a descent direction for ${L_\sigma }(\bar {\boldsymbol{B}})$ to solve $\bar {\boldsymbol{B}}$, that is becausei.e., $\zeta (\bar {\boldsymbol{B}}) - \bar {\boldsymbol{B}} = - {\left[\dfrac{{{\boldsymbol{W}}(\bar {\boldsymbol{B}})}}{{{\sigma ^2}}} + 2\lambda {\bar {\boldsymbol{A}}^\text{H}}\bar {\boldsymbol{A}}\right]^{ - 1}}G(\bar {\boldsymbol{B}})$.

The specific proof process can be found in Proof 2.

In addition, define the matrixdue to $\dfrac{{{\boldsymbol{W}}(\bar {\boldsymbol{B}})}}{{{\sigma ^2}}} + 2\lambda {\bar {\boldsymbol{A}}^\text{H}}\bar {\boldsymbol{A}}$ is a positive definite matrix [29], so

This indicates that the inner product of $\zeta (\bar {\boldsymbol{B}}) - \bar {\boldsymbol{B}}$ and ${\boldsymbol{G}}(\bar {\boldsymbol{B}})$ is less than 0, so the direction of $\zeta (\bar {\boldsymbol{B}}) - \bar {\boldsymbol{B}}$ is opposite to the direction of the gradient ${\boldsymbol{G}}(\bar {\boldsymbol{B}})$. Therefore, as it changes along the direction $\zeta (\bar {\boldsymbol{B}}) - \bar {\boldsymbol{B}}$, the value of function $L(\bar {\boldsymbol{B}})$ gradually decreases until it satisfies the threshold for terminating the iteration. □

Proof 2　According to (27), we can obtainso that

Subtract $\left[ {\dfrac{{{\boldsymbol{W}}(\bar {\boldsymbol{B}})}}{{{\sigma ^2}}} + 2\lambda {{\bar {\boldsymbol{A}}}^\text{H}}\bar {\boldsymbol{A}}} \right]\bar {\boldsymbol{B}}$ from both sides of the equation at the same time, then there is

According to the expression of ${\boldsymbol{G}}(\bar {\boldsymbol{B}})$, i.e., (42), then there is

So□

Therefore, by using a backtracking algorithm and designing corresponding step sizes and parameters, the required solution $\bar {\boldsymbol{B}}$ can be obtained. The process of Newton-like method is recorded in Method 1.

The descent direction $\zeta (\bar {\boldsymbol{B}}) - \bar {\boldsymbol{B}}$ is very similar to the descent direction of the Newton method becausewhere ${{{\boldsymbol{W}}(\bar {\boldsymbol{B}})}}/{{{\sigma ^2}}} + 2\lambda {\bar {\boldsymbol{A}}^\text{H}}\bar {\boldsymbol{A}}$ is a fully-rank positive definite matrix, which can be seen as a Hessian-like matrix. The specific derivation process of (45) can be found in Proof 2. In reality, when $\bar {\boldsymbol{B}}$ is a real-valued column vector, ${{{\boldsymbol{W}}(\bar {\boldsymbol{B}})}}/{{{\sigma ^2}}} + 2\lambda {\bar {\boldsymbol{A}}^\text{H}}\bar {\boldsymbol{A}}$ is the Hessian matrix of $ L_\sigma(\bar{\boldsymbol{B}}) $’s positive definite part [30].

Remark　The computational burden of the proposed method comes from $\zeta (\bar {\boldsymbol{B}})$, and the computational complexity of $\zeta (\bar {\boldsymbol{B}})$ mainly comes from the inverse matrix ${\left[ {\dfrac{{{\boldsymbol{W}}(\bar {\boldsymbol{B}})}}{{{\sigma ^2}}} + 2\lambda {{\bar {\boldsymbol{A}}}^\text{H}}\bar {\boldsymbol{A}}} \right]^{ - 1}}$, whose computational complexity is ${{O}}({(N\tilde Q)^3})$, this computational complexity is enormous.

To reduce it, due to ${{\boldsymbol{U}}_s} = \bar {\boldsymbol{A}}\bar {\boldsymbol{B}}$, (45) can be represented as

Define $\tilde {\boldsymbol{Z}} = \left[ {\dfrac{{{\boldsymbol{W}}(\bar {\boldsymbol{B}})}}{{{\sigma ^2}}} + 2\lambda {{\bar {\boldsymbol{A}}}^\text{H}}\bar {\boldsymbol{A}}} \right]$, according to the inverse matrix lemma, ${\tilde {\boldsymbol{Z}}^{ - 1}}$ can be represented as

Therefore, $\zeta (\bar {\boldsymbol{B}}) - \bar {\boldsymbol{B}}$ can be represented aswhere ${\boldsymbol{Z}} = \left[ {{{{{\boldsymbol{I}}_M}}}/({{2\lambda {\sigma ^2}}}) + \bar {\boldsymbol{A}}{{\boldsymbol{W}}^{ - 1}}(\bar {\boldsymbol{B}}){{\bar {\boldsymbol{A}}}^\text{H}}} \right] \in {{{\bf{C}}}^{M \times M}}$, so the computational complexity of ${{\boldsymbol{Z}}^{ - 1}}$ is ${{O}}({M^3})$, the computational complexity of ${{{\bar {\boldsymbol{A}}}}^{\rm{H}}}{{\boldsymbol{Z}}^{ - 1}}$ is ${{O}}((N\tilde Q){M^2})$; the computational complexity of ${{\boldsymbol{Z}}^{ - 1}}{{\boldsymbol{U}}_s}$ is ${{O}}({M^2}Q)$; the computational complexity of ${{{\bar {\boldsymbol{A}}}}^{\rm{H}}}{{\boldsymbol{Z}}^{ - 1}}$ is ${{O}}((N\tilde Q){M^2})$; the computational complexity of $\bar {\boldsymbol{A}}{{\boldsymbol{W}}^{ - 1}}(\bar {\boldsymbol{B}}){\bar {\boldsymbol{A}}^\text{H}}$ is ${{O}}((N\tilde Q){M^2})$; the computational complexity of ${{{\bar {\boldsymbol{A}}}}^{\rm{H}}}{{\boldsymbol{Z}}^{ - 1}}{{\boldsymbol{U}}_s}$ is ${{O}}((N\tilde Q)MQ)$. Due to the significantly higher computational complexity of ${{O}}((N\tilde Q){M^2})$ and ${{O}}((N\tilde Q)MQ)$ compared to other terms, thus the proposed method’s computational complexity is ${{O}}((N\tilde Q){M^2})$+${{O}}((N\tilde Q)MQ)$.

The primary mission of this section is to verify the performance of the proposed method and compare it with other methods. The selected methods for comparison include the method of sparse representation (SR) [21], the method of block sparse representation (BSR) [22], the method of robust weighted subspace (RWS) [24], focal underdetermined system solver (FOCUSS) [31], the maximum number of iterations for the FOCUSS method is 30, with a termination threshold of 10^-3 and a norm factor of 1, weighted block sparse DOA estimation based on signal subspace (WBSS) [14], and Cramer-Rao lower bound (CRLB) [32]. In addition, root mean square error (RMSE) is widely used in parameter estimation, and the DOA estimation accuracy of different methods can be compared using RMSE, which is defined aswhere $\tilde M$ represents the number of Monte Carlo experiments, ${\hat \theta _{\tilde m,q}}$ denotes the calculated ${\theta _q}$ by the $\tilde m{\mathrm{th}}$ Monte Carlo experiment. Probability of resolution (PR) can serve as an evaluation criterion for the DOA estimation’s stability:where the ${P_{\tilde m}}(\theta )$ indicates whether successful estimation of ${\theta _q}$ in the $\tilde m{\mathrm{th}}$ Monte Carlo experiment. If the DOA estimation of any signal is within the standard range, ${P_{\tilde m}}(\theta ) = 1$; otherwise, ${P_{\tilde m}}(\theta ) = 0$. The scenario for ${P_{\tilde m}}(\theta ) = 1$, i.e., successful DOA estimation, is

The 1st experiment compares the RMSE of all methods under different SNRs. The parameter settings are as follows: assuming two narrowband signals are incident on the array, i.e., Q=2. The incident angles of the signal are ${\theta _1} = - 13.1^\circ $ and ${\theta _2} = 43.1^\circ $. The number of elements in the array is M=10, and the mutual coupling coefficients is c = [1, 0.1545−0.4776i], i.e., P=2. The input SNR ranges from −8 dB to 8 dB, with an interval of 2 dB, and the number of snapshots T=200. The relevant parameters in Method 1 are as follows: $\lambda = 2$, ${\sigma _{\min }} = {10^{ - 4}}$, $\gamma = 0.3$, $\eta = 0.3$, $\rho = 0.5$. In the following experiments, unless otherwise specified, all experimental parameters are consistent with the first experiment. As shown in Fig. 2, the RMSE of our proposed method reduces gradually in the range of −8−8 dB and achieves optimal DOA estimation at 4 dB. SR, BSR, and RWS belong to the ${l_1}$ norm method, FOCUSS belongs to the reweighted ${l_2}$ norm method, and the method proposed in this paper belongs to the ${l_0}$ method. After preprocessing (eigenvalue decomposition), the influence of noise is significantly reduced. The simulation results show that in an environment with very little noise, the ${l_0}$ method can achieve better performance. Furthermore, the dynamic selection of the step size β ensures that the proposed method exhibits more stable and better performance. WBSS also belongs to the ${l_1}$ method, but its performance is slightly better than the proposed method when obtaining prior information, indicating that prior information can effectively assist sparse recovery. However, the computation time of WBSS is very long, as shown in Table 1. It should be noted that the 0.1° RMSE of the DOA is inevitable due to the 0.1° deviation between the incident signal and the set grid points, so the RMSE curve is no longer decreasing and there is a big gap between the CRLB and the RMSE of DOA, and with the increase of SNR, the gap becomes bigger. If we want to achieve more ideal RMSE, we need to set denser grid points. However, the calculation time will typically increase significantly as the grid points become denser because the computational complexity of sparse DOA methods is closely linked to the number of grid points. Please refer to the 8th experiment for the specific calculation time. It should be noted that in the proposed method, the estimation of the number of sources is extremely important because the acquisition of the signal subspace is based on the number of sources. If the estimation of the number of sources is inaccurate, it means that the signal subspace obtained is also inaccurate, and the mathematical model of the signal sparsity model is also incorrect, which has a significant negative impact on DOA estimation.

The 2nd experiment compares the RMSE of all methods under the different snapshots. Assuming that the number of snapshots is 100−1000, which interval is 100, the input SNR is 1 dB. As shown in Fig. 3, our proposed method performs equally well as the WBSS method and outperforms other algorithms, which achieves optimal DOA estimation within this grid setting at the number of snapshots about 300. As the number of snapshots increases, the RMSE of the DOA estimation of the proposed method becomes better and better. This is because as the number of snapshots increases, the sampling covariance matrix becomes closer to the theoretical covariance matrix, and ${{\boldsymbol{U}}_s}$ becomes closer to the theoretical signal subspace. Similar to the first experiment, due to the 0.1° deviation between the incident signal and the set grid points, so the RMSE curve is no longer decreasing when K=300.

The 3rd experiment compares the PR of all methods under the different SNRs. Please refer to the 1st experiment for specific experimental parameters. As shown in Fig. 4, the PR performance of the proposed method improves with increasing input SNR, and in all SNR ranges, except for WBSS, the PR of the proposed method is higher than other methods, this is because the proposed method has the highest PR. In addition, the PR curve of the RWS method exhibits a decreasing trend under high SNR, indirectly suggesting that there is potential for enhancing the stability of this method. Our proposed method achieves a probability of successful DOA estimation of 1 at 4 dB. These phenomena indicate that the proposed DOA estimation method under mutual coupling has better stability under the different SNRs.

The 4th experiment compares the PR of all methods under the different snapshots. Please refer to the 2nd experiment for specific experimental parameters. As shown in Fig. 5, the PR performance of the proposed method improves with the increase in the number of snapshots; moreover, except for WBSS, under any number of snapshots, the proposed DOA estimation method under mutual coupling outperforms other methods in PR, highlighting the DOA estimation stability of the proposed method in different snapshots.

The 5th experiment aims to verify the effectiveness of the proposed method under various numbers of mutual coupling coefficients. The number of mutual coupling coefficients is set to P=1, i.e., c = [1]; P=2, i.e., c = [1, 0.1545−0.4776i]; P=3, i.e., c = [1, 0.1545−0.4776i, 0.1345+0.1070i], other parameter settings remain the same as in the 1st experiment. As shown in Fig. 6, when P=1, the performance is relatively optimal compared to P=2 and P=3; When P=3, the performance is relatively poor compared to P=2 and P=1. When SNR=2 dB, the number of mutual coupling coefficients reaches the optimal RMSE at P=1 and P=2. When SNR=8 dB, the number of mutual coupling coefficients approaches the optimal RMSE at P=3.

The 6th experiment compares the estimation performance under different DOAs. Set ${\theta _1}$ = 5.1°, ${\theta _2}$ ranging from 14.1° to 26.1°, with an interval of 3°, SNR = 5 dB, and K=200. As shown in Fig. 7, when ${\theta _2}$=20.1°, the RMSE of ${\theta _1}$ and ${\theta _2}$ are both close to 1°; when ${\theta _2}$=23.1°, the RMSE of ${\theta _1}$ and ${\theta _2}$ reaches the optimal RMSE. The simulation shows that the proposed method is more suitable for big DOA intervals.

The 7th experiment compares the computation time of all algorithms. Except for a fixed SNR of 4 dB, all other experimental conditions are consistent with the 1st experiment. The computer settings used are as follows: Intel Core i5-12500, 3.0 GHz processor, 16 GB RAM. As shown in Table 1, the proposed method and FOCUSS are significantly ahead of SR, BSR, and RWS in computational time. This is because SR, BSR, RWS, FOCUSS, and our proposed method require multiple iterations, but our method and FOCUSS have lower computational complexity, resulting in shorter calculation times. The computational complexity of FOCUSS is comparable to that of the proposed method, but the computational speed of FOCUSS is slightly slower than that of the proposed algorithm. This is because the calculation of FOCUSS requires setting corresponding thresholds and maximum iteration times. This experiment shows that under the current parameter settings, the proposed algorithm has higher computational efficiency.

The 8th experiment compares the computation time of different methods under different grid sizes. Fixed SNR=4 dB, the grid sizes are 0.2°, 0.3°, 0.4°, 0.5°, and 1°, respectively. The other simulation conditions are the same as the 1st experiment. As shown in Fig. 8, even with different grids, the computation time of the proposed method is still comparable to that of FOCUSS, and the proposed method is slightly faster than FOCUSS. Compared to the SR, BSR, RWS, FOCUSS, and WBSS, the proposed method is the fastest.

In this paper, we propose a Newton-like method for recovering sparse signals and achieving DOA estimation under unknown mutual coupling. In the proposed method, a block sparse model about array mutual coupling is established based on the signal subspace. Secondly, a new gradient direction, similar to the Newton’s method, is designed based on the mapping $\zeta (\bar {\boldsymbol{B}})$ and uses matrix inverse lemma to reduce its computational complexity. Finally, the advantage of the proposed method in DOA estimation accuracy is demonstrated through simulation.