Loading [MathJax]/extensions/TeX/boldsymbol.js

基于张量结构的快速三维稀疏贝叶斯学习STAP方法

崔宁 行坤 段克清 喻忠军

王晓戈, 李槟槟, 陈辉, 等. 基于脉内频率编码联合调频斜率捷变波形的ISRJ对抗方法[J]. 雷达学报(中英文), 2024, 13(5): 1019–1036. doi: 10.12000/JR24046
引用本文: 崔宁, 行坤, 段克清, 等. 基于张量结构的快速三维稀疏贝叶斯学习STAP方法[J]. 雷达学报, 2021, 10(6): 919–928. doi: 10.12000/JR21140
WANG Xiaoge, LI Binbin, CHEN Hui, et al. Anti-ISRJ method based on intrapulse frequency-coded joint frequency modulation slope agile radar waveform[J]. Journal of Radars, 2024, 13(5): 1019–1036. doi: 10.12000/JR24046
Citation: CUI Ning, XING Kun, DUAN Keqing, et al. Fast tensor-based three-dimensional sparse Bayesian learning space-time adaptive processing method[J]. Journal of Radars, 2021, 10(6): 919–928. doi: 10.12000/JR21140

基于张量结构的快速三维稀疏贝叶斯学习STAP方法

DOI: 10.12000/JR21140
基金项目: 国家自然科学基金(61871397)
详细信息
    作者简介:

    崔 宁(1995–),男,辽宁抚顺人,中国科学院大学在读博士研究生,主要研究方向为空时自适应处理和雷达动目标检测等

    行 坤(1980–),男,陕西西安人,中国科学院空天信息创新研究院副研究员,中国科学院大学硕士生导师,主要研究方向为机载多功能雷达、地基监视雷达系统总体设计与性能仿真,地面、海面、空中动目标检测与跟踪技术等

    段克清(1981–),男,河北石家庄人,中山大学副教授,硕士生导师,主要研究方向为阵列信号处理、空时自适应处理、压缩感知、深度学习及其在雷达系统中的应用

    喻忠军(1980–),男,四川成都人,中国科学院空天信息创新研究院研究员,中国科学院大学博士生导师,主要研究方向为微系统微集成、微波毫米波电路模块、先进相控阵天馈技术等

    通讯作者:

    段克清 duankeqing@aliyun.com

    喻忠军 yuzj@ucas.ac.cn

  • 责任主编:冯大政 Corresponding Editor: FENG Dazheng
  • 中图分类号: TN957.51

Fast Tensor-based Three-dimensional Sparse Bayesian Learning Space-Time Adaptive Processing Method

