植被物候是研究植物生命周期中反复发生的事件，包括植物的发芽、展叶、开花、叶绿及枯黄等，是生态系统应对气候变化的敏感性指标^[1]，除此之外，植被物候的年际变化以及年内变化可以反映出陆地植被与大气之间的碳循环、水循环以及能量交换^[2-3]。监测物候有利于量化植被的季节性周期和跟踪不规则的自然事件^[4]。因此准确监测全球及区域尺度的植被物候动态对理解气候—生物圈的相互作用意义重大。

植被物候的监测方法主要包括通过人眼记录物候事件、定期摄影、自动重复摄影和卫星遥感^[4-5]。最初的人工监测不仅监测区域十分有限，并且耗费大量的人力物力^[5]。自1980年代中期以来，遥感技术和近地面物候观测技术的发展，极大提高了区域到全球植被物候的认识^[6]。同时遥感数据中的植被指数产品拥有监测植被变化趋势的高时空分辨率特征，可以很好地描述生态系统的物候变化^[7]。许多研究显示，近几十年来，全球气候在逐渐变暖，在此背景下，植被生长也发生变化，其中较为关键的3个植被物候参数：植被生长季开始日期(SOS)、生长季结束日期(EOS)和生长季长度(GSL)在年际尺度上表现出了明显的变化^[7]。同时在不同的海拔和经纬度，植被物候也表现出不一样的空间分布特征^[8]。在北半球大多数陆地生态系统中，植被物候的差异变化归因于海拔与气候变化，尤其是温度和降水变化^[8]。而Wang等^[9]的研究表明在黄河流域，白天温度与夜间温度的增加对植被的生长有着相反的作用；黄文洁等^[10]、郭少壮等^[11]发现高原植被物候与海拔有着密切关系。李丹等^[12]对中国典型半干旱荒漠区植被进行研究发现，持续湿润的气候对沙地植被生长具有正向促进作用，而干旱则会抑制沙地植被生长，植被生长与降水和蒸散量相关性较强，而和温度的关系较小。

干旱、半干旱土地荒漠化是全球性环境问题。这些地区属于生态脆弱地带，对于气候变化十分敏感^[13]，深入研究该地区植被物候动态以及气候变化对植被生长的影响，对于指导荒漠化土地治理以及退化生态环境的恢复具有积极作用^[14]。毛乌素沙地作为中国典型干旱、半干旱过渡区和农牧交错带^[13]，对该地区植被物候的研究目前仍较少。本研究使用Double Logistic函数拟合法和动态阈值法对2000—2019年的NDVI数据进行拟合，提取植被物候关键参数，探究毛乌素沙地的植被物候的变化规律，结合温度和降水量数据，应用空间自相关分析对其与温度和降水的关系进行分析，旨在定量揭示毛乌素沙地的植被生态系统与气候变化之间的响应关系。

1 研究区概况与数据来源 1.1 研究区概况

毛乌素沙地位于36°48′—40°12′N，106°10′—111°53′E，行政区划上跨越3省12县(市区)，包括内蒙古自治区鄂托克前旗、鄂托克旗、乌审旗、伊金霍洛旗，陕西省榆林市、定边县、靖边县、横山区、神木市及宁夏回族自治区的灵武市、陶乐县、盐池县。毛乌素沙地研究区的沙地面积约为4.22×10⁴ km²，海拔高度255~1 948 m，整体海拔自西北向东南倾斜。毛乌素沙地主要植被类型为草地，研究区属于荒漠草原—干草原—森林草原的过渡带，属于典型的半干旱大陆气候^[14]。

1.2 数据来源

