一种基于微波光子超高分辨雷达机翼振动参数估计方法

符吉祥 邢孟道* 徐 丹 王安乐(西安电子科技大学雷达信号处理国家重点实验室 西安 710071)

②(西安电子科技大学信息感知技术协同创新中心 西安 710071)

③(空军预警学院 武汉 430019)

1 引言

得益于全天时全天候和超远距离探测的优势，雷达被广泛应用到军事和民事的应用中，如战场监视和遥感测绘等。逆合成孔径雷达，因为其能对非合作目标进行探测和成像，获得目标的形状、尺寸、姿态和运动信息，进而可以对非合作目标进行识别，在国防安全领域有着重要应用，受到各个国家的高度重视和广泛研究。随着雷达技术和需求的不断提高，雷达向高分辨、多功能和多任务方向发展。然而，传统微波雷达受限于微波器件的“电子瓶颈”[1]，无法产生大带宽信号，使得雷达分辨率受限。近年来，微波光子技术的提出打破了这一瓶颈。微波光子技术将光子技术应用到雷达信号的产生和调制阶段，充分利用了光子的大带宽低损耗和抗电磁干扰的优势，使得雷达能够发射超大带宽的高质量信号。美国自然杂志评价其为照亮雷达未来的技术[2,3]。目前，已经研制成功的微波光子雷达所能产生的最大信号带宽为12 GHz[4]，这大大超过传统的微波雷达。超高的分辨率，使得雷达能够对目标进行更加精细的探测和感知，如目标微振动的观测和提取，多频段同时多任务观测等。然而，分辨率的大幅度提升也给雷达成像和参数估计带来了新的挑战。一是高分辨率会带来严重的越距离单元徙动[5]；二是运动误差补偿的精度提高，为此需要研究新的雷达成像和参数估计方法。

图1 微波光子实测飞机目标1维距离像Fig.1 One dimensional range profile of airplane measured by microwave photonic radar

在对10 G带宽的微波光子雷达实测民航飞机的数据处理过程中，我们发现，飞机目标在转向过程中，因为目标状态的改变，机翼与飞机的主体存在运动不一致的现象，如图1所示的飞机目标的高分辨1维距离像，图1(a)是飞机整体包络对齐后的1维距离像包络，机身的包络如图1(b)所示，可以发现，机身散射点包络变化较为平稳，部分机翼的包络如图1(c)所示，可以发现机翼回波的包络线存在明显的抖动现象。机翼的这种非平稳抖动使得传统的成像方法无法对其进行聚焦成像。机翼的振动相对于平动以及机身相对于雷达的转动来说，其具有幅度小频率高的特点，所以可以将其视为飞机的微动特征。利用微动特征可以评估飞机的运动状态以及飞机的安全状态，进一步为飞机目标的识别和监测作为依据。目前微动参数提取的相关方法多是应用在弹道类目标上，因为弹头散射结构较为简单，散射点较少，所以其可以利用基于时频分析的方法进行微动参数估计[6,7]。具体做法是，首先利用包络对齐平动补偿[8,9]提取单个散射点回波，再通过对其进行时频分析获得其时间频率分布曲线，最后通过对时频曲线进行相关处理即可得到弹头的微动参数。但是对于飞机目标来说，一方面，机翼结构较弹头要复杂得多，纯净的单个散射点回波较难获取，而其他散射点的回波会对参数估计造成严重干扰；另一方面，因为机翼回波的信噪比较低，直接采用时频分析很难从其时频分布图中提取时频曲线。基于此，本文提出一种基于子孔径序列图的振动参数估计方法。首先将机翼的振动建模为单摆运动，然后通过分析信号波数谱表达式，得到短时间机翼散射点聚焦位置和单摆运动参数之间的关系，进一步通过机翼子孔径成像结果序列图中散射点位置的变化来反演单摆运动的参数。最后提出一种修正的极坐标格式算法(Modified Polar Format Algorithm , MPFA)，通过非均匀傅里叶变换和估计的振动参数对振动的机翼进行高效的聚焦成像，并通过构造成像结果熵值优化函数对振动参数进行精估计。因为是利用图像序列图来进行摆动参数的反演，这充分地利用了2维相干积累增益，使得在低信噪比下的参数估计成为可能。另外相比于传统的时频分析方法，一方面，序列图能够提供距离和多普勒两维的信息；另一方面，序列图中各个散射点之间是解耦的，没有交叉项干扰的影响。

2 机翼振动建模分析和振动参数估计

2.1 运动模型分析

根据力学分析可知，机翼的振动与单摆运动存在一致性，所以将其建模为单摆运动。如图2所示，以飞机整体所在平面作为XOY平面，X轴指向为在参考时刻沿着机身方向指向机头的方向，Z轴为飞机平面的法线方向。雷达视线与Z轴的夹角为φ0，其在XOY平面上的投影向量与Y轴之间的夹角为β0。机翼沿着Z轴振动，振动过程中，机翼与XOY平面之间的夹角为ϑ，最大夹角为ϑ0，根据单摆转角公式可知，夹角随慢时间正弦变化，即

