2000~2020年黑龙江省陆地植被NPP驱动因子分析

程春香 ，于 敏，吴门新，薄 宇，毛子军

（1.黑龙江省生态气象中心，哈尔滨 150030；2.东北林业大学森林植物生态学教育部重点实验室，哈尔滨 150040；3.国家气象中心，北京 100081）

人类活动与气候变化影响生态环境[1]。为应对全球环境问题我国提出“2030年实现碳达峰、2060年实现碳中和”的战略目标[2]。净初级生产力（Net primary productivity，NPP）是陆地生态系统生产力重要指标之一，作为植被净固碳量，用于表征植被生长抵消人为碳排放量的潜力[3]，对维持全球碳平衡、减缓气候变化及完成双碳战略目标均具有重要意义。

NPP时空演变及驱动因子研究是评估陆地生态系统固碳能力基础，近年广受学者关注。程春晓等研究表明2001~2010年中国东北不同地表覆盖类型NPP变化趋势差异显著[4]。刘晓光等发现黑龙江省三江平原北部NPP 整体呈线性增加趋势，先缓慢减少后波动增加[5]。Tripathi等利用时滞相关分析法研究气候变量对印度热带落叶林净初级生产力的贡献，结果显示降水是影响生态系统净初级生产力主要变量[6]。刘旻霞等发现气温、降水和土地利用/覆被变化（LUCC）等是影响2000~2014 年青海省NPP重要因子[7]。以往研究大多集中于NPP时空演变分析，黑龙江省驱动因子研究十分有限。黑龙江省作为生态资源大省，是我国生态保护修复工程中重点保护修复区[8-9]。

地理探测器是探测数据分异性，揭示数据背后驱动力的一种新兴方法。可定量分析各影响因子权重，与回归分析等传统统计方法相比，更客观、科学，且具有较高精度[10]。在全球变化和我国双碳战略背景下，本文利用已验证TEC模型[11]计算植被NPP，采用Sen趋势度估计、土地利用动态度等方法，针对植被资源丰富的黑龙江省，开展2000~2020年生态保护修复关键实施期植被NPP及其影响因子时空变化特征研究，并利用地理探测器和相关分析法开展NPP驱动因子定量化分析，对生态保护修复以及生态系统碳循环演变理论研究均具有重大意义，为我国未来生态环境管理和政策制定以及“碳达峰、碳中和”目标实现提供科学参考。

1 研究区域概况与研究方法

1.1 黑龙江省概况

黑龙江省区域辽阔，地形与气候条件复杂，气候变化较为显著，自然资源丰富，森林、农田、草地、湿地、聚落等是黑龙江省典型陆地生态系统中重要组成。黑龙江省森林覆盖率位列全国首位，省内林区主要分布于大兴安岭、小兴安岭以及牡丹江东部；此外黑龙江省作为国家重点湿地省份，在三江平原、松嫩平原、大兴安岭、小兴安岭和东部地区均存在大片湿地区域，约占全国天然湿地总面积1/5；黑龙江省草地资源主要分布在三江平原、松嫩平原、大兴安岭和小兴安岭，黑龙江省作为农业大省，粮食种植面积大、产量高，是我国重要商品粮生产基地[12]。黑龙江省在全球区域碳循环过程中发挥重要作用，且对于气候条件变化敏感，在陆地生态系统碳循环过程中极易受气候变化影响。

1.2 数据资料来源

土地覆盖数据集来源于全球地表覆盖Globe-Land30数据集（国家基础地理信息中心全球地表覆盖数据产品服务网站DOI：10.11769），空间分辨率30 m，已发布2000 年、2010 年和2020 年数据，包括耕地、林地、草地、灌木地、湿地、水体、苔原、人造地表、裸地、冰川和永久积雪等10 个类型。根据黑龙江省实际情况和研究需要，简化合并为7个类型，分别代表森林生态系统、农田生态系统、草地生态系统、湿地生态系统、聚落生态系统、其他生态系统和水体，文中生态系统尺度分析均采用2020年GlobeLand30数据作为分区统计的基础数据。

