基于MODIS数据的新疆草地物候提取方法及变化趋势分析

张仁平，郭靖，马晓芳，郭伟勇

（1.新疆大学资源与环境科学学院，新疆 乌鲁木齐 830046；2.新疆大学绿洲生态教育部重点实验室，新疆 乌鲁木齐 830046；3.新疆林业科学院，新疆 乌鲁木齐 830000；4.中国科学院西北生态环境资源研究院，甘肃 兰州 730000；5.察布查尔县农业农村局，新疆 伊犁 835300）

植被物候是指受气候和环境因素的影响而出现的以年为单位的周期现象，它是指示气候与自然环境变化的重要指标，对于研究全球生态环境变化至关重要［1—3］。以往的植被物候观测是通过定点定时的野外目视观察法，虽然能够准确地获取观测点周围的植被物候信息，但无法应用推广到大区域［4］。遥感监测因具有长时间序列、覆盖范围广等特点，目前成为物候观测的主要手段［5—6］。归一化植被指数（normalized difference vegetation index，NDVI）能够反映植被光合有效辐射的变化情况，被广泛应用于多尺度下植被动态变化监测、土地覆被变化、植物生物物理参数反演及区域环境变化等方面［7—9］，但由于传感器本身及其受云、大气气溶胶、冰雪等因素干扰，从卫星数据中获得的NDVI时间序列数据不可避免地受到不同程度的噪声污染，无法真实的反映植被生长状况，因而需要对其进行重建［10—12］。NDVI时序数据重构方法有很多，其中双逻辑斯蒂函数拟合法（double logistic，D-L）仅依赖NDVI曲线的变化特征，即可客观地反映植被生长过程的细节特征，且该方法对噪声的处理及物候特征提取能力较为突出［13］。S-G滤波法（Savitzky-Golay）通过NDVI数据的上包络线来拟合NDVI时序数列的变化趋势，保留了NDVI曲线的特征，通过一个迭代的过程可消除偏离正常生长趋势线的噪声，而使平滑达到最好的效果。同时该方法不仅不受数据时间空间尺度和传感器的限制，还能清晰地描述NDVI时间序列的长期变化趋势以及局部的突变信息［14］。非对称性高斯函数拟合法（asymmetric gaussian，A-G）能够较好地描述NDVI时序数据的总体变化趋势和全局特征，在噪声点优化及保留原始数据质量方面，具有更好的拟合效果［15］。基于此，本研究选用A-G、D-L和S-G算法重构NDVI时序数据。美国NASA根据Zhang等［16］基于Logistic函数对美国东北地区的MODIS EVI（enhanced vegetation index）进行拟合，形成物候产品数据（MCD12Q2）。因为此种算法不需要界定植物物候的阈值，普遍适用于大范围的物候监测［17—18］。因此，选择最佳的遥感物候提取方式，对于特定区域植被物候监测及准确进行物候监测具有一定的意义。

目前，已有大量关于草地物候时空变化的研究。在国外，Delbart等［19］利用Spot-vegetation（VGT）和NOAA advance very high resolution radiometer（AVHRR）数据分析了1982—2004年欧亚大陆植物物候变化趋势，结果表明该地区植物返青日期有所提前。Zeng等［20］利用MODIS NDVI数据对2000—2010年北半球的植物物候研究显示，植被返青仍然保持提前的趋势。在国内，王宏等［21］基于NOAA NDVI和MSAVI数据对中国北方植物物候变化趋势的分析表明，荒漠草原和典型草原的返青期（start of growing season，SOS）提前，枯黄期推迟。倪璐等［22］基于1986—2015年的GIMMS NDVI数据，定量分析了中国天然草地物候的时空变化特征，结果表明1998年以前全国草地SOS平均提前速率为0.37 d·a—1，1998年以后草地物候变化趋势相反，草地SOS平均推迟速率为0.28 d·a—1。何宝忠等［23］通过MODIS NDVI数据提取了新疆地区2001—2016年典型植被物候期，发现新疆植被物候具有明显的纬向分布和垂直地带性分布特征，且与全球大背景下典型植被物候特征变化趋势相反，新疆植被SOS呈推迟趋势，推迟幅度为1.9 d·10 a—1。然而，张仁平等［24］基于MCD12Q2数据，对新疆地区2001—2014年草地物候的时空变化进行了研究，前者的研究表明新疆草地返青期整体呈现提前趋势，提前速率为0.11 d·a—1。通过以上研究可以看出，不同区域和不同时间段物候的变化区域存在着一定的差异，进一步研究新疆地区草地物候变化趋势，有助于更深刻地理解气候变化带来的影响。

