閆文杰 劉圣軍 劉新儒 彭謙 胡婭敏

(1. 中南大學數學與統計學院,長沙,410075?2. 中國人民大學統計學院,北京,100872?3. 廣東省氣候中心,廣州,510080)

1 引言

華南地區(qū)地處中國南部,屬熱帶、亞熱帶季風氣候區(qū),具有非常顯著的季風氣候特征,形成暴雨的天然條件位居我國前列[1],故暴雨發(fā)生較為頻繁,其中初夏為暴雨頻數最多的時段之一. 研究華南地區(qū)暴雨日數的時空分布特征,并選取合適模型進行暴雨日數的預測,在極端天氣時間預測方面具有重要的應用價值.

對于暴雨日數預測問題,前人已進行過許多相關研究. 在因子選取方面,簡茂球等[2]在2013 年利用相關系數法指出海南持續(xù)性暴雨天氣與熱帶大氣的準雙周振蕩有關,并將海溫納入分析范疇,發(fā)現2010 年在赤道中東太平洋發(fā)生的La Ni?a 現象有利于產生大氣準雙周低頻振蕩. 在抽樣方法方面,為解決暴雨預測中類別不平衡問題,楊艷等[3]將減抽樣的思想應用于AdaBoost 綜合學習算法[4]中. 實驗結果表明,相對于傳統的AdaBoost 算法[5],該算法對銅川暴雨日數的預報準確率顯著提高. 在預測方法方面,劉綠柳等[6]應用一步法和兩步法兩種統計降尺度方法預測暴雨日數,并提出了用標準差作為評估暴雨日數異常等級的評分標準.文中交叉檢驗的結果顯示這兩種方法對于月尺度降水與暴雨日數的預測都具有較好的效果.

除了以上傳統的方法, 近年來人工神經網絡在氣象領域也得到了廣泛應用. 葛彩蓮等[7]在2010 年應用BP 神經網絡對降雨量進行預測, 發(fā)現使用前5 年降雨量來預測后1 年的降雨量的效果與用氣象資料預測降雨量的效果相近,讓這種新預測模式應用于降水預測成為可能. 黎玥君等[8]在2017 年基于BP 神經網絡預測浙北夏季降水,發(fā)現神經網絡隱藏層節(jié)點數量為2 時,擬合降水量的效果最佳. 關鵬洲等[9]在2017 年提出了基于梯度提升決策樹(GBDT)特征選擇的改進AdaBoost 回歸模型、基于多個強回歸器組合的stacking 回歸模型,和Inception 卷積神經網絡對數據集進行分類和回歸,獲得了對短期降雨較好的預測效果.

前人的研究成果從數據處理到模型構建,都為本文提供了寶貴的思路,但也存在亟待改善的地方. 大多數研究工作采用EOF 分解、獨立分量分析等傳統機器學習方法提取特征，選用常用的平均絕對誤差、相對偏差[8]、Ts 評分[3]等指標進行模型檢驗.然而，傳統方法算法簡單，更適合于線性關系建模，且常用檢驗指標不能全面有效地對結果進行評估.為改善上述問題,本文引入新的特征提取方法與檢驗方法.

2 研究方法

2.1 數據

降水數據來源于中國氣象局國家氣象信息中心和廣東省氣候中心數據庫. 本文選用1981 年至2018 年華南地區(qū)192 個站點的數據. 選取的5 個環(huán)流要素場來自于NCEP/NCAR 再分析資料的月數據，分別為: 高度場(hgt),經向風場(vwnd),緯向風場(uwnd),比濕(shum),海平面氣壓(slp),其中水平分辨率為2.5°×2.5°?海溫數據選用NOAA ERSST V3b 海表溫度(SST),水平分辨率為2.0°×2.0°?海冰數據來自NOAA 北半球逐月海冰資料?積雪逐月數據由Rutgers University Global Snow Lab 提供(包括北美洲、北半球、歐亞大陸積雪數據). 數據如表1所示.

