基于数据同化的飞机尾流行为预测

沈淳 李健兵 高航 陈柏纬 韩启光 王雪松

沈淳, 李健兵, 高航, 等. 基于数据同化的飞机尾流行为预测[J]. 雷达学报, 待出版. doi:  10.12000/JR21007
引用本文: 沈淳, 李健兵, 高航, 等. 基于数据同化的飞机尾流行为预测[J]. 雷达学报, 待出版. doi:  10.12000/JR21007
SHEN Chun, LI Jianbing, GAO Hang, et al. Aircraft wake vortex behavior prediction based on data assimilation[J]. Journal of Radars, in press. doi:  10.12000/JR21007
Citation: SHEN Chun, LI Jianbing, GAO Hang, et al. Aircraft wake vortex behavior prediction based on data assimilation[J]. Journal of Radars, in press. doi:  10.12000/JR21007

基于数据同化的飞机尾流行为预测

doi: 10.12000/JR21007
基金项目: 国家自然科学基金(61490649, 61771479, 61625108),湖南省杰出青年基金(2018JJ1030)
详细信息
    作者简介:

    沈 淳(1985–),男,福建漳州人,博士生,工程师,研究方向为空间信息获取与处理

    李健兵(1979–),男,湖南邵东人,博士,国防科技大学电子科学学院教授,主要研究方向为新体制雷达、雷达信号处理

    王雪松(1972–),男,内蒙古人,博士,国防科技大学电子科学学院教授,主要研究方向为极化信息处理、新体制雷达技术、电子对抗

    通讯作者:

    李健兵 jianbingli@nudt.edu.cn

  • 责任主编:夏海云 Corresponding Editor: XIA Haiyun
  • 中图分类号: TN955+.1

Aircraft Wake Vortex Behavior Prediction Based on Data Assimilation

