CN106154236B - 一种基于cemd抑制高频地波雷达射频干扰的方法 - Google Patents

Abstract

本发明涉及雷达信号处理领域，特别涉及一种基于CEMD抑制高频地波雷达射频干扰的方法。本发明在原有高频地波雷达信号处理的基础上，通过利用中位数设置阈值，对高于阈值的原始信号判定存在射频干扰，将原始信号搬移至高频带，利用复经验模态分解算法将原始信号分解出多个本征模态函数，在每个本征模态函数中通过固定阈值法检测定位干扰，并对干扰进行置零处理，对处理后的原始信号搬移后原频，得到重构信号。本发明可同时抑制瞬态和非瞬态干扰，在实际数据处理中，应用效果很好，在降低干扰的同时无损回波信号，抑制效果显著，增加了雷达探测距离，提高了雷达探测精度。

Description

一种基于CEMD抑制高频地波雷达射频干扰的方法

技术领域

本发明属于雷达信号处理领域，特别涉及到高频地波雷达射频干扰抑制的方法，具体是一种基于CEMD抑制高频地波雷达射频干扰的方法。

背景技术

高频地波雷达是一种可以进行大面积海洋环境监测的重要探测设备，并且由于其可全天候、大范围、不受雨雾影响地连续数月甚至数年工作的特点，因而日益受到人们的青睐。我国海岸线全长超过3.2万公里，因此高频地波雷达对国民经济和国防建设具有十分重要的意义。

高频地波雷达一般工作在3-30MHz，而在此频段会有大量射频干扰，射频干扰主要来自于雷达工作频段内密集的短波通讯信号、广播电台信号以及工业干扰。由于这些干扰均是有源的，其功率与雷达回波信号相比十分强大，当其进入接收机后，会极大地降低雷达的数据质量。特别是在傍晚17时以后，由于电离层D层的消失，射频干扰纷纷进入雷达接收机，严重影响系统工作性能，甚至使雷达无法正常工作。

依据持续时间长短，射频干扰可分为瞬态干扰和非瞬态干扰。对于瞬态干扰，主要利用时域算法抑制，通过小波变换、特征分解等技术检测射频干扰，并对干扰位置置零，最后使用数据插值或预测方法重构丢失段数据。但这些算法在抑制非瞬态干扰时性能不佳。针对非瞬态干扰，主要采用自适应波束形成算法。但该算法要求射频干扰保持空间平稳性，这在整个相干积累时间内难以满足，并且该算法无法抑制主瓣干扰。针对雷达信号在整个相干积累时间内同时受瞬态干扰和非瞬态干扰影响的现状，目前尚没有一种很好的算法可以同时抑制上述两种干扰。

发明内容

本发明针对上述问题，提出一种基于高频地波雷达的同时抑制瞬态干扰和非瞬态干扰的方法，避免了现有算法的不利影响。本发明在抑制干扰的同时完整保留回波信号，为现有的高频雷达系统提供一种快速有效的抗射频干扰方法，从而能够增加雷达探测距离，提高雷达探测精度。

本发明的技术方案是一种基于CEMD抑制高频地波雷达射频干扰的方法，包括以下步骤：

步骤1，检测射频干扰。求出原始信号的中位数幅度A_median，设置k*A_median为阈值；k为常数，一般取10-20dB；由于射频干扰功率远高于回波信号，因此原始信号的平均值不能真实反映回波信号平均幅度，用中位数替代，设置数倍中位数为阈值；若原始信号不高于阈值，则判定无干扰，不进行后续处理；否则存在干扰，转到步骤2。

步骤2，对原始信号乘以一个给定的高频信号hf(t)，将信号搬移至高频带。这个步骤可以表达为：

x₁(t)＝x(t)×hf(t) (1)

式中，x₁(t)表示搬移后的原始信号，x(t)为原始信号。

步骤3，对搬移后的信号x₁(t)进行复经验模态分解CEMD(Complex EmpiricalMode Decomposition)得到多个本征模态函数IMF(Intrinsic Mode Function)和一个趋势函数r(t)，如下式：

式中，IMF_p(t)表示第p个本征模态函数。

步骤4，选择第一个IMF，求出该IMF的中位数幅度I_median，设置q*I_median为阈值，检测是否有高于该阈值的部分，若存在，对高于阈值的位置进行定位。接着对该位置进行置零处理。

步骤5，重复步骤4，直至循环完所有IMF和趋势函数。对经过处理后的所有IMF和趋势函数相加，得到抑制干扰后的原始信号。

步骤6，将抑制干扰后的原始信号搬移回原频带，得到抑制射频干扰后的重构信号。

