王勇，李韜，項(xiàng)建弘

哈爾濱工程大學(xué) 信息與通信工程學(xué)院，黑龍江 哈爾濱 150001

波達(dá)方向（direction of arrival，DOA）估計(jì)作為陣列信號(hào)處理領(lǐng)域中的重要研究方向，在聲吶、雷達(dá)和地震勘測(cè)等領(lǐng)域具有很大的應(yīng)用空間[1-3]。其中比較經(jīng)典的有基于子空間理論的子空間分解類算法，即多重信號(hào)分類（multiple signal classification method，MUSIC）算法[4]和信號(hào)參數(shù)旋轉(zhuǎn)不變（eestimating signal parameter via rotational invariance techniques，ESPRIT）算法[5]。以上2 類算法均需得到信號(hào)協(xié)方差矩陣，再對(duì)其進(jìn)行分解，在良好條件下，可實(shí)現(xiàn)高分辨率DOA 估計(jì)，但是需要大量的快拍數(shù)據(jù)作為估計(jì)精度保障，且無法求解相干源信號(hào)。盡管以上算法可通過平滑處理的方法來解決相干信號(hào)的問題，但會(huì)增加計(jì)算量，且會(huì)損失一定的陣列孔徑，使得分辨率下降。另外，最大似然估計(jì)（maximum likelihood，ML）方法也具有較好的估計(jì)性能，但是需要非線性多維搜索，計(jì)算量較大[6]。近年來，壓縮感知（compressed sensing，CS）理論的提出和廣泛應(yīng)用[7-9]，為DOA估計(jì)提供了新的研究途徑，一些學(xué)者利用信源在空域的稀疏性，直接將壓縮感知的經(jīng)典算法運(yùn)用到波達(dá)方向估計(jì)中。由于計(jì)算量小，以正交匹配追蹤（orthogonal matching pursuit，OMP）為代表的貪婪類算法[10-11]常被用來解決該問題，該類算法的關(guān)鍵在于經(jīng)過反復(fù)迭代，找到測(cè)量矩陣中提供較大能量的原子，從而估計(jì)出信號(hào)方向。但在低信噪比情況下通常找不到最優(yōu)解，容易使得DOA估計(jì)精度下降。另一類方法是將l0范數(shù)最小化問題轉(zhuǎn)化為求解l1范數(shù)最小化問題，如L1-SVD 算法[12]估計(jì)精度高，需要采用線性規(guī)劃進(jìn)行求解，計(jì)算量較大，實(shí)時(shí)處理困難。文獻(xiàn)[13]提出平滑l0范數(shù)重構(gòu)算法，其主要思路是通過構(gòu)造一個(gè)帶參的高斯平滑函數(shù)族來近似信號(hào)的l0范數(shù)，通過最速下降法求解平滑函數(shù)的最優(yōu)解。由于其計(jì)算簡(jiǎn)單，目前已取得了許多的研究成果，如文獻(xiàn)[14]將高斯函數(shù)替代為逼近度更高的雙曲正切函數(shù)去擬合l0范數(shù)，并采用修正牛頓法提高了收斂速度。文獻(xiàn)[15]和文獻(xiàn)[16]中分別提出用一組復(fù)合三角函數(shù)和修正近似雙曲正切函數(shù)來逼近l0范數(shù)，雖然重構(gòu)效果有所提高，但函數(shù)模型比較復(fù)雜，使得迭代計(jì)算復(fù)雜度較高。

基于平滑l0范數(shù)類算法的核心是將求解l0范數(shù)問題轉(zhuǎn)化為求解平滑函數(shù)的極值問題，但是該類算法通常對(duì)噪聲敏感，且多用于單快拍下DOA 估計(jì)。為了進(jìn)一步優(yōu)化算法的適應(yīng)能力，引入了多觀測(cè)矢量(multiple measurement vectors，MMV)模型，采用奇異值分解，先對(duì)信號(hào)降維并提取信號(hào)子空間，一定程度上降低了計(jì)算量。為了使得算法重構(gòu)精度更高，提出擬合程度更高的復(fù)合優(yōu)化函數(shù)去逼近l0范數(shù)，并采用加權(quán)機(jī)制，加速稀疏解的獲取，最后利用最速下降法快速求解。通過實(shí)驗(yàn)論證，所提算法在低信噪比、少快拍下具有更優(yōu)的DOA 估計(jì)誤差。

1 壓縮感知框架下的DOA 估計(jì)模型

