王鼎，高衛(wèi)港，吳志東

（信息工程大學(xué)信息系統(tǒng)工程學(xué)院，河南 鄭州 450001）

0 引言

陣列信號(hào)處理技術(shù)自面世以來(lái)，在雷達(dá)[1]、通信[2]、水聲[3]、醫(yī)學(xué)[4]等領(lǐng)域得到了廣泛的應(yīng)用。同時(shí)，隨著陣列處理理論的不斷完善發(fā)展，空間譜估計(jì)技術(shù)因其在遠(yuǎn)距離探測(cè)、低截獲信號(hào)處理等方面的優(yōu)勢(shì)，受到了國(guó)內(nèi)外學(xué)者的關(guān)注。

縱觀空間譜估計(jì)技術(shù)的整個(gè)發(fā)展歷程，從最開(kāi)始的傳統(tǒng)空間譜估計(jì)，即經(jīng)典傅式譜分析在空域上的直接推廣，這種方法因受限于陣列孔徑無(wú)法突破“瑞利限”[5]，到后來(lái)的現(xiàn)代空間譜估計(jì)技術(shù)，其中最具代表性的便是以子空間技術(shù)為基礎(chǔ)的多重信號(hào)分類（MUSIC,multi signal classification）算法[6]。與傳統(tǒng)空間譜估計(jì)技術(shù)相比，MUSIC 算法的估計(jì)精度和角度分辨率都達(dá)到了新的高度。但是，隨著子空間算法理論的不斷發(fā)展，此類算法的局限性也逐漸暴露出來(lái)，對(duì)陣列輸出協(xié)方差矩陣以及信號(hào)/噪聲子空間的準(zhǔn)確度有高度依賴性導(dǎo)致子空間算法在低信噪比及小快拍情況下效果并不理想。此外，在發(fā)展過(guò)程中，還出現(xiàn)了極大似然類測(cè)向方法[7]，這種方法通過(guò)借助特定維數(shù)模型與觀測(cè)數(shù)據(jù)之間的最佳擬合來(lái)實(shí)現(xiàn)空間譜估計(jì)，在一定程度上提高了對(duì)低信噪比和小樣本情況的適應(yīng)能力，但其又因?yàn)檩^低的計(jì)算效率及對(duì)信號(hào)個(gè)數(shù)先驗(yàn)的依賴無(wú)法很好地發(fā)揮性能優(yōu)勢(shì)。

無(wú)論是低信噪比還是小樣本情況，陣列信號(hào)處理問(wèn)題本質(zhì)上都屬于從有限的觀測(cè)數(shù)據(jù)中獲得相應(yīng)的參數(shù)，而近幾年興起的稀疏重構(gòu)技術(shù)在解決信號(hào)缺失條件下信號(hào)恢復(fù)問(wèn)題的過(guò)程中展現(xiàn)出了極大優(yōu)勢(shì)，目前在壓縮感知[8]、信號(hào)分析[9]、圖像處理[10]等領(lǐng)域已有了一定規(guī)模的應(yīng)用，并取得了不錯(cuò)的成果。稀疏重構(gòu)是指利用信號(hào)在特定變換域（如空域）上的稀疏性，從包含信號(hào)所有可能分量的完備或超完備集合中選取少量幾個(gè)基函數(shù)重構(gòu)原始信號(hào)的過(guò)程。目前，稀疏重構(gòu)算法主要包括匹配追蹤[11]、Lp范數(shù)[12]和稀疏貝葉斯[13]三類。其中，匹配追蹤類算法對(duì)觀測(cè)模型的稀疏度有很高的要求，在信號(hào)不同分量之間相關(guān)性較強(qiáng)時(shí)性能不好；Lp范數(shù)類算法的局限性則在于正則化因子參數(shù)的選取問(wèn)題，具體的觀測(cè)模型、信號(hào)環(huán)境等因素都會(huì)對(duì)該參數(shù)的選取產(chǎn)生極大的影響；稀疏貝葉斯學(xué)習(xí)方法是由Tipping[14]提出的利用貝葉斯原理與觀測(cè)模型先驗(yàn)信息相結(jié)合的一類重構(gòu)方法。目前，稀疏貝葉斯重構(gòu)算法在圖像處理領(lǐng)域的發(fā)展較為完善，在其他領(lǐng)域則相對(duì)滯后。

