陳羽璇， 楊勤科,2?， 劉寶元， 黃晨璐， 王春梅,2， 龐國偉,2

(1.西北大學城市與環(huán)境學院，710127，西安；2.旱區(qū)生態(tài)水文與災(zāi)害防治國家林業(yè)局重點實驗室，710127，西安；3.中國科學院 水利部 水土保持研究所，712100，陜西楊凌；4.北京師范大學地理學院，100875，北京)

土壤侵蝕已嚴重威脅世界許多地區(qū)的社會經(jīng)濟發(fā)展[1]，對其治理受到廣泛關(guān)注[2]，土壤侵蝕的治理需以區(qū)域土壤侵蝕調(diào)查制圖為數(shù)據(jù)基礎(chǔ)。自20世紀90年代以來，國內(nèi)外學者對區(qū)域土壤侵蝕評價進行一系列研究[2]。Lu等[3]和Teng等[4]先后用RUSLE對澳洲土壤侵蝕進行評價，Borrelli等[5]在全球尺度上進行土壤侵蝕定量分析，劉寶元等[6]基于抽樣調(diào)查方法進行全國土壤水力侵蝕普查制圖，楊志成等[7]在抽樣調(diào)查基礎(chǔ)上進行縣域和省域土壤侵蝕評價。目前已初步形成對坡面和小流域尺度土壤侵蝕速率的精度評價方法和對區(qū)域尺度侵蝕危險性評價方法，但區(qū)域尺度上實現(xiàn)土壤侵蝕速率的調(diào)查和制圖方法有待研究。

珠江是中國第三大河流，珠江三角洲是我國乃至世界最具活力的經(jīng)濟區(qū)之一。由于地處南亞熱帶，降水量大、土層薄，加上人類活動干擾，珠江流域土壤侵蝕愈發(fā)嚴峻[8-9]。自1990年以來，研究者對珠江流域土壤侵蝕及石漠化等進行一些研究[9-13]，但目前的研究對流域土壤侵蝕空間分異特征和主控因子的認識不足，與珠江三角洲繁榮的經(jīng)濟不相適應(yīng)，不能有效支撐珠江流域水土保持規(guī)劃和決策。為此，迫切需要利用現(xiàn)代科學技術(shù)手段，對該區(qū)土壤侵蝕進行系統(tǒng)研究。

筆者基于抽樣調(diào)查單元獲取的相關(guān)信息和區(qū)域土壤侵蝕因子數(shù)據(jù)，結(jié)合CSLE模型[14]與GIS技術(shù)，采用地圖代數(shù)和空間插值2種方法對研究區(qū)土壤侵蝕速率制圖并進行土壤侵蝕分析，以期為該地區(qū)生態(tài)環(huán)境治理及水土保持工作實施提供可靠的科學依據(jù)，同時促進我國亞熱帶地區(qū)土壤侵蝕的研究。

1 研究區(qū)概況

研究區(qū)包括珠江流域及海南省等相鄰區(qū)(圖1)，總面積約59.14萬km2。珠江流域地貌類型復雜，山地、丘陵占總面積94.4%，平原面積小且分散；西江中上游地區(qū)以石灰?guī)r山地為主，東江流域以花崗巖丘陵為主，北江中上游地區(qū)干流西側(cè)為石灰?guī)r山地、東側(cè)為紅巖盆地、下游以花崗巖丘陵為主。流域地處南亞熱帶，年平均氣溫在14～22 ℃之間；流域內(nèi)年平均降水量1 470 mm，自然植被以常綠闊葉林為主，主要土地利用類型為林地、草地和耕地。主要土壤類型有赤紅壤、紅壤、黃壤和水稻土等。

圖1 研究區(qū)位置圖Fig.1 Location map of the study area

本研究的數(shù)據(jù)主要包括2部分，一是土壤侵蝕抽樣調(diào)查數(shù)據(jù)，研究區(qū)(包括陸地部分55 km緩沖區(qū))共布設(shè)抽樣調(diào)查單元1 695個(圖2)，該數(shù)據(jù)基于高分辨率影像解譯而來，主要包含土地利用類型和水土保持措施2部分解譯內(nèi)容，精度優(yōu)于1∶5 000[15]；二是用于計算區(qū)域土壤侵蝕因子專題數(shù)據(jù)。數(shù)據(jù)及其來源詳見表1。

圖2 抽樣單元分布圖Fig.2 Distribution of sampling units

表1 基礎(chǔ)數(shù)據(jù)一覽表Tab.1 List of basic data

2 研究方法

2.1 土壤侵蝕因子計算

筆者使用CSLE模型[14]計算土壤侵蝕速率。所用參數(shù)為全球降雨侵蝕力因子(R因子)、中國土壤可蝕性因子(K因子)、全球L坡長因子與S坡度因子(LS因子)。這3個因子的計算方法詳見第一次全國水利普查[17]和相關(guān)文獻，本文主要介紹水土保持措施因子的計算。

