基于DNDC的茶园土壤碳素空间分异及其影响因素*——以福建省安溪县为例

税 伟，吴朝玮，付 银

（福州大学环境与安全工程学院，福建福州 350108）

0 引言

茶产业是福建省泉州市安溪县的特色产业，茶叶种植产业化、专业化、集群化是安溪城市化发展的重要驱动力。城市化的快速推进也为安溪县茶产业链形成、品牌打造、茶叶国际销售网络渠道开拓提供了资金技术支持。近年来，地方政府号召积极开展低碳减排行动，因此减少茶园温室气体排放、建设生态茶园成为地方推动城市化和农业绿色发展的重要途径。目前，全球碳库的主体是陆地土壤碳库[1,2]，而城市建设的扩张蔓延带来的土地城镇化过程已经改变了原有土地利用的类型、结构、功能[3]及相应的土壤碳储量时空分布格局。对于以茶叶为经济作物的安溪县，依托低碳茶叶经济技术，建设立体人工复合生态茶园可有效提高土壤有机碳储量和增加茶园碳汇量[4]，这对于保持安溪全县茶叶产量稳定和茶产业链转型升级发展意义重大。另一方面，土地城镇化进程也改变了土壤的理化性质，影响土壤碳释放的强度，表现为CO2排放显著增加[5,6]，而CO2浓度的增加将进一步加剧温室效应，产生一系列危害人类健康与生态的连锁反应[7]。因此，在城市化和推动茶产业绿色专业化发展的背景下，探究安溪县茶园土壤碳素循环过程，明确茶园土壤碳储量状况及其空间分布格局，对于采取区域差异化技术模式建设生态茶园和增强茶园土壤碳汇能力具有重要意义。但以往茶园土壤碳素动态研究多以传统野外采样实测为主，时间、人力等消耗大，且由于受到海拔、坡度、植被覆盖类型、土壤属性以及施肥管理等多种因素的综合作用影响[8]，茶园土壤碳储量和温室气体排放量在空间与方向上表现出较大差异。于是进一步揭示各碳素指标的空间分异的影响因素，将对合理配置茶园土壤营养元素、改善安溪茶园施肥管理措施及最终实现茶园绿色生态发展提供指引。

目前，由李长生等[9]提出的反硝化分解（denitrification decomposition,DNDC）模型已经在模拟农业、草地和森林生态系统的土壤碳氮动态变化及收支平衡状况的研究中表现出良好的适用性[10]。在亚洲包括中国在内的多个国家和地区已将DNDC模型用于水稻[11]、玉米[12]、稻田—小麦轮作[13]等作物的碳氮动态循环研究，其能够契合中国区域的土壤、气候、作物生长状况，研究也已证实该模型在温室气体排放等方面的模拟精度较好。但目前基于DNDC模型的土壤碳氮模拟研究对象以农田和流域居多，关注的作物多为水稻、玉米等粮食作物，缺乏对经济作物特别是茶树的相关研究，对产业化、专业化、集群化发展的安溪茶区这一特定对象的土壤碳素的模拟与分析略显不足。而土壤初始属性、地形等多种因子对专业化茶区土壤各碳素指标的综合作用影响研究也较为缺乏。

因此，文章基于DNDC模型开展安溪县专业化茶区土壤碳素循环模拟，结合GIS空间统计方法（空间自相关分析和冷热点分析）揭示各碳素指标的空间分异特征。在对比分析普通最小二乘法（Ordinary Least Squares Method,OLS）和地理加权回归模型（Geographical Weighted Regression,GWR）等地统计回归模型的基础[14]上，识别各解释因子对目标变量的影响和作用效果差异。研究对专业化茶区的生态建设具有重要意义，以期为优化安溪专业化茶区茶叶种植结构和减少茶园温室气体排放等提供科学决策依据。

1 研究区概况与方法

1.1 研究区概况

