基于遙感與抽樣單元調查的縣域尺度水土流失推算方法

段 倩,齊 斐,羅夢琦,劉 霞*,唐 俊,張春強,吳鎮宇,李 想,姚孝友

1.南京林業大學南方現代林業協同創新中心,江蘇省水土保持與生態修復重點實驗室,江蘇南京210037

2.淮河水利委員會淮河流域水土保持監測中心站,安徽蚌埠233001

水土流失是全球普遍關注的生態環境問題。作為水土流失最嚴重的國家之一，我國已為治理水土流失投入了大量財力、物力和人力。為評價水土流失治理效果，開展水土流失動態監測具有重要的意義。目前，在國外水土流失調查計算過程中，美國主要采用抽樣調查和通用土壤流失方程（USLE），澳大利亞主要采用基于GIS 和USLE 的水土流失網格估算，歐洲主要采用基于專家法和USLE 模型法[1-4]等。在國內，水土流失狀況評價方法主要有定性評價法（綜合評判法）以及中國土壤流失方程（Chinese Soil Loss Equation,CSLE）模型法。前者在我國80 年代、90 年代的水土流失遙感調查中廣泛應用[5]；后者首次與抽樣調查相結合大規模應用于2010～2012 年開展的第一次全國水利普查水土保持情況普查中[6-9]，并在2018 年開始的全國水土流失動態監測用于水力侵蝕計算。在CSLE 模型應用中發現，抽樣密度會影響小尺度尤其是縣域尺度的水土流失狀況推算結果，但是影響大小具有地域性差異；而且不同的水土流失推算方法也會影響最終的結果。趙維軍、張巖等[10-12]以陜西吳起縣為例，對比分析不同密度對水土流失因子精度影響，發現1%密度和4%密度均能很好地代表吳起縣水土流失狀況，并在1%抽樣密度基礎上對比CSLE 模型和綜合評判法估算成果，認為前者更有優越性，但是二者以抽樣單元內分析為主，均未涉及到不同抽樣密度、不同推算方法對縣域水土流失狀況估算的影響；鄒叢榮等[13]在山東蒙陰縣研究1%、4%抽樣密度與三種水土流失推算方法對水土流失狀況的影響，發現4%抽樣密度適用于直接外推法和插值外推法，采用1%抽樣密度適用于全覆蓋計算法（柵格計算法）。故本研究以沂蒙山泰山國家級重點治理區沂水縣為例，分別探索不同抽樣密度、不同推算方法對縣域水土流失狀況的影響，對比分析差異和原因，進一步探討適用于縣域尺度水土流失調查與評價的抽樣密度和推算方法，為縣域水土流失動態監測方法提供數據支撐。

1 材料與方法

1.1 研究區概況

沂水縣地處沂山南麓，隸屬山東省臨沂市，面積2420.08 km2，屬北方土石山區。地貌主要以低山丘陵為主，地勢自西北向東南傾斜，山丘區面積比94.3%。土壤包括粗骨土、棕壤、褐土、紅粘土和潮土五大類，以粗骨土為主。屬暖溫帶季風型大陸性氣候，多年平均氣溫12.3 ℃，多年年均降水量629 mm。水資源豐富，境內沂河為山東省第一大河，以及山東省第三大水庫跋山水庫。植被屬暖溫帶落葉闊葉林區，主要喬木樹種有刺槐（Robinia pseudoacaciaL.）、側柏（Platycladus orientalis（Linn.）Franco）、花椒（Zanthoxylum bungeanumMaxim.）、油松（Pinus tabulaeformisCarr.）等，自然灌木與草本植物主要有黃荊（Vitex negundoL.）、三裂繡線菊（Spiraea trilobataL.）、胡枝子（Lespedeza bicolorsTurcz.）等。地理位置見圖1。

1.2 野外抽樣單元數據采集與處理

