李國強，方海燕

（1杭州市水文水資源監(jiān)測總站，浙江 杭州310016；2中國科學院地理科學與資源研究所 陸地水循環(huán)及地表過程重點實驗室，北京100101）

土壤侵蝕是世界性的環(huán)境問題，它不但引起土壤及其土壤養(yǎng)分的流失，侵蝕掉的泥沙進入河道后還會引起湖泊、水庫的淤積以及環(huán)境污染等問題［1］。東北黑土區(qū)是我國重要的商品糧基地，黑土有機質含量豐富、肥力高、性狀好，土地農(nóng)作物產(chǎn)量高。自20世紀50年代以來，隨著人口的增加和開墾強度的增大，黑土層厚度從開墾初期的60～70cm下降到目前的20～30cm，有些地方已出現(xiàn)破皮黃［2］，嚴重威脅到國家的糧食安全。

土壤侵蝕模型是土壤侵蝕和水土保持研究的重要手段。然而，由于東北黑土區(qū)土壤侵蝕研究相對較晚，小流域徑流、產(chǎn)沙觀測相對較少，該區(qū)土壤侵蝕模擬研究多采用 （R）USLE（（Revised）Universal Soil Loss Equation）模型［3－5］，但（R）USLE 模型計算侵蝕量時忽略了泥沙的沉積和空間的變異性［6］。其他的一些物理機制模型如 WEPP和LISEM等對數(shù)據(jù)要求十分嚴格，在黑土區(qū)流域尺度上尚無應用。比利時魯汶大學研發(fā)的 WaTEM／SEDEM（Water and Tillage Erosion Model and Sediment Delivery Model）模型所需數(shù)據(jù)少，在土壤、地形、土地利用和降雨等基礎數(shù)據(jù)具備的情況下，即可運行［7］。此外，該模型結構簡單，所需校正參數(shù)少，同時又能模擬土壤侵蝕產(chǎn)沙及分布特征；除了產(chǎn)沙觀測數(shù)據(jù)外，湖泊、水庫等的泥沙淤積信息也可對該模型進行校正和驗證。因此，自 WaTEM／SEDEM模型提出以來，在國外許多國家［8－10］得到了很好的應用，在我國的紅壤區(qū)小流域［11－12］也利用該模型開展了工作。

截至目前，國內(nèi)外多應用出口產(chǎn)沙校正和驗證WaTEM／SEDEM模型，使得不同的模型參數(shù)組合都能很好地模擬流域出口產(chǎn)沙量，表明流域內(nèi)侵蝕沉積的多種組合均會有相同的流域出口產(chǎn)沙［13］。然而，利用流域內(nèi)部侵蝕沉積強度數(shù)據(jù)校正模型，將可能真實地再現(xiàn)流域侵蝕沉積特征，同時也能夠模擬流域出口產(chǎn)沙量。

東北黑土區(qū)土壤侵蝕研究開展較晚，一些地區(qū)缺乏水沙觀測數(shù)據(jù)，但該區(qū)有許多的湖泊和水庫泥沙淤積信息數(shù)據(jù)，從而為該模型的應用提供了很好數(shù)據(jù)支持。此外，核示蹤流域內(nèi)侵蝕強度在黑土區(qū)已成功開展，利用核示蹤得到的土壤侵蝕強度可校正模型［13］。因而，在黑土區(qū)聯(lián)合應用 WaTEM／SEDEM模型及核示蹤技術開展土壤侵蝕研究，將是彌補該區(qū)侵蝕產(chǎn)沙數(shù)據(jù)較少、模型模擬研究不足的重要舉措。

1 材料與方法

1.1 研究區(qū)概況

