土壤通过入渗、蒸散发、地表径流和地下水补给等途径调节水分平衡，从而对区域和全球地表水分和能量平衡产生重大影响^[1-2]。土壤饱和导水率(soil saturated hydraulic conductivity，K_s)是研究地下水产流机制、径流发生以及土体可蚀性的重要参数。传统获取K_s值是通过原位测量或野外取样进行室内测定的方式获得，所获数据准确性高，但该方法成本高、耗时耗力，同时由于K_s较大的空间变异性，所测数据很难应用到大区域范围内^[3]。1989年荷兰土壤学家Bouma^[4]首次提出传递函数(pedo-transfer functions，PTFs)的概念，旨在通过大量数据建立土壤易获取性质(有机质、质地、容重等)与土壤水力学参数之间的关系，间接计算土壤水力学参数。传递函数的提出弥补了大量土壤水力学参数的缺失值，大大促进了土壤水分运动、溶质迁移和土壤侵蚀的研究^[5]。近半个世纪，尤其是近30年来，科研工作者围绕K_s传递函数做了大量的研究工作，取得了丰硕的研究成果，也提出了基于不同因子的传递函数。Campbell和Shiozawa^[6]根据英国大量的土壤样本数据，将土壤容重和颗粒大小分布作为因子构建土壤传递函数用于预测K_s；Wang等^[7]利用青藏高原382个土壤数据，将土壤粒径分布、容重、有机碳含量以及首次加入纬度作为输入因子建立了K_s传递函数；同样，Zhang和Schaap^[8]根据1 306个土壤数据，采用神经网络方法分别建立了将土壤颗粒分布、容重、田间持水量作为输入因子的多个传递函数并用于预测K_s。然而，影响K_s的因素极其复杂，其中土壤类型、质地、结构以及土壤大孔隙是影响K_s的主要因素^[9-11]，同时K_s还受区域气候以及土地利用类型的影响^[12-13]，导致K_s呈现极大的时空变异性。因此，在采用土壤传递函数预测K_s时，需要进行区域验证。

横断山地区位于四川盆地西部、云贵高原西北部和青藏高原东部，处于中国第一、二级地势阶梯的过渡部位，属高山峡谷地貌。区内地形地貌、气象水文、地质构造等地质环境条件复杂，土壤侵蚀严重，出露的岩层主要有三叠系上统灰白色砂板岩夹少许灰岩，第四纪松散沉积物质组成为古河床沉积形成的砂、卵石层，岩体较为破碎，风化强烈^[14]，在暴雨季节，崩塌、滑坡、泥石流爆发频繁，严重威胁居民生命安全及重大工程建设^[15-16]。K_s是横断山地区土壤水分运动、自然灾害预测的重要参数^[17]。然而，关于横断山地区K_s的相关研究鲜有报道，且目前在横断山地区水力学参数指标的可用数据非常有限^[18-20]。战海霞^[21]的研究发现，横断山地区土壤多为林地土壤，土体结构松散，根系发达且石砾含量高，土壤大孔隙结构丰富，由此导致的土体渗透系数高达25 mm/min。因此，现有的K_s传递函数在该类土壤K_s参数估算上有待验证。基于此，本研究结合研究区土壤含石率偏高的特点，增加土壤石砾含量(粒径 > 2 mm)作为输入变量，同时选取土壤容重、有机质含量和颗粒分布3种土壤基本特性参数作为输入变量，拟构建本区域K_s传递函数模型。

1 材料与方法 1.1 研究区概况

本研究以贡嘎山海螺沟流域(29°20′ ~ 30°20′ N，101°30′ ~ 102°15′ E)为研究区。贡嘎山位于青藏高原东缘，主峰海拔7 556 m，是横断山地区的最高峰，区域内地势总体呈东西两侧和北面高、中部和东南部低的趋势，岭谷高差一般1 500 ~ 3 000 m，各支沟狭窄，多呈“V”型谷，谷坡坡度多35° ~ 55°；在支沟和主沟交汇处相对开阔，坡度较支沟两岸处低缓，为25° ~ 35°^[22]。该区域属于亚热带季风气候区，年均气温4.1 ℃，年均降水量1 903 mm，雪线和林线分别位于海拔4 900 m和3 700 m附近^[23]。海螺沟流域属于典型山地生态系统，垂直气候带完整，随海拔升高，依次为亚热带、暖温带、寒温带、亚寒带、寒带等不同气候类型^[24]。同时随着气候类型不同，形成明显的植被带，其中，1 900 ~ 2 200 m是阔叶林带，年降水量约2 000 mm，土壤类型为棕壤，植被以苞槲柯(Lihocarpus clestocarous)为优势种；海拔2 200 ~ 2 800 m为针阔混交林带，年降水量2 500 mm，土壤类型为暗棕壤，植被有铁杉(Tsuga chinensis)、槭树(Acer)、桦树(Betula)等组成；海拔2 800 ~ 3 600 m为针叶林带，年降水量2 800 ~ 3 100 mm，土壤类型为山地灰化土，主要树种以麦吊杉(Picea brachytyla)和峨眉冷杉(Abies fabri)为主。研究区土壤均具有粗骨性强、土薄、石砾含量丰富的基本特征^[25-26]。

