基于CSLE和抽樣調查的泰國土壤水蝕評價

杜朝正, 楊勤科, 張加瓊, 馬 波, 王春梅,4, 龐國偉,4

(1.西北大學 城市與環境學院, 西安 710127; 2.菏澤學院 城市建設學院, 山東 菏澤 274015; 3.中國科學院水利部 水土保持研究所, 陜西 楊陵 712100; 4.旱區生態水文與災害防治國家林業和草原局重點實驗室, 西安 710127)

土壤是人類生存與發展的重要自然資源之一[1]，土壤可以為人類提供食物、能源、水源，在維持生態系統服務、生物多樣性方面也具有重要作用[2]。但是，全球大部分地區的土壤受到水蝕的威脅[3]。土壤水蝕是糧食生產和可持續農業面臨的最嚴重威脅之一[4-5]。土壤水蝕不僅會導致土地退化，生態服務水平降低，還會導致下游河流淤積，從而制約區域的可持續發展[6]。

治理土壤水蝕必須以掌握區域土壤水蝕狀況為前提[7-8]，一般可歸納為2種方法：第1種是基于區域各侵蝕因子圖層相乘，計算區域土壤水蝕速率專題圖，可稱為地圖代數法。該方法計算結果具有較好的空間分布，得到廣泛應用，代表性研究有Borrelli[9]、Panagos[10]、查良松[11]等。第2種是抽樣調查的方法，在研究區布設基本抽樣調查單元，結合土壤水蝕模型，計算各抽樣單元水蝕速率，進而推算出研究區土壤水蝕狀況，可稱為抽樣調查法，代表性研究有Nusser[12]、USDA[13]、劉寶元等[7]、Liu等[14]、Yin等[15]。但如何將地圖代數制圖結果與抽樣調查結果相結合，實現對區域土壤侵蝕的定量評價則鮮見報道。

泰國屬于發展中國家，農作區占全國的40%以上，大力發展農業導致農地不斷擴張，而森林不斷萎縮，加之粗放的耕作和管理方式，土壤侵蝕日益加劇，造成農業地區土地嚴重退化[16]。在農作區，農民習慣將坡面裸露，這使土壤完全暴露于暴雨及徑流中，土壤流失現象普遍存在[17]。在山區和高原，大規模砍伐森林、陡坡耕種(45°～60°)加之游移耕作，進一步加劇了土壤侵蝕[18-19]。目前，雖然已有學者對泰國進行土壤水蝕評價，如Rangsiwanichpong等[16]，但數據分辨率較粗，地形因子和水保措施因子計算有待改進。需進一步提高泰國土壤水蝕定量評價的準確性，為泰國土壤水蝕治理提供更為準確的數據支撐。

本研究以泰國為例，在前期工作的基礎上[20]，采用中國土壤流失方程(Chinese Soil Loss Equation，CSLE)模型，基于各侵蝕因子基礎數據和抽樣調查數據，定量計算泰國土壤水蝕速率專題圖，以期全面準確掌握泰國土壤水蝕基本狀況，為該國土壤水蝕治理提供數據支撐。

1 研究區域與基礎數據

泰國(5°30′—21°N，97°30′—105°30′E)位于東南亞的中南半島中部，總面積513 115 km2。泰國地勢西北高東南低。根據地形特征，可將泰國分為北部(山地為主)、東北部(高原為主)、中部(平原為主)、南部(丘陵地帶)4個自然區域帶。泰國地處熱帶地區，主要受熱帶季風氣候影響，年均降水量約1 150 mm。由于降雨量大，盲目開墾耕地，并且存在不合理的耕作方式，水土流失問題較為嚴重[21]。泰國是重要的農產品出口國，土壤侵蝕問題已嚴重制約該國的可持續農業發展[22-23]。

本研究所用基礎數據主要包括2個部分，一是用于地圖代數法制圖整體計算的各水蝕因子數據(表1)；二是土壤水蝕抽樣調查數據，研究區共計175個抽樣調查單元(圖1)，所有抽樣調查單元數據均來自野外調查解譯并通過了質量檢驗[8,24]。

表1 基礎數據一覽表

圖1 抽樣單元分布圖

2 研究方法

本研究首先采用地圖代數法計算研究區土壤水蝕速率專題圖，再以175個抽樣調查單元計算結果為參照，對地圖代數結果進行直方圖匹配，將匹配之后的制圖結果作為研究區最終的土壤水蝕速率專題圖。

2.1 土壤水蝕評價模型

