接壩區(qū)不同面積小流域泥沙來源分析

許晶晶,杜鵬飛,屈麗琴,趙 瑩,黃東浩,周瑞鵬

(1.中國水利水電科學研究院 國際泥沙研究培訓中心,北京 100048;2.沈陽農(nóng)業(yè)大學 水利學院,沈陽 110866)

土壤侵蝕會影響土壤健康和生產(chǎn)力[1]，引發(fā)泥石流和山體滑坡，導致水庫河道淤積、水體面源污染，造成土壤生物多樣性喪失[2]，從而減少緩解和適應氣候變化的能力[3]。由土壤侵蝕形成的泥沙因攜帶了包括土壤養(yǎng)分、重金屬、痕量金屬和氧化物等在內(nèi)的源地屬性[4]，因此，不同類型物源區(qū)的各屬性含量之間可能會存在差異。對沉積泥沙的溯源不僅能夠為明晰土壤侵蝕、搬運、沉積過程提供依據(jù)，而且可以助力水土保持措施的合理布設[5-6]。

接壩區(qū)作為內(nèi)蒙古高原向河北北部山地的緩沖過渡帶，是京津上游地區(qū)的重要水源地和環(huán)境脆弱區(qū)，區(qū)內(nèi)地貌溝深、坡陡，與緊鄰其北部的壩上高原形成了鮮明對比[7]。由于過度開墾和放牧，接壩區(qū)的生態(tài)環(huán)境受到了嚴重破壞，成為21世紀初京津風沙的主要源區(qū)之一，曾一度威脅著首都地區(qū)的生態(tài)安全。位于接壩區(qū)中部的小壩子鄉(xiāng)，地處半干旱農(nóng)牧交錯帶，生態(tài)環(huán)境十分脆弱，曾是京津浮塵和揚沙天氣的重要沙源地[8]。近年來，在大規(guī)模退耕還林還草治理下，區(qū)內(nèi)風沙危害顯著降低，生產(chǎn)生活環(huán)境明顯改善，反而由水力侵蝕引發(fā)的耕層土壤變薄、侵蝕溝擴張、地表破碎化等有加劇的趨勢，河道出現(xiàn)了不同程度的泥沙淤積。小壩子鄉(xiāng)小壩子村修建于2013年的四級谷坊，歷經(jīng)短短7 a，其中兩級已基本淤滿，其嚴重程度可見一斑。

已有研究顯示，接壩區(qū)不同面積小流域的土壤侵蝕存在很大差異[9]。在小流域尺度上，丘陵和低山頂部以風蝕為主，侵蝕溝溝谷以水蝕為主，位于中下坡位的耕地則受風水兩相侵蝕的影響；在微小流域尺度上，坡頂是以風蝕和水力面蝕為主的復合侵蝕，坡面中下部則以風蝕和水力溝蝕較為常見。為了對比這兩個不同面積小流域的泥沙來源，本研究以整個小壩子鄉(xiāng)小流域和小壩子村的一處微小流域為例，對沉積泥沙的主要源地及其貢獻比例進行分析。研究有助于加深對接壩區(qū)土壤侵蝕和泥沙輸移過程的理解與把握，對基于侵蝕預防的水土保持措施布設和土地利用類型調(diào)整具有非常重要的參考意義。

1 研究區(qū)概況

小壩子鄉(xiāng)位于內(nèi)蒙古高原和燕山山地連接的接壩地帶，土地總面積為325.40 km2，流域土壤以棕壤為主[10]。作為一個自然地理單元，小壩子鄉(xiāng)流域(116°12′49″—116°29′30″E，41°22′8″—41°34′6″N)(簡稱鄉(xiāng)流域，圖1A)內(nèi)最高點海拔2 019 m，最低點海拔874 m，地形起伏較大。行政區(qū)域上，南與窟窿山鄉(xiāng)接壤，東與土城鎮(zhèn)相連，北面和西北面分別為壩上的四岔口鄉(xiāng)和大灘鎮(zhèn)。該區(qū)屬于典型的大陸性季風型半干旱氣候，分屬接壩冷涼氣候類型區(qū)。年均氣溫4℃，年均無霜期105 d，多年年均降水量409.30 mm。區(qū)內(nèi)年均降雨量時空分布差異較大，6—9月份的降雨占全年總量的80%，且多以暴雨形式出現(xiàn)。

作為小壩子鄉(xiāng)流域內(nèi)的一個微小流域，小壩子村流域(簡稱村流域，圖1B)位于小壩子村西口，面積為0.164 km2。該流域坡地曾分布有農(nóng)田，因沙化嚴重于20世紀90年代末退耕。目前水力侵蝕嚴重，主溝道總長為242 m，最寬處達29 m，最深處32 m，侵蝕溝擴張迅速，溯源侵蝕顯著。

