基于MODIS EVI的大兴安岭多年冻土区植被物候研究

王雪梅， 温理想， 李佳诺， 郭 蒙

（1. 东北师范大学 地理科学学院 长白山地理过程与生态安全教育部重点实验室，吉林 长春 130024；2. 北京师范大学 环境学院 水环境模拟国家重点实验室，北京 100875）

0 引言

植被物候是植物生命周期中季节变化的节律，由于其对气候变化敏感以及对生态系统功能的重要性，在过去几十年中受到国内外学者们的广泛关注［1］。IPCC 第六次评估报告［2］指出，自1850 至1900年以来，全球地表平均温度已上升约1 ℃，未来20年全球温升将达到或超过1.5 ℃。气候变暖能够引起植物发芽、展叶、开花、落叶和叶颜色等特征的变化，直接影响植被生产力、生物多样性和全球碳循环等［3-4］。作为冰冻圈的重要组成部分，多年冻土对气候变化极为敏感。准确监测多年冻土区植被物候特征对于研究寒区生态系统对气候变化的响应具有重要意义。

传统的物候观测方法是以物候观测网络为主的地面观测，这种基于植物个体尺度的观测难以满足生态系统、区域乃至全球尺度的物候研究［5］。随着遥感技术的发展，基于植被指数（VIs）的监测方法已广泛应用于植被物候研究［6-8］。目前常用的植被指数包括归一化植被指数（NDVI）和增强型植被指数（EVI）等［9］。然而，地面阴影、雪和云等背景的干扰以及恶劣的大气观测条件等因素都会降低NDVI和EVI 监测植被物候的准确性。例如，Hmimina等［10］使用中分辨率成像光谱仪（MODIS）提取的NDVI监测了不同类型植被的季节变化，发现NDVI无法准确反演常绿森林的物候模式。此外，落叶林中乔木通常比草本和灌木的返青时间晚，这意味着由植被指数确定的生长季开始日期可能反映的是草本和灌木，而非优势树种的物候信息［11］。作为植被指数的补充，太阳诱导叶绿素荧光（SIF）为卫星监测植物功能提供了新的可能性［12-13］。SIF 是光合作用的副产品，包含代谢、生化等与植被生长密切相关的信息，与传统的植被指数相比，SIF 具有更强的植物生理基础［14］。

植被指数的时间序列可以表征季节或年际尺度上的植物生长发育的周期变化［15］。然而受气溶胶、水汽、传感器退化等因素的影响，卫星遥感获取的VIs 时间序列含有明显的数据噪声，难以进行物候信息提取和趋势分析［16］，因此需要采用去噪和平滑等处理来重构时序数据。常用的重构方法包括Savitzky-Golay（S-G）滤波法［17］、双逻辑斯蒂函数（Double-logistic）［18］、非对称高斯函数（Asymmetric Gaussian）拟合法［19］、傅里叶变换（Fourier transformation）［20-21］和小波变换（Wavelet transformation）［22］等。预处理后的VIs 时序曲线采用阈值法、拐点法和曲线斜率法等算法提取植被物候参数［23］。目前遥感时序数据预处理的方法较多，通常根据研究区植被生长特点和数据源质量采用最适用的方法［24］。

相关研究表明，随着全球温度的升高，北半球高纬度地区的植被物候显示出春季提前和秋季延迟的趋势［25-26］。然而，由于温度、降水和土壤湿度等非生物因素的影响，植被物候的变化呈现出复杂的时空格局［27］。大兴安岭多年冻土区位于欧亚大陆多年冻土区的南部边缘，是中国唯一的高纬度多年冻土区，对全球气候变化高度敏感［28-29］。为了探究气候变化对大兴安岭多年冻土区植被的影响，本文基于MODIS EVI 时间序列数据，评估了研究区近20 年植被物候的时空变化特征及其对气候变化的响应，以期丰富寒区生态系统植被物候的研究。

1 数据与方法

1.1 研究区概况

