水土保持功能是森林生态系统服务功能中最重要的一种。亚热带红壤区尽管森林覆盖率较高, 但由于地形复杂、雨量多而集中和人地矛盾突出^[1], 且过去几十年过度追求林业经营效益, 植物群落结构单一, 林下植被多样性丧失, 导致水土流失和森林加速退化^[2]。森林自调节周期较长, 治理初期敏感脆弱, 且随时可能发生反弹, 由此造成的一系列生态环境问题将持续制约当地社会经济发展^[3]。

受多重因素的综合影响, 森林水土保持功能在不同植被条件、土壤和地形情况下研究结果差别很大^[4]。植被主要通过垂直结构如林冠层、灌草层和枯枝落叶层截持调蓄降雨、降低雨滴动能, 从而影响林下水文过程和土壤侵蚀过程^[5]。不同植被类型的水土保持能力波动大, 稳定性低^[6]。热带、亚热带的阔叶林对降雨的调节能力较强, 达100 mm以上, 而暖温带、温带森林及亚热带针叶林的均低于100 mm^[7]。土壤层是森林存蓄水分的主要载体, 占林分截持水量的90%以上^[8], 其截留和持水能力主要受土壤结构、质地和孔隙度等因素的影响^[9-10], 而不同抚育措施也会使森林土壤和结构产生差异^[11], 针阔混交林的土壤团聚体含量、直径和有机质均优于针叶林^[12]。地形可影响降雨再分配过程和雨滴动能, 坡度坡长等因子是水土流失的重要影响因素之一^[13]。有研究表明, 存在临界坡度和坡长值, 水土流失量随着坡度和坡长和增加, 呈先增后减趋势^[14], 还有研究表明土壤养分流失与坡度和坡长存在正相关关系^[15]。综上, 对森林水土保持功能的定量研究已从传统的“植被覆盖度”单一指标向“乔灌草”多层次、枯落物、地形、土壤、气候等多因子多指标逐步深入发展, 但研究更多还是聚焦于各因子与水土流失现象间的关系, 较少关注各因子间的交互作用。

我国在水土保持功能的量化研究从黄土高原地区逐步发展而来^[16], 主要观测方法为“径流小区法”, 该方法将径流量和产沙量与各影响因子联系起来, 通过长时间的地面观测研究自然和人为因素对坡面水土流失的影响规律和人工干扰措施效用^[17]。目前国内观测积累对比国外仍有差距^[18], 模型研究方面主要利用RUSEL、LISEM、USEL等通用水土流失模型进行土壤侵蚀预测模拟以及敏感区域判定^[19-21], 降雨侵蚀力因子、土壤可侵蚀性因子和植被人工措施因子是主要参数, 但侵蚀模型应用对地形要求较高, 不同类型地区缺失的背景值如土壤年流失量、降雨侵蚀力因子和土壤可侵蚀性因子仍需要通过连续8 a以上的径流小区观测得到^[22]。所以, 径流小区仍是进行水土保持功能量化研究的重要手段, 但现今对观测调查数据的背后联系挖掘还不够, 对各因素间交互影响的量化研究还不够深入。

目前, 国内外研究水土保持影响因素间关系主要采用变异函数、相关分析、多元回归分析、主成分分析和灰色关联分析等统计学方法^[23-24], 此类方法对数据要求较高, 只能对连续的数值量进行分析且要求样本单元多于30个, 且其结果表示的是自变量对因变量的单一影响, 无法处理类型变量和分析变量间的交互作用。地理探测器(GeoDetector)是研究空间分层异质性以量化度量现象背后驱动因子影响力并揭示其影响机制的空间统计学模型, 它可以通过量化因变量与自变量之间空间分布相似性程度来表征自变量(X)对因变量(Y)产生的影响力, 适合分析类型变量也可以分析经过离散化的数值量, 既能分析地理空间数据也能分析属性空间数据, 且在探测自变量交互作用上存在优势^[25]。地理探测器的出现正好补充了经典统计学和传统地统计学方法存在的短板, 近年来得到了广泛使用。