其中，ϑ0表示振动角，fv表示振动频率，φ0表示振动初相，表示方位慢时间。在飞机沿着图2中曲线轨迹运动的时候，飞机整体还有相对于雷达的平动和转动。对于机身上的散射点p1，因为其不存在振动，所以其相对于雷达的运动可以表示为

其中，R0表示飞机旋转中心O到雷达的初始距离，vr和ar分别表示飞机相对于雷达运动的径向速度和径向加速度。转动分量是随散射点位置发生空变的，其形式为(平面波近似)

其中，(xp1，yp1)表示散射点p1参考时刻的坐标位置，θ表示飞机相对于雷达转动的角度，其是慢时间的函数，即θ(tm)=ωtm,ω表示转速。对于机翼上的散射点，其不仅包含以上的运动，还包含机翼振动对其的运动调制。通过运动分解可以得到机翼上散射点p2各个时刻的坐标位置如式(5)

图2 飞机运动几何模型Fig.2 Geometry model of airplane

其中，(xp2，yp2)表示散射点p2参考时刻的坐标位置。所以散射点p2到雷达的斜距变化为

假设雷达发射线性调频信号，使用去斜方式接收回波，经过剩余视频项补偿之后，机翼散射点雷达回波信号的波数谱表达式为[10]

2.2 振动参数估计方法

首先通过粗成像对机翼机身回波进行分离。因为机翼存在周期性振动而机身只存在平稳的运动，所以在粗成像结果中，机翼和机身上的散射点散焦情况存在明显的差异，即机翼散射点散焦成线性分布而机身散射点散焦呈现块状分布。另外，从粗成像的轮廓中也可以轻易地将机翼和机身进行区分。如图4所示，图4(a)是飞机整体粗成像结果，可以发现，尽管飞机散焦严重，但是仍能够从轮廓中轻易区分机身和机翼，图4(b)是机身散射点块状散焦结果，图4(c)是机翼散射点线性散焦结果。最后通过加窗滤波将成像结果中的机身回波滤出，从而实现机身机翼回波分离。

对机身回波使用keystone算法[13]校正其线性距离走动，得到

其中，sinc(x)=sin(x)/(x),表示径向波数中心。观察式(9)可以发现，keystone变换之后，包络已经对齐，然而散射点回波存在和其距离坐标yp线性相关的方位向二次相位。通过对不同距离单元信号进行调频率估计并对调频率进行线性拟合即可得到转速ω，然后再构造调频率补偿函数

图3 转台模型Fig.3 Turntable model

其中，r表示式(9)中各个距离单元的位置。对式(9)补偿之后直接方位向FFT即可得到聚焦图像

从式(11)中可以发现，成像结果中散射点的聚焦位置相对于其在XOY中的坐标发生了旋转和缩放，旋转角度为方位角β0，缩放因子为sinφ0。所以可以从定标图像中，通过计算机身相对水平方向的夹角来估计β0，通过测量定标图像中机身长度与宽度与飞机真实的尺寸的比值来估计φ0。

对分离后的机翼部分回波进行平动补偿，将其旋转中心和振动中心补偿为机翼中心点O′，如图3(b)所示，再对其进行子孔径成像，子图像的距离和多普勒如式(12)

图4 飞机粗成像结果Fig.4 Coarse imaging results of airplane

其中，Δyq的值可由以下方程组解得

观察机翼散射点回波，发现其距离和方位是解耦的，利用极坐标格式算法的思想，本文提出一种修正的极坐标格式算法，即

将式(17)代入到式(8)中的机翼回波后得到

对式(18)进行2维傅里叶变换即可得到距离的机翼图像。为了提高算法的效率，使用快速非均匀傅里叶变换(Non−Uniform Fast Fourier Transform,NUFFT)[14,15]高效地实现插值成像的操作。

其中，E[·]表示计算熵值操作，表示图像矩阵，矩阵大小为M×N,si为 其第i个 元素。表示图像的总能量。所以构造的最小熵优化函数如式(20)

其中，wing(kr，tm)表示机翼波数域回波，表示利用参数对进行非均匀快速傅里叶变换。使用粗估计的振动参数作为初值，利用式(20)进行参数搜索，得到振动参数的精确估计。这里采用的搜索策略是利用坐标下降法的思想，每次只对1个参数进行搜索，其他参数固定不变，循环迭代，直至相邻两次迭代的图像熵值变化小于门限停止迭代。算法的流程图如图5所示。

3 仿真实验及实测数据验证

本文使用Ka波段雷达参数对飞机散射点模型进行参数估计，为了简化仿真，将飞机两侧机翼振动参数设置为相同的，仿真参数如表1所示。飞机模型长35 m，宽30 m，如图6所示。雷达视线方位角为80°，俯仰角为85°。

