双流机场雷暴天气预报方法研究

李典南，徐 海，许东蓓

(1.贵州省气象台，贵州 贵阳 550002；2.中国民用航空西南地区空中交通管理局气象中心，四川 成都 610202；3.成都信息工程大学大气科学学院高原大气与环境四川省重点实验室，四川 成都 610225)

0 引 言

雷暴是一种常见的带有雷电现象的局地性强对流天气[1]，其水平发展的范围为几千米至几十千米，垂直发展的高度可达8～15 km，持续时间为几分钟到几小时不等[2]。雷暴的发生通常伴有雷暴大风、短时强降水、冰雹甚至龙卷等一种或多种天气现象[3]，给航空飞行带来了十分恶劣的影响。据统计，全球每天约有4.4万个雷暴发生，并且在任意一个时刻都有0.2～0.4万个雷暴在活动，其活动面积占比全球面积高达1%[4]，严重威胁航空飞行的安全。近年来国内外相关统计显示，每年由于气象因素造成的飞行事故超过事故总数的1/3[5]，而在这部分由于气象因素造成的飞行事故中，又有至少1/3是受到雷暴天气的影响[6-8]。由此可见，提高雷暴天气的预报水平，可在一定程度上减少航空飞行事故的发生以及降低飞行事故带来的损失。

雷暴的形成通常需要3个必要条件：深厚且明显不稳定的气层、充沛的水汽以及合适的触发条件[9]。在雷暴发生发展过程中，大气各气象要素在不断发生变化，当某些气象要素的强度达到一定程度时，雷暴便会发生。为了表征某个时刻大气的性质特征，可选用某些参数的数值来衡量，这就引入了对流参数的概念：将基于气块法分析对流条件而得出的参数，统称为对流参数。在实际预报中，可依据对流参数的数值特征来判断大气层结的不稳定性、水汽条件和动力条件等的强弱变化情况，并基于此得出未来强对流天气的发生时间、落区和强度等信息[10]。利用对流参数进行预报的方法有许多，比如指标叠套法[11]，即选取多个因子并计算其指标，通过对这些指标进行分析叠套来进行潜势预报。还可将多个因子相结合为一个新的预报因子，比如将垂直螺旋度与对流有效位能相结合而成的“垂直能量螺旋度指数(VEH)”[12-13]，VEH能有效体现强对流天气的动力因子和热力因子的共同贡献，杨晶轶等[14]通过实际验证，发现VEH对双流机场的雷暴天气的预报具有一定指示性。此外，建立潜势预报模型也是一种常用的预报方法，比如：主分量旋转法、二级逻辑回归法、Bayes判别法、Logistic回归判别法、神经网络法等[15-20]。其大体思路是：筛选出与雷暴发生关系最密切的几个对流参数或物理量，再将这些参数或物理量利用不同方法来建立预报模型。

成都双流国际机场位于103°57′02″E、30°34′47″N，地处四川盆地西部的平原腹地，海拔约为504.3 m，是中国大陆的八大航空枢纽之一，是西南地区重要的客货集散地。双流机场地势西高东低，自西北向东南倾斜[21]，受青藏高原大地形和盆地气候特征的影响，容易在近地面形成较强的垂直风切变以及气流辐合线、低涡等中小尺度天气系统[22]。在夏季，双流机场常处于太平洋副热带高压的西南边缘，加之机场东面河流纵横、水汽充沛[23-25]，使得雷暴天气频发。

为研究双流机场的雷暴天气并为其空中管制提供技术支持，本文参考了常用预报方法的思路来建立潜势预报模型，并在其基础上做了一定改进。全文综合利用了多种气象资料，首先运用统计学方法分析双流机场雷暴天气的时间变化特征，再利用相关性分析筛选出对双流机场雷暴天气指示性较好的物理量因子以及主要雷达特征量，最后在此基础上基于二级逻辑回归法建立潜势预报模型。

1 资料与方法

1.1 资料数据