假設(shè)空間中有L個(gè)來波信號(hào)方向?yàn)棣萯,i∈{1,2,···,L}的遠(yuǎn)場(chǎng)窄帶信號(hào)si(t),i∈{1,2,···,L}，入射到陣元數(shù)為M的均勻線陣中，其接收范圍為 (-0.5π,0.5π)，陣元間距為d。為了獲得稀疏信號(hào)，將整個(gè)空間劃分為 2N+1份，則整個(gè)空間可表示為 {θ1,θ2,···,θ2N+1}，其中 θk=-0.5π+(k-1)π/2N，k∈{1,2,···,2N+1}。從信號(hào)能量角度分析，由于整個(gè)空間中只有L(L?2N+1)個(gè) 遠(yuǎn)場(chǎng)窄帶信號(hào)，即至多只有L個(gè)不同來波方向上才有真實(shí)的信號(hào)能量分布，因此信號(hào)的真實(shí)方向相對(duì)于整個(gè)信號(hào)空間則是稀疏的?；诳臻g等角度均勻劃分后所構(gòu)造的陣列流型矩陣A是一個(gè)M×(2N+1)維矩陣，在壓縮感知理論中，相當(dāng)于一個(gè)過完備冗余字典，可作為測(cè)量矩陣使用。若陣列接收到L個(gè)遠(yuǎn)場(chǎng)窄帶信號(hào)各自在T個(gè)不同時(shí)刻的快拍數(shù)據(jù)，則信號(hào)的MMV 接收模型可表示為

式中：Y為M×T維陣列接收信號(hào)矩陣；S為(2N+1)×T維遠(yuǎn)場(chǎng)窄帶信號(hào)矩陣，其中只有L(L?2N+1)行不為0，為稀疏信號(hào)矩陣。對(duì)于S中的列向量s而言，s中只有L(L?2N+1)個(gè)位置不為0，為稀疏向量；G為M×T維的高斯白噪聲矩陣；A=[a(θ1)a(θ2) ···a(θ2N+1)]，為M×(2N+1)維空間等角度網(wǎng)格化構(gòu)造的陣列流型矩陣，其中a(θk)=為導(dǎo)向矢量，k∈{1,2,···,2N+1}，λ 為信號(hào)的波長(zhǎng)，d為陣元間距。

壓縮感知框架下的DOA 估計(jì)通常是在獲得稀疏網(wǎng)格化構(gòu)造的陣列流型矩陣A和陣列接收信號(hào)矩陣Y之后，利用重構(gòu)算法對(duì)入射信號(hào)矩陣S進(jìn)行恢復(fù)，從而找出S中對(duì)應(yīng)的非零位置集合并將其映射到劃分的網(wǎng)格空間，最終求解出的θj,j∈{1,2,···,L}即為實(shí)際的信號(hào)來波方向。

2 改進(jìn)平滑l0 范數(shù)的DOA 估計(jì)算法

針對(duì)式（1）中信號(hào)接收矩陣Y的列向量y，以及稀疏矩陣S中的稀疏列向量s，在空間存在噪聲情況時(shí)，可轉(zhuǎn)化為下述優(yōu)化問題進(jìn)行求解：

圖1 和圖2 分別在 ρ=0.2 和 ρ=0.1處直觀地觀察高斯函數(shù)fρ(si)與所提復(fù)合優(yōu)化函數(shù)gρ(si)對(duì)于l0范數(shù)的逼近程度，圖中H(si)代表函數(shù)值。

圖1 ρ=0.2時(shí)函數(shù)對(duì)于 l0范數(shù)的逼近程度

圖2 ρ=0.1時(shí) 函數(shù)對(duì)于 l0范數(shù)的逼近程度

可知在 ρ →0這一不斷變小的過程中，高斯函數(shù)和所提優(yōu)化函數(shù)對(duì)于l0范數(shù)的逼近程度均有所提高，并且當(dāng)si=0 時(shí)2 個(gè)函數(shù)的值均為0；在si不為0 的附近，由于所提函數(shù)對(duì)應(yīng)的值能更快地趨向1，相對(duì)于高斯函數(shù)的圖形更為“陡峭”，即所提優(yōu)化函數(shù)逼近l0范數(shù)程度要高于高斯函數(shù)，這也驗(yàn)證了所提優(yōu)化函數(shù)用于逼近l0范數(shù)的優(yōu)越性。因此式（2）可轉(zhuǎn)化為式（3）的最優(yōu)化問題：

為了加速獲得稀疏解，提高算法的收斂性，對(duì)所提函數(shù)采用加權(quán)機(jī)制，通常在迭代過程中，對(duì)于較大的元素si給與一個(gè)較小的權(quán)值，而較小的元素si給定一個(gè)較大的權(quán)值，本文采用的權(quán)值函數(shù)表形式為

可知此時(shí)當(dāng)si較大時(shí)，權(quán)值wi較小，符合較大的元素對(duì)應(yīng)較小的權(quán)值這一要求。

令