Funds: The National Natural Science Foundation of China (61490649, 61771479, 61625108), Hunan Natural Science Foundation for Distinguished Young Scholars (2018JJ1030)
More Information
  • 摘要: 飞机尾流是飞机飞行时在其后方产生的一对反向旋转的强烈湍流,对后续飞机飞行以及机场安全起降影响极大,其演化趋势的预测已成为空中交通安全管制的瓶颈,亟需发展基于实时探测数据的飞机尾流行为预测技术。在雷达探测反演得到的尾流涡心位置和速度环量等特征参数基础上,开展飞机尾流行为预测分析,能够预知飞机尾流危害区域,为机场安全起降动态间隔标准制定提供技术支撑。该文结合风场线性切变和最小二乘拟合方法构建了参数化尾流行为预测模型,解决了经典尾流预测模型气象环境参数未随时间演化实时调整的问题。该文根据复杂风场非线性演化特点,设计了基于无迹卡尔曼滤波的数据同化模型,利用雷达探测数据对尾流行为预测进行实时修正。数值仿真验证和实测数据验证结果表明,基于数据同化的飞机尾流行为预测方法能够根据实时探测数据对尾流行为预测轨迹进行修正,得到更加贴近实测的飞机尾流行为预测轨迹。
  • 图  1  飞机尾流形成示意图

    Figure  1.  Illustration of aircraft wake vortex

    图  2  激光雷达探测飞机尾流场景设置

    Figure  2.  Geometry configuration for Lidar detection of wake vortex

    图  3  飞机尾流多普勒速度RHI回波分布图

    Figure  3.  Doppler velocity distribution of wake vortex in an RHI

    图  4  飞机尾流行为预测流程图

    Figure  4.  Wake vortex behavior prediction flow chart

    图  5  估计背景风场的非尾流区域

    Figure  5.  Regions free of wake vortex that was used to estimate the background wind

    图  6  切向速度与多普勒速度的几何关系

    Figure  6.  Relationship between the tangential velocity and Doppler velocity

    图  7  数据融合方法流程

    Figure  7.  Flowchart of data fusion method

    图  8  飞机尾流行为预测方法对比(水平方向轨迹)

    Figure  8.  Comparison between different wake vortex behavior prediction (Horizontal trajectories)

    图  9  飞机尾流行为预测方法对比(垂直方向轨迹)

    Figure  9.  Comparison between different wake vortex behavior prediction (Vertical trajectories)

    图  10  飞机尾流行为预测方法对比(速度环量)

    Figure  10.  Comparison between different wake vortex behavior prediction (Circulation)

    图  11  香港机场北跑道激光雷达探测场景设置2014

    Figure  11.  Geometry setup of the observation in north runway of Hong Kong international airport, 2014

    图  12  飞机尾流行为预测方法对比

    Figure  12.  Comparison between different wake vortex behavior prediction

    图  13  香港机场南跑道激光雷达探测场景设置2018

    Figure  13.  Geometry setup of the observation in south runway of Hong Kong international airport, 2018

    图  14  飞机尾流行为预测方法对比

    Figure  14.  Comparison between different wake vortex behavior prediction

    表  1  激光雷达探测参数设置

    Table  1.   Detection parameters of the Lidar

    主要参数量值
    雷达波长1.55 μm
    脉冲宽度170 ns
    时间窗长度120 ns
    距离门宽度21 m
    采样率50 MHz
    脉冲积累数1500
    信号噪声比–5 dB
    扫描速度2°/s
    扫描范围0~15°
    下载: 导出CSV

    表  2  飞机尾流预测位置相对误差

    Table  2.   Relative error of predict trajectories

    主要参数DS method (%)DA method (%)
    横向风
    $V_{\rm{c}}^0$ (${\rm{m/s}}$)
    切变率$\beta $纵向风
    $V_{\rm{c}}^{\rm{y}}$ (${\rm{m/s}}$)
    ${E_{r1}}$${E_{r2}}$${E_{r1}}$${E_{r2}}$
    –70.05–0.317.7419.631.121.35
    –100.05–0.326.8328.522.642.41
    –50.01–0.33.423.141.981.76
    –50.10–0.34.164.712.862.17
    –50.05-0.53.593.211.231.01
    –50.05–0.86.876.521.361.71
    下载: 导出CSV

    表  3  飞机尾流速度环量相对误差

    Table  3.   Relative error of wake vortex circulation

    主要参数
    EDR (${{\rm{m}}^2}/{{\rm{s}}^3}$)
    DS method (%)DA method (%)
    ${E_{r1}}$${E_{r2}}$${E_{r1}}$${E_{r2}}$
    0.013.243.581.121.31
    0.033.263.751.251.02
    0.083.042.951.171.16
    0.103.413.251.081.29
    下载: 导出CSV

    表  4  香港机场激光雷达探测参数设置

    Table  4.   Detection parameters of the Lidar in Hong Kong field campaigns

    主要参数量值
    雷达波长(μm)1.54
    距离门宽度(m)25
    脉冲宽度(ns)200
    脉冲重复频率(kHz)20
    探测距离(m)50~6000
    下载: 导出CSV
  • [1] 李健兵, 高航, 王涛, 等. 飞机尾流的散射特性与探测技术综述[J]. 雷达学报, 2017, 6(6): 660–672. doi:  10.12000/JR17068

    LI Jianbing, GAO Hang, WANG Tao, et al. A survey of the scattering characteristics and detection of aircraft wake vortices[J]. Journal of Radars, 2017, 6(6): 660–672. doi:  10.12000/JR17068
    [2] GERZ T, HOLZÄPFEL F, BRYANT W, et al. Research towards a wake-vortex advisory system for optimal aircraft spacing[J]. Comptes Rendus Physique, 2005, 6(4/5): 501–523. doi:  10.1016/j.crhy.2005.06.002
    [3] THOBOIS Ludovic and CARIOU Jean-Pierre. Next generation scanning LIDAR systems for optimizing wake turbulence separation minima[J]. Journal of Radars, 2017, 6(6): 689–698. doi: 10.12000/JR17056.
    [4] HON Kaikwong and CHAN Pakwai. Aircraft wake fortex observations in Hong Kong[J]. Journal of Radars, 2017, 6(6): 709–718. doi: 10.12000/JR17072.
    [5] Vortex State-of-the-Art & Research Needs. Project report under EC contract 2134622015[R]. 2015. doi: 10.17874/BFAEB7154B0.
    [6] CHENG J, TITTSWORTH J, GALLO W, et al. The development of wake turbulence recategorization in the United States[C]. 8th AIAA Atmospheric and Space Environments Conference, Washington, USA, 2016: 1–12.
    [7] HOLZÄPFEL F. Probabilistic two-phase wake vortex decay and transport model[J]. Journal of Aircraft, 2003, 40(2): 323–331. doi:  10.2514/2.3096
    [8] HOLZÄPFEL F. Sensitivity analysis of the effects of aircraft and environmental parameters on aircraft wake vortex trajectories and lifetimes[C]. 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Dallas, USA, 2013: 7–10.
    [9] DE VISSCHER I, WINCKELMANS G, LONFILS T, et al. The WAKE4D simulation platform for predicting aircraft wake vortex transport and decay: Description and examples of application[C]. AIAA Atmospheric and Space Environments Conference, Toronto, Canada, 2010: 7994.
    [10] SARPKAYA T, ROBINS R E, and DELISI D P. Wake-vortex eddy-dissipation model predictions compared with observations[J]. Journal of Aircraft, 2001, 38(4): 687–692. doi:  10.2514/2.2820
    [11] DELISI D P. Development of the VIPER fast-time wake vortex model (development, assumptions, examples, and plans)[C]. WakeNet3-Europe Specific Workshop on "Operational Wake Vortex Models", Belgium, 2011.
    [12] PROCTOR F, HAMILTON D, and SWITZER G. TASS driven algorithms for wake prediction[C]. 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, USA, 2006.
    [13] HOLZÄPFEL F. On the maturity of wake vortex observation, prediction, and validation[C]. WakeNet3-Eurooe 1st Workshop on "Wake Turbulence Safety in Future Aircraft Operations", Paris, France, 2009.
    [14] DE VISSCHER I, WINCKELMANS G, and BRICTEUX L. Some reflections on the achievable quality of operational WV prediction using operational meteorological and aircraft inputs[C]. WakeNet3-Europe 1st Workshop on "Wake Turbulence Safety in Future Aircraft Operations", Paris, France, 2009.
    [15] PRUIS M J. Development of a new probabilistic wake vortex prediction model[C]. WakeNet3-Europe Specific Workshop on "Operational Wake Vortex Models", Belgium, 2011.
    [16] SCHÖNHALS S, STEEN M, and HECKER P. Wake vortex prediction and detection utilising advanced fusion filter technologies[J]. The Aeronautical Journal, 2011, 115(1166): 221–228. doi:  10.1017/S0001924000005674
    [17] LI Jianbing, CHAN P W, WANG Tao, et al. Circulation retrieval of wake vortex with a side-looking scanning Lidar[C]. CIE International Conference on Radar, Guangzhou, China, 2016: 1–4.
    [18] JULIER S, UHLMANN J, and DURRANT-WHYTE H F. A new method for the nonlinear transformation of means and covariances in filters and estimators[J]. IEEE Transactions on Automatic Control, 2000, 45(3): 477–482. doi:  10.1109/9.847726
    [19] JULIER S J and UHLMANN J K. Unscented filtering and nonlinear estimation[J]. Proceedings of the IEEE, 2004, 92(3): 401–422. doi:  10.1109/JPROC.2003.823141
    [20] GERZ T, HOLZÄPFEL F, and DARRACQ D. Commercial aircraft wake vortices[J]. Progress in Aerospace Sciences, 2002, 38(3): 181–208. doi:  10.1016/S0376-0421(02)00004-0
    [21] 屈龙海. 晴空和湿性大气中飞机尾流雷达散射特性的研究[D]. [博士论文], 国防科技大学, 2015: 29–31.

    QU Longhai. Study on the radar scattering characteristics of aircraft wake vortex in clear air and moist air[D]. [Ph. D. dissertation], National University of Defense Technology, 2015: 29–31.
    [22] AHMAD N N and PROCTOR F. Review of idealized aircraft wake vortex models[C]. 52nd Aerospace Sciences Meeting, National Harbor, USA, 2014.
    [23] BURNHAM D. Chicago monostatic acoustic vortex sensing system, Volume I: Data collection and reduction[R]. FAA-RD-79-103, 1979.
    [24] SHEN Chun, LI Jianbing, ZHANG Fulin, et al. Two-step locating method for aircraft wake vortices based on Gabor filter and velocity range distribution[J]. IET Radar, Sonar & Navigation, 2020, 14(12): 1958–1967. doi:  10.1049/iet-rsn.2020.0319
    [25] LI Jianbing, SHEN Chun, GAO Hang, et al. Path Integration (PI) method for the parameter-retrieval of aircraft wake vortex by Lidar[J]. Optics Express, 2020, 28(3): 4286–4306. doi:  10.1364/OE.382968
    [26] 焦云涛. 低空风切变与飞行安全[J]. 民航经济与技术, 1994, (11): 13–14.

    JIAO Yuntao. Windshear at low altitude and flight safety[J]. Civil Aviation Economics and Technology, 1994, (11): 13–14.
    [27] WILSON D K, OSTASHEV V E, GOEDECKE G H, et al. Quasi-wavelet calculations of sound scattering behind barriers[J]. Applied Acoustics, 2004, 65(6): 605–627. doi:  10.1016/j.apacoust.2003.11.009
    [28] 张宏昇. 大气湍流基础[M]. 北京: 北京大学出版社, 2014: 161–165.

    ZHANG Hongsheng. Atmospheric Turbulence Foundation[M]. Beijing: Peking University Press, 2014: 161–165.
    [29] 李金梁. 箔条干扰的特性与雷达抗箔条技术研究[D]. [博士论文], 国防科学技术大学, 2010: 57–58.

    LI Jinliang. Study on characteristics of chaff jamming and anti-chaff technology for radar[D]. [Ph.D. dissertation], National University of Defense Technology, 2010: 57–58.
    [30] SMALIKHO I N, BANAKH V A, HOLZÄPFEL F, et al. Method of radial velocities for the estimation of aircraft wake vortex parameters from data measured by coherent Doppler Lidar[J]. Optics Express, 2015, 23(19): A1194–A1207. doi:  10.1364/OE.23.0A1194
    [31] 沈淳, 高航, 王雪松, 等. 基于激光雷达探测的飞机尾流特征参数反演系统[J]. 雷达学报, 2020, 9(6): 1032–1044. doi:  10.12000/JR20046

    SHEN Chun, GAO Hang, WANG Xuesong, et al. Aircraft wake vortex parameter-retrieval system based on Lidar[J]. Journal of Radars, 2020, 9(6): 1032–1044. doi:  10.12000/JR20046
    [32] GAO Hang, LI Jianbing, CHAN P W, et al. Parameter-retrieval of dry-air wake vortices with a scanning Doppler Lidar[J]. Optics Express, 2018, 26(13): 16377–16392. doi:  10.1364/OE.26.016377
  • [1] 薛策文, 冯晅, 李晓天, 梁文婧, 周皓秋, 王颖.  全极化探地雷达多极化数据融合分析研究 . 雷达学报, doi: 10.12000/JR20104
    [2] 刘宁波, 丁昊, 黄勇, 董云龙, 王国庆, 董凯.  X波段雷达对海探测试验与数据获取年度进展 . 雷达学报, doi: 10.12000/JR21011
    [3] 李健兵, 王雪松.  分布式软目标雷达特性与感知技术述评 . 雷达学报, doi: 10.12000/JR20052
    [4] 施龙飞, 全源, 范金涛, 马佳智.  通信化雷达探测技术 . 雷达学报, doi: 10.12000/JR20088
    [5] 罗迎, 倪嘉成, 张群.  基于“数据驱动+智能学习”的合成孔径雷达学习成像 . 雷达学报, doi: 10.12000/JR19103
    [6] 沈淳, 高航, 王雪松, 李健兵.  基于激光雷达探测的飞机尾流特征参数反演系统 . 雷达学报, doi: 10.12000/JR20046
    [7] 刘宁波, 董云龙, 王国庆, 丁昊, 黄勇, 关键, 陈小龙, 何友.  X波段雷达对海探测试验与数据获取 . 雷达学报, doi: 10.12000/JR19089
    [8] 李道京, 胡 烜.  合成孔径激光雷达光学系统和作用距离分析 . 雷达学报, doi: 10.12000/JR18017
    [9] 张珂殊, 潘洁, 王然, 李光祚, 王宁, 吴一戎.  大幅宽激光合成孔径雷达成像技术研究 . 雷达学报, doi: 10.12000/JR16152
    [10] Hon Kaikwong, Chan Pakwai.  Aircraft Wake Vortex Observations in Hong Kong . 雷达学报, doi: 10.12000/JR17072
    [11] 刘俊凯, 李健兵, 马梁, 陈忠宽, 蔡益朝.  基于矩阵信息几何的飞机尾流目标检测方法 . 雷达学报, doi: 10.12000/JR17058
    [12] 胡程, 刘长江, 曾涛.  双基地前向散射雷达探测与成像 . 雷达学报, doi: 10.12000/JR16058
    [13] 高云泽, 董泽华, 方广有, 纪奕才, 周斌.  嫦娥三号测月雷达第一通道数据处理与分析 . 雷达学报, doi: 10.12000/JR15030
    [14] 闫召爱, 胡雄, 郭商勇, 程永强, 郭文杰, 潘艺升.  星载钠荧光多普勒激光雷达性能分析 . 雷达学报, doi: 10.12000/JR14140
    [15] 王福友, 罗钉, 刘宏伟.  低分辨机载雷达飞机目标分类识别技术研究 . 雷达学报, doi: 10.3724/SP.J.1300.2014.14075
    [16] 马萌, 李道京, 杜剑波.  振动条件下机载合成孔径激光雷达成像处理 . 雷达学报, doi: 10.3724/SP.J.1300.2014.13132
    [17] 杜剑波, 李道京, 马萌.  机载合成孔径激光雷达相位调制信号性能分析和成像处理 . 雷达学报, doi: 10.3724/SP.J.1300.2014.13094
    [18] 李道京, 张清娟, 刘波, 杨宏, 潘洁.  机载合成孔径激光雷达关键技术和实现方案分析 . 雷达学报, doi: 10.3724/SP.J.1300.2013.13021
    [19] 李伟, 王兴亮, 邹鲲, 徐艺盟, 张群.  基于数据融合和陷波滤波的MIMO 雷达抗欺骗干扰算法 (英文) . 雷达学报, doi: 10.3724/SP.J.1300.2012.20060
    [20] 吴谨.  关于合成孔径激光雷达成像研究 . 雷达学报, doi: 10.3724/SP.J.1300.2012.20076
  • 加载中