圖1 研究區(qū)概況和樣品采樣點分布

2 材料與方法

2.1 樣品采集

在夏秋降雨，尤其是個別場次暴雨的作用下，面上表土流失和溝道擴張發(fā)育同時出現(xiàn)，是本區(qū)土壤水力侵蝕的一個重要特征。對鄉(xiāng)流域而言，由于京津風沙源治理工程的持續(xù)推進，當?shù)卣e極進行退耕還林還草，使得研究區(qū)林地、草地分布最廣，面積占比高達九成。流域上游坡頂為草地，坡中為成片的樺木林，總體覆蓋度較高。中下游林地和草地呈塊狀分布，但總體覆蓋度較低，中坡位和下坡位有農(nóng)田分布，表土流失嚴重。加之溝道內(nèi)坍塌的堆積物在夏秋雨季被沖刷，使得底土流失也不容小覷。為全面反映泥沙來源情況，將農(nóng)田表土、林草地表土及侵蝕溝底層物質(zhì)確定為該流域的潛在泥沙源區(qū)。對村流域而言，坡頂為沙化較為嚴重的草地，坡中及坡腳為不同沙化程度的退耕農(nóng)田。由于沙化嚴重，粗骨化較為明顯，在暴雨沖刷下，村流域內(nèi)溝道縱橫，主溝道垮塌現(xiàn)象顯著。在這種情況下，退耕農(nóng)田表土、草地表土和侵蝕溝底層物質(zhì)就成為村流域潛在的泥沙源區(qū)。根據(jù)上述分類，于2018年4月進行了樣品采集。對鄉(xiāng)流域，通過遙感影像研判解譯和實地勘查，選擇了有代表性的源區(qū)17處，其中包括農(nóng)田5處、林草地9處、侵蝕溝3條。為降低各源區(qū)指紋因子的空間變異，增加樣本的代表性，采樣使用多點混合法，以20 m為半徑，重復采集樣品3個進行充分混合。對村流域，使用網(wǎng)格法采集潛在泥沙源區(qū)樣品，共采集草地表層樣11個，沙地表層樣10個，侵蝕溝底層樣13個。上述所有表層樣品的采集深度均為5 cm，侵蝕溝底層樣品的采集部位均在溝壁距表層高度超過30 cm的亞表層和母質(zhì)層。沉積樣品均在流域出口處采集，沿主河道每隔50 m(鄉(xiāng)流域)和5 m(村流域)采集沉積表層樣，采集個數(shù)均為4個，采集深度約為5 cm?？紤]到泥沙輸移的連通性，樣品的采集區(qū)相對靠近溝道，這些位置發(fā)生侵蝕后，泥沙更易搬運至河道直至沉積區(qū)。具體樣品采集信息見表1。

表1 兩個流域樣品采集信息

2.2 樣品分析

樣品采集后，經(jīng)自然風干、手工剔除雜質(zhì)、研磨后，分別過2 mm和0.063 mm的篩子用于樣品測試。測試工作在中科院寒區(qū)旱區(qū)環(huán)境與工程研究所完成。其中，小于2 mm的部分用于測定土壤粒徑分布，測試儀器為激光粒度儀Mastersizer 3000；小于0.063 mm的部分用于測試樣品中的指紋因子，測試儀器為X射線熒光光譜分析儀，測試的指紋因子包括30個，分別是：P，Mn，Cr，Co，Cu，Pb，Zn，Ti，Ni，V，Sr，Ga，Rb，Nb，Ba，Zr，As，Cl，Br，La，Y，Nd，Ce，Al2O3，SiO2，F(xiàn)e2O3，K2O，Na2O，MgO，CaO。

2.3 研究方法

本研究采用基于指紋因子測定的指紋識別技術來計算不同源區(qū)的泥沙來源及其貢獻比。指紋識別技術通過對比物源和泥沙的指紋因子關系來實現(xiàn)對不同源地泥沙貢獻的定量估算，已在我國黃土高原[11-15]、長江上游[16-19]、黔中喀斯特高原[20]、東北黑土區(qū)[21]、南方紅壤區(qū)[22-23]等多地得到了成功應用。這一方法篩選指紋因子時主要包括以下3個步驟：

(1)雙邊范圍檢驗。通過比較源區(qū)的樣品和沉積區(qū)樣品的元素濃度進行雙邊范圍檢驗，當沉積區(qū)樣品元素濃度大于源區(qū)樣品元素濃度的最大值或者小于最小值時，將其剔除，剩下的即為保存性指紋因子。

