基于时频谱图分解的SAR时变窄带干扰抑制方法

摘要

基于时频谱图分解的SAR时变窄带干扰抑制方法，包括以下步骤：获取回波数据；判断回波数据是否存在窄带干扰，若存在干扰执行下一步；根据存在窄带干扰的回波数据的列向量获得混合信号时频谱图，并进行奇异值分解，获得对应的待分离信号矩阵；对待分离信号矩阵进行分离，获得干扰信号矩阵；根据最小二乘准则利用干扰信号矩阵中的干扰信号分量重构干扰信号时频谱图：从混合信号时频谱图中剔除干扰信号得到有用信号时频谱图，并对有用信号时频谱图进行逆短时傅里叶变换，得到剔除干扰后的回波数据的列向量；判断是否遍历完所有回波数据的列向量，若遍历完则得到剔除干扰后的回波数据。本发明具有运算量较小，效率较高并且没有数据缺失的优点。

说明书

技术领域

本发明属于信号处理技术领域，尤其涉及一种雷达信号处理技术中的基于 时频谱图分解的SAR时变窄带干扰抑制方法。

背景技术

窄带干扰是针对合成孔径雷达系统（Synthetic Aperture Radar,简称 SAR）的主要干扰方式之一。所谓窄带干扰是指干扰信号带宽与合成孔径雷达 发射信号带宽之比很小（如1%）的信号。窄带干扰可分为非人为干扰和人为干 扰两类。非人为干扰主要包括与SAR处于同一频段的通讯设备、电视网和其它 低频波段辐射设备；常见的人为干扰主要包括强脉冲干扰和频率引导“瞄准 式”干扰，人为窄带干扰具有干扰功率大、针对性强、干扰实现简单、形式多 变等特点。这些干扰的存在降低了SAR回波数据的信噪比，尤其是在干扰功率 较大的情况下，会使得图像中出现亮线并且使图像变模糊，从而严重降低了 SAR图像质量。因此，如何对窄带干扰形式进行有效的识别和抑制是目前SAR 干扰处理技术的重要课题。

Feng Zhou等人在文献“Narrow-Band Interference Suppression for  SAR Based on Independent Component Analysis”（IEEE Trans.Geosci. Remote Sens.，2013）中提出采用独立分量分析方法进行干扰抑制。该方法利 用有用回波信号与干扰信号的统计特性差异，采用独立分量分析方法直接在信 号数据域提取窄带干扰，并根据最小二乘原则获得对干扰基信号的复包络的精 确估计后，重构干扰信号，从原始回波中消去干扰信号，实现时变窄带干扰抑 制。但该方法存在以下不足：其利用时间平滑增加观测通道，减小了样本的数 目，致使估计出的干扰信号长度缩短，降低了独立性准则的计算精度，使得干 扰信号分离不彻底，从而影响了分离抑制的效果。

申请号为201210001678.2、发明名称为基于复数经验模态分解的时变窄带 干扰抑制方法的中国发明专利申请公开了一种利用复数经验模态分解的时变窄 带干扰抑制方法。该方法将信号分解成一系列本征模态函数，在抑制窄带干扰 的同时充分保留了信号本身的非线性和非平稳特征，但是该方法对局部信号分 析处理未考虑干扰信号的全局统计特性，甚至会分解出虚假的信号，从而不能 有效地刻画干扰信号的时频特性，导致干扰信号分离不彻底。

发明内容

针对现有技术存在的不足，本发明的目的是提供一种基于时频谱图分解的 SAR时变窄带干扰抑制方法，该方法可以准确提取窄带干扰信号的时频特征， 在保留有用信息的同时对时变窄带干扰进行抑制，从而实现对雷达回波信号清 晰成像。

为了实现上述目的，本发明采取如下的技术解决方案：

基于时频谱图分解的SAR时变窄带干扰抑制方法，包括以下步骤：

步骤1、获取回波数据；

雷达以脉冲重复频率发射并接收脉冲，得到以距离为行向量、以方位为列 向量的回波数据；

步骤2、判断回波数据是否存在窄带干扰，若存在干扰执行步骤3，否则 执行步骤7；

