基于RUSLE 和GIS 的博羅縣土壤侵蝕特征分析

陳燕奎，張 豫，朱猷亮

（1.嘉應學院 地理科學與旅游學院，廣東 梅州 514015;2.深圳航天宏圖信息技術有限公司,廣東 深圳 518000）

為了全面認識土壤侵蝕的過程，理解全球變化與區域土壤侵蝕關系[1]，Wischmeier 和Smish 分析了全美各地徑流泥沙觀測資料，提出通用土壤流失方程（Universal Soil Loss Equation,USLE）[2].1985 年美國科學家對USLE 進行修正，研發出修正通用土壤流失方程（Revised Universal Soil Loss Equation,RUSLE）[3].由文獻可知，學者利用RUSLE 方程研究水土流失土壤侵蝕量估算是有意義的，主要圍繞RUSLE 方程的研究進程主要體現在側重土壤侵蝕量估算[4-6]、方程因子改進與增加方程參數估算[7]與GIS、RS 等技術結合做土壤侵蝕的時空分析[8-9]以及在全球不同地區、不同尺度、不同時序等方面深度剖析土壤侵蝕方程的應用效益[10-12]；區域研究多以省、市層面，而南方山區縣級層面的研究還比較少.本文以廣東省博羅縣域為例，引入基于修正的通用土壤流失方程（RUSLE），綜合運用GIS 和RS 技術，試探索定量評價博羅縣流域土壤侵蝕現狀及空間分布，并分別估算不同坡度及土地利用類型下的土壤侵蝕量，為本區域的生態環境保護和政府決策提供參考依據.

1 區域概況與數據源

1.1 區域概況

博羅縣位于廣東省珠江三角洲東北部，地理坐標為北緯23°03′50″～23°43′20″，東經113°49′50″～114°45′50″，全縣總面積約2 858 平方千米，屬東江水系，地貌類型為北部山地丘陵、間有山谷平原，中部丘陵臺地，南部沿東江自東向西的沖積平原等3 個地帶，海拔約為208～368 m，屬亞熱帶季風氣候，年平均氣溫22.1℃，多年平均降水量約為1 918 mm，日照時數約為1 871.5 h.全域地形起伏比較大，海拔-10～1 239 m，土壤則以粉砂土為主，粉質黏粘土為輔，另外森林覆蓋率較高，土地利用類型多樣，植被物種豐富.

1.2 數據源

所用數據來源如表1 所示，在ArcMap10.7 和PIE-Basic6.0 等軟件功能的支持下，建立博羅縣域的土壤流失方程各類因子空間數據庫，并將所有因子圖層設立統一坐標系（CGCS2000 坐標系）和劃分評價單元：30 m×30 m 的柵格單元.此外，土地利用數據類型按照用地性質初步劃分為耕地、有林地、水域、建設用地、草地、裸地以及未利用地7 大類等.

表1 博羅縣土壤侵蝕所需基礎數據

2 數據處理與方法

2.1 土壤侵蝕方程

基于國內外現有成果和博羅縣土壤侵蝕發生的特點，建立了坡面土壤侵蝕模型，即博羅土壤流失方程BLSLE（BoLuo Soil Loss Equation），其基本形式為

式中：A為土壤流失量，t·(hm-2·a-1)；R為降雨侵蝕力因子，MJ·mm·(hm-2·a-1)，1MJ=106J；K為土壤可蝕性因子，t·hm2·h/(hm2·MJ·mm)；L為坡長因子，無量綱；S為坡度因子，無量綱；C為地表植被覆蓋與管理因子，無量綱；P為水土保持措施因子，無量綱.

根據式（1），博羅縣土壤流失方程構建關鍵在于6 個區域性參數取值的確定.

2.2 RUSLE 模型中各因子值的確定

2.2.1 降雨侵蝕力因子（R）