(2)Kruskal-Wallis H檢驗。將通過雙邊范圍檢驗的指紋因子進行Kruskal-Wallis H檢驗，評價單個指紋因子在不同泥沙源區(qū)中是否存在顯著差異。p<0.05表示單個指紋因子在不同物源區(qū)間差異顯著，該指紋因子通過Kruskal-Wallis H檢驗，可用于判定泥沙來源。

(3)多元逐步判別分析。將上一步驟中篩選出的指紋因子進行逐步判別分析，每一步進入判別函數(shù)的都是Wilks′ λ的計算值最小的指紋因子。將所有進入判別函數(shù)的指紋因子進行組合，即為最佳指紋因子組合。

基于指紋識別計算泥沙來源的模型主要包括多元線性混合模型和貝葉斯模型。本研究采用的是由Collins等[24]提出的多元線性混合模型，其計算可靠性和準確性已得到了有效驗證[25-26]。

Walling-Collins多元線性混合模型如下：

(1)

式中:n是泥沙源區(qū)個數(shù);Ci是沉積泥沙中指紋因子i的濃度;m是指紋因子個數(shù);Ps是泥沙源區(qū)s的泥沙貢獻率;Ssi是泥沙源區(qū)s中指紋因子i的平均濃度。

利用擬合優(yōu)度算法檢驗模型計算結果的不確定性，見公式(2)：

(2)

式中:MAF是平均絕對擬合度值,一般認為MAF>0.8時,可以接受多元混合線性模型的結果。

3 結果與分析

3.1 不同面積小流域泥沙粒徑特征

在鄉(xiāng)、村兩個典型流域中，不同土地利用類型的泥沙顆粒組成差異明顯。由圖2A可知，在鄉(xiāng)流域，與農(nóng)田樣品和侵蝕溝樣品相比，不同林草地樣品的粒徑變化較大，采集于地表覆蓋厚達15 cm枯枝落葉層的樺木林表土中的砂粒含量最低，僅為20.27%，采集于河灘、沙化嚴重、基本沒有植被覆蓋地表的次生楊樹林表土樣品砂粒含量最高，達93.14%。農(nóng)田樣品粉粒和黏粒含量是所有源區(qū)樣品中最高的，為40.93%，8.63%，較侵蝕溝和林草地分別高了16.14%，27.63%，2.25%，5.90%，這可能反映了農(nóng)田在長期耕作過程中下翻細粒土壤融入耕層的顯著影響。

由圖2B可知，村流域物源區(qū)樣品砂粒含量均超過70%，退耕農(nóng)田砂粒含量最高，達96.81%，粗?；蠲黠@。對比圖2A和2B，村流域退耕農(nóng)田與侵蝕溝砂粒含量分別為96.81%，91.10%，比鄉(xiāng)流域農(nóng)田和侵蝕溝分別高出46.37%，7.14%。而沿主溝道向流域外延的沉積樣品砂粒含量逐漸降低，黏粒和粉粒含量逐漸增加，說明在搬運過程中，粗顆粒泥沙先沉積，而細顆粒泥沙則能夠被搬運得更遠。

圖2 兩個流域泥沙粒徑百分比

由圖3可見，兩個流域物源區(qū)的中值粒徑中，林草地基本一致，侵蝕溝略有差異，農(nóng)田和退耕農(nóng)田差異較大。村流域退耕農(nóng)田的中值粒徑為190.40 μm，是鄉(xiāng)流域農(nóng)田中值粒徑的7.40倍。沉積區(qū)泥沙的中值粒徑則基本相等。在鄉(xiāng)流域內(nèi)，中值粒徑由沉積區(qū)、侵蝕溝、林草地和農(nóng)田的順序依次減小，在村流域內(nèi)，則略有不同，按沉積區(qū)、退耕農(nóng)田、侵蝕溝和草地的順序依次減小。不論是鄉(xiāng)流域，還是村流域的中值粒徑差異均較為明顯地反映出不同土地利用對土壤或泥沙顆粒組成的影響。沉積區(qū)泥沙中值粒徑最大，這可能與泥沙輸移的分選性有關，即被搬運至流域外的細粒泥沙比例較高而在流域內(nèi)沉積的粗粒泥沙較多。也可能與泥沙連通性有關，本研究結果則直觀反映了沉積區(qū)與較大中值粒徑源區(qū)距離對沉積的影響，即在鄉(xiāng)流域，沉積區(qū)與侵蝕溝的距離最近，而在村流域，沉積區(qū)與退耕農(nóng)田的距離最近。但同時也應該認識到，由于土壤侵蝕的發(fā)生具有顯著的時空變異特征，泥沙連通性理論上也會出現(xiàn)類似的變化，加之其在反映土壤侵蝕嚴重程度上存在滯后性，所以不論對于單次侵蝕過程，還是一定時段的累積過程而言，準確刻畫泥沙連通性的變化都較為困難。

