基于HHT的高頻雷達機動目標參數(shù)估計

張冰瑞, 雷志勇, 周海峰, 周 毅

(中國電子科技集團公司第十四研究所, 江蘇 南京 210039)

0 引 言

高頻雷達工作頻段在3～30 MHz,具有超視距、覆蓋范圍大、性價比高的獨特優(yōu)勢,在反走私稽查、電離層環(huán)境監(jiān)測、早期戰(zhàn)略預警、隱形飛機探測等多個領(lǐng)域具有重要的民用和軍事價值[1-3]。高頻雷達通過長時間相干積累增加目標信噪比,對于高速機動的飛機或?qū)椖繕?多次回波中其運動狀態(tài)變化較大,導致目標多普勒擴展影響能量聚集[4-6]。對這類機動目標的檢測,需要通過一定手段進行運動補償來提高檢測效率,如何估計目標運動參數(shù)是其中的關(guān)鍵技術(shù),逐漸成為該領(lǐng)域的研究熱點。

目前的研究中,機動目標參數(shù)估計主要根據(jù)回波信號的相位信息,通過建立相應數(shù)學模型,利用不同數(shù)學工具對其進行解析和估計,從而獲得目標速度和加速度等信息。常用算法有多通道補償法[7]、多項式相位法[8-10]和基于時頻分析技術(shù)[11-13]等算法。多通道補償法將目標近似為勻加速運動,取若干加速度值分別對回波進行相位補償,根據(jù)不同補償通道的信噪比改善效果估計出目標加速度。文獻[7]提出的多通道補償算法綜合應用了頻率校正、通道合成和峰值檢測,具有計算復雜度和虛警率低的優(yōu)勢。然而此類算法近似為窮舉法,受實際通道數(shù)的限制難以實現(xiàn)目標運動參數(shù)的高精度估計。多項式相位法利用多項式為目標相位建模,通過高階模糊函數(shù)估計出多項式的各階系數(shù),進而獲得運動參數(shù)估計。文獻[8]提出的方法通過高階模糊函數(shù)估計出目標的高階運動參數(shù),并在回波中補償?shù)粼撾A分量,循環(huán)執(zhí)行直到目標所有運動分量都補償?shù)?以此來估計各階運動參數(shù)。該類方法計算量小、精度高,但對信噪比要求較高,在強雜波和多目標的情況下穩(wěn)定性差。基于時頻分析的方法將目標信號變換到時頻域,通過時頻域的能量積累進行目標參數(shù)估計。文獻[11]提出的時頻分析技術(shù)通過構(gòu)造時間-頻率變化率的聯(lián)合域,并沿時間軸的不同直線路徑進行積分來估計目標參數(shù),在低信噪比下取得了較好效果。由于高頻雷達地海雜波能量較大(通常比目標大10～50 dB),常規(guī)時頻分析方法均需先利用海雜波循環(huán)對消法[14]和奇異值分解法[15-16]等抑制雜波,不僅增加了算法的運算量而且參數(shù)估計精度受制于雜波抑制效果的影響。同時,高頻雷達距離分辨率低,經(jīng)常會遇到同一距離單元存在多個目標的場景,常規(guī)算法需要逐個目標進行參數(shù)估計、加速度補償和目標剔除的循環(huán)操作,導致運算量大且算法穩(wěn)定性差。

針對現(xiàn)有算法存在的問題,本文研究基于希爾伯特黃變換(Hilbert-Huang transform,HHT)的機動目標運動參數(shù)估計方法,HHT是文獻[17]提出的時頻分析方法,該算法是數(shù)據(jù)驅(qū)動的自適應信號分析方法,能夠很好地處理非平穩(wěn)和非線性信號。經(jīng)驗模態(tài)分解(empirical mode decomposition,EMD)是其核心算法,通過將信號從時域分解為若干本征模態(tài)函數(shù)(intrinsic mode function,IMF)分量之和,每個IMF分量具有不同的時頻特性,并利用希爾伯特變換計算出每個IMF分量的瞬時頻率,從而獲得信號的時頻分布特性。為便于分析雷達信號,將EMD擴展到復數(shù)域,提出了復數(shù)經(jīng)驗模態(tài)分解(complex empirical mode decomposition, CEMD)方法[18-20]。本文將其用于高頻雷達的機動目標參數(shù)估計,運用CEMD將目標慢時間維信號分解為多個IMF分量,根據(jù)目標和雜波的時頻特性將其分解為不同的IMF分量,利用目標瞬時多普勒頻率的變化特性,通過線性擬合估計出目標運動參數(shù)。與常規(guī)機動目標參數(shù)估計方法相比,該方法無需抑制雜波,可同時對多目標進行參數(shù)估計,數(shù)據(jù)分析結(jié)果表明該方法可以有效應用于目標參數(shù)估計。

