程增飛 趙永波 水鵬朗 徐保慶

(西安電子科技大學(xué)雷達(dá)信號處理國家重點實驗室 西安 710071)

1 引言

作為陣列信號處理的重要分支，波達(dá)方向(Direction of Arrival,DOA)估計在雷達(dá)、聲吶及無線通信等領(lǐng)域有著重要應(yīng)用[1,2]。在入射信源為點信源的假設(shè)下，經(jīng)典的DOA估計算法(MUSIC,ESPRIT等)具有良好的DOA估計性能。但是，在移動無線通信、水下聲吶及無源雷達(dá)等應(yīng)用中，由于在信源及接收天線間不存在視距傳播，陣列天線接收到的信號為多條多徑散射信號的疊加，此時需要用分布式信源(Distributed Source,DS)模型對入射信源進(jìn)行描述，DS模型一般采用信號的中心入射角度和角度擴(kuò)展來描述信號的波達(dá)方向。由于信號模型的失配，經(jīng)典的DOA估計算法不能夠?qū)S的波達(dá)方向進(jìn)行有效估計[1]。

為了解決這一問題，20世紀(jì)90年代以來，通過對DS進(jìn)行建模，相關(guān)領(lǐng)域的學(xué)者從不同角度提出了多種有效的DS參數(shù)估計方法[2?11]。DS可以視為多個散射特性滿足一定統(tǒng)計規(guī)律的散射體的集合，根據(jù)散射體散射特性的不同，DS可以分為相干分布式信源、不相關(guān)分布式信源(Incoherently Distributed Source,IDS)和部分相關(guān)分布式信源，本文主要研究IDS的參數(shù)估計問題。文獻(xiàn)[3]和文獻(xiàn)[4]提出了一種DS的參數(shù)化描述方法，并據(jù)此給出了一種基于廣義MUSIC的DS參數(shù)估計算法(DISPARE)。在經(jīng)典最大似然算法的基礎(chǔ)上，文獻(xiàn)[2]給出了DS參數(shù)估計的最大似然算法，并針對最大似然算法的運算量問題，提出了計算復(fù)雜度較低的加權(quán)子空間擬合算法。在經(jīng)典的DOA估計算法中，Capon算法具有重要地位，根據(jù)IDS的信號模型特點，文獻(xiàn)[5]提出了基于廣義Capon譜估計器的IDS參數(shù)估計算法。和經(jīng)典的DOA估計算法相比，上述算法都具有良好的IDS參數(shù)估計性能，但是，由于其在進(jìn)行參數(shù)估計時需要進(jìn)行2維或高維搜索，算法的計算復(fù)雜度不能滿足實際系統(tǒng)的需要。為了降低算法運算量，文獻(xiàn)[6]提出了一種分布式信源參數(shù)的解耦估計算法。利用信號協(xié)方差矩陣的非結(jié)構(gòu)化模型，該算法首先對入射信源的角度擴(kuò)展進(jìn)行估計，然后基于估計出的角度擴(kuò)展值對入射信源的中心入射角度進(jìn)行計算。雖然該算法的運算量較2維搜素類算法的運算量低，但是其僅適用于單個入射信源的情況。針對上述算法的缺陷，文獻(xiàn)[7]中提出了基于多項式求根的低復(fù)雜度IDS參數(shù)估計算法，如擴(kuò)展求根MUSIC算法等。該類算法不僅計算復(fù)雜度大大降低，而且能夠有效用于多入射信源情況，但是其參數(shù)估計性能存在較大損失。

隨著稀疏表示理論的發(fā)展和完善，其在陣列信號處理領(lǐng)域中的應(yīng)用越來越廣泛[12?15]。本文結(jié)合IDS信源的模型特點，提出一種基于稀疏表示的IDS參數(shù)估計方法。該方法根據(jù)IDS采樣協(xié)方差矩陣的Toeplitz性質(zhì)，采用凸優(yōu)化方法對IDS的采樣協(xié)方差矩陣進(jìn)行擬合；然后，基于IDS協(xié)方差矩陣的兩點近似模型，采用稀疏表示的方法對IDS的角度擴(kuò)展參數(shù)進(jìn)行估計；利用估計出的角度擴(kuò)展參數(shù)，在IDS的角度分布形式已知的基礎(chǔ)上，采用IDS的Jacobi-Anger級數(shù)展開模型構(gòu)建稀疏基矩陣，并據(jù)此估計出IDS的中心入射角度。本文算法的運算量較搜索類算法低，在低信噪比和小快拍條件下具有良好的參數(shù)估計性能。

