接坝区不同面积小流域泥沙来源分析

许晶晶,杜鹏飞,屈丽琴,赵 莹,黄东浩,周瑞鹏

(1.中国水利水电科学研究院 国际泥沙研究培训中心,北京 100048;2.沈阳农业大学 水利学院,沈阳 110866)

土壤侵蚀会影响土壤健康和生产力[1]，引发泥石流和山体滑坡，导致水库河道淤积、水体面源污染，造成土壤生物多样性丧失[2]，从而减少缓解和适应气候变化的能力[3]。由土壤侵蚀形成的泥沙因携带了包括土壤养分、重金属、痕量金属和氧化物等在内的源地属性[4]，因此，不同类型物源区的各属性含量之间可能会存在差异。对沉积泥沙的溯源不仅能够为明晰土壤侵蚀、搬运、沉积过程提供依据，而且可以助力水土保持措施的合理布设[5-6]。

接坝区作为内蒙古高原向河北北部山地的缓冲过渡带，是京津上游地区的重要水源地和环境脆弱区，区内地貌沟深、坡陡，与紧邻其北部的坝上高原形成了鲜明对比[7]。由于过度开垦和放牧，接坝区的生态环境受到了严重破坏，成为21世纪初京津风沙的主要源区之一，曾一度威胁着首都地区的生态安全。位于接坝区中部的小坝子乡，地处半干旱农牧交错带，生态环境十分脆弱，曾是京津浮尘和扬沙天气的重要沙源地[8]。近年来，在大规模退耕还林还草治理下，区内风沙危害显著降低，生产生活环境明显改善，反而由水力侵蚀引发的耕层土壤变薄、侵蚀沟扩张、地表破碎化等有加剧的趋势，河道出现了不同程度的泥沙淤积。小坝子乡小坝子村修建于2013年的四级谷坊，历经短短7 a，其中两级已基本淤满，其严重程度可见一斑。

已有研究显示，接坝区不同面积小流域的土壤侵蚀存在很大差异[9]。在小流域尺度上，丘陵和低山顶部以风蚀为主，侵蚀沟沟谷以水蚀为主，位于中下坡位的耕地则受风水两相侵蚀的影响；在微小流域尺度上，坡顶是以风蚀和水力面蚀为主的复合侵蚀，坡面中下部则以风蚀和水力沟蚀较为常见。为了对比这两个不同面积小流域的泥沙来源，本研究以整个小坝子乡小流域和小坝子村的一处微小流域为例，对沉积泥沙的主要源地及其贡献比例进行分析。研究有助于加深对接坝区土壤侵蚀和泥沙输移过程的理解与把握，对基于侵蚀预防的水土保持措施布设和土地利用类型调整具有非常重要的参考意义。

1 研究区概况

小坝子乡位于内蒙古高原和燕山山地连接的接坝地带，土地总面积为325.40 km2，流域土壤以棕壤为主[10]。作为一个自然地理单元，小坝子乡流域(116°12′49″—116°29′30″E，41°22′8″—41°34′6″N)(简称乡流域，图1A)内最高点海拔2 019 m，最低点海拔874 m，地形起伏较大。行政区域上，南与窟窿山乡接壤，东与土城镇相连，北面和西北面分别为坝上的四岔口乡和大滩镇。该区属于典型的大陆性季风型半干旱气候，分属接坝冷凉气候类型区。年均气温4℃，年均无霜期105 d，多年年均降水量409.30 mm。区内年均降雨量时空分布差异较大，6—9月份的降雨占全年总量的80%，且多以暴雨形式出现。

作为小坝子乡流域内的一个微小流域，小坝子村流域(简称村流域，图1B)位于小坝子村西口，面积为0.164 km2。该流域坡地曾分布有农田，因沙化严重于20世纪90年代末退耕。目前水力侵蚀严重，主沟道总长为242 m，最宽处达29 m，最深处32 m，侵蚀沟扩张迅速，溯源侵蚀显著。

图1 研究区概况和样品采样点分布

2 材料与方法

2.1 样品采集

