气候变化已影响到陆地生态系统的生物、物理和化学过程[-], 植被作为陆地生态系统的重要组成部分, 对气候变化比较敏感[-], 气候变化引起的植被变化已直接影响到区域陆—气交互过程的物质与能量平衡[]。因此, 开展区域及全球尺度的植被变化研究对认识气候变化与陆地生态系统的相互作用有重要意义[]。大尺度植被变化的连续观测对传统生态学研究技术提出了挑战, 随着全球变化生态学的发展[], 遥感技术逐渐成为大尺度生态问题连续观测的重要手段。卫星传感器通过探测陆表的光谱反射信息来识别植被的变化, 利用植被在可见光与近红外波段的反射特征构建的归一化植被指数(Normalized Difference Vegetation Index, NDVI), 常作为植被生长状况及覆盖度变化的最佳指示因子[], 已在不同尺度植被变化研究中广泛应用[-], 特别是时间序列的NDVI能较好的反映陆表植被活动的时间演化和空间变异特征[]。Myneni等[]利用1982—1991年的AVHRR-NDVI数据对全球植被活动进行了研究, 提出了北半球中高纬地区植被活动增强的观点；Piao等[]利用Global Inventory Modeling and Mapping Studies (GIMMS)第二代NDVI产品研究了1982—2006年亚欧大陆生长季的植被活动, 指出NDVI存在两个明显的变化阶段, 1997之前显著上升, 1997年之后NDVI上升趋势减缓；Park等[]同样利用该GIMMS产品研究了东亚地区的植被活动, 得出NDVI增加趋势在20世纪90年代中期发生明显转折；杜家强[]等利用AVHRR和MODIS两个传感器构建的连续NDVI数据对新疆近30年的植被活动展开研究, 认为NDVI的变化趋势在20世纪90年代末存在拐点。此外, 诸多学者关于不同地理单元[-]乃至全国范围[-]的植被活动研究认为, 研究区的植被活动总体增强。
然而, 以上研究存在几方面的问题, 一是大多研究基于上一代GIMMS NDVI[](1982—2006)、SPOT-VGT NDVI(1998—)或MODIS NDVI(2000—至今)等数据, 往往研究时间序列较短；尽管有学者延续不同传感器的NDVI数据开展植被变化研究[], 但不同传感器的性能差异又增加了NDVI数据的不确定性。二是传统的线性分析方法对时序较短的非平稳NDVI序列开展变化趋势研究, 会忽略非平稳数据的非线性波动特征, 致使研究结论和规律存在局限性, 而集合经验模态分解(Ensemble Empirical Mode Decomposition, EEMD)因其良好的自适应性可弥补该缺陷。三是前人研究多从自然单元和行政区划的角度探讨植被活动特征, 忽视了不同陆地生态系统内部的植被活动规律, 而我国类型多样的陆地生态系统往往镶嵌分布, 不同陆地生态系统的植被活动对气候变化的响应存在明显差异[]。为此, 本研究将利用1981—2012年的第三代GIMMS NDVI数据(GIMMS NDVI3g)和集合经验模态分解等非线性分析方法, 从生态系统的视角出发, 开展近30年中国陆地生态系统NDVI时空变化特征研究, 以期揭示气候变化背景下各生态系统的响应, 为我国制定应对气候变化的策略和生态文明建设提供科学依据。
美国国家航空航天局(National Aeronautics and Space Administration, NASA)基于NOAA-AVHRR系列数据开发出的GIMMS NDVI3g产品是一种长时序的全球植被指数数据集[](1981—2012年), 空间分辨率为8km, 时间分辨率是15d。该数据集消除了大气水汽、火山爆发、太阳高度角和传感器灵敏度的变化等因素的影响[], 精度得到有效提高, 良好的数据质量及连续的时间序列使其在大尺度生态过程模拟研究中得到广泛应用。本研究从NASA网站(http://ecocast.arc.nasa.gov/)获取了1982—2012年覆盖全球的GIMMS NDVI3g数据, 运用IDL软件编程对遥感影像进行读取、配准、格式转换以及研究区裁剪等预处理。
中国陆地生态系统类型宏观结构数据源于中国科学院资源环境科学数据中心(http://www.resdc.cn/)。基于遥感解译获取的1:10万土地利用/土地覆盖数据对各生态系统类型进行辨识、研究和重新分类, 形成中国陆地生态系统类型数据集(分别为农田生态系统、森林生态系统、草地生态系统、水体与湿地生态系统、荒漠生态系统、聚落生态系统和其他生态系统)。本研究采用2010年的数据集, 提取不同生态系统类型, 因聚落生态系统和其他生态系统类型所占比例较低, 同时考虑到聚落生态系统的碎小斑块很难在8km分辨率的GIMMS NDVI3g数据中反映出来, 其他生态系统主要为裸土地或裸岩石质地, 其NDVI指标无实际生态学意义, 故研究中未对这两类生态系统类型进行分析。
多变量厄尔尼诺-南方涛动指数(The Multivariate ENSO Index, MEI)来源于美国国家海洋与大气管理局地球系统研究实验室(http://www.cdc.noaa.gov/people/klaus.wolter/MEI/), 集合了海平面气压、海表温度、海面气温、地面纬向风和径向风、总云量等参数, 能够全面地监测和反映厄尔尼诺-南方涛动(El Niño-Southern Oscillation, ENSO)事件[]。太阳黑子数据采用比利时皇家天文台太阳活动数据分析中心(http://www.sidc.be/silso/datafiles)发布的月太阳黑子数序列, 该数据已采用13个月滑动窗口平滑处理。本研究截取1982—2012年的时间序列, 并用集合经验模态分解提取不同时间尺度的周期分量。
中国植被类型丰富多样, 同一生态系统往往兼有多个自然带的不同植被型, 且部分植被型存在明显的季相特征。为更全面描述年际植被活动特征, 本研究合成了年平均归一化植被指数(ANDVI)和年最大归一化植被指数(MNDVI)两种指标[], 分析1982—2012年中国陆地生态系统植被活动的空间变化。其中, ANDVI反映年内植被生长的平均状态, MNDVI反映当年最优水热组合下的植被最佳生长状态, 即最佳植被长势。
采用最小二乘法逐像元拟合1981—2012年间的GIMMS NDVI3g数据序列, 求取各像元近30年的变化斜率。一元线性回归斜率值的正、负分别表示NDVI升高或降低的趋势, 斜率的大小则反映了NDVI值上升或下降的程度, 回归趋势的显著性采用F检验判定。将NDVI变化趋势和不同显著性水平(P=0.05和P=0.01)叠加分析, 可得到极显著上升、显著上升、上升(但不显著)、下降(但不显著)、显著性下降和极显著性下降6种变化趋势组合。
集合经验模态分解(Ensemble Empirical Mode Decomposition, EEMD)[-]是近年来发展的一种适应于非线性、非平稳序列的信号分析方法, 该法在EMD(Empirical Mode Decomposition)的基础上, 引入辅助白噪声, 继承了EMD的自适性, 有效地解决了尺度混合问题；它将原始信号中不同时间尺度的周期特征分解开, 形成具有不同尺度周期的本征模态函数分量(Intrinsic Mode Function, IMF)和一个残差分量(Residual Standard Error, RSE)。其中, 所添加白噪声对原始信号的影响将遵循如下统计规律[]：
式中, σ为标准偏差, 即输入信号与相应IMF分量的偏差, ε为白噪声幅度, N为集合样本数。
经蒙特卡罗法验证, 当自由度趋于无穷, 添加白噪声的序列经过分解后, IMF的能量谱平均密度与平均周期具有如下关系：
式中, En为n个IMF的能量谱平均密度, Tn为相应平均周期。理论上, 白噪声IMF分量的平均周期与谱密度在二维空间中分布在斜率为-1的直线上, 但实际偏差使得连续概率分布下白噪声能量谱密度与平均周期有如下关系：
式中, En为平均能量, Tn为平均周期, df为自由度, 当k = 2.326和1.645分别可以得到99%和95%的置信区间。分布在不同显著水平曲线之上的IMF分量表示通过了相应水平的显著性检验, 即该IMF分量具有显著的物理意义。
EEMD分解15d合成的GIMMS NDVI3g序列及标准差序列时, 辅助白噪声与原始信号的标准差比率为0.2, 集合样本数为1000, 并借助白噪声的集合扰动对IMF分量进行显著性检验, 获取不同IMF的信度水平。其中, IMF分量表征从高频到低频不同时间尺度的NDVI周期波动特征, RSE反映了NDVI时间序列的残余变化趋势。
1982—2012年, 中国ANDVI整体呈上升趋势(0.0029/10a, P < 0.05)。从看出, 全国有59.31%的国土区域ANDVI呈上升趋势, 其中仅44.69%的区域超过5%显著性水平；全国有40.69%的国土区域ANDVI呈下降趋势, 其中有46.43%的区域超过5%显著性水平。从空间上看, 大兴安岭、嫩江流域上游、三江平原、长白山、长三角地区和珠三角地区的ANDVI显著下降, 降幅达到0.0080/10a以上(P < 0.01), 说明这些区域的年平均植被活动在减弱；此外, 西北地区的准噶尔盆地、塔里木盆地、柴达木盆地、甘肃西部和内蒙古西部广大地区的ANDVI降幅不大, 但下降趋势极为明显, 这与上述地区NDVI本底值较低有关。华北平原、黄土高原、秦巴山区、台湾山脉、东南丘陵局部、海南岛局部和横断山区的ANDVI极显著上升, 增幅达到0.0080/10a以上(P < 0.01), 局部地区增幅远超0.0240/10a, 说明这些地区的年平均植被活动在增强；此外, 藏北高原及天山局部地区的ANDVI虽然增幅不高, 但上升趋势极为明显。
近30年全国MNDVI整体也呈上升趋势, 且上升幅度比ANDVI更强(0.0076/10a, P < 0.01)。如所示, 占比65.24%的国土区域MNDVI呈增加趋势, 其中46.85%的区域超过5%显著性水平；占比34.76%的国土区域MNDVI呈下降趋势, 其中超过5%显著性水平的区域仅占39.27%。从空间上看, 极显著下降的区域主要分布在以荒漠生态系统为代表的内蒙古高原西部、准噶尔盆地和塔里木盆地, 说明这些地区年最佳植被长势在逐年减弱；而东北平原、华北平原、黄土高原、天山南北两侧绿洲区、大巴山脉、两广丘陵及云贵高原局部的MNDVI呈极显著上升趋势, 增幅达到0.0240/10a以上(P < 0.01), 说明这些地区年最佳植被长势在逐年增强。
由于MNDVI和ANDVI代表的植被活动意义不同, 导致两者的变化幅度差异较大, 尤其在东部季风区, MNDVI的变化幅度明显高于ANDVI。同时, 二者的变化趋势在一些区域存在差异, 甚至截然相反, 其中东南沿海的浙闵丘陵、台湾岛和海南岛地区的年平均植被活动在增强, 但年最佳植被长势却在变差, 人类活动持续增强导致森林植被破坏可能是重要原因[]；东北松嫩平原的年平均植被活动在减弱, 这与东北地区农田NDVI下降趋势一致, 但年最佳植被长势却在变好, 这种分异可能是暖干化趋势下松嫩平原农田生态系统的种植结构和作物熟制变化引起的。结合ANDVI和MNDVI来看, 全国植被活动显著增强的区域主要是分布在东部季风区的农田和森林生态系统, 而西北地区的荒漠生态系统和东北森林区的植被活动则显著减弱。
中国各生态系统面积占比由高到低依次是草地、森林、农田、荒漠和水体与湿地生态系统；ANDVI和MNDVI最高的是森林生态系统, 分别为0.5297和0.8437, 荒漠生态系统最低, 分别为0.0789和0.1543。线性趋势()表明, 除荒漠生态系统外, 各生态系统的植被活动主要表现出上升趋势。农田生态系统的植被活动显著增强, ANDVI的增幅为0.0061/10a(P < 0.01), MNDVI达到0.0162/10a(P < 0.05), 这与模型模拟的中国农业植被净初级生产力增加趋势一致[], 近30年中国农业种植结构调整、化肥用量增加、大气CO2浓度升高和气候变暖[]。森林生态系统的ANDVI、MNDVI均呈不同幅度的上升趋势, 但ANDVI表现出了较强的空间差异, 东北地区森林生态系统的ANDVI显著下降, 而南方亚热带森林生态系统的ANDVI则显著上升, 这与森林生态系统的碳源/汇空间分布格局较为一致[]；MNDVI极显著上升, 年最好植被长势显著增强可能与气候变化有关, 气温和大气CO2浓度的升高促进了森林植被年内最茂盛时期的长势。草地生态系统的ANDVI显著上升(0.0022/10a, P < 0.05), MNDVI极显著上升(0.0067/10a, P < 0.01), 说明我国草地地上生物量显著增加, 这与沈海花等[]的研究结果一致, 尽管草地生物量在近30年的变化趋势存在两个拐点, 但总体依然呈增加趋势。水体与湿地生态系统的NDVI值是二者的混合信息, NDVI变化主要源于湿地植被变化, 由于近30年我国湿地大面积萎缩[], 特别是水体面积的萎缩和湿地植物由水生向陆生的演替, 均会造成该生态系统植被盖度和生物量的增加, 导致NDVI升高。荒漠生态系统与前4类生态系统截然相反, ANDVI和MNDVI都处于下降趋势, 这与中国西北日益严重的荒漠化现状比较吻合[]。
表明, 农田、森林和草地生态系统的ANDVI均以上升为主, 农田生态系统上升像元比例达70.69%, 森林和草地生态系统分别为68.69%和63.49%。然而, 上升像元中, 农田、森林和草地生态系统分别仅有39.86%、27.31%和27.19%的像元通过了5%显著性检验, 可见, 各生态系统ANDVI近30年的上升态势并非稳定的线性趋势, 而是在上升过程中伴随着较大的波动。近30年, 荒漠生态系统有73.61%的像元呈下降趋势, 其中有54.47%的像元通过了5%显著性检验, 荒漠生态系统的植被存在明显的退化过程。水体与湿地生态系统的上升与下降像元比例几乎相当。各生态系统MNDVI近30年的趋势分析结果基本与ANDVI的结果一致, 即农田、森林和草地生态系统的植被呈波动上升态势, 且显著性更强, 荒漠生态系统呈现较明显的退化趋势, 区别在于水体与湿地生态系统的MNDVI有64.29%的像元表现为上升趋势。
对近30年中国各陆地生态系统的NDVI时间序列进行集合经验模态分解, 均可得到8个本征模态函数分量(IMF)及1个残差分量(RSE), 高频分量(IMF1—IMF3/4)主要为模型添加的白噪声和植被年内季节性生长波动信息, 低频分量(IMF4/5—IMF8)为年际、年代际的植被波动周期()。同时, 利用白噪声的集合扰动, 检验不同波动周期IMF分量的物理意义。
由于大气与海洋环流的脉动和人类活动间接作用于陆地生态系统, 导致植被活动可能存在周期性特征。各生态系统低频分量的平均周期()和显著性检验()表明, 农田生态系统的IMF分量除了IMF8外, 均通过了5%显著性检验, 可见农田生态系统的NDVI具有从年际到年代际的多重波动周期, 但随着波动周期的增长, 其显著性逐渐下降, 农业生产对市场规律的适应和长周期的气候波动都可能是影响因素。森林生态系统1.9a、5.3a和6.9a的波动周期通过了1%显著性检验, 但长周期的年代际波动显著性很低, 可能由于森林生态系统的群落环境相对比较优越, 漫长的群落演替过程受气候波动影响较小。由于草地生态系统包含了北方温带草地、青藏高原高寒草地及南方和东部次生草地, 区域气候条件差异大, 植被活动的气候波动痕迹不明显, 该生态系统仅6.8a和22.2a的波动周期通过了5%的显著性检验。水体与湿地生态系统的植被活动更易受到气候波动的影响, 主要表现出两个短周期的波动。荒漠生态系统的NDVI可能包含了大量非植被的背景信息, 仅9.6a周期通过了5%显著性检验, 可以认为荒漠生态系统的周期特性不具有明确的物理意义。上述分析表明, 以上5类生态系统主要存在1.9—7.6a的年际波动周期特征和9.6—11.7a、22.2a的年代际波动周期特征；其中年际波动周期特征普遍比较显著, 年代际波动周期特征显著性相对较差。由此推论, 各生态系统植被活动的周期性变化可能与海洋和大气环流的脉动有一定联系。
残差分量RSE反映NDVI序列的内在发展趋势, 可以看出, 近30年来, 虽然农田、森林、草地和水体与湿地生态系统的NDVI呈总体上升趋势, 荒漠生态系统的NDVI总体下降, 但内在包含了“上升—下降”两个阶段, 由之前的明显上升趋势转变为不同程度的下降趋势, 这与传统线性趋势显然不同[, ]。暖干化趋势下, 尤其是生长季降水减少引起的干旱胁迫是植被活动发生转折的重要原因[-]。农田、森林、草地、水体与湿地和荒漠生态系统NDVI趋势转折的节点有所差异, 主要集中在20世纪90年代和本世纪初, 植被活动转折期的差异在相关文献[, -]有所印证, 这可能是气候变化响应的差异和生态系统分类精度共同作用的结果。
对于单一生态系统的NDVI图像, 各像元NDVI之间的标准差越大, 则表明该生态系统植被长势的空间异质性越强, 植被活动的地域差异越大。从看出, 各生态系统空间异质性也存在2.1—7.1a的年际波动周期特征和10.7a、30.2a的年代际波动周期特征, 但各生态系统空间异质性的年代际周期波动普遍不具有物理意义, 主要表现为2.1—7.1a的年际周期节律, 这可能反映了短周期气候波动对生态系统植被活动的影响。
标准差序列的残差分量(RSE)表明(), 近30年中国各陆地生态系统植被长势的空间异质性总体增强, 植被活动的区域差异越来越大。荒漠生态系统NDVI的空间异质性持续增强, 直到本世纪后才逐渐减弱, 这与亚洲沙漠边缘植被活动变化的趋势比较一致[], 沙漠边缘有效水分的变化直接影响到植被活动, 进而反映到空间异质性的变化趋势。森林生态系统植被长势的空间异质性在20世纪90年代触底后持续增强, 农田、草地和水体与湿地生态系统的空间异质性一直处于增强趋势, 区域差异加大, Peng等[]指出, 近30年中国南、北方的植被长势存在明显差异, 暖干化趋势对中国北方植被长势的季节性抑制可能是中国陆地生态系统植被活动空间异质性增强的重要原因。
前述分解的周期与ENSO事件2—7a的变化周期、太阳活动有关的黑子周期(11—12a)、海尔周期(22a)、准35a的布吕克纳周期在时间上存在某种契合, 可能表明不同陆地生态系统植被活动与引起海洋和大气环流脉动的ENSO事件、太阳活动和海气相互作用等存在一定联系。Propastin[]和Nagai[]等指出陆表植被活动与ENSO事件显著相关, Lü等利用GIMMS数据研究发现中国各陆地生态系统有32.34%的区域的植被对ENSO事件响应敏感[]。本研究将MEI指数进行EEMD分解, 并与各主要生态系统NDVI具有相近波动周期的IMF分量进行相关分析, 发现ENSO事件存在2.1a、5.3a和7.1a的波动周期, 生态系统NDVI的波动周期与ENSO事件存在不同程度的相关性()。农事生产通过有目的的水分调节适应农时, 明显掩盖了自然气候系统对农田生态系统的周期性影响, 仅IMF6分量与ENSO事件的7.1a周期极显著相关；森林生态系统受人类活动干扰较少, IMF4、IMF5和IMF6分量分别与ENSO事件的2.1a、5.3a和7.1a周期极显著相关, 有着明显的气候短周期脉动痕迹；草地生态系统的周期性特征不明显, 震荡周期与ENSO事件无关；水体与湿地生态系统的IMF6分量与ENSO事件7.1a周期极显著相关, 湿地植被活动存在明显振荡周期；荒漠生态系统植被信息不足, 波动周期与ENSO事件无显著相关性。
此外, 太阳活动与全球气候变化关系密切[], 尽管太阳活动峰、谷值年的太阳总辐射变化极小, 但施加大气系统的影响可以通过热带地区海-气耦合作用得到强化, Hadley/Walker环流随着类La Niña事件的出现而减弱[], 这对东亚地区大气环流异常和中国南北地区反向降水变化影响显著[], 而降水、气温和辐射等气候变量则通过光合作用等生理过程影响植被的生长状况和群落演替。太阳活动11a周期与不同生态系统的NDVI变化存在不同程度的极显著负相关关系(), 这可能与太阳活动峰值年我国南方和东北夏季降水偏少有关, 而生长季尤其是夏季降水对植被生长极为重要, 显然植被活动11a左右的震荡周期含有明显的太阳活动信息。上述分析表明, 气候系统的脉动会引起陆地生态系统的NDVI波动, 海洋与大气环流的脉动与各生态系统植被活动的周期性节律有着一定关联, ENSO事件和太阳活动是推动植被活动周期性振荡的因素之一。
基于1982—2012年的GIMMS NDVI3g和中国陆地生态系统分布数据, 利用一元线性回归、集合经验模态分解和相关分析等手段, 研究了近30年中国各陆地生态系统NDVI的时空变化特征, 初步分析了NDVI变化与气候事件的关系, 研究取得如下几点结论：
(1) 我国近30年来NDVI主要以上升为主, ANDVI的上升幅度为0.0029/10a(P < 0.05), MNDVI的上升幅度为0.0076/10a(P < 0.01)。植被活动显著增强的区域主要是分布在东部季风区的农田和森林生态系统, 显著下降的区域主要是分布于西北的荒漠生态系统和东北森林区。
(2) 从线性特征来看, 尽管各ANDVI和MNDVI变化趋势的显著性有所差异, 但农田、森林、草地和水体与湿地生态系统的NDVI总体呈现非稳定的上升趋势, 上升过程中伴随着较大波动, 荒漠生态系统的NDVI呈下降趋势, 植被退化显著。
(3) 5类生态系统的植被活动主要存在1.9—7.6a的年际周期波动特征和10.7a、22.2a的年代际周期波动特征, 其中1.9—7.6a的年际周期特征普遍比较显著, 年代际周期的显著性相对较差。与线性趋势不同, 各生态系统植被活动的残差趋势包含“上升—下降”两个阶段, 并相继于20世纪90年代到本世纪初发生转折。
(4) 近30年来, 我国各生态系统植被活动的空间异质性存在2.1—7.1a的年际周期节律；随着NDVI的上升, 农田、森林、草地和水体与湿地生态系统的NDVI空间异质性持续增强, 生态系统内部植被活动的区域差异逐渐加大, 荒漠生态系统植被活动的空间异质性增强的趋势于本世纪初后逐渐减弱。
(5) 气候系统的脉动会引起陆地生态系统的NDVI波动, 海洋与大气环流的短周期脉动与各生态系统植被活动的周期性节律有一定关联, ENSO事件和太阳活动可能是推动植被活动周期性振荡的因素之一。