1.   引言

当目标沿雷达视线(Light Of Sight, LOS)方向运动时，其回波信号的载频将发生偏移，即产生多普勒现象。除目标整体运动外，若目标或目标上的某些结构还存在独立的振动或旋转，则称其为微动。微动会在目标主体运动对应的主多普勒谱周围产生边带，即产生微多普勒效应^[1]。直升机、飞机旋转叶片、小型卫星和空间碎片等航空航天目标的典型微动包括自旋、进动和章动^[2]等。

对于空间微动目标，其高分辨雷达回波^[3,4]蕴含着散射中心2维或3维分布等结构信息，同时包含着自旋频率、进动频率及进动角等运动信息，上述信息为准确的目标分类、识别提供了重要支撑。目前，典型空间微动目标的高分辨雷达成像与微动参数估计方法研究^[5–9]已受到雷达成像与雷达自动目标识别领域的广泛关注。

空间微动目标的高分辨成像方法包括参数化方法^[10–14]与非参数化^[15–21]方法两类。其中，参数化成像方法首先建立各种微动形式的参数化模型，进而采用基于模型的参数估计方法实现高分辨成像。主要包括基于固定散射中心模型的成像方法^[10,22–24]及基于滑动散射中心模型^[11,13]的进动目标成像方法。对于章动等复杂微动形式，需要建立非常复杂的参数化模型，并实现大量未知参数的准确求解，由于目标的散射中心坐标与微动参数耦合，因此求解运算量很大。非参数化成像方法则主要包括自适应时频分析^[25]与散射中心航迹关联^{[16,17,26–28]}成像两类。与参数化成像方法相比，非参数化成像方法具有各种微动形式具有鲁棒性，能够避免由于模型失配而引起的较大误差，计算效率较高。

对于非参数化方法，基于航迹矩阵分解的成像方法^[17]可实现自旋、进动、章动等微动目标的高分辨成像。该类方法的关键步骤之一是在距离-慢时间域实现散射中心航迹的精确估计和关联。现有方法采用卡尔曼滤波器和最小欧氏距离准则，实现基于1维斜距信息的航迹关联^[17,29]，当散射中心回波包络交叉点较多、相距较近时容易产生较大的关联误差。此外，获取高质量的高分辨距离像(High Resolution Range Profile, HRRP)也是非参数化成像的关键步骤。

为了解决上述问题，本文提出一种基于距离-瞬时多普勒(Range-Instantaneous Doppler, RID)像序列的微动目标高分辨3维成像新方法。该方法充分利用散射中心在距离-瞬时多普勒域2维分布比距离-慢时间域1维分布可分性更强等特性，提出基于RID像序列的散射中心航迹关联方法，提高了航迹交叉点散射中心的可分性。进而通过带约束条件的矩阵分解求得散射中心3维分布和等效雷达视线矩阵，实现空间微动目标高分辨3维成像。最后，仿真数据证明了算法的有效性。

本文结构如下：第2节介绍了RID序列的生成方法；第3节研究了基于RID序列的航迹矩阵关联方法，以及基于现代谱估计的航迹矩阵精估计方法；第4节研究了基于航迹矩阵分解的微动目标高分辨3维成像方法；第5节以锥体章动目标为例，给出目标航迹关联及3维成像结果；最后一节进行了总结。

2.   距离-瞬时多普勒序列

对于信号 $s\left( t \right)$ ，其短时傅里叶变换(STFT, Short-Time Fourier Transform)满足^[15,30]：

其中， $\omega$ 表示角频率， $\tau$ 表示时延， $w\left( \cdot \right)$ 为窗函数。