1 希爾伯特黃變換

HHT包括2個階段:第一階段為 EMD,將信號按照一定規(guī)則自適應地分解成若干IMF分量;第二階段為Hilbert譜分析,即獲得每個IMF分量的瞬時頻率和能量隨時間的分布。CEMD方法借助了“旋轉(zhuǎn)”這一概念,將復數(shù)信號描述為快速旋轉(zhuǎn)IMF分量和慢速旋轉(zhuǎn)IMF分量的疊加,通過計算復信號s(t)在不同旋轉(zhuǎn)方向φn上的投影s'(t)[18],即

s'(t)=Re[s(t)e-iφnt]

(1)

φn=2nπ/N,n=1,2,…,N

(2)

將各旋轉(zhuǎn)方向上投影包絡(luò)的均值作為復信號包絡(luò)(即慢速旋轉(zhuǎn)分量),通過多次迭代從而將快速旋轉(zhuǎn)分量提取出來。CEMD算法的詳細流程如圖1所示,其中復數(shù)IMF分量需滿足2個條件:①在任意旋轉(zhuǎn)方向φn上,復信號投影分量實部的零點與極值點相等或至多相差一個;②任意時刻,復信號在所有旋轉(zhuǎn)方向上的投影分量確定的切線均值為零。

圖1 CEMD算法流程圖Fig.1 Flowchart of CEMD algorithm

傳統(tǒng)傅里葉變換中,頻率的定義主要針對平穩(wěn)信號而言,即在時域上具有恒定幅度和頻率的正弦或余弦函數(shù)。頻率被定義為整個信號長度中具有恒定幅度的正弦或者余弦函數(shù)。而對于非平穩(wěn)信號,頻率是隨時間不斷變換的,因此用瞬時頻率來描述更為合適。復數(shù)信號經(jīng)過CEMD分解為若干單頻(或窄帶)瞬態(tài)信號A(t)ejω(t),瞬時頻率定義為

(3)

其中

(4)

此時原信號s(t)可表示為時間、頻率和幅度的函數(shù)，即若干IMF分量之和的形式為

(5)

式中，IMFi為第i個IMF分量；r為分解余量；fi(t)和Ai(t)分別為第i個IMF分量的瞬時頻率和幅度。將其表示在時頻平面上便得到s(t)的HHT時頻譜H(f,t)。

2 機動目標模型及參數(shù)估計方法

高頻雷達慢時間維回波信號s(t)的模型可表示為

s(t)=x(t)+c(t)+w(t)

(6)

式中,x(t)為目標信號;c(t)為地海雜波;w(t)為噪聲。當高頻雷達照射海上目標時,回波中就含有大量海雜波,回波頻譜中海雜波在一階Bragg峰fb處有很強的能量,其模型可以表示為[1]

(7)

式中,f0為雷達載頻；θ為電磁波相對于海平面的入射角。機動目標速度模型可用n階多項式表示為

v(t)=v0+a1t+…+antn

(8)

式中,v0為目標初始速度;an為n階加速度。雷達機動目標檢測中通常只需補償一階加速度[6],此時回波信號模型可表示為

x(t)=ej4π/λ(v0t+0.5a1t2)

(9)

高頻雷達中根據(jù)雜波和目標的不同頻譜特性,可采用CEMD算法將目標回波從雜波中分解出來,進而計算出目標時頻譜,根據(jù)目標瞬時多普勒頻率與時間的線性關(guān)系擬合出運動參數(shù)。具體步驟如下。