拜泉縣（E125°30′～126°31′，N47°20′～47°55′）位于黑龍江省齊齊哈爾市東北部，地處小興安嶺余脈與松嫩平原的過渡地帶，屬于東北典型黑土區(qū)，地形特征為漫川漫崗。拜泉縣為中溫帶大陸性氣候，多年平均氣溫為1.28℃，全年平均降水量500mm左右，多集中在7、8、9三個月，降雨歷時短，強度大，易形成徑流侵蝕。土壤以典型黑土和草甸黑土為主。流經(jīng)縣境內(nèi)的河流有通肯河、雙陽河、潤津河。為了發(fā)展農(nóng)業(yè)，拜泉縣從20世紀60年代開始修建水庫塘壩等蓄水工程，沉積了大量泥沙。本文選取該縣境內(nèi)的24個水庫，水庫控制流域土地利用較為單一，以耕地為主（表1）。

1.2 數(shù)據(jù)來源

1.2.1137Cs數(shù)據(jù)采集與處理 開展研究的24個水庫中，齊心水庫控制流域地形適中，土地利用類型較為全面，具有一定的代表性。因此，2013年10月，在齊心水庫控制流域，根據(jù)土地利用類型、地形特點，利用5cm直徑的土鉆采集土壤樣品，采樣深度視黑土層厚度而定，坡面侵蝕區(qū)采集深度40cm左右，沉積區(qū)100cm左右，以保證采樣點所有的137Cs含量都能夠采集到。采集的25個樣點土壤經(jīng)晾干、研磨后過2mm篩，剔除草根和石塊等雜物，稱300g左右供核元素測試。137Cs放射性活度用美國CANBERRA公司配備的高純鍺探頭λ能譜儀（測定范圍為3～3×103keV）測定。在662keV處求算137Cs放射性活度（Bq／kg），核素測定時間為 80 000s左右，測定誤差 ＜8%。

表1 研究區(qū)24座水庫的控制流域面積、耕地比例和年均產(chǎn)沙量Tab.1 The areas，percentages of cultivated lands and sediment yields of the 24reservoir catchments

其中：A（t）為單位面積累積137Cs活度（Bq／m2，由采樣點核素放射性活度與土鉆橫截面積得到）；R為侵蝕強度（kg／（m2·a））；t為137Cs釋放以來經(jīng)歷的時間；d為累積質量深度（kg／m2），λ為137Cs衰變常數(shù)；I（t）為137Cs年沉降通量（Bq／（m2·a））；Γ為土壤侵蝕新沉降的137Cs在耕作之前被侵蝕移動的比例，P為顆粒校正系數(shù)。

在沉積樣點（土壤樣品核素活度大于背景值區(qū)核素活度）：

1.2.2 土壤侵蝕強度計算 MBM2模型在黑土區(qū)核素含量和土壤侵蝕轉換中已得到成功運用［14］。在侵蝕樣點（土壤樣品核素活度小于背景值區(qū)核素活度）：

其中：Aex為測量得到的137Cs活度與其背景值Aref之差；Cd（t′）為沉積采樣點土壤的137Cs活度（Bq／kg）。

1.3 WaTEM／SEDEM 模型的運行

1.3.1 模型簡介 WaTEM／SEDEM是在RUSLE基礎上構建的分布式模型，它分為侵蝕模塊、泥沙輸移能力模塊和泥沙運移模塊［15］。模型年均輸沙能力Ct（kg／（m·a））用來反映地表植被覆蓋對輸沙能力的影響，模型僅需要校正泥沙輸移能力系數(shù)Ktc。地表植被覆蓋較好，Ktc值較低；地表植被覆蓋較差，則Ktc較高。

1.3.2 模型數(shù)據(jù)來源及處理

（1）數(shù)字高程模型（DEM，分辨率為20m）由國家基礎地理信息中心1991年航攝的1∶5萬地形圖數(shù)字化得到；（2）土地利用數(shù)據(jù)利用2000年TM影像（30m分辨率）采用人工目視解譯的方式獲取，并進行了野外考察校正；（3）降水數(shù)據(jù)由拜泉縣氣象局及中國氣象科學數(shù)據(jù)共享網(wǎng)獲得；（4）拜泉縣土壤類型空間分布圖（1∶10萬）由中國科學院東北地理與農(nóng)業(yè)生態(tài)研究所獲得；（5）水庫控制流域年均產(chǎn)沙量（SY）由黑龍江水文圖集中44個觀測站36年實測產(chǎn)沙利用水文比擬法［16］獲取，并經(jīng)三道鎮(zhèn)水庫和齊心水庫高精度GPS差分法和實地采樣得到了驗證。