分析气象数据来源于黑龙江省气象数据中心。数据集主要包括2000~2020年黑龙江省84个国家气象观测站点年降水量（R）、年平均气温（T）、年最高气温（Tmax）和年最低气温（Tmin）数据。根据研究需要对气象数据进行空间插值，获得空间分辨率为1 km的气象要素数据集。空间插值方法选取协同克里金法，空间插值过程中考虑海拔高度因素影响，以提高降水和气温要素插值结果准确性。

黑龙江地图底图数据来源于自然资源部网站（http://bzdt.ch.mnr.gov.cn/），本研究使用的是黑龙江省标准地图——政区简图版（省级轮廓）32 开，审图号：黑S（2018）029号。

1.3 植被NPP计算

本研究植被NPP数据采用TEC模型计算。TEC模型由延昊在2015 年根据LUE 模型改进提出，是一种新的光能利用率模型[13-14]。该模型充分考虑水分胁迫、辐射以及C3/C4植物光合作用的差异。模型产品已完成地面通量观测数据验证，精度优于BESS和MOD17产品[11]。模型公式如下：

式中，NPP、GPP、Rg和Rm分别表示陆地植被净初级生产力、总初级生产力、生长呼吸消耗量以及维持呼吸消耗量（gC·m-2）。ε*为植被最大光能利用率，C3 植物ε*=1.8 gC·MJ-1，C4 植物ε*=2.76 gC·MJ-1。Tε和Wε分别为温度胁迫系数和水分胁迫系数，FPAR为光合有效辐射吸收率，PAR为光合有效辐射（MJ·m-2），Ta为月平均气温（℃）。

1.4 趋势分析

采用一元线性回归法分析研究要素区域整体趋势，可直接计算变化率和显著性检验。分析要素像元级趋势采用Sen趋势度估计法，利用Mann-Kendall（简称MK）趋势检验法作显著性检验。两种方法相结合可降低噪声干扰，在一定程度上提高空间像元趋势准确性[15]。

Sen趋势度是经过计算序列的中值，计算公式：

式中，xj和xi为时间序列数据。如计算出β值为正数，说明研究时间序列呈增加趋势，如计算出β值为负数，说明研究时间序列呈减少趋势。

MK 法对异常值不敏感，且对序列分布无要求。检验中对于序列X=（x1，……，xn），先确定所有对偶值（xj，xk，k＞j）中xj和xk大小关系。研究假设，H0：序列数据排列随机，无显著趋势，H1：序列有上升/下降单调趋势。检验统计变量S计算如下式：

其中

S为正态分布，其均值为0，方差Var（S）按如下公式计算：

当n＞10 时，标准正态统计变量按如下公式计算：

在双边趋势检验中，在给定的0.05 置信水平上，如果|Z|≥1.96，则原假设不成立，时间序列数据在0.05置信水平上存在明显上升或下降趋势。

1.5 协同克里金插值

协同克里金插值法是在普通克里金插值法基础上实现区域化变量最佳估算的方法，具有两个以上的协同区域化属性[16]。在进行气象要素协同克里金插值时，可选择与要素具有相关性的辅助因子，如海拔高度、坡度、坡向、风向等，研究选择海拔高度作为辅助因子进行协同克里金插值。公式如下：

式中，Z(x)为插值点估计值；Zi为第i个样本点气象要素值；λ、λi为权重系数；y(x)为海拔高度；my、mz为海拔高度和气象要素全局平均值。

1.6 相关分析

针对NPP 与气象要素计算相关系数，进行相关分析，可明确不同因子相关关系[17]，相关系数R计算公式如下：

式（11）中，R表示相关系数；n表示研究数据年限；xi为第i年NPP数据；yi为第i年气象要素数据。

通过t检验进行相关性系数显著性检验，公式具体如下：

为定性分析相关关系显著性，按照阈值划分相关系数显著性水平，根据R值与P值大小，R＞0说明自变量与因变量呈负相关，反之自变量与因变量呈正相关。P＜0.01 时，说明相关关系达极显著程度；0.01＜P＜0.05，说明相关关系达显著程度；P＞0.05，说明相关关系不显著。

1.7 土地利用动态度

土地利用动态度是定量研究区域土地覆盖类型变化速率指标，反映土地覆盖类型变化转移速度和剧烈程度。通过计算土地利用动态度可确定土地覆盖类型变率最明显地区[18]。计算公式如下：

