基于復數據經驗模式分解的天波超視距雷達瞬態干擾抑制

周忠根 水鵬朗

(西安電子科技大學雷達信號處理國家重點實驗室 西安 710071)

1 引言

天波超視距雷達(OTHR)工作在高頻(HF)頻段(3～30 MHz)，利用發射的高頻電磁波遇到電離層折射、后向返回散射路徑下視傳播的原理來實現對視距以外遠程目標的探測。它已廣泛應用在軍事和民用領域[1,2]。然而，天波超視距雷達工作的電磁環境非常復雜，不僅存在極強的海(地)雜波，而且存在短波無線電通訊[3]、雷電[4]、流星余跡[5]以及各種人工干擾等多種干擾信號，同時還有大氣噪聲、宇宙噪聲以及各種人為噪聲等。在OTHR中，把存在時間遠小于相干積累時間的干擾稱為瞬態干擾，主要包括無線電通訊、雷電、流星余跡等干擾信號。瞬態干擾在目標檢測多普勒譜域內有很寬的譜，它會影響到幾乎所有或大部分多普勒單元內目標信號的檢測。因此，必須對瞬態干擾進行抑制以提高雷達檢測性能。

目前己有多種時域檢測方法被提出，其中大多數方法需預先抑制海(地)雜波[4,6-10]。文獻[4,6,7]直接在多普勒譜中利用頻域濾波的方式抑制海(地)雜波，該方法簡單直觀、易于工程實現。然而，由于使用的是一個理想高通濾波器，信號變換回時域時會出現拖尾效應。如果有一個瞬態干擾的幅度明顯大于其它干擾的幅度，其拖尾有可能擴展到幾乎整個相干積累時間，那么其它干擾信號就會被遮蔽。這樣低幅度的瞬態干擾就不能被有效的檢測到。文獻[8]用相鄰距離、方位單元上的海(地)雜波在短時間內的高度相關性來進行雜波對消，即對來自于相鄰距離、方位單元的時域回波作協方差矩陣，然后利用特征分解和主元分析法獲得海(地)雜波的估計，最后用原始數據減去該估計值，達到抑制海(地)雜波的目的。文獻[9,10]將雷達回波信號分段構造成矩陣并進行奇異值分解，然后利用奇異值分解的正交性實現雷達回波的正交分解，使海(地)雜波和瞬態干擾分離開來。這些抑制海(地)雜波的方法不涉及頻域濾波，但運算量相對較大。海(地)雜波抑制完后，瞬態干擾就在時域突顯出來，就可確定瞬態干擾所在的時域位置。然后，該時域位置的數據將被挖除并重構。對于數據的重構，一般采用AR模型[8-11]來恢復挖掉的數據。AR模型有兩個缺陷：一方面，為保證預測誤差達到最小，需精心選擇AR模型的階數，計算復雜度較高。另一方面，當一個相干積累時間內存在多個瞬態干擾時，如果將所有的瞬態干擾都挖除，那么剩余的數據量可能太少。這時用AR模型作數據重構時，沒有充分利用全部數據，使得重構效果不佳，進而影響雷達檢測性能。

本文主要研究利用復數據經驗模式分解(Complex Empirical Mode Decomposition,CEMD)來抑制瞬態干擾。1998年，Huang等人[12]提出的希爾伯特-黃變換(HHT)是一種新的時頻分析方法，被認為是近年來對線性平穩信號分析的一個重大突破，能夠很好地處理非平穩和非線性信號。經驗模式分解(Empirical Mode Decomposition,EMD)是它的核心部分，其目的就是把復雜信號分解成多個基本模式分量的和，這些基本模式分量被稱作固有模態函數(IMF)。傳統的EMD只適用于處理實數據序列，而雷達信號處理領域需要處理的都是復數據。因此，有多位學者把EMD擴展到復數域[13-15]，從而使EMD能應用于雷達信號處理領域[16-18]。本文采用文獻[15]的第2種CEMD算法來處理天波超視距雷達信號，提出了一種基于CEMD的天波超視距雷達瞬態干擾抑制方法。該方法首先利用CEMD把雷達回波信號分解成多個IMF，然后由第1個IMF來確定瞬態干擾的位置，再利用相關系數來判斷挖除數據的IMF分量個數，最后進行數據重構。實測數據處理結果表明該方法可以有效地抑制瞬態干擾。

