基于HASM的中国植被NPP时空变化特征及其与气候的关系

赵苗苗 ，刘熠 ，杨吉林, ，刘羽 ，赵娜 ，岳天祥*

1. 中国科学院地理科学与资源研究所/资源与环境信息系统国家重点实验室，北京 100101；2. 中国科学院大学，北京 100049；3. 中国科学院地理科学与资源研究所/陆地表层格局与模拟院重点实验室，北京 100101

在过去的半个世纪里，气候变化已经引起了人们的广泛关注。在自然和人为因素的共同影响下，全球气候和大气组成发生了明显的变化，显著改变了陆地生态系统的过程、结构和功能。大气中温室气体浓度富集导致过去100年地球表面气温平均增加了 0.74 ℃，预计在未来 100年还将升高 1.4－5.8 ℃（IPCC，2013）；同时气候变化也改变了全球降水模式，增加时空变异，伴随着干旱、洪涝等极端气候事件频发（Folland et al.，2001），而这些改变均会对陆地植被净初级生产力（net primary productivity，NPP）产生影响，并由之综合表现出来（Nemani et al.，2003）。

植被净第一性生产力指绿色植物在单位时间、单位面积上所积累的有机干物质总量，作为陆地生态系统碳循环的重要组成部分，NPP直接反映绿色植物在自然环境条件下的生产力和地表质量状况，成为判定生态系统碳源/汇和调节生态过程的主要因子（Field et al.，1998），对于研究全球气候变化具有重要意义。NPP主要受自然因素（气候、土壤、植被等）和人类活动（土地利用变化等）等因素的影响（Potter et al.，1993）。国外学者在NPP时空特征变化及其与气候因子的相关关系方面已经做了很多工作（Melillo et al.，1993；Knapp et al.，2001），国内有关这方面的研究也有很多重要的成果（方精云等，2003；高志强等，2004；Gao et al.，2009；王琳等，2010）。研究表明，过去几十年的NPP受气候变化的影响呈增加趋势（Fang et al.，2001；朱文泉等，2007）。

NPP估算方法中，由于气候观测站点数量有限，传统的基于采样点数据估算NPP的方法易受空间尺度限制，较大尺度的生态系统与气候因子关系研究以模型模拟和遥感观测为主要的研究方法。过程模型理论框架完整、结构严谨、机理清晰，能较好地模拟生态系统物质循环和能量流动，是全球和区域尺度上植被生产力评估的重要手段之一（Lu，2000；冯险峰等，2014）。已有不少过程模型对区域尺度的净初级生产力进行了模拟，以生物地球化学模型和遥感模型为主（朱文泉等，2007；仲晓春等，2016）。近年来，动态植被模型以其机制合理、过程全面引起了广泛重视，应用动态植被模型模拟分析气候对自然植被的可能影响，已经成为全球变化生态学相关研究的一种发展趋势（傅伯杰，2010）。与此同时，利用MOD17A3 NPP遥感数据产品研究植被生产力的时空分布特征及与气候要素的关系已经得到广泛应用和验证（王轶虹等，2017；李登科等，2018）。将NPP动态植被模型模拟结果与遥感数据产品结合，既可以解决遥感观测因天气、云量造成的数据缺失和反演过程中的数据偏差问题，又可以弥补动态植被模型缺乏人类活动影响的不足，在一定程度上提高了结果的可靠性。

为更好地了解近年来中国陆地植被NPP对气候变化响应的空间分异规律，本研究利用动态植被LPJ-GUESS模型模拟 2001－2015年中国陆地生态系统植被NPP的时空分布，并基于HASM降尺度方法与MOD17A3 NPP数据融合。根据气候要素和土壤质地九大分区，分析NPP动态变化趋势及其与气候因子间的相关性，并与其他研究结果进行比较，为评价中国陆地植被生产力在全球碳循环中的地位和作用提供科学依据，有利于认识气候变化对中国陆地生态系统植被的影响，对中国不同陆地生态系统的生产力评估、对气候变化的响应和适应和生态系统的环境建设等工作具有重要的参考价值。