1 材料与方法

1.1 研究区概况

新疆维吾尔自治区地处我国西北隅（34°22′—49°33′N，73°22′—96°21′E），总面积为1.66×106km2，约占全国总面积的1/6。新疆具有“三山夹两盆”地理环境特征，自南向北依次为昆仑山脉地—塔里木盆地—天山山脉—准格尔盆地—阿尔泰山脉，平均海拔2000 m以上。由于远离海洋，且四周高山阻隔，该区形成了典型的温带大陆性气候，具有气温温差较大、降水量少，光热资源充足等特点。新疆草地资源丰富，主要包括低地草甸、高寒草原、高寒草甸、温性荒漠、温性草原、高寒荒漠等6种草地类型，草地总面积约57.26万km2，居全国第3位。该区草地资源不仅是当地农牧民重要的生产资料，而且草地生态系统关系着当地和周边地区的生态与环境状况以及国家生态安全（图1）。

图1 研究区草地类型及地面观测点分布Fig.1 Distributions of the grasslands and ground observation points in Xinjiang

1.2 数据源

NDVI数据来源于美国地质勘探局（https：//earthdata.nasa.gov/eosdis/daacs/lpdaac）的MOD13Q1 V 006，该数据是由16 d合成的250 m×250 m空间分辨率的植被指数产品，本研究下载了2001—2019年共计437景影像为数据源，对新疆地区的草地返青期进行提取和分析。物候产品源于美国国家航空航天局（https：//lpdaac.usgs.gov/）的MODIS土地覆盖动态产品（MCD12Q2），时间为2001—2019年，空间分辨率为500 m×500 m。为使遥感反演物候数据与后期地面物候观测资料的对比分析结果较为可靠，本研究将MOD13Q1 NDVI重采样至500 m。

地面草地返青数据采用2012—2014年2—5月的物候观测数据。该数据共计123个样地，样地大小为500 m×500 m，每个样地设5～8个样方，每个样方大小为1 m×1 m，每隔5～7 d对采样点完成一次观测。当样地中植物开始展叶并正常发育的株（丛）或返青盖度占总株（丛）的百分率为50%时，即认为该样地返青［6］。原始地面物候观测资料为日期型格式，需将其转换为数值型格式以便后期研究分析。

1.3 研究方法

1.3.1NDVI时间序列重构 1）Savitzky-Golay滤波法（S-G）是一种基于最小二乘法的卷积算法，由Savitzky等［25］于1964年最早提出。其基本原理：用一定长度的过滤器和待处理数据作卷积，对待处理的数据作加权多项式拟合，拟合的目标是求得最小均方根误差，而一些远离大多数点的边沿点不会参与拟合。这样拟合时过于偏离正常生长趋势线的噪声部分会被丢弃。其公式可表示为：

式中：Y*j为平滑后的NDVI数据；Y j+i代表原始NDVI数据；Ci为拟合系数，表示从滤波器首部开始第i个NDVI值对应的权重；N为滤波器的长度，其值越大滤波结果越平滑，但也会遗漏能反映真实NDVI的细节数据，反之易产生大量的重复数据。因此本研究经过反复的实验，将Savitzky-Golay滤波窗口的值设为4。