在夏秋降雨，尤其是个别场次暴雨的作用下，面上表土流失和沟道扩张发育同时出现，是本区土壤水力侵蚀的一个重要特征。对乡流域而言，由于京津风沙源治理工程的持续推进，当地政府积极进行退耕还林还草，使得研究区林地、草地分布最广，面积占比高达九成。流域上游坡顶为草地，坡中为成片的桦木林，总体覆盖度较高。中下游林地和草地呈块状分布，但总体覆盖度较低，中坡位和下坡位有农田分布，表土流失严重。加之沟道内坍塌的堆积物在夏秋雨季被冲刷，使得底土流失也不容小觑。为全面反映泥沙来源情况，将农田表土、林草地表土及侵蚀沟底层物质确定为该流域的潜在泥沙源区。对村流域而言，坡顶为沙化较为严重的草地，坡中及坡脚为不同沙化程度的退耕农田。由于沙化严重，粗骨化较为明显，在暴雨冲刷下，村流域内沟道纵横，主沟道垮塌现象显著。在这种情况下，退耕农田表土、草地表土和侵蚀沟底层物质就成为村流域潜在的泥沙源区。根据上述分类，于2018年4月进行了样品采集。对乡流域，通过遥感影像研判解译和实地勘查，选择了有代表性的源区17处，其中包括农田5处、林草地9处、侵蚀沟3条。为降低各源区指纹因子的空间变异，增加样本的代表性，采样使用多点混合法，以20 m为半径，重复采集样品3个进行充分混合。对村流域，使用网格法采集潜在泥沙源区样品，共采集草地表层样11个，沙地表层样10个，侵蚀沟底层样13个。上述所有表层样品的采集深度均为5 cm，侵蚀沟底层样品的采集部位均在沟壁距表层高度超过30 cm的亚表层和母质层。沉积样品均在流域出口处采集，沿主河道每隔50 m(乡流域)和5 m(村流域)采集沉积表层样，采集个数均为4个，采集深度约为5 cm。考虑到泥沙输移的连通性，样品的采集区相对靠近沟道，这些位置发生侵蚀后，泥沙更易搬运至河道直至沉积区。具体样品采集信息见表1。

表1 两个流域样品采集信息

2.2 样品分析

样品采集后，经自然风干、手工剔除杂质、研磨后，分别过2 mm和0.063 mm的筛子用于样品测试。测试工作在中科院寒区旱区环境与工程研究所完成。其中，小于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，Fe2O3，K2O，Na2O，MgO，CaO。

2.3 研究方法

本研究采用基于指纹因子测定的指纹识别技术来计算不同源区的泥沙来源及其贡献比。指纹识别技术通过对比物源和泥沙的指纹因子关系来实现对不同源地泥沙贡献的定量估算，已在我国黄土高原[11-15]、长江上游[16-19]、黔中喀斯特高原[20]、东北黑土区[21]、南方红壤区[22-23]等多地得到了成功应用。这一方法筛选指纹因子时主要包括以下3个步骤：

(1)双边范围检验。通过比较源区的样品和沉积区样品的元素浓度进行双边范围检验，当沉积区样品元素浓度大于源区样品元素浓度的最大值或者小于最小值时，将其剔除，剩下的即为保存性指纹因子。

(2)Kruskal-Wallis H检验。将通过双边范围检验的指纹因子进行Kruskal-Wallis H检验，评价单个指纹因子在不同泥沙源区中是否存在显著差异。p<0.05表示单个指纹因子在不同物源区间差异显著，该指纹因子通过Kruskal-Wallis H检验，可用于判定泥沙来源。

(3)多元逐步判别分析。将上一步骤中筛选出的指纹因子进行逐步判别分析，每一步进入判别函数的都是Wilks′ λ的计算值最小的指纹因子。将所有进入判别函数的指纹因子进行组合，即为最佳指纹因子组合。

基于指纹识别计算泥沙来源的模型主要包括多元线性混合模型和贝叶斯模型。本研究采用的是由Collins等[24]提出的多元线性混合模型，其计算可靠性和准确性已得到了有效验证[25-26]。

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

(1)

式中:n是泥沙源区个数;Ci是沉积泥沙中指纹因子i的浓度;m是指纹因子个数;Ps是泥沙源区s的泥沙贡献率;Ssi是泥沙源区s中指纹因子i的平均浓度。

利用拟合优度算法检验模型计算结果的不确定性，见公式(2)：

(2)

式中:MAF是平均绝对拟合度值,一般认为MAF>0.8时,可以接受多元混合线性模型的结果。

3 结果与分析

3.1 不同面积小流域泥沙粒径特征

在乡、村两个典型流域中，不同土地利用类型的泥沙颗粒组成差异明显。由图2A可知，在乡流域，与农田样品和侵蚀沟样品相比，不同林草地样品的粒径变化较大，采集于地表覆盖厚达15 cm枯枝落叶层的桦木林表土中的砂粒含量最低，仅为20.27%，采集于河滩、沙化严重、基本没有植被覆盖地表的次生杨树林表土样品砂粒含量最高，达93.14%。农田样品粉粒和黏粒含量是所有源区样品中最高的，为40.93%，8.63%，较侵蚀沟和林草地分别高了16.14%，27.63%，2.25%，5.90%，这可能反映了农田在长期耕作过程中下翻细粒土壤融入耕层的显著影响。

