1 引言

雷暴在华南夏季经常发生, 常伴随着强降水、风切变、冰雹、闪电、低能见度等天气现象, 对航空飞行正常运行易造成很大的影响。因此, 雷暴潜势预报对航空业务来说具有重要意义。很多学者建立对流指数来预测雷暴潜势^[1-7]。雷暴的产生机制较复杂, 具有非线性特征。神经网络算法具有非线性映射能力和并行性、适应性、容错性及学习能力, 被许多学者应用于雷暴潜势预报。Agostino^[8]采用神经网络对探空资料和闪电定位资料进行训练, 来预报雷暴的发生以及闪电密度。陈勇伟等^[9]基于BP神经网络采用南京地区2008年的闪电定位资料和探空资料预报了南京地区2009年6-8月的雷暴潜势。杨仲江等^[10]利用闪电定位资料和探空资料, 采用双隐层BP神经网络对太原地区的雷暴潜势进行预报, 发现双隐层BP神经网络比BP神经网络、多元统计回归法预报效果更好。周明薇等^[11]根据2008-2010年夏季邵阳地区的NCEP全球再分析资料和闪电定位资料, 采用支持向量机SVM分类算法建立了6小时雷暴潜势预报模型。但是BP神经网络算法为一种局部搜索的优化方法, 采用了梯度下降法, 学习速度慢, 网络训练失败概率大。Fernández-Delgado等^[12]应用了179个分类器对UCI数据库的相应数据集进行分类预测, 发现随机森林分类器的各项性能表现排在前5名, 优于BP神经网络。

随机森林算法由Breiman Leo^[13]和Adele Cutler等^[14]提出, 该算法结合了Breiman的"Bootstrap aggregating"(自举汇聚法)思想和Ho^[15]的"Random Subspace"(随机子空间)方法。与BP神经网络算法相比, 随机森林算法具有以下优点: (1) 高准确度的分类器; (2) 可处理大量的输入变量; (3) 可评估变量的重要性; (4) 可在内部对于一般化后的误差产生不偏差的估计; (5) 可有效地估计缺失的数据并保持精度; (6) 对于不平衡的分类资料集来说, 可平衡误差; (7) 可计算个例中的亲近度, 对于数据挖掘、侦测离群点和将资料可视化非常有用; (8) 可实现无监督聚类、数据视图和离群点检测; (9) 可检测可变相互作用。因此, 随机森林算法在遥感^[16]、环境^[17]、商业^[18]、医学^[19]等方面得到广泛应用。

广州白云机场每天要制作两次未来24小时重要天气概率预报, 其中雷暴是发生频率最高的重要天气, 雷暴潜势的预报十分重要。平时预报员分析预报雷暴的资料基本来源于天气图、数值模式, 空间分辨率较大, 都超过100 km。而机场预报要求精细化, 100 km(广州白云机场终端区范围) 内的雷暴潜势预报基本靠天气系统分析和经验判断, 因此期望能有可靠的计算方法得到广州白云机场终端区范围内的雷暴潜势预报。本文使用2014-2015年夏季雷达MAX图(最大反射率图) 资料来获得雷暴发生有无数据, 结合探空数据建立广州白云机场终端区内3类区域的12小时随机森林分类模型, 对不同区域的雷暴潜势进行预报对比, 以期获得更高的预报准确率。

2 数据来源及处理 2.1 探空数据选取

选用2014-2015年夏季(6、7、8月)的探空数据和雷达MAX图资料。其中探空数据来自美国怀俄明州立大学提供的清远探空站和香港探空站, 一日2次, 分别是世界时00:00 UTC和12:00 UTC, 共选取20个探空资料对流指数因子(表 1)。表 2为表 1公式中出现的参数的释意。

2.2 雷暴样本选取

雷暴样本在白云机场终端区内选取，以塔台为中心，分别选取离塔台8 km(区域1)、50 km(区域2)、100km(区域3)的3个区域。雷暴发生有无根据雷达MAX图(图 1)来判定，判别条件如下。

