基于调制谱超分辨重构的微多普勒参数估计

赵庆媛, 赵志强, 叶春茂, 鲁耀兵

(北京无线电测量研究所, 北京 100854)

0 引 言

众多研究表明,喷气发动机调制(jet engine modulation, JEM)特征对于窄带雷达气动目标分类是有效的[1-3]。但常规雷达在实际应用中受雷达波形限制,如短驻留、脉冲重复周期(pulse repetition time, PRT)捷变、脉冲重复周期过长等,不能有效地提取稳健的JEM特征,进而影响识别准确率。因此,在雷达资源调度受限的场景下研究准确的JEM特征提取意义重大。

调制谱间隔,又称基频,是气动目标识别的核心特征,取决于微动部件本身的物理特性,为桨叶个数与转速之积,不受观测角度以及雷达载频等参数的影响。文献[2]用复信号自回归双谱法对调制谱间隔进行估计。文献[3]用自相关的方式实现微多普勒周期估计。文献[4-5]将雷达回波转化为时频图,并通过图像边缘提取的方法实现较高精度的微动参数估计。文献[6]在时频变换的基础上结合倒谱变换实现微动部件转动频率的估计。以上的调制谱周期估计方法,都要求相参积累PRT均匀、多普勒频率不折叠,且估计精度受噪声水平影响较大,很难自动进行特征提取。

近年来,空域谱估计算法[7-13]研究成果被进一步扩展到时频域谱估计,取得很好的效果[14-16]。文献[17]用多重信号分类(multiple signal classification, MUSIC)算法实现了雷达时分体制下调制谱的超分辨估计。文献[18]提出幅度相位估计 (amplitude and phase estimation, APES)实现非参数谱估计。文献[19]在APES的基础上提出缺损APES (gapped APES, GAPES)进行缺失数据情况下的谱估计。文献[20]提出用自适应迭代算法(iterative adaptive algorithm, IAA)进行频率估计,显著减少了谱泄露并提高了频谱分辨率。文献[21]将APES与IAA相结合改善了谱泄露问题。以IAA为基础,文献[22-23]在谱稀疏的假设条件下分别提出通过迭代最小化的稀疏学习(sparse learning via iterative minimization，SLIM)算法和基于稀疏迭代协方差的估计(sparse iterative covariance-based estimation, SPICE),进一步提高了谱分辨能力。SLIM算法[22]基于最小二乘算法进行谱估计,而SPICE算法[23-24]则基于加权协方差拟合准则进行谱估计,精度更高,抗噪能力更强,且对于谱的稀疏性要求更低[23],有助于获得更高的调制谱重建精度。

本文基于气动目标调制谱特性对SPICE算法进行了针对性设计,并根据基频形成原理获得基频谱,进而通过峰值搜索基频谱自动估计微多普勒参数。仿真和实测数据分析表明,本文方法在短驻留、PRT捷变及调制谱折叠场景下均可实现稳健精确的基频估计。

1 气动目标微动模型与分析

1.1 理想旋翼雷达回波模型

假设飞机在雷达远场,有一个微动部件,共有N个桨叶,桨叶大小一致,均匀分布。桨叶根部到旋翼中心的距离为L1,桨叶顶部到旋转中心的距离为L2,桨叶角为φ,转速为wr,初始相位角为θ0,桨叶相对于雷达平台方位角为α,俯仰角为β，雷达波长为λ。则飞机旋翼基带回波[25]表示如下:

(1)

式中：θn=θ0+2πn/N+2πwrt-α, 另外:

(2)

(3)

(4)

1.2 调制谱分析

根据第1.1节中的理想回波模型,频域可表示为

(5)

式中:Nf表示单边调制谱谱线个数；fT表示基频；cn表示谱线的复幅度；由λ,N,β,α,φ,L1,L2,θ0及贝塞尔函数共同决定。旋翼回波在频域的状态是离散的谱线,因此调制谱谱线在频域是稀疏的,满足稀疏谱估计的假设条件。

多普勒离散谱线以主体谱线为中心,以固定间隔fT向正负频率对称平移。基频fT为桨叶个数与转动频率之积:

fT=Nwr

(6)

单边调制谱谱线个数Nf与微多普勒单边扩展和基频有关:

(7)

(8)

式中:vtip为桨叶叶尖的转速,根据空气动力学相关理论可知,叶尖的转速一般不会超过声速。

此外,只有当fmax>fT的情况下才能观测到调制谱谱线。巡航状态喷气式飞机的桨叶转动很快,其Nf为零,因此低频雷达不能观测到其调制谱谱线。

1.3 调制谱折叠分析

以上分析的假设条件是脉冲重复频率大于微动调制双边展宽。在脉冲重复频率较高,微动多普勒频率存在折叠的情况下,式(5)可写为

(9)

