2004-2015年贺兰山自然保护区植被NPP时空变化与气候响应

李婷婷, 马 超,2, 郭增长

(1.河南理工大学 测绘与国土信息工程学院, 河南 焦作 454003; 2.河南理工大学 自然资源部矿山时空信息与生态修复重点实验室, 河南 焦作 454003; 3.河南测绘职业学院, 郑州 451464)

联合国大会于2015年9月25日一致通过了题为“改变我们的世界：2030年可持续发展议程”的第70/1号决议，将保护、恢复和促进生态系统可持续发展(目标15)列入全球17项可持续发展目标[1]。目标15具体侧重于可持续地管理森林、恢复退化的土地和成功防治荒漠化、减少退化的自然生境和结束生物多样性丧失。基于遥感技术的长时序变化监测和评估研究，可明确生态变化趋势，提升地方政府实施、监督和评估可持续发展目标的能力，为“中国落实2030年可持续发展议程进展报告”提供数据参考[2]。

在研究生态系统对全球变化的影响、响应和对策中，植被净初级生产力成为一项不可或缺的指标和核心内容[3]。植被净初级生产力(net primary productivity,NPP)是绿色植物总初级生产力(在单位时间和单位面积上通过光合作用产生的全部有机物同化量)扣除自养呼吸后的剩余部分[4]。NPP表示植物光合作用产物固定和转化的效率，直接反映植物在自然环境条件下的生产能力，是评价生态系统结构与功能协调及生物圈人口承载力的状况指标[5]。祝萍等[6]将植被NPP作为评价指标对中国典型国家级自然保护区生境状况时空变化特征进行研究。张镱锂等[7]将青藏高原自然保护区的植被NPP变化过程数据与其相邻等面积区域的植被NPP变化差异进行对比分析，评估自然保护区的保护成效。王培娟等[8]针对长白山自然保护区植被NPP主要影响因子的敏感性进行研究。可见在对自然保护区生态系统的变化检测和评估研究方面，植被NPP已成为一项重要指标。

中国境内山脉众多，其中不乏名山，但能以大尺度地理界线闻名者寥寥无几，贺兰山幸为其一[2]。贺兰山是中国北方干旱地理带的一片绿洲，是宁夏平原的气候屏障和区域生态环境变化的指示器，其生态走向关乎濒危物种的延续与北方干旱荒漠带宁夏段的生态平衡。植被作为干旱区敏感的环境因素，植被覆盖变化的研究已成为揭示干旱区自然环境变化及其规律的重要手段[9]。周梦云等[10]对宁夏贺兰山国家级自然保护区建立前后区域生态脆弱性时空格局变化进行研究。刘胜涛等[11]对宁夏贺兰山自然保护区森林生态系统净化大气的环境功能进行研究。顾延生等[12]对贺兰山中段植被类型及其覆盖变化进行研究。通过文献调研发现，以往的研究多集中于对贺兰山植被进行某些特定生理特征的研究，而没有对贺兰山地区的植被生物量时空变化同气候变化相结合，获得区域生态演化规律。

由于贺兰山自然保护区的局部山地效应，使对该区气候与植被的相关研究变得意义非凡。研究贺兰山自然保护区植被净初级生产力时空变化及其与气候因子的响应关系，可为地方政府提供区域植被生理状况和生产能力的相关信息，对荒漠化生态治理、自然保护区的管理、生态状况评价和促进社会经济可持续发展提供参考依据。

1 研究区概况

宁夏贺兰山为西北为数不多的森林分布区之一，既是荒漠区与荒漠草原区的分界线，又是干旱区与半干旱区的分界线，具有涵养水源和调节气候的作用[11]。由于其地理位置的特殊性和生态系统的脆弱性，早在20世纪80年代就已被列为国家级生态保护区。(贺兰山的管理机构始建于1950年。1982年7月1日，宁夏人大划定贺兰山为省级自然保护区，1988年5月国务院批准宁夏贺兰山自然保护区晋升为国家级自然保护区。)

以1 500 m等高线为基准，并向外延伸1 km来确定研究区的边界范围(图1A)。研究区地处宁夏贺兰山自然保护区核心区，近南北走向，西邻阿拉善高原，东接银川平原。地理坐标为北纬38°07′—39°32′，东经105°17′—106°40′，南北长157.63 km，东西宽120.40 km，面积5 769.25 km2，海拔为1 168～3 497 m。

