京津冀植被净初级生产力时空分布及自然驱动因子分析

王金杰，赵安周 ，胡小枫

1. 河北工程大学地球科学与工程学院，河北 邯郸 056038；2. 河北工程大学矿业与测绘工程学院，河北 邯郸 056038；3. 中国科学院地理科学与资源研究所/资源与环境信息系统国家重点实验室，北京 100101

植被净初级生产力（Net Primary Productivity，NPP）是指植被在单位时间、单位面积内通过光合作用产生的有机物总量减去自养呼吸后的残余部分，是植被与外界环境因子相互作用下的产物（Lieth et al.，1975）。该参数可以直接反映植被在自然环境下的生产能力，同时也是陆地生态系统中物质、能量循环的重要组成部分（Cramer et al.，1999；Ito，2011；陈晓玲等，2016）。因此，分析NPP的时空分布及其影响机制对全球陆地生态系统的变化及碳平衡有着重要的意义（Ruimy et al.，1994；Field et al.，1998）。

早期的 NPP研究主要基于站点观察数据和统计模型模拟（Lieth et al.，1975），这类方法易受空间尺度限制，不适用于大区域 NPP的时空演变的研究（张笑鹤，2011）。随着遥感技术的发展，利用遥感数据和模型估算开展 NPP监测和时空演变分析已经成为研究 NPP的有效手段之一（朱文泉等，2005）。遥感数据时间序列长、覆盖面积大，同时可以提供模型估算所需的植被及环境参数，为大尺度、高精度的 NPP估算研究提供了有力的支撑（赵国帅等，2011）。在众多估算模型中，生态遥感耦合模型利用遥感数据简化了估算过程，同时又对过程模型和光能利用率模型进行了改进，能大大提高NPP的估算精度，是目前估算NPP的主要方法（张锐，2015）。这其中，MOD17A3HGF NPP数据来自给定年份 8天净光合作用合成产品（MOD17A2H）的所有数据之和，是利用BIOME-BGC过程模型和 PSN光能利用率模型的耦合模型模拟得到的。同时该数据产品对每个像素清除了8天LAI/FPAR中质量较差的输入，并通过线性插值法估算其值，改进了MOD17A3数据的质量，目前已在全球不同区域的研究中得到验证（Ge et al.，2021；Venter et al.，2021；Wang et al.，2021）。

NPP与自然环境因素之间的影响机制一直是国内外学者研究的热点问题之一（彭少麟等，2000），诸多学者已从不同区域尺度对 NPP时空分布及其影响因素进行了分析（Mao et al.，2014；Liu et al.，2019；Ge et al.，2021；陈福军等，2011；朱士华等，2017）。目前京津冀地区对于 NPP驱动因素的研究主要采用趋势分析、相关系数等方法（张莎，2015；朱利欣，2018；孙涛等，2020），即多集中在单一的动态变化的角度，对静态的空间分异上的相关性方面研究较少，尤其对于多自然要素因子间交互作用的定量归因相对薄弱（贺倩等，2020）。地理探测器（Wang et al.，2010）可以定量地探测单一因子对 NPP的解释力以及因子间的交互作用，能够较好地解决上述研究中的局限性，现已广泛应用于自然环境（Chen et al.，2021；左丽媛等，2020）、人文社会（Zhang et al.，2020a）、医学（Zhang et al.，2020b）等多个领域。已有学者基于地理探测器分析了 NPP的驱动因素（李金珂等，2019；潘洪义等，2019；贺倩等，2020），但并未考虑到植被类型、土壤类型等因素的影响，使得区域尺度 NPP时空分布的驱动机制尚未完全清晰（李金珂等，2019）。

