楊潔，曹正，杜宇，武凱軍，李雅

(1.中國民用航空中南空中交通管理局氣象中心，廣東 廣州510405；2.中山大學大氣科學學院，廣東 珠海519082；3.中國民航大學航空氣象系，天津300300)

1 引 言

雷暴在華南夏季經常發(fā)生，常伴隨著強降水、風切變、冰雹、閃電、低能見度等天氣現(xiàn)象，對航空飛行正常運行易造成很大的影響。因此，雷暴潛勢預報對航空業(yè)務來說具有重要意義。很多學者建立對流指數(shù)來預測雷暴潛勢[1-7]。雷暴的產生機制較復雜，具有非線性特征。神經網絡算法具有非線性映射能力和并行性、適應性、容錯性及學習能力，被許多學者應用于雷暴潛勢預報。Agostino[8]采用神經網絡對探空資料和閃電定位資料進行訓練，來預報雷暴的發(fā)生以及閃電密度。陳勇偉等[9]基于BP神經網絡采用南京地區(qū)2008年的閃電定位資料和探空資料預報了南京地區(qū)2009年6—8月的雷暴潛勢。楊仲江等[10]利用閃電定位資料和探空資料，采用雙隱層BP神經網絡對太原地區(qū)的雷暴潛勢進行預報，發(fā)現(xiàn)雙隱層BP神經網絡比BP神經網絡、多元統(tǒng)計回歸法預報效果更好。周明薇等[11]根據(jù)2008—2010年夏季邵陽地區(qū)的NCEP全球再分析資料和閃電定位資料，采用支持向量機SVM分類算法建立了6小時雷暴潛勢預報模型。但是BP神經網絡算法為一種局部搜索的優(yōu)化方法，采用了梯度下降法，學習速度慢，網絡訓練失敗概率大。Fernández-Delgado等[12]應用了179個分類器對UCI數(shù)據(jù)庫的相應數(shù)據(jù)集進行分類預測，發(fā)現(xiàn)隨機森林分類器的各項性能表現(xiàn)排在前5名，優(yōu)于BP神經網絡。

隨機森林算法由Breiman Leo[13]和Adele Cutler等[14]提出，該算法結合了Breiman的“Bootstrap aggregating”(自舉匯聚法)思想和Ho[15]的“Random Subspace”(隨機子空間)方法。與BP神經網絡算法相比，隨機森林算法具有以下優(yōu)點：(1)高準確度的分類器；(2)可處理大量的輸入變量；(3)可評估變量的重要性；(4)可在內部對于一般化后的誤差產生不偏差的估計；(5)可有效地估計缺失的數(shù)據(jù)并保持精度；(6)對于不平衡的分類資料集來說，可平衡誤差；(7)可計算個例中的親近度，對于數(shù)據(jù)挖掘、偵測離群點和將資料可視化非常有用；(8)可實現(xiàn)無監(jiān)督聚類、數(shù)據(jù)視圖和離群點檢測；(9)可檢測可變相互作用。因此，隨機森林算法在遙感[16]、環(huán)境[17]、商業(yè)[18]、醫(yī)學[19]等方面得到廣泛應用。

廣州白云機場每天要制作兩次未來24小時重要天氣概率預報，其中雷暴是發(fā)生頻率最高的重要天氣，雷暴潛勢的預報十分重要。平時預報員分析預報雷暴的資料基本來源于天氣圖、數(shù)值模式，空間分辨率較大，都超過100 km。而機場預報要求精細化，100 km(廣州白云機場終端區(qū)范圍)內的雷暴潛勢預報基本靠天氣系統(tǒng)分析和經驗判斷，因此期望能有可靠的計算方法得到廣州白云機場終端區(qū)范圍內的雷暴潛勢預報。本文使用2014—2015年夏季雷達MAX圖(最大反射率圖)資料來獲得雷暴發(fā)生有無數(shù)據(jù)，結合探空數(shù)據(jù)建立廣州白云機場終端區(qū)內3類區(qū)域的12小時隨機森林分類模型，對不同區(qū)域的雷暴潛勢進行預報對比，以期獲得更高的預報準確率。