表1 詳細數據表格

根據廣東省氣候中心暴雨閾值標準,24 小時單站降水量超過50mm 的日數定義為暴雨日數.本文主要研究提前一個月的要素場對華南地區(qū)初夏暴雨日數的影響,因此選取當年5 月的要素作為預報因子來預測華南地區(qū)初夏暴雨日數.

對于部分預測因子,由于缺測或無法測量(如海溫數據只能在海洋區(qū)域測出,陸地區(qū)域無法測量)等原因,會在某些格點或時間點出現異常數據. 本文對異常值直接剔除. 為了消除預報因子之間的量綱影響,本文采用z score 方法對數據進行標準化處理:

其中,zij為對每個格點(第i行,第j列)數據求標準化后的值,xij為初始值, ˉxi為第i行數據的平均值,si為第i行數據的標準差.

2.2 統計降尺度方法框架

統計降尺度法是利用大氣環(huán)流的觀測資料建立大尺度氣候要素和區(qū)域氣候要素之間的統計關系,并把這種關系應用于大氣環(huán)流模式中輸出大尺度氣候信息,進而預測區(qū)域未來氣候變化的一種常用方法[10]. 其基本思路是建立大尺度預報因子和區(qū)域氣候預報量之間的統計關系函數[11]:

其中,y是區(qū)域氣候預報量(氣溫,降水等),x是大尺度預報因子,F(·)為建立的統計關系函數,使用最廣泛的是轉換函數. 統計降尺度基本流程如圖1所示: 在暴雨日數預報問題中,由于選取的預報因子較多,且環(huán)流預報因子為二維格點數據,樣本數量有限,傳統機器學習模型會存在嚴重的過擬合現象. 緩解過擬合的方法主要有兩種[12]:

圖1 統計降尺度方法框架圖

(1)降低預報因子維度. 可通過篩選顯著相關因子和選取高相關區(qū)域格點數據并取區(qū)域平均值的方法來減少維度.

(2)將模型正則化,減少特征參數的數量級.

本文采用高相關區(qū)域選取的方法來提取特征,使用最小角回歸的方法獲取預報因子與預報量之間的轉換函數. 接下來分別詳細討論這兩個方法.

2.3 高相關區(qū)域提取方法

再分析資料數據是包含時間和空間信息的多維格點數據. 預測過程中我們只關注與華南初夏暴雨有相關性的格點數據,因此通常會將某個時間點的空間數據重組成一維數據. 例如,由于本文中環(huán)流場的經緯度范圍分別為0°–175°E,20°S–55°N,分辨率為2.5×2.5,某一時間點的200hPa 高度場數據經度方向71 個數據，緯度方向31 個數據,進行一維重組后,成為維度為2201 的一維向量.但由于相鄰格點數據相關性很高,必然導致模型產生嚴重的多重共線性及過擬合現象. 因此,我們需要對數據進行預處理. 提取高相關區(qū)域是一種有效的數據降維方法. 本文使用預報因子各格點與單站點時間序列數據的線性相關系數來提取預報因子的高相關區(qū)域,其選取操作流程如圖2所示.

圖2 高相關區(qū)域選擇操作流程

以高度場數據與廣東佛岡站點的初夏暴雨日數的相關性分析為例,我們選擇時間范圍為1981年5 月至2018 年5 月,空間范圍為0°–175°E,20°S–55°N 的500hPa 高度場格點數據與廣東佛岡站點初夏暴雨日數進行相關性分析,得到空間相關情況如圖3所示. 設定顯著性水平為0.05，可得高相關區(qū)域數為6，即有6 個區(qū)域的格點數據與因變量的相關系數較大. 由此可得，自變量的維數由原來的2201 維變成6 維.

圖3 500hPa 高度場與暴雨日數相關情況(以佛岡站點為例)

2.4 最小角回歸預測模型

