一種不依賴超參數(shù)的稀疏信號單快拍DOA估計方法

虞 飛, 余 赟, 周利輝, 彭春光

(中國人民解放軍92578部隊, 北京 100071)

0 引 言

信號波達方向(direction of arrival,DOA)估計是陣列信號處理中一項極為重要的研究分支,已在聲納、雷達、通信、導(dǎo)航、地震探測等領(lǐng)域取得廣泛應(yīng)用[1]。以MUSIC[2-3]和ESPRIT[4]算法為代表的傳統(tǒng)子空間類DOA估計算法在高信噪比(signal to noise ratio, SNR)、非強相關(guān)源和多采樣快拍數(shù)的情況下可獲得超分辨估計性能。但在實際應(yīng)用中,高SNR、非強相關(guān)源和多采樣快拍數(shù)都是比較理想的環(huán)境因素,只要有一個因素不滿足,這些算法的DOA估計精度將嚴重下降甚至估計失敗。為此,學(xué)者們提出了很多改進算法,如適用于單快拍的陣列測向方法[5-7],適用于相干情形的空間平滑類算法[8-10]等。然而,這些改進算法都只適用于某一個非理想環(huán)境因素下的高性能估計,當同時存在至少兩個非理想環(huán)境因素時,這類算法依然面臨著失效問題。

近年來,基于稀疏表示框架的陣列參數(shù)估計方法由于在低SNR、有限采樣快拍數(shù)據(jù)、相干信號、信號DOA角度間隔小等非理想條件下同時具有良好的估計性能而引起了相關(guān)學(xué)者的廣泛關(guān)注,并取得了大量高質(zhì)量的研究成果[11-15]。考慮到在實際應(yīng)用中,通常目標信號個數(shù)遠小于傳感器陣元數(shù),目標信號DOA相對于空間來說也是稀疏的,故可將傳統(tǒng)的傳感器陣列輸出模型進行稀疏化表示,得到陣列輸出數(shù)據(jù)的稀疏表示模型。于是,稀疏表示理論中的很多方法很自然地被引入到陣列信號處理問題中,而傳統(tǒng)的陣列參數(shù)估計問題轉(zhuǎn)化為稀疏參數(shù)求解問題?；谙∈璞硎镜年嚵袦y向方法在非理想條件下整體性能優(yōu)于傳統(tǒng)的子空間類測向方法,但分析表明,前者的計算復(fù)雜度明顯高于后者,導(dǎo)致算法對目標信號DOA估計的實時性變得很差。同時,這類方法在求最優(yōu)解過程中,還面臨著一個甚至多個超參數(shù)(或者正則化參數(shù))選取的問題,超參數(shù)的選擇是否合適直接關(guān)系到最終的稀疏恢復(fù)性能。因此,還需要對超參數(shù)進行精細選取,然而目前還沒有明確的方法用來指導(dǎo)超參數(shù)的選擇。

本文依托陣列輸出數(shù)據(jù)的稀疏表示模型,提出了一種基于加權(quán)最小二乘(weighted least squares,WLS)準則的單快拍DOA稀疏迭代自適應(yīng)估計(簡稱為WLS-IAE)算法。該算法不僅保持了稀疏表示框架在非理想條件下的良好估計性能,而且避免了傳統(tǒng)稀疏參數(shù)求解方法的高復(fù)雜度和經(jīng)驗式的超參數(shù)選取問題,僅采用單快拍即獲得高精度DOA估計,因此非常適用于快變目標信號DOA的實時跟蹤測量,具有潛在的工程實用價值。

1 陣列輸出模型

考慮K個遠場窄帶信號入射到由M個無方向性(全向)陣元構(gòu)成的均勻線性陣列(uniform linear array,ULA)中,并假設(shè)K個信號與ULA在同一平面內(nèi)(在實際應(yīng)用場景中,如水下聲納基陣測向,我們經(jīng)常只關(guān)心某一平面內(nèi)信號的入射方位,或者入射信號在該平面內(nèi)的投影,因此這一假設(shè)可以得到保證)[16]。將陣元由1到M進行編號,并以陣元1作為基準或參考陣元。設(shè)參考陣元處的任一接收信號變換到基帶后有如下形式:

s(t)=a(t)ej[ω0t+φ(t)]

(1)

那么,在t時刻,整個陣列的M×1維輸出數(shù)據(jù)模型為

x(t)=As(t)+n(t)

(2)

式中,s(t)=[s1(t),s2(t),…,sK(t)]T為變換到基帶后的參考陣元處接收到的K個信號構(gòu)成的列向量;n(t)=[n1(t),n2(t),…,nM(t)]T為陣列的零均值加性復(fù)高斯白噪聲向量;A=[a(θ1),a(θ2),…,a(θK)]為M×K維導(dǎo)向矢量矩陣,且對于ULA,導(dǎo)向矢量a(θk)[17]可定義為

(3)

式中,d為相鄰陣元間距;λ為信號波長;θk為第k個信號的DOA,通常定義為該信號入射方向與ULA法線方向的夾角,則有θk∈[-π/2,π/2]。