2 信號模型

考慮K個不相關(guān)分布式信源s(θ, θk, σθk),k=1,2,…,K，入射到由N個陣元構(gòu)成的等距線陣列(Uniform linear array,ULA)上，則接收信號可表示為

當(dāng)角自相關(guān)核函數(shù)p(θ, θk, σθk)為單峰對稱函數(shù)且IDS的角度擴(kuò)展較小時，記ωk=2π Δsin(θk)為對應(yīng)角度 θk方向上的空間頻率， σωk=2πΔcos(θk) σθk,Rk可近似表示為

其中λk表示第k個信源的能量，式(3)表明矩陣Rk的秩近似等于2，稱為IDS協(xié)方差矩陣的兩點近似模型[2]，其表明在分布式信源角度擴(kuò)展不是很大時，一個IDS的協(xié)方差矩陣可由兩個理想點信源的協(xié)方差矩陣近似表示。

式(3)是在小角度擴(kuò)展時對IDS協(xié)方差矩陣的近似表示，利用Jacobi-Anger級數(shù)可以得到IDS協(xié)方差矩陣的精確表達(dá)式[2]。記Rk(m,n)表示矩陣Rk的第m行第n列元素，可得

又根據(jù)Jocabi-Anger級數(shù)展開式，式(4)可展開為

其 中 ，z=2π(m?n)Δ， 當(dāng)g為 偶 數(shù) 時 αg=當(dāng)g為奇數(shù)時αg=表示g階第1類Bessel函數(shù)，式(5)稱為IDS協(xié)方差矩陣的Jacobi-Anger級數(shù)展開模型。

3 本文算法

由接收陣列為均勻線陣，容易證明矩陣Rk為Toeplitz矩陣，因此，由有限個Toeplitz矩陣的和仍為Toeplitz矩陣可知，矩陣Rs也是Toeplitz矩陣?；诰仃嘡s的Toeplitz性質(zhì)，為了降低噪聲對參數(shù)估計性能的影響，下面首先采用凸優(yōu)化方法對協(xié)方差矩陣R進(jìn)行擬合，然后采用稀疏表示方法對分布式信源參數(shù)進(jìn)行估計。

3.1 協(xié)方差矩陣擬合

在實際應(yīng)用中，協(xié)方差矩陣R一般由有限次快拍估計得到，即其中Q表示快拍個數(shù)，當(dāng)矩陣R和都可逆，即Q≥N時，考慮式(6)所示的協(xié)方差矩陣擬合準(zhǔn)則：

用rs表示矩陣Rs的第1列，則上述準(zhǔn)則可等價轉(zhuǎn)化為式(7)所示的凸優(yōu)化問題[16]：

其中Toep(·)表示Toeplitz操作，則由式(7)得到rs的估計值s后，可由s唯一確定出Rs的擬合值s。需要說明的是，由于優(yōu)化問題式(7)為一半定規(guī)劃問題[16]，因此由式(7)可以得到唯一的s?；跀M合出的信號協(xié)方差矩陣s及其Toeplitz性質(zhì)，下面首先采用稀疏表示的方法對IDS的角度擴(kuò)展進(jìn)行估計。

3.2 IDS角度擴(kuò)展估計

根據(jù)IDS協(xié)方差矩陣的兩點近似模型，由Rs=，則Rs可近似表示為

其中e1為擬合和稀疏表示誤差，u=[u1,… ,ul,… ,uL]T， 且 當(dāng) γl=ωk±σωk時 ，ul=λk/2， 否 則ul=0。求解式(10)，由向量中非零元素的位置可以得到ωk±σωk的估計值，根據(jù)ωk±σωk的估計值，可計算出分布式信源的角度擴(kuò)展值θk。需要說明的是，根據(jù)ωk±σωk的估計值也可以估計出IDS的中心入射角度，但是，仿真實驗發(fā)現(xiàn)根據(jù)這種方法對中心入射角度的估計精度很差。下面結(jié)合IDS協(xié)方差矩陣的Jocabi-Anger級數(shù)展開近似模型，在估計出的角度擴(kuò)展的基礎(chǔ)上，給出一種基于稀疏表示的IDS中心入射角度的估計算法。

3.3 IDS中心入射角度估計

