0 引言

直升机的旋翼旋转是引起微多普勒效应的典型运动，旋翼微动信息中包含着旋翼的转速信息，利用旋翼转速可以对直升机一类目标进行识别，并且通过旋翼转速估计可以实现旋翼类自旋目标的ISAR成像。因此转速估计已引起了众多学者的关注和研究[1-4]。而宽带线性调频连续波(Wideband Linear Frequency Modulation Continuous Wave, WLFMCW)雷达具有体积小、重量轻、成本低、功耗低等优势，利用它对直升机目标探测可以大大降低雷达设备的复杂度，有利于预警装备的快速部署和使用。

WLFMCW雷达虽然有体制上的优势，但是很多方面异于脉冲体制。首先，WLFMCW雷达在整个周期内都发射信号，旋翼目标不能忽略脉冲内的运动，因此回波中引入了随快时间变化的相位项。当前脉冲体制雷达对旋翼转速估计的主要方法有正交匹配追踪算法、高阶矩函数分析法、Hough变换法[1]、函数构建法[2]和自相关法[3-4]。其中，正交匹配追踪算法、高阶矩函数分析法需要用到回波的相位信息，因此无法直接应用到WLFMCW雷达中；Hough变换法存在着计算量大的问题，无法快速进行转速估计；函数构建法虽然计算量较小，但是在旋翼散射点均匀分布时会失效，难以准确得到转速值；自相关法必须在叶片个数已知的情况下才能准确估计旋翼转速。

综上所述，需要研究一种新的用于WLFMCW雷达旋翼转速快速估计方法。本文从特征曲线图像域出发，充分利用旋翼叶片个数的有限性，首先在自相关基础上得到图像周期，其次通过搜索叶片个数的方法求取相应的转速值，最后利用实数逆Radon变换(Real-Valued Inverse Radon Transform, RIRT)[5]对旋翼进行搜索成像并记录熵值，利用图像最小熵思想来估计旋翼转速。该方法运算量小、方法简便，对于目标模型普适性强。搜索算法在硬件上可以实现并行处理，可实现旋翼转速快速估计，为目标的快速识别奠定基础。

1 WLFMCW雷达的直升机旋翼回波建模与特征分析

宽带线性调频连续波的微动特征曲线会在频率-慢时间域成“正弦”走动，通过对“正弦”曲线周期的分析可以提取出旋翼的转速。下面先对直升机旋翼的回波进行建模和特征分析。

1.1 直升机旋翼回波建模

在远场条件下，直升机回波平动补偿后，此时直升机可等效为悬停状态，如图1(a)所示，设雷达到旋翼中心的距离为RC，方位角和俯仰角分别为α和β，旋翼以恒定的角速度ω=2πfrot绕旋翼中心旋转(frot为旋翼旋转频率)。为了便于分析，考虑α=0°，β=0°的情况，即此时假设雷达、旋翼都投影到一个平面内，如图1(b)所示。旋翼某一个叶片上的散射点P到旋翼中心的距离为r(0≤r≤l，l为叶片长度)，P点到雷达的距离为RP，P点以角速度ω绕C旋转，初始时刻P点的初始旋转角为θ。目标的旋转中心C点位于x轴上，其初始坐标为(xC,0)。选取C为参考点，那么Rref=RC。

若雷达发射锯齿波调频连续波信号，则WLFMCW雷达回波信号为

