Direction of arrival estimation with off-grid target based on sparse reconstruction in impulsive noise

Abstract

This paper proposes a direction of arrival estimation based on sparse signal reconstruction in the presence of alpha noise by the off-grid orthogonal matching pursuit algorithm. Assuming Gaussian distribution as the noise model, the previous sparse reconstruction enhances robustness by utilizing the least square criterion-based direction of arrival estimation algorithms. However, they severely degrade, even invalid, when there is thick trailing and large impulse in the noise molded by $\mathbf{\alpha }$ stable distribution. In addition, due to the discretization of the potential angle space, the accuracy of these methods will be reduced when the target is not completely on the divided mesh exactly. Increasing the grid density to improve estimation effect will increase the computation burden. The compressed sensing signal model is reconstructed by reshaping the fractional lower order covariance matrix of the sensor array received signal. Based on the reshaped signal model, the novel reshaped orthogonal matching pursuit algorithm and reshaped off-grid orthogonal matching pursuit algorithm are derived. Compared to the least square criterion method with Gaussian distribution assumption, the proposed algorithm obtains high-resolution direction of arrival estimation in $\mathbf{\alpha }$ noise. Moreover, the reshaped off-grid orthogonal matching pursuit algorithm improves the direction of arrival estimation accuracy with off-grid target. Numerical simulation results demonstrate the effectiveness of proposed method in direction of arrival estimation with off-grid targets in $\mathbf{\alpha }$ noise.

Introduction

Direction of arrival (DOA) estimation of multi-source signals is an important issue in the field of array signal processing, which is widely used in radar, sonar, wireless communications, and other fields. Many methods have been proposed in previous studies, but most of these methods are based on the assumption that the observed signal is mixed with additive noise, which takes Gaussian distribution as mathematical model. When the statistical properties of Gaussian noise are known, methods based on second-order moments and other algorithms based on the least square (LS) criterion can be used. These kinds of algorithms can effectively suppress Gaussian noise. The algorithm description is intuitive, and the process of calculation is simple.

Unfortunately, the assumption that the additive noise obeys Gaussian distribution is not universal, for there are plenty of other types of noise in nature, and impulsive noise is an important component of them. Research shows that impulsive noise is widely present in artificial audio environment, underwater environment, and various types of electromagnetic environments.^1–3 Furthermore, in many applications, the non-Gaussian noises are expressed as $\alpha$ stable distribution, such as electrical equipment and digital circuit switching noise, a variety of mechanical processing equipment’s electric spark noise, and the wireless communication harass signal noise. So signal processing technology under $\alpha$ noise has practical significance. However, it is restricted in signal processing under $\alpha$ noise because of the following properties:

There is no closed probability density function (PDF) of $\alpha$ stable distribution, so the PDF-based algorithms are unavailable, such as maximum-likelihood method ⁴ and Bayesian estimation algorithm.

Because there are no second-order moment and higher-order moment in $\alpha$ stable distribution, and the second-order moment-based algorithm, ⁵ high-order cumulant-based algorithm ⁶ and the LS criterion-based algorithm cannot be used to deal with Gaussian noise.

In order to solve the above problems, the method of fractional low-order moment is proposed in Ma and Nikias. ⁷ Due to the contribution of Ma and Nikias, ⁷ the existing subspace-based algorithm, such as Multiple Signal Classification (MUSIC) ⁸ and Signal Parameters via Rotational Invariance Techniques (ESPRIT), ⁹ can be used in the presence of $\alpha$ noise. In Zhong et al., ¹⁰ a novel DOA estimation method in impulsive noise is presented by applying particle filtering. In Zeng et al., ¹¹ a $l_p$ -MUSIC algorithm robust to $\alpha$ noise is proposed based on the $l_p$ norm criterion. The DOA estimation approach on spatial smoothing for coherent signals in impulsive noise is proposed in Li and Lin. ¹² However, the subspace-based algorithms are computationally expensive. For the last few years, many scholars have noticed the application of sparse reconstruction model in DOA estimation. Assuming the true DOAs exactly lie on the predefined grid points, though dividing the spatial angle with a predefined set of grid points densely enough, still not all is rational. So true DOAs of signal source in spatial angle are sparse, and sparse reconstruction technique to address DOA estimation problem can be adopted.

Sparse signal reconstruction, as one of the important research fields of Compressed Sensing (CS), ¹³ attracts the attention of many scholars in the last decade. Several approaches are presented in articles. Original Sparse reconstruction model generates a $l_0$ norm optimization problem, which is a non-deterministic polynomial hard (NP-Hard) problem. An approach to avoid NP-Hard problem by converting the $l_0$ norm optimization to $l_1$ norm optimization is proposed in Donoho. ¹⁴ It is a prominent contribution to provide a convex optimization approach for sparse reconstruction. But the $l_1$ norm optimization is a suboptimal approach, which approximates globally optimal solution, when it satisfies certain conditions. In Emadi et al., ¹⁵ orthogonal matching pursuit (OMP)-based DOA estimation is proposed as a solution to the $l_0$ norm optimization problem. In addition, sparse reconstruction with noise is another important topic. Based on the assumption of Gaussian noise, many approaches designed to enhance robustness are proposed (e.g. Basis Pursuit De-Noising (BPDN) algorithm ¹⁶ and Danzig Selector algorithm ¹⁷ ). Based on the assumption of $\alpha$ noise, Zeng et al. ¹⁸ devised a new definition of correlation called $l_p$ -correlation and developed $l_p$ -MP, weak $l_p$ -MP, and $l_p$ -OMP algorithms for robust sparse approximation. In Dai and So, ¹⁹ an off-grid Bayes-optimal algorithm based on posterior probability is proposed for robust DOA estimation with $\alpha$ noise.