采用中國土壤流失方程(Chinese Soil Loss Equation，CSLE)，結合GIS技術，計算研究區土壤水蝕速率。CSLE基本形式為：

A=R×K×L×S×B×E×T

(1)

式中：A為土壤水蝕速率[t/(hm2·a)]；R降雨侵蝕力因子[MJ·mm/(hm2·a·h)]；K為土壤可蝕性因子[t·hm2·h/(hm2·MJ·mm)]；L為坡長因子，無量綱；S為坡度因子，無量綱；B為植被覆蓋與生物措施因子，無量綱；E為工程措施因子，無量綱；T為耕作措施因子，無量綱。

2.2 水蝕因子數據處理

(1)R，K，LS因子數據處理：根據項目組提供的全球降雨侵蝕力因子R、泛第三極土壤可蝕性因子K，通過裁剪、重采樣處理，獲得研究區30 m分辨率的R，K因子柵格數據。將30 m分辨率的多幅分塊LS因子基礎數據進行拼接，并裁剪出研究區范圍。參考中國水利部的關于水土流失動態監測的有關技術文件和Brychta等[25]的研究，計算LS因子上限閾值(坡度30°、坡長100 m時，LS上限=21.2)，最終獲得研究區LS因子專題層數據。

(2)B因子計算：植被覆蓋與生物措施是一種重要的水保措施[26]，參考Borrelli等[9]的研究方法，計算非農業土地類型B因子值。

下載MOD44 B數據，并通過MRT(MODIS Reprojection Tool，MRT)工具提取TC(Tree Cover，TC)、NTV(NonTree Vegetation，NTV)圖層。計算各土地利用類型B因子值采用如下公式。

林地：

BNA=MINB+(MAXB-MINB)×(1-TC)

(2)

非林地：

BP=MINB+(MAXB-MINB)×(1-NTV)

(3)

式中:BNA為林地B因子值；TC為林地覆蓋度,取值為0～1；BP為非林地(即:草地、灌木、裸地等)的B因子值,NTV為植被(除林地外)覆蓋度,取值為0～1。maxB,minB為對應土地利用類型B因子值的最大值和最小值,取值如表2所示。

根據土地利用類型數據，合并各地類B因子值，最終獲得研究區B因子圖層。

表2 各土地利用類型B因子取值范圍

(3) E因子計算：工程措施是用人工或機械挖掘和修筑完成[26]，在前期工作的基礎上，依據抽樣單元布設原則，將研究區劃分為規則矢量網格，確保每個矢量網格內只包括1個抽樣調查單元，將抽樣單元的E因子平均值作為矢量網格的E因子值，再將矢量網格數據轉換為柵格數據，最終獲得研究區的E因子專題圖。

(4) T因子計算：根據水利部水土保持監測中心編寫的《水土流失動態監測技術指南》(2020年8月版)，參考中國各農業熟制區T值表的規則，為Fischer等[27]制作的全球熟制區賦予T因子值。再依據研究區土地利用類型數據，將非耕地區域的T值設置為1，將耕地區域的T值取值為全球T值圖的T值，最終獲得研究區的T因子專題圖。

2.3 土壤水蝕制圖方法

2.3.1 地圖代數法制圖 地圖代數法制圖是將事先準備好的每個土壤侵蝕因子柵格數據相乘，獲得研究區的土壤水蝕速率專題圖。為了確保計算精度及便于顯示，本研究采用WGS84坐標系，最終獲得的區域R，K，LS，B，E，T，A專題數據分辨率均設定為30m。

2.3.2 抽樣單元水蝕速率計算 研究區內共解譯了175個均勻分布的抽樣調查單元，根據解譯單元點位置的地形特征，確定調查單元的形狀為小流域(在山區)或矩形(在平原)，面積大小為0.2～3km2。解譯每個抽樣單元內的土地利用和詳細水土保持措施，其中水土保持措施包括生物措施、工程措施、耕作措施。在研究區范圍內，CSLE模型適用于計算抽樣單元水蝕速率[28]。

對區域R因子、區域K因子、區域LS因子進行多重采樣和裁剪，獲得每個抽樣單元R，K，LS因子專題數據。對各土地利用類型的區域B值圖(250m分辨率)進行多重采樣和裁剪，再根據解譯的土地利用類型和生物措施，計算每個抽樣單元的B因子專題數據。根據解譯的工程措施，為每個圖斑賦E因子值，再將矢量數據轉換為柵格數據，獲得每個抽樣單元的E因子專題數據。根據解譯的土地利用類型和水保措施信息，計算每個抽樣單元的T因子專題數據。再采用CSLE模型計算每個抽樣單元的土壤水蝕速率(針對抽樣調查單元的每個圖斑)，借助python語言實現數據批量處理，最終獲得175個以各抽樣單元為邊界的土壤水蝕速率專題圖(2.5m分辨率)。

