DOA Estimation Based on Sparse Representation of the Fractional Lower Order Statistics in Impulsive Noise
  IEEE/CAA Journal of Automatica Sinica  2018, Vol. 5 Issue(4): 860-868   PDF    
DOA Estimation Based on Sparse Representation of the Fractional Lower Order Statistics in Impulsive Noise
Sen Li, Rongxi He, Bin Lin, Fei Su     
Department of Information Science and Technology, Dalian Maritime University, Dalian 116026, China
Abstract: This paper is mainly to deal with the problem of direction of arrival (DOA) estimations of multiple narrow-band sources impinging on a uniform linear array under impulsive noise environments. By modeling the impulsive noise as fistable distribution, new methods which combine the sparse signal representation technique and fractional lower order statistics theory are proposed. In the new algorithms, the fractional lower order statistics vectors of the array output signal are sparsely represented on an overcomplete basis and the DOAs can be effectively estimated by searching the sparsest coefficients. To enhance the robustness performance of the proposed algorithms, the improved algorithms are advanced by eliminating the fractional lower order statistics of the noise from the fractional lower order statistics vector of the array output through a linear transformation. Simulation results have shown the effectiveness of the proposed methods for a wide range of highly impulsive environments.
Key words: α-stable distribution     direction of arrival (DOA)     fractional lower-order statistics     impulsive noise     sparse representation    
Ⅰ. INTRODUCTION

Direction of arrival (DOA) estimation of multiple emitting sources is an important issue in array processing and has various applications in military, radar, sonar, wireless communications and source localization [1], [2]. A large number of solutions have been proposed to solve this problem during the past years. Usually, these solutions can be categorized into three groups: time-delay based methods, beamforming methods and signal subspace methods. However, the majority of DOA estimation algorithms are developed under certain assumptions: the source signal needs to be statistically stationary and uncorrelated, the number of snapshots is sufficient, and the signal-noise ratio (SNR) is moderately high. Practically, these conditions are barely satisfied, thus these methods achieve the limited estimation accuracy. In order to increase the DOA estimation accuracy, the well-known subspace-based method of multiple signal classification (MUSIC) algorithm and estimation method of signal parameters via rotational invariance techniques (ESPRIT) have been widely used due to their high estimation accuracies but at the price of high complexity.

Recently, the sparse representation technique of signal has been applied in many areas, such as image processing, wireless channel estimation and biomedical signal processing, which also provides a new idea for DOA estimation based on the fact that the number of sources is in general much smaller than that of potential source points when implementing the array processing algorithms. Several DOA estimation methods based on sparse representation have been proposed [3]-[16]. In [3], [4], a whiten sparse covariance-based representation model is first presented for the source parameter estimation by applying the global matched filter (GMF). In [5] the most representative sparse recovery algorithm for DOA estimation ($l_1$-SVD) is proposed, which can effectively estimate DOA with single measurement. By using the singular value decomposition (SVD) of received data matrix, it not only can work in multiple measurement cases but also can reduce the computational complexity. Although the $l_1$-norm minimization is a convex problem and the global minima can be guaranteed easily, the weakness is their undemocratic penalization for larger coefficients, which results in the degradation of signal recovery performance. To conquer this problem, the iterative reweighted $l_1$ minimization is designed [6], [7], where the large weights can be used to discourage nonzero entries in the recovered signals. To improve the convergence rate and estimation accuracy of the $l_{2, 1}$-norm minimization approach, Wei et al. develop a novel greedy block coordinate descent (GBCD) algorithm by using a greedy strategy for choosing descent directions [8]. In [9], a mixed $l_{2, 0}$-norm based joint sparse approximation technique is introduced into DOA estimation where the $l_0$ norm constraint is approached by a set of convex functions, and a method called JLZA-DOA is proposed. Algorithms in [5]-[9] address the DOA estimation problem by directly representing the array output in time domain with an overcomplete basis from the array response vector.

To make use of the second order statistics of the array output, a sparse iterative covariance-based estimation (SPICE) approach for array signal processing by the minimization of a covariance matrix fitting criterion which can be used in both single and multiple measurement cases is proposed in [10]. Another method called $l_1$-SRACV in [11] is also proposed for DOA estimation by using the array covariance matrix sparse representation which exhibits some merits of increased resolution. Because of recovering a joint-sparse inverse problem from multiple measurement vectors, the $l_1$-SRACV algorithm suffers from a high computational cost. Then a new DOA estimation method is proposed in [12], [13] which uses the combination of the Khatri-Rao product and sparse representation to estimate the DOAs of signals by recovering a sparse covariance vector of only a single measurement vector, thereby implying lower computational complexity than the $l_1$-SRACV algorithm. The authors of [14] firstly transform the multiple measurement vectors problem to the virtual single measurement vector (VSMV) problem in the sparse signal representation framework, and then exploit a surrogate truncated $l_1$ function to approximate $l_0$-norm, and successively demonstrate how the nonconvex minimization problem is treated by the difference of convex functions decomposition and the iterative approach. The study of [15] demonstrates why the multiple parameters can be exactly obtained by solving a weighted ``group lasso'' problem in the second-order statistics using a cross-dipole array. In [16], the DOA estimation of the wideband signal has been studied by the sparse representation of the covariance matrix.

All the above mentioned sparse representation based DOA estimation algorithms assume that the ambient noise is Gaussian distributed. However, the noise in practice often exhibits non-Gaussian properties, sometimes accompanied by strong impulsiveness [17]. For example, atmospheric noise (thunderstorms), car ignitions, microwave ovens, office equipments, and other types of naturally occurring or man-made signal sources can result in aggregating noise components that may exhibit high amplitudes for small time intervals. Under investigation, it is found that $\alpha$-stable distribution ($0 < \alpha\leq 2$) is a suitable noise model to describe this type of noise [18]. It can be considered as the greatest potential distribution to characterize various impulsive noises as different characteristic exponent parameters are selected.

