海南岛近20 年GPP 变化格局及驱动因素分析

雷济舟，崔 嵬，朱济帅，张润卿，赵俊福，章 杰，张 翔，孙仲益,5

（1. 海南大学 生态与环境学院，海口 570228; 2. 国家林业和草原局发展研究中心，北京 100714;3. 海南长光卫星信息技术有限公司，海口 570311; 4. 海南省生态环境监测中心，海口 571126;5. 海南省农林环境过程与生态调控重点实验室，海口 570228）

陆地生态系统（以下称为陆生系统）总初级生产力（Gross Primary Production，GPP）是指大气中的CO2通过植被光合作用进入陆生系统并转化为有机碳的总量[1]，在一定程度上决定着碳汇时空分布格局[2]，影响着生态系统物质及能量流动[3]；为实现区域碳平衡，减缓气候变化起到关键作用[4]，具有重要意义[5]。但受到生态系统及其与环境要素间相互作用的影响，GPP 的时空分布具有较大的变异性与异质性特征[6]，因此，对GPP 分布格局的驱动因素定量分析将会有助于提升区域碳循环的理解，也更有助于揭示区域生态系统对人类活动与区域气候变化的响应与反馈。气候变化以及人类活动作为影响热带地区GPP 分布格局的主要驱动因素，相关研究已取得丰富成果，对于前者主要集中于水分[7]、温度[8-9]、CO2浓度[10]以及太阳辐射[11]等气象因素对GPP 时空格局的影响[12]，已通过机理模型[13]、遥感模型[14]及统计方法[15]进行了量化分析与机理解释。而后者，相关研究也证实了诸如城市建设[16]、退耕还林还草[17]、刀耕火种[18]等人类活动会显著改变区域GPP 分布[19]，同时土地利用及覆被变化（Land Use and Cover Change，LUCC）作为人类活动的直接体现[20]，是最为经典的表征指标[21]。热带地区水热条件充足，各种环境因素都会对GPP 产生影响，许多研究对此进行了讨论，其中在亚马逊热带地区通过机理模型揭示出辐射是生产力变化的主导因素之一[22]；热带地区水热条件虽然相对充足，但时空间的不均匀分配使得年际降水对GPP 变化产生较为关键影响[23]，同时GPP 对温度的变化也有着较高的敏感性[24]。在全球尺度的研究中已有结果表明，温度、辐射以及水分对于GPP 的影响程度分别为13.07%、-7.24%和11.74%，具体表现为低纬度地区水分对于GPP 变化起到主导作用，中高纬度地区温度起到主导作用[6]。随着人类在热带地区的活动愈发频繁[25]，越来越多的研究者发现开荒等农业生产行为致使的土地利用变化才是导致GPP 改变的主要驱动因素，这一观点得到在亚马逊雨林、东南亚热带雨林等地的相关研究支持[25-26]。综上所述，热带地区植被结构复杂，环境因素多变，虽然科学界对GPP 时空分布格局的驱动因素及响应机理已基本达成共识，但主导因素及其相对贡献大小依旧为争论焦点。

海南岛是我国第一大热带岛屿，作为相对独立的地理单元是进行热带地区GPP 研究的理想切入点；并且近年来随着海南自由贸易港等政策的落实[27-28]，海南岛建设程度加深，人才引进以及城市扩张等一系列活动导致土地覆盖变化剧烈，为探究人类活动与气候变化对GPP 的影响提供了契机。因此，本研究选择海南岛作为研究靶区，对其近20 年GPP 动态变化趋势进行归因分析，具体目标为：探究2001—2019 年间海南岛GPP 时空变化格局；量化气象要素与人类活动对于年际间GPP 变化的贡献率，揭示主导因素。旨在为解明热带地区生态系统响应区域气候变化提供科学参考也为海南生态文明试验区的政策制定提供理论依据。

1 材料与方法