式（13）中，K为土地利用动态度，Ua和Ub分别为研究阶段初期和末期土地覆盖类型面积，T为样本时间长度。

1.8 地理探测器

地理探测器可有效分析要素之间因果关系。假设变量Y，因子X，公式如下：

式中，q为地理探测器统计值，值域为[0，1]。h为变量Y或因子X分层数或分类数；Nh和N分别为层h和全区单元数；σh2和σ2分别是层q和全区Y值方差。如果分层由变量Y生成，代表变量Y数据本身分异性，值越大数据分异性越明显；如果分层由因子X生成，则q统计值越大表示因子X对变量Y解释力越强，反之则越弱。当q=1，变量Y空间分布完全受因子X控制；当q=0，因子X与变量Y无任何关系。

本研究根据文献分析，最终选取R、T、Tmax、Tmin、土地覆盖类型（2020 年）、海拔高度、坡度、坡向8 个自然因子指标进行NPP 驱动因子分析，部分指标多年区域特征见表1。

表1 2000~2020年NPP影响因子区域统计分析Table 1 Regional statistical analysis of impact factors of NPP from 2000 to 2020

在黑龙江省区域创建100 个随机矢量样点，提取2000~2020 年黑龙江省NPP 和8个因子指标随机样点数据各2 100 个，并按照等间距法进行等级划分获得分级数据。利用NPP与NPP分级数据，运用地理探测器中分异及因子探测模块计算q统计量，分析年NPP 数据分异性。利用NPP 分别与R、T、Tmax、Tmin、土地覆盖类型（2020年）、海拔高度、坡度、坡向等8个因子分级数据，运用地理探测器中分异及因子探测计算q统计量，分析探测2000~2020 年黑龙江省NPP 驱动因子，辅助揭示气候因子、地理环境因子以及土地覆盖类型对NPP 的影响。本研究因子探测采取时空全方位随机采样法（多年、空间随机样点采集，合并分析），实现真正的影响因子探测。q值表示影响因子解释q×100%的NPP[10]。

2 结果与分析

2.1 气象要素变化趋势

2000~2020年黑龙江省年平均气温变化整体呈升高趋势，升高区域高值区主要分布于大兴安岭高海拔区、小兴安岭西北部以及三江平原北部局部，主要为森林和农田，变化率大于0.05 ℃·a-1，并在高值区部分区域，升高趋势显著，占全省面积5.2%。年平均气温变化在小兴安岭东部局部降低，变化率小于-0.02 ℃·a-1，但未达到显著水平（见图1a、1d）。2000~2020年黑龙江省年降水量变化全省均呈增加趋势，90.4%降水量达到显著水平以上。降水量为“中部高-南北低”，除大兴安岭黑龙江区域高海拔区和东南部山区的大青山附近外，增加趋势均达到显著以上水平。降水增加趋势的高值区主要为三江平原中部、小兴安岭东南部及松嫩平原北部局部地区，变化率大于14 mm·a-1（见图1b、1e）。

图1 2000~2020年黑龙江省气象要素变化趋势和地理环境Fig.1 Variation trends and geographical environment of meteorological elements in Heilongjiang Province from 2000 to 2020

2.2 土地覆盖类型变化

计算2000 年和2020 年土地覆盖类型转移矩阵，黑龙江省有550 64.8 km2土地覆盖类型发生改变，其中森林转入量和转出量分别达到18 622.5和15 452.7 km2，转入量主要来自草地和农田，有14 769.6 km2草地转森林和3 412.8 km2农田转森林。黑龙江省2000~2020 年间森林、聚落和水体面积大幅增加，分别增加3 169.8、2 287.3 和972.9 km2，土地利用动态度分别为0.1、1.3 和0.7%·a-1，其中森林增加面积最大，聚落动态变化最为明显。而草地、湿地、农田和其他面积减少，其中草地减少最为明显，减少5 597.1 km2，土地利用动态度为-0.4%·a-1；湿地减少面积位居其次，减少550.0 km2，土地利用动态度为-0.3%·a-1；农田和其他类别，减少面积分别为269.3 和13.5 km2（见表2）。