京津冀是中国北方经济最发达的地区。《京津冀协同发展规划纲要》指出生态环境保护是京津冀一体化协同发展的基础和3个率先突破之一（方创琳，2017）。针对上述研究存在的问题，本文基于MOD17A3HGF数据分析了 2000—2019年京津冀NPP时空演变特征，并利用地理探测器定量探测京津冀 NPP的空间异质性以及气候、地形、植被类型、土壤类型等自然因子的驱动机制，进而为京津冀生态环境建设与改善工作提供科学依据。

1 材料与方法

1.1 研究区概况

京津冀地处 113°27′—119°50′E，36°03′—42°40′N，位于华北平原北部，包括北京、天津、河北省，总面积约为2.18×105km2。该地区东临渤海，西倚太行，北靠燕山山脉，南部多平原，地势自西北向东南倾斜，形成了十分明显的垂直地带特征（刘广，2008）。囊括了平原、台地、丘陵、山地以及高原等多种地貌类型，是中国地貌的缩影（吕国旭等，2017）。气候类型属温带大陆性季风气候，日照充沛，寒暑温差较大，年均气温 10—11 ℃，年降水量为400—800 mm，多集中于夏季。土壤类型丰富，以半淋溶土、半水成土为主。植被类型主要包括针叶林、阔叶林、草地、灌丛以及耕地，其中又以耕地、灌丛为主。

1.2 数据来源与处理

2000—2019年NPP数据来源于美国国家航空航天局（National Aeronautics and Space Administration，NASA）的MOD17A3HGF产品，时间分辨率为 1 a，空间分辨率为 500 m。DEM（Digital Elevation Model）数据来源于地理空间数据云（http://www.gscloud.cn/），空间分辨率为90 m，用于提取研究区的海拔、坡度、坡向等数据。植被类型数据与土壤类型数据分别来源于中国科学院资源环境科学数据中心（http://www.resdc.cn/）的《1꞉1000000中国植被图集》与《1꞉100万中华人民共和国土壤图》（图1）。气象数据则包括从中国气象数据网（http://data.cma.cn/）获取的2000—2019年京津冀及其周边地区34个气象站点的月平均气温、月降水量、相对湿度、平均风速以及日照时数数据，月太阳总辐射量数据利用日照时数数据计算获得，采用反距离权重法（Inverse Distance Weighted，IDW）（Shepard，1968）将各气象数据插值为空间分辨率为500 m的栅格数据。最后，采用双线性内插法（尤玉虎等，2005）将地形、植被类型、土壤类型等数据重采样为500 m分辨率，使所有数据像元大小一致。

图1 京津冀地形图及气象站点分布（a）、植被类型图（b）、土壤类型图（c）Fig. 1 Topographic and meteorological stations distribution map (a), vegetation type map (b) and soil type map (c)of Beijing-Tianjin-Hebei

1.3 研究方法

1.3.1 指标构建与信息提取

已有研究表明，气候和地形都在一定程度上影响着植被的生长（毛德华等，2012）。鉴于此，本文选取平均气温、年均降水量、平均风速、相对湿度、太阳总辐射量、海拔、坡度、坡向、植被类型、土壤类型共 10种涵盖气候、地形、植被、土壤 4个方面的自然因子，探测自然因子对京津冀 NPP的影响机制（表1）。

为使层内平均离散方差最小，而层间平均离散方差最大，采用自然断点法将平均气温、年均降水量、相对湿度、平均风速、太阳总辐射量、海拔、坡度等可量化数据分为9类，将坡向数据按平坡、北坡、东北坡、东坡、东南坡、南坡、西南坡、西坡、西北坡等分为9类。同时将植被类型数据分为常绿针叶林、落叶阔叶林、耕地、草地、灌丛以及非植被6类，将土壤类型数据按土纲分为淋溶土、半淋溶土、水成土、半水成土等12类（图1b和1c）。利用ArcGIS软件渔网点工具生成等间距（0.05°）的采样点9010个，并对NPP数据及各因子的栅格数据数据进行信息提取。

1.3.2 趋势分析