大兴安岭多年冻土区位于我国最北部（49°01′～53°54′ N，119°12′～125°29′ E），海拔在220～1 487 m之间（图1）。该区属于寒温带大陆性季风气候，年平均气温介于-5～2 ℃之间，年平均降水量为460 mm。大兴安岭多年冻土区可分为大片连续多年冻土、岛状融区多年冻土和岛状多年冻土三部分［30-31］，本文以大片连续多年冻土区和岛状融区多年冻土区为研究区，面积约为11.2×104km2，其中大片连续多年冻土区面积约为6.8×104km2，65%～75%的面积分布有多年冻土；岛状融区多年冻土区面积约为4.4×104km2，50%～60%的面积分布有多年冻土。大兴安岭多年冻土区主要林型包括针叶林、阔叶林以及针阔混交林等森林类型［32］，其中兴安落叶松（Larix gmelinii）分布范围最广，属于顶级演替树种。主要的阔叶树种白桦（Betula platyphylla）是林火、砍伐等干扰后的先锋树种，分布范围仅次于落叶松。研究区还以斑块状散布着少量草原、草甸以及沼泽等覆被类型。

图1 研究区地理位置［30］（a）及植被类型［32］（b）Fig. 1 Location and vegetation cover of the study area

1.2 数据来源及预处理

1.2.1 MODIS NDVI和EVI数据

由Terra 和Aqua 两颗卫星携带的MODIS 传感器是目前观测全球生态过程和环境变化的重要数据源。美国国家航空航天局（NASA）的官方网站（https：//search. earthdata. nasa. gov）上发布了多种不同时空尺度的MODIS 产品，本研究使用MODIS13Q1 数据集。MOD13Q1 是使用正弦投影方法的MODIS 三级网格数据产品，已经根据大气和地形的影响对其进行了系统纠正。从MOD13Q1中提取的NDVI 和EVI 具有250 m 的空间分辨率和16 d的时间分辨率。对每期数据执行镶嵌、掩膜和投影变换等操作，得到研究区2000—2019 年NDVI 和EVI的时间序列数据。

1.2.2 SIF数据

OCO-2 是NASA 于2014 年7 月发射的专门用于监测大气中CO2浓度的卫星，重访周期为16 d，赤道过境时间为13：30（当地时间），每天可获取全球尺度的不连续点文件，每个点文件的空间分辨率为1.3×2.25 km2。本研究SIF 数据采用OCO-2 卫星的二级产品OCO-2_L2_Lite_SIF. 8r，它以NetCDF 格式存储，将其转换为Shapefile 格式。为了保证与NDVI、EVI时间分辨率的一致性，按照MODIS 数据合成规则将16 d 的SIF 数据合为一期，每年23 期数据。计算研究区SIF 的平均值，得到2015—2019 年SIF的时间序列数据。

1.2.3 气象数据

气温和降水数据来源于中国气象数据网（http：//data. cma. cn/）提供的大兴安岭多年冻土区及其周围的18个气象站点（图1）的日值数据集。将该数据集进行统计得到月均温和月累积降水量，利用ArcGIS 10.6 软件的克里金插值法进行空间插值得到250 m分辨率的栅格数据。

1.3 研究方法

1.3.1 时间序列重建

原始数据不可避免地存在不同程度的异常值、噪声或空值，这可能导致物候指标提取中的错误［33-34］。因此在物候提取之前，有必要对时间序列数据进行降噪滤波处理。本文采用Savitzky-Golay（S-G）滤波法［35］进行时间序列重建，S-G滤波法是一种基于局部多项式最小二乘法拟合的滤波方法，其计算公式如下：

式中：yj*为滤波后的数据；yj+i代表原始时间序列数据；wi为滤波系数，表示滤波器开始处理的第i个值的权重；m为滤波窗口的大小；N为滑动窗口，其值为2m+1。

1.3.2 物候参数提取

基于TIMESAT 3.3软件，采用动态阈值法来提取SOS、EOS 和LOS 等关键物候指标。由于地表覆被类型、研究方法以及研究区的差异，目前没有统一的阈值标准。软件的开发者Jönsson 等［36］综合多地研究建议将SOS 和EOS 的阈值设定为20%，而Zhao 等［37］在我国东北地区植被物候研究中采用30%阈值取得了很好的效果，同样，Fu 等［38］和Tang等［39］在研究大兴安岭地区植被物候时也将阈值设为30%。参考已有研究，本文使用30%作为植被物候提取的阈值。结果采用年序日（DOY）的形式表示，即从每年第一天算起的实际日数。