1 研究数据和方法

1.1 数据来源

本研究所利用的气象数据由中国气象数据网（http://data.cma.gov.cn/）提供，包括各气象站点的经度、纬度、月均温、月降水量、月平均辐射、CO2浓度以及土壤类型等，共752个地面观测站，时间跨度为2001年1月－2015年12月。利用HASM高精度曲面建模方法对气象站点数据进行插值，叠加中国边界图裁剪获得各气象要素栅格图，空间分辨率为 0.1°×0.1°。

2001－2015 年MOD17A3 NPP数据来自美国NASA网站（https://ladsweb.nascom.nasa.gov/data/search.html），空间分辨率为1 km×1 km。利用MRT（Modis Reprojection Tool）软件对下载的NPP数据进行重投影、拼接处理，为保证数据的匹配精度，所有数据均采用Albers等积圆锥投影。

根据气温、降水等气象要素和土壤质地等将中国分为九大地区，分别为内蒙及长城区、甘新区、东北区、黄土高原区、青藏区、西南区、华南区、黄淮海区、长江中下游区（Yue et al.，2013）。

1.2 研究方法

1.2.1 HASM

高精度曲面建模方法（HASM）是近年来针对困扰曲面建模的误差问题和多尺度问题发展起来的一种基于微分几何学曲面理论的曲面建模方法，已经被成功运用于土壤插值领域、气候模拟和建立DEM 等领域（Yue et al.，2007；Shi et al.，2009；Yue et al.，2010）。

根据曲面的基本定理，曲面 z=f(x, y)是由其第一个和第二个基本变量（Henderson，1998）所决定的。第一个基本系数 E、F、G描述了关于几何性质的信息，如切向量的角度、曲线的长度、区域的面积等等。E、F、G公式如下：

第二个基本变量L、M、N反映了表面的局部扭曲（Liseikin，2004；Toponogov，2006），公式如下：

曲面的基本变量 E、F、G、L、M、N需满足下列高斯方程（Liseikin，2004；Somasundaram，2005；Toponogov，2006；Yue，2011）：

利用标准的有限差分离散化方案（Zhao et al.，2014），可将高斯方程转变为：

式中，h是步长；f是模拟值；i、j表示空间位置；n表示迭代次数。

分别用 A1、A2、A3表示公式（5）中方程左端构成的系数矩阵，用d、p、q表示方程组右端的常数向量，并结合采样点信息，将HASM转化为如下最小二乘问题：

式中，Zn+1为研究区域在空间中离散后的格网点；k为采样信息；S为表示采样位置信息的系数矩阵。

通过引入权重参数 λ，上述约束最小二乘问题可以转化为：

该最优化问题等价于：

在实际应用HASM模型时，需要提供的输入参数为驱动场和精度控制点，其中驱动场可以是其他模型输出的或者根据其他变量计算得到的精度较低的空间变量分布，而精度控制条件则是实际测量、精度较高的数据。在本研究中，将LPJ-GUESS模型模拟输出的NPP的模拟值重采样为1 km分辨率数据，采用 10折交叉验证方法随机抽取 10%MOD17A3 NPP数据为优化控制条件，剩余的90%作为真值验证重采样的降尺度数据精度，不断重复这个步骤，直至模拟精度在误差范围内。运用HASM降尺度方法既保留了曲面整体趋势，又增加了曲面细节信息。

1.2.2 LPJ-GUESS模型

LPJ-GUESS模型是从机理上模拟植物生理生态过程（光合、呼吸、蒸散发等）、生物物理过程（物质和能量交换）以及生物地球化学过程（碳水循环）的动态植被模型（Smith et al.，2001）。通过模拟物种的生长、死亡和更新过程反映植被动态，很大程度上受当前环境状况、现有物种组成以及物种生活规律的影响（Smith et al.，2011）。模型的输入数据包括月均温、月降水量、月平均辐射、CO2浓度以及土壤类型代码等。输出数据包括植被的NPP、生物量、叶面积指数（LAI）、土壤碳储量、蒸散发、径流以及土壤含水量等。

