陣列互耦情況下基于稀疏貝葉斯學習的離網格DOA估計

王緒虎，白浩東，張群飛，田 雨

(1.青島理工大學 信息與控制工程學院，青島 266520；2.西北工業大學 航海學院，西安 710072)

波達方向(direction of arrival,DOA)估計問題(在很多領域受到極大的關注，例如，雷達[1]、聲吶[2-4]，生物醫學[5]等。在過去幾十年已經提出了許多高分辨的DOA估計算法，例如，多重信號分類算法[6](multiple signal classification,MUSIC)、求根MUSIC算法[7](Root-MUSIC)、旋轉不變子空間算法[8](estimating signal parameters via rotational invariance techniques,ESPRIT)等，但上述高分辨算法在低信噪比、少快拍數的情況下，估計性能會嚴重下降，甚至失效。

隨著壓縮感知(compress sensing,CS)理論和稀疏重構技術的不斷發展和成熟，許多學者將其與DOA估計聯系起來，使DOA估計技術進入了新的發展階段。相較于傳統的高分辨算法，基于CS和稀疏重構技術的DOA估計方法，在低信噪比、少快拍的情況下表現出良好的估計性能。該種方法大致分為三類：凸優化方法、貪婪方法[9]、稀疏貝葉斯學習方法[10-11]。例如，文獻[12]研究稀疏信號的表示以及DOA估計問題，提出L1范數奇異值分解(L1norm-singular value decomposition,L1-SVD)方法，把測向問題轉化為求解L1范數的問題，但是其優化問題受到模型正則化參數的影響，影響估計的精度。王偉東等提出的基于稀疏信號功率迭代補償的DOA估計算法，該算法基于補償原理，使信號的稀疏表示近似于L0范數，把DOA估計問題轉化為求解L0范數問題。相比于L1-SVD算法，該算法不需要設置正則化參數，對非相干源具有較高的估計精度。然而，凸優化方法的計算效率限制了其進一步發展。

隨著CS理論研究的不斷深入，稀疏貝葉斯學習方法(sparse Bayesian learning,SBL)被認為與凸優化方法具有相同的全局收斂性，而它計算效率又遠遠優于后者。Ji等將SBL引入CS領域中，提出貝葉斯壓縮感知算法。該算法通常需要滿足聲源的來波方向位于網格點上，才能實現較高的方位估計精度。當入射聲源的波達方向偏離預先劃分好的網格，就導致方位估計精度的降低。因此，Yang等[13]提出了離網格的稀疏貝葉斯學習算法 (off-grid sparse Bayesian learning,OGSBL)，在離網格模型中，在信號真實到達角處采用一階泰勒展開式近似表示，使得估計性能進一步提高。文獻[14]中，提出了求根離格稀疏貝葉斯學習算法(root-off grid sparse Bayesian learning,Root-OGSBL)，降低了OGSBL方法的計算復雜性，在網格間距較小的情況下，提高了計算效率的同時保證了估計的精度。文獻[15]中，使用平均場理論中的變分理論來估計參數的后驗分布，并且提出了三層信號先驗分布，這種分層先驗進一步提升了信號的稀疏性，提升了方位估計的精度。文獻[16]中，提出了稀疏信號新的先驗分布框架—分層合成Lasso先驗，相比于假設伽馬先驗分布，具有更高的稀疏性、更低的重建誤差，提升了方位估計精度。

在實際的聲吶測向系統中，不可避免的存在陣列誤差。例如，水聽器陣元間的耦合、水聽器位置的偏差、水聽器陣元的幅相通道不一致。導致現有常規算法的方位估計精度下降，甚至失效。在文獻[17-18]中，提出了陣元間存在未知互耦下波達方向估計算法，解決了陣元間存在未知互耦情況下，稀疏信號的恢復問題。針對稀疏模型中入射信號方向與離散網格存在誤差，且陣元間的未知互耦沒有被同時考慮的問題，本文提出了一種水聽器陣元間存在不確定互耦，基于稀疏貝葉斯學習的離網格DOA估計方法。

1 信號模型

考慮K個不同方向的遠場窄帶信號入射到如圖1所示的均勻線列陣，該陣列包含M個全向陣元，陣元之間存在互耦，陣元間距為d，因此，該接收模型可以表示為

圖1 M個水聽器組成的均勻線列陣Fig.1 Uniform linear array composed of M-element hydrophones

n(t)

(1)

A=[a(θ1),a(θ2),…,a(θK)]

(2)

式中，a(θk)=[a1(θk),…,am(θk),…,aM(θk)]T為第k個入射信號的導向矢量，am(θk)=ej2π(m-1)d/λsin θk，λ為信號的波長。

式(1)為單快拍觀測矢量模型，當水聽器陣元接收多快拍數據時，多快拍觀測矢量模型可以表示為

Y=CAS+N

(3)

式中，Y=[y(1),y(2),…,y(T)]∈CM×T為多快拍的接收數據矩陣，S=[s(1),s(2),…,s(T)]∈CK×T為多快拍的入射信號數據矩陣，N=[n(1),n(2),…,n(T)]∈CM×T為多快拍的陣列接收噪聲矩陣，T為快拍數。

在本文中，我們研究水聽器陣元間存在未知耦合情況下的DOA估計問題。由文獻[19]可知，水聽器陣元間的互耦可以用矩陣C表示

(4)

該矩陣為對稱的托普利茲矩陣。根據文獻[20]的引理三，我們將上述模型進行處理。首先，定義向量c=[1,c1,c2,…,cM-1]T，且滿足C=toeplitz(c)，toeplitz(c)為使用向量c生成的循環對稱矩陣 ，在t時刻的接收信號可以寫為

yt=CAst+nt=Q(st?c)+nt

(5)

則對應水聽器陣列的多快拍數據可表示為

Y=Q(S?c)+N

(6)

2 稀疏貝葉斯學習模型

2.1 稀疏信號模型

(7)

式(7)中的系統模型就可以轉化為過完備的稀疏模型

(8)

式中，矩陣X是矩陣S的稀疏擴展

X=[x1,x2,…,xT]∈CN×T

(9)

圖2 信號X結構稀疏化模型Fig.2 The structure of sparse matrix X

(10)

在實際的陣列測向系統中，接收信號并不是落在預設的網格點上，此時，采用離網格信號模型，使用離接收信號最近的網格矢量一階泰勒展開式近似表示

(11)

圖3 離格陣列空域模型Fig.3 Off-grid array spatial domain model

因此，由式(8)和(11)可以表示陣元間存在互耦下的離網格模型

Y≈Ω(δ)(X?c)+N

(12)

最終，水聽器陣元間存在未知耦合的離網格稀疏模型建立如式(12)，接下來通過重建稀疏矩陣X來估計入射信號的波達方向，重建稀疏矩陣X中的非零行代表入射信號的實際波達方向。

2.2 模型參數分布假設

2.2.1 噪聲及其精度分布

假設接收的噪聲為獨立同分布的復對稱高斯白噪聲，其分布為

(13)

p(αn)=Γ(αn|a,b)

(14)

稀疏模型的似然函數可表示為

(15)

2.2.2 信號及其精度分布

假設稀疏信號X的多快拍樣本相互獨立，各列服從均值為0、方差為，

(16)

(17)

其中，c和d為伽馬分布的超參數。

2.2.3 互耦矢量及其精度分布

(18)

(19)

式中，e和f為伽馬分布的超參數。

2.2.4 離網格矢量分布

p(δ)=U(-r/2,r/2)

(20)

通過組合式(14)～(20)，可得模型的聯合概率密度

p(X,Y,δ,c,αx,αc,αn)=

p(Y|X,δ,c,αn)p(X|αx)p(c|αc)p(αn)

p(αx)p(αc)p(δ)

(21)

3 稀疏貝葉斯參數學習

利用上述的分布假設，我們提出了一種水聽器陣元間存在互耦，基于稀疏貝葉斯學習的離網格波達方向估計算法。該算法中，我們使用EM算法學習更新離網格誤差矢量和互耦矢量，直至收斂。

首先，利用上述接收的數據Y和參數先驗分布，可以得出信號X的后驗分布

p(X|Y,δ,c,αx,αc,αn)=

αn)p(X|αx)