Funds: The National Natural Science Foundation of China (61871397)
More Information
  • 摘要: 当机载雷达处于非正侧视工作模式时,非平稳杂波会对运动目标检测造成严重干扰。传统三维空时自适应处理(3D-STAP)方法通过构造俯仰-方位-多普勒三维自适应滤波器,可有效抑制非平稳杂波,然而巨大的系统自由度导致其在非均匀杂波环境下训练样本严重不足。虽然稀疏恢复(SR)技术可有效改善样本需求,但庞大的运算开销又使得该技术难以应用于实际。针对上述问题,该文结合机载雷达回3阶张量结构提出一种新的快速三维稀疏贝叶斯学习STAP方法,通过采用运算开销更低的张量处理将大规模矩阵求解拆分为多个小规模矩阵计算,从而大幅降低运算复杂度。详尽的数值实验验证了所提张量基SR-STAP方法可在维持SR-STAP小样本处理性能不变的基础上,将运行时间直接降低数个量级,因此是一种更适用于实际工程的SR-STAP处理方式。

     

  • 相参干扰是现代电子对抗的重要趋势之一[1,2]。相参干扰能够同目标回波一样获得匹配滤波和相参积累增益,因此能够以较小的功率代价在目标周围形成能量高的假目标,严重影响了雷达对目标的检测能力。其中,基于数字射频存储技术[3]的间歇采样转发干扰(Interrupted Sampling Repeater Jamming, ISRJ),克服了高速采样和收发天线高度隔离的难题,能被应用于导弹这类使用收发分时天线的小型干扰机平台[4],在电子对抗领域得到了广泛的应用。

    当前,关于雷达的目标检测[5,6]和抗干扰技术层出不穷[79],然而ISRJ是依据干扰机对雷达信号间歇采样和转发形成的一种脉内的相参干扰。因此,传统脉间频率捷变[7],脉组频率捷变,相位调制[8]等技术难以对抗ISRJ。此外,ISRJ通常被用于自卫式干扰[1],这导致传统的空域技术,例如旁瓣对消和旁瓣匿影技术[9]也失效。目前,针对ISRJ的抗干扰技术大致可分为两个方面:被动抗干扰和主被动联合抗干扰。被动抗干扰是指在干扰进入接收机后通过信号处理的各种手段对干扰进行抑制。主被动联合抗干扰是指发射端通过波形捷变等方式降低雷达信号被干扰机截获概率,同时结合接收端信号处理达到抗干扰目的。

    在被动抗干扰方面,雷达发射的信号为线性调频(Linear Frequency Modulation, LFM)信号。干扰抑制方法主要包含滤波和干扰重构。滤波方法主要依据干扰与目标在能量上的差异对干扰进行抑制。张建中等人[10]将LFM信号划分为多个子段,依据含干扰子段与无干扰子段脉压后能量不同,提取无干扰子段。万鹏程等人[11]利用发射信号调频斜率这个先验信息,将回波信号变换到分数阶域对目标进行提取,然而该方法适用于较低的干信比(Jamming-to-Signal Ratio, JSR),当JSR较大时,提取的目标中有大量的干扰剩余,脉压后会产生距离旁瓣干扰雷达检测。文献[1214]中,通过短时傅里叶变换(Short-Time Fourier Transform, STFT)将回波变换到时频域并根据能量的分布构造带通滤波器对干扰进行抑制,但回波信号的完整性会受到影响,且目标能量有较大损失。盖季妤等人[15]利用目标回波和干扰信号在差分特征空间的差异,抑制干扰的同时实现目标检测,相较于时频滤波算法提高了输出信噪比(Signal-to-Noise Ratio, SNR)。在干扰重构方面,Zhou等人[16]利用回波脉压后的时频分布估计采样次数和转发次数,利用解卷积估计切片宽度进而重构干扰并对消,然而STFT分辨率不高,降低了干扰参数的估计精度。刘孟斐等人[17]利用分段信息熵对ISRJ的关键参数进行估计并判断干扰类型,但是该方法需要对回波信号进行大量分段,计算量大。Lu等人[18]则构造不同长度的匹配滤波器与回波信号进行匹配滤波,依据输出结果估计干扰参数并重构干扰,该方法在低SNR下有着较好的抑制效果,但算法复杂度高,且十分依赖干扰参数的准确估计。

    在主被动联合抗干扰方面,张建中等人[19]提出了基于脉内载频步进LFM信号的抗ISRJ方法,利用子脉冲载频步进使干扰采样后无法干扰相邻子脉冲,从而提取无干扰信号段并积累,在对抗噪声调制干扰时有明显优势。Zhou等人[20,21]则采用发射波形和滤波器联合设计的方法来抑制ISRJ,但该方法依赖间歇采样信号的周期和干扰占空比的精确估计。类似的,Wang等人[22]研究了互补序列和接收机滤波器联合设计。西安电子科技大学的吴耀君等人[2325]基于频率编码信号提出了众多ISRJ对抗方法,其中文献[23]采用聚类算法识别受干扰回波切片,之后利用目标与干扰在时频维的差异进行频率维窄带滤波提取目标,最后分段脉压恢复信号。但干扰和目标在频率维均有较大展宽,在干扰切片较多时,目标与干扰在部分切片中难以分离。文献[24]则利用干扰与目标中心频率的差异在分数阶域对干扰进行抑制,缺点是当中心频率差异较小和JSR较大时,会有较多干扰剩余。类似的,文献[25]从回波信号时频分布出发,将时频分布矩阵沿频率维叠加后,提取回波中无干扰段构造时域滤波器,然后直接对回波脉压输出进行滤波。不足之处是经过处理后一些距离单元的强干扰仍有剩余,当JSR较高时,仍可能在能量上对真实目标进行压制。此外,牛闯等人[26]在脉内频率编码基础上,引入脉内时延捷变,后续采用多级滤波抑制干扰。虽然脉内时延捷变提高了目标与干扰在时频域的区分度,但也降低了发射信号的能量,减小了雷达的探测距离。

    总的来说,被动抗干扰面临着干扰抑制后目标信号完整度不高或参数估计精度要求高,计算量大的问题。而主被动联合抗干扰参数变化灵活,雷达波形被干扰机准确截获概率低,在对抗脉内干扰时有着独特的波形对抗优势。

    受文献[2326]启发,本文从主被动联合抗干扰的角度出发,针对脉内频率编码信号在时频域能量集中度不高,在干扰机采样次数和转发次数均较大时,部分回波切片中存在干扰与目标在时频域耦合导致干扰难以剔除的问题,提出了一种基于脉内频率编码联合调频斜率捷变波形的抗ISRJ方法,利用子脉冲中心频率和调频斜率双参数捷变实现子脉冲相互掩护,利用分数阶傅里叶变换(Fractional Fourier Transform, FrFT)对LFM信号的能量聚焦性和ISRJ的时域不连续性对回波信号进行分数阶域和时域级联滤波。本文的创新点和贡献如下:

    (1) 设计了一种频率编码联合调频斜率捷变的发射波形。该波形提高了目标与干扰在分数阶域的可分性,同时双参数捷变提高了子脉冲间的相互掩护能力并降低了脉压后的旁瓣幅度。

    (2) 干扰抑制过程采用分数阶域-时域联合滤波。克服了单域信号处理难以兼顾较小目标能量损失和较低距离旁瓣的问题。

    (3) 所提方法在应用场景上对干扰机同步采样和非同步采样形成的ISRJ均具有良好的抑制效果。

    干扰机采用收发共用天线用以截获雷达信号并产生ISRJ,故而经过干扰机交替采样、转发形成的干扰波形总是滞后于当前目标回波子脉冲的波形。对此,本文从波形域的角度出发,设计了一种脉内频率编码联合调频斜率捷变的雷达发射波形来对抗ISRJ,其时频示意图如图1所示。

    图  1  脉内频率编码联合调频斜率捷变信号时频图
    Figure  1.  Time-frequency diagram of intrapulse frequency-coded joint FM slope agile waveform

    我们定义子脉冲间的相互掩护能力是干扰机采样转发前一个子脉冲后,后面的子脉冲不受干扰影响的能力。若子脉冲间相关性弱,则子脉冲间的相互掩护能力强。若子脉冲间相关性强,则子脉冲间的相互掩护能力弱。

    图1中,横坐标表示时间,纵坐标表示信号频率,雷达发射信号由多个子脉冲拼接而成。子脉冲均为LFM信号且具有相同的脉宽。不同之处是各子脉冲的中心频率和带宽均不同。通过控制带宽实现子脉冲调频斜率捷变。此外,相邻子脉冲采用正负调频斜率调制,进一步提高了子脉冲间的相互掩护能力。

    假设雷达发射信号的脉冲宽度为Tp,带宽范围为BminBmax。将发射脉冲均匀划分为N个子脉冲,则子脉冲宽度为Tsub=Tp/N,第i个子脉冲的带宽为Bi=Bmin+(a(i)1)ΔB, ΔB=(BmaxBmin)/(N1)为子脉冲最小带宽间隔。a(i){1,2,,N}为子脉冲带宽编码序列,且ij时,a(i)a(j)。相邻子脉冲采用正负调频斜率调制,则第i个子脉冲的调频斜率为ki=(1)i1Bi/Tsub。第i个子脉冲相对初始中心频率的移频量为fi=(b(i)1)ΔfB0/2, Δf=B0/(N1)为子脉冲最小频率间隔,B0为子脉冲载频最大间隔且B0的取值应较大,旨在减小相邻脉冲的相关性。b(i){1,2,,N}为子脉冲频率编码序列。

    雷达发射的脉内频率编码联合调频斜率捷变波形的基带信号为

    s(t)=Ni=1si(t)=Ni=1rect[t(2i1)Tsub/2Tsub]exp{j2πfi[t(2i1)Tsub/2]+jπki[t(2i1)Tsub/2]2} (1)

    其中,si(t)为发射的第i个子脉冲,矩形窗函数为

    rect(u)={1,|u|1/20,

    干扰机形成ISRJ的过程如图2所示。

    图  2  ISRJ机理
    Figure  2.  Mechanism of ISRJ

    图2中,根据转发次数和转发方式的不同,ISRJ干扰可以分为3种类型[18]:间歇采样直接转发干扰(Interrupted Sampling Direct Repeater Jamming, ISDRJ),间歇采样重复转发干扰(Interrupted Sampling Periodic Repeater Jamming, ISPRJ)和间歇采样循环转发干扰(Interrupted Sampling Cyclic Repeater Jamming, ISCRJ),其中,ISDRJ和ISPRJ的主要区别在于切片的转发次数。ISCRJ与前两种干扰的不同之处在于切片的转发方式,ISCRJ转发当前和过去存储的切片。

    设干扰机的采样脉冲串p(t)是一串包络为矩形的脉冲串,其定义如下[4]

    p(t)=rect(tτdτj/2τj)L1i=0δ(tiTs) (2)

    其中,τd是干扰机采样时延,τj是采样脉冲的脉宽,Ts是采样脉冲的重复周期,L是采样脉冲的个数,代表卷积。

    τd=0时,干扰机同步采样,即干扰机采样时刻与雷达信号的起始时刻相同。当τd>0时,干扰机非同步采样,即干扰机采样时刻滞后于雷达信号的起始时刻。

    干扰机对雷达发射信号进行间歇采样相当于雷达发射信号与采样脉冲串的乘积,因此间歇采样后的信号可以表示为

    ss(t)=s(t)p(t) (3)

    进而干扰机通过控制间歇采样后的信号的转发次数和转发方式,可得到3种典型的ISRJ如下[27]

    jISDRJ(t)=Ajss(tτj) (4)
    jISPRJ(t)=AjMm=1ss(tmτj) (5)
    jISCRJ(t)=AjR1r=0ss[tτjr(τj+Ts)] (6)

    其中,jISDRJ(t), jISPRJ(t)jISCRJ(t)分别表示ISDRJ, ISPRJ和ISCRJ。Aj表示干扰幅度,M为ISPRJ的重复转发次数,且M=Ts/τj1R为ISCRJ的循环转发次数,且R=min(M,L),L=Tp/Tsmin()表示取最小值,表示向下取整。

    假设雷达观测场景中有一点目标,与雷达之间的距离为R0,对应的时延为t0=2R0/c,则目标回波信号为

    star(t)=σs(tt0) (7)

    其中,σ表示目标幅度。

    将式(3)代入式(4)—式(6)中,可以分别得到ISDRJ, ISPRJ和ISCRJ的具体表达式。为便于分析,将不同的ISRJ统一用jISRJ(t)表示,假设干扰机的起始距离与目标距离R0相等,则雷达接收到的干扰机产生的ISRJ可以表示为jISRJ(tt0)

    雷达接收的回波信号可以表示为

    sr(t)=star(t)+jISRJ(tt0)+n(t) (8)

    其中,n(t)表示高斯白噪声。

    对于ISRJ,本文采用主被动联合抗干扰策略。下面从变换域可分性和子脉冲相关性两个方面对本文脉内频率编码联合调频斜率捷变波形进行特性分析。

    基于脉内频率编码信号的抗ISRJ方法通常利用时频变换识别受干扰切片中的目标和干扰。STFT是一种广泛使用的时频变换,通过时间上对信号分段处理和傅里叶变换得到信号的时频分布。我们知道LFM信号的频谱近似为矩形,因此脉内频率编码信号的STFT在频率维上有较大的展宽。这导致了在干扰采样和转发次数较大情况下,部分切片中的干扰与目标难以在时频域分离。而LFM信号在最优阶次下的FrFT输出近似为匹配滤波,输出包络近似为sinc函数[28],表现出良好的能量聚焦性。利用这一性质,本文设计的脉内频率编码联合调频斜率捷变波形将回波切片变换到分数阶域,使受干扰切片中的目标和干扰在不同最优阶次下的分数阶域分别能量聚焦,提升了干扰和目标的可分性。下面通过仿真实验进一步验证:

    图3图4分别仿真了脉内频率编码信号受到ISRJ干扰时的时频分布和脉内调频斜率捷变信号受到ISRJ时在分数阶域的最优变换阶次分布。

    图  3  脉内频率编码回波信号时频分布
    Figure  3.  Time-frequency distribution of the intrapulse frequency coded echo signal
    图  4  脉内调频斜率捷变信号的最优变换阶次分布
    Figure  4.  Optimal fractional order distribution of the intrapulse FM slope agile signal

    仿真参数为:Tp=100 μs, N=20,脉内频率编码信号带宽5 MHz,脉内调频斜率捷变信号带宽范围2~10 MHz,采样频率20 MHz, τd=0, τj=Tsub, JSR=10 dB。

    图3中,由于干扰能量远大于目标,因此干扰信号在图3(b)图3(c)图3(d)中表现为黄色的“亮带”,而目标信号表现为“暗带”。脉内频率编码信号通过子脉冲中心频率跳变的方式,提高了目标信号与干扰信号在频率维的区分度,但是当干扰机间歇采样次数和重复转发次数较多时,如图中红圈标注所示,存在部分子脉冲与干扰切片在频率维高度混叠的现象。究其原因,是因为STFT对LFM信号的能量聚焦性较差,干扰信号和目标信号经过变换后在频率维均有较大的展宽,造成目标和干扰难以分离。相较之下,图4(b)图4(c)图4(d)中,脉内调频斜率捷变信号利用受干扰子脉冲中的目标和干扰的调频斜率差异,以及FrFT对LFM信号的能量聚焦性[11,28],使得干扰与目标在分数阶域的不同变换阶次下分别能量聚焦(较亮的点表示干扰,较暗的点表示目标),提高了目标与干扰的可分性。

    3.1节从变换域可分性的角度分析了脉内调频斜率捷变波形在分数阶域对目标与干扰进行分割的优势,然而单参数(即调频斜率)捷变会导致相邻子脉冲相关性较强,回波信号脉压后有较大的旁瓣,不利于目标的检测。为此,在脉内调频斜率捷变的基础上加上子脉冲中心频率捷变。下面结合发射波形自相关输出进行具体分析。

    仿真参数同上,图5对比了脉内调频斜率(单参数)捷变波形与脉内频率编码联合调频斜率(双参数)捷变波形的自相关输出结果。从图中可以看出,脉内调频斜率捷变波形由于相邻子脉冲的相关性较强,导致自相关结果中,主瓣前后有许多较高幅度的距离旁瓣。相较之下,脉内频率编码联合调频斜率捷变信号的旁瓣幅度明显降低。图6仿真了脉内频率编码联合调频斜率捷变波形分别在子脉冲数为4, 10和20时的自相关结果。从图中可知,旁瓣幅度随着子脉冲数的增大而降低。综上,双参数捷变能减小相邻子脉冲的相关性,并且增加子脉冲数能进一步减小发射波形自相关输出的旁瓣。

    图  5  单参数捷变波形与双参数捷变波形的自相关结果对比
    Figure  5.  Comparison of autocorrelation results between single-parameter agile waveform and dual-parameter agile waveform
    图  6  不同子脉冲数下的自相关结果
    Figure  6.  Autocorrelation results of different sub-pulse numbers

    主被动联合抗ISRJ方法中,通常在某一域对干扰进行抑制,例如时域[19,25]、频域[23]和分数阶域[24]。然而,单域信号处理往往难以兼顾干扰抑制后较小的目标能量损失和较少的干扰剩余。针对这一问题,本节在设计的发射波形基础上提出了分数阶域-时域联合抗ISRJ方法,旨在同时减小目标能量损失和降低由剩余干扰引起的距离旁瓣。首先通过FCM算法对回波切片进行分类。之后,对受干扰回波切片进行分数阶域带通滤波提取目标,减少目标能量损失。接着提取无干扰回波切片与发射信号匹配滤波并归一化处理构造时域滤波器,最后对分数阶滤波后回波信号的匹配滤波输出做进一步时域滤波,保护目标所在距离单元信号能量,并降低其他距离单元信号幅度。图7给出了ISRJ抑制方法的流程图。

    图  7  ISRJ抑制方法流程图
    Figure  7.  Flow chart of the proposed ISRJ suppression method
    4.1.1   FCM算法

    FCM[23,29,30]算法是一种软聚类方法,被广泛地应用于数据挖掘、机器学习和图像处理等领域。FCM算法能对样本进行模糊分类并将其划分到多个类别,但各自的隶属度不同。隶属度是指样本属于某类别的概率。隶属度取值在[0,1]之间,每个样本都有属于每个类别的隶属度,并且每个样本的所有隶属度之和为1。隶属度越高表明样本属于某类别的可能性越高。

    假设待分类的样本集{\boldsymbol{X}} = \left\{ {{x_1}{\text{, }}{x_2}{\text{, }} \cdots ,{\text{ }}{x_Z}} \right\}Z为待分类的样本个数,{x_i}为第i个样本,{x_{ij}}为样本{x_i}的第j个特征,每个样本有D个特征。FCM算法将样本集划分成K类,{C_j}\left( {j = 1,2, \cdots ,K} \right)为第j类的聚类中心。{u_{ij}}为样本{x_i}属于第j类的隶属度。

    FCM算法的目标函数和约束条件如下:

    {\boldsymbol{\eta}} ={\displaystyle \sum _{j}^{K}{\displaystyle \sum _{i}^{Z}{u}_{ij}^{\ell }}}{d}_{ij}^{2} (9)
    \sum\limits_{j = 1}^K {{u_{ij}}} = 1,{\text{ }}i = 1,2, \cdots ,Z (10)

    其中,\ell 为模糊加权指数,通常取值为2。{d_{ij}}为样本{X_i}到第j类的聚类中心{C_j}的欧几里得距离。

    FCM算法即是求目标函数{\boldsymbol{\eta}} 在约束条件下的最小值。这一问题可通过拉格朗日数乘法求解并最终推导出{u_{ij}}{C_j}应满足[29]

    {u_{ij}} = \frac{1}{{\displaystyle\sum\limits_{j' = 1}^K {{{\left( {\frac{{{d_{ij}}}}{{{d_{ij'}}}}} \right)}^{\textstyle \frac{2}{{\ell - 1}}}}} }} (11)
    {C_j} = \frac{{\displaystyle\sum\limits_{i = 1}^Z {u_{ij}^\ell {X_i}} }}{{\displaystyle\sum\limits_{i = 1}^Z {u_{ij}^\ell } }}\qquad\;\; (12)

    通过对式(11)和式(12)的迭代计算来更新隶属度{u_{ij}}和聚类中心{C_j}。当达到最大迭代次数或者前后两次迭代得到的隶属度变化小于阈值\varepsilon 时,停止迭代。此时可以根据隶属度{u_{ij}}对每个样本进行分类,{X_i}被划分的类别为

    {\mathrm{Label}}\left( {{x_i}} \right) = \mathop {\arg \max }\limits_{j \in \left\{ {1,2, \cdots ,K} \right\}} \left( {{u_{ij}}} \right) (13)
    4.1.2   回波切片分类

    干扰机的间歇采样、转发导致目标信号只有部分子脉冲受到ISRJ干扰。因此,若将回波信号划分为多个切片,回波切片可分为受干扰的回波切片和未受干扰的回波切片。此外,干扰信号为了达到欺骗和压制的效果,干扰信号能量通常远大于目标信号能量,这导致受干扰回波切片的能量和包络起伏程度要大于未受干扰回波切片的能量和包络起伏程度。考虑到均值和方差可以较好地衡量信号的能量和包络起伏特征,因此本文基于回波切片时域的均值和方差特征,利用FCM算法对回波切片进行分类。详细步骤如下:

    步骤1 将回波信号以发射信号子脉冲宽度{T_{{\mathrm{sub}}}}划分为N个切片。然后计算每个回波切片的均值和方差完成特征提取,D = 2{x_i} = \left( {\text{mean}}\left( {{s_{{\mathrm{r}}\_i}}} \right), {{\mathrm{var}}} \left( {{s_{{\mathrm{r}}\_i}}} \right) \right),{\text{ }}i = 1,2, \cdots ,N{s_{{\mathrm{r}}\_i}}表示第i个回波切片,{\text{mean}}\left( \cdot \right){{\mathrm{var}}} \left( \cdot \right)分别表示求均值和方差操作。回波切片的特征向量构成的样本集为X

    步骤2 依据回波切片的信号成分,设置回波切片的类别数为K。从\left[ {0,1} \right]中随机取值构成N \times K维隶属度矩阵U,并对每一行进行归一化处理。U的第i行第j列的元素为{u_{ij}}

    步骤3 利用FCM算法对回波切片进行分类,算法的伪代码如算法1所示。

    1  模糊C均值算法
    1.  Fuzzy C-mean algorithm
     输入:X:样本集;N:样本数;K:类别数;
     输出: {\mathrm{Label}}\left( {{x_i}} \right),i = 1,2, \cdots ,N :样本所属类别
     初始化:U:隶属度矩阵;
       \ell = 2:模糊加权指数;
       \varepsilon = {10^{ - 5}}:阈值;
       h = 1:迭代次数
     1: repeat
     2: 根据式(12)计算各聚类中心
     3: 根据式(11)更新隶属度矩阵U
     4: until \left\| {{{\boldsymbol{U}}_h} - {{\boldsymbol{U}}_{h - 1}}} \right\| \le \varepsilon
     5: else
     6: h = h + 1
     7: return {\mathrm{Label}}\left( {{x_i}} \right),i = 1,2, \cdots ,N
    下载: 导出CSV 
    | 显示表格

    步骤4 计算每个类别的聚类中心到原点的距离(2-范数)。由于未受干扰回波切片的能量较低,包络起伏程度较小,因此未受干扰回波切片对应类别的聚类中心的2-范数较小,根据每个类别的聚类中心的2-范数的大小关系可以确定各类别是属于受干扰回波切片还是未受干扰回波切片。最后从FCM算法返回值中可知每个切片被划分到的类别,进而判断回波切片是受干扰回波切片还是未受干扰回波切片。

    FrFT是一种广义的傅里叶变换,通过旋转时频平面将信号变换到介于时域与频域的中间域,在处理LFM信号上有着显著的优势[31]

    假设LFM信号x\left( t \right) = {\text{rect}}\left( {t/T} \right)\exp \left( {{\text{j}}\pi k{t^2}} \right)T为脉宽,k为调频斜率,x\left( t \right)的分数阶傅里叶变换为

    {X_\alpha }\left( u \right) = \int\limits_{ - \infty }^\infty {{K_\alpha }} \left( {t,u} \right)x\left( t \right){\text{d}}t (14)

    其中,\alpha 为变换角。 {K_\alpha }\left( {t,u} \right) 为变换核函数,其表达式[31]

    {K}_{\alpha }\left(t,u\right)=\left\{\begin{array}{ll}\sqrt{1-\text{j}{\cot}\,\alpha }\cdot \mathrm{exp}\left\{\text{j2}\pi \left(\dfrac{{t}^{2}+{u}^{2}}{2}\right){\cot}\,\alpha -\text{j}2\pi ut{\csc}\,\alpha \right\}, & \alpha \ne n\pi, \text{ }n为整数 \\ \delta \left(t-u\right), & \alpha =\text{2}n\pi \\ \delta \left(t+u\right), & \alpha +\pi =\text{2}n\pi \end{array} \right. (15)

    k + \cot \alpha = 0时,\alpha 被称作最优变换角{\alpha _{{\mathrm{opt}}}}{\alpha _{{\mathrm{opt}}}}满足[28]

    {\alpha _{{\mathrm{opt}}}} = - {\text{arc}} \cot \left( k \right) (16)

    {\alpha _{{\mathrm{opt}}}}下,式(14)可以进一步写为

    \begin{split} {X_{{\alpha _{{\mathrm{opt}}}}}}\left( u \right) = \,& \int\limits_{ - \infty }^\infty {{K_\alpha }} \left( {t,u} \right)x\left( t \right){\text{d}}t = \sqrt {1 - {\text{j}}\cot {\alpha _{{\mathrm{opt}}}}} \\ & \cdot {{\text{e}}^{{\text{j}}\pi {u^2}\cot {\alpha _{{\mathrm{opt}}}}}} \cdot T\frac{{\sin \left[ {\pi \left( {u\csc {\alpha _{{\mathrm{opt}}}}} \right)T} \right]}}{{\pi \left( {u\csc {\alpha _{{\mathrm{opt}}}}} \right)T}}\\[-1pt] \end{split} (17)

    通常k \gg 1,式(16)可得到近似关系[28]\csc {\alpha _{{\mathrm{opt}}}} \approx k, \left| {\sqrt {1 - {\text{j}}\cot {\alpha _{{\mathrm{opt}}}}} } \right| \approx \sqrt k 。将其代入式(17)得到

    \left| {{X_{{\alpha _{{\mathrm{opt}}}}}}\left( u \right)} \right| \approx \sqrt {BT} \cdot \left| {\sin {\text{c}}\left( {\pi Bu} \right)} \right| (18)

    由式(18)可知,FrFT对LFM信号近似为一个匹配滤波器。在 {\alpha _{{\mathrm{opt}}}} 下,LFM信号经过FrFT后得到的波形包络近似为sinc函数。

    依据FrFT对LFM信号良好的能量聚焦性,可在分数阶域提取目标信号的同时抑制绝大部分干扰,之后可构造时域滤波器对剩余干扰进一步抑制。分数阶域-时域联合滤波的具体步骤如下:

    步骤1 依据先验信息,将各回波切片变换到对应发射信号子脉冲的最优变换角下的分数阶域。

    假设第i个回波切片{s_{{\mathrm{r}}\_i}}\left( t \right)为受干扰回波切片,{s_{{\mathrm{r}}\_i}}\left( t \right)的表达式为

    {s_{{\mathrm{r}}\_i}}\left( t \right) = \sigma {s_i}\left( {t - {t_0}} \right) + {j_{{\mathrm{ISRJ}}\_i}}\left( {t - {t_0}} \right) + {n_i}\left( t \right) (19)

    其中,\sigma {s_i}\left( {t - {t_0}} \right), {j_{{\mathrm{ISRJ}}\_i}}\left( {t - {t_0}} \right){n_i}\left( t \right)分别表示第i个切片中的目标、ISRJ和噪声。

    {\alpha _i}(i = 1,2,···,N)为第i个发射子脉冲的最优变换角。第i个回波切片{s_{{\mathrm{r}}\_i}}\left( t \right){\alpha _i}下的FrFT为

    \begin{split} {F_{{\alpha _i}}}\left[ {{s_{{\mathrm{r}}\_i}}\left( t \right)} \right] \approx \,& \sigma \sqrt {{B_i}{T_{{\mathrm{sub}}}}} \cdot {{\mathrm{e}}^{{\mathrm{j}}{\varphi _i}}} \cdot \sin {\text{c}}\left\{ \pi {B_i}\left[ u -\right.\right. \\ & {u_0} - {f_i}\sin {\alpha _i} - \left( {2i - 1} \right)\\ & \cdot\left.\left.{T_{{\mathrm{sub}}}}\cos {\alpha _i}/2 \right] \right\} + {X_{{\mathrm{j}}\_i,{\alpha _i}}}\left( {u - {u_0}} \right) \\ & + {N_{i,{\alpha _i}}}\left( u \right) \\[-1pt] \end{split} (20)

    其中,{\varphi _i}为第i个切片经过FrFT后的目标相位。{X_{{\mathrm{j}}\_i,{\alpha _i}}}\left( {u - {u_0}} \right){N_{i,{\alpha _i}}}\left( u \right)分别表示第i个切片中经过FrFT后的干扰和噪声。{u_0} = {t_0}\cos {\alpha _i}

    未受干扰回波切片的处理同上,这里不再赘述。由式(20)可知,第i个切片中的目标在{\alpha _i}下的分数阶域形成峰值,而干扰信号与目标信号的调频斜率不同导致干扰信号无法积累。目标的峰值位置为{u_{{\mathrm{peak}}}} = {f_i}\sin {\alpha _i} + \left[ {{t_0} + \left( {2i - 1} \right){T_{{\mathrm{sub}}}}/2} \right]\cos {\alpha _i},据此可构造带通滤波器。

    步骤2 带通滤波器构造。目标的包络为sinc函数波形,其能量主要集中在主瓣之内。因此以其第一零点到原点的距离1/{B_i}为半窗长构造的第i个回波切片的分数阶域带通滤波器如下:

    {{\mathrm{Win}}_i}\left( u \right) = \left\{ {\begin{array}{*{20}{c}} {0,}&{\left| {u - {u_{{\mathrm{peak}}}}} \right| > 1/{B_i}} \\ {1,}&{\left| {u - {u_{{\mathrm{peak}}}}} \right| \le 1/{B_i}} \end{array}} \right. (21)

    步骤3 将FrFT后的回波切片通过分数阶域带通滤波提取目标的主瓣。然后,将提取的信号通过逆FrFT得到回波切片经过分数阶域滤波后的时域表示。所有回波切片经过相同的处理后,得到分数阶域滤波后的回波信号{\tilde s_{\mathrm{r}}}\left( t \right)

    步骤4 时域滤波器构造。FCM算法对回波切片进行分类后,将未受干扰的回波切片提取出来并记为{s_{{\mathrm{free}}}}\left( t \right)。将{s_{{\mathrm{free}}}}\left( t \right)与发射信号卷积,由此构造出归一化的时域滤波器[25]

    H\left( t \right) = \frac{{{s_{{\mathrm{free}}}}\left( t \right) \otimes {s^ * }\left( { - t} \right)}}{{\max \left[ {{s_{{\mathrm{free}}}}\left( t \right) \otimes {s^ * }\left( { - t} \right)} \right]}} (22)

    步骤5 将{\tilde s_{\mathrm{r}}}\left( t \right)与发射信号进行匹配滤波得到{\tilde s_{\mathrm{r}}}\left( t \right)的脉冲压缩输出为

    {\tilde y_{\mathrm{r}}}\left( t \right) = {\tilde s_{\mathrm{r}}}\left( t \right) \otimes {s^ * }\left( { - t} \right) = {\tilde y_{{\mathrm{tar}}}}\left( t \right) + {\tilde y_{{\mathrm{ISRJ}}}}\left( t \right) + {\tilde y_{\mathrm{n}}}\left( t \right) (23)

    其中, {\tilde y_{{\mathrm{tar}}}}\left( t \right) , {\tilde y_{{\mathrm{ISRJ}}}}\left( t \right) {\tilde y_{\mathrm{n}}}\left( t \right) 分别表示经过滤波和脉压后的目标、剩余干扰和噪声。

    步骤6 时域滤波。步骤4中构造的时域滤波器实质是未受干扰回波切片与雷达发射信号匹配滤波后,幅度归一化的结果。因为未受干扰回波切片中包含目标信息,所以其匹配滤波输出为峰值在目标距离单元处的窄主瓣波形。因此,对 {\tilde y_{\mathrm{r}}}\left( t \right) 进一步时域滤波可理解为对 {\tilde y_{\mathrm{r}}}\left( t \right) 进行时域加权,保留 {\tilde y_{\mathrm{r}}}\left( t \right) 在目标距离单元处的能量,并抑制其他距离单元的幅度,降低剩余干扰引起的距离旁瓣。

    {\tilde y_{\mathrm{r}}}\left( t \right) 进一步通过 H\left( t \right) 得到最终时域滤波后的输出 {y'_{\mathrm{r}}}\left( t \right)

    {y'_{\mathrm{r}}}\left( t \right) = {\tilde y_{\mathrm{r}}}\left( t \right) \cdot \left| {H\left( t \right)} \right| (24)

    为验证本文所提基于脉内频率编码联合调频斜率捷变波形的抗干扰方法的有效性,本节仿真实验主要从两个方面进行。一是干扰抑制效果仿真:针对3种典型的ISRJ进行干扰抑制以及抑制后脉压结果对比。二是方法性能分析,包括SNR和JSR分别对FCM算法分类的准确率、干扰抑制后目标检测概率的影响,以及干扰机采样时延和采样脉冲宽度对干扰抑制的影响。仿真中雷达发射波形和干扰相关参数如表1所示。

    表  1  仿真实验的参数设置
    Table  1.  Parameter settings for simulation experiments
    参数 数值
    雷达发射信号脉宽{T_{\rm p}} 100 μs
    子脉冲脉宽{T_{{\mathrm{sub}}}} 5 μs
    子脉冲数N 20
    子脉冲最小带宽{B_{\min }} 2 MHz
    子脉冲最大带宽{B_{\max }} 10 MHz
    子脉冲载频最大间隔{B_0} 10 MHz
    采样频率{f_{\mathrm{s}}} 20 MHz
    干扰机采样脉冲宽度{\tau _{\mathrm{j}}} 5 μs
    间歇采样重复周期{T_{\mathrm{s}}} 25 μs
    干信比(JSR) 20 dB
    信噪比(SNR) 0 dB
    下载: 导出CSV 
    | 显示表格

    为确保仿真参数的一致性,固定带宽编码序列a=[a(1),a(2),···,a(N)]=[8,10,18,15,4,3,20,12,17,11,13,6,19,9,5,7,14,2,16,1],载频编码序列b=[b(1),b(2),···,b(N)]=[6,15,5,12,7,13,4,18,3,16,2,19,9,11,8,20,10,17,1,14]。雷达发射的脉内频率编码联合调频斜率捷变信号的时频图如图8所示。

    图  8  脉内频率编码联合调频斜率捷变信号的时频图
    Figure  8.  Time-frequency diagram of the intrapulse frequency coded joint FM slope agile signal

    图8中,雷达发射信号由多个子脉冲构成,各个子脉冲均为LFM信号。各子脉冲的载频和调频斜率不一致。此外,相邻子脉冲还采取正负调频斜率调制,进一步提高了子脉冲间的相互掩护能力。

    假设干扰机对雷达发射信号进行同步采样,即{\tau _{\mathrm{d}}} = 0图9图10图11分别仿真了不同方法对单个ISDRJ, ISPRJ和ISCRJ的抗干扰效果。

    图  9  ISDRJ抑制
    Figure  9.  Interference suppression for ISDRJ
    图  10  ISPRJ抑制
    Figure  10.  Interference suppression for ISPRJ
    图  11  ISCRJ抑制
    Figure  11.  Interference suppression for ISCRJ

    图9(a)图10(a)图11(a)分别仿真了回波信号受到单个ISDRJ, ISPRJ和ISCRJ干扰时的时域波形,3种场景下的干扰占空比依次为0.25,0.80和0.50,目标在采样点2000处。由图9(b)图10(b)图11(b)可知,3种典型ISRJ经过脉压后均能在目标前后生成大量假目标,其中ISPRJ经过脉压后产生的假目标数目最多,ISDRJ产生的假目标数目最少。对比干扰抑制前后的脉压输出结果可知,本文所提方法能有效抑制ISRJ,脉压后的干扰幅度被极大程度的削弱,同时能检测到目标信号。图9(c)图10(c)图11(c)分别仿真了不同方法抗ISDRJ, ISPRJ和ISCRJ后的脉压输出的对比结果。其中,文献[10]采用分段脉压的方式对受干扰回波切片进行识别和剔除,剔除过程会不可避免造成目标能量较大损失,且随着干扰占空比的增大,脉压后目标峰值能量损失也增大。文献[11]采用分数阶域带通滤波的方式,相较文献[10]减少了目标能量损失。然而滤波过程中,滤波器通带内的干扰信号也被保留下来。同时随着干扰占空比增大,剩余干扰的能量也增大,从而在目标周围形成了较大幅度的旁瓣。文献[25]根据回波信号时频分布提取未受干扰回波切片构造时域窄带滤波器,随后直接对回波信号的脉压输出进行滤波处理。从结果看,落于时域窄带内的干扰增强了目标的峰值能量,但对于目标位置外的一些距离门内的强假目标,经过抑制后仍能形成高能量的尖峰,影响雷达检测。相较之下,本文采用分数阶域联合时域滤波的方法,利用分数阶域带通滤波减小目标能量损失并抑制绝大部分干扰,然后对滤波后的脉压输出进行时域滤波进一步减小了干扰剩余引起的距离旁瓣,从而达到较好的干扰抑制和检测效果。

    5.2.1   FCM算法的分类准确率

    回波切片的正确分类是实现干扰抑制的关键步骤。为评价FCM算法的分类性能,定义分类准确率(Classification Accuracy, CA)如下

    {\mathrm{CA}} = \frac{{100}}{{{\mathrm{MCN}}}}\sum\limits_{i = 1}^{{\mathrm{MCN}}} {\sum\limits_{j = 1}^{{N_{{\mathrm{total}}}}} {\frac{{{E_{ij}}}}{{{N_{{\mathrm{total}}}}}}} } (25)

    其中,{\mathrm{MCN}}表示蒙特卡罗实验的次数,{N_{{\mathrm{total}}}}表示回波中受干扰回波切片的总数,{E_{ij}}表示第i次蒙特卡罗实验中第j个受干扰回波切片的分类结果。当第j个受干扰回波切片被正确分类时,{E_{ij}}的值为1,否则为0。

    设回波信号脉压前SNR为–10~10 dB,JSR分别为5 dB, 10 dB, 15 dB和20 dB。蒙特卡罗实验次数设为200。文献[24]方法利用ISRJ的时域不连续性和能量压制性,在时域对受干扰回波切片进行识别。本文采用的FCM方法与文献[24]方法对受干扰回波切片的分类准确率和分类数目的对比结果如图12所示。

    图  12  FCM方法和文献[24]方法在不同SNR和JSR下对受干扰回波切片的分类准确率和分类数目
    Figure  12.  CA and identification number of the FCM method and Ref. [24] method for the interfered echo slice at different JSRs and SNRs

    图12可以看出,当SNR和JSR较小时,例如{\text{SNR}} = - 10{\text{ dB}}, {\text{JSR}} = 5{\text{ dB}}, FCM方法和文献[24]方法对受干扰回波切片的分类准确率均低于100%,这是因为SNR和JSR较低时,回波切片的主要成分是噪声,噪声的能量较强导致受干扰回波切片中干扰的特征不明显,影响了分类的效果。对照分类数目的曲线可知,低SNR和JSR下,FCM方法对受ISDRJ干扰的回波切片的分类数目大于实际数目,同时分类准确率小于100%,这说明既存在受干扰回波切片未被正确分类的情况也存在未受干扰回波切片被错误分类的情况,此外FCM方法和文献[24]方法对受ISPRJ和ISCRJ干扰回波切片的分类数目均小于实际数目,表明均有部分受干扰回波切片未被识别。随着SNR和JSR的增大,受干扰回波切片中的干扰能量逐渐增强,噪声的影响逐渐降低,使得受干扰回波切片中干扰特征逐渐明显,两种方法对回波中受干扰回波切片的分类性能也随之提升。对比两种方法的CA曲线和分类数目曲线可知,在相同JSR下,FCM方法的分类准确率达到100%时对SNR的要求更低。因此,FCM方法相较于文献[24]方法在低SNR下更加稳健。

    5.2.2   干扰抑制后目标检测概率

    抑制干扰的最终目的是检测目标。因此,本文将目标检测概率(Target Detection Probability, TDP)[32]作为干扰抑制性能的评价标准,TDP定义如下

    {\mathrm{TDP}} = \frac{{100}}{{{\mathrm{MCN}}}}\sum\limits_{i = 1}^{{\mathrm{MCN}}} {{E_i}} (26)

    其中,{E_i}表示第i次蒙特卡罗实验中对目标的检测结果。当目标被正确检测到时,{E_i}等于1,否则为0。

    本文所采用的分数阶域-时域联合滤波方法旨在通过多域联合滤波减小目标能量损失并降低剩余干扰带来的距离旁瓣。因此,性能分析方面与单域信号处理进行对比。设置JSR为0~50 dB,SNR分别为–10 dB, –5 dB和0 dB,蒙特卡罗实验次数为500。在本文设计的发射信号基础上,分别采用分数阶域滤波方法[11]、时域滤波方法[25]以及本文所提的分数阶域-时域联合滤波方法对3种典型干扰进行干扰抑制,抑制后的目标检测概率曲线如图13图15所示。

    图  13  不同方法抗ISDRJ后的目标检测概率
    Figure  13.  The TDP of echo signal containing ISDRJ after interference suppression by different methods
    图  14  不同方法抗ISPRJ后的目标检测概率
    Figure  14.  The TDP of echo signal containing ISPRJ after interference suppression by different methods
    图  15  不同方法抗ISCRJ后的目标检测概率
    Figure  15.  The TDP of echo signal containing ISCRJ after interference suppression by different methods

    图13图15可知,在较低的JSR下,例如JSR为0~20 dB,3种方法对干扰抑制后的目标检测概率均能达到100%,说明回波信号经过干扰抑制后的脉压输出能在目标处形成明显的峰值,通过峰值搜索能准确地检测到真实目标的位置。但当JSR较高时,3种方法抗干扰后对目标的检测性能均出现不同程度的降低。将TDP≥90%时,不同方法所能允许的最大JSR值作为指标来衡量不同方法的抗干扰能力,不同方法在不同SNR下的JSR容限如表2所示。

    表  2  不同方法在不同SNR下的JSR容限(dB)
    Table  2.  The JSR tolerance of different methods at different SNRs (dB)
    方法 –10 dB –5 dB 0 dB
    ISDRJ ISPRJ ISCRJ ISDRJ ISPRJ ISCRJ ISDRJ ISPRJ ISCRJ
    文献[11]方法 40 32.5 37.5 40 32.5 37.5 40 32.5 37.5
    文献[25]方法 30 22.5 27.5 40 27.5 32.5 50 35.0 37.5
    所提方法 50 50.0 50.0 50 50.0 50.0 50 50.0 50.0
    下载: 导出CSV 
    | 显示表格

    表2可以看出,针对3种ISRJ,所提方法在相同SNR下的JSR容限要高于文献[11]方法和文献[25]方法,文献[11]和文献[25]方法在高JSR下的抗干扰能力有所欠缺。其中,文献[25]直接对受干扰回波的脉压输出进行距离维上的加权来降低干扰的幅度,尽管假目标的幅度被极大抑制,但一些距离门内的强假目标仍能形成尖峰。尤其是在高JSR的情况下,这些距离门内的假目标仍能在能量上对目标进行压制,从而影响目标检测。而文献[11]在分数阶域带通滤波的过程中也会保留一部分干扰信号,在高JSR条件下,这些剩余干扰会产生幅度较大的距离旁瓣影响目标检测。此外,文献[25]和文献[11]方法的JSR容限与干扰占空比呈负相关关系,这是因为干扰占空比越大,回波中受干扰回波切片的数目就越多,导致回波脉压后会产生能量更强的假目标以及滤波后会有更多的干扰剩余,进而降低了高JSR下雷达发现目标的能力。对比之下,本文所提方法对SNR和干扰占空比不敏感,在JSR=50 dB时,抗干扰后的目标检测概率>90%,具有更好的干扰抑制性能。

    5.2.3   采样时延对干扰抑制的影响

    上述仿真实验均是针对干扰机同步采样的情形。下面针对干扰机非同步采样的情况,即{\tau _{\mathrm{d}}} \ne 0,对所提抗干扰方法的性能进行仿真分析。由于干扰机延时采样,回波切片大致可分为下列4种情形:无干扰回波切片,被少部分干扰的回波切片,被大部分干扰的回波切片以及被完全干扰的回波切片。因此将FCM算法的待分类的类别数K设为4。假设{\tau _{\mathrm{d}}} = 1{\text{ }} {\text{μs}}, {\text{SNR}} = 0{\text{ dB}}, {\text{JSR}} = 20{\text{ dB}} ,所提方法对3种ISRJ抑制后的脉压结果如图16所示。图17进一步仿真了{\tau _{\mathrm{d}}}/{T_{{\mathrm{sub}}}}的取值在0~0.5,蒙特卡罗实验次数为200,其他参数不变情况下,所提方法抗ISRJ后的脉压输出的目标幅度随{\tau _{\mathrm{d}}}/{T_{{\mathrm{sub}}}}的变化曲线。

    图  16  干扰抑制后脉压结果({\tau _{\mathrm{d}}} = 1{\text{ }} {\text{μs}})
    Figure  16.  Pulse compression results after interference suppression ({\tau _{\mathrm{d}}} = 1{\text{ }} {\text{μs}})
    图  17  采样延时与干扰抑制后目标归一化幅度关系
    Figure  17.  Relationship curve between sampling delay and normalized amplitude of target after interference suppression

    图16中,在干扰机非同步采样的情况下,所提方法均能有效对抗3种ISRJ,较好保留目标能量的同时也降低了距离旁瓣的影响。图17中,对于非同步ISDRJ和非同步ISCRJ,回波信号中的干扰切片数目不多,FCM算法能够较好地完成分类,因此干扰抑制后的目标幅度较为平稳,所提方法对非同步ISDRJ和非同步ISCRJ有着稳定的抗干扰性能。相较之下,对于非同步ISPRJ,干扰抑制后的目标幅度有着明显的起伏。这是因为ISPRJ的干扰占空比较大,回波中干扰切片多,FCM算法难以对所有回波切片进行准确分类,导致根据未受干扰回波切片构造的时域滤波器的主瓣位置偏离目标所在距离单元。因此在某些干扰机采样时延下,滤波后的目标能量有相对较大的衰减,导致了ISPRJ曲线的较大起伏。总体而言,干扰抑制后的目标归一化幅度>0.6,所提方法能有效抑制非同步采样下的ISRJ,并较好地保留了目标能量。

    5.2.4   采样脉冲宽度对干扰抑制的影响

    脉内波形捷变抗ISRJ的思想同传统脉间波形捷变抗全脉冲转发干扰的思想一致,即雷达波形变化速度快于干扰机复制转发干扰的速度。若干扰机采样脉冲宽度小于发射信号子脉冲宽度,则会造成转发的干扰信号会有部分处在当前子脉冲内,由于参数高度相似而难以被抑制,并且在高JSR和干扰切片数多场景下,脉内波形捷变方法失效。因此,在实际对抗场景中,雷达需要预先进行环境感知,对干扰参数进行粗略估计,使发射波形捷变速度快于干扰机复制转发干扰速度。

    上述实验是在干扰机采样脉冲宽度等于雷达发射信号子脉冲宽度的特例下进行的。下面仿真分析干扰机采样脉冲宽度大于雷达发射信号子脉冲宽度的情况。当干扰机采样脉冲宽度{\tau _{\mathrm{j}}}不是子脉冲宽度{T_{{\mathrm{sub}}}}整数倍时,干扰机对子脉冲有的部分采样,有的全部采样,因此将FCM算法的类别数设为4。当{\tau _{\mathrm{j}}}/{T_{{\mathrm{sub}}}} = 1.6, {\text{SNR}} = 0{\text{ dB}}, {\text{JSR}} = 20{\text{ dB}} ,干扰机采样延时为0时,所提方法对3种ISRJ抑制后的脉压结果如图18所示。进一步设置{\tau _{\mathrm{j}}}/{T_{{\mathrm{sub}}}}的取值为1~2,其他参数不变,蒙特卡罗实验次数为200,本文方法抗ISRJ后脉压输出的目标幅度随{\tau _{\mathrm{j}}}/{T_{{\mathrm{sub}}}}的变化曲线如图19所示。

    图  18  干扰抑制后脉压结果({\tau _{\mathrm{j}}}/{T_{{\mathrm{sub}}}} = 1.6)
    Figure  18.  Pulse compression results after interference suppression ({\tau _{\mathrm{j}}}/{T_{{\mathrm{sub}}}} = 1.6)
    图  19  采样脉冲宽度与干扰抑制后目标归一化幅度关系
    Figure  19.  Relationship curve between sampling pulse duration and normalized amplitude of target after interference suppression

    图18可知,所提方法在干扰机采样脉冲宽度大于发射信号子脉冲宽度时,仍能有效抑制干扰。图19中,3种干扰抑制后的目标幅度曲线随采样脉冲宽度的变化均有较大起伏。主要有两点原因。一是抑制前一些采样脉冲宽度下的干扰形成的假目标落在目标所在距离单元,导致滤波后的剩余干扰提高了目标所在距离单元的幅度。二是由于采样脉冲宽度大于子脉冲宽度,回波切片分类困难,使得一些采样脉冲宽度下构造的滤波器的主瓣偏离目标所在距离单元,在时域滤波时,削弱了目标所在距离单元的幅度。综合两种因素,造成了3种干扰抑制后目标幅度曲线的较大起伏。

    在自卫式干扰背景下,本文针对间歇采样转发式干扰提出了一种基于脉内频率编码联合调频斜率捷变波形的抗干扰方法。该方法利用干扰机收发分时特点对回波切片进行分类,然后通过联合滤波可有效抑制干扰。基于方法过程和实验结果,该方法有如下优点:

    (1) 雷达发射波形采用脉内双参数捷变,脉压后旁瓣幅度低,子脉冲间的相互掩护能力强;

    (2) FCM算法在低SNR下对单个干扰机同步采样形成的ISRJ的回波切片分类效果较好。当{\text{SNR}} \ge - 2.5{\text{ dB}} {\text{JSR}} \ge 5{\text{ dB}} 时,对受干扰回波切片的分类准确率可达100%。此外,所提联合滤波方法兼顾了较小的目标能量损失和较低的距离旁瓣。同时在高JSR下,所提方法的目标检测概率≥90%。

    (3) 所提方法不仅能有效对抗干扰机同步采样场景下的3种典型的ISRJ,还能在干扰机非同步采样场景下对干扰进行有效抑制。

    然而需要说明的是,本文方法需要检测回波前沿来划分回波切片,因此所提方法适用于雷达近距离探测场景,在该场景下,目标信噪比本身比较大,然后可借助脉压前相参积累、慢时间特征量和雷达跟踪状态下滑窗检测等手段检测回波前沿。而当雷达远距离探测或目标回波能量本身很小时,则由于无法检测前沿导致方法失效。此外,在多个时间上错位的复合干扰或强噪声干扰的背景下,FCM算法无法对回波中切片进行准确分类,使得干扰抑制性能不佳。这些问题都将在未来的工作中进一步研究。

  • 图  1  平面阵机载雷达几何照射图

    Figure  1.  The geometry of planar phased-array airborne radar

    图  2  张量基处理流程

    Figure  2.  The flow of tensor-based processing

    图  3  两种不同结构对比

    Figure  3.  The comparison of two different structures

    图  4  三维杂波谱估计结果

    Figure  4.  The estimation of three-dimensional clutter spectrum

    图  5  不同SR-STAP方法的SCNR损失结果

    Figure  5.  The SCNR loss results of different SR-STAP methodss

    图  6  传统STAP方法的SCNR损失结果

    Figure  6.  The SCNR loss results of conventional STAP methods

    图  7  平均SCNR损失结果

    Figure  7.  The result of average SCNR loss

    图  8  不同数据大小下实际运行时间

    Figure  8.  The running time of different input data size

    表  1  TMSBL算法

    Table  1.   TMSBL algorithm

     输入:方位字典{{\boldsymbol{S}}_{\rm{a}}},俯仰字典{{\boldsymbol{S}}_{\rm{e}}},多普勒字典{{\boldsymbol{S}}_{\rm{d}}},训练样本集合{\boldsymbol{X}}
     输出:稀疏系数{\boldsymbol{\varXi}}
     初始化:过完备字典{\boldsymbol{S}} = {{\boldsymbol{S}}_{\rm{e}}} \otimes {{\boldsymbol{S}}_{\rm{a}}} \otimes {{\boldsymbol{S}}_{\rm{d}}},稀疏控制系数
         {{\boldsymbol{\gamma}} _0} = {{\boldsymbol{e}}_{ {N_{\rm{a} } }{M_{\rm{e} } }{K_{\rm{d} } } } },均值{Y_0} = 1,噪声方差\sigma _0^2,最大迭代
         次数{i_{\max} },收敛阈值\mu
     1:对于 i = 1:{i_{\max} }执行
     2:   对于j = 1:{N_{\rm{a}}}{M_{\rm{e}}}{K_{\rm{d}}}执行
     3:     计算方位索引{ {\rm{loc} }_{\rm{a}}},俯仰索引{{\rm{loc}}_{\rm{e}}}和多普勒索引{{\rm{loc}}_{\rm{d}}}
     4:     {\mathcal{T}_{:,:,:,j} } =
           {\mathcal{F}_{N,M,K} }\left\{ { {\mathcal{V}_{NK} }\left\{ { {\gamma _j}{\boldsymbol{S} }_{:,{ {\rm{loc} }_{\rm{d} } } }^{\rm{d} }{ {\left( { {\boldsymbol{S} }_{:,{ {\rm{loc} }_{\rm{a} } } }^{\rm{a} } } \right)}^{\rm{T} } } } \right\}{ {\left( { {\boldsymbol{S} }_{:,{ {\rm{loc} }_{\rm{e} } } }^{\rm{e} } } \right)}^{\rm{T} } } } \right\}
     5:   结束循环
     6:   {\boldsymbol{C} } = {\mathcal{M}_{NMK,NMK} }\left\{ {\mathcal{T}{ \times _2}{\boldsymbol{S} }_{\rm{d} }^{\rm{H} }{ \times _3}{\boldsymbol{S} }_{\rm{a} }^{\rm{H} }{ \times _4}{\boldsymbol{S} }_{\rm{e} }^{\rm{H} } } \right\} + {\sigma ^2}{\boldsymbol{I}}
     7:   {\boldsymbol{Y}} = {\mathcal{M}_{ {N_{\rm{a}}}{M_{\rm{e}}}{K_{\rm{d}}},NMK} }
           \cdot\left\{ {\mathcal{D} \odot \left( { {\mathcal{C}^{ - 1} }{ \times _1}{\boldsymbol{S} }_{\rm{d} }^{\rm{H} }{ \times _2}{\boldsymbol{S} }_{\rm{a} }^{\rm{H} }{ \times _3}{\boldsymbol{S} }_{\rm{e} }^{\rm{H} } } \right)} \right\}{\boldsymbol{X} }
     8:   对于j = 1:{N_{\rm{a}}}{M_{\rm{e}}}{K_{\rm{d}}}执行
     9:      {\varSigma _{j,j} } = {\gamma _j} - { {\boldsymbol{Q} }_{j,:} }{ {\boldsymbol{T} }_{:,j} }
     10:     {\gamma _{j + 1} } = \dfrac{1}{L}\left\| { { {\boldsymbol{Y} }_{j,:} } } \right\|_2^2 + {\varSigma _{j,j} }
     11:     {D_j} = \dfrac{ { {\varSigma _{j,j} } } }{ { {\gamma _j} } }
     12:   结束循环
     13:   {\sigma ^2} = \dfrac{1}{ {NMKL} }\left\| {{\boldsymbol{X}} - {\boldsymbol{SY}}} \right\|_{\rm{F} }^2

            + \dfrac{ {\sigma _{i - 1}^2} }{ {NMK} }\mathop \sum \limits_{j = 1}^{ {N_{\rm{a} } }{M_{\rm{e} } }{K_{\rm{d} } } } \left( {1 - {D_j} } \right)
     14:   如果 \left\| {{\boldsymbol{Y}} - {{\boldsymbol{Y}}_{i - 1} } } \right\|_{\rm{F} }^2/\left\| {\boldsymbol{Y}} \right\|_{\rm{F} }^2 \le \mu 跳出循环
     15:结束循环
    下载: 导出CSV

    表  2  计算复杂度

    Table  2.   Computational complexity

    方法复乘法次数计算复杂度
    OMP\left( {NMK{N_{\rm a}}{M_{\rm e}}{K_{\rm d}} + r_{\rm s}^3 + NMKr_{\rm s}^2 + 2NMK{r_{\rm s}}} \right)L{K_{\rm OMP}}O\left( {NMK{N_{\rm a}}{M_{\rm e}}{K_{\rm d}}L{K_{\rm OMP}}} \right)
    MIAA\left( {2{N_{\rm a}}{M_{\rm e}}{K_{\rm d}}{{\left( {NMK} \right)}^2} + {{\left( {NMK} \right)}^3} + \left( {L + 1} \right){N_{\rm a}}{M_{\rm e}}{K_{\rm d}}NMK + L{N_{\rm a}}{M_{\rm e}}{K_{\rm d}}} \right){K_{\rm MIAA}}O\left( {{N_{\rm a}}{M_{\rm e}}{K_{\rm d}}{{\left( {NMK} \right)}^2}{K_{\rm MIAA}}} \right)
    MFOCUSS\left( {NMK{N_{\rm a}}{M_{\rm e}}{K_{\rm d}}L + {{\left( {NMK} \right)}^3} + 2{{\left( {NMK} \right)}^2}{N_{\rm a}}{M_{\rm e}}{K_{\rm d}} + NMK{{\left( {{N_{\rm a}}{M_{\rm e}}{K_{\rm d}}} \right)}^2}} \right){K_{\rm MFOC}}O\left( {NMK{{\left( {{N_{\rm a}}{M_{\rm e}}{K_{\rm d}}} \right)}^2}{K_{\rm MFOC}}} \right)
    MSBL \begin{gathered} \left( {{{\left( {{N_{\rm a}}{M_{\rm e}}{K_{\rm d}}} \right)}^3} + 4NMK{{\left( {{N_{\rm a}}{M_{\rm e}}{K_{\rm d}}} \right)}^2} + \left( {3{{\left( {NMK} \right)}^2} + {L^2} + 2NMKL} \right){N_{\rm a}}{M_{\rm e}}{K_{\rm d}}} \right. \hfill \\ \left. { + 4{{\left( {NMK} \right)}^3}/3 + {{\left( {NMK} \right)}^2} + NMKL} \right){K_{\rm MSBL}} \hfill \\ \end{gathered} O\left( {{{\left( {{N_{\rm a}}{M_{\rm e}}{K_{\rm d}}} \right)}^3}{K_{\rm MSBL}}} \right)
    MFCSBL\left( { {N_{\rm a} }{M_{\rm e} }{K_{\rm d} }\left( {5{ {\left( {NMK} \right)}^2} + 2L{{NMK} } + 4{{NMK} } + 3} \right) + { {\left( {NMK} \right)}^3} + L{{NMK} } + L} \right){K_{ {\rm{MFCSBL} } } }O\left( { {N_{\rm a} }{M_{\rm e} }{K_{\rm d} }{ {\left( {NMK} \right)}^2}{K_{\rm MFCSBL} } } \right)
    TMSBL \begin{gathered} \left( {{N_{\rm a}}{M_{\rm e}}{K_{\rm d}}(3NMK + NM + N + {L^2} + 2NMKL + NM{K^2} + NK{M^2})} \right. + 2{(NMK)^3}/3 \hfill \\ \left. { + {{(NMK)}^2}(1 + {K_{\rm d}} + {M_{\rm e}}) + NMK({N_{\rm a}}{M_{\rm e}}KN + NM{N_{\rm a}}{K_{\rm d}} + L)} \right){K_{\rm TMSBL}} \hfill \\ \end{gathered} O\left( {{N_{\rm a}}{M_{\rm e}}{K_{\rm d}}NM{K^2}{K_{\rm TMSBL}}} \right)
    下载: 导出CSV

    表  3  雷达系统参数

    Table  3.   Radar system parameters

    参数符号数值
    载机速度{v_{\rm{p}}}150 m/s
    载机高度H8000 m
    行/列阵元数M/N6/8
    脉冲数K8
    雷达波长\lambda 0.1 m
    脉冲重复频率{f_{\rm{r}}}8100 Hz
    阵元间距d0.05 m
    主波束指向\theta /\varphi –90°/0°
    杂噪比60 dB
    下载: 导出CSV
  • [1] REED I S, MALLETT J D, and BRENNAN L E. Rapid convergence rate in adaptive arrays[J]. IEEE Transactions on Aerospace and Electronic Systems, 1974, AES-10(6): 853–863. doi: 10.1109/TAES.1974.307893
    [2] 谢文冲, 段克清, 王永良. 机载雷达空时自适应处理技术研究综述[J]. 雷达学报, 2017, 6(6): 575–586. doi: 10.12000/JR17073

    XIE Wenchong, DUAN Keqing, and WANG Yongliang. Space time adaptive processing technique for airborne radar: An overview of its development and prospects[J]. Journal of Radars, 2017, 6(6): 575–586. doi: 10.12000/JR17073
    [3] MENG Xiangdong, WANG Tong, WU Jianxin, et al. Short-range clutter suppression for airborne radar by utilizing prefiltering in elevation[J]. IEEE Geoscience and Remote Sensing Letters, 2009, 6(2): 268–272. doi: 10.1109/LGRS.2008.2012126
    [4] DUAN Keqing, XU Hong, YUAN Huadong, et al. Reduced-DOF three-dimensional STAP via subarray synthesis for nonsidelooking planar array airborne radar[J]. IEEE Transactions on Aerospace and Electronic Systems, 2020, 56(4): 3311–3325. doi: 10.1109/TAES.2019.2958174
    [5] SUN Ke, MENG Huadong, WANG Yongliang, et al. Direct data domain STAP using sparse representation of clutter spectrum[J]. Signal Processing, 2011, 91(9): 2222–2236. doi: 10.1016/j.sigpro.2011.04.006
    [6] YANG Zhaocheng, LI Xiang, WANG Hongqiang, et al. On clutter sparsity analysis in space-time adaptive processing airborne radar[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 10(5): 1214–1218. doi: 10.1109/LGRS.2012.2236639
    [7] DUAN Keqing, XU Hong, YUAN Huadong, et al. Three-dimensional sparse recovery space-time adaptive processing for airborne radar[J]. The Journal of Engineering, 2019, 2019(19): 5478–5482. doi: 10.1049/joe.2019.0343
    [8] GUO Yiduo, LIAO Guisheng, and FENG Weike. Sparse representation based algorithm for airborne radar in beam-space post-Doppler reduced-dimension space-time adaptive processing[J]. IEEE Access, 2017, 5: 5896–5903. doi: 10.1109/ACCESS.2017.2689325
    [9] HAN Sudan, FAN Chongyi, and HUANG Xiaotao. A novel STAP based on spectrum-aided reduced-dimension clutter sparse recovery[J]. IEEE Geoscience and Remote Sensing Letters, 2017, 14(2): 213–217. doi: 10.1109/LGRS.2016.2635104
    [10] WANG Zetao, XIE Wenchong, DUAN Keqing, et al. Clutter suppression algorithm based on fast converging sparse Bayesian learning for airborne radar[J]. Signal Processing, 2017, 130: 159–168. doi: 10.1016/j.sigpro.2016.06.023
    [11] QIU Wei, ZHOU Jianxiong, ZHAO Hongzhong, et al. Three-dimensional sparse turntable microwave imaging based on compressive sensing[J]. IEEE Geoscience and Remote Sensing Letters, 2015, 12(4): 826–830. doi: 10.1109/LGRS.2014.2363238
    [12] SIDIROPOULOS N D, DE LATHAUWER L, FU Xiao, et al. Tensor decomposition for signal processing and machine learning[J]. IEEE Transactions on Signal Processing, 2017, 65(13): 3551–3582. doi: 10.1109/TSP.2017.2690524
    [13] ZHAO Rongqiang, WANG Qiang, FU Jun, et al. Exploiting block-sparsity for hyperspectral Kronecker compressive sensing: A tensor-based Bayesian method[J]. IEEE Transactions on Image Processing, 2020, 29: 1654–1668. doi: 10.1109/TIP.2019.2944722
    [14] 姜磊, 王彤. 机载雷达自适应对角加载参数估计方法[J]. 电子与信息学报, 2016, 38(7): 1752–1757. doi: 10.11999/JEIT151003

    JIANG Lei and WANG Tong. An adaptive estimation method of diagonal loading parameter for airborne radar[J]. Journal of Electronics &Information Technology, 2016, 38(7): 1752–1757. doi: 10.11999/JEIT151003
    [15] MA Zeqiang, LIU Yimin, MENG Huadong, et al. Jointly sparse recovery of multiple snapshots in STAP[C]. 2013 IEEE Radar Conference, Ottawa, Canada, 2013: 1–4. doi: 10.1109/RADAR.2013.6586083
    [16] DUAN Keqing, WANG Zetao, XIE Wenchong, et al. Sparsity-based STAP algorithm with multiple measurement vectors via sparse Bayesian learning strategy for airborne radar[J]. IET Signal Processing, 2017, 11(5): 544–553. doi: 10.1049/iet-spr.2016.0183
    [17] TROPP J A and GILBERT A C. Signal recovery from random measurements via orthogonal matching pursuit[J]. IEEE Transactions on Information Theory, 2007, 53(12): 4655–4666. doi: 10.1109/TIT.2007.909108
    [18] YANG Zhaocheng, LI Xiang, WANG Hongqiang, et al. Adaptive clutter suppression based on iterative adaptive approach for airborne radar[J]. Signal Processing, 2013, 93(12): 3567–3577. doi: 10.1016/j.sigpro.2013.03.033
    [19] YANG Shuyuan, LI Bin, WANG Min, et al. Compressive direction-of-arrival estimation via regularized multiple measurement FOCUSS algorithm[C]. 2014 International Joint Conference on Neural Networks, New York, USA, 2014: 2800–2803. doi: 10.1109/IJCNN.2014.6889967.
    [20] BADER B W and KOLDA T G. Matlab tensor toolbox version 2.6[EB/OL]. http://www.sandia.gov/tgkolda/TensorToolbox. 2015.
    [21] GUERCI J R. Space-Time Adaptive Processing for Radar[M]. Boston: Artech House, 2014: 52–73.
    [22] HALE T B, TEMPLE M A, RAQUET J F, et al. Localized three-dimensional adaptive spatial-temporal processing for airborne radar[C]. 2002 International Radar Conference, Edinburgh, UK, 2002: 191–195. doi: 10.1049/cp:20020275.
    [23] 洪玺, 王文杰, 殷勤业. 基于多级维纳滤波器的空时自适应信号处理及其在无线通信系统中的应用[J]. 信号处理, 2017, 33(3): 430–436. doi: 10.16798/j.issn.1003-0530.2017.03.026

    HONG Xi, WANG Wenjie, and YIN Qinye. Multistage wiener filter based space and time adaptive signal processing and its application in wireless communication system[J]. Journal of Signal Processing, 2017, 33(3): 430–436. doi: 10.16798/j.issn.1003-0530.2017.03.026
    [24] 王璐, 吴仁彪. 直接数据域空时自适应单脉冲方法[J]. 系统工程与电子技术, 2016, 38(12): 2738–2744. doi: 10.3969/j.issn.1001-506X.2016.12.09

    WANG Lu and WU Renbiao. Direct data domain space-time adaptive monopulse method[J]. Systems Engineering and Electronics, 2016, 38(12): 2738–2744. doi: 10.3969/j.issn.1001-506X.2016.12.09
  • 加载中
图(8) / 表(3)
计量
  • 文章访问数: 1986
  • HTML全文浏览量: 983
  • PDF下载量: 169
  • 被引次数: 0
出版历程
  • 收稿日期:  2021-09-26
  • 修回日期:  2021-12-03
  • 网络出版日期:  2021-12-23
  • 刊出日期:  2021-12-28

目录

/

返回文章
返回