植被物候是植被整个生长发育过程的周期性现象，通过发芽、展叶、开花、落叶等关键节点的变化，反映对气候变化的响应。植被物候在全球生态系统中扮演关键角色，其动态变化与气候因子密切相关，因此，物候变化特征及其对气候变化的响应机制，常被视为探究全球生态系统演变规律的参考依据[1]。近年来，极端气候事件的发生频次、强度等均有所增加，而东北地区是极端气候事件变化较明显的地区，其植被类型丰富，对气候变化较敏感[2-3]。因此，对东北三省的植被物候变化进行监测，并探究其对极端气候事件的响应机制，对生态系统的科学研究、环境资源的管理与保护具有重要意义。
目前，利用遥感手段获取的数据资源被广泛用于大尺度长时间序列的植被物候监测，有效弥补了人工观测提取物候指标的局限性[4]。屈鹏等[5]利用MODIS数据提取了青藏高原草地物候参数；缪丽娟等[6]利用MODIS NDVI数据提取了中国城市植被物候参数，二者均表明了通过遥感手段获取的物候指标具有一定的可靠性。随着极端气候事件的频发，针对东北地区极端气候事件的相关研究，逐渐成为学界广泛关注的热点领域。刘俊男等[7]研究了东北地区极端降水的时空变化特征；辛玉洁等[8]研究了东北地区极端气温的时空演变特征。当前关于东北地区植被物候对气候变化的响应研究多集中在植被物候和平均气候间的关系，对于植被春季物候（Start of the growing season，SOS）对极端温度事件的响应研究较少，而与平均气候相比，极端气候对植被物候的影响似乎更为明显[9]。因此，为充分认识气候变化事件对东北三省植被的影响，探究植被春季物候对极端温度事件的响应机制至关重要。
本研究以植被SOS为研究植被响应极端气候事件的重要物候指标，基于2001—2020年研究区MODIS NDVI数据，运用动态阈值法提取物候指标，利用ERA5-Land再分析数据集提取极端温度事件指标，研究植被春季物候的变化特征，通过Pearson相关分析法对植被春季物候与极端温度事件的关系进行分析，为探析东北三省植被物候对极端温度事件的响应机制提供参考。
1 材料与方法
1.1 研究区基本情况
研究区属温带湿润、半湿润大陆性季风气候区，地貌类型主要有平原、丘陵和山地，土地利用类型以耕地、林地、草地、水域和未利用地为主[10]。
1.2 数据来源与预处理
植被指数数据来自MODIS的MOD13A2产品中的归一化植被指数（Normalized difference vegetation index，NDVI），空间分辨率为1 km，时间分辨率为16 d，数据时间覆盖2001—2020年。气象数据来自ERA5-Land Daily Aggregated数据集，选取日最高气温、日最低气温以及日降水量作为气象因子，空间分辨率为0.1°×0.1°，时间分辨率为1 d，数据时间覆盖1981—2020年。土地利用数据来自MODIS的MCD12C1 v061产品，空间分辨率为0.05°×0.05°，时间分辨率为1 a，数据时间覆盖2001—2020年。
为消除人为因素和植被类型自然变化的影响，利用MCD12C1数据进行逐年土地利用类型分析，筛除粮食、果树等栽培植被类型的分布区，仅保留研究时间内植被类型未发生变化的自然植被覆盖区域数据，并以此作为研究区植被覆盖区域[11]。利用ArcGIS 10.8软件基于植被分布矢量图层对上述数据进行掩膜提取，并通过重采样将数据空间分辨率统一至1 km。
1.3 研究方法
1.3.1 植被物候参数的提取
采用Savitzky-Golay滤波拟合法对研究区2001—2020年MODIS NDVI数据逐像元进行时间序列曲线拟合，再利用动态阈值法逐像元提取植被物候参数，其计算如式（1） 。
1.3.2 极端气候事件提取
表1 极端气候事件指标 |
| 类型 | 指标名称 | 代码 | 定义 | 单位 |
|---|---|---|---|---|
| 高温 指标 | 重度高温天数 | TX90P | 每月日最高温>90%分位值的天数 | d |
| 日最高气温极高值 | TXx | 每月日最高温的最大值 | ℃ | |
| 重度高温最长持续天数 | MaxDurTX90P | 每月日最高温>90%分位值的最长连续天数 | d | |
| 低温 指标 | 重度低温天数 | TN10P | 每月日最低温<10%分位值的天数 | d |
| 日最低气温极低值 | TNn | 每月日最高温的最小值 | ℃ | |
| 重度低温最长持续天数 | MaxDurTN10P | 每月日最低温<10%分位值的最长持续天数 | d | |
1.3.3 Theil-Sen斜率估计与Mann-Kendall检验
采用Theil-Sen斜率估计计算植被SOS变化趋势，可去除异常值的影响，斜率计算如式（2） [3]。
采用Mann-Kendall检验对结果显著性进行检验以判断其显著变化程度，统计量 计算如式（3） 。
为正态分布，方差计算如式（5） 。
标准化检验统计量Z计算如式（6） 。
设置显著性水平 ，若 ，则表示时间序列数据发生了显著变化。结合Theil-Sen斜率估计和Mann-Kendall检验方法获取发生显著变化的区域， 且 表示数据显著升高； 且 表示数据显著下降。
1.3.4 相关性分析
2 结果与分析
2.1 植被SOS的变化趋势
由图1可知，2001—2020年，植被SOS像元均值多数在一年中的第100～130天，约占研究区面积的80%以上，其中分布在第110～120天之间的区域占比最大，为45.62%，其次为分布在第100～110天之间的植被SOS区域，其占整体像元的27.25%，分布在第120～130天之间的区域较少，仅占9.84%，其余均值区间占比均不超过5%。由图2可知，85.19%的区域呈现物候提前的变化趋势，其中大多分布在-0.5～0 d/a，占研究区面积的50.32%，其次分布在-1.0～-0.5 d/a，占研究区面积的25.42%。Theil-Sen斜率估计与Mann-Kendall检验表明，物候呈显著提前趋势的区域占整体像元的12.97%。由图3可知，研究区植被SOS整体以-0.25 d/a的速率提前（P>0.05）。综上，研究区在2001—2020年植被SOS主要呈提前变化趋势。
2.2 植被SOS对极端温度事件的响应
由表2可知，研究区植被SOS与极端高温指标之间以负相关关系为主，3—5月TX90P与植被SOS呈负相关的区域占比先增加后减少，其中通过显著性检验的像元占整体像元比例分别为22.65%、27.42%和2.90%；而植被SOS与TX90P呈现正相关的区域占比较少，且呈显著正相关的像元比例在1%左右。上述数据说明，发生在3月、4月TX90P事件对植被SOS的影响更为明显，且使大部分区域植被SOS呈提前趋势。3—5月植被SOS与TXx存在相关关系的区域占比和变化趋势与TX90P相似，其中呈现显著正相关关系的区域不足1%，而通过显著性检验的负相关区域占比分别为14.98%、18.85%和4.29%，3月、4月与TX90P相比比例减少，说明3月、4月植被SOS对TX90P的响应更加显著。3—5月植被SOS与MaxDurTX90P呈负相关的像元比例逐渐减少，其中呈显著负相关关系的像元占比分别为25.08%、13.38%和2.76%，说明3月、4月MaxDurTX90P对植被SOS的影响较大，而呈正相关的区域虽然有所提升，但其中表现为显著正相关的像元占比仅为1%，可见MaxDurTX90P与植被SOS主要存在负相关关系。整体而言，3月、4月植被SOS与极端高温指标的负相关关系更为显著，极端高温事件增多在一定程度上会促使研究区植被SOS提前。
表2 研究区植被SOS与极端高温指标相关性检验统计 单位：% |
| 代码 | 月份 | 相关性 | |||
|---|---|---|---|---|---|
| 正相关 | 显著正相关 | 负相关 | 显著负相关 | ||
| TX90P | 3 | 21.36 | 1.14 | 78.64 | 22.65 |
| 4 | 18.40 | 0.36 | 81.60 | 27.42 | |
| 5 | 39.58 | 0.47 | 60.42 | 2.90 | |
| TXx | 3 | 23.84 | 0.43 | 76.16 | 14.98 |
| 4 | 19.58 | 0.43 | 80.42 | 18.85 | |
| 5 | 36.37 | 0.74 | 63.63 | 4.29 | |
| MaxDurTX90P | 3 | 16.99 | 0.18 | 83.01 | 25.08 |
| 4 | 26.51 | 0.42 | 73.49 | 13.38 | |
| 5 | 46.09 | 1.15 | 53.91 | 2.76 | |
由表3可知，植被SOS与极端低温指数相关关系较复杂，各指标在不同月份表现不同，且存在正负相关关系的区域占比差异较小。3月植被SOS与TN10P的正负相关关系像元占比相近，在4月、5月呈正相关的区域比例更大，且呈显著相关的像元比例略高，说明4月、5月TN10P对植被SOS的影响较大，在一定程度上TN10P的增加可能会使植被SOS推迟。植被SOS与TNn的相关系数在3月、5月正负比例接近，而在4月呈负相关的区域较多，其通过显著性检验的负相关区域像元占比达10.03%，说明植被SOS对4月的TNn响应最为明显。植被SOS与MaxDurTN10P的相关关系比例基本相同，具有显著相关关系的区域大多在5%以下，且呈显著负相关的比例略高于呈显著正相关比例，说明植被SOS对MaxDurTN10P的响应并不显著。整体而言，植被SOS对极端低温事件的响应并不明显。
表3 研究区植被SOS与极端低温指标相关性检验统计 单位：% |
| 代码 | 月份 | 相关性 | |||
|---|---|---|---|---|---|
| 正相关 | 显著正相关 | 负相关 | 显著负相关 | ||
| TN10P | 3 | 49.48 | 1.92 | 50.52 | 5.34 |
| 4 | 70.25 | 8.99 | 29.75 | 1.87 | |
| 5 | 61.76 | 5.96 | 38.24 | 2.20 | |
| TNn | 3 | 55.90 | 4.38 | 44.09 | 2.03 |
| 4 | 29.98 | 1.63 | 70.02 | 10.03 | |
| 5 | 42.40 | 1.82 | 57.60 | 2.06 | |
| MaxDurTN10P | 3 | 44.91 | 2.02 | 55.09 | 5.94 |
| 4 | 54.51 | 1.46 | 45.49 | 2.20 | |
| 5 | 51.59 | 2.88 | 48.41 | 3.36 | |
3 结论
植被通过调节自身的形态结构和生长来适应所处的生长环境，从而最大化地利用环境资源。然而，当极端气候事件发生时，植被对环境的适应平衡将会被打破，这间接影响了植被的发育机制、生长态势以及分布格局[14]。本研究基于2001—2020年MODIS NDVI数据提取了研究区植被春季物候，分析多年变化特征，并结合ERA5-Land Daily Aggregated数据集中的日最高气温和日最低气温数据，利用Pearson相关分析法探究了研究区植被春季物候对极端温度事件的响应特征。具体结论如下。
（1）研究区2001—2020年植被SOS多集中在每年的第100～130天，其面积占研究区面积的80%以上，其中分布在第110～120天之间的区域占比最大。
（2）2001—2020年，植被SOS整体以提前的变化趋势为主，在像元尺度上，变化速率集中在-1.5～0.5 d/a，85.2%的区域表现为植被SOS提前，12.97%的植被SOS为显著提前变化。
（3）植被SOS与极端高温指标主要呈负相关，其对3月、4月极端高温指标的响应更为显著，与极端低温指标的相关关系较为模糊且对其响应并不明显。