青藏高原草地物候时空变异特征分析

Analysis of Spatial and Temporal Variability of Grassland Phenology on the Qinghai-Tibet Plateau

摘要: 气候变暖引起的极端自然灾害频发，对全球生态系统的结构和功能造成巨大挑战。青藏高原地理条件特殊，生态环境脆弱，表现出更优于其他地区的气候敏感性，而草地作为青藏高原主要的植被类型，对整个青藏高原生态系统有着不可忽略的影响，因此研究青藏高原草地物候变化具有重要价值。本文利用GOSIF数据集构建SIF (Sun-Induced Fluorescence)指数，同时采用Sen斜率估计对青藏高原草地物候时空演变格局进行探究，所得结论如下：1) 青藏高原草地返青期主要集中于第110~140天之间，其空间分布具有明显差异性，大体呈自东南向西北逐渐推迟的规律；枯黄期主要集中于第260~280天之间，其空间分布上并无显著规律；生长季长度主要集中于第130~189天之间，空间分布规律也较为明显，自东南向西北生长季长度逐渐变短。2) 在2001~2020年间，青藏高原草地返青期年际变化趋势为整体提前，提前约0.23 d/a左右；青藏高原草地枯黄期年际变化趋势表现为整体推迟，推迟约1 d/a左右；青藏高原草地生长季长度同样表现为整体延长趋势，延长约1 d/a左右。

Abstract: The frequent occurrence of extreme natural disasters caused by climate warming poses great challenges to the structure and function of global ecosystems. With special geographical conditions and fragile ecological environment, the Qinghai-Tibet Plateau is more sensitive to climate than other regions. As the main vegetation type on the Qinghai-Tibet Plateau, grassland has a significant impact on the entire ecosystem of the Plateau. Therefore, it is of great value to study the phenological changes of grassland on the Plateau. In this paper, the GOSIF data set is used to construct the index of SIF (Sun-Induced Fluorescence), while Sen slope estimation is used to explore the spatial-temporal evolution pattern of grassland phenology over the Tibetan Plateau, and the conclusions are as follows: 1) The greening period of grassland on the Qinghai-Tibet Plateau was mainly concentrated between the 110th and 140th day, and its spatial distribution was significantly different, and it gradually delayed from southeast to northwest. The yellowing stage was mainly concentrated between the 260th and 280th day, and there was no significant pattern in its spatial distribution. The growing season length mainly concentrated in 130~189 days, and the spatial distribution pattern was also obvious. The growing season length gradually decreased from southeast to northwest. 2) From 2001 to 2020, the interannual change trend of grassland greening stage on the Qinghai-Tibet Plateau was advanced by 0.23 d/a. The interannual variation trend of the Tibetan Plateau grassland yellowing stage was delayed as a whole by about 1 d/a. The length of grassland growing season on the Qinghai-Tibet Plateau also showed an overall lengthening trend, extending by about 1 d/a.

1. 引言

全球气候变化是当前科学界的热门话题，也是社会界和决策界共同关注的焦点。IPCC (Intergovernmental Panel on Climate Change)相关报告指出：自1950年以来，全球气候系统发生了过去几十年间从未有过的变化；1880~2012年间，在不同时期，全球绝大多数地区都经历了一定的升温过程，平均升高约0.85℃；1990~2010的30年间，地表温度以远超于1850年以来任何时期的速率增温，全球气候变暖现已成为不再争论的事实 [1] 。气候变暖引起海平面上涨、干旱暴雨等极端自然灾害频发，对全球生态系统结构和功能造成巨大挑战。地表植被作为生态系统的初级生产者，是陆地生态系统的重要组成要素，也是联系水文–土壤–大气和其他生态要素的重要桥梁。因此，植被动态变化在讨论有关气候变化方面发挥着不可或缺的作用，植被物候也作为反映植被动态变化的有效指示因子，成为当下研究的热点之一 [2] 。物候学是一门以研究动植物的生活周期以及这种周期怎样响应季节和年际环境要素的学科，在诸多领域中都具有非常重要的意义。植被物候是指植物因气候和生活环境等要素的影响而发生的以年为周期的规律性生长变化，不但能反映植物自身特征，也能直接体现植被乃至生态系统对气候变化的响应程度。近几年，伴随全球变暖加剧，植被自身的生长发育过程也出现了一定程度的改变，作为大自然的最佳指示器，植被物候的时空变异特征也逐渐受到相关学者的关注 [3] 。

