间歇采样转发干扰的时频域识别方法
2021-12-21余康林匡华星王超宇
余康林,匡华星,王超宇
(中国船舶集团第七二四研究所,南京 211153)
引 言
以数字射频存储技术(Digital radio frequency memory,DRFM)为基础的间歇采样转发干扰,其干扰策略最早由Wang等[1]和Sparrow等[2]分别提出,目前典型的间歇采样转发干扰样式主要有间歇采样直接转发干扰、间歇采样重复转发干扰和间歇采样循环转发干扰。由于间歇采样转发干扰与目标回波在时域和频域上存在较强相干性,匹配滤波后形成大量逼真假目标,传统干扰识别方法难以辨识间歇采样转发干扰,使得雷达无法有效检测到目标。针对间歇采样转发干扰识别率较低的问题,亟需研究具有较强识别能力的干扰识别方法。
目前间歇采样转发干扰识别方法的研究热点主要集中于Wigner⁃Ville变换、分数阶傅里叶变换[3]、模糊函数[4]和短时傅里叶变换(Short⁃time Fourier transform,STFT)等时频分析方法。文献[5]提出了一种基于脉冲压缩和STFT的特征提取方法,该方法先对雷达回波进行脉冲压缩,再进行STFT,提取距离维谱峰个数作为特征参数区分间歇采样转发干扰,该方法在较高干噪比时能够辨识几类间歇采样转发干扰,但在较低干噪比时识别能力不足。文献[6]提出了一种基于短时分数阶博里叶变换(Short⁃time fractional Fourier transform,STFRFT)的间歇采样转发干扰识别方法,该方法在对雷达回波进行STFRFT后,在STFRFT域提取切片个数、切片宽度、中心频率和调频斜率作为特征参数进行干扰识别,但该方法搜索复杂度较高。文献[7]提出了基于频谱相像系数的间歇采样转发干扰识别方法,该方法提取了频谱矩形相像系数和三角相像系数作为干扰识别的特征参数,但需提前训练分类器,实时性较差。文献[8]提出了一种基于奇异值分解的间歇采样转发干扰识别方法,该方法在对干扰进行奇异值分解后,提取奇异谱熵函数作为干扰识别的特征参数,能在较低干噪比下仍具有较好的识别效果,但对滑窗宽度的选取较为敏感。
针对现有识别方法存在的问题,本文在干扰建模的基础上,通过推导和分析,论证了几类间歇采样转发干扰二维时频分布的差异,在此基础上提出了一种基于STFRFT降维的间歇采样转发干扰识别方法,详述了特征参数提取方法和干扰识别流程,并将其时频分布特征总结为先验知识,无需训练分类器即可对几类间歇采样转发干扰进行分类识别。最后还分析了当间歇采样转发干扰切片宽度发生变化时,所提方法识别率的变化情况。
1 间歇采样转发干扰时域模型
间歇采样转发干扰(Interrupted⁃sampling and repeater jamming,ISRJ)的工作原理是干扰机截取、储存雷达发射信号的某一片段后,对该片段进行一次或多次转发,重复截取⁃储存⁃转发过程,直至脉冲结束。上文提及的3种间歇采样转发干扰的区别在于转发策略不同,工作原理图可参考文献[9]。间歇采样直接转发干扰与间歇采样重复转发干扰的差异在于干扰机对当前截取的雷达发射信号片段的转发次数不同,间歇采样直接转发干扰对当前截取的信号片段只转发一次,而间歇采样重复转发干扰多次转发当前截取的信号片段;间歇采样循环转发干扰与两者的不同在于:间歇采样循环转发干扰在转发完当前截取的信号片段后,还会逆序转发之前截取的全部信号片段。自卫式干扰场景下,雷达回波可描述为
式中:s(t)、j(t)和n(t)分别为真实目标回波、干扰和噪声。
对于脉冲压缩雷达,不考虑载频,假设雷达发射信号时域模型为
式中:rect(⋅)为矩形窗函数;T为脉宽;k为雷达发射信号调频斜率
3类间歇采样转发干扰对应时域模型如下:
(1)间歇采样直接转发干扰(Interrupted⁃sampling and direct repeater jamming,ISDJ),其工作原理图如图1所示,对应时域模型为
图1 间歇采样直接转发干扰Fig.1 Interrupted-sampling and direct repeater jamming
式中:M为切片个数,TJ为切片宽度。
(2)间歇采样重复转发干扰(Interrupted⁃sampling and periodic repeater jamming,ISPJ),其工作原理图如图2所示,对应时域模型为
图2 间歇采样重复转发干扰Fig.2 Interrupted-sampling and periodic repeat⁃er jamming
式中:N为每个切片的转发次数,α(m,n)=(m-1)(N+1)+n为第m次切片进行第n次转发对应的时延系数。
(3)间歇采样循环转发干扰(Interrupted⁃sampling and cyclic repeater jamming,ISCJ),其工作原理图如图3所示,对应时域模型为
图3 间歇采样循环转发干扰Fig.3 Interrupted-sampling and cyclic repeater jamming
式中:α(m)=m(m+1)/2-1为第m次切片对应的时延系数,β(m,n)=n(n+1)/2+m(n-1)为第m个切片进行第n次转发对应的时延系数。
2 间歇采样转发干扰时频特征分析
对于间歇采样转发干扰,时延相同的切片在FRFT域对应同一个峰值点。由上述干扰信号时域模型可知,ISDJ每次转发的切片时延相同,因此其在FRFT域仅有一个峰值;而ISPJ对每个切片进行多次转发,因此其在FRFT域具有多个峰值,并且峰值个数等于截取次数;同样ISCJ在FRFT域也具有多个峰值,其沿第一个峰值点所在时间轴切片的脉冲个数等于截取次数。下面通过FRFT推导来验证上述分析。
式中k为雷达发射信号调频斜率。
下面对ISDJ、ISPJ和ISCJ进行FRFT推导,以间歇采样直接转发干扰为例,其FRFT为
由式(9)可得间歇采样直接转发干扰FRFT系数归一化幅度谱近似为|Xp(u)|≈|sinc(πTJ(ucscα+k TJ))|,ISDJ在u0=-k TJsinα时取得最大值,在FRFT域存在一个峰值,并且主瓣宽度为2|sinα|/TJ。
同理,间歇采样重复转发干扰FRFT系数归一化幅度谱近似为
其在FRFT域的峰值点为u0=-n k TJsinα,观察式(9)和式(10)可以发现间歇采样重复转发干扰的FRFT是对间歇采样直接转发干扰FRFT的结果在FRFT域进行多次频移得到,其在FRFT域具有多个峰值点。
间歇采样循环转发干扰FRFT系数归一化幅度谱近似为
其在FRFT域的峰值点为u0=-β(m,n)k TJsinα,在FRFT域也存在多个峰值点。
由式(9~11)可以看出,对间歇采样转发干扰FRFT后,ISDJ在FRFT域表现为仅有一个峰值,而ISPJ和ISCJ在FRFT域有多个峰值,因此可以将峰值个数作为识别ISDJ的特征参数。但由于3类间歇采样转发干扰FRFT的结果为sinc函数,会产生能量较大的旁瓣,在进行谱峰搜索时,会造成严重的虚警,并且仅从FRFT域提取特征参数无法对ISPJ和ISCJ进行有效区分。
为了克服在FRFT域提取特征参数存在的问题,同时进一步区分ISPJ和ISCJ,本文采取了基于短时滑窗的FRFT分析方法,在对间歇采样转发干扰采取时频域截断的同时,通过加窗处理,可以保证间歇采样转发干扰在STFRFT域有较强能量聚集性的同时具有较低的旁瓣,并且能够得到间歇采样转发干扰的局部时频信息,为间歇采样转发干扰的识别提供更多信息支持。
工程应用中是对采样后的信号进行离散短时分数阶傅里叶变换(Discrete short⁃time fractional Fou⁃rier transform,DSTFRFT),而DSTFRFT实际上是对滑窗截取的每个信号片段进行离散分数阶傅里叶变换(Discrete fractional Fourier transform,DFRFT)。下面介绍DSTFRFT的算法实现流程和基于STFRFT域的间歇采样转发干扰识别方法。
3 基于STFRFT域的间歇采样转发干扰识别
3.1 DSTFRFT的算法实现流程
对于函数x(t),其p阶STFRFT定义为
式中:g(t)为窗函数,核函数Kp(t,u)的定义如式(7)所示。在实际使用中需要对STFRFT进行离散化处理,对于离散信号x(l),其DSTFRFT定义为
式中:x(n)为滑窗截取的信号片段,n'=(i-1)Δ,N为滑窗长度,Δ为滑窗间隔。
核函数Kp(n,m)表达式为
式中:Δu为FRFT域采样间隔,Δt为时域采样间隔。由STFRFT的定义可知,对于信号x(t)的ST⁃FRFT实际上是将滑窗截取的每个信号片段进行FRFT,工程应用上,为了实现DFRFT的快速计算,一般将DFRFT转换为快速傅里叶变换(Fast Fourier transform,FFT),其算法实现主要包括以下几个步骤:
(1)窗函数、滑窗长度和滑窗间隔的确定。由文献[10]可知,高斯窗函数在短时分数阶傅里叶变换域仍具有最小的时宽带宽积,因此为了具有较好的时频分辨率,本文窗函数选取高斯窗。为了便于推导,这里以单分量LFM信号为例,分析高斯窗对其频率分辨率的影响。加窗LFM信号的FRFT变换为
式中σ为标准差。可得加窗LFM信号的FRFT变换等效于高斯函数的傅里叶变换,因此加窗LFM信号FRFT的幅度为
可以看出加窗LFM信号FRFT后的幅度分布仍为高斯函数,其主瓣宽度受σ控制,因此可以通过调节σ来改善频率分辨率,并且其主瓣宽度与LFM信号调频率无关,因此具有较好的鲁棒性。
为了在有限的信号观测时间范围内更好地描述信号的时频特性[11],参考LFM信号STFT的滑窗长度选取准则[12],对于连续多项式函数φ(tm),假设其(n+2)阶可导,对于STFRFT,滑窗最佳时宽为
由于对信号进行加窗处理,会导致信号产生能量损失,因此需要使用重叠相加法来降低信号的能量损失,同时,为了降低DSTFRFT的算法复杂度,滑窗间隔一般选取为滑窗长度的0.1~0.3倍[13]。
(2)离散傅里叶变换算子构造。规定Δt=1/fs为时域采样间隔,fs为采样率,为了构造离散傅里叶变换算子,则Δu与Δt应满足ΔuΔt=|sinα|/M,M为离散傅里叶变换的区间长度,则FRFT域采样间隔Δu=fs⋅|sinα|/M。将Δu与Δt代入式(13),可得
因此通过构造离散傅里叶变换算子,DSTFRFT就转换为对x'(n)进行离散傅里叶变换。
(3)时域chirp乘法运算。将滑窗截取的信号片段x(n)乘以exp(iπn2Δt2cotα)。
(4)根据sinα的符号进行相应处理。若sinα>0,则对x'(n)进行傅里叶变换;若sinα<0,则进行逆傅里叶变换。
(5)频域chirp乘法运算。将滑窗截取的信号片段变换后的结果乘以Aαexp(iπm2Δu2cotα),遍历所有信号片段,即得到离散信号x(l)的DSTFRFT。
3.2 基于STFRFT域的干扰识别流程
对于间歇采样转发干扰,时延相同的切片FRFT后在FRFT域对应同一个峰值。通过上述分析可知,ISDJ在FRFT域仅有一个峰值,沿峰值点所在时间轴的脉冲个数等于截取次数;ISPJ在FRFT域具有多个峰值,峰值个数等于转发次数,每个峰值点在时间轴的脉冲个数等于截取次数;对于ISCJ,其在FRFT域也具有多个峰值,除了第一个峰值点沿时间轴的切片包含多个脉冲,其余峰值点沿时间轴的切片脉冲个数仅有一个,并且第一个峰值点对应的脉冲个数等于截取次数,因此可统计谱峰个数和脉冲个数作为识别间歇采样转发干扰的特征参数。特征提取与识别流程如图4所示。
图4 间歇采样转发干扰识别流程图Fig.4 Recognition flow chart of interrupted-sampling and repeater jamming
包括以下几个步骤:
(1)对STFRFT后的时频图进行降维处理。降维方法如下:对于FRFT域波形,每个频点的幅值对应着该频点所在时间轴的最大值,遍历所有频点,即得到该干扰FRFT域波形;对于时域波形,在获得干扰FRFT域峰值点位置后,沿峰值点所在时间轴对二维时频图进行切片,即得到该峰值点在时域的脉冲分布情况。
(2)对降维后的数据进行标准化和归一化。标准化和归一化可以降低数据中较大噪声和异常值的影响,提高干扰与噪声的区分度。
(3)自适应门限设置。由于干扰在FRFT域具有较高的能量聚集性,在FRFT域形成很窄的尖峰,因此干扰与噪声存在较大的区分度。假设归一化后FRFT域的幅度为X1,X2,…,XN,考虑3种门限,分别为FRFT域幅度谱加权平均值、均值和最大值向下取3 dB,对比3种门限下干扰识别方法识别率的变化。其中加权平均值定义为
对于时域波形,在干噪比较高的情况下,每个峰值点所在时间轴的脉冲幅度较为接近,因此可将每个峰值点幅度向下取3 d B作为门限统计脉冲个数。
(4)间歇采样转发干扰在FRFT域与时域的特征参数统计。在FRFT域对干扰进行过门限谱峰搜索,统计谱峰个数,若峰值个数为1,则为间歇采样直接转发干扰,否则为另两种干扰;相比于间歇采样重复转发干扰,间歇采样循环转发干扰第一个峰值对应的脉冲个数与剩余峰值对应的脉冲个数不相等,为了降低搜索复杂度,本文取干扰的第一个峰值和最后一个峰值进行脉冲个数统计,若两峰值点对应的脉冲个数相等,则为间歇采样重复转发干扰,否则为间歇采样循环转发干扰。
最后,对本文干扰识别方法特征参数提取时的运算及搜索复杂度进行分析。假设回波时频变换后对应时频矩阵大小为M×N,则二维时频变换需要N Mlog2M次运算。将二维时频矩阵降维至FRFT域需要M(N-1)次最大值搜索,FRFT域谱峰个数提取需要M-2次搜索,时域脉冲个数提取需要2N次搜索,总搜索次数约M N+2N-2;文献[6]中时频矩阵二值化处理需要2M N次运算,其中最大值搜索需M N次,二值化运算需M N次,从二值化后的时频矩阵中提取切片宽度、切片个数等特征参数时需要至少M N+16N次搜索,总搜索次数为3M N+16N。因此本文干扰识别方法搜索次数与文献[6]相比搜索次数可减少约70%。
为了验证所提方法的有效性,下面给出了几组仿真实验对所提方法进行验证。
4 仿真结果及分析
本节通过几组仿真实验对上述间歇采样转发干扰识别方法进行验证,并且通过对比分析了识别方法的性能。表1给出了仿真中使用的主要参数,采样率设为150 MHz,干信比(Jamming⁃to⁃sig⁃nal ratio,JSR)设为0~20 dB,干噪比(Jamming⁃to⁃noise ratio,JNR)设 为-10~10 d B,信 噪 比(Signal⁃to⁃noise ratio,SNR)设为-25~10 d B。由于雷达接收机在采样得到雷达回波后,一般先下变频将高频信号解调为基带信号,再对雷达回波进行后续处理,因此为了方便分析,这里载频设为0 MHz。
表1 仿真参数设置Table 1 Simulation parameter setting
图5~7分别给出了切片宽度2.5μs、干噪比-3 dB时,上述仿真条件下3种间歇采样转发干扰二维时频分布降维到FRFT域和时域对应的幅度图,此时FRFT对应的sinα为4×10-13。图5为3类间歇采样转发干扰时频分布降维后的FRFT域幅度谱,从图中可以看出间歇采样直接转发干扰FRFT域幅度谱具有一个谱峰,而间歇采样重复转发干扰和间歇采样循环转发干扰FRFT域幅度谱具有多个谱峰,与上述分析一致,可将谱峰个数作为区分间歇采样直接转发干扰和另两类间歇采样转发干扰的特征参数。图6和图7分别为间歇采样重复转发干扰和间歇采样循环转发干扰沿峰值点所在时间轴切片后得到的脉冲分布幅度谱,图6、7中可以看出间歇采样重复转发干扰沿每个峰值点所在时间轴切片后对应的脉冲个数相等,而间歇采样循环转发干扰沿第一个峰值点切片后对应脉冲个数与剩余峰值点对应脉冲个数不等,因此可将该特征作为区分间歇采样重复转发干扰和间歇采样循环转发干扰的特征参数。由图5~7可以看出,3类间歇采样转发干扰在时频域分布特征不同,因此可从时频域区分三者。下面通过仿真实验对上述干扰识别方法进行验证。
图5 间歇采样转发干扰FRFT幅度图Fig.5 FRFT amplitude spectrum of interrupted-sampling and repeater jamming
图6 ISPJ两峰值点沿时间轴的脉冲分布幅度图Fig.6 Time domain amplitude of two peaks of ISPJ
图7 ISCJ两峰值点沿时间轴的脉冲分布幅度图Fig.7 Time domain amplitude of two peaks of ISCJ
仿真实验1
令干信比为15 d B,干噪比为-10~10 dB,干扰参数设置如表1所示,仿真步骤如下:
(1)设定切片宽度为2.5μs,对某一固定的干噪比,利用干扰模型生成对应的数据样本;
(2)利用上述干扰识别算法对生成的数据样本进行处理,并进行干扰识别,若识别结果与干扰类型一致,则表示此次识别结果正确;
(3)重复步骤1~2,对每种干扰进行200次Monte Carlo实验,记录该干噪比下每种干扰的识别率;
(4)改变干噪比,重复步骤1~3,绘制每种干扰不同干噪比下的识别率变化曲线;
(5)改变切片宽度,重复步骤1~4,绘制不同切片宽度下每种干扰对应的识别率变化曲线,最后绘制每种干扰不同干噪比下的平均识别率。
图8为间歇采样转发干扰在上述干扰识别方法下的识别结果,其中图8(a、b、c)为FRFT域阈值取FRFT域幅度谱加权平均值对应结果。由于间歇采样直接转发干扰降维后的FRFT域幅度谱在不同切片宽度下都只包含一个谱峰,相同干噪比下,降维后的FRFT域幅度分布较为接近,由图8(a)可以看出,不同切片宽度下,间歇采样直接转发干扰的识别率虽然存在波动,但是波动范围较小,因此切片宽度的变化对识别间歇采样直接转发干扰的影响较小,在干噪比大于-4 dB时,识别率能够达到90%以上;对于间歇采样重复转发干扰,转发次数不变时,切片宽度减小将会导致其沿每个峰值点所在时间轴切片包含的脉冲个数增加,受噪声的起伏影响,低干噪比下,在统计每个峰值点对应脉冲个数时正确概率将降低,因此切片宽度2.5μs时的识别率大于切片宽度2μs和1μs时的识别率,并且在干噪比大于-6 dB时,识别率可达到90%以上;对于间歇采样循环转发干扰,切片宽度减小将会导致第一个峰值点与剩余峰值点对应脉冲个数差距变大,因此在低干噪比下,切片宽度为1μs时的识别率高于切片宽度2.5μs和2μs时的识别率,并且在干噪比大于-8 d B时识别率能达到90%以上。
同时,对降维后的FRFT域幅度谱取不同阈值时几类干扰平均识别率的变化进行了分析,并与文献[6]算法进行对比。识别结果如图8(d)所示,其中阈值1为取FRFT域幅度谱加权平均值,阈值2为取FRFT域幅度谱最大值向下取3 d B,阈值3为取FRFT域幅度谱均值。由图8(d)可知,当阈值设为FRFT域幅度谱最大值向下3 dB时识别效果最好,在干噪比为-6 d B时,识别率接近100%,并且在低干噪比下,该识别方法识别率高于文献[6]算法。这是因为STFRFT对LFM信号具有能量聚集性,而噪声在二维时频面上是随机分布的,因此在较低干噪比下干扰的峰值幅度也远高于噪声幅度,干扰与噪声之间具有较高的区分度,通过选择合适的阈值就能区分干扰与噪声。但阈值设置过低时,将无法有效区分干扰和噪声;阈值设置过高时,将会导致干扰信息丢失。低干噪比下,阈值取为最大值向下3 d B时的值将高于FRFT幅度谱平均值和FRFT幅度谱加权平均值,并且不会导致干扰信息丢失,因此识别效果最好。
图8 间歇采样转发干扰不同干噪比下识别率Fig.8 Recognition rate of ISRJ at different JNRs
仿真实验2
令干信比为15 dB,信噪比为-25~10 d B,干扰仿真参数如表1所示,几类干扰不同信噪比下的识别率如图9所示。由图9可知,当FRFT域阈值为最大值向下取3 d B时,在低信噪比下识别率高于文献[6]算法,并高于另两种阈值;但是当阈值取FRFT域幅度谱均值时无法对几类干扰进行有效识别。
图9 间歇采样转发干扰不同信噪比下平均识别率Fig.9 Recognition rate of ISRJ at different SNRs
令信噪比为0 dB,干信比为0~20 dB,干扰仿真参数如表1所示,几类干扰不同干信比下识别率如图10所示。由图10可知,当FRFT域阈值取最大值向下3 d B时,干扰识别率高于另两种阈值,在低干信比下,识别率高于文献[6]算法。
图10 间歇采样转发干扰不同干信比下平均识别率Fig.10 Recognition rate of ISRJ at different JSRs
综上所述,3类间歇采样转发干扰的识别率主要与干噪比和干信比相关,并且3种干扰对应识别率随着干噪比和干信比的升高逐渐提高;当FRFT域阈值取FRFT域幅度谱最大值向下3 d B时,干扰识别效果最好,在低干信比和低干噪比下,识别率高于文献[6]算法,但取均值时无法有效识别几类干扰。
5 结束语
该文针对间歇采样转发干扰识别率较低的问题,提出了一种基于STFRFT的间歇采样转发干扰识别方法,该方法结合了FRFT对LFM信号的强能量聚集性和STFT加窗处理对旁瓣的抑制以及对频率分辨率改善的优势。同时为了降低特征提取的复杂度,本文通过二维时频降维处理,在降维后的FRFT域和时域统计特征参数,达到了降低搜索复杂度的目的,并且通过对特征参数的分析还能获知间歇采样转发干扰的干扰策略。当切片宽度发生变化时,在谱峰可分的情况下,上述算法同样有效。仿真结果表明,所提干扰识别方法对噪声不敏感,在低干噪比下也能较好辨识几类间歇采样转发干扰。并且由于干扰与真实目标回波的时频能量分布在不同的条带中,在干扰样式识别的基础上,可以通过剔除干扰对应时频能量点达到抑制干扰的目的,再通过逆STFRFT得到干扰抑制后的雷达回波。