2 復數據經驗模式分解(CEMD)

EMD是一種基于信號局部極值特征的自適應信號處理方法，它將任一時間信號分解成一系列的IMF，每個IMF同時滿足兩個條件：(1)信號極值的數目和零點的數目必須相等或最多相差一個；(2)在任意一點處，極大值包絡和極小值包絡的平均值為零。EMD算法可參考文獻[12]。傳統的EMD只適用于處理實數據序列，Rilling等人[15]把它擴展到復數域，它的基本思想是：信號是由快速旋轉分量(fast rotations)層疊在慢速旋轉分量(slower rotations)之上組成的。下面簡要描述其中一種算法：假設原始復信號為x(t)，

(1)初始化：r0=x(t),imf0=x(t),i=1；

(2)抽取第i個固有模態函數(IMF)：

(a)初始化：imfi,0=imfj-1,j=1；

(b)給定一系列方向

對imfi,j-1在方向φk上投影：

(c)抽取pφk(t)的極大值里插值得到方向φk的切線(t)；

(d)計算所有方向上的切線均值：

(e)求得分量

(f)如果滿足Rilling等人[19]提出的中止條件，則imfi=imfi,j。否則j=j+1轉到步驟(b)；

(3)令ri=ri-1-imfi；

(4)如果ri的極值點個數多于或等于η,η是預先給定的殘差包含極值點數的門限(一般η取值為3)。設i=i+1，轉到步驟(2)；否則，轉到步驟(5)；

3 基于CEMD的瞬態干擾抑制方法

3.1 雷達回波模型及其CEMD分解

CEMD分解實質上是一組濾波過程，它將信號中不同尺度的波動或趨勢逐級分解開來，產生一系列具有不同頻率尺度的IMF。每個IMF體現了蘊含在原信號中不同頻率尺度的振蕩特性。它是從高頻到低頻逐次分解的基于信號本身的自動的、自適應時變濾波過程。含高頻成分的IMF分量總是最先被分解出來，即局部特征時間尺度短的波動總是最先分解出來。與傳統方法相比,CEMD時空濾波的一個顯著的優點就是濾波后的結果能夠充分保留原始信號本身所固有的非平穩特征。天波超視距雷達回波信號是典型的非平穩信號，適合用CEMD來分析。假設x(t)是某一距離-方位分辨單元在一個相干積累時間內收到的時域回波信號，則x(t)可以用一個線性組合信號來表示。

其中c(t)是海(地)雜波，s(t)是目標信號，i(t)是瞬態干擾，in(t)表示其它干擾和噪聲。在雷達回波信號中，瞬態干擾、雜波和目標回波在時域以及頻域上有非常明顯的差異：一方面在時域上瞬態干擾能量集中，持續時間短，而雜波和目標信號能量分散，持續時間長。即瞬態干擾的局部特征時間尺度比雜波和目標信號的短，可認為是快速旋轉分量，先被CEMD抽取出來。另一方面在頻域上瞬態干擾能量分散，占有很寬的譜帶，有時甚至會占滿整個頻譜，而雜波和目標信號則能量集中，并且海(地)雜波的多普勒頻譜的主譜線在零頻率附近。因CEMD的濾波特點是從高頻到低頻逐次分解，所以可認為海(地)雜波是慢速旋轉分量，它們在瞬態干擾后被抽取出來。對實測的OTHR數據采用CEMD分解某一距離-方位分辨單元的回波信號，如圖1所示。圖中畫出的是回波信號、各個IMF分量(用imf1,imf2,…,imfq表示)及剩余分量res的實部波形，幅度均進行了歸一化處理。從圖中發現原始數據中明顯含有兩處瞬態干擾，其中瞬態干擾被分解到前面3個分量。這表明瞬態干擾確實是先被分解出來，并且很顯然第1個分量imf1比其它分量更能準確地反映瞬態干擾出現的位置。