颅内静脉窦血栓属于特殊性的脑静脉系统缺血性脑血管病，患者的主要临床症状包括意识障碍、视乳头水肿以及局灶性神经功能障碍等[1]。颅内静脉窦血栓具有发病形势多样、发病原因复杂以及临床症状、体征缺乏特异性等特点，且患者病情进展迅速，增加了误诊以及漏诊率[2-3]。因此，本研究通过分析诊断时间与颅内静脉窦血栓患者临床特征及预后的关系，旨在为颅内静脉窦血栓的早期诊断提供参考依据。现报道如下。

2）非对称性高斯函数拟合法（A-G）是基于不对称高斯函数的非线性最小二乘拟合算法，由Jonsson等［26］于2002年提出。其原理是对植被生长过程用分段高斯函数模拟，然后通过对各高斯拟合曲线进行平滑连接实现时间序列重建，公式可表示如下：

式中：f（t）为局部拟合函数；t为时间；c1和c2控制曲线的基准和幅度；a1决定峰值和谷值的位置；a2，a3，a4，a5控制曲线的宽度和陡峭度；g（a1，…，a5）为时间t的高斯函数；局部拟合后通过整体函数F（t）连接。

式中：F（t）为整体函数；tl，tr为时间序列数据中未拟合部分的区间；（t）为［0，1］之间的剪切数；fl（t）、fc（t）、fr（t）分别为区间内左边最小值、中间最高值和右边最小值的局部拟合函数。

3）双逻辑斯蒂函数拟合法（D-L）是一种半局部拟合的经典方法，由Beck等［27］在2006年提出。其原理是首先将整个时间序列中时间点对应值按极大、极小分成多个区间，然后分别对该区间进行D-L函数局部拟合，其局部拟合方式与非对称高斯拟合方法（A-G）类似。该方法具有无需确定阈值，仅依赖NDVI曲线的变化特征，即可客观地反映植被生长过程的特征。其公式表示为：

式中：a1，a2，a3，a4控制左、右半部分的宽度和陡峭度。其整体连接函数F（t）与非对称性高斯函数拟合法（A-G）相同，见公式4。

1.3.2物候期的提取 植被物候信息提取的方法较多，常见的有拟合法、最大斜率法、主成分分析法、滑动平均法、阈值法等［28］，虽然以上方法都是基于植物的生理特征和遥感机理，但对于不同的应用区域、数据质量等方面有所差异。本研究基于前人研究经验及新疆地区实际情况［23，29］，采用Jonsson等［26］提出的动态阈值法，将NDVI增长达到当年NDVI振幅20%的时刻定义为返青期的阈值，基于TIMESAT软件逐年提取新疆地区的草地返青期。

MCD12Q2产品是以MODIS增强型植被指数（enhanced vegetation index，EVI）为数据源，先对EVI时序数据进行预处理，消除气溶胶、云和观测角度等影响，然后用滑动窗口（连续5个时相EVI值）来判断持续升高和持续降低区间，当所在区间极大值和振幅满足给定阈值条件后，判断该升高—降低区间为一个生长周期［30］，最后对于每一个生长周期采用分段Logistic函数拟合，并利用其曲线变化率确定植被生长始期、成熟期、末期等。其中，将EVI增长的节点日期定义为返青期。

1.3.3对比验证 地面实测数据是植被物候遥感监测精度评价的最有效手段。基于新疆地区2012—2014年共计123个样地观测数据作为真值，采用皮尔森相关系数（R），偏差（Bias）和均方根误差（root mean square error，RMSE）对4种方法提取的返青期进行对比分析。

皮尔森相关系数［31］是用来描述两个随机变量线性相关的统计量，其值介于—1与1之间。

式中：n为用于验证的样本数据；x i为遥感提取的返青期；xˉ为由遥感提取的2012—2014年返青期平均值；yi为地面观测的返青期；为由地面观测的2012—2014年返青期平均值。当R＞0时，表明遥感提取的返青期与地面观测的返青期呈现正相关；当R＜0时，结果相反；当R＝0时，表示遥感提取的返青期与地面观测的返青期不存在线性相关关系。