安溪地处福建东南，经纬度范围117°36′E～118°17′E，24°50′N～25°26′N。安溪划归泉州市管辖，与南安、永春等市县相邻；县域内地貌以低山、中低山和丘陵为主，平原地区面积占比少；安溪县西北部是中低山区域，丘陵在安溪县境内广泛分布，平原主要分布在安溪东部的凤城镇、城厢镇等乡镇[15]。安溪县属南、中亚热带海洋性季风气候，地形地貌的差异使得东部外安溪区域形成南亚热带气候，该区年均温19～21℃，年降水量1 600 mm；西部内安溪区域由于受北方气流影响较大，坡谷地形成复杂气候状况，年均温16～18℃，年降水量1 800 mm，四季分明，因此该区域也是乌龙茶优良产区[16]。

安溪县位居中国重点产茶县第一位，全县茶树品种繁多，其中国家级茶树良种有6个，茶叶品质、产量均位于全国前列。整体上，安溪县茶叶种植专业化水平较高。税伟等[14]基于综合区位商法对安溪县茶叶种植专业化水平进行区域划分，发现其具有“西高东低”的空间格局，其中位于县城中心的凤城镇和城厢镇专业化水平较低，湖头镇、魁斗镇、尚卿乡、蓬莱镇、官桥镇等乡镇专业化水平处于中等地位，祥华乡、感德镇等安溪县西部、西北部乡镇专业化水平较高。

1.2 数据获取与预处理

该研究采用的模型是DNDC 9.5版本，在区域尺度上构建DNDC模型所需的数据详细如下。土壤类型数据来源于资源环境数据云平台，研究区安溪县域内主要有红壤、黄红壤、黄色赤红壤等9种类型土壤。土壤属性数据来源于中国土壤数据库的实测数据，按照类型进行查询匹配。气象数据来源于中国气象数据共享网，将2017年安溪县气象观测站的历史观测数据通过ArcGIS插值计算得到。作物参数数据来源于DNDC 9.5版软件提供的一套茶叶种植的基本生理和物候学参数，通过对相应参数进行率定修正以适应安溪茶园茶树种植模拟。农田管理数据主要包括施肥量、施肥日期以及作物播种收割日期等。茶树的种植日期参考黄东方[17]对安溪县铁观音的研究，铁观音在10—11月或早春2—3月定植，一般每年除草4～5次。施肥量主要包括无机肥和有机肥的数量和种类，其中无机肥核算是以安溪县农业局提供的各行政乡镇的无机肥输入量为依据；有机肥数据以安溪县统计局提供的各区域同年的猪、牛、羊和禽类等年末存栏头数为基准，折算系数参考蔡明等[18]研究结果，依据各作物所需比例，将无机肥和有机肥施用量分别按照公式（1）折算出安溪县茶园有机肥的施用量。各施肥日期及施肥比例参照安溪县农业农村局提供的测土配方施肥数据，按照春茶、夏茶与暑茶采摘前2∶2∶3比例来分配。

式（1）中，Fi是第i种作物单位面积施肥量（kg/hm2），Fsum是研究区施肥总量（kg），Sj是研究区第j作物的播种面积（hm2），Pj是第j种作物专家推荐的施肥量（kg/hm2），n是研究区作物种类。

1.3 DNDC模型模拟与验证

根据DNDC模型的模拟数据格式要求，参考田展等[19]的研究将整个安溪县域划分为多个格网单元，并假设每个格网属性是均一的。在ArcGIS10.2平台上将气候、土壤和地形数据与茶园面矢量文件进行空间连接，按照一致性与相似性原则，分割茶园矢量面数据，最终获得581个有效模拟单元。每个模拟单元分别以初始土壤属性的最大、最小值进行模拟，并取平均值代表该单元的值。该研究主要记录并分析茶园10～20cm土壤有机碳（SOC10～20cm）含量、凋落物输入碳量、土壤二氧化碳（CO2）排放量、土壤甲烷（CH4）排放量的变化规律。