所用的资料主要为：①2013—2018年双流机场的逐小时观测资料。该资料由双流机场例行每小时记录一次，要素包括风向、风速、气温、露点、相对湿度、修正海平面气压、云、能见度、天气现象等等。②欧洲中期天气预报中心(ECMWF)2013—2018年的ERA-interim逐6 h再分析资料(水平分辨率为0.125°×0.125°)，主要为高空资料。③成都市气象局多普勒天气雷达S波段SC型号的雷达产品资料，主要包含回波顶高、1.5°仰角基本反射率、3.4°仰角基本反射率、反射率最大回波出现的高度、组合反射率、垂直累积液态水含量和3 h累积降水量，所选取的数值是以双流机场为中心、半径20 km范围内的最大值。

1.2 研究方法

1.2.1 相关系数的计算 相关系数是最早由统计学家卡尔·皮尔逊设计的统计指标，是研究变量之间线性相关程度的量，一般用字母r表示。记某一天雷暴过程的发生与否为Y，发生为Y=1，不发生为Y=0；记某因子为X，X的样本序列为Xi(i=1,2,3,…)。为判断X因子与雷暴发生与否Y的相关性，可运用点双序列公式[26]来计算相关系数：

(1)

1.2.2 多元回归方程的建立 利用二级逻辑回归的思路建立多元回归方程组(第一级预报方程、第二级消空方程)[27]，该方法的基本思路是把一个预报事件与多个预报因子均看成随机事件，将随机事件Y的出现用“1”记录，未出现用“0”记录，将事件Y与n个预报因子之间的关系看成是n个因子已经出现的条件下事件Y出现与否的关系。

在逻辑回归中，X和Y均为“0”“1”化后的矩阵，具体方法为：X表示因子与阈值之间的关系，当因子的值大于其阈值，则记为X=1，反之记为X=0；Y表示当日有无雷暴过程，若当日有雷暴过程发生则记为Y=1，反之记为Y=0。通过拟合计算，可建立如下形式的多元回归方程(b0,b1,b2, ……,bn为回归参数)：

Y=b0+b1×X1+b2×X2+b3×X3…+bn×Xn

(2)

1.2.3 预报因子的检验方法 在确定好预报因子和预报模型之后，将一份样本数据带入检验，以评估其预报效果，具体步骤为：

首先，将预报因子X1、X2、X3……Xn的指标(阈值)带入对应的第一级预报方程，用k1代表计算出的Y值；再将样本数据带入同一方程，用Y(预报)代表利用样本数据计算出的Y值。若有Y(预报)≥k1，则预报为有雷暴过程，反之为无雷暴过程。

第二步，将预报因子Xn+1、Xn+2、Xn+3……X2n的指标(阈值)带入对应的第二级消空方程，用k2代表计算出的Y值，同样再将样本数据带入第二步的方程得出Y(消空)，预报有无雷暴的方法与第一步相同。

第三步，结合第一级预报方程、第二级消空方程的预报结论来综合判定：若某日Y(预报)预报为“无”，则预报结论为“无”。若某日Y(预报)预报为“有”，则进一步利用第二级消空方程进行判定——若Y(消空)也预报为“有”，则结论为“有”；若Y(消空)预报为“无”，则结论为“无”。

第四步，通过计算和统计，将预报结论带入王名才等[28]公式来分析预报效果：

(3)

(4)

(5)

其中，N为实况发生雷暴的总次数，n为预报有雷暴的总次数，n1为正确预报出雷暴的次数，n′为预报有雷暴但实况无雷暴的次数。

2 双流机场雷暴变化特征

根据双流机场的逐小时天气现象资料，将每一次雷暴天气自开始时间至结束时间的整个过程记为1次(除1.2节外均按此标准统计)。通过统计(表1)，2013—2018年双流机场共有221 d出现雷暴过程，累计发生283次。其中伴随降水过程的雷暴共发生了187 d，累计218次，占雷暴总次数的77.03%；其余为干雷暴过程，共发生了62 d，累计65次，占雷暴总次数的22.97%。下面对这6 a的雷暴天气进行月变化和日变化的特征讨论。