WaTEM／SEDEM模型輸入數(shù)據(jù)有DEM、降雨侵蝕力因子R、土壤可蝕性因子K、土地管理因子C和水土保持措施因子P，河流水系圖等。河流水系由DEM通過ArcGIS水文模塊提取，匯水面積大于1 km2即可視為河道［9］；研究區(qū)C和P值參考黑土區(qū)已有研究［5，17］及實地調(diào)查得到。土壤可蝕性因子K根據(jù)拜泉縣1984年土壤普查資料及黑土區(qū)有關研究成果［17－18］。研究區(qū)及周邊地區(qū)有些站點缺少逐日降水數(shù)據(jù)，利用年均降雨量方法獲取降雨侵蝕力R 值［19］：

其中：R為多年平均降雨侵蝕力（MJ·mm／（h·hm2·a）），P為多年平均降雨量（mm）。在ArcGIS9.3中通過克里金插值方法得到各水庫控制流域降雨侵蝕力值。

1.3.3 模型校正 利用137Cs示蹤法，根據(jù)25個樣點的土壤侵蝕強度對模型參數(shù)進行校正。先初步根據(jù)已發(fā)表的論文確定Ktc的范圍，以減少模型運行次數(shù)，之后根據(jù)Ktc組合，得到流域侵蝕沉積特征，提取各采樣點的侵蝕沉積強度，并與137Cs示蹤得到的侵蝕沉積速率對比；根據(jù)Ktc高值和低值的多次組合運行模型。采用模型有效系數(shù)（cNS）確定Ktc最佳參數(shù)組合。cNS計算方法如下：

其中n是指觀測的數(shù)目；Oi是指觀測值；Omean是指觀測平均值；Pi是指預測值。cNS的取值范圍是（－∞，1），cNS值越接近于1，模擬效果越好。

2 結果分析

2.1 采樣點侵蝕沉積強度

有研究［20］已經(jīng)測量了研究區(qū)的137Cs背景值，因而，根據(jù)背景值和25個采樣點核素面積活度及MBM2（公式1和2），得到了采樣點侵蝕強度。侵蝕強度受土地利用類型和地形的影響，坡耕地侵蝕強度平均值為－24.5t／（hm2·a），最大侵蝕強度達到－45.7 t／（hm2·a），最大沉積速率為15.4t／（hm2·a）。林地平均侵蝕強度為－1.1t／（hm2·a），最大侵蝕強度為－4.5t／（hm2·a），最小值（沉積速率）為0.6t／（hm2·a）。灌木林地侵蝕強度介于二者之間，平均侵蝕強度為－5.6t／（hm2·a）。研究區(qū)草地平均侵蝕強度為－9.5t／（hm2·a）。整個齊心流域采樣點侵蝕速率變異很大，平均侵蝕強度為－19.8t／（hm2·a）。

2.2 模型校正和驗證

利用以上25個采樣點的侵蝕沉積強度數(shù)據(jù)對模型進行校正。根據(jù)地表徑流泥沙輸移敏感性的不同，對兩種截然不同的土地利用類型賦予不同的Ktc值。林草地泥沙輸移能力系數(shù)賦予低值Ktc－low，耕地則賦予高值 Ktc－h(huán)igh，Ktc變化范圍在0～12之間。Ktc－low和Ktc－h(huán)igh不同組合下得到采樣點模擬侵蝕強度，通過137Cs示蹤和模擬侵蝕強度公式（4）計算 NS變化值。圖1顯示最佳 Ktc－low和 Ktc－h(huán)igh值組合為4m和8m，此時cNS為0.62。

