高效的多跳頻信號2D-DOA估計算法

于欣永, 郭 英, 張坤峰, 眭 萍, 李 雷, 李紅光, 孟 濤

(1. 空軍工程大學信息與導航學院, 陜西 西安 710077; 2. 通信網信息傳輸與分發技術重點實驗室,河北 石家莊 050081; 3. 空軍通信士官學校, 遼寧 大連 116100)

0 引 言

跳頻通信因其保密性好、抗干擾能力強、組網能力強等特點,廣泛應用于航空航天通信領域,如Link16、美軍軍事衛星等,同時也向通信偵察提出了嚴峻的考驗[1-3]。網臺分選是跳頻信號偵查的一個關鍵環節,它可以從復雜的電磁環境中將不同跳頻電臺分選出來,實現敵我識別,對敵方重要的通信電臺實施跟蹤、干擾、打擊、破壞,為軍事通信對抗情報生成、干擾引導及長期情報偵查提供依據和支撐。

針對跳頻信號的偵查處理目前大多側重于信號檢測和時頻域的參數估計,而對跳頻信號的空域信息研究較少[4-6]。而跳頻信號的空域信息如波達方向等在網臺分選、信號跟蹤和干擾中具有重要作用。文獻[7-8]利用構造空時頻模型進行跳頻信號的波達方向(direction of arrival,DOA)估計,但是文中算法只能適用于超定情況,而實際復雜電磁環境中,大多情況處于欠定條件下;文獻[9]運用空時頻方法進行線性調頻信號DOA估計及信號分離,取得了良好的效果,但是只適用于線性調頻信號;文獻[10]結合空時頻模型和多重信號分類(multiple signal classification,MUSIC)算法進行跳頻信號DOA估計,為跳頻信號DOA估計提出了新的思路,這種方法可以適用于適定與欠定條件,估計精度較理想,但是算法復雜度高;文獻[11-12]用root-MUSIC算法代替了文獻[10]中的MUSIC算法,降低了計算復雜度,但是文中算法對陣列結構要求苛刻,實際工程應用中難以滿足條件。文獻[7-12]所提算法都只是關注的一維DOA信息,沒有考慮二維DOA(two dimensional-DOA,2D-DOA)問題,無法定位三維的空間目標。針對跳頻信號的2D-DOA估計問題,現有的研究較少,文獻[13-14]將空間極化時頻分布思想引入到跳頻信號參數估計中,并結合旋轉不變子空間(estimation of signal parameters via rotation invariant technique,ESPRIT)算法實現跳頻信號的2D-DOA估計,可以在平行線陣和“L型”線陣下高精度估計出多跳頻信號二維波達方向,但是這種算法要求信噪比高,而且ESPRIT需要進行參數配對,增加了計算量,降低了算法的估計性能。文獻[15-16]提出了一種MUSIC對稱壓縮譜(MUSIC symmetrical compressed spectrum,MSCS)算法估計連續信號波達方向算法,這種算法在MUSIC算法基礎上減少了譜峰搜索的次數,算法復雜度較小,但是只適用于連續信號。

為了利用跳頻信號的空2D-DOA信息輔助多跳頻信號的網臺分選,本文提出一種空時頻分析與MSCS相結合的多跳頻信號2D-DOA估計的算法。首先在時頻圖上提取跳頻信號的自相時頻點,構造每個hop的空時頻矩陣(spatial time-frequency distribution,STFD),進而得到時頻域的協方差矩陣;然后將共軛子空間的思想引入到傳統MUSIC算法中,通過對噪聲子空間及其共軛的交集進行奇異值分解,實現噪聲子空間的降維,進而構造新的降維空間譜函數;最終根據新的降維空間譜,只需在半譜內峰值搜索,實現跳頻信號的2D-DOA估計。最后,對估計出來的2D-DOA信息進行聚類實現跳頻信號的網臺分選。另外,本文在時頻圖處理時采用形態學濾波進行去噪,提高了算法在低信噪比時的魯棒性。仿真實驗表明所提算法與MUSIC算法相比,在保證估計性能的同時將計算復雜度降低了一半。