(22)

式中，p(X|αx)的概率密度函數為

(23)

X后驗分布也為高斯分布，滿足

(24)

式中，μt為均值向量，Σx為協方差矩陣

(25)

Σx=[αnΥH(δ,c)Υ(δ,c)+diag(αx)]-1

(26)

其中，Υ(δ,c)=Ω(δ)(IN?c)，另外，定義所有快拍下均值矢量組成的矩陣

μ[μ1,…,μt,…,μT]

(27)

其次，我們對超參數進行學習，超參數的學習過程等價于最大化所有待估計參數的后驗概率分布

(28)

最大化待估計參數的后驗概率分布又等價于最大化似然函數

L(δ,c,αx,αc,αn)=

ε{lnp(Y|X,δ,c,αn)p(X|αx)p(c|αc)

p(αn)p(αx)p(αc)p(δ)}

(29)

式中，ε{·}表示條件后驗期望EX|Y,δ,c,αx,αc,αn{·}。接下來，詳細推導超參數迭代更新的表達式，下述推導中涉及到矩陣對向量變量的求導，參考附錄B。

忽略式(29)中與互耦矢量無關的概率項，可得到關于互耦矢量的似然函數

L(c)=ε{lnp(Y|X,δ,c,αn)p(c|αc)}=

(30)

