杜朝正,楊勤科,王春梅,龐國偉

(1 西北大學 城市與環(huán)境學院, 陜西 西安 710127;2 菏澤學院 城市建設學院, 山東 菏澤 274015)

土壤侵蝕是全球性生態(tài)環(huán)境問題之一[1]，被侵蝕的土壤有機碳是溫室氣體碳成分的重要來源之一[2]，土壤侵蝕速率大小直接影響水土流失區(qū)域的土壤侵蝕經(jīng)濟損失狀況[3]。治理土壤侵蝕必須以土壤侵蝕定量評價和制圖為基礎[4-6]，即通常采用土壤侵蝕模型定量計算土壤侵蝕速率[7]。逐步建立的土壤侵蝕模型有美國的WEPP模型[8]、USLE[9]和RUSLE[10]，歐洲的EROSEM模型[11]，荷蘭的LISEM模型[12]以及中國的CSLE[13]等。國外較多研究采用USLE/RUSLE計算區(qū)域土壤水蝕速率，如Panagos等[14]編制歐盟土壤流失圖，Teng等[15]編制澳大利亞土壤侵蝕圖。國內(nèi)學者多采用CSLE計算區(qū)域土壤水蝕速率，如章文波等[16]采用CSLE預測坡面和田間土壤侵蝕，程琳等[17]采用CSLE對陜西省土壤侵蝕進行定量評價。

一般進行區(qū)域土壤侵蝕定量評價方法可歸納為2種：第一種是基于侵蝕因子圖層整體運算的方法，可稱為地圖代數(shù)法，代表性研究有程琳等[17]、Borrelli等[18]、Yang等[19]、周來等[20]；第二種是基于抽樣調(diào)查結果和插值的方法，在研究區(qū)布設基本抽樣調(diào)查單元，結合土壤侵蝕模型計算抽樣單元土壤侵蝕速率，使用空間插值模型來完成區(qū)域土壤侵蝕評價制圖，簡稱為空間插值法，代表性研究有劉寶元等[21]、Yin等[22]。地圖代數(shù)法和空間插值法均已被廣泛應用，但這2種方法有何異同還鮮有報道。

本研究以青藏高原典型樣區(qū)為例，以現(xiàn)有土壤侵蝕因子數(shù)據(jù)和抽樣調(diào)查單元數(shù)據(jù)為數(shù)據(jù)源，分別采用基于侵蝕因子柵格圖層的地圖代數(shù)法和基于抽樣調(diào)查單元的空間插值法，計算研究區(qū)土壤水蝕速率，并對計算結果進行對比分析，以期為土壤侵蝕定量評價與制圖提供方法參考。

1 材料與方法

1.1 研究區(qū)域與基礎數(shù)據(jù)

研究區(qū)位于青藏高原典型樣區(qū)，總面積170萬km2，海拔為331～7 182 m，平均海拔超過4 000 m[23]，年均降水量分布不均，由西北向東南呈逐漸遞增的形式分布，是多條重要河流的水源地，也被稱為“亞洲水塔”，為超過14億人提供水源。目前該區(qū)域土壤侵蝕研究相對滯后[24]，迫切需要對土壤侵蝕進行試驗觀測和制圖[25]。本研究所用青藏高原典型樣區(qū)的基礎數(shù)據(jù)見表1。

表1 青藏高原典型樣區(qū)的基礎數(shù)據(jù)Table 1 Basic data of typical sample areas of the Tibetan Plateau

1.2 空間插值法

1.2.1 抽樣單元布設與土壤水蝕速率的計算 在研究區(qū)外圍創(chuàng)建55 km緩沖區(qū)，在研究區(qū)和緩沖區(qū)內(nèi)，按照0.25緯度×0.25經(jīng)度布設樣點，共計2 631個抽樣單元，再根據(jù)《泛第三極土壤侵蝕抽樣調(diào)查單元遙感解譯工作大綱》，利用高分辨率遙感數(shù)據(jù)，解譯抽樣單元內(nèi)土地利用類型和水保措施信息，用于支持抽樣單元植被覆蓋與生物措施因子(B)、工程措施因子(E)、耕作措施因子(T)的計算。再結合降雨侵蝕力因子(R)、土壤可蝕性因子(K)、坡度坡長因子(LS)等數(shù)據(jù)，根據(jù)水利部水土保持監(jiān)測中心編寫的《水土流失普查技術規(guī)定》和文獻[26]，采用CSLE模型計算各抽樣單元的土壤水蝕速率。