步骤3、根据存在窄带干扰的回波数据的列向量获得混合信号时频谱图 S，并对混合信号时频谱图S进行奇异值分解，获得对应的待分离信号矩阵Y；

对回波数据的列向量进行短时傅里叶变换，得到混合信号时频谱图S，并 按照下式对混合信号时频谱图S进行奇异值分解：

S＝UΛV^H

上式中的U为混合信号时频谱图的左特征向量矩阵，Λ为奇异值对角矩 阵，V为混合信号时频谱图的右特征向量矩阵，(·)^H表示矩阵共轭转置操作；

将分解得到的混合信号时频谱图的左特征向量矩阵U的前L个特征向量作 为行向量，构成待分离信号矩阵Y，其中的L为主分量个数，主分量个数L根 据最小描述长度准则计算得到：

$L=\arg \min \{-m(n-L)\log \left((n-L){\left(\underset{p=L+1}{\overset{n}{\Pi }}{\sigma }_p\right)}^{\frac{1}{n-L}}{\left(\underset{p=L+1}{\overset{n}{\Sigma }}{\sigma }_p\right)}^{-1}\right)+0.5L(2n-L)\log \left(m\right)\}$

上式中的argmin(·)表示计算使目标函数取最小值时对应变量值的运算操 作，m为混合信号时频谱图的左特征向量矩阵U的行数，n为混合信号时频谱 图的左特征向量矩阵U的列数，log(·)表示求对数运算操作，(·)表示对第 L+1至第n项的求乘积运算操作，σ_p表示奇异值对角矩阵Λ的第p个奇异值， (·)表示对第L+1至第n项的求和运算操作；

步骤4，对待分离信号矩阵Y进行分离，获得干扰信号矩阵I，包括以下子 步骤；

子步骤4a、随机产生大小为L×L的解混矩阵W；

子步骤4b、计算迭代矩阵z； 迭代矩阵z＝[z₁，z₂，…，z_p，…，z_L]，其中，z_p为迭代矩阵z的第p个列向 量，$z_p=E\left\{Y{\left(w_p^{H}Y\right)}^{*}{(0.1+{\left|w_p^{H}Y\right|}^2)}^{-1}\right\}-E\left\{0.1{(0.1+{\left|w_p^{H}Y\right|}^2)}^{-2}\right\}{w}_p,$E(·)表示求期望运算操作，Y为待分离信号矩阵，w_p为解混矩阵W的第p 个列向量，(·)^H表示矩阵共轭转置操作，(·)^*表示共轭运算操作，|·|表示求 绝对值运算操作；

子步骤4c、对解混矩阵和迭代矩阵的所有列向量依次判断是否满足 若所有列向量都满足该条件，则执行子步骤4d，否则，对不 满足条件的w_p值进行更新，令所述执行子步骤4b；其中，||·||表 示求模值运算操作，z_p为迭代矩阵z的第p个列向量，w_p为解混矩阵w的 第p个列向量，ε表示设定的阈值；

子步骤4d、将解混矩阵w各列向量去相关得到最终解混矩阵 $\tilde{W}=W{\left(W^{H}W\right)}^{-1/2};$

子步骤4e、根据最终解混矩阵计算干扰信号矩阵

步骤5、根据最小二乘准则，将干扰信号矩阵I中的干扰信号分量按照下 式重构干扰信号时频谱图：

$S_k=i_k\left(i_k^{+}S\right)$

其中，S_k表示第k个干扰信号分量的时频谱图，i_k表示干扰信号矩阵I的 第k个行向量，+表示求矩阵的广义逆运算操作，S为混合信号时频谱图；

步骤6、从混合信号时频谱图S中剔除干扰信号得到有用信号时频谱图S′， 并对有用信号时频谱图S′进行逆短时傅里叶变换，得到剔除干扰后的回波数据 的列向量；

有用信号时频谱图其中，|·|表示求绝对值运算操 作，表示对第1项至第L项的求和运算操作，S_k为第k个干扰信号的时频 谱图，exp(·)表示求指数运算操作，j为虚数符号，Θ为对混合信号时频谱图S 求得的相位角；

步骤7，判断是否遍历完所有回波数据的列向量，若遍历完，则执行步骤 8，否则执行步骤2；

