全球气候变化及其带来的影响已受到广泛关注, IPCC第五次评估报告指出, 2003—2012年平均温度比1850—1900年上升了0.78℃, 而降水在地区和季节上的差异也有所增加, 未来极端性天气气候事件的发生概率可能将进一步增加^[1]。中国的气候变化与全球的变化趋势是一致的, 平均增温速度高于全球或北半球同期, 而北方和青藏高原尤为明显^[2-3]；降水趋势变化存在明显的区域差异^[4-5], 其中东北地区^[6]和西部地区^[7-8]的降水出现下降趋势。

气候作为生态系统演变的重要控制因子, 它的任何变化都会对生态系统产生影响^[9-10], 进而对地带性植被的建群种分布产生影响^[11-12], 所谓地带性植被是指在分布上与气候带界线大体相符的大面积植被类型, 它综合反映了各自的大气候特点和一定的生态地理空间^[13-14]。有学者沿用Kira的温暖指数(WI)和寒冷指数(CI)划分地带性植被^[15], 有学者以WI为基础进行拓展^[16], 或在此基础上引入湿度指数(HI)^[17], 这些划分方法均是建立在气候是地带性植被类型分布的主要影响因子基础上的, 因此对东北地区的划分基本一致, 即东北地区的自然植被属欧亚森林、草原植物亚区和中国-日本森林植物亚区, 主要地带性植被为寒温带针叶林、温带针叶混交林和暖温带落叶阔叶林。东北区地处全球气候变化的敏感区域, 地域辽阔, 地势复杂, 气候变率大^[18-19], 因此许多学者就气候变化对该区域植被分布产生的影响做了一些研究, 如冷文芳等^[20]在模拟100a后气温和降水变化的基础上, 使用logistic回归模型预测东北森林主要建群树种的空间分布变化, 结果表明多数树种的覆盖率可能降低, 只有红松、长白落叶松和蒙古栎有可能增加。程肖侠^[21]也在预测模拟未来100a气候情景下东北地区森林类型的变化, 认为气候变暖不利于东北主要森林类型的生长, 主要的针叶树种比例可能会减少, 阔叶树种比例会增加。以上研究均是建立在气候预测基础上的, 但目前研究表明东北地区的气候变化已成不争的事实^[22], 这种变化对该区域森林主要建群树种的分布区域是否已经产生了影响, 影响程度如何还少见报道, 本研究正是针对这一问题, 在总结前人研究成果的基础上, 收集近十几年新增加的科研资料, 根据东北地带性建群树种和常见种温暖指数范围, 从生态气候学的角度探讨气候变化对东北三省森林植被地理分布的影响。该研究不仅有利于区域森林植被的恢复与重建, 也可为该区域自然资源的合理利用、土地退化防治、生物多样性的保护等区域可持续发展战略和规划的制定提供理论参考。

1 研究区概况

东北三省包括黑龙江省、吉林省和辽宁省, 位于118.83°—135.09°E, 53.56°—38.72°N之间, 南邻黄海和渤海, 东面和北面有鸭绿江、图们江、乌苏里江和黑龙江环绕, 西面与内蒙古接壤, 总土地面积78.7万km², 境内有大、小兴安岭、长白山系和东北平原(包括松辽平原、辽河平原和三江平原), 海拔高度在2—2667m之间。该地区属温带季风气候, 其中北部大兴安地区属寒温带气候, 其它地区属温带气候, 年平均气温-4.15—10.99℃, 年降水量366.90—1080.47mm, 日平均日照时数5.22—8.10h。

2 数据来源与处理 2.1 地理信息数据