1.2 样品采集与分析

于2019年8月在海螺沟流域开展野外调查，在常绿阔叶林、针阔混交林、针叶林3种植被带各选取3块代表性样地，除去地表凋落物后，每个样地挖一个1 m深的土壤剖面。每个剖面分四层(0 ~ 10、10 ~ 20、20 ~ 50和50 ~ 100 cm)，用环刀(高度5 cm，直径5 cm，体积100 cm³)取原状土样用于测定土壤容重、石砾含量和K_s；并在环刀周围0 ~ 20 cm范围内采集土样，带回实验室自然风干后研磨分别过0.25 mm和2 mm筛孔用于土壤有机质和颗粒组成的测定。K_s利用定水头下的马氏瓶渗透仪进行测定^[27]，容重采用环刀烘干法测定^[28]，石砾含量以环刀内土样的全部石砾(粒径 > 2 mm)与环刀内全部土样的体积百分比计算，土壤有机质采用重铬酸钾外加热法测定^[29]，土壤颗粒组成采用吸管法测定^[30]，土壤颗粒粒级按照美国制土壤颗粒组成分组，分为3级：0.05 ~ 2 mm砂粒，0.002 ~ 0.05 mm粉粒和 < 0.002 mm黏粒。

1.3 研究方法

本研究利用Excel 2010及SPSS 22.0对研究区所测得土壤参数进行基本统计分析，并将研究区实测数据分别代入不同土壤传递函数模型中进行计算，比较估算值与实测值的误差，采用平均误差(ME)、均方根误差(RMSE)、决定系数(R²)评价传统土壤传递函数的预测精度。在此基础上，选取土壤容重、石砾含量、有机质含量、土壤粒径组成作为输入参数，利用相关分析、主成分分析和多元逐步回归分析研究K_s的主要影响因素，建立横断山地区森林土壤K_s传递函数。

本研究所选取的9种经典传递函数包括Weynants模型^[31]、Cosby模型^[32]、Wosten模型^[33]、Campbell模型^[6]、Pucket模型^[34]、Vereecken模型^[35]、Saxton模型^[36]、Julià模型^[37]、Ahuja模型^[38]，其具体的函数表达式见表 1。

新建传递函数的具体构建步骤如下：①建立多元回归模型Y=a + b×x₁ + c×x₂ + d×x₃(a、b、c分别为各自变量x₁、x₂、x₃的回归系数)；②回归系数显著性检验，将不显著变量剔除，再次进行回归；回归方程显著性F检验，确保P < 0.05；③模型检验：进行模型多重共线性检验，如果变量之间存在共线性，则剔除该变量；自相关检验，消除自相关带来的误差；④重复以上步骤① ~ ③，直至模型中各自变量与因变量达到统计学意义上的显著相关。

模型评价指标的计算：

$ {\text{ME = }}\frac{{\sum {{{{\text{(}}{K_{\text{s}}} - {K_{\text{e}}}{\text{)}}}^{\text{2}}}} }}{n} $

(1)

$ {\text{RMSE = }}\sqrt {\frac{{\sum {{\text{(}}{K_{\text{s}}} - {K_{\text{e}}}{\text{)}}} }}{n}} $

(2)

$ {R^{\text{2}}}{\text{ = }}\frac{{{\text{cov(}}{K_{\text{s}}}{\text{, }}{k_{\text{e}}}{\text{)}}}}{{{\text{var(}}{k_{\text{s}}}{\text{)var(}}{k_{\text{e}}}{\text{)}}}} $

(3)

式中：K_s为实际测量的土壤饱和导水率(mm/min)；K_e为模型预测的土壤饱和导水率(mm/min)；n为测量的数据点的数量；ME反映了估算值对实测值的平均偏离程度，ME值越接近0，模型拟合效果越好；RMSE反映估算值与实测值的吻合程度，RMSE值越小，则越接近于实测值；R²反映了估算值与实测值在变化趋势上的相似性，R²越接近1，说明估算值与实测值在变化趋势上越相似^[39]。

2 结果与讨论 2.1 研究区土壤理化性质

