考慮降水時間相關性的地面觀測-雷達-衛(wèi)星遙感逐時降水融合方法研究

阮惠華，張鈞民 ，許劍輝，戴曉愛

(1.廣東省氣象探測數據中心，廣東 廣州 510641；2.廣東省科學院廣州地理研究所/廣東省遙感與地理信息應用重點實驗室/廣東省地理空間信息技術與應用公共實驗室，廣東 廣州 510070；3.成都理工大學地球科學學院，四川 成都 610059)

1 引 言

隨著氣象觀測系統(tǒng)的迅猛發(fā)展，通過地面自動氣象站、雷達、衛(wèi)星遙感等手段獲取的不同時空分辨率的降水觀測數據越來越多[1]；各行業(yè)對高質量的時空連續(xù)降水數據產品要求越來越高，特別是針對極端天氣情況下高時空分辨率的降水數據，對洪澇災害及滑坡監(jiān)測等具有極其重要的意義[2-4]。

綜合氣象站-雷達-衛(wèi)星遙感-地表輔助參量等多種資料，多源數據融合與數據同化技術已成為獲取高精度、高質量、時空連續(xù)的多時空尺度降水融合網格數據產品的有效手段[1,5]。根據地面觀測-雷達、地面觀測-衛(wèi)星遙感、雷達-衛(wèi)星遙感、地面觀測-雷達-衛(wèi)星遙感等不同源的降水數據組合，相關學者開發(fā)了包括地理加權回歸方法(Geographically Weighted Regression，GWR)[6-7]、回歸克里金插值方法[8-9]、貝葉斯模式平均方法(Bayesian Model Averaging，BMA)[10-12]、最 優(yōu) 插值[13-14]、多尺度模式融合方法[15]、隨機森林與XGBoost 機器學習方法[16-19]等不同的多源降水數據融合模型，并在子午河盆地、淮河流域、揚子江流域、華南地區(qū)等區(qū)域進行應用分析。為獲取更高時間分辨率的降水數據產品，相關學者利用“概率密度函數+貝葉斯模型平均(BMA)+最優(yōu)插值”組合方法有效地融合地面觀測降水-雷達定量降水估測產品-衛(wèi)星遙感降水產品，研制了高精度的1 km 分辨率的逐時降水融合產品[1,20-21]，并對這些逐時降水融合產品在四川、黃土高原、長江流域、華南地區(qū)等區(qū)域的適用性進行對比評估[22-24]。結果表明，在華南區(qū)域，融合降水產品與實際降水量的偏差較小，對實際降水量的再現能力強。

盡管現有的逐時降水數據融合研究取得一定進展，針對龍舟水、臺風水等不同類型的暴雨過程，逐時多源降水融合方法的適用性仍需進一步研究。此外，現有的逐時降水融合方法較少考慮暴雨過程降水的時間相關性。受到熱帶與中緯度天氣系統(tǒng)的共同影響，華南地區(qū)成為我國暴雨最頻繁的地區(qū)之一。基于此，本文利用廣東省北部山區(qū)2018 年汛期強降水時段4 月23—28 日(北京時間，下同)龍舟水、5 月7—11 日龍舟水和9 月16—17 日臺風“山竹”3 次暴雨過程250 個氣象站逐時降水數據，以及雷達降水、衛(wèi)星遙感降水、地形與海岸線距離等輔助變量，在對不同類型的逐時降水數據進行質量控制處理基礎上，基于機器學習算法與地統(tǒng)計學理論，開展空間分辨率為1 km逐時氣象站點-雷達-衛(wèi)星遙感降水的融合試驗，充分考慮相鄰時刻降水的時間相關性，以期獲取復雜地形區(qū)域暴雨過程高精度的降水融合數據，并對逐時降水融合試驗結果進行分析評估。

2 數據及預處理

2.1 研究區(qū)概況

研究區(qū)主要位于粵北山區(qū)(圖1)。研究區(qū)海拔最低6 m，最高1 427 m，高海拔區(qū)域主要集中在東西兩邊，中間區(qū)域海拔整體相對較低，主要集中在清遠東南部的英德、清新，清城境內的北江河谷。研究區(qū)屬亞熱帶季風氣候，雨水資源豐富，平均年降水量1 631.4～2 149.3 mm，年平均降水日(日降水量≥0.1 mm/d)為160～173 天，是廣東省一個典型的“雨窩”[25]。

