造林后表层土壤有机碳变化及影响因素分析

王 艺 杰，于 丽 君，张 稳，孙 文 娟*

(1.中国科学院植物研究所植被与环境变化国家重点实验室，北京 100093;2.中国科学院大学，北京 100049;3.中国科学院大气物理研究所大气边界层物理和大气化学国家重点实验室,北京100029)

0 引言

大气CO2浓度升高[1]是全球变化的重要问题，树木通过光合作用固碳是缓解全球CO2浓度上升的有效策略之一[2],因此，造林和再造林被认为是增加陆地碳固存、缓解大气CO2浓度上升的有效管理方式[2]。陆地生态系统的大部分碳储存在土壤中，土壤碳储量大于植被和大气碳储量之和[3]，其碳输入与输出平衡对于大气CO2浓度有重要影响[4],尤其是土壤有机碳比无机碳参与地球生物化学循环的活跃性更高[5]，因此，合理评估大范围造林后土壤有机碳密度变化，可为国家实施森林管理、规划减排计划提供可靠的数据支持，且有助于国际社会明确人工林土壤有机碳库在抵消碳排放中的作用。

由于缺乏可靠的造林后土壤有机碳密度估算方法，目前主要采用测定数据外推或基于土壤有机碳密度变化与环境要素关系的简单假设进行估算。例如：方精云等基于土壤固碳占整个人工林生态系统固碳量约1/3的比例，结合森林清查数据折算出1981-2000年中国土壤碳汇[6]；Nave等将时间作为唯一自变量估算美国造林后0～10 cm的土壤有机碳密度变化量[7]。这些简单估算忽略了影响造林后土壤变化的诸多因素及其相互作用的复杂性，导致结果不确定性较大。有研究指出，造林后土壤有机碳密度变化与造林年限存在典型的时间函数关系[8,9]。Nave等进一步提出，如果土壤有机碳密度与植被的结构、生产力以及气候和土壤的物理特性等因素存在定量响应关系，则可以通过植被生产力等数据估计土壤有机碳密度[10]。因此，本文通过收集全球造林后土壤有机碳变化的测定数据，分析影响人工林表层(0～20 cm)土壤有机碳密度变化的因素，并建立人工林表层土壤有机碳指标与植被、气候、土壤和土地利用因素关系的定量模型，以此估算我国不同区域造林后土壤有机碳变化。

1 数据及研究方法

1.1 数据来源与筛选

本文2009-2020年全球造林数据来源于Web of Science和中国知网等文献数据库，2009年前数据根据已发表的meta分析或综述文献[8,9,11-15]获得，对初步获取的340篇中英文文献进一步筛选：文献记述中包含造林前土地利用类型、土壤容重、造林树种、林龄信息，土壤采样深度大于20 cm，且必须包含对照组；若文献中采用时空代替法进行长时序研究，则将造林年限小于3年的土壤作为对照组。经过筛选，共有104篇文献满足条件。

样点数据在土壤、气候环境、树种和造林年限等方面的代表性，是建立区域适用的造林后土壤有机碳指标估算模型的必要前提。经上述筛选后的数据共包含123个研究样点(图1)，390组成对数据，造林年限跨度从6个月到91年，包含针叶树23种(油松、马尾松、杉木、赤松、黑松、辐射松、湿地松、云杉、冷杉等)、阔叶树31种(杨属、桉属、合欢属、栎属等)。样点的土壤质地包括粉壤土(粉粒含量<60%)、黏土(黏粒含量>60%)、砂土(砂粒含量>90%)等，涵盖从低纬度热带到中高纬度寒温带等主要陆地生态系统气候带，造林前土壤有机碳含量(有机碳重量百分比)范围为0.03%～6.15%，整理和计算的土壤数据均为0～20 cm深度。

注：审图号为GS (2016)1665号，底图无修改。

1.2 自变量选取与缺失值补充

依据已有研究[10-12,16]，本文先验性地选取植被类型、造林前土地利用类型、造林年限、乔木地上/地下有机碳密度、土壤初始有机碳密度、气候湿润度、土壤pH值、土壤黏粒+粉粒含量(clay+silt)作为自变量，对土壤有机碳变化进行分析。其中植被类型、造林前土地利用类型和造林年限信息在文献的样点描述中均有记述，个别文献中未明确记述的自变量通过以下公式进行填补。

ABC=SAB×TD×0.5/1 000

(1)

UBC=SUB×TD×0.5/1 000

(2)