LPJ-GUESS模型适合于多尺度生态系统过程的动态模拟，并已在全球范围内得到广泛的应用，如 Smith et al.（2001）和 Hickler et al.（2004）模拟了欧洲和美国东北的植被动态，Tagesson et al.（2009）结合遥感数据模拟了瑞典森林生态系统的净初级生产力（NPP），Piao et al.（2013）利用LPJ-GUESS模拟了GPP、NPP、NBP和温度敏感性对气候变化和 CO2趋势的响应。LPJ-GUESS模型也被用于模拟中国陆地生态系统及北京山区森林的碳循环（Piao et al.，2010；Liu et al.，2009）。

在应用LPJ-GUESS模型进行模拟时，先用2001－2015年的年均气象观测数据运行1000年，使植被生态系统和土壤结构达到平衡状态，各生态系统变量趋于稳定，且NPP和异养呼吸量（HR）、凋落物量（LT）相等。LPJ-GUESS模型平衡以后，用2001－2015年的观测数据进行动态模拟。

1.2.3 NPP年际变化

2001－2015 年每个像元NPP值年际变化趋势分析采用一元线性回归趋势线法，回归直线的斜率采用最小二乘法求得（王轶虹等，2017）。

式中，n为拥有相同坐标的栅格数即年份，n=15；Yi为第i年的像元NPP值；θslope表示一元线性回归线的斜率，θslope＞0说明NPP在2001－2015年间的变化趋势是增加的，反之减少。

1.2.4 NPP与气候相关分析

基于2001－2015年每个像元的NPP与对应的年均温和降水的相关系数计算，每个栅格包括 15对数据（龙慧灵等，2010）。

式中，变量i为年序号；n取值为15；Yij为第i年第j个像元的NPP数据；Cij为第i年第j个像元的气象因子数据；RYCj为第j个像元NPP与气候因子的相关系数。

1.2.5 假设检验

计算检验统计量tr的公式如下：

若tr在0.05水平上大于临界值，则说明两变量间有显著关系（仲晓春等，2016）。

运用Excel 2016和ArcGIS 10.2进行绘图。

2 结果与分析

2.1 植被NPP的时空分布

2.1.1 近15年全国植被NPP的空间分布

由图1可知，近15年全国植被平均NPP空间分布地域性差异明显，总体上沿水热梯度呈由西北向东南逐渐增加的趋势。其中，西北沙漠地区和青藏高原腹地的植被NPP值最小，大都在100 g·m-2（以 C计，下同）以下，面积约占8.95%。其次是其周边地区，植被NPP值在100－300 g·m-2之间，面积约占25.89%。从西北内陆向东南沿海过渡带，植被NPP值在300－600 g·m-2之间，区域范围较大，约占50%以上。四川盆地大部分地区及海南省和东北平原局部植被NPP值在700 g·m-2以上。

图1 2001—2015年全国植被平均NPP分布图Fig. 1 Spatial distribution of the annual mean NPP in China from 2001 to 2015

表1 2001－2015年全国九大区域的面积及植被NPP总量占全国的比例Table 1 The total NPP and area proportion of the nine regions in China in 2001-2015

