1.   引言

逆合成孔径雷达(Inverse Synthetic Aperture Radar, ISAR)可以全天候、全天时、远距离获取非合作目标的高分辨图像，在军事和民用上都具有非常重要的应用价值^{[1, 2]}。高效的距离多普勒(Range-Doppler, RD)算法或时频分析法可以得到非合作目标的高分辨ISAR图像；然而，所得的图像位于距离多普勒域，不能反映出目标的真实尺寸。为了后续更好地对目标进行识别，需要对目标ISAR图像进行方位定标(cross-range scaling)；而估计转动角速度(Rotation Angle Velocity, RAV)是方位定标的一个难点和研究热点。

针对该问题，近年来国内外学者先后提出多种ISAR方位定标算法^[3—11]。文献[3]提出一种Radon-Ambiguity方法，通过在不同距离单元的平面上检测线性调频信号的调制率，进而估计出RAV；但在计算过程中需要2维搜索，计算量很大。文献[4－7]方法通过选择不同距离单元内的特显点，提取出其相位系数来估计RAV；然而在实际ISAR成像中很难找到理想的散射点，并且特显点的信噪比会影响RAV的估计精度，同时在不同距离单元估计出的RAV值起伏较大。文献[8]在转动速度、转动加速度的2维平面上以图像信息熵或对比度为准则进行搜索，当ISAR图像质量最好时，对应参数即为转动速度与转动加速度的估计值；然而，该方法需要进行多维搜索，无法兼顾算法的精度和运算效率。文献[9]提出一种基于压缩感知的方法，利用目标在成像空间内具有稀疏特性，实现目标高分辨成像和方位定标；但存在对噪声敏感且计算复杂度高等问题。文献[10]利用两幅RD图像之间的旋转相关特性估计RAV，此方法是在假设旋转中心已知的情况进行的，然而在实际图像处理中旋转中心很难获取，限制了此方法的实用性。文献[11]利用2维快速傅里叶变换(2-D FFT)以及极坐标映射估计RAV，该方法将ISAR图像的旋转变为极坐标域中沿极角的平移，进而估计RAV。虽然避免了旋转中心的估计，但是由2维频域向极坐标转换过程中的插值操作会带来RAV估计精度损失和计算复杂度高的问题。因此，需要进一步研究高效的ISAR图像方位定标方法。

本文提出一种基于伪逆极坐标快速傅里叶变换(Pseudo Polar Fast Fourier Transform, PPFFT)的高效ISAR方位定标方法。首先利用高效的PPFFT把两幅不同时刻ISAR图像的旋转转化为伪极坐标域中沿极角方向的平移；然后定义了一种新的代价函数-积分相关函数，来粗估目标的RAV；最后，采用二分法快速得到最优的RAV，完成方位定标。相比于现有ISAR定标算法，本文算法避免了插值操作带来的精度损失和高计算复杂度问题。

2.   ISAR图像RD域与频域关系建模及分析

图1给出了经过平动补偿后ISAR成像的几何模型，其中目标以均匀角速度${\omega _0}$旋转，雷达与目标旋转中心之间的距离为${r_a}$；则对旋转目标上任一散射点$P\left( {{r_0},{\theta _0}} \right)$在${t_k}$时刻与雷达之间的瞬时斜距可写为

图 1  ISAR成像的几何模型示意图

下载: 全尺寸图片 幻灯片

由于雷达到目标的距离远大于目标尺寸，在远场条件下照射到目标的电磁波可用平面波近似；因此，斜距表达式(1)可近似写为

其中，${x_0} = {r_0}\cos {\theta _0}\;;\;{y_0} = {r_0}\sin {\theta _0}$，根据距离时间延迟和多普勒频率变化特性，则该散射中心$P$在距离多普勒域上的坐标可表示为^[10,11]

其中，${X_0} = - {{{r_a}} / {{\eta _r}}}$为常数项，${f_s}$是距离向采样频率，${\rm PRF}$是脉冲重复频率，$\rm c$是光速，$\lambda$是波长，$M$是累积脉冲数，${\eta _r} = {c / {\left( {2{f_s}} \right)}}$和${\eta _a} = {{\lambda {\rm PRF}} / {\left( {2M{\omega _0}} \right)}}$分别是距离向尺度因子和方位向尺度因子。