Notably, the above methods do not preprocess the measured data to suppress noise. Although these algorithms can estimate targets in the presence of noise, it leads to reduced accuracy, when the signal-noise ratio (SNR) is at a lower level. In order to improve the accuracy, Ma et al. ²⁰ suppress Gaussian noise by second-order statistics and developed Khatri–Rao (KR) subspace approach, meanwhile Li et al. ²¹ suppresses $\alpha$ noise by fractional lower order statistics and reconstructed the CS signal model based on KR subspace approach. The convex optimization method, namely a suboptimal approach based on $l_1$ norm optimization, is used to solve the sparse reconstruction problem in Li et al. ²¹

In this paper, the CS signal model is reconstructed by reshaping the fractional lower order covariance (FLOC) matrix based on KR subspace approach and devises two OMP algorithms, namely reshaped orthogonal matching pursuit (R-OMP) and reshaped off-grid orthogonal matching pursuit (R-OGOMP), to address DOA estimation problem with $\alpha$ noise. Particularly, the R-OGOMP algorithm has three notable advantages. First, noise suppression performed by FLOC improves the accuracy when SNR is low. Second, OMP is an approach for $l_0$ norm optimization, which is a NP-Hard problem, superior to a suboptimal approach for $l_1$ norm optimization in precision. Third, off-grid framework formed by first-order Taylor expansion effectively deals with the targets, which are not exactly on the predefined grids without significant increase in calculations. In addition, to simplify the numerical calculation, the differential property of KR product is summarized and proved.

The symbols mentioned in the text are explained as follows: the $|•|$ indicates the absolute value of scalar; the $\Vert •\Vert$ indicates the norm of a vector; the superscripts ${(•)}^T$ , ${(•)}^{*}$ , and ${(•)}^H$ indicate transpose, conjugate, and conjugate transpose, respectively; the $\odot$ , $\otimes$ , and ${\otimes }_{KR}$ indicate Hadamard product, Kronecker product, and KR product.

$\alpha$ Stationary distribution

The PDF for the $\alpha$ stable distribution does not have closed-form. Therefore, the $\alpha$ stable distribution is usually represented by characteristic function as follows ¹

$$\varphi (t)=\exp \left\{\mathit{jat}-\gamma {\left|t\right|}^{\alpha }\left[1+j\beta \mathrm{sgn}(t)\omega \left(t,\alpha \right)\right]\right\}$$

(1)where