作為權(quán)值函數(shù)的對(duì)角矩陣，在迭代尋優(yōu)過程中，將權(quán)值函數(shù)加權(quán)到目標(biāo)稀疏信號(hào)

上，可得：

在后續(xù)的優(yōu)化迭代過程中，較小的加權(quán)值能夠一定程度減緩目標(biāo)信號(hào)中較大元素的下降，而較大的加權(quán)值能加速目標(biāo)信號(hào)中較小元素的下降，因此針對(duì)稀疏信號(hào)s而言，目標(biāo)信號(hào)的大多數(shù)元素能更快的趨近于0。由于在噪聲環(huán)境下，來波信號(hào)主要集中在元素值較大的位置，因此引入加權(quán)機(jī)制在一定程度上提升了算法的收斂性。

通過加權(quán)函數(shù)的引入，式（3）可轉(zhuǎn)化為式（4）優(yōu)化問題：

對(duì)于式（4），本文采用最速下降法進(jìn)行求解，稀疏信號(hào)向量s的迭代更新方向記為

式中：μ是步長(zhǎng)，為一個(gè)較小的常數(shù)；d為最速下降法的下降方向。

根據(jù)目標(biāo)函數(shù)gρ(si)可知其一階導(dǎo)數(shù)為

根據(jù)式（4）優(yōu)化的方向進(jìn)行優(yōu)化迭代，并將每次迭代的結(jié)果進(jìn)行梯度投影，其投影方式為

綜上所提算法主要包括2 個(gè)循環(huán)，首先確定一個(gè)遞減序列 ρ1,ρ2,···,ρj，外循環(huán)控制 ρ的取值由ρ1變化到 ρj，內(nèi)循環(huán)則是針對(duì)每一個(gè)Ysv，計(jì)算函數(shù)的優(yōu)化方向d并對(duì)信號(hào)進(jìn)行迭代更新和進(jìn)行梯度投影，從而最終實(shí)現(xiàn)信號(hào)的重構(gòu)。

具體的計(jì)算步驟：

本文針對(duì)的是多快拍模型，由于上述流程只是針對(duì)Ysv的某一列向量ysv求得的重構(gòu)解，因此需要對(duì)Ysv所有的列向量分別重復(fù)上述步驟求出各自重構(gòu)解，然后再求平均值，得到多快拍模型下的稀疏重構(gòu)解，最終再映射到劃分網(wǎng)格上，獲得DOA 估計(jì)值。

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

在Matlab 平臺(tái)上進(jìn)行實(shí)驗(yàn)仿真驗(yàn)證本文所提算法的可行性，并與SL0、OMP、L1-SVD 算法進(jìn)行對(duì)比，分析所提算法在低信噪比、少快拍下的DOA 估計(jì)性能優(yōu)勢(shì)。本實(shí)驗(yàn)中均采用間距為半波長(zhǎng)的均勻線陣，并將空間從 (-90?,90?)等角度均勻劃分1 81份，每份間隔1°，噪聲環(huán)境為高斯白噪聲。估計(jì)指標(biāo)采用均方根誤差和估計(jì)成功率，其中DOA 估計(jì)的均方根誤差定義如下：

式中：L為空間目標(biāo)信號(hào)源個(gè)數(shù)，K為蒙特卡洛實(shí)驗(yàn)次數(shù)，為第k次蒙特卡洛實(shí)驗(yàn)中第l個(gè)信源的DOA 估計(jì)值，θl為第l個(gè)信源的真實(shí)來波方向。

DOA 估計(jì)成功率定義為

式中F(θl)為 入射方向?yàn)?θl的信源在滿足估計(jì)誤差RMSE下，能夠被估計(jì)出來的總次數(shù)。本文實(shí)驗(yàn)假設(shè)估計(jì)值與真實(shí)值的RMSE小于或等于1°時(shí)，則認(rèn)為估計(jì)成功。

實(shí)驗(yàn)1驗(yàn)證所提算法在不同快拍條件下的DOA 估計(jì)結(jié)果有效性。假設(shè)陣元數(shù)M=16，信噪比 為5 dB，有3 個(gè)目標(biāo)源 信號(hào)，分別從-30?、0?、20?入射，其中-30?和 0?為相干信號(hào)源。分別在單快拍和快拍數(shù)為10 的條件下進(jìn)行仿真，并對(duì)空間譜進(jìn)行了歸一化處理，最終DOA 估計(jì)結(jié)果如圖3和圖4 所示。

圖3 單快拍下DOA 估計(jì)結(jié)果

圖4 快拍數(shù)為10 的 DOA 估計(jì)結(jié)果