把式(3)和式(4)表示为矩阵形式，可得散射点$P$在不同时刻${t_{k1}}$和${t_{k2}}$得到的两幅RD图像${I_1}$和${I_2}$中的位置具有式(5)关系^[10,11]：

式中，${\theta _d} = {\omega _0}\left( {{t_{k2}} - {t_{k1}}} \right) = {\omega _0}\Delta {t_k}$是两幅RD图像之间的旋转角度，${C_r} = {{{\eta _a}} / {{\eta _r}}}$为相对尺度因子。从式(5)可知，第2幅RD图像${I_2}$是由第1幅RD图像${I_1}$经过旋转-拉伸变换得到的。下面就不同时刻得到的两幅RD图像在频域的关系进行分析；对给定的$N \times N$大小的RD图像$I\left( {u,v} \right)$，其2维傅里叶变换$\hat I\left( {{\omega _x},{\omega _y}} \right)$可表示为

根据傅里叶变换性质可知，对频率变量$\left( {{\omega _x},{\omega _y}} \right)$在不同坐标下，可以得到不同的2维傅里叶变换形式；如2维频率变量在卡迪尔坐标采样网格下，$\left( {{\omega _x},{\omega _y}} \right) = \left( {k,l} \right)$, $k,l = ( - {M / 2},{M / 2} - 1)$, $M$为采样点数，其2维傅里叶变换频域形式可写为

(7)

从式(7)可知，由于在时域和频域都是均匀采样，因此可用2维快速傅里叶变换(2-D FFT)实现。如果2维频率变量$\left( {{\omega _x},{\omega _y}} \right)$是在极坐标网格下采样，则在极坐标网格中：${\omega _x} = {r_k}\cos {\theta _l},{\omega _y} = {r_k}\sin {\theta _l}$, ${r_k} = k,{\theta _k} = {{2{{π}} l} / L}$, $k,l = ( - {M / 2},{M / 2} - 1)$；其中，$M$和$L$分别表示沿极径${r_k}$和极角${\theta _l}$方向的采样点数，在极坐标网格下其2维傅里叶变换形式为

(8)

根据式(5)和式(8)，对不同时刻得到的两幅RD图像${I_1}$和${I_2}$，在极坐标采样网格下做2维傅里叶变换，其在极坐标域的关系式为

从式(9)可知，两幅ISAR图像${I_1}$和${I_2}$在RD域的旋转和拉伸分别对应于极坐标中极角的平移和极径的伸缩，然而当2维频率变量$\left( {{\omega _x},{\omega _y}} \right)$处于极坐标网格时，其2维傅里叶变换没有相应的快速变换，计算复杂度急剧增加，给数据的实时处理带来困难。因此需要进一步研究高效的ISAR图像定标方法。

3.   本文提出的ISAR方位定标方法

根据第2节分析可知，不同时刻得到的两幅ISAR图像在RD域旋转对应为极坐标傅里叶变换域沿极角的平移；在此基础上，本节提出一种基于PPFFT的快速方位定标方法。

3.1   PPFFT原理

伪逆极坐标傅里叶变换(PPFT)是一种在伪逆极坐标采样网络下的2维傅里叶变换，其2维频率变量在伪逆极坐标网格下采样形式为^{[12, 13]}

式中，${P_1}$和${P_2}$满足的条件为：${P_1} \triangleq \left\{ \left( { - \displaystyle\frac{{2l}}{N}k,k} \right)\biggr|\right.$$\left.- \displaystyle\frac{N}{2} \le l \le \displaystyle\frac{N}{2}, - N \le k \le N \right\}$, ${P_2} \triangleq \left\{ \left( {k, - \displaystyle\frac{{2l}}{N}k} \right)\biggr|\right.$$\left. - \displaystyle\frac{N}{2} \le l \le \displaystyle\frac{N}{2}, - N \le k \le N \right\}$，其中$k$表示“伪极径”，$l$表示“伪角”。图2(a)和图2(b)分别表示式(10)中的伪逆极坐标网格下采样区域${P_1}$和${P_2}$，完整的伪逆极坐标采样网格如图2(c)所示。根据对应的关系伪逆极坐标采样网格用$\left( {r,\theta } \right)$表示可写为

