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

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

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

1.3 水土流失推算方法

1.3.1 水土流失計(jì)算模型本研究采用中國(guó)水土流失通用模型（CSLE）[11]，根據(jù)《北方土石山區(qū)水土流失綜合治理技術(shù)標(biāo)準(zhǔn)》（SL 665-2014）進(jìn)行水土流失狀況統(tǒng)計(jì)分析。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為坡長(zhǎng)、坡度因子（無(wú)量綱）；B、E、T為水土保持措施因子（無(wú)量綱）。

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

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

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

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

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

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

全覆蓋計(jì)算法是基于詳細(xì)的土地利用矢量數(shù)據(jù)（根據(jù)SPOT 6 獲?。M(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)平均法計(jì)算并賦值。

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

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

2 結(jié)果與分析

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

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

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

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

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

圖4 1%密度下插值外推法和全覆蓋計(jì)算法土壤侵蝕分布圖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%；全覆蓋計(jì)算法水土流失面積1074.17 km2，占44.39%。土壤侵蝕強(qiáng)度均是以微度侵蝕為主，其次為輕度侵蝕。見(jiàn)圖5。

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

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

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

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

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

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

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

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

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

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

3 討論

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

4 結(jié)論

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

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

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

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