研究区属于典型高山森林系统，不同土层的基本理化性质如表 2所示，可见，不同土壤深度土壤理化性质变化较为明显，随土壤深度增加，容重和土壤黏粒含量逐渐增大，表层土壤容重范围介于0.69 ~ 0.98 g/cm³，这是由于森林表层土壤因根系丰富，土壤较疏松，随着土壤深度的增加，土壤结构越为紧实，导致较深层次土壤容重增加^[40]。表层(0 ~ 10 cm)土壤凋落物丰富，有机质含量平均含量为46.59 g/kg，平均含石率为2.77%；随土壤深度增加，含石率逐渐升高，而有机质含量逐渐降低，50 ~ 100 cm土层有机质含量仅为15.33 g/kg，而含石率增加到19.33%。从土壤颗粒组成看，砂粒是主要的粒级，占比均介于77.81% ~ 97.27%，且表现出随土层深度的增加逐渐增加的趋势。变异系数可反映实测数据与平均值之间的差异性，可表示数据的离散程度，用实测数据标准差与平均值的比来表示。变异系数 < 0.15属于弱变异性，0.16 ~ 0.35属于中等变异性，> 0.36表示较强的变异性^[41]。由表 2可知，不同土层容重变异系数 < 0.35，在空间变异上属于弱变异性；表层土壤有机质含量变异程度较小，深层土壤有机质变异系数 > 0.36，变异程度较大；土壤表层含石率变异系数为0.80，变异程度较大；全土层土壤砂粒变异系数范围介于0.02 ~ 0.07，变异程度最小。

2.2 经典土壤饱和导水率传递函数在研究区的适用性

本文选用表 1所列举的9种经典土壤传递函数，对研究区土壤K_s值进行估算。表 3为应用9种K_s传递函数计算所得的估算值与实测值的RMSE、ME和R²。结果发现，9种传递函数ME值均大于10，表明估算值大大低估了研究区实测K_s值，无法达到预测效果，Ahuja模型预测效果优于其他8种传递函数，但仍然达不到预测精度；9种传递函数RMSE值介于16.60 ~ 24.46 mm/min，预测精度均较差，Ahuja函数RMSE值最小为16.60，预测精度高于其他8种传递函数，但预测效果仍属较差。从估算值与实测值的R²看，9种经典传递函数预测精度均很差，即使相对略好的Saxton模型和Ahuja模型的R²值分别为0.51和0.31，仍无法实现准确预测。从图 1也可以看出，Weynants、Cosby、Wosten、Campbell、Pucket、Saxton和Julià模型的估算结果均位于1∶1线底部，严重低估实测值，无法达到预测效果。Vereecken和Ahuja模型结果虽分布于1∶1线两侧，但较1∶1线偏离程度较大。

从输入变量的构成来看，Cosby、Pucket和Julià传递函数仅将土壤颗粒组成作为输入变量来估算K_s，适用于低海岸平原地区黏粒含量大于20% 且大孔隙不发达的土壤^{[32，34，37]}；Campbell传递函数在考虑土壤粒径分布的基础上加入容重作为输入变量，构建K_s传递函数，适用于美国北部土壤类型^[6]；Rawls等^[42]指出土壤有机质含量是影响K_s的重要因素，Weynants、Wosten和Vereecken传递函数在考虑土壤粒径分布和容重的基础上，将有机质含量作为输入变量构建传递函数；Saxton模型则加入土壤饱和含水量作为因子；Ahuja模型是以孔隙度为唯一变量构建的模型。然而，目前以传统影响因子构建的传递函数均大大低估了本研究区土壤K_s。究其原因，有两个方面，其一，横断山地区高山森林生态系统植物根系发达，容易形成土壤大孔隙，为土壤优先流提供通道^[43]，大大改变土壤K_s，然而，目前经典传递函数均没有考虑土壤大孔隙对土壤K_s的影响；其二，由坡积物发育的研究区土壤的石砾含量高，而石砾是形成大孔隙的重要因素之一^[44]，因此，此类土壤大孔隙极为丰富。土壤大孔隙是影响K_s的重要因素，关于土壤大孔隙与K_s之间的关系已有大量报道，王金悦等^[45]研究发现林地土壤大孔隙的总数量决定了K_s 66% 的变异，大孔隙的平均体积决定了K_s 79% 的变异；刘目兴等^[46]发现半径 > 1 mm的大孔隙数量决定了K_s 85% 的变异。显然，大孔隙的数量与体积对土壤水分的运移与分布有着重要的影响，大孔隙的数量越多、体积越大，土壤K_s值越大。综上所述，本研究区植物根系发达，土壤石砾含量高，导致土壤大孔隙发达，土壤导水能力极强，故而采用传统经典传递函数会大大低估本区的K_s值。因此，本研究将土壤石砾含量作为新增的输入变量，考虑由石砾作用所造成的土壤优先流对土壤K_s的影响，从而构建适合本研究区的K_s传递函数，用于预测土壤K_s值。

