基于复数据经验模式分解的天波超视距雷达瞬态干扰抑制

周忠根 水鹏朗

(西安电子科技大学雷达信号处理国家重点实验室 西安 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.