雖然通過Jocabi-Anger級數(shù)展開可以得到IDS協(xié)方差矩陣的精確表達(dá)式，在實際應(yīng)用中，一般用有限階Jocabi-Anger級數(shù)展開式對IDS的協(xié)方差矩陣進(jìn)行近似。若記Jocabi-Anger級數(shù)展開階數(shù)為G，則Rk(m,n)可近似表示為

3.4 稀疏表示問題的求解

對于式(10)和式(14)中的稀疏表示問題，可以采用多種方法進(jìn)行求解，本文采用循環(huán)加權(quán)最小二乘(Iteratively Reweighted Least Square,IRLS)算法對其進(jìn)行求解[17,18]。由于式(10)和式(14)的形式完全一樣，為方便起見，用式(15)對其進(jìn)行統(tǒng)一表示。

其中0≤p＜1為范數(shù)約束值，ε表示允許的稀疏表示誤差。由于0≤p＜1，上述優(yōu)化問題的目標(biāo)函數(shù)為非凸的，導(dǎo)致式(16)不存在閉合形式的最優(yōu)解。IRLS算法通過在每一次迭代過程中求解一個加權(quán)形式的最小二乘問題來對上述非凸問題進(jìn)行求解，由于加權(quán)最小二乘問題具有閉合形式的最優(yōu)解，因此，該算法在每次循環(huán)過程中都能夠得到最優(yōu)形式的閉合解，運算量較小。具體而言，若設(shè)IRLS算法在第i次循環(huán)中得到的最優(yōu)解矢量為ti，則在第i+1次循環(huán)過程中IRLS算法需求解式(17)所示的無約束優(yōu)化問題：

其中μi+1是第i+1次循環(huán)中的正則化參數(shù)，γi+1=為第i+1次循環(huán)中的加權(quán)矢量，⊙表示Hadamard積。由此，得到基于稀疏表示的IDS參數(shù)估計算法步驟如下所述：

步驟 3 設(shè)定正則化參數(shù) λ0=λspread，初始解向量范數(shù)約束值p=0.8，采用IRLS算法對式(10)進(jìn)行求解，得到向量u的估計值，并由其峰值位置確定出角度擴(kuò)展的估計值θk,k=1,2, …,K；

步驟 5 設(shè)定正則化參數(shù) λ0=λcentral，初始解向量范數(shù)約束值p=0.8，采用IRLS算法對式(14)進(jìn)行求解，得到向量v的估計值，并由其峰值位置確定出中心角度k,k=1,2,…,K。

在算法運算量方面，本文算法的運算量主要集中在采用式(7)對協(xié)方差矩陣的擬合及采用RILS算法對式(10)和式(14)的求解。對于式(7)中的凸優(yōu)化問題，采用內(nèi)點法對其進(jìn)行求解時，其運算量為O N6.5()，而采用IRLS算法對(10)和式(14)進(jìn)行求解時，若記算法循環(huán)迭代總次數(shù)為I，其運算量為，因此本文算法的運算量為。而對于DISPARE和廣義Capon算法，若在進(jìn)行2維搜索時，其在角度擴(kuò)展及中心入射角度維劃分的網(wǎng)格點數(shù)分別為H和L，則兩算法的計算復(fù)雜度分別為O(HLKN2+N3)和O(HLN3)。由于一般情況下網(wǎng)格點數(shù)H和L遠(yuǎn)遠(yuǎn)大于陣元數(shù)N、信源數(shù)K及迭代次數(shù)I，因此，本文算法的計算量較DISPARE和廣義Capon算法小。

4 計算機(jī)仿真

本節(jié)對上述基于稀疏表示的IDS參數(shù)估計算法進(jìn)行分析，以驗證算法的有效性。為了對比分析本文算法的性能，同時對DISPARE算法、廣義Capon算法及擴(kuò)展求根MUSIC算法進(jìn)行了計算機(jī)仿真。實驗中，設(shè)定陣列天線為由N=10個陣元構(gòu)成的ULA，且Δ=0.5,IDS的角自相關(guān)核函數(shù)為高斯型函數(shù)，Jocabi-Anger級數(shù)展開階數(shù)G=30，分別對各算法性能隨接收信號信噪比、快拍數(shù)及IDS角度擴(kuò)展的變化情況進(jìn)行分析。

仿真1 矩陣擬合對采樣協(xié)方差矩陣的影響