在USLE 方程中，降雨侵蝕力因子（R）被定義為降雨動能和最大30 分鐘降雨強度的乘積.引用文獻[13]、[14]基于日降雨量提出計算廣東地區降雨侵蝕力的簡易算法模型：

式中，iR表示第i個半月時段的侵蝕力值（MJ·mm·hm-2·-h·-a）；k表示該半月時段內的天數；pj表示半月時段內第j天的日雨量，要求日雨量≥12 mm，否則以0 計算，12 mm 與侵蝕性降雨標準對應；α和β為模型待定參數.模型參數α和β值的計算方程為

式中，Pd12表示日降雨量≥12 mm 的日平均雨量（mm），Pd12表示日降雨量≥12 mm 的年平均雨量.

文獻[13]指出，在年降雨量較小的地區，模型的決定系數較小、估算多年平均侵蝕力的相對誤差較大，而在年降雨量比較豐富的地區相反，模型的決定系數較大、估算多年平均侵蝕力的相對誤差較小.而本研究區年降雨量比較豐富，20 年日降雨量≧12 mm 的年平均雨量為1 615 mm，所以該公式在本地區可行.收集到博羅縣及其附近增城、龍門、河源、惠陽、惠東共6 個站點2000～2021 年的日降雨量數據計算博羅縣降雨侵蝕力的空間分布，具體地理位置見圖1.小部分缺失數據，通過最鄰近法空間插值方法，保持數據序列的完整和連續性.

圖1 氣象站點位置

首先，根據式（3）、（4）和收集到的各氣象站點日降雨數據，計算各氣象站點的模型參數α和β值（見表2），再根據式（2）計算各氣象站點半月時段的降雨侵蝕力值并累加得出各年降雨侵蝕力值，取各站點2000～2021 年平均降雨侵蝕力值為該站點的降雨侵蝕力值，應用樣條函數法插值生成研究區降雨侵蝕力分布如圖2 所示.

圖2 博羅縣降雨侵蝕力分布

表2 各站點模型參數與年平均降雨侵蝕力值

2.2.2 土壤可蝕性因子（K）

結合博羅縣的具體調查確定土壤機械顆粒組成，選擇Williams 等提出的EPIC 模型，計算土壤可蝕性因子的值.按USDA 土壤類型分類標準對1∶100 萬土壤類型進行分類，將博羅縣區域USDA 土壤分類為粉質粘土、粉質粘壤土、粉土、壤土和砂土，各種土壤在區域所占比例如圖3 所示.

圖3 博羅縣區域USDA土壤分類

根據博羅縣USDA 土壤類型分類圖，以及各種土壤類型中粘粒、粉粒、砂粒和有機碳含量，運用侵蝕-生產力影響模型EPIC（Erosion Productivity Impact Calculation）中土壤可蝕性因子估算公式（5）計算各種土壤類型的土壤可蝕性[15]，土壤可蝕性值計算結果如表3 和圖4 所示.

圖4 土壤可蝕性值分布

圖5 L·S 乘值的空間分布特征

表3 土壤組成及土壤可蝕性值

式中，KEPIC為EPIC 模型計算的土壤可蝕性因子（t·hm2·h·hm-2·MJ-1·mm-1），SN、SIL和CLA分別為美制土壤粒級分類標準中的砂粒、粉粒和黏粒含量，SN1=1-SN/100，Cα為有機碳含量（%）.

2.2.3 坡度因子（S）與坡長因子（L）

坡度和坡長均是土壤侵蝕量變化的加速因子，主要體現在坡度高低和坡面長度對土壤沙土侵蝕速度、搬運量和長度的影響.同等條件下，坡度越大，坡面的沖刷能力越強，而且搬運量將隨著波長（L）的增大而增大.對于波長和坡度研究，小尺度研究可根據實測值代入分析，對于大區域的估算，可從DEM 數據提取坡度和波長.修正土壤流失方程（RULSE）的坡長坡度因子計算方法為：