2001－2015 年全国植被NPP变化范围为1.9－1131 g·m-2·a-1，全国植被 NPP 平均为 376 g·m-2·a-1。从气候分区来看，每个区植被 NPP均值、面积和NPP总量占全国的比例不同（表1）：青藏区面积和植被NPP总量占全国的比例最大，均为23%，植被NPP 均值为 385.47 g·m-2·a-1；甘新区植被 NPP 均值最小，为164.98 g·m-2·a-1，面积占全国的23%，年均NPP总量只占全国的10%；西南区植被NPP均值最大为 531.81 g·m-2·a-1，面积占全国的 11%，年均NPP总量却占全国的16%；东北区植被NPP均值为 513.10 g·m-2·a-1，面积占比 10%，年均 NPP 总量占比13%；长江中下游植被NPP均值为498.12 g·m-2·a-1，面积占比 11%，年均 NPP 总量占比 14%。其他地区如黄土高原区（459.89 g·m-2·a-1）、内蒙及长 城 区 （ 337.17 g·m-2·a-1）、 华 南 区 （ 480.49 g·m-2·a-1）、黄淮海区（441.54 g·m-2·a-1）面积和年均NPP总量占全国的比例相当。

2.1.2 近15年全国植被NPP的时间变化

对中国2001－2015年年均NPP值逐像元进行趋势变化分析（图 2）。统计发现，69.25%的像元NPP 呈增加趋势，其中最大趋势为 72.64 g·m-2·a-1，以青藏高原局部、海南东部和台湾北部地区增幅最大，全国 NPP平均变化趋势为 2.29 g·m-2·a-1，有84%的像元 NPP 年际变化范围为-3－5 g·m-2·a-1。

由图3可知，近15年全国植被NPP总体上呈逐年增加趋势，变化范围在3314.101－3847.70 Tg之间，植被平均NPP值为3519.59 Tg。从NPP增加趋势看，近 15年全国 NPP的斜率约为 23.41 Tg·a-1，即每年全国增加的NPP约为23.41 Tg。相关系数为0.392，未达显著水平。

近15年，由于气温逐渐增加，降水分布不均，全国区域植被 NPP总量变化波动明显，除西南区（-2.48 Tg·a-1）和内蒙及长城区（-0.15 Tg·a-1）呈减少趋势外，其余各区都呈增加趋势。2015年全国植被NPP最大，比最小年份（2002年）多736 Tg。占全国植被 NPP总量比例最大的青藏区总体波动较大，2002年和2006年都有小幅的减少，总体变化范围在710.41－928.52 Tg之间。西南区和长江中下游区变化范围和总体走势较相似，均在2008年、2010年和2014年出现了极小值。甘新区和内蒙及长城区则在2008年和2013年出现了极大值。

图2 2001－2015年全国植被NPP总量线性斜率分布图Fig. 2 Spatial distribution of the annual total NPP linear trend in China from 2001 to 2015

图3 2001－2015年全国九大区域植被NPP总量变化Fig. 3 Change pattern of annual total NPP in the nine regions in China from 2001 to 2015

2.2 植被NPP与气候的相关性

图4 2001－2015年全国植被NPP与年均温的相关系数分布（a）及其t检验结果（b）Fig. 4 Correlation relationship between average NPP and temperature in China during the years from 2001 to 2015

植被净第一生产力受各种自然因素（气温、降水、太阳辐射、土壤类型和地形等）和土地利用方式改变等人类活动因素的影响。为了说明气候对全国植被NPP的影响，对2001－2015年全国植被NPP与 15年来的气温和降水进行相关性分析。同时查阅相关系数的显著性资料得到，当n=15时，r的绝对值大于0.59为极显著相关，大于0.44为显著相关，其余为一般相关，由此生成相关系数空间分布图（图4a，图5a）。并对其进行t检验分析，在α=0.05（v=13，临界值为1.77）水平上得到t统计变量分布图（图4b，图5b）。

图5 2001－2015年全国植被NPP与年降雨量的相关系数分布（a）及其t检验结果（b）Fig. 5 Correlation relationship between average NPP and precipitation in China during the years from 2001 to 2015

近 15年来，全国只有小部分地区的植被年平均 NPP与年平均温度之间呈正相关，比例为32.18%，其中8.71%为显著或极显著正相关。大部分地区为负相关，其中62.39%为一般负相关。而相关性在α=0.05置信水平上，只有22.18%通过了检验，这说明 2001－2015年全国大部分地区年平均NPP与年平均温度相关性不显著。