圖1 研究區(qū)域及氣象站點分布 a.廣東省邊界范圍；b.研究區(qū)域及氣象站點分布。

2.2 研究數據

本文以廣東省氣象局提供的逐時降水數據作為地面基準數據，250個氣象站點的空間分布如圖1a 所示。隨機選擇空間均勻分布的200 個氣象站點(占80%)作為訓練集站點，50 個氣象站點(占20%)為測試集站點并用來評價不同降水融合模型的精度(圖1b)。選取2018 年4 月23—27 日(過程I)、5 月7—10 日(過程II)、9 月16—17 日(過程III)3次暴雨過程的降水數據進行逐時降水數據融合試驗，相關數據介紹如表1所示。

表1 相關數據來源及介紹

雷達降水：由1 km 分辨率的逐6 min 雷達定量降水估測產品計算得到，具體的計算過程參考文獻[25]。

與海岸線距離、DEM：這兩個地表輔助參量的處理過程主要參考文獻[25]。

IMERG 和GSMaP 降 水：IMERG(Integrated multi-satellite retrievals for GPM) 和 GSMaP(Global satellite mapping of precipitation)是新一代全球降雨觀測計劃(Global precipitation measurement, GPM)廣泛應用的衛(wèi)星降水產品，具有覆蓋范圍更廣(擴展到全球)、時空分辨率更高的優(yōu)勢(IMERG：空間分辨率為0.1 °、時間分辨率為0.5 h；GSMaP：空間分辨率為0.1 °、時間分辨率為1 h)[26-28]。本文利用面到點克里金插值方法[27]分別對IMERG 和GSMaP 降水產品進行空間降尺度，得到空間分辨率為1 km 的逐時IMERG 和GSMaP降水產品。

3 方 法

3.1 考慮降水時間相關性的XGBoost 逐時降水數據融合方法

本文提出的考慮時間相關性的XGBoost逐時降水融合方法技術路線如圖2所示。

圖2 考慮時間相關性的多源逐時降水融合方法技術流程圖

(1) 利用普通克里金插值算法對前兩個時刻的地面氣象站逐時降水進行插值，得到空間分辨率為1 km的逐時網格降水數據(PreSiteOK)；前兩個時刻的降水插值結果(PreSiteOK)作為降水時間相關性引入到XGBoost模型中。

(2) 采用地統(tǒng)計學的面到點克里金插值(Area-to-point Kriging，ATPK)[27-29]分別對空間分辨率為0.1 °的逐時IMERG 和GSMaP 降水產品進行空間降尺度。面到點克里金插值方法是利用已知面對未知點進行插值估計的空間降尺度方法，其原理與普通克里金類似，即未知點值為其所在面以及附近面數據的線性加權求和[27]：

權重λx的求解如下：

式中，x為高空間分辨率(1 km)的待插值網格點，K為低空間分辨率(0.1 °)的網格數目，vi為低空間分辨率的網格，表示面與點協(xié)方差函數，表示面與面協(xié)方差函數為空間分辨率為0.1 °網格的IMERG 和GSMaP 降水產品,λx為權重，μx為拉格朗日乘子,為降尺度的空間分辨率為1 km的IMERG和GSMaP降水產品。

(3) XGBoost 是一種基于從GBDT(梯度提升樹)算法改進和擴展而來的經過優(yōu)化的集成學習算法[30]，本文擬采用XGBoost 算法建立當前時刻氣象站-前兩個時刻網格插值-雷達-衛(wèi)星遙感多源降水預測模型：

通過步驟(3)獲取每時刻的XGBoost降水融合模型以及對應的逐時降水殘差。將空間分辨率為1 km 的逐時雷達降水數據、前兩個時刻網格降水數據、空間降尺度的衛(wèi)星遙感降水產品和地表輔助變量作為輸入，利用構建的XGBoost 逐時降水融合模型預測1 km分辨率的逐時降水。