式中，S為坡度因子；θ為DEM 提取的坡度值；L為坡長因子；λ為DEM 提取的坡長值：22.13 m 為標準小區的坡長；α為坡度坡長因子.

該模式是學者在美國耕地坡度現狀基礎上建立的，因各地區的坡度坡長變化不用，該模式并非適合直接套用到其他研究地區，只能作為基礎模型范本，由此，實際使用該模式需要根據地方實際進行模型參數修正.

（1）坡度因子修正.美國耕地坡度普遍小于20%（11.3°），博羅縣大部分耕地分布于山地丘陵的坡地上，坡度值大于15°，坡地利用在博羅縣丘陵山地較為普遍，因此借鑒文獻[16]對坡度在9%～55%的陡坡侵蝕的研究對坡度因子進行修正，坡度小于14°沿用原有的坡度因子計算方式，其余坡度值，重新提出坡度因子計算方式：

因此將修正后的坡度因子進行合并，則得到

（2）波長因子修正.坡長因子在ArcMap10.7 軟件下進行修正計算，通過PIE-Basic6.0 等軟件DEM 提取坡度數值圖，可以根據坡度值推算出α坡度坡長因子，公式如下：

式中β為細溝侵蝕和面蝕的比值.

2.2.4 地表植被覆蓋與管理因子（C）

植被覆蓋率與土壤水土保持具有很強的相關性，一般而言植被覆蓋度越大，其植被葉、樹干對雨水沖刷地表起到降速與阻隔作用和植被根系對土壤起到固定作用就越大，從而使地表土壤在受到雨水作用時越不輕易被剝蝕.學者們對C值研究，定義了C值取值介于0～1 之間，植被覆蓋度（?）越高，其C取值越小，表示植被對土壤侵蝕的抑制作用越大.采用蔡崇法等[17]的分析方法，用于分析植被覆蓋度（f）與對應C值的關系.

（1）計算植被覆蓋度（f）.利用PIE-Basic6.0 軟件將2019-2021 年的1 月，4 月，7 月，10 月的TM-L8 數據集，提取涵蓋博羅縣范圍的影像數據，利用紅外與近紅外波段，運用ENVI 軟件的波段運算功能計算博羅縣植被歸一化指數（NDVI）年度NDVI 平均值：

式中，f為植被覆蓋度（%）；NDVI為歸一化植被指數，其中NDVImax和NDVImin為歸一化植被指數中的最大值和最小值.

（2）計算植被覆蓋與管理因子（C），公式如下：

根據文獻[17]，當植被覆蓋度為0 時，最易發生土壤侵蝕，C取值為1，覆蓋度大于78.3%時，發生土壤侵蝕概率小，C取值為0.在PIE-Basic6.0 軟件的支持下，利用柵格計算器功能，得出研究區植被覆蓋與管理因子空間分布圖，如圖6 所示.

圖7 水土保持措施因子空間分布

2.2.5 水土保持措施因子（P）

水土保持措施因子（P）主要以采取水土保持措施后的土壤侵蝕量與對應的未采取水土保持措施的情況下順坡種植時的土壤侵蝕量的比值.取值范圍0～1 之間，丘陵山區的隨著地表利用方式的不同，其P值也有差異，從而土壤侵蝕的作用也有不同，P值越小，表示水土保持能力越好，土壤越不易被侵蝕.參考文獻[18]、[19]對不同的土地利用類型的P進行賦值，以博羅縣2021 年的土地利用現狀數據為基礎（見表4），運用ArcMap 分析工具得到P 值柵格圖層.

表4 不同土地利用類型P值

3 結果與分析

3.1 研究區土壤侵蝕空間分布特征

將上述各因子的初始數據進行坐標匹配，CGC2000 坐標系，劃分柵格單元格大小30 m×30 m，利用ArcMap10.7 軟件的柵格計算功能將博羅縣RUSLE 模型各因子計算結果進行柵格疊加計算，得到博羅縣的土壤侵蝕模數分布圖（見圖8），博羅縣年平均土壤侵蝕模數為9721.24 t·km-2·a-1，年侵蝕量892.61×104t.

