引言

逆合成孔径雷达(Inverse synthetic aperture radar，ISAR)可以对目标进行高分辨率二维成像^[1]。ISAR为了在距离向得到高的分辨率，普遍采用的方法是发射宽带线性调频信号，接收时采用脉冲压缩技术，距离向分辨率${\rho _{\rm{r}}} = \frac{c}{{2B}} $；在横向为了得到高的分辨率，主要依靠运动目标相对于雷达的转动得到，横向分辨率${\rho _{\rm{a}}}{\rm{ = }}\frac{\lambda }{{2\Delta \theta }}{\rm{ = }}\frac{\lambda }{{2\mathit{\Omega }T}} $。当ISAR系统参数已知时，距离向分辨率根据公式可以很容易计算。由于ISAR主要用于对非合作目标成像。目标的运动轨迹、成像积累时间内相对雷达的转动角度大小并不知道，所以ISAR的横向定标相对较难。根据横向分辨率公式可以看出，要确定横向分辨率就需要估计目标相对雷达的转角，因此ISAR图像的横向定标可以等效于成像积累时间内目标的转角估计问题。

目前ISAR目标转角估计的方法主要有：轨迹拟合法^[2]、多天线干涉法^[3]、层析成像法^[4]、互相关法^[5]、回波参数法^[6]及特征线斜率法^[7]。轨迹拟合法利用雷达窄带跟踪数据，对目标轨迹进行滤波、拟合等处理，从而计算目标相对于雷达的转角；但由于窄带信号跟踪数据精度较低，同时对于目标自身转动的估计无能为力，对复杂运动目标的估计误差更大。多天线干涉法利用多天线相干处理技术对两幅天线接收到的回波信号进行相位干涉，估计转动矢量大小，进而完成横向定标，但系统太复杂，成本较高。层析成像法利用卷积反投影方法获得在各种角度参数情况下的目标图像，然后根据峰值最大准则估计最优转角，但每次层析成像的计算量比较大。互相关法通过搜索使两幅ISAR图像的相关性最大的等效转动角速度来完成横向定标。该方法假设目标的等效旋转中心和几何中心一致，但在实际成像中这种假设并不总是正确。回波信号参数法是对一类方法的统称，通常在慢时间域对回波信号进行建模分析，认为回波信号是线性调频信号，或者具有二次相位项，并且其中的参数同转动角速度存在确定关系，通过数学工具估计这些参数，从而计算旋转角速度，进而完成横向定标。该方法定标效果受制于运动补偿效果和参数估计的精度，尤其对于等效转速度较低的目标，慢时间域线性调频信号调频率很低，较难估计。特征线斜率法利用目标俯视图的两条特征线即可完成定标；但该方法需要对图像进行二值化处理，然后进行Radon变换。在进行二值化处理时，不同阈值对于后期估计精度有很大影响，且阈值难以通过计算得到；同时在进行Radon变换时可能会产生模糊峰值。

分数阶傅里叶变换作为傅里叶变换的广义形式，是一种统一的时频变换，同时反映了信号在频域和时域的信息，由于其独有的特性，近年来在众多领域得到了广泛应用^[8-11]。本文基于回波信号参数法，给出了一种基于分数阶傅里叶变换(Fractional Fourier transform, FrFT)的ISAR横向定标方法，此方法假设目标在相干积累时间内运动平稳。在此前提下，本文首先对ISAR回波信号进行分析建模，确定回波信号在慢时间域为线性调频信号，然后利用FrFT估计出信号的调频率，求得相干积累时间内目标等效转动角速度，从而实现横向定标。

1 回波信号分析

相干积累时间运动平稳的目标，经过运动补偿以后，转化为转台运动。转台成像示意图如图 1所示。图 1描述的是一个匀速转动的三维目标投影到二维成像平面上。r₀为目标的等效旋转中心到雷达的距离，Ω为目标绕自身o点的等效旋转角速度。在起始时刻(即t=0)，目标上的任意一点A(x_a, y_a, z_a)到雷达的距离为

$ r = {\left[ {r_0^2 + r_{\rm{a}}^2 + 2{r_0}{r_{\rm{a}}}\sin \left( {\theta + \mathit{\Omega }t} \right)} \right]^{\frac{1}{2}}} $

(1)

式中:$ {r_{\rm{a}}} = \sqrt {x_{\rm{a}}^2 + y_{\rm{a}}^2} $。目标的几何尺寸要远小于目标到雷达的距离，那么式(1)可近似表示为

$ r \approx {r_0} + {x_{\rm{a}}}\sin \mathit{\Omega }t + {y_{\rm{a}}}\cos \mathit{\Omega }t $

(2)

目标上该点的多普勒频率可以表示为

$ {f_{\rm{d}}} = - \frac{2}{\lambda }\frac{r}{{{\rm{d}}t}} = - \frac{{2{x_{\rm{a}}}\mathit{\Omega }}}{\lambda }\cos \mathit{\Omega }t + \frac{{2{y_{\rm{a}}}\mathit{\Omega }}}{\lambda }\sin \mathit{\Omega }t $

