生态系统服务是指生态系统与生态过程所形成及维持人类赖以生存的自然环境条件和效用^[1]。森林生态系统有明显的土壤保持作用, 林冠层能够拦截降水且减弱雨滴对土壤表层的直接冲击和侵蚀, 根系可以固持土壤, 阻滞土壤颗粒流失, 改良土壤结构^[2-12]。Borrelli等研究认为减少全球4.1%的森林面积会使全球增加52%的土壤侵蚀量^[3]。目前森林生态系统土壤保持功能的研究重点集中在区域土壤侵蚀量、侵蚀强度^[5-10]或土壤保持功能价值^[6-17]等方面, 鲜见对不同类型森林生态系统土壤保持功能差异的研究^[6]。中国生态系统研究网络(Chinese Ecosystem Research Network, CERN)的森林生态系统台站覆盖北半球热带、亚热带和温带主要的森林生态系统类型, 是研究不同类型森林生态系统土壤保持功能差异的天然实验场；此外, CERN台站经过长期封山育林, 人类活动干扰较少, 均为保存完好的地带性植被, 是相同水热条件下自然本底较好的生态系统, 估算其土壤保持功能可为相同类型生态系统土壤保持功能评估提供参照^[18], 在一定程度上解决土壤保持功能评估区域可比性差的问题。

森林生态系统结构的空间异质性会导致森林生态系统服务功能的空间异质性^[19]。根据生态系统结构-过程-功能关系, 森林生态系统类型和其植被覆盖度都密切影响着土壤保持功能的强弱^{[13, 20]}。对不同森林生态系统类型土壤保持功能而言, 混交林优于纯林^[21], 天然林优于人工林^[22-25]；植被覆盖度上升则能有效遏制土壤侵蚀的发生^{[11, 20]}。森林生态系统实际发挥的土壤保持功能还会因为地形^{[8, 12]}、气候^[25]等环境条件以及林龄^[26]的不同而产生差异。因此, 对森林生态系统服务功能的动态变化与其环境影响因素的研究已成为目前生态系统服务功能的研究热点^[16], 但仍缺乏分区分类型的森林生态系统环境因素如何影响其土壤保持功能的研究。

森林生态系统是重要的陆地生态系统, 森林覆盖面积占陆地面积的20.31%^[27]。CERN森林生态系统台站包括落叶针叶林、常绿针叶林、落叶阔叶林、常绿阔叶林、人工纯林等各类型森林生态系统, 具有一定的生态系统代表性, 适合开展不同类型森林生态系统的土壤保持功能研究。本研究选用CERN的10个典型森林生态系统为研究对象, 模拟分析了森林生态系统2005—2015年土壤保持功能的时空分异规律, 并探讨了其影响因素。

1 数据与方法 1.1 研究区概况

本研究选择CERN的10个典型森林生态系统的综合观测场样地, 涵盖热带(西双版纳)、亚热带(鼎湖山、鹤山、哀牢山、会同、神农架、贡嘎山)和温带(茂县、北京、长白山)。样地纬度范围为21.96°N—42.4°N, 海拔变化范围为77.3—3160 m, 样地生态系统类型、地理位置、土壤类型等详细信息见表 1。

1.2 数据来源