圖3 兩個流域中值粒徑對比

3.2 最佳指紋因子組合確定

對測定的30個指紋識別因子進行分析。第一步，利用雙邊范圍檢驗篩選出兩個流域中分別有25個和29個保存性指紋因子；第二步，通過非參數(shù)Kruskal-Wallis H檢驗篩選在不同物源區(qū)間有顯著差異的指紋因子(p<0.05)，兩個研究區(qū)分別有7個和25個潛在指紋因子通過檢驗；第三步，通過逐步判別分析確定最佳指紋因子組合。各級指紋因子篩選結果見表2。由于村流域的地形更簡單，泥沙在搬運、沉積過程中更為穩(wěn)定，所以三級篩選中，村流域通過檢驗的指紋因子個數(shù)均比鄉(xiāng)流域多。

表2 兩個流域各級篩選后指紋因子對比

兩個不同面積流域的逐步判別分析結果見表3。鄉(xiāng)流域以K2O，Cl，Ce為最佳指紋因子組合，能以85.7%準確度判別該研究區(qū)的農(nóng)田表土、林草地表土和侵蝕溝3個泥沙源區(qū)。每步累積判別率為61.9%，66.7%，85.7%，單步驟判別率分別為61.9%，33.3%，57.1%。村流域以Ba，P，CaO，SiO2這4個因子作為最佳指紋因子組合，能以89.5%準確度判別該研究區(qū)的退耕農(nóng)田表土、草地表土和侵蝕溝3個泥沙源區(qū)。每步累積判別率為60.5%，63.2%，78.9%，89.5%，單步驟判別率分別為60.5%，65.8%，57.9%，63.2%。

表3 兩個流域的逐步判別分析結果

3.3 不同面積小流域泥沙來源分析

3.3.1 不同源區(qū)的泥沙相對貢獻率 根據(jù)最佳指紋因子組合中各因子的物源區(qū)和沉積區(qū)泥沙樣品中的濃度，利用Walling-Collins模型計算不同物源區(qū)的泥沙相對貢獻率。農(nóng)田表土和退耕農(nóng)田表土分別是鄉(xiāng)流域和村流域主要的泥沙源區(qū)。在鄉(xiāng)流域，農(nóng)田表土對沉積泥沙的貢獻比高達71.6%，是第二源區(qū)林草地表土(15.9%)的4.5倍，第三源區(qū)侵蝕溝僅貢獻了12.5%的泥沙。在村流域，退耕農(nóng)田表土、侵蝕溝和草地表土的泥沙貢獻比差距相對小一些，分別占比44.4%，34.2%，21.4%。

為了檢驗計算結果的穩(wěn)定性，使用平均絕對擬合度值(mean absolute fit，MAF)評估模型(公式4)對計算結果進行了評估，鄉(xiāng)流域和村流域兩個研究區(qū)的泥沙來源的MAF分別為0.97，0.98，遠高于0.8的臨界可接受值，表明Walling-Collins模型對兩個流域泥沙貢獻率的計算結果是可靠的。

3.3.2 不同源區(qū)的泥沙面積比 為明確不同土地利用類型在流域內(nèi)的面積占比與泥沙貢獻比的關系，利用樣品采集當年(2018年)的遙感影像對流域土地利用類型進行了解譯。由于小壩子鄉(xiāng)土地利用類型比較復雜，解譯過程中參照《區(qū)域水土流失動態(tài)監(jiān)測技術規(guī)定》，依據(jù)土地利用二級分類標準進行，解譯結果見圖1和表4。

由圖1和表4可知，鄉(xiāng)流域土地利用類型以林地和草地為主，占總面積的90.21%，其次為農(nóng)田，占比6.47%，侵蝕溝占比2.13%，剩余為建設用地等其他用地。村流域的退耕農(nóng)田和草地分別占比41.38%，57.26%。根據(jù)現(xiàn)場調(diào)查，結合無人機影像，確定的侵蝕溝面積為0.014 2 km2，面積占比8.64%。

表4 兩個流域的土地利用類型面積及占比