2.3 研究区土壤饱和导水率与基本理化性质的相关性

土壤K_s与理化性质的相关性分析结果如图 2所示，可见，K_s与含石率(GF)在P < 0.01水平下显著负相关，与容重(BD)在P < 0.05水平下显著负相关，表明土壤含石率和容重是影响横断山区森林土壤K_s的主要限制因子。K_s与有机质含量(SOM)、黏粒含量(clay)的相关性不大，表明横断山区土壤K_s受有机质含量和土壤颗粒组成影响不显著。

2.4 研究区土壤饱和导水率传递函数的构建

目前，国内外K_s传递函数主要采用容重、有机质含量、土壤孔隙度以及土壤粒级分布作为输入变量^[32-36]，经验证，并不适用于横断山研究区。在石质山区，只针对土壤部分展开研究而不考虑石砾的作用是不全面的，横断山地区土壤属森林土壤，具有较高的石砾含量、发达的根系以及丰富的大孔隙，这些特点对土壤入渗、蒸发、径流等水文过程产生重要影响^[47]。一方面，王慧芳等^[48]的研究表明，土壤中存在的石砾有显著增加K_s的作用，且碎石含量越高，增幅越明显，导致传统K_s传递函数对含石砾土壤的适用性很差^[49]。另一方面，当土壤中石砾含量较高时，石砾与土壤交界处更容易形成大孔隙，而土壤中大孔隙的数量及连通性是影响K_s的主要因素，当石砾含量增加到某一范围时，大孔隙之间的连通性会显著提高，形成有效的大孔隙通道，发生优先流，水分移动速度加快，K_s发生明显的跃升^[50]，Xu等^[51]等研究发现，当石砾含量 > 15% 时，K_s随着石砾含量的增加呈现显著增加趋势。本文在研究石砾含量对K_s的影响时，重点关注土壤中 > 2 mm的石砾^[52]，研究成果可为石砾含量丰富的森林土壤K_s研究提供帮助。针对此类林地土壤，本研究选择土壤含石率(GF)、容重(BD)、有机质含量(SOM)、黏粒含量(clay)以及粉粒含量(silt)作为初始输入变量，利用多元逐步线性回归分析法分析不同变量组合构建模型的AIC值，AIC鼓励数据拟合的优良性同时尽量避免出现过度拟合的情况。所以优先考虑的模型应该是AIC值最小的组合。由表 4可知，其他变量不变的情况下，减少土壤容重作为输入变量增加模型AIC值，同时降低R²；加入土壤含石率明显降低模型AIC值，同时提高模型R²。选择土壤容重、含石率、黏粒含量和粉粒含量作为输入变量，构建传递函数，模型AIC值最低同时R²最大，效果最佳。本研究新建的K_s传递函数表示为：K_s = 9.48 + 12.32×BD + 0.29×SOM – 1.94×GF + 2.89×silt –5.34×sand，经检验，新建传递函数ME值为5.13 mm/min，RMSE值为7.04 mm/min，预测误差明显低于9种经典传递函数，R²达0.67(n=40)(表 3)。从图 1可以看出，新建传递函数预测值和实测值分布于1∶1线两侧，预测效果较好。

杜阿朋等^[53]研究表明，六盘山叠叠沟小流域森林土壤石砾中最多的是粒径在2 ~ 4 mm和4 ~ 6 mm的破碎石砾，二者之和占总石砾体积的62.39%，在该区研究石砾对水文的影响时应将重点放在2 ~ 6 mm的石砾上。喻明美等^[54]研究表明，环刀尺度测量K_s在白云山5种不同类型含石砾林地土壤具有一定的代表性。因此，充分兼顾理论和实际可操作性，本研究基于环刀尺度样品数据构建的K_s传递函数具有可靠性和可行性，基本适用于本研究区土壤。

3 结论

研究区土壤石砾含量高，大孔隙丰富，导水能力强，表层土壤K_s值最高可达57.01 mm/min，经验证，现有9种经典传递函数对该区的K_s预测值与实测值偏差极大，不在同一个数量级上，大大低估了实测值，不适用于该区域森林土壤K_s预测。本研究增加土壤石砾含量(粒径 > 2 mm)作为输入变量，新建传递函数：K_s = 9.48 + 12.32×BD + 0.29×SOM – 1.94×GF + 2.89×silt – 5.34×sand，预测值与实测值相关系数为0.67，可获得较为准确的K_s估算值，该模型可作为研究区K_s传递函数使用。