甘肃祁连山自然保护区草地时空变化及其对气候的响应

宋伟宏，王莉娜，张金龙

（1.甘肃省环境监测中心站，甘肃 兰州 730020；2.中国科学院西北生态环境资源研究院/甘肃省遥感重点实验室，甘肃 兰州 730000）

祁连山是由森林、草地交错分布组成的复合生态系统[1-3]。草地作为祁连山生态系统的重要组成部分，在维持区域生态平衡，阻止腾格里、巴丹吉林和库姆塔格沙漠南侵，保障黄河和内陆河径流补给等方面具有不可忽视的作用[4-5]。同时，草地也是祁连山地区畜牧业发展的重要生产生活资料，与牧民的家庭经济收益及区域经济发展密不可分[6-7]。及时掌握祁连山草地植被生长状况、演变趋势及其与气候的关系，对区域复合生态系统功能的科学评估具有参考价值，对指导地方生态环境系统修复治理具有重要借鉴意义。

遥感技术以其速度快、范围广、连续性好和成本低等优点被广泛应用在草地变化的监测中[8-11]。Akiyama等[12]曾提出遥感监测结合地理信息系统在草地监测方面很有前景，而且随着卫星传感器的发展和地理信息系统的推广，精准、实时的草地监测也将会成为可能。Zhao等[13]基于遥感技术，实现了对阴山北麓草地生态功能区植被覆盖率的动态监测。冯双双[8]基于Landsat TM和Landsat 8 OLI TIRS影像建立了草地植被退化遥感监测模型，实现对坝上草原草地退化的动态监测。Eckert等[14]研究发现MODIS NDVI时间序列分析对于探测植被变化区域和识别土地退化和再生是适宜的。葛静等[15]基于UAV和MODIS遥感数据对黄河源东部地区的高寒草地盖度动态进行了监测。李肖娟[16]基于遥感数据对祁连山NPP时空变化特征和趋势做了分析。Shao等[17]基于多期Landsat卫星图像对青藏高原阿坝县草地保护工程前后的草地变化进行了监测，从中得到的草地退化和再生的变化规律有效地辅助青藏高原环境恢复措施的制定。气候变化对植被的影响已经被国内外众多学者关注并展开研究，且主要体现在植被NDVI与气温、降水的关系和响应特征上[18-21]。Gu等[22]基于MOD13Q1 NDVI数据对红河流域的植被覆盖与气候因子的响应进行研究。Duo等[23]对华北平原过去33年的植被覆盖度进行了监测。有研究表明，降水在影响西北地区植被变化的自然因素中起着最主要作用[24]。张禹舜[25]基于MODIS遥感数据研究了近11年祁连山年际、年内、生长季及不同季节植被NPP变化及其与气候因子的响应。Guo等[26]通过重构SPOT_Vegetation NDVI数据集实现对中国北方锡林郭勒草原生长气的遥感监测并分析了与气象因子的关系。王巧玲等[1]对祁连山北麓草原植被NDVI与气候的关系研究发现，天祝县和肃北县的NDVI分别与气温和年降水量显著正相关。

本研究以祁连山自然保护区为研究区，在前人研究的基础上结合多期土地覆被数据、植被指数和气象资料，分析祁连山自然保护区内不同草原类型NDVI空间变化特征及趋势，并从季节尺度分析草地演变对气温和降水的响应特征，以期进一步认识该区域气候变化对草地演变的影响程度。

1 材料与方法

1.1 研究区概括

甘肃祁连山自然保护区位于甘肃省境内祁连山自然保护区中、东段，地跨武威、金昌、张掖3市 8 个县区，地理位置 97°23′34″-103°45′49″ E，36°29′57″-39°43′39″ N，北至祁连山森林分布下线与河西走廊相邻，南沿祁连山主脉与青海省接壤，西至肃南县界，东至天祝县界，平均海拔4 000-4 500 m。由于特殊的外营力作用及复杂的地质发展期，使得区域内的地形与气候形成了东西南北向差异明显的地貌形态，植被和土壤的水平与垂直地带性规律表现明显[27]，自上而下分为5个生物气候带，分别为山地荒漠草原气候带、山地草原气候带、山地森林草原气候带、亚高山灌丛草甸气候带和高山亚冰雪稀疏植被气候带。在涵养水源、维护区域生态平衡方面具有重要地位的祁连山区通过其内部的森林、草地等生态系统构成的复合生态系统调蓄区域降水、冰川和积雪融水。