表1 双流机场2013—2018年雷暴发生天数和次数统计Tab.1 Statistics on the number of days and the number of thunderstorms at Shuangliu Airport from 2013 to 2018

2.1 月变化特征

2013—2018年的初雷在2月、3月各出现1次，4月出现4次；终雷在8月出现2次，9月出现3次，10月出现1次。图1为双流机场2013—2018年雷暴逐月发生次数统计图，由图可见雷暴的发生次数随月变化呈显著的单峰型：1—7月雷暴发生次数随时间增加，7月雷暴发生的次数最多，8月雷暴发生的次数仅次于7月。经计算，每年7、8月雷暴发生的次数之和，均达到各年全年发生总次数的一半以上。从9月开始雷暴发生次数明显减少，11月—次年1月基本无雷暴发生。可见春季(3—5月)和夏季(6—8月)是雷暴的多发季，其中以夏季7、8月发生次数最多。这是因为在春季和夏季，双流机场边界层附近的气流受日照加热抬升较强，并且机场东侧水域丰富，近地面潮湿，尤其在炎热的夏季，双流机场处太平洋副热带高压和青藏高压的过渡区中，常受西南暖湿气流和高原波动槽的影响，所以形成了高温高湿且不稳定的环境，有利于雷暴天气发生[29]。而秋季和冬季四川盆地层结较稳定，水汽条件和抬升条件较弱，雷暴发生条件不足，故发生次数较少。

图1 双流机场2013—2018年雷暴各年的逐月发生次数柱状图Fig.1 Histogram chart of the number of monthly thunderstorms occurrence in Shuangliu Airport from 2013 to 2018

2.2 日变化特征

将每小时作为一个时段，若某个时段有雷暴过程发生或有雷暴过程持续，则记该时段发生1次，统计结果如图2所示。由图可知，雷暴在上午发生的次数最少，6 a里平均每小时仅发生十几次。午后雷暴发生次数逐渐增加，21时—次日06时(北京时，下同)是高发时段，6 a共计发生521次，占比总次数的58.54%，这期间平均每小时发生雷暴超过40次。若以每3 h为1个时段，则雷暴发生次数最多的时间段是21时—00时，共计发生193次；其次是00—03时169次，03—06时159次，18—21时98次。综合来看，雷暴在一天之内的时间分布呈明显的“夜雷多、日雷少”特征。

图2 双流机场2013—2018年雷暴在一天各时段内发生次数折线图Fig.2 Line chart of the number of thunderstorms occurrence in Shuangliu Airport from 2013 to 2018 during various periods of the day

形成该特征的主要原因有：第一，双流机场所处地形较闭塞，近地面潮湿，云雾较多，夜间云层上部辐射冷却速度较快，云层下部由于云雾的保暖作用而冷却较慢，气层上冷下暖趋于不稳定，易形成雷暴天气[30]。第二，夏季太平洋副热带高压伸入我国西南地区且活动频繁，太平洋副热带高压白天多西伸，使双流机场受西南暖湿气流影响，为雷暴的发生发展预备了良好条件；夜间多东退，高原西风槽东移，多在后半夜移动到双流机场附近引起雷暴天气[31]。此外，弱冷锋在经过四川盆地时，由于受到地形抬升和太阳辐射的共同作用，在白天常表现为锋消；而夜间随着地表辐射增强，又逐渐表现为锋生，易形成锋面雷暴[32]。

3 预报因子的筛选

3.1 物理量因子

基于雷暴形成的3个必要条件[33]，初步选取了850 hPa假相当位温、850与500 hPa假相当位温差、相对湿度、比湿等物理量因子以及对流有效位能、K指数、强天气威胁指数、热力总指数、抬升指数、风切变指数等对流参数(以下统称为物理量因子)。借助欧洲中心ERA-interim逐6 h再分析资料，针对双流机场2013—2018年的雷暴天气进行物理量因子的相关系数计算(公式(1))。计算样本为2013—2018年每年的4—8月(2014年8月资料不全，故排除，下同)，共计887 d，其中有208 d发生雷暴天气。物理量因子Xi的取值为每日的算术平均值。