An important characteristic of the $\alpha$-stable distribution is that only moments of order less than $\alpha$ exist. Therefore the performance of the DOA estimation algorithm based on the second order statistics of the array output will severely degrade in the presence of the $\alpha$-stable non-Gaussian noise. One way to alleviate this problem is to introduce new covariance estimates. Authors in [19] propose new subspace DOA estimation methods based on fractional lower order moment (FLOM) matrices, namely FLOM_MUSIC. However, it is limited in the range of $2\geq \alpha > 1$. Authors in [20] introduce a new subspace algorithm based on the phased fractional lower order moment (PFLOM), namely PFLOM_MUSIC, which it is applicable for $0 < \alpha\leq 2$. In [21], a subspace-augmented MUSIC technique for recovering the joint sparse support of a signal ensemble corrupted by additive impulsive noises is introduced. In order to mitigate the performance degradation of the DOA estimation methods based on the sparse representation of the second order statistics of the array output, new algorithms are proposed in this paper by using the sparse representation of the fractional lower order statistics vector of the array output. To enhance the robustness performance of the proposed algorithms, the improved algorithms are advanced by eliminating the fractional lower order statistics of the noise from the fractional lower order statistics vector of the array output through a linear transformation. Computer simulation experiments are presented to illustrate the performance superiority of the proposed methods over the DOA estimation method based on the sparse representation of the second order statistics of the array output under $\alpha$-stable noise environments.

Ⅱ. $\alpha$-STABLE DISTRIBUTION

The $\alpha$-stable distribution's probability density function does not have closed form. It can be conveniently described by its characteristic function as

$ \phi(t)=e^{\{jat-\gamma|t|^{\alpha}[1+j\beta {\rm sgn}(t)\varpi(t, \alpha)]\}} $ (1)

where $\varpi(t, \alpha)=\tan\frac{\pi\alpha}{2}$, if $\alpha \neq 1$; $\varpi(t, \alpha)=\frac{2}{\pi}\log|t|$, if $\alpha$ $=1$, and ${\rm sgn}(t)$ is $|t|$ if $t\neq 0$ and 0 if $t=0$. $\alpha$ is the characteristic exponent, it controls the thickness of the tail in the distribution and is restricted in $0 < \alpha\leq 2$. $\gamma$ is the dispersion parameter and is similar to the variance of the Gaussian distribution. $\beta$ is the symmetry parameter. If $\beta$ $=0$, the distribution is symmetric and the observation is referred to as the S$\alpha$S (symmetry $\alpha$-stable) distribution. $a$ is the location parameter. When $\alpha=2$ and $\beta=0$, the $\alpha$-stable distribution becomes a Gaussian distribution. The tails of stable distribution with characteristic exponent $0 < \alpha < 2$ are significantly thicker than that of the Gaussian distribution and the smaller $\alpha$, the thicker the tails. An important difference between the Gaussian and the $\alpha$-stable distributions ($0 < $ $\alpha$ $ < $ $2$) is that only moments of order less than $\alpha$ exist for the $\alpha$-stable distribution. Due to the non-existence of the second order statistics of the $\alpha$-stable distribution when the characteristic exponent is restricted in $0 < \alpha < 2$, the second order statistics, such as correlation and covariance, do not make sense. Therefore, the fractional lower order statistics (FLOS) have been defined [18], such as the fraction lower order moment (FLOM) in [19] and the phased fractional lower order moment (PFLOM) in [20].

Ⅲ. DOA ESTIMATION BASED ON SPARSE REPRESENTATION OF SECOND ORDER STATISTICS VECTOR

Consider the case of $K$ narrow far-field signals $s_1(t)$, $s_2(t)$, $...$, $s_K(t)$ with different DOAs $\theta_1$, $\theta_2$, $...$, $\theta_K$ arriving at a uniform linear array (ULA) with $M$ sensors in presence of additive noises $n_1(t)$, $n_2(t)$, $...$, $n_M(t)$. Assume that the noise is i.i.d random variable and is not correlated with signals. The received signal vector is given by

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

where

$ \begin{align*} &X_{M\times1}(t)=[x_1(t), x_2(t), ..., x_M(t)]^T\\ &A_{M\times K}(\theta)=[{a}(\theta_1), {a}(\theta_2), ..., {a}(\theta_K)]\\ &S_{K\times 1}(t)=[s_1(t), s_2(t), ..., s_K(t)]^T\\ &N_{M\times1}(t)=[n_1(t), n_2(t), ..., n_M(t)]^T \end{align*} $

where $x_m(t)$, $m=1, 2, ..., M$, is the output of the $m$th array element, ${a}(\theta_n)$, $n=1, 2, ..., K$, is the steering vector which can be expressed as

$ {a}(\theta_n)=\left[1, e^{-j\frac{2\pi}{\lambda}d\sin\theta_n}, ..., e^{-j\frac{2\pi}{\lambda}(M-1)d\sin\theta_n}\right]^T $ (3)

where $\lambda$ is the carrier wavelength of the signal, $d$ is the intersensor spacing.

Assume the noise in (2) is zero-mean Gaussian white noise with the power of $\sigma_n^2$, the second order statistics covariance matrix of the array output can be expressed as

$ R=E(X(t)X^H(t))=A(\theta)R_sA^H(\theta)+\sigma_n^2 I_M $ (4)

where the source covariance matrix $R_s=E(s(t)s^H(t))={\rm diag}\{\sigma_s\}$ is diagonal with source signal power vector $\sigma_s=$ $[\sigma_1^2, $ $..., \sigma_K^2]^T$ and $I_M$ denotes the $M\times M$ identity matrix. $\sigma_i^2$, $i=1, ..., K$, is the source signal power.

Applying the vectorization operator in (4), we have [22]

$ y=\mbox{vect}(R)=B(\theta)\sigma_s+\sigma_n^2\mbox{vect}(I_M) $ (5)
$ B(\theta)=[a^\ast(\theta_1)\otimes a(\theta_1), ..., a^\ast(\theta_K)\otimes a(\theta_K)] $ (6)