1.1 研究区概况海南岛（18°10′—20°10′N，108°37′—111°03′E），地处热带北缘，属热带季风海洋性气候，温差较低，全年高温，年平均气温22.5～25.6 ℃。年平均降雨约为1 640 mm，降水充足但分配不均，雨旱两季分明。海南岛地势为四周低平，中间高耸，呈穹隆山地形，中部山区以五指山、鹦哥岭为隆起核心，向外围逐级下降，由山地、丘陵、台地、平原构成环形层状地貌，梯级结构明显，土地覆盖类型丰富，以森林、农田以及草原为主，其中，森林生态系统包括常绿阔叶林、常绿针叶林、落叶阔叶林以及混交林等。

1.2 研究数据

1.2.1 土地利用类型数据海南岛土地利用类型数据集来源于MODIS（Moderate-resolution Imaging Spectroradiometer）Land Cover 产品（MCD12Q1），时间分辨率为1 a，空间分辨率为500 m，本研究时间范围为2001—2019 年[29]。本研究采用了IGBP（International Geosphere Biosphere Programme）分类方案所确定的17 个土地覆盖方案。结合相关研究[30]以及本研究实际情况，对土地覆盖类型进行重分类，重分类情况见表1。

表1 根据IGBP 分类系统对土地利用类型的重分类

1.2.2 GPP 数据2001—2019 年海南岛GPP 数据集为MODIS 产品MOD17A2HV006。该产品为4 级标准产品[31]，主要通过光能利用率模型进行计算，逐8 天合成空间分辨率为500 m 的GPP 产品，从区域尺度到全球尺度已经被广泛应用[32]。本研究将GPP 数据求和处理成月值尺度后转化为年尺度进行后续分析。

1.2.3 气象数据本研究所使用2001—2019 年逐月Ta、PAR 以及相对湿度（Relative Humidity，RH）气象数据，其中Ta、RH 数据均源自于国家地球系统科学数据中心（http://www.geodata.cn/），其空间分辨率均为1 km，是利用全国基本气象站观测数据进行空间插值并综合地形数据所绘制的空间分布数据[33-34]。根据Ta 与RH 数据并结合Tetens 经验公式[35]对海南岛逐月VPD(饱和水汽压差)的数据进行计算；PAR(光合有效辐射)的数据为Global Land Surface Satellite（GLASS）产品（http://www.glass.umd.edu/index.html），其空间分辨率为0.01°，时间分辨率为1 d 是基于多源遥感数据和地面实测数据进一步反演得到的长时间序列、高精度的全球地表遥感产品[36]。

1.3 研究方法

1.3.1 技术路线首先以海南岛年际间用地类型未改变的栅格作为数据样本，利用其多年VPD、Ta 及PAR 年距平值作为解释变量，GPP 年距平均值作为目标变量，通过高斯过程回归（Gaussian Process Regression， GPR）与随机森林（Random Forest，RF）两种机器学习算法构建预测模型，GPR 是常用的监督分类学习方法，是一种基于贝叶斯方法的非参数概率模型，回归的目的是通过学习样本，经过训练得到输入变量与输出变量之间的函数关系，常常用于小样本回归分析。RF 是经典的基于分类和回归树的集成学习算法，可以解释若干自变量对因变量的作用；其次，选择较好的机器学习算法后，利用年际间VPD、PAR 以及Ta 的差异作为输入数据，对其年际间气象要素差异所引起的GPP 差异进行估算，完成气象因素对GPP 影响的相对贡献率的计算；最后，以年际间用地类型发生改变的栅格作为研究目标，计算相邻年间LUCC 改变的栅格其GPP 的变化值，进而完成LUCC 对GPP 影响的相对贡献率的计算；最终进行海南岛GPP 变化主导因素的探讨（图1）。

图1 研究流程图