2 模型的稀疏表示與稀疏矩陣方程的求解

x(t)=Φγ(t)+n(t)

(4)

min‖γ(t)‖0s.t.‖x(t)-Φγ(t)‖2≤β

(5)

式中,l0-范數(shù)‖γ(t)‖0表示向量γ(t)中非零元素的個數(shù)；‖·‖2表示向量的2-范數(shù)；正則化參數(shù)β用來指定允許的噪聲水平。

直接求解優(yōu)化問題式(5),必須篩選出向量γ(t)中所有可能的非零元素,由于搜索空間過于龐大,故此方法是非確定性多項式時間(non-deterministic polynomial-time, NP)難題。

考慮到l1-范數(shù)是最接近于l0-范數(shù)的凸目標函數(shù),目前,使用最廣泛的求解方法是將式(5)的l0-范數(shù)最小化問題轉(zhuǎn)化為凸松弛的l1-范數(shù)最小化問題(簡稱l1-min算法)[18],即

min‖γ(t)‖1s.t.‖x(t)-Φγ(t)‖2≤β

(6)

式中,‖·‖1表示向量的l1-范數(shù)。該式可以很容易轉(zhuǎn)化為二階錐規(guī)劃(second-order cone programming,SOCP)問題的一般形式,從而在二階錐規(guī)劃的框架下求解。但在實際應(yīng)用時,目標可能來向的角度范圍θ在網(wǎng)格劃分時其數(shù)量級將達到104,這將明顯降低二階錐規(guī)劃問題的求解速度,導(dǎo)致算法對DOA估計的實時性變得很差。

3 基于WLS的稀疏迭代自適應(yīng)估計方法

假設(shè)稀疏信號矢量γ(t)與噪聲矢量n(t)之間相互統(tǒng)計獨立,則由式(4)可得

ΦPΦH+Rn

(7)

P=diag[p1,p2,…,pn,…,pN]

(8)

定義干擾(即稀疏信號矢量γ(t)中除了第n個信號γn(t)之外的其他所有信號)的協(xié)方差矩陣為

(9)

(10)

為求得式(10)的最優(yōu)解,將該式展開并處理得如下形式:

(11)

式中,

故關(guān)于γn(t)的加權(quán)最小二乘估計為

(12)

將式(9)代入式(12),根據(jù)矩陣求逆引理可得

(13)

(14)

(15)

相應(yīng)地有

(16)

重復(fù)計算:

forn=1,2,…,N

end

直至收斂。

(17)

(18)

(19)

式中,向量em表示M×M維單位陣IM=[e1,e2,…,eM]的第m列。

綜上,WLS-IAE算法最終形式歸納如下。

初始化:

重復(fù)計算:

forn=1,2,…,N

end

直至收斂。

4 算法計算量分析

l1-min算法和本文提出的WLS-IAE算法都是基于稀疏表示框架的單快拍DOA估計方法,其中,l1-min算法作為依賴于超參數(shù)的稀疏估計類方法具有一定代表性,在其基礎(chǔ)上發(fā)展了很多改進算法。下面詳細分析WLS-IAE算法的計算量,并與l1-min算法進行比較。

令MDN表示復(fù)數(shù)乘除法運算的次數(shù)。首先分析本文算法初始化的計算量,再分析算法迭代過程的計算量。

(20)

結(jié)合WLS-IAE算法初始化過程的計算量,則WLS-IAE算法實現(xiàn)所需的MDN總次數(shù)約為

(21)

l1-min算法一般是將式(6)的凸松弛的約束優(yōu)化問題轉(zhuǎn)化為二階錐規(guī)劃問題的一般形式,并在二階錐規(guī)劃的框架下采用內(nèi)點法求解,其計算復(fù)雜度為O(N3)[18]。

考慮到在實際應(yīng)用中,M?N,nitr≈10,則nitr?N,且一般有nitr?M成立,則WLS-IAE算法的計算復(fù)雜度為O(M2N)。

因此,雖然同為基于稀疏表示框架的DOA估計算法,WLS-IAE算法的計算復(fù)雜度遠低于l1-min算法,具有更好的實時性。

5 仿真實驗

在以下仿真實驗中,考慮由15個陣元構(gòu)成的均勻線性傳感器陣列,相鄰兩個陣元間距為信號波長的一半,即d=λ/2。假設(shè)兩個單位功率的窄帶信號分別從不同方向入射到上述傳感器陣列,實驗中取陣列對信號的單快拍采樣數(shù)據(jù)。SNR定義為

下列實驗中,如無特別說明,取SNR=10 dB。

(1)算法收斂特性

考慮兩個相干信號DOA參數(shù)分別以θ1=-10°,θ2=60°入射到上述均勻線陣。圖1給出了采用WLS-IAE算法分別在初始化、迭代1次、迭代5次、迭代10次過程中對DOA估計的歸一化稀疏功率譜圖。其中算法初始化采用的是CBF算法。兩個相干信號之間的關(guān)系表示為s2(t)=s1(t)?ejη,其中,?ejη為相關(guān)系數(shù),反映了兩個信號之間的數(shù)值關(guān)系。實驗中取?=1,η=π/6。