本文采用一元线性回归分析法逐像元地分析2000—2019年京津冀 NPP的变化趋势，并采用Mann-Kendall趋势检验法来判断其变化趋势的显著性，计算公式为：

式中，Slope表示NPP线性回归系数，Slope>0时NPP呈增长趋势，反之则呈减少趋势。n为研究年数，变量i为年序号，xi表示第i年的NPP值。

1.3.3 变异系数

本文利用变异系数（Coefficient of Variation，CV）来分析 2000—2019年京津冀 NPP的波动程度，计算公式为：

式中，CV为变异系数；n为研究年数；xi为第i年NPP值；¯为京津冀研究时段平均NPP。CV值越大，表明NPP变化程度越大；反之，表明NPP变化程度越小。

1.3.4 地理探测器

地理探测器是由Wang et al.（2010）提出的探测空间分异性，以及揭示其背后驱动因子与影响机制的一组统计学方法，该方法通过比较目标因子与影响因子空间分布的一致性来判定影响因子的解释力与贡献力（王欢等，2018）。其核心思想是：如果某自变量对因变量有重要影响，那么该自变量与因变量的空间分布具有相似性。地理探测器包括因子探测器、交互作用探测器、风险探测器以及生态探测器4种探测器。

（1）因子探测器：因子探测器主要通过因子解释力的大小来定量地衡量各自然因子对京津冀NPP空间分异的影响程度。其中，解释力q值表达公式为：

式中，h=1,…,L为影响因子的分层数；Nh与N分别表示影响因子在层h及整个京津冀的样本数；σh与σ分别表示影响因子在层h及整个京津冀的NPP值的方差。q的值域为[0, 1]，q值越大，表明该影响因子对京津冀NPP空间分布的解释力越强，q=0表示因子与NPP无任何关系，q=1则表示因子可以完全解释NPP的空间分布。

（2）交互作用探测器：交互作用探测器用来探测两影响因子交互作用时是否会增强或减弱对NPP的解释力，或者这些因子对NPP的影响是相互独立的。其判断依据如表2。

表2 影响因子交互作用类型Table 2 Types of the interaction between two influencing factors

（3）风险区探测器：风险区探测器用来判断两个因子在子区域之间的属性均值是否有显著差异，用于搜索NPP均值高的区域。该探测器用t统计量来检验：

式中，表示影响因子层h处NPP的均值；nh表示影响因子在层h的样本数；Var表示方差。

（4）生态探测器：生态探测器用于比较两影响因子X1与X2对京津冀NPP的空间分布的影响是否有显著差异，以F统计量来衡量：

式中，NX1与NX2分别表示影响因子X1与X2的样本数；L1与L2分别表示影响因子X1与X2的分层数；SSWX1与SSWX2分别表示影响因子X1与X2所分层的层内方差之和。

2 结果

2.1 NPP时空分布

2.1.1 NPP年际变化

2000—2019年京津冀植被NPP整体呈显著增长趋势，其增长速率为 5.85 g·(m2·a)−1（以 C 计，下同）。NPP 年均值在 235.05—386.36 g·(m2·a)−1之间波动，平均值为 323.99 g·(m2·a)−1。其中，最大值出现在2016年，最小值出现在2000年（图2）。

图2 2000—2019年京津冀年均NPP变化趋势Fig. 2 Trend of annual average NPP in Beijing-Tianjin-Hebei from 2000 to 2019

2.1.2 NPP空间变化及稳定性

2000—2019年京津冀年均NPP空间分布如图3a。从全区范围看，京津冀NPP介于25.97—749.29 g·(m2·a)−1之间；高值区主要集中在北部以及西部的燕山和太行山脉，其 NPP 值高于 400 g·(m2·a)−1；低值区主要集中在西北部坝上高原地区以及东南部平原地区，其 NPP 值普遍低于 300 g·(m2·a)−1。

