基于近似l0范數(shù)的實(shí)數(shù)化DOA估計(jì)

吉愛(ài)國(guó) 劉偉平 劉志強(qiáng)

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

0 引 言

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

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

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

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

1 陣列信號(hào)模型

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

圖1 陣列信號(hào)DOA估計(jì)模型

第t次快拍時(shí)，陣列接收到的信號(hào)為:

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

(1)

信號(hào)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積變換

對(duì)于兩個(gè)矩陣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積。對(duì)于兩個(gè)向量a=[a1,a2,…,an]T和b=[b1,b2,…,bm]T，其Kronecker積可以表示為:

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

(4)

假設(shè)空間信號(hào)之間互不相干，則式(2)中的RS為對(duì)角陣，即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維陣列流型矩陣，稱(chēng)其為虛擬陣列流型矩陣，y對(duì)應(yīng)虛擬陣元，p對(duì)應(yīng)虛擬信源。因此通過(guò)KR積變換，將陣列多測(cè)量矢量(Multiple Measurement Vectors，MMV)模型轉(zhuǎn)換為虛擬陣列單測(cè)量矢量(Single Measurement Vectors，SMV)模型。

將式(5)改寫(xiě)為：

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可對(duì)y進(jìn)行降維[14]：

(8)

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

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

(9)

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

(10)

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

(11)

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

(12)

3 DOA估計(jì)算法

(13)

(14)

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

(15)

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

(16)

存在噪聲時(shí)，有：

(17)

(18)

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

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

(19)

(20)

(21)

由于牛頓方向中的Hesse矩陣(即▽2Fσ(x))并非正定矩陣，無(wú)法保證牛頓方向?yàn)橄陆捣较颍市鑼?duì)其進(jìn)行修正，構(gòu)造一個(gè)新矩陣G=▽2Fσ(x)+εkI，其中:εk是一個(gè)合適的正數(shù)，I是單位矩陣。取：

(22)

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

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

y=vec(R)

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

(3) 使用迭代算法求解DOA估計(jì)問(wèn)題。

Forj=1,2,…,J

σ=σj

Forl=1,2,3,4,5

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

修正牛頓法x=x+d

end

σ=ρσj

end

4 仿真實(shí)驗(yàn)

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

(23)

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

圖2 RMSE隨SNR的變化

圖3 DOA估計(jì)成功率隨SNR的變化

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

圖4 DOA估計(jì)成功率隨N的變化

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

表1 不同算法計(jì)算時(shí)間的比較

5 結(jié) 語(yǔ)

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