1.   引言

逆合成孔径雷达(Inverse Synthetic Aperture Radar, ISAR)成像技术凭借独特的全天时、全天候以及高分辨的优势，被广泛应用于战场侦察、空中预警、天文观测等领域^[1−3]。对ISAR而言，其方位分辨能力来源于多普勒效应，目标与雷达的径向相对运动引起二者距离的变化，从而影响雷达回波的时延，改变回波的频率和相位信息，多普勒带宽的大小决定了方位分辨率的高低。

然而，多普勒效应并非一定有益于ISAR成像。在实际工程实践中，目标部件相对于主体往往叠加有局部运动，例如无人机和直升机旋翼的转动、各类军舰及装甲车上天线的转动，弹道导弹弹头的进动等，这些局部运动统称为微动^[4,5]。主体运动会产生多普勒频率，微动也会产生多普勒频率，该频率会叠加到主体运动产生的多普勒频率上，产生一定的多普勒展宽，这就是微多普勒效应^[6,7]。

在早期的低分辨ISAR中，目标部件的旋转、振动等物理现象常被粗略补偿或者直接忽视。但随着雷达技术的不断发展，特别是微波光子技术的出现，目前ISAR系统的带宽已经能够达到数GHz至10 GHz^[8,9]，再通过增大合成孔径转角使得二维分辨匹配，现代超高分辨ISAR系统已经具备厘米级成像能力。如图1所示，超大带宽对距离分辨率有着巨大提升，尤其对小型目标的探测意义显著。然而，分辨率的提升也对处理精度提出了更高的要求，使得目标上的微动干扰很难再被轻视，微多普勒效应直接影响ISAR成像质量，在经过方位向压缩之后，微动部分形成方位向散焦条带，导致图像局部模糊^[10,11]，因此，微多普勒效应的抑制对于超高分辨ISAR成像有着显著的意义。

图  1  几种典型目标所占距离门数目随信号带宽的变化

Figure  1.  Variation of the number of range cells for several typical targets with bandwidth

下载: 全尺寸图片 幻灯片

微多普勒信号是典型的时变非平稳信号，而且其带宽通常远大于主体回波，常规的动目标指示以及杂波对消类算法对此效果不佳^[12]。目前，微多普勒抑制方法主要分为参数化和非参数化两大类^[13]。参数化方法需要针对具体的微动形式进行精细建模，分析微动的真实特性，通过估计微动的频率、幅度、相位等信息对回波进行补偿^[14]，对模型的准确度要求非常高。文献[15]将微多普勒信号建模为正弦调频(Sinusoidal Frequency Modulated, SFM)信号，提出了正弦调频傅里叶贝塞尔变换(Sinusoidal Frequency Modulation Fourier-Bessel Transform, SFMFBT)算法，对多成分的微多普勒分量进行了估计；文献[16]提出了一种联合雷达视线角以及距离多普勒曲线的振动参数粗估计方法，并在此基础上通过构造最小熵优化函数进行精估计，以实现微动补偿功能。相比于参数化算法，非参数化算法则通过预先分离主体回波与微多普勒信号来实现微动的抑制^[17,18]，这类方法对模型的依赖程度较小，能够更好地适应复杂真实场景。文献[19]提出一种基于独立成分分析(Independent Component Analysis, ICA)的分布式雷达组网环境下的微动消除算法，该方法可以准确而稳健地改善图像质量，但是不适用于单基ISAR系统；文献[20]提出一种基于经验模态分解(Empirical Mode Decomposition, EMD)的微多普勒抑制算法，EMD通过将信号中的极大极小值连线，自适应地把实信号分解为各个本征模态函数(Intrinsic Mode Function, IMF)，将非平稳信号平稳化^[21]。基于EMD的微多普勒抑制算法将EMD扩展到了复数域，实现对ISAR回波的分解，然后通过过零检测法筛选主体回波，达到抑制微多普勒的目的。该方法简单实用，但也有着缺陷。首先EMD本身存在模态混叠现象^[22]，进而影响后续的成像模态选取，使得微动抑制不彻底或者造成主体回波的损失；其次，过零检测法虽然去除了高频部分，但微多普勒分量在低频的残余仍然会造成图像质量下降。文献[23]提出一种基于局部均值分解(Local Mean Decomposition, LMD)的快速旋转目标微多普勒分离技术，成功分离了涡轮螺旋桨飞机实测数据中的微动信号，但是LMD同样面临模态混叠问题，而且该算法在低分辨前提下未对分解结果加以筛选，只保留了残余项作为主体回波，未能避免成像主体的损失。