步骤8，获取剔除干扰后的回波数据；

前述步骤中1≤p≤L，1≤k≤L。

进一步的，所述步骤2中通过一维距离像的峭度来判断回波数据是否存在 干扰，判断过程为：对步骤1的回波数据的列向量进行傅里叶变换，得到一维 距离像，判断一维距离像的峭度κ是否大于等于设定的幅度阈值，若是，则认 为回波数据中存在窄带干扰；

一维距离像的峭度

上式中的E()表示求期望运算操作，y为一维距离像的幅度值，μ为一维距 离像幅值的均值，d为一维距离像幅值的标准差。

进一步的，所述幅度阈值为10～200。

本发明方法通过对回波数据进行短时傅里叶变换增加观测通道，并对时频 谱图进行奇异值分解，降低待处理数据的维数，克服了现有技术中独立分量分 析方法通过时间平滑增加观测通道导致样本缺失的不足，使得本发明运算量 小，效率高并且没有数据缺失；其次，本发明方法充分利用回波信号和干扰信 号的时频特性差异，通过盲源分离方法将干扰信号与回波信号在时频谱图域独 立的分离出来，能够很好的刻画干扰信号的时频特性，克服了现有技术中基于 复数经验模态分解方法仅考虑局部统计特性而造成干扰抑制不充分的缺点，使 得本发明具有信号损失小，抑制效果好的优点。

附图说明

图1为本发明方法的流程图；

图2为仿真实验对混合信号进行傅里叶变换的距离维幅频特性图；

图3为仿真实验对混合信号进行短时傅里叶变换的二维时频谱图；

图4为仿真实验对混合信号进行短时傅里叶变换的三维时频谱图；

图5为仿真实验对提取的干扰信号进行短时傅里叶变换的二维时频谱图；

图6为仿真实验对提取的干扰信号进行短时傅里叶变换的三维时频谱图；

图7为仿真实验对剔除干扰信号后的有用信号进行短时傅里叶变换的二维 时频谱图；

图8为仿真实验对剔除干扰信号后的有用信号进行短时傅里叶变换的三维 时频谱图。

具体实施方式

为了让本发明的上述和其它目的、特征及优点能更明显，下文特举本发明 实施例，并配合所附图示，做详细说明如下。

本发明方法的基本思路是：获得的回波数据后，首先，对回波数据进行窄 带干扰定性识别；其次，对存在干扰的回波数据进行短时傅里叶变换得到相应 的二、三维的混合信号时频谱图，并对其进行奇异值分解，获得对应的待分离 信号矩阵；接着，利用窄带干扰信号与有用信号在时频域的差异性和独立性， 采用盲源分离方法对得到的存在干扰的待分离信号矩阵进行分离，得到干扰信 号矩阵；之后，采用最小二乘方法重构对应的二、三维的干扰信号分量的时频 谱图，并与混合信号时频谱图进行相减，得到有用回波信号时频谱图；最后， 对有用回波信号时频谱图进行短时傅里叶变换得到消去干扰信号后的有用回波 信号，实现SAR时变窄带干扰抑制。

以上是本发明的核心思想，下面将结合本发明实施例中的附图，对本发明 实施例的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发 明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通 技术人员在没有做出创造性劳动前提下获得的所有其他实施例，都属于本发明 保护的范围。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明 还可以采用其它不同于在此描述的其它方式来实施，本领域技术人员可以在不 违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例 的限制。

参照图2，图2为本发明方法的流程图，本发明方法的步骤如下：

步骤1、获取回波数据；

雷达以脉冲重复频率发射并接收脉冲，得到以距离为行向量、以方位为列 向量的回波数据；

步骤2、判断回波数据是否存在窄带干扰；

本发明可优选通过一维距离像的峭度来判断回波数据是否存在干扰，判断 过程为：对回波数据的列向量进行傅里叶变换，得到一维距离像，判断一维距 离像的峭度κ是否大于等于设定的幅度阈值，若是，则认为回波数据的列向量 中存在窄带干扰，执行步骤3，否则执行步骤7；优选的，前述幅度阈值为 10～200；

一维距离像的峭度

