在夏秋降雨,尤其是个别场次暴雨的作用下,面上表土流失和沟道扩张发育同时出现,是本区土壤水力侵蚀的一个重要特征。对乡流域而言,由于京津风沙源治理工程的持续推进,当地政府积极进行退耕还林还草,使得研究区林地、草地分布最广,面积占比高达九成。流域上游坡顶为草地,坡中为成片的桦木林,总体覆盖度较高。中下游林地和草地呈块状分布,但总体覆盖度较低,中坡位和下坡位有农田分布,表土流失严重。加之沟道内坍塌的堆积物在夏秋雨季被冲刷,使得底土流失也不容小觑。为全面反映泥沙来源情况,将农田表土、林草地表土及侵蚀沟底层物质确定为该流域的潜在泥沙源区。对村流域而言,坡顶为沙化较为严重的草地,坡中及坡脚为不同沙化程度的退耕农田。由于沙化严重,粗骨化较为明显,在暴雨冲刷下,村流域内沟道纵横,主沟道垮塌现象显著。在这种情况下,退耕农田表土、草地表土和侵蚀沟底层物质就成为村流域潜在的泥沙源区。根据上述分类,于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 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,Al₂O₃,SiO₂,Fe₂O₃,K₂O,Na₂O,MgO,CaO。

本研究采用基于指纹因子测定的指纹识别技术来计算不同源区的泥沙来源及其贡献比。指纹识别技术通过对比物源和泥沙的指纹因子关系来实现对不同源地泥沙贡献的定量估算,已在我国黄土高原^[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多元线性混合模型如下:

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

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

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

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

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

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

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

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

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