图 2  伪逆极坐标采样网格

下载: 全尺寸图片 幻灯片

其中，$r_k^1,r_k^2,\theta _l^1,\theta _l^2$满足的条件为：$r_k^1 \!=\! k\sqrt {4{{\left(\! {\displaystyle\frac{l}{N}} \!\right)}^2} \!\!\!+\! 1} ,$$r_k^2 = k\sqrt {4{{\left( {\displaystyle\frac{l}{N}} \right)}^2} + 1}$, $\theta _l^1 = {{{π}} / 2} - \arctan \left( {\displaystyle\frac{{2l}}{N}} \right),$$\theta _l^2 = \arctan \left( {\displaystyle\frac{{2l}}{N}} \right)$。

类似于式(7)和式(8)定义的迪卡尔坐标和极坐标下2维傅里叶变换，ISAR图像在式(10)和图2表示的伪逆极坐标采样网格中，其2维傅里叶变换形式分别为^[12,13]

(12a)

幸运的是伪逆极坐标采样网格下的2维傅里叶变换能够快速实现。图3描述了式(12a)快速实现的示意图；首先对原始图像沿y 轴方向补零，接着做快速傅里叶变换(FFT)；然后，沿x 轴方向做Chirp-Z变换^[13]，得到了2维频率在伪逆极坐标采样网格下对应的区域${P_1}$；同理，可以得到对应的区域${P_2}$。将${P_1}$和${P_2}$叠加就可得到整个图像的在伪逆极坐标采样网格下2维傅里叶变换结果。实现过程仅需进行2次1维FFT和2次1维Chirp-Z变换^[12,13]，避免了2维插值操作。因此PPFFT具有高的精度和较低的计算复杂度。

图 3  伪逆极坐标采样区域${P_1}$的2维傅里叶变换快速实现示意图

下载: 全尺寸图片 幻灯片

3.2   基于PPFFT的RAV快速估计

根据3.1节的PPFFT的原理，首先对不同时刻${t_{k1}}$和${t_{k2}}$得到的两幅RD图像${I_1}$和${I_2}$进行PPFFT，用${P_{k1}}$和${P_{k2}}$分别表示两幅RD图像PPFFT后的幅度信息，则其具有式(13)关系：

其中，$\left( {{\omega _{x0}},{\omega _{y0}}} \right)$为旋转中心，其位于2维频域零频处；式(13)表明两幅RD图像旋转对应伪逆极坐标系下沿伪极角$\theta$的平移。因此，只要估计出${t_{k2}} - {t_{k1}}$时间段内两幅图像的旋转角度$\Delta \theta$，就可以得到目标的RAV。为了估计两幅ISAR图像的旋转角度$\Delta \theta$；本文定义了一种新的1维代价函数-径向积分相关函数，其表达式为

其中，${\varphi _m}\left( {{\theta _i}} \right) = \displaystyle\sum\nolimits_{j = 1}^N {{P_{k,m}}\left( {{r_j},{\theta _i}} \right)}, \;\;m = 1,2$为两幅RD图像经PPFFT后沿径向进行积分的结果；根据式(14)，两幅RD图像的旋转角度$\Delta \theta$可通过式(14)取最大值时的位置$k$获得。为了简化计算，利用傅里叶变换的卷积特性，最大值对应的位置$k$可由式(15)获得。

式中，${F_\theta }\left( \cdot \right)$表示傅里叶变换，${F_\theta }^{ - 1}\left( \cdot \right)$表示傅里叶逆变换，下标$\theta$表示沿极角方向。根据获得的两幅RD图像的旋转角度$\Delta \hat \theta$，就可得到RAV的估计值$\hat \omega = {{\Delta \theta } / {\left( {{t_2} - {t_1}} \right)}}$以及相对尺度因子${C_{r0}}$。

