應(yīng)用雙線性模型的離格波達(dá)角估計方法*

王友順，高童迪，袁正道,2，丁永春

(1.河南開放大學(xué) 人工智能工程研究中心，鄭州450002；2.鄭州大學(xué) 信息工程學(xué)院，鄭州450001；3.中國船舶集團(tuán)有限公司第七一三研究所，鄭州450001)

0 引 言

利用傳感器陣列估計多個遠(yuǎn)場信號的波達(dá)角(Direction of Arrival,DOA)是陣列信號處理的基礎(chǔ)問題，在雷達(dá)、聲吶和地震預(yù)測等領(lǐng)域得到了廣泛的應(yīng)用[1-2]。針對DOA估計問題，國內(nèi)外研究者提出了多種不同類型的方法，如子空間分解、支持向量機(jī)、矩陣特征空間分類等[3-4]。

隨著壓縮感知技術(shù)的出現(xiàn)[5]，研究者們將DOA估計歸納為稀疏恢復(fù)問題[6]，并提出了多種不同的估計算法。例如，基于正交匹配追蹤算法[7]，將波達(dá)角范圍劃分為均勻網(wǎng)格，假定實際的波達(dá)角準(zhǔn)確位于某個網(wǎng)格上。相比于子空間劃分等其他方法，壓縮感知類方法通常能夠獲得更好的性能，特別是在快拍數(shù)量較少的情況下。稀疏貝葉斯學(xué)習(xí)[8]是一類經(jīng)典稀疏估計算法，通過給變量增加稀疏引導(dǎo)先驗，能夠以較少的觀測實現(xiàn)更優(yōu)的估計性能。但是當(dāng)真實來波方向并不精確位于網(wǎng)格上時，上述方案會導(dǎo)致較嚴(yán)重的估計偏差。更精細(xì)的網(wǎng)格劃分一方面會導(dǎo)致較高的運(yùn)算復(fù)雜度，另一方面會使感知矩陣的列向量之間具有較高的相關(guān)度[9]，導(dǎo)致估計算法收斂性變差。

近年來國內(nèi)外多個團(tuán)隊提出了基于離格信號的稀疏重構(gòu)算法。文獻(xiàn)[9]提出了一種基于稀疏最小二乘的離格算法，假定真實角度不必限于網(wǎng)格之上，能夠計算出誤差并對估計值進(jìn)行修正。文獻(xiàn)[10]則提出了一種求根稀疏貝葉斯學(xué)習(xí)方法，在線性陣列場景下能夠取得很好的估計性能。由于復(fù)雜度較低，經(jīng)典的近似消息傳遞算法(Approximate Message Passing,AMP)在稀疏貝葉斯學(xué)習(xí)框架下得到了廣泛應(yīng)用[11]，但是普通的AMP算法在遇到較“壞”的觀測矩陣時(如病態(tài)條件數(shù)、列相關(guān)、低秩等)會導(dǎo)致迭代過程中出現(xiàn)發(fā)散[12]。特別是在DOA場景下，陣列流形具有較強(qiáng)的相關(guān)性，普通的AMP算法不穩(wěn)定。為了解決上述問題，多種不同的改進(jìn)算法相繼被提出：文獻(xiàn)[13]提出了一種向量化AMP算法(Vector-AMP)，通過兩次標(biāo)量化-期望傳播方法，大大提高了魯棒性；在此基礎(chǔ)上，文獻(xiàn)[14]歸納了雙線性VAMP算法(Bilinear Adaptive-VAMP,BAd-VAMP)，將VAMP算法擴(kuò)展到了雙線性場景，在字典矩陣學(xué)習(xí)、矩陣補(bǔ)全等多個場景下得到了應(yīng)用。

