杜 妍

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

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

1.3 降雨侵蝕力計(jì)算方法

研究表明，水土流失領(lǐng)域著名專家Xie等[6]提出的降雨侵蝕力模型在我國(guó)具有良好應(yīng)用性，故而利用該模型計(jì)算站點(diǎn)尺度降雨侵蝕力，公式如下：

(1)

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

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

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

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

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

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

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

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

(4)基于訓(xùn)練好的XGBoost模型對(duì)柵格變量集進(jìn)行空間預(yù)測(cè)，生成廣東省1km空間分辨率的降雨侵蝕力因子?xùn)鸥衩妗?/p>

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

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

2.1 廣東省站點(diǎn)觀測(cè)降雨侵蝕力統(tǒng)計(jì)特征

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

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

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

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

2.2 廣東省降雨侵蝕力與GPM降水之間關(guān)系特征

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

圖4 廣東省GPM降水量與站點(diǎn)降雨侵蝕力之間相關(guān)性散點(diǎn)圖

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

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

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

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

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

圖6 廣東省降雨侵蝕力提取精度散點(diǎn)圖

3 結(jié)論

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

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

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