所述步骤1和步骤4中，采用固定阈值法检测射频干扰，由于射频干扰功率远高于回波信号，因此原始信号的平均值不能真实反映回波信号平均幅度，用中位数替代，设置数倍中位数为阈值，对不高于阈值的原始信号不进行后续处理，并且针对射频干扰被分解到多个本征模态函数的特点，相比于步骤1，在步骤4中设置了一个较低的阈值。该方法简单快捷，节省了算法时间，满足雷达毫秒级的实时处理要求。

所述步骤2中，对原始信号乘以高频信号，将其搬移至高频带。在大多数的传统算法和本文的后续步骤中，均采用对射频干扰直接置零来抑制干扰，该处理会在低频段带来新的噪声，而低频带正是回波信号聚集的频段。在步骤2中先将原始信号搬移至高频带，完成射频干扰抑制后再将重构信号搬移回原频带，噪声则被搬移到了负频，避免了新产生的低频噪声对回波信号的影响。

本发明的优点在于：

1.在高频地波雷达中采用CEMD算法进行射频干扰抑制处理，在处理过程中无需区分瞬态干扰和非瞬态干扰，可同时抑制这两种射频干扰。

2.CEMD算法相比于传统的EMD算法，可以适用于复信号处理，在高频地波雷达中，原始信号的实部、虚部存在相关性，因此采用CEMD算法。通过CEMD将原始信号按频率高低分解出不同的IMF，对每个IMF进行处理，从而在抑制射频干扰的同时，最大程度地保留了有用信号特征。

3.采用数倍中位数作为阈值，原始信号中包含了高功率的射频干扰，此时的平均值不能真实反映回波信号的平均幅度，因此替代为中位数。并且本算法首先检测是否含有射频干扰，对不含有射频干扰的原始信号不进行处理，节省了雷达实时处理时间。

4.将原始信号搬移至高频带，避免了在IMF中置零带来新的噪声。若原始信号不进行频谱搬移，直接在IMF中对干扰位置进行置零处理，会给重构信号低频段带来新的噪声，而低频段是回波信号所处的位置。本发明首先采用频谱搬移的方式将原始信号搬移至高频带，接着将重构后的信号搬回原频，避免了低频噪声的影响。

5.对IMF中干扰位置进行置零处理，不会损失其他IMF中该位置的回波信号，达到了抑制干扰的同时无损回波信号，而且该处理速度快捷，满足雷达实时工作要求。

本发明通过以上创新点提出了一种新型高频地波雷达射频干扰抑制方法。本发明所述方法基于射频干扰和回波信号的功率、频率差异性，进而抑制射频干扰，因此可同时抑制瞬态和非瞬态射频干扰，雷达信号在相干积累时间内通常被这两种干扰同时影响，因此本发明所述方法对雷达信号的实际应用具有重要价值。此方法在实际数据处理中应用效果很好。在降低干扰的同时无损回波信号，抑制效果显著，增加了雷达探测距离，提高了雷达探测精度。

附图说明

图1，高频地波雷达的工作原理框图；

图2，射频干扰经系统处理的过程；

其中，图(a)为射频干扰与本振信号混频图；图(b)为射频干扰经过低通滤波器图；

B—扫频带宽，T—扫频周期，t₁—干扰起始时间，t₂—干扰起始时间，f_i—干扰频率，f₀—工作频率，b—滤波器带宽；

图3，本发明的算法流程图；

图4，含有瞬态干扰的原始信号抑制前后对比(实线代表原始信号，点虚线表示重构信号)；

图5，瞬态干扰抑制前后的距离多普勒谱；其中，图5(a)为原始的距离多普勒；图5(b)为抑制瞬态干扰后的距离多普勒谱；

图6，含有非瞬态干扰的原始信号抑制前后对比(实线代表原始信号，点虚线表示重构信号)；

图7，瞬态干扰抑制前后的距离多普勒谱，其中，图7(a)为原始的距离多普勒；图7(b)为抑制瞬态干扰后的距离多普勒谱。

具体实施方式

下面以具体实例结合附图，对本发明作进一步说明。

线性调频体制的高频地波雷达的工作原理及射频干扰的特征介绍如下：

高频地波雷达的工作原理框图如图1所示(其中1为接收天线，2为线性扫频的本振信号，3为低通滤波，4为距离解压缩变换，5为多普勒变换)。雷达采用线性调频波形体制，在每个扫频周期内，回波信号经混频、低通滤波、采样和快时域的离散傅里叶变换得到该扫频周期的距离谱。距离谱中的每个谱点对应于一个距离元的采样点，在多个扫频周期的相干积累时间内，对多个距离谱的采样序列进行慢时域的离散傅里叶变换得到该距离元的多普勒谱。对探测范围内的所有距离元都进行慢时域的离散傅里叶变换，则得到距离-多普勒二维谱。