2010~2020 年与2000~2010 年，两阶段土地覆盖类型变化对比显示，森林在2000~2010年大幅增加后出现减少，农田减少后2010~2020年又大面积增加。草地在两阶段均大面积减少。与2000~2010年相比，2010~2020 年湿地面积减少速度大幅降低。聚落在2000~2010 年小幅增加，在2010~2020年期间出现大幅度增加。2000~2020年土地覆盖类型转移面积小于2000~2010 年和2010~2020 年转移面积之和，表明黑龙江省土地覆盖类型变化在部分区域内存在一定波动反复（见表2）。

表2 2000~2020年黑龙江省土地覆盖变化面积和变化率Table 2 Area and rate of land cover change in Heilongjiang Province from 2000 to 2010

2.3 NPP时空分布与趋势

2000~2020 年黑龙江省年均NPP 均值（简称NPP平均值）为602.7 gC·m-2。不同生态系统NPP平均值，森林＞草地＞农田＞湿地＞聚落，分别为668.4、586.1、560.5、523.4 和490.1 gC·m-2。森林、农田和草地是2000~2020年黑龙江省生态系统年均NPP总量（植被年均总净固碳量）前三名，分别为130.7、105.8 和34.3 TgC（1 Tg=1012g），占全省46.0%、37.2%和12.1%。湿地和聚落受生产力和分布面积限制，年均NPP 总量较小，仅为5.0 和5.3 TgC，占全省1.8%和1.9%（见表3）。

表3 2000~2020年黑龙江省年均NPP不同类型生态系统的统计分析Table 3 Statistical analysis of annual average NPP of different ecosystems in Heilongjiang Province from 2000 to 2020

2000~2020年黑龙江省年NPP均值呈极显著增加趋势（P＜0.01，R2=0.81），变化率为8.55 gC·m-2·a-1（见图2）。2010年是我国生态保护修复重大工程及解决三农农业问题的重要时间节点[8]，阶段性分析2000~2010 年（P＜0.01，R2=0.53）和2010~2020 年（P＜0.01，R2=0.64）NPP 同样呈增加趋势。2010~2020 年NPP 均值变化率（10.85 gC·m-2·a-1）高于2000~2010 年（7.69 gC·m-2·a-1）。

图2 2000~2020年黑龙江省年NPP均值年际变化Fig.2 Interannual variation of average annual NPP in Heilongjiang Province from 2000 to 2020

2000~2020年不同生态系统年NPP均值极显著增加，农田增加最快，为10.2 gC·m-2·a-1，其次是聚落、草地、森林和湿地，分别为9.3、8.3、7.1和6.4 gC·m-2·a-1。2000~2010 年和2010~2020 年不同生态系统年NPP 均值同样呈增加趋势，大部分变化率达到显著水平，但2000~2010年森林变化率和2010~2020年湿地变化率未达到显著水平，即年NPP 均值增长趋势不显著。2010~2020年不同生态系统年NPP均值变化率除湿地外均高于2000~2010年，变化差值聚落最大，为4.6 gC·m-2·a-1，农田和森林其次，为3.7和2.6 gC·m-2·a-1（见表4）。

表4 2000~2020年不同类型生态系统年NPP均值统计分析Table 4 Statistical analysis of average annual NPP of different ecosystems from 2000 to 2020

2.4 NPP驱动因子分析

NPP 分异性探测结果显示，2000~2020 年黑龙江省NPP 具有较高数据分异性，q统计量为0.89。因子探测各因子对NPP 的q统计量，R＞坡度＞海拔高度＞土地覆盖类型＞Tmax＞T＞坡向＞Tmin，降水量对NPP 影响最大，解释力为25%；地理环境因子坡度、海拔高度和土地覆盖类型对NPP 也具有重要影响，解释力分别为20%、18%和14%，各因子对NPP 影响均达到显著以上水平（见表5）。

表5 2000~2020年NPP分异及影响因子统计分析Table 5 Statistical analysis of NPP differentiation and impact factor from 2000 to 2020

2.5 NPP与气象要素相关分析

2000~2020年黑龙江省植被NPP与降水量和平均气温均呈正相关，即随降水量增加或平均气温升高，黑龙江植被NPP 也增加，相关系数分别为0.76 和0.24，其中降水量相关系数达到极显著（见表6）。