Bias是指遥感返青期与地面实测值之间的差别，其描述的是模型本身的准确度［32］。其计算公式如下：

RMSE同样用来表示遥感返青期与地面实测值的偏差，从而反映不同方法提取返青期的模拟精度［33］。

式中：N代表地面观测点的个数；Ei为第i个遥感提取的返青期；Oi为第i个地面观测点的返青期。Bias与RMSE越接近0，表明遥感提取的返青期与地面观测的返青期越接近。

1.3.4趋势分析 本研究采用线性回归模型，从像元尺度上来定量描述新疆地区草地返青期（start of growing season，SOS）的变化趋势，以此来分析研究区草地物候的时空变化特征［34］。其公式表示如下：

式中：b代表斜率，其值大于0，说明变化趋势呈现推迟的趋势，相反则是提前的趋势；n为年数；i为年序号，Pi表示第i年对应的植被返青期。

同时，采用F检验对变化趋势的显著性进行评价。

式中：k＝1，n＝19。通过F检验表查不同水平对应的R值，得到F0.01（1，18）＝8.53，F0.05（1，18）＝4.94，F0.1（1，18）＝3.05。当|R|≥0.59时，通过置信度0.01水平的检验（变化趋势极显著）；当0.49＜|R|＜0.59时，通过置信度0.05水平的检验（变化趋势显著）；当0.40＜|R|＜0.49时，通过置信度0.1水平的检验（年际变化趋势略显著）。

2 结果与分析

2.1 对比验证

研究中用于4种方法对比分析一致的地面实测点，均为123个。受MCD12Q2产品数据缺失较多的影响，研究区部分地面实测点对应的MCD12Q2产品返青期为空值，最终可用于分析的仅有26个地面实测点。

从图2可以看出，2012—2014年新疆地面实测点返青期主要集中在每年第68～122天，A-G、D-L、S-G反演的返青期主要集中在第60～120天，而MCD12Q2产品的返青期在第100～130天。可见，A-G、D-L和S-G反演的返青期比地面实测返青期提前，而MCD12Q2产品返青期推迟。从返青期提取结果与地面观测数据的相关系数R来看，A-G、D-L返青结果和地面实测数据的相关性较高，分别达到0.879、0.880，而S-G、MCD12Q2产品与地面实测数据的相关性较低，分别为0.870、0.626。从返青期提取结果与地面观测数据的误差指标RMSE来看，MCD12Q2最小，RMSE为12.692 d，其次是A-G和D-L监测结果，RMSE分别为16.395 d、18.754 d，而S-G方法监测结果与实测返青期误差量较大，RMSE大于20 d。从返青期提取结果与地面观测数据的拟合度Bias来看，AG和D-L提取结果与地面观测数据的拟合度较优，分别为—0.058、—0.059，MCD12Q2提取结果与地面观测数据的拟合度较差，为—0.163。上述分析结果表明，4种方法中，基于A-G方法返青期提取结果精度较高。

图2 返青期观测数据与遥感反演数据的对比分析Fig.2 The comparison between phenology observed in the ground stations and phenology extracted by asymmetric gaussian，double logistic，Savitzky-Golay and MCD12Q2 retrieval

由于遥感提取的物候参数与地面实测存在尺度不统一等问题［35］，因而进一步从空间层面上探究4种方法提取新疆草地返青的最优模型。标准差衡量的是样本数据的离散程度。一般来说，标准差越小，表明数据越聚集；标准差越大，则表明数据越分散。

从图3可以看出，4种方法反演的返青期在北疆地区的标准差总体上较低，而南疆地区返青期的标准差较高。结合表1，4种方法中返青期标准差（＜15 d）和标准差（15～30 d）的面积比例较大。A-G、D-L、S-G及MCD12Q2产品的返青期标准差（＜30 d）的面积比例分别为82.19%、76.53%、67.98%、96.81%，以A-G方法提取的返青期标准差（＜30 d）的面积比例最多。MCD12Q2产品的返青期标准差（＜30 d）的面积比例最高，但是其反演的返青期在空间上空值较多。