数字高程模型(DEM)为SRTM地形产品V4.1版本数据, 空间分辨率为90m, 数据来源于中国科学院计算机网络信息中心国际科学数据镜像网站(http://www.gscloud.cn)。

研究所需的县级和省级行政区划数据来自于中国气象局下发的1:25万基础地理信息, 对数据进行拓补检查, 去除省界和县界的间隙。地面观测站位置数据为中国气象局下发的矢量数据。

2.2 土地覆盖类型数据

土地类型覆盖数据(MCD12Q1) 来源于美国NASA LPDAAC(The Land Processes Distributed Active Archive Center)EOS数据中心, 该数据的空间分辨率为500m, 是利用当年Terra和Aqua星观测所得的数据描述土地覆盖类型, 信息提取主要技术是监督决策树分类, 产品包含5种不同土地覆盖分类方案的分类结果和对应的质量信息数据集, 本研究选用“IGBP的全球植被分类方案”的分类结果, 该分类方案根据国际地圈生物圈计划(IGBP), 本文所用数据的时间为2001年和2013年, 空间范围为东北地区, 利用MODIS Reprojection Tool software(MRT)V4.0, 进行影像的拼接、重采样(统一空间分辨率为500m)和投影转换(由Sinusoidal投影转换为经纬度投影)等预处理。

2.3 气象数据

气象数据来源于黑龙江省、吉林省和辽宁省气象局的1961—2013年的196个气象台站逐日整编资料, 包括日平均气温和日降水量, 站点的地理位置和海拔高度如图 1。

2.3.1 气象数据预处理

气候变化研究必须建立在可靠的资料基础上, 而由于气象台站迁移、观测数据误差和缺失、城市热岛效应等原因, 会产生数据序列的不均一性, 因此在进行分析之前需要对数据进行预处理, 最大限度地降低因数据误差而影响分析结果的准确性。

(1) 质量控制

由于观测仪器、资料处理和传输过程可能引起误差, 资料中很难避免存在一定的错误, 最大限度地减少这些错误的可能影响是非常必要的^[23], 因此本研究采用翟盘茂^[24]的方法对气象数据进行质量控制, 删除错误数据。另外Karl^[25]也指出, 把与时间有关的资料缺测假设为零时, 会在气候变化趋势研究中引起误差, 本研究中监测缺测数据的方法是当一年内缺测超过20d(约大于一年内实际观测日数的5%)时, 该年的统计量被设为缺测, 当45a序列中缺测总量超过5a时, 则该测站被排除使用。

(2) 正态分布检验

多数气候诊断方法和预测模型是在气候变量呈正态分布假定前提下进行的, 所以对气候变量是否呈正态分布形态的检验是十分必要的。本文利用SPSS软件对各气候区逐日气温和相对湿度数据进行正态检验, 显著性水平sig.有98.7%均≥0.05, 说明各气候区逐日气温数据绝大多数近似正态分布, 分析时可不对其进行标准化处理。与气温数据不同, 日降水量的变化并不具有逐步的、连续的特征, 因而判断降水数据是否服从正态分布是以月为自变量进行检验, 检验结果绝大多数服从正态分布, 因而分析时不作标准化处理。

(3) 热岛效应分析

IPCC第3次评估报告认为, 从全球平均来看, 城市化对地表平均气温记录的影响微乎其微, 其实际作用要比观测到的增温值小一个数量级^[26], 而Wang^[27]和赵宗慈^[28]却认为, 热岛效应对长期气温序列研究的影响不容忽视, 于是本研究根据2015年统计年鉴中人口数据, 删除人口大于等于100万人和人口大于等于50万人的城市测站的气温观测数据, 计算未删除、删除人口大于等于100万人城市和删除人口大于等于50万人的城市的年平均气温变化速率, 结果分别为0.0149、0.0126、0.0145℃/a, 它们之间的差别比增温值(21世纪以来比19世纪60年代高0.16℃)小一个数量级, 因此本研究不考虑城市热岛对气温序列变化的影响。

(4) 数据插补

本文对缺测测站的气象数据采用Newton法进行插补。

2.3.2 气象数据突变检测

通过对东北地区1961—2013年气温和降水突变点的检验, 为探寻气候变化对研究区主要地带性植被类型分布可能带来的影响提供时间节点。

(1) Mann-Kendall检验

Mann-Kendall(简称MK)是气象学/气候学中经常用来进行突变检验的一种方法, 具体方法如下：

设原始时间序列为y₁, y₂, …, y_n, m_i表示第i个样本y_i大于y_j(1≤j≤i)的累积数, 定义统计量：

在原序列随机独立等假设下, d_k的均值和方差分别为：

将上面公式的d_k标准化, 得：

UF_k组成一条UF曲线, 通过信度检验可得出其是否有明显的变化趋势。把此方法引用到反序列中, 计算得到另一条曲线UB, 则两条曲线在置信区间内的交点确定为突变点。

(2) 滑动T检验

滑动t检验法, 它是用来检验两随机样本平均值的显著性差异^[29]。为此, 把一个长度为n的连续随机变量x分成两个样本子集x₁和x₂, 让μ_i、S_i2和n_i分别代表x_i的平均值、样本方差和样本长度(i=1, 2)。两样本子集始终间隔一个样本, 这样检验的就是某一年后n₂年和前n₁年均值的显著性差异。

原假设H₀:μ₁-μ₂=0。定义t统计量为：

这里, S_P是联合样本方差,

为σ₂的无偏估计(E[S_P²]=σ²), 显然t₀-t(n₁+n₂-2) 分布, 给定信度α, 得到临界值t_α, 计算t₀后在H₀下比较t₀与t_α, 当t₀≥t_α时, 否定原假设H₀, 即说明其存在显著性差异。当|t₀| < t_α时, 则接受原假设H₀。

用MTT来检验某一个时间序列, 可以得到相应的t₀统计量序列, 峰值和谷值所对应的年份称之为突变年份、突变点或者过渡年份。

(3) 累积距平检验

累积距平也是一种常用的、由曲线判断变化趋势的方法, 对于序列x, 其某一时刻t的累积距平表示为：

其中

将n个时刻的累积距平值全部算出, 即可绘出累积距平曲线进行趋势分析。

3 结果与分析 3.1 气象数据突变检验 3.1.1 研究区划分

东北三省跨越寒温带、中温带和暖温带, 且地形复杂, 山地、丘陵和平原相间, 海拔高度跨越2665m, 气候变化趋势和突变点可能会存在差异, 因此本研究参照《中国生态区划研究》^[30]提供的中国生态地理分区数据(来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)提供的矢量数据), 分别检验东北三省的7个生态地理分区(图 1)的年平均气温和年降水量的突变点。