2.3.3 直方圖匹配制圖 以均勻分布的175個抽樣調查單元為基礎，抽樣單元計算結果可全面表達水土保持措施的影響，可真實反映研究區的統計特征；地圖代數法制圖結果具有較好的空間分布特征。為了結合兩者的優點，嘗試以175個抽樣單元水蝕速率專題圖結果(針對每個抽樣單元的每個柵格值統計)的累積頻率為參照，對地圖代數制圖結果進行累積頻率直方圖匹配，最終獲得研究區新的土壤水蝕速率專題圖。

3 結果與分析

3.1 土壤水蝕因子分析

按照上述方法獲得的研究區各水蝕因子專題層如圖2所示。研究區R因子空間分布上呈現出由東南向西北逐漸減弱的趨勢，主要原因為泰國屬于熱帶季風氣候，降雨量由東南沿海向西北內陸逐漸減小。K因子空間分布呈現出與土地利用類型相關的空間格局，K因子值高的區域與耕地分布區域相對應，K因子值低的區域主要對應林地區域。LS因子呈現由西北向東南逐漸減小的分布格局，主要原因為泰國西北部為山地，地形起伏較大，坡度較陡，中部為平原區，地形較為平緩。B因子值的空間分布與土地利用類型的分布有關，按照林地、草地、耕地依次增大。E因子的分布與國家政策和當地耕作習慣有關，在泰國范圍內工程措施較少。T因子只針對農地地塊，T因子值的大小與輪作方式、管理措施相關，輪作方式、管理措施又存在區域間差異，在南部丘陵地帶T因子值較大，在東北部T因子值較小。

3.2 泰國水蝕速率特征

3.2.1 空間分布特征分析 泰國水蝕速率專題圖顯示水蝕速率高值區主要集中于泰國北部和東北部，水蝕速率低值區主要分布于中部平原區(圖3)。這與泰國的地形和耕地分布密切相關，泰國北部屬于山地，地形起伏較大，降雨容易引起土壤流失；東北部地區存在大量耕地區域，耕地區域B因子值較大，并且該區域的K，R因子值也較大。中部水蝕速率較小，主要原因為中部地形以平原為主，并且R因子相對較小。

3.2.2 統計特征分析 參考Bahadur等[29]、Borrelli等[9]的分級方案，結合泰國區域的實際情況，對泰國土壤侵蝕速率的侵蝕等級劃分結果顯示：泰國侵蝕速率以<500t/(km2·a)為主(74.2%)，隨著侵蝕等級升高，侵蝕面積比率呈減少趨勢。土壤侵蝕速率>15 000t/(km2·a)的面積僅占研究區總面積的0.6%，但年侵蝕量約占研究區侵蝕總量的21.5%，可見研究區局部侵蝕劇烈。泰國平均水蝕速率為687.9t/(km2·a)，是全球平均土壤水蝕速率[287.2t/(km2·a)]的2.4倍，個別地區達到1 000t/(km2·a)以上(占面積13.2%，占侵蝕總量72.0%)(表3)。可見，與全球平均水蝕速率相對比，研究區土壤水蝕較為嚴重。

3.2.3 熱點地區分析 參考Borrelli等[9]的研究，結合本研究分級方案，根據水蝕治理的緊迫程度，將水蝕速率以500，2 500t/(km2·a)為閾值設置為熱點地區。>500t/(km2·a)的熱點地區沿泰國周邊分布，該熱點地區地形LS因子值(均值為7.5)是泰國區域LS因子值(均值為3.2)的2.3倍，降雨侵蝕力也較大[均值為6 419.2MJ·mm/(hm2·a·h)]。>2 500t/(km2·a)的熱點地區主要分布于泰國北部和東北部，如此分布主要是地形LS因子、植被覆蓋與生物措施B因子、土地利用方式等因素共同作用的結果。該熱點地區LS因子值(均值為7.3)、B因子值(均值為0.849)是泰國區域的LS因子值(均值為3.2)、B因子值(均值為0.507)的2.3倍和1.7倍，降雨侵蝕力也較大[均值為6 587.6MJ·mm/(hm2·a·h)]，這導致該熱點地區水蝕較為嚴重(圖4)。對該熱點地區進一步分析發現，熱點地區84.1%為耕地，平均坡度為13.4°，因此，坡耕地進一步加劇了局部土壤水蝕。>500t/(km2·a)的熱點地區面積占研究區面積的25.8%，年侵蝕量占研究區侵蝕總量的84.3%，該熱點地區需制定長期治理規劃，逐步減少土壤水蝕的影響；>2 500t/(km2·a)的熱點地區面積占研究區總面積的4.7%，但年侵蝕量占研究區侵蝕總量的53.1%，該熱點地區迫切需要盡快治理，并且治理效果顯著。

