一种基于全频带谱分析的雷达压制式干扰检测方法与流程

本发明涉及雷达信号处理技术领域，尤其涉及一种基于全频带谱分析的雷达压制式干扰检测方法。

背景技术：

现代战争中电磁环境日益复杂，雷达抗干扰面临越来越多的挑战，压制式干扰是对雷达威胁较大的一类干扰；压制式干扰是通过释放强噪声干扰，其从天线旁瓣进入雷达主通道，与目标回波叠加在一起，使雷达接收机过载、饱和，严重影响雷达的正常检测性能。

雷达抗干扰能力已成为雷达系统设计的一项核心战术指标，干扰检测性能也已成为雷达抗干扰能力的重要组成部分。良好的干扰检测性能使得雷达能够根据干扰检测结果及时自动选取最优的抗干扰措施，进而使得雷达的整体抗干扰性能得到充分发挥。

雷达干扰检测功能实现的方法很多，但多数方法较复杂，不利于工程实现，具有很大的应用局限性；对于传统机械扫描雷达来说，受到自身体制局限性，通常采用在主通道外增加一路辅助通道进行干扰检测，用此方法检测得到的信息往往存在虚警，不够准确；此外，由于主通道和辅助通道工作在当前雷达工作频点下，所以只有干扰处在雷达当前工作频点时才能够被有效检测到，大大影响了干扰检测功能的完备性。

技术实现要素：

针对以上现有技术存在的不足，本发明的目的在于提出一种基于全频带谱分析的雷达压制式干扰检测方法，该方法能够实现雷达所有工作频带的频谱功率分析，方法简单，实用性强，识别率高，易于工程实现，并提高了雷达的抗干扰能力。

为达到上述目的，本发明的实施例采用如下技术方案予以实现。

一种基于全频带谱分析的雷达压制式干扰检测方法，所述方法包括如下步骤：

步骤1，获取雷达发射第l个脉冲后的回波射频信号，对所述回波射频信号依次进行模数转换和正交采样，得到第l个脉冲覆盖雷达工作频带的基带数据；从所述基带数据中均匀截取M段子基带数据，对每段子基带数据依次进行时域加窗和N点快速傅里叶变换，得到该段子基带数据的N个频点数据，从而得到第l个脉冲M段子基带数据的M×N个频点数据；×表示乘号，其中，l初值为1，l＝1，...，L，L为雷达一个扫描周期内包含的脉冲个数，且为自然数，M≥1，N≥1；

步骤2，分别计算第l个脉冲M段子基带数据的M×N个频点数据的模值平方，将M段子基带数据中每段子基带数据对应的第n个频点的模值平方进行累加，得到第n个频点对应的频谱值，令n依次取1，…，N，从而得到第l个脉冲N个频点分别对应的频谱值；将所述N个频点分别对应的频谱值进行累加，得到第l个脉冲的频谱累加值；M和N为自然数

步骤3，令l加1，并依次重复步骤1至步骤2，从而得到雷达一个扫描周期内L个脉冲分别对应的频谱累加值；

步骤4，获取加权中值滤波线性窗，并根据雷达一个扫描周期内L个脉冲分别对应的频谱累加值计算每个脉冲的频谱累加值的中间值、频谱累加值的平均值，然后根据每个脉冲的频谱累加值的中间值、频谱累加值的平均值得到该脉冲对应的加权中值滤波后的值，从而得到L个脉冲对应的加权中值滤波后的值，并组成雷达一个扫描周期内L个脉冲的加权中值滤波结果；

步骤5，根据所述雷达一个扫描周期内L个脉冲的加权中值滤波结果，确定雷达一个扫描周期内的干扰区域集合，所述干扰区域集合包含J个干扰区域，J为自然数；

步骤6，获取J个干扰区域内包含的脉冲分别对应的方位码，根据J个干扰区域内包含的脉冲分别对应的方位码、J个干扰区域内包含的脉冲分别对应的加权中值滤波后的值，确定J个干扰区域内包含的脉冲分别对应的干扰方位、干扰强度、干扰类型以及干扰频点。

