CN102520396A - 基于复数经验模态分解的时变窄带干扰抑制方法 - Google Patents

Abstract

本发明公开了一种基于复数经验模态分解的时变窄带干扰抑制方法。其步骤为：(1)获取SAR回波(2)判断是否存在时变窄带干扰(3)复数经验模态CEMD分解(4)距离脉压及走动校正(5)参数估计(6)运动补偿(7)方位脉压(8)多视处理。本发明学习的窄带干扰抑制方法可适用于距离和方位二维时变的NBI。本发明充分利用复数经验模态分解方法，对阈值抑制强干扰分量进行设置，较好的抑制强干扰能力，并最大可能的保留了信号，SAR成像更加清晰。

Description

基于复数经验模态分解的时变窄带干扰抑制方法

技术领域

本发明属于信号处理技术领域，更进一步涉及雷达信号处理领域中的基于复数经验模态分解的时变窄带干扰抑制方法。本发明可用于抑制雷达回波信号中的窄带干扰，减少有用信号的损失，实现清晰的SAR成像。

背景技术

窄带干扰是合成孔径雷达(SAR)系统中一种常见的干扰形式。它主要包括射频干扰、自激干扰和人为释放干扰等。而射频干扰主要指SAR雷达(特别是低波段SAR雷达)受到同一波段的电视网、通讯网等其它辐射设备的非相参电磁信号的干扰。这些干扰的存在降低了SAR图像的信噪比，在干扰功率较大的情况下，会在图像中出现亮线和图像模糊；自激干扰是雷达本身信号泄漏所形成的干扰，该干扰形式的功率较大，会在SAR图像中形成干扰条带，覆盖目标；人为释放干扰这种强干扰信号能够遮盖或淹没有用信号，影响合成孔径雷达的成像性能，增加图像识别和后处理的难度。因此，抑制窄带干扰是SAR成像面临的关键问题。

张双喜等人在文献“基于小波变换的合成孔径雷达窄带干扰抑制算法”(西安电子科技大学学报，第38卷第5期，2011年10月)提出了一种基于小波变换的窄带干扰抑制方法，该方法首先将回波信号的频谱变换到小波域中，再把窄带干扰所对应的小波系数检测出来并通过设计具体的时频域滤波器把这些干扰系数去除掉，尽量减小信号的损失，实现窄带干扰的抑制。但是该方法的不足是，在不同的尺度下均需要对选定提取的窄带干扰的阈值进行设置和恒虚警概率检测，导致对小波变换及其每一尺度下的小波系数检测需要消耗大量的时间。

任超等人在文献“一种改进的空时自适应处理干扰抑制算法”(兵工学报，第31卷第12期，2010年12月)提出了采用改进的空时自适应方法进行干扰抑制的方法，该方法是将级联IIR陷波滤波器与空时FIR滤波器组相结合，级联IIR滤波器可利用有限的阶数在频域实现很窄的零陷，有效抑制窄带干扰。但是该方法的不足是，在调节陷波器参数时采用的是LMS自适应算法，它收敛速度较慢，且很难估计窄带干扰的频率。

发明内容

本发明的目的在于克服上述现有技术的不足，提出一种基于复数经验模态分解的时变窄带干扰抑制方法。本发明弥补了小波变换中多次阈值设置和恒虚警检测的耗时多，实时处理能力有限，改进的空时自适应方法收敛速度慢，需要精确估计窄带干扰的频率。而本发明所涉及的复数经验模态分解方法窄带干扰抑制方法充分利用窄带干扰的识别，并在抑制窄带干扰的同时充分保留了信号本身的非线性和非平稳特征。

实现本发明的基本思路是：将回波信号进行复数经验模态分解，分解为不同的本征模态函数，对分解得到的本征模态函数进行分析识别存在窄带干扰，自适应地抑制干扰对应的本征模态函数，接着设置阈值干扰的强干扰分量，再通过数据重构得到无干扰的数据，然后进行对无干扰的信号进行距离脉压和参数估计，根据估计的运动误差进行运动补偿，再对补偿后的信号进行方位脉压，最后对方位脉压信号进行多视处理，得到清晰的SAR图像。

为实现上述目的，本发明包括以下步骤：

(1)获取SAR回波

1a)雷达获取SAR回波数据；

1b)对回波数据在每一个距离单元上进行傅立叶变换，得到距离域信号；

(2)判断是否存在时变窄带干扰