式中：ABC、UBC分别为乔木地上、地下有机碳密度(Mg C hm-2)；SAB、SUB分别为单株乔木地上、地下生物量(kg/株)；TD为乔木密度 (株/hm2)，利用树木异速生长方程计算[17]；常数0.5为林木生物量含碳量[18]。

SOCDbf=BD×D×SOC

(3)

(4)

式中：SOCDbf为土壤初始有机碳密度(Mg C hm-2)；BD为土壤容重(g/cm3)，缺失值通过Adams方程计算[19]；D为土壤深度(cm)；SOC为土壤有机碳含量(%)。

IdM=P/(T+10)

(5)

式中：IdM为气候湿润度[20]；P为年降水量 (mm),T为年均气温(℃),二者缺失值采用AgMERRA样点区的1980-2010年气候平均值代替[21]。

本研究从GSDE(Global Soil Dataset for use in Earth System Models)数据库[22]获取样点坐标对应的土壤pH值及其机械组成数据，为使GSDE数据库中的土壤深度分层与本文关注的土壤深度一致，需对GSDE数据进行土壤深度转化，计算公式为：

(6)

(7)

Material_content0～xi=Material_massxi/Soil_massxi

(8)

式中：Soil_mass为单位面积的土壤质量；Material_mass为单位面积的物质质量；Material_content为物质百分含量(%)；0～xi表示从0 cm到xi深度的值，i为土壤层序数。

1.3 土壤有机碳指标计算

本文表征土壤有机碳变化的指标包括土壤有机碳密度占生态系统有机碳密度的比例(PSOCD)和土壤有机碳密度变化量(ΔSOCD)。因不同土地利用方式下土壤容重不同，相同深度下土壤质量也不同[23]，因此需要先进行土壤等质量法校正(Equivalent Soil Mass，ESM)，使造林前后的土壤在相同质量下进行有机碳密度对比[9,12,23]，计算公式为：

SOCafn×100

(9)

ΔSOCD=SOCDcorr-SOCDbf

(10)

PSOCD=SOCDcorr/(SOCDcorr+ABC+UBC)

(11)

式中：SOCDcorr表示校正后的土壤有机碳密度(Mg C hm-2)；SOCDaf表示造林后土壤有机碳密度(Mg C hm-2)，其计算同式(3)；i为土壤层序数；n为最底层土壤序数。

1.4 模型建立与退耕还林后土壤固碳情景模拟

研究全球造林后土壤有机碳指标变化的文章以meta分析为主[8,11]，但该方法只能研究分类变量，当变量是数值型(特别是连续型)数值时，需先划分为分类变量[8,11]，该过程会造成大量信息损失，且无法定量分析影响因素的总效应。本文根据前人研究结果，选取对土壤有机碳指标有显著或潜在影响的自变量[8,10-13]，采用逐步线性回归模型，选择AIC(Akaike′s Information Criterion)最小的方程形式，分析其对土壤碳指标的影响。鉴于回归模型易受数据极值影响，要求分析数据满足方差齐性和正态性等[24]，本文同时使用随机森林模型[25]分析土壤有机碳指标变化的影响因素。首先将单个决策树的特征空间最大值设为3(即自变量总数的1/3)，根据遍历调优的结果，将随机森林建立子树的数量设为1 100；然后使用交叉检验方法，以均方根误差(RMSE)(式(12))检验模拟值和观测值之间的偏差[26,27]，选择相应变量代入随机森林模型。

(12)

式中：Pi和Oi分别为样本i的模拟值和观测值；n为样本量。

退耕还林后土壤有机碳指标变化的情景模拟使用全国尺度的土壤理化数据、年均温和年降水数据，造林区域、造林密度和树种数据来源于中国《造林技术规程》《主要树种龄级与龄组划分》[28,29]。中国《造林技术规程》主要以年积温和年降水量为标准划分造林气候区[28]，得到9个造林气候区，除极干旱区不适宜乔木生长外，其余造林气候区分别为寒温带区、中温带区、暖温带区、亚热带区、热带区、半干旱区、干旱区和高寒区。

2 结果分析

2.1 土壤有机碳指标变化影响因素定量分析