本发明的有益效果：第一，本发明方法结构清晰，易于硬件实现，实用性高；第二，本发明方法无需通过切换雷达工作频点，就能够实现雷达工作频带的全频带干扰检测，应用范围大大增强；第三，中值滤波方法的使用，有利于减少脉冲噪声的带来的检测虚警，自适应门限的应用，提高应用的灵活性，提高系统适应环境变化的能力；第四，本发明参数可变，能够根据使用过程中的外界环境变化和应用需求进行灵活调节，有效提高方法使用的灵活性、复用性和可移植性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明一种基于全频带谱分析的雷达压制式干扰检测方法流程示意图；

图2a为切比雪夫窗时域波形示意图；

图2b为切比雪夫窗频域波形示意图；

图3为加切比雪夫窗与不加切比雪夫窗分别进行快速傅里叶变换后得到的频谱示意图；

图4为加权中值滤波算法的流程示意图；

图5为利用自适应门限划分干扰区域流程示意图；

图6a为中值滤波前频谱累加和的曲线示意图；

图6b为中值滤波后频谱累加和的曲线示意图；

图7a为方位融合后的干扰标志位示意图；

图7b为更新后的自适应门限示意图；

图8为某脉冲32频点功率分布示意图。

具体实施方式

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

本发明实施例提供一种基于全频带谱分析的雷达压制式干扰检测方法，参照图1，所述方法包括如下步骤：

步骤1，获取雷达发射第l个脉冲后的回波射频信号，对所述回波射频信号依次进行模数转换和正交采样，得到第l个脉冲覆盖雷达工作频带的基带数据；从所述基带数据中均匀截取M段子基带数据，对每段子基带数据依次进行时域加窗和N点快速傅里叶变换，得到该段子基带数据的N个频点数据，从而得到第l个脉冲M段子基带数据的M×N个频点数据；×表示乘号，其中，l初值为1，l＝1，...，L，L为雷达一个扫描周期内包含的脉冲个数，且为自然数，M≥1，N≥1。

对所述第l个脉冲的基带数据中截取的M段数据依次进行时域加窗和快速傅里叶变换，依次得到加窗后第l个脉冲的M段数据y_l1，…，y_lM和第l个脉冲的M段N点频域结果Y_ll，…，Y_lM；其中，M≥1，N≥1，每段数据包含N个频点，则第l个脉冲的基带数据中截取的M段数据包含M×N个频点。

步骤1的子步骤为：

(1a)i∈{1，…M}，对第l个脉冲的第1段N个频点数据x_l1利用切比雪夫窗进行时域加窗，x_l1＝[x_l11，…，x_l1N]，得到加窗后第l个脉冲的第1段数据y_l1，y_l1＝[y_l11，…，y_l1N]，然后对所述加窗后第l个脉冲的第1段数据y_l1做N点快速傅里叶变换，得到第l个脉冲的第1段N点频域数据Y_l1，Y_l1＝[Y_l11，…，Y_l1N]。

(1b)对第l个脉冲的第i段N个频点数据x_li利用切比雪夫窗进行时域加窗，x_li＝[x_li1，…，x_liN]，得到加窗后第l个脉冲的第i段数据y_li，y_li＝[y_li1，…，y_liN]，然后对所述加窗后第l个脉冲的第i段数据y_li做N点快速傅里叶变换，得到第l个脉冲的第i段N点频域数据Y_li，Y_li＝[Y_li1，…，Y_liN]。

(1c)令i加1，重复执行子步骤(1b)，直到依次得到加窗后第l个脉冲的第M段数据y_lM和第l个脉冲的第M段N点频域数据Y_lM，y_lM＝[y_lM1，…，y_lMN]，Y_lM＝[Y_lMl，…，Y_lMN]，此时得到加窗后第l个脉冲的M段数据y_l1，…，y_lM和第l个脉冲的M段N点频域数据Y_l1，…，Y_lM。