注：耕地區域的B因子值為1；無工程措施區域的E因子值為1；非耕地區域的T因子值為1。

3.3 各土地利用類型的土壤水蝕分析

通過比較不同土地利用類型的土壤水蝕狀況(表4)，我們發現從耕地到草地、林地的平均水蝕速率明顯降低。耕地平均水蝕速率是草地、灌叢、林地平均水蝕速率的2.0倍、2.2倍、2.8倍，說明相對于耕地，林草地具有更好的水土保持作用。在泰國區域，耕地面積占研究區面積的比率最大(49.6%)，平均水蝕速率最高[1 020.2t/(km2·a)]，進而耕地侵蝕量占研究區總侵蝕量的比率最大(73.6%)，這是研究區平均水蝕速率較高的重要原因。

圖3 泰國土壤水蝕速率

表3 不同侵蝕等級的土壤水蝕面積及水蝕量

進一步統計后我們發現，在耕地區域中，>5 000t/(km2·a)的區域占耕地總面積的比率高達4.6%。草地、灌叢、林地在>5 000t/(km2·a)的區域占該地類面積的1.6%，0.6%，0.1%，明顯比耕地小。林地、灌叢、草地具有較好的儲水固土能力，因此，在有條件的區域，將高侵蝕級別的耕地改變為林地、灌叢、草地，對于減少水土流失具有重要意義。

3.4 主控因子分析

分類決策樹狀圖顯示了水蝕因子的重要性，終端節點映射出每個像元值歸因于分類決策樹中最頻繁的節點類(圖5—7)。結合因子節點可視化圖以及泰國土地利用類型圖可知，節點11，12，14，15主要受R因子(降雨侵蝕力因子)影響，主要分布于林地區域，在該區域由于降雨量較大而容易發生侵蝕。節點7，8主要受K因子(土壤可蝕性因子)影響，主要分布于東北部高原區的耕作區域，該區域常年耕作(包括一年多熟)，并且不合理的耕作方式和管理措施而引起土壤可蝕性較強，進而容易引發土壤侵蝕。節點4，5主要受LS因子(坡度坡長因子)影響，主要分布于中部平原區的耕作區和南部丘陵區的耕作區，在中部平原區，耕地的坡長較長；在南部丘陵區的地面起伏較大，容易引發土壤侵蝕。

4 討 論

4.1 計算結果合理性評估

本研究在UpperNamWa流域、NamChun流域、KhlongYai流域計算的土壤水蝕速率均值[2 154.6，706.5，223.1t/(km2·a)]與Bahadur[29]、Shrestha等[30]、Shrestha等[31]計算的均值[2 127，600，200t/(km2·a)]相近。本研究計算泰國平均土壤水蝕速率的結果為687.9t/(km2·a)，與Borrelli等[9]的結果類似[713.4t/(km2·a)]而與Rangsiwanichpong等[16]采用RUSLE模型的結果[812.5t/(km2·a)]存在差異，主要與采用不同分辨率數據導致LS因子準確度不同有關。Rangsiwanichpong等在計算水蝕速率過程中，采用1km分辨率的DEM計算LS因子，依據1km的DEM提取坡度，并對坡耕地賦水保措施值，由于分辨率較粗，這使計算的LS因子、水保措施因子誤差較大。而本研究采用30m分辨率的DEM數據，計算結果更貼近真實值。通過對比可知，本研究計算結果較為合理，具有一定的可信度。

圖5 土壤侵蝕因子重要性分類決策樹

與抽樣單元結果對比，泰國區域內175個抽樣單元土壤侵蝕速率的總體均值為706.8t/(km2·a)，研究區計算結果的175個相應區域的總體均值為642.7t/(km2·a)，兩者較為接近，總體計算精度為90.9%，進一步驗證了計算結果具有一定的可信度。