为实现散射中心航迹的准确关联，需要获得其距离-瞬时多普勒像序列。假设雷达发射大时宽-带宽积脉冲信号，距离脉压后回波共包含 ${N_{\rm{r}}}$ 个距离单元，则对存在回波的距离单元 ${r\!_n}$ ( $n \in \left[ {{N_1},{N_2}} \right]$ , $\left[ {{N_1},{N_2}} \right]$ 为存在回波的距离单元区间)分别做STFT以得到其时频图 ${{{I}}\!_n}\left( {{f\!_{\rm{d}}},{t_{\rm{m}}}} \right)$ ，其中 ${f\!_{\rm{d}}}$ 表示多普勒， ${t_{\rm{m}}}$ 表示慢时间。随后，将时频图堆成3维矩阵 ${{Q}}({r\!_n},{f\!_{\rm{d}}},{t_{\rm{m}}})$ 。最后，沿时间轴 ${t_{\rm{m}}} = {t_i}$ ( $i \in \left[ {1,{N_{\rm{a}}}} \right]$ , ${N_{\rm{a}}}$ 为方位单元数)取出2维矩阵切片，即得到 ${t_i}$ 时刻的RID像。连续变换 ${t_{{i}}}$ 即可获得RID像序列。该过程示意图如图1所示。

图  1  RID像生成过程示意图

Figure  1.  The process of RID image series generation

下载: 全尺寸图片 幻灯片

3.   航迹矩阵关联及航迹矩阵精估计

3.1   基于分水岭法的RID像散射中心提取

为了利用RID像序列实现散射中心航迹关联，需要提取每幅图像中散射中心的2维坐标。分水岭(watershed)算法^[31]能够精确定位图像中的微弱边缘，并获得封闭且连续的分割曲线，因此适用于提取RID像中的散射中心支撑域。此外，相比于基于统计学的图像分割算法，该算法计算量小且分割较为准确，适用于图像数据的实时处理。因此，本文首先使用采用分水岭方法对RID像进行图像分割以获得每个散射中心对应的支撑域，然后计算每个支撑域对应的散射中心质心，并将此质心作为散射中心2维坐标的粗估计。基于watershed方法的RID像分割过程实现方法如下：

Step1：将原始图像归一化后，通过设定门限值将其转化为二值图像；

Step2：计算二值图像中每个像素点到其最近非零点的距离(如果像素本身非零，则其本身为最近的非零点，因此距离为0)，用于替代该像素点的像素值，得到矩阵 ${{D}}$ ；

Step3：令 ${{D}} = - {{D}}$ ，得到梯度图像；

Step4：采用分水岭方法对Step3中得到的梯度图像进行分割^[31]。

散射中心分割完成后，提取分割后每个散射中心的支撑域。具体步骤为：首先将分割后图像的1值和0值点赋为0，并将其他点赋为255；然后求二值图像的连通域；最后取其边界得到微动目标每个散射中心的支撑域。最后，将每个散射中心对应支持区的质心作为RID图像中每个散射中心2维坐标的估计。其中，质心计算方法如下：

其中， $f\left( {u,v} \right)$ 表示 $\left( {u,v} \right)$ 点处的像素值，u, v分别表示像素点的横坐标和纵坐标。

将每一时刻的RID像都做上述处理，则可获得散射中心在各个时刻对应的坐标。

3.2   航迹矩阵关联

由于微动目标具有惯性，因此认为相邻两幅距离-瞬时多普勒图像中同一散射中心的坐标连续变化，从而基于最近邻法实现航迹关联。设第 $i$ ( $i \in \left[ {1,{N_{\rm{a}}}} \right]$ )幅RID图像中的第 $j$ ( $j \in \left[ {1,{P}} \right]$ , ${P}\,$ 为散射中心个数)个散射中心的坐标向量为 ${{{a}}_{ij}}$ ，计算该点与第 $i + 1$ 幅RID图像中各个散射中心坐标向量的欧氏距离，选取与其欧氏距离最小的散射中心作为与该点关联的散射中心，即计算式(3)：

依次计算RID序列中相邻两幅图像中各散射中心的关联点，从而实现RID图像中各散射中心的关联，并得到微动目标航迹矩阵 ${{W}}$ 的粗估计。该矩阵的每一列对应一个散射中心在观测时间内的瞬时斜距。具体而言，基于RID序列的航迹关联实现方法如下：

Step1：初始化航迹矩阵 ${{{W}}\!_{{{{N}}_{\rm{a}}} \times {{P}}}}\,$ ，令所有元素都为0；

Step2：将 $P\,$ 个散射中心的初始时刻瞬时斜距写入 ${{{W}}\!_{{{{N}}_{\rm{a}}} \times {{P}}}}\,$ 的第1行；