where $\otimes$ denotes Kronecker product. It is interesting to see that (5), similar to (2), can be taken as the array output of single snapshot where $B(\theta)$, $\sigma_s$ and $\mbox{vect}(I_M)$ are the virtual manifold matrix with dimension $M^2\times K$, equivalent source vector, and equivalent noise vector, respectively. The new signal vector $y$ can be sparsely represented in a redundant basis. Define a set $\hat{\theta}=[\hat{\theta}_1, \hat{\theta}_2, ..., \hat{\theta}_Q]$ which denotes potential source location of interest and assume that the true DOAs are exactly on this set. The number of the potential source locations $Q$ should be much greater than the number of actual sources $K$ and the number of virtual array sensors $M^2$. Define the overcomplete basis $B(\hat{\theta})=[a^\ast(\hat{\theta}_1)\otimes a(\hat{\theta}_1), ..., a^\ast(\hat{\theta}_Q)\otimes a(\hat{\theta}_Q)]$ and the signal power vector $\nu=[\nu_1, \nu_2, ..., \nu_Q]^T$ where $a^\ast(\hat{\theta}_i)\otimes a(\hat{\theta}_i)$ denotes the steering vector of the virtual array and the elements of vector $\nu$ have $K$ nonzeros, that is, $\nu_j$ $=$ $\sigma_i^2$ if $\hat{\theta}_j=\theta_i$, $i=1, ..., K$. As a result, $y$ can be rewritten as the following form:

$ y=B(\hat{\theta})\nu+\sigma_n^2\mbox{vect}(I_M). $ (7)

Hence the DOA estimation can be reduced to the detection of the nonzero elements of $\nu$. In practice, the unknown $y$ is estimated from the $N$ snapshots, let $\hat{y}$ be the estimation of $y$, then $\hat{y}=\mbox{vect}(\hat{R})$, where $\hat{R}=\frac{1}{N}\sum_{t=1}^NX(t)X^H(t)$. Define $\Delta y$ as the estimation error, then $\Delta y=\hat{y}-y$. Let $\hat{\nu}$ be the estimate of $\nu$, the DOA estimation problem can be further converted into the following convex optimization problem [12]:

$ \min\|\hat{\nu}\|_1, \quad {\rm s.t.}~\|\hat{y}-B(\hat{\theta})\hat{\nu}-\sigma_n^2\mbox{vect}(I_M)\|_2\leq \varepsilon $ (8)

where $\varepsilon$ is a parameter which means how much of the error we wish to allow and plays an important role in the algorithm performance. It can be known that the error $\Delta y$ satisfies asymptotically normal (AsN) distribution [23],

$ \Delta y=\hat{y}-y=\mbox{vect}(\hat{R}-R)\sim \mbox{AsN}\left(0_{M^2, 1}, \frac{1}{N}R^{T}\otimes R\right). $ (9)

Define the weighted matrix $W^{-\frac{1}{2}}=\sqrt{N}R^{-\frac{T}{2}}\otimes R^{-\frac{1}{2}}$, then

$ W^{-\frac{1}{2}}\Delta y\sim \mbox{AsN}\left(0_{M^2, 1}, I_{M^2}\right). $ (10)

Then from (7), we can further get

$ \left\|W^{-\frac{1}{2}}\left[\hat{y}-B(\hat{\theta})\hat{\nu} -\sigma_n^2\mbox{vect}(I_M)\right]\right\|_2^2 \sim \mbox{As}\chi^2(M^2) $ (11)

where $\mbox{As}\chi^2(M^2)$ is the asymptotic chi-square distribution with $M^2$ degrees of freedom. Therefore, the parameter $\varepsilon$ should be introduced such that

$ \left\|W^{-\frac{1}{2}}\left[\hat{y}-B(\hat{\theta})\hat{\nu} -\sigma_n^2\mbox{vect}(I_M)\right]\right\|_2^2\leq \varepsilon^2 $

with a high probability $pc$, that is $\varepsilon={\rm sqrt}({\chi^2_{pc}(M^2)})$. Let $\hat{W}^{-\frac{1}{2}}$ $= \sqrt{N}\hat{R}^{-\frac{T}{2}}\otimes \hat{R}^{-\frac{1}{2}}$ be the estimate of the weighted matrix $W^{-\frac{1}{2}}$ and $\hat{\sigma}_n^2$ be the estimate of $\sigma_n^2$ by the average of $M$ - $K$ smallest eigenvalues of the eigenvalue decomposition (EVD) of the estimate covariance matrix $\hat{R}$, then the statistically robust and tractable formula for DOA estimation can be reduced as follows:

$ \min\|\hat{\nu}\|_1, \quad {\rm s.t.}~ \left\|\hat{W}^{-\frac{1}{2}}\left[\hat{y}-B(\hat{\theta})\hat{\nu} -\hat{\sigma}_n^2\mbox{vect}(I_M)\right]\right\|_2\leq \varepsilon. $ (12)

This DOA estimation algorithm based on the sparse representation of the second order statistics covariance vector can be named as SS_SOSCV algorithm.

Ⅳ. DOA ESTIMATION BASED ON SPARSE REPRESENTATION OF FRACTIONAL LOWER ORDER STATISTICS VECTOR

When the noise in (2) is $\alpha$-stable impulsive noise with a characteristic exponent $0 < \alpha < 2$, the performance of the SS_SOSCV algorithm will degrade since the covariance matrix is not defined for $0 < \alpha < 2$. In this case, introducing a modified covariance matrix instead of the covariance matrix can alleviate the problem. In this paper, we introduce two DOA estimation methods based on sparse representation of fractional lower order statistics vector, and the improved algorithms which can enhance the robustness of the proposed algorithms are further studied.

A. SS_FLOMV Algorithm

The FLOM matrix $C$ which is suitable for $\alpha$-stable distribution noise environments can be used to replace the covariance matrix $R$ in (4). The $(i, k)$ element of matrix $C$ can be defined as

$ C_{i, k}={\rm E}\{x_i(t)|x_k(t)|^{p-2}x_k^\ast(t)\} $ (13)