為反映各源區(qū)泥沙輸移沉積的嚴重程度，引入泥沙面積比來定量表達各侵蝕源區(qū)泥沙搬運至流域出口沉積區(qū)的比例和其所占流域面積的比例，定義為這兩個比例的比值。當比值結果大于1時，表明流域內(nèi)該類型物源區(qū)流失土壤的比例偏高，相應地，土壤侵蝕強度也越高；反之，表明流域內(nèi)該類型物源區(qū)流失土壤的比例不是很嚴重，值越接近于0，土壤流失后搬運至沉積區(qū)的比例越少。

由圖4可見，除鄉(xiāng)流域的林草地和村流域的草地表土泥沙面積比小于1外，其他4個侵蝕源區(qū)的泥沙面積比均高于1。尤其是鄉(xiāng)流域的農(nóng)田表土泥沙面積比高達11.07，其土壤侵蝕、搬運的比例遠高于村流域退耕農(nóng)田表土泥沙面積比的1.07。一方面，這可能與農(nóng)田在人為擾動后，表土變得疏松，從而更易于流失有關。每年的雨季也恰是農(nóng)田的除草季，相比于退耕農(nóng)田，農(nóng)田的人類活動干擾更加強烈，由此導致表土更容易被侵蝕掉。另一方面，從農(nóng)田在鄉(xiāng)流域的分布來看，大多集中在下坡位，距離溝道較近，受水系連通性強的影響，表土侵蝕后更容易被水流搬運至河道中，從而其在流域出口沉積的概率也更高。

圖4 兩個流域的源地泥沙面積比

鄉(xiāng)流域和村流域的侵蝕溝泥沙面積比分別為5.87，3.96，表明該研究區(qū)域的侵蝕溝不斷擴張發(fā)育導致的底層土壤流失是不容小覷的。侵蝕溝的快速擴張使坡面徑流和泥沙更加容易在溝道內(nèi)匯聚，從而增大了流域內(nèi)泥沙輸移的連通性，間接加劇了土壤侵蝕后的輸移比例。鄉(xiāng)流域的林草地和村流域的草地泥沙面積比分別只有0.18，0.43，與林草地植被覆蓋度總體較高，能夠增加泥沙攔截，從而降低侵蝕搬運量有關。

從上述結果來看，村流域后續(xù)的水土流失治理重點應著眼于侵蝕溝，鄉(xiāng)流域除了侵蝕溝外，還要重視農(nóng)田的水土流失防治，實施以工程措施、生物措施和耕作措施相結合的綜合治理模式。具體來說，在主溝道中修建淤地壩，支溝中修建梯級谷坊，進行分段控制，并在溝頭植樹種草減緩溯源侵蝕。坡耕地應減少雨季擾動的頻次和強度，通過免耕、等高耕作等措施降低水力侵蝕的危害。

4 結論與討論

本研究在接壩區(qū)選取了小壩子鄉(xiāng)和小壩子村兩個不同面積的小流域定量計算了不同源區(qū)的泥沙來源比例。對研究區(qū)樣品的粒徑分析發(fā)現(xiàn)，同流域內(nèi)不同源區(qū)間泥沙粒徑差異明顯，不同流域同一源區(qū)的泥沙粒徑差異也很大，但兩個流域沉積泥沙粒徑較為接近。鄉(xiāng)流域農(nóng)田表土、林草地表土和侵蝕溝3種源區(qū)泥沙來源的相對貢獻率分別為71.6%，15.9%，12.5%；村流域農(nóng)田表土、草地表土和侵蝕溝3種源區(qū)泥沙來源的相對貢獻率分別為44.4%，21.4%，34.2%。鄉(xiāng)流域泥沙面積比最大的是農(nóng)田，比值達11.07，其次為侵蝕溝5.87，最小的林草地僅為0.18；村流域中泥沙面積比由大到小分別為侵蝕溝3.96，退耕農(nóng)田1.07，草地0.43。研究結果表明，侵蝕溝和農(nóng)田是未來水土流失防治的重點區(qū)域。

本研究明晰了接壩區(qū)不同面積小流域泥沙來源的基本特點，相關結果可為該區(qū)水土流失治理和土地利用規(guī)劃編制提供參考。但同時應該注意到，由于接壩區(qū)內(nèi)海拔高差大，受降雨、植被、泥沙連通性及人為擾動強度的影響，不同流域的侵蝕產(chǎn)沙和輸移沉積也會存在一些差異。在今后的研究中應綜合考慮這些因素的綜合影響，以便深入比較不同條件下各尺度流域的泥沙來源情況，分析不同空間尺度逐級遞進對計算結果的影響，從而為泥沙來源跨尺度銜接計算提供參考。