1.2.2 空間插值計算 以計算的抽樣單元土壤水蝕速率平均值、GLC30中的土地利用類型數(shù)據(jù)為基礎，通過空間插值法計算研究區(qū)土壤水蝕速率，本研究采用地統(tǒng)計協(xié)同克里金插值，該插值方法能夠容納多個解釋變量，較好地反映多個環(huán)境因子對目標變量的影響[27]。與普通克里金方法相比，協(xié)同克里金插值可以輸入多個協(xié)變量作為影響因素，使區(qū)域土壤水蝕速率的計算結果更加準確。具體步驟為：(1)統(tǒng)一土地利用分類系統(tǒng)，使解譯土地利用類型與GLC30中的數(shù)據(jù)相對應；(2)按照一級土地利用類型，統(tǒng)計每個抽樣單元的平均土壤水蝕速率；(3)按照一級土地利用類型進行空間插值，依據(jù)統(tǒng)計結果，并以降雨侵蝕力因子(R)、土壤可蝕性因子(K)、地形因子(LS)為協(xié)變量，按照不同土地利用類型依次進行地統(tǒng)計協(xié)同克里金插值，獲得耕地、林地、草地、居住及建設用地、裸土等的插值面；(4)空間插值后處理：按照GLC30中的一級土地利用類型分布，將各土地利用類型插值結果圖層進行合成，最終獲得區(qū)域土壤水蝕速率柵格數(shù)據(jù)層，其中湖泊、坑塘、永久性冰川積雪、戈壁、裸巖等賦值為0。

為了確保計算精度及方便對比結果，本研究插值結果采用WGS84坐標系，區(qū)域水蝕速率專題層分辨率設定為1 km。

1.3 地圖代數(shù)法

在1.2節(jié)中由于采用了CSLE模型計算抽樣單元的土壤水蝕速率，為了便于對比，在地圖代數(shù)法中也采用CSLE模型計算研究區(qū)土壤水蝕速率專題層，其基本形式如下：

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

(1)

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

1.3.1 降雨侵蝕力因子R、土壤可蝕性因子K、地形因子LS數(shù)據(jù)處理[26]根據(jù)項目組提供的大范圍降雨侵蝕力R、中國范圍土壤可蝕性因子K，通過裁剪處理，并重采樣為1 km分辨率柵格數(shù)據(jù)，獲得研究區(qū)的降雨侵蝕力因子R專題層和土壤可蝕性因子K專題層；將項目組提供的多幅地形因子基礎數(shù)據(jù)進行拼接并裁剪出研究區(qū)域，再根據(jù)水利部水土保持監(jiān)測中心編寫的《水土流失普查技術規(guī)定》計算坡度因子上限閾值(坡度30°時坡度因子值，S上限=10.0)和坡長因子上限閾值(坡長100 m時坡長因子值，L上限=2.1)，進而計算得出LS上限閾值(LS上限= 21.2)，并重采樣為1 km分辨率柵格數(shù)據(jù)，最終獲得研究區(qū)LS數(shù)據(jù)專題層。

1.3.2 植被覆蓋與生物措施因子B計算 根據(jù)水利部水土保持監(jiān)測中心編寫的《水土流失普查技術規(guī)定》和文獻[26]中的方法計算B值。本研究中，取2014-2016年NPV均值作為喬木類林下地表植被蓋度數(shù)據(jù)。草地、灌木林地、喬木林地B值采用下式計算：

草地：

(2)

灌木林地：

(3)

喬木林地：

SLRi=0.444 68×e(-3.200 96×GD3)-0.040 99×
e(YBD-GD3×YBD)+0.025，

(4)

(5)

式中：SLRi是第i個半月植被覆蓋與生物措施因子值，無量綱；GD1、GD2、GD3分別為草地、灌木林地、喬木林地的地表植被蓋度，取值0～1；YBD是郁閉度，指喬木林地植被覆蓋度，取值0～1；WRi為第i個半月降雨侵蝕力占全年降雨侵蝕力的比例，取值0～1。