因為EM算法是不斷收斂遞增的，所以估計參數就等價于求似然函數的極值，有如下關系

(31)

因此，對式(31)求導，可得互耦矢量的迭代表達式(32)

(32)

令式(32)等于零，可得

(33)

式中，yt表示陣列接收數據Y的第t列，μt表示稀疏信號的均值矩陣μ的第t列。

同理，忽略與信號精度無關的概率分布，可得到如下似然函數表達式

L(αx)=ε{lnp(X|αx)p(αx)}

(34)

對似然函數求導，令其等于零，可得信號精度的迭代表達式

(35)

Σx(n,n)表示協方差矩陣的第n行n列的元素，μn,t表示稀疏信號的均值矩陣μ的第n行t列的元素。

同理，對于噪聲精度，我們推導其似然函數和迭代表達式，如下

L(αn)=ε{lnp(Y|X,δ,c,αn)p(αn)}

(36)

(37)

對于互耦精度，給出似然函數和迭代表達式

L(αc)=ε{p(c|αc)p(αc)}

(38)

(39)

對于離網格間距，給出似然函數和迭代表達式

L(δ)=ε{lnp(Y|X,δ,c,αn)p(δ)}

(40)

(41)

其中，P∈RN×N，v∈RN×1下式(42)(43)給出

(42)

(43)

最后，基于上述分析算法總結如下：

第一步：輸入水聽器陣元接收數據Y；

第二步：初始化需要更新的參數X,δ,c,αx,αc,αn=mean(var(Y))，超參數a,b,c,d,e,f的設置[15]；

第三步：根據(25)(26)更新后驗期望矢量和協方差矩陣；

第四步：根據(35)更新信號精度αx；

第五步：根據(37)更新噪聲精度αn；

第六步：根據(33)和(39)更新互耦矢量及其精度；

第七步：根據(41)更新離網格間隔矢量；

4 仿真分析

4.1 空間譜分析

根據文獻[21]，水聽器間的耦合效應可以表示為

(44)

其中，互耦系數中的幅度ξ服從均勻分布，即ξ～U[-0.05,0.05]，互耦系數的相位φ也服從均勻分布，即φ～U[0,2π]，參數αc為耦合系數，表征相鄰陣元之間的耦合效應。理論上，相鄰陣元間的耦合效應是相同的，仿真參數的設置如表1所示。

表1 仿真參數Tab.1 Simulation parameters

陣元間的耦合系數αc為-15 dB，正則化參數為2時，四種方法歸一化空間譜如圖4所示；陣元間的耦合系數αc為-15 dB，正則化參數為3時，四種方法歸一化空間譜如圖5所示；陣元間的耦合系數αc為-5 dB，正則化參數為3時，四種方法歸一化空間譜如圖6所示。

從圖4至圖6可以看出MUSIC方法的空間譜主瓣寬度較寬，旁瓣較高，當增大耦合系數時，兩個主瓣的峰值和真實的角度值相差較大；對比圖4和圖5，可以看出當L1-SVD方法正則化參數為2時，主瓣寬度較窄，旁瓣較低；改變正則化參數為3時，主瓣寬度變寬，且出現了許多旁瓣；對比圖5和圖6可以看出，當耦合系數增大時，OGSBL方法旁瓣數量增加，且旁瓣譜級升高；對比圖4～圖6，本文提出的方法主瓣寬度和旁瓣譜級基本不變。

圖4 空間譜(正則化參數為2，αc=-15 dB)Fig.4 Spatial spectrum(Regularization parameters 2, αc=-15 dB)

圖5 空間譜(正則化參數為3，αc=-15 dB)Fig.5 Spatial spectrum(Regularization parameters 3， αc=-15 dB)

圖6 空間譜(正則化參數為3，αc=-5 dB)Fig.6 Spatial spectrum(Regularization parameters 3， αc=-5 dB)