2 數(shù)據(jù)來源及處理

2.1 探空數(shù)據(jù)選取

選用2014—2015年夏季(6、7、8月)的探空數(shù)據(jù)和雷達MAX圖資料。其中探空數(shù)據(jù)來自美國懷俄明州立大學提供的清遠探空站和香港探空站，一日2次，分別是世界時00:00 UTC和12:00 UTC，共選取20個探空資料對流指數(shù)因子(表1)。表2為表1公式中出現(xiàn)的參數(shù)的釋意。

表1 懷俄明大學探空資料對流指數(shù)及其物理意義

表2 表1中公式中出現(xiàn)的參數(shù)的釋意

2.2 雷暴樣本選取

雷暴樣本在白云機場終端區(qū)內選取，以塔臺為中心，分別選取離塔臺8 km(區(qū)域1)、50 km(區(qū)域2)、100 km(區(qū)域3)的3個區(qū)域。雷暴發(fā)生有無根據(jù)雷達MAX圖(圖1)來判定，判別條件如下。

(1)雷達反射率因子大于等于35 dBZ。研究表明35 dBZ與成熟積雨云的降雨強度和閃電有很好的相關關系，而且常被選為識別和追蹤雷暴算法所選的閾值[20-23]。徐慧[24]、郭巍等[25]均利用雷達反射率因子大于等于35 dBZ作為標準來研究中國區(qū)域的對流初生，證明了這一標準在中國的適用性。同時在實際飛行中，機組多次反映當機載雷達探測到雷達反射率因子大于等于35 dBZ及其以上的云團就實施繞飛。

(2)雷達反射率因子大于等于35 dBZ區(qū)域在10 km以上的連續(xù)區(qū)域。對應的MAX圖中大于等于35 dBZ的像素點要超過23個。

(3)雷達最大反射率對應高度大于等于6 km，這是為了與低云區(qū)別開來。

出現(xiàn)滿足以上3個條件的情況，則判定為有雷暴發(fā)生。沒有雷暴發(fā)生的樣本記為1，有雷暴發(fā)生的樣本則記為2。如圖1所示，2014年6月23日05:40 UTC，區(qū)域1、2、3范圍內均有雷暴云團活動。

圖1 2014年6月23日05:40 UTC廣州白云機場125 km的MAX雷達回波圖 橫向和縱向1.0～17.0 km子圖為橫向和縱向雷達回波最大反射率對應的高度。三個圓圈分別表示以塔臺為中心，離塔臺8 km、50 km、100 km區(qū)域。右側文字說明為雷達的各種參數(shù)：時間(世界時)，雷達雜波濾波方法(離散傅里葉變換)，時間采樣，脈沖重復頻率(1 000 Hz/666 Hz)，范圍(125 km)，水平分辨率(0.417 km/pixel×0.417 km/pixel)，垂直分辨率(0.107 km/pixel)，高度范圍(1.0～17.0 km)。

2.3 樣本處理

將清遠探空站和香港探空站的探空數(shù)據(jù)進行距離線性插值法得到白云機場的探空數(shù)據(jù)。距離線性插值法如下：

其中，LAT和LNG(本場經緯度)、LAT1和LNG1(清遠探空站經緯度)、LAT2和LNG2(香港探空站經緯度)，a為清遠探空站到白云機場的距離，b為香港探空站到白云機場的距離，x為清遠探空站的探空數(shù)據(jù)，y為香港探空站對應清遠探空站的同一天同一時次的探空數(shù)據(jù)，Z為插值后的白云機場探空數(shù)據(jù)。