為解決離格波達(dá)角估計中復(fù)雜度和魯棒性的矛盾，本文提出將DOA估計模型轉(zhuǎn)換為雙線性問題，利用文獻(xiàn)[14]中所提BAd-VAMP算法解決該問題。本文的創(chuàng)新點如下：第一，提出將離格DOA估計模型轉(zhuǎn)換為雙線性問題，并利用BAd-VAMP算法對雙線性模型進(jìn)行推導(dǎo)和運(yùn)算；第二，充分利用離格DOA模型中待估計向量的相同稀疏特征減小迭代中矩陣的維度，降低了迭代算法的復(fù)雜度；第三，提出一種基于殘差的冷/熱重啟的迭代機(jī)制，避免了雙線性算法容易陷入鞍點的問題。在仿真部分，本文建立DOA估計模型對所提算法與文獻(xiàn)中已有方法進(jìn)行對比，結(jié)果表明所提算法在機(jī)載雷達(dá)等快拍數(shù)較少的場景下具有明顯的性能優(yōu)勢。

本文符號說明：Y∈N×L表示維度為N×L的復(fù)矩陣，IN表示N維單位陣，(·)T和(·)H分別代表轉(zhuǎn)置和共軛轉(zhuǎn)置，CN(x;μ,ν)表示均值為μ方差為ν的復(fù)高斯分布，x∈U(a,b)表示變量x服從區(qū)間(a,b)的均勻分布，對向量x取平均表示為〈x〉，矩陣或向量的2范數(shù)表示為‖‖2，tr(A)表示矩陣的求跡運(yùn)算，diag(a)表示利用向量a構(gòu)建對角陣。

1 離格DOA模型

(1)

式中：yt∈M×1為第t個快拍下M個陣元的接收向量；K×1表示第t個快拍下的K個來波信號；w表示均值為0、方差為1/λ的高斯白噪聲，其中λ表示噪聲精度。由于噪聲方差不隨時間變化，本文中省略w的上標(biāo)t。陣列流形表示為

A(θ)=[a(θ1),a(θ2),…,a(θN)]∈M×N。

要提高估計精度，最直接的方法是將信號空間劃分為更小的網(wǎng)格，即加大N。但更大的N會導(dǎo)致復(fù)雜度提升，且陣列流形A(θ)的列相關(guān)性變大。相關(guān)性高的陣列流形A(θ)會導(dǎo)致估計性能變差，對于消息傳遞等迭代類算法甚至?xí)鸢l(fā)散，導(dǎo)致估計失敗。

(2)

式中：

e(θn)=a′(θn)為誤差矢量，其中a′(θn)表示對導(dǎo)向矢量a(θn)依θn求導(dǎo)。式(2)所示的接收模型可以重寫為

yt=(A(θ)+E(θ)diag(β))st+w。

(3)

式中：誤差矩陣和誤差向量分別定義為E(θ)=[e(θ1),e(θ2),…,e(θN)]和β=[β1,β2,…βN]T。當(dāng)DOA估計模型中有L個快拍時，上式可以寫為多重觀測形式：

Y=(A(θ)+E(θ)diag(β))S+W。

(4)

式中：矩陣Y=[y1,y2,…,yL]和S=[s1,s2,…,sL]分別表示觀測矩陣和待估計矩陣，且S具有相同列稀疏特征，即每列非零元素位置相同。由于上式中向量θ為固定值，則可簡記為

Y=(A+Ediag(β))S+W。

(5)

值得注意的是，上述模型中未知參數(shù)有S和β，在信號處理領(lǐng)域，通常將利用一組觀測同時估計兩組獨(dú)立參數(shù)的模型稱為雙線性模型，在下一節(jié)中將式(5)重新建模為雙線性形式。

2 雙線性系統(tǒng)建模

通用雙線性模型可以表示為[14-15]

Y=Φ{β}C+W。

(6)

式中：Y=[y1,y2,…,yL]和C=[c1,c2,…,cN]分別表示觀測和待估計矩陣。目前針對上述雙線性系統(tǒng)的方法有俄亥俄州立大學(xué)Schniter教授提出的PBiGAMP[15]和BAd-VAMP[14]兩種算法。根據(jù)文獻(xiàn)中的仿真結(jié)果可知BAd-VAMP從復(fù)雜度和性能上全面超越PBiGAMP算法，所以本文選擇BAd-VAMP算法作為算法框架。