该研究使用安溪县农业农村局提供的测土配方施肥数据中的土壤有机质含量（SOM）指标对DNDC模型模拟的结果进行精度验证，依据Guo和Gifford提出的转化公式[20]，将模型模拟的SOC含量数据转换为SOM数据进行对比验证。同时，模型的模拟值与实测值的偏差使用均方根误差（RMSE）进行表征，并进一步采用归一化均方根误差（NRMSE）对模拟值和实测值的拟合程度进行描述[21]，NRMSE<25%说明两者一致性程度较高。采用线性拟合度系数R2表征拟合效果，其值越大，模型模拟效果越好。

式（2）至（4）中，Po、Qm分别代表模拟值和实测值，N为验证数据的总个数，代表实测值均值，和为线性拟合的模拟值和模拟值的均值。

1.4 空间统计分析

该研究使用ArcGIS10.2平台进行空间统计分析。

空间自相关分析：用于描述各观测单元的空间邻近程度，该研究采用Moran′s I指数进行测度，其值范围在-1～1，小于0表示负向自相关，大于0值表示正向自相关[22]。

空间冷热点分析：核算观测单元的Getis-OrdGi统计量，从而辨识观测单元明显不同的冷热点。其原理是将观测单元及其邻近单元之和与区域单元总和进行对比[23]，计算统计检验值Zscore。以Zscore的正负、大小来区分观测单元属性在空间上聚类的情况。

回归分析：在各目标变量与解释变量之间进行关系建模分析。一般采用OLS进行最优线性拟合分析，但OLS是在全局意义上估算参数，忽视了局部的空间变化。GWR模型弥补了OLS模型缺陷，通过将地理事物的空间位置差异反映为拟合方程的回归系数随空间位置变化[24]，考虑空间局部非平稳性特征，可解析各解释变量对目标变量空间分布的影响程度。使用AICc准则和R2adj对GWR模型的拟合优度进行检验。

2 结果与讨论

2.1 模型参数率定与验证

DNDC模型默认的作物参数均为北美地区观测的结果，因此需要对模型中的内部参数或过程进行调整。乔帅帅等[25]的研究指出最高生物量、积温、需水量等对DNDC模型模拟结果影响较大。参考《中国茶叶大词典》[26]获知，茶叶生产1 kg干物质耗水量为385 kg；茶叶的经济产量为每公顷可产干茶1 500～3 000 kg，依据国际上通用的植物碳转化系数0.4[27]，核算茶叶的最大籽粒产量600～1 200 kg/hm2；通过历史日气温数据以及参考杨浩湘[28]对安溪铁观音生长积温的研究，安溪县茶叶生长期积温核算结果约为5 500℃；已有研究表明茶叶生长的最适宜气温为25℃[29,30]。通过参数调整以及将模拟值与实测值对比后进行的相应调整，模型模拟效果表现良好。该研究模型运行时间为2017年1月1日至12月31日。

将DNDC模型模拟的SOC含量数据与测土配方施肥数据中的土壤有机质数据进行对比，从而评价模型的模拟精度和安溪县茶园土壤质量状况。由于在DNDC的模拟格网单元内包含多个测土配方施肥样点数据，因此取在格网单元内部的样本点的SOM均值进行验证。经过核算，测土配方施肥数据中SOM的区间值为2.41～43.40 g/kg，均值为20.63 g/kg；模型模拟SOM的区间值为5.33～50.65 g/kg，均值为22.70 g/kg。进一步，模型模拟的SOM值与测土配方施肥数据的SOM实测值间RMSE为2.99，NRMSE为14.50%，线性方程y=0.943 7x+3.135 2，R2为0.88，上述数据表明DNDC模型模拟效果较好。

