鄭超昊,尹志偉,曾鋼鋒,許月萍,周鵬,劉莉
(1.浙江大學 建筑工程學院,浙江 杭州 310058;2.臺州市水利局 浙江 臺州 318000;3.臺州市水利水電勘測設計院有限公司 浙江 臺州 318000)
降水是水循環(huán)的重要組成部分,也是地表水文過程的基本驅動元素[1].降水預報是進行流域水文預報、水資源規(guī)劃以及防洪防汛等研究的基礎資料.降水時空差異性大,因此如何獲取高精度、高分辨率的降水預報是水文領域一大難題[2].確定性降水預報無法完全模擬大氣過程和初始場資料觀測誤差的影響,降雨預報的模式逐漸從確定性預測轉向多成員集合預報,以集合的形式來描述未來降雨發(fā)生的概率情況[3].由于初值擾動、模式設計、集合成員數(shù)等方面的不完善因素,導致集合預報產(chǎn)品存在預報誤差、集合成員離散度低的問題[4]而無法直接使用,必須運用后處理方法提升集合預報的可靠性和準確性,才能保證高質量的降水預報信息供給[5].
在對數(shù)值降水預報的統(tǒng)計后處理方法研究上,代刊等[6]系統(tǒng)介紹了定量降水集合預報的統(tǒng)計學后處理技術;Zhao等[7]評估了分位數(shù)映射方法對季節(jié)性降雨預測的校正;Andrew等[8-9]發(fā)現(xiàn),參數(shù)化后處理方法貝葉斯聯(lián)合概率建模的校正效果較全面.深度學習方法的學習和泛化能力出色[10],具有非線性關系數(shù)據(jù)處理能力和數(shù)據(jù)稀疏地區(qū)較好表現(xiàn)的獨特優(yōu)勢,在偏差糾正后處理研究中被廣泛應用[11]:Qiu等[12]以多任務卷積神經(jīng)網(wǎng)絡考慮時間序列和站點之間的相關關系來進行降水預報,Krasnopolsky等[13]使用人工神經(jīng)網(wǎng)絡(artificial neural network, ANN)對降水集合預報進行后處理,皇甫江等[14]利用深度學習開展偏振雷達定量降水估測研究.在降水預報后處理的應用研究中,CNN和LSTM模型鮮有提及.
本研究提出結合CNN和LSTM的時空深度后處理模型(CNN-LSTM).該模型具有對時間序列數(shù)據(jù)進行分析的能力[15],能夠有效提取數(shù)據(jù)的空間信息特征特點,還能夠同時考慮預報降水數(shù)據(jù)空間特征和時間相關性.本研究以浙江省椒江流域為對象,將不同預報時效的歐洲中期天氣預報中心(European centre for medium-range weather forecasts,ECMWF)數(shù)值預報網(wǎng)格數(shù)據(jù)作為輸入數(shù)據(jù),以CMA-CMORPH實測雨量數(shù)據(jù)作為參照,對后處理前后降水資料進行全面評估.
椒江流域地處浙江中東部,東臨東海,介于東經(jīng)120°21′~121°47′,北緯28°30′~29°31′之間,主要位于臺州市境內(nèi),屬于亞熱帶季風氣候,多年平均氣溫為17.1 ℃,多年平均降雨量為1 634.9 mm.該區(qū)域的降水具有明顯的季節(jié)性,全年降水主要集中在4—6月梅雨季和7—10月臺風期;徑流量年內(nèi)分配主要集中在4—10月.永安溪和始豐溪是位于椒江上游的2條最主要河流[16].永安溪為椒江的源頭,長144 km,流域面積為2 704 km2.始豐溪是椒江流域二級支流,為椒江水系最大支流,主流長132.7 km,流域面積為1 610 km2.如圖1為椒江流域所處浙江省位置分布圖,以及永安溪、始豐溪流域位置以及水系分布圖.
圖1 椒江流域的流域分布、水系、降水格點數(shù)據(jù)示意圖Fig.1 Schematic diagram of basin distribution, water system and precipitation grid data in Jiaojiang river basin
1.2.1 ECMWF降水預報數(shù)據(jù) 使用的預報數(shù)據(jù)為ECMWF數(shù)值降水預報產(chǎn)品,從ECMWF的官網(wǎng)(www.ecmwf.int)下載得到.該數(shù)據(jù)由51個預報成員組成,空間分辨率為0.5°×0.5°,選取協(xié)調(diào)世界時零點(北京時間08:00)起預報的未來0~72 h逐6 h的降水預報進行后處理.考慮到降水預報可能存在的落區(qū)誤差,選取的經(jīng)緯度范圍為東經(jīng)120°~122°,北緯27°30′~29°30′,該范圍完全覆蓋椒江流域及其周邊區(qū)域.以2008—2018年汛期4—10月逐6 h降水資料作為輸入數(shù)據(jù),椒江流域在2012、2015和2016年遭受典型臺風事件,并且這些年份的降水數(shù)據(jù)完整性高,無需插補,故選定2012、2015和2016年作為測試期,其余時間的數(shù)據(jù)作為率定期.
1.2.2 CMA-CMORPH融合數(shù)據(jù) CMA-CMOPRH融合降水數(shù)據(jù)作為驗證數(shù)據(jù),從國家氣象信息中心下載得到,為中國30 000個自動站與CMORPH降水產(chǎn)品融合的逐時降水量網(wǎng)格數(shù)據(jù)集,擁有2008—2018年分辨率為0.1°的逐小時降水數(shù)據(jù).融合數(shù)據(jù)能夠有效減少系統(tǒng)偏差和隨機誤差,CMA-CMOPRH融合降水數(shù)據(jù)總體誤差水平不超過10%,優(yōu)于同類型產(chǎn)品在研究區(qū)域的精度[17].椒江流域47個實測雨量基本站點逐小時降水數(shù)據(jù)不連續(xù)且不充分,不足以支持深度學習模型的建模,而CMA-CMOPRH融合降水數(shù)據(jù)精度較高,數(shù)據(jù)較長,因此利用該降水產(chǎn)品替代實測數(shù)據(jù)對ECMWF數(shù)據(jù)進行后處理.
在CMA-CMORPH融合數(shù)據(jù)中,選取東經(jīng)120°30′~122°,北緯28°~29°30′,2008—2018年汛期4—10月逐小時降水格點數(shù)據(jù).將數(shù)據(jù)進行6 h累計計算(北京時間8:00點作為起始),篩出數(shù)據(jù)中的異常值和缺失值,利用min-max標準化方法進行數(shù)據(jù)歸一化,將6 h累計數(shù)據(jù)按照日期排序.將ECMWF數(shù)值預報降水數(shù)據(jù)轉化成0~72 h預報時效的6 h累計降水預報.篩出各預報時效內(nèi)數(shù)據(jù)的異常值和缺失值,利用min-max標準化對51個成員的各數(shù)據(jù)進行歸一化處理,再將處理后的數(shù)據(jù)按照日期排序,以保證數(shù)據(jù)的一致性.單個成員無法很好地代表ECMWF數(shù)據(jù)的整體性,而集合平均(等權重算術平均)數(shù)據(jù)可以涵蓋所有成員的預報能力,因此采用51個成員集合平均數(shù)據(jù)作為指標評判.將集合平均數(shù)據(jù)進行降尺度,使其分辨率為0.1°.利用實測降水資料的空間分布進行降尺度可以更好地捕捉降水的空間特征,但實際預報應用中不存在實測降水資料,因此本研究選用反距離權重插值法(IDW)進行數(shù)據(jù)降尺度,得到與CMA-CMORPH數(shù)據(jù)具有相同空間尺度的ECMWF數(shù)據(jù)用于后續(xù)評估.
將CMA-CMORPH柵格數(shù)據(jù)和ECMWF降水預報數(shù)據(jù)作為輸入數(shù)據(jù),CNN提取分析空間特性,LSTM獲取時間相關性[18],CNN-LSTM進行預報降水數(shù)據(jù)的時空誤差后處理.
2.2.1 模型基本原理 CNN基礎構成包含3種類型的層級:卷積層、池化層以及全連接層.卷積層通過卷積運算來提取輸入數(shù)據(jù)中的特征[18].在池化層中,濾波器將卷積層提取的特征信息進行壓縮,并且保留提取到的最顯著特征.CNN可以通過多層網(wǎng)絡學習不同感知域對目標網(wǎng)格和周圍區(qū)域的空間依賴性[19].LSTM是特殊的循環(huán)神經(jīng)網(wǎng)絡(RNN),它時序前饋連接的特性,使過去時刻的信息可以對當前時刻的輸出產(chǎn)生影響.相較于傳統(tǒng)ANN,LSTM具有4種結構:更新門、輸入門、遺忘門和輸出門.更新門用于儲存過去信息的累計;遺忘門可以選擇將過去信息中影響較小的因素“遺忘”,使模型累計的系統(tǒng)誤差減少,得到比傳統(tǒng)RNN更優(yōu)的結果[19].
如圖2所示為CNN-LSTM的模型結構,該模型是時間序列數(shù)據(jù)集從輸入到結果輸出的降水后處理過程.
圖2 CNN-LSTM的模型結構Fig.2 Model structure of CNN-LSTM
2.2.2 模型結構 CNN-LSTM包括CNN空間特征提取、LSTM時間相關性分析2個模塊.每個空間位置的降水不孤立存在,而是與周圍空間相關聯(lián),利用CNN模塊提取目標周圍區(qū)域的降水信息,為當前目標點的修正提供參考.過去時間的降水信息會影響當前時刻降水,為此,利用LSTM模塊時序前饋的特點提取降水前后方向的時間相關性.由此可知,CNN-LSTM能夠兼顧時間、空間,可以運用于降水后處理中.
從CNN模塊來看,由于預報數(shù)據(jù)的不確定性會導致單獨成員的預報結果有較大的偏差,將51個預報成員的數(shù)據(jù)融合,生成數(shù)據(jù)集{51, 5, 5},并將率定期所有數(shù)據(jù)組連接成連續(xù)時間序列數(shù)據(jù)集作為耦合模型的輸入數(shù)據(jù),即該數(shù)據(jù)集包括不同成員的預報信息要素,還能包含空間和時間信息.輸入數(shù)據(jù)先通過二維卷積層Conv2D,再通過最大值池化層MaxPool提取最大值特征,利用Dropout方法防止模型過擬合.先將數(shù)據(jù)根據(jù)設定參數(shù)進行卷積、最大值池化,再向量平鋪成一維向量.將CNN模塊輸出得到的向量與神經(jīng)節(jié)點為256個的全連接層相連,再將全連接層轉化成(16×16)的中間結果輸入LSTM模塊.從LSTM模塊來看,16×16的二維中間向量,分別經(jīng)過神經(jīng)元結點為64、32、16的雙向LSTM模型,每一步都通過Dropout方法,再由激活函數(shù)Relu輸出16×16的二維向量與實測數(shù)據(jù)進行訓練.
本研究模型的優(yōu)化器采用Adam,訓練迭代次數(shù)為1 000,每組訓練的數(shù)據(jù)量為64,耦合模型的學習率最初為0.001,每迭代400次其學習率變?yōu)樵瓉淼?.1,模型損失函數(shù)為均方誤差MSE、兩幅圖像相似度的指標SSIM兩者加權的結果,權重設置分別為0.7、0.3.
式中:x、y均為對比圖像,μx、μy分別為x、y像素的平均值,σx、σy分別為x、y像素方差,σxy為x、y的像素協(xié)方差,c1、c2均為穩(wěn)定常數(shù).
以下4個性能指標用于定量評估CNN-LSTM的后處理能力:均方根誤差RMSE、相對偏差RB、相關系數(shù)CC、平均絕對誤差MAE.指標的計算式[20]分別為
式中:n為各個網(wǎng)格點觀測值的數(shù)量,為網(wǎng)格預測降水,為網(wǎng)格實測觀察降水,為網(wǎng)格預測降水的平均值,為網(wǎng)格實測觀察降水的平均值.
為了進行分級檢驗評估,在ECMWF降水預報資料中提取24 h累計降水,計算永安溪、始豐溪、平原區(qū)域面雨量.根據(jù)《中短期天氣預報質量檢驗辦法》降水等級劃分如表1所示,將降水預報和降水實況按照24 h累計降水量劃分成對應等級.對各降水分級進行檢驗,各等級TS(threat score)評分為
表1 降水等級劃分表Tab.1 Classification of precipitation levelsmm
式中:k為1~4級,分別代表小雨、中雨、大雨、暴雨及以上降水預報;NAk為k級降水預報正確次數(shù);NBk為k級降水空報次數(shù);NCk為k級降水漏報次數(shù).ECMWF預報數(shù)據(jù)采用集合平均的數(shù)據(jù)進行TS評分計算.
評估驗證期內(nèi)不同預報時效的預報后處理數(shù)據(jù),對比分析不同預報時效內(nèi)永安溪、始豐溪、平原區(qū)域面雨量;選取驗證期內(nèi)的多場典型暴雨進行分析,著重分析永安溪、始豐溪、平原區(qū)域面雨量.
計算12個預報時效內(nèi)每個網(wǎng)格點的數(shù)據(jù),對比后處理前、后預報數(shù)據(jù)的指標來評估CNNLSTM.流域全部格點不同預報時效的各指標計算結果如圖3(小提琴圖)所示,其中虛線為25%、50%、75%分位線,t為預報時效.原始預報RMSE普遍超過3 mm,最大值超過6 mm;經(jīng)CNNLSTM后處理,各個預報時效的RMSE均有所降低,基本不超過2 mm,最大值為3 mm.經(jīng)CNNLSTM后處理,RB絕大部分落在-0.25附近區(qū)域,而原始預報的RB跨度大.可以看出,原始預報并不是預報時效越短預報技能越高,且預報技能也不隨著預報時效增長呈基本的單調(diào)增減,這是氣象預報領域較為常見的“預報不一致性”現(xiàn)象,即在連續(xù)多次預報中,前后兩次預報的預報結果差異較大.該現(xiàn)象被認為主要與初始誤差和模式誤差有關,相關分析可參考文獻[21].
圖3 模型評估指標在所有格點不同預報時效后處理前、后的對比圖Fig.3 Comparison of model assessment indexes before and after post-processing for all grids with different lead times
如圖4所示為不同預報時效后處理前、后各指標的空間分布圖.可以看出,不同預報時效后處理后的RMSE、RB、MAE、CC在全流域均有很大的提升.其中RMSE從2.6~6.0 mm下降為1.0~3.0 mm,后處理前、后RMSE下降的格點比例為100%;RB從-0.6~1.0穩(wěn)定至-0.35~-0.15;MAE從0.7~1.8 mm下降為0.3~0.7 mm,CC從0.2~0.5上升至0.70~0.90;后處理前、后MAE、CC指標改善的格點比例為100%.綜合所述,CNN-LSTM對于數(shù)值預報降水產(chǎn)品不同預報時效的預報精度均有顯著提升,并且對于整個椒江流域空間精度提高具有普遍性.
圖4 不同預報時效后處理前、后模型評估指標的空間分布圖Fig.4 Spatial distribution of model assessment indexes before and after post-processing for different lead times
如圖5所示,永安溪、始豐溪流域以及平原區(qū)域不同預見期面雨量結果相似.圖中,下標數(shù)字0代表后處理前,下標數(shù)字1代表處理后.各個預報時效的RMSE均落在基準線下側,說明后處理后數(shù)據(jù)的RMSE均小于后處理前數(shù)據(jù)的RMSE,且后處理前RMSE為2.2~4.0 mm,處理后RMSE為1.0~2.0 mm.各個預報時效的CC均落在基準線上方,后處理前的CC為0.0~0.4,處理后的CC>0.8.各預報時效指標的MAE均落在基準線下方,后處理前的MAE為0.75~1.75 mm,處理后的MAE為0.25~0.75 mm.各個預報時效的RB分布在基準線兩側,其中后處理后的RB穩(wěn)定在-0.25附近,后處理前的RB為-0.75~0.25,可見處理后的RB穩(wěn)定在更合理的范圍.綜上所述,在永安溪、始豐溪流域,4個模型評估指標經(jīng)過CNN-LSTM模型后處理后的預報面雨量數(shù)據(jù)較于后處理之前在數(shù)據(jù)精度方面有很大程度的提升.
圖5 模型評估指標在永安溪、始豐溪流域面雨量不同預報時效后處理前、后的對比圖Fig.5 Comparison of model assessment indexes of areal rainfall before and after post-processing over Yonganxi and Shifengxi watershed for different lead times
如表2所示為計算所得永安溪、始豐溪、平原區(qū)域的各量級降水出現(xiàn)的累計頻次N.如表3所示為ECMWF降水預報后處理前、后永安溪、始豐溪、平原區(qū)域的TS評分.由表可知,經(jīng)過后處理3個區(qū)域的TS分別為0.900、0.906、0.894,均高于后處理之前的0.754、0.755、0.752,分別提升了19.4%、20%和18.9%.結果證實經(jīng)CNN-LSTM后處理的雨季預報正確率有顯著提升.
表2 永安溪、始豐溪、平原區(qū)域降水分級累計出現(xiàn)頻次Tab.2 Cumulative occurrence frequency of precipitation classification in Yongan, Shifeng and plain regions
表3 ECMWF預報產(chǎn)品后處理前、后雨季3個區(qū)域的TS評分項Tab.3 TS scores items of ECMWF forecast products before and after post-processing at three regions during rainy season
如圖6所示為各區(qū)域不同降水級別的TS評分圖.后處理后的各等級累計降水TS評分較之后處理前的均有很大程度提升,尤其是小雨級別的TS,數(shù)值由0.77提升至0.91,中雨級別的TS從0.48提升至0.6,大雨級別的TS從0.3提升至0.4,暴雨級別的TS評分也有所提升.在后處理前、后的ECMWF預報產(chǎn)品在降水等級評估中,小雨的預報正確率最高,中雨次之,暴雨的預報正確率較低,這可能是暴雨事件較少導致的.
圖6 ECMWF預報產(chǎn)品后處理前、后的等級累計降水TS評分對比圖Fig.6 Comparison of TS scores of ECMWF forecast products before and after post-processing
選取2012年8月“??迸_風、2015年8月“蘇迪羅”臺風作為典型暴雨,分析永安溪流域、始豐溪流域、平原區(qū)域未來3 d累計6 h面雨量.選取2012年8月6日8:00起未來72 h數(shù)據(jù)作為“海葵”臺風的預報數(shù)據(jù),并且利用分位線代表集合預報中95%置信區(qū)間.由圖7可知(圖中,p為降水量),本場臺風實際暴雨中心位于流域以北區(qū)域,流域內(nèi)最大實測降水量為362.8 mm,原始預報降水中心偏離,最大雨量僅為226.6 mm,小于實測降水量;后處理后預報暴雨中心與實測降水一致,流域內(nèi)累計降水量為332.6 mm.選取2015年8月7日8:00起未來72 h數(shù)據(jù)作為“蘇迪羅”臺風的預報數(shù)據(jù).本場臺風的暴雨中心落在椒江流域,實測累計降水的降水值為338.1 mm,原始預報降水最大雨量為308.6 mm,后處理后的預報降水最大值約為310.4 mm.由對比可知,在2場臺風期間,后處理后的預報數(shù)據(jù)不管是降雨累計,還是暴雨中心都比處理前的更為精確.后處理后的預報數(shù)據(jù)對比原數(shù)據(jù)有相當程度的提升.
圖7 2場臺風72 h累計降水分布圖Fig.7 Distribution of 72h cumulative precipitation for two typhoons
如圖8所示為臺風期間永安溪、始豐溪流域和平原區(qū)域逐6 h面雨量情況,為了進一步評估后處理前、后預報的差異,圖中增加原始集合預報的95%的置信區(qū)間結果.在“海葵”臺風期間,3個區(qū)域最大累計實測降水量分別為35、40、35 mm,預報后處理后累計降水分別為32、36、33 mm;三者預報后處理后的累計降水誤差均不超過10%.在“蘇迪羅”臺風期間,3個區(qū)域最大累計實測降水量分別為50、60、70 mm,后處理后累計降水分別為45、54、64 mm.三者預報后處理后的累計降水誤差均不超過10%,后處理前的累計降水誤差較大.
圖8 2場臺風在3個區(qū)域的逐6 h面雨量對比圖Fig.8 Comparison of 6 h areal rainfall of two typhoons at three regions
對2場臺風預報進行評估,結果如表4所示.在“??迸_風期間,3個區(qū)域后處理后的RMSE下降均超過90%,MAE在后處理后也有大幅度下降, CC在后處理后穩(wěn)定在0.99,RB在后處理后提升效果均超過89%.在“蘇迪羅”臺風期間,3個區(qū)塊后處理后的RMSE分別下降了91.3%、90.7%、91.6%;MAE分別下降了91.3%、90.2%、91.7%;后處理前的CC不足0.55,后處理后的均穩(wěn)定在0.99;RB分別提升了91.3%、90.1%、91.5%.結果表明CNN-LSTM后處理的數(shù)據(jù)模擬典型臺風的效果良好,面雨量模擬精度大幅度超過未經(jīng)后處理的數(shù)據(jù).結果也表明原始數(shù)值預報降水產(chǎn)品對臺風降水預報能力較弱,通過耦合模型后處理后的降水產(chǎn)品對臺風降水預報能力較強.
(1)相較于原始數(shù)值天氣預報數(shù)據(jù),后處理后的RMSE的最大值從6 mm降至3 mm;CC的最大值從0.6升至0.9,RMSE、CC的改善格點比例為100%.說明CNN-LSTM能夠有效提高降水預報精度.
(2)后處理后永安溪、始豐溪流域和平原區(qū)域不同預報時效面雨量的指標有顯著提升:RMSE全面提升,最大值從4 mm降至2 mm;后處理后的CC均超過0.8;MAE從0.75~1.75 mm降至0.25~0.75 mm;RB后處理后的值優(yōu)于后處理前的,并穩(wěn)定在-0.25.
(3)后處理后永安溪、始豐溪以及平原區(qū)域的雨季降水預報的TS評分均大于0.89,高于后處理前的TS評分;后處理后的各降水等級TS評分均好于后處理前的,其中小雨TS評分從不足0.8提升至0.91,中雨的TS評分從不足0.5提升至0.6,大雨TS評分從不足0.4提升至0.4.
(4)2場典型臺風降水后處理后的累計降水量和暴雨中心評估準確性高.對于永安溪、始豐溪流域、平原區(qū)域逐6 h面雨量降水峰值的預報數(shù)據(jù)結果誤差均不超過10%,各類指標評提升均超過85%,說明經(jīng)過CNN-LSTM后處理的ECMWF降水產(chǎn)品可以用于臺風降水預報.
(5)CNN-LSTM適用于數(shù)值降水預報后處理,能夠提高定量降水預報的準確性,對防汛減災都具有十分重要的意義.本研究僅針對南方濕潤地區(qū)得出結論,后續(xù)將針對北方地區(qū)開展相關研究.