基于GPM數(shù)據(jù)和XGBoost算法的降雨侵蝕研究

杜 妍

(韶關市防洪管理中心，廣東 韶關 512026)

土壤侵蝕是自然和人類活動對土地資源的破壞過程，導致了生態(tài)環(huán)境退化和農(nóng)業(yè)生產(chǎn)受損等問題。持續(xù)性水土保持作為一項重要生態(tài)治理工程，對維系土地生態(tài)健康、水資源涵養(yǎng)、植被多樣性有著不可忽視的作用[1-3]。在土壤侵蝕防治實踐中，需要提取一系列生態(tài)因子，而降雨侵蝕力因子是其中關鍵參數(shù)之一，其表征了降雨能量對地表土體產(chǎn)生的擊濺、剝離、運移能力。傳統(tǒng)研究中對區(qū)域尺度降雨侵蝕力因子提取依賴于稀疏雨量站資料，并結(jié)合空間統(tǒng)計方法進行分析處理，其顯著缺點是無法準確刻畫降雨侵蝕力柵格面，尤其在站點稀疏地區(qū)的估算精度存在較大不確定性[4-5]。近年來隨著衛(wèi)星遙感和機器學習技術(shù)發(fā)展，通過衛(wèi)星降水產(chǎn)品和非線性擬合方法，提出了一些快速準確提取降雨侵蝕力因子的方法，例如陳君等[3]利用GBDT和GPM數(shù)據(jù)生成了華南地區(qū)降雨侵蝕力產(chǎn)品，梁宇靖等[2]利用GLAS數(shù)據(jù)得到了中國降雨侵蝕力信息?；谇叭搜芯炕A上，基于GPM衛(wèi)星降水數(shù)據(jù)和XGBoost算法，探索了一種新的提取降雨侵蝕力因子的方法，并利用地面雨量站數(shù)據(jù)驗證了該方法的有效性。

1 材料與方法

1.1 研究區(qū)概況

廣東地處我國大陸南端，南海北岸、珠江流域中下游，總面積17.8萬km2，海岸線長3368km。區(qū)域?qū)賻X南丘陵、珠江三角洲平原和沿海平原地形，總體地貌崎嶇破碎海拔為0～1902m。受西太平洋季風、印度洋季風、赤道低壓和副熱帶高壓影響，區(qū)域形成南亞熱帶濕潤季風氣候，具有全年溫和、降水充沛特點，年平均氣溫18～24℃、降水量在1500～2000mm、無霜期約270～320d、積溫5800～7500℃、日照時數(shù)1500～2100h。其中降雨量集中于4—9月，其他月份降水較少，常有臺風、暴雨、干旱、洪水、雷電和霜凍等氣象災害影響。由于季風性降水集中性強、強度大，在降水月易形成強烈降雨侵蝕破壞力，對地表產(chǎn)生劇烈沖刷、剝離，因此該地水土流失現(xiàn)象極為普遍。如圖1所示。

1.2 數(shù)據(jù)來源