通過對式(5)和式(6)的類比，可以將式(5)所示的DOA估計轉(zhuǎn)換為雙線性模型。定義式(5)中的流形和誤差矩陣為

Φ{β}=A+∑nβnEn。

式中：En=[0,…,e(θn),…,0]∈M×N。則DOA估計模型可以寫為雙線性形式：

Y=(A+∑nβnEn)S+W=Φ{β}S+W。

(7)

根據(jù)式(7)所示的數(shù)據(jù)之間的約束關(guān)系，可以對系統(tǒng)中所有已知和未知變量的聯(lián)合分布函數(shù)利用貝葉斯公式進(jìn)行因式分解[12-13]：

p(Y,S,β,λ)=p(Y|S,β,λ)p(S|γ)p(γ)p(β)p(λ)=

(8)

3 基于BAd-VAMP的DOA估計算法

BAd-VAMP算法可以分為兩部分，分別針對參數(shù)S和β估計，本節(jié)也按照以上劃分進(jìn)行推導(dǎo)。

3.1 參數(shù)S估計部分

定義

(9)

定義

(10)

由于函數(shù)節(jié)點fδ為等效節(jié)點，則消息

其中均值方差分別為

(11)

然后再利用MF規(guī)則估計超先驗參數(shù)γ為

(12)

其中期望和方差分別計算為

(13)

從而完成S估計部分的迭代計算。

3.2 參數(shù)β估計部分

E[lnp(Y|S2,β,λ)]=

Const-λ(βHHβ-βHα-αHβ)。

(14)

式中：H∈N×N，α∈N×1，Const表示常數(shù)。式(14)中矩陣H和向量α中每個元素計算為

(15)

(16)

其中：i,j∈[1:N]，

(17)

(18)

式中：β為前一次迭代中求得的估計值，據(jù)此可得到Φ{β}。對式(14)所示期望依β求梯度，并令其等于零，可得到使期望取最大值的β值，即

βmax=H-1α。

3.3 算法改進(jìn)

(1)迭代阻尼

(2)復(fù)雜度降低

(3)局部極值點判決

由于雙線性問題的“非凸”特性，使得迭代算法可能會陷入局部極值點(鞍點)，給整體的估計精度帶來較大損失。鞍點在雙線性問題中容易識別，即判定收斂后的殘差r=‖Y-Φ{β}S‖2是否大于閾值，如果是則可判定為鞍點。此外為提升算法收斂速度，本文提出了一種冷/熱重啟的解決方案，實現(xiàn)方式是最初迭代時僅迭代較少次數(shù)，選擇初始化NOuter=10，每次迭代進(jìn)行收斂性判斷，當(dāng)識別為收斂即非鞍點時，采用熱重啟方式，將本次所估計的矩陣S作為入口參數(shù)進(jìn)入重啟迭代；反之識別為發(fā)散或鞍點，則采用冷重啟方式，即將矩陣S設(shè)置為隨機(jī)。

4 數(shù)值仿真和復(fù)雜度分析

本文選擇均勻線性陣列，陣元個數(shù)M=30，則虛擬網(wǎng)格個數(shù)N=30，即網(wǎng)格寬度為6°。設(shè)置快拍數(shù)為L=2～12，設(shè)置閾值γmax=103。為衡量本文所提算法的有效性，本節(jié)選擇文獻(xiàn)中已有的網(wǎng)格化稀疏貝葉斯學(xué)習(xí)(Off-Grid Sparse Bayesian Learning,OGSBL)、求根SBL(Root SBL,RTSBL)和高斯-賽德爾根(Gauss-Seidel Root,GSROOT)作為對比，此外還選擇克拉美羅界(Cramer-Rao Bound,CRB)作為參考。

4.1 性能分析

