段 倩,齊 斐,羅夢琦,劉 霞*,唐 俊,張春強(qiáng),吳鎮(zhèn)宇,李 想,姚孝友

1.南京林業(yè)大學(xué)南方現(xiàn)代林業(yè)協(xié)同創(chuàng)新中心,江蘇省水土保持與生態(tài)修復(fù)重點實驗室,江蘇南京210037

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

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

1 材料與方法

1.1 研究區(qū)概況

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

1.2 野外抽樣單元數(shù)據(jù)采集與處理

1.2.1 野外抽樣單元布設(shè)采用分層不等概抽樣法[7,9,14,15]，以國家第一次水利普查布設(shè)的野外抽樣單元為基礎(chǔ)，在山丘區(qū)增加抽樣單元，勾繪0.2～3 km2小流域，使其抽樣密度增加至4%，平原區(qū)維持不變。全縣1%抽樣密度調(diào)查單元24 個，單元控制區(qū)域100 km2；4%調(diào)查單元96 個，單元控制區(qū)域25 km2。野外調(diào)查點分布見圖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 野外抽樣單元數(shù)據(jù)采集與處理基于野外調(diào)查底圖和現(xiàn)場實際狀況進(jìn)行地塊邊界勾繪，填寫野外調(diào)查信息表，詳細(xì)記錄單元內(nèi)各地塊土地利用、水土保持措施等信息[7]。調(diào)查結(jié)束后，將野外調(diào)查地圖配準(zhǔn)矢量化，逐地塊填寫地塊信息。

1.3 水土流失推算方法

1.3.1 水土流失計算模型本研究采用中國水土流失通用模型（CSLE）[11]，根據(jù)《北方土石山區(qū)水土流失綜合治理技術(shù)標(biāo)準(zhǔn)》（SL 665-2014）進(jìn)行水土流失狀況統(tǒng)計分析。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）：基于沂蒙山區(qū)30 年88 個雨量站點的日降雨數(shù)據(jù)，采用章文波[16]公式計算R，利用普通克里金插值法獲取30 m 區(qū)域柵格圖層[17]，通過裁剪獲取沂水縣R 柵格圖層。

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

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

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

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

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

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

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

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

2 結(jié)果與分析

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

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

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

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

圖3 1%抽樣密度下三種推算方法各級土壤侵蝕強(qiáng)度面積對比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%抽樣密度下水土流失數(shù)量及空間差異4%密度下，直接外推法水土流失面積990.03 km2，占40.91%；插值外推法水土流失面積998.04 km2，占41.24%；全覆蓋計算法水土流失面積1074.17 km2，占44.39%。土壤侵蝕強(qiáng)度均是以微度侵蝕為主，其次為輕度侵蝕。見圖5。

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

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

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

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

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

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

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

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

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

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

3 討論

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

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

4 結(jié)論

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

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

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

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