该生态地理分区数据是基于温度指标把全国划分为寒温带、中温带、暖温带、北亚热带等11个区域, 根据干湿状况指标把全国划分为潮湿、湿润、半湿润、半干旱等6个区域, 经上述2种指标组合后把全国划分成46个生态地理分区, 东北三省占其中的7个生态地理分区(图 1)。

3.1.2 突变点检测

从年平均气温MK检验结果可见(图 2), 东北三省7个生态地理分区的年平均气温均存在明显的上升趋势, 其中生态地理区Ⅰ的UF和UB在1982年处交叉, 表明从1982年开始该区的气温可能发生突变, 同时滑动T检验结果和通过计算气温序列的累积距均显示1982年为该区气温突变年份；用同样的方法检测生态地理区Ⅱ—Ⅶ的气温突变年为1988年, 该结论与^[31]一致。

年降水量MK检验结果与年平均气温不同, 生态地理区Ⅱ、Ⅲ、Ⅴ和Ⅵ的MK降水量统计值UF在[-1.96, 1.96]波动, 说明这几个生态地理区在研究时间段内年降水量无显著变化；生态地理区Ⅳ从1976年开始显著下降, 1983年开始显著上升；生态地理区Ⅶ从1963年开始下降, 在2003年出现显著下降趋势。结合MK和滑动T检验结果, 东北地区的7个生态地理区年降水量没有明显的突变点, 因此本文只研究气温突变对东北地区地带性植被分布适宜区的影响。

3.2 气候变化对东北地带性植被分布的影响 3.2.1 植被分布热量指标WI

按照徐文铎^[17]的方法计算水热指标——温暖指数(WI)：

式中, t为5℃以上的月平均气温。

再根据树种的地理分布资料, 按半峰宽(PWH)计算法^[32]确定每个树种的热量分布范围：

式中, X为热量指标的平均值, S为标准差。这种确定方法能使树种约有78%分布在最适热量范围之内^[17]。表 1为6个东北地带性植被建群种和常见种植被类型的WI范围。

3.2.2 WI插值

把第3.1.2节气候突变点检测结果作为分界点, 分别计算突变点前、后各个气象观测站点的WI平均值, 即利用1961—1981年生态地理区Ⅰ和1961—1987年生态地理区Ⅱ—Ⅶ的数据计算突变点前每个站点WI的年平均值(WI₁)；突变点后的WI的平均值(WI₂)所用数据的时间范围为1982—2013年生态地理区Ⅰ和1988—2013年生态地理区Ⅱ—Ⅶ的数据, 计算方法同上。