如表2所示，将通过α=0.05显著性水平检验的因子按相关系数的绝对值从大到小排列，相关系数越大说明因子对雷暴的预报指示性越好。由表2可知，对流有效位能、K指数、850与500 hPa假相当位温差、850 hPa比湿、850 hPa假相当位温、强天气威胁指数对雷暴有较高的指示意义，相关系数绝对值均高于0.4。除此之外，700 hPa比湿、热力总指数、抬升指数与雷暴的相关性也较好，相关系数绝对值均高于0.3。

表2 物理量因子与雷暴发生情况的相关系数Tab.2 The correlation coefficient between the physical factors and thunderstorms

统计学中，箱线图(又称箱图)能以简单的组合图形直观地反映较大容量样本的值分布情况和数据批中隐含的结构信息，可通过分析数据批的分布形状、排除异常数据点，来进行数据批之间的比较[34]。利用相关系数绝对值大于0.3的前9个物理量因子的日均值，绘制雷暴日、非雷暴日的箱线图(图3)。由图可知，每组因子箱线图的箱体之间存在一定交集(值域交集)，计算各因子箱体交叉部分占雷暴日箱体的比例，得出的结论按从小到达排列为：对流有效位能0%；K指数0%；抬升指数6.56%；850与500 hPa的假相当位温差10.49%；850 hPa比湿15.27%；850 hPa假相当位温16.18%；强天气威胁指数20.21%；热力总指数29.07%；700 hPa比湿40.94%。当某因子在雷暴日和非雷暴日的值域交集越少，说明该因子对雷暴的预报指示性越好，即对流有效位能、K指数、850与500 hPa的假相当位温差和抬升指数对雷暴的预报更有利。

图3 2013—2018年887个样本有无雷暴情况下各物理量因子箱线图(上下*号：第99百分位数和第1百分位数；上下短实线：第90百分位数和第10百分位数；矩形内实线：中位数；矩形上下边线：第75百分位数和第25百分位数)Fig.3 Box plot of each physical factors in 887 samples with or without thunderstorms from 2013 to 2018(Upper and lower *: the 99th percentile and the 1st percentile; upper and lower short solid lines: the 90th percentile and the 10th percentile; solid line within the rectangle: median; upper and lower edges of the rectangle: 75th percentile and 25th percentile)

在样本容量固定的情况下，因子过多不仅对拟合方程起的贡献不大，因子本身存在的随机因素还容易对方程的稳定性产生影响[19]。基于雷暴发生的3个基本条件以及相关系数和箱线图所显示结果综合考虑，最后确定出如下4个物理量预报因子：对流有效位能、K指数、850 hPa比湿、850与500 hPa的假相当位温差。

3.2 雷达因子

多普勒天气雷达是监测实况天气的重要手段，广泛运用于强对流天气的监测和预警。许多雷达产品都可用于反映强对流天气的发展状态，比如回波顶高、垂直累积液态水含量常被用来评估风暴强度[35]，雷达回波强度、回波顶高等常用来反映闪电和雷暴云发展的高度[36-37]等等。

结合实际业务经验[37-42]，初步选取了与强对流天气联系较紧密的7个雷达因子进行相关系数的计算(表3)。计算样本为2013—2018年每年4—8月的雷暴日，共计208 d，雷达因子值取逐小时整点时刻的值，共计4 992 h，其中发生雷暴的时间有873 h。表3中雷达因子已按相关系数的绝对值从大到小排列，且均通过α=0.05的显著性水平检验。

表3 多普勒天气雷达因子与雷暴发生情况的相关系数Tab.3 The correlation coefficient between the doppler weather radar factors and thunderstorms