統(tǒng)計區(qū)域1、區(qū)域2范圍內、區(qū)域3范圍內雷暴發(fā)生與否數(shù)據(jù)共得三組數(shù)據(jù)。將同一天的00:00 UTC的探空數(shù)據(jù)與00:00—11:59 UTC雷暴發(fā)生與否數(shù)據(jù)匹配，12:00 UTC的探空數(shù)據(jù)與12:00—23:59 UTC雷暴發(fā)生與否數(shù)據(jù)匹配，因此每組樣本包含20個探空指數(shù)因子和1個雷暴發(fā)生與否數(shù)據(jù)。除去缺失資料的樣本，2014年的數(shù)據(jù)為183組，2015年的數(shù)據(jù)為182組，共計365組。其中區(qū)域1樣本組雷暴概率為135/365≈37%；區(qū)域2樣本雷暴概率為237/365≈65%；區(qū)域3樣本雷暴概率為277/365≈76%。

2.4 預報因子的選取

在篩選合適的輸入因子時，由于對流指數(shù)是連續(xù)型因子，雷暴發(fā)生與否為1、2非連續(xù)型變量，所以不能采用皮爾森相關系數(shù)一類的線性相關系數(shù)來衡量兩者的關系?；疑P聯(lián)度分析是研究因子間的幾何對應關系(即灰色關聯(lián)度)，序列曲線的幾何形狀的相似程度與灰色關聯(lián)度成正比[26-27]，它計算量小，適用于任何類型包括無規(guī)律的數(shù)據(jù)，在多因素，非線性領域得到廣泛應用[28-31]。

灰色關聯(lián)分析的步驟如下。

(1)確定雷暴發(fā)生與否數(shù)列和對流指數(shù)因子數(shù)列。設雷暴發(fā)生與否數(shù)列為X0′={X0′(k)|k=1，2，……，n}；對流指數(shù)因子數(shù)列Yi′={Yi′(k)|k=1，2，……，n}，i=1，2，……，m。這里m為指數(shù)因子個數(shù)20，n為樣本組數(shù)365。

(2)對指標數(shù)據(jù)進行無量綱化。由于各因子的物理意義不同，導致量綱不同，為了便于比較，采用均值化法(式(2))進行數(shù)據(jù)的無量綱處理。

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

(3)計算每個對流指數(shù)因子數(shù)列與雷暴發(fā)生與否數(shù)列對應元素的絕對差值。即|X0(k)-Xi(k)|。

(5)計算每個對流指數(shù)因子數(shù)列與雷暴發(fā)生與否數(shù)列對應元素的灰色關聯(lián)系數(shù)ξi(k)：

其中μ為分辨系數(shù)，此處取μ=0.5。

(6)計算灰色關聯(lián)度r0。分別計算各個指數(shù)與雷暴發(fā)生與否序列對應元素的均值即灰色關聯(lián)度r0i：

分別將20個對流指數(shù)因子與預報量(雷暴發(fā)生與否)作灰色關聯(lián)度分析，經過計算得到表3，表3為白云機場終端區(qū)區(qū)域1、區(qū)域2、區(qū)域3范圍內，20個對流指數(shù)與雷暴發(fā)生有無的灰色關聯(lián)度r0。區(qū)域1范圍內的r0平均值為0.968 9，大于平均值的指數(shù)由大到小依次為PLCL、SWI、BRCV、CT、KI、LFCV、TT、LFCT、BRN、VT、LCVT、CAPV、CAPE、LI。由于LI和LFCV、CAPE和CAPV物理意義近似，所以保留兩者中的大值LFCV、CAPV，刪去LI、CAPE。最終區(qū)域1選取了10個預報因子為PLCL、SWI、BRCV、CT、KI、LFCV、TT、LFCT、BRN、VT。同理，區(qū)域2的r0平均值為0.971 4，共選取了10個預報因子為PLCL、KI、CT、BRCV、SWI、LFCV、TT、VT、LCVT、CAPV。區(qū)域3的r0平均值為0.973 5，共選取了10個預報因子為

表3 白云機場終端區(qū)(區(qū)域1、區(qū)域2、區(qū)域3)20個對流指數(shù)與雷暴發(fā)生有無的灰色關聯(lián)度r0

