基于近似l0范數的實數化DOA估計

吉愛國 劉偉平 劉志強

(青島理工大學信息與控制工程學院 山東 青島 266520)

0 引 言

波達方向估計作為陣列信號處理的重要研究內容，在雷達、射電天文、聲吶、測向等領域起著重要的作用[1]。常見的DOA估計方法可分為三類[2-4]：波束形成法(Beam Forming，BF)，子空間法和最大似然法(Maximum Likelihood，ML)。常規波束形成法計算速度快，但存在“瑞利限”，其估計性能受到陣列孔徑的約束。最小方差無畸變響應(Minimum Variance Distortionless Response，MVDR)算法提高了估計分辨率，但在低信噪比下估計效果較差[2]。最大似然法估計性能較好，但其估計過程中需要非線性多維搜索，計算量較大[4]。以多重信號分類(Multiple Signal Classification，MUSIC)算法為代表的子空間類算法，利用信號子空間和噪聲子空間的正交性來構造空間譜函數，通過譜峰搜索來檢測信號的DOA，具有較高的分辨率，但在低信噪比和快拍數較少時性能會下降[3]。

壓縮感知理論為DOA估計問題帶來了新的解決方法，使得利用少量快拍數據進行DOA估計有了可能。對信源所在的空間進行網格劃分，使信源個數遠小于網格數，便可使用壓縮感知中的稀疏重構方法進行DOA估計。壓縮感知恢復過程的本質是求解l0范數(即非零元素的個數)，但是求解l0范數NP難，因此Chen等[5]提出基追蹤(Basis Pursuit，BP)算法，利用l1范數代替l0范數進行稀疏重構。Malioutov等[6]將l1范數用于DOA估計中，提出l1-SVD算法，先對陣列接收的數據矩陣進行奇異值分解，再通過求解二階錐規劃問題進行DOA估計，然而該算法容易出現偽峰，影響估計性能。韓樹楠等[7]提出了一種加權l1-SVD方法，利用接收數據傅里葉變換后的空間譜來獲取權值，解決了l1-SVD算法的偽峰問題，但該算法穩定性較差，需要多次采樣進行平滑處理。蔡晶晶等[8]將Khatri-Rao(KR)積理論、空域平滑理論和壓縮感知理論相結合，實現了相干信號的DOA估計。

利用l1范數求解稀疏重構反問題精度高，但是計算量較大。近年來，一些學者提出使用光滑函數逼近l0范數，將其轉化為最優化問題并通過迭代算法求解，提高了計算速度[9-13]。Mohimani等[9]提出SL0算法，采用帶參數的高斯函數來逼近l0范數，將求解l0范數轉化為光滑函數的最值問題，重建速度比基追蹤算法快2～3倍。趙瑞珍等[10]提出NSL0算法，使用雙曲正切函數逼近l0范數，選用修正牛頓方法求解最優化問題，提高了運算效率。孫娜等[11]利用最速下降法和擬牛頓法進行稀疏重構，用BFGS公式計算Hesse矩陣的近似值，使得重構精度和峰值信噪比等得到了提高。

然而上述近似l0范數算法僅能用于單快拍、無噪、實數情況下的稀疏重構問題，本文提出了一種基于近似l0范數的實數化DOA估計算法(AL0-DOA)，使得近似l0范數算法可用于多快拍、存在噪聲、復數運算的DOA估計中。先利用Khatri-Rao(KR)積變換，將陣列多測量矢量模型轉換為虛擬陣列單測量矢量模型，再通過降維和實數化進一步降低計算量，同時抑制了噪聲，最后引入平滑函數來近似l0范數，通過修正牛頓算法求解最優化問題得到DOA的估計值。仿真結果表明本文算法計算快，精度較高，可對DOA進行有效估計。

1 陣列信號模型

考慮一M元各向同性的等距線陣(Uniform Linear Array，ULA)，陣元間距d為λ/2，λ是信號波長。假設K個遠場窄帶信號源以平面波入射到線陣上，信源方向分別為θ1,θ2,…,θK，建立如圖1所示的陣列模型。