為降低噪聲對參數(shù)估計性能的影響，本文方法在進(jìn)行參數(shù)估計前，首先對采樣協(xié)方差矩陣進(jìn)行了擬合處理，本實驗就擬合處理對采樣協(xié)方差矩陣的影響進(jìn)行仿真分析。仿真中通過計算原始采樣協(xié)方差矩陣及擬合得到的信號協(xié)方差矩陣同真實信號協(xié)方差矩陣Rs間差值的Frobemius范數(shù)來說明擬合的必要性。設(shè)定采樣快拍數(shù)為100，單個IDS的中心入射角度θ=0°，角度擴(kuò)展 σθ=4°，接收信號信噪比由-5 dB變化到25 dB，在每一接收信噪比獨立進(jìn)行300次實驗，得到及同Rs間差值的Frobemius范數(shù)隨信噪比的變化如圖1所示。

圖1采樣協(xié)方差矩陣及擬合得到的信號協(xié)方差矩陣同真實信號協(xié)方差矩陣間的差值隨信噪比的變化

分析圖1可知，在不同信噪比下，s同Rs間差值的Frobemius范數(shù)一直小于同Rs間差值的Frobemius范數(shù)，這說明擬合過程降低了噪聲對信號協(xié)方差矩陣的影響，能夠提高參數(shù)估計精度。

仿真2 算法性能隨接收信號信噪比的變化情況

本實驗中設(shè)定采樣快拍數(shù)為100，單個IDS的中心入射角度θ=0°，角度擴(kuò)展 σθ=4°，接收信號信噪比由-5 dB變化到25 dB，正則化參數(shù)初始值λspread和λcentral分別為0.1和0.2，在每一接收信號信噪比獨立進(jìn)行300次實驗，得到IDS的中心入射角度及角度擴(kuò)展的測量均方根誤差隨接收信號信噪比的變化情況如圖2所示。

圖2 各算法性能隨接收信號信噪比的變化曲線

由圖2可知，隨信噪比的增大，各算法的參數(shù)估計性能逐漸提高，且本文算法具有較好的性能。在信噪比較低時，求根類算法不能夠?qū)DS的參數(shù)進(jìn)行有效測量，DISPARE算法和廣義Capon算法具有較擴(kuò)展求根MUSIC算法更好的性能，但略遜于本文算法。

仿真3 算法性能隨快拍數(shù)的變化情況

本實驗中設(shè)定接收信號信噪比為10 dB，快拍數(shù)由10變化到500，其它參數(shù)同仿真1中相同，在每一快拍數(shù)下進(jìn)行300次獨立實驗，得到IDS的中心入射角度及角度擴(kuò)展的測量均方根誤差隨快拍數(shù)的變化情況分別如圖3(a)和圖3(b)所示。

圖3 各算法性能隨快拍數(shù)的變化曲線

分析圖3可知，雖然各算法的性能隨快拍數(shù)的增大逐漸提高，但是在快拍數(shù)大于200后，隨著快拍數(shù)的增加，算法性能受快拍數(shù)的影響不再明顯，同時，本文算法對中心角度的估計性能明顯優(yōu)于現(xiàn)有算法。

仿真4 算法性能隨角度擴(kuò)展的變化情況

本實驗中設(shè)定接收信號信噪比為10 dB，快拍數(shù)由100，單個IDS的中心入射角度θ=0°，角度擴(kuò)展σθ由1°變化到10°，其它參數(shù)同仿真1中相同，在每一角度擴(kuò)展值下進(jìn)行300次獨立實驗，得到IDS的中心入射角度及角度擴(kuò)展的測量均方根誤差隨分布式信源角度擴(kuò)展的變化情況分別如圖4(a)和圖4(b)所示。

圖4 各算法性能隨角度擴(kuò)展的變化曲線

由圖4可知，擴(kuò)展求根MUSIC及本文算法對中心入射角度的估計性能隨角度擴(kuò)展的增加并不是單調(diào)變化的，其原因在于擴(kuò)展求根MUSIC及本文算法對中心角度的測量是基于IDS協(xié)方差矩陣的兩點近似模型的，當(dāng)入射信號的角度擴(kuò)展很小時，算法不能夠有效測量出兩個近似點信源的入射角度值，導(dǎo)致算法性能下降。

5 結(jié)論