具体地，由于近区常存在强杂波而且接收机滤波器带外抑制有限，对真实干扰的分析产生不利影响，所以选取样本时应避开近区强杂波距离段，在第l个脉冲的基带数据中截取M段数据；因为干扰会出现在任意距离段，所以截取的M段数据应尽量均匀分布在第l个脉冲的基带数据中；对第l个脉冲的基带数据中截取的M段数据先要进行时域加窗以抑制频谱泄露，同时降低旁瓣干扰；参照图2a和图2b，图2a为切比雪夫窗时域波形图，图2b为切比雪夫窗频域波形图；切比雪夫窗被称作最大振幅比时窗函数，旁瓣具有等波动性，所有旁瓣具有相同的高度。

快速傅里叶变换由于物理意义明确，计算简单，是常用的分析工具；参照图3，为加切比雪夫窗与不加切比雪夫窗分别进行快速傅里叶变换后得到的频谱示意图；从图3可以看出，加切比雪夫窗后再进行快速傅里叶变换，能够起到非相干积累的效果，提升了干噪比，有利于后续分析；其中，N值的大小为雷达所有工作频带每个数据段中要分析的频点个数，M值的大小取决于雷达的最小检测干燥比，所述干燥比为干扰与噪声的功率比值。

步骤2，分别计算第l个脉冲M段子基带数据的M×N个频点数据的模值平方，将M段子基带数据中每段子基带数据对应的第n个频点的模值平方进行累加，得到第n个频点对应的频谱值，令n依次取1，…，N，从而得到第l个脉冲N个频点分别对应的频谱值；将所述N个频点分别对应的频谱值进行累加，得到第l个脉冲的频谱累加值；M和N为自然数。

步骤2中具体为：

(2a)令n＝1，…，N，n的初值为1；

(2b)将第l个脉冲M段子基带数据中每段子基带数据对应的第n个频点的模值平方进行累加，得到第n个频点对应的频谱值Z_ln；

(2c)令n的值加1，并重复执行子步骤(2b)，从而得到第l个脉冲N个频点分别对应的频谱值；

(2d)将所述N个频点分别对应的频谱值进行累加，得到第l个脉冲的频谱累加值S_l。

步骤3，令l加1，并依次重复步骤1至步骤2，从而得到雷达一个扫描周期内L个脉冲分别对应的频谱累加值。

步骤4，获取加权中值滤波线性窗，并根据雷达一个扫描周期内L个脉冲分别对应的频谱累加值计算每个脉冲的频谱累加值的中间值、频谱累加值的平均值，然后根据每个脉冲的频谱累加值的中间值、频谱累加值的平均值得到该脉冲对应的加权中值滤波后的值，从而得到L个脉冲对应的加权中值滤波后的值，并组成雷达一个扫描周期内L个脉冲的加权中值滤波结果。

参照图4，为加权中值滤波算法的流程示意图，步骤4具体包括如下子步骤：

(4a)获取长度为W的加权中值滤波线性窗，W＜＜L，＜＜表示远小于，W为奇数；

(4b)将雷达一个扫描周期内L个脉冲分别对应的频谱累加值依次记为S₁，...，S_l，...，S_L；

(4c)令l＝1，...，L，l初值为1；

(4d)令第l个脉冲对应的频谱累加值位于加权中值滤波线性窗的中间位置，得到加权中值滤波线性窗内包含的W个脉冲分别对应的频谱累加值，将加权中值滤波线性窗内包含的W个脉冲分别对应的频谱累加值进行排序，获取排序后的中间值，并将其作为第l个脉冲的频谱中间值G_l；并对加权中值滤波线性窗内包含的W个脉冲分别对应的频谱累加值求平均，将得到的平均值作为第l个脉冲的频谱累加值的平均值μ_l；

(4e)预设权重因子k1，得到第l个脉冲加权中值滤波后的值Q_l，Q_l＝k₁×G_l+(1-k1)×μ_l；

(4f)令l加1，并依次重复子步骤(4d)至子步骤(4e)，从而得到L个脉冲对应的加权中值滤波后的值，并组成雷达一个扫描周期内L个脉冲的加权中值滤波结果。

其中，k1∈[0，1]，k1由消除脉冲噪声和保持干扰边缘性能要求综合评估得到；特别的，当k1＝0时，干扰边缘保持性能达到最优；当k1＝1时，脉冲噪声可被最大程度消除，此为标准中值滤波。