上式中的E(·)表示求期望运算操作，y为一维距离像的幅度值，μ为一维距 离像幅值的均值，d为一维距离像幅值的标准差；

步骤3、根据存在窄带干扰的回波数据的列向量获得混合信号时频谱图 S，并对混合信号时频谱图S进行奇异值分解，获得对应的待分离信号矩阵Y；

对回波数据的列向量进行短时傅里叶变换，得到混合信号时频谱图S，并 按照下式对混合信号时频谱图S进行奇异值分解：

S＝UΛV^H

上式中的U为混合信号时频谱图的左特征向量矩阵，Λ为奇异值对角矩 阵，V为混合信号时频谱图的右特征向量矩阵，(·)^H表示矩阵共轭转置操作；

将分解得到的混合信号时频谱图的左特征向量矩阵U的前L个特征向量作 为行向量，构成待分离信号矩阵Y，其中的L为主分量个数，主分量个数L根 据最小描述长度准则计算得到：

$L=\arg \min \{-m(n-L)\log \left((n-L){\left(\underset{p=L+1}{\overset{n}{\Pi }}{\sigma }_p\right)}^{\frac{1}{n-L}}{\left(\underset{p=L+1}{\overset{n}{\Sigma }}{\sigma }_p\right)}^{-1}\right)+0.5L(2n-L)\log \left(m\right)\}$

上式中的argmin(·)表示计算使目标函数取最小值时对应变量值的运算操 作，m为混合信号时频谱图的左特征向量矩阵U的行数，n为混合信号时频谱 图的左特征向量矩阵U的列数，log(·)表示求对数运算操作，表示对第 L+1至第n项的求乘积运算操作，σ_p表示奇异值对角矩阵Λ的第p个奇异值， 表示对第L+1至第n项的求和运算操作；

步骤4，对待分离信号矩阵Y进行分离，获得干扰信号矩阵I；

窄带干扰信号与有用信号在时频域具有差异性和独立性，利用该差异性和 独立性可采用盲源分离法进行分离，得到干扰信号矩阵I，其包括以下子步 骤；

子步骤4a、随机产生大小为L×L的解混矩阵W；

子步骤4b、计算迭代矩阵z；

迭代矩阵z＝[z₁，z₂，…，z_p，…，z_L]，其中，z_p为迭代矩阵z的第p个列向 量，$z_p=E\left\{Y{\left(w_p^{H}Y\right)}^{*}{(0.1+{\left|w_p^{H}Y\right|}^2)}^{-1}\right\}-E\left\{0.1{(0.1+{\left|w_p^{H}Y\right|}^2)}^{-2}\right\}{w}_p,$E(·)表示求期望运算操作，Y为待分离信号矩阵，w_p为解混矩阵W的第p 个列向量，(·)^Ｈ表示矩阵共轭转置操作，(·)^*表示共轭运算操作，|·|表示求 绝对值运算操作；

子步骤4c、对解混矩阵和迭代矩阵的所有列向量依次判断是否满足 若所有列向量都满足该条件，则执行子步骤4d，否则，对不 满足条件的w_p值进行更新，令所述执行子步骤4b；其中，||·||表 示求模值运算操作，z_p为迭代矩阵z的第p个列向量，w_p为解混矩阵w的 第p个列向量，ε表示设定的阈值，一般取1e-3；

子步骤4d、将解混矩阵w各列向量去相关得到最终解混矩阵 $\tilde{W}=W{\left(W^{H}W\right)}^{-1/2};$

子步骤4e、根据最终解混矩阵计算干扰信号矩阵

步骤5、根据最小二乘准则，将干扰信号矩阵I中的干扰信号分量按照下 式重构干扰信号时频谱图：

$S_k=i_k\left(i_k^{+}S\right)$

其中，S_k表示第k个干扰信号分量的时频谱图，i_k表示干扰信号矩阵I的 第k个行向量，+表示求矩阵的广义逆运算操作，S为混合信号时频谱图；

步骤6、从混合信号时频谱图S中剔除干扰信号得到有用信号时频谱图S′， 并对有用信号时频谱图S′进行逆短时傅里叶变换，得到剔除干扰后的回波数据 的列向量；