步驟1雷達在某一方位和距離單元上的慢時間維回波信號為S={s(1),s(2),…,s(N)},N為相干積累點數(shù)。采用CEMD對回波S進行分解,可得到若干IMF分量。

步驟2為抑制CEMD分解中的端點效應,即信號兩端極值點的不確定導致的信號包絡(luò)線擬合偏差,這里采用鏡像延拓法[21-22]以靠近端點處的極值點為對稱軸,將原始信號映射成一個周期信號,從而盡可能消除端點的影響。

步驟3合理的CEMD分解停止原則也是算法的關(guān)鍵。本文采用仿柯西收斂準則控制EMD的分解次數(shù),即兩個連續(xù)IMF的標準差SD小于一定閾值時判定分解過程結(jié)束[23]。

(10)

步驟4計算IMF分量的瞬時頻率和幅度,獲得目標HHT的時頻譜,采用最小二乘法對瞬時頻率隨時間的變化曲線進行多項式擬合,獲得目標初速度和加速度的估計。

3 算法性能分析

3.1 單目標場景

對高頻雷達某一距離單元的慢時間維數(shù)據(jù)進行HHT算法處理,距離-多普勒譜如圖2所示。

圖2 回波信號距離-多普勒譜Fig.2 Distance-doppler spectrum of the echo signal

由圖2可知，數(shù)據(jù)中加入一初始速度v0=100 m/s、加速度a1=30 m/s2的目標,目標信噪比20 dB,雜噪比47 dB。

信號時域圖和頻域圖如圖3所示,雖然目標能量較強,但由于目標加速度大,頻譜展寬嚴重,采用常規(guī)CFAR檢測困難。利用HHT算法對信號進行處理。首先對原始信號進行CEMD計算,形成2個IMF分量和1個剩余項，如圖3所示。

圖3 單目標CEMD計算結(jié)果Fig.3 CEMD results of the single target

由圖3(a)中信號實部的時域波形可知,IMF分量的波動速率從大到小逐步降低,對這些信號做FFT;由圖3(b)頻譜圖可知，原始回波是海雜波、目標和噪聲成分的疊加。CEMD處理后第1個IMF分量主要為目標分量,第2個IMF分量和剩余信號主要為海雜波分量,而噪聲存在于所有信號成分中,可見CEMD有效地將目標和雜波分離開來。

計算第1個IMF分量的瞬時頻率和幅度,繪制成HHT時頻譜如圖4所示。圖中用顏色變化表示瞬時頻率隨時間的分布情況,紅色曲線為采用最小二乘法線性擬合的結(jié)果。

圖4 單目標HHT時頻譜Fig.4 HHT Spectrum of the single target

圖5 單目標Dechirp計算結(jié)果Fig.5 Dechirp results of the single target

為研究算法魯棒性,分析信噪比大小對目標速度和加速度估計精度的影響,在100個不同距離單元回波中加入不同信噪比的目標信號,計算HHT算法對目標速度和加速度的估計誤差,并對100次計算結(jié)果取平均,獲得統(tǒng)計結(jié)果如圖6所示。

圖6 不同信噪比下單目標參數(shù)的估計誤差Fig.6 Estimation errors of single target parameters under different SNRs

由圖6可知，Dechirp算法受制于SVD抑制海雜波的性能,低信噪比下速度和加速估計誤差較大;采用HHT算法時,隨著目標信噪比的增大,速度和加速估計誤差呈快速減小的趨勢,當信噪比大于13 dB時,兩種方法的參數(shù)估計誤差均小于1%。

3.2 多目標場景

考慮到高頻雷達距離分辨率低,經(jīng)常會遇到同一距離單元存在多個目標的情況,下面對多目標場景下的算法性能進行分析。在回波數(shù)據(jù)中加入初始速度100 m/s、加速度30 m/s2的目標1,同時加入速度80 m/s、加速度0 m/s2的目標2,兩目標信噪比均為20 dB,對該場景下的數(shù)據(jù)進行HHT分析,結(jié)果如圖7所示。

圖7 多目標CEMD計算結(jié)果Fig.7 CEMD result of multi-targets

