雪宝顶及其周边区域NPP时空格局演变及气候响应

邱 豪，王建华,2， 董廷旭,2，张雪茂， 严 霜

(1.绵阳师范学院资源环境工程学院，四川绵阳 621006；2.生态安全与保护四川省重点实验室，四川绵阳 621006)

0 引言

图1 雪宝顶区域地势图Fig.1 The area of Xuebaoding

植被净初级生产力(NPP)是指绿色植物在一定面积和时间内通过光合作用所产生的有机物总量减去呼吸作用消耗后的积累部分，是能够反应区域气候变化、区域碳积蓄与陆地生态系统响应的重要参数[1-3].NPP是物质与能量转化研究的基础，同时能够作为评价陆地生态系统可持续发展的一个重要生态指标[4-5].近几十年，NPP估算的模型陆续被提出，其中，光能利用率模型认为能够利用任何影响植被生长的要素估算NPP[6].1993年，学者Potter、Felid依次建立起以遥感数据作为估算模型数据源的CASA模型.Potter依据植被NDVI指数以及气象因子实现了基于光能利用率模型的全球陆地生态系统NPP的估算模拟[7].Field随后对CASA模型中光能利用率参数的取值给出了具体的界定[8]，使得CASA模型估算结果的精度和可靠度得到提高.CASA模型所需数据主要来源于遥感卫星数据，覆盖广泛且时间分辨率较高，且所需参数相对较少、误差小，CASA模型已应用于区域植被NPP长期变化的监测和NPP研究中[9-13].朱玉果等通过CASA模型分析宁夏省不同草地2001-2010年NPP的时空分布特征，并探索了不同草地NPP与气象因子的相关性[14];王耀斌等利用CASA模型结合遥感数据，分析了2000-2015年秦巴山区植被净初级生产力，探索了NPP的时空变化及其趋动因子[15].通过遥感手段研究区域植被净初级生产力已经成为如今一种重要手段.

雪宝顶区域位于四川省西北部，包含平武县、松潘县、北川县、茂县和黑水县，地处四川盆地与青藏高原的过渡带，是川西北地区重要的生态屏障，在维系川西北地区生态环境安全中有着重要作用.但是由于雪宝顶区域地形差异大，站点较少导致无法完整掌握雪宝顶区域NPP时空分布.因此，尚未见雪宝顶区域NPP相关研究，探究雪宝顶区域NPP对气候响应的报道也相对较少.本文利用2008-2017年的遥感影像和气象站点资料，通过CASA模型对近10a雪宝顶区域的NPP进行估算，并在此基础上分析雪宝顶区域NPP时空分布特征，探讨与气象因子的相关性，以期填补该区域NPP研究空白，为生态修复提供科学依据.

1 数据与方法

1.1 数据来源及预处理

