植物物候指植物受气候和其他环境因子(地形、水文和土壤等)影响而出现的、以年为周期的自然现象^[1], 对农事预报、指导农牧业生产和适应气候变化的生态保护政策制定具有重要意义^[2-3]。政府间气候变化专门委员会(IPCC)第5次评估报告报出, 全球气候系统的变暖是明确无疑的, 自1901年以来全球地表气温已经升高了0.89 ℃^[4]。全球变暖导致我国植物春季物候期发生显著变化^[5]。自20世纪80年代以来, 我国大部分植被的春季物候期主要表现为植物生长季始期呈现提前趋势, 特别是东北、华北和青藏高原等地区^[6], 而西南地区东部和长江中游等地区的物候期有所推迟^[7]。

青藏高原作为世界的第3极, 幅员辽阔, 由于气候条件和植被类型复杂多样、区域水热条件差异性大等因素影响, 不同学者对其植物物候变化持有不同观点。YU等^[8]发现气候变暖使青藏高原高寒草甸和高寒草原生长季推后, 生长期缩短; DING等^[2]认为青藏高原高寒草地物候年际变化在不同的海拔和自然带上分异显著, 植物生长季长度呈现中、西部缩短, 东部延长的趋势; ZHANG等^[9]发现植物生长季始期逐渐提前; SHEN等^[10]则认为没有证据表明青藏高原植物生长季始期已经提前。雅鲁藏布江流域位于青藏高原南部, 面积为24.2万km², 平均海拔4 600 m, 高差达7 000多m, 是典型的高海拔大流域^[11]。复杂多样的气候条件和丰富的生物多样性造就了雅鲁藏布江流域极其显著的植被垂直地带性^[12], 使其成为研究青藏高原气候变化生态响应的典型代表性区域^[13]。相关研究表明, 全球变暖导致雅鲁藏布江流域和大峡谷自然保护区植被绿度显著增加, 但植被生长率与升温速率随海拔梯度的变化趋势不一致; 降水变化是造成这种现象的一个重要因素, 但基于数据统计分析的解释并不能完全回答植被生长率随海拔变化的差异性, 亟待查明流域不同海拔梯度植物物候变化规律^[14]。

为此, 基于1999—2013年雅鲁藏布江流域植被指数(NDVI)和数字高程模型(DEM)数据, 分析植物返青期、枯黄期和生长季长度的时空变化规律, 研究流域不同海拔梯度的植物物候变化特征, 试图探析全球变暖影响下青藏高原高海拔地区植被生长与海拔梯度的响应关系。

1 研究区概况

雅鲁藏布江发源于我国喜马拉雅山脉北麓海拔约5 750 m的杰马央宗冰川, 自西藏自治区巴昔卡(海拔约150 m)出境后流入印度。流域海拔跨度约为150~7 000 m, 地势西部高、东部低, 其北侧为冈底斯山、念青唐古拉山和唐古拉山, 南侧为喜马拉雅山。流域内上、下游气候条件差异巨大, 上游地区受环流形势和水汽条件的限制, 为温带草原气候区, 无暴雨, 年降水量 < 300 mm, 发育着高寒草原、草甸和灌丛等。中游地区在米拉山的屏障作用下, 潮湿的水汽难以抵达西部地区, 属温带森林草原气候区, 年降水量介于300~600 mm之间, 发育着灌丛草原、亚高山灌丛草甸和亚高山草甸。下游地区属于高山峡谷区, 雅鲁藏布江在南迎巴瓦峰(海拔7 782 m)附近骤然由东流折向南流、再转向西南流, 形成世界上罕见的马蹄形大河, 该区为亚热带湿润气候区, 最大降水量达4 000 mm, 分布着少量的亚热带常绿阔叶林、山地热带雨林和季雨林等森林类型, 近下游为过渡性的针阔混交或常绿阔叶林^{[12, 14]}。