特别的，如果加权中值滤波线性窗内包含的脉冲数小于W，则将加权中值滤波线性窗内空白窗口对应的频谱累加值设为零，再将加权中值滤波线性窗内包含的W个频谱累加值求平均，将得到的平均值作为位于加权中值滤波线性窗的中间位置的脉冲的频谱累加值的平均值。

中值滤波器对处理脉冲噪声非常有效，连续个数小于窗口宽度一半的离散脉冲将被滤除，而对于斜升或斜降的信号不产生影响，窗口可以根据冲激噪声空间密度进行调整，这样能够有效降低干扰指示的虚警率；另外加权中值滤波算法较标准中值滤波具有更优良的滤波性能，能更好地保留原始数据细节和边缘，并且在降低虚警率的同时也减少真实干扰的漏警。

步骤5，根据所述雷达一个扫描周期内L个脉冲的加权中值滤波结果，确定雷达一个扫描周期内的干扰区域集合，所述干扰区域集合包含J个干扰区域，J为自然数。

参照图5，为利用自适应门限划分干扰区域流程示意图，步骤5具体包括如下子步骤：

(5a)根据雷达一个扫描周期内L个脉冲分别对应的加权中值滤波后的值Q₁，...，Q_l，...Q_L，计算雷达一个扫描周期内L个脉冲分别对应的加权中值滤波后的值的均值μ和雷达一个扫描周期内L个脉冲分别对应的加权中值滤波后的值的标准差σ，其表达式分别为：

(5b)根据加权中值滤波后的值的均值μ和加权中值滤波后的值的标准差σ，计算雷达一个扫描周期内L个脉冲的干扰检测自适应门限T，T＝μ+k2×σ，并将其作为第0次迭代后的干扰检测自适应门限T₀，k2表示预设的系数，并设定干扰检测固定门限T2；

(5c)令l＝1，...，L，l初值为1；

(5d)将第l个脉冲加权中值滤波后的值分别与干扰检测自适应门限T和干扰检测固定门限T2比较，若所述第l个脉冲加权中值滤波后的值大于干扰检测自适应门限T，且大于干扰检测固定门限T2，则认为第l个脉冲所处方位存在干扰，并将第l个脉冲所处方位的干扰标志位置为1，否则将第l个脉冲所处方位的干扰标志位置为0；

(5e)令l的值加1，并重复执行子步骤(5d)，直到将第L个脉冲所处方位的干扰标志位置为1或将第L个脉冲所处方位的干扰标志位置为0；从而得到雷达一个扫描周期内L个脉冲对应的干扰标志位0/1序列；

(5f)预设干扰个数门限k3，令迭代次数m＝0，且m＝0，...，MAX，MAX为预设的最大迭代次数；

将雷达一个扫描周期内L个脉冲对应的干扰标志位0/1序列在方位上进行结果融合，目的是将“孤立点”剔除：

(5g)将雷达一个扫描周期内L个脉冲对应的干扰标志位0/1序列中连续0序列包含的0的个数小于预设个数门限k3的0序列全都置为1，进而得到融合后雷达一个扫描周期内L个脉冲的干扰标志位序列；将融合后雷达一个扫描周期内L个脉冲的干扰标志位序列中的各个连续1序列作为干扰区域，进而得到第m次迭代和融合后雷达一个扫描周期内L个脉冲的干扰标志位序列包含的D′个干扰区域；

(5h)设置干扰个数上限k4，所述干扰个数上限为干扰主瓣达到最大角度时包含的脉冲数，分别统计D′个干扰区域中每个干扰区域对应的脉冲个数，并将脉冲个数大于干扰个数上限k4的干扰区域归到第一集合D1，将脉冲个数小于干扰个数上限k4的干扰区域归到第二集合D2；

(5i)若第一集合D1为空集合，或者迭代次数m＝MAX，则执行子步骤(5p)；

(5j)设第一集合D1包含H个干扰区域，d∈{1，…，H}，d的初值为1；

