基于稀疏貝葉斯的多跳頻信號二維波達方向估計

李紅光，郭 英，眭 萍，蔡 斌, 2，蘇令華

(1. 空軍工程大學 信息與導航學院， 西安 710077； 2. 北京信息技術研究所， 北京 10094)

跳頻(FH)通信作為一種重要的擴頻通信類型，在軍事通信系統中得到了廣泛應用[1-3].FH信號的空域信息是一項重要參數，能夠為FH通信網臺分選和干擾引導等任務提供重要依據.Lin等[4]通過構造空時頻矩陣，實現FH信號波達方向(DOA)估計，但該類算法需要滿足信源數小于陣元數的超定條件，在實際應用過程中受限.陳利虎[5]將空時頻分析與多重信號分類(MUSIC)算法相結合，實現了FH信號在超定和欠定條件下的DOA估計，但是該算法復雜度較高.張東偉等[6]通過建立FH信號的極化敏感陣列觀察模型，結合空時頻分布矩陣和多項式求根方法實現DOA估計，該類算法的計算復雜度較低，但對接收陣列結構要求較為嚴格.Fu等[7-8]利用盲源分離的方法進行單源點檢測，根據估計出的載頻信息和混合矩陣實現FH信號的DOA估計，能夠實現欠定條件下的多跳頻信號DOA估計，但該算法只適用于正交組網方式的FH通信，在低信噪比條件下，單源點檢測容易受到噪聲干擾，算法魯棒性較差.上述算法均只能完成一維DOA估計，而對于FH網臺的分選和定位，二維DOA的信息特征辨識度更高.現有文獻對于FH信號的二維DOA估計較少，于欣永等[9]首先構造每一跳信號的空時頻矩陣，將共軛子空間的思想引入到MUSIC算法中，通過半普搜索實現多跳頻信號的二維DOA估計，該算法在一定程度上降低了MUSIC的計算復雜度，但只適用于不相關信源.張東偉等[10-11]通過設計陣元時頻點選取策略，構建了每跳信號空間極化時頻矩陣，并利用最小二乘旋轉不變子空間(ESPRIT)算法完成多跳信號的二維DOA估計，但該算法必須預先知道信源個數，運算量較大.

隨著壓縮感知技術不斷發展，基于稀疏重構理論的DOA估計算法[12-13]相繼出現，這類算法不易受信源數量和相干性的影響，對接收天線陣列誤差和噪聲敏感度低.目前，基于壓縮感知的DOA算法主要應用于非FH信源.韓樹楠等[14]對近似0范數稀疏類算法進行改進，利用奇異值分解估計出DOA，提高了算法在低信噪比、低快拍數條件下的DOA估計性能.但是，該算法在求解凸優化問題的過程中計算復雜度較高，而且容易受空間譜偽峰和正則化參數變化的影響，導致DOA估計不準.Lei等[15-16]采用稀疏貝葉斯學習(SBL)方法完成了DOA估計，此類算法保持了L1范數稀疏重構算法的估計性能，具有更強的可拓展性，計算復雜度也有所降低.

為了實現對FH信號二維DOA估計，更好地輔助于FH信號網臺分選，本文利用L型陣列的結構特點，將FH信號方位向和俯仰向的二維DOA估計問題轉化為兩個一維DOA估計問題，再構建空間頻率的過完備字典，通過SBL重構方法求得方位向和俯仰向的空間頻率，再由Capon[17-18]空間頻率配對算法完成FH信號頻率和空間頻率的正確配對，最終解算出FH信號跳頻率、方位角和俯仰角三維信息.為了降低算法計算復雜度，利用奇異值分解方法對兩個子陣的接收矩陣進行降維預處理，降低矩陣維數，減少運算時間.

1 數學模型

1.1 FH信號L型陣列接收模型