裸土地類B值：當24期半月植被蓋度數(shù)據(jù)中的第14期半月植被蓋度數(shù)據(jù)(the 14th fractional vegetation coverage，fvc14)>0.05時，B裸土=B草地；當fvc14≤0.05時，B裸土=0.5。耕地、居住及建設用地、湖泊、坑塘、永久性冰川積雪、戈壁、裸巖等土地類型B值按《水土流失普查技術規(guī)定》要求賦值計算。

1.3.3 工程措施因子E計算 依據(jù)抽樣單元布設原則，按照0.25緯度×0.25經(jīng)度間隔，繪制研究區(qū)內(nèi)矢量網(wǎng)格，確保每個網(wǎng)格內(nèi)只包括一個抽樣調(diào)查單元，將抽樣單元內(nèi)E的平均值作為矢量網(wǎng)格的E值，再將矢量數(shù)據(jù)轉換為1 km分辨率的柵格數(shù)據(jù)，獲得研究區(qū)的E因子專題圖。

1.3.4 耕作措施因子T計算 首先將劉巽浩等[28]編制的中國耕作制度區(qū)劃數(shù)據(jù)轉換為矢量圖形數(shù)據(jù)，并與《水土流失普查技術規(guī)定》中的輪作分區(qū)的T因子值表進行數(shù)據(jù)關聯(lián)，獲得全國耕作區(qū)劃T值圖，對其進行裁剪，并將矢量數(shù)據(jù)轉換為分辨率1 km的柵格數(shù)據(jù)，獲得研究區(qū)的區(qū)域耕作區(qū)劃T因子圖；然后依據(jù)GLC30數(shù)據(jù)，將非耕地區(qū)域的T值設置為1，將耕地區(qū)域的T值設置為耕作區(qū)劃T值，最終獲得研究區(qū)的T因子專題圖。

為了確保計算精度及方便對比結果，本研究計算采用WGS84坐標系，最終獲得的R、K、LS、B、E、T、A專題層分辨率均設定為1 km。

2 結果與分析

2.1 2種方法計算的土壤水蝕速率的統(tǒng)計特征對比

對2種方法計算的土壤水蝕速率的基本統(tǒng)計信息進行對比，以累積頻率為99%時的土壤水蝕速率為極大值進行統(tǒng)計，結果如表2所示。

表2 基于2種方法的青藏高原典型樣區(qū)土壤水蝕速率的基本統(tǒng)計結果比較Table 2 Comparison of basic statistical results of soil water erosion rates based on two methods in the typical sample areas of Tibet Plateau t/(km2·年)

從表2可以看出，從土壤水蝕速率平均值來看，空間插值法的結果(593.5 t/(km2·年))與地圖代數(shù)法(549.4 t/(km2·年))較為接近，表明2種方法計算的研究區(qū)土壤侵蝕狀況基本一致。標準差可以衡量計算結果的離散程度。從土壤水蝕速率標準差來看，地圖代數(shù)法的結果是空間插值法的1.7倍，說明空間插值法的結果更平滑，更能反映研究區(qū)土壤侵蝕的整體狀況。從土壤水蝕速率極大值來看，地圖代數(shù)法的結果(5 000 t/(km2·年))是空間插值法(2 800 t/(km2·年))的1.8倍，所以地圖代數(shù)法能更好地反映研究區(qū)土壤侵蝕的極大值狀況。

為了進一步比較2種方法計算結果的差別，以土壤水蝕速率計算結果為基礎，設定25為步長，統(tǒng)計各步長頻率，繪制頻率曲線和累積頻率曲線，結果如圖1所示。

圖1 基于2種方法的青藏高原典型樣區(qū)土壤水蝕速率頻率曲線(A)和累積頻率曲線(B)的對比Fig.1 Comparison of soil water erosion rate frequency curve (A) and cumulative frequency curve(B) based on two methods in the typical sample areas of Tibet Plateau