(5k)获取第一集合D1包含的H个干扰区域中的第d个干扰区域，设第d个干扰区域内包含F_d个脉冲，获取所述F_d个脉冲各自对应的加权中值滤波后的值中的最大值，将所述F_d个脉冲各自对应的加权中值滤波后的值中的最大值与第m次迭代后的干扰检测自适应门限T_m求均值；令m的值加1，将该均值作为第m次迭代后的干扰检测自适应门限T_m；

(51)将第d个干扰区域内对应的F_d个脉冲各自加权中值滤波后的值分别与第m次迭代后的干扰检测自适应门限T_m和设置的干扰检测固定门限T2进行比较；

(5m)若第d个干扰区域内第k个脉冲加权中值滤波后的值大于第m次迭代后的干扰检测自适应门限T_m，且大于干扰检测固定门限T2，则认为第d个干扰区域内第k个脉冲所处方位存在干扰，并将第d个干扰区域内第k个脉冲所处方位的干扰标志位置为1；否则将第d个干扰区域内第k个脉冲所处方位的干扰标志位置为0；其中，k∈{1，…，F_d}，k的初值为1，F_d表示第d个干扰区域内对应的脉冲个数；

(5n)令k加1，重复执行子步骤(5m)，直到得到第d个干扰区域内第F_d个脉冲所处方位的干扰标志位置1或第d个干扰区域内第F_d个脉冲所处方位的干扰标志位置0；进而得到第d个干扰区域内的干扰标志位0/1序列；

(5o)令d加1，并重复执行子步骤(5k)至(5n)，直到得到第H个干扰区域内的干扰标志位0/1序列，并将此时得到的第1个干扰区域内的干扰标志位0/1序列到第H个干扰区域内的干扰标志位0/1序列，作为第m次迭代后雷达一个扫描周期内L个脉冲对应的干扰标志位0/1序列，并返回子步骤(5g)；

(5p)将第二集合D2包含的J个干扰区域，作为雷达一个扫描周期内的干扰区域集合。

步骤6，获取J个干扰区域内包含的脉冲分别对应的方位码，根据J个干扰区域内包含的脉冲分别对应的方位码、J个干扰区域内包含的脉冲分别对应的加权中值滤波后的值，确定J个干扰区域内包含的脉冲分别对应的干扰方位、干扰强度、干扰类型以及干扰频点。

需要补充的是，机械扫描雷达天线不断转动，每个脉冲对应转到某个方位，会有一个对应的方位码，该方位码指示了当前的方位。

步骤6具体包括如下子步骤：

(6a)设定门限T3，记雷达一个扫描周期内的干扰区域集合包含J个干扰区域，将第x个干扰区域内包含的脉冲个数记为L_x，x∈{1，…，J}，x的初值为1；

(6b)将第x个干扰区域的第c个脉冲N个频点分别对应的频谱值依次与设定的门限T3比较，如果第c个脉冲第n个频点对应的频谱值大于设定的门限T3，则认为第c个脉冲的第n个频点为第c个脉冲的干扰频点；c的初值为1，c＝1，…，L_x，n的初值为1，n依次取值1，…，N；

(6c)令c的值加1，并重复执行子步骤(6b)，进而依次得到第1个脉冲的干扰频点到第L_x个脉冲的干扰频点，并将第1个脉冲的干扰频点到第L_x个脉冲的干扰频点，作为第x个干扰区域的干扰频点C_x；

(6d)设定窄带干扰判决门限T4，计算第x个干扰区域的干扰频点C_x和第x个干扰区域对应的脉冲个数L_x的比值C_x/L_x，若比值C_x/L_x大于窄带干扰判决门限T4，则认为第x个干扰区域内存在宽带干扰，否则存在窄带干扰；

(6e)令x的值加1，并依次重复子步骤(6b)至(6d)，进而分别得到J个干扰区域各自存在宽带干扰或窄带干扰；

其中，J个干扰区域内的干扰频点功率最大值对应的频点，为窄带干扰频点；J个干扰区域内的干扰频点功率最小值和干扰频率功率最大值组成的区间范围，为宽带干扰频点。