研究数据包括气象数据、土壤数据和植被数据。气象数据包括来自CERN的2005—2015年人工观测日降雨量数据(daily precipitation)、由中国陆地生态系统通量观测研究网络(ChinaFLUX)提供的BNF、DHF和CBF的2005—2010年半小时降雨量数据(30 minute precipitation)；土壤数据包括土壤机械组成数据(soil texture)和土壤有机碳数据(soil organic carbon, SOC), 均为CERN提供的长期地面监测数据, 监测层次为剖面0—100 cm；归一化植被指数数据(normalized difference vegetation index, NDVI)来自于美国国家航空航天局(NASA)网站(http://glovis.usgs.gov)下载的MOD13Q1数据产品, 空间分辨率为250 m, 时间分辨率为16 d；地表径流量数据(surface runoff)为CERN提供的各时段观测数据计算的年均值。

1.3 研究方法 1.3.1 RUSLE模型

本研究所选典型森林生态系统主要土壤侵蚀类型为水力侵蚀, 因此研究重点为水蚀情况。采用目前最为成熟且应用最为广泛的修正的通用土壤流失方程(revised universal soil loss equation, RUSLE)^[28-30]评估典型森林生态系统的土壤保持功能, 其计算公式为：

(1)

(2)

(3)

(4)

式中, A为单位面积土壤保持量(t hm^-2 a^-1), A_p为单位面积潜在土壤侵蚀量(t hm^-2 a^-1), A_r为单位面积实际土壤侵蚀量(t hm^-2 a^-1), Â为土壤保持率(%)。M为降雨侵蚀力因子；K为土壤可蚀性因子；LS为坡长坡度因子；C为地表覆盖因子；P为土壤保持措施因子。

(1) 降雨侵蚀力因子(M)

降雨侵蚀力因子是降雨引发土壤侵蚀的潜在能力^[30]。本研究采用日雨量侵蚀力模型计算降雨侵蚀力^[31-33]：

(5)

式中, M为月降雨侵蚀力(MJ mm hm^-2 h^-1 mon^-1)；D_i为第i天的侵蚀性日雨量(mm)(要求日雨量≥12 mm, 否则以0计算)；k表示天数；α、β是模型待定参数。

(2) 土壤可蚀性因子(K)

土壤可蚀性因子是评价土壤遭受降水侵蚀难易程度的重要指标, 与土壤机械组成和土壤有机碳含量密切相关^[28]。本研究K值计算公式^[30]为：

(6)

式中, K为土壤可蚀性因子(t hm² h MJ^-1 hm^-2 mm^-1)；SAN、SIL、CLA为砂粒(0.05—2 mm)、粉粒(0.002—0.05 mm)、粘粒(< 0.002 mm)的含量(%)；C为有机碳含量(%)。

(3) 坡长坡度因子(LS)

坡长坡度因子反映坡长和坡度对坡面土壤侵蚀的影响。LS值计算公式^[30]为：

(7)

式中, λ为坡长(m)；θ为坡度(°)；m为坡长指数；n为坡度指数。

通过分析国内研究成果^[34], 坡长指数m随坡度变化的取值范围为：

(8)

全国坡度指数n值主要集中于1.3—1.4之间, 本研究取1.35^[34]。因CERN样地大小为100m×100m, 因此本研究λ值为100 m。

(4) 植被覆盖因子(C)

植被覆盖因子反映了不同地面植被覆盖状况对土壤侵蚀的影响。本研究根据NDVI与植被盖度的经验关系计算月植被覆盖度F_c^[35]；利用蔡崇法等^[36]和江忠善等^[34]的方法计算植被覆盖度因子C, 公式为：

(9)

(10)

式中, F_c为植被覆盖度(%), e为自然对数底值。

(5) 水土保持因子(P)

水土保持因子是指在由一定水土保持措施的作用下, 水土流失面积与标准状况下土壤流失面积之比^[30], 其值介于0—1之间。本研究所选典型森林生态系统较少受到人为干扰, 因此本研究中所有站点P因子均设置为1^[37]。

1.3.2 降雨侵蚀力参数优化方法

参数敏感性分析有助于明确日雨量侵蚀力模型参数对降雨侵蚀力估算的影响程度, 从而更精确地进行参数优化。本研究采用OAT(one at a time)参数敏感性分析方法^[38]量化降雨侵蚀力对模型参数的敏感性, 将日雨量侵蚀力模型参数增加和减少10%时降雨侵蚀力的变化百分率均值作为对应参数的敏感性度量。

Kateikord等^[39]实验发现, 影响降雨侵蚀的两个重要参数是降雨持续时间和降雨强度, 降雨持续时间影响地表径流量产生时间；降雨强度影响降雨动能, 雨滴直径增大导致其动能增加, 使雨滴对土壤颗粒的降解力越大, 增加土壤侵蚀量。Wischmeier等^[30]也认为降雨侵蚀力与降雨强度关系密切, 并提出降雨侵蚀力应以次降雨动能(E)与30分钟最大雨强(I₃₀)的乘积来衡量。王万中等^[40]分析了中国降雨侵蚀力资料后认为EI₃₀是一个物理性质的函数, 可以反映雨滴对土壤的剥离和运输能力, 且兼顾我国大多数地区的降雨特性, 对降雨侵蚀力模拟效果好。但由于长时间序列、大范围的次降雨过程数据很难获得, 因此本研究以分布在热带(BNF)、亚热带(DHF)和温带(CBF)3个台站的半小时降雨量数据估算的降雨侵蚀力精确值(公式11—13)对日雨量侵蚀力模型进行参数优化。通过模型决定系数(R²)评价优化结果。

(11)

(12)

(13)

式中, M_i为降雨侵蚀力精确值(MJ mm hm^-2 h^-1 mon^-1)；e_r为降雨动能(MJ mm^-1 hm^-2 h^-1)；ir为雨强(mm/h)；E为日降雨总动能(MJ/mm)；P_r为侵蚀性日雨量(mm)；n为降雨天数。

2 结果与分析 2.1 降雨侵蚀力参数优化与验证结果

通过对日雨量侵蚀力模型的参数进行敏感性分析, 发现参数α和参数β敏感性相近, 分别为10%和12%；参数α变化取值范围在0.4—4.44之间, 参数β的变化取值仅在2.51—3.25之间, 与参数α相较, 参数β区域空间差异较小, 这与章文波的研究结果一致^[33]。因此本研究利用长白山、鼎湖山和西双版纳3个台站30 min观测降雨数据计算降雨侵蚀力为精确值对参数α进行优化, 用优化后的日雨量模型重新模拟了降雨侵蚀力, 并与优化前的降雨侵蚀力模拟效果进行比较(图 1), 发现未进行参数优化的日雨量侵蚀力模型对3个台站的降雨侵蚀力都存在高估, 其中对温带阔叶红松林(CBF)和亚热带季风常绿阔叶林(DHF)的高估十分明显。参数化后对温带和亚热带森林生态系统的模拟效果有显著提高, R²分别由0.5和0.47提高到0.61和0.69, 并且拟合线更接近于1:1线(图 1)；对热带季节雨林的降雨侵蚀力模拟效果提升较小, R²由0.7提升至0.71(图 1)。

2.2 森林生态系统土壤保持功能空间格局分析

根据参数优化后的RUSLE模型的计算结果, 通过土壤保持量和土壤保持率两个指征对各森林生态系统土壤保持功能进行评估。表 2显示了CERN2005—2015年各典型森林生态系统的土壤保持量、土壤保持率和土壤侵蚀量。2005—2015年期间土壤保持量呈现自北至南显著上升(R²=0.65^***)的空间格局(图 2)。土壤保持量变化范围为4.44—891.67 t hm^-2 a^-1, 其中鼎湖山、西双版纳和鹤山土壤保持量较高, 大于500 t hm^-2 a^-1, 其余样地的土壤保持量除长白山低于10 t hm^-2 a^-1, 均处于100—500 t hm^-2 a^-1范围内。实际土壤侵蚀量变化范围为0.07—11.51 t hm^-2 a^-1, 其中贡嘎山实际土壤侵蚀量高于10 t hm^-2 a^-1, 长白山、茂县、神农架、哀牢山和西双版纳实际土壤侵蚀量均低于1 t hm^-2 a^-1。根据《土壤侵蚀分类分级标准》(SL190—2007)^[41]的容许土壤流失量判断, 贡嘎山、鹤山和北京的实际土壤侵蚀量超过了容许土壤流失量, 属轻度侵蚀, 其余各站均未超出容许范围, 属微度侵蚀。

由于不同森林生态系统的土壤保持量受森林结构以及各自所处的气候条件、地形条件等环境因素的影响, 土壤保持量存在明显的量级差异。因此本研究采用土壤保持率反映森林生态系统对土壤侵蚀的调控速率。从表 2可以看出各森林生态系统的土壤保持率均达到97%以上。饶恩明提出2010年全国生态系统平均土壤保持率为95.71%^[42], 孙文义等指出2010年黄土高原草地和林地的土壤保持率分别为88%和97%^[6]。与之比较, CERN森林生态系统各样地土壤保持率均高于全国均值, 可达到97%以上, 说明各森林生态系统土壤保持功能都得到充分发挥。本研究利用OAT敏感性分析方法对RUSLE模型中的环境因子进行了敏感性分析, 选取敏感性>2%的4种环境因子并分析其与土壤保持功能的相关性(表 3), 各环境因子取值均为2005—2015年CERN森林生态系统样地观测数据多年均值。降雨量的空间格局对潜在土壤侵蚀量、土壤保持量的空间格局产生了显著影响(R²=0.52^*, R²=0.52^*), 这是由于水力侵蚀是研究区主要侵蚀类型, 降雨雨滴动能对土壤的打击作用破坏土壤的团粒结构, 随着降雨持续时间的增加产生了地表径流, 且降雨对径流的紊动作用又使径流分散, 对泥沙的搬运能力增强^[43-44], 所以降雨是导致土壤径流产生和土壤溅蚀的基础；土壤自身性质与土壤侵蚀营力的敏感性也存在密切关系, 实际土壤侵蚀量与土壤砂粒含量呈现显著正相关(R²=0.41^*), 本研究根据CERN提供的地表径流量数据发现土壤砂粒含量超过40%的GGF和DHF地表径流量最高, 均大于400mm, 所以考虑其实际土壤侵蚀量可能是由于土壤砂粒含量较高导致地表径流量增加而造成的, 但由于CERN各森林生态系统2010—2015年植被盖度均达到50%以上, 植被冠层拦截雨滴从而消解雨滴动能造成的土壤团粒结构的破坏^[2], 因此森林植被覆盖有效的减少了实际土壤侵蚀量, NDVI与实际土壤侵蚀量存在显著负相关(R²=0.64^**)。

2.3 森林生态系统土壤保持功能年际动态分析

2005—2015年期间, CERN各森林生态系统土壤保持率均呈现上升趋势(图 3), 茂县、鹤山和贡嘎山土壤保持率上升趋势显著, 上升速率分别为0.08%/a、0.11%/a和0.22%/a, 贡嘎山土壤保持率增长最快。其中茂县和鹤山林龄仅为30年左右, 贡嘎山为65年左右。进一步分析土壤保持率年际增长率与林龄的关系(图 4), 发现两者显著正相关, 呈现幂函数关系(y=9.37x^-1.4, R²=0.56^**)。因此中幼龄林的土壤保持率年际增长率较高, 而随着林龄的增加, 土壤保持率年际变化趋于稳定。

3 讨论

本研究基于RUSLE模型估算的CERN各典型森林生态系统土壤保持量与已有研究结果基本一致。裴夏根据2001—2006年CERN的样地长期监测数据估算的长白山阔叶红松林和西双版纳热带季节雨林的土壤保持量分别为9.68、607.64 t hm^-2 a^-1[45], 与本研究结果相近。王斌利用USLE模型计算CERN的9个森林生态系统2005年土壤保持量介于16.87—808.94 t hm^-2 a^-1[46], 本研究结果在长白山、茂县、哀牢山略低于该研究结果, 其余站点略高, 这可能是由研究时段差异导致的。与南方红壤丘陵区、东北地区等相似区域研究结果比较, 本研究结果基本位于其结果区间内^{[8, 44-47]}。

降雨是土壤水蚀的主要动力, 降雨侵蚀力是土壤保持功能评估中最基本也是最重要的因子^[31-33], 因此对降雨侵蚀力的准确估算能有效提高土壤保持功能评估的精度。但由于长期次降雨过程资料较难获取, 所以前人在研究中提出了以降雨量为基础资料的降雨侵蚀力估算方法^{[43, 48]}。本研究利用位于3个气候带的样地的30 min降雨量数据估算的降雨侵蚀力精确值对日雨量侵蚀力模型进行参数优化, 发现该模型在降雨丰富的地区模拟效果明显优于降雨较少的地区, 这与章文波等^[33]对日雨量侵蚀力模型的研究结果一致。王万中^[40]、钟莉娜等^[25]也提出, 在降雨较少的北方地区, 降雨多表现出历时短、雨量集中的特点。这时的土壤侵蚀基本上是与雨滴直径大小和速度相关的雨滴动能通过对土壤颗粒的剥离和地表径流的输移而产生的, 与降雨强度关系密切；而在降雨较丰富的南方地区, 暴雨特性差异则并不显著。而日雨量侵蚀力模型以雨量估算侵蚀力导致其在降雨较少的区域模拟结果偏差较大。因此, 应以不同区域的降雨特征为依据进行降雨侵蚀力估算, 才能提供更为准确的土壤保持功能评估。

虽然森林生态系统的土壤保持功能显著高于其他类型的生态系统^{[2, 22]}, 但本研究结果显示各森林生态系统台站的土壤保持量、潜在土壤侵蚀量和实际土壤侵蚀量的量级并不一致。这是因为土壤保持功能的空间格局是由引起土壤侵蚀的各类气候因素、地形因素, 以及保护土壤的植被共同作用而形成的。例如, 在西双版纳和鼎湖山, 尽管降雨丰沛导致潜在土壤侵蚀量较大, 但由于较高的植被覆盖度(>90%)对降雨形成阻隔, 实际土壤侵蚀量是较低的。而在贡嘎山, 由于坡度较大造成了实际土壤侵蚀量较大, 而且海拔较高(3160 m)导致水热条件差抑制了植被生长, 植被覆盖度为十站中最低仅有54%, 较低的植被覆盖度导致植被无法充分发挥其土壤保持功能, 这表明坡度会对实际土壤侵蚀量造成直接影响, 而海拔会通过影响植被生长进而造成实际土壤侵蚀量的变化, 这与陈龙^[12]、查良松^[11]、陈思旭^[8]等人研究结果相一致。植被覆盖度较低(57%)的北京和植被构成主要为人工纯林的鹤山实际土壤侵蚀量高于相似水热条件的长白山和鼎湖山, 张彪^[49]、郑华^[22]、赖仕嶂^[48]、康惠惠^[20]、杨吉华^[23]、黄承标^[24]等人提出植被覆盖度较低或者植被为人工造林时, 产生的实际土壤侵蚀量会大于植被覆盖度高或原生森林植被的实际土壤侵蚀量。

对于森林生态系统土壤保持功能的研究, 必须结合当地自然条件, 客观认识森林生态系统的功能特征, 但我国当前的森林土壤保持功能研究还集中在单一植被覆盖度指标^[50], 而忽视了其他植被特征对土壤保持功能的影响。本研究发现林龄为30a左右的茂县和鹤山, 植被覆盖度年际变化速率达到0.43%/a和0.88%/a且显著上升, 说明此时植被冠层可能处于逐渐郁闭, 植被尚处于空间结构的形成阶段, 对降雨动能的拦截作用随着植被冠层生长而显著增加, 因此其土壤保持率增长趋势显著。林龄大于100a的西双版纳、哀牢山和鼎湖山的植被覆盖度没有明显的年际变化, 其土壤保持率则没有表现出显著的增长趋势, 但当林龄较大时地面枯落物层累积可降低地表径流流速, 提高土壤下渗作用^[51], 其平均土壤保持率可达到99.4%以上, 本研究认为此时森林生态系统可能处于土壤保持功能得到充分发挥的功能显著阶段, 这与吴钦孝^[21]的研究结果相似。总体来说, 林木生长过程中对土壤以及森林生态系统环境的改变会影响森林生态系统土壤保持功能的发挥^{[21, 52-53]}。

本研究采用RUSLE模型分析了中国典型森林生态系统土壤保持功能时空变化特征及其影响因素, 为评估森林生态系统对土壤侵蚀的调控作用提供了依据。但CERN观测场样地选址多为平缓的坡中或山地下坡区域导致无法判断是否存在短距离沉积, 且由于RUSLE模型仅考虑水蚀动力的固定参数, 未解释森林生态系统结构特征与土壤保持功能的内在联系, 因此今后还需要结合长期和系统的观测实验, 获取长时间序列且精度更高的相关实测数据, 综合分析森林植被对土壤侵蚀的控制作用, 研究土壤保持功能各个因素间的相关性, 从而能更客观全面的评估森林生态系统的土壤保持功能。

4 结论

1) 日雨量侵蚀力模型在降雨较丰富的热带地区模拟效果比在亚热带和温带地区模拟效果好。对日雨量侵蚀力模型进行参数校正后热带地区模拟效果提升较小, 亚热带和温带地区的模拟效果有显著提升。

2) 2005—2015年期间, CERN各典型森林生态系统土壤保持量呈现北低南高的空间格局(R²=0.65^***)。本研究中贡嘎山、鹤山和北京为轻度侵蚀, 其余各站均属微度侵蚀。所有样地土壤保持率均达到97%以上, 说明各森林生态系统发挥了较好的土壤保持功能。

3) 降雨量对潜在土壤侵蚀量和土壤保持量的影响最大, 其决定系数0.52；NDVI与实际土壤侵蚀量呈显著负相关(R²=0.64^**), 土壤砂粒含量与实际土壤侵蚀量呈显著正相关(R²=0.41^*)；植被林龄与土壤保持率年际增长率呈幂函数关系, 说明林龄越小土壤保持率年际变化越快。