在常用逐步回歸方法中,參數選擇可采用前向選擇算法和前向梯度算法. 然而,這兩種算法比較暴力,效率較低. Bradly Efron 在2004 年發(fā)表的文章中提出了一種新的算法——最小角回歸(LAR)算法[13]. 該算法保留了前向梯度算法的精確性,同時簡化了迭代的過程,步驟如下[13,14]:

Step 1: 設有n個經過了標準化的自變量xk(k= 1,2,··· ,n),中心化的因變量y. 計算所有自變量與y的相關系數并排序,選出相關系數最大的一個自變量,不妨設為x1,滿足

其中r(x,y)表示x與y之間的相關系數. 此時將x1加入逼近y的特征集合中.

Step 2: 在x1方向上用x1逼近y,選擇步長θ1,得到下列回歸方程

其中θ1表示回歸系數, ˉy表示y在x1方向上的斜向投影,且

定義殘差

選取θ1使得存在一個未被選取的自變量,不妨設為x2,滿足

此時,殘差yres位于x1和x2的角平分線上. 將x2加入逼近y的特征集合中.

Step 3: 在Step 2 得到的角平分線方向上前進步長θ2,按上述方法更新殘差yres,使得存在另一個未被選取的自變量,不妨設為x3,滿足

Step 4: 循環(huán)上述步驟, 直到殘差‖yres‖2小于給定的值? 或者已經遍歷了所有自變量xk(k=1,2,··· ,n),算法停止.

從幾何角度分析,每一步選擇的前進路徑必須保證已選入模型的變量xk與殘差yres的角度最小(即“最小角”). 因此算法每一次都選擇原路徑與新變量夾角的角平分線方向作為新的前進路徑方向.

最小角回歸算法是一個適用于高維數據的回歸算法. 該算法的最壞時間復雜度和最小二乘法類似,但計算速度卻能與前向選擇算法一樣快. 同時,該方法可以產生分段線性結果的完整路徑,在模型的交叉驗證中極為方便.

2.5 評價檢驗

2.5.1 單站點預測結果檢驗評估指標

氣象學中有許多檢驗評估指標,包括時間距平相關系數(TCC)、空間距平相關系數(ACC)、相對操作特征(ROC)等. 其中TCC原理簡單,運用廣泛,且適用于本文研究問題的單站點輸出結果分析,因此可以作為主要的檢驗指標. 符號一致率,也稱同號率(SS),能夠把握距平值的變化趨勢(升高或者降低). 訓練集決定系數()可以檢驗模型訓練集的擬合情況.測試集決定系數()與訓練集決定系數原理類似,可以檢驗測試集的擬合情況,同時驗證TCC指標的合理性,但標準較為嚴苛.

時間距平相關系數的公式為其中,n1為模型測試集的樣本數量,TESTi為測試集中第i個測試數據,OBSi為與測試數據對應的第i個真實數據,為測試集數據的平均值,為與測試集對應的真實數據的平均值.TCC取值范圍是[?1,1],越接近1,說明預測效果越好.

同號率公式為

其中,I(·)表示示性函數.

測試集決定系數為

訓練集決定系數為

其中,n2為模型訓練集的樣本數量,TRAINj為訓練集中第j個測試數據.

2.5.2 區(qū)域預測結果檢驗評估指標

趨勢異常綜合評分(Ps)是國內氣象領域常用的評價標準,它綜合考慮了同號率、擬合能力以及異常值(偏離平均值較多的觀測值)預報能力. 本文在保留Ps評分部分標準的同時,針對華南地區(qū)暴雨日數預測問題,對評分標準進行調整,將不連續(xù)的分段常數函數調整為連續(xù)函數,得到調整Ps評分(APs).Ps評分標準基于距平百分率進行操作,APs 評分沿用這種模式. 本小節(jié)涉及的變量均為對應數據的距平百分率.

設共有m個站點, 每個站點都有n個測試樣本, 即n年的月平均暴雨日數數據. 當OBS(i,k)*TEST(i,k)< 0, 即第i年第k個站點的觀測值與對應預測值符號相反時, 定義第i年,第k個站點的得分為零