1.3.3 趋势分析

Mann-Kendall（MK）趋势检验法和Sen 斜率法的结合已广泛用于植被变化趋势的分析，是长时间序列数据分析的重要方法［40］。MK 检验法可以有效地检测时间序列数据的变化趋势，而无需遵循特定分布的样本，也不会受到少数异常值的干扰［41］，Sen斜率估计（Theil-Sen Median）是由Sen 提出和开发的一种非参数统计的趋势计算方法［42］，计算公式如下：

式中：Q为Sen 的斜率；Xj和Xi分别为时间j和i处的序列值，中位数是根据时间序列中所有的观测值对计算得出的，最终的Q由N来决定，具体如下：

本研究采用MK趋势检验法识别出通过显著性检验（P<0.05）的像元，同时采用Sen 斜率法量化物候的变化趋势，Sen>0 时表明物候期的趋势是推迟或延长；Sen<0 时，则表明物候期趋势是提前或缩短。

1.3.4 偏相关分析

偏相关分析是指当两个变量同时与第三个变量相关时，将第三个变量的影响剔除，只分析另外两个变量之间相关程度的过程。当控制变量个数为一时，偏相关系数称为一阶偏相关系数；控制变量个数为二时，偏相关系数称为二阶相关系数［43］。本研究采用一阶偏相关系数分析物候与气温、降水的关系，计算公式如下：

式中：Rxy·z是控制变量z后x和y的相关系数；rxy、ryz、rxz分别是变量x和y、y和z、x和z的相关系数。

通过MATLAB 软件将SOS 与3—5月平均气温和降水、EOS 和8—10 月平均气温和降水分别进行偏相关分析以及显著性检验，其中显著正相关为r>0且P<0.05；显著负相关为r<0且P<0.05。

2 结果分析

2.1 NDVI、EVI和SIF时间序列特征分析

2015—2019 年NDVI、EVI 和SIF 的滤波结果及其平均值如图2 所示。整体来看，可以发现三条曲线整体变化特征基本一致，生长季数值都远高于非生长季且均表现出明显的季节性变化规律。NDVI与EVI 波形相近，但EVI 进入生长阶段的时间滞后于NDVI，而结束生长阶段的时间早于NDVI。EVI整体分布更加集中，能更加准确反映植被生长季的变化特征。相比于NDVI，SIF 进入生长阶段的时间滞后于NDVI，而结束生长阶段的时间早于NDVI。SIF 与EVI 几乎同一时间开始快速增长，SIF 到达峰值后快速下降，比EVI 先结束生长季，所以SIF 整体分布更加集中。与NDVI和EVI相比，SIF 具有植物生理基础，并且与植被的生化过程直接相关。但是目前能获取到的SIF 数据时空分辨率较低，无法满足长时序和小区域植被物候的研究需求。SIF 和NDVI 的时间序列特征差异很大，而和EVI 的基本一致。已有研究表明NDVI在高植被覆盖区域会出现过饱和现象［44］，因此，本研究采用EVI开展大兴安岭多年冻土区的植被物候研究。

图2 2015—2019年大兴安岭多年冻土区NDVI、EVI和SIF的时间序列曲线（a）及其均值（b）［（b）图中的点代表植被指数的快速升高点和降低点］Fig. 2 Time series curves （a） and average value （b） of NDVI、EVI and SIF from 2015 to 2019 in permafrost regions of Greater Khingan Mountains ［The points in （b） figure indicate rapid increases and decreases in the vegetation index］

2.2 SOS、EOS和LOS空间分布格局