where $p$ is the order of the moments. Setting $p=2$ reduces (13) to an appropriate covariance matrix under the condition of Gaussianity. However, as we deviate from this condition $p$ should be set to a lower value and it must satisfy the inequality 1 < $p < \alpha\leq2$ so that $C_{i, k}$ is bounded. It can be proved that the FLOM matrix $C$ can be expressed as [19]

$ C=A(\theta)\Lambda_sA^H(\theta)+\xi I_M $ (14)

where the diagonal matrix $\Lambda_s={\rm diag}\{\gamma_s\}$ can be interpreted as the FLOM matrix of the source signals and $\xi$ can be interpreted as the FLOM of the $\alpha$-stable additive noise level. $\gamma_s$ $=[\gamma_1, ..., \gamma_K]^T$, $\gamma_i$, $i=1, ..., K$, is the fractional lower order power of the signals. Applying the vectorization operator on (14), we have

$ y_{\rm FLOM}=\mbox{vect}(C)=B(\theta)\gamma_s+\xi\mbox{vect}(I_M). $ (15)

The vector $y_{\rm FLOM}$ can be sparsely represented in the overcomplete basis $B(\hat{\theta})$ as the following form:

$ y_{\rm FLOM}=B(\hat{\theta})\nu_{\rm FLOM}+\xi\mbox{vect}(I_M). $ (16)

As with the SS_SOSCV algorithm, the DOA estimation can be resolved by the following convex optimization problem:

$ \begin{align} &\min\|\hat{\nu}_{\rm FLOM}\|_1 \\ & {\rm s.t.}\quad \left\|\hat{W}^{-\frac{1}{2}}_{\rm FLOM}\left[ \hat{y}_{\rm FLOM}-B(\hat{\theta})\hat{\nu}_{\rm FLOM}-\hat{\xi}\mbox{vect}(I_M)\right]\right\|_2 \leq\varepsilon \end{align} $ (17)

where the weighted matrix $\hat{W}^{-\frac{1}{2}}_{\rm FLOM}$ can be defined as

$ \hat{W}^{-\frac{1}{2}}_{\rm FLOM}=\sqrt{N}\hat{C}^{-\frac{T}{2}}\otimes\hat{C}^{-\frac{1}{2}} $ (18)

$\hat{C}$ is the estimate of the FLOM matrix $C$ and the $(i, k)$ element of matrix $\hat{C}$ can be defined as

$ \hat{C}_{i, k}=\frac{1}{N}\sum\limits_{t=1}^N\left\{x_i(t)|x_k(t)|^{p-2}x_k^\ast(t)\right\} $ (19)

$\hat{\nu}_{\rm FLOM}$ and $\hat{y}_{\rm FLOM}$ are the estimation of $\nu_{\rm FLOM}$ and $y_{\rm FLOM}$, $\hat{\xi}$ is the estimation of $\xi$ by the average of $M-K$ smallest eigenvalues of the EVD of the matrix $\hat{C}$. This DOA estimation method based on the sparse representation of the FLOM vector can be named as SS_FLOMV algorithm.

B. SS_PFLOMV Algorithm

From the FLOM definition, it can be seen that it is limited in range of $2\geq \alpha>1$, so the SS_FLOMV algorithm is not applicable under the $\alpha$-stable noise with characteristic exponent $0 < \alpha\leq 1$. In [20] a new class of robust bounded covariance matrices based on phased fraction lower order moment (PFLOM) which is applicable for $0 < \alpha\leq 2$ are used. The $(i, k)$ element of PFLOM matrix $\Gamma$ can be defined as

$ \begin{equation} \Gamma_{i, k}={\rm E}\left\{x_i^{\langle b\rangle}(t)x_k^{\langle -b\rangle}(t)\right\}, \quad 0 <b <\frac{\alpha}{2} \end{equation} $ (20)

where the PFLOM operation on a complex number $z$ is

$ \begin{align} z^{\langle b\rangle}= \begin{cases} \frac{|z|^{b+1}}{z^\ast}, &z\neq 0\\ 0, &z=0 \end{cases} \end{align} $ (21)

and the conjugate of the $b$th PFLOM of $z$ as $z^{-\langle b\rangle}=(z^\ast)^{\langle b\rangle} =(z^{\langle b\rangle})^\ast$. It can be proved that the matrix $\Gamma$ can be expressed as [20]

$ \Gamma=A(\theta)\Phi_s A^H(\theta)+\kappa I_M $ (22)

where the diagonal matrix $\Phi_s={\rm diag}\{\varphi_s\}$ can be interpreted as the PFLOM matrix of the source signals and $\kappa$ can be interpreted as the PFLOM of the $\alpha$-stable additive noise level. $\varphi_s=[\varphi_1, ..., \varphi_K]^T$, $\varphi_i$, $i=1, ..., K$, is the phased fractional lower order power of the signals. Applying the vectorization operator on (22) and then sparse representation in the overcomplete basis $B(\hat{\theta})$, we have

$ y_{\rm PFLOM}=\mbox{vect}(\Gamma)=B(\theta)\varphi_s+\kappa \mbox{vect}(I_M) $ (23)
$ y_{\rm PFLOM}=B(\theta)\nu_{\rm PFLOM}+\kappa \mbox{vect}(I_M). $ (24)

Likewise, the DOAs can be estimated by solving the following optimization problem:

$ \begin{align} &\min\|\hat{\nu}_{\rm PFLOM}\|_1 \\ & {\rm s.t.}\ \Big\|\hat{W}_{\rm PFLOM}^{-\frac{1}{2}}\Big[ \hat{y}_{\rm PFLOM}-B(\hat{\theta})\hat{\nu}_{\rm PFLOM} \\ &\qquad-\hat{\kappa} \mbox{vect}(I_M)\Big]\Big\|_2 \leq \varepsilon \end{align} $ (25)

where the weighted matrix $\hat{W}_{\rm PFLOM}^{-\frac{1}{2}}$ can be defined as

$ \hat{W}_{\rm PFLOM}^{-\frac{1}{2}}=\sqrt{N}\hat{\Gamma}^{-\frac{T}{2}}\otimes \hat{\Gamma}^{-\frac{1}{2}} $ (26)

