基于广义keystone和频率变标的微波光子ISAR高分辨实时成像算法

杨利超 高悦欣 邢孟道 盛佳恋(西安电子科技大学雷达信号处理国家重点实验室 西安 710071)

②(上海无线电设备研究所 上海 201109)

1 引言

逆合成孔径雷达成像(Inverse Synthetic Aperture Radar, ISAR)的距离分辨率与方位分辨率分别通过大带宽信号和目标与雷达间的等效转角获得[1]，并广泛应用于军事和民用探测领域。然而传统微波雷达由于受到微波器件的限制难以产生超大带宽信号，这一技术瓶颈限制了ISAR成像分辨率的提升。微波光子雷达是一种新体制雷达，它将微波理论与电子学结合起来[2—7]，能够产生大带宽和高载频的线性调频信号，进而实现更高分辨率的ISAR成像。

但针对微波光子雷达实时成像仍存在如下问题：在成像算法方面，大带宽和高载频的信号使得空变的距离弯曲问题不能忽略，相关法[8,9]等包络对齐方法无法完成距离徙动校正。并且常用的自聚焦算法[10,11](相位梯度自聚焦，最小熵自聚焦等)无法补偿距离弯曲项产生的空变二次相位误差，因此传统的距离多普勒域实时成像算法不再适用；在计算量方面，大带宽信号导致距离向采样点数增大。为保证方位多普勒谱不发生混叠，高载频信号需设置较大的脉冲重复频率，导致方位向成像需积累更多脉冲数。因此虽然波束域ISAR成像算法[11]可解决空变的距离弯曲问题，但由于计算量过大无法应用于微波光子ISAR实时成像。

综上所述，在信号具有大带宽和高载频特点的微波光子雷达ISAR实时成像处理中，需要设计兼具高效率和精确补偿的成像算法。针对实时成像处理中的难点，本文提出了一种基于广义keystone[12,13]和频率变标的微波光子ISAR高分辨实时成像算法。本算法首先利用广义keystone变换去除距离弯曲对包络徙动的影响，并且减小距离走动和保留了方位相位，然后利用滑窗和最小方差准则提取特显点相位，其次配合相位二次拟合粗估计和图像偏置法(Map Drift, MD)精确估计反演目标横向速度，最后利用速度估计结果结合频率变标(Frequency Scaling, FS)[14,15]算法完成空变的距离弯曲校正和方位匹配滤波成像，从而高效地实现微波光子雷达2维高分辨成像。

2 成像模型

微波光子雷达ISAR成像系统发射大带宽线性调频信号，一般采用解线频调(Dechirp)模式接收。接收信号回波表达式为

其中，t为距离向快时间(相对于距离中心)，tm为方位向慢时间，rect(·)为窗函数，Tp和Ta分别为距离和方位的窗长，γ为线性调频信号调频率，fc为信号载频，c为光速，q为散射点总数，¾i为目标散射点i的后向散射系数，Ri(tm)为散射点i的瞬时斜距，Rref为Dechirp接收的参考距离。将式(1)沿距离向作傅里叶变换即可得到目标1维距离向

其中，f为距离频率，目标散射点i的距离向聚焦位置为f=-2γRΔi(tm)/c，各个散射点的距离徙动曲线可用斜距历程RΔi(tm)描述。

微波光子雷达ISAR实时成像中假设短时间内目标平稳运动，此时目标运动与成像几何如图1所示。XOY和xoy分别为雷达坐标系和目标本地坐标系，OY轴垂直于目标航线，OX轴平行于目标航线，为目标某一散射点，Rbi为i距OX轴的垂直距离，xi为i在本地坐标系内的横向位置，目标运动速度为v，设目标中心点o与雷达距离最短的时刻为tm的0时刻，散射点i的斜距历程和方位多普勒表示为

图1 目标运动示意图Fig.1 The geometry of target movement

将式(5)代入式(4)可得目标距离徙动在多普勒域的响应表达式

将式(4)在tm=0点进行泰勒展开得距离徙动在时域的响应表达式

因xiRbi，式(7)可以简化为

式(8)中关于tm的1次项为散射点i的多普勒中心，对应i的聚焦位置。2次项为距离弯曲项，该项系数与目标距离向位置相关，另外式(6)中的响应函数曲率与目标位置相关，由此可知回波包络存在空变的距离弯曲现象。回波包络距离徙动的响应曲线示意图如图2所示。

由图2可知，空变的距离弯曲项造成各个散射点的徙动响应曲线的曲率不同，进而导致不同时刻的散射点距离徙动差不同(ΔR1ΔR2)。距离弯曲的影响表现为其徙动值与距离分辨率的比值，即常规的平动补偿方法有效的前提为目标距离徙动的空变部分小于1/8个距离单元，具有超高距离分辨率的微波光子雷达往往不满足该前提条件。而传统ISAR实时处理往往忽略距离弯曲的空变性，采用平动补偿和相位校正结合的算法进行连续成像。并且在微波光子雷达高载频的情况下，空变的距离弯曲项引起的2次相位误差将造成严重的图像模糊，常规的自聚焦算法无法校正空变误差相位。