植被覆蓋與生物措施因子(B因子)：基于Borrelli等的[5]C因子計算方法，根據(jù)研究區(qū)實際情況改進。對農(nóng)業(yè)用地C因子賦值為1；對非農(nóng)業(yè)用地C因子，結(jié)合GLC30土地利用分類，定義相應(yīng)土地利用C因子值的約束范圍。

CNA=MINCF+(MAXCF-MINCF)(1-TC)

(1)

式中：CNA為林地、灌木覆蓋C因子值；MINCF=0.000 1、MAXCF=0.003；TC為林地覆蓋的比例。

CP=MINC+(MAXC-MINC)·(1-NVS)

(2)

式中：Cp為非林地覆蓋(草地覆蓋等)C因子值；MINC=0.01、MAXC=0.05；NVS為非林地(草地等)覆蓋的比例。

裸土按照式(2)計算，MINC、MAXC取值為0.1、0.5，其計算結(jié)果分辨率為250 m，通過重采樣得到30 m分辨率數(shù)據(jù)。最后基于GLC30土地利用分類合并CNA、Cp、裸土C和耕地C，得到B因子圖層。

工程措施因子(E因子)：利用抽樣單元土地利用數(shù)據(jù)，根據(jù)《區(qū)域水土流失動態(tài)監(jiān)測技術(shù)規(guī)定》中計算方法編寫程序計算得到抽樣單元E因子均值，再對E因子均值的點數(shù)據(jù)進行反距離權(quán)重空間插值得到E因子圖層。

耕作措施因子(T因子)：根據(jù)《區(qū)域水土流失動態(tài)監(jiān)測技術(shù)規(guī)定》，對抽樣單元內(nèi)耕地的輪作措施屬性字段賦值，非耕地賦值1，并基于Borrelli等[5]對農(nóng)業(yè)中水田的處理方法，將水田賦值為0.15，最終得到T因子圖層。

2.2 流域土壤侵蝕速率圖制作

筆者采用2種方法：空間插值法和地圖代數(shù)法。空間插值法是基于抽樣單元的侵蝕速率圖[15]，用地統(tǒng)計協(xié)同克里金插值方法，以R、K、B因子作為協(xié)變量，對耕、林、草、灌木及裸地分別插值生成不同地類圖，并將水體、濕地、不透水面設(shè)置為0，最終得到土壤侵蝕速率專題圖。地圖代數(shù)法是根據(jù)準備好的各因子圖層數(shù)據(jù)，直接進行各因子相乘得到土壤侵蝕圖，考慮到LS、K因子分辨率均為30 m，本研究的運算將分辨率設(shè)置為30 m。

2.3 主控性因子分析方法

基于地理探測器進行土壤侵蝕主控因子分析。地理探測器是一種通過探測要素的空間分層異質(zhì)性及揭示背后驅(qū)動力的一種統(tǒng)計學方法[18]，其大小由q來衡量，q表達式如下：

(3)

(4)

3 結(jié)果與分析

3.1 土壤侵蝕影響因子分布特征

影響土壤侵蝕速率各因子中(圖3)，R因子總體上呈從東至西遞減的規(guī)律，其值域為1 786.81～16 324.18 MJ·mm/(hm2·h·a)，平均值為6 813.56 MJ·mm/(hm2·h·a)。K因子值域為0～0.012 5 t·hm2·h/(MJ·mm·hm2)，均值為0.004 8 t·hm2·h/(MJ·mm·hm2)，有從東到西遞增趨勢。LS因子值域為0.04～21.24，平均值為8.70，高值見于云貴高原、東部山地和海南五指山等地。R、LS因子值普遍較高，是土壤侵蝕的主要誘發(fā)因子。B因子、E因子和T因子，值域分布范圍依次為0～1、0.15～1、0.15～1，其均值依次為0.21、0.91、0.87。B因子值較低，而E和T因子值較高，研究區(qū)植被覆蓋度較好，是土壤侵蝕主要抑制因子。

圖3 CSLE各因子空間分布圖Fig.3 Spatial distribution of CSLE factors

3.2 2種計算方法下土壤侵蝕的對比分析

3.2.1 空間分布特征對比分析 在水利部批準SL190—2007《土壤侵蝕分類分級標準》[19]的基礎(chǔ)上進行細分，編制了土壤侵蝕分級圖(圖4和圖5)，研究區(qū)土壤侵蝕均主要集中分布在研究內(nèi)的貴州省及云南省、廣西中部和廣東省沿海區(qū)，其中強烈和極強烈侵蝕區(qū)域范圍較小且分布在較零散的坡耕地上。2張圖上土壤侵蝕速率較高的地區(qū)與土地利用圖的耕地分布較吻合，而林、草地的侵蝕則相對較弱。從微觀角度看，2圖在細節(jié)紋理結(jié)構(gòu)上略有不同。基于地圖代數(shù)法的結(jié)果局地變化比較明顯，空間插值結(jié)果整體較為平滑，主要原因在于插值所用的侵蝕速率為抽樣單元均值，且插值過程中有平滑作用。

