1 引言
低空风切变一般是指在距离地面600 m高度内， 风速或风向发生剧烈变化， 在近地层存在显著风矢量差异的现象（周建华， 2011）， 具有尺度小、 时间短、 突发性强等特点， 是对飞行安全危害较大的天气现象之一（Deshpande et al， 1991）。飞机在起飞和着陆阶段遭遇严重的低空风切变时， 剧烈的风向与风速变化将导致飞机飞行升力发生突然改变， 若不能及时调整可能引发飞机失速坠地、 冲出跑道等灾难性事故， 根据国际民用航空组织（International Civil Aviation Organization， ICAO）2023年的安全报告， 2022年共有64起飞行事故， 其中25起由湍流、 风切变或雷暴等天气现象引起， 约占事故总数的39%。兰州中川国际机场位于兰州市中心西北方向约75 km处， 是西北地区的重要机场之一， 拥有一条南北向的跑道。由于机场平均海拔大约2 km， 属于高原机场， 并且位于秦王川盆地中部， 周围被3 km高山环绕。机场北面和东面地势平坦， 西面和南面沟壑纵横， 导致风场变化迅速， 对兰州机场的航班运行构成一定的威胁（党冰等， 2012）。随着航班数量的不断增加， 低空风切变对飞行安全的影响愈发显著， 精准预报低空风切变的需求不断增加， 对保障飞行与人民生命财产安全具有重要意义。
各大机场逐步配备了相关探测设备以应对风切变， 主要包括超声波测风仪、 多普勒天气雷达、 风廓线雷达和激光测风雷达。虽然这些先进设备在风场探测方面具有优势， 但对低空风切变的精准观测仍有难度， 同时低空风切变的预报也是一个难点。目前， 风切变的准确预报已成为中小尺度气象学研究的热点课题。随着观测技术和计算条件的改善， 除了使用多普勒天气雷达和风廓线雷达观测低空风切变外， 数值模拟也被广泛应用于机场风切变的研究中。Hon（2020）针对某国际机场在连续两年间内利用数值天气预报（Numerical Weather Prediction， NWP）模式预测低空风切变的表现进行了统计评估， 证明数值模式具备模拟机场低空风切变的性能。Basu and Storm（2010）将WRF模式应用于对美国大平原（The United States Great Plains， USGP）上任意位置处风切变指数的估算任务中。结果表明， WRF模型能定性地捕捉部分低空风切变的特征， 但模式对于大气边界层下部的物理参数化方案有待改进。翁雪玲等（2020）利用WRF模式对大连机场3个由偏北大风造成的典型低空风切变事件进行了模拟。结果显示， WRF模式模拟得到的风速数值小于观测数据， 不能表现出风切变发生时的风速和风向突变特征， 但可以通过计算风速差来预测风切变的发生时刻。
前人研究已经证明数值模式对风切变研究的可行性， 但目前国内外在风切变模拟应用方面的研究较少， 且多侧重于探讨风切变的形成机理。利用中尺度数值模式可以研究低空风切变的形成机理， 但受模拟精度及模式分辨率的限制， 无法对低空风切变发生过程中的小尺度风场进行准确模拟。因此， 许多学者采用计算流体动力学模型结合中尺度数值模式对小尺度风场进行模拟。李磊等（2010）使用CFD模式软件FLUENT对复杂地形风场进行了精细模拟研究， 结果显示， FLUENT较中尺度模式能更精确地刻画小范围复杂地形特征， 能够在小尺度范围内得到更高分辨率的模拟结果。Xiong et al（2022）将CFD对高层建筑周围典型风场的模拟结果与风洞试验的仿真结果进行比对， 验证基于多面体网格的CFD仿真精度， 发现通过改变k-ε湍流模型能有效改善CFD在不同风速条件下的仿真效果。张嘉荣和程雪玲（2020）改进了CFD模型， 采用大涡模拟（Large Eddy Simulation， LES）湍流方案， 并在大涡模拟中加入了温度层结的浮力模型， 发现WRF模式与Fluent耦合使得模式对风向与风速的模拟效果得到了改善。
大涡模拟因能够细致描述湍流结构而广泛应用于小尺度过程的研究中。该方法可以通过WRF-CFD耦合方案或在WRF模式中使用子网格模块实现。孙学金等（2017）将WRF模式最内层设置为百米级多层嵌套网格， 在新疆博斯腾湖周围区域模拟研究时段内风场变化， 结果表明精细化地形数据对WRF-LES模拟结果有重要影响， 而地表粗糙度对WRF-LES近地面风速模拟的影响机理需要进一步研究。张珊等（2023）采用WRF-LES大涡模拟试验， 研究2022年北京冬奥会张家口崇礼地区附近山地复杂地形， 通过分析研究区域风场， 发现WRF-LES能够很好地模拟出局地风场的显著变化特征。尽管已经有许多针对WRF-CFD和WRF-LES耦合预报方案的研究， 但目前对这两种方案得到的小尺度精细化风场结果尚未展开系统比较。二者都可以模拟小尺度风场， 但对风切变模拟结果的差异性研究较少。
传统数值天气预报模式由于计算效率低、 消耗时间长、 模型适应性差等问题， 未能实际应用于机场风切变的预报预警中。近年来， 随着大数据及人工智能技术的发展， 基于数据驱动的人工智能技术作为改进传统预报方式的手段， 逐渐应用于各类气象要素的预测任务中。长短期记忆网络是一种时间循环神经网络， 被广泛应用于回归和分类分析任务中。已经有部分研究将LSTM应用于风速预测中， 如Huang et al（2024）采用基于信息熵的熵集成经验模式分解（entropy ensemble empirical mode decomposition， eEEMD）方法与LSTM耦合的方式， 针对具有复杂时间序列特征的风电数据展开风电预测的研究。但现有的研究大多只针对某一固定区域， 缺乏对模型在不同地形和环流背景条件下风速数据泛化能力的讨论。Irrgang et al（2020）也指出了这一不足， 同时缺乏对WRF模式与LSTM耦合预报方式的研究。
基于以上论述， 本文在使用WRF模式对机场风场进行中尺度模拟的基础上， 通过比较WRF-LES和WRF-CFD耦合方案的结果， 探讨两种耦合方案的优劣。同时构建长短期记忆网络开展单变量与多变量时间序列风速预测， 以讨论LSTM在风速预测中的优势， 期望为提高风场模拟精度， 缩短精细风场模拟时间提供新的视角和方法。
2 案例及实验设计介绍
2.1 案例分析
2021年4月15日06:31（世界时， 下同）， 兰州中川机场飞机使用36号跑道着陆时， 在距36号跑道入口水平1.0 km垂直91.44 m处遭遇低空风切变， 触发低空风切变警告； 16日03:30， 兰州中川机场飞机使用18号跑道起飞时， 在全跑道至18号跑道向北3.0 km范围内152.40~752.80 m高度遭遇低空风切变。两次低空风切变具体信息如表1所示， 连续两天的两次低空风切变， 对飞行业务与飞行安全造成重大威胁。
表1 低空风切变个例基本信息Table 1 Basic Information of Low-level Windshear Cases |
| 日期（年-月-日） | 世界时 | 位置 | 高度/m | 起降阶段 |
|---|---|---|---|---|
| 2021-04-15 | 06:31 | 距36号跑道入口1.0 km | 91.44 | 着陆 |
| 2021-04-16 | 03: 30 | 全跑道至18号跑道向北3.0 km | 152.40~752.80 | 起飞 |
利用ERA5再分析数据给出2021年4月15日和16日低空风切变发生时500 hPa高空和地面天气形势， 如图1所示［该图及文中所涉及的地图均是基于中华人民共和国自然资源部审图号为GS（2019）1822号的标准地图制作， 底图无修改］。4月15日06:00， 中川机场位于内蒙古高原上空高空冷槽槽后， 机场高空为偏西风， 风速达到22 m·s-1以上， 而机场地面主要受冷锋后西侧高压控制， 风向以偏北风为主， 锋区活动导致气温、 风速及风向的迅速变化为风切变的产生创造有利的环境条件； 而16日03:00， 高空冷槽加深发展为切断低压并向东移动， 机场上空风速约为18 m·s-1， 风向仍为偏西风， 机场周围风速切变更加明显， 地面图显示机场西侧高压移动到机场正北方， 高压势力加强， 机场东南方等压线稀疏， 风速较小。总体来说， 在两次低空风切变发生过程中， 机场上空主要受北方高空槽后西北风的影响， 风向以偏西风为主且风速较大， 地面受北部高压控制以偏北风为主， 风速较小， 该天气形势有利于高空动量下传， 为机场低空风切变的发生提供有利条件。
图1 2021年4月15日06:00（a， c）和16日03:00（b， d）500 hPa（a， b）与地面（c， d）天气图蓝色实线为等位势高度线， 红色线为等温线， 黑色风羽指代风向风速， 红色五角星指代中川机场位置， 红色虚线框内为影响机场的主要天气系统， 棕色线表示槽线 Fig.1 Weather charts at 500 hPa （a， b） and the ground （c， d） for 06:00 on 15 （a， c） and 03:00 on 16 （b， d） April， 2021.The blue solid lines representing contours of geopotential height， red lines representing isotherms， black wind barbs indicating wind direction and speed， a red five-pointed star indicating the location of Zhongchuan Airport， and the major weather system affecting the airport enclosed within the red dashed box， The brown line represents the trough line |
2.2 数据来源和研究区概况
兰州中川机场（36.52°N， 103.62°E）位于中国甘肃省兰州市皋兰县中川镇， 距兰州市中心约75 km， 位于青藏高原东北侧与黄土高原交界处， 是中国西北地区主干机场之一。利用先进星载热发射和反射辐射仪全球数字高程模型（Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model， ASTER GDEM）绘制机场周围地形［图2（a）］， 可以看出机场位于秦王川盆地东南部， 周围有3000 m高山环绕， 地势主要为北高南低， 北侧和东侧地形平坦开阔， 西侧和南侧地形沟壑纵横。图2（b）显示了机场的跑道分布， 中川机场有一条南北走向（18/36）的跑道， 跑道长4 km， 宽45 m。由于机场周围地形复杂， 导致其风场复杂多变； 同时该机场平均海拔为1947.2 m， 属于高原机场， 受太阳辐射影响大， 升温较快， 上冷下热容易形成对流， 导致风场的变化， 有利于低空风切变的发生。
图2 原始地形tif文件研究区域范围（a）、 地形（彩色区， 单位： m）及观测站点分布状况（b）黑色断线框表示机场跑道范围， 红色圆点为观测站点位置， 五角星为中川机场中心位置 Fig.2 Original terrain tif file about the scope of the study area （a） and the topographic （color area， unit： m） and observation station distribution in the study area （b）.The black dashed box indicates the airport runway area， black dots indicate the locations of observation stations， and the five-pointed star represents the central location of Zhongchuan Airport |
为测试模型在不同地形及环流背景下的预测精度， 本文选取机场附近自东向西北绕城2号、 MID1、 二次雷达站、 18及36五个观测站点在2021年4月15日06:30与4月16日03:30的观测数据进行分析， 观测站点的位置分布和经纬度坐标如表2所示。
表2 观测站点名称与经纬度位置Table 2 Names and latitude and longitude locations of observation stations |
| 序号 | 站点名称 | 经纬度坐标 |
|---|---|---|
| 1 | 北绕城2号 | 103.647°E， 36.553°N |
| 2 | MID1 | 103.622 °E， 36.513°N |
| 3 | 二次雷达站 | 103.596°E， 36.500°N |
| 4 | 18 | 103.622°E， 36.528°N |
| 5 | 36 | 103.622°E， 36.496°N |
利用机场的观测站点数据来对两次低空风切变发生前后的风向风速变化进行分析， 如图3所示， 15日低空风切变发生时， 可以看到MID1站点风速发生突变， 风速由2 m·s-1增加到8 m·s-1， 风向在15日05:00前主要为偏西风， 随后风向角度逐渐减小， 在06:00左右转变为偏北风， 最后于低空风切变发生时转变为偏北风； 16日低空风切变发生时， 风速也有突变， 由2 m·s-1增大为6 m·s-1， 相较于15日低空风切变， 风速变化较小， 但是总体上的风速增加较大， 同时其风向突变更加明显， 由偏东风快速突变为偏北风。两次低空风切变的发生存在一些相同特征： 在低空风切变发生时均存在风速的迅速增加， 同时风向也发生改变， 由偏东风转变为偏北风。
图3 2021年4月15日00:00至17日00:00 MID1站点的风速（a）、 风向（b）变化黑色实线为观测数据， 红色虚线指代低空风切变发生时间 Fig.3 Variations in wind speed （a） and direction （b） at MID1 Station from 00:00 on April 15 to 00:00 on April 17， 2021.Solid black lines represent observed data， dashed red lines indicate the occurrence time of low-level windshear |
2.3 模式参数设置及试验设计
2.3.1 WRF模式
本文选用WRF-ARW V4.2模式来对中尺度风场进行模拟， WRF模式的优势之一是使用嵌套网格来提高模拟精度与计算效率， 本文使用四层嵌套网格对研究区域进行模拟， 网格间分辨率比率设置为3， 即网格分辨率分别为27 km、 9 km、 3 km和1 km， 模拟时段为2021年4月14日00:00至4月17日00:00， 覆盖两次低空风切变的发生时间， 时间积分步长与嵌套网格相对应， 分别为108 s、 36 s、 12 s和4 s， 以时间间隔6 h的NCEP FNL数据作为边界条件输入； 此外， 为了比较嵌套网格分辨率提升对模式模拟结果的影响， 对每一层嵌套网格格点数都设置为100×100， 并且在垂直方向设置为36层， 四层嵌套网格中心点重合， 中心点设置为（36.512°N， 103.620°E）， 即中川机场中心点位置， 网格设置如图4所示。
由于模式网格分辨率有限， WRF模式对于小于网格点的物理过程无法模拟， 但可以利用多种物理过程参数化方案来表示次网格物理过程。物理过程方案设置如下： 微物理过程采用Thompson方案； 积云对流过程采用Kain-Fritsch方案对云内水汽上升下降过程、 相对粗糙的微物理过程进行描述； 近地面层为与边界层方案配合而选用经修订的 MM5 表层方案（Revised MM5 Monin-Obukhov）； 陆面过程采用Noah方案， 如表3所示。结合前人关于低空风切变特征的研究（Sukanta and Brandon， 2010）， 中川机场风切变主要受到复杂地形及太阳辐射的影响， 因此本文主要针对边界层方案和辐射方案进行参数化方案选取。王腾蛟等（2013）的研究表明， 对于黄土高原丘陵地形条件下大气边界层进行模拟时， MYJ方案在稳定边界层条件下模拟效果良好， 而YSU方案则能更好地实现对流边界层结构特征的模拟， 考虑到中川机场周围具有相似的地形地貌特征， 因此本文选择两种边界层参数化方案YSU（Hong et al， 2006）、 MYJ（Mellor and Yamada， 1982）， 以及两种辐射参数化方案RRTM（Mlawer et al， 1997）和Goddard方案（Chou and Suarez， 1999）分别进行组合， 形成四个参数化方案组合开展参数化方案敏感性试验， 如表4所示。
表3 WRF模式参数化方案设置Table 3 WRF mode parameterization scheme Settings |
| 物理过程 | 参数化方案 | 备注 |
|---|---|---|
| 微物理过程（mp_physics） | Thompson （Thompson et al， 2008） | 包含冰、 雪和霰过程的适用于高分辨率模拟的方案。 |
| 积云对流过程（cu_physics） | Kain-Fritsch （Kain， 2004） | 使用有下沉气流的深对流和浅对流次网格方案。 |
| 近地面层过程（sf_sfclay_physics） | Revised MM5 Monin-Obukhov （Jiménez et al， 2012） | 基于含有Carslon-Boland粘性底层和来自查表中的标准相似函数的Monin-Obukhov方案 |
| 陆面过程（sf_surface_physics） | Noah （Chen and Dudhia， 2001） | 统一NCEP / NCAR / AFWA方案， 含有4个不同层次上土壤温度和湿度、 积雪覆盖面积和冻土物理过程。 |
表4 WRF模式物理过程参数化方案组合Table 4 Physical process parameterization schemes in WRF mode |
| WRF参数化 方案组合 | 辐射方案 | 边界层 方案 | |
|---|---|---|---|
| 长波辐射 | 短波辐射 | ||
| S1 | RRTM | Dudhai | MYJ |
| S2 | RRTM | Dudhai | YSU |
| S3 | Goddard | Goddard | MYJ |
| S4 | Goddard | Goddard | YSU |
为客观评价数值天气预报模式WRF在模拟风向与风速方面的效果， 采用气象监测站点的实测资料对模拟结果进行精度评价。本文选择四个指标参数来对中尺度风场模拟结果进行评价： 均方根误差（Root Mean Squared Error， RMSE）， 平均绝对误差（Mean Absolute Error， MAE）， 一致性指数（Index of Agreement， IA）和相关系数（Correlation Coefficient， R）。
均方根误差是一种常用的衡量模型预测风速值与实际观测风速值之间差异的指标， 它用于评估风速模拟的误差； 平均绝对误差是所有单个观测值与算术平均值的偏差的绝对值的平均， 平均绝对误差可以避免误差相互抵消的问题， 因此能够准确反映实际预测误差的大小； 一致性指数和相关系数用于评价模式的模拟能力。四个参数指标的计算公式如表5所示， 其中S指代模拟数据； O指代观测数据； 下标i指代时间序列值； N指代总体时间序列样本量。
表5 评价指标Table 5 Evaluation indicators |
| 名称 | 公式 |
|---|---|
| 均方根误差（RMSE） | |
| 平均绝对误差（MAE） | |
| 一致性指数（IA） | |
| 相关系数（R） |
2.3.2 CFD模式
计算流体动力学是一种利用数值分析和算法来解决涉及流体流动的问题的科学， CFD 模式指的是在计算流体动力学中使用的数值方法和算法模型， 用于模拟和预测流体在既定初始条件下的运动特征（李绍辉， 2020）。现有的部分研究结果表明， 采用 WRF 模式的输出风场作为 CFD 模式的边界条件， 可以很好地模拟实际大气中的复杂地形区域的微尺度风场（Temel et al， 2018）， 而目前 WRF 模式也可以调用模式自带的LES模块来实现大涡模拟（李叶晴， 2023）， 从而实现微尺度风场的模拟， 因此本文针对两种降尺度风场耦合预报方案展开对比研究。
对于小尺度风场， 本文采用CFD模式的开源软件OpenFOAM与前期WRF模式四层嵌套网格结果进行耦合模拟。WRF-CFD耦合模拟使用高精度的ASTER GDEM地形数据构建CFD模式中的研究区域， 以机场中心向外6 km选择一个12 km×12 km的区域（图2）， 利用建模软件GlobalMapper、 SketchUp与OpenFOAM软件前处理工具SnappyHexmesh得到三维地形网格， 然后利用OpenFOAM的blockMesh工具建立体网格， 网格范围为12 km×12 km×1 km， 网格分辨率设置为100 m。此外， CFD模式的湍流求解方案设置为LES大涡模拟并且选择WALE方案。以WRF最优参数化组合对2021年4月15日00:00至4月17日00:00风场每隔5 min的模拟结果作为边界条件输入CFD模式。由于两次低空风切变发生前后机场上空主要为东北风， 因此选用东面、 北面和顶面作为入流面， 而西面和南面作为出流面， 地面选择为壁面无滑移条件。设置模式运行时间与WRF数据的输入时间一致（24 h）， 模式运行文件按照30 min输出一次。
2.3.3 LSTM风速预测模型
传统的循环神经网络（Recurrent Neural Network， RNN）主要用于处理基于时间序列的数据， 但随着时间推移会对先前的信息有所遗忘， 容易出现梯度消失或梯度爆炸等问题（赵敬娇等， 2021）。LSTM的基本单元较RNN更为复杂， 该模式单元结构包含： 输入门、 遗忘门、 一个或多个记忆单元和输出门（虞浩跃等， 2019）。门控机制表达式如下：
式中： 下标i， f，o分别指代输入门、 遗忘门、 输出门； t表示时间步； 为激活函数； 表示输入值； 指前一时间步的隐状态； ， ， 指输入与门之间的权重参数； ， ， 指隐藏单元与门之间权重参数； ， ， 为偏置参数； 表示前一时间步的细胞状态（Cell State）； 细胞状态与隐状态的计算方法分别如公式（4 ）、 （5 ）所示：
式中： 下标c表示记忆细胞； 为输入与记忆细胞间的权重参数； 指隐藏单元与记忆细胞间的权重参数。
通过引入记忆单元（Memory Cell）和门控机制（Gate Mechanism）， LSTM能够更好地捕捉序列数据中的长期依赖关系， 解决了RNN存在的梯度消失和梯度爆炸问题。同时， LSTM具备自动学习及处理非线性关系的能力， 在应对数据量足够大的复杂模式时， 表现出良好的预测性能。这使得LSTM在时间序列预测、 信号处理等任务中具有优势（胡庆芳等， 2020； Xiang et al， 2020）。考虑到风速数据受到多种环境因素的影响而呈现高度非线性的特征， 本文选用LSTM模型展开对于人工智能在低空风切变预报中应用的研究。
实际观测数据受到观测设备、 观测场地及保障需求等条件的限制， 得到的观测数据集中仅包含风速与风向数据， 但观测数据集采样间隔为15 s， 风速数据十分精细且数据量大。WRF模式输出数据囊括不同高度处位温、 气压、 湿度等多种不同种类的气象要素信息， 但因数据集时间间隔为5 min， 同时间段内数据量远少于观测数据。基于两种数据集的不同特征， 本研究首先将观测数据用于开展单变量时间序列风速预测， 以分析构建的LSTM模型在关于风速预测方面的可行性与准确性； 若LSTM模型在数据量较大的单变量预测中具有良好的表现， 再利用数据量较少的WRF模式输出数据集开展多变量风速预测的尝试， 即利用与风速相关的气象要素构建多特征预测算法对风速进行预测， 探讨基于长短期神经网络与WRF耦合的精细化风场预报方案的可行性。
单变量时间序列风速预测所用数据来自兰州中川国际机场周边二次雷达站、 北绕城二号、 MID1及18观测站点的实测风速数据， 利用变化趋势不同的风速数据表示不同地形及环流背景条件， 据此对模型的泛化性能进行评估。选取时间段为2021年4月15日00:00:00 -23:59:10， 采样间隔为15 s， 试验所用数据集为满足数据中风速序列具备波动性和周期性的特征以便LSTM模型进行预测时取得更好的效果， 研究中进行数据清洗后四个数据集可用数据依次为3881、 3904、 5760及5760， 考虑到4月15日低空风切变发生时次为06:31， 试验采用数据的前30%作为测试集， 剩余数据作为训练集。风速数据如图5所示。
考虑到在传统的数值天气预报模式中， 斜压原始方程组内包含的比湿、 气压、 温度等气象要素与风速之间存在一定的动力或热力关系； 同时， 由于位温在干绝热过程中为保守量， 与温度相比位温是一种稳定的示踪物， 常用于追溯气块或气流的源地及研究演变规律（盛裴轩等， 2003）。综合考虑上述两点原因， 本研究选用比湿、 气压、 地面2 m高度温度及扰动位温数据构建特征集用于开展基于LSTM多变量时间序列预测的研究。
多变量时间序列风速预测所用数据来自机场MID1跑道中部观测站点处三个嵌套网格（d02， d03， d04）下的WRF模式输出数据， WRF模式以FNL数据作为输入， 选取时间段为2021年4月14日00:00:00至4月16日23:55:12， 采样间隔为5 min， 数据集包含100 m高度处比湿、 气压、 扰动位温、 风速及地面2 m高度温度资料， 共含866组数据。因WRF数据输出时间间隔较大， 数据量较少， 为保证模型的准确性而提高训练集占比， 试验采用数据的前80%作为训练集， 剩余数据作为测试集。具体数据分布如图12（a）， （c）， （e）所示。
图6 WRF模式4种参数化方案组合风速模拟结果Fig.6 Combined wind speed simulation results of four parameterized schemes in WRF model |
图7 WRF模式四种参数化方案组合结果评价（d02， d03， d04指代嵌套网格层数）（a） RMSE， （b） MAE， （c） IA， （d） R Fig.7 Evaluation of the Combination Results of Four Parameterization Schemes in the WRF Model （d02， d03， d04 refer to the number）.（a） RMSE， （b） MAE， （c） IA， （d） R |
图8 WRF-LES和WRF-CFD在三个站点的10 m高度风速模拟结果对比（a）18号站点， （b）MID1站点， （c）36号站点 Fig.8 Comparison of wind speed simulation results of WRF-LES and WRF-CFD at 10 m height at three stations.（a） Site 18， （b） MID1 site， （c） Site 36 |
图9 2021年4月15日06:30和4月16日03:30两个时刻的WRF-CFD和WRF-LES的10 m高度水平风场对比红色五角星指代中川机场位置， 红色箭头指代观测站点风速； （a， b） WRF-CFD， （c， d） WRF-LES Fig.9 Comparison of WRF-CFD and WRF-LES horizontal wind fields at 10 m height at 06:30 on April 15 and 03:30 on April 16， 2021.Red five-pointed star refers to the location of Nakagawa Airport， red arrow refers to the wind speed at the observation station.（a， b） WRF-CFD， （c， d） WRF-LES） |
图10 2021年4月15日06:30 （a， c）和2021年4月16日03:30 （b， d）的WRF-CFD （a， b）和WRF-LES （c， d）的飞机起飞或着陆过程中的风速差对比虚线指代飞机飞行过程中的高度变化， 实线指代风速变化， 红色虚线代表遭遇低空风切变位置 Fig.10 Comparison of wind speed differences during aircraft takeoff or landing between WRF-CFD （a， b） and WRF-LES （c， d） for the times of 06:30 on April 15， 2021 （a， c）， and 03:30 on April 16， 2021 （b， d）.The dashed lines represent changes in altitude during the flight， solid lines represent changes in wind speed， red dashed lines indicate the positions where low-altitude wind shear is encountered |
图11 单变量时间序列风速预测结果与绝对误差（a~d）依次为二次雷达站、 北绕城二号、 MID1及18号站点 Fig.11 Wind speed prediction results and absolute errors of unit time series.（a~d） are the Secondary Radar Station， North Ring Road No.2， MID1， and site 18， respectively |
建立模型前为加快模型收敛速度， 同时使WRF模式输出数据不同特征具有可比性， 消除特征间存在的量纲差异， 采用线性归一化的方法对上述两种数据进行数据预处理， 将数据值映射到［0， 1］区间内。
式中： 为归一化后数据； 为原数据； 与 为原数据集最大值和最小值。
本文采用的隐藏层结构为1个LSTM层与1个Dense全连接层， 在第一隐层中定义50个神经元。模型学习算法采用自适应矩估计（Adam）和自适应学习率调整算法（Adadelta）， 模型训练迭代次数（Epochs）、 批处理大小（Batch size）和模型学习算法等超参数调优过程需要结合具体数据确定， 由于本研究中采用的数据格式相近， 因此采用网格搜索（Grid Search）的方式进行超参数选择。模型训练迭代次数选择范围为［20， 120］， 间隔10； 批处理大小选择范围为［4， 30）， 间隔2。
本研究采用平均绝对误差、 均方根误差、 一致性指数、 相关系数和平均绝对百分比误差（Symmetric Mean Absolute Percentage Error， SMAPE）五个指标评价模型性能， 其中MAE和RMSE用于反映预测结果的绝对误差， 从图5中可以观察到部分时段风速值较小， 因此选用SMAPE观察预测结果的相对误差可有效规避因真实值过小而导致误差较大的问题， 表达式如下：
式中： n为样本数量； 为模拟值； 为观测值。
3 模拟结果分析
3.1 WRF模拟效果分析
首先将四种参数化方案组合模拟出的风速数据和与机场跑道处实测数据进行对比， 如图6所示。从结果上看， 不同方案都模拟出了风速变化的基本趋势， 但各方案对于风速变化的响应与观测数据相比具有滞后性； 对于4月15日风速模拟结果， 不同方案间模拟结果差别较大， 方案S2模拟结果显著优于其他方案， 在风速极值的处理上更加接近观测数据； 对于4月16日的低空风切变， 不同方案模拟结果差异并不明显， 都模拟出了风速变化的趋势， 相对而言， S2对于风速的模拟结果更加接近观测数据， 但在16:00 - 24:00的时间段内效果较S1和S4方案差。
为了对模拟结果进行量化比较， 使用前文所述的四个指标参数对风速模拟结果进行评价， 比较不同模式设计方案对风速的模拟效果， 如图7所示。其中d02、 d03和d04指代不同嵌套网格层， 由于d01嵌套网格分辨率过低， 未对d01网格进行分析。从结果来看， 随着嵌套网格分辨率的增加， 模拟精度也随之提高， 最内层网格均方根误差和平均绝对误差最小， 而一致性指数和相关系数最高， 因此从模拟结果和观测数据的四个评价指标来看， 最内层网格模拟效果最好。其次可以比较相同辐射方案之间的模拟效果来对比边界层方案模拟效果， 从观测数据的评价结果来看， S1方案模拟效果较好， 该方案的模拟结果与观测数据相比相关性较强， 同时误差也小于其他方案； S2方案虽然能更好的处理风速极值， 但受到模拟结果滞后于观测数据的影响， S2方案整体的误差偏大。总体而言， 随着模式分辨率的提升， 风速模拟效果显著改善； 方案S1模拟效果最好， 模拟结果相关性强而且误差最小。
3.2 WRF-LES与WRF-CFD模拟效果分析
对机场跑道上的三个站点10 m高度处WRF-LES和WRF-CFD的风速模拟结果进行对比。图8显示， 从单个站点的风速变化模拟结果上看， WRF-CFD模拟结果较WRF模拟结果偏大， 这可能是因为WRF模式中使用的参数化方案及输入变量有很多， 而WRF-CFD耦合模式中输入到CFD模式中的边界条件变量仅为WRF模式中输出的风速风向数据， 未将其他变量对风速的影响考虑在内； 而WRF-LES模拟结果优于其他前两种方案， 因为WRF-LES在小于1 km的网格中使用了LES湍流参数化方案， 提升了小尺度湍流的模拟精度， 并且结合WRF模式自身的优势， 可以对辐射等物理过程中使用参数化方案， 因此WRF-LES模拟结果最好， 并且和观测数据风速一致性较好。总结来说， 对于单站点的风速模拟结果， WRF-LES由于考虑了多种物理过程及湍流参数化方案， 模拟效果优于仅考虑物理过程参数化方案的WRF模式以及在小尺度网格仅考虑湍流参数化方案的WRF-CFD耦合模式。
图9结果显示， 在机场的北面及机场附近区域， WRF-CFD模拟结果主要为偏北风， 在西南方风向偏转为东北风， 这些模拟结果与观测数据都有很好的对应， 在机场的南方及东南方向， 风速的模拟结果和观测数据有较大不同。上述结果可能由跑道南北两侧影响边界层风场的主要因素不同所导致： 跑道北侧临近山地， 风场主要受山脉影响； 南侧地势相对平坦且建筑物分布相对密集， 因此边界层风场主要受建筑物影响。彭珍和胡非（2006）在对边界层风场的研究中指出， 气象观测站周围近地面风场会受到周围建筑物的影响， 导致观测风向变得紊乱， 而本文构建的三维地形网格分辨率为100 m， 忽视了大部分建筑物的作用， 所以导致在机场南方及东南方的模拟效果较差。
对于WRF-LES大涡模拟结果， 虽然在机场跑道南方的模拟效果优于WRF-CFD， 但是整体上风场杂乱且对于跑道北部的模拟效果较差， 不能展现出风场的基本状况。通过对两次低空风切变过程的风场比较， 可以得出结论： WRF-CFD对风场的模拟效果较好， 与观测数据一致性较高， 但是由于CFD仅将WRF结果中的风速作为边界条件， 导致该方案风速模拟结果比较均一； 而WRF-LES可以考虑多种物理过程参数化方案， 对风场的模拟更加接近真实风场特征， 但是对比观测数据没有WRF-CFD模拟效果好， 因此， 主要对WRF-CFD模式的风场模拟结果进行分析。在第一次低空风切变发生时刻， 机场附近有风的逆转， 即由偏北风逆时针偏转为西北风， 机场南部风速较小， 且主要为东北风； 对于第二次低空风切变发生时刻， 机场周围主要是偏北风， 在机场东南方偏转为东北风， 通过该时刻的WRF-CFD模拟， 可以看到机场西部由于地形海拔高度的迅速降低， 风速较小的偏北风在经过时风速迅速增大， 于山脚附近产生局地强风， 该结果进一步说明了地形对于机场周围风速的重要影响。
接下来针对WRF-CFD与WRF-LES两种耦合方案对低空风切变的模拟效果进行分析（图10）， 本文中所讨论的两次低空风切变均是飞机在起飞或着陆过程中遭遇的， 因此模拟飞机起飞或着陆的飞行过程中的风速变化。表1给出了两次低空风切变过程的具体信息， 对于第一次低空风切变过程， 两种方案模拟结果都成功模拟出了低空风切变， 通过计算飞机飞行过程的逆风、 侧风和垂直风风速变化， 发现模拟结果都是侧风风速变化最大， 并且均存在侧风风速先增大后减小最后再增大的过程； 两种方案结果在风速模拟数值和低空风切变发生范围上存在较为明显的差异， WRF-CFD模式模拟出的最大风速差较小， 但更接近于飞机遭遇低空风切变的位置， 而WRF-LES模拟出的侧风风速差大， 但影响范围较小， 且远离飞机遭遇低空风切变的位置。
对于第二次低空风切变过程， 该次记录中飞机处于起飞过程中， 并且此次低空风切变对飞机的影响范围较大， 在此期间， 飞机高度在152.4~752.8 m都受到低空风切变影响， 但是从两次低空风切变的模拟效果上来看， 最大风速差都是飞机飞行过程中的逆风风速变化过程中， 对于影响范围更大的第二次低空风切变两个模式模拟出的风速差都达不到7.7 m·s-1的阈值， 即都没有模拟出低空风切变， 这一结果可能由多因素共同作用所导致。从上文WRF模拟结果上看， 第一次低空风切变发生时WRF模式模拟出的风速值大于第二次低空风切变， WRF模拟风速值作为模式输入数据会影响到两种方案的模拟效果， 使得第二次风速差的模拟小于第一次低空风切变； 模式的设计也会对模拟效果产生影响， WRF-LES通过选择不同物理参数化方案提升模拟效果， 而WRF-CFD模式的模拟效果更大程度上受到湍流方案及模式的入流面选择的影响， 因此后续可以通过提升两种模式的模拟效果进一步比较对于低空风切变的模拟效果。结合两次低空风切变过程的模拟结果， 可以认为WRF-CFD对于低空风切变的模拟效果优于WRF-LES方案， 模拟得到的低空风切变更加接近飞机遭遇低空风切变的位置， 但是两种模式对于第二次低空风切变的模拟效果都不好， 模拟得到的风速差都小于低空风切变的阈值， 仍有很大的改进空间。
3.3 LSTM风速预测结果分析
图11对4月16日12:00:00 -23:55:12的风速模拟结果表明， 单变量风速预测的结果均能在一定程度上反映出实际风速变化， 并且在处理风速值低于6 m·s-1的情况时表现良好。表6中模型对于四种不同地形及环流背景条件下的模拟结果误差分布状况大致相同， 证明在低风速（0~6 m·s-1）情景下具备相对较高的精度及较好的稳定性； 然而， LSTM在风速值大于6 m·s-1时预测误差较大， 如图11（b）中于第200个样本点附近真实值为测试集中的最大值， 但预测结果为邻近区域的极小值， 这可能与风速值大于6 m·s-1的风速样本在训练集中所占比例过少相关， 模型不能很好地捕捉较大风速值情况下风速的变化趋势， 这也造成了SMAPE总体偏大。上述结果表明， 该模型较好地把握了不同地形与环流背景条件下风速变化的非线性关系， 较好地处理了传统气象模式在模型适应性上存在的问题。
表6 单变量风速预测结果指标Table 6 Indicators of unit wind speed prediction |
| 二次雷达站 | 18 | 北绕城2号 | MID1 | |
|---|---|---|---|---|
| RMSE/（m·s-1） | 0.6709 | 0.9035 | 1.5274 | 0.8804 |
| MAE/（m·s-1） | 0.4099 | 0.5902 | 0.942 | 0.6622 |
| IA | 0.9496 | 0.9367 | 0.9187 | 0.946 |
| R | 0.9323 | 0.9634 | 0.9159 | 0.9632 |
| SMAPE | 13.4010 | 17.2266 | 20.1793 | 22.9070 |
国内外部分研究采用WRF与LSTM耦合的方式开展河水流量或降水的预测， 如Wang et al（2023）将WRF二维高分辨率（1 km）模式输出数据输入CNN-LSTM模型中对径流流量进行模拟， 证明CNN-LSTM方法与WRF耦合开展洪水流量模拟的可行性。但关于短期风速预测方面的研究多采用实际观测结果构建模型， 缺乏与WRF模式的结合（林海涛， 2010）。基于LSTM在单变量风速预测模拟结果能在一定程度上反映出实际风速状况， 本研究采用WRF-LSTM耦合方法展开多变量融合的预测尝试， 探讨基于LSTM与WRF耦合的精细化风场预报方案的可行性。对于4月16日09:30:00 - 23:55:10模拟结果显示， 长短期神经网络预测的风速曲线与测试集原始风速数据曲线具有高度的一致性， 表7中绝对误差基本维持在0.25 m·s-1以下且相对稳定， 在风速大于6 m·s-1的条件下该模型仍具有很好的表现， 如图12（c）中第100~175的样本点， 图12（e）中第0~50的样本点。同时， 融合多变量的风速预测模型能够很好地处理短时间内风速值突变的情况， 如图12（c）第75~100的样本点处风速值由5 m·s-1迅速增加至11 m·s-1， 而模型误差仅在初始变化时有上升趋势， 自样本点87处达到峰值0.5 m·s-1后下降。多变量风速预测模型在保证预测精度的情况下， 以更高的计算效率实现了对风速的预测。但从图12原始数据图像中可以观察到15日与16日的要素曲线具有一定的相似性， 模型模拟结果精度较高是否同这一原因相关有待进一步讨论。
表7 多变量风速预测结果指标Table 7 Result indicators of multiple wind speed prediction |
| d02 | d03 | d04 | |
|---|---|---|---|
| RMSE/（m·s-1） | 0.0078 | 0.018 | 0.0145 |
| MSE/（m·s-1） | 0.0065 | 0.0149 | 0.0108 |
| IA | 0.9989 | 0.9982 | 0.9988 |
| R | 0.9979 | 0.9975 | 0.9981 |
| SMAPE | 5.028 | 2.2717 | 6.6015 |
图13结果显示多变量风速预测结果与观测数据间存在较大差距， 但同WRF输出数据拟合较好， WRF模式输出的数据同实际观测资料的偏差及变量选择方式存在问题， 导致模型预测风速不能很好地表现出实际风场变化。综上所述， 基于LSTM的多变量时间序列风速预测模型在测试集上具有良好的表现， 但受观测数据中观测项目有限和WRF模式模拟效果欠佳的限制， 多变量融合模型现阶段仅可使用WRF模式输出数据进行预测且结果同实际风速间存在不小差距。通过改进WRF模式模拟方案或增加机场跑道附近观测项目， 多变量时间序列风速预测模型高精度和较传统大气数值预报模式预报所需时间短的实际意义才能得到很好的体现。
图13 WRF输出结果、 多变量风速预测结果与观测数据比较（a）为WRF输出结果、 多变量风速预测结果与观测数据曲线， （b）为WRF输出结果与观测数据间的误差及平滑曲线 Fig.13 Comparison of WRF output results， multivariate wind speed prediction results and observed data.（a） Curves of WRF output results， multivariate wind speed prediction results and observed data， （b） errors and smooth curves between WRF output results and observed data |
4 结论与讨论
本研究利用中尺度大气模式结合大涡模拟和Openfoam， 并尝试融合LSTM机器学习预报技术， 对2021年4月15日06:31飞机着陆过程和2021年4月16日03:30飞机起飞过程遭遇的两次低空风切变进行预测模拟分析。首先利用再分析数据对两次低空风切变过程的天气形势进行分析， 然后以中川机场为中心进行了多组WRF中尺度数值模拟， 得到最优物理过程参数化方案组合； 然后使用ASTER GDEM高分辨率地形数据提高地形精度， 采用WRF-LES大涡模拟和WRF-CFD耦合模式两种方法对机场跑道周围小尺度风场进行模拟， 并对两种模式的模拟结果进行对比； 最后利用MID1观测数据和WRF模式输出数据构建LSTM风速预测模型， 展开单变量与多变量风速预测研究， 得到以下主要结论：
（1） 低空风切变的发生与大尺度环流背景有关。在两次低空风切变发生过程中， 机场上空主要受北方高空槽后西北风的影响， 风向以偏西风为主且风速较大， 地面受北部高压控制以偏北风为主， 风速较小， 该天气形势有利于高空动量下传， 为机场低空风切变的发生提供有利条件。
（2） 随着嵌套网格分辨率的增加， WRF模式模拟结果均方根误差和平均绝对误差减小， 一致性系数和相关系数增大， 模拟效果随着分辨率的提高而变好； 通过调整辐射和边界层参数化方案， 得到模拟效果较好的方案组合（RRTM方案和MYJ方案）。
（3） 对单个站点风速模拟来说， WRF-LES模拟结果比较贴合观测数据， 而WRF-CFD模拟结果则偏大； 对于近地面水平风场来说， 由于考虑多种物理过程参数化方案， WRF-LES的模拟结果更类似真实风场， 但是模拟结果与观测站点差别较大， 而WRF-CFD耦合模式风速模拟结果比较均一， 但是模拟区域的模拟效果较好， 模拟风场与观测站点风速风向基本一致， 并且能够模拟出明显的小尺度涡旋特征。
（4） 对于飞机降落过程中遭遇的两次低空风切变的模拟， WRF-LES和WRF-CFD两种模式都可以模拟出第一次低空风切变， 并且都是侧风风速差最大； 而第二次低空风切变两种模式都是逆风风速差最大， 但是两种模式风速差都没有达到阈值， 都没有模拟出第二次低空风切变， 原因可能是因为传入两种模式的WRF风速数据值较小， 需要在后续工作中进一步验证。
（5） 基于LSTM的单变量风速预测模型能在一定程度上反映出实际风速的变化， 在风速值低于6 m·s-1的情况时表现良好， 并且具备处理不同地形及环流背景条件下风速数据的能力， 但不能很好地捕捉较大风速值情况下风速的变化趋势。多变量风速预测模型较单变量风速预测模型结果精度更高且更为稳定， 能很好地处理风速大值与风速快速变化的情况， 但受观测数据中观测项目有限和WRF模式模拟效果欠佳的限制， 该模型的实际价值有待进一步讨论。
（6） LSTM风速预测模型较传统大气数值预报模式在没有进行物理过程参数化设置的情况下， 该模型较好地把握了不同地形与环流背景条件下风速变化的非线性关系， 一定程度上解决了传统预报模式在模型适应性上存在的问题。同时该模型在保证预测精度的情况下， 以更高的计算效率实现了对风速的预测， 实际应用场景更加广泛。
本文着重分析WRF-CFD与WRF-LES耦合方案在风场及低空风切变预报中的差异性、 基于LSTM进行风速预测的可行性与准确性， 但还有一些科学问题需要深入探讨， 包括： WRF-LES和WRF-CFD两种耦合方案进行预测时误差的主要来源， 在对风场进行模拟时考虑更多的物理过程对结果是否有正面影响和LSTM风速预测模型应用于实际低空风切变预报的方法等等， 这些问题需要进一步深入研究。