为了分析研究区植被物候的空间特征，逐像元计算2000—2019 年大兴安岭多年冻土区植被SOS、EOS 和LOS 等物候参数的均值（图3）。结果显示，SOS 主要集中在96～144 d，即研究区的植被在4 月上旬到5 月下旬进入生长季，占总面积的97.67%。植被SOS 的空间分布差异比较明显，高值区域主要集中在岛状融区多年冻土区，低值则主要分布在大片连续多年冻土区。研究区SOS 的平均值为129.46 d，大片连续多年冻土区平均值为127.29 d，岛状融区多年冻土区平均值为132.78 d，大片连续多年冻土区的SOS 平均值小于岛状融区多年冻土区（表1）。EOS 主要集中在272～320 d，研究区生长季结束的范围为10 月初到11 月中旬，占总面积的93.77%。研究区EOS 的平均值为295.11 d，大片连续多年冻土区平均值为297.96 d，岛状融区多年冻土区平均值为290.73 d。大片连续多年冻土区的EOS 平均值大于岛状融区多年冻土区（表1）。研究区LOS 绝大部分都集中在128～224 d 范围内，占到总面积的98.60%。研究区全区LOS 的平均值为165.65 d，大片连续多年冻土区平均值为170.67 d，岛状融区多年冻土区平均值为157.95 d（表1）。大片连续多年冻土区植被LOS 大于岛状融区多年冻土区。

表1 2000—2019年植被物候参数空间分布的平均值Table 1 Average value of vegetation phenology parameters from 2000 to 2019

图3 2000—2019年大兴安岭多年冻土区植被物候参数的空间分布图Fig. 3 The spatial distribution of phenological parameters in permafrost regions of the Greater Khingan Mountains from 2000 to 2019

2.3 SOS、EOS和LOS变化趋势

2.3.1 SOS、EOS和LOS年际变化趋势

大兴安岭多年冻土区2000—2019 年SOS、EOS、LOS 整体上年际变化（图4）的趋势不明显，且研究全区和大片连续多年冻土区、岛状融区多年冻土区的变化规律基本一致，即同步升高和降低。SOS 在2002 年和2014 年有两次明显低值，在2004年有一处明显的高值；EOS 在2002 年有一处低值，在2012 年和2015 年有两处明显高值；LOS 在2015年左右有一处明显的高值，其他年份则无明显变化幅度。

图4 2000—2019年大兴安岭多年冻土区植被物候参数的年际变化Fig. 4 Inter-annual variation of vegetation phenology parameters in the permafrost zone of Greater Khingan Mountains from 2000 to 2019

2.3.2 SOS、EOS和LOS变化趋势的空间格局

结合MK 检验和Sen斜率法计算2000—2019年SOS、EOS 和LOS 在像元尺度的空间变化趋势，将不具有显著性（P>0.05）和变化趋势为0 的像元排除，如图5所示，只有少部分的像元通过了显著性检验。研究区SOS 显著变化的像元数为4.5×104个，仅占全部像元的2.24%。其中位于大片连续多年冻土区的像元有2.54×104个，位于岛状融区多年冻土区的像元有1.96×104个。20 年来SOS 变化趋势的范围从提前6.67 d 和推迟5.33 d 不等，变化趋势的平均值为-1.23 d·（20a）-1（表2），整体SOS 呈提前趋势。分区来看，大片连续多年冻土区的变化趋势范围为-6.67～3.56 d，平均值为-1.19 d·（20a）-1，岛状融区多年冻土区的变化趋势范围为-6.58～5.33 d，平均值为-1.28 d·（20a）-1，二者的SOS 均呈提前趋势。

表2 2000—2019年研究区物候参数的变化趋势Table 2 The change trend of phenology parameters in the study area from 2000 to 2019

研究区EOS 显著变化的像元数为12.1×104个，占全部像元的6.04%，其中位于大片连续多年冻土区的像元有8.12×104个，占显著变化像元的67.12%，岛状融区多年冻土区的像元有3.98×104个，占显著变化像元的32.88%。20 年来EOS 趋势变化的范围为-10.00～9.33 d，变化趋势的平均值为-0.46 d·（20a）-1（表2），整体EOS 变化趋势呈提前趋势。分区来看，大片连续多年冻土区像元的变化趋势范围为-9.14～9.33 d，平均值为-0.32 d·（20a）-1，岛状融区多年冻土区像元的变化趋势范围为-10.0～8.73 d，平均值为-0.76 d·（20a）-1。研究区LOS变化趋势的平均值为2.39 d·（20a）-1，大片连续多年冻土区为1.91 d·（20a）-1，岛状融区多年冻土区为2.78 d·（20a）-1，过去20 年大兴安岭多年冻土区植被的LOS呈延长的趋势。