1.2.1 野外抽樣單元布設采用分層不等概抽樣法[7,9,14,15]，以國家第一次水利普查布設的野外抽樣單元為基礎，在山丘區增加抽樣單元，勾繪0.2～3 km2小流域，使其抽樣密度增加至4%，平原區維持不變。全縣1%抽樣密度調查單元24 個，單元控制區域100 km2；4%調查單元96 個，單元控制區域25 km2。野外調查點分布見圖2。

圖1 沂水縣地理位置Fig.1 The location of Yishui County

圖2 沂水縣1%和4%抽樣單元分布Fig.2 The location of field units in Yishui County

1.2.2 野外抽樣單元數據采集與處理基于野外調查底圖和現場實際狀況進行地塊邊界勾繪，填寫野外調查信息表，詳細記錄單元內各地塊土地利用、水土保持措施等信息[7]。調查結束后，將野外調查地圖配準矢量化，逐地塊填寫地塊信息。

1.3 水土流失推算方法

1.3.1 水土流失計算模型本研究采用中國水土流失通用模型（CSLE）[11]，根據《北方土石山區水土流失綜合治理技術標準》（SL 665-2014）進行水土流失狀況統計分析。CSLE 模型基本形式為：

式中：A為土壤流失量（t·hm-2a-1）；R為降雨侵蝕力因子（MJ·mm·hm-2·h-1·a-1）；K為土壤可蝕性因子（t·hm2·h·hm-2·MJ-1·mm-1）；L、S為坡長、坡度因子（無量綱）；B、E、T為水土保持措施因子（無量綱）。

降雨侵蝕力因子（R）：基于沂蒙山區30 年88 個雨量站點的日降雨數據，采用章文波[16]公式計算R，利用普通克里金插值法獲取30 m 區域柵格圖層[17]，通過裁剪獲取沂水縣R 柵格圖層。

土壤可蝕性因子（K）：根據Williams 模型[18]和徑流小區觀測資料計算K值，并基于全縣土壤類型圖，獲取全縣土壤可蝕性因子圖層。

坡長、坡度因子（L、S）：基于1:1 萬地形圖，經投影變換后，生成10 m 柵格大小DEM，采用劉寶元[19]修正算法，提取坡度坡長因子。

水土保持措施因子（BET）：基于徑流小區數據和水利普查中提供的措施因子參考值，對地塊進行賦值[14,16]。

1.3.2 水土流失推算方法主要采用直接外推法、插值外推法和全覆蓋計算法[13,14]。

直接外推法和插值外推法均是以抽樣單元內水土流失狀況推算縣域水土流失狀況。前者分別將1%和4%密度抽樣單元內水土流失狀況按照各級土壤侵蝕強度面積百分比推算至其控制區域，并匯總，此法無法直接成圖；后者則對輕度及其以上土壤侵蝕強度面積比和水土流失面積比，采用普通克里金插值，經平衡計算后獲取縣域水土流失狀況，采用累積百分比分級法獲取縣域水土流失狀況空間分布示意圖。

全覆蓋計算法是基于詳細的土地利用矢量數據（根據SPOT 6 獲取），進行BET 因子賦值，疊加縣域R、K、L、S、BET 因子圖層，獲取縣域水土流失狀況及空間分布。其中，B因子是根據縣域林草植被覆蓋度結合野外調查數據進行賦值；E、T因子分別根據野外調查數據采用面積加權平均法計算并賦值。

1.3.3 相對誤差分析采用測量學相對誤差的概念，以某一密度或某一推算方法結果為基準，來反映不同密度不同推算方法的差異程度。

式中：δi為第i級土壤侵蝕強度相對誤差（%）；Δi為第i級土壤侵蝕強度絕對誤差，此文中為兩種對應結果之差，在此為方便后續計算取其絕對值；Li為真值，此文中同種密度下不同推算方法結果比較分析時以全覆蓋計算法結果為真值，同種推算方法下不同密度結果比較分析時以4%抽樣密度結果為真值；δˉ為平均相對誤差（%）；Pi為第i級土壤侵蝕強度面積百分比（無量綱）；i=1,2,…,6，土壤侵蝕強度等級。