PLCL、KI、BRCV、CT、TT、LFCV、SWI、VT、LCVT、CAPV。在白云機場終端區(qū)三個區(qū)域內，PLCL、KI、CT、BRCV、SWI、LFCV、TT、VT、LCVT、CAPV這10個對流參數(shù)與雷暴發(fā)生有無的灰色關聯(lián)度較大，其中PLCL、KI、CT、BRCV這4個指數(shù)的關聯(lián)度一直排在前五里，表明雷暴產生跟抬升凝結高度氣壓、850～700 hPa的溫度和露點溫度、風切變的關系最緊密。

3 隨機森林算法

隨機森林算法是一種包含多個無關聯(lián)的決策樹的分類器，并取所有決策樹中分類結果最多的那類為最終結果，其算法原理如下。

3.1 Bootstrap法重采樣

設集合Z為Z={Z(k)|k=1，2，……，n}，若每次有放回地從集合S中抽取一個樣本，一共抽取n次，形成新的集合Z*，則集合S*中不包含某個樣本Z(i)(i=1，2，……，n)的概率為：

因此，通過Bootstrap法重采樣，除去重復的樣本，新集合Z*約占原集合Z中(1-0.368)×100%=63.2%的樣本。

3.2 隨機森林算法流程

(1)基于Bootstrap方法隨機產生R個訓練集Z1、Z2、……、ZR。

(2)每個訓練集生成對應的決策樹C1、C2、……、CR；設樣本的屬性個數(shù)為A，從A個屬性中隨機抽取aa(aa為大于0且小于A的整數(shù))個屬性用最好的分裂方式對當前節(jié)點進行分裂，得到分裂屬性集。

(3)所有決策樹不剪枝。

(4)利用每個決策樹對測試集樣本X進行測試，得到對應的測試結果C1(X)、C2(X)、……、C R

(X)。

(5)將R個決策樹中輸出最多的分類類別作為測試集樣本X所屬的分類類別。

3.3 輸入輸出數(shù)據(jù)處理

為了對比分析，設立3個樣本輸入組：區(qū)域1樣本輸入組、區(qū)域2范圍內樣本輸入組、區(qū)域3范圍內樣本輸入組。其中區(qū)域1樣本輸入組包括365組10個預報因子，區(qū)域2范圍內樣本輸入組包括365組10個預報因子，區(qū)域3范圍內樣本輸入組包括365組10個預報因子。為了簡化計算，縮小量值，加快網絡收斂，將樣本輸入組進行[0.1，0.9]之間的歸一化：

其中，P為歸一化前的樣本輸入數(shù)據(jù)，Pmin為P的最小值，Pmax為P的最大值，Q為歸一化之后的矩陣。

樣本輸出組為365組雷暴發(fā)生有無數(shù)據(jù)。365組輸入輸出樣本中，將183組用來訓練網絡，182組用于檢驗。

4 結果分析

4.1 區(qū)域1預報結果及泛化分析

本文采用randomforest-matlab開源工具箱進行計算。區(qū)域1共有10個預報因子，因此樣本的屬性個數(shù)A為10，分裂屬性個數(shù)aa為的取整數(shù)值3。為了減少隨機性的影響，決策樹顆數(shù)T分別取{50，100，150，200，……，950，1 000}共20種情況，每種情況運行100次模型，取其準確值(即測試結果里預報值1，實況值1和預報值2，實況值2的總和除以182)的最高值作為當前決策樹顆數(shù)下的分類正確率。圖2為區(qū)域1樣本輸入隨機森林模型中，決策樹顆數(shù)對性能的影響，可看到當決策樹顆數(shù)為400時，分類正確率最高。

圖2 區(qū)域1樣本輸入隨機森林模型中，決策樹顆數(shù)對性能的影響

設立T=400，aa=3，運行隨機森林模型，將結果用兩變量預報驗證列聯(lián)表[32]表示得到表4。采用4個指數(shù)(如式組(8)所示)來評估預報結果，分別為：臨界成功指數(shù)CSI[33]、預報準確率AF[34]、虛假報警率FAR[32]和探測概率POD[32]。