1.3.2 机器学习模型构建及评价本研究通过GPR 及RF 机器学习算法构建模型以分析气象因素对GPP 的相对贡献率。如表2 所示，海南岛近20 年间各用地类型均有不同程度改变，裸地与湿地类型未变栅格所占比例较低，用于训练模型的样本数量较少，因此选择GPR 算法；而林地、草地、耕地等其他类型用地样本数量充足，利用RF 能够取得较好效果。所有用地利用的模型训练均采用五折交叉验证进行精度的评估，以决定系数（Coefficient Of Determination，R2）、均方根误差（Root Mean Square Error，RMSE，记作SE）作为模型精度评价指标。

表2 模型评价精度表

1.3.3 归因分析方法为探究LUCC 以及气象因素对于海南岛GPP 的影响程度，本研究中LUCC以及气象因素相对贡献度计算公式如下：

式中，Con(LUCC)(i,i+1)为第i年至第i+1 年间LUCC对于GPP 变化的相对贡献率，GPPi为i年海南岛GPP 总值，GPP(LUCC)(i,i+1)为i至i+1 年间LUCC所驱动的GPP 变化值；Con(Climate)(i,i+1)为i年至第i+1 年间气象因素对GPP 变化的相对贡献率，f为1.3.1 所介绍的机器学习模型。VPD(i,i+1)、Ta(i,i+1)以及PAR(i,i+1)分别为第i年至第i+1 年间VPD、Ta 以及PAR 的变化量，其中，i为年份；n为各土地利用类型，n=1，2……7。

2 研究结果

2.1 海南岛GPP 变化情况海南岛时间上，近20 年期间GPP（以C 表示，下同）时间变化如图2-a所示，GPP 年际变化明显，总体上研究时期GPP 呈现上升趋势（0.44 Tg·a-1）。2005 年海南岛GPP 值最低仅为56.46 Tg C，2017 年后GPP 值持续增加，在2019 年达到峰值74.12 Tg C。其GPP 主要由3 种土地利用类型构成，林地为GPP 值最高的土地利用类型，约占到全年GPP 总值的85%，19 年期间GPP 值以0.46 Tg·a-1增加；耕地作为GPP 值仅次于林地的土地利用类型，其年GPP 值最大占比达到10.7%，耕地19 年期间变化幅度不明显，其年GPP 值以0.07 Tg·a-1增加；草地GPP 值在19 年期间持续降低（-0.11 Tg·a-1）；空间上，近20 年间海南岛87.8%的面积呈现出GPP 增加趋势(图2-b)；其中，海南岛北部以及东北部是GPP 极显著增加的集中区，中部山区部分地区GPP 无明显变化，而海口市与三亚市周边小部分区域的GPP 表现为减少趋势，且碎片化程度较高，占比为8.9%。

图2 海南岛GPP 时空间变化趋势

2.2 海南岛LUCC 变化情况对近20 年海南岛土地利用变化数据进行分析（表3），研究时期内转移面积小于1 km2的分支不进行展示。结果显示，2001—2005 年，海南岛土地利用转移总面积为2 736.94 km2。其中，林地的转出面积最大，净转出面积为309.13 km2，主要转化为耕地；2005—2010 年土地转移总面积为4 164.04 km2，草地转化成林地以及耕地转化为林地为最主要的转出方式；2010—2015 年期间，海南岛的土地覆盖变化主要发生在林地、草地以及耕地的地区面积分别为1 188.10 km2、1 074.56 km2以及1 070.85 km2；2015—2019 年最主要的转移方式是林地转化为耕地和草地转化为林地，其转化净面积分别为541.25 km2和214.65 km2。

表3 海南岛近20 年间土地利用每5 年的转移情况