2a)若存在，则对距离域信号作傅立叶逆变换，得到时域信号，执行步骤(3)；

2b)若不存在，执行步骤(4c)；

(3)复数经验模态CEMD分解

3a)将时域信号作为初始信号，确定时域信号发射方向个数；

3b)采用矢量相乘方法计算每一个发射方向上的投影矢量；

3c)绘制投影矢量波形图，在投影矢量波形图上找出多个极大值，通过三次样条插值得到投影矢量的包络；

3d)将原点与投影矢量极大值点的均值位置连接的矢量作为平均旋转矢量，获得一个本征模态函数；

3e)判断本征模态函数是否满足经验条件，若满足，则记录本征模态函数，执行步骤3f)；否则，执行步骤3b)；

3f)将初始信号减去本征模态函数作为初始信号，执行步骤3a)，直到获得所有本征模态函数；

3g)将所有本征模态函数和残余分量叠加得到分解后信号；

3h)将4倍的时域信号赋值为窄带干扰的阈值抑制强干扰分量，去除分解后信号中高于阈值抑制强干扰分量的信号部分，得到干扰抑制后信号；

(4)距离脉压及走动校正

4a)将干扰抑制后信号在方位向平均分块，方位子块个数与采样点数相同；

4b)对每个方位子块在距离上进行傅里叶变换，得到干扰抑制后距离域信号；

4c)对干扰抑制后距离域信号进行基于距离-方位维的距离线性走动校正，得到校正后频域信号；

4d)将校正后频域信号乘以距离脉压参考函数，并进行傅里叶逆变换，得到距离脉压后时域信号；

(5)参数估计

5a)采用相邻脉冲相关法进行多普勒中心估计；

5b)采用图像偏置法对每个方位子块进行瞬时多普勒调频率估计；

5c)采用基于瞬时多普勒调频率估计的运动误差分离与估计方法实时估计出天线相位中心的运动误差；

(6)运动补偿

6a)将步骤(4)得到的距离脉压后时域信号乘以运动误差对应的相位项，得到运动误差相位补偿后信号；

6b)将运动误差相位补偿后信号进行傅里叶变换，得到相位补偿后频域信号；

6c)相位补偿后频域信号乘以线性校正相位，并进行傅里叶逆变换，得到包络补偿后时域信号；

(7)方位脉压

7a)将方位窗函数乘以多普勒中心偏移复共轭相位项与多普勒调制复共轭相位项作为方位压缩参考函数；

7b)将包络补偿后时域信号与方位压缩参考函数作卷积，得到方位脉压后信号；

(8)对方位脉压后信号进行非相干叠加，完成多视处理，得到清晰的SAR图像。

本发明与现有技术相比具有以下优点：

第一，本发明采用了复数经验模态分解方法，将干扰信号分解为不同的本征模态函数，自适应地抑制干扰对应的本征模态函数，通过数据重构得到无干扰的数据，克服了现有技术小波变换中多次阈值设置和恒虚警检测的耗时多导致对目标的成像缺乏实时性，使得本发明具有了实现简单，效率高，实时性高的优点。

第二，本发明采用了对本征模态函数进行分析识别存在窄带干扰，克服了现有技术中改进的空时自适应方法需要对窄带干扰频率精确估计的不足，使得本发明具有了收敛速度快，无需对窄带干扰频率进行精确估计，并能充分保留信号本身的非线性和非平稳特征的优点。

附图说明

图1为本发明的流程图；

图2为本发明成像结果对比图。

具体实施方式

下面结合附图1对本发明的步骤做进一步详细描述。

步骤1：获取SAR回波

雷达以一定的脉冲重复频率发射并接收脉冲，得到以距离为行向量以方位为列向量的SAR回波，对回波数据在每一个距离单元上根据下式进行傅立叶变换，得到距离域信号：

$S\left(k\right)=\underset{K=0}{\overset{N-1}{\mathrm{\Σ}}}s\left(n\right)\exp \left(-j\frac{2\mathrm{\π}}{N}\mathrm{kn}\right)$

其中，S(k)为距离域信号，∑为求和符号，k为第k个距离单元，N为距离采样点数，s(n)为回波时域信号，exp表示指数函数的底，j为虚数单位，n为第n个距离采样点。

步骤2：观察距离域信号频谱，判断是否存在时变窄带干扰

