2003—2018年米仓山地区植被物候时空变化及对气候的响应

邵周玲,周文佐,李 凤,周新尧,杨 帆

西南大学地理科学学院, 遥感大数据应用重庆市工程研究中心,重庆 400715

物候是生物长期适应自然环境周期性变化,形成与此相适应的生长发育节律的现象,对研究植物和气候关系具有重要意义[1- 3]。第五次IPCC报告指出在过去的130年间全球温度上升超过了0.85℃[4],全球气候变暖会对陆地生态系统在季节和年际尺度上产生影响,如植物在不同时期出现发芽、展叶、开花、结果、枯叶等物候现象,这些现象对环境变化敏感且容易观测,因此物候通常作为气候变化的生物感应器[5- 6]。植物生长季的变化对全球碳的固定、CO2浓度变化以及全球水气循环都有重要影响,揭示植物物候变化成为物候学当前的热点研究问题[7]。

在过去数十年中,受气候变暖影响全球普遍出现不同程度的春季物候提前、秋季物候推迟以及生长季延长的现象,不同区域物候变化受自然环境差异影响存在明显的差异,气候变化对植物物候沿纬向和垂直方向上的分布变化都有显著影响[8-10]。我国学者对于植被物候的研究主要集中在中国的北温带以及高海拔地区[11- 12],其中东北地区[13- 15]、青藏高原[16- 19]以及内蒙古[20]是主要的研究区域。秦巴山区作为我国大陆南北气候过渡带和重要生态廊道,同时也是气候变化敏感区和生态脆弱区,一直以来都是物候研究的重点区域。目前,已有学者在秦岭做了相关研究并取得一定成果。王钊等[21]、Xia等[22]对秦岭物候的研究发现生长季始期和末期整体上分别出现提前和推迟,并随着海拔升高生长季始期、末期和长度分别出现推迟、提前和延长趋势；邓晨晖等[23]发现秦岭北坡的生长季始期较南坡早,生长季末期南坡则较北坡晚,在海拔>2700m的高海拔区,生长季始期、末期及生长季长度随海拔上升的波动幅度更强；张晓东等[24]研究结果显示2000—2015年秦岭东部伏牛山地区生长季始期和生长季末期总体呈推迟变化。气温和降水被认为是影响物候变化最直接的因子[25]。已有研究结果表明气温是影响秦岭物候变化最重要的因素,特别是同期气温升高对物候始期的提前以及末期的推迟具有主导控制作用[21,26- 27]。

虽然目前对秦巴山区的物候研究已有初步的结论,但是由于区域植被类型丰富,垂直地带性差异显著,气候变化对于秦巴山地不同垂直带的植被物候的影响还缺乏深入的研究。本文基于2003—2018年MODIS NDVI数据分析秦巴山区典型山地-米仓山地区的植被物候时空格局变化,并且对不同垂直带、不同植被类型的物候时空特征进行研究,最后结合气象资料分析物候变化与主要气候因子的关系,为评价全球变化背景下米仓山地区植被对于气候变化的响应提供科学依据。

1 研究区概况

米仓山地区(30°50′—33°16′N,105°50′—109°46′E)位于秦巴山区的南部,地处陕西省、四川省和重庆市交界处,属于大巴山的西段,呈西北-东南走向,海拔在200—2500m(图1)。地处中国的北亚热带,受东南季风和西南季风的共同影响,年均温15℃左右,年降水量达1200mm以上,气候温和湿润,四季分明。米仓山地区植被类型丰富,主要有常绿阔叶林、落叶阔叶林、常绿针叶林、落叶阔叶灌木林、常绿阔叶灌木林、针阔混交林以及草地。

图1 研究区地形与主要植被类型分布Fig.1 The landform and main vegetation types of study area

2 数据来源与研究方法

2.1 数据来源与预处理