2.3 LUCC 对GPP 变化的影响研究时期，海南岛LUCC 所驱动的GPP 变化主要发生在林地、草地以及耕地地区（图3）LUCC 所驱动的GPP 年际变化明显。其中，林地地区LUCC 所驱动的GPP呈增加趋势，在2018—2019 年林地地区LUCC 所驱动的GPP 达到最大值8.41×10-2Tg C，占比57%；草地地区19 年期间呈显著微弱增加趋势，2005—2006 年草地地区达到LUCC 所驱动的GPP峰值7.56×10-2Tg C，2004—2005 年草地地区LUCC 所驱动的GPP 仅为-5.99×10-2Tg C；耕地地区在近20 年内也是呈现出微弱增加趋势，在2008—2009 年LUCC 所驱动的GPP 值最大。

图3 各土地利用类型在LUCC 所驱动GPP 变化

2.4 相对贡献率变化LUCC 以及气象因素对于海南岛GPP 相对贡献率的年际变化如图4 所示，年际间LUCC 与气象要素的相对贡献率同频率较高，但程度相差较大，气象因素主导海南岛年际间GPP 变化。LUCC 的相对贡献率近20 年期间波动较为剧烈，其中2008—2009 年期间相对贡献率达0.31%，而在2004—2005 年LUCC 对于GPP 的相对贡献率最小为-0.24%；气象因素对年际间GPP 变化的相对贡献率逐年差异明显，2005—2006 年其正向相对贡献率达到最大值10.79%，而在2010—2011 年气象要素对于GPP 的影响程度为历年来最低值-14.99%。

图4 相对贡献率变化图

3 讨 论

3.1 LUCC 对GPP 的影响LUCC 一直是陆生系统GPP 的主要驱动因素。林地、草地以及耕地地区是海南岛LUCC 发生的主要区域。研究时期LUCC 共变化15 528.40 km2。其中，植树造林6 749.33 km2，累积GPP 达到0.23 Tg C，这与相关研究结论[37]基本一致即植树造林对于GPP 增长有着较为良好的促进作用。这一现象的原因可能是椰林[38]、橡胶林[39]等林种具有较高的经济价值。此外，近52%的新林地资源主要来源于草地，改善了海南岛植被覆盖，强化了林地生态系统的防线保护作用；耕地也在此过程中净增加679.41 km2，GPP 在此过程中增加1.63×10-2Tg C。耕地的快速扩张一定程度改变原有的生态环境，由单一经济作物所构成的生态系统抵抗力稳定性降低可能会在未来引发水土流失，所以如何合理规划耕地资源是要深入考虑的。

LUCC 对GPP 年际间变化的相对贡献率远低于气象因素，其一在于，海南岛林地面积高（79.1%），并可占全岛GPP 的84.8%，并且主要集中在中部山区，由于天然林保护等政策[40-41]，中部山区人类活动强度较低，所以在近20 年来，林地并未发生较大的面积变化（近20 年变化不足2%）；二是由于MCD12Q1 产品的将海南岛土地利用与覆盖分为14 类，而为方便分析本研究将其重分类为7 类，不同分类系统所造成的相对贡献率影响(图5-a)最高可相差3 倍（2004—2005 年）。尺度效应也是重要因素，如图5-b 所示，将海南岛近20 年逐年逐市县进行拆解，并使用IGBP 原始分类系统，LUCC 在各市县中对年际间GPP 变化的相对贡献率普遍高于全岛平均值，LUCC 相较于气象因素而言作用范围小，因此在进行全岛LUCC 相对贡献率计算时不同地区间差异相互抵消，因此也造成LUCC 的相对贡献率较低。LUCC 相较于气象因素，对于小尺度范围的GPP年际间变化起到主要作用。LUCC 与气象因素对于年际间GPP 变化的相对贡献率具有较高的同频性，其原因在于相对贡献率的算法以前一年GPP 为基准，而非年际间GPP 差异为基准，好处在于避免出现年际间GPP 相差极小而导致的极端贡献；但负效应变为同频效应以2005—2006 年为例，极端事件第18 号热带风暴“达维”对全岛GPP 产生显著负面影响[42]，GPP 为近20 年最低值，因此2005—2006 年会产生极高的LUCC 与气象因素的正向相对贡献，各市县贡献率受整体变化所影响。