若距离域信号频谱中出现窄带尖峰，尖峰幅度为回波信号幅度的几倍或几十倍，则判断存在窄带干扰，再对距离域信号根据下式进行傅立叶逆变换，得到时域信号，执行步骤3：

$s\left(n\right)=\frac{1}{N}\underset{N}{\overset{N-1}{\mathrm{\Σ}}}S\left(k\right)\exp \left(j\frac{2\mathrm{\π}}{N}\mathrm{kn}\right)$

其中，∑为求和符号，s(n)为回波时域信号，N为距离采样点数，∑为求和符号，n为第n个距离采样点，S(k)为距离域信号，k为第k个距离单元，exp表示指数函数的底，j为虚数单位。

若未出现窄带尖峰，则判断不存在窄带干扰，执行步骤4。

步骤3：复数经验模态CEMD分解

将时域信号作为初始信号，确定时域信号发射方向个数，方向个数为32～64，其中个数需要根据软硬件平台要求进行确定，发射方向满足

1≤n≤N

将初始信号和代表方向的相位项

相乘并取它的实部，采用矢量相乘方法，实现初始信号在发射方向上的投影，计算每一个发射方向上的投影矢量

其中

第n个发射方向，N为发射方向个数，j为虚数单位。

绘制投影矢量波形图，在投影矢量波形图上找出多个极大值，通过三次样条插值得到投影矢量的包络，将原点与投影矢量极大值点的均值位置连接的矢量作为平均旋转矢量，通过一个发射方向的时域信号减去平均旋转矢量，得到一个本征模态函数：

h_k(t)＝x_k(t)-m_k(t)

其中h_k(t)为本征模态函数，k为第k次循环，x_k(t)为初始信号，m_k(t)为平均旋转矢量。

根据下式判断本征模态函数是否满足经验条件：

$\mathrm{SD}=\mathrm{\&Sigma;}\frac{{|h_{i+1}\left(t\right)-h_i\left(t\right)|}^2}{{h^2}_{i+1}\left(t\right)}<0.3$

其中，SD为循环中止条件参数，∑为求和符号，h_i(t)为本征模态函数，下标i和i+1表示相邻两个发射方向的间距，若满足，则记录本征模态函数，将初始信号与本征模态函数相减作为初始信号，重复步骤3，若不满足，则进行下一个循环的判定，直到获得所有本征模态函数，确定最后一次分解后残余分量，初始信号与本征模态函数相减得到，将所有的本征模态函数和残余信号叠加得到分解后信号，将4倍的时域信号赋值为窄带干扰的阈值抑制强干扰分量，观察分解后信号的幅度，去除掉其中幅度高于阈值抑制强干扰分量的部分，得到干扰抑制后信号。

步骤4：距离脉压及走动校正

将干扰抑制后信号在方位向平均分块，方位子块个数与采样点数相同，这里取为512，其中个数需要根据仿真时的软硬件条件进行选择，对每个方位子块在距离上进行傅里叶变换，得到干扰抑制后距离域信号，再对干扰抑制后距离域信号乘以线性频率相位因子

其中exp为指数函数的底，j为虚数单位，ΔR为距离走动量，c为光速，f_r为发射信号载频，进行基于距离-方位维的距离线性走动校正，得到校正后信号。

将校正后信号乘以距离脉压参考函数

其中

为矩形窗函数，exp为指数函数的底，

为快时间，γ为发射信号的调频率，得到距离脉压后频域信号。

对距离脉压后频域信号进行傅里叶逆变换，得到距离脉压后时域信号。

步骤5：参数估计

采用相邻脉冲相关法，对相邻脉冲在对应距离点上的相位共轭相乘再求平均，得到多普勒中心的估计值。

对第一个脉冲依次在其距离点上的相位取共轭，与第二个脉冲对应的距离点上的相位共轭相乘，再将所有距离点对应相乘的结果进行作和再求其平均，第二个脉冲和第三个脉冲作同样的操作，以此类推，直到所有的脉冲做完该操作，从而获得多普勒中心的估计值。

采用图像偏置法对每个方位子块进行瞬时多普勒调频率估计。

第一步，将距离脉压后时域信号平均分为两段信号，前后半段信号分别乘以

其中exp为指数函数的底，j为虚数单位，γ_e0为采用惯导系统对多普勒调频率的粗估计值，t_m为慢时间，得到前后孔径信号；

第二步，对前后孔径信号分别进行方位傅里叶变换，得到前后孔径信号的两视图像；