式中：φ(nfT)表示调制谱谱线折叠到(-fr/2,fr/2]区间的谱线集合,谱线个数为2Nf+1。φ(nfT) 可表示为

(10)

式中：fr表示脉冲重复频率。

2 基于SPICE的调制谱间隔估计

2.1 算法流程

由第1.2节分析可知,旋翼回波在频域是离散谱线,调制谱是稀疏的,满足SPICE谱估计的假设条件。

根据调制谱的稀疏性和固定间隔特性,本文提出SPICE谱估计与基频组功率累积相结合的微多普勒参数估计方法。第1步是对雷达时域回波运用SPICE谱估计获得高分辨调制谱,第2步基于SPICE谱进行基频组功率累积,并结合先验基频范围,得到基频谱(fundamental frequency spectrum, FFS),第3步根据FFS特性进行基频估计,如图1所示。

图1 基频估计算法示意图Fig.1 Schematic diagram of fundamental frequency estimation algorithm

2.2 SPICE谱估计

雷达时域回波表示为y=[y1,y2,…,yn,…,yN]T,其中,yn是目标在第n个PRT目标回波脉冲压缩后的复幅度,yn=y(tn) (n=1,2,…,N),tn对应第n个PRT,N表示相参积累脉冲个数。

SPICE谱估计模型可表示如下:

(11)

(12)

式中：‖·‖2表示矩阵的Frobenius范数；R-1/2为R-1的厄尔米特正定平方根。

通过对上式进行简化和变换[23],最小化代价函数f可转换为下面的线性约束最小化问题:

(13)

SPICE算法通过迭代运算实现全局最优化,包括初始化和迭代步骤。

初始化:

(14)

重复下面迭代过程直至收敛:

(15)

式中：上标i表示迭代次数。本文设定收敛标准为‖Pi+1-Pi‖/‖Pi‖达到极小值。

值得注意的是,本文中的SPICE算法涉及超参数K值的选取,K值设置与谱分辨需求有关。谱分辨率为脉冲重复频率除以K值。高重波形工作条件下,为实现较高的基频估计精度,K值需相应增大。目前还有一些研究尝试用无网格设置的方式实现SPICE估计,但都面临抗噪性差、需要较多快拍数且无法适应非均匀阵列等问题[27-30]。

2.3 基频组功率累积及基频估计

根据式(6)～式(8)和式(10)可得到基频fT对应的基频组φ(nfT)。其中,vtip由空气动力学原理可设置为340 m/s,基频范围根据目标先验信息进行设置,仰角β可根据航迹信息得到。

将基频组对应的谱功率进行累加即为基频组功率累积。定义FFS为基频组功率累积占比。fT对应的FFS表示如下:

(16)

倒谱[31]方法也适用于求解基频,但其在调制谱存在折叠的情况下效果不理想。

值得注意的是FFS有多个峰值,其中峰值对应的最大基频为实际调制谱间隔,其他峰值对应微动基频的公约数。

为了进一步稳健提取调制谱间隔,对FFS进一步处理得到FFS组功率占比平均值FFS(fT)/(2Nf+1)。并设定FFS组功率占比和FFS组功率占比平均值均为峰值的基频对基频估计值。至此,实现了自动参数估计。

3 验证与分析

3.1 仿真验证与分析

为了验证本文算法估计基频的有效性,利用式(1)～式(4)进行雷达回波仿真,并分析驻留时间、PRT参差、调制谱折叠以及噪声水平对算法的影响。

雷达载频设定为250 MHz。俯仰角、方位角和桨叶角均设置为0°,初始相位随机设置。根据式(8)计算得到fmax。目标微动参数见表1,目标1为直升机Z9,目标2为直升机M171,目标3为螺旋桨飞机Y5,目标4为喷气式飞机。选取两类直升机作为仿真目标的目的有两个, 一是研究奇偶桨叶数对本文方法是否有影响;二是验证基频估计对于气动目标型号识别的有效性。

其中目标4的基频大于调制谱单边展宽fmax,因此无法观测到微动对应的调制谱线,在仿真以及实测数据中此结论得到验证。

本文结合目标和雷达相关参数将基频范围设定为10～200 Hz。

3.1.1 短驻留场景下的仿真

在本节中,我们通过仿真验证短驻留场景下本文方法对于气动目标基频估计的有效性,并分析相参积累时间对基频估计的影响。

设置PRT为1 ms,相参脉冲个数为70,相参积累时间为70 ms。加入复正态白噪声,均值为零,噪声标准方差为0.1。SPICE算法中的K设置为1 000,频域分辨率为1 Hz。

设定SPICE谱估计和IAA谱估计收敛标准均为‖Pi+1-Pi‖/‖Pi‖<10-4。