就变化趋势而言（图 3b），2000—2019年京津冀植被整体恢复显著。NPP呈增加与减少趋势的面积分别占 97.62%和 2.38%，其中有 80.29%和0.31%的区域呈显著增加与显著减少趋势。京津冀北部地区NPP增加显著，NPP下降的区域主要零星分布于南部及沿海地区。

就变异系数而言（图 3c），2000—2019年京津冀植被 NPP整体情况较稳定，平均变异系数为17.25%。其中，变异系数低于15%的稳定区域占京津冀总面积的 36.31%，主要分布在西南与东北部部分地区；变异系数介于15%—40%的较稳定区域占京津冀总面积的 63.11%，主要分布于京津冀西北、中部以及东南部部分区域；变异系数大于40%的极不稳定区域主要分布在燕山山脉及太行山脉等高海拔山地地区。

图3 2000—2019年京津冀NPP均值（a）、变化趋势（b）及变异系数（c）空间分布Fig. 3 Spatial distribution of the mean (a), change trend (b), and coefficient of variation (c) of NPP in Beijing-Tianjin-Hebei from 2000 to 2019

2.2 NPP自然影响因子地理探测

2.2.1 单因子影响力探测

地理探测器因子探测结果表明，各影响因子对京津冀 NPP空间分布的影响程度有明显差异（表3）。各影响因子对应的q值大小排序为：平均气温>海拔>土壤类型>坡度>平均湿度>太阳总辐射量>植被类型>平均风速>降水>坡向。从影响因子解释力的角度来看，平均气温、海拔、坡度、土壤类型的解释力均在30%以上，是影响京津冀NPP空间分布的最主要因素；降水、平均风速、平均湿度、太阳总辐射量以及植被类型的解释力在 10%—20%之间，是影响京津冀 NPP空间分布的次要因素；坡向的解释力只有 1.12%，说明京津冀植被NPP的空间分布受坡向影响较小。

表3 单影响因子q值Table 3 q values of single impact factor

从各影响因子q值变化情况来看（图4），除平均风速外，2000—2019年其它各影响因子q值均呈上升趋势，表明除平均风速外的各影响因子对NPP空间分布的解释力逐步增强。分时间段来看，2000—2005年，除降水量和平均风速外，其余影响因子的q值均呈上升趋势，其中平均气温q值增长幅度最大；2005—2010年，降水量、平均风速及坡向解释力有所提升，且平均风速解释力已接近降水量，其余影响因子q值均呈下降趋势，其中相对湿度q值下降幅度最大；2010—2015年，相对湿度和太阳辐射量的解释力增加，其他影响因子q值呈下降趋势；2015—2019年，各个影响因子的q值均呈上升趋势，其中相对湿度上升幅度最大，且逐渐逼近平均气温的q值，成为对京津冀NPP空间分异解释力第二高的影响因子。

图4 2000—2019年各影响因子q值变化Fig. 4 Changes in q values of impact factors in 2000-2019

2.2.2 自然因子显著性差异分析

各自然因子对京津冀NPP空间分布的影响的差异显著性如表4所示。从表中可以看出，平均气温及海拔与其他影响因子对 NPP的影响均存在显著性差异，除坡度外，土壤类型与其它影响因子存在显著性差异，降水与平均风速对 NPP的影响无显著性差异。以上分析进一步表明，平均气温、海拔、土壤类型及坡度对京津冀 NPP空间分布的影响最大。

表4 影响因子显著性差异（置信水平95%）Table 4 Significant difference of influencing factors (confidence level 95%)

2.2.3 因子交互作用探测

因子交互作用探测结果表明（表5），京津冀NPP空间分异特征的形成并非只受单一因子影响，而是各因子共同作用的结果。同时也可以看出，各影响因子间的交互作用对京津冀 NPP的影响并非相加等简单的叠加过程，而是呈非线性增强及双因子增强关系。就解释力而言，平均气温与其它因子交互作用的解释力均达到45%以上。平均气温与土壤类型的交互作用最为明显（q=0.6112），其次是平均气温与降水（q=0.6028）、平均气温与太阳总辐射量（q=0.5899）、降水与海拔（q=0.5826）；交互作用最不明显的为降水与坡向，解释力只有12.3%。

