直升机除质心平动外伴随的旋翼转动会对雷达回波相位产生周期性调制，这种调制现象称为微多普勒效应(micro-Doppler effect，MDE)[1]．微多普勒效应能反映出目标的几何构成和运动特性，可以用来确定目标的属性，为目标识别提供新的途径[2]，因此利用微多普勒效应对直升机进行识别是目前的一个研究热点．直升机旋翼建模方法主要有RCS模型[3-4]、积分模型[5-7]和散射点模型[8-9]三种．值得指出的是：当散射点模型中散射点均匀分布且间隔小于雷达波长的1/2时，可等效为积分模型，此时目标回波信号中会出现微动导致的闪烁现象[10]．若仅用少量散射点对目标进行等效，则回波信号中不会出现闪烁，此时微动特征的提取问题就转化为多分量周期性微多普勒信号的参数估计问题．关于多分量周期性微多普勒信号参数估计的解决方法主要可分为四类．a．基于信号分解的方法．这类方法能够将多分量信号分解成单分量信号，如通过chirplet变换[11]、AM-LFM分解[12]和经验模态分解[13]等算法将微多普勒分量从回波中分离出来，但是由于回波中各微多普勒分量交叠严重，当两个分量频率相接近时，难以进行有效分离和微动参数提取．b．基于变换域的方法．目标振动和旋转引起的微多普勒信号均表现为正弦调频(sinusoidal frequency modulation，SFM)信号形式，因此可通过构建匹配算子建立具有独特运算定义的正弦调频信号域，将回波分解在正交三角函数基上获得回波的微多普勒频谱[14]，这类方法具有较好的估计精度和抗噪性能，但是当信号分量多于3个时交叉项严重，难以根据微多普勒频谱提取微动特征．c．基于图像域的方法．这类方法运用较多，先通过不受交叉项干扰的Gabor变换[8]进行时频分析获取时频图像，再利用Hough变换[15]、iRadon变换[16]等图像处理方法将图像空间中的曲线检测问题转化为参数空间中的峰值检测问题来提取目标的微动特征，其优势在于只要能推导出目标微动在时频域中的曲线方程，就可以实现目标微动特征的提取，但这类方法严重依赖时频分析方法的性能，当回波中存在多个微动信号分量时，时频聚集性差，且运算量随参数空间维数的增加成指数增长．d．基于稀疏重构的方法．由于散射点模型的目标回波可视为少数强散射中心回波的叠加，天然具有稀疏特性，因此可以采用稀疏重构的方法分析微多普勒信号，这类方法基于微多普勒信号的正弦调制特性，构建相应的稀疏字典，通过正交匹配追踪(orthogonal matching pursuit，OMP)[17]等算法对微多普勒分量进行反演，在欠采样条件下获得了较好的估计效果，但是当搜索网格失配时鲁棒性不高．针对以上问题，本研究提出一种基于谱聚集度指标的旋翼目标微多普勒参数估计方法，并利用奇异值比谱法(singular value ratio，SVR)[18-19]对搜索过程进行优化，仿真结果表明该方法能够有效提取目标的微动特征，为直升机目标识别提供了新的技术途径．1 目标回波模型直升机的雷达回波主要是机身回波和旋转部件回波的矢量和，机身部分只有平动分量，而旋转部件除平动分量外还有微动分量，因此这里主要给出旋转部件的雷达回波模型．目标与雷达的位置关系如图1所示，以雷达位置O为原点建立雷达坐标系O-UVW，以直升机旋翼中心Q为原点，XY面保持水平，X轴平行于U轴建立参考坐标系Q-XYZ；同时，以Q为原点，旋翼所在平面为xy面建立目标坐标系Q-xyz．设O与Q的距离为R0，其俯仰角和方位角分别为α和β，0≤α≤π/2，为方便分析且不失一般性，假设雷达波束在β=0°时照射目标．10.13245/j.hust.238753.F001图1雷达与目标位置关系设雷达发射信号为窄带线性调频信号，记为st(t^,tm)=Arect(t^/TP)exp[2πj(fct+μt^2/2)]，(1)式中：A为发射信号幅度；Tp为脉冲宽度；fc为载频；μ为调频斜率；t^为快时间；tm=mTr为慢时间，其中，m为第m个回波，Tr为脉冲重复周期；t=t^+tm为总时间．经脉冲压缩和杂波抑制后旋翼的回波信号[8]可表示为sr(t^,tm)=∑k=1K∑i=1IσikATpsinc[B(t^-2Rik(tm)/c)]⋅exp[-4πjRik(tm)/λ]，(2)式中：K为叶片数目；I为目标单个叶片上的散射点数目；σik为散射系数；B为发射信号带宽；λ为雷达波长；函数项sinc(⋅)为距离项，包含了目标的位置信息和走动信息；函数项exp(⋅)为多普勒项，它包含了目标的多普勒信息；Rik(tm)为旋翼上散射点到雷达的距离．由图1位置关系可得出目标上散射点到雷达的距离为Rik(tm)=R0+vt+Rinitr0，(3)式中：R0为雷达坐标系与参考坐标系的距离；v为目标平动速度矢量；r0为散射点P在目标坐标系中的位置，r0=[rikcos(ωtm+θik),riksin(ωtm+θik),0]T；Rinit为初始欧拉角确定的初始旋转矩阵．按照z-x-z协定，将目标坐标系分别围绕z轴旋转ϕe，围绕x轴旋转θe，再围绕z轴旋转φe，从而可转化成参考坐标系Q-xyz，则散射点P在参考坐标系中的坐标为Rinitr0=a11rikcosΩ+a12riksinΩa21rikcosΩ+a22riksinΩa31rikcosΩ+a32riksinΩ，(4)式中：Ω=ωtm+θik，其中，ω为旋转角速度，θik为初相；0≤rik≤l为散射点距旋翼中心的距离．此时，目标的瞬时频率为fd=12πdΦ(tm)dtm=2λ[v+ddtm(Rinitr0)]Τn．(5)当目标处于远场条件时，R0≫vtm+Rinitr0，n可以近似为n=R0/R0，这就是雷达的LOS方向，式(5)中第一项是由目标平动引起的多普勒频率fd=2vΤn/λ，第二项则是由旋转引起的微多普勒频率fmD=(2/λ)d(Rinitr0)Τn/dtm，代入可得fmD=(2ωrik/λ)m2+n2sin(Ω-arctan(n/m))，(6)式中：m=a11cosα+a31sinα；n=a12cosα+a32sinα．当目标处于水平飞行或悬停状态时，即ϕe,θe,φe均为0°，此时式(6)改写为fmD=(2ωrik/λ)cosαsinΩ．由此可见：旋翼对回波呈现正弦调制特性，且以平动多普勒频移为中心分布在两侧，边带宽度与目标尺寸、转速、俯仰角和目标姿态有关．2 微动特征提取方法微动特征提取主要是通过分析回波的微多普勒效应，从中获取能反映目标属性的特征量，实现对目标结构、尺寸等参数的估计．通过构建多个解调算子来计算相应的谱聚集度，由于不同解调算子与微多普勒信号相乘得到的谱聚集度差异性很大，因此可利用此特点找到解调算子与目标微动参数之间的映射关系，估计出相应的微动参数．但是该过程须要在多维特征空间进行遍历式的搜索，运算量随特征空间的维度增加而爆发式增长．因此，引入奇异值比谱法来对搜索过程进行优化，通过检测信号的主周期来降低搜索空间的维度，为确保算法实时性提供基础．当目标仅由叶尖散射点生成时，式(2)中回波慢时间维信号由多个微多普勒信号分量组成，即s(tm)=∑k=1Ksk(tm)，分量个数与叶片数量一致，根据其相位项规律定义一个解调算子Φ(ȃ,ω̑,φ̑)=exp(ja⌢cos(ω⌢tm+φ⌢))，式中：a⌢为叶片尺寸算子，取值范围为4~10 m；ω⌢为叶片转速算子，取值范围为3~7 Hz；φ⌢为叶片初相算子，取值范围为0~2π．同时，定义信号的离散傅里叶谱聚集度(concentration of spectrum，CS)为γCS=E1N2∑k=1K∑n=0N-1Xk(n)4，(7)式中：E(⋅)为均值函数；Xk(n)为离散形式下信号sk(tm)第n个频点的傅里叶系数；N为信号长度．此时单一分量信号与解调算子乘积的频谱为X(f)=∫sk(tm)Φ(a,ω,φ)e-j2πfdtm．根据Parseval定理[20]有∫X(f)2df=∫sk(tm)2dtm．即引入解调算子前后微多普勒信号的能量不会发生变化，且回波中K个微多普勒分量是正交的，因此有∫s(tm)2dtm=∫∑k=1Ksk(tm)2dtm=∑k=1K∫sk(tm)2dtm，即回波信号s(tm)的能量等于各微多普勒分量能量之和．当解调算子与微多普勒分量均不匹配时，按照式(7)可计算出此时的谱聚集度指标为γCS1=E1N2∑k=1K∑l=0N-1Xk(l)4，式中Xk(l)为微多普勒信号与解调算子乘积在离散形式下第l个频点的傅里叶系数．当解调算子与回波中第k'个微多普勒分量匹配而与其他K-1个分量不匹配时，i'表示不与解调算子匹配的微动分量，在第l'个频点的傅里叶系数为Xi''(l')，那么第k'个分量在频域中仅在零频处有值，记为X'(0)．由于解调算子并不影响信号总能量大小，因此根据Parseval定理有1NX'(0)2=1N∑l'=0N-1Xk''(l')2．那么根据式(7)可得γCS2=E1N2∑i'=1,i'≠k'K-1∑l=0N-1Xi''(l')4+1N2X'(0)4≫E1N2∑i'=1,i'≠k'K-1∑l=0N'-1Xi''(l')4+1N2∑l=0N-1Xk''(l)4=E1N2∑k'=1K∑l=0N-1Xk''(l)4=γCS1．可见解调算子与第k'个微动分量信号匹配时，γCS2≫γCS1．解调后微多普勒信号的频谱如图2所示．图2(a)为微多普勒信号的原始频谱，为验证谱聚集度指标的可行性，假设解调算子中的a,φ参数已成功匹配，通过改变ω来模拟解调算子不匹配、准匹配和完全匹配的情形．图2(b)、(c)和(d)分别为ω=10.2π,10.02π,10.00π解调后微多普勒信号的频谱．图2(b)中解调算子不匹配，频谱杂项较多，严重影响了极值点的判断；图2(c)中解调算子参数与真值误差较小，属于准匹配场景，此时已经可观察到显著的极值点；图2(d)中解调算子完全匹配，零频处出现强峰值，与不匹配分量在频谱中对应的带宽和幅值具有较大差异．10.13245/j.hust.238753.F002图2解调后微多普勒信号频谱由此可见：当微动回波信号s(tm)中某一个微多普勒分量与解调算子匹配时，其谱聚集度远大于解调算子与任意微动分量均失配时的谱聚集度，且锐化程度很高．此时根据解调算子的参数即可估计出目标叶片数、叶片长度、旋转频率及叶片初相．但是若直接在三维空间对微动参数进行遍历搜索，运算量过大，因此还引入了一种优化算法对搜索空间进行降维处理．SVR算法已广泛用于信号主周期分量的估计，本研究通过该算法估计出目标旋转速度ω^，并将其作为先验信息，此时只要在Φ{a,φ}空间中找出所有γCS的极大值，便可得出对应的a^和φ^．由于每个微多普勒分量相位项中初相角有差异，必然存在多组符合要求的估计值，其数量即为目标叶片数．基于谱聚集度指标的参数估计结果精度很高，可对a所在数据行逐列搜索，并在aoφ平面中逐次消除所有极大值，依次估计出旋翼微动参数，算法具体步骤如下．步骤1 依据旋翼转动频率的先验信息，以步长Δf设定频率的变化范围[fmin,fmax]，对范围内的每一个频率fr依据L=round(fs/fr)计算得到(fmax-fmin)/Δf+1个L的值．步骤2 将旋翼目标所在距离单元内的长度为Ns的慢时间维微动回波数据转换成大小为GL阶的二维数据矩阵，其中G=floor(Ns/L)-1．步骤3 对每个二维数据矩阵作奇异值分解，求得各个矩阵的奇异值比，其中最大的奇异值比谱值对应的频率即为旋翼的转动频率，将对应的目标旋转速度记为ω^．步骤4 以步长Δa和Δφ分别设定旋翼长度、叶片初相的搜索范围为[amin,amax]和[φmin,φmax]，结合ω^可以构建[(amax-amin)/Δa+1]×[(φmax-φmin)/Δφ+1]个解调算子exp(jacos(ω^tm+φ))．步骤5 根据谱聚集度的定义计算不同解调算子下旋翼目标微多普勒信号的谱聚集度．步骤6 初始化参数P=2，首先找出aoφ空间中最大的两个γCS，得出旋翼叶片长度估计值a^1和a^2，以及初相位估计值φ^1和φ^2，以(a^1,φ^1)，(a^2,φ^2)坐标为圆心，将半径3个单元内的γCS置零，再次计算此时aoφ平面中最大的γCS，记录此时叶片长度估计值a^3及初相估计值φ^3，由于旋翼叶片连续两个叶片之间的相位差相等，若φ^1，φ^2和φ^3之间的差值不满足该条件，则叶片数的估计为K^=P；若φ^1，φ⌢2和φ⌢3之间的差值满足该条件，则进一步执行置零操作，更新P=P+1，计算a^P，φ^P及连续叶片之间差值，直到搜索出所有符合条件的微动参数，搜索终止时的P值即为叶片数估计值．3 仿真分析为分析本文方法对目标参数的估计性能，设置仿真目标参数为叶片数3个、转动频率5.15 Hz、叶片长度6.1 m、初相为(30°,150°,270°)．参考某窄带雷达确定仿真雷达参数，载频1 GHz，重复频率4 kHz，脉宽100 μs，带宽1 MHz．仿真中信噪比(signal to noise ratio，SNR)均为脉冲压缩后的数值．仿真电脑搭载64 bit Windows 10系统，Intel Core i7-8550U CPU@1.80 GHz处理器，运行内存16 GiB，仿真软件版本为Matlab R2016a．本文方法先利用奇异值比谱法对仿真回波慢时间维数据展开分析，设置转速搜索范围为[8π,14π]，以步进0.02π对范围内全部采样点进行遍历式搜索，其目标转速搜索结果如图3所示．在图3中，仅在f^0=5.15 Hz处出现奇异值比极大值，与其他采样点上的奇异值比差异明显，能准确估计出目标微多普勒信号中的主周期分量，即目标转速．10.13245/j.hust.238753.F003图3转速估计结果在准确估计出转速的前提下，对三维搜索空间进行降维处理，图4为二维搜索空间中叶片长度和初相的估计结果．图4(a)是本文方法的估计结果，图4(b)和(c)分别是OMP和Hough变换(HT)的估计结果．为对比上述3种方法的性能，将CS，OMP和HT方法中叶片长度搜索范围均设置为4~10 m，步进为0.1 m，初相搜索范围均设置为[0°,360°]，步进为1°．10.13245/j.hust.238753.F004图4叶片数、叶片长度和初相估计结果由图4可见：上述3种方法的估计结果中均能看到3个强峰值，即出现了3组符合要求的估计值，这意味着回波信号中存在3个不同的初相角，即叶片数为3，与预设值一致．在搜索参数设置相同的情况下，本文、OMP和HT方法对目标微动参数的估计值保持一致，叶片长度均为6.1 m，初相为(30°,150°,270°)，与预设值完全一致，说明在搜索网格设置合理的情况下，这3种方法均能完成微动特征提取，但本文方法和OMP方法的锐化程度更高，更易于分辨．通常而言，利用时频分析方法进行直升机旋翼微动特征提取时，大都须提前区分叶片数量的奇偶性，这是由于当叶片数为奇数时，微多普勒频率在正频部分和负频部分的正弦曲线不对称；而当叶片数为偶数时，若某一个叶片出现正频最大值，与之对称的叶片必定出现负频最大值，即正频部分和负频部分的正弦曲线对称．这表明叶片数奇偶性不同时，呈现出的规律不尽相同，须区分讨论．但是从理论分析可知：构建的解调算子并不关注叶片的奇偶性，只要慢时间维回波与解调算子构成匹配关系即可估计出目标的属性，且有几个叶片就能构成几组匹配关系，即aoφ空间出现几个峰值，因此这里不对叶片数为偶数的情况进行赘述．为从多个角度说明本文方法的优越性，进一步讨论了上述3种算法的运算时间、信噪比适应性．通过100次蒙特卡罗实验，SVR算法的平均耗时为7.1 s，但是上述方法均采用了SVR来进行优化，因此不考虑SVR算法的耗时．由于HT属于图像处理方法，须要获取回波慢时间维的时频信息，因此在讨论这种方法的运算时间时还须加上时频分析方法的耗时，上述算法的平均运算时长分别为5.06，89.12和7.10 s．从算法耗时来看，OMP，HT和本文方法依次递减，即所提出的方法方法运算速度最快．常规窄带雷达的回波在没有干扰的情况下，脉冲压缩后的信噪比通常为5~8 dB，但是实际复杂电磁环境下会低得多，因此这里仅讨论脉压后信噪比小于0 dB的场景．图5为不同信噪比条件下的微多普勒参数估计误差，图5(a)为叶片长度估计误差，图5(b)为初相估计误差．必须指出的是：SVR算法在信噪比低于-10 dB时失效，因此在SNR小于-10 dB时，假设转速估计准确．从图5可见：OMP在信噪比低于-19 dB时失效，本文方法在信噪比低于-17 dB时失效，HT在信噪比低于-16 dB时失效，但是在-13~-16 dB范围内存在不小的误差，因此从信噪比适应性能来看，OMP、本文和HT方法依次减弱．10.13245/j.hust.238753.F005图5估计误差为验证网格失配情况下上述方法的性能，通过将叶片长度和初相分别设置为6.13 m和30.7°，让搜索网格失配．仿真发现本文方法和OMP方法性能均下降了5~6 dB，其估计误差与偏离网格程度直接相关．而HT方法基本不受网格影响，这得益于该方法利用的是时频图像信息，时频曲线具有一定展宽，在参数变化不明显的情况下，其时频结果基本相同．从仿真分析结果可以看出：SVR算法能够准确估计出微多普勒信号中的主周期，即目标旋转速度，本文方法相比于OMP和HT方法具有更快的运算速度，相比HT方法具有更好的信噪比适应性能，但略逊于OMP方法．在搜索网格不失配情况下，估计误差为0，与OMP和HT方法一致，而当搜索网格失配时，信噪比适应性下降5~6 dB，与OMP方法相当，估计误差由偏离网格的程度决定．可以看出本文方法相比于传统方法具有一定的综合对比优势．