1 跳頻信號的陣列快拍模型

假設,在觀測時間Δt內共有K個跳頻信號sn(t),其跳周期為Tn,第k跳的載頻、初相為ωnk、φnk,起始跳的時長為Δt0n,則sn(t)的表達形式[17]為

(1)

式中,t′=t-(k-1)Tn-Δt0n;υn是信號sn(t)的基帶復包絡;rect表示單位矩形脈沖。

假設N個跳頻信號被如圖1所示的陣列接收,圖中兩均勻線陣(uniform linear array, ULA)ULA1、ULA2的陣元數為M,陣元間距為d1,陣列的間距為d2,且滿足max(d1,d2)

圖1 雙平行均勻線陣Fig.1 Double parallel uniform linear array

跳頻信號的載頻隨著偽隨機碼隨機跳變,是一個寬帶信號,但是只分析某一跳時,可以當作分段窄帶信號進行處理。假設跳頻信號S以方位角θ∈[0,2π]和俯仰角φ∈[0,π/2]從遠場入射到接收陣列上。陣元增益為1,則子陣ULA1對入射波S的導向矢量aS1(θ,φ)為

aS1(θ,φ) =[1,e-j2πd1sin θcos φ/λ,…,e-j2π(M-1)d1sin θcos φ/λ]T

(2)

式中,λ為信號的波長。ULA2對入射波S的導向矢量為aS2(θ,φ),相比較aS1(θ,φ)差別是由陣列距離d2引起,則aS2(θ,φ)為

aS2(θ,φ) =[1,e-j2π(d1sin θcos φ+d2sin θsin φ)/λ,…,

e-j2π[(M-1)d1sin θcos φ+d2sin θsin φ]/λ]T

(3)

則陣列對入射波S的導向矢量a(θ,φ)為

(4)

由導向矢量組成的陣列流型矩陣為

A=[a1(θ,φ),a2(θ,φ),…,aN(θ,φ)]

(5)

則跳頻信號的陣列快拍模型為

X(t)=Y(t)+N(t)=AS(t)+N(t)

(6)

式中,X(t)=[x1(t),x2(t),…,xM(t)]T為陣列接收到的數據矢量;S(t)=[s1(t),s2(t),…,sN(t)]T為跳頻信源數據矢量;N(t)=[n1(t),n2(t),…,nM(t)]T為噪聲數據矢量。

2 跳頻信號STFD構造

2.1 STFD基本原理

信號xi(t)的Cohen類離散時頻分布為

(7)

式中,φ(l,τ)為核函數。信號xi(t)與xj(t)離散形式的Cohen類互時頻分布為

(8)

定義信號x(t)的空時頻分布為

(9)

式中,[DXX(t,f)]ij=Dxixj(t,f)表示各陣列輸出信號之間的互時頻分布。根據式(6)和式(9)得到接收陣列的時頻域協方差矩陣為

E[DXX(t,f)]=E[DYY(t,f)]+E[DNN(t,f)]=

ADSS(t,f)AH+E[DNN(t,f)]

(10)

2.2 高精度STFD構造

選取自項時頻點構造STFD是2D-DOA估計性能的前提,特別在低信噪比情況下,短時傅里葉變換得到的時頻圖效果不佳。本文利用形態學濾波的方法對時頻圖進行修正,在修正的時頻圖上完成各個hop的提取。

假設時頻圖上任意一時頻點頻率的模值為F(t,f)∈[Fmin(t,f),Fmax(t,f)],對短時傅里葉變換得到的跳頻信號時頻圖做閾值分割處理轉化為二值圖像,處理后的時頻圖矩陣為

(11)

本文閾值ε通過自適應信噪比迭代算法來確定，算法流程如下。