表5 影响因子交互作用q值Table 5 q values of interaction factors

2.2.4 NPP高值区识别及因子层间NPP差异性分析

京津冀 NPP高值区以及各影响因子的适宜范围或类型如表6所示。结果表明，对于气候因子，随降水量及相对湿度的增加，京津冀 NPP逐渐增加，分别在602.46—655.64 mm、63.39%—67.50%时NPP达到最大值；随平均风速的增加NPP呈减少趋势，风速 1.56—1.82 m·s−1时 NPP 值最高；NPP随平均气温及太阳总辐射量呈先增加后减少的变化趋势，分别在8.37—9.31 ℃、5415.28—5489.42 MJ·m−2时 NPP值最大。随高程、坡度及坡向的不同，NPP呈波动变化趋势，分别在1670—2803 m、5.13°—7.41°、西南坡向处NPP达到最大值；常绿针叶林及淋溶土类型处的NPP值最高。

就有显著差异的分层组合百分比而言（表6），各影响因子的百分比在 44.44%—97.22%之间。各气候因子的层间差异相对较大，具有显著差异的百分比均在 90%以上，其中相对湿度的层间差异最大，百分比达到 97.22%；坡度和坡向的层间差异较小，有显著差异的百分比在49%以下，其中层间差异最小的是坡度，百分比只有44.44%。

表6 影响因子适宜范围或类型（置信水平95%）Table 6 Suitable ranges or types of impact factors (confidence level 95%)

3 讨论

3.1 NPP数据精度验证及时空变化特征

本文研究2000—2019年京津冀植被NPP时空分布特征发现，时间上，京津冀植被 NPP整体呈增长趋势，与以往研究结论一致（张莎，2015；朱利欣，2018；孙涛等，2020）。究其主要原因在于2000年退耕还林、还草生态工程实施后，京津冀生态质量有所改善，同时，诸多生态安全屏障的打造也对京津冀的生态环境起到了一定的保护作用（崔海宁，2011）。

空间分布上，京津冀 NPP均值由西北向东南呈低—高—低的阶梯式分布，燕山山脉及太行山脉，植被多为林地，NPP较高；坝上高原、东南部平原地区植被主要为草地和耕地，NPP值较低，这与张莎（2015）、吕国旭等（2017）、朱利欣等（2019）研究结果一致。京津冀植被整体处于稳定恢复阶段，西北坝上高原及山地地区人口密度小，是生态工程建设重点实施的区域，植被恢复显著，中部及沿海部分地区受城市扩张、农田灌溉等人类活动影响较大，植被出现退化现象。对比张莎（2015）、朱利欣等（2019）的研究结果，本文中京津冀中部地区植被退化趋势有所减缓，造成这种差异的原因可能与研究时间、数据来源及质量差异有关。

3.2 影响因子

从影响因素来看，平均气温是对京津冀植被NPP空间分布影响最大的因素，它可以影响植被生长环境的地表温度及空气温度，还可以通过影响光合作用、蒸腾作用等植被生理特征来直接影响植被的存活与生长。降水对京津冀植被 NPP空间分布的影响力较小，这与侯英雨等（2007）的研究结果稍有差异，表明静态角度下 NPP空间分异的主要影响因子与动态角度下 NPP变化的主要影响因子存在差异性（贺倩等，2020），主要原因可能在于京津冀地区耕地面积较大，农田灌溉等人类活动使降水对NPP空间分异的影响下降。土壤类型、海拔和坡度也是影响NPP空间分布的主要因素，土壤类型可以直接影响植被的生长，还可以限制植被的类型和分布（杨淑萍等，2019）；海拔直接影响植被类型的分布和气候状况，而坡度通过坡面侵蚀强度间接影响植被生长，可见，地形也是通过控制区域的水热和土壤条件来影响其他环境变化，进而对植被格局产生影响（李金珂等，2019；潘洪义等，2019）。此外，研究发现风速对京津冀NPP的影响呈下降趋势，这说明京津冀防风固沙等工作成效显著。