图(16) / 表 (4)
计量
  • 文章访问数:  73
  • HTML全文浏览量:  69
  • PDF下载量:  23
  • 被引次数: 0
出版历程
  • 收稿日期:  2021-01-22
  • 修回日期:  2021-03-16
  • 网络出版日期:  2021-04-02

基于数据同化的飞机尾流行为预测

doi: 10.12000/JR21007
    基金项目:  国家自然科学基金(61490649, 61771479, 61625108),湖南省杰出青年基金(2018JJ1030)
    作者简介:

    沈 淳(1985–),男,福建漳州人,博士生,工程师,研究方向为空间信息获取与处理

    李健兵(1979–),男,湖南邵东人,博士,国防科技大学电子科学学院教授,主要研究方向为新体制雷达、雷达信号处理

    王雪松(1972–),男,内蒙古人,博士,国防科技大学电子科学学院教授,主要研究方向为极化信息处理、新体制雷达技术、电子对抗

    通讯作者: 李健兵 jianbingli@nudt.edu.cn
  • 责任主编:夏海云 Corresponding Editor: XIA Haiyun
  • 中图分类号: TN955+.1

摘要: 飞机尾流是飞机飞行时在其后方产生的一对反向旋转的强烈湍流,对后续飞机飞行以及机场安全起降影响极大,其演化趋势的预测已成为空中交通安全管制的瓶颈,亟需发展基于实时探测数据的飞机尾流行为预测技术。在雷达探测反演得到的尾流涡心位置和速度环量等特征参数基础上,开展飞机尾流行为预测分析,能够预知飞机尾流危害区域,为机场安全起降动态间隔标准制定提供技术支撑。该文结合风场线性切变和最小二乘拟合方法构建了参数化尾流行为预测模型,解决了经典尾流预测模型气象环境参数未随时间演化实时调整的问题。该文根据复杂风场非线性演化特点,设计了基于无迹卡尔曼滤波的数据同化模型,利用雷达探测数据对尾流行为预测进行实时修正。数值仿真验证和实测数据验证结果表明,基于数据同化的飞机尾流行为预测方法能够根据实时探测数据对尾流行为预测轨迹进行修正,得到更加贴近实测的飞机尾流行为预测轨迹。

注释:
1)  责任主编:夏海云 Corresponding Editor: XIA Haiyun

English Abstract