步驟1假設初始閾值為ε0=[Fmin(t,f)+Fmax(t,f)]/2;

步驟2按照式(11)用ε0截取短時傅里葉得到的時頻圖,把時頻圖劃分為兩個部分;

步驟3對劃分后的兩部分分別計算其平均頻率模值,并取其和的一半作為新的閾值ε1;

步驟4用ε1按照式(11)重新進行截取,按照步驟3一直迭代至第N+1次的閾值與第N次的閾值相等,即該閾值為最佳分割閾值ε。

fΘB=min{f(x-x′,y-y′)-b(x′,y′)|(x′,y′)∈Db}

(12)

f⊕B=max{f(x-x′,y-y′)+

b(x′,y′)|(x′,y′)∈Db}

(13)

式中,b(x′,y′)為結構元素;Db為b(x′,y′)的定義域。腐蝕能夠消除影響對某一跳提取的散點,同時會使邊界收縮;膨脹可以填補信號中的空洞,而且會將邊界擴張。因此,通過腐蝕和膨脹處理可以消除時頻圖中的干擾噪聲點以及散點而對圖像中的其他點沒有影響,進而得到一個清晰穩健的時頻圖,再通過坐標變換將該時頻圖從二值圖像轉換為原始圖像。

圖2(a)、圖2(b)分別為3個同步跳頻信號在信噪比為-5 dB下只經過短時傅里葉變換的時頻圖和采用形態學濾波修正之后的時頻圖。可見,短時傅里葉變換得到的時頻圖背景噪聲大,時頻聚焦性差;采用形態學濾波去噪之后的時頻圖能量更集中,更有利于自相時頻點的提取。

圖2 多跳頻信號時頻分布圖Fig.2 Time frequency distribution of multi-frequency-hopping signals

3 基于MSCS的2D-DOA估計

3.1 MUSIC算法原理

假設陣列同時接收到跳頻信號的個數為L,則時頻域協方差矩陣E[DXX(t,f)]經過特征值分解可表示為

(14)

式中,US表示信號子空間;UN表示噪聲子空間;Σ為E[DXX(t,f)]特征值組成的對角矩陣。由噪聲子空間UN構造的MUSIC算法空間譜函數PMUSIC(θ,φ)為

(15)

根據PMUSIC(θ,φ)在(θ,φ)二維域進行譜峰搜索,使得PMUSIC(θ,φ)取得極值的(θ,φ)即為所求2D-DOA,由于傳統MUSIC算法需在全場范圍內進行譜峰搜索,使得算法復雜度過高。

3.2 MSCS算法原理

(16)

對式(2)和式(3)的導向矢量兩端取共軛可得

(17)

(18)

由式(17)和式(18)可得

(19)

對式(16)兩端取共軛,利用式(19)得

(20)

(21)

由式(21)可知,當φ∈[0,π],有

(22)

同理可得,當φ∈[π,2π]時,有

PMSCS(θ,φ-π)=PMSCS(θ,φ)

(23)

令噪聲子空間UN=[UN1,UN2,…,UNM-L],則

(24)

將式(24)代入式(21)并取倒數可得

(25)

由式(16)和式(25)可得

(26)

3.3 2D-DOA估計

對R進行奇異值分解得

R=UΛVH

(27)

式中,Λ為對角陣,則

(28)

(29)

4 算法流程及復雜度分析

4.1 算法流程

基于以上分析論述,本文所提的基于MSCS的多跳頻信號2D-DOA估計算法具體步驟如下。

步驟1根據式(11)采用自適應信噪比作為閾值對短時傅里葉得到的跳頻信號時頻圖二值化處理,然后用式(12)、式(13)對處理后的圖像進行形態學濾波修正,得到清晰穩健的時頻圖;

步驟2在修正后的時頻圖上提取每個hop的自項時頻點根據式(9)構建STFD,并根據式(10)求得其時頻域的協方差矩陣;

