基于RUSLE模型的梅河口市土壤侵蝕動態分析

劉 盼， 任春穎， 王巖松， 劉建祥

(1.中國科學院 東北地理與農業生態研究所濕地生態與環境重點實驗室, 長春 吉林， 130102;2.中國科學院大學, 北京 100049; 3.松遼水利委員會松遼流域水土保持監測中心站, 長春 吉林， 130021)

東北黑土區是中國主要的糧食生產基地[1]，也是東北地區重要的生態屏障。長期以來，自然因素和不合理的人類活動導致東北黑土區出現了侵蝕溝發育嚴重[2]、土層變薄[3]等水土流失問題，進而導致農業生產能力下降[1]，阻礙了黑土區的農業生產和經濟發展[1,3]。土壤侵蝕定量分析與評價是水土流失治理的前提，而掌握區域土壤侵蝕動態變化對評價該地區的水土保持治理效果、指導水土保持措施優化配置、水土資源的保護和可持續利用具有重要意義[4]。

土壤侵蝕定量評價方法有侵蝕模型法、數字高程法和核示蹤法[5]。常用的侵蝕模型包括基于經驗統計和基于物理過程的侵蝕模型兩大類[5]。基于物理過程的侵蝕模型[6-7]通常需要輸入大量參數，參數不易獲得且輸入參數的可變性會導致結果的不確定性增加。因此，基于物理的模型應用范圍受到限制[8]。修正通用土壤流失方程(RUSLE)[9]是目前國內外應用最為廣泛的一種經驗模型，它具有模型結構簡單、參數獲取容易、計算快速等優點。該模型與GIS/RS技術相結合，可以實現柵格尺度土壤侵蝕模擬與預測，被廣泛用于土壤侵蝕量估算[10-11]、動態變化分析[4]、土壤侵蝕空間分異特征分析[12]以及水土流失敏感性分析[13]等。基于RUSLE模型在東北黑土區的研究主要集中在流域尺度的土壤侵蝕及動態變化研究[14-16]。水土保持項目的實施以流域和縣域尺度居多[17]，而目前在東北黑土區開展的基于縣域尺度的土壤侵蝕動態研究較少深入開展影響因素方面的分析。

吉林省梅河口市是國家農業綜合開發黑土區水土流失重點治理項目縣之一，截至2017年，水利部對該區域的水土流失進行了多年治理，治理措施是否有效、效果是否顯著，是國家以及相關部門亟需了解的問題。本研究擬結合3S技術，運用RUSLE模型計算梅河口市2010年、2017年的土壤侵蝕模數，并從土地利用類型、坡度和水土保持措施3個方面深入分析土壤侵蝕變化特征及其影響因素，以期為區域水土保持治理工作提供科學依據。

1 研究區概況

梅河口市位于吉林省東南部，地理位置為42°08′—43°02′N，125°15′—126°03′E，面積2 174.60 km2。地形南北斜長、東西較窄。地勢西南和東北兩端稍高，中部較低，地面高程為230～440 m。梅河口市屬北溫帶大陸性季風氣候，年平均氣溫為4.6 ℃，年平均降水量為708.3 mm。土壤類型主要有草甸土、白漿土、水稻土、沖積土、棕壤等。全市分為3個綜合區域，西南地區種植多為旱地；中部平原是以水稻、玉米為主的糧食產區；東北部是以玉米、大豆為主的低山丘陵農牧經濟區。全市主要有輝發河、一統河、大沙河三大河流，海龍水庫和堿水水庫兩大灌區。

2 數據與方法

2.1 數據來源與遙感影像處理

本研究選用的遙感數據為2017年Landsat OLI數據和2010年Landsat TM數據，軌道號為11830，11831，空間分辨率為30 m，影像可從地理空間數據云(http:∥www.gscloud.cn/)免費獲取。輔以1∶5萬地形圖對Landsat影像進行幾何校正，校正使用二次多項式和最近鄰像元法且將均方根誤差控制在0.5個像元內；其次對影像進行大氣校正。對兩期影像進行鑲嵌并使用研究區邊界進行裁剪得到研究區遙感影像數據。經過預處理的影像可作為分類的數據源。

非遙感數據包括：①土壤數據。土壤理化性質源自全國第二次土壤普查成果，土壤類型分布源自1∶100萬全國土壤數據。②1∶5萬數字化地形圖數據。包含等高線、高程點、水系等要素。③降雨數據。源自梅河口市氣象站點30 a逐日侵蝕性降雨量資料。④文字資料。包含水土流失治理項目可行性報告和工程措施設計圖等，由松遼水利委員會提供。