假設有K個遠場FH信號s(t)=[s1(t)s2(t) …sK(t)]T入射到一個L型陣列，如圖1所示，圖中αi和βi分別為FH信號si(t) (i=1,2,…,K)入射方向與子陣x及y的夾角.該L型陣列位于x-y平面，子陣x和y互成90°，每個子陣分別由M(M>K)個陣元組成，陣元間距均為d且d

圖1 FH信號L型陣列接收示意圖Fig.1 Schematic diagram of FH signal L-array reception

在觀測時間內，第i個FH信號可表示為

(1)

式中：νi為si(t)的幅度；Bi為在觀測時間內si(t)的總跳數；fi,bi為si(t)第bi跳的載頻；Ti為si(t)的跳頻周期；t′為觀測時長，t′=t-(bi-1)Ti.

假設接收陣列中各陣元是各向同性的，則第i個FH信號在子陣x和y的導向矢量分別為

(2)

(3)

式中：τi,m=(M-m)fi,bidcosθisinφi/c；υi,m=(M-m)fi,bidsinθisinφi/c，m=1,2,…,M.

由圖1所示的空間幾何關系可得

將式(4)，(5)代入式(2)，(3)得

(6)

(7)

由式(6)和(7)可知，FH信號在子陣x和y的方位角、俯仰角和跳頻率三維信息可由空間角和跳頻率二維信息表示.為了進一步降低算法求解難度，現定義第i個FH信號在頻率fi,b處的方位向空間頻率為Hi,b，俯仰向空間頻率為Vi,b，其表達式為

由式(8)和(9)可知，子陣x和y的導向矢量為

(10)

(11)

由式(10)和(11)可知，FH信號在子陣x和y只有空間頻率一維信息，則接收陣列流型矩陣可分別表示為

(12)

(13)

從而得到子陣x和y接收的FH信號分別為

1.2 基于空間頻率的FH信號稀疏表示

考慮到子陣x和y接收的FH信號空間頻率具有稀疏性，分別將方位向空間頻率和俯仰向空間頻率進行均勻離散劃分，得到超完備的空間頻率集合，分別為F={F1,F2,…,FN}和G={G1,G2,…,GN}，Fn和Gn分別表示空間頻率集合F和G中的第n個空間頻率元素，n=1,2,…,N，N為劃分網格數且N?M>K.則由式(14)和(15)可得基于空間頻率的子陣x和y單測量陣列輸出模型：

在多快拍數條件下，將子陣x和y每次接收的快拍數據L作為1個快拍幀進行運算，則基于空間頻率的稀疏表示為

(20)

同理，對于子陣y接收的FH信號也可以估計出其俯仰向空間頻率，即

(21)

2 算法原理

2.1 輸出矩陣降維預處理

在實際應用過程中，為了提高FH信號DOA的估計精度，一般需要增加FH信號的快拍數L，但是隨著快拍數L的增加，輸出矩陣的維數增多，這必然會導致計算量大幅增加，算法的時效性降低.為了在保證算法估計精度的基礎上減少算法運算時間，采用奇異值分解(SVD)的方法分別對輸出矩陣X和Y進行降維處理.根據文獻[19-20]，矩陣X和Y的SVD表達式為

式中：Ux和Uy分別為X和Y的M×M維左特征正交矩陣；Σx和Σy為M×L維對角陣；Wx=[Wx1Wx2]和Wy=[Wy1Wy2]分別為X和Y的L×L維右特征正交矩陣，Wx1和Wy1分別為Wx和Wy的前K列.Wx1和Wy1中分別包含了X和Y經過SVD分解后的特征值所對應的所有特征向量，因此可得X和Y經SVD降維后的矩陣Xsvd和Ysvd分別為

Xsvd和Ysvd包含了子陣x和y接收的所有FH信號空間頻率信息，可將式(18)和(19)轉化為