根據(jù)137Cs示蹤法校正后的模型，利用表1中24個水庫（不包括齊心水庫）控制流域產(chǎn)沙量（YS）與產(chǎn)沙模數(shù)YSS（YS與流域面積的比值）進行模型驗證。圖2表明，模型模擬的YSS和YS值和觀測值相近，其相關系數(shù)分別達到0.86和0.7，模型對于年產(chǎn)沙的模擬效果優(yōu)于產(chǎn)沙模數(shù)。

表2 核素137Cs示蹤得到的采樣點土壤侵蝕（負數(shù)）沉積（正數(shù)）強度Tab.2 Activities of 137Cs and its derived soil erosion rates for the sampling points

3 討論

在東北黑土區(qū)，已成功利用核示蹤技術開展了土壤侵蝕研究。方海燕等［14］研究發(fā)現(xiàn)，在東北鶴山農(nóng)場面積為28.5hm2的農(nóng)耕地小流域上，土壤凈侵蝕強度為2.0t／（hm2·a），遠小于本研究中耕地平均土壤侵蝕強度；而流域內(nèi)侵蝕區(qū)的侵蝕強度均值為－15.0t／（hm2·a），略小于本研究得到的土壤侵蝕強度。在拜泉縣小流域坡面上，于寒青等［20］發(fā)現(xiàn)等高種植、梯田和順坡坡面上土壤平均侵蝕強度為分別為－23、－25和－35t／（hm2·a），大于本研究平均值，這是因為他們的采樣點位于坡度較大的坡面侵蝕區(qū)。從以上比較來看，研究中得到的侵蝕速 率數(shù)據(jù)可信。

圖1 利用采樣點的侵蝕速率對泥沙輸移能力系數(shù)Ktc－low（a）和Ktc－h(huán)igh（b）進行校正Fig.1 Validations of the Ktc－low（a）and Ktc－h(huán)igh（b）values using 137Cs－derived erosion rates of the sampling points

圖2 水庫控制流域產(chǎn)沙觀測值與WaTEM／SEDEM模擬值：（a）年均產(chǎn)沙量，（b）年均產(chǎn)沙模數(shù)Fig.2 Modeled sediment yield per year（a）and specific sediment yield（b）versus the counterparts from observed values for the 24reservoir catchments

研究發(fā)現(xiàn)，利用流域出口產(chǎn)沙數(shù)據(jù)校正模型時，雖然有時流域出口產(chǎn)沙的模擬效果很好，但流域內(nèi)土壤侵蝕沉積模式可能有多個。如在西班牙Barasona水庫流域，Alatorre等［13］發(fā)現(xiàn)有一系列的Ktc高值和低值組合，均能得到相同的cNS值。137Cs之所以能夠校正WaTEM／SEDEM模型，除了能夠反映空間上土壤侵蝕沉積模式外，另外一個重要特點即是二者的時間尺度相同或相似。137Cs示蹤得到的是50年以來土壤的平均侵蝕強度；WaTEM／SEDEM也是一個模擬多年時間尺度上土壤平均侵蝕強度的模型，只要模型模擬時將二者的時間跨度調(diào)為一致，即可消除時間尺度的影響。目前，利用137Cs示蹤得到的土壤侵蝕強度校正 WaTEM／SEDEM模型已成功應用［13］。

本研究中，利用核示蹤法得到的最佳Ktc組合為4m和8m，該值小于Shi等人［11］在王家橋流域Ktc組合值（15m 和55m），也小于 Van等人［8］利用西班牙21個流域產(chǎn)沙值得到的Ktc值組合（75m和250m）。在西班牙小流域，Alatorre等人［13］利用137Cs方法得到的Ktc最佳組合值為2.84m和9.05 m，與本研究較為接近。在數(shù)據(jù)不充分時，有些研究甚至將 Ktc－low和 Ktc－h(huán)igh固定為一個比值［11］。有研究發(fā)現(xiàn)［21］，不同侵蝕類型區(qū)和不同數(shù)字高程模型分辨率情況下，校正得到的參數(shù) Ktc－low和 Ktc－h(huán)igh值差別會很大。齊心水庫控制面積為9.9km2，研究中25個土壤采樣點偏少，但利用采樣點獲得的侵蝕沉積速率對模型校正后，流域出口產(chǎn)沙的模擬效果很好（圖2）。值得注意的是，YS較YSS的模擬精度高，這是因為YS值變化幅度大，因而對于YS的模擬效果優(yōu)于YSS模擬能力。

