李文龍,袁 利,魏文杰,劉卓昊,高睿瑜,牛 勇,趙傳普，張榮華?

(1.山東泰山森林生態(tài)系統(tǒng)國家定位研究站，山東農(nóng)業(yè)大學(xué)林學(xué)院，271018，山東泰安；2.淮河水利委員會淮河流域水土保持監(jiān)測中心站，233001，安徽蚌埠)

及時準確掌握區(qū)域水土流失情況、推動水土流失科學(xué)治理，已成為生態(tài)環(huán)境建設(shè)的迫切需要。隨著科技的發(fā)展，水土流失監(jiān)測的方法和體系正逐步完善[1]。區(qū)域尺度水土流失監(jiān)測包括基于土壤侵蝕分類分級標準的3因子定性分析法、基于監(jiān)測評價多因子指標體系的綜合指數(shù)法、基于水土流失調(diào)查單元的抽樣統(tǒng)計分析法、遙感和地統(tǒng)計法和基于物理經(jīng)驗?zāi)Ｐ偷哪Ｐ头╗2-4]等5類評價方法。

2018年，水利部實施《全國水土流失動態(tài)監(jiān)測規(guī)劃(2018—2022年)》，以縣為單元開展水土流失動態(tài)監(jiān)測，掌握全國水土流失情況[5]。2020年，水利部水土保持監(jiān)測中心頒布《2020年度水土流失動態(tài)監(jiān)測技術(shù)指南》(簡稱“技術(shù)指南”)[6]，明確水力侵蝕計算基于CSLE模型，結(jié)合野外調(diào)查，利用柵格統(tǒng)計的方法，計算土壤侵蝕模數(shù)。2021年，水利部水土保持監(jiān)測中心印發(fā)《水土流失圖斑落地工作技術(shù)細則(試行)》(水保監(jiān)〔2021〕41號)，以年度水土流失動態(tài)監(jiān)測成果為基礎(chǔ)，通過水土流失圖斑確定與判斷軟件，將柵格形式表達的土壤侵蝕強度圖概化為以矢量形式表達的土壤侵蝕強度圖，為水土保持規(guī)劃設(shè)計和綜合治理等服務(wù)。

目前已有學(xué)者開展基于不同空間數(shù)據(jù)結(jié)構(gòu)的土壤侵蝕評價：吳思穎等[7]、陳鋒等[8]基于柵格數(shù)據(jù)結(jié)構(gòu)對安溪縣、元陽縣等進行土壤侵蝕計算、統(tǒng)計。鄒海天等[9]從計算精度、空間分布、計算過程等方面對比分析外井溝小流域柵格、矢量2種數(shù)據(jù)結(jié)構(gòu)土壤侵蝕計算結(jié)果的差異。已有研究發(fā)現(xiàn)，基于柵格數(shù)據(jù)結(jié)構(gòu)進行土壤侵蝕計算、統(tǒng)計，單個矢量地塊會存在多個柵格和多種侵蝕強度，難以直接應(yīng)用于以矢量地塊為對象的水土流失治理和規(guī)劃，而直接采用矢量進行縣域尺度評價的相關(guān)研究較少，分析不同空間數(shù)據(jù)結(jié)構(gòu)對縣域尺度土壤侵蝕評價的影響還缺乏深入的研究。

鑒于此，筆者以河南省汝州市為研究區(qū)域，采用基于柵格數(shù)據(jù)結(jié)構(gòu)的3種方法和基于矢量數(shù)據(jù)結(jié)構(gòu)的2種方法，采用CSLE模型，分析縣域尺度不同計算方法和統(tǒng)計方法的土壤侵蝕差異，并探討不同方法的適用性，為縣域尺度水力侵蝕評價提供參考。

1 研究區(qū)概況

汝州市位于河南省中西部， E 112.30°～113.05°,N 33.94°～34.24°(圖1)。地勢呈南北高、中間低的特點，山地丘陵面積占78.1%。屬暖溫帶大陸性季風(fēng)氣候，四季分明。年均氣溫15 ℃，年均降水量606.5 mm，汛期集中在7月上旬至8月下旬。土壤類型主要有褐土、潮土為主，兼有部分砂漿黑土。植被為暖溫帶闊葉落葉林帶。總面積1 573 km2，土地利用以耕地、林地為主。屬伏牛山中條山國家級水土流失重點治理區(qū)，2020年水力侵蝕面積403.51 km2，占土地面積的25.65%。