青藏高原地处亚洲大陆中南部，海拔平均4000 m以上，是世界上最高、最年轻的高原，有着水平和垂直地带性的完美体现，气候、土壤和植被十分丰富，因其自身特殊的地理条件和脆弱的生态环境，成为全球变化的敏感区域。1980s中期以来，青藏高原表现出的增温速率是同纬度其他地区的两倍，该区域的平均温度更是每10年上升约0.4℃ [4] ，显现出高山地区更优于其他地区的气候敏感性，被称为北半球甚至是全球气候变化的放大器，研究其植被物候变化具有重要价值。草地资源作为我国的重要自然资源，占有国土面积的40%以上，在生态系统的调节和修复等方面发挥着毋庸置疑的作用 [5] ，且青藏高原天然草地作为青藏高原主要的植被类型，占地面积在60%左右，对整个青藏高原生态系统有着不可忽略的影响 [6] ，因此，研究青藏高原草地物候时空变异，有利于认识当地植被的空间分布格局，及时把握当地生态系统的变化情况，对于青藏高原生态环境保护及农牧业生产发展具有实际指导意义，是模型模拟和评估全球气候变化下生态系统变化研究的一个前提条件 [7] 。

2. 研究区概况

2.1. 地理位置

青藏高原南始喜马拉雅山南部，北达昆仑山、阿尔金山和祁连山北部，介于25˚~40˚N，68˚~105˚E之间 [8] [9] [10] ，大部分地区位于我国境内，但也有少数边缘地区在印度、巴基斯坦、尼泊尔等多个国家，总面积约2.5 × 10⁶ km²，大约占我国陆地总面积的三分之一左右 [11] [12] [13] 。青藏高原海拔范围在3000~5000米之间，发育了多类高山生态系统，是国内大江大河的主要发源地，其水资源的丰富储量使其拥有“亚洲水塔”的美誉，对我国自然生态环境和人类活动产生深刻影响 [14] [15] [16] 。

2.2. 气候特征

青藏高原巨大的海拔高度决定了当地的太阳辐射十分强烈，稀薄的空气使其能享受到比同纬度地区更多的太阳光照，但也正因如此，当地的气温普遍偏低，日较差大，从降雨来看，季风和暖湿气流由于受山脉阻挡的影响，其降雨自东南地区向西北地区逐级递减 [17] ，且以夜雨较多。因此，青藏高原的总体气候特点可以概括为：冬季寒冷漫长，夏季凉爽多雨，干湿分明，多大风 [18] 。

2.3. 植物资源

青藏高原是我国主要放牧区，天然草地类型众多，主要可以分为高山草甸、高山灌丛草甸、亚高山疏林灌木草甸、高山草原、高山山地荒漠、高寒沼泽、山地灌丛和高山稀疏及垫状草地八种，其中高寒草地分布面积最广，所占面积最大 [19] 。青藏高原也是我国重要的原始林区，以乔松、高山松、云南松等树种最为常见，尽管种类繁多，但整体覆盖面较小，所占比例较低 [12] 。

3. 数据来源与研究方法

3.1. 植被类型数据