参照安溪县测土配方施肥数据的土壤养分丰缺指标分级标准，发现安溪县域内茶园SOM含量较为丰腴，说明其土壤肥力较好。该结论与杨如兴[31]等通过实验测定安溪县土壤有机质得出的安溪县土壤有机质含量丰富的结论较为一致。一般来说，优质茶园SOM在20 g/kg以上[32]，安溪县茶园SOM均值大于20 g/kg，这与安溪县是著名的专业化茶叶产区有关。与其他专业化茶区相比，安溪县茶园SOM含量低于四川蒙顶山茶园[33]SOM含量，其均值为27.22 g/kg，但高于河南信阳毛尖茶区[34]SOM含量，其均值为16.16 g/kg。

表1 模型率定参数

2.2 土壤碳素转化与循环结果

基于DNDC模型模拟茶园土壤碳素转化与循环结果表明（表2），各碳素输入指标中土壤有机碳SOC含量是茶园土壤最大的碳源，均值为11 007.90kg/hm2，其次为凋落物输入碳量，但凋落物输入碳量相对较少，均值为505.52 kg/hm2。各碳素输出指标中，土壤CO2排放量占比最高，均值为883.78 kg/hm2，CH4排放最低，均值为111.68 kg/hm2。这表明土壤通过呼吸作用产生的CO2排放是茶园土壤有机碳的主要损失途径，这一结论与王峰[35]等的研究结论一致。与前人研究相比，该研究模拟的茶园CH4排放量低于水稻田[36]、稻—稻—绿肥、稻—稻—油[37]轮作模式，这与茶园有机肥施用量相对较少有关。该研究模拟的茶园土壤有机碳含量年度变化（dSOC）为负值，表明安溪茶园土壤肥力呈现下降的趋势，可能与茶树枯枝落叶、修剪物以及施肥输入减少密切相关。

表2 基于DNDC区域尺度的各碳素指标模拟值kg/hm2

2.3 各模拟碳素指标空间分异

表3列出了空间自相关分析的结果。各碳素指标的Moran′s I均大于0，表明安溪县土壤各碳素指标具有明显的空间自相关性。各碳素指标均表现为正相关的聚集分布，且符合地理学第一定理的规律，其中SOC10～20cm的正相关性最强，其后依次为凋落物输入碳量、CH4排放量、CO2排放量。

表3 各碳素指标描述性统计

各碳素指标的空间分布结果表明，安溪县茶园土壤碳素储量丰富，各碳素空间上主要分布在安溪的长坑乡、感德镇、祥华乡等中低山和低山地区。首先，安溪县SOC10～20cm空间分布呈现较大差异，SOC10～20cm主要分布在感德镇、祥华乡、长坑乡、大坪乡等中低山和低山地区（图1-1）。实地调研发现这些区域茶叶种植专业化水平较高，多采用茶园集中管理方式开展科学施肥，土壤肥力较好，且有机肥施肥较多，因此SOC10～20cm含量较高。实地调研还发现剑斗镇、白濑乡、金谷镇等茶叶种植专业化水平也较高，但该区域茶园坡度较低，且位于晋江流域上游区域，土壤水分含量较高，土壤透气性相对较差，因此土壤微生物活性较低，导致有机质降解缓慢[38]；另一方面，该区域受水流作用的冲刷作用较强，造成土壤有机质易流失，因此SOC含量较低。总体上，安溪县SOC10～20cm空间分异的重要特征是专业化茶区较非专业化茶区更高。其次，凋落物输入碳的空间分布与SOC10～20cm较为一致，主要位于安溪县西部和南部低山和中低山地区（图1-2），茶叶等凋落物或茶树修剪残渣作为肥料返园，经微生物作用可提高土壤有机碳含量，增强土壤肥力，改善土壤结构，减少土壤侵蚀[39]。CO2和CH4排放量的高值区与凋落物输入碳高值区的分布较为一致（图1-3，图1-4），其原因在于茶叶凋落物经过腐殖质化和微生物分解等过程最终产生CO2，但在厌氧环境中会释放出CH4[40]。

图1-4 甲烷排放量