由于东北三省的气象观测站点少且分布不均, 尤其在地形复杂的山区气象观测站分布尤其稀疏, 这些有限的气象观测站无法代表相距较远地区的气候状况, 基于此数据计算的WI的代表性也会较差, 为解决这一问题, 本研究对单点的WI数据进行插值, 由于WI是基于气温计算而来, 而气温与经度、纬度和海拔高度又具有较好的相关性^[33-34], 因此本研究以经度、纬度和海拔高度为自变量, 分别以WI₁和WI₂为因变量拟合多元线性回归模型, 拟合结果见表 2。

利用回归方程推算出各个观测站的WI₁和WI₂的拟合值, 用拟合值减去实际计算值得到的差值分别记为WI_1a和WI_2a, 再利用IDW插值方法把这些点的差值转换为面, 用这个面数据作为订正系数来订正利用回归方程拟合出的WI₁和WI₂面数据, 得到订正后的WI₁和WI₂的面数据, 再从中提取各个观测点订正后的点数据, 并减去实际计算值得到的差值分别记为WI_1b和WI_2b, 分别统计对比WI_1a与WI_1b和WI_2a与WI_2b的标准误差、标准差和数值范围, 结果表明, WI_1b标准误差、标准差和数值范围分别比WI_1a小0.13、1.77、6.82cm, WI_2b标准误差、标准差和数值范围分别比WI_2a小0.12、1.57、7.83cm, 说明经过订正后的WI值与实际计算值更为接近。

3.2.3 植被类型分布区域变化

分析图 3可见, 气温突变点前适宜在东北三省分布面积最大的是温带针阔混交林, 为55.89万km², 主要分布于松嫩平原、辽河平原、三江平原、小兴安岭和长白山海拔小于800m左右的区域, 在气温突变点后该植被类型的适宜分布面积为47.22万km², 总体适宜面积减少了8.67万km², 其中主要收缩区域为辽河平原, 面积为11.71km², 扩张区域为大兴安岭东部和南部以及小兴安岭海拔在400—500m和长白山海拔在800—1100m的区域, 面积为3.04万km²。

气温突变点后适宜分布面积增加的植被类型为暖温带落叶阔叶林和温带草原, 其中暖温带落阔叶林在气温突变点前的适宜分布面积为55.39万km², 主要分布在除大兴安岭、小兴安岭和长白山海拔低于600m以下的其它区域, 在气温突变点后适宜该植被类型分布的面积为60.90万km², 总体适宜面积增加了5.51万km², 成为东北地区适宜分布最广泛的植被类型, 其中主要收缩的区域在东北南部沿海区域, 面积为1.62万km², 扩张的区域在松嫩平原北部、小兴安岭海拔在300—400m和长白山海拔在600—800m之间的区域, 面积为7.13万km²。

气温突变点后适宜分布面积减小最多的植被类型为寒温带针叶林, 也是气温突变点前后适宜面积变化最大的植被类型。在气温突变点前适宜分布的面积为46.54万km², 主要分布在大兴安岭、小兴安岭、长白山、松嫩平原北部和三江平原大部, 气温突变点后面积为30.91万km², 总体适宜面积减少了15.63万km², 其中主要收缩的区域在松嫩平原北部和三江平原大部, 面积为15.83万km², 扩张的区域为大兴安岭海拔在800m以上的区域, 面积仅有0.20万km²。

气温突变点后适宜分布面积增加最大的植被类型为温带草原, 在气温突变点前适宜分布的面积为50.47万km², 主要分布在松嫩平原大部、辽河平原中部和北部、三江平原、长白山海拔低于600m以下的区域, 在气温突变点后适宜该植被类型分布的面积为53.46万km², 总体适宜面积增加了2.99万km², 其中收缩的区域在辽河平原西部和北部、长白山南部, 面积为4.14万km², 扩张的区域在小兴安岭海拔在300—400m和长白山海拔在800—1100m的区域, 面积为7.13万km²。

高山冻原在东北地区适宜分布面积最小, 主要分布在大兴安岭海拔在900m以上和长白山海拔在1800m以上的区域, 气温突变点前的适宜分布面积为104.79km², 突变点后减小了27.91km²；亚高山矮曲林的分布面积也较小, 主要分布区域为大兴安岭海拔在700—900m, 气温突变点前的适宜分布面积为3.03万km², 突变点后面积减小为2.79万km²。