為進一步探討模型模擬的精確性，利用校正好的模型得到了各個小流域總侵蝕量，進而根據(jù)小流域總產(chǎn)沙量與總侵蝕量的比值，得到了泥沙輸移比RSD（圖3），結果表明研究區(qū)小流域泥沙輸移比變化于0.16和0.40之間，且隨著流域面積的增加，泥沙輸移比呈現(xiàn)減小的趨勢。這符合傳統(tǒng)的產(chǎn)沙模數(shù)－流域面積關系［22］。此外，在鶴山農(nóng)場，Dong等［23］研究發(fā)現(xiàn)面積為27.29km2的小流域泥沙輸移比為0.1；簡金世［24］在東北漫川漫崗區(qū)三個水文站得到平均泥沙輸移比為0.145；Fang等［14］發(fā)現(xiàn)黑土區(qū)28.5hm2的小流域RSD大于0.5，這可能與流域面積小有關。此外，野外考察發(fā)現(xiàn)，河道上游的溝道表現(xiàn)為侵蝕，而河道中下游多有泥沙沉積。這一方面說明圖3的合理性，另一方面也表明與面積有關的YSS模擬效果要低于YS，這也從側面證實了圖2的準確性。

以上表明，利用137Cs示蹤得到的土壤侵蝕沉積速率能夠很好校正WaTEM／SEDEM模型。因此，在水沙監(jiān)測資料稀少地區(qū)，核素示蹤和模型 Wa－TEM／SEDEM的結合運用，對于該區(qū)及類似研究區(qū)的土壤侵蝕研究具有很好的應用前景。

圖3 水庫控制小流域泥沙輸移比隨流域面積的變化特征Fig.3 Changing pattern of sediment delivery ratio with catchment areas

4 結論

WaTEM／SEDEM是一個基于RUSLE的模型，模型參數(shù)少，數(shù)據(jù)需求量低，水文站實際觀測數(shù)據(jù)和淤積信息及137Cs法都可對模型進行校正和驗證。利用24個水庫流域產(chǎn)沙量驗證后，發(fā)現(xiàn)模型模擬的效果較好。取得的結論如下：

（1）采樣點坡耕地侵蝕強度為－45.7～－15.4 t／（hm2·a），平均為－24.5t／（hm2·a），林地侵蝕強度較小，僅為－1.1t／（hm2·a）。灌木林地侵蝕強度介于二者之間，平均侵蝕強度為－5.6t／（hm2·a）。整個齊心流域采樣點侵蝕速率變異很大，平均侵蝕強度為－19.8t／（hm2·a）。

（2）研究區(qū)模型泥沙輸移系數(shù)Ktc最佳組合為4 m和8m，此時模型有效系數(shù)cNS為0.62。通過24個水庫控制流域驗證得到產(chǎn)沙量YS和產(chǎn)沙強度YSS的模擬值與實測結果相近，模擬效果很好。

（3）核示蹤技術和 WaTEM／SEDEM模型聯(lián)合運用，能夠很好地開展黑土區(qū)土壤侵蝕和水土保持研究，也可為該區(qū)的水土保持工作提供科學依據(jù)和理論指導。

［1］Novotny V，Chesters G.Delivery of sediments and pollutants from nonpoint sources：A water quality perspective［J］.Journal of Soil Water Conservation，1989，44：568－576.

［2］范昊明，蔡強國，王紅閃.中國東北黑土區(qū)土壤侵蝕環(huán)境［J］.水土保持學報，2004，18（2）：66－70.

［3］張憲奎，許靖華，盧秀琴，等.黑龍江省土壤流失方程的研究［J］.水土保持通報，1992，12（4）：1－9.

