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

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

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

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

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

1 材料與方法

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

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

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

1.2 空間插值法

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

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

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

1.3 地圖代數(shù)法

在1.2節(jié)中由于采用了CSLE模型計(jì)算抽樣單元的土壤水蝕速率，為了便于對比，在地圖代數(shù)法中也采用CSLE模型計(jì)算研究區(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ù)項(xiàng)目組提供的大范圍降雨侵蝕力R、中國范圍土壤可蝕性因子K，通過裁剪處理，并重采樣為1 km分辨率柵格數(shù)據(jù)，獲得研究區(qū)的降雨侵蝕力因子R專題層和土壤可蝕性因子K專題層；將項(xiàng)目組提供的多幅地形因子基礎(chǔ)數(shù)據(jù)進(jìn)行拼接并裁剪出研究區(qū)域，再根據(jù)水利部水土保持監(jiān)測中心編寫的《水土流失普查技術(shù)規(guī)定》計(jì)算坡度因子上限閾值(坡度30°時(shí)坡度因子值，S上限=10.0)和坡長因子上限閾值(坡長100 m時(shí)坡長因子值，L上限=2.1)，進(jìn)而計(jì)算得出LS上限閾值(LS上限= 21.2)，并重采樣為1 km分辨率柵格數(shù)據(jù)，最終獲得研究區(qū)LS數(shù)據(jù)專題層。

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

草地：

(2)

灌木林地：

(3)

喬木林地：

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

(4)

(5)

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

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

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

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

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

2 結(jié)果與分析

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

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

表2 基于2種方法的青藏高原典型樣區(qū)土壤水蝕速率的基本統(tǒng)計(jì)結(jié)果比較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可以看出，從土壤水蝕速率平均值來看，空間插值法的結(jié)果(593.5 t/(km2·年))與地圖代數(shù)法(549.4 t/(km2·年))較為接近，表明2種方法計(jì)算的研究區(qū)土壤侵蝕狀況基本一致。標(biāo)準(zhǔn)差可以衡量計(jì)算結(jié)果的離散程度。從土壤水蝕速率標(biāo)準(zhǔn)差來看，地圖代數(shù)法的結(jié)果是空間插值法的1.7倍，說明空間插值法的結(jié)果更平滑，更能反映研究區(qū)土壤侵蝕的整體狀況。從土壤水蝕速率極大值來看，地圖代數(shù)法的結(jié)果(5 000 t/(km2·年))是空間插值法(2 800 t/(km2·年))的1.8倍，所以地圖代數(shù)法能更好地反映研究區(qū)土壤侵蝕的極大值狀況。

為了進(jìn)一步比較2種方法計(jì)算結(jié)果的差別，以土壤水蝕速率計(jì)算結(jié)果為基礎(chǔ)，設(shè)定25為步長，統(tǒng)計(jì)各步長頻率，繪制頻率曲線和累積頻率曲線，結(jié)果如圖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ù)法的結(jié)果為單峰曲線，空間插值法的結(jié)果為雙峰曲線。從土壤水蝕速率累積頻率曲線來看，在地圖代數(shù)法的結(jié)果中，土壤水蝕速率小值累積頻率百分比較大且土壤水蝕速率大值累積頻率曲線更長，可知地圖代數(shù)法結(jié)果更突出土壤水蝕速率小值和大值，空間插值法結(jié)果更突出均值附近土壤水蝕速率值。出現(xiàn)這種狀況的原因如下：(1)在地圖代數(shù)法的的工程措施因子E計(jì)算過程中，以一個(gè)抽樣單元的E均值作為其周圍0.25緯度×0.25經(jīng)度范圍的E因子值，由于0.25緯度×0.25經(jīng)度范圍內(nèi)的E因子值不是均勻分布，這導(dǎo)致0.25緯度×0.25經(jīng)度范圍內(nèi)的E因子值要么被縮小要么被放大，進(jìn)而使該范圍內(nèi)土壤水蝕速率小值和土壤水蝕速率大值的百分比均增大。(2)在空間插值法的計(jì)算過程中，以抽樣單元的土壤水蝕速率均值為輸入數(shù)據(jù)，這使得空間插值法的結(jié)果主要反映每個(gè)抽樣單元的均值狀況，不能反映出每個(gè)抽樣單元土壤水蝕速率最小值和最大值狀況。所以地圖代數(shù)法的結(jié)果中土壤水蝕速率小值和大值所占比例均大于空間插值法結(jié)果，空間插值法所得的結(jié)果中土壤水蝕速率在均值附近出現(xiàn)第二個(gè)波峰。由于地圖代數(shù)法的結(jié)果中土壤水蝕速率小值占比較大，同時(shí)空間插值法的結(jié)果側(cè)重反映抽樣單元土壤水蝕速率均值狀況，抽樣單元的土壤水蝕速率最小值和最大值均未參與插值運(yùn)算，從而導(dǎo)致當(dāng)土壤水蝕速率較小時(shí)，地圖代數(shù)法累積頻率高于空間插值法；當(dāng)土壤水蝕速率較大時(shí)，地圖代數(shù)法的累積頻率曲線更長。

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

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

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

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

圖2顯示，從微觀來看，空間插值法的結(jié)果較為平滑，土壤水蝕速率較大值(>3 000 t/(km2·年))的分布較為集中；而地圖代數(shù)法的結(jié)果較為粗糙，土壤水蝕速率較大值(>3 000 t/(km2·年))的分布較為分散。以2種結(jié)果為基礎(chǔ)，計(jì)算局部方差(local variance，LV)，判斷2種結(jié)果的局部變異情況[30]。LV平均值越大，表明數(shù)據(jù)表面的局部變化越大、表面越粗糙。經(jīng)計(jì)算，地圖代數(shù)法的LV均值(696 138.3)大于空間插值法(107 957.6)，說明地圖代數(shù)法的結(jié)果能較好地反映局部區(qū)域變化情況。這是因?yàn)樵诘貓D代數(shù)法中，獲得準(zhǔn)確的面狀水保措施較為困難，難以表達(dá)真實(shí)的土壤水蝕速率；但該方法采用的是逐柵格單元運(yùn)算，每個(gè)柵格單元都具有侵蝕因子值，所以地圖代數(shù)法的計(jì)算結(jié)果能夠更好地反映局部區(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種方法計(jì)算的土壤侵蝕強(qiáng)度面積所占比例的對比

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

圖3 基于2種方法的青藏高原典型樣區(qū)不同土壤水蝕速率等級(jí)面積所占比例的對比 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ù)法的結(jié)果均能反映研究區(qū)的土壤侵蝕強(qiáng)度級(jí)別狀況，其中空間插值法的結(jié)果相對更合理。

2.4 本研究結(jié)果與其他研究結(jié)果的對比

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

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

3 討 論

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

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

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

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

4 結(jié) 論

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

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