2.4 气温和降水对植被物候的影响

2.4.1 气温和降水的变化趋势

图6 为研究区2000—2019 年气温和降水量的变化情况。全年和3—5 月气温呈波动上升，8—10月气温波动变化呈微弱下降趋势，总体上无显著趋势。全年和3—5 月、8—9 月降水量呈波动变化，呈微弱上升趋势。

图6 2000—2019年气温和降水量的变化Fig. 6 Changes in temperature （a） and precipitation （b） from 2000 to 2019

2.4.2 气温和降水对SOS的影响

由于研究区SOS多年平均值为129.46 d（表1），植被开始生长的时间与3—5 月气温和降水密切相关。本文将SOS 与3—5 月的平均气温和总降水量进行逐像元偏相关分析（图7）。总体上看，植被SOS 受气温的影响较降水明显（面积比例分别为15.51%和3.15%）。大兴安岭多年冻土区植被SOS与3—5 月气温的显著偏正负相关的面积占比分别为0.48%和15.03%，表明大部分区域随着气温升高，SOS 呈提前的趋势。SOS 与3—5 月降水的显著偏正负相关的面积占比相对较少，分别为2.20%和0.95%。降水对SOS 的影响有很强的空间异质性，呈正相关部分多分布于研究区北部，大部分区域的降水增多会导致SOS的推迟。

图7 植被SOS与3—5月气温和降水的偏相关分析结果Fig. 7 Skewed correlation analysis of vegetation SOS and March—May temperature （a） and precipitation （b）

2.4.3 气温和降水对EOS的影响

研究区多年平均EOS 值为295.11 d（表1），植被结束生长的时间与8—10 月气温和降水密切相关。本文将EOS与8—10月的平均气温与总降水量进行逐像元偏相关分析（图8）。可以看出，植被EOS 受气温和降水的影响较大，且受气温影响较降水显著（显著相关的面积比例分别为7.35%和3.96%）。EOS 与8—10 月的气温呈明显的正相关性，其中显著偏正负相关的面积占比分别是6.81%和0.54%，即降水一定的情况下，气温升高，EOS 将推迟。EOS 与8—10月降水量的显著偏正负相关的面积占比分别为3.46%和0.50%，表明研究区的EOS 与降水量主要表现为正相关，即温度一定时，降水增加，EOS推迟。

图8 植被EOS与8—10月气温和降水的偏相关分析结果Fig. 8 Skewed correlation analysis of vegetation EOS and August—October temperature （a） and precipitation （b）

3 讨论

本研究中EVI 比NDVI 更适合大兴安岭多年冻土区的物候研究，这是由于NDVI 算法仅使用红光和近红外波段，当植被覆盖率很高时，红光波段会迅速饱和，从而产生饱和效应［45］。EVI 在NDVI 的基础上改进了算法，使用了蓝光波段并改进了残留气溶胶的处理方法，减少了大气和土壤背景的影响，避免了植被覆盖率高的地区出现饱和的问题［46-49］。所以EVI在大兴安岭地区的性能优于NDVI。目前，已有一些学者采用不同的数据和方法开展东北地区植被物候研究。Tang 等［39］估算了1982—2012 年大兴安岭地区的SOS 和EOS，发现SOS 主要分布在第90～150 d 之间，而EOS 的范围则是第245～305 d之间。Yu 等［50］计算了1982—2015 年中国东北地区的SOS，范围从一年中的第100～140 d 不等，而EOS的范围为280～320 d。Liu等［51］计算了中国温带植被的平均EOS 值，其结果主要分布在第270～310 d。本文SOS 分布在第96～144 d，EOS 的范围为第272～320 d，与已有结果基本一致。像元空间分辨率越高，越能够减少混合像素对物候信息提取的干扰。本文使用的EVI 数据像元大小为250 m，高于大多数已有研究的空间分辨率，其物候提取的精度也会有所提高。