2 結果與分析

2.1 同一抽樣密度不同推算方法水土流失特征

2.1.1 1%抽樣密度下水土流失數量及空間差異在1%抽樣密度下，沂水縣直接外推法水土流失面積832.93 km2，占34.42%；插值外推法水土流失面積811.94 km2，占33.55%；全覆蓋計算法水土流失面積1053.06 km2，占43.51%。土壤侵蝕強度均是以微度侵蝕為主，其次為輕度侵蝕。見圖3。

通過不同推算方法結果對比發現，直接外推法、插值外推法各級土壤侵蝕強度面積和水土流失面積差異較小，水土流失面積僅相差20.99 km2，占全縣0.87%。但二者與全覆蓋計算法的結果存在差異，水土流失面積分別偏低220.13 km2、241.12 km2；輕度及其以上侵蝕相對誤差較高為12.18～250.33%。在1%密度下，以全覆蓋計算法水土流失面積為基準，直接外推法和插值外推法相對誤差分別為20.90%和22.90%，各級土壤侵蝕強度平均相對誤差分別為22.69%和24.84%。見圖3。

從空間分布來看，1%抽樣密度下插值外推法和全覆蓋計算法存在明顯差異。插值外推法水土流失區域集中分布在西部和西北部，而全覆蓋計算法分布相對零散。見圖4。

圖3 1%抽樣密度下三種推算方法各級土壤侵蝕強度面積對比Fig.3 The results of three estimation methods based on 1%sampling ratio

圖4 1%密度下插值外推法和全覆蓋計算法土壤侵蝕分布圖Fig.4 Soil erosion intensity by different calculation methods based on 1%sampling density

2.1.2 4%抽樣密度下水土流失數量及空間差異4%密度下，直接外推法水土流失面積990.03 km2，占40.91%；插值外推法水土流失面積998.04 km2，占41.24%；全覆蓋計算法水土流失面積1074.17 km2，占44.39%。土壤侵蝕強度均是以微度侵蝕為主，其次為輕度侵蝕。見圖5。

通過對比發現，直接外推法和插值外推法計算結果差異極小，水土流失面積僅相差8.01 km2，占0.33%，各級土壤侵蝕強度面積相差15 km2以下。且二者與全覆蓋計算法的計算結果差異減小，水土流失面積分別低84.14 km2和76.13 km2，各占3.48%和3.15%，但中度及其以上侵蝕仍有較大差異，相對誤差為17.47～31.71%。4%密度下，以全覆蓋計算法水土流失面積為基準，直接外推法和插值外推法相對誤差分別為7.09%和7.83%，各級土壤侵蝕強度平均相對誤差分別為10.04%和9.66%。見圖5。

圖5 4%抽樣密度下三種推算方法各級土壤侵蝕強度面積對比Fig.5 The results of three estimation methods based on 4%sampling ratio in Yishui

從空間分布來看，4%抽樣密度下，插值外推法和全覆蓋計算法水土流失空間分布差異較大，前者水土流失區域片狀分布，中度及其以上侵蝕主要分布在西部，東部各鄉鎮幾乎均為微度侵蝕，而后者呈均勻異質性零散分布，在全縣各個鄉鎮廣泛存在。見圖6。

圖6 4%密度下插值外推法和全覆蓋計算法土壤侵蝕分布圖Fig.6 Soil erosion intensity by different calculation methods of 4%sampling density

2.2 同一推算方法下不同抽樣密度水土流失特征與分布

2.2.1 1%和4%密度下直接外推法結果差異對比1%、4%抽樣密度，沂水縣水土流失面積隨抽樣密度增大而升高，相對誤差15.87%；隨著侵蝕強度的增加，各級土壤侵蝕強度相對誤差變大，其中劇烈侵蝕相對誤差達96.00%。直接外推法只能根據抽樣單元水土流失狀況，大體估算水土流失分布區域。

