吉愛國 劉偉平 劉志強

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

0 引 言

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

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

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

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

1 陣列信號模型

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

圖1 陣列信號DOA估計模型

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

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

(1)

信號X(t)的協(xié)方差矩陣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)

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

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

(5)

式中:(·)*表示矩陣的共軛。比較式(5)與式(1)，可以發(fā)現(xiàn)式(5)中A*⊙A為新的M2×K維陣列流型矩陣，稱其為虛擬陣列流型矩陣，y對應(yīng)虛擬陣元，p對應(yīng)虛擬信源。因此通過KR積變換，將陣列多測量矢量(Multiple Measurement Vectors，MMV)模型轉(zhuǎn)換為虛擬陣列單測量矢量(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進(jìn)行降維[14]：

(8)

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

2.2 實數(shù)化的虛擬陣列

(9)

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

(10)

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

(11)

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

(12)

3 DOA估計算法

(13)

(14)

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

(15)

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

(16)

存在噪聲時，有：

(17)

(18)

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

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

(19)

(20)

(21)

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

(22)

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

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

y=vec(R)

(2) 對y進(jìn)行降維、實數(shù)化和去噪。

(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 仿真實驗

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

(23)

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

圖2 RMSE隨SNR的變化

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

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

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

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

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

5 結(jié) 語

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