(3)

式中:λ为雷达波长。由于相干积累时间很短，此时目标的旋转角度Δθ=Ωt很小，则式(3)可以近似成

$ {f_{\rm{d}}} = - \frac{{2{x_{\rm{a}}}\mathit{\Omega }}}{\lambda } + \frac{{2{y_{\rm{a}}}{\mathit{\Omega }^2}}}{\lambda }t $

(4)

由式(4)可知，各散射点回波信号在慢时间上表现为线性调频信号，目标的等效旋转角速度与该信号的调频率存在确定关系，其调频率K和多普勒初始频率f₀分别为

$ {f_{\rm{0}}} = - \frac{{2{x_{\rm{a}}}\mathit{\Omega }}}{\lambda } $

(5)

$ K = \frac{{2{y_{\rm{a}}}{\mathit{\Omega }^2}}}{\lambda } $

(6)

式(6)中，由于距离向定标很容易完成，则y_a可知，λ已知，所以根据回波信号的调频率K，就能计算出旋转角速度Ω。可见旋转角速度的估计问题转换为估计信号的调频率K。若雷达发射脉冲信号的脉冲周期用T表示，则发射时刻t_m=mT(m=0, 1, 2, …)称为慢时间。本文的研究基于慢时间域，后文中没有特别强调时间的情况下的时间均指慢时间。

2 分数阶傅里叶变换原理

通常，函数f(u)的p阶分数阶傅里叶变换可以表示为f_p(u) ^[12]。分数阶傅里叶变换的基本表达式为

$ {f_p}\left( u \right) = \int_{ - \infty }^{ - \infty } {{K_p}\left( {u,t} \right)f\left( t \right){\rm{d}}t} $

(7)

式中K_p(u, t)为分数阶Fourier变换的核函数，则

$ {K_p}\left( {u,t} \right) = \left\{ \begin{array}{l} {A_\alpha }\exp \left( {{\varphi _\alpha }} \right)\;\;\;\;\;\;\;\alpha \ne n{\rm{ \mathsf{ π} }}\\ \delta \left( {u - t} \right)\;\;\;\;\;\;\;\;\;\;\;\;\alpha = 2n{\rm{ \mathsf{ π} }}\\ \delta \left( {u + t} \right)\;\;\;\;\;\;\;\;\;\;\;\alpha = \left( {2n \pm 1} \right){\rm{ \mathsf{ π} }}\; \end{array} \right. $

(8)

式中:${A_\alpha } = \sqrt {1-{\rm{j}}\cot \alpha }, \alpha = \frac{{\rm{ \mathit{ π} }}}{2} $，p≠2n，n为整数，且φ_α=jπ(u²cotα-2utcscα+t²cotα)。

当分数阶次p=1时，有$ \alpha = \frac{{\rm{ \mathit{ π} }}}{2}$，A_α=1，由式(7)得

$ {f_1}\left( u \right) = \int_{ - \infty }^{ - \infty } {{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}ut}}f\left( t \right){\rm{d}}t} $

(9)

由上可见f₁(u)就是f(t)的普通傅里叶变换，同理f_-1(u)就是f(t)的普通傅里叶逆变换。由此可认为分数阶傅里叶变换是一种广义的傅里叶变换。因为参数$ \alpha = \frac{{p{\rm{ \mathit{ π} }}}}{2}$仅出现在核函数中三角函数的参数位置，所以参数p的定义以4为周期，因此只需考虑区间p∈-2, 2即可。当p=0时，f₀(u)=f(u)；当p=±2时，f_±2(u)=f(-u)。

3 有效转动角速度的估计

线性调频(Linear frequency modulation，LFM)信号在不同的分数阶傅里叶域上呈现不同的能量聚集性^{[13, 14]}。LFM信号的分数阶傅里叶变换相当于将信号的时频分布投影在不同旋转角度的频率轴上，在特定的旋转角度，LFM信号能量会得到高度集中。所以检测LFM信号的基本方法如下：在变量旋转角α的某个取值区间内，利用FrFT对信号进行处理，形成能量在参数(α, u)平面上的二维分布，并按阈值进行搜索即可实现。

单个LFM信号的时频分布及其FrFT如图 2所示，LFM信号与时间轴的夹角为β，斜率即为其调频率，旋转轴与时间轴的夹角为α，当α=β时即为该信号的傅里叶变换。当α与β正交时，该LFM信号在分数阶傅里叶域上的投影聚集在一点上，此时能量最强^[15]，该LFM信号的调频率为

$ K = - \cot \left( {{\alpha _0}} \right) $

(10)

在ISAR图像中，同一个距离单元内往往包含多个散射点。各散射点引入的多普勒频率的时频分布如图 3(a)所示。式(5)决定了这些信号在频率f轴上的截距，式(6)决定了这些信号的调频率，即图 2中的斜率。

