第72卷第1期
2017年1月地理学报ACTAGEOGRAPHICASINICAVol.72,No.1January,20171982-2013年青藏高原植被物候变化及气象因素影响
孔冬冬1，张强2,3,4，黄文琳1，顾西辉1
（1.中山大学水资源与环境系，广州510275；2.北京师范大学环境演变与自然灾害教育部重点实验室，
北京100875；3.北京师范大学地表过程与资源生态国家重点实验室，北京100875；
4.北京师范大学减灾与应急管理研究院，北京100875）
摘要：根据NDVI3g数据，本文定义了18种植被物候指标研究植被物候变化情况。根据1:100
万植被区划，把青藏高原划分为8个植被区分。对物候变化比较显著的区域，采用最高温度、最
低温度、平均温度、降水、太阳辐射数据，运用偏最小二乘法回归（PLS）研究物候变化的气候成
因。结果表明：①青藏高原生长季初期物候指标，转折发生在1997-2000年，转折前初期物候指
标平均提前2~3d/10a；青藏高原末期物候指标转折发生在2004-2007年左右，生长季长度物候
指标突变发生在2005年左右，转折前末期物候指标平均延迟1~2d/10a、生长季长度平均延长
1~2d/10a；转折之后生长季初期物候指标推迟趋势的显著性水平仅为0.1，生长季末期物候指
标、生长季长度指标趋势不显著。②高寒草甸与高寒灌木草甸是青藏高原物候变化最剧烈的
植被分区。高寒草甸区生长季长度的延长主要是由生长季初期物候指标提前导致的。高寒灌
木草甸区生长季长度的延长主要是由于初期物候指标的提前，以及末期物候指标的推迟共同
作用导致的。③采用PLS进一步分析气象因素对高寒草甸与高寒灌木草甸物候剧烈变化的影
响。表明，温度对物候的影响占主导地位，两植被分区均显示上年秋季、冬初温度对生长季初
期物候具有正的影响，该时段温度一方面会导致上年末期物候指标推迟，间接推迟生长季开始
时间；另一方面高温不利用冬季休眠。除夏季外，其余月份最小温度对植被物候的影响与平均
温度、最高温度的影响类似。降水对植被物候的影响不同月份波动较大，上年秋冬季节降水对
初期物候指标具有负的影响，春初降水对初期物候指标具有正的影响。8月份限制植被生长季
的主要因素是降水，此时降水与末期物候指标模型系数为正。太阳辐射对植被物候的影响主
要在夏季与秋初。PLS方法在物候变化研究中具有较好的效果，本文研究结果将会对植被物候
模型改进，提供有力的科学依据。
关键词：物候；青藏高原；NDVI3g；PLS
DOI:10.11821/dlxb201701004
1引言
植物物候是指植物受气候和其他环境因子的影响而出现的以年为周期的自然现象[1]，是植物长期适应环境的季节性变化而形成的生长发育节律[2]。物候是植被的一个敏感性和关键性特征，它能够反应植被生长变化与气候变化，帮助理解植被生态过程、量化气象收稿日期：2016-09-18;修订日期：2016-11-28
基金项目：国家杰出青年科学基金项目(51425903);国家基金委创新群体项目(41621061)[Foundation:National
ScienceFoundationforDistinguishedYoungScholarsofChina,No.51425903;FundforCreativeResearch
GroupsofNationalNaturalScienceFoundationofChina,No.41621061]
作者简介：孔冬冬(1993-),男,河南周口人,博士生,主要从事生态水文研究。E-mail:kongdd@mail2.sysu.edu.cn
通讯作者：张强(1974-),男,山东沂水人,博士,教授,博导,主要从事流域气象水文学研究、旱涝灾害机理、流域地表
水文过程及其对气候变化的响应机制以及流域生态需水等研究。E-mail:zhangq68@bnu.edu.cn
39-52页
变化对陆地生态系统的影响[3]。青藏高原作为世界的第三极，高寒草地约152.15×104km2（占青藏高原总面积257.24×104km2的59.15%），是中国乃至亚洲的重要牧区之一，青藏高原的生态过程对保障中国乃至东亚生态安全具有独特的屏障作用[4]，因此，研究青藏高原物候变化具有重要意义。
虽然关于青藏高原的植被物候已经很多研究，但是较多关注于生长季开始时间的变化及成因[3,5-7]。尽管如此，关于生长季开始时间的变化和成因依然具有较大争议[3,5-7]。不仅仅针对生长季开始时间，本文定义18个物候指标研究生长季各个阶段物候指标的变化情况。不同区域，不同植被类型，植被物候也存在较大的差异，简单地对青藏高原物候指标求整体平均会覆盖物候指标的空间差异性，因此本文对青藏高原植被物候按照植被分区进行研究。对物候变化比较显著的区域，本文使用PLS判断气象变量的影响。在分析气候变化对植被物候影响时，植被对气象变化的影响具有滞后效应[8]。沈妙根等根据气象变量前n天值进行叠加，判断生长季前期气象变量的影响[9-10]。但是气候变量对植被物候的影响可能不是连续的，简单的进行前期叠加可能会掩盖一些现象。而PLS可以识别出是哪段时间气象因子对植被物候造成的影响，且兼具主成分分析和多元回归的优点。本文主要围绕以下两个问题进行研究：①青藏高原物候到底发生了怎样的变化；②气象因子对青藏高原物候变化的影响。本文根据NDVI3g数据，定义了18种植被物候指标研究植被物候变化情况。根据1:100万植被区划，把青藏高原划分为8个植被分区；对物候变化比较显著的区域，采用最高温度、最低温度、平均温度、降水、太阳辐射数据等，运用PLS研究物候变化的气候成因。
2数据
2.1NDVI数据
本文采用归一化植被指数（NormalizedDifferenceVegetationIndex,NDVI）分析植被物候变化趋势。NDVI与植被初级生产力、叶覆盖、生物量[11]具有很好的相关性，已经被广泛应用于量化植被生长趋势与生长过程[11]。本文使用的NDVI是由全球监测与模型组利用NOAA系列卫星合成分辨率为1/12°的半月NDVI第三代数据集NDVI3g（时间跨度为1981年7月-2013年12月，http://ecocast.arc.nasa.gov/data/pub/gimms/）。
为减少云、大气对NDVI的干扰，数据首先采用谐波时间序列分析法（HarmonicAnalysisofNDVITimeSeries,HANTS）对NDVI数据进行平滑处理。HANTS能够解析的代表生长季NDVI，被广泛的应用于物候指标提取[4,12]。由于物候提取方法很容易受到非生长季NDVI影响，而积雪覆盖会显著影响非生长季NDVI[13]，因此NDVI数据在提取物候指标时需要进行除雪处理[7]。本文利用日平均温度（半月中含有5天连续低于0°）确定可能存在的积雪覆盖[13]。在选择研究区域时，把植被较为稀疏的荒漠区和NDVI季节性变化不明显的区域进行截除。超过10年以上NDVI年内变化小于0.1的区域判定为该地区NDVI季节性变化比较微弱。同时除去含有多个生长季与生长季跨年的区域，截除方法参考GARONNA等[12]。
2.2气象数据
本文使用的气象数据来源于中国气象驱动数据集，该数据集是中科院青藏高原研究所开发的一套近地面气象与环境要素再分析数据集[14]，空间精度为0.1°×0.1°。本文使用的温度数据为该数据集3小时时间尺度温度数据，据此计算逐日最高温度、最低温度、平均温度；降水和太阳辐射数据采用日时间尺度数据。
根据1:100万植被图集[15]，本文把青藏高原划分为9个植被分区[13]，其中热带雨林区
。域植被季节性变化不明显，在选择研究区域过程中已截去，余下8个植被分区[16]（图1）
图1青藏高原植被类型分区
Fig.1VegetationclustersintheTibetanPlateau3方法
3.1植被物候指标的计算方法
NDVI数据经过除噪、除雪处理之后，本文定义了18种植被物候指数研究植被物候变化情况。NDVI数据计算物候之前，首先采用spline样条函数把NDVI数据插值到逐日，和参数函数不同，spline可以最大程度保留NDVI数据的变化形状[12,17]。
（1）NDVI阈值法（TRS）
NDVI-NDVIminNDVIratio=（1）maxmin
式中：NDVIratio第一次超过阈值的日期定义为生长季开始日期；NDVIratio第一次低于阈值
的日期定义为生长季结束日期[12,17]。如沈妙根等[18-19]，本文计算阈值取0.2与0.5，下文分别简称TRS2、TRS5。
（2）NDVI最大斜率法（DES）NDVI的斜率的峰值和谷值分别对应生长季开始与生长季结束时间，同样生长季开始和结束时间之间的间隔为生长季长度[12,18]。同时NDVI最大值对应的时间定义为峰值时间。
（3）曲率曲线法对于S型曲线，物候转移日期都可以根据局部曲率k′变化极值进行估计[20]。曲率曲线k'的两个局部最大值定义为返青期和成熟期，两个局部最小值定义为衰老期和休眠期。曲率曲线k'计算公式如下：
f″(t)k′=（2）(1+f′(t))
式中：f′(t)是NDVI序列的一阶导数；f″(t)是NDVI序列的二阶导数。
（4）Gu法f′(t)的最大值和最小值用来定义NDVI曲线在恢复期和衰老期的切线。恢复期切线与基线和最值线的交点对应的日期，定义为上升期和稳定期；衰老期切线与最值线和基线的交点对应的日期定义为下降期和衰退期[21]。考虑到中间季节的NDVI下降，Filippa等[17]定义了一条稳定线，来拟合稳定期与下降期之前的NDVI
，稳定线与衰退
线的交点对应的日期重新定义为下降期。
3.2分段回归
分段回归在植被NDVI量级，及植被物候转折[3]方面有较为广泛的应用。本文采用分段回归分析植被物候指标突变转折点。
β+βt+εt≤αy=01（3）β0+β1t+β2(t-α)+εt>α{
式中：t为时间（年）；y为响应变量（植被物候指标）；β0、β1、β2是回归系数；β0为截距，β1与β1+β2分别为转折前后趋势；ε为误差项；α为响应变量转折点的位置。t检验用来检测原假设H0：β2是否显著不为0（即突变前后趋势有无显著差别）[3]。当检验显著性水平达到0.1时，转折点达到了变异程度，α称为变异点。为保证趋势分析具有一定的数据长度，令转折点的范围在1986-2009年间。同时，采用Mann-Kendall（M-K）法检测转折前后物候指标变化趋势。
3.3偏最小二乘法回归
偏最小二乘法回归（Partialleastregression,PLS），兼具主成分分析和多元回归的优点，克服了预测变量相关导致的多元共线性[22]，PLS在植被变化归因[22]以及近年来植被物候变化研究中均有应用[6]。采用VIP指数（VariableimportanceIndex）挑选显著影响变量。VIP＞0.8，说明该变量具一定的解释意义；VIP≥1，说明该变量具有明显解释意义。本文采用较为严格的限制条件，VIP≥1时认为该变量影响显著。有关PLS的具体计算过程请参考文献[23]。本文采用PLS研究植被物候对气候的响应，对于物候发生显著变化的植被分区，首先把物候指标和气象变量根据植被分区计算区域平均。对于某一物候指标，该物候指标的区域平均值作为响应变量，对应的前期1-12月的区域平均气象变量（月平均最高温度、月平均最低温度、月平均平均温度、月累积降水、月平均太阳辐射共12×5=60个气象变量）作为独立变量。PLS模型系数类似于多元回归斜率，反映了物候指标对气候变量的敏感性，VIP≥1时说明自变量对因变量具有显著的解释意义。4结果分析
4.1不同植被分区物候均值
图2显示了青藏高原15种物候指标在各个植被分区的多年均值情况（未包含生长季长度指标），青藏高原生长季主要在4底到10月初，青藏高原东北地区物候明显早于西南地区。不同方法定义的物候指标均值分布存在明显差异。其中返青期、上升期、TRS2生长季开始时间对应的时间比较接近，大致在4月底5月初。成熟期主要分布在7月底到9月初，不同分区差异较大，温带草原物候偏早，温性草原偏晚，相差将近1个半月。NDVI峰值时间主要集中在8月初。衰老期、下降期分布相对集中，且较为吻合，平均出现在8月底9月初。在9月中下旬进入生长季末期。TRS2生长季结束时间不同分区波动较大，尤其是在常绿阔叶林和温性草原地区，可能是因为这些地区植被覆盖度比较高，优势物种衰退之后，底层的草本植被仍然可以保持一定的NDVI值[24]，因此这些地区TRS2生长季结束时间波动较大，阈值偏小（0.2）。Yu等[6]认为在使用阈值法提取青藏高原物候指标时，生长季结束时间取阈值0.6拟合情况较为接近真实。同时，TRS5生长季结束时间与DES生长季结束时间分布集中且较为相似，说明采用TRS5与DES推求的生长季结束时间较为接近真实。10月中下旬植被逐渐进入衰退期，11月中旬植被进入休眠期。衰退期与休眠期明显晚于采用阈值法和最大斜率法定义的生长季结束日期。
注：TRS2表示NDVI阈值为0.2；TRS5表示NDVI阈值为0.5；DES表示最大斜率法；下同。
图2青藏高原不同植被分区物候指标多年平均分布情况
Fig.2Multi-annual(1982-2013)averageofphenologicalmetricsinTibetanPlateauindifferentvegetationclusters4.2物候指标变化
青藏高原生长季初期物候指标，转折发生在1997-2000年（图4），很多研究认为1998年左右青藏高原生长季初期物候发生转折[25-26]，这与本文研究非常吻合。对于具有显著突变点的植被物候指标，上升期、TRS2、TRS5、DES生长季开始时间，转折前均呈现显著下降趋势（显著性水平达0.05）（图3a、3b），相应的下降趋势分别为2.3d/10a、3.5d/10a、3.7d/10a、3.3d/10a，生长季初期物候指标平均提前2~3d/10a。之前研究认为1982-1999年，区域平均生长季开始时间提前速率达4.5~10.2d/10a[3,19,25]，可能是本文采用MK计算线性趋势剔除异常值的影响，使本文得出的提前速率略小于前人研究。转折之后，上升期、TRS5、DES生长季开始时间等指标，呈推迟现象，延迟幅度在7~10d/10a，但显著性水平只达到0.1。其中TRS2生长季开始时间，转折之后推迟趋势达0.05显著性水平，但推迟幅度较小，仅4.4d/10a。众多生长季初期物候指标只有TRS2在转折之后趋势较为显著，因此转折之后初期物候指标推迟可信度较低。关于青藏高原生长季初期物候指标转折之后趋势具有较大争议[3,6-7]，但是青藏高原西南地区初期物候推迟
。对于青藏被普遍接受[4,19]，只有西南地区的温带草原和北部的荒漠呈推迟显著（图3b）
高原生长季末期物候指标，转折发生在2004-2007年左右（图4）。休眠期、衰退期、DES、TRS5生长季结束时间，转折前显著上升，这些物候指标平均延迟1~2d/10a。转折后生长季末期物候指标具有提前趋势，衰退期对应的提前趋势为-17.6d/10a（显著性水平0.1）；TRS5生长季结束提前趋势为-7.8d/10a（显著性水平达0.05）。Che等[27]认为青藏高原末期物候指标趋势不明显，1999-2011年趋势为0.96d/10a（p=0.3）
。同时，除衰
注：TRS2表示NDVI阈值为0.2；TRS5表示NDVI阈值为0.5；DES表示最大斜率法；下同。
图2青藏高原不同植被分区物候指标多年平均分布情况
Fig.2Multi-annual(1982-2013)averageofphenologicalmetricsinTibetanPlateauindifferentvegetationclusters4.2物候指标变化
青藏高原生长季初期物候指标，转折发生在1997-2000年（图4），很多研究认为1998年左右青藏高原生长季初期物候发生转折[25-26]，这与本文研究非常吻合。对于具有显著突变点的植被物候指标，上升期、TRS2、TRS5、DES生长季开始时间，转折前均呈现显著下降趋势（显著性水平达0.05）（图3a、3b），相应的下降趋势分别为2.3d/10a、3.5d/10a、3.7d/10a、3.3d/10a，生长季初期物候指标平均提前2~3d/10a。之前研究认为1982-1999年，区域平均生长季开始时间提前速率达4.5~10.2d/10a[3,19,25]，可能是本文采用MK计算线性趋势剔除异常值的影响，使本文得出的提前速率略小于前人研究。转折之后，上升期、TRS5、DES生长季开始时间等指标，呈推迟现象，延迟幅度在7~10d/10a，但显著性水平只达到0.1。其中TRS2生长季开始时间，转折之后推迟趋势达0.05显著性水平，但推迟幅度较小，仅4.4d/10a。众多生长季初期物候指标只有TRS2在转折之后趋势较为显著，因此转折之后初期物候指标推迟可信度较低。关于青藏高原生长季初期物候指标转折之后趋势具有较大争议[3,6-7]，但是青藏高原西南地区初期物候推迟
。对于青藏被普遍接受[4,19]，只有西南地区的温带草原和北部的荒漠呈推迟显著（图3b）
高原生长季末期物候指标，转折发生在2004-2007年左右（图4）。休眠期、衰退期、DES、TRS5生长季结束时间，转折前显著上升，这些物候指标平均延迟1~2d/10a。转折后生长季末期物候指标具有提前趋势，衰退期对应的提前趋势为-17.6d/10a（显著性水平0.1）；TRS5生长季结束提前趋势为-7.8d/10a（显著性水平达0.05）。Che等[27]认为青藏高原末期物候指标趋势不明显，1999-2011年趋势为0.96d/10a（p=0.3）
。同时，除衰
图3青藏高原物候指标突变前后趋势及全部年份整体趋势
Fig.3phenologicalmetrics'trendsbeforeandafterturningpointsandtrendsofthewholeresearch
period
注：2表示1982年，20表示2000年，30表示2010年；其中空白表示突变前或者突变后序列长度小于5年，
灰色表示90%置信水平下突变不显著。
图4青藏高原分段回归物候指标突变年份
Fig.4Theturningpointsofphenologicalmetricsdetectedbypiecewise退期、TRS5生长季结束时间，其余末期物候指标趋势不显著。对于青藏高原生长季长度物候指标，其突变时间与生长季末期物候指标突变时间相近，约在2005年。转折之前TRS5生长季长度平均延长2.4d/10a，DES生长季长度延长1.3d/10a。转折之后生长季长
度变化趋势不显著。
对于具有显著突变点的植被物候指标的全年趋势，大多数分区植被物候指标并未有明显的推迟或提前趋势，这可能由突变前后变化趋势完全相反导致的（图3a、3b、3c）。尽管如此，高寒草甸和高寒灌木草甸分区依然有物候指标趋势显著性水平达0.05。全年而言青藏高原的高寒草甸与高寒灌木草甸是物候指标变化最剧烈的区域（图3d、3e），可能是这些地区草地对气候变化比较敏感。在高寒草甸区，返青期、上升期提前幅度较大，达到2~5d/10a（显著性水平达0.05）；TRS5、DES生长季长度显著延长，而生长季结束时间等物候指标与休眠期、衰退期等趋势均不显著，这说明高寒草甸区生长季长度的延长主要是由生长季初期物候指标提前导致的。在高寒灌木草甸区，TRS2、TRS5、DES生长季长度指数均显著延长（显著性水平均达0.05）；返青期、TRS2生长季开始时间均显著提前，但提前幅度较小，分别为1.5d/10a、1.1d/10a；休眠期、衰退期推迟显著，平均推迟1~2d/10a。这说明高寒灌木草甸区生长季长度的延长可能是由于生长季前物候指标的提前，以及生长季后期物候指标的推迟共同作用导致的。高寒草甸与高寒灌木草甸是物候指标变化最剧烈的区域。因此下文，对这两个植被分区趋势变化显著的8个物候指标（返青期、上升期、休眠期、衰退期、DES、TRS2、TRS5生长季开始时间与生长季长度），采用PLS进一步探究其原因。
4.3气象因素对植被物候的影响
对于最高温度，全部植被分区、全部月份均呈现显著上升趋势，除10月份高寒草甸与常绿阔叶林地区显著性水平为0.1，其余均是0.05显著性水平显著上升（图5）。11月-2月增温幅度达到0.15℃/a，冬季和初春升温幅度最高，升温幅度明显高于全国。最低温度的升温幅度不如最高温度，且2-6月、10-11
月很多分区升温不显著。和最高温度类图51982-2013年青藏高原气象因素平均温度、最高温度、最低温度、降水、太阳辐射变化趋势
Fig.5Multi-annual(1982-2013)trendsofaveragetemperature,maximumtemperature,minimum
temperature,precipitationandsolarradiation
似，平均温度在冬季升温幅度最大，夏季和秋季升温幅度较小，且一些地区升温不显著。对于降水，大多数地区呈现上升趋势，降水显著上升主要集中在5-8月，而同时期太阳辐射则呈显著下降趋势。
4.3.1高寒草甸区
（1）生长季初期物候指标（返青期、上升期、DES，TRS2，TRS5生长季开始时间）温度对植被物候指标的影响强于降水和太阳辐射，平均温度、最高温度、最低温度对生长季初期物候指标的影响大致相同（图6a~6e）。上年9月-次年2月平均温度与返青期、上升期、DES，TRS5生长季开始时间的PLS模型VIP均大于1，具有显著解释意义。其中上年9-12月温度与生长季初期物候指标正相关，该时段温度上升（图5a~5c）导致初期物候指标延迟；而当年1-2月、4月、6月温度升高导致植被物候提前（图6a、6c、6e）。温度升高会使植被生长发育所需的有效积温提前达到；同时温度升高会提升酶的活性，加快植被物候进程[28]。上年秋末冬初降水与初期物候指标模型系数为负（图6a~6e），尤其在上年10、12月份，高寒草甸10月份降水显著上升（显著性水平达0.05），12月份降水趋势不显著。3、5月份降水对DES、TRS5生长季开始时间的影响，逐渐由正转负（图6c、6e）。5月份降水与上升期、TRS2生长季的PLS模型系数为正，可能是上升期、TRS2生长季开始时间定义的生长季开始时间偏早（图2）。上年夏末秋初太阳辐射与生长季初期物候指标正相关（图6a~6e），而当年的2-3月的太阳辐射对TRS2、TRS5生长季开始时间具有负的影响。可能是因为上年夏末秋初高太阳辐射会导致生长季结束时间推迟[13]，从而导致当年生长季开始时间推迟；而2-3月太阳辐射升高说明降雪减少，有利于有效积温的积累。同时太阳辐射升高，意味着光照时间延长、光周期缩短，刺激植被
返青。
注：方框内部系数表示VIP≥1；红色虚线表示年分割线，下同。
图6青藏高原高寒草甸均值方差标准化物候指标与气象变量PLS系数
Fig.6PLSmodelcoefficientsofthecenteredandscaledphenologicalandmeteorologicaldatainalpinemeadows
（2）生长季末期物候指标（休眠期，衰退期）
对于末期物候指标，温度对植被物候指标的影响依然占主导地位，相对于最高温度与平均温度，最低温度对末期物候指标的影响减弱（图6f、6g）。温度对末期物候指标的影响主要集中在当年夏末与秋季。上年7月-当年8月，温度与末期物候指标的模型系数接近于0；当年秋季温度对生长季末期物候指标的影响最大，模型系数接近于0.1。夏末最高温度、太阳辐射与末期物候指标负相关，最低温度为正相关。秋季温度上升是生长季末期物候指标延长的主要因素。降水与末期物候指标的关系比较复杂，模型系数正负变化较大。6月、8月太阳辐射对末期物候指标具有负的影响。
（3）生长季长度物候指标（DES、TRS2、TRS5生长季长度）
温度对生长季长度物候指标的影响占主导地位，平均温度与最高温度PLS模型系数类似（图6h~6j）。除当年8月份，其余时间最高温度对生长季长度的影响强于最低温度。上年9月至当年4月温度，对应的生长季长度与初期物候指标的PLS模型系数刚好吻合，说明该时段温度影响初期物候指标；夏末与秋季温度主要影响末期物候指标，从而调节生长季长度。夏末与秋季最小温度与最高温度对生长季长度指标的影响差别较大。8月份，最低温度、降水对物候有正的影响，太阳辐射为负的影响。可能是该时段太阳辐射升高，增加蒸腾作用，使植物可利用水分减少，该时段主要限制因素为降水。对于生长季长度物候指标，上年12月份降水影响为正，当年4月份降水为负，而8月份影响又转为正。可能是因为上年12月份与当年8月份降水增多，可以增加植被的可利用水资源，而4月份降水增多，意味着降雪增多，天气变冷不利于植被返青。而太阳辐射对生长季长度物候指标的影响主要体现在上年和当年8月份。8月份太阳辐射上升会提前生长季结束时间，上年生长季结束时间推迟可能会对当年生长季开始时间带来负的影响。
4.3.2高寒灌木草甸区从生长季长度PLS模型系数来看，TRS2与DES、TRS5生长季长度指标PLS模型系数差别较大，在夏、秋季与降水、太阳辐射的模型系数甚至和DES、TRS5对应的完全相反，这可能是因为TRS2在计算生长季结束时间时误差较大，从而生长季长度误差也很大（图7h~7j）。
（1）生长季初期物候指标（返青期、上升期、DES，TRS2，TRS5生长季开始时间）温度对高寒灌木草甸初期物候指标的影响不如高寒草甸（图7a~7f）。返青期PLS模型系数与其他初期物候指标差别较大，可能是采用GU法推求的返青期在高寒灌木草甸明显偏早导致的（图2）。而上升期、DES、TRS2、TRS5生长季开始时间PLS模型结果较为相似。与高寒草甸现象类似，上年温度（尤其是9、10月份）对初期物候指标具有正的影响。年后平均温度、最高温度对春季物候指标的影响则转为负，尤其是4月、6月（图7c~7e）。与高寒草甸不同的是，冬季平均温度、最高温度、最低温度，对初期物候指标不具显著解释意义。降水对初期物候指标的影响强于温度和太阳辐射。上年7月至当年4-6月，降水与初期物候指标PLS模型系数从正转负，年后又由负转正，尤其是在当年4月份。
（2）生长季末期物候指标（休眠期，衰退期）
对于生长季末期物候指标，温度的影响强于降水与太阳辐射（图7f、7g）。生长季末期物候指标的前3个月，平均温度、最高温度、最低温度对末期物候指标的影响较为一致，温度升高有利于末期物候指标推迟。可能此时升温会增强植被光合作用，减少叶绿素退化，延迟进入生长季末期[13]。5-6月份最低温度的PLS模型系数为负，此时最高温度对末期物候指标的影响不显著，可能5-6月份温度首先影响初期物候指标，最小温度升高使初期物候指标提前，从而使末期物候指标提前。降水、太阳辐射对末期物候指标的影
图7青藏高原高寒灌木草甸均值方差标准化物候指标与气象变量PLS系数
Fig.7PLSmodelcoefficientsofthecenteredandscaledphenologicalandmeteorologicaldatainalpineshrubmeadows响主要集中在当年8月份，此时降水、最低温度对末期物候指标具有正的影响，太阳辐射为负的影响，平均温度与最高温度对末期物候指标的影响比较微弱。可能此时限制植被生长的主要因素是最低温度与降水，而太阳辐射、最高温度增加会导致植被可利用水分减少[10]。
（3）高寒灌木草甸生长季长度（DES、TRS2、TRS5生长季长度）
与高寒草甸类似，上年温度对生长季长度有负的影响。可能的原因为，一方面上年末期物候指标推迟，会导致当年生长季初期物候指标相应的延迟；另一方面，气温升高不利于植被冬季休眠[6,29]，从而对生长季长度造成负的影响。但是高寒灌木草甸年前温度对生长季长度的影响主要集中在夏末秋初，且影响相对微弱（图7h、7j）。当年4-9月份平均温度、最高温度对生长季长度的影响都比较大，PLS模型系数为正，温度升高导致生长季长度延长；6月份，最高温度、太阳辐射对生长季长度具有正的影响，而最低温度对生长季长度具有负的影响。降水对生长季长度的影响集中在上年的秋冬季节，以及当年的秋季。上年冬末，降水对生长季长度有正的影响，该时段降水增多导致植被返青期可利用水分增多，诱导提前进入生长季。上年及当年秋季降水对生长季长度具有负的影响，可能是因为秋季限制植被生长发育的主要因素主要是最高温度与太阳辐射[28]，而降水增多，可能会导致最高温度和太阳辐射降低。太阳辐射对生长季长度物候指标的影响主要集中在当年6月、9月份，且具有正的影响（图7h、7j）。
5讨论与结论
5.1讨论
数据上，本文对NDVI3g
数据进行了精细的除雪处理。同时，为避免单一方法提取
物候指标的误差，定义了18种植被物候指标研究植被物候变化情况。考虑到区域的差异性，根据1:100万植被区划，把青藏高原划分为8个植被区分。方法上，PLS兼具主成分分析和多元回归的优点，克服了预测变量相关导致的多元共线性[22]，其得到了广泛的认可，并应用于诸多自然科学领域。对青藏高原物候变化最剧烈的两个植被分区（高寒草甸与高寒灌木草甸），综合考虑最高温度、最低温度、平均温度、降水、太阳辐射等，采用PLS进行成因分析。
物候变化是温度、降水、太阳辐射等因素综合作用的结果，温度对物候的影响占主导地位。研究发现，上年秋季、冬初温度对生长季初期物候具有正的影响。该结论与之前许多研究相符：张福春等采用积分回归法研究北京春季树木开花期实测物候资料，得出冬季气温偏高不利于休眠，反而使开花期推迟[29]。Guo[30]采用实测物候资料研究中国杏树物候变化、Eike等[31]采用实测物候指标研究加利福尼亚胡桃物候、Yu等[6]采用GIMMS-NDVI研究青藏高原物候变化，均发现冬季升温导致生长季初期物候指标推迟。夏季温度对植被物候的影响较为特殊，最小温度与最高温度对植被物候的影响存在明显差异：在高寒草甸，8月份最高温度、太阳辐射与末期物候指标负相关、最低温度正相关，高寒草甸多年平均最高温度与最大降水均在7月份，8月份降水明显减少，而太阳辐射、温度依然很高，可能是该时段较高的太阳辐射、温度导致植被蒸腾耗水增加，可利用水分减少，使植被加快进入枯黄。同时8月份降水与末期物候指标正相关也佐证了这一观点（图6f、6g）；在高寒灌木草甸，6月份最高温度、太阳辐射对生长季长度具有正的影响，而最低温度对生长季长度具有负的影响。可能是6月份最低温度的升高会导致植被夜间呼吸作用消耗的有机物增多，不利于植被生长。而最高温度、太阳辐射增加有利于植被进行光合作用，叶绿素增多。此外，秋季温度对末期物候指标具有正的影响，相对于夏季，秋季植被需水减少，温度升高会增强光合作用酶的活性，减缓叶绿素退化，从而推迟末期物候指标[13]。
降水对物候的影响不同月份波动较大。在高寒草甸，上年秋末冬初降水与初期物候指标模型系数为负，尤其在上年10、12月份；在高寒灌木草甸，上年冬季降水与初期物候指标模型系数为负。这与沈妙根等人的研究一致，生长季前期降水增多，会使青藏高原大多数地区春季物候提前[32]。反之，缺水则会推迟物候[28]。春初，降水与初期物候指标模型系数由负转正，可能是因为3、4月份降水主要以降雪的形式表现[26]，降雪增多会致使温度、太阳辐射偏低，不利于植被返青；在高寒草甸5月份降水与初期物候指标模型系数为负，可能5月气温回暖，降水的主要形式是雨水，降水增多有利于提高植被可利用水分，加速植被返青。8月份限制植被生长的主要因素为降水，降水增多促使生长季末期物候指标推迟，这与Che等人[27]研究结果一致。
相对于温度与降水，太阳辐射对初期物候指标影响较小，主要集中在6、8、9月份。这与Liu等[13]人结论一致，这可能是因为太阳辐射主要因纬度而不同，对气候变化不敏感（除多云的条件下）。在高寒草甸与高寒灌木草甸，6月份、8月份太阳辐射与末期物候指标负相关。6月份太阳辐射升高导致春季物候提前，间接导致生长季末期物候指标提前。8月份可能是夏季太阳辐射过高，太阳辐射下降反而有利于末期物候指标推迟[27]。在高寒灌木草甸，9月份太阳辐射对生长季长度指标具有正的影响，此时高太阳辐射可以加强植被光合作用能力，延缓脱落酸的堆积[13]，因此对生长季长度具有正的影响。
5.2结论
（1）青藏高原初期物候指标，转折发生在1997-2000年，转折前初期物候指标平均提前2~3d/10a；青藏高原末期物候指标转折发生在2004-2007年左右，生长季长度物候指
标突变发生在2005年左右，转折前末期物候指标平均延迟1~2d/10a、生长季长度平均延长1~2d/10a；转折之后生长季初期物候指标推迟趋势的显著性水平仅为0.1，生长季末期物候指标、生长季长度指标趋势不显著。
（2）高寒草甸与高寒灌木草甸是青藏高原物候变化最剧烈的植被分区。高寒草甸区生长季长度的延长主要是由生长季初期物候指标提前导致的。在高寒灌木草甸区生长季长度的延长主要是由于初期物候指标的提前，以及末期物候指标的推迟共同作用导致的。
（3）采用PLS进一步分析气象因素对高寒草甸与高寒灌木草甸物候剧烈变化的成因。两植被分区均显示上年秋季、冬初温度对生长季初期物候具有正的影响，该时段温度一方面会导致上年末期物候指标推迟，间接推迟生长季开始时间；另一方面高温不利用冬季休眠，从而秋季及冬初温度升高对初期物候具有推迟作用。除夏季外，其余月份最小温度对植被物候的影响与平均温度、最高温度的影响类似。降水对植被物候的影响不同月份波动较大，上年秋冬季节降水对初期物候指标具有负的影响；当年春初降水对初期物候指标具有正的影响。8月份限制植被生长季的主要因素是降水，此时降水与末期物候指标模型系数为正。太阳辐射对植被物候的影响主要在6、8、9月份。
PLS方法在物候变化研究应用中具有较好的效果，结果符合实际，大多数现象能在他人研究结果中的得到印证。本文研究结果将会对植被物候模型改进，提供有力的科学依据。
参考文献(References)
[1]ZhangXuexia,GeQuanshen,ZhengJingyun.RelationshipsbetweenclimatechangeandvegetationinBeijingusing
remotesenseddataandphenologicaldata.ActaPhytoecologicaSinica,2004,28(4):499-506.[张学霞,葛全胜,郑景云.北京地区气候变化和植被的关系:基于遥感数据和物候资料的分析.植物生态学报,2004,28(4):499-506.]
[2]LuPeiling,YuQiang,HeQingtang.Responsesofplantphenologytoclimaticchange.ActaEcologicaSinica,2006,26
(3):923-929.[陆佩玲,于强,贺庆棠.植物物候对气候变化的响应.生态学报,2006,26(3):923-929.]
[3]ZhangG,ZhangY,DongJ,etal.Green-updatesintheTibetanPlateauhavecontinuouslyadvancedfrom1982to2011.
ProceedingsoftheNationalAcademyofSciences,2013,110(11):4309-4314.
[4]DingMingjun,ZhangYili,SunXiaomin,etal.SpatiotemporalvariationinalpinegrasslandphenologyintheQinghai-
TibetanPlateaufrom1999to2009.ChineseScienceBulletin,2012,57(33):3185-3194.[丁明军,张镱锂,孙晓敏,等.近10年青藏高原高寒草地物候时空变化特征分析.科学通报,2012,57(33):3185-3194.]
[5]ShenM.SpringphenologywasnotconsistentlyrelatedtowinterwarmingontheTibetanPlateau.Proceedingsofthe
NationalAcademyofSciences,2011,108(19):E91-E92.
[6]YuH,LuedelingE,XuJ.WinterandspringwarmingresultindelayedspringphenologyontheTibetanPlateau.
ProceedingsoftheNationalAcademyofSciences,2010,107(51):22151-22156.
[7]ShenM,SunZ,WangS,etal.Noevidenceofcontinuouslyadvancedgreen-updatesintheTibetanPlateauoverthelast
decade.ProceedingsoftheNationalAcademyofSciences,2013,110(26):E2329.
[8]CongN,WangT,NanH,etal.Changesinsatellite-derivedspringvegetationgreen-updateanditslinkagetoclimatein
Chinafrom1982to2010:Amultimethodanalysis.GlobalChangeBiology,2013,19(3):881-891.
[9]ShenM,TangY,ChenJ,etal.Influencesoftemperatureandprecipitationbeforethegrowingseasononspring
phenologyingrasslandsofthecentralandeasternQinghai-TibetanPlateau.AgriculturalandForestMeteorology,2011,151(12):1711-1722.
[10]ShenM,PiaoS,ChenX,etal.Strongimpactsofdailyminimumtemperatureonthegreen-updateandsummer
greennessoftheTibetanPlateau.GlobalChangeBiology,2016,22(9):3057-3066.
[11]BuitenwerfR,RoseL,HigginsSI.Threedecadesofmulti-dimensionalchangeingloballeafphenology.Nature
ClimateChange,2015,5(4):364-368.
[12]GaronnaI,deJongR,deWitAJW,etal.Strongcontributionofautumnphenologytochangesinsatellite-derived
growingseasonlengthestimatesacrossEurope(1982-2011).GlobalChangeBiology,2014,20(11):3457-3470.
[13]LiuQ,FuYH,ZengZ,etal.Temperature,precipitation,andinsolationeffectsonautumnvegetationphenologyin
temperateChina.GlobalChangeBiology,2016,22(2):644-655.
[14]HeJie,YangKun.ChinaMeteorologicalForcingDataset.Lanzhou:ColdandAridRegionsScienceDataCenterat
Lanzhou,2011.[何杰,阳坤.中国区域高时空分辨率地面气象要素驱动数据集.兰州:寒区旱区科学数据中心,2011.]
[15]EditorialBoardoftheVegetationAtlasofChina,ChineseAcademyofSciences.VegetationAtlasofChina.Beijing:
SciencePress,2001.[中国科学院中国植被图编辑委员会.中国植被图集.北京:科学出版社,2001.]
[16]HouXueyu.VegetationGeographyandtheChemicalCompositionofItsDominantPlants.Beijing:SciencePress,1982.
[侯学煜.中国植被地理及优势植物化学成分.北京:科学出版社,1982.]
[17]FilippaG,CremoneseE,MigliavaccaM,etal.Phenopix:ARpackageforimage-basedvegetationphenology.
AgriculturalandForestMeteorology,2016,220:141-150.
[18]WhiteMA,deBeursKM,DidanK,etal.Intercomparison,interpretation,andassessmentofspringphenologyinNorth
Americaestimatedfromremotesensingfor1982-2006.GlobalChangeBiology,2009,15(10):2335-2359.
[19]ShenM,ZhangG,CongN,etal.Increasingaltitudinalgradientofspringvegetationphenologyduringthelastdecade
ontheQinghai-TibetanPlateau.AgriculturalandForestMeteorology,2014,189/190:71-80.
[20]KlineM.Calculus:AnIntuitiveandPhysicalApproach.2nded.DoverPublications,1998.
[21]GuL,PostWM,BaldocchiDD,etal.CharacterizingtheSeasonalDynamicsofPlantCommunityPhotosynthesis
acrossaRangeofVegetationTypes.SpringerNewYork,2009:35-58.
[22]HouMeiting,HuWei,QiaoHailong,etal.ApplicationofPartialLeastSquares(PLS)regressionmethodinattribution
ofvegetationchangeineasternChina.JournalofNaturalResources,2015,30(3):409-422.[侯美亭,胡伟,乔海龙,等.偏最小二乘(PLS)回归方法在中国东部植被变化归因研究中的应用.自然资源学报,2015,30(3):409-422.]
[23]WoldS,SjostromM,ErikssonL.PLS-regression:Abasictoolofchemometrics.ChemometricsandIntelligent
LaboratorySystems,2001,58(2):109-130.
[24]FridleyJD.Extendedleafphenologyandtheautumnnicheindeciduousforestinvasions.Nature,2012,485(7398):359-
362.
[25]PiaoS,CuiM,ChenA,etal.Altitudeandtemperaturedependenceofchangeinthespringvegetationgreen-update
from1982to2006intheQinghai-XizangPlateau.AgriculturalandForestMeteorology,2011,151(12):1599-1608.
[26]ShenM,PiaoS,DorjiT,etal.PlantphenologicalresponsestoclimatechangeontheTibetanPlateau:Researchstatus
andchallenges.NationalScienceReview,2016,2(4):454-467.
[27]CheM,ChenB,InnesJL,etal.Spatialandtemporalvariationsintheenddateofthevegetationgrowingseason
throughouttheQinghai-TibetanPlateaufrom1982to2011.AgriculturalandForestMeteorology,2014,189/190:81-90.
[28]WangLianxi,ChenHuailing,LiQi,etal.Researchadvancesinplantphenologyandclimate.ActaEcologicaSinica,
2010,30(2):447-454.[王连喜,陈怀亮,李琪,等.植物物候与气候研究进展.生态学报.2010,30(2):447-454.]
[29]ZhangFuchun.EffectsofglobalwarmingonplantphonologicaleventsinChina.ActaGeographicaSinica,1995,50(5):
402-410.[张福春.气候变化对中国木本植物物候的可能影响.地理学报,1995,50(5):402-410.]
[30]GuoL,DaiJ,WangM,etal.Responsesofspringphenologyintemperatezonetreestoclimatewarming:Acasestudy
ofapricotfloweringinChina.AgriculturalandForestMeteorology,2015,201:1-7.
[31]LuedelingE,GassnerA.PartialLeastSquaresRegressionforanalyzingwalnutphenologyinCalifornia.Agricultural
andForestMeteorology,2012,158/159:43-52.
[32]ShenM,PiaoS,CongN,etal.PrecipitationimpactsonvegetationspringphenologyontheTibetanPlateau.Global
ChangeBiology,2015,21(10):3647-3656.
VegetationphenologychangeinTibetanPlateaufrom1982to
2013anditsrelatedmeteorologicalfactors
KONGDongdong1,ZHANGQiang1,2,3,HUANGWenlin1,GUXihui1
(1.DepartmentofWaterResourcesandEnvironment,SunYat-senUniversity,Guangzhou510275,China;
2.KeyLaboratoryofEnvironmentalChangeandNaturalDisaster,BeijingNormalUniversity,Beijing100875,
China;3.StateKeyLaboratoryofEarthSurfaceProcessesandResourceEcology,BeijingNormalUniversity,
Beijing100875,China;4.AcademyofDisasterReductionandEmergencyManagement,
BeijingNormalUniversity,Beijing100875,China)
Abstract:UsingNDVI3gvegetationindex,wedefined18phenologicalmetricstoinvestigatephenologychangeintheTibetanPlateau(TP).Consideringheterogeneityofvegetationphenology,wedividedTPinto8vegetationclustersaccordingto1:1000000vegetationclustermap.Usingpartialleastregression(PLS)method,weinvestigatedimpactsofclimatevariablessuchastemperature,precipitationandsolarradiationonvegetationphenology.Resultsindicatedthat:(1)Turningpointsofthedateofthestartofgrowingseason(SOS)metricsaremainlyobservedduring1997-2000,beforewhichSOSadvanced2-3d/a.Turningpointsofthedateoftheendofgrowingseason(EOS)andlengthofgrowingseason(LOS)metricsarefoundduring2005and2004-2007,respectively.Beforetheturningpoint,EOShasadelayedtendencyof1-2d/10a,andLOShasalengtheningtendencyof1-2d/10a.Aftertheturningpoint,thetendencyofSOSandEOSmetricsisquestionable.Meanwhile,lengtheningofLOSisnotstatisticallysignificant;(2)Alpinemeadowsandalpineshrubmeadowsaresubjecttothemostremarkablechanges.LengtheningLOSofalpinemeadowismainlyduetoadvancedSOSanddelayedEOS.Nevertheless,lengtheningLOSofalpineshrubmeadowisattributedmainlytoadvancedSOS;(3)UsingPLSmethod,wequantifiedimpactsofmeteorologicalvariablessuchastemperature,precipitationandsolarradiationonphenologychangesofalpinemeadowsandalpineshrubmeadows,indicatingthattemperatureisthedominantmeteorologicalfactoraffectingvegetationphenology.Inthesetworegions,autumnoflastyearandearlywintertemperatureoflastyearhaveapositiveeffectonSOS.Firstly,increasedtemperatureinthisperiodwouldpostponelastyear'sEOS,andhenceindirectlydelaySOSofthecurrentyear;Secondly,warmingautumnandearlywinterhavethepotentialtonegativelyimpactfulfilmentofchillingrequirements,leadingtodelayofSOS.Exceptsummer,minimumtemperaturehasasimilareffectonvegetationphenology,whencomparedtoaverageandmaximumtemperature.Furthermore,precipitationeffectsonphenologyfluctuatewidelyacrossdifferentmonths.Precipitationoftheautumnandwinter/springofthelastyearhasanegative/positiveeffectonSOS.Besides,precipitationactsasthekeydriverconstrainingvegetationgrowthinAugust,duringwhichprecipitationhasapositiveimpactonEOS.Therefore,solarradiationcanexertimpactsonvegetationphenologymainlyduringsummerandearlyfall.Ourresearchwillprovideascientificsupportfortheimprovementofvegetationphenologymodel.
Keywords:phenology;TibetanPlateau;NDVI3g;PLS