3.2 基于CEMD的瞬態干擾抑制方法

觀察圖1中的雷達回波信號經CEMD分解后第1個IMF分量imf1的波形特征，發現在瞬態干擾出現的時段，波形變化范圍大，而其它時段變化比較平緩。因此，可用一滑窗在imf1上滑動，通過比較滑窗內數據的局部方差和imf1的方差來判斷該時段內是否有瞬態干擾。如果局部方差大于imf1的方差，則認為該時段內出現了瞬態干擾。瞬態干擾的位置一旦確定，就可以進行數據挖除和信號重構。為此，提出了一種基于CEMD的瞬態干擾抑制方法，具體算法如下：

圖1 CEMD分解某一OTHR實測數據的實部波形

(1)對某一個距離單元的回波 信號X={x(1),x(2),…,x(N)},N為數據的采樣點數，對其進行CEMD分解得

其中殘余量res代表原始數據序列的一種趨勢。

(2)確定干擾位置 設窗函數

其中L為窗函數的長度。計算imf1的標準方差σ，把σ作為檢測瞬態干擾的門限。然后，讓矩形窗在第1個分量信號imf1上滑動，計算每一個窗內所含數據的標準方差。當＞σ時就判斷該段回波數據含有瞬態干擾，這樣就確定了干擾出現的位置，同時把所有干擾位置記錄下來。

(3)抑制瞬態干擾 計算X與各個imfk(k=1,2,…,q)之間的相關系數ρk,k=1,2,…,q，取最大相關系數對應的信號分量，不妨設為imfh,1≤h≤q。把前h-1個信號分量imfj(j=1,2,…,h-1)對應于干擾出現的位置處的數據置為零，記為(j=1,2,…,h-1)。

(4)重構信號 把步驟(3)中處理后的各信號分量與其它信號分量相加即可得到去除瞬態干擾的重構信號

步驟(2)中，干擾的位置確定后，也可挖除原始回波信號對應干擾時段的數據來抑制瞬態干擾，然后通過AR模型預測恢復挖除的數據。這里不采用這種方法，主要是因為用AR模型來重構數據不僅增加了計算量，還需精心選擇模型的階數。另外，如果原始數據包含多處瞬態干擾，那么用來重構信號的數據比較少，重構效果不好，影響雷達目標檢測性能。步驟(3)通過計算相關系數來確定信號的主分量，這是考慮到天波超視距雷達的回波信號中海(地)雜波能量最強，是回波信號的主要部分。因此，可計算回波信號與各IMF分量之間的相關系數來判定哪個IMF分量是主雜波。由于CEMD是從高頻到低頻分解信號，主雜波以后的信號分量可認為都是雜波分量。

4 實測天波雷達數據處理

對實測的OTHR數據，采用基于CEMD的干擾抑制方法處理其某一距離-方位分辨單元的各次回波，發現該數據中含有瞬態干擾。CEMD分解結果如圖1所示，實驗中人為添加了一個歸一化頻率為0.1436，信雜比為-40dB的復正弦弱信號目標。窗函數w(n)的長度L=16，數據長度N=1024。圖2(a)是瞬態干擾抑制前原始信號的時頻圖，可以明顯看到有兩處瞬態干擾，其強度比雜波強度弱一些。還有一處不是很明顯，在靠近第1個采樣點附近。圖2(b)是瞬態干擾抑制后重構信號的時頻圖，表明基于CEMD的干擾抑制方法處理后，瞬態干擾完全得到了抑制。圖2(c)是抑制掉的瞬態干擾信號的時頻圖，從圖中可知能準確地檢測到瞬態干擾信號并確定它們的位置。該實測數據包含了3處瞬態干擾信號，一處大約在第130到第145個采樣點之間，另一處大約在第625個到第655個采樣點之間，還有一處在靠近第1個采樣點附近。這也可從圖2(a)和2(d)看出來。圖2(d)是瞬態干擾抑制前、后原始數據和重構數據的實部波形對比圖(幅度進行了歸一化處理)，可見信號突變部分變得很平滑。