逐步线性回归模型分别定量描述了PSOCD和ΔSOCD与植被、气候、土壤和土地利用因素的关系，如表1所示:对PSOCD有正向作用的自变量为土壤黏粒+粉粒含量(P<0.001)与土壤初始有机碳密度(P<0.001)，其回归系数分别为0.15和0.36；有负向作用的为土壤pH值(P<0.01)、乔木地上(P<0.001)和地下(P<0.05)有机碳密度、造林年限(P<0.01)，其系数分别为-1.70、-0.28、-0.33和-0.15；造林前土地利用类型对PSOCD的影响排序为农田>草地>荒地>林地(P<0.001)。对ΔSOCD有正向作用的是土壤黏粒+粉粒含量(P<0.001)、乔木地下有机碳密度(P<0.01)、造林年限(P<0.01)和气候湿润度(P<0.05)，其回归系数分别为0.26、0.36、0.16和0.07；有负向作用的为土壤初始有机碳密度(P<0.001)、土壤pH值(P<0.1)和乔木地上有机碳密度(P<0.05)，其系数分别为-0.15、-0.87和-0.07；造林前土地利用类型对ΔSOCD的影响排序为农田>荒地>草地>林地(P<0.05)，植被类型对其影响排序为针阔混交林>常绿阔叶林>常绿针叶林>落叶阔叶林>落叶针叶林(P<0.05)。

在随机森林模型中，自变量与因变量的关系没有解析表达式，仅能通过偏相关曲线进行判断，因此将随机森林模型的偏相关曲线与逐步线性回归各自变量对因变量的作用方向相结合(图2、图3)，便于理解随机森林模型的解析效果。由图2可知，随机森林模型与逐步线性回归模型结果基本一致，即PSOCD随土壤初始有机碳密度和土壤黏粒+粉粒含量增加而增加，随乔木地上/地下有机碳密度、土壤pH值和造林年限增加而减少。其中，土壤黏粒+粉粒含量高有利于土壤形成团聚体，吸附并保护碳，因此有利于土壤有机碳密度的增加[11,12]。土壤pH值影响树木生长和土壤碳过程，低pH值会降低树木的生长速度，从而减少土壤的碳输入，同时会降低土壤有机质的分解速率，有利于土壤有机碳积累[30]；反之，高pH值地区树木生长速度更快，土壤有机碳分解也更快，综合结果是碱性土壤的PSOCD随着pH值增加迅速降低(图2c)。PSOCD随造林年限增加而下降，主要是由于土壤有机碳密度的增加速率小于乔木有机碳密度的增加速率。已有研究通常将气候带作为分类变量分析温度和降水对土壤有机碳的影响[12,15]，认为从温带地区到亚热带地区，土壤有机碳密度有增大趋势[11,12,19]，这与本研究结果一致(图2d)。

注：正负号表示变量在逐步线性回归模型中的符号方向，未出现则表示变量在逐步线性回归中被剔除， *、**、***分别表示P<0.05、P<0.01、P<0.001，下同。

图3 不同自变量与 ΔSOCD的随机森林模型偏相关曲线及其在逐步线性回归模型中的作用方向Fig.3 Partial correlation curves of ΔSOCD and independent variables based on random forest model and the action direction of independent variables in stepwise linear regression model

由图3可知，除土壤pH值、气候湿润度和乔木地上有机碳密度外，其他自变量在随机森林模型中的偏相关趋势与在逐步线性回归中的作用方向一致，说明ΔSOCD随土壤黏粒+粉粒含量(小于80%时)、乔木地下有机碳密度和造林年限(小于60年时)增加而增加，随土壤初始有机碳密度增加而减少，而与土壤pH值的关系更接近二次曲线(图3c)。Rowley等指出酸性土壤通过Fe3+和Al3+的吸附作用保护土壤有机碳，碱性土壤主要通过Ca2+保护土壤有机碳，只有中性土壤是土壤微生物利用土壤有机碳的“机会窗口”[31]，该结论与本研究曲线形状表现的意义一致且偏相关曲线的拐点也与“机会窗口”的pH值范围基本一致。ΔSOCD与气候湿润度的关系也表现为二次曲线，当温度与降水协同更有利于土壤呼吸增加时，土壤有机碳密度也随之降低，表现为ΔSOCD处于低值，从“三基点”的角度解释，达到最大土壤呼吸的气候湿润度小于达到最大土壤有机碳输入的气候湿润度，ΔSOCD在该点之前随着气候湿润度增加而减少，之后随着气候湿润度增加而增加，因此，ΔSOCD与气候湿润度的曲线趋势与文献[15]中气候带对土壤有机碳密度变化速率的影响趋势相似。

