青藏高原特殊的地理环境对气候变化具有放大效应(冯松 等，1998)，所以青藏高原成为了全球气候变化的脆弱区和敏感区(李林 等，2010)。青藏高原冻土分布广泛，气候变化显著影响其发育和分布(王绍令 等，1996)。近地表土壤的冻融循环过程与冻土的时空分布、陆面/水文过程以及生态系统有紧密的关联(Li 等，2012)，土壤的冻融作用会改变土壤剖面的水热状态和生化特性，进而影响大气与地表间的能水交换、径流分布和植被生理状态(杜子银 等，2014)。随着全球气温的升高，近地表冻融循环已经受到影响。在1967年—2013年间，中国的近地表年内冻融循环次数以0.39±0.05 cycles/10a的速度减少，并且循环次数的降低与气温的升高有显著的负相关关系(Peng 等，2016)。1978年—2008年，中国地表土壤冻结时间减少了34.3±16.5 days(Jin 等，2015)，其中，青藏高原的变化最为明显，近地表土壤冻结时间以16.8 days/10a的速度减少(1988—2007年)(Li 等，2012)。更多永久冻土的融化，增加了大气中温室气体的含量，促进气候变暖的趋势，同时也会加速寒区水资源的蒸发和流失，导致沙漠化的发展(Wang 等，2000)。

近年来，微波遥感技术的发展为监测地表冻融状态和冻融循环过程提供了有效的途径。液态水的复介电常数明显不同于与冰等其他物质，因此微波信号可以有效地观测土壤中液态水含量的变化(谢燕梅 等，2013)。例如，当气温降低时，土壤中的液态水发生相变，转化为固态冰，此时土壤复介电常数(包括实部和虚部)均明显减小，使地表后向散射系数和亮度温度发生显著变化，据此实现对地表冻融状态的监测(曹梅盛 等，2006)。微波的波长长，不仅可以穿透地表，探测土壤内部信息，而且受大气中云等物质的干扰较小，还能够全天时地进行观测，所以微波遥感在监测地表冻融状态方面相比其他波段具有一定的优势；被动微波遥感能够获得逐日的影像数据，可以利用其长期稳定地观测较大区域的地表冻融状态(晋锐 等，2009)。目前常用的被动微波遥感影像判别近地表冻融状态的算法主要有3种：双指标算法(Judge 等，1997；Zhang和Armstrong，2001；Zhang 等，2003；Zuerndorfer 等，1990；Zuerndorfer和England，1992)、决策树算法(Jin 等，2009)和时间序列变化检测算法(Smith 等，2004)，它们本质上都是根据冻融循环的微波辐射和散射特性提出和发展的。

春季地表温度回升，冻土融化，植被生长并进入新的物候循环，它们的物候信息也能够敏感地指示气候变化。返青期是指植物生长季的开始，它的变化不仅能够指示气候变化，还能够显著影响青藏高原草地生态系统的生产能力、结构组成，进而影响地表和大气的能量、水分和碳循环，反作用于气候变化(Piao 等，2008)。传统的返青期观测方法主要依靠地面实测，这种方法多用于个体和物种水平上的物候观测，在大范围整体植被观测上有一定的局限。遥感技术的发展为植被的返青期等物候观测提供了新的方法，而且这种观测是在区域的群落和生态系统层面上，具有更强的宏观性(Lüdeke 等，1996)。遥感观测方法是基于遥感影像时间序列植被指数来提取返青期的，利用时间序列遥感影像识别植被返青期主要有两个步骤：(1)时间序列植被指数的拟合去噪与重建；(2)基于植被指数时间变化曲线的返青期提取(范德芹 等，2014)。由于近几十年来的全球气候变暖，许多学者对青藏高原植被的物候期进行了长时间序列的观测和趋势分析。有学者(Zhang 等，2013；丁明军 等，2011；范德芹 等，2014)认为近30年青藏高原植被的返青期具有提前的趋势，并且伴随着植被生长季延长，高原陆地生态系统碳积累呈增多的现象；而Piao等人(2011)认为，自1982年至今青藏高原植被返青期并没有明显的变化趋势，而是在不同阶段返青期的提前和推迟有所不同。所以，近几十年来青藏高原植被返青期随着气候变暖的变化情况目前仍存在一定分歧。