图5 算法流程图Fig.5 Flow chart of the proposed algorithm

表1 仿真雷达参数Tab.1 Simulation radar parameters

图6 飞机散射点模型Fig.6 Scatter model of airplane

图7 机身聚焦成像定标结果Fig.7 Focused and scaled result of airplane body

图8 不同时刻机翼子孔径成像结果Fig.8 Sub−aperture imaging results of airplane wing at different time

图7是机身聚焦成像定标结果，从图7中可以看出机身主轴和Y轴存在夹角，通过之前的分析可知，其夹角大小和雷达视线的方位角是相等的，通过机头标记点A的坐标位置，可得雷达视线的方位，通过测量飞机宽度BC与实际中机身宽度之间的比例，可以得到对俯仰角的估计，即，其中W0为飞机的真实宽度。

图8是机翼不同子孔径成像结果，可以发现其方位多普勒存在周期性的变化。以图8(a)中用圆圈标记的散射点1,2,3为例，分别提取其不同子孔径方位多普勒和距离的变化，提取的多普勒变化曲线如图9(a)所示，可以发现曲线存在明显的阶梯跳跃，这是因为子孔径图像多普勒分辨率有限，不能反映多普勒的连续变化，为了降低跳变对参数估计的影响，对提取的多普勒曲线进行平滑滤波处理，得到的曲线如图9(b)所示。为了进行振动角的估计，还需要对曲线进行去均值处理，处理结果如图9(c)所示。通过测量多普勒曲线波峰之间的时间差即可估计出振动频率，即。提取的散射点1, 2, 3的距离变化曲线如。通过式(16)可以得到初相的估计值。最后对振动角和初相参数进行精搜索，机翼回波熵值优化函数的代价曲面如图11所示，可以发现，熵值是平滑变化的，且是凸的，所以使用本文所提搜索策略能够快速得到参数的精确估计结果。最终振动参数的估计结果如表2所示。机翼区域精细化处理前后的对比图像如图12所示。

图9 提取的不同散射点的多普勒变化曲线Fig.9 Extracted Doppler curves of different scatterers

图10 提取的不同散射点的距离变化曲线Fig.10 Extracted range curves of different scatterers

图11 图像熵值优化函数代价曲面Fig.11 Cost surface of entropy optimization function

表2 仿真数据机翼振动参数估计结果Tab.2 Vibration parameters estimation result of simulation data

为了进一步验证本文方法的有效性，下面对实测数据进行处理，实测机身成像定标结果如图13所示，通过测量机身轴线与Y轴之间的夹角，可以得到方位角的估计，即。观测飞机的机型为A320，通过公开资料可知其机身宽度为3.95 m，通过测量定标图像中机身宽度，可以得到俯仰角的估计。图14是振动较为明显的机翼不同子孔径的距离多普勒成像结果，通过提取图中用圆圈标记的散射点图10所示，通过在距离变化曲线上选取两个时间点，联立式(14)和式(15)求解得到振动角的估计值1的距离和多普勒变化，如图15所示。最小熵优化函数的代价平面如图16所示，振动参数估计结果如表3所示。机翼区域精细化处理前后的对比图像如图17所示，从图中用圆圈标记的散射点的聚焦情况可以看出，使用本文方法对机翼区域精细化处理之后，散射点的聚焦质量有所改善。

图12 机翼区域精细化处理前后结果对比Fig.12 Before and after accurate processing of airplane wings

图13 实测数据机身聚焦成像定标结果Fig.13 Focused and scaled airplane body result of the measured data

图15 提取的散射点的距离和多普勒变化曲线Fig.15 Range and Doppler curves of extracted scatter

图16 图像熵值优化函数代价曲面Fig.16 Cost surface of image entropy optimization function

表3 实测数据机翼振动参数估计结果Tab.3 Vibration parameters estimation result of measured data

4 结束语

本文针对机翼在飞机非平稳运动的状态下振动的现象，提出一种利用微波光子超高分辨雷达估计机翼振动参数的方法。首先通过对平稳运动的机身进行聚焦成像和定标处理，从中估计雷达视线角。然后通过对机翼振动进行建模，推导了机翼回波的解析式，并分析了振动对成像的影响。针对振动导致的机翼成像结果散焦的问题，本文首先通过子孔径序列成像提取散射点距离和多普勒曲线，再联合雷达视线角参数进行振动参数的初步估计，然后提出了修正的极坐标格式算法，构造了最小熵优化函数对振动参数进行精确估计。利用估计的振动参数，修正的极坐标格式算法能够对机翼进行聚焦成像。本文对仿真数据和飞机实测数据进行了试验，试验结果验证了本文方法的有效性。

图17 机翼区域精细化处理前后结果对比Fig.17 Before and after accurate processing of airplane wings