线性调频本振信号可表示为

其中，f₀为雷达载频，K＝B/T为扫频速率，B为扫频带宽，T为扫频时宽。

不失一般性，考虑射频干扰为单频信号，

I(t)＝A(t)exp(j2πf_it) (4)

其中，A(t)为射频干扰的幅度，f_i为射频干扰的频率。

经混频、低通滤波后输出的干扰信号为

式中，Δf＝f_i-f₀为射频干扰与雷达工作频率的频差。

从上式可以看出，一个单频干扰经解调、滤波后的输出是一个带宽受限的线性调频信号。该信号的带宽为低通滤波器的带宽b。

图2说明了频率为f_i的射频干扰经系统处理的过程。与雷达工作频率同频段的射频干扰与线性调频本振信号混频后将产生新的线性调频信号，经过低通滤波器后在时域上明显具有短时性，其持续时间而在频域(即距离谱)上极大地增加了基底噪声，通常会淹没Bragg峰，从而严重影响雷达的探测性能。

本发明正是利用射频干扰解调后时间短、强度大的特性，对原始信号进行复经验模态分解得到多个IMF，并分别进行干扰抑制，最大程度地保留了回波信号。具体实施步骤参见图3：

步骤1，检测射频干扰。求出原始信号的中位数幅度A_median，设置k*A_median(k为常数，一般取10-20dB)为阈值。若原始信号不高于阈值，则判定无干扰，不进行后续处理；否则存在干扰，转到步骤2。

步骤2，对原始信号乘以一个给定的高频信号hf(t)，将信号搬移至高频带，高频信号的频率通常为一半的采样频率。这个步骤可以表达为：

x₁(t)＝x(t)×hf(t) (6)

式中，x₁(t)表示搬移后的原始信号，x(t)为原始信号。

步骤3，对搬移后的信号x₁(t)进行复经验模态分解CEMD(Complex EmpiricalMode Decomposition)得到多个本征模态函数IMF(Intrinsic Mode Function)和一个趋势函数r(t)，如下式：

式中，IMF_p(t)表示第p个本征模态函数。有关CEMD的分解步骤可参见相关文献。

步骤4，选择第一个IMF，求出该IMF的中位数幅度I_median，设置q*I_median(q为常数，取8-12dB)为阈值，检测是否存在高于阈值的部分，若不存在，跳过该IMF。否则对高于阈值的位置进行定位，接着对该位置进行置零处理。

步骤5，重复步骤4，直至循环完所有IMF和趋势函数。对经过处理后的所有IMF和趋势函数相加，得到抑制干扰后的原始信号。

步骤6，将抑制干扰后的原始信号搬移回原频带，这与步骤2为正反变换。最终，得到抑制射频干扰后的重构信号。

图4和图5为对瞬态干扰进行抑制前后的原始信号、距离多普勒谱对比，瞬态干扰得到了极大抑制，信噪比得到提升。

图6和图7为对非瞬态干扰进行抑制前后的原始信号、距离多普勒谱对比，非瞬态干扰被显著抑制，信噪比有所提升。

Claims (1)

1.一种基于CEMD抑制高频地波雷达射频干扰的方法，其特征在于：包括以下步骤：

步骤1，检测射频干扰：求出原始信号的中位数幅度A_median，设置k*A_median为阈值，k为常数，其范围取10-20dB；若原始信号的幅度中存在高于阈值的部分，则存在干扰，转到步骤2；

步骤2，对原始信号乘以一个给定的高频信号hf(t)，将信号搬移至高频带，高频信号的频率为采样频率的二分之一；此步骤表达为：

x₁(t)＝x(t)×hf(t) (1)

式中，x₁(t)表示搬移后的原始信号，x(t)为原始信号；

步骤3，对搬移后的信号x₁(t)进行复经验模态分解CEMD得到多个本征模态函数IMF和一个趋势函数r(t)，如下式：

式中，IMF_p(t)表示第p个本征模态函数；

步骤4，选择第一个IMF，求出该IMF的中位数幅度I_median，设置q*I_median为阈值，q为常数，取8-12dB；检测是否存在高于阈值的部分，若不存在，跳过该IMF；否则对高于阈值的位置进行定位，接着对该位置进行置零处理；

步骤5，重复步骤4，直至循环完所有IMF；对经过处理后的所有IMF和趋势函数相加，得到抑制干扰后的原始信号；

步骤6，将抑制干扰后的原始信号搬移回原频带，这与步骤2为正反变换；最终，得到抑制射频干扰后的重构信号。