與實地考察對比，2018年11月，張加瓊等[20]對泰國北部的土壤侵蝕狀況進行了實地考察，并沿考察路線布設了18個抽樣單元。其18個抽樣單元水蝕速率的總體均值為1 767.7t/(km2·a)，本研究在相應18個區域計算的水蝕速率總體均值為1 267.3t/(km2·a)，準確度為71.7%，差異較大，主要原因為：計算方法不同。張加瓊等在計算考察的18個抽樣單元的B因子值時，采用NPV數據作為林下植被蓋度，這將低估生物措施的作用，使計算的B因子值偏高，進而使計算的水蝕速率偏高。本研究采用Borrelli等方法計算B因子值，結果更加合理。

圖6 主控因子節點可視化圖

圖7 土地利用類型

4.2 治理建議

(1) 進一步加強土壤水蝕的基礎理論研究。雖然泰國政府意識到土壤侵蝕的危害，但目前的土壤水蝕觀測和調查資料匱乏[20]，難以滿足土壤侵蝕治理的需求，現有資料也難以支撐準確評價水土保持效益，因此，建議進一步加強土壤侵蝕的基礎理論研究，加大人才、科技、經費支持力度，構建土壤侵蝕的學科理論體系，為制定合理的水保措施和防治土壤流失提供理論支撐。

(2) 因地制宜制定土壤侵蝕治理措施。土地利用方式及變化對土壤侵蝕有重要影響[32]，經過多年耕作，土壤抗侵蝕能力明顯下降[33]，合理的水保措施能夠提高土壤抗蝕性[34]。離散分布于林地內部的耕地大多屬于熱點地區，侵蝕較為強烈，建議在林地區域內部禁止毀林開墾，并進行退耕還林草。東北部高原區的耕作區域主控因子為土壤可蝕性因子，針對該區域內的耕地，建議制定合理的耕作措施，選擇合適的耕作措施類型和輪作制度，減少該區域的土壤流失。南部丘陵區的耕作區主控因子為坡度坡長因子，建議將坡耕地改為林草地或免耕。

(3) 提高作物單位面積產量。耕地的水蝕速率最大，而泰國49.6%區域屬于耕地，建議控制耕地面積數量，提高單位面積作物產量。泰國作物單產具有較大提升潛力。如清萊府的稻米單產3 500kg/hm2，而土壤、氣候條件類似的云南稻米單產為12 500kg/hm2[35]。

(4) 加強國際交流，引進切實可行的水土保持措施。從其他國家引進水保措施，并因地制宜加以改進，提高水土保持效益。如在坡面上覆蓋的植被及其殘留物，能明顯減少地表徑流和土壤流失，雜草覆蓋是一種保護梯田的最便宜、最有效的方法[19]，這些方法也適用于泰國局部區域。

5 結 論

(1) 在泰國區域進行土壤水蝕定量評價，應優先選擇直方圖匹配方法制圖。地圖代數法制圖結果具有較好的空間分布特征，但由于受到空間數據分辨率的限制，不能精確表達詳細水保措施；在抽樣單元尺度上，由于包括詳細水保措施，計算的抽樣單元水蝕速率較為準確，具有較為準確的統計特征，但難以精確表達空間分布。直方圖匹配制圖結果兼具兩者的優點，因此，優先選擇直方圖匹配方法在泰國區域進行土壤水蝕定量評價制圖。

(2) 雖然已有文獻(如文獻[9]、文獻[16])對泰國進行土壤水蝕定量評價，但采用分辨率較粗，而本研究采用30m分辨率數據計算泰國土壤水蝕速率，并且通過與175個抽樣調查結果對比、與實地考察結果相對比可知，本研究計算結果較為理想。文獻[16]未對研究區進行主控因子分析，而本研究通過主控因子分析發現，東北部高原區的耕作區、中部平原區和南部丘陵區的耕作區、林地區域的主控因子分別為土壤可蝕性因子、坡度坡長因子、降雨侵蝕力因子，針對不同區域的主控因子，可采取更改土地利用方式、制定合理的耕作措施等不同的治理方案，治理研究區土壤流失問題。

(3) 泰國土壤水蝕較為嚴重，局部侵蝕較為劇烈。在各土地利用類型中，耕地水蝕最為嚴重，林草地水蝕速率較小，應因地制宜地制定不同區域的土壤流失治理措施，建議控制耕地面積數量，提高耕地單位面積產量。水蝕速率>2 500t/(km2·a)的熱點地區主要分布于泰國北部和東北部，該熱點地區84.1%區域為耕地，在有條件的區域，將該熱點地區的耕地改變為林草地，對于減少水土流失具有顯著效果。