圖8 博羅縣土壤侵蝕模數空間分布

從空間分布特征分析，將博羅縣土壤侵蝕分布按照鎮級行政區劃進行統計分析，得到的各行政區土壤侵蝕情況統計如表5 所示，按照南方紅壤丘陵侵蝕等級劃分，石灣鎮及園洲鎮屬微度，土壤侵蝕較為嚴重的行政區排名前五的是羅浮山、橫河鎮、公莊鎮、柏塘鎮及泰美鎮，土壤侵蝕總體發生區域分布在西北部及東中部山地丘陵區，該區域海拔相對較高，山地多，降雨充沛，可以看出雨水沖刷和地形起伏對土壤侵蝕的影響比較顯著，為土壤侵蝕提供了有利條件.

表5 博羅縣各行政單元土壤侵蝕量

博羅縣屬于南方紅壤丘陵區，分為6 個侵蝕等級，其大部分地區年平均土壤侵蝕模數屬于輕微度水平，按照該標準分類，得到博羅縣各等級土壤侵蝕面積如表6 所示，各等級土壤侵蝕空間分布如圖9 所示.各侵蝕強度的面積大小依次是微度＞劇烈＞輕度＞強度＞極強度＞中度，分別占47.45%、15.36%、14.28%、10.81%、7.77%、4.34%.雖然博羅縣土壤侵蝕等級極強度以上比重只占23.13%，但土壤侵蝕量卻占到82.51%，其中劇烈等級的平均侵蝕模數高達46 759.54 t·km-2·a-1，遠高于當地容許土壤流失標準500 t·km-2·a-1，微弱侵蝕面積僅為47.45%，而劇烈侵蝕年土壤侵蝕量卻為微度侵蝕年土壤侵蝕量的63 倍.

圖9 博羅縣土壤侵蝕分級

表6 博羅縣土壤侵蝕分級

3.2 土壤侵蝕與坡度空間分布特征分析

地形是土壤侵蝕的主要影響因素，利用博羅縣DEM 數據生成坡度圖，按照水利部土壤侵蝕分類分級標準將坡度范圍劃分為6 個等級區，統計得出各等級下的土壤侵蝕分析如表7 所示.土壤侵蝕面積最大的為坡度0°～5°，占土壤侵蝕面積比例的51.12%，其土壤侵蝕量的比例僅占總量的10.24%，其平均土壤侵蝕模數為3 265.26 t/（km2·a），處于中度以下侵蝕級別，也是6 個坡度等級中侵蝕模數最低的.8°～15°和15°～25°兩個坡度等級是僅次于0°～5°土壤侵蝕面積，分別占到土壤侵蝕面積的19.13%和15.56%；兩者的土壤侵蝕量比例分別為25.11%和36.81%，為6 個坡度等級中侵蝕量比例最大的兩個值，另外25°以上兩個坡度等級，雖然其侵蝕面積較少，僅占土地總面積的2.82%，但其土壤侵蝕量卻占到侵蝕總量的19%，而且是6 個坡度等級中平均侵蝕模數最大的兩個，遠超過南方紅壤侵蝕最高等級的基數（15 000 t/（km2·a））.

表7 博羅縣不同坡度下土壤侵蝕特征分析

3.3 土壤侵蝕與土地利用空間分布特征

利用土地利用類型圖和土壤侵蝕強度等級圖進行疊加分析分析如表8 所示.各土地利用類型下的土壤侵蝕量來看，有林地范圍內的土壤侵蝕量占到總侵蝕量的45.08%，耕地和園地范圍的土壤侵蝕量占到總侵蝕量的25.16%和13.43%，另外雖然裸地范圍僅占土地總面積的0.45%，但其土壤侵蝕量卻占到總土壤侵蝕量的6.41%，屬于激烈侵蝕等級.

表8 博羅縣不同土地利用下土壤侵蝕特征分析