宁夏贺兰山是我国东西部气候的分界，属季风气候向大陆性气候的过渡，也是我国200 mm等降水量线所经过的地区[12]。年均气温6～8℃，寒暑变化强烈，年均降水176.4～313.0 mm，全年干旱少雨。

贺兰山因海拔高，地势东缓西陡，导致自然气候产生显著垂直分异，构成独特的干旱区山地植被垂直带，孕育出荒漠草原，耐旱乔木、灌木，油松、山杨林和青海云杉等植被类型[13]，利用中国1∶100万植被类型空间分布图将研究区的植被类型进行合并，得到宁夏贺兰山植被类型空间分布(图1B)。

图1 宁夏贺兰山的地形与植被分布

2 数据与方法

2.1 研究数据

2.1.1 数据来源与处理 表达植物在自然环境下的生产能力的植被净初级生产力NPP数据，源于美国一级大气存档和分布式数据发布中心(Level-1 and Atmosphere Archive & Distribution System Distributed Active Archive Center, LAADS DAAC)，研究选取2004—2015年的MODIS (Moderate Resolution Imaging Spectro radiometer)数据产品：MOD17A3H NPP数据(https:∥ladsweb.modaps.eosdis.nasa.gov)，其空间分辨率为500 m，主要用于研究NPP的变化与气候响应的关系。该NPP数据是基于MODIS/TERRA卫星的遥感参数，通过利用BIOME-BGC模型与光能利用率模型建立的NPP估算模型模拟得到全球陆地植被净初级生产力年际变化数据集，目前已在全球不同区域的植被生长状况、生物量估算、环境变化检测等研究中得到验证和广泛应用[14-15]。利用MODIS数据重投影工具MRT (Modis Reprojection Tool)将NPP数据定义为UTM投影，WGS_84椭球的地理数据，用于分析2004—2015年贺兰山植被NPP变化特征。

高程数据源于美国太空总署(NASA)和国防部国家测绘局(NIMA)联合测量的航天飞机雷达地形测绘使命数字高程模型(SRTM3 DEM,http:∥srtm.csi.cgiar.org/SELECTION/inputCoord.asp)，水平分辨率为90 m，主要用于以1 500 m等高线来划定研究区范围。

气象数据源于资源环境数据云平台(http:∥www.resdc.cn/Default.aspx)，分别为2004—2015年年降水量空间插值数据集，用于分析降水量对植被NPP的影响；2004—2015年年平均气温空间插值数据集，用于分析气温对植被NPP的影响。

植被类型空间分布数据源于资源环境数据云平台(http:∥www.resdc.cn/Default.aspx)1∶100万中国植被图，通过ArcGIS对研究区植被类型空间分布图进行投影变换，将植被类型合并为阔叶林、针叶林、灌丛、草地和荒漠这5种，为植被NPP阈值分割提供基础参考。

2.1.2 数据质量验证 研究区属国家级自然保护区，大范围实地调查难以实施，加之研究时间序列长(12 a)，生物量分布的时间和空间异质性大，通过有限的实地调查所获得的平均生物量来推算整个研究区的生物量也会存在较大误差。针对类似情况，已有研究者[16]将NPP估算结果与土地利用分类数据产品进行对比，通过局部土地利用变化趋势来验证整体精度，具有片面性。本文在研究区均匀选择120个样本点，利用遥感数据产品MOD13A1，通过将样本点12 a的平均NPP与平均NDVI做对比来验证精度，检验结果显示二者具有非常高的拟合度(r=0.944，p<0.01)，MOD17A3H NPP数据质量可靠。

2.2 研究方法

2.2.1 偏差分析 平均NPP是基于像元对多年NPP数据进行计算求得的平均值，用于得到研究区多年的平均NPP空间分布图，其计算公式如下：

(1)

偏差反映某时间段内NPP偏离多年NPP均值的程度，即某段时间内NPP的盈亏情况[17]。其计算公式如下：

(2)

2.2.2 相关性分析 相关性分析[18]可以反映要素之间的相关程度和相关方向，本文采用Pearson相关系数法对植被年NPP与气温、降水量的响应关系进行研究，其计算公式如下：

(3)

2.2.3 趋势分析 一元线性回归法[17]可以对研究数据进行趋势分析，得到其变化速率，本文采用一元线性回归分析法对研究区总NPP、气温和降水量进行研究，其计算公式如下：

Yi=β0Xi+β1

(4)

式中：Yi为因变量；Xi为自变量；β0为斜率；β1为一随机变量。

(5)