四種方法的蒙特卡洛仿真結果如圖7所示，其中L1-SVD方法的正則化參數設置為3。圖7(a)是耦合系數為-15 dB時、圖7(b)是耦合系數為-10 dB時、圖7(c)是耦合系數為-5 dB時50次蒙特卡洛仿真的統計結果。從圖7可以看出，水聽器陣元間存在耦合時，本文方法的估計性能和穩健性是最優的。

(a) 耦合系數αc=-15 dB

L1-SVD方法在不同正則化參數時的蒙特卡洛仿真結果如圖8所示。圖8(a)為正則化參數為2時、圖8(b)為正則化參數為3時，50次蒙特卡洛仿真統計結果。從圖8可以看出，改變正則化參數對L1-SVD方法的估計性能影響較大，因而該方法的實用性較差。

(a) 正則化參數為2

綜上所述，當水聽器陣元間存在互耦時，相比于MUSIC、L1-SVD、OGSBL方法，本文提出方法的估計性能最優，穩健性最高。

陣元間耦合系數為-5 dB，正則化參數為3時，四種方法信號估計結果如表2所示，表中的估計結果為50次蒙特卡洛仿真試驗的統計均值。從表中可以看出，本文方法的估計結果與真實角度的偏差最小，表明本文方法的估計性能優于MUSIC、L1-SVD、OGSBL方法。

表2 真實角度與本文方法估計DOA結果對比Tab.2 Comparison of the real angle and the DOA estimation result of the algorithm proposed in this paper

4.2 分辨力分析

不同方法的分辨能力不同，本小節從信噪比的高低、陣元間的耦合系數的大小研究本文提出方法的角度分辨能力。當滿足下述的條件時，則認為能夠正確分辨兩個角度[22]。

(45)

仿真試驗條件如表1所示，設定耦合系數αc=-7 dB，圖9給出了分辨成功率隨角度間隔的變化規律，其中(a)和(b)子圖分別為信噪比為0和5 dB時的統計結果。從圖中可以看出，信噪比為0時，三種方法成功估計出兩個角度的間隔分別為9°(本文方法)，10°(OGSBL)，12°(MUSIC)。本文方法的角度分辨力優于OGSBL和MUSIC方法。從圖9中可以看出，信噪比為0和5 dB時，可成功估計兩個角度的間隔為9°和7°。隨著信噪比的增加，本文方法的分辨力逐漸增強。

(a) SNR=0 dB 分辨成功率和角度間隔的關系

仿真試驗條件如表1所示，設定信噪比為5 dB，圖10給出了分辨成功率隨角度間隔的變化規律，其中(a)和(b)子圖分別為耦合系數為-12 dB和-3 dB時的統計結果。從圖中可以看出耦合系數為-12 dB時，三種方法成功估計出兩個角度的間隔分別為6°(本文方法)，7°(OGSBL)，8°(MUSIC)。本文方法的角度分辨力優于OGSBL和MUSIC方法。從圖10中可以看出，耦合系數為-12 dB和-3 dB時，可成功估計兩個角度的間隔為6°和11°。隨著耦合系數的增加，本文方法的分辨力逐漸降低。

(a) 耦合系數αc=-12 dB，分辨成功率和角度間隔的關系

仿真試驗條件如表1所示，設定入射角度間隔為8°和耦合系數為-7 dB時，圖11給出了分辨成功率隨信噪比的變化規律；設定入射角度間隔為11°和信噪比為5 dB時，圖12給出了分辨成功率隨耦合系數的變化規律。

圖11 分辨成功率和信噪比的關系Fig.11 The relationship between resolution probability and SNR

從圖11可以看出隨著信噪比的增加，三種方法的分辨成功率逐漸增加。從圖中可以看出，角度間隔為8°時，本文方法和OGSBL方法可成功估計出兩個角度的信噪比分別為1 dB和4 dB，而MUSIC算法無法完全正確估計出兩個真實的入射角度。本文方法的可成功估計兩個信號的信噪比低于OGSBL方法，本文方法的分辨力優于OGSBL方法和MUSIC方法。

從圖12可以看出隨著耦合系數的增加，三種方法的分辨成功率逐漸降低。從圖中可以看出，角度間隔為11°時，三種方法成功估計出兩個角度的耦合系數分別為-3 dB(本文方法)，-7 dB(OGSBL)，-7 dB(MUSIC)。本文方法的可成功估計兩個信號的耦合系數高于OGSBL方法和MUSIC方法，本文方法的分辨力優于OGSBL方法和MUSIC方法。

圖12 分辨成功率和耦合系數的關系Fig.12 The relationship between resolution probability and coupling coefficient