(4) 對XGBoost 降水融合模型的逐時降水殘差進行普通克里金插值，得到1 km 分辨率的逐時降水殘差，并與模型預測結果進行相加得到1 km分辨率的高質量逐時降水融合結果。

3.2 融合試驗設計

為了進一步驗證不同逐時降水融合模型的性能，設計了“考慮降水時間相關性的XGBoost逐時降水數據融合(本文提出方法，記為方案I)”、“不考慮降水時間相關性的XGBoost逐時降水數據融合(XGBoost 方法，記為方案II)”、“不考慮降水時間相關性的隨機森林逐時降水數據融合(RF方法，記為方案III)”三組逐時降水融合對比試驗方案。具體的試驗設計為：方案I，利用普通克里金插值算法對前兩個時刻的地面氣象站逐時降水進行插值處理，并作為輔助變量，采用XGBoost算法進行逐時氣象站點降水-前兩個時刻逐時網格降水-雷達降水-衛(wèi)星遙感降水融合分析；方案II，不考慮前兩個時刻逐時網格降水，采用XGBoost 算法進行逐時氣象站點降水-雷達降水-衛(wèi)星遙感降水融合分析；方案III，不考慮前兩個時刻逐時網格降水，采用RF算法進行逐時氣象站點降水-雷達降水-衛(wèi)星遙感降水融合分析。

3.3 評價方法

結合檢驗站點的降水觀測數據，采用主要的評價指標對逐時降水融合結果進行定量評價，包括均方根誤差(RMSE)、平均絕對誤差(MAE)和決定系數(R2)，其計算公式如下：

式中，n表示用于檢驗的站點數，Ok表示第k個氣象站點降水觀測數據，表示第k個氣象站點對應的逐時降水融合結果，Oˉ表示氣象站點降水觀測數據的平均值。

4 結果與討論

結合研究區(qū)的50個地面站點逐小時降水觀測數據，分別利用決定系數(R2)、均方根誤差(RMSE)和平均絕對誤差(MAE)3 個統(tǒng)計指標評價2018 年4 月23—28 日(龍舟水)、5 月7—11 日(龍舟水)和9月16—17 日(臺風“山竹”)3 次暴雨過程的小時尺度多源降水數據融合結果的精度。

圖3 對應2018 年3 次暴雨過程的方案I、方案II和方案III逐時降水融合結果。結果顯示，3次暴雨過程中，方案I融合結果比方案II和方案III融合結果的精度高，R2均大于0.609。9 月暴雨過程方案I 融合結果的RMSE 和MAE 最小，分別為2.198和1.065 mm；方案I 降水融合的R2達到了0.78，高于方案II和方案III降水融合的R2。這可能是因為方案I不僅考慮了當前時刻地面降水與雷達降水、衛(wèi)星遙感降水、地表輔助變量間的非線性關系，還將前兩個時刻地面降水作為自變量引入到模型中，考慮了降水數據的時間相關性。此外，圖3 也顯示了3 次暴雨過程的方案II 逐時降水融合都優(yōu)于方案III 逐時降水融合，3 次暴雨過程的最大逐時降水量逐漸增加，暴雨過程II 的最大逐時降水量超過50 mm；暴雨過程I的逐時降水融合結果比過程II 和過程III 的逐時降水融合結果的精度低。結果表明，在降水量相對較大的情況下，三種方案都能獲取較高精度的逐時降水融合數據；然而，在降水量相對較小的情況下，三種方案的降水融合結果精度相對較低。這可能是因為在降水量相對較小的情況下，衛(wèi)星遙感降水產品的估測存在不足[31]，逐時降水融合模型引入衛(wèi)星遙感降水產品導致降水融合結果存在一定的誤差。

圖3 2018年過程I(a)、過程II(b)和過程III(c)三種方案的逐時降水融合與站點降水觀測比較

