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

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

0 引 言

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

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

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

1 希爾伯特黃變換

HHT包括2個階段:第一階段為 EMD,將信號按照一定規(guī)則自適應(yīng)地分解成若干IMF分量;第二階段為Hilbert譜分析,即獲得每個IMF分量的瞬時(shí)頻率和能量隨時(shí)間的分布。CEMD方法借助了“旋轉(zhuǎn)”這一概念,將復(fù)數(shù)信號描述為快速旋轉(zhuǎn)IMF分量和慢速旋轉(zhuǎn)IMF分量的疊加,通過計(jì)算復(fù)信號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ò)的均值作為復(fù)信號包絡(luò)(即慢速旋轉(zhuǎn)分量),通過多次迭代從而將快速旋轉(zhuǎn)分量提取出來。CEMD算法的詳細(xì)流程如圖1所示,其中復(fù)數(shù)IMF分量需滿足2個條件:①在任意旋轉(zhuǎn)方向φn上,復(fù)信號投影分量實(shí)部的零點(diǎn)與極值點(diǎn)相等或至多相差一個;②任意時(shí)刻,復(fù)信號在所有旋轉(zhuǎn)方向上的投影分量確定的切線均值為零。

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

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

(3)

其中

(4)

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

(5)

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

2 機(jī)動目標(biāo)模型及參數(shù)估計(jì)方法

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

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

(6)

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

(7)

式中,f0為雷達(dá)載頻；θ為電磁波相對于海平面的入射角。機(jī)動目標(biāo)速度模型可用n階多項(xiàng)式表示為

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

(8)

式中,v0為目標(biāo)初始速度;an為n階加速度。雷達(dá)機(jī)動目標(biāo)檢測中通常只需補(bǔ)償一階加速度[6],此時(shí)回波信號模型可表示為

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

(9)

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

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

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

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

(10)

步驟4計(jì)算IMF分量的瞬時(shí)頻率和幅度,獲得目標(biāo)HHT的時(shí)頻譜,采用最小二乘法對瞬時(shí)頻率隨時(shí)間的變化曲線進(jìn)行多項(xiàng)式擬合,獲得目標(biāo)初速度和加速度的估計(jì)。

3 算法性能分析

3.1 單目標(biāo)場景

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

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

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

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

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

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

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

圖4 單目標(biāo)HHT時(shí)頻譜Fig.4 HHT Spectrum of the single target

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

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

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

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

3.2 多目標(biāo)場景

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

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

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

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

圖8 目標(biāo)1的HHT時(shí)頻譜Fig.8 HHT spectrum of target 1

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

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

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

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

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

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

4 結(jié) 論

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

[1] 周萬幸.天波超視距雷達(dá)發(fā)展綜述[J].電子學(xué)報(bào),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ìn)的通道補(bǔ)償法檢測高頻雷達(dá)機(jī)動目標(biāo)[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ìn)峰,李萬閣,艾慧,等.基于改進(jìn)時(shí)頻分析方法的高頻雷達(dá)機(jī)動目標(biāo)檢測算法研究[J].電子與信息學(xué)報(bào),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]. 電波科學(xué)學(xué)報(bào), 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)驗(yàn)?zāi)B(tài)分解邊緣效應(yīng)抑制方法綜述[J]. 儀器儀表學(xué)報(bào), 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.