由图2B可知，村流域物源区样品砂粒含量均超过70%，退耕农田砂粒含量最高，达96.81%，粗粒化最明显。对比图2A和2B，村流域退耕农田与侵蚀沟砂粒含量分别为96.81%，91.10%，比乡流域农田和侵蚀沟分别高出46.37%，7.14%。而沿主沟道向流域外延的沉积样品砂粒含量逐渐降低，黏粒和粉粒含量逐渐增加，说明在搬运过程中，粗颗粒泥沙先沉积，而细颗粒泥沙则能够被搬运得更远。

图2 两个流域泥沙粒径百分比

由图3可见，两个流域物源区的中值粒径中，林草地基本一致，侵蚀沟略有差异，农田和退耕农田差异较大。村流域退耕农田的中值粒径为190.40 μm，是乡流域农田中值粒径的7.40倍。沉积区泥沙的中值粒径则基本相等。在乡流域内，中值粒径由沉积区、侵蚀沟、林草地和农田的顺序依次减小，在村流域内，则略有不同，按沉积区、退耕农田、侵蚀沟和草地的顺序依次减小。不论是乡流域，还是村流域的中值粒径差异均较为明显地反映出不同土地利用对土壤或泥沙颗粒组成的影响。沉积区泥沙中值粒径最大，这可能与泥沙输移的分选性有关，即被搬运至流域外的细粒泥沙比例较高而在流域内沉积的粗粒泥沙较多。也可能与泥沙连通性有关，本研究结果则直观反映了沉积区与较大中值粒径源区距离对沉积的影响，即在乡流域，沉积区与侵蚀沟的距离最近，而在村流域，沉积区与退耕农田的距离最近。但同时也应该认识到，由于土壤侵蚀的发生具有显著的时空变异特征，泥沙连通性理论上也会出现类似的变化，加之其在反映土壤侵蚀严重程度上存在滞后性，所以不论对于单次侵蚀过程，还是一定时段的累积过程而言，准确刻画泥沙连通性的变化都较为困难。

图3 两个流域中值粒径对比

3.2 最佳指纹因子组合确定

对测定的30个指纹识别因子进行分析。第一步，利用双边范围检验筛选出两个流域中分别有25个和29个保存性指纹因子；第二步，通过非参数Kruskal-Wallis H检验筛选在不同物源区间有显著差异的指纹因子(p<0.05)，两个研究区分别有7个和25个潜在指纹因子通过检验；第三步，通过逐步判别分析确定最佳指纹因子组合。各级指纹因子筛选结果见表2。由于村流域的地形更简单，泥沙在搬运、沉积过程中更为稳定，所以三级筛选中，村流域通过检验的指纹因子个数均比乡流域多。

表2 两个流域各级筛选后指纹因子对比

两个不同面积流域的逐步判别分析结果见表3。乡流域以K2O，Cl，Ce为最佳指纹因子组合，能以85.7%准确度判别该研究区的农田表土、林草地表土和侵蚀沟3个泥沙源区。每步累积判别率为61.9%，66.7%，85.7%，单步骤判别率分别为61.9%，33.3%，57.1%。村流域以Ba，P，CaO，SiO2这4个因子作为最佳指纹因子组合，能以89.5%准确度判别该研究区的退耕农田表土、草地表土和侵蚀沟3个泥沙源区。每步累积判别率为60.5%，63.2%，78.9%，89.5%，单步骤判别率分别为60.5%，65.8%，57.9%，63.2%。

表3 两个流域的逐步判别分析结果

3.3 不同面积小流域泥沙来源分析

3.3.1 不同源区的泥沙相对贡献率 根据最佳指纹因子组合中各因子的物源区和沉积区泥沙样品中的浓度，利用Walling-Collins模型计算不同物源区的泥沙相对贡献率。农田表土和退耕农田表土分别是乡流域和村流域主要的泥沙源区。在乡流域，农田表土对沉积泥沙的贡献比高达71.6%，是第二源区林草地表土(15.9%)的4.5倍，第三源区侵蚀沟仅贡献了12.5%的泥沙。在村流域，退耕农田表土、侵蚀沟和草地表土的泥沙贡献比差距相对小一些，分别占比44.4%，34.2%，21.4%。

为了检验计算结果的稳定性，使用平均绝对拟合度值(mean absolute fit，MAF)评估模型(公式4)对计算结果进行了评估，乡流域和村流域两个研究区的泥沙来源的MAF分别为0.97，0.98，远高于0.8的临界可接受值，表明Walling-Collins模型对两个流域泥沙贡献率的计算结果是可靠的。