$\hat{\Gamma}$ is the estimation of the PFLOM matrix $\Gamma$ and the $(i, k)$ element of matrix $\hat{\Gamma}$ can be defined as

$ \hat{\Gamma}_{i, k}=\frac{1}{N}\sum\limits_{t=1}^N\left\{x_i^{\langle b\rangle}(t)x_k^{\langle -b\rangle}(t)\right\} $ (27)

$\hat{\nu}_{\rm PFLOM}$ and $\hat{y}_{\rm PFLOM}$ are the estimation of $\nu_{\rm PFLOM}$ and $y_{\rm PFLOM}$, $\hat{\kappa}$ is the estimation of $\kappa$ by the average of $M-K$ smallest eigenvalues of the EVD of the matrix $\hat{\Gamma}$. This DOA estimation method based on the sparse representation of the PFLOM vector can be named as SS_PFLOMV algorithm.

C. Improved Algorithms

Equations (16) and (24) can be unified and expressed as

$ y_{\rm FLOS}=B(\hat{\theta})\nu_{\rm FLOS}+(\xi|\kappa)\mbox{vect}(I_M). $ (28)

Notice that the vector $(\xi|\kappa)\mbox{vect}(I_M)$ has only $M$ nonzero elements, then these elements of $y_{\rm FLOS}$ corresponding to the positions of nonzero elements in $(\xi|\kappa)\mbox{vect}(I_M)$ can be removed and the rest $M(M-1)$ entries of $y_{\rm FLOS}$ corresponding to the positions of zeros elements in $(\xi|\kappa)\mbox{vect}(I_M)$ can be preserved. Mathematically, this operation can be formulated as

$ \begin{align} y_{\rm IFLOS}=&\ Jy_{\rm FLOS}\\ =&\ J\left\{B(\hat{\theta})\nu_{\rm FLOS}+(\xi|\kappa)\mbox{vect}(I_M)\right\}\\ =&\ D(\hat{\theta})\nu_{\rm FLOS} \end{align} $ (29)

where $J$ is an $M(M-1)\times M^2$ selecting matrix and can be represented as

$ J=[J_1, J_2, ..., J_{M-1}]^T $ (30)

where

$ \begin{align} J_m=&\ [\pmb{e}_{(m-1)(M+1)+2}, \pmb{e}_{(m-1)(M+1)+3}, \\ &\ ..., \pmb{e}_{m(M+1)}] \in\mathbb{R}^{M^2\times M}, \\ &\qquad\qquad\qquad\qquad m=1, ..., M-1 \end{align} $ (31)

$\pmb{e}_i$, $i=(m-1)(M+1)+2, ..., m(M+1)$, is an $M^2\times 1$ column vector with 1 at the $i$th position and 0 elsewhere. $D(\hat{\theta})$ $=$ $JB(\hat{\theta})\in C^{M(M-1)\times Q}$ is the new steering matrix. This elimination operation avoids the estimation of the fractional lower order statistics of the impulsive noise and further reduces the effect of the impulsive noise. Hence the DOA estimation can be obtained by the following minimization:

$ \min\|\nu_{\rm FLOS}\|_1, \quad {\rm s.t.}~\left\|y_{\rm IFLOS}-D(\hat{\theta})\nu_{\rm FLOS}\right\|_2 <\varepsilon_I. $ (32)

Letting $\hat{y}_{\rm IFLOS}$ be the estimation of $y_{\rm IFLOS}$ and $\Delta y_{\rm IFLOS}=\hat{y}_{\rm IFLOS}-y_{\rm IFLOS}$ be the estimation error, we can get $\Delta y_{\rm IFLOS}$ $=$ $ J\Delta y_{\rm FLOS}$. From (9), we know that $\Delta y_{\rm IFLOS}$ satisfies

$ \Delta y_{\rm IFLOS}\sim \mbox{AsN}\left(0_{M(M-1), 1}, \frac{1}{N} J\left((C|\Gamma)^T\otimes (C|\Gamma)\right)J^T\right). $ (33)

Define the weighted matrix ${W}_{\rm IFLOS}^{-\frac{1}{2}}$ as

$ {W}_{\rm IFLOS}^{-\frac{1}{2}}=\sqrt{N} J^{-\frac{1}{2}}\left((C|\Gamma)^{-\frac{T}{2}}\otimes (C|\Gamma)^{-\frac{1}{2}}\right)J^{-\frac{T}{2}} $ (34)

and its estimation as $\hat{W}_{\rm IFLOS}^{-\frac{1}{2}}$, then

$ W_{\rm IFLOS}^{-\frac{1}{2}}\Delta y_{\rm IFLOS}\sim \mbox{AsN}\left(0_{M(M-1), 1}, I_{M(M-1), 1}\right) $ (35)
$ \left\|W_{\rm IFLOS}^{-\frac{1}{2}}\left[\hat{y}_{\rm IFLOS} -D(\hat{\theta})\hat{\nu}_{\rm FLOS}\right]\right\|_2^2 \sim \mbox{As}\chi^2[M(M-1)]. $ (36)

Therefore, a parameter $\varepsilon_I$ should be selected such that

$ \left\|W_{\rm IFLOS}^{-\frac{1}{2}}\left[\hat{y}_{\rm IFLOS} -D(\hat{\theta})\hat{\nu}_{\rm FLOS}\right]\right\|_2^2\leq \varepsilon_I^2 $

with a high probability $pc$, that is

$ \varepsilon_I=\sqrt{\chi^2_{pc}(M(M-1))} . $

Likewise, the DOAs can be estimated by solving the following optimization problem:

$ \begin{align} &\min\|\hat{\nu}_{\rm FLOS}\|_1\\ & {\rm s.t.}\quad \Big\|\hat{W}_{\rm IFLOS}^{-\frac{1}{2}}\big[\hat{y}_{\rm IFLOS} -D(\hat{\theta})\hat{\nu}_{\rm FLOS}\big]\Big\|_2\leq \varepsilon_I. \end{align} $ (37)

The DOA estimation methods based on (37) by using the FLOM matrix $C$ and the PFLOM matrix $\Gamma$ can be named as SS_IFLOMV and SS_IPFLOMV, respectively.