式中：β0为变化速率；n为年数(时间序列为2004—2015年，即n=12)；k为时间序列；Xk为第k年的自变量。

3 结果与分析

3.1 MODIS NPP分类处理

通过对2004—2015年的NPP数据进行计算，得到贺兰山12 a的平均NPP空间分布图(图2)。依据贺兰山自然保护区主要植被类型的空间分布情况和NPP数值的对应关系，将研究区划分荒漠、草原、灌丛、针叶林、阔叶林5个等级，其分别对应NPP阈值为[0,65]，(65,85]，(85,150]，(150,250]，(250,500]，单位为g C/(m2·a)[13]。

12年的NPP均值分割图显示出贺兰山自然保护区生态系统的景观结构，植被NPP由外向内逐渐递增，即从低山浅山的荒漠草原带，到中低山区的灌丛林带，再到高海拔山区的针阔混交林带，植物按海拔垂直分布的层次明显，构成了当地典型的山地生态景观。其景观辨识度优于1∶100万中国植被覆盖专题图。

图2 2004-2015年平均NPP密度分割图

3.2 MODIS NPP时间变化分析

研究对贺兰山自然保护区2004—2015年植被初级净生产力(NPP)进行了统计分析(表1)。

表1 植被净初级生产力(NPP)统计

NPP区间在[0,65]的荒漠面积占比总体呈波动性减少的趋势，荒漠面积在2005年、2011年和2015年荒漠面积达3 500 km2以上，占研究区总面积的60%左右。12 a间荒漠面积占比的最大变幅为35.13%；NPP区间在(65,85]的草原面积占比总体呈波动性增加的趋势，草原面积在2010年达到961 km2，占研究区面积的16.66%。12 a间草原面积占比的最大变幅达9.65%；NPP区间在(85,150]的灌丛面积占比总体呈波动增加的趋势，灌丛面积在2012年达到最大值1 832.75 km2，占研究区面积的31.77%。12 a间灌丛面积占比的最大变幅达19.99%；NPP区间在(150,250]的针叶林面积占比较为平稳，约768.75 km2，占研究区面积的13.32%。12 a间针叶林面积占比的最大变幅为5.82%；NPP区间在(250,500]的阔叶林面积占比最少，且波动变化较大，整体呈增加趋势。阔叶林面积达225 km2以上，占研究区总面积的4%左右。12 a间阔叶林面积占比的最大变幅为4.33%。

2004—2015年研究区植被NPP均值为71.22 g C/(m2·a)，同期全国植被年平均NPP值273.5 g C/(m2·a)[19]，研究区植被NPP均值仅为同期全国植被年平均NPP值的26.04%；植被总NPP变化范围为0.324 4～0.501 5 Tg C/a，波动幅度达154.6%，总体在变化区间0.177 1 Tg C/a以内呈缓慢增加趋势；植被总NPP增长速率为0.281 6 g C/(m2·a)，远低于同期全国植被NPP增长率1.141 5 g C/(m2·a)[19]。偏差分析得到总NPP偏差值在12 a中有7 a低于平均水平，出现生态亏缺的年份多于生态盈余。

3.3 MODIS NPP空间变化分析

3.3.1 苏峪口国家森林公园剖面水平地带性变化 苏峪口国家森林公园是贺兰山自然保护区植被覆盖较好的区域，区内设有贺兰山森林生态系统定位观测研究站，研究其生态剖面对认识整个自然保护区生态演化具有代表性。以苏峪口国家森林公园核心区中点(东经105°54′60.00″，北纬38°42′53.96″)为原点，设置W—E和N—S两个剖面，分析核心区两个地理方位的植被生产力变化，可以形成二维的时空变化晕渲图(图3)，图中黑色实横线位置为苏峪口国家森林公园核心区中点位置所在，在W—E剖面上像素点的位置为28 pixel，在N—S剖面上像素点的位置为77 pixel。

图3 苏峪口国家森林公园NPP剖面晕渲图

在W—E剖面时序变化图中(图3A)，苏峪口国家森林公园植被净初级生产力在水平地带性分布上东部要优于西部，其中W—E剖面东部NPP总量为0.014 5 Tg C，占比66.55%，西部NPP总量为0.007 3 Tg C，占比33.45%。由于地理分界的原因，西部为大面积的荒漠和草原所占据，东部浅山地带有部分草原和灌丛区。从时间序列的整体来看，2004年、2007年、2012—2014年(5 a)植被净初级生产力整体要优于2005—2006，2008—2011年、2015年(7 a)，这与总NPP变化的偏差分析结果一致。在植被NPP较好的5 a里，荒漠面积明显减少，草原和灌丛面积相应明显增加，针叶林和阔叶林有小范围面积的增加。