(1) 雷达反射率因子大于等于35 dBZ。研究表明35 dBZ与成熟积雨云的降雨强度和闪电有很好的相关关系，而且常被选为识别和追踪雷暴算法所选的阈值^[20-23]。徐慧^[24]、郭巍等^[25]均利用雷达反射率因子大于等于35 dBZ作为标准来研究中国区域的对流初生，证明了这一标准在中国的适用性。同时在实际飞行中，机组多次反映当机载雷达探测到雷达反射率因子大于等于35 dBZ及其以上的云团就实施绕飞。

(2) 雷达反射率因子大于等于35 dBZ区域在10 km以上的连续区域。对应的MAX图中大于等于35 dBZ的像素点要超过23个。

(3) 雷达最大反射率对应高度大于等于6 km，这是为了与低云区别开来。

出现满足以上3个条件的情况，则判定为有雷暴发生。没有雷暴发生的样本记为1，有雷暴发生的样本则记为2。如图 1所示，2014年6月23日05:40 UTC，区域1、2、3范围内均有雷暴云团活动。

2.3 样本处理

将清远探空站和香港探空站的探空数据进行距离线性插值法得到白云机场的探空数据。距离线性插值法如下：

$ \begin{array}{*{20}{l}} \;\;\;\;a=\sqrt{(\mathrm{LAT} 1-\mathrm{LAT})^{2}+(\mathrm{LNG} 1-\mathrm{LNG})^{2}} \times \pi \times \\6370 / 180 \\ \;\;\;\;b=\sqrt{\left(\mathrm{LAT} 2-\mathrm{LAT}^{2}+(\mathrm{LNG} 2-\mathrm{LNG})^{2}\right.} \times \pi \times \\6370 / 180 \\ \;\;\;\;Z=(a \times y+b \times x) /(a+b) \end{array} $

(1)

其中，LAT和LNG(本场经纬度)、LAT1和LNG1 (清远探空站经纬度)、LAT2和LNG2(香港探空站经纬度)，a为清远探空站到白云机场的距离，b为香港探空站到白云机场的距离，x为清远探空站的探空数据，为香港探空站对应清远探空站的同一天同一时次的探空数据，Z为插值后的白云机场探空数据。

统计区域1、区域2范围内、区域3范围内雷暴发生与否数据共得三组数据。将同一天的00:00 UTC的探空数据与00:00—11:59 UTC雷暴发生与否数据匹配，12:00 UTC的探空数据与12:00— 23:59 UTC雷暴发生与否数据匹配，因此每组样本包含20个探空指数因子和1个雷暴发生与否数据。除去缺失资料的样本，2014年的数据为183组，2015年的数据为182组，共计365组。其中区域1样本组雷暴概率为135/365≈37°%; 区域2样本雷暴概率为237/365≈65°%; 区域3样本雷暴概率为277/365=76%。

2.4 预报因子的选取

在筛选合适的输人因子时，由于对流指数是连续型因子，雷暴发生与否为1、非连续型变量，所以不能采用皮尔森相关系数一类的线性相关系数来衡量两者的关系。灰色关联度分析是研究因子间的几何对应关系(即灰色关联度)，序列曲线的几何形状的相似程度与灰色关联度成正比^[26-27]，它计算量小，适用于任何类型包括无规律的数据，在多因素，非线性领域得到广泛应用^[28-31]。

灰色关联分析的步骤如下。

(1) 确定雷暴发生与否数列和对流指数因子数列。设雷暴发生与否数列为$X_{0}^{\prime}=\left\{X_{0}^{\prime}(k) \mid k=\right.$1，2，……，n}；对流指数因子数列$Y_{i}^{\prime}=\left\{Y_{i}^{\prime}(k) \mid k=\right.$ 1，2，……，n}，i = 1，2，……，m。这里m为指数因子个数20，n为样本组数365。

(2) 对指标数据进行无量纲化。由于各因子的物理意义不同，导致量纲不同，为了便于比较，采用均值化法(式(2))进行数据的无量纲处理。

$ X_{j}(k)=\frac{Y_{j}^{\prime}(k)}{\frac{1}{n} \sum_{k=1}^{n} Y_{j}^{\prime}(k)} $

(2)

其中，j=0，1，……，20，k=1，2，……，365。

(3) 计算每个对流指数因子数列与雷暴发生与否数列对应元素的绝对差值。即$\text { | } X_{0}(k)-\left.X_{i}(k)\right|$。

(4) 确定$\min\limits_{1 \leq i \leq \mathrm{m}} \min\limits_{1 \leq k \leq n}\left|X_{0}(k)-X_{i}(k)\right|$与$\max\limits_{1 \leq i \leq m} \max\limits_{1 \leq k \leq n}\left|X_{0}(k)-X_{i}(k)\right|$。

(5) 计算每个对流指数因子数列与雷暴发生与否数列对应元素的灰色关联系数ξ_i(k):

$ \xi_{i}(k)=\frac{\min\limits_{i} \min\limits_{k}\left|X_{0}(k)-X_{i}(k)\right|+\mu \times \max\limits_{i} \max\limits_{k}\left|X_{0}(k)-X_{i}(k)\right|}{\left|X_{0}(k)-X_{i}(k)\right|+\mu \times \max\limits_{i} \max\limits_{k}\left|X_{0}(k)-X_{i}(k)\right|} $

(3)

其中为分辨系数，此处取μ=0.5。

(6) 计算灰色关联度r₀。分别计算各个指数与雷暴发生与否序列对应元素的均值即灰色关联度r_0i：

$ r_{0 i}=\frac{1}{n} \sum_{k=1}^{n} \xi i(k) $

(4)

分别将20个对流指数因子与预报量(雷暴发生与否)作灰色关联度分析, 经过计算得到表 3, 表 3为白云机场终端区区域1、区域2、区域3范围内, 20个对流指数与雷暴发生有无的灰色关联度r0。区域1范围内的r₀平均值为0.968 9, 大于平均值的指数由大到小依次为PLCL、SWI、BRCV、CT、KI、LFCV、TT、LFCT、BRN、VT、LCVT、CAPV、CAPE、LI。由于LI和LFCV、CAPE和CAPV物理意义近似, 所以保留两者中的大值LFCV、CAPV, 删去LI、CAPE。最终区域1选取了10个预报因子为PLCL、SWI、BRCV、CT、KI、LFCV、TT、LFCT、BRN、VT。同理, 区域2的r₀平均值为0.971 4, 共选取了10个预报因子为PLCL、KI、CT、BRCV、SWI、LFCV、TT、VT、LCVT、CAPV。区域3的r₀平均值为0.973 5, 共选取了10个预报因子为PLCL、KI、BRCV、CT、TT、LFCV、SWI、VT、LCVT、CAPV。在白云机场终端区三个区域内, PLCL、KI、CT、BRCV、SWI、LFCV、TT、VT、LCVT、CAPV这10个对流参数与雷暴发生有无的灰色关联度较大, 其中PLCL、KI、CT、BRCV这4个指数的关联度一直排在前五里, 表明雷暴产生跟抬升凝结高度气压、850~700 hPa的温度和露点温度、风切变的关系最紧密。

3 随机森林算法

随机森林算法是一种包含多个无关联的决策树的分类器, 并取所有决策树中分类结果最多的那类为最终结果, 其算法原理如下。

3.1 Bootstrap法重采样

设集合Z为Z={Z(k)|k=1, 2, ……, n}, 若每次有放回地从集合S中抽取一个样本, 一共抽取n次, 形成新的集合Z*, 则集合S*中不包含某个样本Z(i)(i=1, 2, ……, n)的概率为:

$ p=\left(1-\frac{1}{n}\right)^{n} $

(5)

当n→∞时, 有,

$ \lim\limits_{n \rightarrow \infty} p=\lim\limits_{n \rightarrow \infty}\left(1-\frac{1}{n}\right)^{n}=\mathrm{e}^{-1} \approx 0.368 $

(6)

因此, 通过Bootstrap法重采样, 除去重复的样本, 新集合Z* 约占原集合Z中(1-0.368) × 100%= 63.2%的样本。

3.2 随机森林算法流程

(1) 基于Bootstrap方法随机产生R个训练集Z₁、Z₂、……、Z_R。

(2) 每个训练集生成对应的决策树C₁、C₂、……、C_R; 设样本的属性个数为A, 从A个属性中随机抽取aa(aa为大于0且小于A的整数)个属性用最好的分裂方式对当前节点进行分裂, 得到分裂属性集。

(3) 所有决策树不剪枝。

(4) 利用每个决策树对测试集样本X进行测试, 得到对应的测试结果C₁(X)、C₂(X)、……、C_R (X)。

(5) 将R个决策树中输出最多的分类类别作为测试集样本X所属的分类类别。

3.3 输入输出数据处理

为了对比分析, 设立3个样本输入组: 区域1样本输入组、区域2范围内样本输入组、区域3范围内样本输入组。其中区域1样本输入组包括365组10个预报因子, 区域2范围内样本输入组包括365组10个预报因子, 区域3范围内样本输入组包括365组10个预报因子。为了简化计算, 缩小量值, 加快网络收敛, 将样本输入组进行[0.1, 0.9] 之间的归一化:

$ \mathrm{Q}=0.1+\frac{P-P_{\min }}{P_{\max }-P_{\min }} \times 0.8 $

(7)

其中, P为归一化前的样本输入数据, P_min为P的最小值, P_max为P的最大值, Q为归一化之后的矩阵。

样本输出组为365组雷暴发生有无数据。365组输入输出样本中, 将183组用来训练网络, 182组用于检验。

4 结果分析 4.1 区域1预报结果及泛化分析

本文采用randomforest-matlab开源工具箱进行计算。区域1共有10个预报因子, 因此样本的属性个数A为10, 分裂属性个数aa为10的取整数值3。为了减少随机性的影响, 决策树颗数T分别取{50, 100, 150, 200, ……, 950, 1 000}共20种情况, 每种情况运行100次模型, 取其准确值(即测试结果里预报值1, 实况值1和预报值2, 实况值2的总和除以182)的最高值作为当前决策树颗数下的分类正确率。图 2为区域1样本输入随机森林模型中, 决策树颗数对性能的影响, 可看到当决策树颗数为400时, 分类正确率最高。

设立T=400, aa=3, 运行随机森林模型, 将结果用两变量预报验证列联表^[32]表示得到表 4。采用4个指数(如式组(8)所示)来评估预报结果, 分别为: 临界成功指数CSI^[33]、预报准确率AF^[34]、虚假报警率FAR^[32]和探测概率POD^[32]。

$ \left\{\begin{aligned} \mathrm{CSI} &=\frac{e}{e+f+g}=35.23 \% \\ \mathrm{AF} &=\frac{e+h}{e+f+g+h}=68.68 \% \\ \mathrm{FAR} &=\frac{f}{e+f}=45.61 \% \\ \mathrm{POD} &=\frac{e}{e+g}=50 \% \end{aligned}\right. $

(8)

随机森林有一个重要的特点: 没有必要对它进行交叉验证或用一个独立的测试集来获得误差的一个无偏估计。它可在内部进行评估, 即在生成过程中就对误差建立了一个无偏估计。袋外错误率oob是随机森林泛化误差的一个无偏估计, 它的结果近似于需要大量计算的k折交叉验证。如式(9)所示, 区域1的袋外错误率oob为31.32%, 表明超过2/3的样本能正确分类, 泛化性能较好。

$ \mathrm{oob}=\frac{f+g}{e+f+g+h}=31.32 \% $

(9)

4.2 区域2分析结果

区域2共有10个预报因子, 因此样本的属性个数A为10, 分裂属性个数aa为$\sqrt{10}$的取整数值3。图 3为区域2范围内样本输入随机森林模型中, 决策树颗数对性能的影响, 可看到当决策树颗数为250时, 分类正确率最高。

设立T=250, aa=3, 运行随机森林模型得到的结果为表 5和式组(10)。对比区域1预报结果, CSI、AF、POD值都要高, 而虚假报警率FAR值偏低, 表明雷暴预报效果更好了。当雷暴未来24小时发生可能性大于70% 时, 航空重要天气MDRS通报业务将启动。预报准确率AF和探测概率POD超过70%, 开始对业务有指示作用。区域2范围内的袋外错误率为28.57%, 泛化性能有提升。

$ \left\{\begin{aligned} \mathrm{CSI} &=\frac{e}{e+f+g}=63.64 \% \\ \mathrm{AF} &=\frac{e+h}{e+f+g+h}=71.43 \% \\ \mathrm{FAR} &=\frac{f}{e+f}=27.78 \% \\ \mathrm{POD} &=\frac{e}{e+g}=84.26 \% \\ \text { oob } &=\frac{f+g}{e+f+g+h}=28.57 \% \end{aligned}\right. $

(10)

4.3 区域3范围内分析结果

区域3共有10个预报因子, 因此样本的属性个数A为10, 分裂属性个数aa为$\sqrt{10}$的取整数值3。图 4为区域3样本输入随机森林模型中, 决策树颗数对性能的影响, 可看到当决策树颗数为100时, 分类正确率最高。

设立T=100, aa为$\sqrt{10}$的取整数值3, 运行随机森林模型得到的结果为表 6和式组(11)。三组预报结果里, 区域3的预报结果CSI、AF、POD值最高, 而虚假报警率FAR值最低, 表明雷暴预报效果最好。区域3的袋外错误率为24.73%, 泛化性能不错。

$ \left\{\begin{aligned} \mathrm{CSI} &=\frac{e}{e+f+g}=71.34 \% \\ \mathrm{AF} &=\frac{e+h}{e+f+g+h}=75.27 \% \\ \mathrm{FAR} &=\frac{f}{e+f}=22.22 \% \\ \mathrm{POD} &=\frac{e}{e+g}=89.6 \% \\ \mathrm{oob} &=\frac{f+g}{e+f+g+h}=24.73 \% \end{aligned}\right. $

(11)

5 结论

本文将2014-2015年夏季雷达MAX图资料获取雷暴发生有无数据, 结合美国怀俄明州立大学探空数据建立广州白云机场终端区3个区域的12小时随机森林分类模型, 通过对比可得到以下结论。

(1) 白云机场终端区区域1、区域2、区域3内, PLCL、KI、CT、BRCV、SWI、LFCV、TT、VT、LCVT、CAPV这10个对流参数与雷暴发生有无的灰色关联度较大。雷暴的产生跟抬升凝结高度气压、850~700 hPa的温度和露点温度、风切变的关系最紧密。

(2) 终端区区域面积越大, 雷暴发生样本比例越高, 临界成功指数CSI、预报准确率AF、探测概率POD越来越高, 虚假报警率FAR越来越低, 表明预报出来的准确率越来越高。区域2、区域3的预报准确率AF和探测概率POD超过70%, 给航空天气预报员分析预报雷暴提供了参考, 对航空重要天气MDRS通报业务有指示作用。

(3) 终端区三个区域的测试样本超过2/3能得到正确分类, 表明随机森林算法泛化性能较好。区域范围越大, 袋外错误率越低, 泛化性能越好。