然而从式(5)可知，两幅图像之间除了旋转以外还存在尺度伸缩；尺度伸缩会对RAV的估计精度产生影响，所以为了更精确地估计RAV，需要粗估RAV对原始图像做旋转-拉伸补偿，从而精确估计出旋转角度。这里使用二分法求解最优的RAV；二分法是一种能够重新分配搜索区间并选择根所在的子区间进行2次搜索，进一步寻找真实根的经典优化算法^[14]。本文利用二分法来得到最优的旋转角度，选择合适的搜索区间$\left[ {a,b} \right]$，在小区间内实现迭代搜索，在完成多次迭代搜索就可得到目标精确RAV，从而完成ISAR图像方位定标。整个算法的实现过程如图4所示。

图 4  提出的ISAR方位定标算法实现流程图

下载: 全尺寸图片 幻灯片

3.3   计算复杂度分析

本小节分析提出方法的计算复杂度，并与文献[11]提出的基于2维快速傅里叶变换和极坐标映射估计RAV的方法进行比较。假设距离维采样点数为${N_r}$，方位维采样点数为${N_a}$；且长度为$N$的信号做1次复数相乘需要的计算量为$O\left( N \right)$，做1次FFT计算复杂度为$O\left( {N{{\log }_2}N} \right)$。

文献[11]提出的估计RAV方法主要有以下几步构成：首先，对不同时刻得到的两幅RD图像进行2维FFT操作变到2维频域；然后，利用插值操作从2维频域映射到极坐标域；最后，在极坐标域定义2维相关函数来粗估RAV，并利用优化算法来迭代求解最优RAV。其总的计算复杂度大约为

其中，${N_{\ker }}$为插值核长度，$Q$为算法优化过程需要的迭代次数。

根据图4提出的算法实现流程图，提出的ISAR方位定标算法主要有以下几步构成；首先对不同时刻得到的两幅RD图像进行PPFFT变换到伪逆极坐标域，$N$点Chirp-Z变换需要的计算复杂度为$O\left( {N + 2N{{\log }_2}N} \right)$^[15]；然后，在伪逆极坐标域定义1维积分相关函数来粗估RAV，并利用二分法来迭代求解最优RAV。提出的算法总的计算复杂度大约为

其中，$M$为二分法求解过程所需要的迭代次数。

根据式(16)和式(17)可知，文献[11]的方法由于需要使用插值操作和2维代价函数，导致其计算复杂度偏高。而本文提出的方位定标算法由于避免了插值操作，且定义了新的1维代价函数-积分相关函数；因此，其计算复杂度远低于文献[11]的方法。

4.   仿真结果分析及实测数据验证

4.1   计算机仿真结果及其分析

为了验证本文算法的有效性，现把所提方法应用到飞机模型的ISAR方位定标中。雷达工作参数和飞机目标运动参数见表1。图5为经过平动补偿后在${t_{k1}}$和${t_{k2}}$两个不同时刻进行RD成像所得到的ISAR图像，很明显第2幅RD图像是第1幅RD图像的旋转-拉伸结果。

表 1  雷达参数和目标运动模型

参数	数值
载波频率	5.6 GHz
波长	0.0536 m
传输信号带宽	400 MHz
距离采样频率	512 MHz
脉冲重复频率	150 Hz
有效回波脉冲	600
旋转角速度	0.0436 rad/s

下载: 导出CSV 

| 显示表格

图 5  两个不同时刻得到的RD图像

下载: 全尺寸图片 幻灯片

对${t_{k1}}$和${t_{k2}}$两个时刻的RD图像进行PPFFT，得到图像在伪逆极坐标下的变换结果，将图像的旋转转化为沿伪极角方向的平移；接着对伪逆极坐标域的图像沿径向积分；图6(a)为利用式(14)定义的代价函数得到的相关曲线，根据峰值检测可得旋转目标的RAV为${\omega _0} = 0.0455$ rad/s，最后根据式(5)对RD图像进行旋转-拉伸补偿并完成方位定标。