有用信号时频谱图其中，|·|表示求绝对值运算操 作，表示对第1项至第L项的求和运算操作，S_k为第k个干扰信号分量的 时频谱图，exp(·)表示求指数运算操作，j为虚数符号，Θ为对混合信号时频谱 图S求得的相位角；

步骤7，判断是否遍历完所有回波数据的列向量，若遍历完则执行步骤 8，否则执行步骤2；

步骤8，获得剔除干扰后的回波数据；

前述步骤中1≤p≤L，1≤k≤L。

本发明方法可以弥补独立分量分析法导致样本缺失、复数经验模态分解法 对干扰信号分离不彻底的不足，充分利用雷达信号与干扰信号时频特性的差异 性，利用时频谱图分解方法直接在信号时频域提取窄带干扰并进行抑制，能够 保留有用回波信号，有效地抑制时变的窄带干扰。

本发明的效果可以通过以下的仿真实验进一步说明，仿真时采用MATLAB （R2010b）软件进行仿真。

仿真条件如下：

由于雷达回波信号可以近似为线性调频信号，因此，仿真观测信号由一个 线性调频信号,一个窄带干扰信号和噪声混合组成，线性调频信号的调频率为 7×10¹²Hz/s，干扰信号载频f_c为8MHz，采样频率为150MHz，观测时间为[- 3.4,3.4]μs，信噪比为10dB。

仿真内容

如图2所示，对混合信号进行傅里叶变换后得到图2所示的距离维幅频特 性图，其中横坐标为频率单元，纵坐标为归一化幅度值。从图2可看出，由于 窄带干扰的存在，线性调频信号频谱中存在明显的带内冲激干扰，在一个小频 率段中幅值出现了突变。

如图3和图4所示，对混合信号进行短时傅里叶变换后，得到图3所示的 混合信号二维时频谱图和图4所示的混合信号三维时频谱图。图3的横坐标为 时间单元，纵坐标为频率单元，由图3可看出，斜线为线性调频信号的时频分 布，水平线为干扰分量的时频分布，干扰分量的存在影响了对线性调频信号的 时频特性的辨识；图4的X坐标为频率单元，Y坐标为时间单元，Z坐标为幅 度值，由图4可看出，突起为干扰分量的时频分布，底部隆起的尖峰为线性调 频信号的时频分布，干扰分量能量很强，压制了线性调频信号的时频特征，从 而影响了对线性调频信号的时频特性的辨识。

如图5和图6所示，采用时频谱图分解方法提取窄带干扰分量得到重构的 如图5所示的干扰信号二维时频谱图和如图6所示的干扰信号三维时频谱图。 图5的横坐标为时间单元，纵坐标为频率单元，由图5可看出，水平线为干扰 分量的时频分布，与图3相比，提取的干扰信号，只保留了干扰信号的时频信 息,没有混杂线性调频信号的时频信息；图6的X坐标为频率单元，Y坐标为时 间单元，Z坐标为幅度值，由图6可看出，突起为高能量干扰信号的时频分 布，与图4相比，提取的干扰信号，只保留了干扰信号的时频信息,并很好的 刻画了时频谱图的幅度强度和时变特性。

如图7和图8所示，将重构的干扰信号从混合信号中剔除后得到图7所示 的有用信号二维时频谱图和图8所示的有用信号三维时频谱图。图7的横坐标 为时间单元，纵坐标为频率单元，由图7可看出，斜线为线性调频信号的时频 分布，被恢复出来的信号在时频域中已没有与窄带干扰的交叉项，只保留了原 线性调频信号的时频信息，由于高能量的干扰信号被去除，噪声的时频特征也 可清晰可见；图8的X坐标为频率单元，Y坐标为时间单元，Z坐标为幅度 值，由图8可看出，突出的斜线表示恢复的线性调频信号的时频变化信息，底 部的随机起伏表示噪声的时频特征，高能量的干扰信号被清除，恢复的线性调 频信号时频信息被很好的保留。通过仿真实验说明了本方法的有效性。

以上所述，仅是本发明的较佳实施例而已，并非对本发明做任何形式上的 限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何 熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示 的技术内容做出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发 明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修 改、等同变化与修饰，均仍属于本发明技术方案的范围内。