步驟5根據式(29)構建MSCS空間譜函數,在(θ,φ)二維域進行半譜峰搜索,完成2D-DOA的高效估計。

4.2 算法復雜度分析

假設根據空間譜函數進行譜峰搜索時,搜索步長為ζ,在MUSIC算法中,根據式(15)構造空間譜函數時算法的復雜度為3M2-2ML,則MUSIC算法估計多跳頻信號2D-DOA的計算量G為

G=π(3M2-2ML)/ζ

(30)

(31)

(32)

可以看出,在進行跳頻信號2D-DOA估計時,相比較MUSIC算法,MSCS的算法復雜度降低一半左右。

5 仿真與分析

假設ULA1、ULA2的陣元間距d1為2.5 m,陣列的間距d2為3 m;空間中存在5個跳頻信號FH1～FH5的二維DOA參數(θ,φ)分別為(20°,30°),(40°,50°),(60°,70°),(70°,60°),(80°,40°);跳周期均為10 μs,歸一化載頻在0～0.5隨機跳變,采樣率為100 MHz。

實驗2～實驗5均為經過形態學濾波處理的條件下做的仿真,每次實驗進行200次蒙特卡羅計算,選用均方根誤差(root mean square error, RMSE)和估計成功率作為衡量算法性能好壞的標準。RMSE的定義為

(33)

η=N1/N

(34)

式中,N1代表200次蒙特卡羅實驗中方位角(俯仰角)估計誤差小于2°的實驗次數,N代表總實驗次數,即N=200。

5.1 實驗1

為了驗證形態學濾波對本文所提算法進行跳頻信號2D-DOA估計性能的影響,假設陣列同時接收到的3個跳頻信號FH1、FH2和FH3,各跳的快拍數為2 000,陣列ULA1和ULA2的陣元數均為4,信噪比以2 dB為步進從-10 dB遞增到30 dB,有無經過形態學濾波修正的MSCS算法估計性能如圖3和圖4所示。

圖3 實驗1參數估計成功率Fig.3 Estimated success rate of parameters in experiment 1

圖4 實驗1參數的RMSEFig.4 RMSE of parameters in experiment 1

實驗1的結果表明,不管有無經過形態學濾波處理,MSCS算法方位角及俯仰角的估計成功率都隨著信噪比的增加而逐漸增大,RMSE都隨著信噪比的增加而逐漸減小;當信噪比小于10 dB時,經形態學濾波處理的估計成功率和RMSE均優于未經形態學濾波處理的算法;當信噪比大于10 dB時,兩種情況估計性能相當。

5.2 實驗2

假設空間中存在3個跳頻信號FH1、FH2和FH3,各跳的快拍數為2 000,ULA1和ULA2的陣元數均為4,信噪比以2 dB為步進從-10 dB遞增到30 dB,MSCS算法和MUSIC算法的2D-DOA估計性能曲線如圖5和圖6所示。

實驗2的結果表明,兩種算法對跳頻信號的方位角及俯仰角的估計成功率隨著信噪比的增加而逐漸增大,而RMSE隨著信噪比的增加而逐漸減小;當信噪比達到10 dB左右時,兩種算法的參數估計成功率均達到100%;整體上MSCS算法的估計成功率要略高于MUSIC算法;信噪比在-10～10 dB時,MSCS算法的估計精度要略差于MUSIC算法,但是當信噪比大于10 dB時兩種算法的估計精度相當。

圖5 實驗2參數估計成功率Fig.5 Estimated success rate of parameters in experiment 2

圖6 實驗2參數的RMSEFig.6 RMSE of parameters in experiment 2

5.3 實驗3

為了驗證陣元數對算法的影響,假設空間中存在3個跳頻信號分別FH1、FH2和FH3,各跳的快拍數為2 000,ULA1、ULA2的陣元數以2為步進從4遞增到16,在信噪比分別0 dB和20 dB時,MSCS算法和MUSIC算法的估計性能如圖7和圖8所示。