圖1 研究區(qū)地理位置圖Fig.1 Geographic location map of the study area

2 材料與方法

2.1 數(shù)據(jù)及處理

1)遙感影像：2 m分辨率GF- 1和ZY- 3遙感影像(時相2020年4—6月)、250 m分辨率MODIS- NDVI產(chǎn)品(2017—2019年)。

2)土地利用/水土保持措施解譯：對影像進行投影變換、配準、鑲嵌、裁剪等處理，借助ArcGIS，按技術(shù)指南[6]，采用人機交互解譯方法獲取土地利用矢量圖斑。

3)土壤/植被覆蓋度/地形數(shù)據(jù)：土壤來源于河南省土壤類型圖(1∶20萬)；植被覆蓋度基于MODIS- NDVI產(chǎn)品，采用參數(shù)修正法[6]計算；地形采用30 m分辨率DEM圖(地理空間數(shù)據(jù)云)，通過LS提取工具[10]計算坡長(L)、坡度(S)因子。數(shù)據(jù)均重采樣為10 m分辨率，并矢量化，用于細化土地利用矢量圖斑。

4)CSLE模型因子數(shù)據(jù)：降雨侵蝕力因子(R)、土壤可蝕性因子(K)由淮河流域水土保持監(jiān)測中心站提供，耕作措施因子(T)、水土保持工程措施因子(E)、生物措施因子(B)根據(jù)技術(shù)指南[6]直接賦值或計算。數(shù)據(jù)均重采樣為10 m分辨率。

2.2 基于柵格的土壤侵蝕計算與統(tǒng)計

基于中國土壤流失方程(CSLE)，按SL190—2007《土壤侵蝕分類分級標準》[11]和技術(shù)指南[6]判定土壤侵蝕強度，利用ArcGIS軟件生成縣域土壤侵蝕強度柵格圖層。

在統(tǒng)計侵蝕結(jié)果時，采用以下3種方法：

1)柵格統(tǒng)計法：基于土壤侵蝕強度柵格，統(tǒng)計柵格數(shù)量，獲得不同土壤侵蝕強度面積。

2)技術(shù)指南土壤侵蝕地塊水土流失評價方法(簡稱“地塊評價法”)：基于土壤侵蝕強度柵格，以土地利用矢量地塊為單元，地塊內(nèi)輕度及以上侵蝕強度柵格數(shù)量超過地塊面積50%時，則判定該地塊為水土流失地塊，統(tǒng)計獲得土壤侵蝕面積。

3)水土流失圖斑確定與判斷軟件統(tǒng)計法(簡稱“軟件判斷法”)：基于土壤侵蝕強度柵格，以土地利用矢量地塊為單元，人為判定人為水土流失地塊強度，利用北京師范大學(xué)水土流失圖斑確定與判斷軟件，對其他地類以地塊侵蝕模數(shù)、高強度指數(shù)2個指標，按大小排序，以柵格統(tǒng)計的土地利用二級類型土壤侵蝕面積為目標值進行面積累加，至與目標值接近時，判定為水土流失地塊，統(tǒng)計不同土壤侵蝕強度面積。高強度指數(shù)：地塊內(nèi)強烈(含)以上高強度侵蝕面積占地塊內(nèi)水蝕或風(fēng)蝕面積之比。

軟件判定原則：以土地利用二級類為比較標準，基于圖斑計算的水蝕面積與基于柵格統(tǒng)計水蝕面積相對誤差≤0.5%，水蝕各強度級圖斑累計面積與基于柵格統(tǒng)計的水蝕各強度級面積相對誤差≤1%,基于水土流失圖斑計算的水蝕面積與基于柵格統(tǒng)計的水蝕面積相對誤差≤1%。

2.3 基于矢量的土壤侵蝕計算與統(tǒng)計