本文植被类型数据采用1:100万全国植被类型数据，来源于中国科学院资源环境科学数据中心(http://www.resdc.cn/)。该数据集由侯学煜院士主编及相关院所、部门及高校的250余位专家共同编制，详细反映了我国多种植被类型的分布情况及水平和垂直地带性分异规律，是研究生态环境变化和生态环境保护的重要参考依据。该数据的空间分辨率为1 km，定义为Albers投影，在ArcGIS中利用掩膜、重分类、按栅格属性提取等工具获得青藏高原4种草地类型(见图1)。

3.2. 遥感数据

太阳诱导叶绿素荧光(Sun-Induced Fluorescence, SIF)作为新兴的遥感指数，是指植被在太阳光的照射下由光合中心直接发出的一种特有的光谱信号，并在红光波段(690 nm)和近红光波段(740 nm)出现两个峰值 [20] 。SIF代表的并非传统意义上的，类似于NDVI等由植被绿度来表征潜在光合作用的植被指数，其代表的更多是植被自身生长发育所进行的实际光合作用，更倾向于生理意义 [21] 。有研究发现，尽管植被类型不同，但GPP与SIF的相关性依然良好，由此说明，SIF在表示植被光合活性时仍有一定的代表意义 [22] [23] [24] [25] [26] 。

Figure 1. Grassland vegetation types and distribution in the Qinghai-Tibet Plateau

图1. 青藏高原草地植被类型及分布

本文所用2001~2020年SIF数据由Li等 [27] 提供，是基于OCO-2 (Orbital Carbon Observatory 2)卫星的全球高时空分辨率(0.05˚，8天) SIF数据集(即GOSIF)。GOSIF采用数据驱动方法，利用OCO-2离散测深数据、MODIS遥感数据和气象再分析数据构建预测SIF模型 [27] ，与原始OCO-2 SIF具有相似且合理的季节周期，但具有较高的时空分辨率、持续的全球覆盖和较长的数据记录，对于理解全球光合作用的长期趋势非常重要。

3.3. 物候参数提取

本文基于TIMESAT 3.3软件平台，使用双logistic曲线(D-L)拟合法对SIF数据的时间序列进行重构，以达到平滑数据的效果，之后选择常用的动态阈值法以5%和20%的阈值分别进行返青期(Start of the Growing Season, SOS)和枯黄期(End of the Growing Season, EOS)等各项物候指标的提取，其中生长季长度(Length of the Growing Season, LOS)通常由EOS与SOS做差得出 [28] 。

3.4. 趋势分析

Theil-Sen Median方法又被称为Sen斜率估计，是一种稳健的非参数统计的趋势计算方法，通常与Mann-Kendall (MK)显著性检验结合来揭示斜率变化的趋势和幅度。该方法计算效率高，受测量误差和离群数据的影响较小，在长时间序列数据的趋势分析应用较广。

Theil-Sen趋势分析：

$\text{Slope}=\text{Median}\left[\frac{\left(x_j-x_i\right)}{\left(j-i\right)}\right],\forall j>i$ (1)

式中，x_i和x_j分别是i和j时的序列值，1 ≤ i ≤ j ≤ n，n为序列长度，Slope的正负表示趋势变化的方向，若Slope大于零，则为正向变化，若Slope小于零，则为逆向变化。

Mann-Kendall趋势检验：

对于序列值 $X_i=\left(x_1,x_2,\cdots ,x_n\right)$ ，统计量S计算如下：

$S={\displaystyle \underset{i=1}{\overset{n-1}{\sum }}{\displaystyle \underset{j=i+1}{\overset{n}{\sum }}\mathrm{sgn}\left(x_j-x_i\right)}}$ (2)

$\mathrm{sgn}\left(x_j-x_i\right)=\{\begin{array}{c}1,x_j-x_j>0\\ 0,x_j-x_j=0\\ -1,x_j-x_j<0\end{array}$ (3)

当n ≥ 10时，统计量S接近标准正态分布，Z值可用于趋势检验：

$Z=\{\begin{array}{c}\frac{S-1}{\sqrt{VAR\left(S\right)}},S>0\\ \text{0,}S=0\\ \frac{S+1}{\sqrt{VAR\left(S\right)}},S<0\end{array}$ (4)

$VAR\left(S\right)=\frac{\left[n\left(n-1\right)\left(2n+5\right)-{\displaystyle \underset{i=1}{\overset{m}{\sum }}{t}_i\left(t_i-1\right)\left(2t_i+5\right)}\right]}{18}$ (5)

式中，n是序列中数据个数，m是序列中重复出现的数据组的数目，t_i是第i组重复数据组中的重复数据个数。标准化统计量Z的正值表示增加或上升的趋势，而负值则表示序列中的减少或下降趋势。为了检验显著性，将计算出来的Z取绝对值与从标准表获得的标准正态累积值Z_1−P/2在p%显著性水平上进行比较，以接受或拒绝零假设。

4. 结果与分析

4.1. 青藏高原草地物候空间分布格局

2001~2020年青藏高原草地生长季开始日期(a)、结束日期(b)和生长季长度(c)的空间分布如图2所示，可以看出青藏高原草地生长季开始日期主要集中在110~140天之间，即4月下旬到5月下旬的时间范围内，且生长季开始日期的分布特征具有较为明显的空间差异性，从图2(a)可以看到青藏高原草地返青大体呈现自东南向西北逐渐推迟的规律，返青期最早出现在东部，即在第99天左右开始返青，最晚主要出现在草地分布的西北部区域，在第146天左右开始返青。与返青期不同，青藏高原草地生长季结束日期在空间分布上并不存在一定规律(图2(b))，生长季结束日期主要集中在260~280天之间，即9月下旬到10月上旬的时间范围内，枯黄期最早出现在东北地区，即在第253天左右开始枯黄，像元所占比重极低；枯黄期最晚主要出现在西南区域，在第296天左右开始枯黄，像元比重也较小。对于生长季长度而言，青藏高原草地生长季长度也有着和生长季开始日期一样较为明显的空间分布规律(图2(c))，自东南向西北生长季长度逐渐变短，且生长季长度主要集中于130~189天之间，生长季长度最短，主要位于青藏高原东北角，在西部和中部地区也有零散分布，生长季长度仅为118天左右；生长季长度最长主要位于青藏高原东部及西南部，在偏北及中部地区也有零星分布，生长季长度长达189天左右。

(a) (b) (c)

Figure 2. Spatial distribution of SOS (a), EOS (b) and LOS (c) mean values in the Qinghai-Tibet Plateau from 2001 to 2020

图2. 2001~2020年青藏高原草地SOS (a)、EOS (b)、LOS (c)均值空间分布

4.2. 青藏高原草地物候年际变化特征

2001~2020年，青藏高原草地生长季开始日期、结束日期及生长季长度均值的年际变化如图3所示。从图3(a)的线性拟合结果来看，青藏高原草地生长季开始日期在2001~2020年间整体呈现提前趋势，但不同年份返青期的表现却各不相同。2004年之前，青藏高原草地返青期一直呈现较为稳定的提前趋势，并在2005年出现了第一次返青期推迟，而后在2008年之前再次持续提前，并在2009年又一次发生返青期推迟，但自2010年青藏高原草地返青期出现小幅提前后，青藏高原草地返青期出现频繁波动，自2015年后才逐步趋于稳定。从图3(b)的线性拟合结果来看，青藏高原草地生长季结束日期在2001~2020年间整体呈现推迟趋势。自2002年青藏高原草地枯黄期出现大幅提前后，青藏高原草地枯黄期在2003~2007年间频繁波动，2008年后才逐步趋于稳定，并在2013年后波动上升。青藏高原草地生长季长度在2001~2020年间的整体趋势可以从图3(c)的线性拟合结果看出，主要呈整体延长趋势。在2008年前，青藏高原草地生长季长度大体呈较为一致的延长趋势，但2009年出现大幅下降，之后上升幅度变小，2019年才出现大幅上升。

4.3. 青藏高原草地物候空间趋势分析

本节利用Theil-Sen Median方法对2001~2020年青藏高原草地SOS、EOS和LOS的空间变化趋势进行分析，结果如图4所示。各物候指标的Slope值表示斜率变化，若slope值为正，则表示该区域的物候变化为上升趋势，若Slope值为负，则表示该区域的物候变化为下降趋势。从图4(a)中可以看出，在2001~2020年间，青藏高原草地SOS的趋势变化范围大致在−1.7 d/a~1.4 d/a之间，青藏高原草地SOS的

Figure 3. Interannual variation characteristics of SOS (a), EOS (b) and LOS (c) in the Tibetan Plateau from 2001 to 2020

图3. 2001~2020年青藏高原草地SOS (a)、EOS (b)、LOS (c)年际变化特征

(a) (b)(c)

Figure 4. Spatial variation trends of SOS (a), EOS (b) and LOS (c) in the Tibetan Plateau from 2001 to 2020

图4. 2001~2020年青藏高原草地SOS (a)、EOS (b)、LOS (c)空间变化趋势

整体空间变化趋势以提前为主，提前约0.23 d/a左右，提前趋势占整个青藏高原草地SOS的76%左右，而推迟趋势仅在24%左右的地区呈现，主要集中在西南区域。青藏高原草地EOS趋势变化范围在图4(b)中体现，大致范围在−1.2 d/a~2.2 d/a之间，整体空间变化趋势以推迟为主，推迟约1 d/a左右，推迟趋势占整个青藏高原草地EOS的69%左右，提前趋势只占青藏高原草地EOS的31%左右，主要集中在东南部地区。图4(c)表示青藏高原草地LOS的趋势变化范围，大致在−1 d/a~2.7 d/a之间，青藏高原草地LOS的整体空间变化趋势以延长为主，延长约1 d/a左右，延长趋势占整个青藏高原草地LOS的88%左右，只有零散的，占整个青藏高原草地LOS 12%左右的地区呈缩短趋势。

同时将Slope结果结合MK显著性检验，对青藏高原草地物候指标SOS、EOS和LOS的空间变化趋势进一步说明，如图5所示。从图5(a)中可以看出，青藏高原草地SOS在大部分地区都以提前不显著趋势为主，占整个青藏高原草地SOS的63%左右，显著提前趋势仅占整个青藏高原草地SOS的13%左右，主要位于在东南、东北和中部地区，推迟不显著趋势占整个青藏高原草地SOS的23%左右，显著推迟趋势更是占整个青藏高原草地SOS的1%不到。在图5(b)中，青藏高原草地EOS大多以推迟不显著趋势为主，同样占整个青藏高原草地EOS的63%左右，而显著推迟趋势却占整个青藏高原草地EOS的6%不到，主要在东北和西南角有零星分布，体现不显著趋势则占整个青藏高原草地EOS的30%左右，显著提前趋势的占比仅接近1%左右。由图5(c)所见，青藏高原草地LOS以延长不显著趋势为主，占整个青藏高原草地LOS的74%左右，而显著延长趋势仅占整个青藏高原草地LOS的14%左右，主要位于东北、东南和中部地区，与SOS的显著提前趋势分布类似，缩短不显著趋势则占整个青藏高原草地LOS的12%左右，分布较为离散，显著缩短趋势甚至占整个青藏高原草地的0.1%不到。

(a) (b)(c)

Figure 5. Spatial variation trends and significance of SOS (a), EOS (b) and LOS (c) in grassland of Qinghai-Tibet Plateau from 2001 to 2020 (P < 0.05)

图5. 2001~2020年青藏高原草地SOS (a)、EOS (b)、LOS (c)空间变化趋势及显著性(P < 0.05)

5. 结论

本文主要从青藏高原草地物候空间分布格局、年际变化特征及变化趋势分析三个方面对青藏高原草地物候时空演变特征展开讨论，主要结论如下：

1) 青藏高原草地生长季开始日期主要集中在第110~140天之间，即4月下旬到5月下旬的时间范围内，且生长季开始日期的分布特征具有较为明显的空间差异性，大体呈现自东南向西北逐渐推迟的规律；青藏高原草地生长季结束日期在空间分布上并无规律，生长季结束日期主要集中在第260~280天之间，即9月下旬到10月上旬的时间范围内；青藏高原草地生长季长度也有着和生长季开始日期一样较为明显的空间分布规律，自东南向西北生长季长度逐渐变短，且生长季长度主要集中于第130~189天之间。

2) 青藏高原草地生长季开始日期在2001~2020年间整体呈现提前趋势，但不同年份返青期的表现却各不相同，2004年之前，青藏高原草地返青期一直呈现较为稳定的提前趋势，但2010年后，青藏高原草地返青期出现频繁波动；生长季结束日期在2001~2020年间整体呈现推迟趋势，并于2003~2007年间频繁波动，2008年后才逐步趋于稳定，且在2013年后波动上升；生长季长度在2001~2020年间整体呈延长趋势，但2009年出现大幅下降，之后上升幅度变小，2019年才出现大幅上升。

3) 在2001~2020年间，青藏高原草地SOS的趋势变化范围大致在−1.7 d/a~1.4 d/a之间，整体空间变化趋势以提前为主，提前约0.23 d/a左右，提前趋势占整个青藏高原草地SOS的76%左右；青藏高原草地EOS趋势变化大致范围在−1.2 d/a~2.2 d/a之间，整体空间变化趋势以推迟为主，推迟约1 d/a左右，推迟趋势占整个青藏高原草地EOS的69%左右；青藏高原草地LOS的趋势变化范围大致在−1 d/a~2.7 d/a以内，整体空间变化趋势以延长为主，延长约1 d/a左右，延长趋势占整个青藏高原草地LOS的88%左右。

4) 将斜率结果结合MK显著性检验，青藏高原草地SOS以提前不显著趋势为主，占整个青藏高原草地SOS的63%左右，显著提前趋势仅占13%左右，推迟不显著趋势占23%左右，显著推迟趋势更是占1%不到；青藏高原草地EOS以推迟不显著趋势为主，同样占整个青藏高原草地EOS的63%左右，显著推迟趋势占6%不到，提前不显著趋势则占30%左右，显著提前趋势的占比仅接近1%左右；青藏高原草地LOS以延长不显著趋势为主，占整个青藏高原草地LOS的74%左右，而显著延长趋势仅占14%左右，与SOS的显著提前趋势分布类似，缩短不显著趋势则占12%左右，显著缩短趋势甚至占0.1%不到。