图5 土地分类对相对贡献率的影响及各市(县)相对贡献率变化

3.2 气象因素对GPP 的影响本研究仅选取Ta、VPD 以及PAR[43]作为影响GPP 年际间变化的主要气象驱动因素，未考虑如大气CO2浓度、氮沉降等环境因素，这为本研究带来一定不确定性，但根据Sun 等人研究表明[6]，环境要素中气象因素是环境因素中主导GPP 年际间变化的关键因素，因此不确定性能够最小程度降低。此外，本研究选取VPD 而非降水作为水分条件，主要由于像海南岛这类低纬度湿润地区，大气的湿润程度是限制植物生长的主要因素[44]。其次，对于气象因素间的交互作用，本研究所采用机器学习算法（GPR 与RF），属于数据驱动模式，能够减弱自变量间的交互关系[15]。

气象因素相对贡献较高，分别选取正负相对贡献率最高的2008—2009 年（2005—2006 年为最高正相对贡献，但主要受极端事件“达维”影响，因此选取2008—2009 年）与2010—2011 年为例；2008—2009 年，PAR、Ta 及VPD 均呈现超距平的转好态势；而2010—2011 年PAR 与Ta 呈现出远超平均年际间变化的降幅，同时水分制约也略有所增加（图6）；可见本研究所提出的基于机器学习的数据驱动算法能够捕捉气象要素对GPP 的影响。除气象要素与LUCC 外，人类的管理与经营活动，如耕地农业管理及政策支持[45-46]、橡胶林下经营[47]、天然林的保护[48]、极端事件以及病虫害等异常事件也会对海南岛GPP 产生影响[49]。

图6 海南岛各气象要素变化图

3.3 不确定性分析首先本研究所使用的GPP与LUCC 数据来源于MODIS 产品，MODIS 作为光学遥感产品应用于热带地区时不可避免精度会降低[50]，这是传统光学遥感产品共同面临问题。MODIS 的LUCC 产品在全球尺度的精度为73.6%[51]，本研究将可能出现异物同谱的相似类进行了重分类，共7 大类，这能够大大降低其分类精度所带来的不确定性(图5-a)，这里不能够忽视。其次，MODIS 的GPP 产品基于光能利用率算法，是现阶段主流大区域尺度的遥感GPP 算法，包括EC-LUE[52]、CASA 以及VPM[6]等多种形式，其中可吸收光合有效辐射比率（fAPAR，fraction of Absorbed PAR）是模型核心部分，无论何种形式模型均将fAPAR 考虑为关于NDVI 等植被指数的函数；长时间序列植被指数主要源于MODIS 产品，因此本研究中直接使用了MODIS 的GPP 产品，避免二次计算时增加不确定性。此外关于计算方案引入的不确定性主要为，当年际间LUCC发生改变时，对GPP 的影响则全部归于LUCC 的贡献；但实则气象要素依旧对该区域产生影响，可此部分的影响并未考虑到本研究中，低估了气象因素造成的总体影响。综上，鉴于本研究的时间范围与尺度，结合各遥感产品的精度，最终确定MCD12Q1 为研究数据，该数据产品在国内的研究已广泛应用于国内相关研究[53]，并取得较为可靠的成果。考虑到海南岛全年较高的云覆盖，产品精度有限。未来要进一步研究LUCC 对GPP 影响等问题时，更高精确度土地利用数据产品是不可或缺的；本研究使用机器学习模型虽然可以在一定程度上将气象因素与人类活动进行解耦并单独分析其对GPP 的影响，实际环境中GPP 的影响是多元驱动的，无法对影响因素的共同作用进行驱动分析。本研究中假设LUCC 发生改变时GPP的变化归因为LUCC 所引起，将气象因素作用归于LUCC 贡献中，这也带来了一定不确定性。