图1-3 二氧化碳排放量

图1-1 10～20cm土壤有机碳量

图1-2 凋落物输入碳量

各碳素指标的空间冷热点分析结果表明，各碳素的冷热点分布存在明显的空间异质性。SOC10～20cm热点区域主要分布在安溪县西部（图2-1），包括大坪乡、感德镇北部、桃舟乡西部、福田乡等低山和中低山区域；冷点区域主要位于安溪县东部和南部，具体为长坑乡南部、蓝田乡东部、祥华乡南部、虎邱镇北部和西南部等区域。总体上SOC10～20cm冷热点分布格局呈现“东冷西热”的特征。凋落物输入碳热点区域分布在长坑乡、祥华乡、大坪乡等低山和中低山区域，安溪县东南部地区有少量分布；冷点区域主要分布在县域东部和南部的丘陵和低山地区（图2-2）。CO2冷热点区域分布格局与SOC10～20cm总体上较为相似，但是祥华乡和长坑乡是CO2排放较为突出的热点区域；剑斗镇西部与感德乡交界处是明显的冷点区（图2-3）。CH4排放量热点区域主要集中在感德乡、大坪乡、祥华乡东部和长坑乡西部等中低山区域，剑斗镇、蓬莱镇、白濑乡、湖头镇、湖上乡、金谷镇等是CH4排放的冷点区（图2-4）。总体上，安溪县西部和南部的低山和中低山区域等专业化水平较高的茶区是各碳素指标分布的热点区域，一方面表明该区域土壤肥力较好，另一方面该区域温室气体排放量亦较大。

图2-4 甲烷排放量冷热点分析

图2-1 10～20cm土壤有机碳量冷热点分析

图2-2 凋落物输入碳量冷热点分析

2.4 碳素空间分异的影响因素

参考前人研究发现气候、地形、土壤初始属性等[10]因素与碳素的空间分布存在显著相关性。因此，该研究首先选择降雨、经纬度、地形因子、初始土壤属性因子（SOCmax，SOCmin，pHmax，pHmin，土壤黏粒含量最大值，土壤黏粒含量最小值，土壤容重最大值，土壤容重最小值）与各碳素指标进行相关性分析，筛选出与各碳素指标具有相关性的变量，再使用OLS对各相关变量进行线性拟合，经过t检验和筛除具有共线性的变量，最后剩余变量作为解释变量输入GWR模型并依据拟合结果进行分析。

土壤碳素中CO2、CH4是温室气体的重要来源，识别其影响因素对茶园固碳减排具有重要意义，因此该研究仅对影响茶园土壤CO2和CH4排放的因素进行探讨。经过相关性分析、OLS建模和共线性筛除，发现初始土壤属性因子对茶园CO2、CH4排放具有明显影响。GWR模型的拟合方程残差结果不具有空间自相关性，表现为随机分布。进一步利用AICc信息准则，得到CO2、CH4的R2值分别为0.73和0.79，R2adj分别为0.71和0.79，上述结果显示GWR模型能够模拟各解释变量的影响。

图2-3 二氧化碳排放量冷热点分析

2.4.1 CO2排放影响因素分析

初始土壤属性中黏粒含量最大值对茶园土壤CO2的排放影响最强，且具有明显的负向作用效果（表4），对剑斗镇、湖头镇、湖上乡、白濑乡、金谷镇等属于丘陵和平原的区域负向影响最为明显（图3）。SOCmax对茶园土壤CO2的排放正向作用效果显著（表4），其中在西坪镇、大坪乡、长坑乡等低山和中低山茶区影响作用较大（图3），考虑该区域多为专业化茶区，土壤肥沃、土壤有机碳含量大，因此通过土壤呼吸作用转换生成CO2的可能性较大。

图3 CO2排放量回归系数空间分布

表4 CO2排放量参数估计

2.4.2 CH4排放影响因素分析