图2(a)～图2(d)表示目标1～目标4的谱估计,对比SPICE、IAA和基于泰勒窗的快速傅里叶变换(fast Fourier transform, FFT),可以看到SPICE算法对频率分辨率最高,IAA次之,FFT最差。对于目标1～目标3,SPICE算法都能够正确估计多个调制谱谱峰。

图2 70 ms驻留情况下的谱估计结果Fig.2 Spectrum estimation result of 70 ms dwell situation

目标4对应喷气式飞机,其SPICE谱估计结果是一根谱线,无法观测到调制谱间隔,这是由于调制谱间隔大于微动多普勒展宽。结合理论分析和喷气式飞机实测微动特性可以推测,即使在无遮挡状态下,低于1 GHz载频的低频雷达也很难观测到正常巡航状态下喷气式飞机调制谱谱峰。

图3(a)～图3(d)表示FFS,横坐标为基频,纵坐标为FFS组功率累积占比。从图3可以看到目标1的FFS最大值对应的基频除了24 Hz还有12 Hz;目标2的FFS最大值对应的基频为15 Hz;目标3的FFS具有4个最大值,对应的基频除了80 Hz外,还有40 Hz、20 Hz和16 Hz;目标4的基频组功率累积占比都很低。通过对FFS组功率占比平均值峰值搜索进一步进行筛选,可知目标1～目标3的基频分别为24 Hz、15 Hz和80 Hz。这与表1中的理论计算结果相一致。

相参积累时间是影响基频估计准确率的一个重要因素。我们将噪声标准方差设定在0.1,PRT设置为1 ms,相参脉冲积累个数范围为5～80。每个相参脉冲积累时间网格进行100次蒙特卡罗仿真得到基频估计准确率。

图3 70 ms驻留情况下的FFSFig.3 FFS of 70 ms dwell situation

表1 气动目标微动部件参数Table 1 Micro-motion parts parameters of pneumatic target

通过图4可知,要实现有效的基频估计,目标1～目标3的积累时间分别应大于43 ms、68 ms和13 ms。通过与表1比对可知在高信噪比情况下,相参积累时间略大于一个微动周期,即一个基频的倒数时,本文方法即可实现有效的基频估计。

图4 驻留时间对基频估计的影响Fig.4 Effect of dwell time on fundamental freguency estimation

在运算速度方面,在70 ms驻留场景下,SPICE和IAA达到收敛条件的迭代次数分别为879和1 993。FFT、SPICE和IAA在Matlab环境下运行时间分别为0.01 s、0.43 s和1.02 s(计算机主频为3.80 GHz,内存为128 G)。计算过程中采用共轭梯度算法[32]降低SPICE和IAA算法的运算量。

3.1.2 PRT参差场景下的仿真

设置PRT1为1 ms,PRT2为2 ms,PRT1和PRT2相交替,相参积累个数均有30个,驻留时间为90 ms,SPICE算法中K值设置为1 000。加入复正态白噪声,均值为零,噪声标准方差为0.1。

从图5和图6可知,在PRT参差的情况下,本文算法具有明显的优势。而IAA的基频谱在实际基频处也有尖峰,但并不明显,FFT的方法在实际基频处未形成尖峰。

图5 PRT参差情况下的谱估计结果Fig.5 Spectrum estimation result of different PRTs

图6 PRT参差情况下的FFSFig.6 FFS of different PRTs

通过进一步仿真验证,本文算法对于PRT脉组参差的场景仍然适用,且在目标微动周期保持稳定的条件下,性能不会随着PRT切换时间延长而恶化。

3.1.3 调制谱折叠对基频估计的影响

以目标1为例,设定PRT为1～5 ms,每种PRT均对应50个相参积累脉冲个数,加入零均值高斯复噪声,标准方差为0.1,进行100次蒙特卡罗仿真。

表2对比了本文方法和SPICE谱估计与倒谱相结合的FFS估计方法,可以看到脉冲重复周期为1 ms,调制谱不存在折叠时,本文基频估计准确率比SPICE估计与倒谱相结合的方法更高。随着PRT的变大,折叠次数的增加,倒谱的性能急剧降低,在PRT大于3 ms时,其基频估计准确率恶化到15%。而本文方法对于折叠次数并不敏感,PRT为4 ms,仍能够得到高准确率的基频估计结果。

表2 PRT对基频估计准确率的影响

倒谱并未利用基频的先验信息,因此即使在微多普勒频率未折叠的情况下,基频估计准确率也低于本文算法。从表2的准确率对比来看,本文方法基频估计性能更优。

SPICE算法在谱稀疏的前提下才有效,因此随着折叠程度的持续增加,谱稀疏度降低,会影响SPICE谱估计,进而影响FFS估计。

3.1.4 噪声对基频估计的影响

噪声是影响基频估计准确率的重要因素,本节分析不同噪声环境下本文方法对目标1～目标3的基频估计准确率。