3.3.2 不同源区的泥沙面积比 为明确不同土地利用类型在流域内的面积占比与泥沙贡献比的关系，利用样品采集当年(2018年)的遥感影像对流域土地利用类型进行了解译。由于小坝子乡土地利用类型比较复杂，解译过程中参照《区域水土流失动态监测技术规定》，依据土地利用二级分类标准进行，解译结果见图1和表4。

由图1和表4可知，乡流域土地利用类型以林地和草地为主，占总面积的90.21%，其次为农田，占比6.47%，侵蚀沟占比2.13%，剩余为建设用地等其他用地。村流域的退耕农田和草地分别占比41.38%，57.26%。根据现场调查，结合无人机影像，确定的侵蚀沟面积为0.014 2 km2，面积占比8.64%。

表4 两个流域的土地利用类型面积及占比

为反映各源区泥沙输移沉积的严重程度，引入泥沙面积比来定量表达各侵蚀源区泥沙搬运至流域出口沉积区的比例和其所占流域面积的比例，定义为这两个比例的比值。当比值结果大于1时，表明流域内该类型物源区流失土壤的比例偏高，相应地，土壤侵蚀强度也越高；反之，表明流域内该类型物源区流失土壤的比例不是很严重，值越接近于0，土壤流失后搬运至沉积区的比例越少。

由图4可见，除乡流域的林草地和村流域的草地表土泥沙面积比小于1外，其他4个侵蚀源区的泥沙面积比均高于1。尤其是乡流域的农田表土泥沙面积比高达11.07，其土壤侵蚀、搬运的比例远高于村流域退耕农田表土泥沙面积比的1.07。一方面，这可能与农田在人为扰动后，表土变得疏松，从而更易于流失有关。每年的雨季也恰是农田的除草季，相比于退耕农田，农田的人类活动干扰更加强烈，由此导致表土更容易被侵蚀掉。另一方面，从农田在乡流域的分布来看，大多集中在下坡位，距离沟道较近，受水系连通性强的影响，表土侵蚀后更容易被水流搬运至河道中，从而其在流域出口沉积的概率也更高。

图4 两个流域的源地泥沙面积比

乡流域和村流域的侵蚀沟泥沙面积比分别为5.87，3.96，表明该研究区域的侵蚀沟不断扩张发育导致的底层土壤流失是不容小觑的。侵蚀沟的快速扩张使坡面径流和泥沙更加容易在沟道内汇聚，从而增大了流域内泥沙输移的连通性，间接加剧了土壤侵蚀后的输移比例。乡流域的林草地和村流域的草地泥沙面积比分别只有0.18，0.43，与林草地植被覆盖度总体较高，能够增加泥沙拦截，从而降低侵蚀搬运量有关。

从上述结果来看，村流域后续的水土流失治理重点应着眼于侵蚀沟，乡流域除了侵蚀沟外，还要重视农田的水土流失防治，实施以工程措施、生物措施和耕作措施相结合的综合治理模式。具体来说，在主沟道中修建淤地坝，支沟中修建梯级谷坊，进行分段控制，并在沟头植树种草减缓溯源侵蚀。坡耕地应减少雨季扰动的频次和强度，通过免耕、等高耕作等措施降低水力侵蚀的危害。

4 结论与讨论

本研究在接坝区选取了小坝子乡和小坝子村两个不同面积的小流域定量计算了不同源区的泥沙来源比例。对研究区样品的粒径分析发现，同流域内不同源区间泥沙粒径差异明显，不同流域同一源区的泥沙粒径差异也很大，但两个流域沉积泥沙粒径较为接近。乡流域农田表土、林草地表土和侵蚀沟3种源区泥沙来源的相对贡献率分别为71.6%，15.9%，12.5%；村流域农田表土、草地表土和侵蚀沟3种源区泥沙来源的相对贡献率分别为44.4%，21.4%，34.2%。乡流域泥沙面积比最大的是农田，比值达11.07，其次为侵蚀沟5.87，最小的林草地仅为0.18；村流域中泥沙面积比由大到小分别为侵蚀沟3.96，退耕农田1.07，草地0.43。研究结果表明，侵蚀沟和农田是未来水土流失防治的重点区域。

本研究明晰了接坝区不同面积小流域泥沙来源的基本特点，相关结果可为该区水土流失治理和土地利用规划编制提供参考。但同时应该注意到，由于接坝区内海拔高差大，受降雨、植被、泥沙连通性及人为扰动强度的影响，不同流域的侵蚀产沙和输移沉积也会存在一些差异。在今后的研究中应综合考虑这些因素的综合影响，以便深入比较不同条件下各尺度流域的泥沙来源情况，分析不同空间尺度逐级递进对计算结果的影响，从而为泥沙来源跨尺度衔接计算提供参考。