图 6  粗估的1维相关曲线

下载: 全尺寸图片 幻灯片

然而，由于受两幅RD图像之间拉伸的影响，一次粗估所得的RAV存在较大的误差，需要根据粗估得到的RAV补偿图像拉伸带来的影响，再进一步精细估计最优的RAV。采用二分法求最优RAV的迭代收敛曲线如图7所示，算法迭代了10次左右就趋于收敛，计算得到最优的RAV为${\omega _0} =$ 0.04358 rad/s。

图 7  迭代收敛曲线

下载: 全尺寸图片 幻灯片

为了更好地验证本文方法优势，本文和文献[11]方法进行比较。在本文仿真参数下，文献[11]方法得到最优的RAV为${\omega _0} = 0.04342$ rad/s，定标后ISAR图像如图8(a)所示。而利用本文方法得到最优的RAV为${\omega _0} = 0.04358$ rad/s，定标后ISAR图像如图8(b)所示。定标结果与目标真实尺寸基本一致；在RAV估计精度方面，文献[11]方法是0.41%，而本文方法是0.05%，说明本文方法能够对图像旋转角进行更为精确的估计，定标性能更好。

图 8  方位定标后的ISAR图像

下载: 全尺寸图片 幻灯片

为了更清晰直观比较计算复杂度，使用Intel Double-core处理器，CPU主频2.4 GHz，内存4 GB, Window 7操作系统，32位台式主机，在MATLAB 2014a环境下，分别运行了文献[11]方法和本文方法，运行时间如表2所示。由表2可知，本文提出的定标方法计算量远小于文献[11]方法。为了验证分析噪声环境下算法的稳健性，图9给出了本文方法和文献[11]方法在不同信噪比(SNR)下RAV估计的均方根误差(RMSE)；其中，RMSE获得是通过100次独立的门特卡罗实验。从图9可知，在SNR较高时，相比文献[11]方法，本文方法具有更好的估计精度。随着SNR的降低，本文方法和文献[11]方法RAV估计精度趋于一致，但其计算复杂度远低于文献[11]的方法，与上述的仿真结果和理论分析相吻合。

表 2  两种方法运行时间对比(s)

方法名称	粗估所用时间	定标总时间
文献[11]方法	107.224	958.659
本文方法	12.203	96.712

下载: 导出CSV 

| 显示表格

图 9  RAV估计的均方根误差随信噪比变化曲线

下载: 全尺寸图片 幻灯片

4.2   实测数据处理

Yak-42飞机模型如图10(a)所示，雷达载频${f_c}$为5.52 GHz，脉冲重复频率${\rm PRF}$为25 Hz，发射信号带宽${B_r}$为400 MHz，方位向脉冲回波个数为128。对Yak-42实测时所得的回波数据进行旋转角度估计，计算得到最优的RAV为${\omega _0} = 0.0092$ rad/s；根据估计所得出的RAV进行方位定标，得到定标后ISAR图像如图10(c)所示。同时，图10(b)给出了利用文献[11]方法方位定标后的ISAR图像。从图10可以看出，文献[11]方法定标得出的Yak-42的长度约为30.65 m，翼展约为31.86 m；本文定标得出的Yak-42的长度约为32.82 m，翼展约为31.24 m，更接近图10(a)中Yak-42飞机的实际尺寸；并且飞机轮廓清晰，几何尺寸明显，有利于后续对运动目标高效识别。

图 10  实测数据图像

下载: 全尺寸图片 幻灯片

5.   结论

本文提出一种基于PPFFT快速ISAR方位定标方法。该方法主要思想是通过PPFFT将RD图像的旋转运动映射为伪逆极坐标域中沿极角方向的平移运动，并通过定义了新的1维代价函数-积分相关函数和经典二分法估计得到精确的RAV。提出的方法避免了插值操作和2维代价函数，具有较低的计算复杂度，在工程实现中具有一定的优势。最后，计算机仿真和实测数据结果验证了本文方法的有效性。