张子凡，熊茂秋，李福杰，刘晓煌，郝玉恒，邢莉圆，王新华，赖 明，袁鹏程
（1.中国地质大学(武汉)地理与信息工程学院, 湖北 武汉 430074；
2.自然资源要素耦合过程与效应重点实验室, 北京 100055；
3.中国地质大学(武汉)地理与资源学院, 湖北 武汉 430074；
4.中国地质调查局乌鲁木齐自然资源综合调查中心, 新疆 乌鲁木齐 830000；
5.中国地质调查局自然资源综合调查指挥中心, 北京 100055；
6.中国地质调查局烟台海岸带地质调查中心, 山东 烟台 264000）
内蒙古地区是我国主要的草地资源分布区，也是我国的畜牧业区和粮食的主要生产地[1]。植被净初级生产力(net primary productivity, NPP)是指绿色植物在单位时间、单位面积内所累积的有机干物质总量[2]。NPP 是生态系统碳汇、碳源的主要参数，是陆地生态系统循环中必不可少的内容，也是评估生态系统质量与健康状况的重要依据[3]。目前研究植被净初级生产力的手段有实测法和模型估算法[4]，实测法主要通过抽样测量实地单位面积和时间内的生物量进行NPP 的调查和计算。随着遥感技术的发展，利用遥感观测数据估算区域内，长时间序列的NPP 值时空变化成为主流的研究趋势[5]。基于遥感数据结合模型估计NPP 的方法主要有CASA模型、BEPS 模型和BIOME-BGC 模型等[6-8]。通过长时间的遥感数据，结合统计学上的方法，对区域内NPP 值的变化趋势和变化可持续性研究已有较多[9-12]。研究NPP 长时间序列时空变化及其影响因子，以及影响因子的时空分异性，对研究环境的变化以及不同区域碳汇的影响因素有重要意义。
目前，利用区域的遥感观测数据，估算NPP 值，结合气候、环境等要素，研究NPP 值变化对气候和人类活动的响应对区域内生态环境监测的预测有重要意义，王耠熠等[13]通过计算NPP 栅格像元点在长时间序列上的变化，利用相关性和偏相关性分析方法，研究了若尔盖高原区NPP 值与降水量和温度之间的关系，同时重点研究了最高温和最低温对NPP 的 响 应。Yin 等[14]利 用Google Earth Engine 平台研究了森林NPP 对植被类型和气候的响应，得出降水对NPP 的影响较强，郭睿妍等[15]利用GEE 遥感数据云平台，采取岭回归分析和冗余分析的方式，研究了黄河流域森林NPP 变化及其影响因子。闫妍等[16]依据研究区内不同流域之间影响因素的不同以及CASA 模型的性质，通过从流域和省域尺度，来构建NPP 值与气候相关要素的多元回归模型，从时间上以月尺度来研究得出湖南省地区温度和水汽压对NPP 的影响较大，而降水相对较小，表明了不同区域气候因子对NPP 值的影响不同。随着机器学习技术的发展，在统计学的角度上，根据实测生物量数据和实测NPP 值结合各个因子等数据采取机器学习算法中的各类训练的方法，对区域内NPP 和生物量进行回归预测，可以以决策树的方式进行非线性预测，同时可以测定各因子数据的贡献率[17-18]。然而这些研究NPP 变化的影响因子大多是以区域每个栅格点的长时间序列的NPP 和气候数据，使用统计学的方法，进行相关性，偏相关性，以及整个区域的最小二乘法、岭回归模型和机器学习中的决策树算法进行回归分析。进行时间序列分析时，容易忽略对于时间变化不显著但地理变化显著的因子，比如地形等因子。所以只通过时间序列数据变化进行回归分析时存在局限性，而以整个区域的栅格点为对象进行全局回归分析，会存在空间局部不平稳性的问题，无法反映出NPP 分布和影响因素的空间异质性。由于位于地理接近的区域，物体空间分布存在一定的自相关性和聚集性，内蒙古NPP 区域自相关性显著[19]，研究区域内NPP 时空变化的驱动因子时，进行区域全局分析多元线性回归OLS 模型以及采用决策树模型进行分线性回归分析存在局部聚集等局部不平稳性问题。而时空地理加权回归模型，通过核函数确定时空距离的权重矩阵，可以较好地反映影响因子时空非平稳性的问题，可以从时间和空间的角度研究驱动因素[20-21]，这样可以分析不同区域内各因子不同的回归系数。
本文研究基于2000-2019 年内蒙古草原自然资源大区的MOD17A3-NPP 数据、landsat7,8 的NDVI数据和国家气象网站的气象站点数据。采用统计分析的方法，研究该区域的NPP 时空变化及其可持续性。利用时空地理加权回归模型，研究不同时空尺度下NPP 值的驱动因子，得出不同时空尺度的回归系数，量化不同区域NPP 变化的影响因素，同时也可为NPP 的预测提供帮助。
研究区域为全国自然资源综合区划结果中的华北草耕自然资源大区的北部区域[22-23]，命名为内蒙古草原资源大区，该大区是我国草地资源分布的主要区域(37°37′12″～51°02′42″ N, 106°28′40″～121°28′03″ E)，主要包括河套平原和内蒙古高原、黄土高原的部分地区，位于我国内蒙古自治区，面积约6.51 × 108hm2，主要的自然资源为草原，约占63.42% (图1)。
图1 研究区概况图Figure 1 Overview of the study area
2.1 数据源以及预处理
本研究主要使用了NPP 数据、气象数据、NDVI数据、高程数据和土地利用数据。NPP 数据来源于2000-2019 年的MODIS 陆地4 级标准数据NPP 产品数据集MOD17A3-NPP 数据，来自美国航空航天局网站(https://ladsweb.modaps.eosdis.nasa.gov/)，其空间分辨率为500 m，坐标系为WGS_1984 坐标系，投影为UTM 投影。气象数据来自于中国气象数据网(http://data.cma.cn)，选取了内蒙古、陕西、山西和河北省123 个气象站点2000-2019 年的降水量、平均气温、相对湿度、光照时长数据，采用克里金插值法进行插值，得到与NPP 数据投影坐标系相同和分辨率相同的年尺度降水、平均气温、相对湿度、光照时长的栅格数据。NDVI 数据来自美国陆地卫星landsat7,8 数 据 的 反 演(https://landsat.gsfc.nasa.gov/)。高程数据和土地利用数据来自中科院资源环境数据中心(https://www.resdc.cn/)，选取1 km 分辨率的土地数据和dem 数据。
2.2 研究方法
2.2.1 Theil-Sen Median 方法和Mann-Kendall 显著性检验
通过Theil-Sen Median 方法植被净初级生产力的趋势变化，利用Mann-Kendall 方法进行变化的显著性检验[24]。该方法可以减少外部带来的误差，同时对趋势的判断受异常值的影响较小，有非常高的适用性，可以更准确反映长时间序列下植被净初级生产力的变化。通过这个方法研究区域内每一个像元的时间序列变化，确定2000-2019 年NPP 的时空变化特征。Theil-Sen Median 方法公式如下：
式中：Z为NPP变化的显著性，n为研究的时间序列，tp为每年重复的NPP值重复的个数。定义α为显著性检验水平，Z(1-α)/2为标准正态方差，当|Z| >Z(1-α)/2时，表示时间序列ET 在α水平上存在显著变化。
2.2.2 Hurst 指数
Hurst 指数是用来评估时间序列信息长时间依赖的方法[25]。在本研究中，用Hurst 指数来研究植被净初级生产力变化的可持续性，可以判断未来变化的趋势。具体公式和计算步骤如下：
定义研究区间NPP 变化的时间序列：
对Hurst 指数进拟合：
Hurst 指数(H)在0.5 ＜H＜ 1 时，表明时间序列是一个可持续序列，即未来变化趋势与现在变化趋势一致，H越接近1，持续性变化越强。H= 0.5 则说明NPP时间序列为随机序列，不具备时间相关性，0 ＜H＜ 0.5 则表明未来变化趋势与现在有相反的趋势，H越接近0，反持续性变化越强。
2.2.3 时空地理加权回归(GTWR)模型
基于最小二乘法的多元线性回归模型可反映全局的变量因子与因变量的关系，但对于局部聚集性比较强的数据，利用全局分析无法反映因子的局部不平稳特征。地理加权回归模型可以反映局部空间的不平稳性，时空地理加权回归模型是在地理加权回归模型基础上，引入了时间维度[26]，与地理加权回归模型相比，能更好地揭示时空非平稳性变化规律，不仅能够反映空间变化趋势，还可以提高对时间数据演变规律的认知，从而增强对时空数据的分析能力[27]。本文选择2000-2019 年的数据，研究NPP值与降水、气温、湿润度、光照时长、高程、植被覆盖度(NDVI)的相互关系，计算回归系数，同时也能根据回归系数依据气候NDVI 等值对未来NPP 值进行预测。时空地理加权回归模型公式如下：
3.1 植被净初级生产力的时空变化分析
统计2000-2019 年各像元点NPP 均值的空间分布特征(图2)，并计算各年份区域内NPP 的均值，2000-2019 年区域NPP 平均值约为208.96 g·(m2·a)-1，NPP 值范围为0～727.0 g·(m2·a)-1。NPP 较高的区域主要位于东部和南部，地表资源主要为森林资源和耕地资源，较低的区域主要位于西北部的荒漠地区(图1、图2)。从时间序列上看，2000-2019 年，内蒙古草原自然资源大区NPP 均值波动较为明显，但整体上呈现较快的上升趋势，平均每年上升约4.06 g·(m2·a)-1，NPP 均 值 在2001 年 最 低，仅 约150.55 g·(m2·a)-1，在2019 年 最 高，NPP 均 值 约 为247.79 g·(m2·a)-1。在2000-2019 年，2001-2002、2007-2008、2011-2012 和2017-2018 年，区域内NPP 均值急剧上升，在2000-2001、2006-2007 和2014-2015年下降较为明显，其余年份变化较为平缓(图3)。
二是松花江黑龙江洪水量级大，持续时间长。松花江发生1998年以来最大流域性较大洪水，其中嫩江上游发生超50年一遇特大洪水，松花江上游发生超20年一遇大洪水。黑龙江发生1984年以来最大流域性大洪水，下游发生超100年一遇特大洪水。嫩江、松花江干流水位超警戒历时46天，黑龙江干流水位超警戒历时58天。
图3 2000-2019 年内蒙古草原自然资源大区植被净初级生产力均值及变化趋势图Figure 3 Mean net primary productivity (NPP) values and their change trend for natural resources in the grassland region of Inner Mongolia from 2000 to 2019
基于2000-2019 年内蒙古草原自然资源大区的NPP 栅格数据，采用Theil-sen 斜率趋势分析方法和Mann-Kendall 显 著 性 检 验 方 法，取|Z| > 1.645、1.96 和2.576，表示研究序列通过了置信度为90%、95%和99%的显著性检验，分别对应为弱显著、显著和极显著。得到NPP 值变化趋势以及变化显著性的空间分布图(图4a)，大部分区域NPP 值呈上升趋势，约占整个区域的94.01%，仅5.99%的区域呈现减少趋势。并且大部分区域NPP 值增加趋势表现为极显著增加，约占总研究区面积的60.31%，其次为显著增加，约占研究区域的19.57%。NPP 呈减少趋势的区域零星分布，其中不显著减少所占的区域最多，约5.88%，极显著、弱显著和显著减少的区域仅占约0.15%，零星分布于建设用地相关的区域，仅有约1.01 × 105hm2。
图4 2000-2019 年植被净初级生产力值变化趋势和变化显著性(a)以及可持续性(b)Figure 4 Net primary productivity change trend and significance (a) and persistence of the change trend (b) from 2000 to 2019
将NPP 的变化趋势与Hurst 指数相叠加得出NPP 值的变化趋势和持续形式，分为持续减少(2000-2019 年呈现减少趋势，H> 0.5)、潜在增加(2000-2019 年呈减少趋势，H＜ 0.5)、无变化趋势(H= 0.5)、潜在减少(2000-2019 年呈现增加趋势，H＜ 0.5)、持续增加(2000-2019 年呈现增加趋势，H> 0.5)。由图4b 可知，在2000-2019 年，尽管内蒙古草原资源大区94%的区域呈上升趋势，但未来潜在减少的区域最多，约占整个区域面积的74.15%，分布于该区域西部，南部的大部分地区，其次为未来持续增长的区域，约占整个区域的19.45%，主要分布于该大区的东部和东北部的有林地地区。处于不确定游离状态的区域仅有约5.83%，主要分布在荒漠和湿地地区。未来NPP 呈现持续减少的区域仅有约0.23%，零星分布于该大区的中部。
3.2 NPP 值时空变化的驱动因子分析
通过时空地理加权回归模型，选取2000、2005、2010、2015、2019 年的降水、年平均气温、光照时长、相对湿度等气候因子以及海拔高度地形因子和NDVI 植被覆盖度因子。通过Arcgis 渔网工具全局选取了1 056 个采样点，研究区域内NPP 值时空变化与各因子之间的关系，然后进行点的插值得出植被净初级生产力值变化和各影响要素回归系数分布(图5、图6)。最终通过时空地理加权回归模型(GTWR)得出的模拟拟合度R2= 0.925，该结果表明回归精度较好，可以进行分析研究。
图6 2019 年植被净初级生产力值变化、各影响要素回归系数分布Figure 6 Changes in NPP of vegetation in 2019, distribution of regression coefficients of various influencing factors
统计区域2000-2019 各影响因子得到平均回归系数，2000-2019 年降水、相对湿度和NDVI 变化和NPP 值相关性较强，主要为正向相关关系，气温、光照时长和蒸散发主要为负向相关关系(表1)。在2000 年，降水对NPP 值变化的回归系数最高值为1.113，平均值为0.368，呈负相关的区域约占5%，在空间分布上，主要分布在西南少部分地区，其他区域呈正相关，影响程度从西南到东北存在明显的空间分异性，从西南到东北地区呈现线性增加的趋势；
2000 年，气温对NPP 影响主要为负相关，回归系数的平均值为-0.027，最高值为0.803，最低值为-0.631，其中正向相关影响的主要分布于该区锡林郭勒盟市的西部，占整个研究区域的20%，其他区域呈负向相关关系，并且存在明显的空间异质性，在研究区域西南影响程度较大，东北部次之，中部最小；
光照时长对NPP 值的影响相对较小，在2000 年影响平均回归系数为-0.061，最高值为-0.051，最低值为-0.066，几乎所有区域呈负相关性影响，影响程度从西南到东北有显著的线性减小趋势；
湿润度对NPP 值的影响在2000 年影响平均回归系数为0.228，最高值为0.233，最低值为0.221，几乎所有区域呈现正向相关，影响程度从西南到东北呈线性增加的趋势；
2000 年蒸散发对NPP 的影响平均回归系数为-0.020，最高值为-0.030，影响程度由西南到东北区域逐渐递减；
NDVI 值与NPP 值主要是正向相关关系，且回归系数较高，其中在西南地区最高，由西南向东北呈递减趋势；
2000 年海拔高度与NPP 值回归系数的平均值为0.080，最低值为-0.032，最高值为0.177，其中海拔高度与NPP 呈现负相关性的区域占37.45%，主要是在该大区的西南部，其他区域呈现正相关，影响程度从该区西南到东北有减小的趋势(图5)。
表1 2000-2019 年植被净初级生产力影响因子回归系数Table 1 Regression coefficients for net primary productivity influencing factors, 2000-2019
图5 2000 年植被净初级生产力值变化、各影响要素回归系数分布Figure 5 Changes in NPP of vegetation in 2000, distribution of regression coefficients of various influencing factors
2019 年，降水、温度、光照时长、湿度、蒸散发等气候因子对NPP 值影响的回归系数的平均值分别为0.257、-0.145、-0.201、0.081 和-0.001，海拔和NDVI 值的回归系数分别为0.078 和0.559。与2000相比，降水、湿度和蒸散的影响程度在减少，光照、温度和NDVI 值的影响程度有所增加，海拔影响基本不变(表1)。2019 年，降水对NPP 值的影响的回归系数最高值为0.715，最低值为-0.158，其中呈负相关的区域，约占10%，主要在呼伦贝尔市西北和巴彦淖尔市东南，呈正相关的区域影响程度由西南向东北呈增加趋势；
气温对NPP 的影响，最高值为0.901，最低值为-1.109，大部分区域呈现负相关影响，约占59.56%，呈正相关的区域主要在该大区东北部和巴彦淖尔市西部；
2019 年，光照时长对NPP 的影响回归系数最高值为-0.197 4，最低值为-0.202 4，平均值为-0.201，所有区域光照时长与NPP 值呈现负相关影响，影响程度从西南到东北呈线性减弱；
相对湿度对NPP 值的回归系数最高为0.077，最低为0.084，影响程度较小，全部区域湿度与NPP 值呈正相关影响，与2000 年不同，影响程度由西南到东北逐渐减弱；
2019 年地表蒸散对NPP 值影响的回归系数最高值为-0.08，最低值为-0.01，影响程度变化不大，由西南到东北递减；
2019 年NDVI值和NPP 值的相关系数最低值为0.08，平均值为0.559，所有区域呈现正相关关系，其中在东北部和西南部的部分地区，NDVI 值对NPP 的影响较小，影响程度存在显著的聚集性；
海拔高度对NPP 值的影响与2000 有所不同，最低值为0.048，最高值为0.118，全区域呈正向关，有显著的空间异质性，由西南向东北递增(图6)。
本研究通过对2000-2019 年NPP 值的变化及其驱动因子研究得出，降水量变化是内蒙古草原自然资源大区NPP 值增加的主要气候因素，NDVI 值与NPP 值密切相关，但NDVI 值越高的地区，NPP受NDVI 变化的影响相对较小。另外内蒙古草原自然资源大区大部分区域温度和光照时长的变化与NPP 值呈负相关关系，同时也有学者研究指出一定程度的气温升高可以提高光合效率和资源利用效率从而增加NPP[29-30]，可能因为在内蒙古半干旱地区，由于水分较少，温度升高会引起水分的减少，从而使植被利用水分减少反而降低了光合效率，地表水分是影响NPP 变化的主要原因。越干旱的地区温度、光照时长和蒸散发对NPP 值的负相关性就越强(图5、图6)。而NPP 与温度呈正相关的区域主要在东北部和锡林郭盟市西部植被覆盖指数较高的地方。各类要素对NPP 值的影响从该大区西南到东北有显著的线性变化，有明显的空间异质性，可以为区域碳储量和自然资源的分区研究提供帮助。
目前本研究还存在许多不足，主要对NPP 值时空变化的趋势进行分析。以及利用气候、海拔、植被覆盖度等因子对NPP 值进行时空地理加权回归分析，通过回归系数研究各因子对NPP 值的影响。主要还是以统计学以及统计学相关的模型进行分析，没有深入研究内部变化机理，所得到的影响系数只是统计学上的相关性。由于2000 年和2019 年的气温和植被类型存在较大差异，有些区域的NPP值的影响回归系数在2000 年与2019 年相差较大，甚至相关趋势都不同，这可能是不同气候对NPP 的影响存在一个阈值，这个阈值需要大量的时空数据分析，可能是以后需要研究的方向。另外NPP 的变化可能与人类的生产活动有很大相关性。综合区域自然因子和工农业、畜牧业、人口密度、GDP 等社会经济因子研究NPP 的变化是今后需要解决的问题。
在2000-2019 年间内蒙古草原自然资源区域NPP 均值呈波动上升趋势，平均每年增加约207.72 g·(m2·a)-1。在2000-2019 年94.01%的区域NPP 值呈上升趋势，79.88%的区域NPP 值呈显著上升(|Z| >1.96)，呈显著下降的区域极少，仅有约0.15%；
利用Hurst 指数来判断每个像元点变化的可持续性，得出尽管2000-2019 年有94.01%的区域NPP 值呈现上升趋势，但在未来会持续增加的区域仅占整个研究区域的17.45%主要分布在研究区域东部和东北部地区，74.15%的区域在未来有潜在下降的趋势。
采用时空地理加权回归模型，研究NPP 值变化的影响因子。得出在2000-2019 年，该研究区域内降水、相对湿度和NDVI 值的影响主要呈正相关，降水和NDVI 值对NPP 的影响较高，在2000 和2019年相关系数分别为0.368、0.257 和0.402、0.259。温度、蒸散发和光照时间主要呈负相关关系。其中温度和光照时间影响相对较大，2000 年和2019 年回归系数平均值分别为-0.027、-0.145 和-0.061、-0.201。蒸散发影响较小，2000 年和2019 年回归系数的平均值分别为-0.02 和-0.001。
各因子对NPP 的影响存在空间异质性。降水与NPP 值的变化主要呈正相关关系，但在2000 年和2019 年仍然有一部分区域呈现负相关，分别占5%和10%，主要在该大区西南和东北的部分区域。温度主要是负相关，2000 和2019 呈正相关的区域主要在锡林郭盟西部；
另外光照时长、相对湿度、蒸散量和高程对NPP 值的影响在不同区域也有所不同，从该区西南到东北有显著的线性变化趋势。