在N—S剖面时序变化图中(图3B)，苏峪口国家森林公园植被净初级生产力在水平地带性分布上北部优于南部，其中N—S剖面北部NPP总量为0.042 5 Tg C，占比57.35%，南部NPP总量为0.031 6 Tg C，占比42.65%。受地形及日照影响，北半球植被NPP的高值部分主要分布在坡面北侧(阴坡)[20]，因沿N—S剖面线存在多个山坡，每座山坡的北侧植被NPP均高于南侧，所以植被NPP的高值部分呈北高南低交替分布且NPP等值线北侧梯度大，南侧梯度小。2004—2015年N—S剖面线北部除2004年、2012年因降水丰富使得灌丛面积增加近20个像素(约10 km)外，其余年份绝大部分均由荒漠和草原覆盖，且面积趋于稳定；2004—2015年N—S剖面线南部由荒漠、草原和灌丛交错分布，2004年、2007年、2012—2014年针叶林边际线向南部扩张了近40个像素(约20 km)，其余年份几乎全部由荒漠、草原和灌丛覆盖；2004—2015年N—S剖面线中部由针叶林和阔叶林交叉分布，且阔叶林林线在2004年和2012年所占范围居多。

3.3.2 贺兰山山脊线生态剖面垂直地带性变化 因研究区呈地质学上的NNE—SSW向狭长分布，为全面表达研究区生态剖面上NPP变化细节，沿贺兰山自然保护区地形山脊线建立生态剖面线，为研究生态剖面上NPP变化细节与DEM之间的对应关系，对研究区时间序列2004—2015年NPP剖面与DEM进行叠加分析(图4)。

图4 山脊线剖面与DEM叠加分析时序变化

由山脊线剖面与DEM叠加分析时序变化(图4)可知，植被NPP垂直地带性分布趋势明显，整体呈针叶林—阔叶林—灌丛—草原—荒漠的结构变化。由图4的研究区生态剖面的NPP与DEM叠加分析可知，植被NPP变化情况与DEM关系密切，越接近山脊线中点(越在研究区中心)，DEM值越大，植被NPP值就越大。海拔在1 500 m左右几乎无植被覆盖，植被NPP值接近0值；海拔在1 500～2 000 m对应的植被NPP值的变化范围为0～100 g C/(m2·a)；海拔在2 000～2 500 m对应的NPP值的变化范围为100～200 g C/(m2·a)，且变化趋势高度一致；海拔在于2 500 m时，植被NPP值基本维持在300 g C/(m2·a)左右。当海拔高于2 500 m以后，植被NPP值开始呈下降趋势，植被NPP值在200～300 g C/(m2·a)变化。

3.4 MODIS NPP的气候响应

3.4.1 气候变化情况 2004—2015年贺兰山区局地气候呈湿冷化(图5)。年均气温呈波动下降趋势，年均温度虽存在波动性但变化较为平稳，变化区间在1.5℃以内，最高值出现在2013年，为8.0℃，最低值出现在2011年，为6.6℃，气温变率为-0.01℃/a；年均降水量呈波动上升趋势，研究区年均降水量变化区间在14 mm以内，在2012年出现的最高年均降水量为312.9 mm，在2005年出现的最低年均降水量为176.4 mm，降水变率为2.77 mm/a。

图5 2014-2015年研究区气候变化趋势

3.4.2 气候与NPP的响应分析 利用SPSS®Pearson相关性分析得到年均气温与总NPP间的相关性系数r=0.050 2，置信度p=0.876，则年均气温与总NPP不相关。通过图6A可直观看出气温与NPP的空间相关性，[0.45，0.95]高度正相关区间主要集中在苏峪口国家森林公园附近区域，面积约为21.75 km2，占研究区面积的0.38%；[0.25，0.45)中度正相关区间主要零星分布在研究区中部，面积约为305.25 km2，占研究区面积的5.29%；93.13%的区域NPP与气温的相关系数均在[-0.25，0.25)低度相关区间范围内；剩余0.12%的区域NPP与气温的相关系数均在[-0.45，-0.95)中度负相关区间范围内，主要分布在于宁夏平原相交的边界区域。[-0.95,-0.45)高度负相关区间所占面积为0。