图2 回波包络响应曲线Fig.2 Echo envelope curve in azimuth

3 频率变标成像

变标类成像算法[14,15]是完成空变距离弯曲校正的有效手段。并且该类算法计算步骤主要包含快速傅里叶变换(Fast Fourier Transform, FFT)与矩阵点乘操作，相比于插值类算法计算量较小，适用于ISAR实时成像。频率变标算法是适用于Dechirp接收模式的距离徙动校正算法，在目标横向运动速度v已知时(本文算法的横向速度估计方法将在第4节介绍)，可有效去除距离弯曲的空变性。在此对频率变标算法去空变的原理作简单介绍，文献[14]指出式(1)中剩余视频相位(Residual Video Phase, RVP)可表示为与原信号时域卷积的形式，令，式(1)可简写为

其中，rectr(·)表示信号的距离窗函数，*表示信号卷积，将式(9)展开为波束域表达式[14]

其中，a(fa)称为变标因子，式(11)中关于Kr的1次项对应空变的距离弯曲项，A(xi)为方位成像所需相位项，Hp(Kr)为泰勒展开的高次项。为去除距离弯曲的空变性，频率变标算法引入的变标函数为

将式(11)与变标函数式(15)相乘，得到信号

式(16)中关于Kr的1次项解除了距离与方位频率的耦合，使得各个散射点的弯曲程度一致。去除距离弯曲空变性后，FS算法仍需进行RVP校正，逆频率变标，距离徙动校正3步完成包络对齐，最后利用A(xi)项进行匹配滤波成像，滤波函数为exp[j4fca(fa)/c]。FS算法的具体操作步骤与变标过程的证明可参考文献[14]。FS操作过程中需进行3次距离FFT, 2次方位FFT与6次矩阵点乘(补偿函数)。假设信号的方位向点数为M，距离向点数为N, FS算法成像的计算量为MN(1.5*log2M+log2N+6)次复数乘法。

4 基于广义楔石变换的横向速度估计方法

本文第3节内容主要分析了频率变标算法在微波光子ISAR实时成像中的可行性，而应用FS算法的前提是目标横向运动速度已知。但横向速度的概念在ISAR成像处理中不常使用，并且常规的ISAR目标运动参数估计方法[16,17]主要用于目标径向速度、加速度等参数的估计和信号平动补偿。因此本文提出了一种基于广义keystone的ISAR横向速度估计方法。为保证参数估计精度，本文算法从目标的相位信息反演其运动参数，具体分两步处理：提取信号方位向相位调频率和反演目标横向运动速度。

由斜距近似表达式(8)可知目标散射点i的方位调频率与其距离向位置Rbi有关，即信号调频率具有空变性，因此准确反演横向速度需提取目标特显点相位历程。利用微波光子雷达距离分辨率高的特点，短时间内可使用广义keystone变换将回波包络线分离，其变换关系为

将式(18)，式(8)代入式(1)中并忽略RVP项，可得广义楔石变换(Generalized Keystone Trans−form, GKT)变换后的信号表达式为

将式(19)在t=0处泰勒展开，并稍加整理得

式(20)中关于t的相位项对应散射点的距离向聚焦位置，其距离弯曲经过GKT变换后被消除，并且距离走动减小为原来的1/2。此时目标各个散射点的包络线很大程度上完成了分离，并且各个散射点短时间内的距离徙动值显著减小，可以认为在一个小窗长特显点的包络处于同一距离单元。而式(20)关于fc仍然保留了横向速度v的2次项，因此可以通过提取特显点方位调频率反演v的精确值。此时可通过滑窗和最小方差准则选取特显点，具体的选点和相位提取流程如下：

步骤1 确定起始位置。对时域回波信号作广义keystone变换，将变换后的回波包络线作为输入，设定滑窗长度L，设定方差门限值Ar。搜索幅度方差最小的距离单元S(n，τm),n表示距离单元坐标。以S(n，τm)为中心取出子数据块进行相位提取；

步骤2 相位提取。设定每次滑窗有重叠部分，在选取的字数据块中以信号幅度最小方差为准则选取每次滑窗的特显点距离单元，并且提取其相位历程，记为[h1，h2，···，hk],k为滑窗次数。在滑窗过程中记录中心点对于斜距Rbi；

步骤3 相位拼接。由于相位解缠绕时可能存在初相差异，因此利用相邻子孔径相位的重叠部分，提取相邻相位历程的初相差异，最后将子孔径相位历程[h1，h2，···，hk]拼接为特显点i的相位历程Hi。

