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

符吉祥 邢孟道 徐丹 王安乐

引用本文:
Citation:

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

    作者简介: 符吉祥(1992–),男,安徽蚌埠人,博士生,研究方向为ISAR运动补偿和成像。E-mail: jixiang_xd@126.com;邢孟道(1975–),男,浙江绍兴人,博士,教授,西安电子科技大学前沿交叉研究院副院长,研究方向为雷达成像、动目标检测、目标识别。E-mail: xmd@xidian.edu.cn;徐   丹(1990–),女,陕西咸阳人,博士生,研究方向为3维ISAR图像重构、电磁成像。 E-mail: 1045759961@qq.com;王安乐(1988–),男,山东临沂人,博士,中国人民解放军空军预警学院讲师,研究方向为微波光子雷达、光生微波基准信号等。E-mail: anlehit@163.com.
    通讯作者: 邢孟道, xmd@xidian.edu.cn
  • 基金项目:

    国家杰出青年自然基金(61825105)

  • 中图分类号: TN957

Vibration-parameters Estimation Method for Airplane Wings Based on Microwave-photonics Ultrahigh-resolution Radar

    Corresponding author: XING Mengdao, xmd@xidian.edu.cn ;
  • Fund Project: The National Science Fund for Distinguished Young Scholars (61825105)

    CLC number: TN957

  • 摘要: 机翼在飞机运动不平稳状态下存在振动,在微波光子超高分辨雷达观测下,这种振动会造成机翼难以聚焦,针对这一难题,该文提出一种基于微波光子超高分辨雷达的机翼振动参数估计方法。首先通过粗成像将机身和机翼回波进行分离,再通过对机身成像和定标结果估计雷达视线角(LOS)。然后对机翼进行子孔径序列成像,提取散射点的距离和多普勒变化曲线,再联合雷达LOS以及距离和多普勒曲线对振动参数进行粗估计,最后通过修正的极坐标格式算法(MPFA)以及构造最小熵优化函数对振动参数进行精估计。该文首次提出了修正的极坐标格式算法,其能够对复杂运动的目标进行距离和方位向的解耦,如震动的机翼和摇摆的舰船等。仿真和实测数据的处理结果验证了该方法的有效性和实用性。
  • 图 1  微波光子实测飞机目标1维距离像

    Figure 1.  One dimensional range profile of airplane measured by microwave photonic radar

    图 2  飞机运动几何模型

    Figure 2.  Geometry model of airplane

    图 3  转台模型

    Figure 3.  Turntable model

    图 4  飞机粗成像结果

    Figure 4.  Coarse imaging results of airplane

    图 5  算法流程图

    Figure 5.  Flow chart of the proposed algorithm

    图 6  飞机散射点模型

    Figure 6.  Scatter model of airplane

    图 7  机身聚焦成像定标结果

    Figure 7.  Focused and scaled result of airplane body

    图 9  提取的不同散射点的多普勒变化曲线

    Figure 9.  Extracted Doppler curves of different scatterers

    图 10  提取的不同散射点的距离变化曲线

    Figure 10.  Extracted range curves of different scatterers

    图 11  图像熵值优化函数代价曲面

    Figure 11.  Cost surface of entropy optimization function

    图 12  机翼区域精细化处理前后结果对比

    Figure 12.  Before and after accurate processing of airplane wings

    图 8  不同时刻机翼子孔径成像结果

    Figure 8.  Sub-aperture imaging results of airplane wing at different time

    图 13  实测数据机身聚焦成像定标结果

    Figure 13.  Focused and scaled airplane body result of the measured data

    图 14  不同子孔径机翼成像结果序列图

    Figure 14.  Airplane wing imaging results of different sub-apertures

    图 15  提取的散射点的距离和多普勒变化曲线

    Figure 15.  Range and Doppler curves of extracted scatter

    图 16  图像熵值优化函数代价曲面

    Figure 16.  Cost surface of image entropy optimization function

    图 17  机翼区域精细化处理前后结果对比

    Figure 17.  Before and after accurate processing of airplane wings

    表 1  仿真雷达参数

    Table 1.  Simulation radar parameters

    信号带宽(GHz)脉冲宽度(μs)载频(GHz)PRF (Hz)采样率(Hz)脉冲数振动角(°)振动频率(Hz)初相(°)
    1015035666.665008192110
    下载: 导出CSV

    表 2  仿真数据机翼振动参数估计结果

    Table 2.  Vibration parameters estimation result of simulation data

    振动频率(Hz)振动角(°)振动初相(°)
    实际值估计值实际值估计值实际值估计值
    10.9910.980–0.18
    下载: 导出CSV

    表 3  实测数据机翼振动参数估计结果

    Table 3.  Vibration parameters estimation result of measured data

    振动频率(Hz)振动角(°)振动初相(°)
    2.480.0712.93
    下载: 导出CSV
  • [1] 潘时龙, 张亚梅. 微波光子雷达及关键技术[J]. 科技导报, 2017, 35(20): 36–52.PAN Shilong and ZHANG Yamei. Microwave photonic radar and key technologies[J]. Science &Technology Review, 2017, 35(20): 36–52.
    [2] MCKINNEY J D. Technology: Photonics illuminates the future of radar[J]. Nature, 2014, 507(7492): 310–311. doi: 10.1038/507310a
    [3] CAPMANY J and NOVAK D. Microwave photonics combines two worlds[J]. Nature Photonics, 2007, 1(6): 319–330. doi: 10.1038/nphoton.2007.89
    [4] YAO Yao, ZHANG Fangzheng, ZHANG Ying, et al. Demonstration of ultra-high-resolution photonics-based Ka-band inverse synthetic aperture radar imaging[C]. Optical Fiber Communication Conference, San Diego, California, USA, 2018: Th3G.5. doi: 10.1364/OFC.2018.Th3G.5.
    [5] CHEN C C and ANDREWS H C. Target-motion-induced radar imaging[J]. IEEE Transactions on Aerospace and Electronic Systems, 1980, AES-16(1): 2–14. doi: 10.1109/TAES.1980.308873
    [6] 李彦兵, 张曦文, 李飞, 等. 一种大加速度机动目标微动参数估计方法[J]. 电子与信息学报, 2017, 39(1): 82–87. doi: 10.11999/JEIT160261LI Yanbing, ZHANG Xiwen, LI Fei, et al. Estimation of micro-motion feature for large accelerated target[J]. Journal of Electronics &Information Technology, 2017, 39(1): 82–87. doi: 10.11999/JEIT160261
    [7] 张翼, 朱玉鹏, 黎湘. 基于微多普勒特征的目标微动参数估计[J]. 信号处理, 2009, 25(7): 1120–1124. doi: 10.3969/j.issn.1003-0530.2009.07.022ZHANG Yi, ZHU Yupeng, and LI Xiang. Micro-motion parameter estimation of ballistic missile target based on micro-Doppler feature[J]. Signal Processing, 2009, 25(7): 1120–1124. doi: 10.3969/j.issn.1003-0530.2009.07.022
    [8] ZHU Daiyin, WANG Ling, YU Yusheng, et al. Robust ISAR range alignment via minimizing the entropy of the average range profile[J]. IEEE Geoscience and Remote Sensing Letter, 2009, 6(4): 204–208. doi: 10.1109/LGRS.2008.2010562
    [9] PENG Shibao, XU Jia, WANG Libao, et al. A new ISAR range alignment method based on particle swarm optimizer[C]. Proceedings of 2009 Asian-Pacific Conference on Synthetic Aperture Radar, Xi’an, Shaanxi, China, 2009: 618–621. doi: 10.1109/APSAR.2009.5374278.
    [10] 保铮, 邢孟道, 王彤. 雷达成像技术[M]. 北京: 电子工业出版社, 2005: 1–336.BAO Zheng, XING Mengdao, and WANG Tong. Radar Imaging Technology[M]. Beijing: Publishing House of Electronics Industry, 2005: 1–336.
    [11] ZHANG Lei, SHENG Jialian, DUAN Jia, et al. Translational motion compensation for ISAR imaging under low SNR by minimum entropy[J]. EURASIP Journal on Advances in Signal Processing, 2013, 2013(1): 33. doi: 10.1186/1687-6180-2013-33
    [12] 符吉祥, 孙光才, 邢孟道. 一种大转角ISAR两维自聚焦平动补偿方法[J]. 电子与信息学报, 2017, 39(12): 2889–2898. doi: 10.11999/JEIT170303FU Jixiang, SUN Guangcai, and XING Mengdao. A two dimensional autofocus translation compensation method for wide-angle ISAR imaging[J]. Journal of Electronics &Information Technology, 2017, 39(12): 2889–2898. doi: 10.11999/JEIT170303
    [13] LI Xiaolong, CUI Guolong, YI Wei, et al. Range migration correction for maneuvering target based on generalized keystone transform[C]. 2015 IEEE Radar Conference, Arlington, VA, USA, 2015: 95–99. doi: 10.1109/RADAR.2015.7130977.
    [14] GAO Jingkun, DENG Bin, QIN Yuliang, et al. Efficient terahertz wide-angle NUFFT-based inverse synthetic aperture imaging considering spherical wavefront[J]. Sensors, 2016, 16(12): E2120. doi: 10.3390/s16122120
    [15] MA Zheng, ZHANG Yong, and ZHOU Zhennan. An improved semi-Lagrangian time splitting spectral method for the semi-classical Schrödinger equation with vector potentials using NUFFT[J]. Applied Numerical Mathematics, 2017, 111: 144–159. doi: 10.1016/j.apnum.2016.08.015
  • [1] 马兵强 . 滑动聚束FMCW-SAR的子孔径波数域成像算法. 雷达学报, 2013, 2(3): 319-325. doi: 10.3724/SP.J.1300.2013.13044
    [2] 赵雨露张群英李超纪奕才方广有 . 视频合成孔径雷达振动误差分析及补偿方案研究. 雷达学报, 2015, 4(2): 230-239. doi: 10.12000/JR14153
    [3] 马萌李道京杜剑波 . 振动条件下机载合成孔径激光雷达成像处理. 雷达学报, 2014, 3(5): 591-602. doi: 10.3724/SP.J.1300.2014.13132
    [4] 田雪梁兴东李焱磊董勇伟 . 基于子孔径包络误差校正的SAR高精度运动补偿方法. 雷达学报, 2014, 3(5): 583-590. doi: 10.3724/SP.J.1300.2014.14068
    [5] 丁振宇谭维贤王彦平洪文吴一戎 . 基于波数域子孔径的机载三维SAR偏航角运动误差补偿. 雷达学报, 2015, 4(4): 467-473. doi: 10.12000/JR15016
    [6] 詹学丽王岩飞王超刘碧丹 . 一种基于脉冲压缩的机载条带SAR重叠子孔径实时成像算法. 雷达学报, 2015, 4(2): 199-208. doi: 10.12000/JR14126
    [7] 王小青孙海清种劲松 . 从SAR 子孔径序列图像中提取海洋动态特征的改进相位相关法(英文). 雷达学报, 2013, 2(1): 30-38. doi: 10.3724/SP.J.1300.2013.13016
    [8] 时晨光汪飞周建江陈军 . 基于低截获概率优化的雷达组网系统最优功率分配算法(英文). 雷达学报, 2014, 3(4): 465-473. doi: 10.3724/SP.J.1300.2014.13140
    [9] 赵娟白霞 . 一种适用于TDOMP算法的测量矩阵优化方法. 雷达学报, 2016, 5(1): 8-15. doi: 10.12000/JR15131
    [10] 赵永科吕晓德 . 一种联合优化的无源雷达实时目标检测算法. 雷达学报, 2014, 3(6): 666-674. doi: 10.12000/JR14005
    [11] 鲁彦希何子述程子扬刘爽利 . 多目标跟踪分布式MIMO雷达收发站联合选择优化算法. 雷达学报, 2017, 6(1): 73-80. doi: 10.12000/JR16106
    [12] 王建峰林赟郭胜龙喻玲娟洪文 . 圆迹SAR的建筑物全方位优化成像方法研究. 雷达学报, 2015, 4(6): 698-707. doi: 10.12000/JR15069
    [13] 郝天铎崔琛龚阳孙从易 . 基于凸优化方法的认知雷达低峰均比波形设计. 雷达学报, 2018, 7(4): 498-506. doi: 10.12000/JR18002
    [14] 吴敏张磊刘松杨邢孟道 . OFDM-ISAR的稀疏优化成像与运动补偿. 雷达学报, 2016, 5(1): 72-81. doi: 10.12000/JR16017
    [15] 王文钦程胜娟邵怀宗 . 基于稀疏矩阵和相关函数联合优化的MIMO-OFDM线性调频波形复用设计与实现方法. 雷达学报, 2015, 4(1): 1-10. doi: 10.12000/JR14148
    [16] 滑文强王爽郭岩河谢雯 . 基于邻域最小生成树的半监督极化SAR图像分类方法. 雷达学报, 2019, (): 1-13. doi: 10.12000/JR18104
    [17] 赵勇胜赵拥军赵闯 . 联合角度和时差的单站无源相干定位加权最小二乘算法. 雷达学报, 2016, 5(3): 302-311. doi: 10.12000/JR15133
    [18] 陈世超罗丰胡冲聂学雅 . 基于多普勒谱非广延熵的海面目标检测方法. 雷达学报, 2019, (): 1-11. doi: 10.12000/JR19012
    [19] 熊丁丁崔国龙孔令讲杨晓波 . 基于互相关熵的非高斯背景下微动参数估计方法. 雷达学报, 2017, 6(3): 300-308. doi: 10.12000/JR17007
    [20] 葛建军李春霞 . 一种基于信息熵的雷达动态自适应选择跟踪方法. 雷达学报, 2017, 6(6): 587-593. doi: 10.12000/JR17081
  • 加载中