乔木地上有机碳密度的作用在两个模型中存在差异(图3e)，乔木地上、地下有机碳密度在逐步线性回归中的非标准化系数分别为-0.07和0.36(表1)。实际情况下乔木地上、地下生长同步，因此，从两者系数的比例看，只有乔木地下与地上有机碳密度比例小于1∶5.14的树种，土壤有机碳密度才会随着乔木有机碳密度增加而减少，可近似认为根冠比小于0.19的树种造林后土壤有机碳密度会下降。由于针叶树和阔叶树的生物量分配策略不同，阔叶树的根系发达[32]，使得地下生物量更高的阔叶树提高了根系对土壤有机碳的输入，故树种的根冠比越大，越有利于土壤有机碳积累。有研究显示，根据文献数据计算出华北落叶松的根冠比是0.26，而野外采样获得的华北落叶松的根冠比是0.17[33]，不同研究中由于样点针叶树种的根冠比存在差异，导致土壤有机碳密度变化的计算结果出现土壤有机碳损失[13,34]、累积或者保持不变[35]。

在随机森林模型中，通常根据IncMSE指标(去除自变量后随机森林预测误差增加的百分数)明确自变量对因变量的影响程度，数值越大说明自变量越重要。基于此，对逐步线性回归模型和随机森林模型的自变量进行重要性排序(图4)。由图4a可知，随机森林模型中，PSOCD自变量的重要性从大到小为土壤初始有机碳密度(76.59%)、乔木地上有机碳密度(73.73%)、乔木地下有机碳密度(42.97%)、气候湿润度(33.22%)、土壤pH值(28.39%)、土壤黏粒+粉粒含量(26.98%)、造林年限(26.12%)、植被类型(18.33%)和土地利用类型(17.01%)，而逐步线性回归模型中自变量排序与之略有不同，但主导因素一致，均为土壤初始有机碳密度、乔木地上有机碳密度、乔木地下有机碳密度，而土壤pH值、土壤黏粒+粉粒含量和造林年限是次要因素。由图4b可知，ΔSOCD的自变量在随机森林模型中的重要性从大到小为土壤黏粒+粉粒含量(52.86%)、土壤初始有机碳密度(52.71%)、造林前土地利用类型(49.94%)、土壤pH值(41.89%)、气候湿润度(36.79%)、植被类型(28.24%)、乔木地下有机碳密度(21.88%)、乔木地上有机碳密度(21.75%)和造林年限(21.39%)。其中，乔木地上、地下有机碳密度和造林年限的回归方程标准化系数排序靠前，但重要性排序靠后；而气候湿润度和土壤pH值的重要性排序靠前，但标准化系数排序靠后。两种模型结果出现差异的原因可能是统计理论方法不同，逐步线性回归模型是对全部数据进行方差分解，而随机森林模型采用有放回抽样形成不同子集，分析每个自变量对因变量预测结果的影响。但两模型中土壤黏粒+粉粒含量和土壤初始有机碳密度均位居前列，而逐步线性回归方程中造林前土地利用类型也达到显著水平。综合考虑可认为主导ΔSOCD的因素为土壤黏粒+粉粒含量、土壤初始有机碳密度和造林前土地利用类型，其他因素可归为次要因素。

注：纵轴旁数字表示逐步线性回归中标准化系数绝对值的排序。

2.2 退耕还林后土壤有机碳指标变化情景模拟

本文将10折交叉检验的均方根误差(RMSE)作为自变量筛选依据，RMSE越小说明预测值与观测值越接近。重复5次检验后得到土壤有机碳指标的RMSE与自变量数量的关系(图5)，可以看出，PSOCD、ΔSOCD的RMSE在自变量数量大于4后，基本不变。结合图4的分析结果，选择图4中IncMSE排序前四的自变量(PSOCD的模拟变量为乔木地上/地下有机碳密度、土壤初始有机碳密度和气候湿润度，ΔSOCD的模拟变量为土壤黏粒+粉粒含量、土壤pH值、土壤初始有机碳密度和造林前土地利用类型)作为随机森林模型的驱动变量，进而模拟我国不同区域退耕还林后土壤有机碳时空变化(表2)。

图5 基于随机森林模型的0～20 cm土壤有机碳指标交叉检验Fig.5 Cross-validation of SOC0～20 estimation based on random forest model

表2 基于随机森林模型估算我国退耕还林后0～20 cm PSOCD 和ΔSOCD变化Table 2 Estimation of PSOCD and ΔSOCD changes in 0～20 cm soil layer under the "Grain-for-Green" Program based on random forest model