表1 不同方法反演的返青期标准差面积比例Table1 Standard deviation area ratio of SOSinverted by different methods（%）

图3 不同方法反演的返青期（SOS）标准差分布Fig.3 Distr ibution of the standar d deviation of SOSin Xinjiang inverted by differ ent methods

综上所述，4种方法中A-G方法反演的返青期的结果最佳。因此，后面关于新疆草地返青期的空间分布及时间变化趋势等的研究分析都是基于A-G提取的返青期来展开的。

2.2 草地物候空间分布特征

图4为2001—2019年新疆草地平均返青期的空间分布格局。新疆地区草地返青期分布具有明显的地域差异，呈现出自北向南逐渐推迟的分布特征。2001—2019年该地区返青期主要集中在第60～140天，多年均值约为第100天。其中，北部准噶尔盆地和伊犁河谷区域的草地返青时间最早，早于第80天，而阿尔泰山、天山中部及昆仑山等区域的草地返青时间最晚，晚于第140天。结合图5得到，新疆不同草地类型返青期存在一定的差异。其中，高寒草甸与高寒草原的返青期最晚，约为第130天，其次是高寒荒漠与低地草甸，分别为第127、124天，再次为温性草原，为第103天，温性荒漠返青时间最早，约为第86天。

图4 基于A-G方法新疆草地返青期（SOS）的空间分布Fig.4 Spatial distribution of SOSin Xinjiang based on A-G method

图5 基于A-G方法不同草地类型返青期（SOS）的多年均值Fig.5 Multi-year mean values of SOS for different grassland types based on the A-G method

2.3 草地物候年际变化特征

基于像元尺度对2001—2019年新疆草地的返青期与年份之间进行了线性拟合，得到过去19年间该区返青期的变化趋势。从图6可以看出，北部准噶尔盆地边缘地带、天山南部及昆仑山东段少数区域的草地返青期呈现提前的趋势，提前的面积比例约为46.93%，提前幅度多为1 d，而草地返青期推迟的区域主要集中于天山东部和昆仑山北部边缘地区，推迟的面积比例约为53.07%，推迟幅度约为1 d。总体上来看，草地返青期通过P＜0.1的有效像元比例约为10.00%，其中2.0%像元的SOS在0.05的显著性水平上呈显著的提前趋势，平均每年提前1～2 d，主要分布在阿尔泰山以南、准噶尔盆地西部边缘区，而2.7%像元的SOS在0.05的显著性水平上呈显著的推迟趋势，平均每年约推迟2 d，主要在塔里木盆地西北部边缘及准噶尔盆地东部少数区域。

图6 基于A-G方法新疆草地返青期（SOS）的空间变化及显著性检验Fig.6 Var iation tr end and significance test of grassland SOSin Xinjiang based on A-G method

由图7可知，低地草甸、温性荒漠和高寒荒漠类型下，A-G方法提取的返青期均呈现推迟的变化趋势，温性草原、高寒草原及高寒草甸的返青期表现为提前。从返青期推迟变化来看，温性荒漠的SOS推迟最为明显，推迟速率为0.456 d·a—1，其次是低地草甸，推迟速率为0.379 d·a—1，而高寒荒漠的返青期稍有推迟，推迟速率为0.313 d·a—1。从返青期提前变化来看，温性草原的返青期明显提前，提前速率为0.384 d·a—1，高寒草原次之，返青期的提前速率为0.196 d·a—1，而高寒草甸的返青期略有提前，提前速率为0.127 d·a—1。分时间段来看，2001—2010年，低地草甸、温性荒漠、温性草原、高寒草原及高寒草甸的返青期均呈现推迟趋势，分别推迟2、13、13、1、11 d，而2010—2019年温性荒漠、温性草原、高寒草原及高寒草甸的返青期变化趋势相反，分别提前8、14、3、5 d。