当前, 探究土壤侵蚀过程和水土保持措施防蚀机理以及利用土壤侵蚀预测模型, 是水土保持科学领域亟需解决的重点问题^[26]。本文采用地理探测器统计分析方法, 结合径流小区关于降雨、地表径流、泥沙、土壤和植被等多种地面观测数据, 以对应侵蚀模型相关参数, 分析了在亚热带红壤区针叶纯林采取补植阔叶树种和施肥生态修复措施的水土保持效应, 找出了影响亚热带红壤区坡面水土保持能力的主要因素, 量化了各影响因素单独和交互作用后对坡面产沙产流的影响力, 为探究亚热带红壤区水土流失防治机理和土壤侵蚀预测模型的参数本地化提供科学基础。

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

研究区位于福建省龙岩市新罗区小池镇龙门溪小流域内(北纬24°47″—25°35″, 东经116°40″—117°20″), 地处九龙江北溪上游, 是龙岩市256万人的饮用水与工农业水源地, 区域生态安全意义重大。研究区以红壤土为主, 属亚热带海洋性季风气候, 平均海拔560 m, 年均温度18.4 ℃, 年均降雨量1768 mm, 日照时数1804 —2060 h, 全年气候温和, 无霜期长, 雨量充沛, 适合亚热带作物和林木的生长^[27]。近几十年来由于当地养殖业、种植业发展和采伐、开矿等人为活动干扰, 导致地表植被及下垫面遭到破坏, 水土流失严重, 生物多样性减少, 生态系统出现功能性失调。2011年遥感调查发现新罗区水土流失总面积达128.8 km², 约占土地总面积的4.8%^[27]。研究区内植被结构简单, 以马尾松(Pinus massoniana)和杉木(Cunninghamia lanceolata)中龄林为主, 水源涵养及水土保持功能低下。

1.2 试验地布设与处理

由于土壤养分含量较低, 在研究区内形成了近20 a林龄的以马尾松、杉木为造林树种的“老头松”和“老头杉”覆盖的针叶纯林。在2014年4月设置了11个坡度25°、中坡位和正北坡向的固定样方(15 m×20 m), 选择枫香(Liquidambar formosana)、深山含笑(Michelia maudiae Dunn)、木荷(Schima superba)、红叶石楠(Photinia x fraseri)、千年桐(Aleurites montana)等乡土阔叶树种, 在试验地补植幼龄(5 a)阔叶树和施肥处理。

对照Ⅳ区(后文简称针叶纯林)为针叶纯林, 对应样地10、11, 不作干涉；处理Ⅰ区(简称补植施肥3∶2林地)为补植阔叶树至针阔混交比例3∶2林地, 对应样地1、2、3, 采用穴施法^[28]对全部乔木施复合肥(N16:P16:K16), 穴长25 cm、宽20 cm、深30 cm, 穴距50 cm, 施肥穴统一选择于树干正上坡距其基部50 cm处, 每株施肥500 g, 随即覆土；处理Ⅱ区(简称补植3∶2林地)为补植阔叶树至针阔混交比例3∶2林地, 对应样地4、5、6；处理Ⅲ区(简称补植7∶2林地)为补植阔叶树至针阔混交比例7∶2林地, 对应样地7、8、9。

本研究采用样方收获法、土壤剖面法、环刀法和烘干法等调查取样方法^[29-30], 测量各样地的植被、土壤等相关数据(表 1)。

1.3 径流小区设置

参考相关标准^[31], 在11个样方中逐一布设标准径流观测小区, 小区投影面积为100 m²(20 m×5 m), 周围用PVC板围护并设有5—10 m的隔离防护带以消除区间干扰, 围埂应互相连接紧密, 埋深牢靠, 地表出露20 cm, 埋入地下20 cm。根据小池镇50 a来最大日降雨量209 mm, 计算日降雨径流量^[32]为8.36 m³(植物措施径流系数确定为0.4), 径流小区的洪峰流量采用经验公式^[33]计算为0.71 m³/h。设集流槽横断面为矩形, 底宽0.20 m, 深0.15 m, 长度与径流小区宽度一致, 槽底部设计为两侧向中间倾斜的形式, 坡度0.03, 便于地表径流及泥沙的集中, 顶部加设盖板。引水管选用管径100 mm PVC管, 分流孔用Φ50 mm分流管, 采用7孔分流法, 用钢管焊接在分流桶外壁, 选择中间一孔将分流水量引入集流桶中, 分流孔孔底高度0.6 m, 集流桶进水孔孔底高度0.6 m。分流箱和集流桶固定于径流小区下方阶梯平台, 分流箱采用1.2 mm铁板焊接而成, 外刷防锈漆, 直径0.8 m, 顶部加盖加锁, 桶内设观测标记。集流桶直径0.8 m, 桶内设观测标记。分流箱、集流桶底部连接钢管, 并安装阀门, 作为排水孔。