上述现象均表明各自然因子对京津冀植被NPP空间分布的影响存在一定的耦合关系。交互作用探测结果也进一步表明，与单一自然因子相比，因子两两之间的协同作用均增强了其对植被 NPP空间分异的解释能力。其中，解释力最高的3种交互作用均为平均气温与另一因子的协同作用。气温和土壤类型交互作用最强，一方面，气温和土壤类型本身是影响京津冀NPP空间分布的最主要因素；另一方面，气温是土壤养分的主要制约因素，土壤类型又直接关系着植被能否有效吸收利用养分（Knorr et al.，2005）。京津冀地处半湿润、半干旱地区的过渡带，平均气温和降水的协同作用对京津冀 NPP时空分异的影响也尤为重要，植被对水热组合条件变化的响应与适应将对植被的生理生态过程、物质积累与分配以及生态系统的结构和功能产生巨大影响（Butler et al.，2012；吕晓敏等，2015）。太阳辐射是植被进行光合作用的重要环境因子，影响着植被对大气CO2的吸收（张弥等，2009）。同时，太阳辐射变化也是造成地表温度和大气温度变化的主要原因，且气温与太阳辐射量呈正相关关系，因此，平均气温与太阳总辐射量的协同作用对京津冀NPP空间分布的影响也不可忽视。

未来京津冀地区生态保护与建设工作应着重考虑气温、土壤及地形对植被的影响，同时也应考虑多种自然因子之间的协同作用，从多视角、多维度探究自然因子对植被的影响（左丽媛等，2020），合理制定生态保护措施，加大生态保护力度。

3.3 不足之处

京津冀地区人口密度极大，人类活动对植被生长的影响也是巨大的，本文仅从气候、地形、植被和土壤等方面定量分析了自然因子对京津冀植被NPP空间分异的影响机制，而人类活动对于研究区NPP时空分异产生的综合影响还有待进一步的研究。同时，地理探测器目前对影响因子的空间分区还没有明确的划分标准，加上因子选取以及数据获取方面的限制，其因子解释力存在一定的主观性（潘洪义等，2019；陶帅等，2020）。

4 结论

本文以 2000—2019年京津冀 NPP数据为基础，探讨了研究区 NPP的时空演变特征，并选取气候、地形、植被、土壤等方面共10个自然因子，运用地理探测器定量探测了各影响因子对京津冀NPP时空分异的驱动机制，主要结论如下：

（1）时间上，2000—2019年京津冀植被 NPP整体呈显著增长趋势，其增长速率为 5.85 g·(m2·a)−1，年平均值为 323.99 g·(m2·a)−1。

（2）空间上，NPP高值区主要集中在北部及西部的燕山和太行山脉，低值区主要集中在西北坝上草原及东南部平原地区；2000—2019年京津冀植被整体恢复显著，NPP呈显著增加与减少趋势的面积分别占80.29%和0.31%；NPP整体情况较为稳定，平均变异系数为 17.25%，其稳定区域占京津冀总面积的36.31%。

（3）各影响因子对植被 NPP空间分异的解释力差异明显，其中平均气温、海拔、土壤类型和坡度是影响植被 NPP空间分异的最主要因素（q>30%，P<0.01）；2000—2019年，除平均风速外的各影响因子对 NPP空间分异的解释力呈上升趋势。

（4）自然因子对京津冀植被 NPP的影响存在交互作用，均呈非线性增强及双因子增强关系，其中平均气温与土壤类型的交互作用最为明显（q=0.6112）。同时，研究揭示了促进植被生长的各自然因子最适宜特征，以进一步了解自然因子对植被NPP的影响机制。