由于调频率由散射点的纵向距离单元y_a和目标相对雷达的转速Ω决定，所以，同一纵向距离单元的各散射点的回波信号的调频率相同。本文只关心信号的调频率，所以同一距离单元的多个LFM信号的FrFT变换对应的幅度，会在同一个确定的变换参数下达到最大值，如图 3(b)所示。

设某纵向距离单元的观测信号为

$ s\left( t \right) = \sum\limits_{i = 0}^{M - 1} {{\sigma _i}\exp \left[ {{\rm{j}}2{\rm{ \mathsf{ π} }}\left( {{f_{ci}}t + \frac{1}{2}K{t^2}} \right)} \right]} $

(11)

式中：M为该纵向距离单元内包含的散射点的总个数；σ_i为第i点的散射强度，且

$ K = \frac{{2{\mathit{\Omega }^2}{y_{\rm{a}}}}}{\lambda } $

(12)

对该信号进行离散化，并做FrFT，采样频率为f_s，离散化后总点数为N。如果幅度最大时的角度为α₀，离散化后回波信号的调频率估计值将由式(10)变为

$ K = \frac{{ - \cot \left( {{\alpha _0}} \right)f_{\rm{s}}^2}}{N} $

(13)

通过式(12)可知，所有散射点具有相同的旋转角速度和波长，所以调频率K与纵向距离y₀成线性关系，可表示为K=ay_a，其中$\alpha = \frac{{2{\mathit{\Omega }^2}}}{\lambda } $。理论上使用单个距离单元的信号就可以估计调频率，进而完成横向定标。但是为了避免单个距离单元估计中可能会造成的误差，通常选择一组不同纵向距离单元进行处理，最后进行最小二乘拟合，估计得到最优的a值，则角速度的估计值为

$ \mathit{\Omega } = {\left| {\frac{{a\lambda }}{2}} \right|^{\frac{1}{2}}} $

(14)

因观测时间T已知，则可得到最终转角估计值：Δθ=ΩT，横向分辨率为$ {\rho _{\rm{a}}}{\rm{ = }}\frac{\lambda }{{(2\Delta \theta )}}$，即可完成横向定标。对基于FrFT的横向定标进行总结，其算法流程如图 4所示。

4 仿真实验

雷达发射线性调频信号，系统参数为：载频f_c=10 GHz；信号带宽B=400 MHz；脉宽T_p=80 μs；脉冲重复频率PRF=250 Hz；采样频率f_s=800 MHz；距离分辨率为0.375 m。仿真目标散射点模型如图 5所示，目标尺寸为24 m×56 m，反射系数均为单位值。设定目标转台成像旋转角速度ω=0.03 rad/s。

图 6为模拟目标未定标的ISAR图像，图像大小为246×278。由图 6可见，图像未做定标，距离维和方位维像素代表的真实尺寸不一致，模型ISAR图像发生了明显变形，需要进行ISAR图像定标，调整图像的横向尺寸后，使图像反映出目标的真实形状，以便估计目标的真实尺寸。

对第195个距离单元的信号进行FrFT，结果在参数(α, u)平面上的二维分布如图 7所示，FrFT阶数的取值范围定为0≤α≤2，FrFT域是指参数FrFT变换后的参数u。按照阈值对其进行二维搜索，得到峰值点所对应的旋转角α=1.578 7 rad/s，利用式(13)得到调频率的估计值为K=1.765 8 Hz/s。

本文选取60个距离单元进行计算，先利用3σ准则剔除野值，再利用最小二乘准则进行拟合。图 8为所选距离单元对应的纵向距离和调频率K值。利用最小二乘准则进行拟合，直线斜率的估计值为0.059 0 Hz/(m·s)^-1，代入式(14)得到角速度的估计值为0.029 7 rad/s，横向分辨率为0.453 4 m。利用计算得到的目标横向分辨率，对目标横向尺寸进行调整，得到定标后的图像如图 9所示。根据图 6和图 9可以计算出目标定标前后尺寸，如表 1所示。

定标前得到的目标尺寸为149.91×56.86，无法得到目标的真实尺寸，不利于后期的目标识别和特征提取。图 8为定标后的结果，给出了模拟目标横向和纵向各两个点的位置，可以进一步计算得到目标横向长度24.45 m，纵向长度57.32 m，与模拟目标的仿真尺寸(24 m×56 m)比较接近。

5 结束语

转角估计是ISAR成像中方位向定标的前提，更是特征提取、目标识别等应用的基础。本文在慢时间域对回波信号进行了分析，确定了转动角速度和回波信号调频率之间的关系，并利用分数阶傅里叶变换估计目标回波信号的调频率，从而得到目标的有效转动角速度，进而完成横向定标。最后对卫星目标的运动进行了仿真成像，使用本文方法对仿真数据进行了处理，实验结果证明了本文方法的有效性和可行性。