区域内全年降水量主要集中在5-9月，年均降水量257.2～389.9 mm，年均降水量随海拔高度升高增多，在达到最大降水高度后，又随海拔高度升高而逐渐减少。祁连山区大部分地区年平均气温1～4 ℃，林区极端最高气温34 ℃，极端最低气温-30.2 ℃。≥ 10 ℃活动积温在海拔2 500 m以下地区为1 500～2 500 ℃·d，海拔3 000 m以上地区为200～500 ℃·d[28]。研究区具体地理位置如图1所示。

1.2 数据源与分析方法

1.2.1 数据源及处理

MODIS NDVI数据来源于美国地质勘探局(https://lpdaac.usgs.gov)的MOD13Q1，该产品是一个采用Sinusoidal投影方式的3级陆地植被数据，其空间分辨率为250 m，时间分辨率为16 d。本研究选用2000-2018年间共414景影像，基于MRT工具对其进行拼接、波段提取及投影转换，并利用祁连山自然保护区矢量边界对其裁剪，再通过ArcGIS软件提取2000-2018年16 d合成的NDVI数据。

土地覆被数据：基于国家土地覆被分类标准，将祁连山自然保护区土地覆被分为耕地、林地、草地、水域、城乡工矿居民用地和未利用地六大类。Landsat影像(表1)来自于美国地质勘探局网站，运用ENVI 5.1对影像数据进行辐射定标、大气校正等预处理，通过ArcGIS10.3中人机交互的方式解译了祁连山自然保护区2000年、2005年、2010年和2015年4期土地覆被数据。经验证，kappa系数达到86.93%，解译精度较高。

草原类型数据：来源于中国科学院1∶100万全国植被类型图，基于1∶100万中国植被图数据将甘肃祁连山草地植被分为温带丛生禾草草原、嵩草杂类草高寒草甸、禾草苔草高寒草原和温带丛生矮禾草矮半灌木荒漠草原四类，文中依次简称典型草原、高寒草甸、高寒草原和荒漠草原。

气候数据：来源于中国气象数据中心逐日气温与降水量数据(图1)，通过ANUSPLIN软件对其进行空间插值。ANUSPLIN是基于局部薄盘样条光滑函数开发的较为成熟的气象要素插值软件，可以同时批量处理多个空间表面，且在长时间序列气象数据的处理中有很大优势，将最终插值结果同植被NDVI做响应分析。

1.2.2 趋势分析法

采用一元线性方程分析研究区生长季NDVI对应像元的时间变化趋势，从而反映祁连山自然保护区整体空间变化规律及不同草原区域格局的变化，其公式表示如下：

图1 祁连山自然保护区地理位置及数字高程模型Figure 1 Location and digital elevation model of the Qilian Mountain National Nature Reserve

表1 祁连山自然保护区2000、2005、2010、2015年Landsat影像信息Table 1 Landsat image information for the Qilian Mountain National Nature Reserve for 2000, 2005, 2010, and 2015

式中：SLOPE代表斜率，其值大于0，说明变化趋势呈现改善的趋势，相反则是退化的趋势；n为年数；i为年序号，NDVIi表示第i年的NDVI最大。参照前人的研究[29]，将SLOPE划分为7个等级，分别为严重退化(≤-0.009 1)、中度退化(-0.009 0～-0.004 6)、轻微退化(-0.004 5～-0.001 0)、基本不变(-0.000 9～0.000 9)、轻微改善(0.001 0～0.004 5)、中度改善(0.004 6～0.009 0)、明显改善(≥ 0.009 1)。

1.2.3 气象插值方法

ANUSPLIN是基于薄盘样条理论的专业插值软件，包括普通薄盘和局部薄盘样条函数。由于局部薄盘光滑样条法除普通的样条自变量外，还允许引入线性协变量子模型，诸如温度和海拔之间的相关关系[30]，因此本研究采用该方法对所需气象站点数据(气温与降水量)进行插值。其原理表述如下：

式中：Zi为i点的因变量；χi为d维插值自变量向量，f(χi)是χi的 未知光滑函数；yi是p维的独立协变量矩阵，b为yi的 未知p维系数矩阵；ei是期望为0、方差为wiσ2的独立变量，wi是i的相对误差方差，σ2为常数，是所有数据点上的误差方差。需注意：当p=0时，即缺少bTyi，模型可简化为普通薄盘光滑样条，而若f(χi)为0时，即没有独立样条变量，ANUSPLIN不可用。本研究以经度、纬度作为独立变量，高程为独立协变量，对研究区气象点数据进行空间插值。

1.2.4 相关关系分析

本研究从栅格尺度上对不同草原植被指数NDVI与气候因素包括年降水量、年平均气温进行相关性运算，从而在空间上表达NDVI与气候因素的相关程度的分布特征。计算两个时间序列xt和yt的相关系数：

式中：σx,σy表示x和y的均方差；和y的均值；Cxy为x和y的协方差；rxy表示x和y的相关系数；n为时间序列长度。

相关性系数表达了两个变量之间的相关性程度，其值的范围[-1，1]，绝对值越大，说明相关性越好。本研究的相关性研究划分标准如表2所列。

表2 相关系数与相关程度分类Table 2 R values and correlation classification

2 结果与分析

2.1 祁连山自然保护区2000-2015年草地面积时空动态变化分析

分析2000年、2005年、2010年和2015年祁连山自然保护区草地的面积变化发现(表3)，草地是祁连山自然保护区内重要的土地覆被类型，其面积占比在47%左右上下波动。总体来说，祁连山自然保护区内草地的面积变化主要表现在前期（2000-2005年）减少，后期2010年和2015年又相对增加；2000-2015年其面积变化总体上呈现为减少，减少了1.98 km2。

表3 祁连山自然保护区2000-2015年草地面积变化Table 3 Grassland area change in the Qilian Mountain National Nature Reserve from 2000 to 2015

在退耕还林还草等一系列生态保护与修复工程实施的背景下，草地在生长过程中受到地表水热组合和人类活动的双重影响，所以祁连山自然保护区内草地面积的变化与气候变化和人类活动关系紧密。进一步分析2000年与2015年草地与其他土地覆被类型间的时空转移及变化特征(表4、图2)，结果表明：草地地类转移最明显的图谱单元为“32”和“23”，即草地转为林地和林地转为草地，转移面积分别为799.01和786.96 km2，转移区域主要分布在祁连山自然保护区的中部与东南部；其次为“63”和“36”，即未利用地转为草地和草地转为未利用地，转移面积分别为500.3和469.69 km2，转移区域主要集中在祁连山自然保护区的北部与西南边缘；再次为“31”和“13”，即草地转为耕地和耕地转为草地，转移面积分别为94.86和88.17 km2，主要分布于自然保护区的中部边缘和东南部边缘区域；再次为“34”和“43”，即草地转为水域和水域转为草地，转移面积分别为34.86和22.71 km2，其空间分布差异显著，主要集中分布在保护区的北部、中部、南部；最后为“35”和“53”，转移面积占比相对较小，草地转出为农村居民用地的面积为5.57 km2，城乡工矿居民用地转入草地的面积为4.12 km2，零星分布在研究区的东南边缘。

表4 祁连山自然保护区2000和2015年草地与其他土地覆被类型间的转移矩阵Table 4 Transitions between grassland and other land-use types in the Qilian Mountain National Nature Reserve between 2000 and 2015 km2

图2 草地地类转移图谱单元空间差异(a)与面积变化(b)Figure 2 Grassland class transition map (a) and area change (b)

综上可得，草地转为林地明显，这说明区域内退耕还林还草等一系列生态保护与修复工程的实施以及生态经济林的建设使得祁连山自然保护区生态功能区的生态功能进一步加强，有助于区域水源涵养功能的发挥。草地转为未利用地和城乡工矿居民用地的面积大于转入面积，数据量仅为保护区草地面积的0.4%，随着近年生态环境整改的进一步落实，数据会进一步减少。

2.2 不同草原类型的遥感监测

2.2.1 季节尺度下不同草原类型的时空变化趋势

计算分析2000-2018年典型草原、高寒草甸、高寒草原和荒漠草原4类草地春季、夏季、秋季和冬季NDVI的变化趋势(图3-图6)。

2000-2018年春季，4种类型草地长势变化差异不明显，以无明显变化为主(图3、表5)。整个自然保护区轻微改善的草地类型中，典型草原的面积最大，占保护区草地总面积的18.54%；其次是高寒草甸，面积占保护区草地总面积的7.58%。此外，高寒草甸的退化面积占自然保护区草地总面积的8.64%，在4个草地类型中退化最为显著，主要分布在自然保护区的中部和东部。

2000-2018年夏季，4种类型草地均以轻微改善和中度改善为主(图4、表6)。相比其他草地类型，高寒草甸轻微改善的面积最大，占到保护区草地总面积的20.06%。中度改善的草地中，典型草原的面积最大，占保护区草地总面积的23.44%，其次是高寒草甸，面积占保护区草地总面积的13.52%，高寒草原中度改善的面积占保护区草地总面积的8.13%，中度改善区域重要分布在保护区的中部和西北部。

2000-2018年秋季，4种类型草地均以轻微改善为主(图5、表7)。其中，典型草原轻微改善的面积最大，占保护区草地总面积的31.51%；其次是高寒草甸，面积占比为28.07%。

图3 2000-2018年祁连山自然保护区春季不同草原类型植被NDVI的空间变化趋势Figure 3 Spatial variation trends in the NDVI for different grassland types in the Qilian Mountain National Nature Reserve in spring of 2000-2018

图4 2000-2018年祁连山自然保护区夏季不同草原类型植被NDVI的空间变化趋势Figure 4 Spatial variation trends in the NDVI of different grassland types in the Qilian Mountain National Nature Reserve in summer of 2000-2018

2000-2018年冬季，4种类型草地同样以轻微改善为主(图6、表8)。相比其他草地类型，典型草原轻度改善的面积最大，占保护区草地总面积的27.28%；其次是高寒草甸，面积占比为24.76%。结合空间分布来看，高寒草甸在张掖市东南部靠近青海省门源县存在少量草地轻微退化，典型草原在张掖市东南部靠近金昌市永昌县同样分布了一些退化草地，其他区域草地普遍轻微改善。

综上，2000-2018年，祁连山自然保护区不同类型草地的长势变化在季节尺度上均存在较明显的差异，同一季节，4种草地类型草地的主要变化趋势相似。整个保护区除了春季，其他季节草地均以改善为主。草地在春季的退化面积要明显大于其他季节，草地在夏季的改善程度要明显大于其他季节，秋季改善程度大于冬季。

2.2.2 不同草原类型年际时空变化趋势

2000-2018年，祁连山自然保护区典型草原、高寒草甸、高寒草原、荒漠草原的NDVI总体呈现明显的增加趋势，增加速率分别为0.005 3·(19a)-1、0.003·(19a)-1、0.005 3·(19a)-1和 0.006 2·(19a)-1(图7)。各个草原类型的NDVI值均在2001年达到最小，而2018年达到最大值，其中，典型草原、高寒草原、高寒草甸和荒漠草原的NDVI值范围分别为0.38～0.53、0.35～0.50、0.39～0.47 和 0.26～0.41，可见高寒草甸的NDVI值比较稳定，而其他3种类型的NDVI值波动较大。

2.2.3 不同草原类型变化对气候的时空响应特征

本研究从栅格尺度上对不同草原植被指数NDVI与气候因素包括年降水量、年平均气温进行相关性运算，得到NDVI与气候因素相关程度的空间分布特征及对应的面积统计结果(图8、表9)。

图5 2000-2018年祁连山自然保护区秋季不同草原类型植被NDVI的空间变化趋势Figure 5 Spatial variation trends in the NDVI for different grassland types in the Qilian Mountain National Nature Reserve in autumn of 2000-2018

从典型草原植被NDVI与年降水的相关性空间分布来看，整体上NDVI与年降水量相关性较气温高。NDVI与年降水的相关性表达以轻度和中度相关为主，相关性表达为显著相关的区域主要分布在自然保护区的西部，占整个区域总面积的23.06%；NDVI与气温的相关性表达以轻度相关为主，相关性表达为中度相关、显著相关的区域主要分布在张掖市东部与金昌市接壤的区域，其中，显著相关的区域仅占整个区域总面积的3.62%。从高寒草甸植被NDVI与年降水的相关性空间分布来看，整体上NDVI与年降水的相关性表达以轻度相关为主，显著相关的区域主要集中在保护区的西南部边缘地区，占区域总面积的17.02%；NDVI与年平均气温的相关性较低，整体上以轻度相关为主，中度相关的区域在自然保护区中部零星分布，显著相关区域仅占区域总面积的3.59%。从高寒草原植被NDVI与年降水的相关性空间分布来看，整体上NDVI与年降水的轻度、中度和显著相关的面积差异较小，显著相关的面积占区域总面积的35.33%，在4类草地类型中面积占比最大。相比年降水，年平均气温与NDVI的相关性较低，整体上以轻度相关为主，轻度相关的比重占区域总面积的82.02%。荒漠草原植被NDVI与年降水的相关性表达以中度相关为主，显著相关占区域总面积的20.33%。与其他3类草地类型一样，年平均气温与NDVI的相关性水平要低于与年降水的相关性水平，显著相关的区域面积占10.16%。

总体而言，甘肃祁连山自然保护区内不同草原类型植被NDVI与年降水、年平均气温的相关性时空响应特征具有一定的差异，但整体来看，草地植被NDVI与年降水的相关性水平均高于年平均气温，与年平均气温的相关性水平主要表达为轻度相关为主。说明区域年降水比年平均气温对草原植被NDVI的影响程度更大。

分析祁连山自然保护区不同草原类型季节NDVI与当季、上一季气温相关关系表明(表10)，4类草原类型的季节NDVI与当季气温的相关性均不显著。其中，高寒草甸植被NDVI与气温的相关性最大，相关性系数为0.311，其余草原类型与植被NDVI的相关系数均小于0.3。不同草原类型植被NDVI与气温的相关性除高寒草甸外其余3个类型草原以夏季最高，而秋季最差。4类草原植被NDVI与上季度气温的相关性大多呈不显著特点，秋季各草原类型NDVI与气温的相关性最高，滞后性最为明显。其中，仅有典型草原秋季NDVI与夏季气温相关关系达显著水平。

图6 2000-2018年祁连山自然保护区冬季不同草原类型植被NDVI的空间变化趋势Figure 6 Spatial variation trends in the NDVI for different grassland types in the Qilian Mountain National Nature Reserve in winter of 2000-2018

表5 祁连山自然保护区不同草原类型春季变化趋势的面积统计Table 5 Spatial trends for different grassland types in the Qilian Mountain National Nature Reserve in spring 2000-2018

表6 祁连山自然保护区不同草原类型夏季变化趋势的面积统计Table 6 Spatial trends of different grassland types in the Qilian Mountain National Nature Reserve in summer of 2000-2018

表7 祁连山自然保护区不同草原类型秋季变化趋势的面积统计Table 7 Spatial trends for different grassland types in the Qilian MountainNational Nature Reserve in autumn of 2000-2018

表8 祁连山自然保护区不同草原类型冬季变化趋势的面积统计Table 8 Spatial trends for different grassland types in the Qilian Mountain National Nature Reserve in winter of 2000-2018

分析祁连山自然保护区不同草原类型季节NDVI与当季、上一季降水的相关关系表明(表11)，不同草原类型植被NDVI与降水量的相关性均高于气温，各草原类型与降水量的相关性大于0.3，其中，荒漠草原植被NDVI与降水量的相关系数为0.612，达到0.01显著水平。在4类草原类型中，夏季与秋季的植被NDVI与降水量的相关性明显高于春季与冬季。除夏季外，4类草原类型植被NDVI与上季降水量的相关系数基本上大于0.3，其中，典型草原春季NDVI与冬季降水量的相关系数为0.553(P＜0.05)；荒漠草原春季NDVI与冬季降水量的相关系数为0.626 (P＜0.01)；荒漠草原秋季NDVI与夏季降水量显著相关，荒漠草原冬季NDVI与冬季降水量显著负相关。春季不同草原植被NDVI与上季度降水量的滞后性最为明显，其次是冬季，而夏季滞后性最弱。

图7 2000-2018年祁连山自然保护区不同草原类型NDVI的变化趋势Figure 7 Variation in the NDVI for different grassland types in the Qilian Mountain National Nature Reserve from 2000 to 2018

3 讨论与结论

本研究在分析2000-2015年祁连山自然保护区草地面积变化的基础上，基于MOD13Q1 NDVI产品，结合草原类型，从季节、年际两个尺度分别分析了祁连山自然保护区草地植被NDVI的时空变化趋势，并讨论了不同草原类型NDVI对气候变化响应特征。结果表明:

1) 2000-2015年草地面积变化中，草地作为祁连山自然保护区内重要的土地覆被类型，其变化主要体现为2000-2005年减少，于2010年和2015年又相对增加。草地转出以转为林地最明显，转出为耕地和水域较为明显，草地转出为林地、耕地和水域的面积分别小于各自转变为草地的面积；草地转出为城乡工矿未利用地和未利用地的面积大于其转入为草地的面积。