不同区域植被物候特征的差异主要取决于植被类型、高程和气候等因素［52］。本研究发现大片连续多年冻土区的SOS 均值小于岛状融区多年冻土区，EOS 的均值大于岛状融区多年冻土区，这与Fu等［38］和Yu 等［50］的研究结果一致。岛状融区多年冻土区西南部为森林与草原的过渡地带，主要植被类型包括草原、阔叶林和农田等，不同植被类型生长习性不同，物候差异较大。草原和农田的LOS 小于森林的LOS，阔叶林的LOS 小于针叶林［53］。不难发现SOS、EOS 和LOS 等物候参数的空间格局与植被类型密不可分。高程差异带来最明显的影响是水热条件不同，从而导致不同海拔的物候差异。通常随着海拔的升高温度降低，而本研究区位于西伯利亚冷高压的边缘，冷高压中心空气下沉使得高压边缘上空空气不足，因此空气绝热增温上升，从而出现低海拔处的温度低于高海拔地区“逆温”现象，这种现象在大兴安岭地区一年四季均存在［54-55］。将图3结果与DEM（图1）叠加发现LOS 高值和高海拔区域基本对应，大片连续多年冻土区的海拔高于岛状融区多年冻土区，大片连续多年冻土区的LOS 大于岛状融区多年冻土区。大兴安岭多年冻土区逆温现象的存在以及植被类型的不同是导致不同类型冻土区植被物候差异的主要原因。

近20 年植被物候参数在研究区尺度上没有显著的变化趋势。在像元尺度上，大部分的像元不具备显著性的变化趋势，只有少部分像元通过了显著性检验。对显著变化的像元进行分析发现，研究区SOS 和EOS 整体上呈提前趋势，EOS 呈微弱的提前趋势，LOS 整体上呈增长趋势。大片连续多年冻土区和岛状融区多年冻土区SOS、EOS、LOS变化趋势同研究区一致。通过对植被物候与气温和降水进行偏相关分析发现大兴安岭多年冻土区植被物候主要受到气温的显著影响。研究区全年温度和降水呈增加趋势，不同月份表现不同的变化趋势，3—5月气温上升，降水下降，都有利于SOS 的提前。而8—10月气温呈下降趋势，降水上升，导致EOS 呈现微弱的提前趋势。与EOS 相比，SOS 对气温的响应更为显著，SOS 提前与气温升高密切相关，这与俎佳星等［56］和丛楠等［57］结果基本一致。

4 结论

本文基于SIF、NDVI 和EVI 等植被指数开展大兴安岭多年冻土区植被物候的研究。首先比较了三种指数提取物候的差异和适应性，发现EVI 最适用于大兴安岭地区的物候研究。采用MODIS EVI提取SOS、EOS 和LOS 等关键植被物候参数，分析大兴安岭多年冻土区植被物候的时空变化及其对气候变化的响应。主要结论如下。

（1）大兴安岭多年冻土区NDVI、EVI 和SIF 的时间序列能够反映植被的季节变化，可用于植被物候信息的提取。EVI 在高植被覆盖区性能优于NDVI，同时EVI 与SIF 曲线更加一致。EVI 最适合大兴安岭多年冻土区植被物候的研究。

（2）由于逆温现象的存在和植被类型的差异，大兴安岭大片连续多年冻土区的SOS 均值小于岛状融区多年冻土区，大片连续多年冻土区的EOS 均值大于岛状融区多年冻土区，大片连续多年冻土区的LOS大于岛状融区多年冻土区。

（3）过去20 年大兴安岭多年冻土区植被的SOS呈提前趋势，EOS 呈微弱提前趋势，SOS 提前的趋势大于EOS，因此研究区LOS 呈延长的趋势。大片连续多年冻土区和岛状融区多年冻土区呈现相同的变化趋势。

（4）大兴安岭多年冻土区气温对植被物候的影响程度高于降雨量，部分区域植被SOS 和EOS 分别受到3—5 月、8—10 月气温的影响。3—5 月气温升高，对SOS 有提前作用；8—10 月气温降低，对EOS有提前作用。