图7 基于A-G方法不同草地类型返青期（SOS）的年际变化趋势Fig.7 The interannual variation trend of SOSin different grassland types based on A-G method

3 讨论

不同的遥感物候提取方法各有优缺点，本研究采用4种不同构建原理的植被物候提取算法，提取了新疆地区草地植被返青期。通过123个地面实测数据对不同方法的提取结果进行比较和分析，发现A-G方法反演的返青期结果和地面实测数据的相关性较高，拟合度较优，误差相对较小。虽然4种数据中MCD12Q2产品的返青期标准差（＜30 d）的面积比较最大，但由于该产品数据缺失较多，地面有效实测点仅有26个，因而该方法验证结果误差相对较大。

新疆属于典型的干旱半干旱地区，生态系统脆弱，研究该区域的植被物候特征变化趋势，对深入认识气候变化对陆地生态系统的影响具有典型性。安佑志［36］的研究发现，受海拔因素的影响，新疆北部山区的返青时间要比同纬度地区植被稍晚些，其中阿尔泰山的返青时间最晚，准噶尔盆地和伊犁河谷地区的返青时间最早。同时得到20世纪90年代之前新疆北部地区的返青期以提前为主，南部部分地区的返青期显著延迟，本研究估算的返青期具有明显的垂向分布差异结论与此基本一致。已有研究表明，气候变暖导致了21世纪前物候返青期的提前，但随着升温程度的减弱，这种趋势在近几十年可能会发生变化［3，22］。Piao等［37］认为，自1982年至今青藏高原植被返青期并没有明显的变化趋势，而是在不同阶段返青期提前和推迟有所不同。张雯等［38］的研究表明，1982—2013年内蒙古地区植被返青期总体呈现推迟趋势，尤以2010—2013年推迟趋势较为明显。从新疆草地返青期整体变化趋势来看，本研究与何宝忠等［23］的研究结果较为一致，均认为新疆草地的返青期呈现推迟趋势。从以上分析可以看出，植被物候的变化与植被类型、研究区域、研究时段以及遥感数据等有关系。

4 结论

本研究通过4种不同物候提取方法，进行新疆草地返青期的提取，并基于地面实测数据，对遥感提取返青期进行精度对比验证，同时基于最优模型分析了新疆草地返青期时空变化。主要结论如下：

1）4种方法中A-G方法反演的新疆草地返青期的结果最佳。基于A-G方法提取的新疆草地返青期和实测点的草地返青期相关性较高（0.879），RMSE相对较小（16.395 d）。同时A-G方法提取的新疆草地返青期标准差（＜30 d）的面积比例最多，达到82.19%。

2）新疆地区草地返青期主要集中在第60～140天，具有自北向南逐渐推迟的明显地域差异。其中，北部准噶尔盆地和伊犁河谷区域的草地返青时间最早，而阿尔泰山、天山中部及昆仑山等区域的草地返青时间最晚。不同草地类型返青期存在一定的差异，以高寒草甸与高寒草原的返青时间（第130天）最晚，而温性荒漠返青时间（第86天）最早。

3）新疆地区的草地返青期总体呈现微弱的推迟趋势，其中阿尔泰山以南、准噶尔盆地西部边缘区（2.0%像元）的SOS在0.05的显著性水平上呈显著的提前趋势，平均每年约提前1～2 d，而塔里木盆地西北部边缘及准噶尔盆地东部少数区域（2.7%像元）的SOS在0.05的显著性水平上呈显著的推迟趋势，平均每年约推迟2 d。不同草地类型中，低地草甸、温性荒漠和高寒荒漠的返青期均呈现推迟的变化趋势，而温性草原、高寒草原及高寒草甸的返青期表现为提前趋势。此外，温性荒漠、温性草原、高寒草原及高寒草甸的返青期以2010年为界，呈现“先推迟后提前”的变化趋势。