D. Algorithm Computational Costs and Steps

The main computational costs of the SS_FLOMV and SS_PFLOMV algorithms include the calculation of the FLOM matrix $C$ or the PFLOM matrix $\Gamma$, the EVD of the matrix $C$ or $\Gamma$ to estimate the parameter $\xi$ and $\kappa$, and solving the optimization problem of (17) and (25), require $O(NM^2)$, $O(M^3)$, and $O(Q^3)$, respectively. As the SS_IFLOMV and SS_IPFLOMV algorithms do not need to estimate the parameter $\xi$ and $\kappa$, so their computational costs are slightly lower than those of the SS_FLOMV and SS_PFLOMV algorithms. But the computational costs of these four algorithms are higher than those of subspace-based FLOM_MUSIC and PFLOM_MUSIC algorithms, where the main complexity of these two algorithms is in calculating the array covariance matrix $R$ and its EVD.

From the above analysis, the SS_FLOMV and SS_ PFLOMV algorithms' steps can be summarized as following:

Step 1: Obtain the FLOM estimate matrix $\hat{C}$ or the PFLOM estimate matrix $\hat{\Gamma}$ using the array received data by (19) or (27). Then apply the vectorization operator on them to get the vector $\hat{y}_{\rm FLOM}$ or $\hat{y}_{\rm PFLOM}$.

Step 2: Get the estimation of the parameter $\hat{\xi}$ or $\hat{\kappa}$ by the average of $M-K$ smallest eigenvalues of the EVD of the matrix $\hat{C}$ or $\hat{\Gamma}$.

Step 3: Calculate the weighted matrix $\hat{W}_{\rm FLOM}^{-\frac{1}{2}}$ or $\hat{W}_{\rm PFLOM}^{-\frac{1}{2}}$ by (18) or (26).

Step 4: Solve the convex optimization problem of (17) or (25) to get the estimation of the vector $\hat{\nu}_{\rm FLOM}$ or $\hat{\nu}_{\rm PFLOM}$.

Step 5: Estimate the DOAs according the location of nonzero elements in the vector $\hat{\nu}_{\rm FLOM}$ or $\hat{\nu}_{\rm PFLOM}$.

The SS_IFLOMV and SS_IPFLOMV algorithms' steps are similar to those of SS_FLOMV and SS_PFLOMV algorithms except that Step 2 is applying the elimination operation (29) on the vector $\hat{y}_{\rm FLOS}$ to get the vector $\hat{y}_{\rm IFLOS}$.

Ⅴ. SIMULATION RESULTS

In this section, a series of numerical experiments under different conditions are conducted to compare the performances of the proposed SS_FLOMV, SS_PFLOMV, SS_IFLOMV and SS_IPFLOMV algorithms with those of the FLOM_MUSIC, PFLOM_MUSIC and SS_SOSCV methods. Throughout this section, the convex optimization problems of (12), (17), (25) and (37) are resolved by using the software package CVX [24], the probability $pc$ in the proposed algorithms is set as 0.999. An $M=8$ element ULA with an intersensor pacing of half a wavelength is used. The direction grid is set to have 181 points sampled from $-90^\circ$ to $90^\circ$ with $1^\circ$ intervals. Two performance criteria are used to assess the performances of the algorithms. The first one is the probability of resolution. The DOAs are considered to be resolved within $1^\circ$ estimate error. 2000 independent Monte Carlo experiments are performed. The experiment number that DOAs can be resolved is denoted by $N_{ok}$, then the probability of resolution is defined as $N_{ok}/2000$. In the case that DOAs can be resolved, set $\overline{\theta}_i(n)$, $i=1, 2, ..., K$, as the estimation of $\theta_i$ for the $n$th Monte Carlo experiment, the average mean square error (RMSE) of the DOAs estimation is defined as

$ { RMSE}=\frac{1}{N_{ok}}\sum\limits_{n=1}^{N_{ok}}\sqrt{\frac{1}{K} \sum\limits_{i=1}^K\left(\bar{\theta}_i(n)-\theta_i(n)\right)^2} . $ (38)

As the characteristic of the $\alpha$-stable distribution makes the use of the standard SNR meaningless, a new SNR measure, generalized signal-to-noise ratio (GSNR), is defined as [18]

$ {GSNR}=10\lg\frac{\delta_s^2}{\gamma} $ (39)

where $\delta_s$ is the variance of the signal, $\gamma$ is the dispersion parameter of the $\alpha$-stable noise.

Example 1: Three sources impinging on array for $-50^\circ$, $0^\circ$, and $50^\circ$ under the condition of $\alpha$ stable distribution noise with characteristic exponent $\alpha=1.5$ are considered. The GSNR is 10 dB and the number of snapshots is fixed at 100. Figs. 1-7 are the normalized spatial spectrums of FLOM_MUSIC, PFLOM_MUSIC, SS_SOSCV, SS_FLOMV, SS_PFLOMV, SS_IFLOMV, and SS_IPFLOMV algorithms, respectively. It can be seen that sparse representation based methods have higher resolution than the subspace based methods, that is, the normalized spatial spectrums in Figs. 3-7 are sharper than those in Figs. 1 and 2. In these sparse representation based methods, the SS_FLOMV and SS_PFLOMV ones proposed in this paper have a better spatial spectrum performance than the SS_SOSCV algorithm. And the improved SS_IFLOMV and SS_IPFLOMV algorithms have the best spatial spectrum performances.

Download:
Fig. 1 Normalized spatial spectrum of FLOM_MUSIC algorithm.
Download:
Fig. 2 Normalized spatial spectrum of PFLOM_MUSIC algorithm.
Download:
Fig. 3 Normalized spatial spectrum of SS_SOSCV algorithm.
Download:
Fig. 4 Normalized spatial spectrum of SS_FLOMV algorithm.
Download:
Fig. 5 Normalized spatial spectrum of SS_PFLOMV algorithm.
Download:
Fig. 6 Normalized spatial spectrum of SS_IFLOMV algorithm.
Download:
Fig. 7 Normalized spatial spectrum of SS_IPFLOMV algorithm.