1.4 指标监测

研究监测了2015年2月至2016年5月期间11个标准径流小区共25组降雨量、地表径流量及泥沙量数据, 监测方法简要介绍如下^[31]：(1)降雨量由设置在径流小区附近的翻斗式数字雨量计观测记录, 实时监测日降雨量、降雨次数及降雨过程等指标；(2)地表径流量为当次降雨雨停后对分流桶水位进行3次重复测量取平均值计算所得；(3)泥沙量测量为雨后对各径流小区集流桶进行充分搅拌后, 采集水样500 mL, 用滤纸过滤后烘干测定。每次雨后测量后对分流箱和集流桶进行排水排沙清理, 确保不产生滞后影响。

1.5 数据处理与分析 1.5.1 数据处理

对试验的实测数据进行初步筛选, 去除离群值。径流量和泥沙量数据采用SPSS 22.0软件One-way ANOVA的LSD法进行差异显著检验。数据分析与绘图在Origin 2017完成。

根据前期文献调研, 将观测数据分为4类影响因子：降雨因子包括总降雨量、降雨天数、平均降雨量、最强雨量和临近降雨量; 土壤因子包括土壤厚度、腐殖质厚度, 土壤含水率和土壤容重; 植被因子包括林分密度、平均树高、胸径、针阔比、灌草层Patrick指数、盖度和生物量；生态修复措施因子包括补植阔叶树至针阔比3∶2、补植至7∶2、施肥和对照4种不同处理。

1.5.2 数据分析方法

地理探测器包括四个探测器, 分别为：(1)用于探测现象背后是否存在分异及各因子对其分异性的解释度的因子探测器；(2)识别和量化不同因素的相互影响及共同效果的交互作用探测器；(3)判断两个影响因素属性均值是否存在显著差别的风险区探测器；(4)比较两因素对现象的空间分布影响是否有显著差异的生态探测器^[25], 本研究使用的为因子探测器和交互作用探测器。采用2016年6月的土壤和植被因子实地调查数据进行地理探测器分析。

(1) 因子探测

因子探测器可量化自变量和因变量在空间分布上的一致性程度, 其原理为：如存在变量X影响变量Y的现象, 那么它们的空间分布趋于耦合, 其分布既可以是线性也可以是非线性。两变量空间分布的耦合程度用可以地理探测器q统计度量, 反映X对Y的决定力, 进而揭示现象背后可能存在的因果关系。

将研究对象划分为h=1, …, L个层, N和σ²分别为总体的样本单元数和方差, SSW和SST为层内方差之和和总方差, 则q值的表达式为：

(1)

① 若分层是对Y的划分, 那么q=0说明Y没有空间分异性；q=1说明Y存在完美的空间分异性, 即Y空间分异程度达到100%；

② 若分层是根据X对Y的划分, 那么q=0说明X对Y没有任何影响力；q =1说明X完全控制了Y的空间分布, 即X对Y影响力达到100%。

(2) 交互作用探测

交互作用探测器可识别不同影响因子间的两两交互作用, 即评估变量X₁和X₂共同作用时是否会增加或减弱对变量Y分布的影响力(或不同变量间不产生交互作用)。其原理为：首先计算两个变量X₁和X₂对Y的因子探测结果, 分别记为q(X₁)和q(X₂), 再计算它们共同作用后对Y分布的决定力, 记为q(X₁∩X₂), 最后比较q(X₁)、q(X₂)与q(X₁∩X₂)数值大小并对其作用关系进行分类(表 2)：