圖1 陣列信號DOA估計模型

第t次快拍時，陣列接收到的信號為:

X(t)=AS(t)+N(t)t=1,2,…,N

(1)

信號X(t)的協方差矩陣R為:

R=E[X(t)XH(t)]=

AE[S(t)SH(t)]AH+E[N(t)NH(t)]=

ARSAH+I

(2)

2 KR積虛擬陣列模型

2.1 KR積變換

對于兩個矩陣A=[a1,a2,…,ak]∈Cn×k和B=[b1,b2,…,bk]∈Cm×k，其KR積定義為:

A⊙B=[a1?b1,a2?b2,…,ak?bk]∈Cnm×k

(3)

式中:“⊙”和“?”分別表示KR積和Kronecker積。對于兩個向量a=[a1,a2,…,an]T和b=[b1,b2,…,bm]T，其Kronecker積可以表示為:

a?b=[a1b,a2b,…,anb]T=vec(baT)

(4)

假設空間信號之間互不相干，則式(2)中的RS為對角陣，即RS=diag(p)，其中p=[p1,p2,…,pK]T，矢量化式(2)可得：

(A*⊙A)p+vec(I)

(5)

式中:(·)*表示矩陣的共軛。比較式(5)與式(1)，可以發現式(5)中A*⊙A為新的M2×K維陣列流型矩陣，稱其為虛擬陣列流型矩陣，y對應虛擬陣元，p對應虛擬信源。因此通過KR積變換，將陣列多測量矢量(Multiple Measurement Vectors，MMV)模型轉換為虛擬陣列單測量矢量(Single Measurement Vectors，SMV)模型。

將式(5)改寫為：

y=(A*⊙A)p+vec(I)=GBp+vec(I)

(6)

式中：A*⊙A=GB，G∈CM2×(2M-1)為列正交矩陣；B=[b(θ1),b(θ2),…,b(θK)]∈C(2M-1)×K，b(θk)=[ej(M-1)ωk,ejMωk,…,ejωk,1,e-jωk,…,e-jMωk,e-j(M-1)ωk]T。

(7)

令W=GTG=diag(1,2,…,M-1,M,M-1,…,2,1)，式(6)左乘W-1/2GT可對y進行降維[14]：

(8)

比較式(8)與式(5)，可以發現式(8)中W1/2B為陣列流型矩陣，其維數由M2×K縮減為(2M-1)×K，減少了冗余信息，降低了運算量。

2.2 實數化的虛擬陣列

(9)

H1W-1/2GTvec(I)=0H2W-1/2GTvec(I)=0

(10)

故將式(8)分別左乘H1和H2得：

(11)

式中：BR和BI均為實數。綜上所述，通過左乘矩陣H可將復數運算轉化為實數運算，同時除去了噪聲，可表示為：

(12)

3 DOA估計算法

(13)

(14)

求解式(14)可利用CVX工具箱，但其計算量較大。文獻[9]提出使用高斯函數來近似l0范數，采用最速下降法求解式(13)，本文稱該方法為KR-SL0算法。為了進一步提高近似效果，可以采用式(15)所示的雙曲正切函數作為近似l0范數。

(15)

式中：σ為一個很小的正數，它控制著函數的光滑程度；xi表示向量x的第i個元素。則有：

(16)

存在噪聲時，有：

(17)

(18)

當σ的取值很小時，函數Fσ(x)存在很多局部最小值點，很難求出其全局最優解；σ的取值越大，函數Fσ(x)越光滑，越不能精確估計l0范數。所以選擇一個遞減的序列σ1,σ2,σ3,…，對σ的每一個取值，求函數Fσ(x)的最優值，直到σ足夠小，避免陷入函數Fσ(x)的局部最小值。

本文使用修正牛頓法求解式(18)的問題[10]。首先計算Fσ(x)的牛頓方向d=-▽2Fσ(x)-1▽Fσ(x)。

(19)

(20)

(21)