圖4 顯示了3 次暴雨過程不同融合方法的逐時降水融合的R2、RMSE 和MAE 箱圖。對比3 次暴雨過程的逐時降水融合，方案I 的3 次暴雨過程降水融合的中位數R2分別為0.415、0.466和0.734，都高于其他兩種方法降水融合的中位數R2；除了5月暴雨過程方案I 降水融合的中位數RMSE 略大于方案II 降水融合的中位數RMSE，方案I 的暴雨過程I 和III 降水融合中位數RMSE 都小于其他兩種方法降水融合的中位數RMSE，暴雨過程I和III降水融合的中位數RMSE 分別為0.418 mm 和1.735 mm；方案I 的3 次暴雨過程降水融合的中位數MAE 都高于其他兩種方法降水融合的中位MAE。結果表明，方案I降水融合結果優(yōu)于方案II和方案III降水融合結果。這進一步表明在逐時降水融合過程中引入多個相鄰時刻的降水信息能有效提高降水融合模型的精度。

圖4 不同融合方法的逐時降水融合的R2(a)、RMSE(b)和MAE(c)箱圖

比較方案II 和方案III 降水融合發(fā)現，暴雨過程I 和過程II 的方案II 降水融合的中位數R2小于方案III 降水融合的中位數R2，但3 次暴雨過程方案II 降水融合的R2均值都大于方案III 降水融合的R2均值；3次暴雨過程的方案II降水融合的中位數RMSE 和MAE 整體上小于方案III 降水融合的中位數RMSE 和MAE。這表明，方案II 降水融合的精度整體上優(yōu)于方案III降水融合。從圖4 也可看出，方案I、方案II 和方案III 三種融合方案在過程III 都能獲取較高精度的逐時降水融合結果，其中位數R2分別達到了0.734、0.722和0.643，遠高于暴雨過程I和過程II逐時降水融合的中位數R2。

表2 和圖5～7 分別顯示了3 次暴雨過程方案I、方案II和方案III的R2表現最好的時刻占總暴雨過程比例情況、逐時降水融合的R2、RMSE和MAE時間序列以及對應時間50個地面站點逐時降水觀測的平均值。從表2 可看出，3 次暴雨過程方案I獲得更高精度的逐時降水融合，優(yōu)于方案II 和方案III 降水融合，R2占比全部大于53.82%，遠高于其他兩種降水融合方案。

表2 不同降水融合方案R2表現最好的時刻占總暴雨過程的比例情況

圖5 暴雨過程I逐時降水融合R2(a)、RMSE(b)、MAE(c)精度評價

從圖5 看出，暴雨過程I 的方案I 降水融合整體上優(yōu)于方案II 和方案III 降水融合(表2)。但也出現方案I 降水融合精度比方案II 和方案III 降水融合精度低的情況，比如暴雨過程I 的23 日23 時方案I 降水融合的R2遠小于方案II 和方案III 降水融合的R2。這可能是因為暴雨過程I的23 日21—22 時降水量比較大，而到了23 時降水量顯著減少，出現有的區(qū)域可能沒有降水的情況，方案I 引入21—22時兩個時刻相對較大降水信息反而增加降水融合誤差。從圖5也可看出，盡管方案II降水融合精度整體上稍微優(yōu)于方案III 降水融合，但也有不少時刻出現方案II降水融合精度差于方案III降水融合。

從圖6看出，暴雨過程II的方案I降水融合R2、RMSE 和MAE 整體上優(yōu)于方案II 和方案III 降水融合，特別是在持續(xù)降水量比較大的暴雨過程II的9 日02 時—10 日12 時這段時間內，方案I 降水融合精度明顯優(yōu)于方案II和方案III降水融合。在這段時間內，部分方案II 降水融合的R2低于方案III 降水融合的R2，但這兩種方法的RMSE 和MAE相差不大。

圖6 暴雨過程II逐時降水融合R2(a)、RMSE(b)、MAE(c)精度評價

圖7 顯示了暴雨過程III 的方案I、方案II 和方案III 降水融合的R2、RMSE 和MAE 時間序列，盡管出現方案II和方案III降水融合的R2略優(yōu)于方案I降水融合的R2，如暴雨過程III的16日10時；但整體上，方案I 降水融合的R2、RMSE 和MAE 明顯優(yōu)于方案II 和方案III 降水融合的R2(表2)、RMSE 和MAE。這表明在這次暴雨過程中，考慮前兩個時刻降水信息能更好地改善逐時降水融合模型精度。此外，圖7 也顯示了，方案II 降水融合精度整體上優(yōu)于方案III 降水融合精度，因為方案II 降水融合具有更高的R2和更低的RMSE和MAE。