鉴于地理探测器方法在分析类型数据方面具备的优势, 故将降雨、土壤、植被和生态修复措施因子进行离散化处理^[34](降雨因子根据《降水强度等级划分标准》分级, 植被和土壤因子根据《GBT26424—2010森林资源规划设计调查技术规程》分级, 其余如灌草多样性和生物量根据个人先验知识等比分级), 以探测各因子与地表径流量、坡面产沙量在时空分布特征上的一致性, 以q值量化各因子对其分异性的解释力, P值检验其显著性, 并探究各因子交互作用对坡面水土保持能力的影响。

2 结果与分析 2.1 补植施肥措施对针叶纯林水土保持能力的影响

径流量受不同补植施肥措施影响显著, 且降雨量越大, 各措施的径流量差异越明显(图 1)。不同措施林地径流量的时间变化趋势与降雨基本相同, 降雨量是影响坡面产生径流的主要因子, 其影响力强于补植施肥措施。观测期内, 在针叶纯林(4240 L)内采取补植阔叶树措施, 可有效减少径流量, 补植施肥3∶2林地、补植3∶2林地和补植7∶2林地径流总量分别为4342 L、2937 L和2272 L, 径流变化幅度为4.8%、-30.04%和-46.01%。同为针阔混交比3∶2的混交林, 施肥林地较不施肥林地径流总量提高了47.83%。方差分析结果显示, 补植3∶2林地和补植7∶2林地径流总量显著小于其余两种林地(P＜0.05)。

补植阔叶树可有效减小针叶纯林土壤流失(图 2)。对比针叶纯林(1561 g), 补植施肥3∶2林地、补植3∶2林地和补植7∶2林地泥沙量为1044 g、767 g和521 g, 分别减少了52.76%、69.71%和76.30%, 显著小于针叶纯林(P＜0.05)。泥沙量变化呈季节性分布, 春夏季(3月—8月)土壤流失远高于秋冬季(9月-来年2月), 泥沙量均值分别为66 g和8 g。针叶纯林春夏季泥沙量显著高于其他措施林地, 而秋冬季各措施样地泥沙量差异较小, 补植7∶2林地泥沙量显著低于其他林地。

2.2 坡面径流量和泥沙量的因子探测结果

林地环境对坡面径流量单因子影响力强于生态修复措施, 影响力大小排序为：降雨因子、植被因子、土壤因子和补植施肥措施因子(图 3)。降雨因子对径流量的单因子影响力最高, q值大于0.5, 且结果显著(P＜0.01)。其余影响因子对径流量均有一定影响, 其中单因子影响力最大的依次是土壤容重、林分密度、灌草盖度、树高和针阔比, q值均在0.08左右(P＜0.05)。生态修复措施因子补植处理和施肥的单因子影响力相对较小, q值在0.03左右。

地表径流量对泥沙量单因子影响力最高, q值达0.84, 远高于降雨因子。降雨因子在土壤流失过程中起主导作用, 除降雨天数为0.14外, 其余q值均高于0.2。林分密度、土壤容重、灌草层盖度、土壤含水率、灌草层总生物量、多样性指数对泥沙量的解释度均在0.12左右, 结果显著, 说明它们是影响坡面土壤流失的重要因素。

2.3 坡面地表径流量和泥沙量的交互作用探测结果

利用地理探测器对19个影响水土支持能力的因子进行交互作用探测(图 4)发现, 各影响因子间主要表现为非线性增强和双因子增强作用, 占比分别超过53%和35%。非线性增强作用主要出现在降雨因子和其他类型因子的交互中, 交互后影响力远大于单因子影响力；双因子增强作用则出现在降雨因子间交互作用和非降雨因子间交互作情况下, 交互影响力较单因子略有提升。

林分密度、灌草层盖度、土壤容重、胸径、针阔比、树高、灌草生物量和土壤含水率等8个因子与降雨因子交互后对地表径流量分布的影响力大于0.9。林分密度、灌草层盖度和土壤容重与降雨总量交互作用后影响力为1, 比单因子作用增强40%, 表明该类强烈的因子间交互作用主导了径流量的产生。

降雨因子与林分密度、灌草层盖度、灌草层生物量、土壤含水率和容重交互后对泥沙量影响力大于0.95, 降雨总量与其他因子交互后q值大于0.98甚至达到1, 表现出对泥沙量分布的强烈影响力。单因子影响较小的因素如抚育施肥、补植处理、土壤深度和腐殖质深度, 可与其他因子产生双因子或非线性的增强作用, 如补植措施对泥沙量单因子影响力为0.05, 与降雨因子交互后影响力高达0.98, 与其余非降雨因子交互后也可以达到0.13。