(6f)获取J个干扰区域内包含的脉冲各自对应的方位码，分别记为F₁，…F_i′，…F_L′，L′为J个干扰区域内包含的脉冲个数，并根据J个干扰区域内包含的脉冲各自对应的加权中值滤波后的值分别为Q₁，…Q_i′，…Q_L′，计算得到J个干扰区域内包含的脉冲各自对应的干扰方位，其中将J个干扰区域中第l′个脉冲的干扰方位记为其表达式为：

其中[]表示取整运算。

根据一个雷达扫描周期的干扰区域集合包含J个干扰区域，计算得到J个干扰区域内的干扰强度为max(Q₁，…，Q_L′)，Q₁表示第1个脉冲加权中值滤波后的值，Q_L′表示第L′个脉冲加权中值滤波后的值；i′∈{1，…，L′}，L′为J个干扰区域内包含的脉冲个数。

下面结合MATLAB仿真实验对本发明效果做进一步的说明和分析

(一)仿真条件：

为方便起见，本发明采用雷达实际采集数据进行仿真，所述实际采集数据是经过专用数据采集器采集得到，并拷贝至PC机硬盘并保存，仿真输入数据长度为1个雷达扫描圈，包含3000个脉冲，每个脉冲数据为32个采样点。

(二)仿真内容：

按照上述步骤及仿真条件，在软件MATLAB 2015b中按照本发明进行仿真实验，简要说明如下：

1、将每个脉冲32点快速傅里叶变换结果累加，得到累加值，这样一个雷达扫描圈得到3000点的累加结果；

2、设置中值滤波线性窗长度W＝15，权重k1＝0.75。对3000点的累加结果进行中值滤波；

3、设置固定门限＝3e7，求取自适应门限，均值＝9.3494e6，标准差＝2.2845e7，所以自适应门限＝均值+标准差＝3.2194e7，固定门限＝3e7，利用中值滤波后的值与自适应门限和固定门限依次比较大小，同时大于两个门限的中值滤波结果其干扰标志位才会被置1。得到干扰标志位序列，不大于2个点的连续的0将会被置1，得到方位融合后的干扰标志位序列；

4、根据干扰标志位序列中连续1所在区域和预设的干扰区域内脉冲个数上限k4＝50，判决干扰区域。

5、在每个干扰区域，设置频点干扰检测门限T3＝1e7，设置宽窄带干扰判决门限T4＝4，进行频点干扰判决，判断干扰类型、干扰频点和干扰方位。

3、仿真结果分析：

1、如图6a，为中值滤波前频谱累加和的曲线图，即一个雷达扫描圈得到3000脉冲的快速傅里叶变换累加结果。从图6a中可以看出，原始累加结果在方位上有较大起伏，而且部分方位存在脉冲噪声干扰带来的“毛刺”；

2、如图6b，为中值滤波后频谱累加和的曲线图；从图6b中可以看出，异常的过大的值和过小的值被滤除，这样在后续干扰区域判决时可以有效减少干扰判决的虚警。

3、如图7a，为方位融合后的干扰标志位示意图；从图7a中可以看出，存在3段连续的连续1组成的序列，分别是437-470、2190-2227、2400-2498，其中437-470是干扰区域1；2190-2227是干扰区域2；2400-2498包含脉冲个数大于预设个数门限50，所以应该包含两个或两个以上的干扰区域，所以更新自适应门限，重置标志位，直到所有干扰区域内包含脉冲数小于个数门限。

4、如图7b，为更新后的自适应门限示意图，其值为1.405e8，进行2400-2498区域内干扰标志位重置，最终划分得到干扰区域为2432-2458、2471-2488；

5、如图8，为某脉冲32频点功率分布示意图，即干扰区域中某脉冲的频点功率值；从图8中可以看出，干扰集中在26-28这3个频点上。以干扰区域437-470为例，包含63个过门限的频点，脉冲数为34，平均每个脉冲包含的干扰频点数为63/34＝1.85＜4，干扰类型为窄带干扰，干扰强度＝6.9272e7，干扰频点为第26频点，干扰方位为454。

综上所述，仿真实验验证了本发明的正确性，有效性和可靠性。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。