1)原始矢量地塊均值計算法(簡稱“原始矢量法”)：以土地利用矢量圖斑為單元，提取圖斑內(nèi)CSLE 7個因子均值，進行乘積運算，獲取圖斑侵蝕模數(shù)，判定土壤侵蝕強度。

2)細化矢量地塊均值計算法(簡稱“細化矢量法”)：將土壤類型、坡度、植被覆蓋度等矢量數(shù)據(jù)，逐步與土地利用矢量圖斑進行疊加，提取細化圖斑內(nèi)CSLE 7個因子均值，計算獲取細化圖斑的土壤侵蝕模數(shù)、強度。依據(jù)技術(shù)指南[6]，最小圖斑控制面積≥400 m2。

2.4 土壤侵蝕計算結(jié)果差異分析

從土壤侵蝕強度，不同土地利用土壤侵蝕，及計算量、圖斑大小、成圖效果等方面，分析不同計算和統(tǒng)計方法結(jié)果差異。

不同侵蝕強度或土地利用侵蝕面積差異率：

C=A-Q；
D=C/Q。

式中：C為不同侵蝕強度或土地利用侵蝕面積差值，km2；A為非基準方法不同土壤侵蝕強度或土壤侵蝕面積，不同方法以下標進行區(qū)分，km2；Q為基準方法不同土壤侵蝕強度或土壤侵蝕面積，不同方法以下標進行區(qū)分，km2；D為侵蝕強度或不同土地利用侵蝕面積差異率，%。

計算量從計算所需數(shù)據(jù)量大小MB和運算時間長短min反映。圖斑大小采用圖斑數(shù)量、圖斑面積標準差、單位面積圖斑數(shù)量等來表達[12]。成圖效果利用Fragstats 4.2軟件提取聚集度指數(shù)(aggregation index)[13]，分析水力侵蝕面積及強度空間分布。聚集度指數(shù)取值越小，越接近于0，圖斑越離散。

3 結(jié)果與分析

3.1 基于柵格計算的3種統(tǒng)計方法土壤侵蝕差異

以動態(tài)監(jiān)測采用的柵格統(tǒng)計結(jié)果為基準，對比其他2種方法差異。從侵蝕強度(表1)看，地塊評價法無法生成不同侵蝕強度面積，且侵蝕面積差值和差異率均較大。而軟件判斷法侵蝕面積相差0.05 km2，差異率僅0.01%，且只在輕度和中度侵蝕面積略有差異。

表1 不同柵格統(tǒng)計方法不同侵蝕強度差異

從不同土地利用土壤侵蝕面積(表2)看，地塊評價法中，林地差異率最小，為0.48%；交通建設(shè)用地差異率最大，達121.9%。軟件判斷法中，園地、其他交通用地侵蝕面積沒有差異；草地差異率最大，為-2.17%，遠小于地塊評價法。

表2 不同柵格統(tǒng)計方法不同土地利用侵蝕面積差異

從計算量看，柵格統(tǒng)計法數(shù)據(jù)量小，運算耗時短，簡便快捷。地塊評價法、軟件判斷法數(shù)據(jù)量與計算步驟增加，計算及統(tǒng)計時間均有不同程度增加，地塊評價法增加時間較多。

從成圖效果(圖2)看，地塊評價法無法生成土壤侵蝕強度空間分布圖。柵格統(tǒng)計法以柵格成圖，聚集度指數(shù)88.89%，離散程度高，連續(xù)性差，空間分布相對分散，成圖效果較差。軟件判斷法以矢量成圖，聚集度指數(shù)96.81%，離散程度低，連續(xù)性好，空間分布相對集中，成圖效果較好。

圖2 柵格計算法土壤侵蝕空間分布Fig.2 Soil erosion distribution under different raster calculation methods

3.2 基于矢量的2種計算方法土壤侵蝕差異

以原始矢量法為基準，對比分析細化矢量法土壤侵蝕差異。從表3可見，輕度侵蝕面積差值最大，為-81.43 km2，但差異率最小，為-17.08%；劇烈侵蝕差值最小，為0.01 km2，但卻是從無到有的變化。