本節利用SVD算法將M×L維輸出矩陣X和Y降低到M×K維，且SVD算法本身對子陣接收的噪聲產生抑制作用.由于噪聲的奇異值相對于FH信號非常小，所以在SVD過程中能夠將FH信號中的大部分噪聲分離出來，提高了算法在低信噪比條件下的估計性能.

2.2 基于SBL的空間頻率估計

本節以子陣x為例，求解方位向空間頻率.根據式(26)，其觀測噪聲矩陣Ex-svd各列均是相互獨立的，且服從均值為0，協方差矩陣為(σ2/L)IM(IM為M×M維單位矩陣)的高斯分布.則Xsvd關于FH信號幅度的概率密度函數可表示為

(28)

本節SBL引入中間參數γ=[γ1γ2…γN]T表示FH信號在空間頻率集Dx上的功率譜分布，各參數間相互獨立，且有

(29)

式中：N(0,Γ)為均值為0，協方差矩陣為Γ的高斯分布，Γ=diag(γ).

(30)

(31)

(32)

式中：US為均值；ΣS為方差.

對式(30)取對數并忽略其常數項，即可得到優化γ的目標函數：

(γ,

(33)

(34)

通過式(34)的j次迭代運算后，可估計出所有FH信號的方位向空間頻率集

則其對應噪聲功率的無偏估計值為

(35)

(36)

式中：λm表示Xsvd的第m個奇異值.

2.3 Capon空間頻率配對方法

根據式(8)，其方位向空間角可表示為

αi=arccos(cHi,b/fi,b)

(37)

可見，只要將FH信號的跳頻率和空間頻率正確配對，即可估計出方位向和俯仰向的空間角，從而解算出方位角、俯仰角和跳頻率三維信息.

將式(14)進行傅里葉變換，可得到子陣x接收的混合信號在fi,b處的頻域輸出為

x(fi,b)=Ax(Hi,b)sx(fi,b)+ex(fi,b)

(38)

式中：sx(fi,b)為子陣x接收的FH信號在fi,b處的頻域輸出；ex(fi,b)為子陣x接收的噪聲在fi,b處的頻域輸出.

根據Capon波束形成思想，對方位向空間頻率上的FH信號功率進行歸一化處理，并最小化子陣x輸出功率，其約束方程為

(39)

式中：ω為最優權值系數，根據文獻[24]有

(40)

Rx(fi,b)為頻率fi,b處的頻域協方差矩陣，

Rx(fi,b)=E[x(fi,b)xH(fi,b)]

(41)

(42)

(43)

(44)

(45)

(46)

2.4 算法流程及復雜度分析

本文算法的具體步驟如下：

(1) 構造方位向和俯仰向空間頻率字典F和G.

(2) 根據式(18)和(19)構造方位向和俯仰向空間頻率的陣列輸出矩陣X和Y，流型矩陣Dx和Dy.

(3) 根據式(24)和(25)降低輸出矩陣X和Y經過SVD分解得到降維后的矩陣Xsvd和Ysvd，減少計算量.

τ=

(47)

圖2 相干FH信號的空間頻率估計和空間角與頻率配對結果Fig.2 Spatial frequency estimation and pairing results of coherent FH signals

3 實驗結果與分析

本節通過實驗驗證，對比分析本文算法和文獻[13]基于MUSIC的對稱壓縮譜算法(MSCS)的二維DOA估計性能.假設有若干個遠場FH信號分別從不同角度入射到L型均勻陣列上，子陣陣元數均為8，陣元間距為1 m，角度范圍 -90°～90°，頻率范圍0～100 MHz，則空間頻率Hi,b和Vi,b的范圍為 -0.333 3～0.333 3，超完備字典網格數N=400.

實驗1驗證本文算法對相干FH信號的空間頻率估計以及與跳頻率正確配對的有效性.假設兩個遠場相干FH信號的入射角分別為(-25°，57°)和(33°，-48°)，FH周期均為100 μs，頻率集均為{41,62,87} MHz，采樣率為200 MHz，采樣時長為1 ms，信噪比為0.經計算，兩個FH信號的空間頻率分別為(0.103 9，-0.048 4)，(0.157 1，-0.073 3)，(0.220 4，-0.102 8)和(-0.085 2，-0.055 3)，(-0.128 8，-0.083 6)，(-0.180 7，-0.117 4).圖2所示為兩個相干FH信號在快拍數為80條件下的空間頻率估計和配對結果.圖中：P為功率；Spf為空間頻率；Spa為空間角.從圖2(a)和2(b)可知，SBL-Capon算法能夠實現對相干FH信號空間頻率的有效估計.從圖2(c)和2(d)可知，在完成空間頻率正確估計后，利用基于FFT變換的頻率估計和Capon空間頻率配對算法， 將空間頻率和跳頻率進行正確配對并得到空間角，根據入射角空間幾何關系，計算出方位角和俯仰角，實現了FH信號跳頻率和二維DOA的聯合估計.由圖2(c)還可以看出，SBL-Capon算法能夠將FH信號頻率集進行分選，即相同空間角度的跳頻率來自于同一個FH信號，當觀測時間足夠長時，即可完成FH信號的跳頻率集估計，能夠有效輔助FH網臺分選.

實驗2對比分析SBL-Capon和 MSCS兩種算法在不同空間頻率間隔下的空間頻率估計精度.由空間頻率定義和空間角幾何關系可知，空間頻率的估計精度對DOA的解算精確性有重要影響.為了方便控制仿真中各FH信號間的空間頻率間隔，由于MSCS算法只適用于非相干信源，所以取入射角相近的兩個非相干FH信號.其入射角分別設為(25°，64°)和(27°，62°)，頻率集均為{71,86,97} MHz，跳周期分別為100和50 μs，采樣率為200 MHz，采樣時長為1 ms，信噪比0 dB，快拍數為80.經計算，兩個FH信號的空間頻率分別為(0.192 8,0.089 9)，(0.233 5,0.108 9)，(0.263 4,0.122 8)和(0.186 2,0.094 9)，(0.225 5,0.114 9)，(0.254 4,0.129 6).圖3為SBL-Capon和 MSCS兩種算法的空間頻率估計結果.從圖3(b)和3(d)可知，當FH信號源的空間頻率間隔較小時(即空間角度相近)，SBL-Capon算法的估計精度優于 MSCS算法的估計精度.

圖3 非相干FH信號的空間頻率估計結果Fig.3 Spatial frequency estimation results of non-coherent FH signals

圖4所示為空間頻率估計的方均根誤差隨最小空間頻率間隔的變化情況，圖中ΔS為最小空間頻率間隔，空間頻率范圍為0.01～0.1，最小間隔0.02，空間頻率的方均根誤差為

RMSEFV=

(48)

從圖4可知，當兩個FH信號的最小空間頻率間隔較小時，SBL-Capon的空間頻率估計性能優于MSCS算法的空間頻率估計性能.隨著空間頻率間隔增加，MSCS算法的空間頻率估計精度逐漸提高，SBL-Capon算法的估計性能基本保持不變，不易受空間頻率間隔影響，當空間頻率間隔大于0.07時，SBL-Capon和MSCS算法的估計性能基本一致.

實驗3對比分析SBL-Capon和 MSCS兩種算法在不同信噪比下DOA方均根誤差和估計成功率.假設：3個非相干遠場FH信號，入射角在 -90°～90° 范圍內隨機產生，角度間隔不小于2°；FH周期分別為150，100和80 μs；頻率集分別為65～95 MHz范圍內的3個隨機頻率，且同一個FH信號的相鄰跳點頻率間隔大于5 MHz；采樣率為200 MHz，采樣時長為1 ms，快拍數為80；信噪比范圍 -10～20 dB，間隔為2 dB.

圖4 空間頻率估計的方均根誤差隨最小空間頻率間隔變化Fig.4 Root mean square error varies of the spacial frequency with minimum spatial frequency interval

DOA方均根誤差RMSEDOA和估計成功率ηDOA分別為

RMSEDOA=

(49)

(50)

圖5 不同信噪比下DOA估計性能Fig.5 DOA estimation performance under different SNR

圖6 不同快拍數下DOA估計性能Fig.6 DOA estimation performance under different snapshots

圖5所示為RMSEDOA和ηDOA隨信噪比SNR增加的變化情況.從圖5可知，隨著信噪比增加，兩種算法的DOA估計方均根誤差逐漸減小，估計成功率逐漸提高.在低信噪比條件下，SBL-Capon算法的估計精度和成功率均高于MSCS算法.主要原因是由于SVD降維預處理和SBL稀疏求解過程都對噪聲起到抑制作用，當信噪比大于10 dB時，兩種算法的估計性能基本一致.

實驗4對比分析SBL-Capon和 MSCS兩種算法在不同快拍數下DOA方均根誤差和估計成功率.假設：3個非相干遠場FH信號，入射角在 -90°～90° 范圍內隨機產生，角度間隔不小于2°；FH周期分別為150，100和80 μs；頻率集分別為 65～95 MHz范圍內的3個隨機頻率，且同一FH信號的相鄰跳點頻率間隔大于5 MHz；采樣率為200 MHz，采樣時長為1 ms，信噪比0 dB；快拍數范圍50～550，間隔50.圖6所示是RMSEDOA和ηDOA隨快拍數變化情況.從圖6可知：隨著快拍數的增加，兩種算法的DOA估計性能逐漸提高；在低快拍數條件下，SBL-Capon的DOA估計精度和成功率均優于MSCS算法；當快拍數大于450后，兩種算法的估計性能基本保持不變.

實驗5對比分析SBL-Capon算法在不同字典網格數N下DOA方均根誤差RMSEDOA和運算時間t.假設：3個非相干遠場FH信號，入射角在 -90°～90° 范圍內隨機產生，角度間隔不小于2°；FH周期分別為150，100和80 μs；頻率集分別為65～95 MHz范圍內的3個隨機頻率，且同一FH信號的相鄰跳點頻率間隔大于5 MHz；采樣率為200 MHz，采樣時長為1 ms，快拍數為80；信噪比范圍 -10～20 dB，間隔2 dB；N取值分別為300，400和500.運算時間均是在MATLAB2014環境下，通過TIC和TOC指令計算平均運算時間t.

圖7所示為不同網格數N條件下RMSEDOA和t隨信噪比增加的變化情況.從圖7可知，在相同信噪比條件下，字典網格數越小，算法運算時間越短.這主要是由于字典網格數減少使得流型矩陣各列之間相關性減弱，從而減小了超參數求解的迭代次數，加快了算法收斂，但是DOA估計精度也有所下降.因此，在實際應用中，要根據FH信號帶寬范圍合理劃分字典間隔.

圖7 不同網格數下DOA估計性能和運算時間Fig.7 DOA estimation performance and operation time under different grid numbers

4 結語

為了更好地輔助FH信號網臺分選，本文提出一種基于稀疏貝葉斯學習的FH信號二維DOA估計算法.該算法首先根據入射角的空間幾何關系，將方位角、俯仰角和跳頻率三維信息轉換為一維空間頻率信息，降低了稀疏重構算法復雜度.然后經過SVD降維處理，減少矩陣運算維數，通過SBL算法估計出空間頻率，利用FFT頻率估計算法和Capon空間頻率配對算法將空間頻率和跳頻率正確配對，解算出方位角和俯仰角.實驗結果表明，在低信噪比或低快拍數條件下，該算法能有效估計出相干或非相干多FH信號的二維DOA信息，而且不易受空間頻率間隔影響.如何進一步提高本算法時效性，以及實現FH信號DOA和極化參數的聯合估計是下一步研究工作.