實(shí)際中，無(wú)論是具有超分辨率的子空間算法，還是稀疏重構(gòu)算法，都無(wú)法達(dá)到理論估計(jì)精度，這是因?yàn)樵趯?shí)際中無(wú)法避免由于各種原因?qū)е碌年嚵姓`差，包括各個(gè)陣元和通道的幅度和相位不具有一致性[15]、天線之間相互影響產(chǎn)生互耦效應(yīng)[16]、陣元安裝位置出現(xiàn)偏差引入的位置誤差[17]等，都會(huì)使陣列流型和理想的陣列流型之間存在較大的偏差，從而使各種算法的性能出現(xiàn)一定程度的惡化，甚至完全失效。針對(duì)幅相誤差和互耦誤差共存的情況，研究人員[18-20]提出了一系列的聯(lián)合校正方法，其中，文獻(xiàn)[18]利用均勻線陣互耦矩陣的特殊結(jié)構(gòu)，通過(guò)交替迭代的方式對(duì)陣列誤差參數(shù)進(jìn)行優(yōu)化校正。文獻(xiàn)[19]利用陣列協(xié)方差矩陣的關(guān)系式構(gòu)造了一種二次代價(jià)函數(shù)，并通過(guò)交替迭代進(jìn)行參數(shù)估計(jì)。文獻(xiàn)[20]提出了基于協(xié)方差匹配技術(shù)的誤差聯(lián)合校正算法。文獻(xiàn)[21]利用MIMO 系統(tǒng)，引入若干個(gè)經(jīng)過(guò)精確校準(zhǔn)的輔助陣元對(duì)誤差及波達(dá)方向進(jìn)行聯(lián)合估計(jì)。

基于對(duì)當(dāng)前研究現(xiàn)狀的分析，本文針對(duì)幅相誤差和互耦誤差共存的情況，提出了基于信號(hào)空域稀疏性的貝葉斯測(cè)向與誤差校正方法，并分析了信號(hào)測(cè)向精度及相應(yīng)誤差校準(zhǔn)的理論下限，仿真結(jié)果驗(yàn)證了方法的有效性。

本文用到的一些數(shù)字符號(hào)說(shuō)明如下。向量和矩陣均用粗斜體表示。xT和xH分別表示向量的轉(zhuǎn)置和共軛轉(zhuǎn)置分別表示L1 范數(shù)、L2 范數(shù)、Frobenous 范數(shù)。表示矩陣A的行列式。tr(A) 表示矩陣A的跡。(A)j,:表示矩陣A的第j行，(A):,j表示矩陣A的第j列，(A)i,j表示矩陣A的第(i,j) 個(gè)元素。xj表示向量x的第j個(gè)元素。diag(A) 表示由A的對(duì)角線元素構(gòu)成的向量；diag(x) 表示對(duì)角矩陣，其對(duì)角向量對(duì)應(yīng)于向量x；toeplitz(x) 表示由x形成的托普利茲矩陣。x'(?)表示x(?)關(guān)于?的導(dǎo)數(shù)。Re(·) 和Im(·) 分別表示對(duì)復(fù)數(shù)進(jìn)行取實(shí)和取虛操作。x⊙y表示x和y的Hadamard 積。表示給定分布概率條件下對(duì)應(yīng)的期望運(yùn)算符。

1 信號(hào)模型

假設(shè)在空域中有K個(gè)遠(yuǎn)場(chǎng)窄帶信號(hào)入射到M元陣列，入射方向?yàn)?，那么M×1維的陣列接收數(shù)據(jù)為

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

在實(shí)際中，式(1)所對(duì)應(yīng)的理想情況是不存在的，更普遍的是幅相誤差與互耦誤差同時(shí)存在，在這種情況下，陣列輸出的計(jì)算式為

采用文獻(xiàn)[22]中使用的離格模型，式(2)可以寫為

為了進(jìn)一步明確陣列輸出與互耦誤差系數(shù)和幅相誤差參數(shù)之間的關(guān)系，對(duì)式(3)做如下變換

通過(guò)對(duì)空間角度域進(jìn)行均勻采樣，將式(3)和式(4)進(jìn)行空域擴(kuò)展，得到陣列輸出的超完備表示形式為

2 模型求解

2.1 參數(shù)求解

引入中間變量γ=[γ1,…,γ L]T表示方 向?1,…,? L處對(duì)應(yīng)的信號(hào)分量幅度的方差，即其中R=diag(γ)。

陣列輸出的概率密度函數(shù)為

結(jié)合式(6)和式(7)，得到給定γ,c,δ2,β條件下的概率為

令式(11)的導(dǎo)數(shù)為零，可得c第q次迭代時(shí)的值為

令式(13)的導(dǎo)數(shù)為零，可得f第q次迭代時(shí)的值為

令式(15)的導(dǎo)數(shù)為零，可得第q次迭代時(shí)δ2的值為

令式(17)的導(dǎo)數(shù)為零，可得第q次迭代時(shí)的值為

其中，constant 是一個(gè)常數(shù)余項(xiàng)，P是一個(gè)半正定矩陣，計(jì)算式為

其中，μt由式(25)得到，實(shí)際上，β與是聯(lián)合稀疏的，對(duì)應(yīng)的K個(gè)源位置處有K個(gè)非零值，只需計(jì)算對(duì)應(yīng)γl的最大K個(gè)位置處的β值，其余值設(shè)為零。因此，β、ν和P的大小可以分別減小為K、K和K×K維。

由式(19)可知

對(duì)于式(12)、式(13)、式(16)和式(18)中給出的期望值取決于條件概率密度函數(shù)并且的形式不是直觀給出的，因此有必要對(duì)上述參數(shù)更新方法進(jìn)行簡(jiǎn)化。

由于在t時(shí)刻，信號(hào)幅度僅與該時(shí)刻的陣列觀測(cè)值x(t) 有關(guān)，根據(jù)式(6)和式(7)可以得到關(guān)于x(t) 的后驗(yàn)概率密度函數(shù)為

將式(26)～式(30)代入?yún)?shù)的更新式便可實(shí)現(xiàn)對(duì)參數(shù)c,δ2,γ,β的更新。而后估計(jì)出信號(hào)幅度的均值與方差，進(jìn)入下一次的EM 迭代。由文獻(xiàn)[23]可知，EM 算法具有嚴(yán)格的收斂性，該迭代過(guò)程逐漸逼近穩(wěn)定的估計(jì)值，當(dāng)?shù)諗繒r(shí)，由各變量的估計(jì)值便可以獲得信號(hào)的空域分布、陣列誤差模型和噪聲功率等。

3 參數(shù)估計(jì)方差的克拉美羅界

N個(gè)陣列觀測(cè)樣本x1,…,xN的概率密度函數(shù)為

對(duì)概率密度函數(shù)（式(31)）取對(duì)數(shù)可以得到如下對(duì)數(shù)似然函數(shù)

其中，?和δ2是實(shí)數(shù)域上的參數(shù)，s(t) 和c是復(fù)向量。因此，將s(t) 和c分解為實(shí)部和虛部之和的形式，分別為其中，分別代表復(fù)向量的實(shí)部與虛部。經(jīng)過(guò)上述分解之后，得到的參數(shù)集為

根據(jù)式(32)的似然函數(shù)，可以得到關(guān)于噪聲方差δ2的克拉美-羅下界（CRLB,Cramer-Rao lower bound）為

在存在陣列互耦或者陣列幅相誤差時(shí)，式(35)給出的CRLB 只是簡(jiǎn)單地將誤差向量分解為實(shí)部與虛部，這并不能很好地反映出相關(guān)參數(shù)的物理屬性，因此，結(jié)合文獻(xiàn)[24]中的結(jié)論，可以將的CRLB 轉(zhuǎn)化為的CRLB。這樣就可以更加清晰地得出復(fù)向量c的估計(jì)精度。

4 仿真實(shí)驗(yàn)與結(jié)果分析

本節(jié)將通過(guò)模擬仿真實(shí)驗(yàn)對(duì)本文提出的多種誤差同時(shí)存在時(shí)的校正方法進(jìn)行驗(yàn)證，下面對(duì)一些固定參數(shù)進(jìn)行說(shuō)明，陣列為8 元均勻直線陣，半徑波長(zhǎng)比為0.5，信號(hào)入射角度為(30.3° 36.8°)，算法最大迭代次數(shù)為1 500 次，停止閾值為10-3。采用均方根誤差作為評(píng)判指標(biāo)，其中，N為蒙特卡羅仿真次數(shù)，在本文中設(shè)置為為第i次迭代的估計(jì)值，?為真實(shí)值?；ヱ钕禂?shù)為c=(0.6 +j0.4,0.2 +j0.1)T。值得一提的是，在這里假設(shè)第一個(gè)陣元完全校準(zhǔn)，因此，幅度誤差分別為(0 0.15 0.30 0.20 0.25 -0.20 -0.20 -0.30)，相位誤差分別為(0° -20°-30° -40° 35° 25° 20° 40°) 。

4.1 校正前后的空間譜圖

本節(jié)驗(yàn)證幅相誤差和互耦誤差對(duì)信號(hào)空間譜的影響，并探究本文方法對(duì)信號(hào)的分辨能力。假設(shè)2 個(gè)不相關(guān)信號(hào)以(30.3° 36.8°)的方向入射至8 元均勻直線陣，信噪比為10 dB，快拍數(shù)為100，其余參數(shù)不變。圖2 給出利用本文方法誤差校正前后的空間譜。其中，虛線輔助線為信號(hào)真實(shí)方向。

圖2 誤差校正前后的空間譜

從圖2 中可以看出，在多種誤差同時(shí)存在的情況下，空間譜的峰值被削弱，譜峰也發(fā)生相應(yīng)的偏移，導(dǎo)致分辨率下降。根據(jù)文獻(xiàn)[25]，在統(tǒng)計(jì)意義下，幅相誤差對(duì)靠近來(lái)波方向的空間譜峰影響較大，其中，幅度誤差只影響空間譜的峰值尖銳程度，不會(huì)改變譜峰的位置，相位誤差則會(huì)導(dǎo)致譜峰發(fā)生偏移，而互耦誤差則對(duì)兩者皆有影響。利用本文方法對(duì)幅相誤差和互耦誤差進(jìn)行聯(lián)合估計(jì)，可以很好地對(duì)來(lái)波信號(hào)方向進(jìn)行估計(jì)，提高空間譜的分辨能力。

4.2 性能分析

本節(jié)比較的方法包括離格稀疏貝葉斯推斷（OGSBI,off-grid sparse Bayesian inference）、求根稀疏貝葉斯學(xué)習(xí)（rootSBL,root sparse Bayesian learning）[26]和L1SVD 算法，分別從信噪比和快拍數(shù)2 個(gè)角度對(duì)本文方法的性能進(jìn)行探究，結(jié)果分別如圖3～圖6 所示。其中，當(dāng)探究信噪比對(duì)算法的影響時(shí)，信噪比的范圍為-10～10 dB，快拍數(shù)固定為100；當(dāng)探究快拍數(shù)對(duì)算法的影響時(shí)，快拍數(shù)范圍為10～100，信噪比固定為0，其余參數(shù)不變。

圖3 信號(hào)方位估計(jì)均方根誤差隨信噪比變化曲線

圖4 信號(hào)方位估計(jì)隨信噪比變化箱線

圖5 信號(hào)方位估計(jì)均方根誤差隨快拍數(shù)變化曲線

圖6 信號(hào)方位估計(jì)隨快拍數(shù)變化箱線

從圖3 和圖5 中可以看出，當(dāng)分別采用信噪比和快拍數(shù)為變量時(shí)，本文方法均能達(dá)到很好的效果，且逼近克拉美羅界。在信噪比較小和快拍數(shù)較低的情況下，稀疏貝葉斯類算法中的rootSBL 和OGSBI 與本文方法表現(xiàn)大致相同，但隨著信噪比和快拍數(shù)的增加，rootSBL 和OGSBI 算法雖然考慮了離散網(wǎng)格對(duì)信號(hào)方位估計(jì)造成的影響，但互耦誤差的存在使估計(jì)的空間譜偽峰較高，一定程度上限制了這2 種算法的性能。L1SVD 算法未考慮網(wǎng)格誤差的影響，效果最差。

從圖4 和圖6 中可以看出，隨著信噪比和快拍數(shù)的提升，箱線圖中的中位線逐漸逼近真實(shí)的來(lái)波方向，整體的波動(dòng)程度迅速下降。2 次實(shí)驗(yàn)的蒙特卡羅次數(shù)均為100 次，直觀來(lái)看，當(dāng)信噪比提升時(shí)，箱線圖的波動(dòng)程度相比于快拍數(shù)增大時(shí)收斂得更快，箱體更小。

將快拍數(shù)設(shè)置為 10～100，信噪比設(shè)置為-10～10 dB，其余條件保持不變，繪制互耦誤差和幅相誤差隨信噪比快拍數(shù)的變化曲面，分別如圖7 和圖8 所示。

圖7 互耦誤差隨信噪比快拍數(shù)變化曲面

圖8 幅相誤差隨信噪比快拍變化曲面

從圖7 和圖8 中可以看出，隨著信噪比和快拍數(shù)的提高，互耦誤差和幅相誤差曲面均能迅速下降，其中，互耦誤差曲面比幅相誤差曲面整體偏低的原因有2 個(gè)：1) 在觀測(cè)模型中，幅相誤差與信號(hào)方向之間的相互耦合效果要比互耦誤差與信號(hào)方向之間的耦合效果強(qiáng)；2) 采用實(shí)驗(yàn)值與真實(shí)值之間差值的二范數(shù)作為評(píng)價(jià)指標(biāo)，本文中互耦誤差系數(shù)只有2 個(gè)，幅相誤差系數(shù)則有8 個(gè)。

表1 和表2 是在信噪比為10 dB、快拍數(shù)為100、其余條件與上文保持一致條件下的互耦系數(shù)和幅相系數(shù)估計(jì)結(jié)果。

表1 互耦系數(shù)估計(jì)結(jié)果

表2 幅相系數(shù)估計(jì)結(jié)果

5 結(jié)束語(yǔ)

本文提出了幅相誤差和互耦誤差同時(shí)存在時(shí)基于稀疏貝葉斯算法的求解方法。該方法共包含2 個(gè)階段，第一階段，在空域稀疏的前提下構(gòu)建了誤差存在時(shí)接收信號(hào)的超完備模型，得到接收信號(hào)的后驗(yàn)概率密度函數(shù)。由于密度函數(shù)具有極強(qiáng)的非線性，直接優(yōu)化十分困難，第二階段采用EM 算法對(duì)該函數(shù)進(jìn)行迭代優(yōu)化求解相應(yīng)參數(shù)。此外，本文還推導(dǎo)了相應(yīng)參數(shù)的克拉美羅界。仿真實(shí)驗(yàn)結(jié)果表明，本文方法測(cè)向性能優(yōu)于OGSBI、rootSB 和L1SVD 中的方法，同時(shí)能夠?qū)Ψ嗾`差和互耦誤差進(jìn)行有效的校正。