2) 2000-2018年祁连山自然保护区除春季外，其他季节均以轻微改善为主，夏季植被NDVI的改善状况最佳。其中，春季高寒草甸的退化面积占自然保护区草地总面积的8.64%，而典型草原轻度改善草地面积占草地总面积的18.54%；夏季轻微改善的草地中高寒草甸的面积最大，占到保护区草地总面积的20.06%，中度改善的草地中典型草原的面积最大，占保护区草地总面积的23.44%；秋季保护区草地改善程度整体优于冬季。

3)祁连山自然保护区不同草原植被NDVI均呈现增加的趋势，且以荒漠草原的增加速率最大0.006 2·(19a)-1，其次是典型草原与高寒草原，均为0.005 3·(19a)-1，而高寒草甸的增加速率最小，仅为 0.003·(19a)-1。戴声佩等[31]基于 SPOT VGT-NDVI数据对祁连山不同草地植被NDVI变化趋势进行研究，发现典型草原和平原草地NDVI增加速率高于高寒草甸，虽本研究采用的数据和方法与其不尽相同，但得到的结果基本一致。

图8 2000-2018年不同草原类型NDVI与气温(左)、降水量(右)的相关性Figure 8 Correlations between NDVI and temperature (left panels), and precipitation (right)for different grassland types from 2000 to 2018