得到特显点i的相位历程Hi后估计分为两步：2次曲线拟合粗估计信号调频率；图像偏置法[1]精估计剩余调频率误差。设2次曲线拟合粗估计的调频率值为k1，对Hi作粗补偿

横向速度估计方法可分解为选点和估计两部分。选点操作中广义keystone计算复杂度较高，其余操作计算量可忽略。广义keystone操作时的插值操作可由线性调频Z变换(Chirp Z−Transform,CZT)变换代替，其中包括4次方位向FFT和2次复数点乘。假设信号的方位向点数为M，距离向点数为N，选点操作的计算量为2MN(log2N+1)次复数点乘。估计操作中图像偏置法估计剩余误差时计算复杂度较高，其余操作计算量可忽略。图像偏置法中包括5次FFT与1次复数点乘，估计操作的计算量为2.5Nlog2N+N次复数点乘。因此速度估计的计算复杂度为2MN(log2N+1)+2.5Nlog2N+N次复数乘法。

5 仿真实验与实测验证

结合横向速度估计方法和频率变标算法，实时成像处理流程图如图3所示。

5.1 仿真实验分析

为验证本文算法的有效性，本文首先进行点阵的速度估计实验和成像仿真实验。仿真参数如表1所示。

图3 算法流程图Fig.3 Flow diagram of proposed algorithm

表1 仿真数据参数Tab.1 Simulation parameters

仿真实验在10 m×10 m的网格中布了3×3的点阵，用以验证算法的可行性。图4为提取特显点相位后的2阶拟合结果，速度粗估计值为65.54 m/s。经粗速补偿后，MD精估计速度为65.08 m/s，与实际值相符。将精确估计的横向速度输入到FS算法中，得到的包络对齐结果为图5(a)，方位匹配滤波成像结果为图5(b)。本文实验采用传统RD算法作为对比，其处理流程包括包络对齐与相位校正两步。图6(a)为信号包络，可以看到包络存在严重的空变距离弯曲现象，与本文第2节图2的描述一致，因此传统相关包络对齐法无法处理。图6(b)为RD算法的成像结果，存在严重的模糊现象。而本文算法完成了空变的距离弯曲校正，且获得了较好的聚焦效果。

图4 2次相位拟合结果Fig.4 Quadratic phase curve fitting

图5 本文算法仿真处理结果Fig.5 Simulation processing results with proposed algorithm

图6 传统RD算法仿真处理结果Fig.6 Simulation processing results with traditional RD process

图7 不同信噪比下速度估计的均方误差曲线Fig.7 MSE curve of velocity estimation

为验证本文所提横向速度估计方法的抗噪性能，本文在不同信噪比下分别进行了100次蒙特卡洛(Monte−Carlo)实验，计算了各个信噪比下的均方误差(Mean Squared Error, MSE)，信噪比变化范围为—5～30 dB。均方误差计算结果如图7所示，在信噪比大于5 dB时速度估计均方误差曲线明显下降并且趋于稳定(在0.2上下波动)。这意味着本文所提速度估计算法对噪声具有较强的鲁棒性。在低信噪比条件下，能够保证目标横向速度的可靠估计。

5.2 实测数据验证

实验参数如表2所示，目标为民航飞机。

图8为提取特显点相位后的2阶拟合结果，两步速度估计后得到的目标横向速度为85.54 m/s。将横向速度输入到FS算法中，得到的包络对齐局部结果为图9(a)，方位匹配滤波成像结果为图9(b)。将本文算法与传统RD算法作对比，图10(a)为传统算法包络对齐结果，空变的距离弯曲现象未能解决。图10(b)为RD算法的成像结果，存在严重的模糊现象。本文处理算法的包络对齐结果与聚焦效果明显优于传统RD算法。尤其在机翼部分，本文算法处理结果目标精细聚焦并且结构清晰。

表2 实测数据参数Tab.2 Measured data parameters

图8 2次相位拟合结果Fig.8 Quadratic phase curve fitting

图9 本文算法实测处理结果Fig.9 Real data processing results with proposed algorithm

图10 传统RD算法实测处理结果Fig.10 Real data processing results with traditional RD process

6 结论

微波光子雷达具有超高2维分辨率，开展关于微波光子雷达实时成像方面的研究对其应用具有重要作用。传统ISAR实时成像算法一般不考虑回波信号的距离弯曲项对包络和相位的影响，无法在微波光子雷达实时成像中应用。本文提出了可解决距离弯曲问题且计算效率可观的实时成像算法。首先利用广义keystone变换减小目标距离徙动对包络的影响，然后滑窗选择特显点以估计其方位调频率和反演目标横向速度，最后在获得精确的横向速度后，本文算法采用FS实现空变距离徙动校正和方位匹配成像。算法的计算效率高并且获得了较好的聚焦效果。