从水文特点来看, 流域可分为上游区、中游区和下游区, 上游为高原宽谷区(日喀则市拉孜县里孜村以上), 中游为宽窄相间河谷区(位于里孜村至林芝地区米林县派镇之间), 下游为峡谷(派镇以下)。上游海拔多在5 000 m以上, 下游多在3 000 m以下, 中游地处藏南宽谷和藏东南高山峡谷区, 海拔基本介于2 800~5 000 m之间。按河谷特征和行政边界, 自江源而下流域可划分为4个宽谷段, 即马泉河宽谷、日喀则宽谷、山南宽谷和米林宽谷, 面积分别占流域总面积的11.0%、33.3%、16.2%和39.5%。其中, 马泉河宽谷段位于上游地区, 日喀则宽谷段和山南宽谷段位于中游地区, 米林宽谷段位于中游和下游地区^{[11, 15]}。

2 数据与方法 2.1 数据

主要包括1999—2013年的NDVI数据和DEM数据。NDVI数据为SPOT-VGT NDVI数据集(来源于http://www.vgt.vito.be/), 空间分辨率为1 km×1 km, 时间分辨率为10 d, 1 a共有36期数据。该NDVI数据集在发布前已进行了几何校正、大气纠正、辐射校正和地形校正等数据预处理, 使用最大值合成法(MVC)消除了非植被因素对NDVI数据的不利影响, 具体见文献[14]。DEM数据为空间分辨率为1 km×1 km的SRTM DEM数据(来源于http://srtm.csi.cgiar.org/)。将DEM数据和NDVI数据集进行几何精校正, 投影变换统一为GCS-WGS-1984。

2.2 物候期的提取方法

基于遥感数据进行植物物候监测通常有6类方法:阈值法、拟合法、滑动平均法、最大斜率法、累积频率法和主成分分析法^[16]。其中, 阈值法是通过设定一组阈值提取植物物候期, 它是最早出现的一种方法, 应用较为广泛, 包括固定阈值法、动态阈值法和多参量阈值法; 动态阈值法采用的阈值不是某个具体的NDVI值, 而是一种动态比值, 该阈值在时间域和空间域上都具有更好的适应性^[17]。拟合法是近几年发展较快的物候提取方法, 它是利用平滑模型函数拟和时间序列遥感数据, 进而提取物候信息, 包括Logistic函数法、非对称高斯函数法(AG拟合法)和谐波函数法。AG拟合法是利用非高斯函数拟合曲线与比值阈值监测植物物候的开始和结束日期, 具有与原始NDVI曲线的整体特征更加吻合、最优保持原始NDVI质量的优点。该研究通过TIMESAT 3.0软件, 采用AG拟合法对时间分辨率为10 d的SPOT-VGT NDVI时间序列数据进行平滑处理^[18-19], 主要步骤如下。

步骤1(区间提取):提取NDVI时间序列数据区间的谷值和峰值。

步骤2(局部拟合):对位于谷值和峰值之间的NDVI时间序列数据进行2次局部拟合, 让最优化拟合函数能较好地描述NDVI的数据曲线, 公式为

$ f\left( t \right) \equiv f\left( {t;c, x} \right) = {c_1} + {c_2}g\left( {t;x} \right), $

(1)

$ \begin{array}{l} g\left( {t;x} \right) = g(t;{x_1}, {x_2}, {x_3}, {x_4}, {x_5})\\ \;\;\;\;\;\;\;\;\;\; = \left\{ \begin{array}{l} \exp \left[{-{{\left( {\frac{{t-{x_1}}}{{{x_2}}}} \right)}^{{x_3}}}} \right], t > {x_1}\\ \exp \left[{-{{\left( {\frac{{{x_1}-t}}{{{x_4}}}} \right)}^{{x_5}}}} \right], t < {x_1} \end{array} \right.。\end{array} $

(2)

式(1)~(2)中, g(t; x)为高斯函数; c₁和c₂为定义基线和振幅的线性参数; x₁为对应时间变量t的峰值或谷值的位置参数; x₂和x₃为决定函数曲线右半部分的宽度和陡峭度的参数; x₄和x₅为决定函数曲线左半部分的宽度和陡峭度的参数。

步骤3(整体连接):在局部拟合函数的基础上进行整体拟合, 描述整个NDVI变化过程, 公式为

$ F\left( t \right) = \left\{ \begin{array}{l} \alpha \left( t \right){f_{\rm{L}}}\left( t \right) + \left[ {1 - \alpha \left( t \right)} \right]{f_{\rm{L}}}\left( t \right),{t_{\rm{L}}} < t < {t_{\rm{C}}}\\ \beta \left( t \right){f_{\rm{C}}}\left( t \right) + \left[ {1 - \alpha \left( t \right)} \right]{f_{\rm{R}}}\left( t \right),{t_{\rm{C}}} < t < {t_{\rm{R}}} \end{array} \right.。$

(3)

式(3)中, [t_L, t_R]为NDVI的总体变化区间, 包含至少1个左谷值、峰值和右谷值的时间区间; f_L(t)、f_C(t)、f_R(t)分别为区间[t_L, t_R]内左谷值、峰值和右谷值所对应的局部拟合函数; α(t)和β(t)为位于[0, 1]之间的剪切系数。

通过式(1)~(3)得到流域植物物候的时间序列曲线, 拟合后曲线更接近真实情况, 避免了整体数据对局部拟合的干扰。反映植物物候的参数包括生长季始期(又称返青期, the start date of growing season, SOG)、生长季末期(又称枯黄期, the end date of the growing season, EOG)和生长季长度(the length of the growing season, LOG), 其中LOG为EOG与SOG之差。采用JÖNSSON等^[19]改进后的阈值法提取物候参数, 即首先利用AG拟合法得到时间序列曲线, 再使用动态阈值法提取SOG、EOG和LOG。SOG的阈值为时间序列曲线的左谷值到峰值之间的值, EOG的阈值为峰值到右谷值之间的值。阈值的设定对SOG、EOG和LOG的提取至关重要, 但动态阈值的设定要综合考虑区域植被类型的差异性。

根据《中华人民共和国植被图》分类系统, 雅鲁藏布江流域包括针叶林(占流域面积的比例约为7%)、阔叶林(4.5%)、灌丛(5.8%)、草原(7%)、草甸(33.5%)和高山植被(19.2%)等11个植被型组、21个植被型。其中, 针叶林主要为亚热带和热带山地针叶林, 阔叶林包括亚热带落叶阔叶林、亚热带硬叶常绿阔叶林和热带雨林等。森林植被面积约占流域总面积的11.5%^[12]。JÖNSSON等^[18]在提出动态阈值法时建议SOG和EOG阈值为NDVI达到年振幅的20%左右。据此, 宋春桥等^[17]在进行藏北高寒草原物候参数提取时, 将SOG阈值设定为10%, 将EOG阈值设定为20%;李明等^[20]在进行森林植被为主地区的物候参数提取时, 将SOG和EOG阈值设定为50%。笔者在综合考虑流域主要植被类型的基础上, 结合相关文献研究成果, 将物候SOG和EOG阈值均设定为NDVI达到年振幅的20%, 以提取流域植物物候参数(图 1)。

2.3 分析方法

基于各年的NDVI时间序列数据, 利用动态阈值法提取植物的SOG、LOG和EOG参数, 采用线性回归分析方法, 计算1999—2013年流域每个像元的植物SOG、LOG和EOG的年际变化趋势, 回归模型为

$ y = at + b。$

(4)

式(4)中, y为SOG、LOG或EOG值; t为年份(1999—2013年); a为变化率; b为截距。

基于1 km×1 km空间分辨率的流域DEM数据, 以500 m为间隔, 将流域高程划分为12个海拔梯度(中心点高程分别为500, 1 000, 1 500, …, 6 000 m), 每个海拔梯度为该中心点高程加减250 m的高程区间(例如海拔梯度500 m即为高程250~750 m, 海拔梯度1 000 m即为高程750~1 250 m, 依次类推)。基于海拔梯度划分, 绘制流域LOG多年平均值, 以及SOG、EOG和LOG随海拔梯度变化的箱线图(the box-and-whisker plot), 以分析不同海拔梯度植物物候的变化特征。其中, 箱体的高低表示四分位距(interquartile range, IQR)大小, 即第3四分位数(Q3)和第1四分位数(Q1)之差, 反映各物候参数(SOG、EOG和LOG)的集中范围。

3 结果与分析 3.1 1999—2013年流域植物物候的年际变化

由图 2可见, 近15 a雅鲁藏布江流域SOG提前的区域占流域总面积的61.29%, 年变化率为>0~5 d·a^-1的区域最多(59.3%), 年变化率为>5~20和>20 d·a^-1的区域分别占2.0%和0.03%。

流域SOG推迟的区域占流域总面积的38.71%, 以年变化率为>-5~0 d·a^-1的区域较多(35.5%), 年变化率为-20~-5和 < -20 d·a^-1的区域分别占3.1%和0.12%。SOG年变化率为>0~5 d·a^-1的区域主要分布在山南宽谷段、日喀则宽谷段以及米林宽谷段较高海拔地区; SOG年变化率为>-5~0 d·a^-1的区域主要分布在马泉河宽谷段以及日喀则宽谷段的较高海拔地区、米林宽谷段的较低海拔地区。整个流域SOG年变化率的像元平均值为0.29 d·a^-1, 返青期略有提前。

流域EOG提前的区域占流域总面积的45.3%, 以年变化率为>0~5 d·a^-1的区域较多(42%), 年变化率为>5~20和> 20 d·a^-1的区域分别占3.3%和0.02%。流域EOG推迟的区域占流域总面积的54.7%, 以年变化率为>-5~0 d·a^-1的区域较多(50.50%), 年变化率为-20~0和 < -20 d·a^-1的区域分别占4.2%和0.05%。EOG年变化率为>0~5和>-5~0 d·a^-1的区域空间分布与其对应的SOG年变化率为>-5~0和>0~5 d·a^-1的区域相似。整个流域EOG年变化率的像元平均值为-0.18 d·a^-1, 枯黄期略有延迟。

流域LOG延长的区域占流域总面积的40.9%, 以年变化率为>0~5 d·a^-1的区域较多(32.5%), 年变化率为>5~20和>20 d·a^-1的区域分别占8.4%和0.03%。流域LOG缩短的区域占流域总面积的59.14%, 以年变化率为>-5~0 d·a^-1的区域较多(48.9%), 年变化率为-20~0和 < -20 d·a^-1的区域分别占10.2%和0.01%。流域LOG年变化率的像元平均值为-0.47 d·a^-1, 植物生长季长度存在缩短现象。这表明流域植物物候变化存在较大的空间异质性现象, 返青期提前的区域枯黄期不一定推迟, 返青期推迟的区域枯黄期也不一定提前。

3.2 植物物候随海拔梯度的变化特征

1999—2013年雅鲁藏布江流域植物LOG的多年平均值总体上呈现随海拔增加而减小的趋势(图 3)。自海拔500 m开始增加, 在海拔1 000 m时达最大, 其后逐渐减小。流域SOG年变化率介于-10~10 d·a^-1之间, 随海拔梯度的增加呈现明显的先减小再增加趋势。从四分位距来看, SOG年变化率自海拔500~1 500 m呈现减小趋势, 在>1 500~3 500 m总体呈增加趋势(在3 000 m处达最大)。就中位数而言, 随海拔梯度的增加呈现先减小再增加后趋于平稳的变化趋势, 总体上表现为海拔梯度 < 3 000 m时中位数为负值, 其中海拔1 000 ~2 000 m之间的中位数较小; 在海拔> 3 000 m时, 中位数为正值。

流域EOG的年变化率基本介于-5~10 d·a^-1之间, 随海拔梯度的增加总体呈减小趋势。从四分位距来看, EOG年变化率随海拔梯度的增加呈略微减的趋势, 在海拔3 500~4 500 m之间稍低。就中位数而言, 随海拔梯度的增加呈现先减小后趋于平稳的变化趋势, 所在海拔梯度的中位数值均大于0, 海拔超过3 500 m以后基本保持不变。

流域LOG年变化率基本介于-10~15 d·a^-1之间, 表现为随海拔梯度增加而先减小再增加的变化趋势。从四分位距来看, LOG年变化率除在海拔3 000 m处较大外, 其余海拔梯度基本保持不变。就中位数而言, 在海拔500 ~2 500 m基本保持不变(其值在5 d·a^-1左右), 自海拔高度达2 500 m以后, 其随海拔梯度增加呈现先减小后平稳的趋势, 海拔>3 500 m后中位数保持稳定, 其值小于0。总体而言, 随着海拔梯度的增加, 雅鲁藏布江流域植物LOG及其年变化率逐渐减小。

4 讨论

有证据表明, 青藏高原的升温速率随海拔增加而增大^[21], 造成不同海拔梯度高山带植物多样性的变化^[11]和植被生长率的差异性^[14]。笔者以雅鲁藏布江流域为例研究发现, 随着海拔梯度的增加, 植物生长季长度和年变化率总体上呈减小趋势, 该结论与DING等^[2]和朴世龙等^[22]指出的海拔在物候的地域分异中扮演重要作用相一致。温度被认为是影响植物物候最重要的气象因子^[23], 地表气温不仅随着海拔的增加而降低〔速率为-0.65 ℃·(100 m)^-1〕, 而且相关研究表明, 高海拔地区的升温速率要大于低海拔地区^[24-25]。雅鲁藏布江流域不同宽谷段植物物候呈现一定的空间差异性, 日喀则宽谷段和山南宽谷段呈现缩短现象, 海拔较高的马泉河宽谷段和海拔较低的米林宽谷段呈现延长现象; 部分地区植物生长季的延长, 可能是气温升高和年降水量变化对特定植被类型综合作用的结果^{[12, 14]}。该研究试图查明高海拔大流域植物物候的变化特征及其在不同海拔梯度的分异性, 然而, 雅鲁藏布江流域植被类型丰富, 垂直地带性显著, 植物物候变化的影响因素极其复杂。笔者虽基于遥感数据进行了一些有益的探索, 仍存在一定的不足之处。一是未分析不同植被类型物候变化的差异性; 二是即使在相同水热条件下, 特别是在高寒地带, 不同植物的萌动发芽和物候特征亦不尽相同。

地面观测是查明植物物候变化的最佳手段^[7]。面对幅员辽阔的青藏高原, 目前仍缺乏足够的实地观测站点, 并且历史观测数据无法获取等, 阻碍了对物候变化遥感监测的有效验证, 这些科学问题有待下一步认真探索。“十三五”期间, 西藏地区将通过设置主要生态系统综合监测大样地和优势种群辅助监测样点, 建成山南、林芝和日喀则等10个生态监测站, 开展生态安全屏障功能监测与评估, 指标涉及返青期、开花期和枯黄期等观测内容, 可以为验证植物物候提取所需的动态阈值设定提供依据。该研究基于遥感手段和野外调查开展分析, 是在生态监测站正式运行之前的一次有益尝试, 以期科学认知高海拔大流域植物物候变化的时空分异特征, 为开展青藏高原关键区域气候变化生态风险预警提供依据。

5 结论

(1) 雅鲁藏布江流域植物生长季长度存在缩短现象, 年变化率的平均值为-0.47 d·a^-1, 植物物候变化存在较大的空间异质性现象。返青期提前的区域, 枯黄期不一定推迟; 返青期推迟的区域, 枯黄期也不一定提前。雅鲁藏布江流域植物的返青期提前的区域占流域的61.3%, 枯黄期提前的区域占流域的45.27%。

(2) 随着海拔梯度的增加, 雅鲁藏布江流域植物生长季长度和年变化率总体上在逐渐减小, 返青期提前和枯黄期推迟的趋势愈发明显, 生长季长度呈现缩短趋势。部分地区植物生长季的延长, 可能是气温升高和年降水量变化对特定植被类型综合作用的结果。

(3) 雅鲁藏布江流域是青藏高原典型的高海拔大流域，基于该区域研究发现的不同海拔梯度植物物候的变化规律，不仅对“十三五”期间青藏高原生态保护与修复工程的实施具有一定的参考价值，而且对基于水土保持与防风固沙、生物多样性维护、水源涵养等不同主导生态功能建设的生态监测站的观测样地布设和监测指标选取方面具有重要的实践价值。