设置PRT为1 ms,相参脉冲个数为70,相参积累时间为70 ms。加入复正态白噪声,均值为零。噪声标准方差范围为0.1～10,每种噪声环境进行100次蒙特卡罗仿真,并统计基频估计准确率。从图7可以看出,与IAA和FFT相比,本文方法具有更好的噪声鲁棒性。

图7 噪声对基频估计的影响Fig.7 Effect of noise on fundamental frequency estimation

3.1.5 分类识别性能分析

本节通过5种场景的仿真分析基频估计对识别算法的影响。

场景1的设定同第3.1.1节场景设定,即PRT为1 ms,脉冲个数为70;场景2的设定同第3.1.2节场景设定;场景3中PRT为2 ms,脉冲个数为50个。高斯复噪声均值为零,方差为0.1。进行1 200次蒙特卡罗仿真生成雷达回波,其中1 000个回波作为训练样本,200个回波作为测试样本。

交叉场景的设定是为了分析目标识别的泛化性。场景1、场景2交叉是指,用场景1的样本训练出分类器,对场景2的测试样本进行分类。

分类器包括随机森林和支持向量机(support vector machine, SVM)。其中,随机森林分类器中决策树个数设置为50个,并采用袋外错误评估特征重要性;SVM算法选用高斯核函数。输入分类器的特征除了基频外,还包括时域幅度谱归一化方差,频域波形熵、峰值比、能量比、谱线个数、中心矩(2～5阶),以及时频域原点矩、中心距[14]。本文方法的时域和时频域特征均与文献[14]一致,频域特征则基于超分辨谱重构结果进行提取得到。

从表3可以看到,场景1～场景3中,在相同的分类器框架下,本文方法对于3类气动目标的分类准确率均高于传统方法;当场景交叉时,传统方法的识别率显著恶化,而本文方法在随机森林分类器框架下性能没有明显变化。对比随机森林和SVM两型分类器,在场景不交叉时,两者和本文方法相结合识别性能基本一致;场景交叉时,随机森林的准确率比SVM更高。

表3 分类准确率比较

随机森林分类器在训练过程中可自动进行特征重要性评估,给予稳健准确的特征高重要性评分,给予其他特征低重要性评分,从而获得较好的分类识别性能。在场景1～3下训练出的随机森林分类器都将调制谱间隔特征重要性排序最高,从而获得最优的目标识别性能。而SVM更适用于构造多维特征的最优投影空间。

本文方法通过SPICE算法能够重建调制谱,从而获取稳健准确的频域特征,有利于改善后续的分类效果。同时,本文的特征提取方法与随机森林分类器相结合可以获得更优的目标识别性能。

3.2 实测数据验证

实测数据由某甚高频波段雷达外场采集得到。脉冲重复周期为1 ms。目标为AS350和Y5两型飞机,AS350型飞机的主桨叶数为3,L1为0 m,L2为5.3 m,wr为7 rps,因此基频为21 Hz。Y5型飞机为表1的目标3,具体参数见表1。两型飞机均向站飞行,仰角约为10°。

雷达原始回波经过主体移动至零频并对主体进行消卷积和滤波[33]等预处理后,再进行SPICE谱估计和基频组功率累积。图8(a)和图9(a)由相参积累个数为55,也就是驻留为55 ms的回波生成。图8(b)和图9(b)由相参积累个数为18,也就是驻留为18 ms的回波生成。其中,图9(b)中基频20 Hz对应的基频组功率累积占比比80 Hz稍高,通过进一步的FFS组功率占比平均值峰值搜索,可得到Y5型飞机的实际基频为80 Hz。因此,通过FFS估计可以得到AS350型飞机和Y5型飞机的基频分别为21 Hz和80 Hz,与目标微动特性一致。

图8 实测数据的谱估计结果Fig.8 Spectrum estimation result based on measured data

图9 实测数据的FFSFig.9 FFS based on measured data

通过第3.1.1节中的仿真和本节的实测数据分析,可以得出结论,在相参积累时间大于1.5倍微动周期的情况下,本文方法可以实现稳定的基频估计。

4 结 论

本文针对微多普勒参数在雷达波形资源受限场景下估计困难的问题,对调制谱的形成机理进行分析,根据其稀疏性的特点引入SPICE谱估计算法,并根据调制谱间隔特性进行基频组功率累积,从而实现基频的稳健准确估计。SPICE谱估计可改善短驻留、脉冲参差、调制谱折叠以及强噪声等场景下的调制谱质量,基频组功率累积则能够解决调制谱折叠的问题,并进一步提高基频参数估计的噪声鲁棒性。

稳健准确的特征提取是精确目标分类的基础,本文方法有利于后续提高三大类气动目标的识别准确率和分类器泛化能力,并能进一步对部分机型进行判别。