從圖1可以看出，地圖代數(shù)法與空間插值法的土壤水蝕速率頻率曲線和累積頻率曲線差異較大。從土壤水蝕速率頻率曲線來看，地圖代數(shù)法的結果為單峰曲線，空間插值法的結果為雙峰曲線。從土壤水蝕速率累積頻率曲線來看，在地圖代數(shù)法的結果中，土壤水蝕速率小值累積頻率百分比較大且土壤水蝕速率大值累積頻率曲線更長，可知地圖代數(shù)法結果更突出土壤水蝕速率小值和大值，空間插值法結果更突出均值附近土壤水蝕速率值。出現(xiàn)這種狀況的原因如下：(1)在地圖代數(shù)法的的工程措施因子E計算過程中，以一個抽樣單元的E均值作為其周圍0.25緯度×0.25經(jīng)度范圍的E因子值，由于0.25緯度×0.25經(jīng)度范圍內(nèi)的E因子值不是均勻分布，這導致0.25緯度×0.25經(jīng)度范圍內(nèi)的E因子值要么被縮小要么被放大，進而使該范圍內(nèi)土壤水蝕速率小值和土壤水蝕速率大值的百分比均增大。(2)在空間插值法的計算過程中，以抽樣單元的土壤水蝕速率均值為輸入數(shù)據(jù)，這使得空間插值法的結果主要反映每個抽樣單元的均值狀況，不能反映出每個抽樣單元土壤水蝕速率最小值和最大值狀況。所以地圖代數(shù)法的結果中土壤水蝕速率小值和大值所占比例均大于空間插值法結果，空間插值法所得的結果中土壤水蝕速率在均值附近出現(xiàn)第二個波峰。由于地圖代數(shù)法的結果中土壤水蝕速率小值占比較大，同時空間插值法的結果側重反映抽樣單元土壤水蝕速率均值狀況，抽樣單元的土壤水蝕速率最小值和最大值均未參與插值運算，從而導致當土壤水蝕速率較小時，地圖代數(shù)法累積頻率高于空間插值法；當土壤水蝕速率較大時，地圖代數(shù)法的累積頻率曲線更長。

以統(tǒng)計的頻率數(shù)據(jù)為基礎，計算直方圖交(histogram intersection，HI)[29]，判斷2種方法計算結果的直方圖相似性，一般HI取值范圍為0～1，HI值越大，說明對比的2種數(shù)據(jù)結果越相似。經(jīng)計算，2種方法計算結果的HI值為0.715，說明這2種方法計算結果具有較高的相似性，即2種方法計算結果反映的研究區(qū)土壤侵蝕狀況基本一致。

綜上可知，在青藏高原樣區(qū)，從統(tǒng)計特征方面來看，地圖代數(shù)法與空間插值法的計算結果均能反映研究區(qū)土壤侵蝕基本狀況，其中空間插值法的結果能更好地反映研究區(qū)土壤水蝕速率均值左右的狀況，且更貼近真實值。

2.2 2種方法計算的土壤水蝕速率表面特征的對比

采用地圖代數(shù)法和空間插值法2種方法獲得的研究區(qū)土壤水蝕速率表面特征分布結果如圖2所示。圖2顯示，從宏觀方面來看，2種方法計算結果均顯示由北向南研究區(qū)土壤水蝕速率呈逐漸增強的趨勢，分布趨勢與實際相符。這可能是因為北部降雨量較少、地勢相對平坦，南部降雨較多、且地形起伏比較大。土壤水蝕速率較高的區(qū)域大多分布于河流兩側，水動力的搬運能力較強，較易產(chǎn)生水力侵蝕。

圖2顯示，從微觀來看，空間插值法的結果較為平滑，土壤水蝕速率較大值(>3 000 t/(km2·年))的分布較為集中；而地圖代數(shù)法的結果較為粗糙，土壤水蝕速率較大值(>3 000 t/(km2·年))的分布較為分散。以2種結果為基礎，計算局部方差(local variance，LV)，判斷2種結果的局部變異情況[30]。LV平均值越大，表明數(shù)據(jù)表面的局部變化越大、表面越粗糙。經(jīng)計算，地圖代數(shù)法的LV均值(696 138.3)大于空間插值法(107 957.6)，說明地圖代數(shù)法的結果能較好地反映局部區(qū)域變化情況。這是因為在地圖代數(shù)法中，獲得準確的面狀水保措施較為困難，難以表達真實的土壤水蝕速率；但該方法采用的是逐柵格單元運算，每個柵格單元都具有侵蝕因子值，所以地圖代數(shù)法的計算結果能夠更好地反映局部區(qū)域水蝕速率變化情況。