图6 气候与NPP相关系数空间分布

利用SPSS®Pearson相关性分析得到两者之年均降水与NPP间的相关性系数r=0.646 8，置信度p=0.023，则年均降水与NPP呈显著正相关。通过图6B可直观表示年级降水与总NPP的空间相关性，[0.45，0.95]高度正相关区间面积约为2 976.50 km2，占研究区面积的51.59%，所占比例最大，主要分布在草原和灌丛等对降水量响应较强的区域；[0.25，0.45)中度正相关区间面积约为787.25 km2，占研究区面积的13.64%，主要分布在针叶林区域；[0.00，0.25)低度正相关区间面积约为1 961.50 km2，占研究区面积的34.00%，主要分布在海拔较高的阔叶林区域；[-0.25，0.00)低度负相关区间面积约为41.25 km2，占研究区面积的0.71%；[-0.45，-0.25)中度负相关区间面积约为3.46 km2，占研究区面积的0.06%；[-0.95，-0.45)高度负相关区间面积约为0。

4 讨论与结论

4.1 讨 论

全球变化对陆地生态系统影响巨大，预计将成为21世纪大量物种濒于灭绝的主要因素[9]。近年来，国际上对气候变化影响生态系统NPP方面开展了较多工作，研究表明生态系统NPP的正遭受气候变化和人类活动的双重影响[17-18]。

贺兰山作为国家级自然保护区，滥砍盗伐得以有效禁止，且地处高寒、远离城市，几乎不受人类活动影响。因此，气候变化是决定研究区植被生态系统的主要因素，直接或间接地影响着植被的净初级生产力。且中国宁夏贺兰山自然保护区位于季风气候向大陆性气候的过渡区，对气候变化的响应极为敏感。研究表明，研究区植被NPP在12 a内总体呈增加趋势，植被NPP与降水的相关性大于与气温的相关性。由于2004年9月—2005年11月在宁夏发生历史性旱灾[21]，大部分区域降水仅达到199.9 mm，为历史上仅次于1982年(192.4 mm)的第二低值年，年降水日数仅为37 d，导致2005年总NPP成为2004—2015年研究期最小值；2006年宁夏地区出现的罕见暴雨天气扭转了历史性干旱的局势[22]，使得2006年总NPP未受到2004—2005年干旱年的影响，NPP总量得以恢复；2012年宁夏贺兰气象站点监测到日降水量创1951年有气象记录以来的最高值[23]，使得2012年NPP值显著增大。说明在干旱与半干旱地区，降水变化是影响植被NPP的主导因素。这一结果与马安青[24]、刘雪佳[25]、焦伟[26]等的研究结果一致。

研究还发现，在海拔3 000 m左右的区域，植被NPP与气温呈中高度正相关(r∈[0.25，0.45])的同时，对降水呈低度正相关(r∈[0.00，0.25])，说明海拔改变了气候对研究区植被NPP的影响规律。但本文并未将高程、坡度、坡向等地形因素与气候因素结合起来研究其与植被NPP的响应关系，这应是下一步的工作重点。

4.2 结 论

(1) 时间变化规律。贺兰山2004—2015年平均NPP为97.91 g C/(m2·a)；总NPP变化范围为0.324 4～0.501 5 Tg C/a；NPP增长速率为0.281 6 g C/(m2·a)。贺兰山自然保护区是以山地针叶林和青海云杉为代表的荒漠区与半荒漠区典型森林生态系统，森林面积占比仅为10%～20%，森林生态状况不容乐观。

(2) 空间分布规律。以苏峪口国家森林公园核心区中点做W—E，N—S剖面，得到长时序水平地带分布规律：东部植被净初级生产力要优于西部，北部植被净初级生产力要优于南部；将贺兰山沿山脊线生态剖面与DEM叠加分析得到垂直地带分布规律：植被NPP值整体呈现出针叶林—阔叶林—灌丛—草原—荒漠的垂直景观结构变化，NPP高值集中分布在海拔2 500 m以上的高海拔地区，在海拔1 500 m左右几乎无植被覆盖。

(3) 气候响应规律。贺兰山区局地气候呈冷湿化趋势，气温变率为-0.01℃/a，降水变率为2.77 mm/a。相较于气温而言，降水量对研究区NPP的影响较大。植被NPP与气温不相关(r=0.050 2，p=0.876)；植被NPP与呈降水量显著正相关(r=0.646 8，p=0.023)，且高度正相关区域主要分布在草原和灌丛区。