圖1 WLS-IAE算法中DOA估計的歸一化稀疏功率譜

從圖1的4組仿真曲線可以看出,本文提出的WLS-IAE算法在相干信號的真實目標方向上形成了譜峰。隨著迭代次數(shù)的逐漸增加,稀疏功率譜的譜峰變得越來越尖銳,而“偽峰”及旁瓣逐漸消失,意味著該算法在迭代10次時對目標信號DOA具有較高的估計精度和分辨率。在改變其他實驗參數(shù)的情況下,算法一般不超過15次迭代即可達到收斂狀態(tài),說明本文算法具有較快的收斂速度。

(2)算法對相干信號DOA估計的分辨率

考慮兩個相關(guān)系數(shù)為ejπ/6的鄰近相干信號參數(shù)θ1=-10°,θ2=-6°入射到上述均勻線陣,采用CBF算法、WLA-IAE算法和l1-min算法對單快拍陣列接收數(shù)據(jù)分別進行10次蒙特卡羅仿真實驗,得到3種算法對DOA估計的歸一化稀疏功率譜圖,分別對應(yīng)如圖2(a)～圖2(c)所示。

圖2 角度間隔較小的兩個相干信號DOA估計的歸一化稀疏功率譜

由圖2的仿真曲線可以看出,當相干目標信號的角度間隔較小時,采用WLS-IAE算法和l1-min算法均可以在真實的目標方向形成尖銳的譜峰,實現(xiàn)了對目標信號DOA的高分辨率估計。相反,CBF算法在兩個真實目標方向上的譜峰合二為一,說明此時CBF算法無法分辨兩個信號的DOA。綜上表明,在單快拍情形下,基于稀疏表示的陣列測向方法對空間鄰近信號仍具有較好的分辨能力。

(3)SNR對DOA估計精度的影響

考慮兩個相關(guān)系數(shù)為ejπ/6的相干信號DOA參數(shù)分別以θ1=-10°,θ2=60°入射到上述均勻線陣,采用WLA-IAE算法、l1-min算法、CBF算法、快速迭代內(nèi)插波束形成(fast iterative interpolated beamformer,FIIB)算法[5]和稀疏協(xié)方差矩陣迭代的單快拍(sparse covariance matrix iteration with a single snapshot,SCMISS)算法[15]分別進行200次蒙特卡羅仿真實驗,得到目標信號DOA估計的均方根誤差(root-mean-square error,RMSE)隨SNR的關(guān)系曲線,如圖3所示。其中,L次蒙特卡羅仿真實驗得到的DOA估計的RMSE定義為

圖3 DOA估計的RMSE隨SNR的關(guān)系曲線

從圖3的仿真曲線可以看出,隨著SNR的逐漸提高,采用上述5種算法估計的DOA的RMSE都是逐漸減小的,說明5種算法對DOA的估計精度均越來越高。其中,WLA-IAE算法對DOA的估計精度最高,尤其在低SNR時估計精度優(yōu)勢更明顯,表明WLA-IAE算法在相干信號、單快拍、低SNR等非理想情形下均具有很好的估計性能。

CBF算法的估計性能因受陣列孔徑的限制,顯然不如WLS-IAE、l1-min和SCMISS超分辨算法，它們都屬于基于稀疏表示框架的陣列測向方法,在低SNR情況下都具有良好的估計性能。但是,l1-min算法的稀疏恢復(fù)性能依賴于正則化參數(shù)β的精細化選取,一旦選取不當,算法存在收斂至局部極值的風(fēng)險。雖然本文在第3節(jié)給出了正則化參數(shù)β的選擇依據(jù),但正如文獻[18]所述,這種選擇依據(jù)只適用于高SNR情況下,在低SNR時對正則化參數(shù)β只能憑經(jīng)驗公式選取,而WLS-IAE算法屬于一種不依賴超參數(shù)的稀疏估計類方法,因此WLS-IAE和l1-min算法在高SNR時性能相當,而在低SNR時前者性能優(yōu)于后者。

6 結(jié) 論

本文依托陣列輸出數(shù)據(jù)的稀疏表示模型,提出了一種WLS-IAE算法。詳細分析了WLS-IAE算法的基本性能和計算復(fù)雜度,并與依賴于超參數(shù)的經(jīng)典稀疏估計類算法l1-min的計算復(fù)雜度進行了比較。仿真結(jié)果表明:WLS-IAE算法在低SNR、單快拍、信號相干、目標信號空間間隔很小等非理想情況下都具有很好的估計精度和分辨率;WLS-IAE算法的計算效率明顯高于一般的稀疏估計類方法。因此,WLS-IAE算法不僅保持了稀疏表示框架在非理想條件下的良好估計性能,而且避免了傳統(tǒng)稀疏參數(shù)求解方法的高復(fù)雜度和經(jīng)驗公式的超參數(shù)選取問題,因此非常適用于快變目標信號DOA的實時跟蹤測量,具有潛在的工程實用價值。