結合遙感影像的色調、形狀、紋理及空間位置等特征，進行野外調查，以此建立土地利用解譯標志。根據解譯標志，在ArcGIS軟件中，采取人機交互的方式對地物進行勾繪、解譯，得到研究區的土地利用類型數據。根據解譯結果和影像預設野外調查樣點，使用GPS獲取野外驗證樣點64個，對土地利用分類結果進行驗證。最后以遙感影像為基礎，結合水土流失治理項目范圍圖和實地調查結果，解譯獲取研究區水土保持措施分布結果。

2.2 RUSLE模型及土壤侵蝕強度分級

2.2.1 RUSLE模型 本研究利用修正的美國通用水土流失方程(RUSLE)[9]對梅河口市的土壤侵蝕進行定量評價。RUSLE模型為：

A=R×K×L×S×C×P

(1)

式中：A——土壤侵蝕模數；R——降雨侵蝕力因子；K——土壤可蝕性因子；L——坡長因子；S——坡度因子；C——植被覆蓋與管理因子；P——水土保持措施因子。下同。

(1) 降雨侵蝕力因子R。降雨侵蝕力因子是描述降水導致土壤侵蝕潛在能力大小的定量指標。本研究的降雨資料為氣象站點近30 a的逐日侵蝕性降雨量資料，因此選擇章文波等[18]提出的降雨侵蝕力計算模型，其模型為：

(2)

式中：M——某半月時段的侵蝕力值；k——半月內的天數(d)；Pj——半月時間內第jd的侵蝕性日雨量(mm)，其日雨量大于12 mm，否則以0計算；半月時段以每月15日為界進行劃分，每月前15 d作為一個半月，剩下的天數作為另一個半月，全年依次劃分為24個時段；α，β——模型待定參數。

β=0.836 3+18.144/Pd12+24.455/Py12

(3)

α=21.586β-7.189 1

(4)

式中：Pd12——日降雨量≥12 mm的日平均降雨量(mm)；Py12——大于侵蝕性日雨量的年平均雨量(mm)。

基于式(3)—(4)，計算研究區逐年各半月的降雨侵蝕力，累加得到年降雨侵蝕力，取多年平均得到年均降雨侵蝕力。

(2) 土壤可蝕性因子K。土壤可蝕性因子是反映土壤性能和土壤被蝕的難易程度的指標。本研究參考第一次全國水利普查——水土保持情況普查中的公式計算公式K因子[19]計算公式如下：

(5)

(3) 坡度和坡長因子。坡度因子S和坡長因子L用來反映地形地貌特征對土壤侵蝕的影響，通常放在一起考慮，視為地形因子。RUSLE模型中的公式是通過大量試驗得到的經驗公式，適合緩坡區域，針對梅河口市坡度在9%～55%陡坡土壤侵蝕，采用國內研究者發展研究的坡度因子S的計算公式[20]為：

(6)

坡長因子L的計算公式為：

L=(λ/22.1)m

(7)

(8)

式中：S——坡度因子；θ——坡度值(°)；L——坡長因子；λ——坡長(m)。

(4) 植被覆蓋與管理因子C。植被覆蓋與管理因子(C)主要反映有關植被覆蓋和變化對土壤侵蝕的綜合作用，是侵蝕動力的抑制因子。本研究根據相關研究成果[21-25]，基于前期獲取的土地利用類型數據，采用直接賦值的方法得到研究區不同土地利用類型的C因子值(表1)。

表1 梅河口市不同土地利用類型C因子值

(5) 水土保持措施因子P。水土保持措施因子是指其他條件相同時，標準小區實施水土保持措施后的土地土壤流失量與標準小區上順坡耕作土壤流失量之比[21]。借鑒相關參考文獻[24,26]，結合本研究實際水保情況，得到研究區不同水土保持措施類型的P因子值(表2)。

2.2.2 土壤侵蝕強度分級標準 本研究將土壤侵蝕模數按照黑土區水土流失綜合防治標準(SL446-2009)[27]土壤侵蝕強度分級標準進行分級(表3)。

表2 梅河口市不同水土保持措施類型P因子值

表3 土壤侵蝕強度分級標準

3 結果與分析

3.1 土壤侵蝕強度分析

將上述各因子的柵格圖，設置統一的投影和柵格分辨率，通過因子疊加運算得到土壤侵蝕模數，并按黑土區的土壤侵蝕標準進行分級，得到研究區2010年、2017年兩個時期的土壤侵蝕強度分級結果(表4)。

表4 梅河口市土壤侵蝕強度分級統計