圖7 實驗3參數估計成功率Fig.7 Estimated success rate of parameters in experiment 3

圖8 實驗3參數的RMSEFig.8 RMSE of parameters in experiment 3

實驗3的結果表明,兩種算法的估計成功率都隨著陣元數的增加而增加,RMSE隨之減小;在相同的信噪比下,MSCS算法參數估計成功率要高于MUSIC算法;陣元數大于10時,MSCS算法的整體估計精度要略低于MUSIC算法,陣元數小于10時,兩種算法估計精度相當。造成這樣現象是由于隨著陣元數的增加,MSCS噪聲子空間維度下降而使得空間譜變得更加的“尖銳”。

5.4 實驗4

為了驗證快拍數對算法性能的影響,假設空間中的3個跳頻信號分別為FH1、FH2和FH3,ULA1、ULA2的陣元數是4,各跳的快拍數從1 000遞增到3 000,在信噪比分別為0 dB和20 dB時,MSCS算法和MUSIC算法的估計性能如圖9和圖10所示。

圖9 實驗4參數估計成功率Fig.9 Estimated success rate of parameters in experiment 4

實驗4的結果表明,兩種算法的DOA估計性能都隨著快拍數的增加而變優;在相同的信噪比下,MSCS算法的估計成功率總體上要略高于MUSIC算法;但是相同快拍數下MSCS算法的估計精度要略低于MUSIC算法。

圖10 實驗4參數的RMSEFig.10 RMSE of parameters in experiment 4

5.5 實驗5

為了驗證MSCS算法相對于MUSIC算法在復雜度上的優勢,假設空間中5個跳頻信號FH1～FH5全部存在,各跳快拍數為3 000,信噪比以4 dB為步進從-8 dB遞增到20 dB時，兩種算法2D-DOA估計所需的時間如表1所示。

表3 2D-DOA估計所需時間

實驗5的結果表明,在不同的信噪比下,MSCS算法2D-DOA估計的時間約為MUSIC算法的估計時間的一半左右;由此可見,相比較傳統MUSIC算法,本文算法的復雜度降低了一半左右。

6 結 論

跳頻信號的2D-DOA信息能夠有效地輔助多跳頻信號網臺分選,多跳頻信號的2D-DOA估計有著重要的意義。本文詳細推導和闡述了基于空時頻分析與MSCS算法相結合的多跳頻信號2D-DOA高效估計算法,為了提高所提算法在低信噪比下的魯棒性,將形態學濾波的思想引入到時頻圖處理中。理論論證和實驗仿真表明,與傳統MUSIC算法相比,本文所提算法在保證估計性能的同時將算法復雜度降低了一半左右。

參考文獻：

[1] SHA Z C. Online hop timing detection and frequency estimation of multiple FH signals[J]. Etri Journal, 2013, 35(5):748-756.

[2] IBRAHIM M, GALAL I. An improved SDR frequency tuning algorithm for frequency hopping systems[J]. Etri Journal, 2016,8(3):89-105.

[3] HANAWAL M, ABDELRAHMAN M, KRUNZ M. Joint adaptation of frequency hopping and transmission rate for anti-jamming wireless systems[J]. IEEE Trans.on Mobile Computing, 2016, 15(9):2247-2259.

[4] 趙新明,金艷,姬紅兵.a穩定分布噪聲下基于Merid濾波的跳頻信號參數估計[J].電子與信息學報,2014,36(8):1878-1883.

ZHAOX M, JIN Y, JI H B. Parameter estimation of frequency-hopping signals based on Merid filter in a stable noise environment[J]. Journal of Electronics &Information Technology, 2014, 36(8):1878-1883.

[5] 張坤峰,郭英,齊子森,等.基于稀疏貝葉斯重構的多跳頻信號參數估計[J].華中科技大學學報:自然科學版,2017,45(1):97-102.