Step3：令 $i = 1$ , $j = 1$ ，计算第 $i$ 幅RID像中第 $j$ 个散射中心与第 $i + 1$ 幅RID像中所有散射中心的欧氏距离，根据式(3)将最小欧氏距离对应的散射中心瞬时斜距写入 ${{{W}}\!_{{{{N}}_{\rm{a}}} \times {{P}}}}\left( {i + 1,j} \right)$ ；

Step4：令 $j = j + 1$ ，重复step3直到 $j = {P}$ ，实现第 $i$ 幅RID像与第 $i + 1$ 幅RID像的2维航迹关联；

Step5：令 $i = i + 1$ ，重复step3—step4直到 $i = {N_{\rm{a}}} - 1$ ，获得矩阵 ${{{W}}\!_{{{{N}}_{\rm{a}}} \times {{P}}}}\,$ 。

3.3   航迹矩阵精估计

由于RID像的距离分辨率为 ${\rho _r} = c/2B$ ，对于X波段雷达而言通常为10^–2 m量级，精度较低；而利用散射中心支撑域的质心对其2维坐标进行近似也会导致较大误差，从而使航迹矩阵产生抖动，影响微动目标3维成像的精度。针对该问题，在获得散射中心2维关联结果的基础上，可以进一步采用Root-MUSIC等谱估计方法^[32]对散射中心的瞬时斜距进行精估计，并对航迹矩阵进行修正，从而提高对微动目标散射中心3维坐标估计的准确性。

微动目标经运动补偿后的回波信号可表示为：

其中， $p \in \left[ {1,{P}} \right]$ 表示散射中心序号， ${A_{{p}}}$ 表示其幅度， $B$ 为带宽， $c$ 为光速， ${f\!_{\rm{c}}}$ 表示载频， ${R_{{p}}}$ 表示第 $p$ 个散射中心与参考点之间的瞬时斜距。若忽略距离窗，则式(4)可被改写为