圖2給出了一次蒙特卡洛仿真中估計值和真實值的對比，仿真中選擇真實角度為[-28.6°,-18.6°,3.5°,15.6°,31.7°]。本仿真中設(shè)置信噪比為30 dB。從圖2中可以看出，本文所提算法能夠精確捕捉到5個來波信號的離格角度和準(zhǔn)確功率，在無來波信號的角度估計結(jié)果趨近于0。

圖2 真實角度/功率與估計值對比

圖3對比了BAd-VAMP與已有算法歸一化均方誤差(Normalized Mean Square Error,NMSE)隨快拍次數(shù)的變化曲線，數(shù)值仿真結(jié)果為500次蒙特卡洛仿真的平均值。每次仿真中波達(dá)方向以[-30°,-18°,6°,18°,30°]為基礎(chǔ)，疊加U(-3°,3°)的隨機(jī)角度偏移產(chǎn)生。圖3中設(shè)置信噪比SNR=20 dB，來波信號個數(shù)K=5。從仿真結(jié)果中可以看出，本文所提算法在快拍數(shù)較多時性能與現(xiàn)有的RTSBL和GSROOT較為接近。但是隨快拍數(shù)的減少，RTSBL和GSROOT方法性能明顯惡化，使得本文所提算法表現(xiàn)出較大的性能增益。從圖中還能看出，由于不考慮網(wǎng)格偏差，OGSBL算法性能較差，從而證明了離格估計算法的必要性。

圖3 估計性能與快拍個數(shù)關(guān)系圖(SNR=20 dB)

圖4對比了在不同快拍個數(shù)下估計性能隨信噪比的變化曲線，從圖中可以得出與圖3相同的結(jié)論，即在快拍個數(shù)較少時本文所提算法具有明顯的性能優(yōu)勢。要達(dá)到相同的估計性能，本文所提算法需要的快拍個數(shù)較少，更適合于機(jī)載雷達(dá)等快時變場景。

圖4 估計性能隨信噪比變化曲線

4.2 復(fù)雜度分析

本文所提算法可以劃分為S估計和β估計兩部分。S估計在式(9)需要矩陣乘法ΦT{β}Φ{β}，具有復(fù)雜度O(MN2)，其余運(yùn)算均為標(biāo)量形式，復(fù)雜度為O(MNL)。β估計部分復(fù)雜度集中于式(15)～(19)，其中式(15)中利用Ei矩陣的特征和矩陣求跡的特點，可知只需要計算i=j情況，所以式(15)具有復(fù)雜度O(MN)，式(17)需要計算L個矩陣求逆，具有極高的復(fù)雜度。為此,文獻(xiàn)[14]提出了一種奇異值分解(Singular Value Decomposition,SVD)的替代方法。通過SVD分解得

Udiag(η)V=SVD(ΦT{β}Φ{β})，

則式(17)可等效計算為

從而將式(17)復(fù)雜度降低為O(MN2)。式(19)需要矩陣求逆操作，具有復(fù)雜度O(N3)，所以本文所提算法整體復(fù)雜度可寫為O(MN2+MNL+N3)。根據(jù)3.2節(jié)所提方法，通過設(shè)置閾值可以隨迭代降低式(19)中的矩陣H維度，從而使得本文所提算法復(fù)雜度進(jìn)一步降低。參考文獻(xiàn)[10]中的復(fù)雜度分析，表1給出了各算法的復(fù)雜度對比，可以看出本文方法與文獻(xiàn)中現(xiàn)有方法保持同階復(fù)雜度。

表1 復(fù)雜度對比

5 結(jié) 論

本文將離格波達(dá)角估計問題建模為典型雙線性問題，通過最新的BAd-VAMP算法對其進(jìn)行估計。利用離格DOA估計模型的共同稀疏性等特點，本文還對BAd-VAMP算法進(jìn)行了改進(jìn)，降低了算法復(fù)雜度，并且引入迭代阻尼提升了魯棒性。仿真結(jié)果表明，相比于現(xiàn)有的RTSBL等方法，本文所提算法在機(jī)載雷達(dá)等快拍數(shù)較少時表現(xiàn)出較大的性能優(yōu)勢，具有較高的應(yīng)用價值。