陳燕奎，張 豫，朱猷亮

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

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

1 區(qū)域概況與數(shù)據(jù)源

1.1 區(qū)域概況

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

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

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

表1 博羅縣土壤侵蝕所需基礎數(shù)據(jù)

2 數(shù)據(jù)處理與方法

2.1 土壤侵蝕方程

基于國內(nèi)外現(xiàn)有成果和博羅縣土壤侵蝕發(fā)生的特點，建立了坡面土壤侵蝕模型，即博羅土壤流失方程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為水土保持措施因子，無量綱.

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

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

2.2.1 降雨侵蝕力因子（R）

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

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

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

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

圖1 氣象站點位置

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

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

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

2.2.2 土壤可蝕性因子（K）

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

圖3 博羅縣區(qū)域USDA土壤分類

根據(jù)博羅縣USDA 土壤類型分類圖，以及各種土壤類型中粘粒、粉粒、砂粒和有機碳含量，運用侵蝕-生產(chǎn)力影響模型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）

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

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

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

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

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

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

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

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 數(shù)據(jù)集，提取涵蓋博羅縣范圍的影像數(shù)據(jù)，利用紅外與近紅外波段，運用ENVI 軟件的波段運算功能計算博羅縣植被歸一化指數(shù)（NDVI）年度NDVI 平均值：

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

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

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

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

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

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

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

3 結果與分析

3.1 研究區(qū)土壤侵蝕空間分布特征

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

圖8 博羅縣土壤侵蝕模數(shù)空間分布

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

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

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

圖9 博羅縣土壤侵蝕分級

表6 博羅縣土壤侵蝕分級

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

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

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

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

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

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

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

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

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

4 結論

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

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

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

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

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