$$\omega \left(t,\alpha \right)=\{\begin{array}{cc}\tan \frac{\alpha \pi }{2}& \alpha \ne 1\\ \frac{2}{\pi }\log \left|t\right|& \alpha =1\end{array}$$

(2) $\mathrm{sgn}\left(t\right)$ is the sign function

$$\mathrm{sgn}\left(t\right)=\{\begin{array}{l}1t>0\hfill \\ 0t=0\hfill \\ -1t<0\hfill \end{array}$$

(3)

$$-\infty <a<\infty ,\gamma >0,0<\alpha \le 2,-1\le \beta \le 1$$

(4)

The important parameter in the equation is the characteristic exponent $\alpha$ . The value of $\alpha$ implies thickness in tails of distribution. The smaller the value, the greater the probability of random variables appearing far away from the center of the distribution, and the larger the noise impulse amplitude. When $\alpha =2$ , the $\alpha$ stable distribution has the same PDF as the PDF of Gaussian distribution; when $\alpha =1$ , $\beta =0$ , the $\alpha$ stable distribution performs the same statistical property as Cauchy distribution. So the $\alpha$ stable distribution is considered as a generalized Gaussian distribution. Both Gaussian distribution and Cauchy distribution are particular cases of the family.

Figure 1 shows the time domain stochastic process in $\alpha$ stable distribution and the probability distribution of $\alpha$ stable distribution with different $\alpha$ .

Figure 1.

$\alpha$ Stable distribution.

$a$ is the location parameter, which implies mean of the distribution when $1<\alpha \le 2$ , or median when $0<\alpha <1$ ; $\beta$ is the symmetry parameter; $\gamma$ is the scale parameter, also called the dispersion. For study convenience, the standard symmetric distribution model of the $\alpha$ stable distribution in signal processing is adopted. The standard symmetric distribution is defined by confirming $a=0$ , $\beta =0$ , $\gamma =1$ , $0<\alpha <2$ , and denoted by $S\alpha S$ . The noise involved in the numerical simulation in this paper is modeled as $S\alpha S$ .

An important property of the $\alpha$ stable distribution is as follows

If $0<\alpha <2$ , then

$E{\left|X\right|}^p=\infty$ , when $p\ge \alpha$

$E{\left|X\right|}^p<\infty$ , when $0\le p<\alpha$

If $\alpha =2$ , then

$E{\left|X\right|}^p<\infty$ , for all $p\ge 0$

The consequence of the formulas above is that the second-order or higher moment of $\alpha$ stable distributions tends to infinite, when $0<\alpha <2$ . Accordingly, the second-order moment method and high-order cumulant method widely used to suppress Gaussian noise are invalid in $\alpha$ noise. The LS criterion method of sparse signal reconstruction in CS is invalid all the same in $\alpha$ noise.

Compressed sensing

The main research directions of CS are as follows

1. Sparse signal compression acquisition

The technology or equipment compresses high-dimensional sparse signals into low-dimensional non-sparse signals during the acquisition process directly.

2. The over-complete dictionary structure

The methods construct a matrix with columns larger than rows according to the principle of CS.

3. Sparse signal reconstruction

The approaches reconstruct high-dimensional signals from low-dimensional received signals and determine non-zero elements and their values.

In this paper, the sparse signal reconstruction algorithms are mainly used to realize DOA estimation. The sparse reconstruction mathematical model is expressed as

$$y=\Phi x$$

(5)where $y\in {ℂ}^{M\times 1}$ is received vector; $\Phi \in {ℂ}^{M\times N}$ ( $N\gg M$ ) is a matrix that obeys a condition known as RIP, ¹³ we name it measurement matrix; $x\in {ℂ}^{N\times 1}$ is input vector. Sparse reconstruction is a problem try to reconstruct the higher-dimensional vector $x$ from the lower-dimensional vector $y$ under the condition that $y$ and $\Phi$ are known. From the point of view of solving mathematical equation, this is an under-determined system, which cannot be solved. But when the sparse conditions are satisfied, formula (5) converts to a solvable constrained optimization problem

$$\min {\Vert x\Vert }_{0}s.t.y=\Phi x$$

(6)

In the presence of noise, the formula is reformulated as

$$\min {\Vert x\Vert }_{0}s.t.{\Vert y-\Phi x\Vert }_2\le \epsilon$$

(7)

However, $l_0$ norm optimization is the NP-Hard problem. To solve the problem, scholars suggest an effective solution named greedy algorithm. The OMP algorithm is the most typical method.

DOA estimation based on FLOC matrix reconstruction

Consider the case that there are K narrow-band far-field signal sources $s_1\left(t\right)$ , $s_2\left(t\right)$ $\dots$ $s_K\left(t\right)$ in the spatial with different DOAs ${\theta }_1$ , ${\theta }_2$ $\dots$ ${\theta }_K$ arriving at a uniform linear array (ULA) with M sensors. The received signals of each array element have additive noise $n_1\left(t\right)$ , $n_2\left(t\right)$ , $\dots$ , $n_M\left(t\right)$ . Without loss of generality, the noise satisfies the assumption that it is independent identically distributed and uncorrelated to the signal. The received signal vector of the array can be expressed as

$$X\left(t\right)=A\left(\theta \right)S\left(t\right)+N\left(t\right)$$

(8)where

$$\begin{array}{l}{X}_{M\times 1}\left(t\right)={\left[x_1\left(t\right),x_2\left(t\right),\dots ,x_M\left(t\right)\right]}^T\\ A_{M\times K}\left(\theta \right)=\left[a\left({\theta }_1\right),a\left({\theta }_2\right),\dots ,a\left({\theta }_K\right)\right]\\ S_{K\times 1}\left(t\right)={\left[s_1\left(t\right),s_2\left(t\right),\dots ,s_K\left(t\right)\right]}^T\\ N_{M\times 1}\left(t\right)={\left[n_1\left(t\right),n_2\left(t\right),\dots ,n_M\left(t\right)\right]}^T\end{array}$$

where $x_m\left(t\right)$ , $m=1,2,\dots ,M$ is the output value of the $m\text{th}$ array element; $a\left({\theta }_k\right)$ , $k=1,2,\dots ,K$ is the $k\text{th}$ steering vector, which can be expressed as

$$a\left({\theta }_k\right)={\left[1,e^{-j\frac{2\pi }{\lambda }d\mathrm{\hspace{0.17em}sin}{\theta }_k},\dots ,e^{-j\frac{2\pi }{\lambda }\left(M-1\right)d\mathrm{\hspace{0.17em}sin}{\theta }_k}\right]}^T$$

(9)where, $\lambda$ is the carrier signal wavelength; $d$ is the intersensor spacing.

If the noise item in formula (8) is $\alpha$ noise and the characteristic exponent $\alpha \in \left(0,2\right)$ , then the FLOC matrix can be expressed as ²¹

$$C=A\left(\theta \right){\Lambda }_s{A}^H\left(\theta \right)+\xi I_M$$

(10)where, $A\left(\theta \right)$ is the steering matrix in formula (8); ${\Lambda }_s=\mathit{diag}\left\{{\delta }_s\right\}$ is a diagonal matrix; the elements of main diagonal ${\delta }_s={\left[{\delta }_1^2,\dots ,{\delta }_K^2\right]}^T$ can be treated as the fractional low-order power of the signal source; $\xi$ can be treated as the power level of the additive $\alpha$ noise and $I_M$ denotes the $M\times M$ identity matrix; $\xi$ is the eigenvalues obtained by eigenvalue decomposition of FLOC matrix $C$ ; so the item $\xi I_M$ can be treated as a known vector.

The $\left(i,k\right)$ element of matrix $C$ can be defined as

$$C_{i,j}=E\left\{x_i\left(t\right){\left|x_k\left(t\right)\right|}^{p-2}{x}_k^{*}\left(t\right)\right\}$$

(11)where, $p$ denotes the fractional lower order exponent, which satisfies $1<p<\alpha <2$ and guarantees $C_{i,j}<\infty$ is bounded. When $p=2$ , the formula (10) is consistent with the covariance matrix in Gaussian noise. Move the item $\xi I_M$ to the left-hand side of formula (10), and then

$$\tilde{C}=A\left(\theta \right){\Lambda }_s{A}^H\left(\theta \right)$$

(12)

Where $\tilde{C}=C-\xi I_M$

Applying the vectorization operator in Ma et al. ²⁰ on formula (10), it is reshaped into CS form as follows

$$y_{\mathit{FLOM}}=\mathit{vec}\left(\tilde{C}\right)=\tilde{A}\left(\theta \right){\delta }_s$$

(13)where, $\mathit{vec}(•)$ represents the expansion of a matrix by columns into a vector; $\tilde{A}\left(\theta \right)$ is the KR product reshaped matrix, which is expressed as

$$\tilde{A}\left(\theta \right)=A^{*}\left(\theta \right)\underset{KR}{\otimes }A\left(\theta \right)=\left[a^{*}\left({\theta }_1\right)\otimes a\left({\theta }_1\right),\dots ,a^{*}\left({\theta }_K\right)\otimes a\left({\theta }_K\right)\right]$$

(14)where, $\tilde{A}\left(\theta \right)$ is $M^2\times K$ reshaped matrix and ${\delta }_s$ is $K\times 1$ equivalent source vector. Formula (13) is reformulated as

$$\tilde{y}=\tilde{A}\left(\theta \right){\delta }_s$$

(15)where $\tilde{y}=\mathit{vec}\left(C\right)$

Define a set $\tilde{\theta }=\left[{\tilde{\theta }}_1,{\tilde{\theta }}_2,\dots ,{\tilde{\theta }}_Q\right]$ obtained by meshing the spatial angle; the number $Q$ of the predefined set of grid should be greater than the number $K$ of the signal sources and the row number $M^2$ of reshaped matrix; define $\tilde{A}\left(\tilde{\theta }\right)=\left[a^{*}\left({\tilde{\theta }}_1\right)\otimes a\left({\tilde{\theta }}_1\right),\dots ,a^{*}\left({\tilde{\theta }}_Q\right)\otimes a\left({\tilde{\theta }}_Q\right)\right]$ is the over complete dictionary; $a^{*}\left({\tilde{\theta }}_1\right)\otimes a\left({\tilde{\theta }}_1\right)$ can be interpreted as the steering vector of the virtual array; $s={\left[s_1,s_2,\dots ,s_Q\right]}^T$ can be interpreted as the sparse signal source power vector in spatial which contains $K$ nonzero elements; the $j\text{th}$ nonzero element $s_j$ equals to ${\delta }_k^2$ , and the predefined angle ${\tilde{\theta }}_j$ equals to DOA of the signal source ${\theta }_i$ , $i=1,\dots ,K$ ; the nonzero element index of vector $s$ is the position index of the sequence of spatial angles and the nonzero element value is the power of the signal source.

From the above definition, the sparse $\tilde{y}$ can be formulated as

$$\tilde{y}=\tilde{A}\left(\tilde{\theta }\right)s$$

(16)

Through contrast, formulas (16) and (5) have the similar sparse structure and the DOA estimation problem can be addressed by applying the algorithms of sparse reconstruction. The problem is converted into the $l_0$ -norm optimization and expressed as

$$\min {\Vert s\Vert }_0,s.t.{\Vert \tilde{y}-\tilde{A}\left(\tilde{\theta }\right)s\Vert }_2\le \epsilon$$

(17)

OMP algorithm can be utilized to solve the $l_0$ -norm optimization. Thus, the DOA of the target signal source is determined according to the nonzero elements index of the vector $s$ .

R-OMP algorithm

The basic idea of OMP algorithm is to treat the observation vector as a linear combination of the basis of the overcomplete dictionary matrix. The nonzero elements in the sparse vector are the coefficients of the basis. The zero elements in the sparse vector mean that the coefficient of the basis is zero, so the vectors in the overcomplete dictionary corresponding to the zero elements have no contributions to the linear combination of the observation vector. Therefore, sparse reconstruction problem can be addressed as long as the basis and corresponding coefficients that contribute the most to the linear combination are found. In OMP algorithm, the contributions of vectors of overcomplete dictionary in the linear combination is determined by calculating the correlation coefficient between the observation vector and the vector of overcomplete dictionary. The vectors of the overcomplete dictionary of the first $K$ ones with the largest contribution are regarded as the basis. The position of its nonzero coefficient in sparse vector can be determined by the index of the basis. The value of the coefficient can be obtained by iterative calculation of residuals. In the process of numerical calculation, the correlation coefficient is transformed into the inner product of the observation vector and the overcomplete dictionary vector under the background of Gaussian noise. The algorithm needs to satisfy several preconditions:

The target number is known.

The inner product operation can be performed between the observed vector and overcomplete dictionary vector.

The target DOA falls exactly on the predefined grid.

The sparse reconstruction problem is addressed based on the reshaped model referred to formula (17) through OMP algorithm and name it R-OMP. The steps of R-OMP algorithm are summarized in Algorithm 1.

Algorithm 1. R-OMP.

Inputs: Observation vector $\tilde{y}$ , reshaped overcomplete dictionary $\tilde{A}$ , sparsity $K$ .

Initialize: Residue error vector $r^0=\tilde{y}$ , support set ${\Omega }^0=\varnothing$ .

Steps:

$\mathit{for}$ $k=1:K$

 Step1: Obtain the location index $i_k=\begin{array}{lll}\arg & \underset{1\le i\le M}{\max }& c_p\left(a_i,r^{k-1}\right)\end{array}$ ;

 Step2: Update support set ${\Omega }^k={\Omega }^{k-1}\cup i_k$ ;

 Step3: Solution of least square

    ${\rho }^k={\left(\tilde{A}{\left(:,\Omega \right)}^H\tilde{A}\left(:,\Omega \right)\right)}^{-1}\tilde{A}{\left(:,\Omega \right)}^H*\tilde{y}$ ;

 Step4: Update residue $r^k=\tilde{y}-\tilde{A}*{\rho }^k$ ;

end

Step5: $\widehat{x}\left(\Omega \right)=\rho$ ;

Outputs: The sparse signal $\widehat{x}$ and the residue error vector $r$ .

R-OGOMP algorithm

R-OMP algorithm mentioned in section “R-OMP Algorithm” requires that the target DOAs exactly fall on the divided grids before DOA can be accurately estimated by sparse reconstruction. Unfortunately, the target has a high probability of not falling on the grids actually. It reduces the estimation accuracy because of grid mismatch.

Therefore, a novel off-grid OMP algorithm is proposed based on FLOC matrix reconstruction named R-OGOMP. The algorithm realizes DOA estimation based on sparse reconstruction in $\alpha$ noise. And it also improves the estimation accuracy by utilizing off-gird approach. The basic idea is expounded as follows.

The first-order Taylor series expansion for $\tilde{A}\left(\theta \right)$ in formula (16) is performed^22,23

$$\tilde{y}=\left(\tilde{A}+\tilde{B}\cdot \Delta \beta \right)s$$

(18)where, $\tilde{A}$ is virtual steering matrix in formula (16); $\tilde{B}$ is $M^2\times Q$ mismatch matrix obtained by taking the partial derivative of $\tilde{A}\left(\theta \right)$ about $\tilde{\theta }$ ; $\beta ={\left[{\beta }_1,{\beta }_2,\dots ,{\beta }_Q\right]}^T$ denotes the mismatch deviation vector where $\beta =\tilde{\theta }-{\theta }_{\mathit{true}}$ ; ${\theta }_{\mathit{true}}$ denotes the true DOA; and $\Delta \beta =\mathit{diag}\left(\beta \right)$ is diagonal matrix operation. According to the literature, ²⁴ $u=\beta \odot s$ , $u$ and $s$ have the same sparse structure. So formula (18) is converted into the optimization problem

$$\begin{array}{l}\underset{s\in {ℂ}^Q,\beta \in {ℂ}^Q}{\min }{\Vert s\Vert }_0+{\Vert u\Vert }_0\\ {\begin{array}{ll}s.t.& \Vert y-\left(\tilde{A}s+\tilde{B}u\right)\Vert \end{array}}_2\le \epsilon \end{array}$$

(19)

As long as the sparse structure of the vector $s$ is determined, the mismatch deviation vector $\beta$ has the same sparse structure. Inspired by Zeng et al. ¹¹ and Liu et al. ²² , we apply alternating direction iterative method to obtain $s^{\left(p\right)}$ and $u^{\left(p\right)}$ . The mathematical model is expressed as followse 

$${\widehat{s}}^{\left(p\right)}=\arg \underset{x^{\left(p\right)}\in {ℂ}^Q}{\min }{\Vert s^{\left(p\right)}\Vert }_{0}s.t.{\Vert \left(y-{\tilde{B}}^{\left(P\right)}{u}^{\left(p\right)}\right)-{\tilde{A}}^{\left(P\right)}{s}^{\left(p\right)}\Vert }_2\le \epsilon$$

(20)

$${\widehat{u}}^{\left(p\right)}=\arg \underset{u^{\left(p\right)}\in {ℂ}^Q}{\min }{\Vert u^{\left(p\right)}\Vert }_{0}s.t.{\Vert \left(y-{\tilde{A}}^{\left(P\right)}{\widehat{s}}^{\left(p\right)}\right)-{\tilde{B}}^{\left(P\right)}{u}^{\left(p\right)}\Vert }_2\le \epsilon$$

(21)where, ${\widehat{s}}^{\left(p\right)}$ is the sparse construction estimation value of $s^{\left(p\right)}$ in the $p\text{th}$ iteration; and ${\widehat{u}}^{\left(p\right)}$ is the estimation value of $u^{\left(p\right)}$ . Through addressing the problem of formulas (20) and (21), the sparse signal ${\widehat{s}}^{\left(p\right)}$ and ${\widehat{u}}^{\left(p\right)}$ can be obtained. And ${\widehat{\beta }}^p$ can be obtained by utilizing the nonzero elements of ${\widehat{s}}^{\left(p\right)}$ to solve the problem ${\widehat{\beta }}_i^p=\mathrm{Re}\left({\widehat{u}}_i^p/{\widehat{s}}_i^p\right)$ for any ${\widehat{s}}_i^{\left(p\right)}\ne 0$ . Then update the following iteration operators

$$s^{\left(p+1\right)}={\widehat{s}}^{\left(p\right)}$$

(22)

$$u^{\left(p+1\right)}={\widehat{u}}^{\left(p\right)}$$

(23)

$${\tilde{\theta }}^{\left(p+1\right)}={\tilde{\theta }}^{\left(p\right)}+{\widehat{\beta }}^{\left(p\right)}$$

(24)

$${\tilde{A}}^{\left(P+1\right)}=\tilde{A}\left({\tilde{\theta }}^{\left(P+1\right)}\right)$$

(25)

$${\tilde{B}}^{\left(P+1\right)}={\tilde{A}}^{\prime \left(P+1\right)}$$

(26)

Because $\tilde{A}$ is given by the KR product, its matrix structure is complex. It is more complex to calculate the partial derivatives of matrix directly, because it takes up more computing resources and slows down computing time. The differential property of KR product is presented and proved.

Property 1: For an arbitrary matrix $A=\left[a_1,a_2,\dots ,a_N\right]$ , $A\in {ℂ}^{M\times N}$ , $a_i\in {ℂ}^{M\times 1}$ , $i=1,2,\dots ,N$ , and matrix $B=\left[b_1,b_2,\dots ,b_N\right]$ , $B\in {ℂ}^{M\times N}$ , $b_i\in {ℂ}^{M\times 1}$ , $i=1,2,\dots ,N$ , there is

$$\left(A\underset{KR}{\otimes }B\right)\prime =A\prime \underset{KR}{\otimes }B+A\underset{KR}{\otimes }B\prime$$

where, $(•)\prime$ denotes differential operation; both $A$ and $B$ are the function matrix of the parameter $\theta$ . Proof

: $a_i\left(\theta \right)$ and $b_j\left(\theta \right)$ is the function vector of parameter $\theta$ in $A$ and B respectively.

$$\begin{array}{l}\left(a_i\otimes b_j\right)\prime =\left({\left[a_{1,i}\cdot b_j^T,a_{2,i}\cdot b_j^T,\cdots ,a_{M,i}\cdot b_j^T\right]}^T\right)\prime \\ ={\left[\left(a_{1,i}\cdot b_j^T\right)\prime ,\left(a_{2,i}\cdot b_j^T\right)\prime ,\cdots ,\left(a_{M,i}\cdot b_j^T\right)\prime \right]}^T\\ ={\left[\left(a_{1,i}^{\prime }\cdot b_j^T+a_{1,i}\cdot b_j^{\prime T}\right),\left(a_{2,i}^{\prime }\cdot b_j^T+a_{2,i}\cdot b_j^{\prime T}\right),\cdots ,\left(a_{M,i}^{\prime }\cdot b_j^T+a_{M,i}\cdot b_j^{\prime T}\right)\right]}^T\\ ={\left[a_{1,i}^{\prime }\cdot b_j^T,a_{2,i}^{\prime }\cdot b_j^T,\cdots ,a_{M,i}^{\prime }\cdot b_j^T\right]}^T+{\left[a_{1,i}\cdot b_j^{\prime T},a_{2,i}\cdot b_j^{\prime T},\cdots ,a_{M,i}\cdot b_j^{\prime T}\right]}^T\\ =a_i^{\prime }\otimes b_j+a_i\otimes b_j^{\prime }\\ \left(A\underset{KR}{\otimes }B\right)\prime =\left[\left(a_1\otimes b_1\right)\prime ,\left(a_2\otimes b_2\right)\prime ,\cdots ,\left(a_N\otimes b_N\right)\prime \right]\\ =\left[\left(a_1^{\prime }\otimes b_1^{\hspace{0.17em}}+a_1^{\hspace{0.17em}}\otimes b_1^{\prime }\right),\left(a_2^{\prime }\otimes b_2+a_2\otimes b_2^{\prime }\right),\cdots ,\left(a_N^{\prime }\otimes b_N+a_N\otimes b_N^{\prime }\right)\right]\\ =\left[a_1^{\prime }\otimes b_1,a_2^{\prime }\otimes b_2,\cdots ,a_N^{\prime }\otimes b_N\right]+\left[a_1\otimes b_1^{\prime },a_2\otimes b_2^{\prime },\cdots ,a_N\otimes b_N^{\prime }\right]\\ =A\prime \underset{KR}{\otimes }B+A\underset{KR}{\otimes }B\prime \end{array}$$

The efficiency of matrix operation is improved by applying Property 1. $\tilde{B}$ in formula (18) derived from property 1 expresses as

$$\tilde{B}=A^{\prime *}\underset{KR}{\otimes }A+A^{*}\underset{KR}{\otimes }{A}^{\prime }$$

(27)

With the analysis above, the steps of R-OGOMP algorithm is summarized in Algorithm 2.

Algorithm 2. R-OGOMP.

Inputs: Fractional lower order exponent $p$ , observation vector $\tilde{y}$ , predefined grid vector $\tilde{\theta }$ , steering matrix $A$ , mismatch matrix $B$ , sparsity $K$ .

Initialize: Residue error vector $r^0={\widehat{r}}^0=\tilde{y}$ ,   support set ${\Omega }^0=\varnothing$ ,

$\widehat{u}=0$ , $\widehat{s}=0$ , $j=0$ , $k=0$ , ${\tilde{A}}^0=A^{*}\underset{KR}{\otimes }A$ and   ${\tilde{B}}^0=A^{\prime *}\underset{KR}{\otimes }A+A^{*}\underset{KR}{\otimes }{A}^{\prime }$

Steps:

  $\mathit{while}$ (stopping criteria is not satisfied)

  $j=j+1$ ;

  $\mathit{for}$ $k=1:K$ ;

  Step 1: Obtain the location index

     $i_k=\begin{array}{lll}\arg & \underset{1\le i\le M}{\max }& c_p\left(a_i,{\widehat{r}}^j\right)\end{array}$ ;

  Step 2: Update support set ${\Omega }^k={\Omega }^{k-1}\cup i_k$ ;

  Step 3: Solution of least square

     ${\rho }^k={\left({\tilde{A}}^j{\left(:,\Omega \right)}^H{\tilde{A}}^j\left(:,\Omega \right)\right)}^{-1}{\tilde{A}}^j{\left(:,\Omega \right)}^H*r^j$ ;

  Step 4: Update residue ${\widehat{r}}^j=r^j-{\tilde{A}}^j\left(:,\Omega \right)*{\rho }^k$ ;

end

Step 5: $\widehat{x}\left(\Omega \right)={\rho }^k$ ;

Step 6: Estimate ${\widehat{u}}^j$ : ${\widehat{u}}^j={\left({\tilde{B}}_{\Omega }^{jH}{\tilde{B}}_{\Omega }^j\right)}^{-1}{\tilde{B}}_{\Omega }^{jH}*{\widehat{r}}^j$ ;

Step 7: Estimate ${\widehat{\beta }}^j$ : ${\widehat{\beta }}_i^j=\mathrm{Re}\left({\widehat{u}}_i^j/{\widehat{s}}_i^j\right)$ ;

Step 8: Update angle: ${\tilde{\theta }}^{\left(j+1\right)}={\tilde{\theta }}^{\left(j\right)}+{\widehat{\beta }}^{\left(j\right)}$ ;

Step 9:Update iteration operators:

$s^{\left(j+1\right)}={\widehat{s}}^{\left(j\right)}$ ,

$u^{\left(j+1\right)}={\widehat{u}}^{\left(j\right)}$ ,

${\tilde{A}}^{\left(j+1\right)}=\tilde{A}\left({\tilde{\theta }}^{\left(j+1\right)}\right)$ ,

${\tilde{B}}^{\left(j+1\right)}={\tilde{A}}^{\prime \left(j+1\right)}$ ;

Step 10: Update $r,\widehat{r}$ : $r^{j+1}={\widehat{r}}^{j+1}=r^j-{\tilde{B}}_{\Omega }^j*{\widehat{u}}^j$ ;

end

Outputs: The sparse signal $s^{\left(j+1\right)}$ and $u^{j+1}$ , estimation of angle ${\tilde{\theta }}^{\left(j+1\right)}$ .

Then discuss the stopping criteria:

Stop the iteration steps after a fixed number $j=J$ , i.e. $j=100$ .

Run until the norm of the residual declines to a level $\epsilon$ . That is, ${\Vert y-{\tilde{A}}^{\left(P\right)}{s}^{\left(p\right)}-{\tilde{B}}^{\left(P\right)}{u}^{\left(p\right)}\Vert }_2\le \epsilon$ .

Halt the algorithm when residual error rate drops below a threshold.

Numerical simulations

In this section, some numerical experiments are provided to assess the performance of the R-OMP and R-OGOMP algorithms. The $\alpha$ stable distribution is taken to model the noise. The general signal-to-noise ratio (GSNR) is defined as

$$\mathit{GSNR}=10\log \left(\frac{1}{\gamma N}{{\displaystyle \sum _{n=1}^N\left|s\left(k\right)\right|}}^2\right)$$

(28)where, $\gamma$ is the scale parameter of $\alpha$ noise; $N$ is the snapshots; and $s\left(k\right)$ is the sampling values of the signal sources. The mean square error（MSE） of the DOAs estimation is defined as

$$\mathit{MSE}={\displaystyle \sum {\hspace{0.17em}}_{n=1}^N}{\displaystyle \sum {\hspace{0.17em}}_{k=1}^K}{\left({\widehat{\theta }}_k^n-{\theta }_k\right)}^2/NK$$

(29)where $K$ is the number of target DOAs; ${\widehat{\theta }}_k^n$ is the $k\text{th}$ DOA estimation on the $n\text{th}$ snapshot; and ${\theta }_k$ is the true DOA.

Figure 2 shows the measurement signal polluted by the impulsive noise and Figure 3 shows the measurement signal polluted by the Gaussian noise. It can be observed that the measurement signal in Figure 2 contains more impulse compared with Figure 3.

Figure 2.

Measurement signal polluted by the impulsive noise.

Figure 3.

Measurement signal polluted by the Gaussian noise.

Determination of the value of p

To determine the value of p, we design two simulations experiments to observe the effect on accuracy and success probability when the p value changes. First, Shao and Nikias ¹ point out that the definition of FLOC is valid when $1\le p<\alpha$ . Second, we consider the noise in the following simulations is $\alpha$ noise with the characteristic exponent $\alpha =1.5$ . So, p value should be taken in intervals $[1,1.5)$ .

The assumption of the experiment is that the numbers of sensors in the array and the snapshots are $M=8$ and $N=500$ , respectively. The intersensor spacing $d$ is half a wavelength. The DOAs of two spatial signal sources ( $K=2$ ) are picked at random in the interval ${\theta }_1\in \left[-20°,-10°\right]$ and ${\theta }_2\in \left[25°,35°\right]$ . The numbers of Monte Carlo trials are $M_c=200$ . The characteristic exponent $\alpha =1.5$ and $\mathit{GSNR}=5dB$ . For R-OGOMP algorithm, the predefined angle is $\mathit{grid}=1°$ .

Figure 4 shows that the performance of proposed algorithms on accuracy has not changed much when p value is 1.1, 1.2, 1.3, and 1.4. Figure 5 shows that the success probability of proposed algorithms is higher when p value is 1.2 and 1.3. Comprehensive consideration, we take $p=1.3$ as the default parameters in the following simulations.

Figure 4.

MSE of DOA estimation of two signal sources in $S\alpha S$ noise with $\alpha =1.5$ versus p.

Figure 5.

Success probability of DOA estimation of two signal sources in $S\alpha S$ noise with a = 1.5 versus p.

Multiple source DOA estimation in α noise

In the simulation, the effectiveness of proposed algorithms in $\alpha$ noise is demonstrated, and the resolution is compared. There are three signal sources in spatial and the DOAs, of them are $\theta =\left[\begin{array}{lll}-36.7°& 3.5°& 37.1°\end{array}\right]$ . The number of the sensors in the ULA is $M=8$ , and the intersensor spacing $d$ is half a wavelength. The snapshots are $N=200$ . For R-OGOMP algorithm, the predefined grid angle is $\mathit{grid}=2°$ . The noise is $\alpha$ noise with the characteristic exponent $\alpha =1.5$ and $\mathit{GSNR}=5dB$ . We take $p=1.3$ for FLOC-based methods. Figures 6 to 9 show the spatial spectra of four independent trials, respectively. Figures 6 and 7 show that the R-OMP and R-OGOMP algorithms can realize DOA estimation based on sparse reconstruction in $\alpha$ noise. The spatial spectrums in Figures 6 and 7 are sharper than those in Figures 8 and 9. It demonstrates that the proposed methods in this paper have higher resolution than $l_p$ -MUSIC and FLOM-MUSIC algorithms.

Figure 6.

Normalized spatial spectrum of R-OMP algorithm is $S\alpha S$ noise with $\alpha =1.5$ and $\mathit{GSNR}=5dB$ . The red circles show the true value of DOAs.

Figure 7.

Normalized spatial spectrum of R-OGOMP algorithm is $S\alpha S$ noise with $\alpha =1.5$ and $\mathit{GSNR}=5dB$ . The off-grid angle is $2^o$ . The red circles show the true value of DOAs.

Figure 8.

Normalized spatial spectrum of $l_p$ -MUSIC algorithm is $S\alpha S$ noise with $\alpha =1.5$ and $\mathit{GSNR}=5dB$ . The step-size is 0.1.

Figure 9.

Normalized spatial spectrum of FLOM-MUSIC algorithm is $S\alpha S$ noise with $\alpha =1.5$ and $\mathit{GSNR}=5dB$ . The step-size is 0.1.

Estimation accuracy comparison

Monte Carlo experiments are carried out to compare the performance of the proposed algorithms with MMV-OMP, FLOM-MUSIC, and $l_p$ -MUSIC algorithms on the accuracy. The numbers of sensors in the array and the snapshots are $M=8$ and $N=500$ , respectively. The intersensor spacing $d$ is half a wavelength. The DOAs of two spatial signal sources ( $K=2$ ) are picked at random in the interval ${\theta }_1\in \left[-20°,-10°\right]$ and ${\theta }_2\in \left[25°,35°\right]$ . $p=1.3$ is taken for R-OGOMP, Fractional Lower Order Moment-MUSIC(FLOM-MUSIC), and $l_p$ -MUSIC algorithms. The numbers of Monte Carlo trials are $M_c=200$ .

Figure 10 illustrates the MSE of the five methods decreases with the increase of GSNR. The characteristic exponent of $\alpha$ noise is $\alpha =1.5$ . For R-OGOMP algorithm, the predefined angle is $\mathit{grid}=2°$ . It can be seen that the conventional sparse reconstruction-based DOA estimation method named Multiple Measurement Vector-OMP (MMV-OMP) algorithm is not robust in $\alpha$ noise. FLOM-MUSIC and $l_p$ -MUSIC algorithms have an improved performance. R-OMP and R-OGOMP algorithms have the best performance compared with other methods in most cases. The R-OGOMP algorithm is superior to R-OMP algorithm, especially when $\mathit{GSNR}\ge -2dB$ .

Figure 10.

MSE of DOA estimation of two signal sources in $S\alpha S$ noise with $\alpha =1.5$ versus GSNR.

Figure 11 shows the trend of MSE of the five methods versus $\alpha \in \left[1.1,2\right]$ . The predefined angle of MMV-OMP, R-OMP, and R-OGOMP algorithms is $\mathit{grid}=1°$ , meanwhile the step size of FLOM-MUSIC and $l_p$ -MUSIC algorithms is $s=1$ . The GSNR is 5 dB. It shows that the MSE decreases with the increase of $\alpha$ and the R-OGOMP algorithm has the best performance on estimation accuracy. The estimation accuracy of MMV-OMP algorithm is lower than other four algorithms. The reason is the MMV-OMP is not robust in $\alpha$ noise. The performance of FLOM-MUSIC algorithm is better than R-OMP algorithm. When $\alpha =2$ , the noise model converts to Gaussian noise and the performance on accuracy of all the five algorithms is similar.

Figure 11.

MSE of DOA estimation of two signal sources in $S\alpha S$ noise versus $\alpha$ at $\mathit{GSNR}=5dB$ .

Success probability comparison

In this section, the Monte Carlo experiments are performed to compare the success probability of MMV-OMP, R-OMP, R-OGOMP, FLOM-MUSIC, and $l_p$ -MUSIC algorithms. The assumption of the experiment is that $M=8$ , $N=500$ , $K=2$ , ${\theta }_1\in \left[-20°,-10°\right]$ , ${\theta }_2\in \left[25°,35°\right]$ , and $p=1.3$ . It can be defined the numerical simulation is success, when the DOA estimated deviation sum of two signal sources is less than 1°.

Figure 12 shows the trend of success probability versus GSNR. Make $\alpha =1.5$ and $\mathit{GSNR}\in \left[-9,21\right]$ . The numerical results indicate conclusions as follows:

Figure 12.

Success probability of DOA estimation of two signal sources in $S\alpha S$ noise with $\alpha =1.5$ versus GSNR.

The MMV-OMP algorithm is not robust in $\alpha$ noise, so its success probability is the lowest at different GSNR. Without noise suppression, the success probability of the MMV-OMP algorithm is zero when $\mathit{GSNR}\le 0dB$ . Because of the $\alpha$ noise suppression, R-OGOMP, R-OMP, and FLOM-MUSIC algorithms have high success probability even when $\mathit{GSNR}\le 0dB$ . The performance of R-OGOMP algorithm on success probability is the best, especially when $\mathit{GSNR}\in \left[-9,0\right]$ .

Figure 13 shows the trend of success probability versus $\alpha$ . Make $\mathit{GSNR}=5dB$ and $\alpha \in \left[1.1,2\right]$ . The result indicates that the performance of R-OGOMP algorithm is the best and the success probability is 1 when $\alpha \ge 1.6$ . When $\alpha =2$ , the performance of the five algorithms is similar.

Figure 13.

Success probability of DOA estimation of two signal sources in $S\alpha S$ noise versus $\alpha$ at $\mathit{GSNR}=5dB$ .

Off-grid performance comparison

The statistical performance of R-OGOMP, R-OMP, FLOM-MUSIC, and $l_p$ -MUSIC algorithms on different predefined grid angles is compared in this section. The assumption of the experiment is that $M=8$ , $N=500$ , $K=2$ , ${\theta }_1\in \left[-20°,-10°\right]$ , ${\theta }_2\in \left[25°,35°\right]$ , and $p=1.3$ . It can be defined success, when sum of DOA estimation deviation is < 1°. The predefined grid angle is equivalent to step size of $l_p$ -MUSIC and FLOM-MUSIC algorithms.

Figure 14 shows that the MSE of R-OGOMP algorithm is significantly better than other algorithms as the predefined grid angle increases.

Figure 14.

MSE of DOA estimation of two signal sources in $S\alpha S$ noise with $\alpha =1.5$ and at $\mathit{grid}\in \left[1,10\right]$ Versus $\mathit{grid}\in \left[1,10\right]$ .

Figure 15 shows that the algorithm with highest success probability versus increasing predefined grid angle is R-OGOMP. The lowest one is $l_p$ -MUSIC algorithm.

Figure 15.

Success probability of DOA estimation of two signal sources in $S\alpha S$ noise with $\alpha =1.5$ and at $\mathit{GSNR}=10dB$ versus $\mathit{grid}\in \left[1,10\right]$ .

Conclusion

This paper proposes a new array manifold by reshaped FLOC matrix and the off-grid model is introduced. The DOA estimation problem is addressed by OMP algorithm based on sparse reconstruction of the reshaped vector, which are named R-OMP and R-OGOMP algorithms. Different from robust-based algorithms, accuracy performance is improved by noise suppression. Because of the off-grid model, the accuracy of DOA estimation is improved when there are off-grid targets in spatial and the number of predefine grid is reduced. In the numerical calculation, the differential property of KR product is summarized to reduce computational complexity and proved. Simulations results demonstrate that the proposed algorithm realizes DOA estimation in presence of $\alpha$ noise. The R-OGOMP algorithm outperforms several existing DOA estimation methods in terms of estimation accuracy and probability of resolution.