Example 2: Three sources impinging on array for $-50^\circ$, $0^\circ$ and $50^\circ$ under the condition of $\alpha$ stable distribution noise with characteristic exponent $\alpha=1.5$ are considered, the number of snapshots is fixed at 100. Figs. 8 and 9 show the comparisons of the probability of resolution and the RMSE with the increase of GSNR between the proposed methods and the SS_SOSCV method, respectively. It can be seen that the probability of resolution and RMSE performance of all methods improve with the increased GSNR. However, the performances of the proposed methods which are based on the sparse representation of the fractional lower order statistics vector are much better than that of the second order statistics based methods. At the same time, the SS_IFLOMV and SS_IPFLOMV methods have better performances than the SS_FLOMV and SS_PFLOMV methods, since the effects of the noises on the algorithms are further reduced by the linear transform on the fractional lower order statistics vector of the array output. It can also be seen that the performances of the methods which are based on the sparse representation of the PFLOM vector are slightly better than those of the methods which are based on the sparse representation of the FLOM vector.

Download:
Fig. 8 Probability of resolution versus GSNR.
Download:
Fig. 9 RMSE of DOA estimation versus GSNR.

Example 3: Three sources impinging on array for $-50^\circ$, $0^\circ$ and $50^\circ$ under the condition of $\alpha$ stable distribution noise with characteristic exponent $\alpha=1.5$ are considered under the condition of $GSNR =$ 4 dB. Figs. 10 and 11 show the simulated performances of five algorithms versus the number of snapshots. It can be seen from Fig. 10 that the proposed SS_IFLOMV and SS_IPFLOMV algorithms have similar probability of resolution performances, which are better than those of the SS_PFLOMV and SS_FLOMV. The SS_SOSCV method has the worst probability of resolution performance compared with the other algorithms. It can be seen from Fig. 11 that the RMSEs of the proposed algorithms decrease monotonically with the number of snapshots, the proposed SS_IFLOMV and SS_IPFLOMV algorithms show a more satisfactory performance than the SS_FLOMV and SS_PFLOMV algorithms, especially when the snapshot is smaller than 400.

Download:
Fig. 10 Probability of resolution versus number of snapshots.
Download:
Fig. 11 RMSE versus number of snapshots.

Example 4: In this example, the performances of the proposed algorithms versus the characteristic exponent $\alpha$ of the noise are assessed. The other simulation conditions are similar to those in Example 1 except that the GSNR is set at 10 dB. Firstly, the situation that the characteristic exponent varies from $\alpha=1$ to $\alpha=2$ is considered. The probability of resolution and RMSE performances of the five methods are displayed in Figs. 12 and 13. It can be seen from these figures that the results are similar to those of the above mentioned examples. As expected, the resolution capability improves and the RMSE decreases with an increased characteristic exponent and the performance of the FLOS based methods outperforms the SOS based methods. The performances of the SS_IFLOMV and SS_IPFLOMV algorithms outperform the SS_FLOMV and SS_PFLOMV algorithms, and, at the same time, the performances of the PFLOM vector based methods outperform the FLOM vector based methods.

Download:
Fig. 12 Probability of resolution versus characteristic exponent $\alpha$.
Download:
Fig. 13 RMSE versus characteristic exponent $\alpha$.

Although the FLOM and PFLOM have the good performances in suppressing the $\alpha$-stable impulsive noise, they are applicable for different impulsive environment. The FLOM is limited in range of $2\geq \alpha>1$ and the PFLOM is applicable for $0 < \alpha\leq 2$. In other words, although FLOM can be calculated by the average in practice, there is no definition for FLOM in theory for $0 < \alpha\leq 1$. So, it can be predicted that the performance of the SS_FLOMV algorithm is inferior to that of the SS_PFLOMV algorithm. The SS_FLOMV algorithm will not even work, when the characteristic exponent of the impulsive noise is in the range of $0 < \alpha\leq 1$. To verify this, Figs. 14 and 15 show the simulated performances of the SS_FLOMV and SS_PFLOMV algorithm under the condition of $0.1\leq\alpha\leq 1$. It can be seen that the probability of resolution of the SS_FLOMV algorithm is zero, which means it does not work, when $0.1 < \alpha < 0.6$. So at this time the RMSE of the SS_FLOMV algorithm does not exist either. However, the SS_PFLOMV algorithm maintains a stable lower probability of resolution and has a fluctuating RMSE when $0.1 < \alpha < 0.6$. And the performance of SS_FLOMV algorithm is much inferior than that of SS_PFLOMV algorithm when $0.6\leq \alpha\leq 1$.

Download:
Fig. 14 Probability of resolution versus characteristic exponent $\alpha$ when $\alpha\leq 1$.}
Download:
Fig. 15 RMSE versus characteristic exponent $\alpha$ when $\alpha\leq 1$.
Ⅵ. CONCLUSION

The new methods based on sparse representation of the fractional lower order statistics vector are proposed for the DOAs estimation under $\alpha$-stable distribution impulsive noise environments. To enhance the performance of the proposed algorithms, the improved algorithms are advanced by a linear transformation on the fractional lower order statistics vector of the array output. Simulation results have shown the effectiveness of the proposed methods for a wide range of highly impulsive environments.