圖7 暴雨過程III的逐時降水融合R2(a)、RMSE(b)、MAE(c)精度評價

圖8 顯示了2018 年3 次暴雨過程不同融合方法預測的累計降水空間分布，在3次暴雨過程方案I、方案II與方案III方法的降水融合結果具有相似的空間特征。

圖8 不同融合方法的累計降水融合結果空間分布圖 a.暴雨過程I；b.暴雨過程II；c.暴雨過程III。

暴雨過程I的累積最大降水量不超過110 mm，暴雨過程II 的累積最大降水量增加到180 mm，到了暴雨過程III，累積最大降水量超過了210 mm。在累積降水量比較高的區(qū)域，方案II 降水融合結果的高值范圍整體上比方案I 和方案III 降水融合結果的高值范圍大；在累積降水量比較低的區(qū)域，方案II 降水融合結果的低值范圍整體上比方案I和方案III 降水融合結果的低值范圍大。此外，方案III 在極大與極小降水融合的性能表現不太好，在降水量相對較高和較低的區(qū)域，方案III 降水融合不能較好地刻畫降水分布的空間細節(jié)。總體上來看，方案I 降水融合結果獲得較滿意的結果，具有更多的空間分布特征。

總的來說，相比方案II 和方案III 逐時降水融合模型，方案I 獲得較高精度的逐時降水融合結果，逐時降水融合具有更豐富的空間分布特征。然而，方案I 逐時降水融合結果仍存在一定的誤差，這可能與引入的空間分辨率為0.1 °的IMERG和GsMaP衛(wèi)星遙感降水產品有關。粗分辨率的衛(wèi)星遙感降水產品在高值和低值的降水估測仍存在一定的誤差[31-33]。未來的研究應該引入更多時間序列暴雨過程的逐時降水信息，通過長短期記憶人工神經網絡(LSTM)挖掘暴雨過程降水的時間特征信息[34]，并以此優(yōu)化逐時降水融合模型，進一步提高逐時降水融合精度。此外，結合遙感數據獲取的云量、水汽等輔助參數[35]，在充分挖掘降水量與云量、水汽之間時空特征信息基礎上，引入到逐時降水數據融合模型中提高多源降水數據融合模型的精度。

5 結 論

本文以廣東省北部山區(qū)2018年汛期強降水時段4 月23—28 日龍舟水、5 月7—11 日龍舟水和9月16—17 日臺風“山竹”3 次典型暴雨過程的逐時降水為研究對象，建立基于XGBoost 與克里金插值算法的地面觀測-雷達-衛(wèi)星遙感多源逐時降水融合模型，充分考慮相鄰時刻降水的時間相關性，開展了空間分辨率為1 km 的逐時降水融合試驗，并與不考慮降水時間相關性的XGBoost和隨機森林(RF)算法的逐時降水融合模型進行對比，得到如下結論。

(1) 考慮降水時間相關性的方案I逐時降水融合模型融合結果精度整體上高于方案II和方案III逐時降水融合結果，與氣象站逐時降水數據最接近，因為方案I通過引入前兩個時刻逐時網格降水數據充分考慮了相鄰時刻降水的時間相關性，使得降水融合結果能夠更好地刻畫空間細節(jié)。在出現降水量極大和極小值的區(qū)域，不考慮降水時間相關性的方案II和方案III逐時降水融合精度有所降低。

(2) 與方案II 和方案III 逐時降水融合結果相比，方案I逐時降水融合結果在不同暴雨過程的準確性均有明顯改進，3次暴雨過程的RMSE分別降低了6.0%和9.7%、6.3%和9.5%、10.5%和30.0%。

(3) 方案II 的精度整體上優(yōu)于方案III，3 次暴雨過程的RMSE 分別降低了3.9%、3.4%和21.8%，表明XGBoost算法在刻畫逐時氣象站降水與雷達降水、衛(wèi)星遙感降水、地表輔助參量間的非線性關系能力方面比RF算法更有優(yōu)勢。