圖2(e)是干擾抑制前、后原始信號和重構信號的多普勒譜對比圖，顯見重構信號雜波兩邊的譜比原始信號的低，并且目標突顯出來了。圖2(f)是多普勒改善效果圖，用原始信號多普勒譜減去重構信號多普勒譜而得到的。多普勒譜基底最大可下降約37 dB，平均下降約8.74 dB。最后，在實測數據確定了干擾位置對其挖除后，嘗試用AR模型來恢復挖除的數據。這時數據重構的效果跟AR模型的階數選擇有關，因此要有一個判斷準則來估計模型的階數，比如可采用由Rissanen提出的最小描述長度(MDL)信息量準則[20]來選擇模型階數，使其一步預測的平均誤差最小。圖3(a)是基于CEMD的干擾抑制方法處理結果，可見目標明顯突出來了。圖3(b)是采用AR模型處理結果，目標沒明顯突出。綜合圖2和圖3說明基于CEMD的抑制方法能更有效地抑制瞬態干擾。

為進一步驗證算法的有效性，對另一距離-方位單元含有應答信號的回波數據處理后發現，原來被掩蓋的應答信號突顯出來了，如圖4所示。以上實驗結果說明基于 CEMD的干擾抑制方法可有效抑制瞬態干擾，能極大地提高雷達目標檢測性能。

圖2 瞬態干擾抑制示意圖

5 結束語

雷達回波信號中瞬態干擾有明顯不同于雜波和目標信號的時頻特征，根據這一特點，本文提出了一種基于CEMD的天波超視距雷達瞬態干擾抑制方法，它在CEMD分解中同時完成瞬態干擾的檢測和數據重構。試驗結果表明，基于CEMD的天波超視距雷達瞬態干擾抑制方法可準確地檢測出瞬態干擾的位置并可有效地恢復數據，抑制瞬態干擾的效果好。該方法避免了雜波抑制和預測重構，運算量小，實用性強，有很高的工程應用價值。特別是它對相干積累時間內瞬態干擾的數量沒有特別要求，即使瞬態干擾的數目較多，仍然可以達到較好的干擾抑制和信號重構效果。

圖3 AR模型與基于CEMD方法抑制效果對比圖

圖4 含有目標信號瞬態干擾抑制效果圖

[1]Gill E,et al..The effect of bistatic scattering angle on the high frequency radar cross sections of the ocean surface[J].IEEE Geoscience and Remote Sensing Letters,2008,5(2)：143-146.

[2]Fabrizio G A.Over the horizon radar[C].IEEE Radar Conference,Rome,May 26-30,2008：1-2.

[3]Fabrizio G A,Gershman A B,and Turley M D.Robust adaptive beamforming for HF surface wave over-the-horizon radar[J].IEEE Transactions on Aerospace and Electronic Systems,2004,40(4)：510-525.

[4]Barnum J R and Simpson E E.Over-the-horizon radar sensitivity enhancement by impulsive noise suppressing[C].IEEE National Radar Conference,Syracuse,New York,USA,May 13-15,1997：252-256.

[5]Wei Min,Li Jun,Liu Tao,et al..Meteor trail interference model in HF environment[C].2006 International Conference on Communication,Circuits and Systems Proceedings,Guilin,China,June 25,2006：626-628.