气象数据来源于国家气象信息中心(http://www.nmic.cn/)选取2008-2017年平武县、松潘县、北川县及其周边共10个气象站的日降水量(mm)、月平均气温(℃)、太阳总辐射(MJ/m2)、站点编号与经纬度坐标.在ARCGIS10.2中，通过空间插值的方法得到研究区的月降水总量(mm)、月平均气温(℃)以及太阳总辐射(MJ/m2)分辨率为250 m*250 m的栅格图.

遥感数据选取来源于美国宇航局(https://ladsweb.modaps.eosdis.nasa.gov/)分辨率为250 m*250 m的全球植被16日合成产品MOD13Q1，对数据进行格式转化，将影像重投影为WGS_1984，利用ENVI软件对重投影后的影响进行研究区裁剪，并对所需波段进行合成处理.土地覆被数据来源清华大学的土地利用/土地覆被数据，并根据NPP估算需求将其分为针叶林、阔叶林等7大类.

1.2 NPP估算方法

本文采用了基于光能利用率的CASA模型对雪宝顶自然保护区进行净初级生产力估算，CASA模型的净初级生产力估算公式如下：

NPP(x,t)=APAR(x，t)×ε(x,t)×0.5

式中APAR(x，t)表示像元x在t月吸收的光合有效辐射， ε表示像元x在t月的实际光能利用率.

APAR(x，t)是根据植被对红外及近红外波段的反射特征实现的，植被吸收的光合有效辐射取决于太阳总辐射和植被本身特征，可由以下公式获得[7]：

APAR(x,t)=SOL(x,t)×FPAR(x,t) ×0.5

式中SOL(x,t)表示t月在像元x处的太阳总辐射量，FPAR(x,t)为植被层对入射光合有效辐射的吸收比例，常数0.5表示植被所能利用的太阳有效辐射占太阳总辐射的比例.

ε光能利用率的估算公式如下[16]：

ε=Tε1×Tε2×Wε×ε*

Tε1表示在低温和高温时植物内在的生化作用对光合作用的限制而降低植被净初级生产力，Tε2表示环境温度从最适温度向高温或低温变化时植物光能利用率逐渐变小的趋势，Wε反映了植物所能利用的有效水分影响光能转化率程度[17].

1.3 NPP时空分异模型

1.3.1 趋势分析 一元线性回归分析能够模拟每个栅格的变化趋势，通过显示每一时间段像元的空间变化特征，综合反映一定时间序列的区域格局演变规律[18].该方法通过拟合不同时间段数据，消减异常因子对NPP的影响，准确地反映区域NPP的变化趋势.Slope即雪宝顶区域的NPP近10a时间序列每个栅格像元的线性变化斜率反应了在10a的时间序列中的变化趋势，公式如下：

式中，Slope为像元NPP回归方程的斜率；NPPi为第i年的NPP值；n为研究时段跨度(2008-2017共10a).

当Slope<0时，表明该时间序列中，区域NPP处于下降趋势，且值越小下降趋势越明显；反之则为增加.趋势显著性通过相关系数进行检验，计算公式如下：

1.3.2 Hurst指数 基于重标极差(R/S)分析方法的Hurst指数最早是由英国水文学家赫斯特在研究尼罗河水库流量和贮存能力的关系时提出，能够有效定量描述时间序列信息自相似性和长相关性[19].

NPP时间序列NPPi，i=1，2，3，…，n，对于任意正整数m，定义该时间序列：

(1)差分序列

ΔNPPi=NPPi-NPPi-1

(2)均值序列

(3)累计离差

(4)极差

对于比值R(m)/S(m)≌R/S，若存在如下关系R/S∝mH，则说明分析的时间序列存在Hurst现象，H称为Hurst指数，该指数值与时间序列持续性关系见表1.

1.3.3 相关性分析 相关系数能够定量表达两个变量之间相关程度，本文针对NPP值与降水、太阳辐射与气温等气象因子的进行相关度分析：

图2 雪宝顶区域年均NPP分布特征Fig.2 Annual average NPP distribution characteristics in the Xuebaoding region

式中，r为X和Y之间的简单相关系数，Xi为第i月的NPP值，Yi月的相关因子值，分别为NP平均值和对应相关因子平均值，i为时间月变量.

2 结果与分析

2.1 NPP空间分布格局

2008-2017年雪宝顶区域的年平均植被净初级生产力如图2所示，研究区近10a的多年NPP均值为398.34 gCm-2a-1，东部山地丘陵地区(北川、平武)多年NPP均值主要集中在400～500 gCm-2a-1，其中雪宝顶自然保护区、千佛山等区域NPP均值超过600 gCm-2a-1；西北部高原草地地区多年NPP均值主要集中在300～400 gCm-2a-1；西南部多年NPP均值分布差异较大，最低不足100 gCm-2a-1，最高超过600 gCm-2a-1；由此，雪宝顶区域的NPP分布存在较强的空间异质性，这与雪宝顶区域的地理位置和气象特征有关.东部的丘陵山地地区处于四川盆地与青藏高原的边缘地区，年均温度在11～13℃左右，降水充沛.且植物群落主要以阔叶林和针叶林为主，植被覆盖度高，净初级生产力强.西北部处于青藏高原东缘，地形多为丘状高原，分布以草地和沼泽为主，植被覆盖度高，净初级生产力较强.研究区西南部多为高山河谷，海拔差异较大，低海拔地区以阔叶林为主，植被覆盖度高，净初级生产力强.研究区内的雪宝顶、三奥雪山因常年覆盖积雪造成植被覆盖度极低，净初级生产力弱.

2.2 雪宝顶区域NPP时间变化特征

雪宝顶区域2008-2017年NPP值变化如图3所示.近10a来，雪宝顶区域NPP处于波动下降趋势，下降趋势明显，线性增长率为-6.862 8 gCm-2a-1.2012年的NPP年均值最低，为318.65 gCm-2a-1，2008年的NPP均值最高，为403.06 gCm-2a-1.将研究区NPP均值划分为6个等级，统计其分布面积占比(图3)，近10a雪宝顶区域NPP值200～300 gCm-2a-1和300～400 gCm-2a-1的面积整体呈现增加的趋势，而400～500 gCm-2a-1和≥500 gCm-2a-1整体则不断下降.NPP值300～400 gCm-2a-1和400～500 gCm-2a-1面积占总面积的60%以上，其中NPP值300～400 gCm-2a-1面积常年超过总面积的50%.研究区NPP值年内动态变化呈现典型单峰的特性(图4)，三月份植被生长期开始的时候，NPP值急速增加，在六月份因为太阳辐射的原因，NPP值增加在一定程度上放缓或者出现下降，到7、8月份达到峰值，9月份以后随着生长季的结束，NPP值迅速下降.

不同地表覆被类型对应的年均NPP值的差异较大，其中，常绿阔叶林的多年NPP均值最高，达到了555.74 gCm-2a-1，主要分布在研究区东部丘陵高山地区以及西南海拔低的地区；其次是农用地、灌丛、高山草甸和高原草甸，其它地表覆被类型多年NPP均值都在300 gCm-2a-1以下，其中城市的多年NPP值最低仅为102.2 gCm-2a-1.

图3 2008-2017多年NPP均值年际变化以及年内变化Fig.3 The interannual and annual changes of NPP from 2008 to 2017

图4 不同地表覆被NPP年际以及年内变化Fig.4 The interannual and annual changes of NPP

2.3 雪宝顶区域NPP的空间变化特征

图5 雪宝顶区域NPP变化趋势Fig.5 Variation trend of NPP in Xuebaoding region

2.3.1 NPP的空间变化特征 本文通过一元线性回归逐像元分析近10a雪宝顶区域NPP变化趋势，从图5中可以看出，研究区95.7%的区域，其线性斜率小于0，仅有4.3%的区域，其线性斜率大于0.从空间变化分析，研究区自东向西NPP下降趋势逐渐减弱，其中趋势分布在-20～-10 gCm-2a-1的区域面积最广，占研究区的43%，而下降趋势最明显的区域集中研究区东部平武县和北川县；NPP呈现增加趋势的地区集中分布在黑河流域和雪宝顶，特别是黑河流域，增长率最高可达55 gCm-2a-1.

图6 NPP变化HURST指数Fig.6 HURST index of NPP

2.3.2 雪宝顶区域NPP变化持续性分析 本研究逐像元计算了雪宝顶区域近10a的Hurst指数(图6)以分析雪宝顶区域NPP变化趋势的持续性.雪宝顶区域的Hurst指数在0.21～0.87≥之间，均值为0.49，Hurst指数小于0.5的区域占雪宝顶区域的56%，而大于Hurst指数大于0.5的面积占雪宝顶区域的44%，表明大部分雪宝顶区域NPP具有较弱的持续性，其变化特征反向特征要高于同向特征.Hurst指数高值主要分布在研究区西北部若尔盖大草原边缘和西南部黑河流域地区，其值在0.75左右；低值最要分布在雪宝顶以东平武县境内.从不同的地表覆被类型分析得出，针叶林、阔叶林、山地草地和农牧地的Hurst均值小于0.5，其中，阔叶林的Hurst均值最低，仅为0.46；灌丛草甸、草甸草原和沼泽草甸Hurst均值大于0.5，沼泽草甸Hurst均值最大，达到了0.56，其次为灌丛.通过将研究区近10aNPP变化斜率与Hurst指数叠加分析获得近10a雪宝顶区域NPP变化持续性特征(图6).从图6中可以得到，NPP变化趋势减少转增加、持续减少、持续增加和增加转减少的面积分别占研究区总面积的54%、40%、4%和2%，即目前NPP呈现减少的大部分地区在未来会出现逆转，未来表现出逆转主要分布在雪宝顶以东的丘陵山地地区，地表覆被类型以阔叶林和针叶林为主；而仍有40%的区域NPP会持续减少，持续减少的区域主要分布在研究区西北地区的丘状高原地带，该区域地表覆被类型以草甸草原为主.

2.4 雪宝顶区域NPP与气象因子相关性分析

2.4.1 年际NPP与气象因子相关性分析 本研究逐像元分析了近10aNPP与同期气象因子(降水、气温和太阳总辐射)的相关关系(图7)，并进行显著性检验.结果显示，年际NPP变化与年际气温与降水波动关系不大，整体呈不明显负相关，这可能是因为该区域近10a年均降水和温度较为稳定，变化波动不大.其中，研究区NPP年际变化与降水波动呈正相关的区域集中分布在东部低海拔区域和黑河流域沿岸，该区域多为河谷气候，降水变化大且不稳定，地表覆被类型多为阔叶林为主；而年际NPP变化与太阳总辐射因子存在着显著的正相关，呈现正相关的面积超过了95%，因此，对太阳辐射量的依赖性较大.通过分析发现雪宝顶区域年际变化尺度上，太阳辐射量是影响NPP年际变化的主要气象因子.

图7 雪宝顶区域NPP与气象因子空间相关性 Tig.7 Spatial correlation between NPP and meteorological factors in the Xuebaoding region

2.4.2 不同地表覆被类型NPP与气象因子相关性分析 不同的地表覆被类型NPP对气象因子的也存在明显不同(表2).通过表2发现，各类地表覆被类型NPP与太阳辐射相关系数均在0.65以上，且通过显著性检验(P<0.05)的面积比均超过50%，说明研究区各类地表覆被NPP与太阳辐射存在显著正相关.研究区各类地表覆被与降水、温度的相关系数均为负值；其中，沼泽草地NPP与温度存在显著负相关，可能与沼泽草地中分布以富养植被为主；然而研究区大多表覆被类型与降水、气温存在不显著负相关，通过显著性检验(P<0.05)的面积比均在7%以内，可能与该区域常年降水稳定，气温波动不大，使得气象环境稳定且适宜，使得该区域大部分植被固碳量相对稳定.

3 讨论

本文采用光能利用率模型(CASA)对雪宝顶区域近10a的植被净初级生产力进行估算，并对该区域NPP的时空分布变化以及与气象因子的相关性进行分析，由于光能利用率模型(CASA)的机理与其它模型不同，其所估算的NPP值必然有所不同，且不可避免存在一定误差，就雪宝顶区域而言，本文对各类地表覆被的NPP值估算与前人对平武县、北川县等区域的模拟值与平均值基本一致[19-21]，且空间分布与美国宇航局植被净初级生产力产品MOD17A3基本一致.

虽然本文研究表明近10a雪宝顶区域NPP以-6.862 8 gCm-2a-1的速率下降，这可能是由于2008年地震以后，大部分区域植被遭到破坏，且期间不断出现各种地质灾害，使得植被恢复受到影响有关；但是未来雪宝顶绝大部分区域呈上升趋势，这可能与近些年各类灾后生态修复工程以及该区域地表覆被类型主要以森林、草地有关.本文受数据源限制，仅选取10a的时间段进行NPP估算，未能将灾前10a的时空分布与变化进行模拟分析，且在估算NPP时未能结合该区域地灾较多特性将地形因子与地灾影响因子加入模型之中，这是今后需要进一步研究的问题.

4 结论

本文采用CASA模型估算雪宝顶区域2008-2017年的NPP，简单分析了其时空变化特征与气候因子的关系，取得以下几点结论：(1) 雪宝顶区域NPP空间分布差异较大，以雪宝顶为中界呈东高西低的分布，东部是地表覆被多为阔叶林的丘陵山地地区区域，西部是地表覆被多为草甸草原的高原丘陵地带；(2) 近10a雪宝顶区域年均NPP值为398.34 gCm-2a-1，呈现波动下降的趋势，以年均-6.862 8g Cm-2a-1的速率减少；(3)雪宝顶近10a呈现整体下降的趋势，下降NPP面积占区域95.7%，主要集中在东部低海拔地区；自东向西下降率递减，且未来大部分区域呈现上升趋势，占区域总面积的58%；(4)在年际尺度上，太阳辐射是影响区域NPP的主要因素；且对于各类地表覆被，太阳辐射也是主要影响因素，仅有灌丛与沼泽草地存在显著负相关.