DOA Estimation Based on Sparse Representation of the Fractional Lower Order Statistics in Impulsive Noise

Sen Li,Rongxi He,Member,IEEE,Bin Lin,Member,IEEE,and Fei Sun

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 αstable 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.

I.INTRODUCTION

D IRECTION 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(l1-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 l1-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 l1minimization 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 l2,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 l2,0-norm based joint sparse approximation technique is introduced into DOA estimation where the l0norm 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 l1-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 l1-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 l1-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 l1function to approximate l0-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 α-stable distribution(0 ＜ α ≤ 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 α-stable distribution is that only moments of order less than α 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 α-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≥α＞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＜α≤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 α-stable noise environments.

II.α-STABLE DISTRIBUTION

The α-stable distribution’s probability density function does not have closed form.It can be conveniently described by its characteristic function as

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

Consider the case of K narrow far- field signals s1(t),s2(t),...,sK(t)with different DOAs θ1,θ2,...,θKarriving at a uniform linear array(ULA)with M sensors in presence of additive noises n1(t),n2(t),...,nM(t).Assume that the noise is i.i.d random variable and is not correlated with signals.The received signal vector is given by

where

where xm(t),m=1,2,...,M,is the output of the m th array element,a(θn),n=1,2,...,K,is the steering vector which can be expressed as

where λ 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,the second order statistics covariance matrix of the array output can be expressed as

where the source covariance matrix Rs=E(s(t)sH(t))=diag{σs}is diagonal with source signal power vector σs=[...,σ2K]Tand IMdenotes the M × M identity matrix.,i=1,...,K,is the source signal power.

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

where⊗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(θ),σsand vect(IM)are the virtual manifold matrix with dimension M2×K,equivalent source vector,and equivalent noise vector,respectively.The new signal vector y can be sparsely represented in a redundant basis.Define a setwhich 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 M2.Define the overcomplete basisand the signal power vector ν=[ν1,ν2,...,νQ]Twheredenotes the steering vector of the virtual array and the elements of vector ν have K nonzeros,that is,νjAs a result,y can be rewritten as the following form:

Hence the DOA estimation can be reduced to the detection of the nonzero elements of ν.In practice,the unknown y is estimated from the N snapshots,letbe the estimation of y,thenwhereDefine Δy as the estimation error,thenLet ˆν be the estimate of ν,the DOA estimation problem can be further converted into the following convex optimization problem[12]:

where ε 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 Δy satisfies asymptotically normal(AsN)distribution[23],

Then from(7),we can further get

where Asχ2(M2)is the asymptotic chi-square distribution with M2degrees of freedom.Therefore,the parameter ε should be introduced such that

with a high probability pc,that isbe the estimate of the weighted

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

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

When the noise in(2)is α-stable impulsive noise with a characteristic exponent 0＜α＜2,the performance of the SS SOSCV algorithm will degrade since the covariance matrix is not defined for 0＜α＜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 α-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

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＜α≤2 so that Ci,kis bounded.It can be proved that the FLOM matrix C can be expressed as[19]

where the diagonal matrix Λs=diag{γs}can be interpreted as the FLOM matrix of the source signals and ξ can be interpreted as the FLOM of the α-stable additive noise level.γs=[γ1,...,γK]T,γi,i=1,...,K,is the fractional lower order power of the signals.Applying the vectorization operator on(14),we have

The vector yFLOMcan be sparsely represented in the overcomplete basis B(ˆθ)as the following form:

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

B.SS PFLOMV Algorithm

From the FLOM definition,it can be seen that it is limited in range of 2≥α＞1,so the SS FLOMV algorithm is not applicable under the α-stable noise with characteristic exponent 0＜α≤1.In[20]a new class of robust bounded covariance matrices based on phased fraction lower order moment(PFLOM)which is applicable for 0＜α≤2 are used.The(i,k)element of PFLOM matrix Γ can be defined as

where the PFLOM operation on a complex number z is

and the conjugate of the b th PFLOM of z as z-〈b〉=(z∗)〈b〉=(z〈b〉)∗.It can be proved that the matrix Γ can be expressed as[20]

where the diagonal matrix Φs=diag{ϕs}can be interpreted as the PFLOM matrix of the source signals and κ can be interpreted as the PFLOM of the α-stable additive noise level.ϕs=[ϕ1,...,ϕK]T,ϕ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(ˆθ),we have

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

C.Improved Algorithms

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

Notice that the vector(ξ|κ)vect(IM)has only M nonzero elements,then these elements of yFLOScorresponding to the positions of nonzero elements in(ξ|κ)vect(IM)can be removed and the rest M(M-1)entries of yFLOScorresponding to the positions of zeros elements in(ξ|κ)vect(IM)can be preserved.Mathematically,this operation can be formulated as

where J is an M(M-1)×M2selecting matrix and can be represented as

where

Therefore,a parameter εIshould be selected such that

with a high probability pc,that is

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

The DOA estimation methods based on(37)by using the FLOM matrix C and the PFLOM matrix Γ 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 Γ,the EVD of the matrix C or Γ to estimate the parameter ξ and κ,and solving the optimization problem of(17)and(25),require O(N M2),O(M3),and O(Q3),respectively.As the SS_IFLOMV and SS_IPFLOMV algorithms do not need to estimate the parameter ξ and κ,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ˆ as following:

Step 2:Get the estimation of the parameter orby the averageof M-K smallest eigenvalues of the EVD of the matrix o

Step 3:Calculate the weighted matrixby(18)or(26).

Step 4:Solve the convex optimization problem of(17)or(25)to get the estimation of the vectoror

Step 5:Estimate the DOAs according the location of nonzero elements in the vectoror

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 vectorto get the vector

V.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°to 90°with 1°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°estimate error.2000 independent Monte Carlo experiments are performed.The experiment number that DOAs can be resolved is denoted by Nok,then the probability of resolution is defined as Nok/2000.In the case that DOAs can be resolved,setθi(n),i=1,2,...,K,as the estimation of θifor the n th Monte Carlo experiment,the average mean square error(RMSE)of the DOAs estimation is defined as

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

where δsis the variance of the signal, γ is the dispersion parameter of the α-stable noise.

Example 1:Three sources impinging on array for-50°,0°,and 50°under the condition of α stable distribution noise with characteristic exponent α=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.

Fig.1.Normalized spatial spectrum of FLOM_MUSIC algorithm.

Fig.2.Normalized spatial spectrum of PFLOM_MUSIC algorithm.

Example 2:Three sources impinging on array for-50°,0°and 50°under the condition of α stable distribution noise with characteristic exponent α =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.

Fig.3.Normalized spatial spectrum of SS_SOSCV algorithm.

Fig.4.Normalized spatial spectrum of SS_FLOMV algorithm.

Fig.5.Normalized spatial spectrum of SS_PFLOMV algorithm.

Fig.6.Normalized spatial spectrum of SS_IFLOMV algorithm.

Fig.7.Normalized spatial spectrum of SS_IPFLOMV algorithm.

Fig.8.Probability of resolution versus GSNR.

Example 3:Three sources impinging on array for-50°,0°and 50°under the condition of α stable distribution noise with characteristic exponent α=1.5 are considered under the condition of GS N R=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.

Fig.9.RMSE of DOA estimation versus GSNR.

Fig.10.Probability of resolution versus number of snapshots.

Fig.11.RMSE versus number of snapshots.

Example 4:In this example,the performances of the proposed algorithms versus the characteristic exponent α 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 α =1 to α =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.

Fig.12.Probability of resolution versus characteristic exponent α.

Fig.13.RMSE versus characteristic exponent α.

Although the FLOM and PFLOM have the good performances in suppressing the α-stable impulsive noise,they are applicable for different impulsive environment.The FLOM is limited in range of 2≥α＞1 and the PFLOM is applicable for 0＜α≤2.In other words,although FLOM can be calculated by the average in practice,there is no definition for FLOM in theory for 0＜α≤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＜α≤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≤α≤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＜α＜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＜α＜0.6.And the performance of SS_FLOMV algorithm is much inferior than that of SS_PFLOMV algorithm when 0.6≤α≤1.

Fig.14.Probability of resolution versus characteristic exponent α when α ≤1.

Fig.15.RMSE versus characteristic exponent α when α ≤ 1.

VI.CONCLUSION

The new methods based on sparse representation of the fractional lower order statistics vector are proposed for the DOAs estimation under α-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.