初始土壤属性中黏粒含量最大值对安溪茶园土壤CH4的排放影响最强，其次是黏粒含量最小值，且两者负向作用效果明显（表5），影响程度空间格局呈现“西北高东南低”的特征（图4）。土壤黏粒含量最大值和最小值影响强度空间分布相似，但初始土壤黏粒含量最小值对长坑乡和祥华乡的负向作用效果更强。SOCmax对安溪茶园CH4的排放具有明显的正向效果影响（表5），并在安溪县西、北部的感德镇、桃舟乡、福田乡等专业化水平较高的茶区表现尤其显著（图4）。

表5 CH 4排放量参数估计

图4 CH 4排放量回归系数空间分布

该研究发现初始土壤属性因子中的土壤黏粒含量和SOC对茶园土壤CO2、CH4排放的影响在空间和作用强度上具有明显的差异性。土壤黏粒含量总体上对CO2、CH4表现为负向作用效果，这与Sass和Fisher[41]的研究表明黏粒含量高的土壤CH4排放较少的结论一致，土壤黏粒含量大将不利于土壤O2和CH4的扩散。另一方面，虽然在其他条件不变的情况下，土壤有机碳含量越多，对CO2和CH4排放的正向作用越大，但是由于土壤有机碳组成结构的差异，土壤有机碳与CH4的关系仍不明晰，部分研究者发现土壤活性有机碳含量对其中的微生物活性具有直接影响，进而影响到CH4的排放。Yagi等[42]研究认为易矿化碳与CH4排放存在线性关系，Vermoesen等[43]研究发现水溶性有机碳含量对CH4排放产生重大影响，而焦燕等[44]研究发现水溶性有机碳含量对土壤CH4排放的影响不大。该研究中土壤初始属性因子中的pH和土壤容重未对CO2、CH4产生影响，不排除由于模型模拟等造成的误差，但相关研究已经证实土壤容重[45]和土壤pH[46]等均会影响土壤CH4和CO2的排放，且排放差异受到诸多影响因子的综合作用。

由于不同耕作管理方式对茶园土壤水分等物质的入渗影响存在差异[47]，进而影响到茶园土壤碳素的动态循环，进一步研究应该关注合理的茶园管理方式来权衡茶园土壤碳素输入和温室气体排放的关系。可开展茶园最佳管理措施模拟评价[48]，以在保证茶叶产量的前提下减少茶园土壤温室气体过度排放。今后还可考虑土壤类型、翻耕、覆膜、灌溉量、肥料种类、施用量、施肥方式、施肥时间以及种植不同茶树间作物等因素对茶园土壤碳循环的影响。

3 结论与建议

3.1 结论

（1）经过参数率定和验证，DNDC模型区域尺度模拟的茶园土壤碳素结果与实测数据具有较高的一致性。

（2）安溪县茶园SOC含量丰富，土壤肥力较好。凋落物输入碳是茶园土壤碳素输入的主要来源，而土壤通过呼吸作用转化为CO2是土壤碳素输出的重要途径。

（3）茶园土壤碳素指标中SOC10～20cm量、凋落物输入碳量、CO2排放量、CH4排放量均表现为空间聚类分布。各碳素指标热点区域均主要分布在安溪县西北部和南部的低山和中低山区域。

（4）初始土壤属性中黏粒含量和SOC对安溪茶园土壤温室气体排放的影响显著，其中土壤黏粒含量最大值对CO2和CH4排放量的空间分布影响最强且表现为负向作用效果，其次是SOCmax，表现为正向作用效果。

3.2 建议

（1）应鼓励安溪县西部低山和中低山区域的专业化茶区将茶树修剪残渣等适量返园，以提高茶园土壤肥力。

（2）进一步科学调控安溪西北部、南部等区域的茶园施肥量、施肥时间等，科学制定施肥管理方式，以权衡茶园碳素输入和温室气体排放的关系，推动安溪茶产业积极响应“双碳”目标，走可持续高质量发展道路。