由圖7可知,原始信號為海雜波、目標1、目標2和噪聲的疊加,經(jīng)CEMD處理信號被分解為3個IMF分量和1個剩余項,第1個IMF分量主要為目標1回波,第2個IMF分量主要為目標2回波,第3個IMF分量和剩余項主要為雜波。

第1個IMF分量的HHT時頻譜如圖8所示。

圖8 目標1的HHT時頻譜Fig.8 HHT spectrum of target 1

嘗試Dechirp算法處理該數(shù)據(jù),結(jié)果如圖9所示。

圖9 多目標Dechirp計算結(jié)果Fig.9 Dechirp results under multi-objectives

由圖9可知，SVD將雜波剔除后,根據(jù)信號頻譜隨補償系數(shù)的變化獲得目標速度和加速度的估計值,分別為100.7 m/s和29.7 m/s2。

采用與第3.1節(jié)相同的方法,分別計算HHT和Dechirp兩種算法下,分析信噪比變化對目標1參數(shù)估計精度的影響(目標2信噪比始終保持為20 dB),結(jié)果如圖10所示。

圖10 不同信噪比下目標1的參數(shù)估計誤差Fig.10 Estimation errors of target 1 parameters under different SNRs

由圖10可知，估計誤差隨目標信噪比的變化趨勢與單目標場景下的統(tǒng)計結(jié)果相似。表明算法在多目標場景下具有較強魯棒性。

4 結(jié) 論

針對高頻雷達強雜波背景下的機動目標檢測問題,給出了基于HHT的機動目標參數(shù)估計算法。該算法利用CEMD從回波中分離出地海雜波和目標信號,進而計算目標不同時刻下的瞬時頻率獲得HHT時頻圖,最終通過線性擬合獲得目標運動參數(shù)的估計。數(shù)據(jù)分析結(jié)果表明,該算法可以有效從強雜波背景中分離出目標信號,直接對分離出的目標信號進行參數(shù)估計,從而避免了常規(guī)算法為消除地海雜波而采取的復雜處理過程;多目標場景下的處理結(jié)果顯示,該算法可以將不同目標信號分離開,進而同時對多目標進行運動參數(shù)估計,對后續(xù)檢測而言可同時對多目標進行速度補償和檢測,有助于簡化機動目標檢測的處理流程;對不同信噪比目標的運動參數(shù)估計結(jié)果顯示,隨著目標信噪比的增大,速度和加速估計誤差呈快速減小的趨勢,信噪比大于13 dB時估計誤差小于1%。

[1] 周萬幸.天波超視距雷達發(fā)展綜述[J].電子學報,2011,39(6):1373-1378.

ZHOU W X. An overview on development of skywave over-the-horizon radar[J].Acta Electronica Sinica,2011,39(6):1373-1378.

[2] FRAZER G J. Application of MIMO radar techniques to over-the-horizon radar[C]∥Proc.of the IEEE International Symposium on Phased Array Systems and Technology, 2016:1-3.

[3] LIU Z, SU H, HU Q. Transient interference suppression for high-frequency skywave over-the-horizon radar[J]. Iet Radar Sonar & Navigation, 2017, 10(8):1508-1515.

[4] XU J, ZHOU X, QIAN L C, et al. Hybrid integration for highly maneuvering radar target detection based on generalized radon-fourier transform[J]. IEEE Trans.on Aerospace & Electronic Systems, 2017, 52(5):2554-2561.

[5] QUAN Y, ZHANG L, LI Y, et al. OTHR spectrum reconstruction of maneuvering target with compressive sensing[J]. International Journal of Antennas & Propagation,2014,2014(10): 1-10.

[6] WANG G, XIA X G, ROOT B T, et al. Moving target detection in over-the-horizon radar using adaptive chirplet transform[J]. Radio Science, 2016, 38(4):1-24.

[7] 雷志勇, 黃銀河, 周海峰, 等. 改進的通道補償法檢測高頻雷達機動目標[J]. 系統(tǒng)工程與電子技術(shù), 2009, 31(4): 822-825.