图(17)表(3)
计量
  • 文章访问数:  176
  • HTML浏览量:  173
  • PDF下载量:  77
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-01-02
  • 录用日期:  2019-03-16
  • 网络出版日期:  1970-01-01
  • 刊出日期:  2019-04-28

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

    通讯作者: 邢孟道, xmd@xidian.edu.cn
    作者简介: 符吉祥(1992–),男,安徽蚌埠人,博士生,研究方向为ISAR运动补偿和成像。E-mail: jixiang_xd@126.com;邢孟道(1975–),男,浙江绍兴人,博士,教授,西安电子科技大学前沿交叉研究院副院长,研究方向为雷达成像、动目标检测、目标识别。E-mail: xmd@xidian.edu.cn;徐   丹(1990–),女,陕西咸阳人,博士生,研究方向为3维ISAR图像重构、电磁成像。 E-mail: 1045759961@qq.com;王安乐(1988–),男,山东临沂人,博士,中国人民解放军空军预警学院讲师,研究方向为微波光子雷达、光生微波基准信号等。E-mail: anlehit@163.com
  • ①. 西安电子科技大学雷达信号处理国家重点实验室   西安   710071
  • ②. 西安电子科技大学信息感知技术协同创新中心   西安   710071
  • ③. 空军预警学院   武汉   430019