以相同方法绘制箱线图(图4)，并计算各雷达因子箱体交叉部分(值域交集)占雷暴日箱体的比例，结论如下：1.5°仰角基本反射率0%；垂直累积液态水含量0%；3 h累积降水量13.52%；3.4°仰角基本反射率15.79%；回波顶高20.00%；组合反射率52.17%；反射率最大回波出现的高度46.03%。综合相关系数和箱线图两方面考虑，最终确定的4个雷达预报因子为：回波顶高、1.5°仰角基本反射率、3.4°仰角基本反射率、垂直累积液态水含量。

图4 2013—2018年4992个样本有无雷暴过程的情况下各雷达因子箱线图(上下*号：第99百分位数和第1百分位数；上下短实线：第90百分位数和第10百分位数；矩形内实线：中位数；矩形上下边线：第75百分位数和第25百分位数)Fig.4 Box plot of each radar factors in 887 samples with or without thunderstorms from 2013 to 2018(Upper and lower *: the 99th percentile and the 1st percentile; upper and lower short solid lines: the 90th percentile and the 10th percentile; solid line within the rectangle: median; upper and lower edges of the rectangle: 75th percentile and 25th percentile)

4 潜势预报模型的建立

4.1 预报因子阈值的计算

利用雷暴日物理量预报因子月平均值以及雷暴发生的整点时刻的雷达预报因子月平均值，绘制如图5所示折线图。图5a显示，物理量预报因子的月平均值随时间变化的趋势基本一致：最小值出现在4、5月，最大值出现在7、8月，具有明显的月变化特征或季节变化特征，说明物理量预报因子阈值的计算需要分春季、夏季考虑。而图5b显示，雷达预报因子的月平均值无明显的时间特征，说明满足雷暴发生条件的雷达预报因子受季节影响较小，故阈值的计算无需分季节考虑。

图5 2013—2018年双流机场雷暴日物理量预报因子(a)、雷达预报因子(b)的月平均值Fig.5 The monthly average of thunderstorm daily physical predictors (a) and radar predictors (b) of Shuangliu Airport from 2013 to 2018

将物理量预报因子按季节分为春夏两组(4、5月为春季，6—8月为夏季)，雷达预报因子不分组；再将每组每个因子的值按从大到小排列，选取每个因子的第一、四分位值作为该因子的预报阈值[17]，即双流机场雷暴天气的预报指标，计算结果如表4所示。

表4 双流机场雷暴天气的预报因子阈值Tab.4 Thresholds of predictors for thunderstorms in Shuangliu Airport

4.2 潜势预报方程的建立

参照1.2.2节，用X1、X2、X3、X4分别表示850 hPa比湿、K指数、850与500 hPa的假相当位温差、对流有效位能，利用这些物理量预报因子建立如公式(2)所示的春季Y1、夏季Y2多元回归预报方程，作为第一级预报方程：

Y1=0.012+0.100×X1+0.081×X2+0.180×X3+0.038×X4

(6)

Y2=0.052+0.142×X1+0.243×X2+0.200×X3+0.031×X4

(7)

若仅使用Y1、Y2方程进行预报，容易造成命中率较高、虚警率也较高的情况(已验证，结果略)，故进一步利用雷达预报因子建立Y3方程作为第二级消空方程：

Y3=0.050+0.064×X5+0.163×X6+0.063×X7+0.093×X8

(8)

其中，X5、X6、X7、X8分别表示回波顶高、1.5°仰角基本反射率、3.4°仰角基本反射率、垂直累积液态水含量。最后再对求得的上述回归方程Y1、Y2、Y3进行F检验[43]，结果表明这3个方程是显著的。

4.3 预报效果的检验

