Analysis of seasonal spatial and temporal variation patterns of forest vegetation NEP and climate drivers in Hainan Island

海南岛位于中国的南部，热带北缘，地理位置介于108°03′～111°03′E和18°10′～20°10′N之间。海南岛的气候类型属热带海岛季风气候，长夏无冬，1年分干湿2季，全年平均温度22.5～25.6 ℃，年均降雨量900～2 500 mm^[26-27]。海南岛拥有较高的森林覆盖度，其中，热带山地雨林是海南岛热带森林植被中面积最大、分布集中的垂直自然地带性的植被类型^[28]，主要分布在吊罗山、五指山、霸王岭、尖峰岭以及鹦哥岭等600 m以上海拔的山地。而低地森林主要分布在600 m以下中低海拔区域，包括红树林、橡胶林等^[29]。

土地利用数据集源自中国研制的30 m空间分辨率全球地表覆盖数据GlobalLand30^[30]，共包含2000、2010、2020年3期的土地利用数据（http://www.globallandcover.com/）。基于海南岛的边界图裁出海南岛的土地利用数据集，然后统计近20 a海南岛的林地区域，绘制出稳定的林地区域图像，即在以上3组土地利用数据中均为林地的区域会被提取出来，形成稳定的林地分布图（图1），再将其分辨率重采样至0.01°以匹配NEP数据的空间分辨率。

本研究使用的NEP是来自全球环境研究中心（Center for Global Environmental Research，CGER），该产品提供了1999—2018年空间分辨率为0.01°，时间分辨率为10 d的全球净生态系统交换（Net Ecosystem Exchange，NEE）数据^[31] （ https://db.cger.nies.go.jp/DL/10.17595/20200227.001.html.en）。初始NEE数据格式被转换为栅格，取相反数即为本研究所使用的NEP数据。 全年共分为36幅栅格图像，再分别将每10 d的NEP均值数据合成月尺度数据和年尺度数据。