基金项目:  国家杰出青年自然基金(61825105)

摘要: 机翼在飞机运动不平稳状态下存在振动,在微波光子超高分辨雷达观测下,这种振动会造成机翼难以聚焦,针对这一难题,该文提出一种基于微波光子超高分辨雷达的机翼振动参数估计方法。首先通过粗成像将机身和机翼回波进行分离,再通过对机身成像和定标结果估计雷达视线角(LOS)。然后对机翼进行子孔径序列成像,提取散射点的距离和多普勒变化曲线,再联合雷达LOS以及距离和多普勒曲线对振动参数进行粗估计,最后通过修正的极坐标格式算法(MPFA)以及构造最小熵优化函数对振动参数进行精估计。该文首次提出了修正的极坐标格式算法,其能够对复杂运动的目标进行距离和方位向的解耦,如震动的机翼和摇摆的舰船等。仿真和实测数据的处理结果验证了该方法的有效性和实用性。

English Abstract

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

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

      图  1  微波光子实测飞机目标1维距离像

      Figure 1.  One dimensional range profile of airplane measured by microwave photonic radar

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

      图  2  飞机运动几何模型

      Figure 2.  Geometry model of airplane

      $ \vartheta \left( {{t_\rm{m}}} \right){\rm{ = }}{\vartheta _0}\sin \left( {2{{π}} {f_\rm{v} }{t_\rm{m}} + {\varphi _0}} \right) $

      其中,${\vartheta _0}$表示振动角,${f_\rm{v} }$表示振动频率,${\varphi _0}$表示振动初相,${t_\rm{m} }$表示方位慢时间。在飞机沿着图2中曲线轨迹运动的时候,飞机整体还有相对于雷达的平动和转动。对于机身上的散射点${\rm{p} _1}$,因为其不存在振动,所以其相对于雷达的运动可以表示为

      $ {R_{{\rm{p} _1}}}\left( {{t_\rm{m} }} \right) = {R^\rm{T} }\left( {{t_\rm{m} }} \right) + R_{{\rm{p} _1}}^{\rm{R}} \left( {{t_\rm{m} }} \right) $

      其中,${R^\rm{T}}\left( {{t_\rm{m} }} \right)$$R_{{\rm{p} _1}}^{\rm{R}} \left( {{t_\rm{m} }} \right)$分别表示散射点${\rm{p} _1}$的平动和转动分量。平动分量对飞机上所有的散射点来说都是相同的,可以将其表示为多项式的形式

      $ {R^\rm{T} }\left( {{t_\rm{m} }} \right) = {R_0} + {v_\rm{r}}{t_\rm{m}}{\rm{ + }}\frac{1}{2}{a_\rm{r}}t_{\rm{m}}^2 +··· $

      其中,${R_0}$表示飞机旋转中心$O\,$到雷达的初始距离,${v_\rm{r}}$${a_\rm{r}}$分别表示飞机相对于雷达运动的径向速度和径向加速度。转动分量是随散射点位置发生空变的,其形式为(平面波近似)

      $ \begin{align} R_{{{\rm{p}}_1}}^{\rm{R}}\left( {{t_{\rm{m}}}} \right){\rm{ = }}&{\left[ {\begin{array}{*{20}{c}} {{x_{{{\rm{p}}_1}}}}\\ {{y_{{{\rm{p}}_1}}}}\\ 0 \end{array}} \right]^{\rm T}}\!\!\underbrace {\left[ {\begin{array}{*{20}{c}} {\cos \theta \left( {{t_{\rm{m}}}} \right)}&{\sin \theta \left( {{t_{\rm{m}}}} \right)}&0\\ { - \sin \theta \left( {{t_{\rm{m}}}} \right)}&{\cos \theta \left( {{t_{\rm{m}}}} \right)}&0\\ 0&0&1 \end{array}} \right]}_{{旋转矩阵}}\\ {\rm{}}& \cdot\underbrace {\left[ {\begin{array}{*{20}{c}} {\sin {\beta _0}\sin {\phi _0}}\\ {\cos {\beta _0}\sin {\phi _0}}\\ {\cos {\phi _0}} \end{array}} \right]}_{{投影向量}} \end{align} $

      其中,$\left( {{x_{{\rm{p}_1}}},{y_{{\rm{p}_1}}}} \right)$表示散射点${\rm{p}_1}$参考时刻的坐标位置,$\theta $表示飞机相对于雷达转动的角度,其是慢时间的函数,即$\theta \left( {{t_\rm{m} }} \right) = \omega {t_\rm{m} }$, $\omega $表示转速。对于机翼上的散射点,其不仅包含以上的运动,还包含机翼振动对其的运动调制。通过运动分解可以得到机翼上散射点${\rm{p}_2}$各个时刻的坐标位置如式(5)

      $ \begin{aligned} \left[ {\begin{array}{*{20}{c}} {{x_{{{\rm{p}}_2}}}\left( {{t_{\rm{m}}}} \right)}\\ {{y_{{{\rm{p}}_2}}}\left( {{t_{\rm{m}}}} \right)}\\ {{z_{{{\rm{p}}_2}}}\left( {{t_{\rm{m}}}} \right)} \end{array}} \right] =& \underbrace {\left[ {\begin{array}{*{20}{c}} {\cos \theta \left( {{t_{\rm{m}}}} \right)}&{\sin \theta \left( {{t_{\rm{m}}}} \right)}&0\\ { - \sin \theta \left( {{t_{\rm{m}}}} \right)}&{\cos \theta \left( {{t_{\rm{m}}}} \right)}&0\\ 0&0&1 \end{array}} \right]}_{{旋转矩阵}} \underbrace {\left[ {\begin{array}{*{20}{c}} 1&0&0\\ 0&{\cos \vartheta \left( {{t_{\rm{m}}}} \right)}&{\sin \vartheta \left( {{t_{\rm{m}}}} \right)}\\ 0&{ - \sin \vartheta \left( {{t_{\rm{m}}}} \right)}&{\cos \vartheta \left( {{t_{\rm{m}}}} \right)} \end{array}} \right]}_{{振动矩阵}}\left[ {\begin{array}{*{20}{c}} {{x_{{{\rm{p}}_2}}}}\\ {{y_{{{\rm{p}}_2}}}}\\ 0 \end{array}} \right] \end{aligned} $

      其中,$\left( {{x_{{\rm{p} _2}}},{y_{{\rm{p}_2}}}} \right)$表示散射点${\rm{p} _2}$参考时刻的坐标位置。所以散射点${\rm{p} _2}$到雷达的斜距变化为

      $ \begin{align} {R_{{\rm{p}_2}}}\left( {{t_{\rm{m}}}} \right){\rm{ = }}&{{{R}}^\rm{T} }\left( {{t_\rm{m} }} \right) + R_{{\rm{p} _2}}^{\rm {R}} \left( {{t_{\rm{m}}}} \right) ={{{R}}^{\rm{T}} }\left( {{t_{\rm{m}}}} \right) + {\left[ {\begin{array}{*{20}{c}} {{x_{{\rm{p} _2}}}\left( {{t_{\rm{m}}}} \right)} \\ {{y_{{\rm{p} _2}}}\left( {{t_{\rm{m}}}} \right)} \\ {{z_{{\rm{p}_2}}}\left( {{t_{\rm{m}}}} \right)} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {\sin {\beta _0}\sin {\phi _0}} \\ {\cos {\beta _0}\sin {\phi _0}} \\ {\cos {\phi _0}} \end{array}} \right] \end{align} $

      其中,$R_{{\rm{p} _2}}^{\rm{R}} \left( {{t_{\rm{m}}}} \right)$表示散射点${\rm{p} _2}$的转动分量。

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

      $ \begin{align} s\left( {{k_{\rm{r}}},{t_{\rm{m}}}} \right) =& \underbrace {\sum\limits_{p = 1}^{{N_1}} {{A_p}\exp \left( { - {\rm{j}}{k_{\rm{r}}}{R_p}\left( {{t_{\rm{m}}}} \right)} \right)} }_{{机身}} {\rm{ + }}\underbrace {\sum\limits_{q = 1}^{{N_2}} {{A_q}\exp \left( { - {\rm{j}}{k_{\rm{r}}}{R_q}\left( {{t_{\rm{m}}}} \right)} \right)} }_{{机翼}} \end{align} $

      其中,${k_{\rm{r}}} = 4{{π}} \left( {\gamma \hat t + {f_{\rm{c}}}} \right)/{\rm{c}}$表示径向波数,$\hat t$表示快时间,$\rm{c} $表示电磁波在自由空间中的速度,${N_1}$为机身上散射点个数,${N_2}$为振动的机翼上的散射点个数。因为平动分量对于所有的散射点都相同,其可以通过平动补偿方法[11,12]进行补偿。补偿之后的模型转换成带振动的转台模型,如图3(a)所示。平动补偿之后的回波表达式可以表示为

      图  3  转台模型

      Figure 3.  Turntable model

      $ \begin{align} s\left( {{k_{\rm{r}}},{t_{\rm{m}}}} \right) =& \underbrace {\sum\limits_{p = 1}^{{N_1}} {{A_p}\exp \bigr( { - {\rm{j}}{k_{\rm{r}}}\left( {{x_p}\sin \left( {{\beta _0} - \theta \left( {{t_{\rm{m}}}} \right)} \right) + {y_p}\cos \left( {{\beta _0} - \theta \left( {{t_{\rm{m}}}} \right)} \right)} \right)\sin {\phi _0}}\bigr)} }_{{机身}}\\ {\rm{}}&{\rm{ + }}\!\!\underbrace {\sum\limits_{q = 1}^{{N_2}} {{A_q}\exp \bigr( { - {\rm{j}}{k_{\rm{r}}}\left( {{x_q}\sin \left( {{\beta _0} - \theta \left( {{t_{\rm{m}}}} \right)} \right)\sin {\phi _0} + {y_q}\left( {\sin {\phi _0}\cos \vartheta \cos \left( {{\beta _0} - \theta \left( {{t_{\rm{m}}}} \right)} \right) - \cos {\phi _0}\sin \vartheta } \right)} \right)} \bigr)} }_{{机翼}}\quad\quad \end{align} $

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

      图  4  飞机粗成像结果

      Figure 4.  Coarse imaging results of airplane

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

      $ \begin{align} {s_{\rm{body} }}\left( {r,{t_\rm{m} }} \right) =&\sum\limits_{p = 1}^{{N_1}} {A_p}{\rm{sinc}} \left( {r - \left( {x_p}\cos {\beta _0} + {y_p}\sin {\beta _0} \right)\!\sin {\phi _0}} \right)\\ {\rm{}}&\cdot\exp \left( { - {\rm{j}}{k_{\rm{rc} }}\left( \left( {{y_p}\cos {\beta _0} - {x_p}\sin {\beta _0}} \right)\omega {t_\rm{m} } + \left( {{x_p}\cos {\beta _0} + {y_p}\sin {\beta _0}} \right)\left( {1 - \frac{{{\omega ^2}t_{\rm{m}}^2}}{2}} \right) \right)\sin {\phi _0}} \right) \quad \end{align} $

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

      $ {H_{\rm{c}}}\left( {r,{t_{\rm{m}}}} \right) = \exp \left( { - {\rm{j}}{k_{{\rm{rc}}}}r\frac{{{\omega ^2}t_{\rm{m}}^2}}{2}} \right) $

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

      $ \begin{align} {s_{\rm{body} }}\left( {r,{f_\rm{d} }} \right) =& \sum\limits_{p = 1}^{{N_1}} {A_p}\;{\rm{sinc}} \bigr( r - \left( {{x_p}\cos {\beta _0} + {y_p}\sin {\beta _0}} \right) \\ \\ {\rm{}}& \cdot\sin {\phi _0} \bigr){\rm{sinc}} \bigr( {f_\rm{d} } - \omega \left( {y_p}\cos {\beta _0}\right. \\ {\rm{}}& \left.- {x_p}\sin {\beta _0} \right)\sin {\phi _0}\bigr) \end{align} $

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

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

      $\left. \begin{aligned} r_q^i \approx &\Delta {x_q}\left( {\sin {\beta _0} - \cos {\beta _0}\omega t_{\rm{c}}^i} \right)\sin {\phi _0} \\ {\rm{}}&+ \Delta {y_q}\sin {\phi _0}\left( {\cos {\beta _0} + \sin {\beta _0}\omega t_{\rm{c}}^i} \right) \\ f_{{\rm{d}} q}^i \approx &- \frac{{{k_{\rm{rc} }}}}{{2{{π}} }}\Bigr( - \Delta {x_q}\sin {\phi _0}\cos {\beta _0}\omega \\ {\rm{}}& - \Delta {y_q}\left( \cos {\phi _0}2{{π}} {f_{\rm{v}}}{\vartheta _0}\cos \left( {2{{π}} {f_\rm{v} }t_{\rm{c}}^i + {\varphi _0}} \right)\right. \\ {\rm{}} & + \sin {\phi _0}\sin {\beta _0}\omega \bigr)\Bigr) \end{aligned} \right\}$

      其中,$r_q^i$$f_{{\rm{d}} q}^i$分别表示第$i$个子孔径图像中第$q$个散射点的距离和多普勒,$t_{\rm c}^i$表示第$i$个子孔径的中心时刻,$\left( {\Delta {x_q},\Delta {y_q}} \right)$表示第$q\,$个散射点相对于等效旋转中心$O'$的坐标。观察式(12)可以发现,散射点的距离与振动参数无关,而散射点的多普勒随着慢时间余弦变化,频率为${f_\rm{v} }$,幅度与散射点的坐标$\Delta {y_q}$以及机翼振动角${\vartheta _0}$成正比。所以通过提取散射点的多普勒变化曲线,测量相邻波峰之间的时间${T_\rm{v}}$,即可得到振动频率的估计,即

      $ {f_{\rm{v}}} = 1/{T_{\rm{v}}} $

      假设测量得到的散射点多普勒波峰和波谷之间的差异为$\Delta f_{{\rm{d}} q}^i$,振动角可由式(14)计算得到

      $ {\vartheta _0} = \frac{{\Delta f_{{\rm{d}} q}^i}}{{2{k_{\rm{rc}}}{f_\rm{v}}\left| {\Delta {y_q}} \right|\cos {\phi _0}}} $

      其中,$\Delta {y_q}$的值可由以下方程组解得

      $ \left. \begin{aligned} r_q^k \approx &\Delta {x_q}\Bigr( {\sin {\beta _0} - \cos {\beta _0}\omega t_{\rm{c}}^k} \Bigr)\sin {\phi _0} \\ {\rm{}}&+ \Delta {y_q}\sin {\phi _0}\Bigr( {\cos {\beta _0} + \sin {\beta _0}\omega t_{\rm{c}}^k} \Bigr) \\ r_q^j \approx & \Delta {x_q}\Bigr( {\sin {\beta _0} - \cos {\beta _0}\omega t_{\rm{c}}^j} \Bigr)\sin {\phi _0} \\ {\rm{}}&+ \Delta {y_q}\sin {\phi _0}\Bigr( {\cos {\beta _0} + \sin {\beta _0}\omega t_{\rm{c}}^j} \Bigr) \end{aligned} \right\} $

      其中,$r_q^j$表示第$j$个子孔径图像中第$q$个散射点的距离,$r_q^k$表示第$k$个子孔径图像中第$q\,$个散射点的距离,$t_{\rm{c}}^j$$t_{\rm{c}}^k$表示第$j$$k$个子孔径的中心时刻。振动初相${\varphi _0}$可由式(16)计算得到

      $ {\varphi _0} = \arccos \left( {\frac{{f_{{\rm{d}} 0q}^i}}{{{k_{\rm{rc} }}\Delta {y_q}\cos {\phi _0}{f_\rm{v} }{\vartheta _0}}}} \right) - 2{{π}} {f_\rm{v} }t_{\rm{c}}^i $

      其中,$f_{{\rm{d}} 0q}^i{\rm{ = }}f_{{\rm{d}} q}^i - {\rm{mean}} \left( {{{{f}}_{{\rm{d}} q}}} \right)$,表示去除多普勒中心后的散射点多普勒频率,$\rm{mean} \left( \cdot \right)$表示取均值操作,${{{f}}_{{\rm{d}} q}}$表示第$q$个散射点的多普勒向量。

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

      $ \left. \begin{aligned} {k_x} =& {k_\rm{r} }\sin \left( {{\beta _0} - \theta } \right)\sin \phi \\ {k_y} =& {k_\rm{r}}\left( {\cos {\phi _0}\sin \vartheta + \sin {\phi _0}\cos \vartheta \cos \left( {{\beta _0} - \theta } \right)} \right) \\ \end{aligned} \right\} $

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

      $ {s_{\rm{wing} }}\left( {{k_x},{k_y}} \right) = \sum\limits_{q = 1}^{{N_2}} {{A_q}\exp \left( { - {\rm{j}}\left( {{k_x}{x_q} + {k_y}{y_q}} \right)} \right)} $

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

      为了能够进一步提高振动参数${{Θ}} {\rm{ = }}\left[ {{\vartheta _0},{\varphi _0}} \right]$的估计精度,构造基于图像熵值的代价函数,通过优化参数使得图像熵值达到最小来对振动参数进行精估计。图像熵能够衡量图像的聚焦质量,同一图像熵值越小,其聚焦性能越好,图像熵定义如式(19)

      $ {\rm{E}} \left[ {{s}} \right] = - \sum\limits_{i = 1}^{MN} {\frac{{{{\left| {{s_i}} \right|}^2}}}{{{P_\rm{s} }}}\ln \frac{{{{\left| {{s_i}} \right|}^2}}}{{{P_\rm{s} }}}} $

      其中,$\rm{E} \left[ \cdot \right]$表示计算熵值操作,${{s}}$表示图像矩阵,矩阵大小为$M \!\times\! N$, ${s_i}$为其第$i$个元素。${P_\rm{s} } \!=\! \displaystyle\sum\nolimits_{i = 1}^{MN} {{{\left| {{s_i}} \right|}^2}} $表示图像的总能量。所以构造的最小熵优化函数如式(20)

      $ \hat {{Θ}} = \arg \!\mathop {\min }\limits_{{Θ}} {\rm{E}} \left( {\Omega \left( {{{{s}}_{\rm{wing} }}\left( {{k_\rm{r} },{t_\rm{m} }} \right),{{Θ}} } \right)} \right) $

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

      图  5  算法流程图

      Figure 5.  Flow chart of the proposed algorithm

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

      信号带宽(GHz)脉冲宽度(μs)载频(GHz)PRF (Hz)采样率(Hz)脉冲数振动角(°)振动频率(Hz)初相(°)
      1015035666.665008192110

      表 1  仿真雷达参数

      Table 1.  Simulation radar parameters

      图  6  飞机散射点模型

      Figure 6.  Scatter model of airplane

      图7是机身聚焦成像定标结果,从图7中可以看出机身主轴和Y轴存在夹角,通过之前的分析可知,其夹角大小和雷达视线的方位角是相等的,通过机头标记点A的坐标位置,可得雷达视线的方位角为${\beta _0} = \arctan \left( {31.39/5.5\;} \right) = {80.0618^ \circ}$,通过测量飞机宽度BC与实际中机身宽度之间的比例,可以得到对俯仰角的估计,即${\phi _0} = \arcsin \left( {\left| {\rm BC}\, \right|/{W_0}} \right) $$ = {85.017^ \circ }$,其中${W_0}$为飞机的真实宽度。

      图  7  机身聚焦成像定标结果

      Figure 7.  Focused and scaled result of airplane body

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

      图  9  提取的不同散射点的多普勒变化曲线

      Figure 9.  Extracted Doppler curves of different scatterers

      图  10  提取的不同散射点的距离变化曲线

      Figure 10.  Extracted range curves of different scatterers

      图  11  图像熵值优化函数代价曲面

      Figure 11.  Cost surface of entropy optimization function

      图  12  机翼区域精细化处理前后结果对比

      Figure 12.  Before and after accurate processing of airplane wings

      振动频率(Hz)振动角(°)振动初相(°)
      实际值估计值实际值估计值实际值估计值
      10.9910.980–0.18

      表 2  仿真数据机翼振动参数估计结果

      Table 2.  Vibration parameters estimation result of simulation data

      图  8  不同时刻机翼子孔径成像结果

      Figure 8.  Sub-aperture imaging results of airplane wing at different time

      为了进一步验证本文方法的有效性,下面对实测数据进行处理,实测机身成像定标结果如图13所示,通过测量机身轴线与Y轴之间的夹角,可以得到方位角的估计,即${\hat \beta _0} = \arctan \left( {6.59/0.82} \right){\rm{ = }}{82.9^ \circ }$。观测飞机的机型为A320,通过公开资料可知其机身宽度为3.95 m,通过测量定标图像中机身宽度${\rm{AB}}{\rm{ = }}3.25\;\rm{m}$,可以得到俯仰角的估计${\hat \phi _0} = {55.36^\circ}$图14是振动较为明显的机翼不同子孔径的距离多普勒成像结果,通过提取图中用圆圈标记的散射点1的距离和多普勒变化,如图15所示。最小熵优化函数的代价平面如图16所示,振动参数估计结果如表3所示。机翼区域精细化处理前后的对比图像如图17所示,从图中用圆圈标记的散射点的聚焦情况可以看出,使用本文方法对机翼区域精细化处理之后,散射点的聚焦质量有所改善。

      图  13  实测数据机身聚焦成像定标结果

      Figure 13.  Focused and scaled airplane body result of the measured data

      图  14  不同子孔径机翼成像结果序列图

      Figure 14.  Airplane wing imaging results of different sub-apertures

      图  15  提取的散射点的距离和多普勒变化曲线

      Figure 15.  Range and Doppler curves of extracted scatter

      图  16  图像熵值优化函数代价曲面

      Figure 16.  Cost surface of image entropy optimization function

      图  17  机翼区域精细化处理前后结果对比

      Figure 17.  Before and after accurate processing of airplane wings

      振动频率(Hz)振动角(°)振动初相(°)
      2.480.0712.93

      表 3  实测数据机翼振动参数估计结果

      Table 3.  Vibration parameters estimation result of measured data

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

参考文献 (15)

目录

    /

    返回文章
    返回