第三步，根据下式，得到估计的前后孔径信号的多普勒差：

${\widehat{\mathrm{\&gamma;}}}_e^{\&prime;}=\frac{4\&times;{\mathrm{PRF}}^2}{N^2}\mathrm{\&Delta;n}$

其中

为多普勒差，PRF为脉冲重复频率，N为采样点数，Δn为两视图像相关估计移动量；

第四步，根据下式得到估计的多普勒调频率

${\widehat{\mathrm{\&gamma;}}}_e={\widehat{\mathrm{\&gamma;}}}_e^{\&prime;}+{\mathrm{\&gamma;}}_{e0}$

其中，

为估计的多普勒调频率，

为多普勒差，γ_e0为采用惯导系统对多普勒调频率的粗估计值。

将瞬时多普勒调频率分解为关于沿航向速度误差项和径向加速度误差项，分离出沿航向速度误差和径向加速度误差，进而获得天线相位中心的运动误差。

步骤6：运动补偿

首先，将步骤4得到的距离脉压后时域信号乘以运动误差对应的相位项

得到运动误差相位补偿后信号，其中exp为指数函数的底，

Δr为运动误差，t为时间，Δa为径向加速的误差，λ回波信号波长。

将运动误差相位补偿后信号进行傅里叶变换，得到相位补偿后频域信号；

相位补偿后频域信号乘以线性校正相位

并进行傅里叶逆变换，得到包络补偿后信号，其中exp为指数函数的底，j为虚数单位，f_r为发射信号载频，r为运动误差，c为光速，得到包络补偿后频域信号。

其次，对沿航向速度误差进行补偿，补偿如下：其中Φ(t)为补偿函数，t为时间，

为多普勒调频率的变化，V为载机的平均速度，Δv为沿航向速度误差，R为目标与雷达航线的最短距离。

步骤7：方位脉压

将方位窗函数a_a(t_m)，乘以多普勒频移相位项

其中exp为指数函数的底，j为虚数单位，f_d表示多普勒频率，t_m表示慢时间，γ表示发射信号的调频率，作为方位压缩参考函数。

将包络补偿后信号与方位压缩参考函数作卷积，得到方位脉压后信号。

步骤8：对方位脉压后信号进行非相干叠加，完成多视处理，得到清晰的SAR图像。

下面结合附图2对本发明的效果做进一步说明。

附图2为根据机载条带式SAR对某地区的田野和村庄录取的实测数据采用不同方法的进行成像的结果对比图。图2(a)为本发明实测数据采用自适应频域置零陷波法的成像结果图，图2(b)为本发明实测数据采用CEMD分解的成像结果图。

实测数据的参数如下：工作波段为X波段，发射脉冲重复周期T_p＝10us，带宽B＝180MHz，成像分辨率为1m×1m，数据受到时变的NBI干扰。回波数据矩阵为4084×1792(距离单元×方位单元)。

图2(a)中横坐标为方位单元，纵坐标为距离单元，可以看到图中目标出现了散焦，这是因为频域陷波法在抑制干扰的同时，造成了频谱裂痕，导致旁瓣增高，另外还有剩余干扰信号的痕迹，这是由于实测数据的信噪比一般不是很高，在成像结果中被干扰覆盖的区域并没有很好地重现出来。

图2(b)中横坐标为方位单元，纵坐标为距离单元。由图2(b)可以看到采用CEMD分解的成像图中有效地抑制窄带干扰，具有良好的聚焦效果，较好地恢复了被干扰遮盖的目标信号，图像的信噪比更高，对比度更好，很好的消除掉干扰痕迹，成像更清晰。

Claims (9)

1.一种基于复数经验模态分解的时变窄带干扰抑制方法，包括以下步骤：

(1)获取SAR回波

1a)雷达获取SAR回波数据；

1b)对回波数据在每一个距离单元上进行傅立叶变换，得到距离域信号；

(2)判断是否存在时变窄带干扰

2a)若存在，则对距离域信号作傅立叶逆变换，得到时域信号，执行步骤(3)；

2b)若不存在，执行步骤(4c)；

(3)复数经验模态CEMD分解

3a)将时域信号作为初始信号，确定时域信号发射方向个数；

3b)采用矢量相乘方法计算每一个发射方向上的投影矢量；

3c)绘制投影矢量波形图，在投影矢量波形图上找出多个极大值，通过三次样条插值得到投影矢量的包络；