由表2可知，我国各气候区针叶树的PSOCD区域均值从造林5年后的80.15%～92.87%下降到造林40年后的27.35%～61.09%，5～40年间降幅从寒温带(52.71%)到热带(61.61%)逐渐增大，从干燥气候区(干旱和半干旱区，29.57%～35.37%)到湿润气候区(52.71%～61.61%)降幅也逐渐增大。阔叶树PSOCD均值从造林5年后的61.21%～85.29%下降到造林40年后的25.67%～53.94%，期间降幅从寒温带(31.40%)到亚热带(41.92%)逐渐增大，从干燥气候区(31.35%～33.73%)到热带和亚热带气候区(34.21%～41.92%)降幅也逐渐增大。在内陆(干旱、半干旱和高寒)地区，阔叶树造林40年后PSOCD均值小于针叶树，而其他气候区阔叶树造林40年后PSOCD均值大于针叶树(表2)，造林40年后暖湿地区的PSOCD均值小于冷干地区。退耕还林后ΔSOCD随造林时间有不同程度增加，在内陆地区针叶树土壤有机碳密度的增量和增幅均小于阔叶树，而其他气候区针叶树土壤有机碳密度的增量和增幅均大于阔叶树。造林后5～40年针叶树土壤有机碳密度增幅最大的是中温带(5.43 Mg C hm-2)，增幅最小的是干旱区(0.23 Mg C hm-2)；阔叶树土壤有机碳密度增幅最大的是高寒区(8.87 Mg C hm-2)，增幅最小的是半干旱区(1.40 Mg C hm-2)，高寒气候区针叶树和阔叶树对土壤有机碳密度增量(造林40年后)和增幅(造林后5～40年间)的影响差异很大。

Deng等研究我国退耕还林土壤有机碳密度变化时发现，在针叶树造林后31～40年内，0～20 cm土壤固碳速率大于阔叶树[36]，人工针叶林0～20 cm土壤的平均固碳速率为0.03 Mg C hm-2a-1，人工阔叶林为0.17 Mg C hm-2a-1[36]。与之相比，本研究预测我国造林40年后针叶林和阔叶林0～20 cm土壤平均固碳速率分别为0.21 Mg C hm-2a-1和0.22 Mg C hm-2a-1(气候区内样本量加权平均)，模型预测的退耕还林土壤固碳作用较高，尤其是针叶林土壤固碳速率。但陈越的研究显示，青海、东北和新疆地区针叶树造林后0～20 cm土壤固碳速率为0.58 Mg C hm-2a-1，大于相同区域中阔叶树造林后土壤平均固碳速率(0.30 Mg C hm-2a-1)[27]，与本研究模拟的造林后5～40年间针叶林ΔSOCD增幅大于阔叶林ΔSOCD增幅结果一致。因此考察针叶树造林后的ΔSOCD时，需考虑更长时间尺度的土壤有机碳积累。

3 结论与展望

本文基于收集的全球造林后土壤有机碳指标变化文献数据，通过逐步线性回归模型和随机森林模型分析影响土壤有机碳变化的因素，并利用统计分析和随机森林模型模拟中国造林后土壤有机碳指标的变化。结论如下：1)土壤有机碳密度占生态系统有机碳密度的比例(PSOCD)主要受土壤初始有机碳密度和乔木地上、地下有机碳密度影响；土壤有机碳密度变化量(ΔSOCD)主要受土壤质地和土壤初始有机碳密度影响。土壤黏粒+粉粒含量越高，越有利于土壤固碳；但土壤初始有机碳密度越高的地区，造林后更易出现土壤有机碳密度下降，因此在农田或荒地上造林更有利于表层土壤有机碳密度增加。土壤有机碳输入随着乔木生长而增加，树种根冠比越大，越有利于土壤有机碳密度增加。但对整个生态系统而言，比土壤有机碳密度增加速度更快的森林碳增加速度使土壤固碳对生态系统固碳的贡献率随时间下降，因此考察人工林生态系统碳汇效应时，需考虑长时间尺度的土壤有机碳积累。2)随机森林模型模拟结果说明，我国造林40年后ΔSOCD自北向南随年均温增大而增大，从内陆向沿海随年降水量增大而增大，但造林后5～40年间，寒温带地区ΔSOCD的增幅更大，说明暖湿地区造林后植被生长速度快，土壤固碳效果更好，而寒冷地区的固碳潜力更大；造林40年后PSOCD自北向南、自内陆向沿海下降且幅度逐渐增大。研究结果可为未来基于大量数据利用随机森林模型预测区域或国家尺度造林后土壤有机碳指标变化提供支持，同时为在管理措施和造林技术方面保护土壤碳库提供指导。