表3 不同矢量計算法不同侵蝕強度差異

從表4可見，細化矢量法不同土地利用侵蝕面積均小于原始矢量法。其中，建設(shè)用地差異率最小，為4.46%；交通運輸用地差異率最大，為50.33%。

表4 不同矢量計算法不同土地利用侵蝕面積差異

從圖斑大小(表5)可見，不同因素疊加后，圖斑數(shù)量增長7.27倍，單位面積數(shù)量增長7.47倍，而平均面積減小0.06 km2，面積標準差縮小7.25倍，圖斑尺寸減小，圖斑破碎程度增大。外井溝研究[9]發(fā)現(xiàn)，因素疊加后，56.93 km2的小流域圖斑數(shù)量增長3.12倍，單位面積數(shù)量擴大8.6倍，地塊平均面積減小0.05 km2?？梢?，無論是小流域還是縣域尺度，細化矢量法對刻畫更為細致，但會使圖斑數(shù)量成倍增加，同時增加破碎、無效圖斑，增加后續(xù)計算量，之后研究中可探索應(yīng)用微流域來表達地形對圖斑細化的影響。

表5 細化矢量法圖斑信息變化

從計算量看，細化后圖斑矢量數(shù)據(jù)增大至337 MB，計算用時增長為原來的1.5 倍，統(tǒng)計時間增加1 min。

如圖3所示，原始矢量法、細化矢量法聚集度指數(shù)分別為93.65%、96.81%，細化矢量法離散程度較低，成圖效果較好。同時，細化矢量法對不同侵蝕強度刻畫更細致，改善原始矢量法對中度及以上強度侵蝕刻畫較弱的情況。

圖3 矢量計算法土壤侵蝕空間分布Fig.3 Soil erosion distribution under different vector calculation methods

3.3 基于柵格計算與矢量計算方法差異

由于地塊評價法無法描述各侵蝕強度面積，筆者對其余4種方法侵蝕結(jié)果進行方差分析。結(jié)果(表6)可見，P<0.01，結(jié)果差異具有高度統(tǒng)計意義。其方差也顯示出各方法計算結(jié)果具有較大的差異性，依據(jù)不同情況選擇不同方法會對最終結(jié)果造成影響。

軟件判斷法與柵格統(tǒng)計法在不同侵蝕強度和土地利用侵蝕面積差異都較小。因此，以動態(tài)監(jiān)測采用的柵格統(tǒng)計法為基準，分析矢量計算法土壤侵蝕差異。

表6 不同方法侵蝕強度結(jié)果方差分析

如表7所示，原始矢量法劇烈侵蝕差異率最大，為-100%；輕度侵蝕差異率最小，為-45.79%。細化矢量法劇烈侵蝕無差異，強烈侵蝕差異率最大，為-278.16%。

表7 柵格統(tǒng)計與矢量計算法不同侵蝕強度差異

如表8所示，原始矢量法、細化矢量法交通運輸用地差異率均最大，分別為-2 539.42%、1 210.95%，建設(shè)用地差異率均最小，分別為-5.97%、1.25%。

表8 柵格統(tǒng)計與矢量計算法不同土地利用侵蝕差異

小流域尺度[9]上，細化矢量法與柵格統(tǒng)計法侵蝕面積差異率為8.47%；侵蝕強度面積上，強烈侵蝕無差異，中度和極強烈侵蝕相差0.01 km2，輕度差值較大也僅為2.66 km2。但在縣域尺度上，計算結(jié)果差值和差異率都較大，較難實現(xiàn)準確的土壤侵蝕評價。

從計算量看，柵格數(shù)據(jù)計算明顯快于矢量數(shù)據(jù)。柵格統(tǒng)計法不涉及矢量數(shù)據(jù)，歷時最短；細化矢量法數(shù)據(jù)量最大，計算、統(tǒng)計歷時最長；軟件判斷法則介于兩者之間。

如圖2和圖3所示，與柵格統(tǒng)計法相比，矢量計算法聚集度指數(shù)較高，連續(xù)性優(yōu)，離散程度較低，成圖效果較好，而與軟件判斷法基本一致。