［4］張樹文，王文娟，李穎，等.近50年來三江平原土壤侵蝕動態(tài)分析［J］.資源科學，2008，30（6）：843－849.

［5］王文娟，張樹文，方海燕.東北典型黑土區(qū)坡溝侵蝕耦合關系［J］.自然資源學報，2012，27（12）：2113－2122.

［6］Verstraeten G.Regional scale modeling of hillslope sediment delivery with SRTM elevation data［J］.Geomorphology，2006，81（1／2）：128－140.

［7］WaTEM／SEDEM manual ［EB／OL］.http：／／www.kuleuven.be／geography／frg／modelling／，2014－11－30.

［8］Van R A，Paolo B，Robert J，et al.Modeling sediment yields in Italian catchments［J］.Geomorphology，2005，65（1／2）：157－169.

［9］Verstraeten G，Prosser lan P，F(xiàn)ogarty P.Predicting the spatial patterns of hillslope sediment delivery to river channels in the Murrumbidgee catchment，Australia［J］.Journal of Hydrology，2007，334（3／4）：440－454.

［10］Boix－Fayos C，De V J，María M M，et al.The impact of land use change and check－dams on catchment sediment yield［J］.Hydrological Process，2008，22（25）：4922－4935.

［11］Shi Z H，Ai L，F(xiàn)ang N F，et al.Modeling the impacts of integrated small watershed management on soil erosion and sediment delivery：A case study in the Three Gorges Area，China［J］.Journal of Hydrology，2012，438／439：156－167.

［12］方海燕，孫莉英，聶彬彬，等.基于 WaTEM／SEDEM模型的雙楓潭流域侵蝕產(chǎn)沙模擬 ［J］.陜西師范大學學報：自然科學版，2014，42（1）：92－97.

［13］Alatorre L C，Begueria S，Lana－Renault N，et al.Soil erosion and sediment delivery in a mountain catchment under scenarios of land use change using a spatially distributed numerical model［J］.Hydrology and Earth System Sciences，2012，16（5）：1321－1334.

［14］Fang Haiyan，Sun Liying，Qi Deli，et al.Using137Cs technique to quantify soil erosion and deposition rates in an agricultural catchment in the black soil region，Northeast China［J］.Geomorphology，2012，169／170：142－150.

［15］盛美玲，方海燕.WaTEM／SEDEM模型及其應用研究進展與展望［J］.地理科學進展，2014，33（1）：85－91.

［16］孫和平.小型水庫泥沙淤積成因分析及淤積量的初步估算［J］.地下水，2005，27（3）：221－222.

［17］魏建兵，肖篤寧，李秀珍，等.東北黑土區(qū)小流域農(nóng)業(yè)景觀結果與土壤侵蝕的關系［J］.生態(tài)學報，2006，26（8）：2608－2615.

［18］翟偉峰，許林書.東北典型黑土區(qū)土壤可蝕性K值研究［J］.土壤通報，2011，45（5）：1209－1213.

［19］章文波，付金生.不同類型雨量資料估算降雨侵蝕力［J］.資源科學，2003，25（1）：35－41.

［20］于寒青，李勇，Nguyen M L，等.基于FRN技術的我國不同地區(qū)典型土壤保持措施的有效性評價［J］.核農(nóng)學報，2012，2（2）：340－347.

［21］Verstraeten G，Poesen J，Gillijins K，et al.The use of riparian vegetated filter strips to reduce river sediment loads：An overestimated control measure？［J］Hydrological Processes，2006，20：4295－4267.

［22］方海燕，蔡強國，李秋艷.產(chǎn)沙模數(shù)與流域面積關系研究進展［J］.地理科學進展，2008，27（6）：63－69.

［23］Dong Y F，Wu Y Q，Zhang T Y，et al.The sediment delivery ratio in a small catchment in the black soil region of Northeast China［J］.International Journal of Sediment Research，2013，28（1）：111－117.

［24］簡金世.松花江流域不同侵蝕類型區(qū)泥沙輸移比的估算［D］.西安：西北農(nóng)林科技大學資源環(huán)境學院，2011.