LEI Z Y, HUANG Y H, ZHOU H F, et al. Improved channel compensation method for maneuvering target detection in OTH radar[J]. Systems Engineering and Electronics, 2009, 31(4): 822-825.

[8] LU K, LIU X Z. Enhanced visibility of maneuvering targets for high-frequency over-the-horizon radar[J]. IEEE Trans.on Antennas and Propagation, 2005, 53(1): 404-411.

[9] XIN Z, LIAO G, YANG Z, et al. A fast ground moving target focusing method based on first-order discrete polynomial-phase transform[J]. Digital Signal Processing, 2017, 60(1):287-295.

[10] KANG B, BAE J, LEE S, et al. Isar rotational motion compensation algorithm using polynomial phase transform[J]. Microwave & Optical Technology Letters,2016,58(7):1551-1557.

[11] 胡進峰,李萬閣,艾慧,等.基于改進時頻分析方法的高頻雷達機動目標檢測算法研究[J].電子與信息學報,2015,37(8):1843-1848.

HU J F, LI W G, AI H, et al. Maneuvering target detection algorithm based on improved Ttime-frequency analysis method in skywave radar[J]. Journal of Electronics & Information Technology, 2015, 37(8): 1843-1848.

[12] XIA X G. Discrete chirp-Fourier transform and its application to chirp rate estimation[J].IEEE Signal Process，2008,48(11): 3122-3133.

[13] SURESH P, THAYAPARAN T, VENKATARAMANIAH K. Fourier-Bessel transform and time-frequency-based approach for detecting manoeuvring air target in sea-clutter[J]. Iet Radar Sonar & Navigation, 2015, 9(5):481-491.

[14] 李雪, 郭曉彤, 王岳松,等. 基于已知傳播模式數(shù)目的海雜波抑制方法研究[J]. 電波科學學報, 2016, 31(4):700-705.

LI X, GUO X T, WANG Y S, et al. Sea clutter suppression algorithm based on ionospheric propagation mode number[J]. Chinese Journal of Radio Science, 2016, 31(4):700-705.

[15] CHEN Z, HE C, ZHAO C, et al. Using SVD-FRFT filtering to suppress first-order sea clutter in HFSWR[J]. IEEE Geoscience & Remote Sensing Letters, 2017, PP(99):1-5.

[16] 薄超, 顧紅, 蘇衛(wèi)民. 基于高階奇異值分解的 OTHR 海雜波抑制算法[J]. 系統(tǒng)工程與電子技術(shù), 2014, 36(5):872-878.

BO C, GU H, SU W M. OTHR sea clutter suppression algorithm based on higher order singular value decomposition[J]. Systems Engineering and Electronics, 2014, 36(5):872-878.

[17] 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.

[18] TANAKA T, MANDIC D P. Complex empirical mode decomposition[J]. IEEE Signal Processing Letters, 2007, 14(2): 101-104.

[19] ALTAF M U, GAUTAMA T, TANAKA T, et al. Rotation invariant empirical mode decomposition[C]∥Proc.of the IEEE International Conference Acoustics, Speech, Signal Processing, 2007, 3: 1009-1012.

[20] GüNTüRKüN U. Bivariate empirical mode decomposition for cognitive radar scene analysis[J]. IEEE Signal Processing Letters, 2015, 22(5):603-607.

[21] 杜陳艷,張榆鋒,楊平,等. 經(jīng)驗模態(tài)分解邊緣效應抑制方法綜述[J]. 儀器儀表學報, 2009, 30 (1): 55-60.

DU C Y, ZHANG Y F, YANG P, et al. Ap-proaches for the end effect restraint of empirical mode decom-position algorithm[J].Chinese Journal of Scientific Instrument,2009,30(1): 55-60.

[22] ZHAO J, ZHAO P, CHEN Y. Using an improved BEMD method to analyse the characteristic scale of aeromagnetic data in the gejiu region of Yunnan, China[J]. Computers & Geosciences, 2016, 88(7):132-141.

[23] ZHENG J, CHENG J, YANG Y. Partly ensemble empirical mode decomposition: an improved noise-assisted method for eliminating mode mixing[J].Signal Processing,2014,96(3):362-374.