3 讨论

在本研究中, 在针阔比例3∶2的混交林采取穴施法施复合肥500 g降低了林地水土保持能力, 与同类研究的结果^[35]不一致, 这可能与土壤贫瘠和穴施法扰动周边环境有关。通过对比2014年4月和2016年6月的每木检尺数据, 发现各林地两年间树高平均增加(1.08±0.33) m, 胸径增加(0.46±0.17) cm, 乔木树干细小, 冠幅稀疏, 生长缓慢, 属于生长受限的“老头树”, 由于土壤养分极度匮乏, 故施肥500 g难以产生效果。还有研究指出贫瘠林地施肥可能导致土壤结皮率的增加, 降低土壤渗透能力, 从而导致径流量和泥沙量增加^[36]。

亚热带红壤区坡面水土保持能力单因子影响力从强到弱排序为：降雨因子＞植被因子＞土壤因子＞补植施肥措施因子。林下径流和泥沙量的产生与降雨密切相关, 姜芃等^[37]使用稀土示踪技术发现降雨是导致南方地区坡面侵蚀的主要因素, 侵蚀量与降雨强度有关, 与本文结果相似。林分密度、灌草盖度、灌草层生物量可有效表征植被覆盖情况, 而植被覆盖率已被证实与坡面径流泥沙流失量高度线性相关^[38]。径流量对坡面土壤流失过程影响强于降雨, 这可能与研究区成土母岩为花岗岩有关, 该类土壤粒径大、砂粒含量高、黏性差、结皮率高, 在降雨后形成地表径流和泥沙的速度快, 导致径流对土壤的侵蚀冲刷更加剧烈, 有研究^[39]表明花岗岩红壤在降雨为45 mm强度下8分钟产流, 50分钟内坡面径流增长迅速, 在产流后泥沙浓度达25 g/L, 高于其他类型土壤, 与本文结果一致。施肥、土壤深度和腐殖质深度的单因子影响力较小, 可能是由于该3类数据差异化不明显导致, 且抚育施肥属于间接性措施, 无法对坡面水文过程产生直接影响。

降雨因子与土壤植被等因子的增强交互作用对径流泥沙的产生具有很强的影响力, 远远高于单一因子的影响。无论是从单因子影响还是从交互作用影响来说, 降雨因子都是导致坡面水土流失的最关键因素^[40], 而林分密度、灌草层盖度和土壤容重与降雨的非线性增强交互作用对亚热带红壤区针叶纯林坡面水土保持功能产生强烈影响, 土壤容重对土壤的入渗持水和抗侵蚀能力影响很大, 且受母岩、成土过程、微生物和周边气候等因素影响存在较强异质性, 李凤鸣等^[41]通过人工模拟降雨发现容重和土壤侵蚀量呈极显著负相关关系, 随着容重变大侵蚀量明显减少, 其相关性不受雨强和坡度影响。林分密度、灌草层盖度等因素容易进行人为控制, 在人工林抚育、森林生态修复和水土流失治理工作中可多加考虑。

4 结论

在亚热带红壤区针叶纯林补植阔树叶可有效提升坡面水土保持能力, 补植的针阔混交比例为7∶2总体效果优于3∶2。在侵蚀退化严重区域, 对“老头树”使用穴施法单木施复合肥500 g难以产生促进生长作用。

森林水土保持能力影响因子交互后主要呈现为非线性增强和双因子增强作用, 除降雨和地表径流量外, 林分密度、灌草层盖度和土壤容重单因子影响力高, 与降雨因子交互作用后产生对坡面径流量和泥沙量的绝对影响力, 在开展森林生态修复、细化亚热带红壤区坡面侵蚀模型参数时可更多关注此类因子。

研究存在一定局限性：一是枯落物层也是决定森林坡面水土保持功能的重要层次, 在研究中缺少对枯落物的实测数据, 未能探讨枯落物的影响力大小及其交互作用；二是地理探测器要求输入数据为类型变量, 研究根据行业标准和先验知识等方式对数值变量进行离散化, 有一定主观性。