研究使用数据为16 d，250 m分辨率的MOD13Q1数据集。该时间序列数据集来源于美国的NASA LP DAAC工作组(http://lpdaac.usgs.gov/lpdaac/products/modis_products_table) (LAAD PROCESSES DISTRIBUTED ACTIVE ARCHIVE CENTER)。数据集包括2000年1月至2019年12月共460期影像。利用MRT工具从MOD13Q1产品的460期影像中提取NDVI数据层，并进行拼接、格式转换和重投影。原始MODI13Q1产品为-3 000~10 000的DN值，为了得到真实NDVI值，将DN值转为NDVI值的关系式为：

$ {\rm{NDVI}} = {\rm{DN}} \times 0.000{\rm{ }}1 $

(1)

然后将多年平均NDVI＜0.05的像元剔除，这些像元可能是荒漠、水体或者其他的非植被覆盖区。

气象数据来源于国家气象科学数据中心(http://data.cma.cn)毛乌素沙地及周围14个气象台站，其时间分辨率为1 d。为了与物候数据进一步进行相关性分析，对温度和降水量进行克里金插值，使气象数据分辨率达到250 m，从而与NDVI数据分辨率保持一致。

2 数据处理方法 2.1 植被物候的提取

植被物候参数的提取包括：①对时间序列进行拟合，降低噪声的影响；②提取物候关键参数。

TIMESAT可对多种类型传感器的时间序列数据进行处理^[15]。包括对时间序列数据进行重新拟合降噪，并得到生长季节的起始点等信息^[15]。对于植被指数时间序列，TIMESAT包括3种自带的拟合方法：非对称高斯函数拟合，Double Logistic(D-L)曲线拟合和Savitzky-Golay滤波法^[16]。本研究采用Double Logistic函数拟合法对植被指数进行拟合。Double Logistic函数是一种半局部拟合法。方法是首先将时间序列中的时间点值按照极大值或者极小值分为多个区间，然后依次对区间进行局部拟合^[17]。

$ \begin{aligned} g(t, {a_1}, \ldots , {a_4}) =& \frac{1}{{1 + {\rm{exp〔}}\frac{{\left( {{a_1} - t} \right)}}{{{a_2}}}〕}} - \\ &\frac{1}{{1 + {\rm{exp〔}}\frac{{\left( {{a_3} - t} \right)}}{{{a_4}}}〕}} \end{aligned} $

(2)

式中：参数a₁, a₂, a₃, a₄为掌握时间序列曲线的宽度和幅度。利用D-L法先对时间序列曲线的局部进行拟合, 再结合完整的时间序列函数特征, 重新构建噪声较小的植被指数时间序列^[18]。

采用Jönsson和Eklundh^[19]改进后的动态阈值法提取植被物候参数，动态阈值法将NDVI增长(下降)幅度为振幅的某个百分比时对应的时间定义为植被SOS(EOS)^[19-20]。该方法的阈值不是具体的植被指数值，而是一个比值形式，即该像元的植被指数值与对应年份时间序列的振幅比^[21-22]。该阈值法具有较强灵活性和适应性^[23]，被广泛用于提取植被物候参数。本研究将同一年内NDVI增长(下降)幅度为振幅的20%定义为植被SOS(EOS)^[23]。

2.2 趋势分析和相关性分析 2.2.1 植被物候趋势分析

本研究利用Theil-Sen median趋势分析法描述物候的变化。对植被时间序列进行计算，以时间为自变量，每年的物候关键期作为因变量，变化趋势计算公式为：

$ \beta = {\rm{mean}}\left( {\frac{{{x_j} - {z_i}}}{{j - i}}} \right), {\forall _j} > i $

(3)

式中：β为趋势线的斜率; x_j，z_i表示第i年物候关键期的值；i为年变量。当β>0时, 表示物候关键期推迟或延长, β < 0时, 表示物候关键期提前或缩短^[24]。

Mann-Kendall检验法可用于检验长时间序列的显著性，其最大的优点是对数据误差有较强的抵抗能力，其显著水平的检验结果较为科学可行^[25]。采用M-K显著性检验方法(95%置信区间)来检验植被生长季开始日期(SOS)、生长季结束日期(EOS)和生长季长度(GSL)的年际变化显著性。

2.2.2 相关性分析

相关性分析用于定量分析自变量降水或气温与因变量物候期的相关性，寻找与物候参数相关性较为显著的控制因子，再逐个像元分析物候参数与气象因子之间的显著相关关系^[26]。其计算公式为：

$ r = \frac{{\sum\limits_{i = 1}^n {\left( {{x_i} - \bar x} \right)\left( {{y_i} - \bar y} \right)} }}{{\sqrt {\sum\limits_{i = 1}^n {{{\left( {{x_i} - \bar x} \right)}^2}} \sum\limits_{i = 1}^n {{{\left( {{y_i} - \bar y} \right)}^2}} } }} $

(4)

式中：y为气候因子某年的一个栅格数值；为气候因子多年均值；x_i为在第i年物候关键期的某个栅格数值；x为物候关键期多年均值。

相关系数r≤0表示负相关，r＞0表示正相关；利用F检验对相关系数的结果做显著性检验。

3 结果 3.1 植被物候时空变化 3.1.1 植被物候时间变化

由图 1可见，毛乌素沙地的植被SOS呈提前趋势。从2000年的111 d(4月中旬)提前到2019年的97 d(4月初)，提前了14 d，年均提前0.7 d；EOS总体变化不明显，呈微弱的提前趋势；生长季长度GSL表现为延长趋势，从2000年的201 d(6.7个月)延长到2019年的214 d(7.1个月)，共延长了13 d，年均延长0.65 d。

3.1.2 植被物候空间变化

由图 2可知，毛乌素沙地植被SOS显著提前(p＜0.05)地区面积为16.78%)，主要分布在毛乌素沙地的东部。大部分地区SOS表现为不显著提前(61.02%)。仅有21.81%的地区表现为不显著推迟趋势，主要分布在毛乌素沙地的西部。EOS的提前/推迟表现不明显。其中，不显著提前的区域占整个植被覆盖区的52.05%，不显著推迟的区域占整个植被覆盖区的44.49%。GSL受SOS和EOS影响，总体呈延长趋势。其中，显著延长(p＜0.05)区域面积为14%，主要分布在毛乌素沙地的东北区域和北部部分区域。大部分区域GSL表现为不显著延长趋势(61.56%)。24.44%的区域表现为不显著缩短趋势，主要分布在毛乌素沙地的南部和西部。

3.2 植被物候空间格局

由图 3可看出，毛乌素沙地的植被物候在空间上具有明显的东西差异。SOS时间范围为第90—160 d，即一年的3—5月上旬。空间上表现为由东到西逐渐推迟的趋势，在研究区的中部及东部，植被生长季开始较早，在研究区的西部，植被生长季开始较迟。SOS主要集中在第90—140 d，即一年中的3—4月下旬。其中，第100—110 d的像元数占总像元数的24.5%，第110—120 d的像元数占总像元数的19.8%。EOS范围为第290—340 d，即一年的9月上旬至11月上旬。

空间上与SOS的东西差异不同，表现为由北向南为逐渐推迟的趋势。北部地区植被生长季结束较早，南部植被生长季结束较迟。EOS主要集中在第300—330 d，即一年的10—11月。其中，第300—310 d的像元数占总像元数的12.5%，第310—330 d的像元数占总像元数的84.9%。GSL范围为第150—250 d，即5至8个月。在空间上，GSL与SOS的空间分布一致，表现为由东到西逐渐缩短的趋势。中部及东部地区植被生长季较长，西部地区植被生长季较短。GSL主要集中在第180—250 d，即6~8.3个月。其中，第200—210 d的像元数占总像元数的18.7%，第210—230 d的像元数占总像元数的22.8%。

3.3 毛乌素沙地植被物候与气候的关系 3.3.1 年均气温和年降水量的时间变化

由图 4可看出，2000—2019年，毛乌素沙地年平均气温呈显著上升趋势(p＜0.05)，年降水量呈微弱增加趋势，总体在250~450 mm之间波动。毛乌素沙地气候总体上表现为年平均气温显著上升，年降水量不显著增加的趋势^[27]。

3.3.2 毛乌素沙地植被物候与气候因子的关系

毛乌素沙地植被SOS、植被EOS与气候因子的相关系数空间分布如图 5所示。

由图 5可知，植被SOS与3月降水量呈较强的负相关关系，其中与3月降水量呈负相关的地区面积为83.85%，其中显著负相关面积为9.4%，主要分布在鄂托克旗、伊金霍洛旗、乌审旗、榆林、神木以及鄂托克前旗的北部等地区。正相关地区面积为16.15%，分布在毛乌素沙地南部的靖边、定边以及鄂托克前旗的南部等地区。

植被SOS与春季季前累计降水量呈负相关关系(90.81%)，其中呈显著负相关的区域面积为6.7%，主要分布在伊金霍洛旗、乌审旗、鄂托克前旗、榆林西北部、盐池、定边、靖边以及鄂托克旗的中部地区。呈正相关地区面积为9.19%，主要分布在鄂托克旗、伊金霍洛旗的北部小部分区域、榆林以及神木南部的小部分区域。植被SOS与3月平均温度呈负相关关系，表现为负相关关系的面积为60.01%，其中显著负相关面积为2.34%，呈正相关的地区面积为39.99%；与4月平均温度相关性较弱，其中呈正相关(51.95%)和负相关(48.05%)的地区面积基本持平；与4月降水量相关性较弱，呈负相关关系地区面积为61.70%，其中显著负相关面积为0.27%；与春季季前平均温度呈正相关关系(58.85%)，其中显著正相关面积为1.29%，呈负相关区域的面积为41.15%。植被EOS与9月降水量呈正相关关系，呈正相关关系面积为63.67%，其中显著正相关面积为4.32%，主要分布在盐池、鄂托克前旗、定边西北部、靖边、乌审旗南部。负相关面积为36.33%，分布在鄂托克旗，神木、伊金霍洛旗以及灵武地区。植被EOS与10月降水量呈负相关关系面积为57.66%，其中显著负相关面积为0.97%；与9月温度呈负相关关系面积为55.3%，其中显著负相关区域面积为3.26%；与10月温度呈正相关区域面积为50.5%，其中显著正相关区域面积为3.56%；与秋季季前平均温度呈正相关区域面积为54.95%；与秋季季前累计降水呈正相关区域面积为54.99%，呈负相关区域面积为45.06%。

4 讨论 4.1 植被物候时空变化

毛乌素沙地植被物候受水热条件控制，植被SOS主要呈提前趋势(77.8%)，其中显著提前区域为整个研究区的16.78%，这与韩红珠等^[24]的结论一致。在毛乌素沙地中部的乌审旗、鄂托克前旗以及东部的神木、榆林、北部的伊金霍洛旗，植被最先进入生长季，在鄂托克旗西部、灵武、定边南部地区等，植被进入生长季较晚，这与朱娅坤等^[25]、王静璞等^[26]的研究结果一致。SOS出现较早的地区多为草地，而出现较晚的地区多为耕地，且海拔较高，人为和地形因素的影响较大^[27]。毛乌素沙地EOS在2000—2019年变化不明显，空间上具有明显的南北差异，北部的植被结束生长季较早，一般在9月下旬至10月上旬。南部植被结束生长季较晚，一般在10月下旬至11月上旬。可能与纬度存在一定的负相关关系。植被GSL受SOS和EOS的影响，空间上表现为由东到西逐渐缩短。SOS的提前是GSL延长的决定因素。在毛乌素沙地东部乌审旗、鄂托克前旗、伊金霍洛旗、横山、神木、榆林和靖边等地的GSL为6~7个月，西部鄂托克旗、灵武、盐池和定边等地的GSL为5~6个月。由于MODIS NDVI数据长度为2000—2019年共20 a，这可能是大部分区域的植被物候变化趋势不明显的原因。而GIMMS NDVI数据有着较长的时间序列(1982—2015年)，是其他数据无法比拟的优点^[25]。缺点是空间分辨率较粗，无法进一步研究细节信息，未来可考虑将MODIS与GIMMS进行数据融合操作，将高时空分辨率的NDVI数据用于物候研究。毛乌素沙地的植被分布稀疏，植被信号较弱^[27]，植被稀疏使得地表裸露，土壤背景反射率较高，可能会影响对植被物候的判断^[28]。未来物候的研究需要地面物候数据的验证来确保结果的准确性^[29]。

4.2 植被物候对气候变化的响应

植被SOS提前受春季季前累计降水量和3月降水量影响较为显著，受春季季前平均温度、3月温度和4月温度影响较小。说明前一年冬季降雪的增多有利于植被生长季开始日期提前^[30]，不同月份的降水量对植被生长季开始日期的影响不同，3月降水量对SOS的影响较强，春季季前累计降水量和4月份降水量则次之。植被EOS提前受秋季季前累计降水量和9月降水量的影响较显著，受季前平均温度以及9月和10月平均温度影响较小。秋季季前累计降水量和9月降水量对植被EOS的推迟呈正相关关系，说明在生长季结束日期前后，降水量的增多有助于延缓植被生长季的结束^[31]，其中植被EOS与9月降水量相关性最强，原因可能是植被GSL主要集中在10月左右，植被对降水的响应具有一定的滞后性^[32]。

研究区的中部及东部地区为毛乌素沙地荒漠化较严重的区域，该区域植被SOS提前，GSL变长，且与春季季前累计降水量和3月降水量影相关性较强。说明在荒漠化区域，降水是影响植被生长的主导因素，这与王涛等^[28]、王旭洋等^[29]和Du Jun等^[30]的观点一致。Du Jun等^[30]发现，在干燥环境下，植被对温度的敏感度通常较低。在这种情况下，虽然适应植被生长的温度已经达到了要求，但是由于降水量低，土壤环境条件不理想，降雨成为植被生长的主要限制因素^{[31, 32]}。

本研究在探究植被物候与气候之间的响应关系时，只选取了平均温度和降水量数据，而实际上最高气温，最低气温和极端气候事件等也会对植被生长产生影响^[33]。另外，除了气象因素外，地形和人为因素的影响也不容忽视^[34]，在高海拔和低海拔，同一植被也会表现出不同的物候现象。大量的人工植树种草，人工灌溉都会对植被生长产生干预^[35]。因此未来研究可进行小区域间植被物候关键影响因素的探究。此外，植物生态系统之间的竞争、资源限制以及植被遗传因素也对植被物候产生了一定的影响^[35-36]。

5 结论

(1) 2000—2019年，毛乌素沙地植被SOS总体呈提前趋势，提前趋势不明显。显著提前面积为16.78%，年均提前0.4 d；植被EOS表现为微弱的提前趋势，年均提前0.3 d；植被GSL总体呈延长趋势，其中显著延长区域面积为14%，年均延长0.75 d。

(2) 空间尺度上，70.8%的毛乌素沙地植被SOS在第90—140 d，表现为由东到西逐渐推迟的趋势。88.9%的植被EOS范围为第290—340 d，由北向南逐渐推迟。72.6%植被GSL范围为第150—250 d，表现为由东到西逐渐缩短。

(3) 毛乌素沙地植被SOS与春季季前累计降水量和3月降水量相关性较为显著，其中与3月降水量的负相关关系最为显著。与4月降水量、春季季前平均温度、3月及4月温度相关性较弱；植被EOS与9月份降水量的正相关关系较强，与秋季季前平均温度、季前累计降水量、9月温度、10月温度、10月降水量的相关性都表现为较弱的相关性。