其中， ${A'\!\!_{{p}}} = {A_{{p}}}\exp \Bigr( {{\rm j}4{{π}} \left( {{f\!_{\rm{c}}} - B/2} \right)\Delta {R_{{p}}}\left( {{t_{\rm{m}}}} \right)/c} \Bigr)$ ，散射中心对应的角频率为 ${\omega _{{p}}} = 4{{π}} \Delta f\Delta {R_{{p}}}\left( {{t_{\rm{m}}}} \right)/c$ , $\Delta f = B/{N_{\rm{r}}}$ , ${N_{\rm{r}}}$ 为距离单元数， $n \in \left[ {1,{N_{\rm{r}}}} \right]$ 。接下来，对每次距离向回波精确估计 ${\omega _{{p}}}$ 以求出 $\Delta {R_{{p}}}\left( {{t_{\rm{m}}}} \right)$ ，从而得到抑制旁瓣和噪声后的高质量高分辨1维距离像(HRRP)。

为了采用Root-MUSIC方法，首先构造距离回波的协方差矩阵：

其中， $m$ 表示窗长，且

通过Z变换找到与单位元距离最近的P个根可以求得角频率 ${\omega _{{p}}}$ 。随后，由 $\Delta {R_{{p}}}\left( {{t_{\rm{m}}}} \right) = {\omega _{{p}}}c/4{{π}} \Delta f$ 得到精估计的瞬时斜距 $\Delta {R_{{p}}}\left( {{t_{\rm{m}}}} \right)$ 。最后，通过最小欧氏距离准则将 ${t_{\rm{m}}}$ 时的瞬时斜距写入 ${{{W}}\!_{{{{N}}_{\rm{a}}} \times {{P\,}}}}$ 的相应行中，即可得到精估计的航迹矩阵 ${{{W}}'\!\!_{{{{N}}_{\rm{a}}} \times {{P}}}}$ 。此外，当回波的信噪比较低时，可以通过构造观测字典，采用噪声稳健的稀疏信号重构方法^[33,34]获得HRRP，并实现航迹矩阵的精估计。

4.   基于航迹矩阵分解的微动目标高分辨3维成像

根据运动的相对性，对于微动目标上的固定散射中心，其在距离-慢时间域的航迹矩阵可以表示为：

其中， $P$ 为散射中心个数， ${N_{\rm{a}}}$ 为方位单元数，矩阵 ${{R}}$ 表示不同时刻的等效雷达视线矩阵， ${{S}}$ 表示目标的散射中心坐标矩阵。根据式(8)可知，从 ${{W}}$ 中重构矩阵 ${{S}}$ 则可得到目标3维散射中心坐标。本文采用基于矩阵奇异值分解的方法重构矩阵 ${{S}}$ ^[17,21]。

利用矩阵奇异值分解法，航迹矩阵可以分解为 ${{{W}}\!_{1{{{N}}_{\rm{a}}} \times {{K}}}} = {{{U}}\!_{{{{N}}_{\rm{a}}} \times {{K}}}}{{{Σ}} _{{{K}} \times {{K}}}}{{V}}_{{{K}} \times {{K}}}^{\rm{T}}$ 。对于3维微动，根据矩阵秩的特性， ${{{Σ}} _{{{K}} \times {{K}}}}$ 的前3个奇异值较大，而其余奇异值趋近于零。因此可做如下近似：

其中，近似后得到 ${{R}}' = {{{U}}\!_1}$ , ${{S}}' = {{{Σ}} _1}{{{V}}\!_1}$ 。并且对于任意可逆矩阵 ${{{A}}_{3 \times 3}}$ , ${{R}}'{{S}}' = \left( {{{R}}'{{A}}} \right)\left( {{{{A}}^{ - 1}}{{S}}'} \right)$ 成立。

根据 ${{R}}$ 的定义可知 ${{R}}$ 各行构成的行向量的模为1。将 ${{R}}$ 用行向量的形式表示为 ${{R}} = {\left[ {{{{l}}_1}\;{{{l}}_2}\; ·\!·\!· \;{{{l}}_{{{{N}}_{\rm{a}}}}}} \right]^{\rm{T}}}$ ，则下列等式成立：

其中， ${{I}}$ 是单位矩阵。估计值 ${\hat {{A}\,}}$ 为式(10)的最小均方解，则 ${{\hat{{A}\,}}^{ - 1}}{{S}}'$ 相当于 ${{S}}$ 的等距映射。

对于任意的正交矩阵 ${{{A}}_1}$ ，满足下列关系：

其中， ${{{R}}'{\hat{{A}\,}}}$ 与 ${{{A}}_1}$ 相乘相当于旋转雷达视线， ${{\hat{{A}\,}}^{ - 1}}{{S}}'$ 与 ${{A}}_1^{\rm{T}}$ 相乘相当于散射中心关于原点旋转。由于满足 ${{{A}}_1}{{A}}_1^{\rm{T}} = {{I}}$ ，根据式(12)计算矩阵 ${{{A}}_1}$ ：

其中， ${{\hat{{l}}}_0}$ 是初始时刻雷达视线方向矢量，令 ${\hat{{A}\,}}{{\hat{{A}\,}}^{\rm{T}}} =$ $\left( {\begin{array}{*{20}{c}} {{{{l}}_{{1}}}}&{{{{l}}_{{2}}}} \\ {{{{l}}_{{2}}}}&{{{{l}}_{{3}}}} \end{array}} \right)$ ，则 ${\hat{{l}}} = \left[ {\begin{array}{*{20}{c}} {{{l}}_{{1}}} \\ {{{l}}_{{2}}} \\ {{{l}}_{{3}}} \\ \end{array}} \right]$ 。结合矩阵奇异值分解所得的 ${{R}}'$ , ${{S}}'$ 以及估计出的矩阵 ${\hat{{A}\,}}$ ，可以得到等效雷达视线矩阵为 ${{R}} = {{{R}}'{\hat{{A}\,}}}{{{A}}_1}$ ，散射中心3维坐标矩阵为 ${{S}} = {{A}}_1^{\rm{T}}{{\hat{{A}\,}}^{ - 1}}{{S}}'$ 。

通过上述航迹矩阵分解方法可以获得微动目标3维散射中心分布，进而实现空间微动目标高分辨3维成像。整体算法流程图如图2所示。

图  2  基于航迹矩阵分解的微动目标高分辨成像算法流程图

Figure  2.  The flow chart for high-resolution imaging of micro-motion targets based on trajectory matrix decomposition

下载: 全尺寸图片 幻灯片

5.   实验与分析

本节采用仿真数据对所提算法进行验证。微动目标散射中心分布如图3(a)所示，该目标由9个散射中心组成。仿真参数为：带宽2 GHz，载频10 GHz，脉冲重复频率PRF=2000 Hz，观测时间为1 s。章动目标自旋角频率为1 Hz，锥旋角频率为0.4 Hz，摆动角频率为0.1 Hz，摆动幅度为5°。回波信号的信噪比为20 dB。

距离脉压后的目标回波如图3(b)所示，其中最底部曲线对应锥顶散射中心。9个散射中心航迹交叉点较多，基于1维距离像关联难度较大。采用watershed方法从图3(c)所示RID像中提取散射中心支撑域的结果如图3(d)所示，进而从中计算出各散射中心坐标，如图3(e)所示，其中蓝色圆圈表示散射中心支撑域轮廓，红色标记表示通过计算得到的散射中心坐标。由图可知，散射中心轮廓清晰，分割效果良好。基于RID像序列的距离-多普勒-慢时间3维关联结果如图3(f)所示，在距离-时间维的关联结果如图3(g)所示，其中不同颜色代表不同散射中心的航迹。由该图可知，该方法能够有效避免交叉点处关联错误等问题，获得准确的散射中心航迹关联结果。

利用Root-MUSIC的谱估计方法对航迹矩阵进行精估计，结果如图3(h)所示。最后，采用航迹矩阵分解法获得微动目标3维散射中心分布的结果如图4(a)所示，其中红色星号表示估计值，蓝色圆圈表示真实值。可以看出，成像结果与真实散射中心分布一致，从而证明了本文所提算法的有效性。等效雷达视线矩阵估计结果如图4(b)所示。

图  4  章动目标3维成像结果

Figure  4.  3D image of the nutation target

下载: 全尺寸图片 幻灯片

下载: 全尺寸图片 幻灯片

为测试所提成像方法的抗噪性能，在保持其他参数不变的条件下，给目标回波中分别加入信噪比为0 dB, 5 dB, 10 dB, 15 dB, 20 dB的高斯白噪声。在每个信噪比下做50次蒙特卡洛实验，并按照式(13)计算均方根误差(Root Mean Square Error, RMSE)：

其中， $n \in \left[ {1,{N_{\rm{m}}}} \right]$ , ${N_{\rm{m}}}$ 为蒙特卡洛实验次数， $p \!\in\! \left[ {1,{P}}\, \right]$ , ${P}$ 为散射中心个数， $\left[ {T_{{{x}},p}^n,T_{{{y}},p}^n,T_{{{z}},p}^n} \right]$ 和 $\left[ {S_{{{x}},p}^n,S_{{{y}},p}^n,S_{{{z}},p}^n} \right]$ 分别表示第 $n$ 次蒙特卡洛实验中第 $p$ 个散射点的真实坐标和估计坐标。最终，不同信噪比下的RMSE曲线如图5所示，可以看出，RMSE随着SNR的增加而降低。

图  5  均方根误差随信噪比的变化曲线

Figure  5.  Variation of the RMSE with SNR

下载: 全尺寸图片 幻灯片

6.   结束语

针对传统参数化成像方法对复杂微动目标建模困难，未知参数求解运算量大等问题，本文提出一种基于RID图像序列的微动目标非参数化高分辨3维成像方法。该方法首先基于watershed法对RID图像进行分割提取散射中心，进而基于最近邻准则对散射中心航迹进行关联，接着通过Root-MUSIC方法实现航迹矩阵的精估计。最终，通过航迹矩阵分解实现微动目标的高分辨3维成像。该方法有效避免了参数化成像方法未知参数求解困难，易产生模型失配等不足。同时，2维关联方法克服了散射中心航迹交叉严重时，传统1维关联方法引起的关联误差，实现了复杂微动目标的高分辨3维成像。

在未来工作中，将研究低信噪比环境下的散射中心关联方法及非参数化微动目标高分辨3维成像方法，并进一步研究基于高分辨图像及等效雷达视线矩阵的微动目标特征提取及识别方法。