圖2 基于2種方法的青藏高原典型樣區(qū)土壤水蝕速率的表面特征分布Fig.2 Surface characteristic distribution of soil water erosion rate based on two methods in the typical sample areas of Tibet Plateau

2.3 2種方法計算的土壤侵蝕強度面積所占比例的對比

對2種方法計算得土壤水蝕速率進行分級，并統(tǒng)計各等級侵蝕強度面積所占比例，結果如圖3所示。圖3顯示，2種方法計算結果中，當土壤水蝕速率為0～25 t/(km2·年)時，其面積所占比例最大，并且不同土壤水蝕速率等級面積所占比例的變化趨勢基本一致，說明2種方法計算結果反映的研究區(qū)土壤侵蝕強度級別狀況基本一致。當土壤水蝕速率為>250～3 000 t/(km2·年)時，空間插值法的面積比例均高于地圖代數(shù)法；當土壤水蝕速率≤250 t/(km2·年)或者>3 000 t/(km2·年)時，地圖代數(shù)法的面積比例均明顯高于空間插值法，說明空間插值法的結果更突出中間土壤侵蝕強度級別，而地圖代數(shù)法的結果更突出較小和較大土壤侵蝕強度級別。表2中，采用空間插值法和地圖代數(shù)法的所得的土壤水蝕速率平均值分別為593.5和549.4 t/(km2·年)，均位于500～750 t/(km2·年)侵蝕強度級別，而在此級別下空間插值法的面積比例高于地圖代數(shù)法，可知空間插值法的結果相對更合理。

圖3 基于2種方法的青藏高原典型樣區(qū)不同土壤水蝕速率等級面積所占比例的對比 Fig.3 Comparison of area ratios of different soil water erosion rate grades based on two methods in the typical sample areas of Tibet Plateau

綜上可知，空間插值法和地圖代數(shù)法的結果均能反映研究區(qū)的土壤侵蝕強度級別狀況，其中空間插值法的結果相對更合理。

2.4 本研究結果與其他研究結果的對比

Wang等[31]基于USLE和數(shù)學模型，綜合評價西藏地區(qū)土壤侵蝕敏感性，判斷西藏地區(qū)土壤侵蝕發(fā)生的可能程度，評價的趨勢結果與本研究計算結果相吻合。此外，已有部分學者對研究區(qū)的平均土壤水蝕速率進行了定量計算，故將本研究采用空間插值法和地圖代數(shù)法對不同研究區(qū)的平均土壤水蝕速率進行了計算，并將所得結果與前人研究結果(文獻值)進行了對比，結果如表3所示。

表3 基于本研究2種方法所得不同研究區(qū)土壤水蝕速率與前人結果的對比Table 3 Comparison of soil water erosion rate in different study areas between the two methods and results from previous studies

表3顯示，康琳琦等[32-33]采用USLE模型計算的青藏高原1980-2013年土壤水蝕速率為400～800 t/(km2·年)，與本研究區(qū)2種方法所得結果相近。Shao等[34]、李俊杰等[37]、宮奎方等[38]計算結果與本研究中2種方法計算結果也均相近，說明2種方法計算結果具有一定的可信度。Borrelli等[18]采用地圖代數(shù)法計算的青藏高原典型樣區(qū)潛在土壤水蝕速率均值為660.9 t/(km2·年)，本研究采用空間插值法和地圖代數(shù)法計算的土壤水蝕速率均值分別為593.5和549.4 t/(km2·年)，比Borrelli等的結果稍小。分析其主要原因是，雖然本研究采用的CSLE與Borrelli等[18]采用的USLE/RUSLE沒有本質差別，但Borrelli等的計算結果未考慮詳細的水保措施，計算結果僅僅為土壤水蝕潛在研究，并不是真實的土壤水蝕速率，本研究2種方法均考慮了詳細的水保措施，計算結果更接近真實值。本研究與文獻[35-36]關于納通河小流域計算結果差異較大，可能原因為文獻[35-36]中的結果為137Cs采樣點的平均值，未計算非采樣點的區(qū)域土壤水蝕速率；而本研究結果為區(qū)域的平均值，這兩種平均值之間可能存在一定差異。本研究與文獻[39]關于拉薩河流域計算結果差異較大，可能原因是文獻[39]中的水保措施較粗放，僅根據(jù)土地利用信息進行了主觀賦值，且賦值偏高，其中文獻[39]將鹽堿地、沙地、冰川/永久積雪等土地類型賦值為1，這與水利部水土保持監(jiān)測中心編寫的《水土流失普查技術規(guī)定》賦值為0差異較大。通過對比可知，本研究的2種方法計算結果較為合理，具有一定的可信度。