表9 不同草原类型NDVI与气温、降水量的面积比例Table 9 Percentage area of different grassland types in the NDVI with temperature and precipitation under different degrees of correlation %

表10 不同草原类型季节NDVI与气温的相关系数Table 10 Correlation coefficients between NDVI and temperature for different grassland types

表11 不同草原类型季节NDVI与降水量的相关系数Table 11 Correlation coefficients between NDVI and precipitation for different grassland types

4)祁连山自然保护区不同草原类型植被NDVI与年降水、年平均气温的相关性呈现出一定的差异，但整体上与年降水的相关性水平均要高于年平均气温，与年平均气温的相关性水平主要表达为轻度相关且鲜有显著相关。说明年降水比年平均气温对草原植被NDVI的影响程度更大。张禹舜[25]研究祁连山净初级生产力NPP变化对气候的响应中也同样指出草甸、草原和荒漠植被对降水的响应程度要明显大于气温。

5)祁连山自然保护区降水量变化对不同草原植被NDVI的影响高于气温，且与草原植被NDVI的滞后性较为明显，其中以春季的滞后性最强，其次是冬季，而夏季滞后性最弱。黄德青等[32]的研究表明祁连山自然保护区生物量与降水量正相关，而与气温未达到显著水平，这与本研究降水量是影响不同草原植被NDVI的主要因素的结论吻合。