REFERENCES
[1] H. C. So, "Source localization: algorithms and analysis, " in Handbook of Position Location: Theory, Practice, and Advances, S. A. Zekavat and R. M. Buehrer, Eds. New York, NY, USA: Wiley-IEEE Press, 2011. http://onlinelibrary.wiley.com/doi/10.1002/9781118104750.ch2/pdf
[2] H. Krim and M. Viberg, "Two decades of array signal processing research: the parametric approach, " IEEE Signal Process. Mag., vol. 13, no. 4, pp. 67-94, Jul. 1996. http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=526899
[3] J. J. Fuchs, "On the use of the global matched filter for DOA estimation in the presence of correlated waveforms, " in Proc. of the 42nd Asilomar Conf. on Signals, Systems and Computers, Pacific Grove, CA, USA, 2008, pp. 269-273. http://ieeexplore.ieee.org/document/5074406/
[4] J. J. Fuchs, "Identification of real sinusoids in noise, the global matched filter approach, " in Proc. of the 15th IFCA Symp. System Identification, Saint-Malo, France, 2009, pp. 1127-1132. http://dx.doi.org/10.3182/20090706-3-FR-2004.00187
[5] D. Malioutov, M. Cetin, and A. S. Willsky, "A sparse signal reconstruction perspective for source localization with sensor arrays, " IEEE Trans. Signal Process., vol. 53, no. 8, pp. 3010-3022, Aug. 2005. https://academic.oup.com/bioinformatics
[6] X. Xu, X. Wei, and Z. Ye, "DOA estimation based on sparse signal recovery utilizing weighted l1-norm penalty, " IEEE Signal Process. Lett., vol. 19, no. 3, pp. 155-158, Mar. 2012. http://ieeexplore.ieee.org/document/6129390/
[7] N. Hu, Z. F. Ye, D. Y. Xu, and S. H. Cao, "A sparse recovery algorithm for DOA estimation using weighted subspace fitting, " Signal Process., vol. 92, no. 10, pp. 2566-2570, Oct. 2012. http://www.sciencedirect.com/science/article/pii/S0165168412001119
[8] X. H. Wei, Y. B. Yuan, and Q. Ling, "DOA estimation using a greedy block coordinate descent algorithm, " IEEE Trans. Signal Process., vol. 60, no. 12, pp. 6382-6394, Dec. 2012. http://dl.acm.org/citation.cfm?id=2710152.2710644&coll=DL&dl=GUIDE
[9] M. M. Hyder and K. Mahata, "Direction of arrival estimation using a mixed l2, 0 norm approximation, " IEEE Trans. Signal Process., vol. 58, no. 9, pp. 4646-4655. http://ieeexplore.ieee.org/document/5466152/
[10] P. Stoica, P. Babu, and J. Li, "SPICE: a sparse covariance-based estimation method for array processing, " IEEE Trans. Signal Process., vol. 59, no. 12, pp. 629-638, Feb. 2011. http://www.ams.org/mathscinet-getitem?mr=2790596
[11] J. H. Yin and T. Q. Chen, "Direction-of-arrival estimation using a sparse representation of array covariance vectors, " IEEE Trans. Signal Process., vol. 59, no. 92, pp. 4489-4493, Sep. 2011. http://dl.acm.org/citation.cfm?id=2333345
[12] Z. Q. He, Q. H. Liu, L. N. Jin, and S. Ouyang, "Low complexity method for DOA estimation using array covariance matrix sparse representation, " Electron. Lett., vol. 49, no. 3, pp. 228-230, Jan. 2013. http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=6457588
[13] Z. Q. He, Z. P. Shi, and L. Huang, "Covariance sparsity-aware DOA estimation for nonuniform noise, " Digital Signal Process., vol. 28, pp. 75-81, May 2014. http://www.sciencedirect.com/science/article/pii/S1051200414000505
[14] Y. Tian, X. Y. Sun, and S. S. Zhao, "DOA and power estimation using a sparse representation of second order statistics vector and l0-norm approximation, " Signal Process., vol. 105, pp. 98-108, Dec. 2014. http://www.sciencedirect.com/science/article/pii/S0165168414002321
[15] Y. Tian, X. Y. Sun, and S. S. Zhao, "Sparse-reconstruction-based direction of arrival, polarisation and power estimation using a cross-dipole array, " IET Radar, Sonar Navig., vol. 9, no. 6, pp. 727-731, Jun. 2015. http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=7127163
[16] Z. H. Sha, Z. M. Liu, Z. T. Huang, and Y. Y. Zhou, "Covariance-based direction-of-arrival estimation of wideband coherent chirp signals via sparse representation, " Sensors, vol. 13, no. 9, pp. 11490-11497, Aug. 2013. http://pubmedcentralcanada.ca/pmcc/articles/PMC3821291/
[17] A. Mahmood, M. Chitre, and M. A. Armand, "PSK communication with passband additive symmetric α-stable noise, " IEEE Trans. Commun., vol. 60, no. 10, pp. 2990-3000, Oct. 2012. http://ieeexplore.ieee.org/document/6253211/
[18] C. L. Nikias and M. Shoa, Signal Processing with Alpha-Stable Distributions and Applications. New York, NY, USA: Wiely, 1995.
[19] T. H. Liu and J. M. Mendel, "A subspace-based direction finding algorithm using fractional lower order statistics, " IEEE Trans. Signal Process., vol. 49, no. 8, pp. 1605-1613, Aug. 2001. http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=934131
[20] H. Belkacemi and S. Marcos, "Robust subspace-based algorithms for joint angle/Doppler estimation in non-Gaussian clutter, " Signal Process., vol. 87, no. 7, pp. 1547-1558, Jul. 2007. http://dl.acm.org/citation.cfm?id=1233044
[21] G. Tzagkarakis, P. Tsakalides, and J. L. Starck, "Covariation-based subspace-augmented MUSIC for joint sparse support recovery in impulsive environments, " Signal Process., vol. 93, no. 5, pp. 1365-1373, May 2013. http://www.sciencedirect.com/science/article/pii/S0165168412004161
[22] W. K. Ma, T. H. Hsieh, and C. Y. Chi, "DOA estimation of quasi-stationary signals with less sensors than sources and unknown spatial noise covariance: a Khatri-Rao subspace approach, " IEEE Trans. Signal Process., vol. 58, no. 4, pp. 2168-2180, Apr. 2010. http://ieeexplore.ieee.org/document/5290056/
[23] B. Ottersten, P. Stoica, and R. Roy, "Covariance matching estimation techniques for array signal processing applications, " Digit. Signal Process., vol. 8, no. 3, pp. 185-210. http://kth.diva-portal.org/smash/record.jsf?pid=diva2:501168
[24] M. Grant, S. Boyd, and Y. Ye, CVX: MATLAB Software for Disciplined Convex Programming. CVX Research, Inc., 2008.