ZHANG K F, GUO Y, QI Z S,et al. Parameter estimation for multiple frequency-hopping signals based on sparse Bayesian reconstruction[J]. Journal of Huazhong University of Science and Technology, 2017,45(1):97-102.

[6] FU K C, CHEN Y F. Blind iterative maximum likelihood-based frequency and transition time estimation for frequency hopping systems[J]. IET Communications, 2013, 7(9):883-892.

[7] LIN C H, FANG W H. Joint angle and delay estimation in frequency hopping systems[J]. IEEE Trans.on Aerospace & Electronic Systems, 2013, 49(2):1042-1056.

[8] LIU X, LI J, MA X. An EM algorithm for blind hop timing estimation of multiple FH signals using an array system with bandwidth mismatch[J]. 2007, 56(5):2545-2554.

[9] ZHANG Y, MA W, AMIN M G. Subspace analysis of spatial time-frequency distribution matrices[J]. IEEE Trans.on Signal Processing, 2001, 49(4):747-759.

[10] 陳利虎, 王永明, 張爾揚. 基于空時頻分析的多FH/DS信號DOA估計[J]. 信號處理, 2009, 25(8):1309-1313.

CHEN L H, WANG Y M, ZHANG E Y. Directions of arrival estimation for multicomponent frequency-hopping/direct sequence spread spectrum signal based on spatial time-frequency analysis[J]. Signal Processing, 2009, 25(8):1309-1313.

[11] ZHANG C, LI L. Parameter estimation of multi frequency hopping signals based on compressive spatial time-frequency joint analysis[J]. Pacific Journal of Mathematics, 2014, 136(1): 85-101.

[12] 張東偉, 郭英, 齊子森,等. 多跳頻信號波達方向與極化狀態聯合估計算法[J]. 電子與信息學報, 2015,37(7):1695-1701.

ZHANG D W, GUO Y, QI Z S, et al. Joint estimation algorithm of direction of arrival and polarization for multiple frequency hopping signals[J]. Journal of Electronics and Information Technology, 2015,37(7):1695-1701.

[13] 張東偉,郭英,張坤峰,等.多跳頻信號頻率跟蹤與二維波達方向實時估計算法[J].電子與信息學報,2016,38(9):2377-2384.

ZHANG D W, GUO Y, ZHANG K F, et al. Online estimation algorithm of 2D-DOA and frequency tracking for multiple frequency-hopping signals[J]. Journal of Electronics & Information Technology, 2016, 38(9): 2377-2384.

[14] 張東偉, 郭英, 齊子森, 等. 跳頻信號2D-DOA與極化參數的欠定估計[J]. 哈爾濱工業大學學報, 2016, 48(4):121-128.

ZHANG D W, GUO Y, QI Z S, et al. Underdetermined estimation of 2D-DOA and polarization for frequencyhopping signals[J].Journal of Harbin Institute of Technology,2016,48(4):121-128.

[15] 閆鋒剛,劉帥,金銘,等.基于MUSIC對稱壓縮譜的快速DOA估計[J].系統工程與電子技術,2012,34(11):2198-2202.

YAN F G, LIU S, JIN M, et al. Fast DOA estimation based on MUSIC symmetrical compressed spectrum[J]. System Engineering and Electronics, 2012, 34(11):2198-2202.

[16] 閆鋒剛, 劉帥, 金銘, 等. 基于降維噪聲子空間的二維陣列DOA估計算法[J]. 電子與信息學報, 2012, 34(4):832-837.

YAN F G, LIU S, JIN M, et al. 2-DDOA Estimation method based on dimension descended noise subspace[J]. Journal of Electronics & Information Technology, 2012, 34(4):832-837.

[17] 馮永新,徐美榮,錢博,等.一種差分跳頻頻率轉移函數算法[J].航空學報,2013,34(3):655-661.

FENG Y X, XU M R, QIAN B, et al. A frequency transform function algorithm for differential frenquency hopping[J]. Chinese Journal of Aeronautics, 2013, 34(3):655-661.