當OBS(i,k)*TEST(i,k)≥0,即第i年第k個站點的觀測值與對應預測值符號相同時(包含其中某項為0 的情況),定義第i年,第k個站點的得分公式為

其中,α1,α2為經驗參數,用于避免出現分子或分母為0 的情況(本文通過多次實驗比較,將α1與α2均設為0.05). 上式分子部分用于對異常值的預報能力,分母部分衡量模型對觀測值的擬合程度.

與Ps評分準則類似,構造第i年所有站點的綜合APs評分

3 實驗結果

在模型訓練過程中,本文將數據按7:3 的比例分為訓練集和測試集. 時間距平相關系數,同號率和測試集決定系數用于評估測試集上的預測效果,訓練集決定系數是訓練集上的檢驗標準. 本小節(jié)涉及到的變量均為距平值.

3.1 單站點預測結果比較

本文分別基于主成分分析法和高相關區(qū)域選取法對數據降維和提取特征,并使用最小角回歸算法對站點暴雨日數進行預測.

我們將提取高相關區(qū)域和主成分分析兩種提取特征的方法進行對比分析. 主成分分析是一種使用最廣泛的數據降維算法,其基本思想是將n維特征映射到k(k<

對自變量進行主成分分析,用協方差矩陣的特征值計算累計方差貢獻率,從而得到累計貢獻率與主成分數量關系,如圖4所示.可以看出,累計貢獻率Ψm與主成分數量m的函數關系為凹函數,累計貢獻率的增長速度隨著主成分數量逐步引入而降低. 當引入第11 個主成分后,累計貢獻率高于0.7?當引入第16 個主成分后,累計貢獻率高于0.8?當引入第23 個主成分后,累計貢獻率高于0.9?累計貢獻率在主成分數量為37 時達到最大值1. 本文選用23 個主成分,樣本貢獻率高于0.9.

圖4 累計樣本貢獻率與主成分數量關系

對預報因子進行高相關區(qū)域提取時,根據最小角回歸算法特性,相關系數檢驗的顯著性水平設為0.05,以增強自變量集的解釋能力. 我們使用1981 年至2006 年的數據進行訓練,將2007 年至2018 年的數據進行測試,再與觀測值比對. 設置殘差上限為2×10?16,非零參數數量設置為15. 實驗結果如圖5所示.

圖5 最小角回歸預測初夏暴雨日數距平

用時間距平相關系數(TCC),同號率,訓練集決定系數和測試集決定系數,分別對基于兩種特征提取方法的最小角回歸算法結果進行檢驗,并將得到的四個結果進行對比分析,結果如表2所示,其中紅色的表示最好結果、藍色的表示次好結果.

由表2 的數據可得,對于兩個站點的模型輸出結果,除廣州站點的訓練集決定系數,使用主成

表2 廣州站點與海口站點初夏暴雨日數預測結果檢驗

分分析的方法略優(yōu)于高相關區(qū)域提取的方法外,其他所有的檢驗指標都是高相關區(qū)域提取方法明顯占優(yōu). 這說明基于高相關區(qū)域提取的最小角回歸算法具有較高的預測精度.

3.2 區(qū)域的分析比較

本節(jié)從整個華南區(qū)域的角度對預測結果進行分析. 為了更好地分析最小角回歸算法的預測效果,與分別采用基于主成分分析和高相關區(qū)域提取兩種特征選擇的多元線性回歸方法進行對比. 為了方便表示,將兩種特征選擇方法下的兩種預測模型分別用代號A,B,C,D 進行表示,具體對應規(guī)則見表3.

表3 代號及對應模型

選取氣象學中常用的檢驗指標TCC,分析比較A,B,C,D 四種模型在測試集上的得分. 將所有站點TCC得分分布情況列出如表4所示. 表中第一列表示TCC的劃分區(qū)間,第二至五列的數據表示TCC得分落在對應區(qū)間的站點數占總站點數的百分比.

表4 四種模型下華南地區(qū)192 站點TCC 得分分布百分比

統計數據表明,模型A,C 的分布百分比隨著TCC的增加,呈現先增后減的分布趨勢?模型A中站點的TCC主要集中于區(qū)間[?0.2,0.8],模型C 中站點的TCC主要集中于區(qū)間[?0.2,0.6],說明基于主成分分析方法選擇特征的模型預測效果整體較差. 模型B,D 的站點百分比隨著TCC增加大致呈遞增趨勢,兩種模型中站點的TCC均集中于效果最佳的[0.6,1]區(qū)間,且占比均為50%左右,不同的是模型B 在[0.8,1]區(qū)間的占比更高,而模型D 在[0.6,0.8]區(qū)間的占比更高. 這表明基于高相關區(qū)域提取方法的兩種預測模型在TCC上的表現接近.

上表從空間尺度上比較了模型的預測效果. 為了進一步從時間尺度上比較分析,對四種模型2011 年到2018 年的測試集進行APs評分,結果如圖6所示.

由圖6可知,整體而言模型B,D 處于更高的評分區(qū)間. 模型D 的APs評分均值最高為0.9359,所有年份均顯著地高于其他模型,且方差最小(1.7×10?4),測試集的預測結果能穩(wěn)定地位于較小的區(qū)間內,整體效果均顯著優(yōu)于其他模型. 模型B 預測效果次之,APs評分均值為0.8871,方差為1.3×10?3. 模型A,C 的測試集評分均值最低,方差最大,說明這兩種模型對華南站點初夏暴雨日數預測精度不足,且穩(wěn)定性較差. 四種模型的APs評分對應均值與方差如表5所示.

表5 四種模型APs 評分對應均值與方差

圖6 四種模型測試集APs 評分對比

4 結論

本文旨在研究前期的大氣環(huán)流因子(高度場(hgt), 經向風場(vwnd), 緯向風場(uwnd), 比濕(shum),海平面氣壓(slp))及海溫、海冰和積雪數據對華南地區(qū)初夏暴雨的影響,從而實現對華南地區(qū)初夏暴雨日數的預測. 本文首先使用高相關區(qū)域提取的方法對數據降維,并與常用的主成分分析方法進行比較. 然后用兩種方法選擇的特征分別構建華南地區(qū)初夏暴雨日數的機器學習預測模型,并運用以時間距平相關系數為主的單站點預測結果檢驗標準和以APs評分為主的區(qū)域性檢驗標準對混合構建的四種模型進行檢驗分析,得到以下兩點結論.

(1)在特征提取方面,從單站點預測結果分析,基于高相關區(qū)域提取方法的模型結果基本反映了真實距平值的變化,基于主成分分析方法的模型結果與真實情況偏差較大. 其他4 種評價指數也驗證了基于高相關區(qū)域提取的最小角回歸預測模型的明顯優(yōu)勢. 從區(qū)域預測結果分析可見,基于高相關區(qū)域提取的兩種模型的時間相關系數分布在區(qū)間[0.6,1],基于主成分分析方法的兩種模型的時間相關系數近似正態(tài)分布,主要位于區(qū)間[0.2,0.6]. 且高相關區(qū)域提取方法的兩種模型的APs評分主要在區(qū)間[0.88,0.96],而基于主成分分析方法的兩種模型的APs評分主要在區(qū)間[0.70,0.89].

(2)在預測模型方面,基于高相關區(qū)域提取方法的最小角回歸算法預測結果的時間相關系數分布在區(qū)間[0.6,0.8],而相應的多元線性回歸模型的時間相關系數相對均勻地分布在區(qū)間[0.6,1],分布相差不明顯. 但高相關區(qū)域提取方法的最小角回歸算法的APs評分主要在區(qū)間[0.91,0.96]上,而對應的多元線性回歸模型的APs評分主要在區(qū)間[0.81,0.93].

綜上,基于高相關區(qū)域提取方法的最小角回歸算法能較好地預測華南初夏暴雨日數.