3.4 植被覆蓋度與土壤侵蝕的影響

利用植被覆蓋度（C）和土壤侵蝕空間分布圖層進行疊加分析，劃分為5 個等級：＜20%、20%～35%、35%～50%、50%～65%、＞65%.統計各覆蓋度等級下的土壤侵蝕情況如表9 所示.植被覆蓋度小于20%和大于65%的區域土壤侵蝕面積高達26.78%和24.57%，其他植被覆蓋度級別下的土壤侵蝕面積變化不大，植被覆蓋度小于20%的平均侵蝕模數最大，為30 045.57 t/（km2·a），雖然植被覆蓋度低的區域為城鎮村建設用地和水域地區，但其還包含裸地、稀疏林等，其雨水沖刷和地形坡度的影響下，使得其平均侵蝕模數較大.而植被覆蓋度高的地區，其平均侵蝕模數較小，均為2 240.67 t/（km2·a），侵蝕量也最小，僅為62.98 t/（a），雖然其地形起伏大，因其林木高而密，林木種質多樣，林下層級多樣，對雨水沖刷起到一定減緩作用.另外20%～65%三個覆蓋度等級，其土壤侵蝕變化不大，但其侵蝕量則占到土壤侵蝕總量的81.23%，這些地區多為山區周圍，地形起伏變化大，而且種植了大量的人工林、稀疏林木、次生林，特別是桉樹類，林分結構較為單一，林下植被較為匱乏，在雨季，特別是大雨、暴風雨等，雨水的沖刷力大，植被對土壤保持起到的作用較弱，其土壤侵蝕情況不容忽視.

表9 博羅縣不同被覆蓋度下土壤侵蝕分析

4 結論

本文利用修正的通用土壤流失方程RUSLE 模型，在GIS 技術支持下，定量估算了博羅縣土壤侵蝕強度，并結合遙感影像、DEM、土地利用等數據，通過空間分析和數理統計方法分析了博羅縣內土壤侵蝕的空間分布特征及影響因素，得出如下結論.

（1）2021 年，博羅縣土壤平均侵蝕模數為9 721.24 t·km-2·a-1，表明水土流失較為嚴重，年侵蝕量892.61×104t.6 個侵蝕強度等級的面積大小依次是微度＞劇烈＞輕度＞強度＞極強度＞中度，其中劇烈等級侵蝕面積雖然只占到15.36%，但土壤侵蝕量卻占到總侵蝕量的54.7%.

（2）石灣鎮及園洲鎮屬微度，而西北部及東中部山地丘陵區的羅浮山、橫河鎮、公莊鎮、柏塘鎮及泰美鎮等區域，該區域海拔相對較高，山地多，地形起伏大、降雨充沛，雨水沖刷大，土壤侵蝕較為嚴重.

（3）博羅縣土壤侵蝕量與坡度存在明顯的正相關，每增加一個等級，平均侵蝕模數呈增大趨勢，而且增大幅度較大；不同的土地利用類型會直接影響土壤侵蝕的量，有林地的面積為52.64%，其土壤侵蝕量比重占45.08%，而裸地雖然面積比重僅占0.45%，但其土壤侵蝕量卻占到總土壤侵蝕量的6.41%，屬于激烈侵蝕等級，因此裸土地是土壤侵蝕易發地區，需要進一步加強水土保持治理力度，做到治理與預防并重.

（4）植被覆蓋度小于20%和大于65%的區域土壤侵蝕面積高達26.78%和24.57%，其他植被覆蓋度級別下的土壤侵蝕面積變化不大；植被覆蓋度20%～65%之間的土壤侵蝕面積為48.64%，但其土壤侵蝕總量比重達81.23%，該區域多為山區丘陵，地形起伏大，分布著大量的人工林、次生林等林分結構單一、林下植被較為匱乏的經濟林木，雨水沖刷嚴重，不利于水土保持，該區域土壤侵蝕情況不容小覷.