气温和降水数据来源于国家科技基础条件平台—国家地球系统科学数据共享服务平台−黄土高原科学数据中心(http://loess.geodata.cn)，该数据集是根据CRU发布的全球0.5°气候数据集以及WorldClim发布的全球高分辨率气候数据集，通过Delta空间降尺度方案在中国地区降尺度生成的，并使用496个独立气象观测点数据进行了数据验证^[32]。

DEM数据集来自中科院资源环境科学与数据中心（https://www.resdc.cn/），该数据集为基于最新的SRTM V4.1数据经整理拼接生成的90 m的分省数据。数据采用WGS84椭球投影。随后将其空间分辨率采样至0.01°，并基于1.2中提取的林地区域进一步提取林地所处的海拔高度。笔者使用ArcGIS中的空间分析工具根据海拔计算出坡度和坡向。

针对存在季节性变化的月度NEP，采用季节性分解，将其分解为趋势因子（Trend Component，TC）、季节性因子（Seasonal Component，SC）和误差因子（Error Component，EC）。TC能准确把握数据的长期变化；SC能捕捉到数据1年内的周期性变化；而EC能反映那些不能被趋势或季节效应解释的变化^[33]。通过相加模型可以表示为：

式中：$NEP_{t}$指每月的NEP，$TC_{t}$、$SC_{t}$、$EC_{t}$分别指趋势、季节性以及随机误差。

线性滤波器是估计时间序列趋势常用的方法，最常见的线性滤波器之一是滑动平均。针对NEP月度数据，选择了12点移动平均法，线性滤波器如下：

式中：$TCt$指去除季节性效应的NEP趋势值；1999—2018共20年，每年12个月，合计240个月，t从第7月开始取值。

在获取趋势因子之后，通过下式可得到季节性因子：

式中：$SC_{t}$为第$t$月时算得的季节性因子。

这个季节性因子的估计也包含了每个t时间的随机误差因子$ECt$，通过计算每个月季节性效应估计的平均值并在所有年份重复此序列来获取整体的季节性效应。

式中：${i} \text { 为年份， } {j} \text { 为月份 } $，j取1~12，$NEP_{ij}$为第i年第j月的平均NEP，$NEP_{j}$为20年的j月平均NEP。

基于像元尺度的趋势分析法能模拟研究区中每个栅格单元的变化趋势，从而反映NEP变化的方向和速率。计算公式为：

式中：Slope是NEP逐像元线性回归方程的斜率；代表年数，n为时间跨度。当Slope>0时，NEP呈增加趋势；当Slope=0时，NEP基本稳定，无明显变化；当Slope<0时，NEP呈减少趋势。

使用BRT来评估1999—2018年气候因子与地形因子对于海南岛森林NEP变化的相对影响程度。BRT分析具有容纳任何数据分布的能力，因此，在分析过程中无需进行数据的转换。在进行BRT分析之前，对森林NEP与气候因子以及地形因子进行皮尔逊相关分析和显著性分析。使用R语言中的GBM包进行BRT分析，将1999—2018年平均NEP的逐像元数据作为响应变量，同一时期的气候因子和地形因子作为解释变量。BRT的本质就是1个加法模型，每次建立模型是在之前建立模型损失函数的梯度下降方向。

基于海南岛森林的分布情况^[29]，以600 m海拔为界限，将海南岛森林分为低海拔森林和高海拔森林分别研究。笔者计算了近20年NEP趋势与气候因子趋势的相关性，如下所示：

式中：$x_{i}$为1999—2018年温度或降雨的逐像元趋势，$y_{i}$为1999—2018年NEP的逐像元趋势，为像元数。随后对NEP趋势、温度趋势以及降雨趋势构建了多元线性回归模型，并以回归系数来估计各因子的贡献度，进而确定不同条件下对NEP影响的主导因素。

对1999—2018年月度NEP数据进行时间序列分析，发现海南岛森林NEP随时间的变化趋势及季节性波动。通过使用滑动平均滤波器，得到了近20年去除季节性影响的较为平稳的NEP变化趋势（图2）。总体来看，NEP在20年间变化起伏波动，在2010年前后出现最大趋势，而在2016年出现最低趋势。笔者进一步分析了NEP年内的季节性效应，发现年内NEP呈现先增后减的季节影响，1−5月，NEP逐渐增大，在5月达到最大值，随后的6−12月，NEP逐渐减小。同理，可以得到降雨的年内季节性效应。1−9月，降雨随月份逐渐增加，在9月达到最大值，随后的10−12月，降雨逐渐减少。基于此可以将海南岛1−3月划分为干季，7−9月划分为湿季。

通过计算每年海南岛森林NEP的均值，得到 1999—2018年逐年NEP的年际变化趋势图（图3-a），整体来看NEP随时间呈现不显著的下降趋势，有机碳的变化率为−0.57 g·m⁻² ·a（P>0.05），分段来看，前10年NEP呈现不显著的增长趋势，有机碳的变化率为3.3 g·m⁻² ·a（P>0.05）；后10年NEP呈现显著下降趋势，有机碳的变化率为−8.4 g·m⁻² ·a（P<0.05）。近20年NEP的均值(有机碳)为483.23 g·m⁻² ·a，其中，有9年的NEP高于均值。在2010年，NEP达到最大值，为534.68 g·m⁻² ·a。在2016年，NEP有最小值，为439.47 g·m⁻² ·a。

同样，在1999—2018年逐年干季和湿季的NEP趋势变化图中（图3-b）可以看到湿季NEP有不显著的下降趋势，而干季NEP有不显著的上升趋势。湿季NEP的均值为136.76 g·m⁻²(有机碳)，其中有11年NEP超过均值，最大值出现在2005年，为156.21 g·m⁻²，最小值出现在2015年，为117.10 g·m⁻²。干季NEP的均值为102.09 g·m⁻²，仅有9年NEP超过均值，最大值出现在2009年，为129.54 g·m⁻²，最小值出现在2005年，为81.24 g·m⁻²。干、湿2季对比，除2009年干季NEP超过湿季NEP，其余年份湿季NEP的值均高于干季。

通过绘制全年、干季以及湿季的NEP空间分布图（图4-a、b、c）来观察NEP在空间方位上的分布。年际的NEP最高值为825.30 g·m⁻² ·a(有机碳)，最低值为−54.64 g·m⁻² ·a，其中NEP的高值集中在海南岛中部和东南部，而在海南岛北部以及海岸线附近多为NEP低值聚集。此外，大部分地区的NEP均大于0。干季的NEP最高值为236.03 g·m⁻²，最低值为−69.89 g·m⁻²，与年际NEP相似，干季NEP的高值多分布在海南岛中部五指山、霸王岭以及尖峰岭一带，NEP低值则聚集在海南岛东北部。湿季NEP的分布较为破碎化，高值和低值相间，其中最高值为211.84 g·m⁻²，最低值为−19.74 g·m⁻²，湿季NEP最高值较干季减少了24.19 g·m⁻²，但最低值较干季增加了50.15 g·m⁻²。

应用一元线性回归对海南岛NEP近20年的年际趋势变化进行分析，进而绘制了NEP的趋势变化空间分布图（图4-d, 图4-e, 图4-f）。年际的趋势变化如图4所示，NEP呈现显著下降趋势的区域占比为18.49%（P<0.05），NEP呈现显著增长趋势的区域占比为17.99%（P<0.05），NEP无显著变化的区域占比为63.52%。从分布来看，NEP显著增长区域主要集中在海南岛东北部，而显著下降区域主要集中在海南岛中部。干季的趋势变化：NEP呈现显著下降趋势的区域占比为6.39%（P<0.05），NEP呈现显著增长趋势的区域占比为23.05%（P<0.05），NEP无显著变化的区域占比为70.56%。在空间分布上，干季的NEP趋势变化与年际的趋势变化相似，但是呈显著下降趋势的区域明显减少。湿季的趋势变化：NEP呈现显著下降趋势的区域占比为23.66%（P<0.05），NEP呈现显著增长趋势的区域占比为12.65（P<0.05），NEP无显著变化的区域占比为63.69%。与年际趋势相比，湿季的NEP显著增长趋势占比降低了5.43%，显著下降趋势占比提高了5.17%；与干季趋势相比，湿季的NEP显著增长趋势占比降低了10.40%，显著下降趋势占比提高了17.27%。

通过对不同季节的NEP与气候因子、地形因子进行皮尔逊相关分析以及显著性检验，发现年际NEP与气温、降雨、海拔以及坡度均有极显著的相关性（P<0.01），相关系数分别为−0.53、−0.03、0.57、0.50，而与坡向无显著相关性（图5-a）；干季NEP与气温、降雨、海拔以及坡度均有极显著的相关性（P<0.01），相关系数分别为−0.54、−0.07、0.73、0.66，而与坡向无显著相关性（图5-b）；湿季NEP与气温、降雨、海拔、坡度以及坡向均有极显著的相关性（P<0.01），相关系数分别为0.07、−0.04、−0.09、−0.10、−0.03（图5-c）。基于以上分析，选择气温、降雨、海拔和坡度作为自变量因子用于BRT分析，以便计算各因素在不同季节对NEP的贡献程度。在进行参数调优后，在R语言中使用GBM包进行BRT分析，得到不同季节不同因子对NEP的相对贡献率。如图6-a所示，对年际NEP影响程度最大的是海拔，其相对贡献率为45.46%，其次分别为降雨、气温和坡度，相对贡献率依次为25.78%、16.17%、12.59%；如图6-b所示，对干季NEP影响程度最大的是海拔，其相对贡献率为40.58%，其次为坡度、降雨和温度，相对贡献率依次为28.13%、16.12%、15.17%；如图6-c所示，对湿季NEP影响程度最大的是降雨，其相对贡献率为30.75%，其次为温度，海拔和坡度，相对贡献率依次为26.77%、21.88%、20.60%。除湿季外，海拔在年际NEP和干季NEP中扮演着重要的角色。

为进一步说明海南岛森林NEP的主要影响因素，分别对高海拔区域和低海拔区域进行了NEP趋势和气候因子（气温、降雨）趋势的相关性分析。结果显示，在海拔600 m以上的森林生态系统中，全年的NEP趋势和湿季的NEP趋势与降雨趋势均有显著的负相关性，相关系数分别为−0.11和−0.19（图7-d, 图7-f），而干季的NEP趋势则与降雨趋势没有显著的相关性（图7-e）。不论是年际的NEP趋势还是干季、湿季的NEP趋势均与温度的趋势没有显著的相关性（图7-a, 图7-b, 图7-c）。而在海拔600 m以下的森林生态系统中，则呈现出十分不同的模式。年际的NEP趋势对于温度趋势有微弱的显著负相关性，相关系数为−0.02（图8-a），对于降雨趋势有极显著的正相关性，相关系数为0.46（图8-d）；干季的NEP趋势与温度趋势没有显著的相关性（图8-b），与降雨趋势有极显著的负相关性，相关系数为−0.06（图8-e）；湿季的NEP趋势与温度趋势、降雨趋势均有极显著的正相关关系，相关系数分别为0.11和0.30（图8-c,f）。如图9所示，多元线性回归结果表明，在干季600 m以下的森林生态系统中，降水趋势对NEP趋势有显著的负贡献，贡献率为−53%。在湿季600 m以下的森林生态系统中，温度趋势对NEP趋势有显著的正贡献，贡献率为90%。

海南岛森林NEP表现出了季节性变化，但是与亚马逊地区不同的是，亚马逊地区在季节性变换时会存在碳动态的转变，如碳源和碳汇的相互转换^{[8, 11-12,34-35]}，而海南岛热带森林不论是从年际角度出发，还是从干季、湿季角度出发，均呈现稳定碳汇状态。所不同的是，海南岛森林NEP在湿季时普遍要比干季高。地理位置差异所带来的气候类型差别可能是海南岛热带森林和亚马逊热带雨林呈现不同碳动态格局的原因。亚马逊地区属热带雨林气候，长年高温多雨，对干旱的敏感性更高^[36]；而海南岛则属热带海岛季风气候，受季风影响分明显的干季和湿季。相对固定的季节模式，使得海南岛森林固碳动态虽有季节性变化，但总体趋势保持稳定。更为重要的一点是，不只是气候变化的影响，农业扩张所带来的森林砍伐，火灾和干旱之间的相互作用也是造成亚马逊森林碳损失的原因^[36-37]。相比而言，海南岛的森林保存的较为完整。此外，本研究采用基于多期土地利用遥感数据提取稳定林地的方法，一定程度上规避了森林变化或者损失而带来的净生态系统生产力的损失。

从趋势的相关性上来看，相比气温，降雨与NEP有着更为显著的相关性。但是在不同海拔下，降雨与NEP的相关性表现出了相反的模式。在海拔600 m以下，年际、湿季降雨趋势分别与年际、湿季NEP趋势呈现出显著的正相关，但是在海拔600 m以上，又呈现出截然相反的显著负相关性结果。群落结构和生态系统过程往往随海拔梯度变化^[24]，一般而言对海拔的响应都是由于温度变化所驱动，但也不绝对，降雨等因素也会随海拔梯度变化^[24,38]。海拔梯度下，海南岛降雨趋势的变化不是十分明显，所以这种结果可能是因为群落结构发生变化。海南岛海拔600 m以上保存较为原始的热带山地雨林，在年际和湿季雨量较为充沛的阶段，随着降雨趋势的增加，附生植物旺盛生长会导致NEP下降。而在低海拔区域，由于人为干扰，林型更为复杂，人工林的固碳能力相较高海拔原始林更为强大。

仅从气候因素考虑，在旱季时，降雨对低海拔森林NEP变化有显著的负向作用，而在湿季时，温度对低海拔森林NEP变化有积极作用。气温和降雨是影响生态系统生产力的主要的因素^[39]，但是随着全球变化加速，陆地生态系统的碳交换对于气候变化的响应和反馈仍旧存在不确定性^[40-41]。尤其在热带森林生态系统，降雨量在年内的波动导致季节性干旱的产生，干旱事件对热带森林生态系统树木的生长以及生态系统的功能有着严重的影响^[42]。低海拔区域更多为人工林，以橡胶林为例，干旱会导致其落叶甚至死亡，光合作用受到抑制^[43]。另一方面，旱季的降雨会增加土壤水分^[44]，进一步影响了土壤呼吸^[45]。雨季时，温度呈现出显著的正向作用。降水的充沛使得水分并不是NEP的主要限制因子，而此时的温度却有着积极的调节能力，与NEP变化呈现出一致性。研究表明，温度的升高是利于生态系统呼吸^[16]、生态系统总产量^[46]以及净碳吸收^[47]。可能在雨季时，由于不再受到水分条件的限制，海南岛森林植被得到较大恢复，处于生长阶段的旺盛期，此时温度对于整个生态系统的光合作用影响远远大于呼吸作用。

近20年海南岛森林NEP的年际变化波动较大，存在明显的季节性。湿季的NEP高于干季的NEP。空间上，NEP呈显著增长趋势主要分布在海南岛东北部，显著下降趋势主要分布在海南岛中部。不同季节NEP与气候因子的关系有所差异，在干季低海拔区域，降雨对NEP有显著的负向贡献，贡献度为−53%（P<0.05）;在湿季低海拔区域，温度对NEP有显著的正向贡献，贡献度为90%（P<0.05）。