[6]Liu Tao and Wang Jie.OTHR impulsive interference detection based on AR model in phase domain[J].WSEAS Transactions on Signal Processing,2009,5(4)：147-156.

[7]Liu Tao,Chen Xiao-xu,Wang Jie,et al..Subspace impulsive interference suppression in OTHR[J].Progress in Electromagnetics Research,2009,C(7)：167-181.

[8]邢孟道,保錚,強勇.天波超視距雷達瞬態干擾抑制[J].電子學報,2002,30(6)：823-826.Xing Meng-dao,Bao Zheng,and Qiang Yong.Transient interference excision in OTHR[J].Acta Electronica Sinica,2002,30(6)：823-826.

[9]黃亮,文必洋,鄧巍.高頻地波雷達抑制瞬態干擾研究[J].電波科學學報,2004,19(2)：166-170.Huang Liang,Wen Bi-yang,and Deng Wei.Suppressing instantaneous interference of high frequency ground wave radar[J].Chinese Journal of Radar Science,2004,19(2)：l66-170.

[10]陳希信,黃銀河.基于矩陣奇異值分解的高頻雷達瞬態干擾抑制[J].電子與信息學報,2005,27(12)：1879-1882.Chen Xi-xin and Huang Yin-he.A SVD-based approach of suppressing transient interference in High-Frequency radar[J].Journal of Electronics＆Information Technology,2005,27(12)：1879-1882.

[11]王階,劉濤,陳曉旭,等.天波超視距雷達流星余跡瞬態干擾抑制方法[J].電子測量與儀器學報,2009,23(10)：67-71.Wang Jie,Liu Tao,Chen Xiao-xu,et al..Method of meteor trail impulsive interference suppression in OTHR[J].Chinese Journal of Radio Science,2009,23(10)：67-71.

[12]Huang N E,Shen Z,Long S R,et al..The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings of the Royal Society A,1998,454(1971)：903-995.

[13]Tanaka T and Mandic D P.Complex empirical mode decomposition[J].IEEE Signal ProcessingLetters,2007,14(2)：101-104.

[14]Altaf M U,Gautama T,Tanaka T,et al..Rotation invariant empirical mode decomposition[C].IEEE International Conference Acoustics,Speech,Signal Processing,Honolulu,HI,2007,3：1009-1012.

[15]Rilling G,Flandrin P,Gon?alves P,et al..Bivariate empirical mode decomposition[J].IEEE Signal Processing Letters,2007,14(12)：936-939.

[16]Elgamel S A and Soraghan J J.Enhanced monopulse radar tracking using empirical mode decomposition[C].Proceedings of the 7th European Radar Conference,September 30-October 1,Paris,France,2010：57-60.

[17]王宏,Narayanan R M,周正歐,等.基于改進EEMD的穿墻雷達動目標微多普勒特性分析[J].電子與信息學報,2010,32(6)：1355-1360.Wang Hong,Narayanan R M,Zhou Zheng-ou,et al..Micro-Doppler character analysis of moving objects using through-wall radar based on improved EEMD[J].Journal of Electronics＆Information Technology,2010,32(6)：1355-1360.

[18]余志斌,金煒東.多分量 LFM 雷達輻射源信號的經驗模式分解[J].西南交通大學學報,2009,44(1)：49-54.Yu Zhi-bin and Jin Wei-dong.EMD for multi-component LFM radar emitter signals[J].Journal of Southwest Jiaotong University,2009,44(1)：49-54.

[19]Rilling G,Flandrin P,and Gon?alves P.On empirical mode decomposition and its algorithms[C].IEEE-EURASIP Workshop on Nonlinear Signal Image Process (NSIP 2003),June,Grado,Italy,2003：8-11.

[20]Rissanen J.Universal Prior for the integers and estimation by minimum description length[J].Annals of Statistics,1983,11(2)：416-431.