本文中植被指数数据采用美国国家航空航天局(https://earthdata.nasa.gov)的16d合成250m MYD13Q1数据产品。采用官方提供的MRT(MODIS Reprojection Tool)对数据进行预处理,最后裁切得到研究区2003—2018年的NDVI数据。气候数据来源于中国气象科学数据中心(http://data.cma.cn),包括2003—2018年米仓山地区及其附近共28个气象站点的平均气温和降水量日数据。先将气温和降水数据整理为月均温、月累计降水量,再采用ANUSPLIN工具将气象数据插值为250m的栅格数据。地形数据为30 m分辨率的ASTER GDEM 数据(http://www.gscloud.cn),将其重采样为250m分辨率的DEM数据。植被数据为来源于欧空局的300m分辨率的陆地覆盖数据(https://www.esa.int),通过重采样将分辨率统一为250m。

2.2 研究方法

2.2.1物候参数提取

利用TIMESAT软件中的S-G滤波(Savitzky-Golay filtering)方法[28]对数据进行滤波平滑,去除时间序列中的随机噪声点,以便于物候参数提取。

(1)

式中,Yj为拟合后的NDVI,Y代表原始序列数据,Ci为第i个NDVI值的滤波系数,N为平滑窗口的大小。

采取动态阈值法对米仓山地区植被物候期进行提取,参考伏牛山、秦岭等附近地区的物候研究[23- 24],将阈值设置为50%,得到三个关键物候参数:生长季始期(Start of the growing season, SOS)、生长季末期(End of the growing season, EOS)、生长季长度(Length of growing season, LOS)。采用年序日(Day of year, DOY)表示物候参数,即距当年1月1日的实际天数。SOS受NDVI的影响很大,植被稀少生长季不可能开始,提取物候时应首先去除数据集中NDVI<0.1的这部分无效数据[29]。

NDVI(SOS,EOS)=(NDVImax-NDVImin)×50%

(2)

式中,NDVImax为一年中植被指数的最大值；NDVImin为植被指数上升和下降阶段的最小值；当NDVI上升或者下降达到变化振幅的50%时分别提取SOS和EOS。

2.2.2趋势分析

采用Theil Sen(T-Sen)斜率分析米仓山地区植被物候变化趋势。T-Sen斜率是一种非参数趋势度计算方法,通过计算长时间序列数据中任意两个时段的斜率,并选取中位数来表征该时段的变化趋势。T-Sen>0,则表明物候变化趋势是推迟或延长；T-Sen<0,则表明物候变化趋势是提前或缩短。T-Sen计算结果结合Mann-Kendall(M-K)统计检验(置信度为90%、95%、99%)将物候变化趋势划分为以下7类:极显著提前或缩短(T-Sen<0,P<0.01)、显著提前或缩短(T-Sen<0,0.010,P<0.01)、显著推迟或延长(T-Sen>0,0.010,0.050.1)。

2.2.3偏相关分析

分析物候与气温、降水的关系时采用偏相关系数来描述。多元相关分析中,简单相关分析可能不能够真实反映两个变量之间的关系,偏相关是在排除其他变量的影响后计算两个变量之间的关系。其表达式如下:

(3)

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

3 结果分析

3.1 植被多年平均物候期

为了研究米仓山地区物候总体状况的空间格局,通过ArcGIS10.2计算出16a的SOS、EOS、LOS均值(图2)。生长季始期主要集中在第80—110d(3月下旬到4月中旬),占总面积的 86.75%；生长季始期在第70—90d的区域主要分布在米仓山东北部,即汉江流域以南低海拔地区；生长季始期为第110—125d的区域主要集中在米仓山南部和东部边缘,并且南部受地形影响呈纵列分布；此外,在一些河谷地带,由于地势较低,少许植被生长季始期早于第70d,而在中部和南部高海拔地带植被生长季始期较晚。生长季末期主要集中在第250—300d(9月中旬至10月下旬),占总面积的95.73%；其中北部和中部的东缘部分区域最早,主要集中在第270d之前；在东北部边缘、中部以及南部边缘低海拔地区生长季末期最晚,主要集中在第290—300d；其他区域集中在第270—290d。生长季长度主要集中在130—210d,占区域面积的96.13%；其中长于195d的区域主要分布在东北部,与生长季始期在第90d以前的分布区域大致相同,主要分布在低海拔河谷平地区域；生长季长度短于130d的区域呈点团状分布,除在南部呈点状零星散布以外,在嘉陵江以东和汉江以南区域也呈团状聚集分布,并且与生长季末期在第250d以前的区域分布大致相同。从水平方向上,由于米仓山地区跨越的纬度范围较小,主要受到地形和气候的影响,其物候变化规律不明显；垂直方向上,河谷、低山地区由于海拔较低、气温较高,植被生长季始期较早,生长季长度较长；高山区域由于海拔较高、气温较低,植被生长季末期较早,生长季长度较短。

图2 米仓山地区2003—2018年物候均值空间分布Fig.2 Spatial pattern of multi-annual mean phenophase in Micang Mountains from 2003 to 2018

米仓山地区物候期与海拔关系密切,利用2013—2018年SOS、EOS、LOS多年均值,选取50m高程为间隔统计得到植被物候随海拔变化趋势(图3)。由于600m以下区域主要为耕地,物候变化受人类活动影响较大,SOS波动幅度较大,无明显规律；600m以上区域SOS随着海拔升高而推迟。SOS随海拔变化的一元线性斜率为0.006,即海拔每上升100m,SOS大约推迟0.6d。EOS与海拔的关系较为复杂, 1000m以下区域随海拔变化物候规律不明显；1000—2000m区域EOS随海拔升高迅速提前；2000m以上呈现波动延后现象。LOS随海拔升高的变化规律与EOS大致相同,存在2000m分界线。

图3 不同海拔上物候参数的变化情况Fig.3 Variations of phenological parameter in different elevations

3.2 物候年际变化趋势分析

利用2003—2018年米仓山地区各年份物候均值,得到SOS、EOS、LOS年际变化趋势(图4)。由图可知,年均物候期随时间均呈波动变化,整体上SOS有提前趋势,提前幅度为0.47d/a,未通过显著性检验(P=0.20),EOS提前0.22d/a,LOS延长0.25d/a,均通过0.01显著性检验。

图4 2003—2018年物候年际变化趋势Fig.4 Inter-annual trends of phenophase from 2003 to 2018

基于T-Sen趋势分析方法在像元尺度上分析米仓山地区2003—2018年植被物候变化趋势,并结合Mann-Kendall统计法对结果进行显著性检验(图5)。SOS拟合斜率值介于-12.56—12.68之间,斜率为负的区域面积占比为74.03%,表明整体上米仓山地区植被生长季始期呈提前趋势。其中,极显著提前、显著提前、弱显著提前像元分别占1.97%、6.27%、6.77%,主要在米仓山地区的重庆市境内以及汉中市南部分散分布。达到显著性(P<0.1)延后的区域较少,极显著延后、显著延后、弱显著延后的像元共占0.85%。从达到显著变化的区域分布来看,显著性提前的区域约是显著性延后区域像元的18倍。EOS拟合斜率值介于-23.18—22.33之间,斜率为负的像元占比为60.67%,表明大部分区域生长季末期呈提前趋势。其中,显著提前区域主要在汉中市、安康市、达州市的中部及北部、重庆市的东北部以及巴中市的中北部。生长季末期显著延后的区域主要分布在达州市、巴中市、重庆市的南部以及广元市。LOS拟合斜率值介于-16.93—15.25之间,斜率为正的像元占比为53.25%,表明米仓山植被生长季长度整体呈延长趋势。其中,极显著延长、显著延长、弱显著延长像元占比分别为2.21%、2.17%、0.59%,主要在重庆市、广元市、巴中市和达州市等地区分布。达到显著性(P<0.1)缩短的区域较少,主要分布在海拔较高的地区,如汉中市和巴中市的交界处,以及重庆市东北部。极显著缩短、显著缩短、弱显著缩短的像元共占2.33%。

图5 2003—2018年米仓山地区植被物候期年际变化空间分布Fig.5 Spatial distribution of inter-annual variation of phenophase in Micang Mountains from 2003 to 2018

3.3 不同植被类型物候空间分布及变化趋势

3.3.1不同植被类型的物候期垂直变化

为了更加深入地探究植被物候的规律,将米仓山地区海拔以500m为间隔,划分为5个梯度带(<500、500—1000、1000—1500、1500—2000、>2000m)。根据米仓山地区植被垂直分布,进一步分析主要植被类型的SOS、EOS、LOS在各垂直带上的物候特征(图6)。米仓山地区植被主要分布在500—1500m,占总面积的56.59%,2000m以上占1.14%。主要的植被类型为常绿阔叶林和落叶阔叶林,占总面积的40.65%,其次是常绿针叶林和落叶阔叶灌木林,分别占总面积的9.38%和7.04%。

图6 各垂直带上不同植被类型面积占比Fig.6 Area percentage of different vegetation types in different vertical belts

不同植被类型在各级垂直带的物候期统计结果如表1。随着海拔上升,不同植被均出现了SOS推迟的现象,但是不同植被变化趋势略有区别。落叶阔叶林SOS最早,早于同级垂直带其他植被类型SOS均值,而且落叶阔叶林SOS随海拔升高推迟幅度最大,从<500m(94.2d)到>2000m(106.9d)的垂直带上推迟了17.9d。草地是SOS最迟的植被类型,迟于同级垂直带其他植被类型SOS均值。常绿型森林植被(常绿阔叶灌木林、常绿阔叶林、常绿针叶林)SOS迟于落叶型森林植被(落叶阔叶灌木林、落叶阔叶林),常绿型森林植被SOS均值在500—1000、1000—1500、1500—2000、>2000m垂直带上晚于落叶型森林分别为9.6、7.3、3.7、2.1d。500m以上区域,随着海拔升高,不同植被间的SOS差距缩小。

表1 不同植被类型在各级垂直带的物候期/d

随着垂直带的海拔升高各植被类型EOS均出现提前,草地提前幅度最小,为16.5d；常绿阔叶灌木林提前幅度最大,为36.6d。常绿型森林植被与落叶型森林植被间的EOS差异不明显；在500m以上区域,不同植被类型间EOS差异增加,分别在500—1000、1000—1500、1500—2000、>2000m垂直带上各植被类型间EOS极差为4.6、7.3、16.3、20.9d,即在2000m以上区域植被间的EOS差异最为明显。

随着垂直带的海拔升高各植被类型LOS缩短,缩短幅度为27.9d—47.6d,即高、低海拔之间的LOS相差一个多月。随着垂直带升高,草地LOS缩短幅度最小(27.9d),常绿阔叶林缩短幅度最大(47.6d)。常绿型森林植被LOS普遍短于落叶型森林,在500—1000、1000—1500、1500—2000、>2000m垂直带上它们之间的差值分别为11.2、7.9、6.2、5.0d,主要由于常绿型森林植被的SOS晚于落叶型森林植被。

3.3.2不同植被类型的物候期时间变化

将米仓山地区2003—2018年不同植被的物候期(SOS、EOS)进行线性回归分析并做显著性检验,得到16a物候期的变化趋势及显著性(图7)。草地、常绿阔叶灌木林的生长季始期提前趋势最明显,提前幅度分别为-0.80、-0.71d/a、其中常绿阔叶灌木林达到显著提前(P<0.05)趋势；除此之外,常绿针叶林和常绿阔叶林也呈提前趋势,提前速度分别为0.49、0.43d/a。米仓山地区植被生长季末期均呈提前趋势,且通过P<0.05显著性检验,其中落叶阔叶林和针阔混交林生长季末期提前明显,年际变率分别为-0.46、-0.57d/a。总体上,生长季始期相比生长季末期提前趋势更明显,生长季长度为生长季末期与始期之差,因此生长季长度延长主要由于生长季始期的提前。

图7 不同植被生长季始期和末期的年际变化Fig.7 Inter-annual change trends of the starting and ending dates of different vegetation growth seasons in Micang Mountains

3.4 物候与气候因子间的相互关系

3.4.1生长季始期对气温、降水的响应

生长季始期与3—5月的平均气温和降水逐像元做偏相关分析(P<0.05),得到其偏相关系数空间分布图(图8)。与气温的偏相关分析结果显示,2003—2018年生长季始期与3—5月平均气温的偏相关系数均多为负值,表明温度升高生长季始期提前。生长季始期与3月平均气温呈负相关的区域占比为77.87%,多分布于汉江以南低海拔地区以及研究区中部；生长季始期与4月平均气温呈负相关的区域占67.59%,多分布于研究区的东部。与降水的偏相关分析结果显示,生长季始期与3月和5月降水负相关系数像元占比分别为55.08%、53.28%,在空间上正相关和负相关分布较为均匀。生长季始期与4月降水呈负相关的区域占65.77%,表明4月降水增加生长季始期提前。在显著性水平上(P<0.05),生长季始期与3月平均气温的相关性显著的像元最多(11.08%),其次是4月降水(9.42%),表明在大部分地区生长季始期主要受3月的平均气温以及4月降水的影响。

图8 生长季始期与3—5月气温/降水的偏相关系数空间分布Fig.8 Spatial distribution of partial correlation coefficients between SOS and temperature/precipitation in March, April and May

3.4.2生长季末期对气温、降水的响应

生长季末期与8—10月的平均气温和降水逐像元做偏相关分析,得到偏相关系数空间分布图(图9)。与气温的偏相关分析结果显示,2003—2018年生长季末期与8月和9月平均气温的正负相关性较为平衡,其中呈正相关的像元占比分别为49.84%、42.58%,生长季末期与8月平均气温呈正相关区域主要分布于米仓山地区南部,而与9月平均气温呈正相关区域主要在米仓山地区北部,到10月生长季末期与平均温度的正相关区域占比有所增加,为59.73%。与降水的偏相关分析结果显示,生长季末期与8月和9月降水的正负相关性相差较小,正相关像元占比分别为51.14%、57.59%。生长季末期与10月降水呈负相关的像元占88.41%,表明大部分区域随着降水增加生长季末期提前。在显著性水平上(P<0.05),生长季末期与8—10月平均气温的相关性显著的像元占比较少,分别为1.29%、2.97%、3.24%。生长季末期与8、9月的降水相关性显著的像元占2.42%和4.24%,而与10月降水相关性显著的像元最多(20.64%),多分布于米仓山地区的北部汉江流域、中部的达州市和巴中市以及东部的安康市。表明在米仓山地区生长季末期主要受10月降水的负向影响。

图9 生长季末期与8—10月气温/降水的偏相关系数空间分布Fig.9 Spatial distribution of partial correlation coefficients between EOS and temperature/precipitation in August, September and October

4 讨论

米仓山地区植被物候特征差异主要取决于植被类型、高程、气候等因素。不同植被类型生长习性不同,物候差异较大。然而,高程差异带来最明显的影响是水热分配不均,从而导致不同海拔物候差异。加之全球变暖的大背景下,中纬度地区温度上升,导致物候始期提前[11]。当前已有研究利用不同的数据和方法对秦巴山区中秦岭的植被物候进行了分析,大多数研究结论表明区域植被物候存在SOS提前、EOS推迟和LOS延长现象,并且在垂直方向上物候期也呈现规律性变化[30-32]。秦岭森林物候期结果表明海拔每升高100m,SOS推迟2d,2200m以下区域EOS提前明显,2200m以上区域变化平缓,LOS也有同样的现象[32]。就整体来看,SOS和LOS变化趋势与秦岭的物候研究保持一致,但是在米仓山地区EOS存在提前趋势。本文发现在米仓山地区物候期对海拔的响应存在2000m界限,在海拔2000m以上区域EOS和LOS变化趋势相对较小,而其下区域随海拔升高出现明显的EOS推迟和LOS缩短现象。探究其原因可能由于温度随着海拔升高以0.65℃/100m下降,降水量在一定高度内随着海拔升高而增加,而到达最大降水高度后降水量减少,造成EOS和LOS变化趋势发生转折。2003—2018年米仓山地区不同植被SOS 提前趋势不同,其中草地提前幅度最大,说明区域内草地对气候变化最敏感,这与丛楠等[33]在北半球中高纬度地区的研究结果相同。

物候始期对气候变化的响应程度比末期更显著,在亚洲中低纬度区域较少受到水分胁迫的影响,春季气温升高成为SOS提前的主要原因[33]。相比于物候始期,温度对末期影响较弱,主要由于EOS受多种因子协同影响。已有研究结果表明,植被的SOS主要受温度控制,而EOS则主要受降水的影响[14,34]。例如,张晓东等[24]研究表明秦岭东部伏牛山植被EOS主要受9月降水的影响,Xie等[35]发现降水增加导致秋季物候不显著提前。本文发现影响米仓山地区SOS提前最主要的因素是3月气温,其次是4月降水；而EOS提前主要受10月降水的影响,其原因可能是当地秋末时节温度较低,降水增加加剧低温效应,从而导致EOS提前。

北半球春季温度升高,生长季始期提前以及生长季延长可能造成植被生产力增加。植被物候能反映气候变化,不同区域的地理条件差异较大,不同的植被类型的物候对区域气候变化的响应也不一样。

5 结论

本文基于2003—2018年MODIS NDVI数据,采用S-G滤波和动态阈值法提取了米仓山地区物候参数,分析了米仓山地区物候参数时空分布和变化趋势,并讨论了物候变化与水热条件的相互关系。2003—2018年米仓山地区植被SOS出现在3月下旬到4月中旬,海拔每上升100m,SOS大约推迟0.6d；EOS主要发生于9月中旬至10月下旬；LOS主要集中在130—210d。除低海拔地区受人类活动影响物候波动较大外,EOS和LOS随海拔变化存在2000m分界线,其下随海拔升高物候EOS明显推迟或LOS缩短,而其上物候变化趋于平缓。米仓山地区16a来SOS呈提前趋势,提前速度为0.47d/a,提前区域占74.03%；EOS提前幅度为 0.22d/a,通过显著性检验(P<0.05)区域仅有4.9%；LOS略有延长,延长幅度为0.26d/a,呈延长趋势区域占53.25%。

米仓山地区不同植被物候期差异明显,常绿型森林植被物候始期晚于落叶型森林植被；草地、常绿阔叶灌木林SOS提前趋势最明显,提前幅度分别为-0.80d/a、-0.71d/a,相对于其他植被类型对气候变化更加敏感；EOS提前趋势最明显的是针阔混交林和落叶阔叶林。米仓山地区SOS主要受3月平均气温和4月降水的影响,3月平均气温升高以及4月降水增加导致物候始期提前；EOS提前主要受10月降水的影响。