2010年、2017年研究區的土壤侵蝕以輕度及以下侵蝕強度為主，占整個土壤侵蝕面積的比例分別為84.62%和85.10%；其次為中度侵蝕、強烈侵蝕和劇烈侵蝕，極強烈侵蝕面積所占比例最少。研究區整體土壤侵蝕較輕，但強烈及以上侵蝕強度仍占有一定比例，說明研究區局部區域水土流失嚴重。2010—2017年，研究區的平均土壤侵蝕模數由698.75 t /(km2·a)下降至678.25 t /(km2·a)，土壤侵蝕狀況有所改善。不同強度等級的土壤侵蝕呈現出“兩增四減”的趨勢，即微度和輕度侵蝕所占面積和比例有所增加，中度及以上強度的土壤侵蝕面積和比例減少。不同強度對應平均侵蝕模數以劇烈侵蝕下降較多，輕度和中度侵蝕略有下降，而其他等級的平均侵蝕模數基本不變。

2010年、2017年微度侵蝕呈東北—西南帶狀分布(圖1)，但2017年微度侵蝕分布范圍較2010年有所擴大，擴大區域主要位于研究區的北部和中部。2010年，輕度和中度侵蝕分布主要分布在研究區的南部和北部；強烈和極強烈侵蝕零散分布在研究區北部；劇烈侵蝕小片狀分布在研究區南部。2017年，輕度侵蝕主要位于研究區西南部，較2010年分布范圍有所縮小；強烈和極強烈的分布區域與2010年基本一致；劇烈侵蝕以東南部為主且分布范圍較2010年有所擴大。

圖1 梅河口市土壤侵蝕強度及變化分布

對2017年、2010年的土壤侵蝕強度分級結果求差，得到研究區土壤侵蝕變化空間分布(圖1)。其中，小于0的值表示侵蝕減輕，0表示侵蝕狀況基本不變，大于0的值表示侵蝕加重。研究區土壤侵蝕狀況為基本不變>侵蝕減輕>侵蝕加重，對應面積和所占比例依次為1 880.99，214.60，79.0 km2，86.50%，9.87%，3.63%。土壤侵蝕減輕的區域主要位于研究區的北部和東北—西南帶狀區域，其他區域零星分布；土壤侵蝕加重的區域主要在研究區的東南部、西南小片狀以及中部小帶狀區域。

3.2 土壤侵蝕與影響因素分析

3.2.1 土地利用與土壤侵蝕 不同土地利用類型對應不同植被覆蓋和人類干擾程度，不合理的人類活動會加劇土壤侵蝕[3]。研究區的土地利用類型包含耕地、林地、草地、居工用地(居民地與工礦交通運輸用地)、水域和其他用地(包含未利用地、裸巖、沙地等)6大類。2010年、2017年研究區土地利用以耕地和林地為主，二者之和分別占土地利用總面積的89.96%和86.53%，其次為居工用地、水域、草地和其他用地。2010—2017年耕地和林地的面積分別減少41.29 km2和33.24 km2；居工用地、水域、草地、其他用地的面積分別增加39.34 km2，1.47 km2，31.99 km2，1.74 km2。其中，居工用地增加的面積主要來自耕地；草地增加的面積主要來自耕地和林地。

不同土地利用類型對應不同侵蝕強度(表5)。2010年、2017年耕地以中度及以下強度侵蝕為主，強烈及以上強度侵蝕所占比例也較大，約10%。林地、居工用地和其他用地以輕度及以下侵蝕強度為主，所占比例超過對應土地利用類型總侵蝕的90%。草地和水域的侵蝕強度極低。2010—2017年，耕地的微度和輕度侵蝕比例增加，其他侵蝕等級所占比例減少。水域的土壤侵蝕情況基本不變。草地、居工用地、其他用地的土壤侵蝕具有類似的變化趨勢，以微度侵蝕所占比例減少，其他強度等級的侵蝕所占比例增加。其中，其他用地強烈以上侵蝕強度所占比例增加最大。

表5 梅河口市不同土地利用類型的土壤侵蝕強度

2010年不同土地利用類型的平均侵蝕模數大小為耕地>林地>其他用地>居工用地>草地>水域，2017年為耕地>其他用地>林地>居工用地>草地>水域。2010—2017年研究區平均土壤侵蝕模數減少了20.50 t /(km2·a)。其中，耕地的土壤侵蝕模數降低了37.49 t /(km2·a)，而其他用地、草地和居工用地的侵蝕模數變化較大，分別增加了82.05，48.34，33.84 t /(km2·a)，林地和水域的土壤侵蝕模數略有增加。