3d)将原点与投影矢量极大值点的均值位置连接的矢量作为平均旋转矢量，获得一个本征模态函数；

3e)判断本征模态函数是否满足经验条件，若满足，则记录本征模态函数，执行步骤3f)；否则，执行步骤3b)；

3f)将初始信号减去本征模态函数作为初始信号，执行步骤3a)，直到获得所有本征模态函数；

3g)将所有本征模态函数和残余分量叠加得到分解后信号；

3h)将4倍的时域信号赋值为窄带干扰的阈值抑制强干扰分量，去除分解后信号中高于阈值抑制强干扰分量的信号部分，得到干扰抑制后信号；

(4)距离脉压及走动校正

4a)将干扰抑制后信号在方位向平均分块，方位子块个数与采样点数相同；

4b)对每个方位子块在距离上进行傅里叶变换，得到干扰抑制后距离域信号；

4c)对干扰抑制后距离域信号进行基于距离-方位维的距离线性走动校正，得到校正后频域信号；

4d)将校正后频域信号乘以距离脉压参考函数，并进行傅里叶逆变换，得到距离脉压后时域信号；

(5)参数估计

5a)采用相邻脉冲相关法进行多普勒中心估计；

5b)采用图像偏置法对每个方位子块进行瞬时多普勒调频率估计；

5c)采用基于瞬时多普勒调频率估计的运动误差分离与估计方法实时估计出天线相位中心的运动误差；

(6)运动补偿

6a)将步骤(4)得到的距离脉压后时域信号乘以运动误差对应的相位项，得到运动误差相位补偿后信号；

6b)将运动误差相位补偿后信号进行傅里叶变换，得到相位补偿后频域信号；

6c)相位补偿后频域信号乘以线性校正相位，并进行傅里叶逆变换，得到包络补偿后时域信号；

(7)方位脉压

7a)将方位窗函数乘以多普勒中心偏移复共轭相位项与多普勒调制复共轭相位项作为方位压缩参考函数；

7b)将包络补偿后时域信号与方位压缩参考函数作卷积，得到方位脉压后信号；

(8)对方位脉压后信号进行非相干叠加，完成多视处理，得到清晰的SAR图像。

2.根据权利要求1所述的基于复数经验模态分解的时变窄带干扰抑制方法，其特征在于：步骤3a)所述的发射方向个数为32～64。

3.根据权利要求1所述的基于复数经验模态分解的时变窄带干扰抑制方法，其特征在于：步骤3b)所述的矢量相乘方法就是将初始信号和代表方向的相位项相乘并取实部，实现初始信号在发射方向上的投影。

4.根据权利要求1所述的基于复数经验模态分解的时变窄带干扰抑制方法，其特征在于：步骤3d)所述的本征模态函数是由一个发射方向的时域信号减去该发射方向上的平均旋转矢量得到。

5.根据权利要求1所述的基于复数经验模态分解的时变窄带干扰抑制方法，其特征在于：步骤3e)所述的经验条件的计算公式为

$\mathrm{SD}=\mathrm{\&Sigma;}\frac{{|h_{i+1}\left(t\right)-h_i\left(t\right)|}^2}{{h^2}_{i+1}\left(t\right)}<0.3$

其中，SD为循环中止条件参数，∑为求和符号，h_i(t)为本征模态函数，下标i和i+1表示两个相邻中间间距。

6.根据权利要求1所述的基于复数经验模态分解的时变窄带干扰抑制方法，其特征在于：步骤3g)所述的残余信号为最后一次分解的剩余信号，即由初始信号与最后一个本征模态函数相减得到。

7.根据权利要求1所述的基于复数经验模态分解的时变窄带干扰抑制方法，其特征在于：步骤4c)所述的距离线性走动校正指的是频域信号乘以线性频率相位因子得到的。

8.根据权利要求1所述的基于复数经验模态分解的时变窄带干扰抑制方法，其特征在于：步骤5a)所述的相邻脉冲相关法是对相邻脉冲在对应距离点上的相位共轭相乘再求平均，得到多普勒中心的估计。

9.根据权利要求1所述的基于复数经验模态分解的时变窄带干扰抑制方法，其特征在于：步骤5c)所述的基于瞬时多普勒调频率估计的运动误差分离与估计方法是从瞬时多普勒调频率中分离出沿航向速度误差和径向加速度误差，获得天线相位中心的运动误差。

Images