利用不相關(guān)分布式信源協(xié)方差矩陣的性質(zhì)，本文給出了一種基于稀疏表示的不相關(guān)分布式信源參數(shù)估計方法。該方法將IDS的2維參數(shù)估計問題轉(zhuǎn)化為了兩個1維稀疏表示問題，因此該算法的計算復(fù)雜度較2維搜索類算法低；同時，所提算法在進(jìn)行參數(shù)估計前首先根據(jù)IDS協(xié)方差矩陣的特點，對采樣協(xié)方差矩陣進(jìn)行了擬合，這能夠有效降低噪聲對算法性能的影響，提高算法的參數(shù)估計性能。

[1] Jantti T P.The influence of extended sources on the theoretical performance of the MUSIC and ESPRIT methods:Narrow-band sources[C]. Proceedings of International Conference on Acoustics,Speech and Signal Processing 1992,San Francisco,1992:429-432.

[2] Trump T and Ottersten B.Estimation of nominal direction of arrival and angular spread using an array of sensors[J].Signal Processing,1996,50(2):57-69

[3] Valaee S,Champagne B,and Kabal P.Parametric location of distribute source[J].IEEE Transactions on Signal Processing,1995,43(9):2144-2153.

[4] Meng Y,Stoica P,and Wong K M.Estimation of the direction of arrival of spatially dispersed signals in array processing[J].IEE Proceedings-Radar,Sonar,and Navigation,1996,143(1):1-9.

[5] Hassanien A,Shahbazpanshi S,and Gershman A B.A generalized Capon estimator for localization of multiple spread sources[J].IEEE Transactions on Signal Processing,2004,52(1):280-283.

[6] Besson O and Stoica P.Decoupled estimation of DOA and angular spread for a spatially distributed source[J].IEEE Transactions on Signal Processing,2000,48(8):2185-2194.

[7] Bengtsson M and Ottersten B.Low-complexity estimator for distributed sources[J].IEEETransactionsonSignal Processing,2000,48(7):1872-1882.

[8] 韓英華,汪晉寬,宋昕.相干分布式信源二維波達(dá)方向估計算法[J].電子與信息學(xué)報,2009,31(2):323-326.Han Ying-hua,Wang Jin-kuan,and Song Xin.2D DOA estimation algorithm for coherently distributed source[J].Journal of Electronics&Information Technology,2009,31(2):323-326.

[9] 楊學(xué)敏,李廣軍,鄭植.基于稀疏表示的相干分布式非圓信號的參數(shù)估計[J].電子與信息學(xué)報,2014,36(1):164-168.Yang Xue-min,Li Guang-jun,and Zheng Zhi.Parameter estimation of coherently distributed non-circular signal based on sparse representation[J].Journal of Electronics&Information Technology,2014,36(1):164-168.

[10] Han K Y and Nehorai A.Distributed source processing with nested linear array[C].IEEE 8th Sensor Array and Multichannel Signal Processing Workshop,Coruna,2014:521-524.

[11] Yang X M,Li G J,Zhang Z,et al..Low-complexity 2D central angle estimation of coherently distributed sources with cross-correlation matrix[J].Electronics Letters,2014,50(16):1118-1120.

[12] Candes E J and Wakin M B.An introduction to compressive sampling[J].IEEE Signal Processing Maganize,2008,25:21-30.

[13] Zheng J and Kaveh M.Sparse spatial spectral estimation:a covariance fitting algorithm,performance and regulation[J].IEEETransactions onSignal Processing,2013,61(11):2767-2777.

[14] Wang J,Sheng W X,Han Y B,et al..Adaptive beamforming with compressed sensoring for sparse receiving array[J].IEEE Transactions on Aerospace and Electronic Systems,2014,50(2):823-833.

[15] Zhao T,Eldar Y C,and Nehorai A.Direction of arrival estimation using co-prime arrays:a super resolution viewpoint[J].IEEE Transactions on Signal Processing,2014,62(21):5565-5576.

[16] Yang Z,Xie L H,and Zhang C S.A discretization-free sparse and parametric approach for linear array processing[J].IEEE Transactions on Signal Processing,2014,62(19):4959-4973.

[17] David W and Srikantan N.Iterative reweighted ?1and ?2methods for finding sparse solutions[J].IEEE Journal of Selected Topics in Signal Processing,2010,4(2):317-329.

[18] Ba D,Babadi B,Purdon P L,et al..Convergence and stability of iteratively re-weighted least squares algorithm[J].IEEETransactions onSignal Processing,2014,62(1):183-195.

[19] Chartrand R.Exact reconstruction of sparse signals via nonconvex minimization[J].IEEE Signal Processing Letters,2007,14(10):707-710.