選擇準確而具有代表性的氣象環(huán)境數(shù)據(jù)是提取大尺度降雨侵蝕力分布信息的重要基礎。文中使用的數(shù)據(jù)主要有：①地面雨量站數(shù)據(jù)，從氣象數(shù)據(jù)共享服務網(wǎng)(https：//urs.earthdata.nasa.gov/)收集到研究區(qū)2020年85個氣象站點日降水資料；②GPM衛(wèi)星降水數(shù)據(jù)，從google earth engine開源平臺獲取“NASA/GPM_L3/IMERG_V06”產(chǎn)品，生成研究區(qū)GPM降水數(shù)據(jù)；Aster DEM數(shù)據(jù)，從地理空間數(shù)據(jù)云下載研究區(qū)DEM產(chǎn)品，用以提取高程、柵格點的經(jīng)度、緯度等信息。

1.3 降雨侵蝕力計算方法

研究表明，水土流失領域著名專家Xie等[6]提出的降雨侵蝕力模型在我國具有良好應用性，故而利用該模型計算站點尺度降雨侵蝕力，公式如下：

(1)

式中，Rday、Pd—日降雨侵蝕力、日侵蝕性降雨量；α—物候參數(shù)，規(guī)定5—9月該參數(shù)取值為0.3937，10—4月取0.3101。

1.4 基于XGBoost算法空間尺度降雨侵蝕力提取流程

XGBoost(extreme gradient boosting)是一種集成學習算法，其在梯度提升基礎上，通過迭代地添加弱學習器來提高模型性能[7-8]。即每次迭代中，XGBoost在上一次迭代的結(jié)果基礎上擬合一個新的決策樹，并通過梯度提升對之前預測結(jié)果的殘差進行擬合，如圖2所示。通過這種方式，XGBoost能夠逐步提升模型的泛化能力，因此具有高效性、精度高、靈活性特點，可以適應不同應用場景。該算法具體數(shù)學過程原理，詳見Chen等[7]人的研究。

圖2 XGBoost算法結(jié)構(gòu)圖

借鑒陳君[3]和梁宇靖[2]等人相關經(jīng)驗，使用XGBoost算法提取廣東省降雨侵蝕力因子簡要流程如下：

(1)利用公式(1)計算并合成全部85站點2020年降雨侵蝕力。

(2)對原GPM衛(wèi)星數(shù)據(jù)、DEM影像等進行重投影、去噪、重采樣、空間統(tǒng)計等預處理，生成1km空間分辨率的環(huán)境變量柵格數(shù)據(jù)集；然后利用全部站點點位提取柵格變量集里各變量數(shù)值信息，進而構(gòu)建樣本集；隨機選擇其中70%的樣本為建模集(n=59)、30%為驗證集(n=25)。

(3)利用訓練集結(jié)合Python3.9程序設計XGBoost回歸算法，因變量為站點測量侵蝕力，自變量為GPM降水量、海拔、經(jīng)度和緯度。為確保模型穩(wěn)健性，利用Grid方法對XGBoost算法進行優(yōu)化。

(4)基于訓練好的XGBoost模型對柵格變量集進行空間預測，生成廣東省1km空間分辨率的降雨侵蝕力因子柵格面。

(5)利用未參與模型訓練的氣象站點空間點位進行獨立驗證，即XGBoost算法預測得到的侵蝕力柵格面對應位置處侵蝕力值y′為參考值，以地面觀測站計算的侵蝕力y為真值，計算決定系數(shù)(R2)、平均絕對誤差(MAE)和均方根誤差(RMSE)，對區(qū)域降雨侵蝕力提取精度進行量化評估。

2 廣東省降雨侵蝕力提取結(jié)果

2.1 廣東省站點觀測降雨侵蝕力統(tǒng)計特征

廣東省2020年降水量最低值為1458.26mm，出現(xiàn)在德清站，最高值為增城站的2411.42mm，最大、最小相差953.16mm；進一步統(tǒng)計表明全部站點降水量平均值為1893.45mm，在全國屬豐水區(qū)，其空間離差系數(shù)為16.52%，表明全境降水量存在明顯不均性。見表1。另利用Xie模型得到的區(qū)域侵蝕力值介于5080.23～17958.04[MJ·mm/(hm2·h·a)]之間，全省平均值為9574.49[MJ·mm/(hm2·h·a)]，其侵蝕強度局全國較高水平。經(jīng)單樣本Kolmogorov-Smirnov檢驗發(fā)現(xiàn)，上述序列變量漸進顯著性sig值均<0.05，說明其不符合正態(tài)分布假設，在后續(xù)分析中須進行對數(shù)變換處理進而減少數(shù)據(jù)噪聲[9]。

表1 廣東省85個氣象站點降水量和降雨侵蝕力統(tǒng)計特征

利用GIS空間可視化功能得到研究區(qū)2020年降雨量分布柵格面圖3。研究區(qū)原GPM降水量數(shù)據(jù)的粗集產(chǎn)品揭示了降水量分布格局，可見該年度降水量中心位于珠三角及江門南部沿海地區(qū)，GPM數(shù)值可達2500mm；降水量低值出現(xiàn)在粵東和粵西丘陵地區(qū)，像素值在1700mm以下；總體呈現(xiàn)自珠三角核心區(qū)向周邊減少的分布特征。

圖3 廣東省GPM降水量空間分布圖

2.2 廣東省降雨侵蝕力與GPM降水之間關系特征

GPM降水量產(chǎn)品作為研究區(qū)降雨侵蝕力因子提取的重要協(xié)變量之一，GPM降水數(shù)值大小及其分布特征對區(qū)域尺度侵蝕力特征有一定影響。經(jīng)相關性分析發(fā)現(xiàn)圖4，二者擬合關系形式為y=7.4689x-4737.6，其R2=0.48，P<0.01，說明衛(wèi)星降水產(chǎn)品與近地表降雨侵蝕力之間存在密切關系，因此可作為區(qū)域降雨侵蝕力空間提取的有效協(xié)變量之一。

圖4 廣東省GPM降水量與站點降雨侵蝕力之間相關性散點圖

2.3 廣東省降雨侵蝕力提取結(jié)果分析

圖5a直觀呈現(xiàn)了廣東省2020年降雨侵蝕力空間細致分布。與圖5b對比可知，區(qū)域降雨侵蝕力與降水量具有高度一致性，這是由于前者是后者的函數(shù)。區(qū)域降雨侵蝕力分布范圍介于5584～15318[MJ·mm/(hm2·h·a)]之間，空間平均值為10864[MJ·mm/(hm2·h·a)]，離差系數(shù)為36.79%，意味著該省域降雨侵蝕力空間分異性顯著。其中降雨侵蝕高值區(qū)聚集于珠三角中北部邊緣、陽江和江門南部一帶，局部可達12000[MJ·mm/(hm2·h·a)]以上；結(jié)合相關觀測資料可知，2020年濕潤季風在該局地停留時間長、降雨強度大，由于這些地區(qū)處于城建核心區(qū)和城鄉(xiāng)結(jié)合帶，強烈降雨侵蝕可能誘發(fā)崩崗、滑坡、泥石流等災害。侵蝕力低值區(qū)出現(xiàn)在粵東、粵西和粵北邊緣地區(qū)，在10000[MJ·mm/(hm2·h·a)]以下，主要由于這些地區(qū)當前降水量較小而強度低；然而這些地區(qū)長期處于低影響開發(fā)、良好植被覆蓋狀態(tài)，因此該地產(chǎn)生的侵蝕風險可忽略不計。

圖5 廣東省降雨侵蝕力和降水量空間分布圖

2.4 降雨侵蝕力提取精度分析

研究區(qū)25個獨立驗證樣本點的XGBoost模型預測降雨侵蝕力與地面站點侵蝕力之間擬合散點圖，如圖6所示。分析發(fā)現(xiàn)二者散點分布于1∶1直線兩側(cè)，其擬合關系形式為y=0.6304x+3894.4，經(jīng)檢驗得到其P值<0.01，表明二者關系具有統(tǒng)計顯著水平。利用公式(5-7)計算得到其R2為0.36，表明侵蝕力真值與預測值具有良好一致性；且MAE和RMSE依次為1451.6、2880.9[MJ·mm/(hm2·h·a)]，意味著該降雨侵蝕力柵格面具有較低不確定性，因此圖4所示柵格面可用作區(qū)域水土流失評價或侵蝕模數(shù)計算的替代品。需指出的是，盡管該提前精度在可接受范圍內(nèi)，但由于區(qū)域內(nèi)侵蝕范圍大，到其RMSE也較大，因此該侵蝕力柵格面仍有改進潛力。

圖6 廣東省降雨侵蝕力提取精度散點圖

3 結(jié)論

(1)GPM降水與降雨侵蝕力之間存在密切線性關系，可成為空間尺度降雨侵蝕力研究新的有效數(shù)據(jù)源。

(2)XGBoost模型基于其非線性擬合優(yōu)勢，生成了精度可靠的侵蝕力柵格面，并具有精細水平表現(xiàn)力。

(3)該改進的方案可為廣東省水土流失評價、生態(tài)建模等貢獻新思路。本研究仍存在不足，例如GPM數(shù)據(jù)本身仍為粗糙產(chǎn)品(空間分辨率為10km)，利用簡單重采樣方法將其轉(zhuǎn)換為1km分辨率可能產(chǎn)生尺度轉(zhuǎn)換誤差。因此，后續(xù)研究應先將GPM數(shù)據(jù)進行降尺度處理，明確其尺度轉(zhuǎn)換誤差后用以輔助生成精細的降雨侵蝕力分布式信息。