圖4 地圖代數(shù)法土壤侵蝕速率圖Fig.4 Soil erosion rate map based on map algebra method

圖5 空間插值法土壤侵蝕速率圖Fig.5 Soil erosion rate map based on spatial interpolation method

3.2.2 統(tǒng)計特征對比分析 采用地圖代數(shù)和空間插值法制圖的土壤侵蝕速率均值分別為791.78 t/(km2·a)和615.37 t/(km2·a)，兩者值較接近。地圖代數(shù)法計算的土壤侵蝕速率均值較高，可能與水保措施難以全面納入計算有關(guān)。筆者計算的土壤侵蝕速率值與前人研究[20-22]對比(表2)，比較接近，可以認為本研究的結(jié)果較為可信。對比研究區(qū)2種方法的主要土地利用類型的土壤侵蝕速率結(jié)果(表3)發(fā)現(xiàn)，土壤侵蝕速率值較為接近?？臻g插值法制圖結(jié)果，耕地、林地的土壤侵蝕速率均值略小，這是由于地圖代數(shù)法是從整個研究區(qū)的計算，分辨率相對較粗，難以獲取詳細的水土保持措施，空間插值法從抽樣單元小流域可計算得到準確的水土保持措施值，計算的土壤侵蝕速率值較小。

表2 土壤侵蝕速率對比參照表Tab.2 Soil erosion rates of this study and extracted from references t/(km2·a)

表3 主要土地利用土壤侵蝕速率對比表Tab.3 Average soil erosion rates of land use types t/(km2·a)

3.2.3 土壤侵蝕強度對比分析 根據(jù)SL190—2007《土壤侵蝕分類分級標準》[19]對土壤侵蝕強度進行統(tǒng)計(表4)，2種方法計算的水土流失面積(輕度及其以上侵蝕)差值為10.17%，空間插值法略大；土壤侵蝕均主要集中于微、輕度侵蝕；但其在微、輕度侵蝕頻率分布略有差異，這可能與空間插值時耕地土壤侵蝕速率值較集中有關(guān)；在中度及以上的侵蝕強度上，土壤侵蝕速率結(jié)果頻率相差不大。

表4 空間插值法與地圖代數(shù)法土壤侵蝕強度統(tǒng)計Tab.4 Soil erosion intensity statistics by spatialinterpolation and map algebra %

總體上看，2種方法計算的土壤侵蝕速率結(jié)果均可表達研究區(qū)土壤侵蝕狀況，考慮到土壤侵蝕防控的需要，優(yōu)先選用空間插值法。

3.3 土壤侵蝕主控因子分析

基于地理探測器的分析結(jié)果(表5)表明，整個研究區(qū)和水土保持區(qū)劃的各個區(qū)，各因子影響程度從大到小的排序基本為土地利用、水保措施(特別是B因子)、地形、土壤和降水。除浙閔山地丘陵區(qū)外，土地利用方式為土壤侵蝕的主控因子，影響因子q值均達到47%以上；其次，B因子對土壤侵蝕速率及其空間分布的影響較大；T因子的影響排第三，華南沿海丘陵臺地區(qū)和海南及南海諸島丘陵區(qū)的q值達到40%以上，可能與該區(qū)復雜的耕作方式有關(guān)，其原因有待進一步探究分析。K因子和R因子的影響較小，這與研究區(qū)以林草植被為主、植被覆蓋普遍較高有關(guān)。

表5 各影響因子q值統(tǒng)計表Tab.5 Statistical table of q values of each impact factor

4 結(jié)論

1)研究區(qū)土壤侵蝕特征環(huán)境是降雨侵蝕力較大、坡度較陡及植被覆蓋度較高，前2項是土壤侵蝕的主要誘發(fā)因子，后者則是土壤侵蝕的主要抑制因子。

2)本文通過地圖代數(shù)和空間插值2種方法制圖得到土壤侵蝕速率均值分別為791.78 t/(km2·a)、615.37 t/(km2·a)，研究區(qū)土壤侵蝕廣泛分布，坡耕地是土壤侵蝕主要的發(fā)生部位。

3)土地利用類型為影響土壤侵蝕速率的主控因素，B因子對土壤侵蝕速率及其空間分布影響較大，總體上講，由于研究區(qū)植被覆蓋度普遍較高，因而土壤和降雨的影響較小。調(diào)整土地利用結(jié)構(gòu)、不斷優(yōu)化植被的水土保持生態(tài)功能是值得注意的問題。

4)地圖代數(shù)和空間插值兩種方法制圖結(jié)果相比，空間插值法的結(jié)果更接近土壤侵蝕的實際情況，可作為區(qū)域土壤侵蝕制圖首選方法，但兩者之間的差異和相互結(jié)合有待討論。

中國科學院水利部水土保持研究所張曉萍研究員、水利部水土保持監(jiān)測中心王愛娟高工、珠江水利委員會水土保持監(jiān)測中心金平偉主任對本研究給予多方支持，一并致謝。