用于检验的样本为2018年4—8月的相关数据，利用1.2.3节的检验方法：首先，将表3中所示的物理量预报因子的阈值按季节带入对应的第一级预报方程(公式(6)或公式(7))，再将雷达预报因子的阈值带入第二级消空方程(公式(8))，用k1、k2、k3分别代表计算出的Y1、Y2、Y3的值，即各方程的阈值：k1=Y1=19.00、k2=Y2=85.32、k3=Y3=7.85。再将2018年4—8月每日的物理量预报因子日均值依季节带入对应的第一级预报方程Y1、Y2，在春季，当Y1(预报)≥k1时，预报为有雷暴过程，反之无雷暴过程；同理，在夏季，当Y2(预报)≥k2时，预报为有雷暴过程，反之无雷暴过程。第三步，运用第二级消空方程，将样本中每日逐小时的雷达预报因子值带入Y3，当Y3(消空)≥k3时，预报为有雷暴过程，反之为无雷暴过程。第四步，结合两级方程的结论来综合判定：若Y1(或Y2)、Y3均预报为“有”，则最终结论为“有”，其他情况下最终结论均为“无”。最后，将最终预报结果带入公式(3)、公式(4)和公式(5)进行计算，结论如表5所示。

表5 利用Y1、Y2、Y3方程预报2018年4—8月雷暴天气Tab.5 The conclusion of using Y1, Y2, Y3 equations to forecast thunderstorm from April to August 2018

由表5可知，预报指标在春季的命中率、虚警率、临界成功指数分别为100%、50.00%、50.00%，命中率虽然很高，但由于空报数较多，使得虚警率较高、临界成功指数较低。指标在夏季的命中率、虚警率、临界成功指数分别为80.49%、8.33%、75.00%，命中率虽不如春季，但因为空报数较少，所以虚警率较低、临界成功指数较高。因此可以认为这套指标以及建立的两级预报方程，对双流机场雷暴天气的预报具有一定指示意义，且综合来看在夏季的预报效果更好。在实际运用中，可将第一级预报方程用作24 h潜势预报，若预报结论显示第2 d有雷暴天气，则加强天气监测，再配合第二级消空方程，以更准确地判断临近时间的天气情况。

5 结论与展望

5.1 结论

本文利用双流国际机场2013—2018年的逐小时气象观测资料、欧洲中心ERA-interim逐6 h再分析资料、多普勒天气雷达产品资料，对双流机场的雷暴天气进行了时间特征分析，并选取预报因子建立潜势预报模型，结果表明：

①双流机场2013—2018年共有221 d发生雷暴天气，累积发生283次，77.03%的雷暴过程伴随降水。雷暴发生次数随时间呈单峰型，春夏为多发季，其中以夏季7、8月发生的次数最多。雷暴多于午后发生，21时—次日06时是高发时段，呈现明显的“夜雷多、日雷少”特征。

②依据雷暴发生的基本条件和实际业务经验选取物理量和雷达参数，通过相关分析，筛选出预报因子：对流有效位能、K指数、850 hPa比湿、850与500 hPa假相当位温差、回波顶高、1.5°仰角基本反射率、3.4°仰角基本反射率、垂直累积液态水含量。再基于二级逻辑回归的思路建立预报方程和消空方程：利用物理量预报因子分春、夏两季建立两个第一级预报方程；利用雷达预报因子建立一个第二级消空方程。通过回代检验，结果表明，夏季的命中率低于春季，但因为夏季预报的空报数较少，所以虚警率和临界成功指数相对于春季都表现较好。由此可认为所建立的潜势预报模型对双流机场雷暴天气的预报具有一定指示性，且综合来看在夏季的预报效果更好。

5.2 反思与展望

论文中仍有许多的内容需要完善和探讨，首先是建立预报方程方面，使用的仅仅是二级逻辑回归法，若选择多种方法(比如神经元方法、SVM方法等)，就可通过对比检验得出最优方案，或形成集成预报。其次是方程的检验方面，使用的是2018年的数据进行回代，若资料充足，选取2019年或之后的资料进行检验，会更具有说服力。最后，在雷达资料选取上，由于双流机场的雷达资料缺失较多，故本文使用的雷达资料是成都市气象局的多普勒天气雷达产品资料，但若采用双流机场的雷达资料进行分析，得出的结论会更加有说服力。