3.2.4 植被类型分布地理中心变化

气温突变点前后各植被类型分布的地理中心均发生了不同程度的移动(图 4), 从移动方向上看, 北部的高山冻原和寒温带针叶林向西北方向移动, 其它植被类型分布的地理中心均向东北方向移动。从移动距离上看, 最大的是南部地区的亚高山矮曲林, 移动了135.44km；其次是温带针阔叶混交林, 移动了71.80km；移动距离最小的是南部的高山冻原, 只移动了1.39km；其它植被类型的移动距离在22.20—71.47km之间。

3.2.5 推算分布区与实际分布区对比

由于本研究缺乏实地调查资料, 因此借助2001年和2013年MCD12Q1数据的分类结果, 与3.2.2中各植被类型的分布区域进行对比, 在MCD12Q1数据的分类结果中落叶针叶林、混交林、落叶阔叶林和草地分别对应3.2.2中植被类型的寒温带针叶林、温带针阔叶混交林、暖温带落叶阔叶林和温带草原。研究中首先利用MCD12Q1数据的分类结果提取2013年各土地覆盖类型与2001年的空间差异, 增加部分的数据分别设为A1、A2、A3和A4, 减少部分的数据分别设为B1、B2、B3和B4；再设3.3.2研究结果中各植被类型在突变点后增加部分的数据为a1、a2、a3和a4, 减少部分的数据分别设为b1、b2、b3和b4, 分别求取A1和a1、…、A4和a4、B1和b1、…、B4和b4的共同区域并计算面积。

结果表明(图 5), 寒温带针叶林无共同减少区域；共同增加的区域分布在大兴安岭海拔在800m以上的地区, 增加面积为0.01万km², 占A1面积的4.82%, 另外82.57%的A1面积分布在气温突变点前后未变的适宜区域。温带针阔叶混交林无共同减少区域；共同增加的区域在大兴安岭东部、小兴安岭海拔在400—500m区域和长白山海拔在800—1100m的区域, 增加面积为1.27万km², 占A2面积的21.46%, 另外61.01%的A2面积分布在气温突变点前后未变的适宜区域；暖温带落叶阔叶林无共同减少区域；共同增加的区域在松嫩平原北部、小兴安岭海拔在300—400m和长白山海拔在600—800m之间的地域, 面积为0.78万km², 占A3面积的20.23%, 另外59.94%的A3面积分布在气温突变点前后未变的适宜区域；温带草原共同增加的区域在小兴安岭海拔在300—400m和长白山海拔在800—1100m的区域, 面积为0.10万km², 占A4面积的8.03%, 另外80.15%的A4面积分布在气温突变点前后未变的适宜区域；共同减小的区域在辽河平原西部和北部以及长白山南部, 占B4面积的18.72%, 另外70.24%的B4面积分布在气温突变点后的不适宜区域。可见气候变化对各植被类型的分布已经产生一定影响。

4 结论与讨论 4.1 结论

以东北三省为研究区, 以1961—2013年的196个气象台站逐日整编资料为基础数据, 以温暖指数的适宜范围为标准, 研究气候变化对东北地带性植被适宜分布区域的影响, 并结合2000年和2013年土地类型覆盖数据(MCD12Q1), 对推算适宜分布区域进行验证, 结果表明：

(1) 近50年来大兴安岭气温突变点为1982年, 其它地区为1988年；东北三省年降水量没有明显的突变点。

(2) 通过建立WI与经度、纬度和海拔高度的线性拟合方程, 对缺少气温观测地区进行数据插补, 再利用差值对插补数据进行订正, 订正后的WI值与实际计算值更为接近。

(3) 气温突变点后适宜东北三省生长的主要植被类型种类无变化, 但各植被类型的适宜分布区域均发生变化, 气温突变点前适宜在东北三省分布面积最大的是温带针阔混交林, 气温突变点后适宜分布面积有所下降, 其中收缩的区域为辽河平原, 扩大的区域为大兴安岭东部和南部以及小兴安岭海拔在400—500m和长白山海拔在800—1100m的区域；气温突变点后适宜在东北三省分布面积最大的是暖温带落阔叶林, 其中的收缩的区域在南部沿海区域, 扩大的区域在松嫩平原北部、小兴安岭海拔在300—400m和长白山海拔在600—800m之间的区域；气温突变点后面积减小最多的是寒温带针叶林, 也是气温突变点前后适宜面积变化最大的植被类型；气温突变点后面积增加最大的是温带草原。