研究區中部小塊條狀區域的土地利用(附圖7—8)由耕地轉換為其他用地導致區域土壤侵蝕加重。其他用地一般為裸巖、沙地和未利用地等，缺乏植被對地面的保護作用，而雨水對地表的強烈沖刷作用使土壤中的顆粒物發生遷移，從而更易發生侵蝕。研究區東南部的土地利用由林地轉換為耕地和草地，西南小片狀區域則由水域轉換為草地，類似相對難以發生侵蝕的土地利用類型向易于發生侵蝕的土地類型的轉換加重了局部區域的土壤侵蝕。總之，土地利用類型的轉換實質為對應的C因子取值發生改變，在同等條件下，C因子取值越大，土壤侵蝕越嚴重。

3.2.2 坡度與土壤侵蝕 坡度是影響土壤侵蝕發生的主要因素之一，侵蝕強度隨著坡度的變化而變化。本研究根據黑土區的坡度侵蝕劃分標準[27]將坡度劃分為9個等級，分析不同坡度下土壤侵蝕特征(圖2)。土壤侵蝕發生在坡度0°～3°的面積為1 458.6 km2，占整個研究區的三分之二；土壤侵蝕分布在坡度3°～5°，5°～15°，15°～25°的比例分別為10.05%，17.88%，4.23%；不足1%的土壤侵蝕發生在坡度大于25°的區域。總體上，研究區土壤侵蝕面積隨著坡度的增加而減少，95%的土壤侵蝕分布在坡度小于15°的區域，侵蝕面積達2 065.87 km2。因此，水土流失治理應集中在坡度低于15°的區域。

2010年、2017年，微度侵蝕主要發生在坡度小于1.5°的區域，侵蝕面積所占比例分別為82.73%和84.72%，隨著坡度的升高所占比例有逐漸減少的趨勢，在3°以上趨于平穩。輕度侵蝕主要分布在0.25°～3°和5°～15°的坡度帶上，兩年所占比例分別為70.75%和68.41%。強烈和極強烈侵蝕分布具有一致性，0°～3°所占比例較低，隨著坡度的增加侵蝕分布比例增加，以5°～8°分布比例最大。劇烈侵蝕則以5°以上分布為主，在8°～15°帶上所占面積比例最大。

2010年、2017年，平均土壤侵蝕模數總體隨著坡度的增加而增大。2010年，0°～4°帶上對應輕度及以下強度侵蝕，4°～25°對應中度侵蝕。2017年，0°～5°對應輕度及以下強度侵蝕，15°～25°是強烈侵蝕，其他坡度帶為中度侵蝕。2017年與2010年相比，<8°的各個坡度帶上的土壤侵蝕模數都有所下降，而8°～25°帶上的侵蝕模數有所增加，以15°～25°坡度帶侵蝕模數增加最多。

圖2 梅河口市坡度-土壤侵蝕關系

3.2.3 水土保持措施與土壤侵蝕 研究區的水土保持措施包含地埂、改壟、梯田、封禁治理、水保林和其他治理措施6大類(表6)。首先，根據研究區5人班小流域2017年水保措施實施面積和遙感解譯措施面積計算水保措施解譯面積的準確率，得到解譯面積總體精度為88.93%，梯田、水保林、封禁治理和其他措施解譯的準確率依次為85.19%，91.49%，87.07%，92.31%。其次，在ArcGIS軟件中以水保措施斑塊為基礎數據，運用隨機點工具獲取水保措施驗證點，并結合谷歌高分影像和專家知識對驗證點的精度進行驗證，獲得水保措施驗證點的精度為81.25%。解譯獲取研究區2011—2013年水保措施實施面積為146.59 km2，水保林所占比例超過一半。2014—2017年實施水保措施的面積為25.84 km2，以地埂、水保林和封禁治理為主。2011—2017年，共實施水保措施面積172.43 km2，占整個研究區面積的7.9%。其中，水保林實施面積最大為83.27 km2，占水土保持措施比例為48.29%；其次為地埂和改壟，實施面積和所占比例依次為36.21 km2，27.54 km2，21.01%，15.97%；封禁治理和其他治理措施所占比例較少。

