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

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

0 引 言

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

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

近年來國內(nèi)外多個(gè)團(tuán)隊(duì)提出了基于離格信號(hào)的稀疏重構(gòu)算法。文獻(xiàn)[9]提出了一種基于稀疏最小二乘的離格算法，假定真實(shí)角度不必限于網(wǎng)格之上，能夠計(jì)算出誤差并對估計(jì)值進(jìn)行修正。文獻(xiàn)[10]則提出了一種求根稀疏貝葉斯學(xué)習(xí)方法，在線性陣列場景下能夠取得很好的估計(jì)性能。由于復(fù)雜度較低，經(jīng)典的近似消息傳遞算法(Approximate Message Passing,AMP)在稀疏貝葉斯學(xué)習(xí)框架下得到了廣泛應(yīng)用[11]，但是普通的AMP算法在遇到較“壞”的觀測矩陣時(shí)(如病態(tài)條件數(shù)、列相關(guān)、低秩等)會(huì)導(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ǔ)全等多個(gè)場景下得到了應(yīng)用。

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

本文符號(hào)說明：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個(gè)快拍下M個(gè)陣元的接收向量；K×1表示第t個(gè)快拍下的K個(gè)來波信號(hào)；w表示均值為0、方差為1/λ的高斯白噪聲，其中λ表示噪聲精度。由于噪聲方差不隨時(shí)間變化，本文中省略w的上標(biāo)t。陣列流形表示為

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

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

(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估計(jì)模型中有L個(gè)快拍時(shí)，上式可以寫為多重觀測形式：

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

(4)

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

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

(5)

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

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

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

Y=Φ{β}C+W。

(6)

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

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

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

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

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估計(jì)算法

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

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

定義

(9)

定義

(10)

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

其中均值方差分別為

(11)

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

(12)

其中期望和方差分別計(jì)算為

(13)

從而完成S估計(jì)部分的迭代計(jì)算。

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

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

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

(14)

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

(15)

(16)

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

(17)

(18)

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

βmax=H-1α。

3.3 算法改進(jìn)

(1)迭代阻尼

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

(3)局部極值點(diǎn)判決

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

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

本文選擇均勻線性陣列，陣元個(gè)數(shù)M=30，則虛擬網(wǎng)格個(gè)數(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給出了一次蒙特卡洛仿真中估計(jì)值和真實(shí)值的對比，仿真中選擇真實(shí)角度為[-28.6°,-18.6°,3.5°,15.6°,31.7°]。本仿真中設(shè)置信噪比為30 dB。從圖2中可以看出，本文所提算法能夠精確捕捉到5個(gè)來波信號(hào)的離格角度和準(zhǔn)確功率，在無來波信號(hào)的角度估計(jì)結(jié)果趨近于0。

圖2 真實(shí)角度/功率與估計(jì)值對比

圖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，來波信號(hào)個(gè)數(shù)K=5。從仿真結(jié)果中可以看出，本文所提算法在快拍數(shù)較多時(shí)性能與現(xiàn)有的RTSBL和GSROOT較為接近。但是隨快拍數(shù)的減少，RTSBL和GSROOT方法性能明顯惡化，使得本文所提算法表現(xiàn)出較大的性能增益。從圖中還能看出，由于不考慮網(wǎng)格偏差，OGSBL算法性能較差，從而證明了離格估計(jì)算法的必要性。

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

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

圖4 估計(jì)性能隨信噪比變化曲線

4.2 復(fù)雜度分析

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

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

則式(17)可等效計(jì)算為

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