由于牛頓方向中的Hesse矩陣(即▽2Fσ(x))并非正定矩陣，無法保證牛頓方向為下降方向，故需對其進行修正，構造一個新矩陣G=▽2Fσ(x)+εkI，其中:εk是一個合適的正數，I是單位矩陣。?。?/p>

(22)

綜上所述，使用AL0-DOA算法進行DOA估計的步驟如下:

(1) 計算數據的協方差矩陣的估計值，并矢量化。

y=vec(R)

(2) 對y進行降維、實數化和去噪。

(3) 使用迭代算法求解DOA估計問題。

Forj=1,2,…,J

σ=σj

Forl=1,2,3,4,5

計算修正牛頓方向d=-G-1▽Fσ(x)

修正牛頓法x=x+d

end

σ=ρσj

end

4 仿真實驗

使用均勻直線陣列，陣元間距為半波長，信源為隨機高斯信號，噪聲為均值為零的高斯白噪聲，將空間角度從-90°至90°以1°為間隔均勻劃分為181個網格。AL0-DOA算法中σ的初始值為2，σ的終止值為10-4，σ的減少步長σ為0.7。均方誤差(RMSE)定義為:

(23)

假設有3個非相干窄帶信源，其入射角度為[-50°, 0°, 55°]，陣元個數為12，快拍數N=500，信噪比SNR從-6～6 dB以2 dB為間隔，每個信噪比下進行500次蒙特卡洛實驗。圖2和圖3分別為KR-SL0、KR-L1和AL0-DOA算法的DOA估計均方誤差RMSE和成功率隨SNR的變化曲線。可以看出，隨著信噪比的增加，幾種算法的均方誤差都逐漸減小，估計成功率逐漸增加。本文提出的AL0-DOA算法的RMSE始終比KR-SL0和KR-L1算法低，DOA估計成功率始終比KR-SL0和KR-L1算法高。信噪比為-2 dB時，本文算法估計成功率達到90%以上；信噪比大于2 dB時，其估計誤差為0，與其他算法相比，本文算法在低信噪比下的表現更好。

圖2 RMSE隨SNR的變化

圖3 DOA估計成功率隨SNR的變化

假設有4個非相干窄帶信源，其入射角度為[-50°, -20°, 20°，55°]，陣元個數為12，信噪比SNR為15 dB。圖4為KR-SL0、KR-L1和AL0-DOA算法的DOA估計成功率隨快拍數N的變化曲線。可以看出，隨著快拍數的增加，幾種算法的估計成功率逐漸增加。本文算法DOA估計成功率始終比KR-SL0和KR-L1算法高，在快拍數為50時，估計成功率達到了90%以上，說明本文算法可以在更少的快拍數下取得更優的估計性能。

圖4 DOA估計成功率隨N的變化

假設有3個非相干窄帶信源，其入射角度為[-50°, 0°, 55°]，陣元個數為12，快拍數N=500，信噪比SNR為15 dB。對KR-SL0、KR-L1、L1-SVD和AL0-DOA每種算法進行10次實驗，CPU主頻為3.2 GHz，仿真軟件為MATLAB 2014a。表1為每種算法10次實驗的總運算時間。結合上述實驗結果分析表中數據可以發現：L1-SVD 算法需要求解二階錐規劃問題，因而所需時間最長；KR-SL0算法計算時間最短，但是其估計性能較差；KR-L1算法運行時間較長，而且其估計性能一般；本文算法具有運算時間短，估計性能高的優點。

表1 不同算法計算時間的比較

5 結 語

本文提出了一種基于近似l0范數的實數化DOA估計算法，將無法直接求解l0范數的問題轉化為平滑函數的最優化問題，并將其用于DOA估計中，通過修正牛頓算法快速求解，解決了DOA估計中的多快拍、含噪、復數運算的問題。利用KR積理論對數據協方差矩陣進行變換，將陣列多測量矢量模型轉換為虛擬陣列單測量矢量模型。降維和實數化過程進一步降低了運算量，去噪過程提高了算法的魯棒性。仿真結果表明，本文算法可對DOA進行有效估計，克服了l1范數計算量大的問題，保留了l0范數運算快的優點。