青藏高原冻土分布广泛，春季气温升高，土壤融化后植被开始返青，土壤的冻融状态与植被生长季的开始和结束有着十分重要的关系。目前学者的研究都是关注土壤的冻融状态或植被物候各自对长时间气候变化的响应而并未涉及它们之间的相互关系，所以为了进一步研究两者间的内在关系，本文的研究内容如下：(1)采用SMMR、SSMI和SSMIS被动微波亮温数据和GIMMS(Global Inventory Modeling and Mapping Studies)NDVI 3g(Third generation GIMMS NDVI from AVHRR sensors)植被数据集获取青藏高原1982年—2013年长时间序列的地表冻融状态和植被返青期；(2)统计分析青藏高原1982年—2013年地表冻土融化首日和植被返青期的变化；(3)分析地表冻土融化首日与植被返青期的相互关系，探讨两者与青藏高原气温变化的关系及其响应特征。

为了长期观测气候变暖条件下的地表冻融循环变化，本文采用美国国家冰雪数据中心(NSIDC)提供的长时间序列逐日被动微波亮温数据，包括Nimbus-7 SMMR Pathfinder Daily EASE-Grid Brightness Temperatures数据(http://nsidc.org/data/nsidc-0071.html/[2016-09-14])和DMSP SSMI-SSMIS Pathfinder Daily EASE-Grid Brightness Temperatures数据(http://nsidc.org/data/NSIDC-0032.html/[2016-09-14])。虽然L波段微波辐射计也已经在轨运行，但是为了保持长时间序列亮温的微波频率的一致性，本文采用上述传感器的亮度温度作为数据源。其中SMMR亮温数据的时间覆盖范围为1982年—1987年，SSM/I亮温数据的时间覆盖范围为1987年—2009年，SSMIS亮温数据的时间覆盖范围为2009年—2016年，表1描述了各微波辐射计的卫星平台和传感器主要参数。

表 1(Table 1)

表 1 被动微波辐射计的主要参数 Table 1 Main parameters of passive microwave remote sensing 主要参数 传感器 SMMR SSM/I SSMIS 平台 NIMBUS-7 DMSP-F08 DMSP-F11 DMSP-F13 DMSP-F17 运行时期 1978-10-26— 1987-8-20 1987-9-7—1991-12-30 1991-12-3—1995-9-30 1995-5-3—2009-4-29 2006-12-14至今 频率/GHz 6.6，10.7，18.0，37.0 19.35，22.24，37.05，85.50 19.35，22.24，37.05，85.50 19.35，22.24，37.05，85.50 19.35，22.24，37.05，91.66 极化方式 V, H V, H V, H V, H V, H 空间分辨率/km 25×25 25×25 25×25 25×25 25×25 入射角/(°) 50.2 53.1 53.1 53.1 53.1 重访周期/d 2 1 1 1 1 扫描宽度/km 780 1400 1400 1400 1400 过赤道时间(UTC) 12:00 (A) 24:00 (D) 6:20 (A) 18:20 (D) 17:17 (A) 5:17 (D) 17:58 (A) 5:58 (D) 17:31 (A) 5:31 (D)

表 1 被动微波辐射计的主要参数 Table 1 Main parameters of passive microwave remote sensing

为保证和微波亮温数据集的空间分辨率和时间覆盖的尽量匹配，本文采用了全球监测与模型研究组的GIMMS NDVI 3g数据集(http://glcf.umd.edu/data/gimms/[2016-11-20])。该数据集是由美国国家航天航空局(NASA)发布的全球植被指数产品，均已经过辐射校正、几何校正和图像增强等影像预处理，空间分辨率为8 km，时间分辨率为15日，覆盖自1982-01—2013-12，可以为中高纬度地区的物候研究提供长时间序列的植被指数数据，在长时间序列叶面积指数(LAI)、总初级生产力(GPP)、光合有效辐射吸收比例(FPAR)等遥感数据产品的制备、物候信息的提取、干旱区荒漠化分析中发挥了重要的作用(Miao 等，2015；Wang 等，2014；Zhu 等，2013)。

为了保持地表冻融和返青期覆盖的时间一致，选取1982年—2013年为本文的研究期。

在双指标算法中，高程是决定地表冻融状态分类阈值的重要因子之一，本研究DEM数据是由西部环境与生态科学数据中心(http://westdc.westgis.ac.cn/[2016-11-25])提供的500 m分辨率DEM数据EASE-GRID投影转换和重采样而成。中国地区综合土地覆盖数据集(冉有华 等，2010)是对中国2000年1∶10万土地覆盖数据进行合并和矢栅转换得到的土地利用数据产品，本研究将其作为双指标算法中针对不同地表类型确定冻融分类阈值的土地分类标准。中国1∶100万植被数据集(侯学煜，2001)是由中国科学院等单位编制出版，由于青藏高原东南部的低海拔地区存在覆盖度较高的常绿乔木林，难以通过NDVI的时间序列提取其返青期，所以利用植被数据集将这部分地区剔除。中国气象局(CMA)在青藏高原布设的90个气象站点，可提供1982年—2006年青藏高原地表0—160 cm的土壤温度逐日观测数据(图1)，用于验证双指标算法的地表冻融状态识别精度，并用于分析植被返青时的日均地表0 cm土壤温度状况。

图 1 青藏高原CMA气象站点分布图 Figure 1 Distribution of CMA meteorology stations on the Qinghai-Tibet Plateau

双指标算法选取卫星降轨37 GHz垂直极化亮度温度( ${T_{{\rm{b}}\_37{\rm{V}}}}$ )和降轨37 GHz与降轨18/19 GHz垂直极化的负亮温谱梯度( $S{G_{\rm{d}}}$ )两个指标来对地表冻融状态进行识别。降轨37 GHz垂直极化亮度温度与地表温度和气温具有高度的相关性，而其发射频率对土壤中的水分含量并不敏感；冻结土壤由于体散射增强而导致亮温“变暗”，即在微波的高频由体散射引起的散射消光比低频强，因而呈现负亮温梯度。所以当双指标均小于各自阈值时，地表为冻结土壤，否则为融化土壤：

式中， ${T_{{\rm{d}}\_{\rm{cutoff}}}}$ 是降轨37 GHz垂直极化亮度温度的阈值， $S{G_{{\rm{d\_cutoff}}}}$ 是降轨37 GHz与降轨18/19 GHz垂直极化亮度温度光谱差的阈值。

由于在传统的双指标算法中没有考虑不同地表类型对亮度温度的影响， ${T_{{\rm{b}}\_37{\rm{V}}}}$ 采用统一的亮温阈值，这会导致地表冻融状态的判别产生一定的误差。因此，本文采用Jin等人(2015)改进后的双指标算法，即根据中国气象局台站的日最低地表温度，针对不同地表类型的 ${T_{{\rm{b}}\_37{\rm{V}}}}$ 阈值分别进行了标定(表2)，标定阈值的地表类型包括水浇地、旱地、有林地、疏林地、灌木林、高覆盖草地、中覆盖草地、低覆盖草地、高原高覆盖草地(h>3500 m)、高原中覆盖草地(h>3500 m)、高原低覆盖草地(h>3500 m)及未利用地(Jin 等，2015)。

表 2(Table 2)

表 2 双指标分类算法中Tb-37 V阈值(Jin 等，2015) Table 2 Thresholds of Tb_37 V in the dual-indices classification algorithm 地表类型 SMMR-Tb_37 V/K SSM/I-Tb_37 V/K 水浇地 258.39 260.47 旱地 259.50 262.82 有林地 260.42 262.24 疏林地 258.16 261.16 灌木林 258.57 260.53 高覆盖草地 250.51 260.94 中覆盖草地 258.38 260.07 低覆盖草地 257.95 259.54 高原高覆盖草地(h>3500 m) 258.74 260.24 高原中覆盖草地(h>3500 m) 254.72 255.51 高原低覆盖草地(h>3500 m) 252.49 254.90 未利用地 251.91 254.57

表 2 双指标分类算法中Tb-37 V阈值(Jin 等，2015) Table 2 Thresholds of Tb_37 V in the dual-indices classification algorithm

在识别逐日地表冻融状态后，以当年7月1日至次年6月30日为一个完整年度，将地表首次为冻结状态的日期作为冻结首日，地表最后一次为冻结状态的日期为冻结终日，其次日为融化首日，而冻结首日和冻结终日之间的天数定义为冻结期。

由于植被生长过程的渐变特征，NDVI的时间序列曲线理论上应该非常平滑(De Jong 等，2011)，但是在NDVI遥感数据产品中，由于时间分辨率和云层、气溶胶或冰雪覆盖等影响而使NDVI值出现异常，因此需要对植被指数遥感产品进行滤波重建(龙鑫 等，2013)。时间序列植被指数数据的拟合去噪方法包含两类：曲线拟合法和基于滤波函数的拟合方法。前者主要包括非对称性高斯函数拟合法(Jonsson和Eklundh，2002)和双Logistic函数拟合法(Beck 等，2006)；后者主要有Savitzky-Golay滤波法(Chen 等，2004)，NDVI时间序列谐波分析法(Jakubauskas 等，2001)，均值迭代滤波法(Ma和Veroustraete，2006)和变权重滤波方法(Zhu 等，2012)等方法。

本文采用Savitzky-Golay滤波法对GIMMS NDVI 3g时序数据进行滤波重建，该方法已成为目前植被指数时序数据平滑降噪的主要方法之一。这种方法在时域内采用局域多项式最小二乘法拟合来进行滤波，与传统的选取常数拟合的方法不同，它使用一个高阶多项式来实现最小二乘拟合滑动窗口近似的功能。相比于其他的滤波方法，它的优势在于既能够降低噪声又能保证曲线的形状不会发生改变(Luo 等，2005)，同时能够更好地逼近曲线的上包络线，捕获到极值周围的波动细节，有利于更加准确地识别植被生长的开始和结束等物候信息。同时，对于沙漠、农作物、森林等不同地表类型的NDVI在时间序列中不会产生极端值(Song 等，2011)。Savitzky-Golay滤波法的一般方程如下：

式中， ${Y_{j + i}}$ 为NDVI的原始值； $Y_j^*$ 为NDVI的拟合值； ${C_i}$ 为滤波窗口中第j个NDVI值的权重系数；m为滤波窗口一半的宽度；N为卷积数，等于滤波窗口的大小(2m+1)。根据宋春桥等人(2011)对青藏高原时序NDVI重建时滤波窗口调试结果，将滤波窗口的大小设置为5。图2展示了青藏高原内某像元2001年—2002年S-G滤波前后对比图。

图 2 青藏高原内某像元2001年—2002年S-G滤波前后对比图 Figure 2 Comparison between original NDVI and filtering by S-G filter in one pixel of the Qinghai-Tibet Plateau from 2001 to 2002

提取返青期的常用方法有阈值法(固定阈值法，最大斜率阈值法，动态阈值法)、曲率法、主成分分析法等，不同方法对不同研究区和植被类型的适用性不同，要根据研究的区域和植被的特征选取最佳提取方法。本文根据目前青藏高原返青期的研究结果(丛楠和沈妙根，2016；Piao 等，2011；范德芹 等，2014)，采用动态阈值法进行返青期的提取。

阈值法是指先设定一个阈值，当时域内某时刻的NDVI达到该阈值时，这个时刻被定义为植被的返青期。动态阈值法是指当某时刻NDVI与时域内最小值的差值达到了时域内最大值与最小值差值的特定比例时，该时刻被定义为植被的返青期(Richardson 等，2006)。相比固定阈值法和最大斜率阈值法，动态阈值法消除了背景值的差异对返青期提取的影响。动态阈值法的一般表达式如下：

青藏高原面积较大，不同地区海拔差异大，地形复杂，因此本文借鉴Liu(2014)在研究青藏高原植被物候变化时采用的方法，选取当春季 ${\rm{NDV}}{{\rm{I}}_{{\rm{thd}}}}$ 的值达到0.3时为青藏高原植被返青期的开始(刘双俞 等，2014)。同时，为了避免在植被稀少地区的土壤变化引起的返青期识别不准问题，本文只对植被生长季4—10月平均NDVI > 0.2的地区进行研究( Zhou 等，2001)。

图3为该实验数据处理与分析的技术路线图。

图 3 数据处理与分析流程 Figure 3 Flow chart of data process and analysis

根据被动微波亮温数据和双指标算法，获得青藏高原1982年—2013年的逐日地表冻融状态，分类结果以日为单元进行文件存储，文件类型为ASCIIGRID，投影方式为等积割圆柱投影，双标准纬线为南北纬30°。采用1982年—2006年90个分布于青藏高原地区的中国气象局站点的日最低地表0 cm土壤温度观测结果进行地面验证，结果表明1982年—2006年地表冻结状态的分类精度为85.06%，融化状态的分类精度为78.49%，总体分类精度为82.15%。图4以2001年7月—2002年6月一个完整年度为例，展示了青藏高原逐月的地表冻结状态的分布范围和发展变化规律，地表冻结最大和最小面积分别为2.44×10⁶ km²和8.01×10⁵ km²，分别占青藏高原面积的94.19%和30.86%。青藏高原西北地区存在一部分终年不融的永久冻土，每年8月开始，冻土的分布从西北向东南扩展，至次年1月达到最大，至此只有高原东南的墨脱县及其周围地区终年不冻，该地区平均海拔约1200 m，年均温约16 ℃，年降水量约2500 mm，植被以热带雨林为主。

图 4 青藏高原2001年7月—2002年6月逐月地表冻融状态范围 Figure 4 Monthly distribution of surface freeze/thaw status on the Qinghai-Tibet Plateau from July, 2001 to June, 2002

根据1982年—2013年被动微波亮温识别的长时间序列地表冻融状态结果，在空间上，青藏高原平均地表冻结时间集中在10月28日至次年3月7日，地表融化首日集中在第132—147日(5月12—27日)，融化首日最早发生在2007年，为第128日(5月8日)；最晚发生在1989年，为第154日(6月3日)。图5展示了青藏高原1982年—2013年平均地表融化首日的空间分布情况。青藏高原的地形和冻土分布复杂，所以融化首日的分布也很复杂：整体上，它受海拔和纬度的影响，东部地区的冻土融化早于西部，南部则早于北部。在青海东部和四川北部，地表冻土融化首日早于第130日(5月10日)；青海中部和西部，地表融化首日通常在第131—150日(5月11—30日)；青海西部以及西藏大部分地区，地表融化首日晚于第151日(6月1日)。图中空白区域为生长季4—10月平均NDVI < 0.2的植被或常绿乔木林植被覆盖的地区。

图 5 青藏高原1982年—2013年地表冻土融化首日的空间分布 Figure 5 Distribution of surface frozen soil thawing date over the Qinghai-Tibet Plateau from 1982 to 2013

根据1982年—2013年青藏高原逐年平均地表融化首日来看，32年间，地表融化首日具有一定的波动性，但在波动中存在一定的变化趋势：1982年—1987年，地表融化首日呈现推迟趋势(1.93±1.81 d/a, R=0.47)；而在之后的长达20年之中呈现明显的提前趋势(–0.67±0.20 d/a, R=0.63)，这与青藏高原过去几十年内表现出的气候变暖、冻土退化、活动层增加等现象一致；在2006年—2013年，融化首日又有一定的推迟(0.97±0.84 d/a, R=0.42)。

经过S-G滤波法的植被指数时间序列重建和动态阈值法的返青期提取，得到1982年—2013年青藏高原地表植被的返青期。结果表明，在空间上，32年间青藏高原植被平均返青期集中在每年的第139—149日(5月19—29日)，平均返青期最早年份在1998年，为第135日(5月15日)，最晚年份在2012年，为第152日(6月1日)。图6展示了青藏高原1998年—2013年植被平均返青期的空间分布情况，与冻结地表融化首日发生的时间相似，由于海拔和纬度地带性而导致的植被类型的差异，青藏高原西部植被返青期比东部普遍偏晚；高原东部和东南部气候较为温暖湿润，植被以山地常绿或落叶阔叶林、热带雨林，山地针叶林、高山草甸为主，植被返青期早于第135日(5月15日)；中部以半干旱的高山草原和山地灌丛为主，植被返青主要发生在第135—145日(5月15—25日)；西部则是高寒半荒漠和荒漠，土壤裸露，植被稀少，返青期晚于第145日(5月25日)，生长期短。图6中高原东南部存在一些植被返青期晚于第155日的区域是由于一些其他的常绿或常绿落叶混交植被所造成的。

图 6 青藏高原1982年—2013年植被返青期的空间分布 Figure 6 Distribution of vegetation green-up date over the Qinghai-Tibet Plateau from 1982 to 2013

青藏高原1982年—2013年植被返青期的变化与地表融化首日的分段性变化趋势相似，在1982年—1987年，返青期推迟(0.28±1.01 d/a, R=0.14)；而在1987年—2006年的20年间，返青期表现出轻微的提前趋势，但其平均速度(–0.13±0.16 d/a, R=0.20)远小于融化首日的提前；在2006年—2013年，返青期显著推迟(1.04±0.52 d/a, R=0.77)。

根据被动微波和NDVI对青藏高原地表融化首日和返青期的提取，1982年—2013年青藏高原植被返青期与融化首日具有显著的相关性(R=0.51，P=0.003)：在1982年—1987年，融化首日和返青期均有推迟(R=0.92，P=0.01)，而在1987年—2006年均出现提前(R=0.59，P=0.006)，2006年—2013年间又均呈现推迟(R=0.64，P=0.086)。1982年—2013年植被返青期通常发生在地表融化首日的3天后(3.94±5.58日)。图7和图8分别展示了青藏高原1982年—2013年地表融化首日和植被返青期变化趋势和相关性。整体上来看，在两者呈推迟趋势时相关性高于两者提前时，这说明返青期的推迟主要是由于地表冻土融化首日推迟导致，而返青期的提前除了融化首日提前外，温度等其他因素也产生了较大的影响。

图 7 青藏高原1982年—2013年地表融化首日与植被返青期 Figure 7 Thawing time and vegetation green-up date over the Qinghai-Tibet Plateau from 1982 to 2013

图 8 青藏高原1982年—2013年植被返青期与地表融化首日的关系 Figure 8 Relation between vegetation green-up date and frozen soil thawing time over the Qinghai-Tibet Plateau from 1982 to 2013

实际上，植被根区土壤的融化本身是其返青的必要条件。土壤融化提供水分，植被根部光合作用与呼吸作用加强，净初级生产量开始积累，植被返青，因此植被返青与地表土壤融化应当具有一定的关系。从本质上来说，土壤融化和植被返青期均受到高原春季气温和地表温度的控制，当温度达到一定条件，两种现象才有可能发生。通过青藏高原气象站的日均地表0 cm土壤温度数据(1982年—2006年)的统计结果来看，土壤开始融化时，地表日平均温度为13.64℃；植被返青时，地表日平均温度为13.84℃，即1982年—2006年植被返青时日均地表0 cm的日平均温度约比地表土壤开始融化时高0.2℃，说明返青期比地表冻土融化的温度阈值要高，所以春季植被返青期普遍晚于地表土壤的融化时间。

但温度并不是决定地表土壤融化和植被返青的唯一因素，降水也是重要的影响因子，它能够影响冻土和植被对温度升高的响应。根据青藏高原部分气象站点的年均春季地表0 cm土壤温度(3—6月)显示，1982年—2006年间气象站点附近的年均春季地表0 cm土壤温度呈显著的上升趋势(图9)，这说明青藏高原的气候是变暖的。虽然地表融化首日和植被返青期大部分时段内呈现出提前的趋势，但少部分时段也存在推迟现象。根据Zhang等人(2015)的研究结果，1982年以来青藏高原不同地区降水的变化趋势不同，温度升高的情况下标准化降水蒸发指数(SPEI)、地面气温梯度以及日照时间均存在差异。总体而言，在温度和降水均升高的情况下，青藏高原相对较为干旱的地区植被返青提前，相对湿润的地区则会延迟；反之在温度升高而降水不变或者减小的情况下，相对干旱的地区植被的返青期推迟，相对湿润的地区则提前。换言之，降水能够影响植被对温度的敏感性。同样，地表冻融状态也会受到降水的影响，降水能够起到绝热作用以减弱气温对地表融化的驱动(Wang 等，2015)。

根据地表融化首日和返青期的结果分析，还发现1982年—2013年返青期的变化速度比融化首日慢。例如，1982年—1987年融化首日的平均变化速度为1.93 d/a，是返青期的7倍；1987年—2006年融化首日的平均变化速度为–0.67 d/a，是返青期的5倍。冻土(非生物)的融化是对高原气候变化的直接响应，而植被(生物)对气候变化的响应是相对间接的。对于这一现象，有Piao等人(2011)认为，青藏高原植被返青期的温度阈值随着气候的变化不是一成不变的，在长时间气候变化过程中，植被返青期的温度阈值会因生物自身的调节机制，而随着气候变暖而提高；而地表土壤融化首日则是与外界环境的温度直接相关的。根据1982年—2006年青藏高原气象站的春季地表0 cm土壤温度与植被返青温度阈值的统计结果得到与Piao等人(2011)一致的观点(图10)，植被返青期的温度阈值与春季平均地表0 cm土壤温度有着显著的正相关关系(R=0.72，P<0.001)，即植被返青期的温度阈值会随着气候变化而变化。

除此之外为了探究植被返青期的提前是否会影响生长季的长度以及植被生长季长度在时间序列上的变化，利用GLIMMS NDVI 3g的S-G滤波重建结果，提取了1982年—2013年青藏高原植被生长季，结果显示青藏高原植被的生长季长度与返青期没有明显的相关关系(R=0.28，P=0.124)，且32年间生长季的长度也没有明显的变化趋势(R=0.33，P=0.062)，即未发现返青期的变化会导致生长季的延长或缩短。因此，在长时间气候变化过程中，青藏高原植被很可能会对气候变化做出适应性调节。

也有类似研究表明，青藏高原气候的普遍变暖并不一定会导致返青期显著提前，如有学者指出，青藏高原植被返青期和夏季(6—8月)季前夜间温度具有明显的负相关关系，这说明这段时间内夜间变暖确实有利于植被返青期提前；但同时又指出返青期与夏季季前白天温度的相关性相对较弱，而且不同地区受温度的影响也并不相同：白天升温对返青期的提前作用主要出现在高原相对湿润的地区，而在干旱区，白天升温对返青期呈现抑制作用(Shen 等，2016)，这也能够解释本文前述的“地表融化首日提前的速度大于植被返青期提前的速度”和“气候变暖的条件下返青期没有持续提前”的现象。

图 9 1982年—2006年青藏高原气象站春季地表温度 Figure 9 Spring surface temperature from the CMA meteorology sites on the Qinghai-Tibet Plateau from 1982 to 2006

图 10 1982年—2006年青藏高原气象站春季地表温度与返青期阈值的关系 Figure 10 Relation between air temperature threshold of vegetation green-up and spring surface temperature across the CMA meteorology sites in the Qinghai-Tibet Plateau from 1982 to 2006

本文通过双指标算法对全球长时间序列被动微波亮度温度数据集进行处理，得到了青藏高原地区年份的地表冻融状态数据集；通过建立一套合理的滤波重建和返青期提取的方法，根据全球长时间序列植被指数数据集提取出青藏高原地区年份的逐年返青期信息；利用以上两个遥感数据产品，分析了青藏高原地表冻融状态和返青期的时空分布特征，研究了两者的变化趋势、相关关系以及它们对温度变化的响应。从中发现，青藏高原平均地表冻结时间集中在10月28日至次年3月7日，地表融化首日集中在5月12—27日，植被平均返青期集中在每年的5月19—29日；青藏高原的地表融化首日和返青期具有相似的分段性变化趋势：即在1982年—1987年，地表融化首日和返青期均有推迟，而在1987年—2006均出现提前，2006年—2013年间又均呈现推迟。根据统计结果表明，1982年—2013年植被返青期通常发生在地表融化首日的3天后(3.94±5.58日)，并且两者具有显著的相关关系；但是在青藏高原持续升温的背景下，两者并未持续提前。因为不论是地表融化还是植被返青，不仅仅受温度的影响，还受降水等其他环境因素的控制；同时我们还发现了植被对于温度升高具有一定的自我调节机制，它可能通过提高返青的温度阈值和控制生长季的长度来对气候变暖带来的结果做出负反馈。

在本研究中使用了布设在青藏高原的超过80个CMA气象站点的日地表0 cm土壤温度观测数据，但是站点的分布海拔和分布密度仍然不能直接代表青藏高原温度的整体变化趋势。同时，由于GIMMS的时间分辨率和数据质量的原因，返青期的提取不可避免存在一定的误差，后续研究可以利用通量塔等进行观测和验证。另外，为了避免低植被覆盖和常绿植被覆盖造成返青期提取的误差，本研究去除了NDVI极低(NDVI在4—10月的均值小于0.2)和极高(植被覆盖为常绿乔木林)的区域，前者多为西部高海拔地区，对气候变化响应最为敏感，这也会在一定程度上影响结果。