4 討論

1)空間數(shù)據(jù)結(jié)構(gòu)對土壤侵蝕計算結(jié)果有一定的影響。柵格計算時，CSLE模型7因子均為10 m分辨率柵格，但E、T因子是采用中心屬性值法(rule of centric cell)柵格化獲得，此過程是一個有損轉(zhuǎn)換過程，且精度損失與柵格大小呈正相關(guān)[14]。可研究采用最大面積值法(rule of maximum area)[15]柵格化，探索不同方法對因子精度的影響，選擇最優(yōu)轉(zhuǎn)化方法。在矢量計算時，采取分區(qū)統(tǒng)計平均值(mean)獲取圖斑因子值，概化了地塊內(nèi)部因素，可采用眾數(shù)(majority)和中值(median)獲取因子值，分析不同方式下結(jié)果差異。

2)圖斑尺寸對矢量計算過程和結(jié)果產(chǎn)生影響。在圖斑解譯過程中，最小圖斑以400 m2控制，但山地丘陵區(qū)林地、耕地、其他交通用地等連續(xù)分布，解譯出超5 km2大圖斑17個。在原始矢量法獲取因子值時，圖斑尺寸越大，概化信息越多，結(jié)果差異越大。而細化矢量法計算過程中，因圖斑內(nèi)部因素復(fù)雜，此17個圖斑細化后數(shù)量增至1萬7 100個，產(chǎn)生大量無效圖斑，實際操作過程中通常需要進行小圖斑綜合處理[16]，增加數(shù)據(jù)量和工作量。因此，在圖斑解譯過程中，應(yīng)探索控制大圖斑面積和細化圖斑方法，減少大圖斑對計算的影響。

3)參考SL 190—2007《土壤侵蝕分類分級標準》[11]，選取81塊經(jīng)過軟件判斷法生成的侵蝕圖斑進行野外侵蝕強度復(fù)核。經(jīng)實地調(diào)查，20塊圖斑侵蝕強度偏高或偏低、地塊邊界不準確，以耕地、林地為主，準確率75.31%。錯誤的圖斑面積均較大，內(nèi)部坡度、坡長多數(shù)變化較大，覆蓋度等因素復(fù)雜交錯，在實際土地利用圖斑勾繪過程中，可根據(jù)遙感影像勾繪完整坡面的圖斑，即可減輕坡度坡長變化帶來的影響，同時又可控制圖斑面積，提高判斷準確率。水土流失圖斑落地與判斷軟件也可依據(jù)不同土地利用方式的不同覆蓋度特點，探索不同土地利用修正系數(shù)，減輕覆蓋度復(fù)雜變化帶來的影響。

5 結(jié)論

1)5種計算、統(tǒng)計方法的土壤侵蝕結(jié)果均有所差異。軟件判斷法和柵格統(tǒng)計法土壤侵蝕差異較小，而地塊評價方法、原始矢量法、細化矢量法較柵格統(tǒng)計法差異較大。

2)基于柵格計算的3種統(tǒng)計方法，柵格統(tǒng)計法數(shù)據(jù)量小，運算時間短，但成圖效果差，較難實現(xiàn)地塊化管理，適用于快速、準確掌握土壤侵蝕面積、強度和不同土地利用土壤侵蝕面積而無需制圖的縣域尺度評價。軟件判斷法較柵格統(tǒng)計法侵蝕差異最小，運算時間適中，成圖效果佳，又可實現(xiàn)侵蝕矢量化表達，適合縣域尺度土壤侵蝕評價與制圖。地塊評價法無法生成不同侵蝕強度面積，且較柵格統(tǒng)計法侵蝕面積和不同土地利用侵蝕面積差異較大，不適合縣域尺度土壤侵蝕評價。

3)基于矢量的2種計算方法，數(shù)據(jù)量大，運算時間長，成圖效果較好，但土壤侵蝕面積、強度及不同土地利用侵蝕面積較柵格計算法差異較大，無法準確反映土壤侵蝕情況，適用于土壤類型、坡度、植被覆蓋度相對一致的平原區(qū)土壤侵蝕評價。