为了解决上述问题，本文提出了一种基于变分模态分解(Variable Mode Decomposition, VMD)^[24]与优选的ISAR微多普勒抑制方法。通过将实数域VMD拓展到复数域，实现了超宽带ISAR方位回波数据的分解；再通过基于图像熵的迭代优化处理，实现了最优成像模态选取，抑制了微多普勒效应对成像的干扰。

2.   回波建模与微多普勒效应分析

2.1   含微动部件目标的ISAR回波建模

ISAR系统通过发射线性调频信号获得高距离分辨力，其回波可表示为

其中，t和$\eta$分别表示快时间和慢时间，${R_i}(\eta )$表示散射点i到雷达的距离历史，c表示光速，${T_{\mathrm{f}}}$表示信号脉宽，${\sigma _i}$表示后向散射系数，${f_{\mathrm{c}}}$和${K_{\mathrm{r}}}$分别表示线性调频信号的载频和调频率。

ISAR信号处理常用转台模型来对距离历史建模^[25]。如图2所示，目标从位置1运动到位置2这一过程实际上可以分成两个独立部分，其等效于从位置1平动到位置3然后转动$\Delta \theta$。由于平动对所有散射点是一致的而且平动补偿易于进行，这里忽略不计，只考虑转动项，对成像平面$xoy$内任一散射点而言，其转动量与坐标$({x_i},{y_i})$的关系可以近似为

图  2  ISAR转台模型

Figure  2.  Turntable model of ISAR

下载: 全尺寸图片 幻灯片

其中，${\omega _{\mathrm{g}}}$表示主体的转动角速度。

实际上，传统的转台模型只适合描述刚体目标，在超高分辨前提下对含有微动部件的目标进行成像时，微动的影响不容忽视。

