基于稀疏貝葉斯的陣列幅相和互耦誤差聯合校正方法

王鼎，高衛港，吳志東

（信息工程大學信息系統工程學院，河南 鄭州 450001）

0 引言

陣列信號處理技術自面世以來，在雷達[1]、通信[2]、水聲[3]、醫學[4]等領域得到了廣泛的應用。同時，隨著陣列處理理論的不斷完善發展，空間譜估計技術因其在遠距離探測、低截獲信號處理等方面的優勢，受到了國內外學者的關注。

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

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

實際中，無論是具有超分辨率的子空間算法，還是稀疏重構算法，都無法達到理論估計精度，這是因為在實際中無法避免由于各種原因導致的陣列誤差，包括各個陣元和通道的幅度和相位不具有一致性[15]、天線之間相互影響產生互耦效應[16]、陣元安裝位置出現偏差引入的位置誤差[17]等，都會使陣列流型和理想的陣列流型之間存在較大的偏差，從而使各種算法的性能出現一定程度的惡化，甚至完全失效。針對幅相誤差和互耦誤差共存的情況，研究人員[18-20]提出了一系列的聯合校正方法，其中，文獻[18]利用均勻線陣互耦矩陣的特殊結構，通過交替迭代的方式對陣列誤差參數進行優化校正。文獻[19]利用陣列協方差矩陣的關系式構造了一種二次代價函數，并通過交替迭代進行參數估計。文獻[20]提出了基于協方差匹配技術的誤差聯合校正算法。文獻[21]利用MIMO 系統，引入若干個經過精確校準的輔助陣元對誤差及波達方向進行聯合估計。

基于對當前研究現狀的分析，本文針對幅相誤差和互耦誤差共存的情況，提出了基于信號空域稀疏性的貝葉斯測向與誤差校正方法，并分析了信號測向精度及相應誤差校準的理論下限，仿真結果驗證了方法的有效性。

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

1 信號模型

假設在空域中有K個遠場窄帶信號入射到M元陣列，入射方向為，那么M×1維的陣列接收數據為

圖1 陣列信號接收模型

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

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

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

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

2 模型求解

2.1 參數求解

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

陣列輸出的概率密度函數為

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

令式(11)的導數為零，可得c第q次迭代時的值為

令式(13)的導數為零，可得f第q次迭代時的值為

令式(15)的導數為零，可得第q次迭代時δ2的值為

令式(17)的導數為零，可得第q次迭代時的值為

其中，constant 是一個常數余項，P是一個半正定矩陣，計算式為

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

由式(19)可知

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

由于在t時刻，信號幅度僅與該時刻的陣列觀測值x(t) 有關，根據式(6)和式(7)可以得到關于x(t) 的后驗概率密度函數為

將式(26)～式(30)代入參數的更新式便可實現對參數c,δ2,γ,β的更新。而后估計出信號幅度的均值與方差，進入下一次的EM 迭代。由文獻[23]可知，EM 算法具有嚴格的收斂性，該迭代過程逐漸逼近穩定的估計值，當迭代收斂時，由各變量的估計值便可以獲得信號的空域分布、陣列誤差模型和噪聲功率等。

3 參數估計方差的克拉美羅界

N個陣列觀測樣本x1,…,xN的概率密度函數為

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

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

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

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

4 仿真實驗與結果分析

本節將通過模擬仿真實驗對本文提出的多種誤差同時存在時的校正方法進行驗證，下面對一些固定參數進行說明，陣列為8 元均勻直線陣，半徑波長比為0.5，信號入射角度為(30.3° 36.8°)，算法最大迭代次數為1 500 次，停止閾值為10-3。采用均方根誤差作為評判指標，其中，N為蒙特卡羅仿真次數，在本文中設置為為第i次迭代的估計值，?為真實值。互耦系數為c=(0.6 +j0.4,0.2 +j0.1)T。值得一提的是，在這里假設第一個陣元完全校準，因此，幅度誤差分別為(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 校正前后的空間譜圖

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

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

從圖2 中可以看出，在多種誤差同時存在的情況下，空間譜的峰值被削弱，譜峰也發生相應的偏移，導致分辨率下降。根據文獻[25]，在統計意義下，幅相誤差對靠近來波方向的空間譜峰影響較大，其中，幅度誤差只影響空間譜的峰值尖銳程度，不會改變譜峰的位置，相位誤差則會導致譜峰發生偏移，而互耦誤差則對兩者皆有影響。利用本文方法對幅相誤差和互耦誤差進行聯合估計，可以很好地對來波信號方向進行估計，提高空間譜的分辨能力。

4.2 性能分析

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

圖3 信號方位估計均方根誤差隨信噪比變化曲線

圖4 信號方位估計隨信噪比變化箱線

圖5 信號方位估計均方根誤差隨快拍數變化曲線

圖6 信號方位估計隨快拍數變化箱線

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

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

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

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

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

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

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

表1 互耦系數估計結果

表2 幅相系數估計結果

5 結束語

本文提出了幅相誤差和互耦誤差同時存在時基于稀疏貝葉斯算法的求解方法。該方法共包含2 個階段，第一階段，在空域稀疏的前提下構建了誤差存在時接收信號的超完備模型，得到接收信號的后驗概率密度函數。由于密度函數具有極強的非線性，直接優化十分困難，第二階段采用EM 算法對該函數進行迭代優化求解相應參數。此外，本文還推導了相應參數的克拉美羅界。仿真實驗結果表明，本文方法測向性能優于OGSBI、rootSB 和L1SVD 中的方法，同時能夠對幅相誤差和互耦誤差進行有效的校正。