近 15年来，全国绝大部分地区（97.15%）的植被年平均 NPP与降雨量之间呈正相关，其中71.23%为显著或极显著，这说明降雨对植被 NPP的影响在全国范围内几乎都为正效应，而且降雨对植被NPP影响的效应大于温度的影响。而NPP与年降雨量之间的相关性在α=0.05水平上有79.26%通过了检验，说明全国大部分地区的植被年平均NPP与年降雨量之间存在显著的正相关关系。

3 讨论

3.1 不同研究的NPP模拟结果

自上世纪90年代开始，大量中国学者应用模型（遥感模型和生态系统过程模型）方法对中国陆地生态系统植被第一性生产力进行了研究。由于所研究的时期、利用的模型和数据空间分辨率不同，得到的NPP模拟结果也不相同（表2）。除BIOME-BGC模型的模拟结果为1.5 Pg·a-1（以C计，下同），其余模型的模拟结果在2.138－4.65 Pg·a-1之间。由表2可知，即使采用同一个模型、同种分辨率对全国陆地生态系统植被NPP进行模拟，如CASA遥感估算模型，研究结果也不相同（高志强等，2008；陈福军等，2011）。有研究表明，中国陆地生态系统NPP 平均为(2.828±0.827) Pg·a-1（高艳妮等，2012），本研究结果与黄玫（2006）的研究结果接近，低于仲晓春等（2016）的估算值，一方面原因是研究的年份不同，仲晓春等人研究的是2001－2010年，黄玫的则是1961－2000年，另一方面，仲晓春等人采用的是MODIS NDVI数据结合CASA模型计量中国植被的NPP，黄玫则是采用耦合了土壤碳氮动力学子模块的AVIM框架的模型作为基础，模拟植被NPP的分布。研究表明，与过程模型相比，遥感产品会高估植被NPP，而过程模型机理清晰，框架完善，但参数众多，参数本地化比较困难（高艳妮等，2012）。本研究利用HASM高精度曲面建模方法，结合陆地生态系统动态植被过程模型和 MOD17A3 NPP遥感数据产品估算中国陆地生态系统植被的NPP，研究结果具有一定的可信度。

3.2 中国陆地生态系统植被NPP空间分布地域性特征及其与气候的关系

中国陆地生态系统植被 NPP空间分布地域性特征明显，总体上沿水热梯度从东向西、从南向北递减，NPP最大值出现在海南和四川盆地地区，新疆和青藏高原则最低。九大分区中甘新区地处西北内陆，气候干燥，年降水量较小，植被以落叶松和云杉为主，生产力较低，面积占比大于年均NPP总量占比；而西南、东北和长江中下游地区处于亚热带湿润、温带半湿润季风气候区，植被丰富，生产力较高。近年来，随着全球气温的不断升高，全国植被NPP总体上在波动中呈逐年增加的趋势，尤以青藏区增幅最明显，而云贵高原和东北平原植被NPP呈逐年减小的趋势，这与谷晓平等（2007）的研究结果一致，主要原因是由于云贵高原降水减少所致，也有研究认为是太阳辐射减少的原因（Wang et al.，2017）。

表2 各模型参数说明及模拟结果Table 2 Simulation results and parameters used in different model

中国陆地植被NPP与气候的相关性分析表明，气温与植被NPP的相关性不显著，只有在青藏高原地区年均温和植被NPP呈显著正相关，其他地区相关性不明显；主要原因是青藏高原属高海拔地区，受人类活动的影响较小，高原腹地年平均温度在0 ℃以下，大片地区最暖月平均温度也不足10 ℃，所以受温度的影响较大。降水是影响陆地生态系统NPP的主要原因，除云贵高原和四川盆地部分地区外，年降水量与植被NPP呈显著正相关，这与前人研究结果（朱文泉等，2007；仲晓春等，2016）略有不同。主要原因是云贵高原和四川盆地地区河流资源丰富，植被生长不受水分胁迫的影响，因此与降水的相关关系不显著。另外，近年来随着气候变暖活动的加剧，降水分布不均且逐年减少，植被受水分较温度的影响偏大。