如图3所示，以雷达为原点O建立坐标系$(X,Y,Z)$，部件转动中心位于点$O'$，其瞬时方位角和俯仰角分别为$\alpha (\eta )$和$\beta (\eta )$，${R_0}(\eta )$表示$O'$离雷达的距离历史，包含了$O'$平动和自身的转动。再以$O'$为原点建立局部坐标系$(X',Y',Z')$，点P表示旋转部件上第p个散射点，在$X'O'Y'$平面内绕$O'Z'$轴以角速度${\omega _{\mathrm{r}}}$匀速转动。构建辅助坐标系$(X'',Y'',Z'')$，其由$(X,Y,Z)$平移得到，且与坐标系$X'O'Y'$拥有共同的原点$O'$。直线$O'Q$为平面$X'O'Y'$与平面$X''O'Y''$的交线，$X'$轴与$O'Q$的夹角为${\psi _0}$，$X''$轴与$O'Q$的夹角为${\theta _0}$，$Z'$轴与$Z''$轴的夹角为${\varphi _0}$。

图  3  含微动部件目标的ISAR几何模型

Figure  3.  ISAR geometric model of targets with micro-motion components

下载: 全尺寸图片 幻灯片

假设点P的转动半径为${a_p}$，初相为${\zeta _p}$，则方位零时刻点P在$(X',Y',Z')$坐标系中的坐标为${{\boldsymbol{r}}_p}(0) = {\left[ {{a_p}\sin {\zeta _p}{\text{ }} - {a_p}\cos {\zeta _p}{\text{ }}0} \right]^{\mathrm{T}}}$，则点P的瞬时坐标为

将坐标系$(X',Y',Z')$变换到$(X'',Y'',Z'')$中，由几何构型可知需要经过3次旋转，变换矩阵分别为

由此得到P点在$(X'',Y'',Z'')$中的瞬时坐标为

又因为$O'$点在$(X,Y,Z)$中坐标为

由平移关系得到，P点在$(X,Y,Z)$中的瞬时坐标为

计算${{\boldsymbol{R}}_p}(\eta )$的模值得到P点离雷达的距离历史为

其中，${A_p} = {a_p}\cos \beta (\eta )\cos ({\theta _0} - \alpha (\eta ))$, ${B_p} = {a_p}\cos \beta (\eta ) \cdot\cos {\varphi _0}\sin ({\theta _0} - \alpha (\eta )) - {a_p}\sin \beta (\eta )\sin {\varphi _0}$。

令${C_p} = \sqrt {A_p^2 + B_p^2}$, ${\xi _p} = {\zeta _p} + {\psi _0} - \dfrac{\pi }{2} + \arctan \dfrac{{{A_p}}}{{{B_p}}}$，则式(8)可进一步简化为

其中，${R_{{\mathrm{md}}}}(\eta )$表示微动，其大小受转动频率、转动半径以及转动平面的空间角度影响。

进一步将式(9)代入式(1)并完成脉冲压缩、距离徙动校正、二次相位补偿后^[25,26]，微动散射点的ISAR回波可以表示为

由式(10)可知，微动散射点对应的包络为正弦形式，其幅值${C_p}$与微动幅度正相关，在超高分辨条件下会跨距离单元；同时微多普勒相位也表现为正弦形式，其特性将在2.2节分析。

2.2   微多普勒特性分析

设$\lambda$为载频对应的波长，对式(10)中相位求导，微动散射点的多普勒频率为

式(11)的第1项是微动散射点随主体运动而产生的，第2项是微动自身产生的，微动给回波附加了一项额外的正弦调频信号。通常情况下，由于主体转动角速度${\omega _{\mathrm{g}}}$远小于部件微动的角速度${\omega _{\mathrm{r}}}$，所以微动带来的多普勒展宽会远远大于主体部分各散射点的多普勒带宽。

以载波频率35 GHz的ISAR系统为例，忽略空间角影响，当其对臂展1 m、旋翼转速10 Hz、叶片长度15 cm的无人机成像时，假设机身二阶转动量为0.1 rad/s，则主体运动产生的最大多普勒频移仅为23 Hz，远小于旋翼微动产生的2.2 kHz微多普勒信号带宽。这里利用短时傅里叶变换(Short Time Fourier Transform, STFT)来展示微多普勒信号的时频特性，结果如图4(a)所示，微多普勒在时频域表现为正弦曲线。

图  4  微多普勒信号的STFT谱和频谱图

Figure  4.  STFT spectrogram and spectrum of the micro-Doppler signal

下载: 全尺寸图片 幻灯片

为了进一步分析微多普勒对成像的影响，对式(10)做贝塞尔函数展开：

忽略常数相位$\xi$，将式(12)沿方位向做${\text{FFT}}$，得到二维图像：

其中，${J}_{k }(\cdot)$表示第k阶类贝塞尔函数，${T_{\mathrm{m}}}$表示方位积累时间，${f_{\mathrm{r}}} = {\omega _{\mathrm{r}}}/(2\pi )$。微多普勒的频谱由如图4(b)所示的间隔${f_{\mathrm{r}}}$的谱线组成，在成像结果中表现为散布在微动中心两侧的干扰条带，导致图像局部模糊。

3.   基于变分模态分解与优选的ISAR微多普勒抑制

根据第2节的分析，目标部件的微动会产生额外的微多普勒带宽，导致方位向出现散焦条带。由式(10)可知，回波中的微多普勒信号本质上是多种正弦调频信号的混合，具有带宽大、时变非平稳的特点。为了将微多普勒信号从回波中去除，传统的EMD算法通过极值点拟合的方式自适应地将非平稳回波分解为多个平稳的模函数，这种分解方式在回波成分较复杂时易受模态混叠影响，不适合超高分辨场景对分解精度的要求。本文引入了VMD算法^[24]来克服这一问题，VMD通过增加正则化项、变分优化等手段，提高了泛化能力，基本克服了模态混叠带来的干扰，比传统的EMD分解结果更为准确。

VMD可以将一个实信号$s(t)$分解成N个中心频率${\omega _k}$从低到高分布的平稳模函数${u_k}(t)$，每个模函数可看成频率单一的窄带信号，将信号分解问题转化成式(14)限制性变分问题：

通过引入惩罚因子${\alpha _{\mathrm{p}}}$和拉格朗日因子${\tilde \lambda _1}$将问题再次转化为无限制性变分问题，其表达式如下：

式(15)可以使用交替方向乘子法(Alternate Direction Method of Multipliers, ADMM)进行求解。在第$n + 1$次迭代中，通过更新$u_k^{n + 1}$, $\omega _k^{n + 1}$和$\lambda _1^{n + 1}$寻找无限制性变分问题的最优解，更新过程主要包括以下3个步骤：

步骤1　在频域更新模函数：

其中，$\omega$为角频率，$\tilde s(\omega )$是$s(t)$的傅里叶变换，$\tilde u_i^n(\omega )$表示${u_i}(t)$的傅里叶变换在迭代n轮后的值，本质是对剩余量$\tilde s(\omega ) -\sum\nolimits_{i \ne k} {{{\tilde u}_i}} (\omega )$进行维纳滤波。

步骤2　更新中心频率：

步骤3　更新拉格朗日因子：

其中，${\tau _1}$是${\tilde \lambda _1}$的更新参数。迭代的结束条件为

其中，$\epsilon$为收敛条件满足的阈值。

在ISAR系统中一般采用正交解调的方式获取回波的同相和正交分量，最终处理的回波数据是包含幅度和相位信息的复数信号。因此，有必要对VMD进行扩展，使其能够应用于复信号分解这一场景。

文献[27]提出一种CVMD算法将复信号实部和虚部分别分解再拼凑起来，可以近似实现VMD的扩展功能，但是通过该方法得到的模函数无法保证频率成分的单一，这是因为实部和虚部分解所得模函数的频段并非一一对应，强行拼接会导致模函数失真。为了获得同实信号相同的分解效果，如图5所示，现将待分解的复数信号$x(t)$按频谱的正负划分成两部分：

图  5  扩展VMD原理示意图

Figure  5.  Schematic diagram of the extended VMD

下载: 全尺寸图片 幻灯片

其中，$\Re (\cdot)$表示取实部，${X_{\mathrm{c}}}(\omega )$表示$x(t)$的傅里叶变换，${\text{F}}^{-1}(\cdot)$表示傅里叶逆变换，${(\cdot)}^{*}$表示取共轭，$u(\omega )$为阶跃函数：

${x_ + }(t)$和${x_ - }(t)$都是实信号且各自包含了原信号的正负频谱信息，可以分别对${x_ + }(t)$和${x_ - }(t)$VMD分解，有

对各模函数做希尔伯特变换搬移频谱可得到复信号$x(t)$的VMD表达式：

其中，$\hslash (\cdot)$表示希尔伯特变换，该方法可以将$x(t)$分解成$2N$个频率成分单一的模函数。理论上EMD得到的分解结果完全由信号本身决定，易于受模态混叠的影响，而VMD可以通过优化惩罚因子${\alpha _{\mathrm{p}}}$和分解层数N得到更精细、准确的分解结果，因此对VMD参数的优化变得尤为关键。

此外，从时频谱上看，主体回波表现为一条平行于时间轴的亮线，能量较为集中，而微动散射点的能量分散，分布跨越多个频率单元。因此对单个模函数而言，主体部分对应的能量比微动部分高，通过将模函数按能量高低排序再设置能量阈值T_e可以实现微多普勒分量的抑制。

设待优化变量集合${\boldsymbol{\rho}} = \left\{ {{\alpha _{\mathrm{p}}},N,T_{\mathrm{e}}} \right\}$，用本节所提方法对式(10)沿方位向分解，将离散化后的${S_{\mathrm{r}}}({f_{\mathrm{t}}},\eta )$记为${\boldsymbol{S}} \in {C^{{N_{\mathrm{r}}} \times {N_{\mathrm{s}}}}}$，其分解后按模函数能量排序所得张量为${{\boldsymbol{S}}_{{\mathrm{vmd}}}}({\alpha _{\mathrm{p}}},N) \in {C^{{N_{\mathrm{r}}} \times {N_{\mathrm{s}}} \times 2N}}$，假设能量阈值T_e对应的选择向量为

则经过微多普勒抑制后的回波矩阵为

对${{\boldsymbol{S}}_{\rm{opt} }}$方位向做FFT得到成像结果${{\boldsymbol{S}}_{\mathrm{I}}}$，这里用图像熵作为评价图像质量的指标，熵值越小表示图像质量越高。设二维图像矩阵为${\boldsymbol{I}}(i,j)$，图像熵的计算表达式如下：

代价函数$f({\boldsymbol{\rho}} ) = {{E}}\left( {{S_{\mathrm{I}}}({\boldsymbol{\rho}} )} \right)$，所以问题转化为

这里采用差分进化算法对式(26)进行求解。差分进化算法是一种适用于求解多维空间全局最优性的算法，其不依赖于梯度信息，不易陷入局部最优值^[28]。

基于差分进化算法的回波VMD分解和成像模态选择过程如下所示。

步骤1　初始化参数，设解空间种群数量为${\text{NP}}$，第G次迭代后的解空间表示如下：

步骤2　从当前种群中随机选择3个个体${X_{{\mathrm{r}}1,G}}, {X_{{\mathrm{r}}2,G}},{X_{{\mathrm{r}}3,G}}$进行差分操作，设${F_{\mathrm{v}}}$为变异因子，则得到变异个体：

步骤3　种群中每个个体和变异个体依概率${\text{CR}}$进行交叉，得到试验个体${{\boldsymbol{U}}_{i,G}} = [u_{i,G}^{{\alpha _{\mathrm{p}}}} {\text{ }}u_{i,G}^N {\text{ }}u_{i,G}^{T_{\mathrm{e}}}]$，其中：

${j_{{\text{rand}}}}$为随机的一个分量，确保交叉后的试验个体至少有一维分量由变异个体提供。

步骤4　将当前个体和试验个体分别代入代价函数计算适应值，选择更优的作为下一代：

步骤5　直到迭代次数G达到${G_{\max }}$，或种群最优解达到预定误差精度时算法结束。

综合以上分析，得到基于VMD和模态优选的ISAR微多普勒抑制算法流程，见图6。该算法可以实现对微多普勒分量的精准抑制，极大程度改善微多普勒对成像的影响。

图  6  所提微多普勒抑制方法流程图

Figure  6.  Flow diagram of the proposed micro-Doppler suppression method

下载: 全尺寸图片 幻灯片

本文算法的核心在于扩展VMD将复数回波分解为多个平稳模函数，并通过进一步优化VMD参数和筛选微多普勒的能量阈值保证主体回波的完整保留和微多普勒的彻底抑制。从原理上讲，本文所提微多普勒抑制方法可以和各种ISAR成像算法结合，适用于各种微动目标场景。

4.   仿真与实测数据验证

本节将通过仿真和实测数据验证上述微多普勒抑制方法的性能，首先通过多点微动目标成像实验，将所提方法与直接成像以及基于PGA^[29], EMD等方法的结果作比较，初步验证了所提方法对微多普勒的抑制效果及噪声环境下的稳定性；然后利用实测数据再次论证了本文所提算法在成像质量上的优势。

4.1   仿真实验

首先，仿照旋翼无人机的结构和姿态进行多点目标ISAR仿真，仿真参数如表1所示，点目标布局如图7(a)所示，图中圈出的是旋翼部件，其余为无人机机身。直接ISAR成像结果如图7(b)，旋翼周围出现了明显的模糊，严重影响图像质量。

表  1  仿真参数表

Table  1.  Simulation parameter table

参数	数值
载波频率${f_{\mathrm{c}}}$	35 GHz
信号带宽${B_{\mathrm{r}}}$	10 GHz
脉冲重复频率${\text{PRF}}$	833 Hz
成像积累时间	2.46 s
无人机运动速度	5 m/s
目标俯仰角	15°
旋翼个数	4
单旋翼叶片个数	2
叶片长度l	0.1 m
单叶片微动散射点个数	5
旋翼转速${\omega _{\mathrm{r}}}$	20$\pi$ rad/s

下载: 导出CSV 

| 显示表格

图  7  多点目标布局及直接成像图

Figure  7.  Multi-point target layout and direct imaging diagram

下载: 全尺寸图片 幻灯片

分别应用PGA^[29], EMD^[20], LMD^[23]以及本文所提方法对图7(b)所示微多普勒干扰进行抑制处理，其中，PGA是一种常用的ISAR图像自聚焦方法，对各种原因导致的图像散焦有着一定的改善作用。各种算法所得成像结果如图8所示，相比于PGA算法，模态分解类算法整体上对微多普勒的抑制更为有效，而在模态分解类算法中，本文所提算法对旋翼和主体之间的低频微多普勒抑制更为彻底，同时主体部分也得到了更为完整的保留。

图  8  各种算法处理后的成像结果对比

Figure  8.  Comparison of the imaging results processed by various algorithms

下载: 全尺寸图片 幻灯片

接下来在时频域观察旋翼的微多普勒抑制效果，在仿真的超高分辨条件下微动会超过一个距离单元，导致单个距离门内的微多普勒信号并不完整，这里将微动所在区域的多个距离门累加起来观察和分析微多普勒效应。实验结果如图9所示，区域1为单个旋翼，此时EMD, LMD及本文方法都较好地抑制了微多普勒并完整保留了主体，而PGA算法未能对微多普勒做出正确补偿。

图  9  各种算法处理后的STFT谱对比

Figure  9.  Comparison of the STFT spectrogram processed by various algorithms

下载: 全尺寸图片 幻灯片

区域2所在距离单元中包含了两个旋翼及多个机身散射点，从几种方法的对比可以看出在多个主体散射点与多个微动点混合的复杂场景下，PGA的效果依然不佳，表明了PGA算法的局限性以及设计微多普勒专业处理算法的必要性。在剩余对比结果中，本文方法比传统基于EMD和LMD的方法对回波分解和微多普勒的筛选更为准确，对低频微多普勒有着更好的压制效果。值得一提的是，由于该仿真条件下微多普勒带宽超过了脉冲重复频率的一半，所以其在时频图中不再是完整的正弦曲线，但是本文所提方法依然有效。

为了对微多普勒抑制结果进行定量分析，定义能量相似比指标${P_{\mathrm{c}}}$：

其中，$\displaystyle\sum\nolimits_m {\displaystyle\sum\nolimits_n {\left| {{S_{\sup }}(m,n)} \right|} }$表示微多普勒抑制后旋翼所在距离门的能量，$\displaystyle\sum\nolimits_m {\displaystyle\sum\nolimits_n {\left| {{S_{{\text{ideal}}}}(m,n)} \right|} }$表示无微多普勒理想情况下旋翼所在距离门的能量。显然${P_{\mathrm{c}}}$越小说明对微多普勒分量的抑制和主体保留的效果越好。

无人机两个区域中旋翼的能量相似比结果如表2所示，除此之外，计算整幅成像结果的图像熵、对比度和锐度指标，各种量化指标均证明了本文方法在仿真中的有效性。

表  2  各种算法处理后的图像质量比较

Table  2.  Comparison of image quality processed by various algorithms

算法	能量相似比(区域1)	能量相似比(区域2)	图像熵	对比度	锐度
直接成像	8.64	3.87	8.84	18.65	1.03E+10
PGA	9.60	4.11	8.75	20.19	9.54E+09
基于EMD	0.22	1.23	8.31	21.44	8.72E+09
基于LMD	0.16	0.84	8.28	20.45	8.58E+09
本文方法	0.15	0.16	8.16	22.71	9.15E+09

下载: 导出CSV 

| 显示表格

接下来，为了比较噪声环境下的各种算法的稳定性，在上述仿真回波中添加了不同信噪比的高斯噪声，所得成像结果如图10所示。

图  10  不同信噪比下的仿真成像结果对比

Figure  10.  Comparison of simulated imaging results at different SNR

下载: 全尺寸图片 幻灯片

从图10可以直观看出，本文所提算法具有较强的鲁棒性，在各个信噪比下对噪声的抗性均优于其他对比算法。这一优势来源于基于差分进化的模态筛选机制，其在理论上对背景噪声同样有一定的抑制效果，因此相比于PGA, EMD等传统方法，本文所提方法抗噪声能力更强，保证了成像质量。

4.2   实测实验

为了进一步验证算法的优越性，对真实六旋翼无人机开展成像研究，回波数据来自空军预警学院的微波光子(MicroWave Photonic, MWP) 超宽带雷达样机^[9]，其距离分辨率最高可达1.68 cm。图11给出了大致成像场景和目标形状示意图，雷达详细参数见表3。

图  11  成像场景及目标示意图

Figure  11.  Radar system and target diagram

下载: 全尺寸图片 幻灯片

表  3  六旋翼无人机实测实验参数表

Table  3.  Experimental parameter table for the hexacopter UAV

参数	数值
载波频率${f_{\mathrm{c}}}$	35 GHz
信号带宽${B_{\mathrm{r}}}$	10 GHz
脉冲重复频率${\text{PRF}}$	833 Hz
成像积累时间	0.72 s
目标俯仰角	15°

下载: 导出CSV 

| 显示表格

对无人机直接成像以及进行各种微多普勒抑制处理后的成像结果如图12所示。为了证明实虚拼接的CVMD^[27]方法对微多普勒抑制的不适用性，本次实验增加了在相同模态筛选机制下与CVMD的对比结果。可以看出，在直接成像图中，旋翼所在距离单元均出现了沿方位向的散焦条带；而经过PGA处理后，散焦条带有一定改善，但仍有明显残留，同时估计相位受到微多普勒的干扰导致机身成像质量下降；基于EMD和LMD的方法在高频处实现了对微动的完全抑制，在低频处由于分解精度不足以及筛选机制的欠缺无法分辨主体附近的微多普勒，导致机身周围的微动干扰依然存在；基于CVMD的方法在本文所提模态筛选机制下对微多普勒起到了一定的抑制效果，但由于其分解结果的频率成分不具备唯一性，导致成像模态部分失真。相比之下，本文所提方法不仅在高频抑制了微多普勒，在低频处凭借VMD的高分解精度和模态筛选上的图像熵迭代策略，极大程度抑制了旋翼对主体以及旋翼与旋翼之间的微动干扰，消除了微多普勒效应带来的图像模糊，显著提升了无人机成像质量。

图  12  各种算法处理无人机数据的成像结果对比

Figure  12.  Comparison of the effects of different algorithms on the processing of measured data

下载: 全尺寸图片 幻灯片

对图12第2行红色圈中旋翼中心点做STFT分析，结果如图13所示，与其他方法相比，所提方法抑制了该距离单元内高频微多普勒分量，在低频处凭借分解精度和筛选机制上的优势，较为完整地保留了旋翼中心对应的单频分量，实现了对全频带微多普勒分量的精准抑制。

图  13  处理后旋翼中心点对应方位STFT谱

Figure  13.  STFT spectrogram of the center point on the rotor after processing

下载: 全尺寸图片 幻灯片

为了进一步衡量微多普勒抑制后的聚焦效果，在图14(a)中，对旋翼中心点做残余相位分析，这里残余相位是指实际相位与线性相位的差值。可以看出本文方法基本抑制了微动带来的残余相位，图14(b)展示了该点的方位剖面对比结果，可以直观地说明采用本文方法处理后的回波成像结果具有更高的聚焦度。

图  14  旋翼中心点残余相位与方位剖面对比图

Figure  14.  Azimuth residual phase and profile comparison of the center point on the rotor

下载: 全尺寸图片 幻灯片

表4采用图像熵、对比度及锐度指标对图12中各幅图像的质量进行定量衡量，其中采用本文方法获得的ISAR图像综合质量最优，由此可见，本文方法在对实测数据进行处理时，仍能很好地抑制微多普勒对成像的干扰。

表  4  六旋翼无人机成像质量比较

Table  4.  Comparison of imaging quality of the hexacopter

算法	图像熵	对比度	锐度
直接成像	8.19	22.12	5.44E+11
PGA	7.98	23.59	2.90E+11
基于EMD	7.63	23.28	4.17E+11
基于LMD	7.59	23.14	4.11E+11
基于CVMD	7.67	23.66	4.05E+11
本文方法	7.05	25.79	4.33E+11

下载: 导出CSV 

| 显示表格

最后，本文对所提方法的计算量进行了分析，设平动和转动补偿后的回波矩阵为${S_{M \times N}}$，M和N分别为方位向和距离向采样点数，所提方法在迭代过程中VMD分解维数的平均值为K，不同算法的计算复杂度见表5。由于本文算法处理精度更高，在运算量上略大于其他算法，但可以更准确地实现超高分辨ISAR微多普勒抑制，获得更高质量成像结果。

表  5  不同算法的运算复杂度

Table  5.  Computational complexity of different algorithms

算法	运算复杂度
本文算法	$O(MNK{L_1}{L_2}{\log _2}M)$，${L_1}$为VMD循环次数，${L_2}$为差分进化迭代次数
基于LMD[23]	$O(MN{L_3}{L_4})$，${L_3}$为调频信号迭代次数，${L_4}$为子信号分离次数
基于EMD[20]	$O(MN{L_5}{L_6})$，${L_5}$为EMD循环次数，${L_6}$为复信号投影维度数
PGA[29]	$O(M{L_7}{\log _2}M + MN{L_7})$，${L_7}$为PGA迭代次数

下载: 导出CSV 

| 显示表格

5.   结语

目标部件微动形成的微多普勒信号具有大带宽、时变非平稳等特点，导致其被抑制非常困难，显著影响超高分辨ISAR的目标成像质量。本文利用微多普勒信号与主体回波的时频分布差异，提出一种基于扩展VMD的信号分解算法实现主体与微多普勒的分离，利用差分进化算法以ISAR图像熵为代价函数迭代优化分解参数和能量阈值，得到最优成像模态，在完整保留主体回波的同时抑制了微多普勒对成像的影响。所提方法凭借更高的回波分解精度和最小熵模态筛选机制，相比传统基于EMD和LMD的方法，能够在不影响主体成像的同时精确抑制微多普勒干扰。最后，通过一系列的仿真和真实超大带宽微波光子ISAR无人机成像实验验证了本文方法的有效性，并与PGA自聚焦、基于EMD, LMD等方法进行了残余相位、剖面以及图像质量指标对比，结果表明，所提方法最为准确地抑制了无人机旋翼周围的微动干扰，相比其他方法更加完整地保留了无人机主体机身，获得更高质量的ISAR图像，充分体现本文方法的优越性。