總體上，2011—2013年水土保持措施分布在研究區的北部和中部，2014—2017年分布在其南部和東南部。其中，97.32%的水土保持措施分布在坡度<15°，高程<435 m的區域，其他區域的水土保持措施分布極少。不同類型的水土保持措施分布區域也不同。水保林主要分布在坡度<15°的北部和中部地區，改壟主要分布于研究區北部，地埂和梯田主要分布在坡度<8°的北部和東南部地區，封禁治理和其他水土保持措施主要分布在坡度5°～25°之間的南部和東南部區域。坡度>25°的區域分布著極少比例的封禁治理和水保林。水土保持措施的分布區域與土壤侵蝕減輕區域基本一致。因此，水土保持措施的實施是該地區土壤侵蝕減輕的主要原因。但部分水保項目附近存在少量侵蝕加重的區域，可能是因為水保措施實施過程中的一些不合理操作造成的新的水土流失。

表6 梅河口市各種水土保持措施面積統計

不同水土保持措施的水保效果存在差別(圖3)，2017表示研究階段不同水保措施實施位置處的平均土壤侵蝕模數，2010年對應無水保措施的土壤侵蝕模數。水土保持措施的實施在不同程度上降低了區域的土壤侵蝕。其中，封禁治理措施效果最明顯，對應土壤侵蝕模數降低了1 573.47 (km2·a)；其次為其他措施、梯田和地埂，平均侵蝕模數分別減少了749.43，451.18，344.06 t/(km2·a)；改壟和水保林實施區域的土壤平均侵蝕模數變化相對較小，分別為100.46 t/(km2·a)和35.51 t/(km2·a)。水保措施的治理效果主要表現在坡耕地以及荒山荒坡和林草稀疏區域的治理上。荒山荒坡和林草稀疏區主要治理措施為封禁治理和水保林，封禁治理可以防止人畜進入林區破壞植被，且與水保林共同提高林地的郁閉度，防止雨滴直接濺擊地表土壤，從而減輕區域水土流失。坡耕地區域，改壟和地埂植物帶的實施可以達到初步控制水土流失的作用，而結合梯田修建等水保措施的實施在很大程度上減緩或控制了區域的土壤侵蝕。

圖3 水土保持措施-土壤侵蝕關系

2010—2017年研究區的土壤侵蝕模數變化較小，僅減少20.5 t/(km2·a)，但治理區域的平均土壤侵蝕模數存在明顯降低，由2010年的506.55 t/(km2·a)下降至2017年的352.47 t/(km2·a)，土壤侵蝕量也由8.73×104t減少至6.08×104t，說明區域水土流失治理工作取得一定成效，但治理后的區域仍存在一定強度的侵蝕。因此，區域水土流失治理工作有待加強。

4 結論與討論

(1) 梅河口市2010年、2017年土壤侵蝕以微度和輕度侵蝕為主，其平均土壤侵蝕模數分別為698.75，678.25 t/(km2·a)，土壤侵蝕狀況有所改善。2010—2017年土壤侵蝕變化趨勢為基本不變>侵蝕減輕>侵蝕加重。其中，土壤侵蝕減輕的區域主要位于研究區的北部和東北—西南帶上，侵蝕加重的區域主要位于研究區的東南部及其他小片狀區域。

(2) 梅河口市的土地利用以耕地和林地為主，而土壤侵蝕強度表現為耕地>其他用地>林地>居工用地>草地>水域。2010—2017年，耕地的侵蝕模數減少了37.49 t/(km2·a)，其他用地的侵蝕模數增加最多為82.05 t/(km2·a)。耕地向居工用地和其他用地的轉換以及林地向耕地和草地的轉換等加重了區域的土壤侵蝕。

(3) 梅河口市土壤侵蝕面積隨著坡度的增加而減少，其中，95%的土壤侵蝕發生在坡度<15°的區域，而土壤侵蝕強度隨著坡度的增加而增大。2010—2017年，坡度<8°的區域土壤侵蝕模數下降，土壤侵蝕狀況得到改善，坡度>8°的區域，土壤侵蝕有所加重。

(4) 梅河口市的水土保持措施以水保林、改壟和地埂為主，占整個水保面積的85.52%。水土保持措施主要分布在研究區北部和東南部，與土壤侵蝕減輕分布區域基本一致。2010—2017年，水保項目實施區域的土壤侵蝕模數減少了154.08 t/(km2·a)，土壤侵蝕量減少了2.65×104t，說明水保措施在改善水土保持功能上取得一定成效，但治理后的土壤仍存在較大強度的侵蝕，水土流失治理工作有待加強。

本研究的不足在于：模型中C因子取值建立在土地利用分類基礎之上，而Landsat影像數據較粗的分辨率和解譯經驗的缺乏都可能導致分類結果的不準確，進而影響C值。其次，p值數據的獲取多以黑龍江賓縣的試驗數據為準，而區域的差異可能導致p值的差異。如何獲取精確的p值以及基于高分影像數據獲取C值都有待進一步研究。