exp[j2π(fc(t-τp)+

(1)

参考信号为

exp[j2π(fc(t-τref)+

k

(2)

式中：Tp为连续波调制周期；fc为载频；k为调频斜率；为快时间；tm为慢时间，tm=mTp(m表示第m个回波，共有M个回波)；t为全时间，且有为雷达回波的时延，为散射点P到雷达的距离；τref为参考时间，τref=2Rref/c。对回波信号Dechirp后，差频信号[6-8]为

exp[Φ1+Φ2+Φ3+Φ4+Φ5+Φ6]

(3)

式中，-j4πkΔRΦ4=j8πkΔRRref/c2，Φ5=j4πkvr(tm)·

通过对式(3)中的相位项快时间微分，可得微多普勒信号的瞬时频率[6-7]：

fr=

sin(ωtm+θ)+4πkRrefvr(tm)/c2

(4)

式中， A=

ψ=arctan[-kλ/(cω)]

式(4)中第一项为旋翼目标旋转引起的走动项，它随慢时间变化，导致特征曲线在频率-慢时间上出现“正弦”变化；第二项包含快时间，为主瓣展宽项；第三项为速度-距离耦合项，它导致主瓣位置的走动。从结果可以看出，WLFMCW雷达的在频率-慢时间域的曲线振幅会被放大，放大倍数为rP=(k2λ2+c2ω2)1/2/(kλ)[6-7]。

1.2 直升机旋翼特征曲线特征分析

通过上述分析可以得知，WLFMCW雷达微动特征曲线在频率-慢时间域以“正弦”规律变化。从式(4)可以看出，特征曲线周期与旋翼旋转周期是相同的，因此可以通过利用微动特征曲线提取出旋翼转速。考虑到旋翼具有k个叶片，相邻叶片之间存在着固定的夹角2π/k，因此相邻叶片的特征曲线也存在2π/k的相位差。由此可知，理想状态下旋翼回波特征曲线可以视为由单个叶片的回波在一个周期内平移k次后经叠加而形成的。图2分别是三叶片和四叶片理论上的微动特征曲线示意图。

由图2可知，图像周期为T=2π/(ωk)，而旋翼旋转周期为Tr=2π/ω。因此真实的旋翼周期与图像周期存在着整数倍的关系：Tr=kT。令ωmax=2π/T，ω=2π/Tr，对于旋翼转速ω可表示为

ω=ωmax/k

(5)

从式(5)可知，旋翼转速的估计需提取微动特征曲线的周期以及估计旋翼叶片的个数。本文提出一种基于该原理的旋翼转速快速估计方法。

2 旋翼转速快速估计方法

该方法共分两步进行，第一步通过图像骨架提取的方法估计出图像周期，第二步在第一步基础上通过最小熵思想，利用旋翼成像后的图像熵信息来估计旋翼叶片个数，最终通过旋翼转速与图像周期和叶片个数的关系估计旋翼转速。

2.1 基于骨架特征曲线的图像周期估计

首先用文献[1,8]中的方法，对回波信号进行预处理，将差频信号在距离向上进行傅里叶变换，得到频率-慢时间域的特征曲线。在特征曲线的图像上进行降噪和平滑处理，消除奇点对估计的影响。然后将回波曲线进行二值化处理，并提取特征曲线骨架[8]。此时的特征曲线是二值化后的灰度图，该图像矩阵中只包含“0”和“1”两个值(“1”表示所提取出的骨架)。将骨架上点的位置全部记录后，找出曲线的极值位置，提取出图像周期T，从而可确定ωmax=2π/T。

2.2 基于图像最小熵的叶片个数估计方法

第二步估计旋翼叶片个数。利用RIRT算法对旋翼进行成像，得到旋翼图像信息并计算图像熵值，利用最小熵原理，搜索图像熵值最小时对应的叶片个数，再通过式(5)计算得到真实的转速。分析得知，如果对转速进行全域搜索，转速估计所需计算量较大，搜索时间也将变长，这不利于目标的快速识别。将转速的估计问题转化为对旋翼叶片个数的估计问题，可以大幅降低估计量。直升机的叶片个数通常为2～8个，搜索次数仅为7次，相比于其他搜索方式，搜索叶片个数的方法可以使计算量大为降低。

从文献[9]可知，逆Radon变换的表达式为

(6)

其中，在旋翼目标回波模型中离散化为θ=ωtm=mωTp，ρ=nΔρ，n为距离向上的第n个采样点，共有N个采样值。对时频分布图像而言，纵轴表示瞬时频率：Δρ=1/Ta(Hz)，Ta为信号积累时间。

成像完成后，假设图像的表达式为P(n,m,k)，它是关于旋翼叶片个数k的函数。图像熵的表达式[2]为

entropy(k)=

(7)

估计最优值为

(8)

因此可以得到旋翼转速的估计表达式：

(9)

通过该方法可以快速估计出旋翼转速并可以直接进行成像，为下一步识别工作奠定基础。快速转速估计方法如图3所示。

3 计算量分析

假设信号预处理结束时得到的维度是N×M的矩阵，其中N是快时间方向上的采样点数，M是慢时间上回波个数。本文转速估计算法总共分为骨架图像周期提取和叶片数估计两步。对于第一步骨架图像周期估计，使用到了寻找极大值点的方法搜索峰值点，此时的计算量为O(M-1)。

对于第二步搜索成像，使用到的成像方法为实数逆Radon变换，文献[10]中分别分析了其传统和改进型算法，本文采用成熟的BP算法(如果需要进一步提升计算速度，可以将其更换为速度更快的算法[10])。在基于BP算法的逆Radon变换中，将图像各个方向进行投影，在45°方向投影时，原本长度为N的向量长度L会变为因此距离向上的维度发生了降维。图中成像区域的维度为Na×Nb，为了得到好的成像效果，对维度进行限定，L=Na=Nb。在此基础上，开始对旋翼叶片个数进行搜索，搜索的叶片数为2～8个，那么BP算法计算量为O(7·M·Na·Nb)，即为O(3.5·N2·M)。相比于第二步，第一步计算量可忽略不计。

分析文献[2]中的计算量，流程如图4所示。

估计共分为两步，第一步是构建函数并用离散傅里叶变换进行转速粗估计。其中函数构建是在频率-慢时间域矩阵的基础上对每一行快时间乘以相应系数并求和，此时计算量为O(2N·M)；接着对构建出的函数进行FFT变换，其计算量为O(Mlog2M)，第一步总计算量为O(2I·M+Mlog2M)。第二步为相位补偿搜索计算，假设需要搜索的相位个数为K，慢时间个数为M，那么对回波矩阵中所有点的位置需要乘以补偿相位，并将对应不同慢时间的矩阵相加，最后进行旋翼转速搜索，此时计算量为O(K·I·M2)。相比于第二步，第一步计算量可忽略不计。

对于Hough变换提取微动参数方法，其计算量较大。假设在频率-慢时间域中存在正弦曲线为

f=rcos(ωt+φ)

(10)

式(10)需要3个参数来加以描述，分别是幅度r、角速度ω和初相φ。使用Hough变换提取微动参数时分别对3个参数进行搜索[10]。其计算量为M·N·Rr·Rω·Rφ,其中Rr，Rω，Rφ分别是3种参数的搜索次数，它们的大小与搜索范围和搜索步长相关。可以看出，Hough变换的计算量远远大于本文方法和文献[6]中的方法。

表1中，对比了这3种算法的运算量。

从表1可以看出，本文相比于文献[6]中方法的第一步运算量很小，几乎可以忽略不计。分析第二步运算量，一般情况下N，M数值在同一数量级上，但是K的值要大于3.5倍，乃至数十倍。因此文献[6]中的方法在第二步运算量上大于本文方法，理论上本文方法计算量更小。与Hough变换法相比，本文方法计算量远远小于Hough变换计算量。并且散射点散射强度是否平稳不会对本文方法产生影响，因此本文方法更具备普适性。由于本文在搜索时搜索次数较少，而且搜索结果互不影响，因此可以实现并行处理，更加提升了运算速度。

4 仿真与分析

本文主要考虑旋翼数未知情况下的转速估计问题，因此主要和同样不需要知道旋翼个数的文献[6]方法进行比较；自相关算法需要已知旋翼个数，本文不予比较。为综合验证本文方法的有效性和适用性，这里设置两组对比实验：一是散射点散射强度平稳时的转速估计，所有散射点的散射强度均为“1”；二是散射点散射强度非平稳时的转速估计，其中散射点的散射强度随机。旋翼散射点模型如图5所示，分别对应三叶片旋翼、四叶片旋翼和五叶片旋翼，同一叶片上的相邻散射点间隔为1 m。

雷达参数如下：雷达载频10 GHz，带宽150 MHz，时宽2 ms，信噪比20 dB，观测时间1 s。

仿真1 散射点强度平稳时旋翼转速估计

各类旋翼回波Dechirp后的频率-慢时间曲线如图6所示，为了便于观察，取前0.5 s曲线进行显示。

图6(a)～图6(c)分别是三叶片、四叶片和五叶片模型在转速为15π rad/s时的频率-慢时间域图像；图6(d)～图6(f)是这3种模型在转速为20π rad/s时的频率-慢时间域图像。

对图6曲线骨架提取后，得到骨架特征曲线如图7所示。

图7(a)～图7(c)分别是三叶片、四叶片和五叶片模型在转速为15π rad/s时的特征骨架曲线；图7(d)～图7(f)这3种模型在转速为20π rad/s时的特征骨架曲线。

在特征曲线骨架基础上，提取出不同叶片数对应的旋翼最大转速。将不同转速代入式(6)中，对图6进行RIRT变换，而后通过式(7)计算相应熵值，得到如图8所示的熵值曲线。

通过搜索图像熵值最低时对应的叶片个数，得到旋翼的叶片个数，最后通过式(8)得到旋翼转速。用该转速进行成像，得到的成像结果如图9所示。

从图9可以看出，在不同转速下旋翼的成像结果所显示的叶片长度有所变化，这与式(4)的分析是一致的。WLFMCW雷达幅度放大倍数为rP=(k2λ2+c2ω2)1/2/(kλ)。所以在转速为15π rad/s的情况下，放大倍数为6.362 3，此时旋翼成像的半径应该为31.811 3 m；在转速为20π rad/s的情况下，放大倍数为8.437 1，此时旋翼成像的半径应该为42.185 3 m。从图中可以看出，理论与仿真结果是一致的。

在散射点散射强度相同的条件下，将本文方法与文献[6]中的方法进行对比，对比结果如表2所示。这里需要说明的是，由于计算机性能的限制，在Hough变换方法中对转速的范围进行了限定，其取值范围是[12π,22π]。

通过仿真可以发现，当旋翼散射点散射强度相同时，文献[6]的估计方法已经失效。但本文方法仍可以准确估计出旋翼转速，说明本文方法在模型的适应能力上更好。

仿真2 散射点散射强度随机时旋翼转速估计

图10(a)～图10(c)分别是三叶片、四叶片和五叶片模型在转速为15π rad/s时的频率-慢时间域图像；图10(d)～图10(f)是这3种模型在转速为20π rad/s时的频率-慢时间域图像。

图11(a)～图11(c)分别是三叶片、四叶片和五叶片模型在转速为15π rad/s时的特征骨架曲线；图11(d)～图11(f)是这3种模型在转速为20π rad/s时的特征骨架曲线。

图12(a)～图12(c)为不同叶片数旋翼的图像值曲线，通过最小熵原理可以估计出旋翼的叶片个数，进而可以估计旋翼转速。图13(a)～图13(c)是旋翼转速为15π rad/s时三叶片、四叶片和五叶片旋翼的成像结果；图13(d)～图13(f)是旋翼转速为20π rad/s时三叶片、四叶片和五叶片旋翼的成像结果。

通过表3可以看出，在散射点散射强度不同的条件下，旋翼转速通过文献[6]的方法可以估计出来。两种方法估计精度相近，但本文方法运行速度更快、所需时间更少。通过以上两种情况的对比，说明了本文方法在保证估计精度的基础上具有更好普适性和小计算量的优势。

通过以上两组仿真实验可以发现，通过本文方法，可以对不同叶片个数的旋翼进行较为精确的转速估计和旋翼成像，在散射点散射强度平稳和非平稳的条件下本文方法都可进行直升机旋翼转速的估计，验证了本文方法的有效性。

5 结束语

本文提出了一种基于骨架提取的旋翼转速快速估计方法，通过提取骨架相邻峰值的相位，得到旋翼转速可能的最大值，接着对不同叶片个数进行搜索成像，利用最小熵原理得到成像效果聚焦性最好的图像，此时所对应的转速就是估计的最佳转速。通过这个方法，可以快速估计出旋翼转速，有效提升目标识别的实时性。但是，该方法是基于骨架提取结果完成的，图像域上的骨架提取容易受到噪声影响，在低信噪比条件下该方法成像效果将下降甚至失效。下一步将就如何在低信噪比条件下快速估计进行研究。

参考文献：

[1] 陈永彬,李少东,杨军,等. 一种旋翼叶片微动特征提取新方法[J]. 雷达科学与技术, 2017，15(1):13-18.

CHEN Yongbin, LI Shaodong, YANG Jun, et al. A New Method for Micro-Motion Signature Extraction of Rotor Blades[J]. Radar Science and Technology, 2017，15(1):13-18. (in Chinese)

[2] 华煜明, 赵华, 郭军海. 高速自旋目标ISAR成像的代数重构方法[J]. 雷达科学与技术, 2016，14(2):134-139.

HUA Yuming, ZHAO Hua, GUO Junhai. Algebraic Iterative Imaging Method for Rapidly Rotating Target in ISAR[J]. Radar Science and Technology, 2016, 14(2):134-139. (in Chinese)

[3] 李靖卿,冯存前,贺思三,等. 利用微动信息矩阵的弹道目标特征提取方法[J]. 信号处理, 2016, 32(4):488-495.

[4] 陈永彬,李少东,杨军,等. 旋翼叶片回波建模与闪烁现象机理分析[J]. 物理学报, 2016, 65(13):281-291.

[5] 张群,罗迎. 雷达目标微多普勒效应[M]. 北京:国防工业出版社, 2013:47-82.

[6] 胡杰民,付耀文,胡志刚,等. 高速旋转目标旋转速度估计方法[J]. 电子与信息学报, 2009, 31(9):2069-2073.

[7] BAI Xueru, XING Mengdao, ZHOU Feng, et al. Imaging of Micromotion Targets with Rotating Parts Based on Empirical-Mode Decomposition[J]. IEEE Trans on Geoscience and Remote Sensing, 2008, 46(11):3514-3523.

[8] 李彦兵,张曦文,李飞,等. 一种大加速度机动目标微动参数估计方法[J]. 电子与信息学报, 2017, 39(1):82-87.

[9] 白雪茹. 空天目标逆合成孔径雷达成像新方法研究[D]. 西安:西安电子科技大学, 2011.

[10] 梁颖,张群,何劲,等. 调频连续波ISAR目标微多普勒特征分析[J]. 现代防御技术, 2011, 39(6):71-77.

[11] 梁颖,田韵,张群,等. 基于多项相位变换的FMCW-ISAR微多普勒特征提取方法[J]. 空军工程大学学报(自然科学版), 2012, 13(2):74-78.

[12] 陈梁栋,李少东,陈文峰,等. LFMCW的旋翼目标微动特征提取方法研究[J]. 空军预警学院学报, 2016, 30(6):391-394.

[13] 李康乐,刘永祥,姜卫东,等. 基于逆Radon变换的微动目标重构研究[J]. 雷达科学与技术, 2010, 8(1):74-79.

LI Kangle, LIU Yongxiang, JIANG Weidong, et al. Reconstruction of Target with Micro-Motions Based on Inverse Radon Transform[J]. Radar Science and Technology, 2010, 8(1):74-79. (in Chinese)

[14] 李传中,邓钰栋,苏卫民,等. 面向像素的并行快速后向投影算法[J]. 南京理工大学学报, 2014, 38(5):651-657.