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

符吉祥 邢孟道 徐丹 王安乐

符吉祥, 邢孟道, 徐丹, 等. 一种基于微波光子超高分辨雷达机翼振动参数估计方法[J]. 雷达学报, 2019, 8(2): 232–242. doi:  10.12000/JR19001
引用本文: 符吉祥, 邢孟道, 徐丹, 等. 一种基于微波光子超高分辨雷达机翼振动参数估计方法[J]. 雷达学报, 2019, 8(2): 232–242. doi:  10.12000/JR19001
FU Jixiang, XING Mengdao, XU Dan, et al. Vibration-parameters estimation method for airplane wings based on microwave-photonics ultrahigh-resolution radar[J]. Journal of Radars, 2019, 8(2): 232–242. doi:  10.12000/JR19001
Citation: FU Jixiang, XING Mengdao, XU Dan, et al. Vibration-parameters estimation method for airplane wings based on microwave-photonics ultrahigh-resolution radar[J]. Journal of Radars, 2019, 8(2): 232–242. doi:  10.12000/JR19001

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

doi: 10.12000/JR19001
基金项目: 国家杰出青年自然基金(61825105)
详细信息
    作者简介:

    符吉祥(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

  • 中图分类号: TN957

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

Funds: The National Science Fund for Distinguished Young Scholars (61825105)
More Information
图(17) / 表 (3)
计量
  • 文章访问数:  1500
  • HTML全文浏览量:  334
  • PDF下载量:  142
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-01-02
  • 修回日期:  2019-03-16
  • 网络出版日期:  2019-04-12
  • 刊出日期:  2019-04-01

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

doi: 10.12000/JR19001
    基金项目:  国家杰出青年自然基金(61825105)
    作者简介:

    符吉祥(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
  • 中图分类号: TN957

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

English Abstract

符吉祥, 邢孟道, 徐丹, 等. 一种基于微波光子超高分辨雷达机翼振动参数估计方法[J]. 雷达学报, 2019, 8(2): 232–242. doi:  10.12000/JR19001
引用本文: 符吉祥, 邢孟道, 徐丹, 等. 一种基于微波光子超高分辨雷达机翼振动参数估计方法[J]. 雷达学报, 2019, 8(2): 232–242. doi:  10.12000/JR19001
FU Jixiang, XING Mengdao, XU Dan, et al. Vibration-parameters estimation method for airplane wings based on microwave-photonics ultrahigh-resolution radar[J]. Journal of Radars, 2019, 8(2): 232–242. doi:  10.12000/JR19001
Citation: FU Jixiang, XING Mengdao, XU Dan, et al. Vibration-parameters estimation method for airplane wings based on microwave-photonics ultrahigh-resolution radar[J]. Journal of Radars, 2019, 8(2): 232–242. doi:  10.12000/JR19001
    • 得益于全天时全天候和超远距离探测的优势,雷达被广泛应用到军事和民事的应用中,如战场监视和遥感测绘等。逆合成孔径雷达,因为其能对非合作目标进行探测和成像,获得目标的形状、尺寸、姿态和运动信息,进而可以对非合作目标进行识别,在国防安全领域有着重要应用,受到各个国家的高度重视和广泛研究。随着雷达技术和需求的不断提高,雷达向高分辨、多功能和多任务方向发展。然而,传统微波雷达受限于微波器件的“电子瓶颈”[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平面上的投影向量${\text{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{\text{π}} {f_\rm{v} }{t_\rm{m}} + {\varphi _0}} \right) $$ (1)

      其中,${\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) $$ (2)

      其中,${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 +··· $$ (3)

      其中,${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]}_{\text{旋转矩阵}}\\ {\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]}_{\text{投影向量}} \end{align} $$ (4)

      其中,$\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]}_{\text{旋转矩阵}} \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]}_{\text{振动矩阵}}\left[ {\begin{array}{*{20}{c}} {{x_{{{\rm{p}}_2}}}}\\ {{y_{{{\rm{p}}_2}}}}\\ 0 \end{array}} \right] \end{aligned} $$ (5)

      其中,$\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} $$ (6)

      其中,$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)} }_{\text{机身}} {\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)} }_{\text{机翼}} \end{align} $$ (7)

      其中,${k_{\rm{r}}} = 4{\text{π}} \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)} }_{\text{机身}}\\ {\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)} }_{\text{机翼}}\quad\quad \end{align} $$ (8)
    • 首先通过粗成像对机翼机身回波进行分离。因为机翼存在周期性振动而机身只存在平稳的运动,所以在粗成像结果中,机翼和机身上的散射点散焦情况存在明显的差异,即机翼散射点散焦成线性分布而机身散射点散焦呈现块状分布。另外,从粗成像的轮廓中也可以轻易地将机翼和机身进行区分。如图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} $$ (9)

      其中,${\rm{sinc}}\left( x \right) = \sin \left( {{\text{π}} x} \right)/\left( {{\text{π}} x} \right)$, ${k_{\rm{rc}}}{\rm{ = }}{{4{\text{π}} {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) $$ (10)

      其中,$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)

      从式(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{\text{π}} }}\Bigr( - \Delta {x_q}\sin {\phi _0}\cos {\beta _0}\omega \\ {\rm{}}& - \Delta {y_q}\left( \cos {\phi _0}2{\text{π}} {f_{\rm{v}}}{\vartheta _0}\cos \left( {2{\text{π}} {f_\rm{v} }t_{\rm{c}}^i + {\varphi _0}} \right)\right. \\ {\rm{}} & + \sin {\phi _0}\sin {\beta _0}\omega \bigr)\Bigr) \end{aligned} \right\}$$ (12)

      其中,$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}}} $$ (13)

      假设测量得到的散射点多普勒波峰和波谷之间的差异为$\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}}} $$ (14)

      其中,$\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\} $$ (15)

      其中,$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{\text{π}} {f_\rm{v} }t_{\rm{c}}^i $$ (16)

      其中,$f_{{\rm{d}} 0q}^i{\rm{ = }}f_{{\rm{d}} q}^i - {\rm{mean}} \left( {{{\text{f}}_{{\rm{d}} q}}} \right)$,表示去除多普勒中心后的散射点多普勒频率,$\rm{mean} \left( \cdot \right)$表示取均值操作,${{\text{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)

      将式(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)

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

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

      $$ {\rm{E}} \left[ {\text{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} }}}} $$ (19)

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

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

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

      图  5  算法流程图

      Figure 5.  Flow chart of the proposed algorithm

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

      表 1  仿真雷达参数

      Table 1.  Simulation radar parameters

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

      图  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

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

      Table 2.  Vibration parameters estimation result of simulation data

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

      图  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

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

      Table 3.  Vibration parameters estimation result of measured data

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

参考文献 (15)

目录

    /

    返回文章
    返回