不同生态系统，2000~2020 年黑龙江省植被NPP与降水量均呈极显著正相关，各生态系统相关系数，聚落＞农田＞草地＞湿地＞森林。其中，聚落和农田相关系数较高，相关系数分别为0.88和0.85（见表6），显著以上正相关区域占聚落和农田总面积比例分别为89.4%和92.2%，表明聚落和农田的植被NPP 受降水影响较为明显。而森林可能因为根系较深、湿地因为水分相对充沛，水分限制低于聚落和农田，与降水相关系数相比略低，相关系数分别为0.56和0.67（见表6），显著以上正相关区域占森林和湿地总面积比例分别为48.1%和59.0%。平均气温与各生态系统植被NPP均呈正相关，即平均气温升高有利于植被生长，但相关均不显著，响应最大的是森林生态系统。

表6 不同生态系统类型NPP与年气温和年降水量相关系数Table 6 Correlation coefficients between ecosystem NPP and annual temperature and precipitation in different ecosystem types

3 讨 论

3.1 NPP及其影响因子变化特征

随人口增长、经济发展和技术进步，大气中温室气体浓度迅速上升，导致全球气温普遍性升高，全球降雨量和降雨频率在时间和空间上发生改变，部分区域降水显著增加[19-20]，与本研究结论一致。2000~2020年黑龙江省平均气温升高、降水量增加，气候总体呈现暖湿化趋势，热量条件和水分条件均得到明显改善，与1961~2003年黑龙江省前期研究结果一致[21]。本研究中2000~2020年间黑龙江省森林面积增加最大，高玉娟等研究发现黑龙江省森林资源出现森林面积、林木蓄积量和森林覆盖率恢复性“三增长”，多项生态保护修复工程取得阶段性成果，表明21 年间黑龙江省植树造林成效显著[22]。同时本研究还发现2000~2020年黑龙江省聚落生态系统土地利用动态度最大，面积动态增加最为明显，并且在2010 年后增速加快，说明黑龙江省21 年间城镇化进程逐渐推进，2010 年后城镇化速度明显加快。此外，本研究中2000~2020年黑龙江省陆地植被NPP呈极显著增加趋势，结论与程春晓等在东北地区研究结果一致[4]，显示出2000~2020 年黑龙江省具有持续增强的植被固碳能力。

3.2 植被NPP驱动因子分析

在全球气候变暖背景下，陆地植被与气象要素相互作用过程是现今地球科学领域关注的重点和热点。本研究中2000~2020年黑龙江省降水量对NPP解释力最大，表明在自然因子中降水是影响黑龙江省NPP主导因子。结果与Ji等在2000~2018年我国5个主要森林生态系统NPP年际变化及其驱动因子研究中发现一致[23]，研究发现森林NPP变化在自然因子中主要是由降水变化引起，其次是温度波动。Tripathi 等研究也发现降水是影响印度热带落叶林净初级生产力的主要气候变量[6]，但Cuo 等在青藏高原NPP 年代际和年际变化研究中发现，温度是影响NPP 可变性和长期变化主要影响因子[24]，与上述结论不同。结合Liebing 最小因子定律分析，推断不同区域影响NPP 的主导因子可能因区域气候条件和地理环境差异而不同。

本研究中2000~2020年黑龙江省大部分地区降水量与植被NPP 呈极显著正相关，植被NPP 对降水量比较敏感，显示黑龙江省降水相对匮乏，大部分地区植被生长受到水分限制。同时2000~2020年黑龙江省年降水量总体呈极显著增加趋势，水分条件明显改善。分析可知，黑龙江省区域水分条件相对匮乏和降水发生极显著变化，可能是降水量成为2000~2020年黑龙江省NPP主导影响因子的原因。随着全球气候变化加剧，在关注黑龙江省气候条件对NPP 正向影响的同时，还需警惕未来极端气候或未来气候条件下行对NPP带来的负向影响[25-26]。此外，因子探测显示土地覆盖类型对NPP也具有重要影响，侧面证明我国多项生态保护修复工程的有效性，并可能已通过改变土地覆盖类型对NPP产生重要影响，表明未来应继续实施适合的生态修复政策增加NPP，进一步提高陆地生态系统固碳潜力。

值得指出的是，NPP仅代表植被净固碳量，无法直接揭示陆地生态系统碳源/汇变化，有待于进一步完善NEP估算算法，并开展深入研究。