2.2.2 1%和4%密度下插值外推法結果差異對比1%、4%抽樣密度，水土流失面積同樣隨抽樣密度增大而升高，相對誤差達18.65%；隨著侵蝕強度的增加，各級土壤侵蝕強度相對誤差變大，其中劇烈侵蝕相對誤差達達129.99%。

插值外推法可插值形成水土流失分布示意圖，但空間分布受抽樣單元影響較大。兩種密度下水土流失區域空間分布基本一致，但是4%抽樣密度東部和北部范圍擴大，且中度侵蝕面積增加。見圖4 和圖6。

2.2.3 1%和4%密度下全覆蓋計算法結果差異對比1%、4%抽樣密度，水土流失面積差異較小，僅相差21.11 km2，相對誤差1.97%；從各級土壤侵蝕強度面積來看，強烈及以上侵蝕面積差異較大，尤其是劇烈侵蝕，由于其面積基數小，相對誤差仍達54.71%。

全覆蓋計算法可成詳細的水土流失圖，兩種密度下，土壤侵蝕強度和水土流失分布規律相似，這與全覆蓋計算法是基于全縣因子圖層，以柵格為單元進行計算有關。見圖4 和圖6。

3 討論

（1）主要受坡度影響，直接推算法和插值推算法結果在兩種抽樣密度下差異較大。在CSLE模型各因子中，4%抽樣單元S 因子均值明顯高于1%抽樣密度結果。根據統計數據分析，沂水縣山丘區面積百分比94.33%，平緩坡占51.80%，而1%密度平緩坡占68.06%，4%抽樣密度平緩坡占52.45%，4%密度平緩坡比例接近全縣比例，并明顯低于1%抽樣密度，從而導致了1%抽樣單元的S 因子明顯偏高。

（2）與鄒從榮等[13]結果發現一致，直接外推法和插值外推法結果受抽樣密度影響大，全覆蓋計算法結果受抽樣密度影響小。鄒叢榮等[13]研究發現蒙陰縣三種推算方法下水土流失面積隨抽樣密度增加而降低，直接外推法和插值外推法結果受抽樣密度影響大，全覆蓋計算法結果受抽樣密度影響小；而本文研究發現沂水縣采用三種推算方法時水土流失面積隨抽樣密度增加而升高，在4%密度下，三種推算方法結果數量特征基本一致。蒙陰縣和沂水縣同屬于沂蒙山區，但是二者土地利用方式和石灰巖山區分布有明顯不同，蒙陰縣土地利用以林地和園地為主，全縣內石灰巖山區成條帶狀分布，沂水縣以林地和耕地為主，石灰巖山區在西部成片狀分布。此外，蒙陰縣1%密度下抽樣單元土壤可蝕性因子、坡度因子和工程措施因子均高于4%密度抽樣單元，而沂水縣1%密度下抽樣單元坡度因子低于4%密度抽樣單元，其余因子差異不明顯。

4 結論

（1）1%抽樣密度下，直接外推法與插值外推法計算數據結果相近，但與全覆蓋計算法結果有差異；4%抽樣密度下，三種推算方法計算數據結果接近，但是水土流失空間分布差異明顯；

（2）直接外推法與插值外推法受抽樣密度影響較大，各級土壤侵蝕強度相對誤差一般大于10%，水土流失面積相對誤差和平均相對誤差均高于15%；全覆蓋計算法受抽樣密度影響較小，水土流失面積相對誤差和平均相對誤差均小于10%，但各級土壤侵蝕強度級別分布存在一定差異；

（3）受坡度影響，1%抽樣密度下直接外推法和插值外推法水土流失面積低于4%抽樣密度結果；

（4）需要根據工作需要選擇合適的抽樣密度和推算方法。若需要縣域詳細的水土流失分布圖，可采用1%抽樣密度下的全覆蓋計算法；若采用直接外推法或插值外推法，宜采用4%抽樣密度布設調查單元。