綜上所述，本文方法的角度分辨力優于OGSBL方法和MUSIC方法。當耦合系數一定時，分辨成功率隨著信噪比的增加逐漸提升；當信噪比一定時，分辨成功率隨著耦合系數的增加逐漸降低。

4.3 估計精度分析

本小節從信噪比、快拍數、網格間隔、陣元間的耦合系數研究不同方法的方位估計精度。把均方根誤差(Root Mean Square Error，RMSE)作為衡量方位估計精度的標準，可表示為

(46)

仿真試驗條件如表1所示，設定信噪比為-2 dB到10 dB變化，步長為2 dB，正則化參數為2，圖13給出了均方根誤差隨信噪比的變化規律。從圖中可以看出，當信噪比為10 dB時，本文方法的估計誤差為0.4°低于0.45°(OGSBL)、0.6°(L1-SVD)、0.7°(MUSIC)。隨著信噪比的增大均方根誤差逐漸減小，且本文方法的估計精度高于OGSBL 方法、L1-SVD方法和MUSIC方法。

圖13 不同信噪比，四種方法的方位估計精度Fig.13 The DOA estimation performance of the four algorithms with different SNR

仿真試驗條件如表1所示，設定快拍數為50到300變化，步長為25，正則化參數為2，圖14給出了均方根誤差隨快拍數的變化規律。從圖中可以看出，當快拍數為300時，本文方法的估計誤差為0.35°低于0.57°(OGSBL)、0.62°(L1-SVD)、0.9°(MUSIC)。隨著快拍數的增大均方根誤差逐漸減小，且本文方法的估計精度高于OGSBL 方法、L1-SVD方法和MUSIC方法。

圖14 不同快拍數下，四種方法的方位估計精度Fig.14 The DOA estimation performance of the four algorithms with different snapshot numbers

仿真試驗條件如表1所示，設定網格間隔為{1°,1.5°,2°,2.5°,3°,4°,5°,6°}變化，圖15給均方根誤差隨網格間隔的變化規律。當網格間隔小于2°時，本文提出的方法與OGSBL方法相比，估計精度相差不大，當網格間隔增大到5°時，OGSBL估計誤差將增大，高于本文方法。隨著網格間隔增大，均方根誤差增大，但本文方法的估計精度要優于OGSBL方法。

圖15 不同信噪比，網格間隔與方位估計精度的關系Fig.15 Comparison of the relationship between grid spacing and DOA estimation performance with different SNR

仿真試驗條件如表1所示，設定耦合系數為-15 dB到-3 dB，步長為2 dB，圖16給出了均方根誤差隨耦合系數變化的規律。當信噪比一定時，增大耦合系數，估計誤差逐漸增大，估計精度逐漸降低；當耦合系數為-3 dB時，增大信噪比為5 dB時，估計誤差減小3°，繼續增大信噪比為10 dB時，估計誤差減小1°。當耦合系數一定時，隨著信噪比的增大，本文方法估計誤差逐漸減小，估計精度逐漸提升。

圖16 不同信噪比條件下，耦合系數與方位估計精度的關系Fig.16 Compare the relationship between the adjacent coupling coefficient and the DOA estimation performance with different SNR

5 結 論

為了解決水聽器陣元間存在未知耦合，現有方法波達方向估計精度較低、角度分辨力差的問題，本文提出了一種基于稀疏貝葉斯學習的DOA離網格估計方法。該方法通過建立新的貝葉斯學習模型，使用EM算法，推導離網格矢量和互耦矢量的分布，重新表示空間譜函數，最后通過譜峰搜索估計出波達方向。仿真試驗結果表明：當水聽器陣元間存在互耦時，本文方法的方位估計性能和穩健性優于MUSIC、OGSBL和L1-SVD方法；本文方法的角度分辨力優于MUSIC和OGSBL方法。鑒于上述試驗使用的是兩個不相關信號，在后續的研究中，將會對相關信號做進一步研究。

附錄A

引理1：矩陣C和向量c定義如正文，對于任意向量a(θk)，我們有：

Ca=Qk(a(θk))c

式中，Qk(a(θk))∈CM×M是兩個矩陣Qk1(a(θk))、Qk2(a(θk))的和：

因此，Qk(a(θk))=Qk1(a(θk))+Qk2(a(θk))，Qk1(a(θk))和Qk2(a(θk))分別為

根據上述引理1可得CAst=Q(st?c)，即

附錄B

有復向量u∈CP×1，v∈CP×1，復矩陣A∈CM×P是復向量x∈CN×1的函數，可得