沈淳, 李健兵, 高航, 等. 基于数据同化的飞机尾流行为预测[J]. 雷达学报, 待出版. doi:  10.12000/JR21007
引用本文: 沈淳, 李健兵, 高航, 等. 基于数据同化的飞机尾流行为预测[J]. 雷达学报, 待出版. doi:  10.12000/JR21007
SHEN Chun, LI Jianbing, GAO Hang, et al. Aircraft wake vortex behavior prediction based on data assimilation[J]. Journal of Radars, in press. doi:  10.12000/JR21007
Citation: SHEN Chun, LI Jianbing, GAO Hang, et al. Aircraft wake vortex behavior prediction based on data assimilation[J]. Journal of Radars, in press. doi:  10.12000/JR21007
    • 飞机尾流是机翼上下表面压力差而在其后方形成的一种反向旋转的强烈气流(如图1所示),是飞机飞行时产生升力的必然产物,具有空间尺度大、持续时间长、旋转强烈等特点。在航空安全管制方面,飞机尾流会对后续进入起飞/降落滑道或预定航线的飞机造成影响,可能发生颠簸、横滚,乃至失去控制。在大气物理方面,飞机尾流的存在会改变大气中水的相态,产生凝结尾迹、管道云等自然现象,甚至可能影响地球热辐射和小范围的天气[1]。飞机尾流的特性、探测和演化规律研究已成为当前空中交通安全管制、大气物理等领域共同关注的前沿科学问题。

      图  1  飞机尾流形成示意图

      Figure 1.  Illustration of aircraft wake vortex

      在空中交通管制领域,为了避免飞机因遭遇前机产生的飞机尾流而发生事故,国际民航组织(International Civil Aviation Organization, ICAO)于20世纪70年代针对不同飞机的起飞重量规定了相邻两架飞机起飞/降落的最小安全距离,但这些规则比较保守,很大程度上限制了机场的起降容量[2-4],难以满足航空业快速发展的需求。欧美等国家在此基础上开展了飞机重分类(Recategorization, RECAT)[5,6]研究并形成适合当地地形和气象环境条件的分类体系。中国民航引入了欧美飞机重分类并开始在广州和深圳机场试运行,但欧美RECAT规则主要针对欧美机场运行需求开发,具有地形和气象环境条件约束,我们对该规则背后物理规律的认识还不充分。我国正处于自主创新、快速发展的历史机遇期,开展飞机尾流短时行为预测,为后续飞机尾流危害评估和构建适合中国民航的尾流动态间隔标准奠定基础,对我国掌握相关核心技术和突破国际垄断,具有重大的应用价值和科学意义。

      飞机尾流行为预测是根据飞机参数(重量、翼展、飞行速度、飞行高度等)以及气象环境参数(横向风、风切变、大气湍流、大气分层等)对其飞行时产生的左右涡旋的演化和消散进行反演,预知飞机尾流的运动路线和强度变化情况,提前确定飞机尾流危害区域,达到有效缩短飞机起降安全间隔,提升航空效率和机场容量的目的。

      基于上述目的,欧美等国家开展了长期研究并建立了部分尾流行为预测模型,欧洲的德国宇航中心(Deutsches Zentrum furLuft-und Raumfahrt, DLR)的Deterministic/Probabilistic Two-Phase wake vortex decay and transport model (D2P/P2P)[7,8]模型,将尾流演化分成初始和快速消散两个阶段,考虑了飞机重量、翼展、飞行速度、位置、飞行航迹倾角等飞行参数和迎风速度、横向风、风切变、湍流消散率等天气参数。比利时鲁汶天主教大学(Universite Catholique de Louvain, UCL)[9]建立了DVM/PVM (Deterministic/Probabilistic wake Vortex Model)模型,该模型基于离散尾涡粒子系统并结合Monte-Carlo仿真生成尾流涡旋和近地镜面涡旋,考虑了飞机位置、翼展、重量、飞行速度等参数。美国的Sparpkaya, Robins, Delisi等研究团队[10]联合构建了APA (AVOSS Prediction Algorithm)模型,该模型考虑了横向风、气温、湍流强度等天气参数并利用镜像涡理论模拟近地面效应。美国西北研究协会(Northwest Research Associates, NWRA)[11]建立了VIPER (Vortex algorithm Including Parameterized Entrainment Results)模型,利用从大涡模拟(Large Eddy Simulation, LES)中参数化学习建立的方程分别对左旋和右旋两个尾涡进行预测。TDAWP (TASS Driven Algorithm for Wake Prediction)[12]模型由美国国家航空航天局(National Aeronautics and Space Administration, NASA)建立,该模型基于控制容积法和流体动量、角动量守恒原理,可在不同天气条件和飞机飞行参数下进行预测。

      上述尾流行为预测模型虽然考虑了多种气象条件因素和飞机参数的影响,但是依赖初始风场气象条件的设置,在预测尾流演化过程中风场演化非常复杂。随着时间的变化,初始气象条件和飞机参数的设定与实际气象条件等存在误差,导致涡心位置和速度环量等趋势预测与尾流实际演化趋势存在偏差,进而影响了飞机尾流行为预测的准确性和鲁棒性[13,14]。在此基础上,Pruis等人[15]提出综合运用多个预测模型,给出尾流演化的概率性预测结果,但这种综合方法仍然依赖于初始气象条件的设置,未能解决气象条件实时变化的问题。Schöenhals等人[16]提出应用数据融合的方法,利用线性卡尔曼滤波或扩展卡尔曼滤波对尾流演化过程中预测和实测数据之间的误差进行估计并对预测数据进行修正,但其仍以经典D2P模型为尾流预测模型,并未根据实际风场演化来实时调整预测模型气象参数,修正鲁棒性和后续预测效果不佳。

      本文设计了非线性卡尔曼滤波数据同化预测模型,对尾流预测模型气象环境参数进行实时拟合修正,能根据实时探测数据对代表飞机尾流演化趋势的涡心位置和速度环量等进行实时修正并重新预测。

      针对经典尾流演化预测模型比较依赖初始气象条件和飞机飞行参数等初始条件的问题,结合复杂风场非线性演化的特点,本文提出了基于无迹卡尔曼滤波(Unscented Kalman Filter, UKF)数据同化的飞机尾流行为预测方法,实现对飞机尾流行为演化的实时修正。针对经典尾流演化预测模型初始气象条件和飞机飞行参数设置后未根据实际风场演化实时调整的问题,应用风场线性切变和最小二乘曲线拟合方法,根据实时探测获得的尾流风场数据计算得到预测模型的气象环境参数,构建了参数化尾流实时演化预测模型。

      本文的组织结构如下:第2节对基于数据同化的飞机尾流行为预测思路进行描述;第3节根据风场线性切变和最小二乘拟合理论建立参数化飞机尾流演化预测模型;第4节设计数据同化预测模型,用于尾流行为预测的修正;第5节通过仿真和实测数据验证基于数据同化的飞机尾流行为预测方法的有效性;第6节给出本文的结论。

    • 在晴空条件下,飞机尾流回波散射主要由尾流内部的浮尘等微粒决定,激光雷达波束的衰减较小,探测距离较远且时间和角度分辨率较高,所以激光雷达是当前探测晴空飞机尾流广泛使用的传感器[17]。激光雷达常用的扫描方式为距离高度指示器(Range-Height Indicator, RHI)扫描,通过上下扫描飞机尾流可以得到扫描平面上各个探测单元在不同扫描时刻的飞机尾流回波信号(如图2所示)。激光雷达探测飞机尾流一般采用侧视方式,雷达设置于机场跑道一侧,在与跑道正交的扫描平面上上下来回扫描飞机尾流,激光波束的仰角范围为$[{\alpha _{\min }} - {\alpha _{\max }}]$。当飞机在跑道上起飞或降落,将会在机场跑道上方的起飞/降落滑道上产生一对旋转方向相反的漩涡,它们的涡心位置分别为${O_{c1}}$${O_{c2}}$

      图  2  激光雷达探测飞机尾流场景设置

      Figure 2.  Geometry configuration for Lidar detection of wake vortex

      激光雷达扫描飞机尾流时,尾流中悬浮的粒子会受到激光的照射并散射能量。激光雷达接收粒子散射回来的光信号后,将其转化为电信号,再通过相干处理,得到尾流区域的多普勒信息(如图3所示)。通过处理多普勒信息,可以反演飞机尾流的涡心位置和速度环量等特征参数,为飞机尾流行为预测提供数据支撑。

      图  3  飞机尾流多普勒速度RHI回波分布图

      Figure 3.  Doppler velocity distribution of wake vortex in an RHI

    • 基于数据同化的飞机尾流行为预测(Data Assimilation, DA)方法如图4所示。应用激光雷达探测数据反演得到的尾流涡心位置和速度环量等特征参数,基于经典尾流预测模型和非线性卡尔曼滤波,构建参数化尾流预测模型进行尾流演化行为预测,同时利用演化过程中实时探测到的尾流特征数据进行非线性卡尔曼滤波,对预测的涡心位置和速度环量等特征参数进行修正并重新进行行为预测,同时利用实测数据对预测模型的气象环境参数进行曲线拟合修正,使尾流行为预测更加贴近真实。具体预测步骤如下:

      图  4  飞机尾流行为预测流程图

      Figure 4.  Wake vortex behavior prediction flow chart

      (1) 根据激光雷达实测获得的尾流初始涡心位置和速度环量等特征参数,以及雷达和风速仪等设备测量反演得到的风速、风切变系数、湍流耗散率等风场数据,应用参数化尾流演化预测模型设计非线性状态传输函数,预测得到飞机尾流在一定时间段内的涡心演化和强度变化趋势。

      (2) 在飞机尾流演化过程中,根据接收到的实时测量数据,将速度环量和涡心水平位置、高度等作为预估量,应用无迹卡尔曼滤波,将非线性状态传输得到的尾流特征参数预测值与实测数据进行数据同化,计算对应的状态估计值、估计方差以及卡尔曼滤波增益系数,对代表尾流演化趋势的速度环量、左右涡心位置等特征参数进行修正。同时,根据实测数据进行最小二乘曲线拟合更新修正速度环量演化模型的气象环境参数,根据实测风场数据进行涡心演化模型的背景切向风场速度和尾涡下降速度等参数的实时估计,利用修正后的气象环境参数和速度参数对数据同化的非状态传递函数进行修正。

      (3) 根据修正后的数据同化模型和尾流涡心位置、速度环量等特征参数对尾流演化趋势进行更新预测得到更加贴近真实的预测数据,演化过程中,一旦得到新的尾流实测数据,就重复第(2)步和第(3)步。

    • 参数化飞机尾流演化预测模型在经典尾流演化模型基础上,根据雷达连续探测反演得到的数个涡心位置和速度环量等特征参数进行最小二乘拟合得到速度环量预测模型参数,根据实时探测到的风场数据并利用线性切变得到涡心演化模型参数,进而对飞机尾流的强度(速度环量)和涡心位置的演化趋势进行估计,并在演化过程中利用数据同化对涡心位置、速度环量以及预测模型的气象环境参数进行不断修正,实现速度环量和涡心位置的实时行为预测。

    • 经典D2P模型[7]将尾流速度环量演化分为稳定段和快速消散段两个阶段,模型如式(1)

      $${\varGamma ^{\rm{*}}}({t^*}) = \left\{ \begin{aligned} &{A - \exp \frac{{ - B}}{{\upsilon _1^*({t^*} - T_1^*)}},}&{{t^*} \le T_2^*} \\ &{A - \exp \frac{{ - B}}{{\upsilon _1^*({t^*} - T_1^*)}} - \exp \frac{{ - B}}{{\upsilon _2^*({t^*} - T_2^*)}},}&{{t^*} > T_2^*} \end{aligned} \right.$$ (1)

      其中,$A$, $\upsilon _1^*$, $T_1^*$, $B$等参数通过大量流体动力学仿真总结得到,在预测起始时设置为固定的常数,典型值为$A = 1.1418$, $\upsilon _1^* = 1.78\times{10^{ - 3}}$, $T_1^* = - 3.48$, $B = 0.0121$$\upsilon _2^*$$T_2^*$由湍流耗散率(Eddy Dissipation Rate, EDR)、浮力频率(Brunt-Vaisala frequency)、大气分层系数(Atmospheric stratification)等气象条件参数决定,具体可参考文献[7]。

      尾流在实际演化过程中,受低空风场各种因素的叠加影响,强度变化趋势复杂,但经典D2P模型的参数设置并未考虑风场演化因素的实时影响,无法真实反映出尾流强度的变化趋势,需根据实时探测反演特征值对速度环量演化模型进行修正。本参数化预测模型将$A$, $\upsilon _1^*$, $T_1^*$, $B$等模型参数作为未知量,通过雷达连续探测反演得到的$N(N \ge 4)$个初始环量值,结合式(1)进行最小二乘曲线拟合求得,同时在演化过程中结合后续探测反演得到的环量值,不断对$A$, $\upsilon _1^*$, $T_1^*$, $B$等模型参数进行最小二乘拟合修正。$\upsilon _2^*$$T_2^*$则由湍流耗散率、浮力频率、大气分层系数等实时探测得到的气象条件参数决定。

    • 飞机尾流受左右涡相互作用${V_{\rm{d}}}$不断下降,同时受背景切向风场${{{\mathit{\boldsymbol{V}}}}_{\rm{c}}}$的影响进行漂移。左右涡心相互作用的下降速度$V_{\rm{d}}^{\rm{L}}$$V_{\rm{d}}^{\rm{R}}$分别为

      $$V_{\rm{d}}^{\rm{L}} = \frac{{{\varGamma _{\rm{R}}}}}{{2{\rm{\pi }}{b_0}}},\;\;V_{\rm{d}}^{\rm{R}} = \frac{{{\varGamma _{\rm{L}}}}}{{2{\rm{\pi }}{b_0}}}$$ (2)

      其中,${b_0}$为左右涡心间距,${\varGamma _{\rm{L}}}$$ {\varGamma }_{\rm{R}}$分别为左涡和右涡速度环量值。

      背景切向风场${{{\mathit{\boldsymbol{V}}}}_{\rm{c}}}$应用线性切变理论进行求解,在一个RHI扫描周期中,背景风场由$x$$y$方向的两部分分量组成:

      $${{{\mathit{\boldsymbol{V}}}}_{\rm{c}}} = V_{\rm{c}}^{\rm{x}} \cdot {{\mathit{\boldsymbol{x}}}} + V_{\rm{c}}^{\rm{y}} \cdot {{\mathit{\boldsymbol{y}}}} $$ (3)

      其中,$x$方向的分量可假设为线性切变:

      $$V_{\rm{c}}^{\rm{x}} = V_{\rm{c}}^0 + \beta \cdot y$$ (4)

      $V_{\rm{c}}^0$$x$方向的分量在$y = 0$处的速度,$\beta $为风切变变化率,$y$为垂直方向的高度。背景切向风场在一个RHI扫描周期可假设为基本不变,因此背景切向风场${{{\mathit{\boldsymbol{V}}}}_{\rm{c}}}$可由不存在尾流的区域估计得到,区域划分如图5所示。当使用$M(M \ge 3)$个风场多普勒速度点时,参数$V_{\rm{c}}^0$, $\beta $$V_{\rm{c}}^{\rm{y}}$可以计算获得。

      图  5  估计背景风场的非尾流区域

      Figure 5.  Regions free of wake vortex that was used to estimate the background wind

      根据求解得到的背景切向风场${{{\mathit{\boldsymbol{V}}}}_{\rm{c}}}$和涡心相互作用下降速度${V_{\rm{d}}}$,垂直方向的尾流涡心轨迹演化模型如式(5)

      $${y^*}({t^*}) = {y_0} + ( - {V_{\rm{d}}} + V_{\rm{c}}^{\rm{y}}){t^*}$$ (5)

      水平方向的尾流涡心轨迹演化模型如式(6)

      $${x^*}({t^*}) = {x_0} + (\beta y + V_{\rm{c}}^0){t^*}$$ (6)

      在尾流演化过程中,使用上述演化模型来反演涡心轨迹垂直方向和水平方向的位置。在数据同化过程中,当探测到新的雷达数据时,实时根据新的风场数据进行$V_{\rm{c}}^0$, $\beta $$V_{\rm{c}}^{\rm{y}}$等尾流背景风场参数的估计并重新进行预测。

    • 数据同化模型基于无迹卡尔曼滤波,利用雷达探测得到的飞机尾流涡心位置和速度环量等特征参数对尾流行为预测进行修正,达到根据实测数据进行实时修正的目的。Julier等人[18]根据近似非线性函数的均值和方差远比近似非线性函数本身更容易的特点,提出了基于确定性采样的UKF算法。

      本文基于无迹卡尔曼滤波,针对尾流演化特点和雷达探测实际,设计了非线性传递函数和量测函数,采用无损(Unscented Transformation, UT)变换[19], 利用比例修正对称采样策略获取的Sigma采样点通过非线性函数的传递进行尾流行为预测,进一步结合实测数据并利用量测函数实时修正行为预测结果。

    • 相比于线性卡尔曼滤波,无迹卡尔曼滤波能够模拟非线性系统变化且估计精度能够达到泰勒级数展开的2阶精度,可应用于复杂风场演化预测。飞机尾流演化过程中受复杂风场影响,其演化规律呈现非线性特点,其非线性演化传递和量测函数构建如式(7)

      $$\left\{ \begin{aligned} &{{{{\mathit{\boldsymbol{x}}}}_{k + 1}} = f({{{\mathit{\boldsymbol{x}}}}_k})} \\ &{{{{\mathit{\boldsymbol{z}}}}_{k + 1}} = h({{{\mathit{\boldsymbol{x}}}}_{k + 1}})} \end{aligned} \right.$$ (7)

      其中,${{{\mathit{\boldsymbol{x}}}}_k} = [{\varGamma _{k,1}},{\varGamma _{k,2}},{x_{k,1}},{x_{k,2}},{y_{k,1}},{y_{k,2}}]$为传递状态矢量,${\varGamma _{k,1}}$${\varGamma _{k,2}}$分别为左涡和右涡速度环量值,${x_{k,1}}$${x_{k,2}}$分别为左右涡心水平位置,${y_{k,1}}$${y_{k,2}}$分别为左右涡心垂直位置,${{{\mathit{\boldsymbol{z}}}}_k} = {{{\mathit{\boldsymbol{v}}}}_{\rm{D}}}$为量测矢量,${{{\mathit{\boldsymbol{v}}}}_{\rm{D}}}$为多普勒速度矢量,$f( \cdot )$$h( \cdot )$分别为系统的非线性状态函数和量测函数。

      无迹卡尔曼数据同化的飞机尾流行为预测方法以比例修正对称采样策略选取Sigma点,利用Sigma点通过UT变换来近似尾流轨迹和速度环量非线性变化函数以及量测函数的均值和方差,其状态演化传递流程如下。

      假设$k$时刻的状态估值${\widehat{{\mathit{\boldsymbol{x}}}}_{k|k}}$和估计方差${{{\mathit{\boldsymbol{P}}}}_{k|k}}$,由比例修正对称采样策略,可得到$2n + 1$个Sigma采样矢量点$\chi _i^{'}{\rm{ }}$,以及对应的权值${(W_i^{\rm{m}})^{'}}$${(W_i^{\rm{c}})^{'}}$,然后将采样点进行非线性状态函数传递得:$\gamma _{k + 1|k}^i = f(\chi _i^{'})$,对应演化传递的均值和方差如式(8)

      $$\left\{ \begin{aligned} &{{{\widehat{{\mathit{\boldsymbol{x}}}}}_{k + 1|k}} = \sum\limits_{i = 0}^{2n} {{{(W_i^{\rm{m}})}^{'}}\gamma _{k + 1|k}^i} } \\ & \begin{split} &{{{\mathit{\boldsymbol{P}}}}_{k + 1|k}} = \sum\limits_{i = 0}^{2n} {{(W_i^{\rm{c}})}^{'}}\left(\gamma _{k + 1|k}^i - {{\widehat{{\mathit{\boldsymbol{x}}}}}_{k + 1|k}}\right)\\ &\quad\quad\quad\quad\cdot{{\left(\gamma _{k + 1|k}^i - {{\widehat{{\mathit{\boldsymbol{x}}}}}_{k + 1|k}}\right)}^{\rm{T}}} \end{split}\end{aligned} \right.{\rm{ }}$$ (8)

      $k + 1$时刻没有雷达测量得到的量测值时,传递状态矢量${\widehat{{\mathit{\boldsymbol{x}}}}_{k + 1|k}}$作为数据同化结果输出为尾流演化预测结果;当$k + 1$时刻存在雷达测量得到的量测值时,根据状态演化得到的${\widehat{{\mathit{\boldsymbol{x}}}}_{k + 1|k}}$${{{\mathit{\boldsymbol{P}}}}_{k + 1|k}}$及采样策略公式可得$2n + 1$个Sigma采样矢量点$\zeta _i^{'}$和相应的权值${(W_i^{\rm{m}})^{'}}$${(W_i^{\rm{c}})^{'}}$,经过非线性量测函数传递得$\zeta _{k + 1|k}^{'} = h(\zeta _i^{'})$,量测变量的均值、方差以及协方差为

      $$\left\{ \begin{aligned} & {{{\widehat{{\mathit{\boldsymbol{z}}}}}_{k + 1|k}} = \sum\limits_{i = 0}^{2n} {{{(W_i^{\rm{m}})}^{'}}\zeta _{k + 1|k}^i} } \\ &{{{\mathit{\boldsymbol{P}}}}_{{\rm{zz}},k + 1|k}} = \sum\limits_{i = 0}^{2n} {{(W_i^{\rm{c}})}^{'}}(\zeta _{k + 1|k}^i - {{\widehat{{\mathit{\boldsymbol{z}}}}}_{k + 1|k}})\\ &\quad\quad\quad\quad\quad\cdot{{(\zeta _{k + 1|k}^i - {{\widehat{{\mathit{\boldsymbol{z}}}}}_{k + 1|k}})}^{\rm{T}}} \\ & {{{\mathit{\boldsymbol{P}}}}_{{\rm{xz}},k + 1|k}} = \sum\limits_{i = 0}^{2n} {{(W_i^{\rm{c}})}^{'}}(\gamma _{k + 1|k}^i - {{\widehat{{\mathit{\boldsymbol{x}}}}}_{k + 1|k}})\\ &\quad\quad\quad\quad\quad\cdot{{(\zeta _{k + 1|k}^i - {{\widehat{{\mathit{\boldsymbol{z}}}}}_{k + 1|k}})}^{\rm{T}}} \end{aligned} \right.$$ (9)

      根据该时刻的实测向量${{{\mathit{\boldsymbol{z}}}}_{k + 1}}$可求出滤波增益${{{\mathit{\boldsymbol{K}}}}_{k + 1}}$以及状态估计和估计方差

      $$\left\{ \begin{aligned} &{{{\widehat{{\mathit{\boldsymbol{x}}}}}_{k + 1|k{\rm{ + }}1}}{\rm{ = }}{{\widehat{{\mathit{\boldsymbol{x}}}}}_{k + 1|k}}{\rm{ + }}{{{\mathit{\boldsymbol{K}}}}_{k + 1}}({{{\mathit{\boldsymbol{z}}}}_{k + 1}} - {{\widehat{{\mathit{\boldsymbol{z}}}}}_{k + 1|k}})} \\ &{{{{\mathit{\boldsymbol{K}}}}_{k + 1}} = {{{\mathit{\boldsymbol{P}}}}_{{\rm{xz}},k + 1|k}}{{({{{\mathit{\boldsymbol{P}}}}_{{\rm{zz}},k + 1|k}})}^{ - 1}}} \\ & {{{{\mathit{\boldsymbol{P}}}}_{k + 1|k + 1}} = {{{\mathit{\boldsymbol{P}}}}_{k + 1|k}} - {{{\mathit{\boldsymbol{K}}}}_{k + 1}}{{{\mathit{\boldsymbol{P}}}}_{{\rm{zz}},k + 1|k}}{{\mathit{\boldsymbol{K}}}}_{k + 1}^{\rm{T}}{\rm{ }}} \end{aligned} \right.{\rm{ }}$$ (10)

      ${\widehat{{\mathit{\boldsymbol{x}}}}_{k + 1|k{\rm{ + }}1}}$作为数据同化结果输出为尾流演化预测结果,同时${\widehat{{\mathit{\boldsymbol{x}}}}_{k + 1|k{\rm{ + }}1}}$${{{\mathit{\boldsymbol{P}}}}_{k + 1|k + 1}}$作为下一时刻状态传输估计的输入。

    • 飞机尾流行为预测主要是对代表飞机尾流强度和位置的速度环量与涡心位置等特征参数演化进行估计,系统非线性状态传递函数$f( \cdot )$为状态向量的步进预测

      $${{{\mathit{\boldsymbol{x}}}}_{k + 1}} = f({{{\mathit{\boldsymbol{x}}}}_k}) = \left[ \begin{aligned} & {{\varGamma _{k,1}} + \frac{{\partial {\varGamma _{k,1}}}}{{\partial {t^*}}}\Delta {t^*}} \\ &{{\varGamma _{k,2}} + \frac{{\partial {\varGamma _{k,2}}}}{{\partial {t^*}}}\Delta {t^*}} \\ &{{x_{k,1}} + \frac{{\partial {x_{k,1}}}}{{\partial {t^*}}}\Delta {t^*}} \\ &{{x_{k,2}} + \frac{{\partial {x_{k,2}}}}{{\partial {t^*}}}\Delta {t^*}} \\ &{{y_{k,1}} + \frac{{\partial {y_{k,1}}}}{{\partial {t^*}}}\Delta {t^*}} \\ & {{y_{k,2}} + \frac{{\partial {y_{k,2}}}}{{\partial {t^*}}}\Delta {t^*}} \end{aligned} \right]{\rm{ }}$$ (11)

      根据尾流预测模型,环量${\varGamma _k}$,涡心水平位置${x_k}$和涡心垂直位置${y_k}$关于时间${t^*}$的导数为

      $$\frac{{\partial {\varGamma _k}}}{{\partial {t^*}}} = \left\{ \begin{aligned} & { - \frac{{{B^2}}}{{\upsilon _1^*{{({t^*} - T_1^*)}^2}}}\exp \frac{{ - {B^2}}}{{\upsilon _1^*({t^*} - T_1^*)}},}&{{t^*} \le T_2^*} \\ &{ - \frac{{{B^2}}}{{\upsilon _1^*{{({t^*} - T_1^*)}^2}}}\exp \frac{{ - {B^2}}}{{\upsilon _1^*({t^*} - T_1^*)}} - \frac{{{B^2}}}{{\upsilon _2^*{{({t^*} - T_2^*)}^2}}}\exp \frac{{ - {B^2}}}{{\upsilon _2^*({t^*} - T_2^*)}},}&{{t^*} > T_2^*} \end{aligned} \right.$$ (12)
      $$\frac{{\partial {x_k}}}{{\partial {t^*}}} = \beta {y_k} + V_{\rm{c}}^0$$ (13)
      $$\frac{{\partial {y_k}}}{{\partial {t^*}}} = - {V_{\rm{d}}} + V_{\rm{c}}^{\rm{y}}$$ (14)
    • 排除背景风的情况下,雷达探测可得到扫描平面的风场多普勒速度分布,该多普勒速度为飞机尾流涡旋切向速度在雷达扫描波束上的投影(如图6所示),使用正弦定理可得多普勒速度与切向速度的关系式为

      图  6  切向速度与多普勒速度的几何关系

      Figure 6.  Relationship between the tangential velocity and Doppler velocity

      $${V_{\rm{D}}}(R,\alpha ) = R_{ci}^0 \cdot \frac{{{V_{\rm{t}}}(r)}}{r} \cdot \sin \left(\alpha - \alpha _{ci}^0\right)$$ (15)

      其中,$(R_{ci}^0,\alpha _{ci}^0)$为涡心位置,$i = 1,2$分别代表左涡心和右涡心,$r$为空间点$(R,\alpha )$到涡心的距离

      $$r = \sqrt {{{(R_{ci}^0)}^2} + {R^2} - 2R_{ci}^0R\cos \left(\alpha - \alpha _{ci}^0\right)} $$ (16)

      经典尾流速度模型有Hallock-Burnham速度模型、Rankine速度模型和Lamb-Oseen速度模型等[20],已证实Hallock-Burnham速度模型与实测数据吻合度较高[21],而Rankine模型和Lamb-Oseen模型存在对环量的低估[22],故采用经典尾流速度模型Hallock-Burnham速度模型[23]来计算尾流的切向速度

      $${V_{\rm{t}}}(r) = \frac{\varGamma }{{2{\rm{\pi }}r}} \cdot \frac{{{r^2}}}{{{r^2} + r_{\rm{c}}^2}}$$ (17)

      其中,${r_{\rm{c}}} = 0.052{b_0}$, b0为涡核大小。

      将切向速度关系式(17)代入多普勒速度关系式(15)中,可得多普勒速度${V_{\rm{D}}}(R,\alpha )$的表达式为

      $${V_{\rm{D}}} = \frac{{\varGamma R_{ci}^0}}{{2{\rm{\pi }}}}\frac{{\sin \left(\alpha - \alpha _{ci}^0\right)}}{{{r^2} + r_{\rm{c}}^2}}$$ (18)

      在雷达实际探测过程中得到的多普勒速度为尾流速度分量和背景风场分量的叠加,设计量测函数时需将背景风场分量去除,因此量测向量${{{\mathit{\boldsymbol{z}}}}_k}$设置为多个空间点去除背景风后的多普勒速度矢量,其中一个空间点$(R,\alpha )$的量测输出为

      $$\begin{split} {{\widehat z}_{k + 1}} =& h({x_{k + 1}}) = \frac{{{\varGamma _{k + 1,1}}{R_{c1}}}}{{2{\rm{\pi }}}}\frac{{\sin (\alpha - {\alpha _{c1}})}}{{{{r_1^2}} + r_{\rm{c}}^2}} \\ &+\frac{{{\varGamma _{k + 1,2}}{R_{c2}}}}{{2{\rm{\pi }}}}\frac{{\sin (\alpha - {\alpha _{c2}})}}{{{r_2^2} + r_{\rm{c}}^2}} - ( - {V_{\rm{d}}} + V_{\rm{c}}^{\rm{y}})\cdot\sin \alpha \\ & -(\beta {y_{k + 1}} + V_{\rm{c}}^0)\cdot\cos \alpha \\[-12pt] \end{split} $$ (19)

      其中,${R_{ci}} = \sqrt {x_{k + 1,i}^2 + y_{k + 1,i}^2} $ , ${\alpha _{ci}} = {\rm{arctan}}({y_{k + 1,i}}/ $$ {x_{k + 1,i}})$, ${r_i}$为空间点$(R,\alpha )$到涡心的距离,$i = 1,2$分别代表左涡心和右涡心。多个不同位置空间点的量测输出值${ {\widehat{z}}_{k + 1}}$即构成量测输出向量${\widehat {{\mathit{\boldsymbol{z}}}}_{k + 1|k}}$。同样为加强算法的鲁棒性,随机选择在RHI扫描平面上多个不同位置的空间点作为实测值并构成实测向量${{{\mathit{\boldsymbol{z}}}}_{k + 1}}$用于数据同化。

    • 本节通过仿真实验和实测数据对基于数据同化的飞机尾流演化趋势预测方法的有效性进行验证分析。首先构建飞机尾流仿真数据模块,随后使用两步定位[24]和路径积分方法[25]对预测需要用到的涡心位置和速度环量等尾流特征参数进行估计,进而对预测方法从仿真数据和实测数据两个方面进行验证分析。

    • 飞机尾流仿真数据模块根据预先设定的飞机参数模拟飞机尾流的形成以及在复杂大气环境中受风场影响演化运动的过程,主要包含复杂风场模拟模型、飞机尾流速度和演化模型以及激光雷达回波模型。

      复杂风场模拟模型包含横向风、风切变和大气湍流等,其中平均横向风是指风速的基准值,通常用来表明特定时间的风速平均值;风切变是指空间中任意两点间风向或风速的突然变化,低空情况(距离地面高度小于600 m)下可划分为水平风的水平切变、水平风的垂直切变和垂直风的切变[26]。大气湍流是介质的不规则运动,其参数极其复杂且无规律变化,可采用基于冯·卡门(Von Karman)频谱的拟小波方法来模拟随机湍流场[27,28]

      飞机尾流稳定段的速度可以简化为左右两个漩涡产生的速度场叠加的合速度,其演化情况受大气温度、湍流能量耗散率、风切变以及尾流的初始强度、初始位置等多方因素影响[8],演化模型采用经典两阶段(Probabilistic Two-Phase, P2P)尾流耗散模型[7]来模拟飞机尾流的强度和位置在复杂大气环境中的演化过程。

      激光雷达接收的是探测单元内大量随机分布的气溶胶的所有后向散射信号,忽略粒子之间的互耦,可认为每个探测单元的后向散射信号为所有气溶胶后向散射信号的相干叠加[29]。激光雷达信号处理把相邻N个采样回波信号补零到M点后再作M点的快速傅里叶变换(Fast Fourier Transform, FFT)得到多普勒频谱。为提升回波的信噪比,可以对多个脉冲重复周期(Pulse Repetition Time, PRT)之间的多普勒频谱进行累积,以提高激光雷达的探测能力[30],进一步用矩估计方法可从频谱的1阶矩和2阶矩得到回波信号的多普勒速度和多普勒谱宽等信息。更多关于飞机尾流仿真模块构建的细节可参考文献[31]。

    • 应用5.1节构建的激光雷达探测飞机尾流仿真数据模块,仿真总时长设置为90 s,左右尾流涡的初始速度环量都设置为${\varGamma _0} = 400$ ${{\rm{m}}^2}/{\rm{s}}$,左右涡心的初始位置分别为${O_{c1}}({t_0}) = (550,107)$ ${\rm{m}}$${O_{c2}}({t_0}) = $$ (610,105)$ ${\rm{m}}$,背景风场包含横向风$V_{\rm{c}}^0 = - 5$ ${\rm{m/s}}$、风切变变化率$\beta = 0.05$、纵向风$V_{\rm{c}}^{\rm{y}} = - 0.3$ ${\rm{m/s}}$和大气湍流${\rm{EDR}} = 0.05$ ${{\rm{m}}^2}/{{\rm{s}}^3}$。激光雷达探测参数设置如表1所示。

      表 1  激光雷达探测参数设置

      Table 1.  Detection parameters of the Lidar

      主要参数量值
      雷达波长1.55 μm
      脉冲宽度170 ns
      时间窗长度120 ns
      距离门宽度21 m
      采样率50 MHz
      脉冲积累数1500
      信号噪声比–5 dB
      扫描速度2°/s
      扫描范围0~15°

      对仿真飞机尾流回波数据应用两部涡心定位方法和路径积分方法进行飞机尾流涡心定位和速度环量的计算,将定位和环量计算结果作为激光雷达探测数据输入,对DA方法进行验证,同时加入数据融合[16](Data Fusion, DS)方法预测结果作为对比。DS方法采用经典D2P模型对尾流进行初步短时演化,同时根据实际探测值应用线性卡尔曼滤波方法,将尾流速度环量、垂直距离和水平距离的估计误差(实际探测值与经典D2P模型预测值之差)作为修正量,建立关于该修正量的传递矩阵并进行传递预测,采用泰勒级数展开方法计算误差转移矩阵,当有测量数据时计算相应的卡尔曼增益系数,通过迭代方法将模型预测数据与实测数据进行融合得到修正量,应用修正量对经典D2P理论的预测结果和误差传递矩阵进行修正,数据融合流程如图7所示,具体方法可见文献[16]。

      图  7  数据融合方法流程

      Figure 7.  Flowchart of data fusion method

      经典飞机尾流预测模型根据初始设置气象环境参数和飞机参数进行尾流行为预测,为检验本文提出的数据同化方法对行为预测气象环境参数估计错误修正的有效性和鲁棒性,分别对经典飞机尾流预测模型初始气象环境参数的影响进行分析。

      (1) 涡心位置水平位置演化分析

      横向风$V_{\rm{c}}^0$由正确的–5 ${\rm{m/s}}$设置为错误的–7 ${\rm{m/s}}$,预测结果如图8所示,因为输入错误的气象环境参数,经典飞机尾流预测模型对尾流涡心演化水平位置的预测随着时间演化误差不断变大,DA方法和DS方法都能够根据实时探测得到的涡心位置进行修正,得到更加贴近真实涡心演化趋势的行为预测,但DA方法得到的行为预测曲线相比数据融合方法更加平缓,更加贴近真实值,而DS方法因为D2P预测误差不断加大导致修正效果不佳。

      (2) 涡心位置垂直位置演化分析

      纵向风$V_{\rm{c}}^{\rm{y}}$由正确的–0.3 ${\rm{m/s}}$设置为错误的–0.5 ${\rm{m/s}}$,预测结果如图9所示,与涡心水平位置的演化趋势一致,因为输入错误的气象环境参数,经典飞机尾流预测模型对尾流涡心演化垂直位置的预测随着时间演化误差不断变大,DA方法和DS方法都能够根据实时探测得到的涡心位置进行修正,得到更加贴近真实涡心和环量演化趋势的行为预测,但DA方法得到的行为预测曲线相比DS方法更加贴近真实值。

      图  8  飞机尾流行为预测方法对比(水平方向轨迹)

      Figure 8.  Comparison between different wake vortex behavior prediction (Horizontal trajectories)

      图  9  飞机尾流行为预测方法对比(垂直方向轨迹)

      Figure 9.  Comparison between different wake vortex behavior prediction (Vertical trajectories)

      (3) 速度环量演化分析

      湍流耗散率EDR由正确的0.05 ${{\rm{m}}^2}/{{\rm{s}}^3}$设置为错误的0.01 ${{\rm{m}}^2}/{{\rm{s}}^3}$,预测结果如图10所示,因为输入错误的气象环境参数,经典飞机尾流预测模型对尾流速度环量的预测存在一个固定的误差,DA方法和DS方法都能够根据实时探测得到的涡心位置和速度环量估计进行修正,但DA方法能够更快建立预测模型,得到的行为预测曲线相比DS方法更加贴近真实值。

      进一步对预测的误差进行对比,定义预测值与真实值之间的相对误差,计算公式如下:

      $${E_{ri}} = \frac{1}{N}\sum\limits_{j = 1}^N {\frac{{|{P_{ci}} - P_{ci}^{\rm{T}}|}}{Q}} \times 100\% ,\quad i = 1,2$$ (20)

      其中,${P_{ci}}$是一个RHI的估计位置或环量值,$P_{ci}^{\rm{T}}$为其对应的理论值,$N$为总的RHI数量,$i = 1,2$分别对应左涡和右涡。当计算涡心位置误差时,$Q = $$ {b_0}$;当计算环量误差时,$Q = {\varGamma _0}$

      图  10  飞机尾流行为预测方法对比(速度环量)

      Figure 10.  Comparison between different wake vortex behavior prediction (Circulation)

      设置不同的错误初始气象环境参数,比较DA方法和DS方法应用实测数据进行实时修正阶段的相对误差,涡心位置相对误差结果如表2所示,可以发现横向风误差对经典尾流预测模型影响较大,随着横向风误差变大,DA方法以及DS方法两种行为预测模型误差也会变大;纵向风误差和湍流耗散率误差及其误差大小的变化对行为预测影响较小,两种行为预测模型都能够达到较好的修正效果;在不同的误差设置条件下,相比DS方法,DA方法能够根据雷达实际探测数据进行修正并取得更好的效果。

      表 2  飞机尾流预测位置相对误差

      Table 2.  Relative error of predict trajectories

      主要参数DS method (%)DA method (%)
      横向风
      $V_{\rm{c}}^0$ (${\rm{m/s}}$)
      切变率$\beta $纵向风
      $V_{\rm{c}}^{\rm{y}}$ (${\rm{m/s}}$)
      ${E_{r1}}$${E_{r2}}$${E_{r1}}$${E_{r2}}$
      –70.05–0.317.7419.631.121.35
      –100.05–0.326.8328.522.642.41
      –50.01–0.33.423.141.981.76
      –50.10–0.34.164.712.862.17
      –50.05-0.53.593.211.231.01
      –50.05–0.86.876.521.361.71

      环量估计相对误差如表3所示,可以发现不同的湍流误差条件下,DA方法和DS方法都能够保持良好的修正效果,这是由于湍流耗散率的误差会导致经典预测模型的误差始终保持在恒定的误差值,经过一段时间的迭代,两种行为预测方法都能够将该恒定误差进行有效的修正;同样可以发现DA方法能够较快修正并得到反映真实反演趋势的传递函数。综上分析,DA方法能够较快建立与尾流真实演化趋势一致的传递函数,行为预测结果鲁棒性较高。

      表 3  飞机尾流速度环量相对误差

      Table 3.  Relative error of wake vortex circulation

      主要参数
      EDR (${{\rm{m}}^2}/{{\rm{s}}^3}$)
      DS method (%)DA method (%)
      ${E_{r1}}$${E_{r2}}$${E_{r1}}$${E_{r2}}$
      0.013.243.581.121.31
      0.033.263.751.251.02
      0.083.042.951.171.16
      0.103.413.251.081.29
    • 香港天文台使用多普勒激光雷达在香港国际机场开展了一系列的飞机尾流探测实验,通过将雷达放置在机场不同位置进行长时间的连续探测,验证了激光雷达探测晴空飞机尾流的有效性。本小节使用在香港国际机场应用激光雷达探测得到的多普勒回波数据对数据同化预测方法的有效性进行分析,探测使用的激光雷达型号为windcube 200s,雷达在垂直于跑道的平面上做俯仰向扫描并接收得到该平面的RHI多普勒回波分布图,探测系统参数设置如表4所示,更多探测场景细节描述可以参考文献[32]。

      表 4  香港机场激光雷达探测参数设置

      Table 4.  Detection parameters of the Lidar in Hong Kong field campaigns

      主要参数量值
      雷达波长(μm)1.54
      距离门宽度(m)25
      脉冲宽度(ns)200
      脉冲重复频率(kHz)20
      探测距离(m)50~6000

      (1) 探测场景1

      雷达放置在香港机场附近的亚洲博览馆天台,距离地面高度约为19 m,与机场北跑道垂直距离约为400 m,雷达波束扫描角度范围0.83°~10.77°,扫描角速度1.87°/s,探测场景如图11所示。

      图  11  香港机场北跑道激光雷达探测场景设置2014

      Figure 11.  Geometry setup of the observation in north runway of Hong Kong international airport, 2014

      应用两步涡心定位和路径积分方法反演得到飞机尾流涡心位置和速度环量等特征参数,并对得到的速度环量与前一相邻时间点的速度环量的变化斜率进行判断,斜率绝对值超过阈值时(大于1)的数据点舍弃不用,以滤除应用路径积分方法出现的反演奇异值点。进而将得到的特征参数用于验证飞机尾流行为预测方法,结果如图12所示,可以发现DA方法和DS方法都能够根据实测数据对尾流演化趋势进行修正,但DS方法基于线性卡尔曼滤波模型,对于非线性雷达探测数据,估计出现非常大的偏差,而DA方法能够根据实测数据进行非线性模型修正,行为预测演化趋势与雷达探测值更加贴近。DS方法修正曲线变化剧烈且偏离雷达探测值,与其使用经典D2P模型进行预测时依赖模型初始气象环境参数设置,与真实风场气象环境差异较大,进而导致模型的传递矩阵变化剧烈,无法利用线性卡尔曼滤波进行有效修正有关。

      图  12  飞机尾流行为预测方法对比

      Figure 12.  Comparison between different wake vortex behavior prediction

      (2) 探测场景2

      雷达放置在机场跑道附近,距离地面高度约为7 m,与机场南跑道垂直距离约为240 m,雷达波束扫描角度范围0~50°,扫描角速度5°/s,探测场景如图13所示。

      图  13  香港机场南跑道激光雷达探测场景设置2018

      Figure 13.  Geometry setup of the observation in south runway of Hong Kong international airport, 2018

      同样应用两步涡心定位和路径积分方法反演得到飞机尾流涡心位置和速度环量等特征参数,将得到的特征参数用于验证飞机尾流行为预测方法,结果如图14所示,同样可以发现DA方法预测效果远远好于DS方法且与雷达探测值接近。与探测场景1得到的结果类似,DS方法速度环量值的估计变化非常激烈且偏离探测值,与预测模型的传递矩阵不能匹配实际气象环境变化,且探测值较少不能通过多次迭代进行有效修正有关。相比DS方法,DA方法能够通过少量的观测值建立稳定的尾流演化趋势模型,特别是在尾流靠近地面时,由于地面效应的影响两个尾涡在垂直方向上出现不同程度的反弹,DA方法能够拟合其变化趋势并修正涡心位置,算法的有效性和鲁棒性更高。

      图  14  飞机尾流行为预测方法对比

      Figure 14.  Comparison between different wake vortex behavior prediction

    • 本文通过数据同化方法研究了飞机尾流行为预测问题,根据非线性卡尔曼滤波模型建立了飞机尾流行为预测模型,通过雷达实时探测数据对飞机尾流行为预测轨迹进行实时修正。结合经典飞机尾流预测模型,根据风场线性切变和最小二乘曲线拟合方法,构建了参数化飞机尾流预测模型,解决了经典模型初始化气象参数未随时间演化进行实时调整的问题。在数值仿真数据和实测数据验证中,建立了飞机尾流雷达回波仿真数据模块,通过误差对比分析,验证了数据同化飞机尾流行为预测方法的有效性和鲁棒性。

      本文的工作是基于数据同化方法实现飞机尾流行为预测的尝试,为针对国内实际地形和气象环境条件制定适合中国民用机场的飞机起降动态间隔标准提供了依据,为提高中国机场航班准点率,节省燃油消耗开支,显著提升社会效益等方面提供了技术支撑和理论基础。未来将进一步验证数据同化模型在不同雷达实际工作条件下的有效性,以实现更加复杂风场条件下的飞机尾流行为预测分析和危害评估。

参考文献 (32)

目录

    /

    返回文章
    返回