3 討 論

空間插值法與地圖代數(shù)法2種方法的差異在于對基礎數(shù)據(jù)的掌握程度以及對地面真實侵蝕狀況的估算程度不同。

土壤侵蝕發(fā)生在微觀尺度上，為了掌握土壤侵蝕真實值，必須在微觀尺度上計算土壤水蝕速率，為此美國農(nóng)業(yè)部的自然資源清查和中國水利部的水利普查(水土流失狀況普查)等工作發(fā)展了土壤侵蝕抽樣調(diào)查方法。為了準確獲取微觀尺度上的各侵蝕因子信息，特別是小流域尺度上的水土保持措施信息，必須進行抽樣單元布設并對高分辨率遙感影像解譯；在編制區(qū)域土壤侵蝕圖時，為了計算無采樣地區(qū)的土壤侵蝕速率，必須借助空間插值方法來實現(xiàn)[21-22]。受到方法自身的限制，空間插值法不一定很好地表達區(qū)域土壤侵蝕的局部變異；但該方法以抽樣單元真實水蝕速率為基礎，以降雨侵蝕力因子R、土壤可蝕性因子K、坡度坡長因子LS為協(xié)變量，計算的區(qū)域水蝕速率結果穩(wěn)定且貼近真實值。

地圖代數(shù)法是基于各侵蝕因子專題層進行整體運算，對于大范圍研究區(qū)而言，各侵蝕因子專題層的分辨率與土壤侵蝕發(fā)生的尺度可能不完全匹配，水保措施信息難以用專題層的方式精確表達。因此，在缺少足夠水保措施信息的情況下，地圖代數(shù)法計算的結果是潛在土壤侵蝕，難以完全反映區(qū)域土壤水蝕速率真實情況；但該方法是逐柵格單元運算，每個柵格單元都具有侵蝕因子值，所以，地圖代數(shù)法計算的區(qū)域水蝕速率結果能夠更好地反映局部變化情況。

雖然Borrelli等[18]已在全球范圍內(nèi)研究數(shù)據(jù)的分辨率可達250 m，但本研究側重計算方法的對比，考慮到本研究區(qū)大部分區(qū)域環(huán)境條件下土壤水蝕速率空間變異較小，暫時將計算結果分辨率設定為1 km。雖然本研究也得到了青藏高原典型樣區(qū)不同土壤水蝕速率計算結果，但需要更多的實測數(shù)據(jù)來驗證。以NPV數(shù)據(jù)作為林下蓋度計算林地水蝕速率，可能使林地水蝕速率有些偏差，如何提升林下蓋度計算精度，并結合更多的實測數(shù)據(jù)進行驗證，還需要進一步研究。

4 結 論

在青藏高原典型樣區(qū)，空間插值法和地圖代數(shù)法均可用于區(qū)域土壤侵蝕定量評價與制圖，2種方法計算結果均能反映研究區(qū)土壤水蝕速率總體分布狀況，2種方法計算結果的宏觀分布趨勢相同，平均值相近，直方圖相交值HI達0.715。地圖代數(shù)結果能較好地反映局部土壤侵蝕變異情況，但面域范圍內(nèi)獲得詳細的水保措施信息較為困難；空間插值法以抽樣單元真實水蝕速率為基礎，其結果能更好地反映研究區(qū)的平均土壤水蝕狀況，結果穩(wěn)定且接近真實值，故在青藏高原典型樣區(qū)可優(yōu)先選擇空間插值法。

志謝：本研究在劉寶元教授指導下完成，章文波老師和殷兵博士參與了計算過程的討論，一并表示感謝。