中国陆地生态系统植被 NPP值与气温降水虽然有一定的相关关系，但受到多种因素共同作用的影响，还与人类活动、土地利用方式、管理措施（陶波等，2003；Mohamed et al.，2004；Matsushita et al.，2004）等因素有关，因此并未表现出显著相关关系。

3.3 模型实际应用中存在的不足

不同研究的NPP模拟结果差异较大，现有模型大多数采用简化的方程和经验公式的形式如光能利用率模型（CASA等），而一些基于生态系统过程的模型（BIOME-BGC等）只考虑垂直方向的通量，空间相互独立，没有考虑物质的水平运移，模型由点向面的扩展应用还存在大量的问题，导致结果差异性较大（Yuan et al.，2015；Liang et al.，2015），另外模型的结构、参数的不确定性和输入数据的精度都会影响模型的不确定性分析（Wu et al.，2006）。LPJ-GUESS动态植被模型充分考虑了植被的生长动态，可用于区域和全球尺度的生态系统过程模拟及其时空分析，但在模型实际应用中还存在一些不足，例如模型缺乏对生态系统碳氮耦合过程的深入描述，没有考虑人类活动如土地利用变化对植被生长过程的影响。

另外，随着温室气体浓度的升高，气候变化加剧，模型模拟植被净初级生产力时只考虑植被对气候变化的响应，而忽略了植被自身的适应性以及对气候的反馈，大多数情况下植被对气候的适应并不是简单的线性问题，而是复杂的非线性问题（Yang et al.，2017）。

4 结论

利用2001－2015年气候数据，运用LPJ-GUESS模型逐像元模拟了中国植被的 NPP，结合MOD17A3 NPP数据，基于HASM高精度曲面建模方法，根据气候要素和土壤质地九大分区，分析NPP时空分布特征及其与气候因子间的相关性关系，主要结论如下：

（1）近15年全国植被平均NPP空间分布地域性明显，总体上呈由西北向东南逐渐增加的趋势。其中西北沙漠地区和青藏高原腹地的植被 NPP值最小，大都在100 g·m-2以下。从西北内陆向东南沿海过渡带植被NPP值在300－600 g·m-2之间。四川盆地大部分地区以及海南省和东北平原局部植被NPP 值在 700 g·m-2以上。

（2）近15年全国植被每个像元的NPP变化范围为 1.9－1131 g·m-2·a-1，全国平均植被 NPP 为 376 g·m-2·a-1。青藏区面积和植被NPP总量占全国的比例最大，均为 23%，植被 NPP均值为 385.47 g·m-2·a-1；甘新区植被 NPP 均值最小，为 164.98 g·m-2·a-1，面积占全国的23%，年均NPP总量只占全国的10%；西南区、东北区和长江中下游区植被丰富，年均NPP总量占比大于面积占比。其他地区面积和年均NPP总量占全国的比例相当。

（3）近15年全国有65.29%的区域植被NPP呈增加趋势，最大值为 72.64 g·m-2·a-1，平均值为 2.29 g·m-2·a-1，以青藏部分、黄土高原区、海南大部和台湾北部地区增幅最大。近 15年除云贵高原和四川盆地部分地区外，有 79.26%的地区年平均 NPP与年平均降水之间呈显著正相关，而年平均NPP与年平均温度的这一比例仅为8.71%，分布在青藏高原地区。

（4）利用HASM高精度曲面建模方法，结合陆地生态系统动态植被过程 LPJ-GUESS模型和MOD17A3 NPP遥感数据产品估算中国陆地生态系统植被的NPP，研究结果具有一定的可信度。