表4 12 h區(qū)域1隨機森林分類器預報結果

隨機森林有一個重要的特點：沒有必要對它進行交叉驗證或用一個獨立的測試集來獲得誤差的一個無偏估計。它可在內部進行評估，即在生成過程中就對誤差建立了一個無偏估計。袋外錯誤率oob是隨機森林泛化誤差的一個無偏估計，它的結果近似于需要大量計算的k折交叉驗證。如式(9)所示，區(qū)域1的袋外錯誤率oob為31.32%，表明超過2/3的樣本能正確分類，泛化性能較好。

4.2 區(qū)域2分析結果

區(qū)域2共有10個預報因子，因此樣本的屬性個數(shù)A為10，分裂屬性個數(shù)aa為的取整數(shù)值3。圖3為區(qū)域2范圍內樣本輸入隨機森林模型中，決策樹顆數(shù)對性能的影響，可看到當決策樹顆數(shù)為250時，分類正確率最高。

圖3 區(qū)域2范圍內樣本輸入隨機森林模型中，決策樹顆數(shù)對性能的影響

設立T=250，aa=3，運行隨機森林模型得到的結果為表5和式組(10)。對比區(qū)域1預報結果，CSI、AF、POD值都要高，而虛假報警率FAR值偏低，表明雷暴預報效果更好了。當雷暴未來24小時發(fā)生可能性大于70%時，航空重要天氣MDRS通報業(yè)務將啟動。預報準確率AF和探測概率POD超過70%，開始對業(yè)務有指示作用。區(qū)域2范圍內的袋外錯誤率為28.57%，泛化性能有提升。

表5 12 h區(qū)域2隨機森林分類器預報結果

4.3 區(qū)域3范圍內分析結果

區(qū)域3共有10個預報因子，因此樣本的屬性個數(shù)A為10，分裂屬性個數(shù)aa為的取整數(shù)值3。圖4為區(qū)域3樣本輸入隨機森林模型中，決策樹顆數(shù)對性能的影響，可看到當決策樹顆數(shù)為100時，分類正確率最高。

圖4 區(qū)域2范圍內樣本輸入隨機森林模型中，決策樹顆數(shù)對性能的影響

設立T=100，aa為的取整數(shù)值3，運行隨機森林模型得到的結果為表6和式組(11)。三組預報結果里，區(qū)域3的預報結果CSI、AF、POD值最高，而虛假報警率FAR值最低，表明雷暴預報效果最好。區(qū)域3的袋外錯誤率為24.73%，泛化性能不錯。

表6 12 h區(qū)域3隨機森林分類器預報結果

5 結 論

本文將2014—2015年夏季雷達MAX圖資料獲取雷暴發(fā)生有無數(shù)據(jù)，結合美國懷俄明州立大學探空數(shù)據(jù)建立廣州白云機場終端區(qū)3個區(qū)域的12小時隨機森林分類模型，通過對比可得到以下結論。

(1)白云機場終端區(qū)區(qū)域1、區(qū)域2、區(qū)域3內，PLCL、KI、CT、BRCV、SWI、LFCV、TT、VT、LCVT、CAPV這10個對流參數(shù)與雷暴發(fā)生有無的灰色關聯(lián)度較大。雷暴的產生跟抬升凝結高度氣壓、850～700 hPa的溫度和露點溫度、風切變的關系最緊密。

(2)終端區(qū)區(qū)域面積越大，雷暴發(fā)生樣本比例越高，臨界成功指數(shù)CSI、預報準確率AF、探測概率POD越來越高，虛假報警率FAR越來越低，表明預報出來的準確率越來越高。區(qū)域2、區(qū)域3的預報準確率AF和探測概率POD超過70%，給航空天氣預報員分析預報雷暴提供了參考，對航空重要天氣MDRS通報業(yè)務有指示作用。

(3)終端區(qū)三個區(qū)域的測試樣本超過2/3能得到正確分類，表明隨機森林算法泛化性能較好。區(qū)域范圍越大，袋外錯誤率越低，泛化性能越好。