由圖3 和圖4 可知，歸一化空間譜在信號(hào)的入射方向均出現(xiàn)了峰值，且隨著快拍數(shù)的增加，歸一化空間譜在非信號(hào)源方向上的波動(dòng)更少，在信源方向上的譜峰相對(duì)而言更尖銳、寬度更窄，更容易區(qū)分信號(hào)的來波方向，并且在前2 個(gè)為相干信號(hào)源的方向上也各自出現(xiàn)了譜峰，因此所提算法在較少快拍，甚至單快拍下均能獲取準(zhǔn)確的DOA 估計(jì)結(jié)果，估計(jì)算法不受相干源信號(hào)影響。

實(shí)驗(yàn)2比較在不同信噪比下各算法的DOA估計(jì)性能。假設(shè)目標(biāo)信源的入射角度為 0?、20?，陣元數(shù)量M=16，快拍采樣數(shù)為10，進(jìn)行600 次蒙特卡洛獨(dú)立重復(fù)試驗(yàn)。同L1-SVD、OMP、SL0 算法進(jìn)行比較，與本文所提算法在不同信噪比下的DOA 估計(jì)誤差和估計(jì)成功率的仿真如圖5 和圖6所示。

圖5 不同信噪比下算法的DOA 估計(jì)誤差

圖6 不同信噪比下算法的DOA 估計(jì)成功率

在不同信噪比條件下，由圖5 可知，高信噪比下，各算法均方根誤差均逐漸趨于穩(wěn)定并收斂到0，當(dāng)信噪比低于-4 dB 時(shí)，所提算法相比其他3 種算法具有明顯優(yōu)勢(shì)，其均方根誤差明顯低于其他3 種算法。由圖6 可以看出，信噪比為-3 dB 時(shí)，各算法估計(jì)成功率均高于90%，所提算法趨于100%，在低信噪比下所提算法具有更優(yōu)的估計(jì)性能。

實(shí)驗(yàn)3比較在不同快拍下各算法的DOA 估計(jì)性能。假設(shè)2 個(gè)信源來波方向分別是 0?、20?，陣元數(shù)量16，信噪比RSN=-3 dB，進(jìn)行600 次蒙特卡洛實(shí)驗(yàn)，各算法在不同快拍下的DOA 估計(jì)結(jié)果如圖7 和圖8 所示。

圖7 不同快拍下算法的DOA 估計(jì)誤差

圖8 不同快拍下算法的DOA 估計(jì)成功率

在不同的快拍條件下，由圖7 可以看出，即使在-3 dB 低信噪比條件下，隨著快拍數(shù)的不斷增加，各算法的估計(jì)誤差也能逐漸降低并趨于穩(wěn)定；在快拍數(shù)少于6 時(shí)，估計(jì)的均方根誤差都不同程度的增大，上述所有算法的性能都開始下降，但所提算法相比其他算法的均方根估計(jì)誤差明顯更小。由圖8 可以看出，隨著快拍數(shù)的增加，各算法的估計(jì)成功率也在緩慢增加并趨于穩(wěn)定；當(dāng)快拍數(shù)為10 時(shí)，所提算法幾乎能以100%的概率估計(jì)出信號(hào)的來波方向；在快拍數(shù)為4 時(shí)，其估計(jì)成功率也在90%左右。相比其他算法，本文提出的算法在低信噪比、少快拍下具有更好的估計(jì)成功率。

4 結(jié)論

本文提出了一種加權(quán)改進(jìn)的平滑l0范數(shù)最小化算法，并將其應(yīng)用到多觀測(cè)矢量DOA 估計(jì)中。

1）為了提高算法重構(gòu)精度，提出了新的復(fù)合優(yōu)化函數(shù)來逼近l0范數(shù)。經(jīng)過分析可知，相比原始的SL0 算法，其高斯函數(shù)逼近程度更高。

2）為了提高算法的穩(wěn)定性，本文算法采用了加權(quán)機(jī)制，使得其在迭代過程中能加速稀疏解的獲取，并且在重構(gòu)前可預(yù)先對(duì)接收數(shù)據(jù)進(jìn)行奇異值分解，提取了信號(hào)子空間，使得算法整體的抗噪性能得到提升，也一定程度上降低了計(jì)算量。

經(jīng)過實(shí)驗(yàn)對(duì)比，本文所提方法相比OMP、原始的SL0 和L1-SVD 算法在低信噪比、少快拍下具有更優(yōu)的DOA 估計(jì)性能。因此在低信噪比、少快拍下的實(shí)際應(yīng)用中，本文所提算法具有更高的可靠性和適用性。但本文只針對(duì)一維DOA 估計(jì)，后續(xù)仍需推廣和應(yīng)用到二維DOA 估計(jì)中，以提高算法的適用性。同時(shí)在后續(xù)研究中也需要考慮陣列流型存在幅相誤差這一問題。