(4) 与主要植被类型实际分布变化相对比, 在寒温带针叶林、温带针阔混交林、暖温带落叶阔叶林和温带草原中, 2013年比2001a增加的区域分别有4.82%、21.46%、20.23%和8.03%的面积落在理论推算的增加区域内, 由此可见从理论推算和实际分布上均证明了气候变化对研究区主要植被类型分布可能已经产生了影响。

4.2 讨论

(1) 本研究在检验近50年来气温突变点的基础上, 推算气候变化对东北三省主要植被类型适宜分布区的影响, 并选取典型植被类型利用现有的遥感分类数据进行对比, 研究结果与已有的研究结论有很多一致之处, 如程肖侠等^[21]认为在维持当前气候条件不变的情况下, 东北森林树种组成基本维持平衡, 但气候增暖使主要针叶树种比例下降, 阔叶树种比例增加, 与本研究结论中气温突变点前后东部地区适宜在东北生长的主要植被类型种类无变化, 但寒温带针叶林和温带针阔混交林在气温突变点后适宜分布区面积下降, 暖温带落叶阔叶林适宜分布区面积增加的结论基本一致；在本研究中气温突变点后温带针阔叶混交林的适宜分布区域已经侵入大兴安岭东部, 这一结论验证了程肖侠等^[35]在研究大兴安岭森林演替中指出的气候变化将可能使大兴安岭会以温带针阔叶混交林为主的预测；冷文芳^[20]在HADCM2SUL和CGCM1模型预测假设100a后全球平均气温和降水量增加的基础上, 认为寒温带针叶林主要建群种(兴安落叶松和云冷衫)的覆盖率有明显下降, 这一结论与本研究结果一致, 但对于本研究中的温带针阔混交林的主要树种红松, 本研究的结论更倾向于CGCM1方案的预测结论；晏寒冰^[36]等认为气候变化使兴安落叶松和白桦的分布向西北漂移, 这与本研究中寒温带针叶林适宜分布区域的地理中心向西北移动的结论一致。

(2) 本研究得出的气候变化对研究区主要植被类型分布可能已经产生了影响的结论, 是在理论推算和遥感分类数据相结合的基础上得出的, 这当中主要误差来源有：一是理论推算方面, 本研究推算气候变化对植被类型分布的影响是建立在气候变化突变点基础上的, 气候变化突变检验的结果是近50年来只有温度因子发生了突变, 因此认为温度因子在其中起了主要作用, 尽管在本研究中已检测到降水无显著变化, 但降水、自然灾害、人为活动等诸多因素对植被类型分布也具有一定的影响, 因此本研究得出的结论具有一定的片面性, 但方精云^[37]和邹春静等^[38]认为, 单因素标准对于确定大区域植被带划分具有较好的指示意义, 而且刘向培^[39]的研究结果表明, 在东北地区温度对地表植被动态变化的影响大于降水, 因此本研究只根据温度因素推算气候变化对植被类型分布的结果仍具有一定的参考意义；二是遥感分类数据的精度问题, 本研究直接采用了美国NASA LPDAACEOS数据中心提供的土地类型覆盖数据(MCD12Q1), 由于缺乏历史地面调查数据, 因此无法对其分类精度进行检验, 分类误差也会影响本研究结果准确程度, 虽然在讨论① 中已有相关研究成果支持了本研究的结论, 但本研究给出的每种植被类型分布区域变化的具体数据还有待进一步验证。

(3) 虽然物种分布与气候因素有密切关系, 但物种的迁移对于气候变化的影响往往有一定的滞后性^[40], 也会受到自身的迁移能力、人为干预和诸多外界条件的影响, 而本研究在推算植被类型分布适宜区域变化时并未考虑上述因素的影响, 这也是把推算结果与植被类型实际分布变化作对比时, 共同增加(减少)面积占实际增加(减少)面积比例较低的原因之一。