随着全球气候变暖和人类活动的加剧, 陆地植被生态系统受到不同程度的干扰。植被是物质循环和能量流动的纽带, 对环境变化反应十分敏感^[1]。大兴安岭林草交错区是我国重要的生态脆弱区, 面临草地退化、土地沙化显著、水土流失严重、调蓄功能下降等生态问题的威胁^[2—3]。为科学把握大兴安岭林草交错区植被的时空变化趋势, 有效开展生态脆弱区生态保护和资源开发, 研究其植被的影响因子相对重要性对生态脆弱区稳定可持续发展提供理论基础。

近年来, 我国生态脆弱区植被动态变化及影响因子进行了大量研究。研究表明生态脆弱区生态系统变化的主导因子时空差异较大^[4]。降水是影响青藏高原脆弱区植被NDVI年际变化的主要气候因子, 气温升高则为植被提供更适宜的湿润环境^[5]；同时降水促进了北方农牧交错区植被生长, 而气温升高则导致植被退化^[6]；西北荒漠绿洲交接生态脆弱区蒸散成为影响植被NDVI的主要气候影响因子^[7]；东北林草交错区寒冷干燥, 植被生长即受气温限制, 又受水分限制存在气候记忆现象导致该地区植被动态变化主要受土壤水分影响^[8]。调查研究结果认为该地区植被动态变化与温度显著相关, 而与降水关系不显著^[9], 1982—2009年长时间序列研究进一步表明春季温度对植被动态变化影响高于降水^[10]。然而研究东北林草交错区的典型草原区和森林草原生态子区2001—2010年植被NDVI结果表明降水是该地区的主要影响因子^[11—12]

此外, 随着社会经济的发展, 人类活动对植被的影响强度不断增加逐渐受到重视, 并且人类活动对植被的影响有提升和抑制两面性^[12—14]。研究气候和人类活动因子对植被影响的学者发现气候因子是植被的主导影响因子^{[8, 10]}, 另有研究认为人类活动是植被的主导影响因子^[15—16]。这一矛盾结论表明不同区域不同植被类型对气候和人类活动的响应存在差异, 并且不同时空尺度对结果也存在影响。针对这一问题, 已有研究主要通过利用气候因子对植被进行多元回归, 结合残差得到气候和人类活动的贡献率^[16—17], 缺乏地面数据实证分析和气候因子间多重共线影响的考虑。

鉴于此, 本文利用遥感数据和地面实测数据, 采用偏最小二乘分析(Partial Least Squares Regression, PLSR)和变量投影重要性指标(variable importance in projection, VIP)等方法在区域和点尺度上讨论了1982—2015年林草交错区植被动态变化及其影响因子, 量化气候因子和人类活动因子对植被动态变化的相对重要性, 为保护和开发生态脆弱区, 实现人与自然和谐发展提供科学依据。

1 材料与方法 1.1 研究区概况

大兴安岭林草交错区地处中国东北部内蒙古自治区, 大兴安岭西麓山地向呼伦贝尔草原过度地带, 包括蒙古根河市、额尔古纳市、牙克石市、陈巴尔虎旗、呼伦贝尔市、鄂温克自治旗、新巴尔虎左旗7个市/旗行政区, 总面积约11.66×10⁴ km², 海拔在700—1700 m范围内(图 1)。大兴安岭林草交错区属于生态功能区划的东部季风生态区向西部干旱生态区交界区, 夏季高温多雨, 冬季寒冷干旱, 气温日较差和年较差很大, 年降水量约为400 mm, 风向和降水均有明显的季节性。大兴安岭林草交错区地势东高西低, 东部为大兴安岭西麓呼伦贝尔林区, 海拔在700—1700 m；而西部为呼伦贝尔大草原, 海拔在500—1000 m。植被类型由西南向东北呈现草地-混交林-落叶针叶林过度的变化趋势, 是中国北方重要的生态屏障和国家重要的生态脆弱区。

1.2 数据来源 1.2.1 遥感数据

气象数据来源于英国East Anglia大学的Climatic Research Unit (CRU, https://crudata.uea.ac.uk/cru/data/)。CRU通过整合已有的气候数据库, 重建了一套空间完整、时间连续的月平均地表气候要素数据集, 空间分辨率为0.5°×0.5°, 网格覆盖全球陆地。为表征研究区气候变化特征, 本研究选取1982—2015年平均温度(TMP)、最高温度(TMX)、最低温度(TMN)、降水量(PRE)、水汽压(VAP)、潜在蒸散量(PET)、云量(CLD)等气候要素(表 1)。

植被指数(NDVI)数据集来源于NASA戈达德航天中心GIMMS NDVI3g, 合成时段为15 d, 空间分辨率为1/12°, 时间跨度是1982—2015年。此外, 2015—2018年的NDVI数据为MOD13C1的合成16 d的最大值, 空间分辨率为0.05°的植被指数数据。由于1982 —2015年的CRU TS v4.01气象数据与GIMMS NDVI3g数据分辨率不同, 采用三次卷积法(Cubic Convolution)重采样为1/12°的数据从而进行相应分析。

1.2.2 地面实证数据

(1) 额尔古纳、陈巴尔虎旗2个旗县的气象台站1957—2018年的逐年、逐月平均气温和降水量数据。(2)额尔古纳的上库力农场和陈巴尔虎旗的陶海牧场1955—2019年开垦面积数据(小麦、油菜和其他作物等)。(3)额尔古纳的上库力农场、苏沁农牧场和陈巴尔虎旗的陶海牧场1955—2019年的家畜密度数据；以及当地40户牧民1980—2019年的平均家畜密度调查数据。(4)额尔古纳、陈巴尔虎旗植被NDVI通过将1982—2015年的GIMMS NDVI和MOD13C1 NDVI数据采用三次卷积法(Cubic Convolution)重采样为1 km×1 km, 同时采用像元级融合方法(主成分变化PCT)将数据进行融合, 然后利用5 × 5的窗口提取上库力农场和陶海牧场NDVI。

1.3 研究方法 1.3.1 趋势分析

在研究植被NDVI或气候因子的变化特征时, 常采用线性回归模型的斜率表征该参数随时间的变化趋势(公式1)。本研究采用最小二乘线性拟合计算变化率：

(1)

式中, 变量j为时间序号(j=1, 2, ……, j), k_j为第j年的变量值。slope表示变量的变化趋势, slope > 0表明变量在研究时间段内是增加的, slope < 0表明变量在研究时段内是降低的。

1.3.2 偏最小二乘回归分析

偏最小二乘回归(Partial Least Squares, PLS)分析是一种通过主成分而分析提取多因变量重要特征信息, 剔除多重相关信息构建新变量对因变量进行归回分析的方法, 该方法具有主成分分析、典型相关分析和线性回归分析等方法的特点, 可以有效解决变量间多重相关性的问题。其解决途径是对变量的信息进行重组, 提取信息时既考虑因变量和自变量的线性关系, 又考虑因变量之间的综合性解释性, 从而消除多重共线问题对模型的稳定性的影响^[18]。

1.3.3 残差分析

残差分析(residual analysis)就是通过回归模型中残差提供的信息, 分析模型的可靠性, 并识别模型预测的干扰程度(公式2)。

(2)

式中, ε_i表示植被NDVI第i年的实测值与预测值之差, NDVI_i表示第i年的实测NDVI值, NDVI′_I表示预测NDVI。

Hurst指数(公式3和4)对于植被退化区域、人类活动影响等具有较强的指示作用^[19—20]。当残差Hurst指数为0.5时, 表明系统无外界干扰；而当0.5 < 残差Hurst指数 < 1时, 表明外界干扰也随之增强。

(3)

(4)

式中, n表示样本个数, Y(n)表示残差累计离差和, S(n)表示样本标准差, C₁为常数。

1.3.4 四参数logistic模型

四参数logistic模型可以有效反映家畜密度和开垦面积增长过程的上升点、拐点、下降点能定量的划分增长的起步、发展、巩固和停滞阶段, 充分反映增长的生命周期。公式为：

(5)

式中, a₁为最低的家畜密度极限值；a₂为最高的家畜密度极限值；t为时间；t₀是曲线下凹和上凸的拐点, 是N/2畜的数据收集年份；p是N/2的斜率。当t的取值为(0, t₁), 曲线单调递增且增速较低, 为初级阶段；当t的取值为(t₁ t₀), 曲线单调递增, 增速较快, 将其划为发展阶段；当t的取值为(t₀, t₂), 曲线增速降低, 将此阶段划为巩固阶段；最后, 当t的取值从t₂过渡到曲线终点。

1.3.5 变量重要性排序

采用变量投影重要性指标(variable importance in projection, VIP, 公式6), 定量分析不同因素对NDVI变化的影响程度, 识别驱动植被NDVI变化的主导因素。

(6)

式中, w_jf是j变量和f分量的权重值, SSY_f是第f个分量解释方差的平方和, J是解释变量的个数, SSY_total是因变量解释方差的平方总和, F是分量的总数。w_jf²给出了每个f组分中第j个变量的重要性, VIP_j表征了j变量的贡献度。

2 结果与分析 2.1 大兴安岭林草交错区植被指数与气候因子时空变化分析 2.1.1 植被指数与气候因子空间分布

1982—2015年东北大兴安岭林草交错区年平均NDVI为0.39, 呈现由西南向东北逐渐增加的趋势(0.15—0.59), 主要原因是由西南向东北为草地向林地逐渐变化引起。此外, 1982—2015年气候因子的年均CLD、PET、PRE、TMN、TMP、TMX和VAP分别为49.93%、1.79 mm、406.66 mm、-9.43 ℃、-2.33 ℃、4.76 ℃和5.51 hPa。年均气候因子的空间分布图(图 2)可知, 年均降水由东到西呈现逐渐增加的趋势, 而平均温度、最低温度、最高温度、水汽压和潜在蒸发呈现逐渐降低的趋势。

综合气候因子和NDVI的空间分布可以看出, 植被空间分布与降水和指示干旱的指标空间分布一致, 大体以年降水量400 mm左右为界, 东部植被类型为林地为主, 西部以草地为主。此外, 气温和潜在蒸散量较高年降水量在200—400 mm的西部地区, 干旱程度较强, 适于草本植物。然而, 年降水量在400—500 mm的东北大兴安岭林草交错区的东部气温和潜在蒸散量较低, 有利于乔木和木本植物生长。

2.1.2 植被指数与气候因子趋势分析

1982—2015年大兴安岭林草交错区植被NDVI与气候因子的变化趋势可知(图 3), NDVI变化趋势总体平均变化趋势为-0.02/10a, 即大兴安岭林草交错区植被呈总体上退化的趋势。大兴安岭林草交错区的年均CLD呈增加趋势(0.42%/10a), 一定范围内增加云量有助于增强大气对地面的逆辐射使得昼夜温差降低, 有利于植被的生长, 但当云量过多时减少到达地面的太阳辐射从而降低温度, 阻碍植被生理过程(光合作用)。此外TMN、TMP和TMX呈增温趋势, 其平均变化趋势分别为0.13℃/10a、0.16℃/10a和0.20℃/10 a。PET、VAP和PRE平均变化趋势分别为0.38 mm/10a、0.29 hPa/10a和-16.3 mm/10a, 即潜在蒸发量和水汽压呈增加趋势, 降水呈减少趋势。

2.2 大兴安岭林草交错区气候因子对植被指数的影响分析 2.2.1 气候因子与植被指数相关分析

气候因子与植被NDVI的相关分析结果如图 4所示, CLD与植被NDVI负相关(R=-0.21), 其中24.98%的CLD与植被NDVI显著负相关, 仅0.49%显著正相关, 随CLD的增加植被NDVI降低显著的像元主要分布在东部森林。PET与植被NDVI的平均R为0.004, PET与植被NDVI显著正相关的面积占总面积的18.60%主要分布在新巴尔虎旗西部、鄂温克族自治区、陈巴尔虎旗与额尔古纳市西部草原地区, 显著负相关占总面积的16.01%主要分布在而根河市和牙克石市及额尔古纳市东部地区的东部林地。与此相反, PRE和VAP与NDVI的平均R分别为0.03和0.18, PRE与植被显著正相关占总面积的19.55%, 显著负相关为5.31%, VAP与植被显著正相关占总面积的21.56%, 显著负相关为0.08%。TMN、TMP和TMX与NDVI的平均R分别为0.01、0.02和0.04, 显著相关像元主要分布在西部草地仅占总面积的1.81%、5.23%和12.74%。TMX对NDVI的影响高于最低温度和平均温度, 可能原因是白天温度和夜晚温度对植被碳同化和消耗是不对称的, 并且多数植被光合作用发生在白天, 因此植被NDVI对最高温度更敏感^[21—22]。总体而言, 西部草地NDVI受指示水分和湿度的指标(PET、PRE和VAP)影响高于东部森林, 而温度对大兴安岭林草交错区的影响低于降水。

2.2.2 气候因子与植被指数的偏最小二乘回归分析

利用偏最小二乘回归对大兴安岭林草交错区植被NDVI进行分析结果如图 5所示。结果表明, 回归的平均残差值和H指数分别为0.014和0.53, 其中H指数> 0.5表明人类活动对该区域植被的影响随着H指数的增加而增强, 并且陈巴尔虎旗、额尔古纳市和牙克石市交界地带残差H指数在0.75—0.97范围内, 表明该地区人类活动对植被的影响显著高于气候因子。植被NDVI的最大影响因子分布可知, 林草交错区植被NDVI主导气候因子的重要性排序为VIP _CLD > VIP _PRE > VIP _PET > VIP _VAP > VIP _TMP > VIP _TMN > VIP _TMX, 其中以云量为主导影响因子占总面积的54.20%, 并且该区域的NDVI较高(林地), 主要原因是云可以吸收和反射太阳辐射降低地面温度, 吸收和放射长波辐射增暖地面, 同时可以改变形态影响水分循环^[23]。

2.3 植被指数影响因子实证分析 2.3.1 气候因子与植被指数的相关分析

1957—2018年额尔古纳和陈巴尔虎旗年均温度分别为-2.4℃和-1.5℃, 呈波动上升趋势, 分别为0.3℃/10a和0.5℃/10a(图 6)通过显著性检验(P < 0.05), 说明其增温趋势显著并且陈巴尔虎旗的增温速率高于额尔古纳。1957—2018年额尔古纳和陈巴尔虎旗年平均降水量为355.2 mm和318.5 mm, 年降水量总体呈不显著的波动下降趋势(P>0.05), 下降速率分别为8.1 mm/10a和2.3 mm/10a。1957—1975年两个地区年降水量呈下降趋势, 数值变化在220—480 mm之间波动。随后在1975—1990年降水量呈上升趋势, 数值变化在220—542.9 mm之间波动, 平均值为371.9 mm和339.1 mm。1990—2018年额尔古纳和陈巴尔虎旗年降水量呈显著下变化趋势, 分别为32.0 mm/10a和48.1 mm/10a, 并且陈巴尔虎旗的降水量下降速率高于额尔古纳。总之, 额尔古纳和陈巴尔古纳温度呈上升趋势, 而降水量呈减少趋势。

如图 7所示, 额尔古纳和陈巴尔虎旗的植被NDVI平均值处于显著下降趋势, 其中, 额尔古纳植被NDVI下降的速率更大。NDVI与年均温度和降水的相关分析结果表明, 植被NDVI随温度的增加而降低, 而植被NDVI随着降水的增加而增加, 并且仅在额尔古纳的年降水量与NDVI显著相关, 温度与NDVI的相关性均未通过显著性检验。温度和降水对植被NDVI的解释率较低在5%—33%范围内, 相较于温度的解释率, 降水的解释率更高。可能原因是干旱对该地区的植被NDVI影响作用高于温度。

2.3.2 人类活动因子对植被指数的影响分析

(1) 土地开垦面积趋势分析

我们以额尔古纳的上库力农场和陈巴尔虎旗的陶海牧场为主要调查对象, 统计分析了自20世纪50年代两个国有农场开垦面积的变化趋势。如图 8所示, 1956年上库力农场成立后, 开垦面积急速增长主要原因是国家处于建国初期, 实施鼓励生育政策, 人口数量剧增, 粮食短缺, 为解决粮食产量问题, 区域内实施以生产粮食为主的土地利用政策, 为此国有农场开始大面积的开垦草原, 导致1968年上库力农场的总开垦面积达到2×10⁴ hm², 随着退耕还林等政策的实施开垦面积增长速度降低甚至出现下降。根据增长曲线拟合结果发现1982年该农场的增长速率达到最大1171.19 hm²/a, 随后增长速率下降, 但仍每年在400 hm²速率以上。近几年区域内严格实施退耕还林、还草政策, 农场的开垦面积呈逐年下降的趋势。

相较于上库力农场, 陶海牧场总开垦面积的波动更为剧烈, 具体来说, 从1959到1961年中国经历三年饥荒, 为应对粮食危机, 政府提出“粮食优先”政策, 导致1961年牧场总开垦面积从800 hm²突增至3000 hm², 随后在1963年又下降至1000 hm²。直到1995年国家开始大量投资农业基础设施以及实施粮食保护价收购政策, 使区域内农业发展水平显著上升, 以致陶海牧场的总开垦面积表现逐年显著波动式增加。根据增长趋势分析发现在1990年该农场的增长速率达到最大270.38 hm²/a, 随后增长速率迅速降低。

(2) 家畜密度趋势分析

1955—2018年的平均家畜密度变化趋势如图 9所示, 国有农牧场和个人牧户总体表现为：1955—2004年平均家畜密度维持在0.74只/hm²标准羊单位, 而在2005年平均家畜密度增加到1.71只/hm²标准羊单位, 此后呈现逐年波动式上升趋势, 直到2019年达到4.41只/hm²标准羊单位。根据自然条件限制和家畜之间的竞争, 利用增长曲线对国有农牧场和个人牧户进行了曲线拟合。研究发现：陶海牧场在2009年达到最大畜牧密度增长速率为0.26只羊单位hm^-2 a^-1, 畜牧密度为3.11只羊单位hm^-2, 上库力农场在2006年达到最大畜牧密度增长速率为0.07只羊单位hm^-2 a^-1, 畜牧密度为1.35只羊单位/hm²。然而, 由于研究时段和增长状态的影响, 苏沁农场和个人农牧户畜牧密度不满足增长模型必要条件, 因此, 不存在最大畜牧密度增长速率。

(3) 人类活动因子与植被指数的相关分析

家畜密度和开垦面积与植被NDVI的关系如图 10所示, 结果表明：额尔古纳家畜密度和开垦面积与植被NDVI显著负相关, 并且两者对植被NDVI的解释率在37%—48%范围内。然而陈巴尔虎旗家畜密度与植被NDVI显著负相关, 与开垦面积弱负相关, 两者对植被NDVI的解释率在27%—42%范围内。由此可知, 家畜密度对NDVI的影响高于开垦面积, 主要原因包括两个方面, 首先, 家畜密度增加影响植被的生长、发育和繁殖能力^[24—25], 从而引发优良牧草群落衰退甚至被其他杂草和毒草代替^[26—27], 导致群落生产力下降和生态系统不稳定性增加^{[24, 28]}。其次, 家畜密度增加导致土壤理化性质发生改变, 如：土壤结构变化降低水分渗入、土壤有机质含量和孔隙度减少加速土壤侵蚀, 进而导致土壤养分流失植被生产力降低甚至出现水土流失现象^{[24, 29]}。

2.3.3 植被指数影响因子的重要性分析

利用变量投影重要性指标对温度、降水、家畜密度和开垦面积对植被NDVI的影响程度进行重要性排序, 发现陈巴尔虎旗和额尔古纳市的家畜密度和开垦面积对植被NDVI的影响均高于温度和降水, 而降水对植被NDVI的影响均高于温度。此外, 额尔古纳的温度、降水、家畜密度和开垦面积的VIP排序为VIP_温度 < VIP_降水 < VIP_家畜密度 < VIP_开垦面积(0.23 < 0.6 < 0.61 < 1.79)；陈巴尔虎旗的VIP_温度 < VIP_降水 < VIP_开垦面积 < VIP_家畜密度(0.22 < 0.54 < 1.04 < 1.52)。额尔古纳VIP_家畜密度 < VIP_开垦面积而R_家畜密度 > R_开垦面积主要原因是相关性强调的是两者对植被NDVI的因果关系(直接或间接), 重要性强调的是影响植被NDVI的必要性侧重点不同因此结果可能存在差异。总之, 人类活动因子(畜牧密度和开垦面积)对大兴安岭林草交错区的植被NDVI影响高于气候因子。气候因子中降水的影响作用高于温度, 降水是两个地区植被NDVI的主要气候影响因子。

3 讨论

东北大兴安岭林草交错区植被显著退化的区域主要分布在草地向林地过度的地区, 气候因子对植被NDVI的影响大小存在差异, 其中PET、PRE和VAP存在明显的地区差异, 在NDVI的低值区域(草地为典型植被)水分的影响作用显著, NDVI随着水分的增加而增加。此外, PRE的重要性指数在草地植被类型中的影响最大(图 5), 并且指示水分的潜在蒸散、降水和水汽压分别呈增加、降低和增加趋势, 这与已有研究草地植被影响因素的结论一致, 如：中国的气候呈现变干旱的趋势^[30], 降水是造成东北地区草地退化的主要原因^[31], 降水对不同草地植被的影响高于温度^[32—33]。与此相反, NDVI较高的区域(针叶林、混交林地为典型植被), NDVI随着水分的增加而降低并且植被呈变绿趋势。可能原因是该地区的温度在该地区年度波动较大, 并且云量和温度的影响作用高于降水, 其中最低温度、最高温度和平均温度与NDVI呈显著正相关, 而温度变量呈现增加的趋势, 该结论与已有的温度升高促进中高纬度植被的生长的结论一致^[34—35], 因此, 东部林地的NDVI呈增加趋势。

温度变量与植被指数的相关系数在干旱地区(大兴安岭林草交错区西部)呈负相关, 而与降水相对较高的东部地区呈正相关, 并且最高温度对植被NDVI的影响高于最低温度和平均温度。可能的原因是在干旱地区限制植被生长的因子是降水, 温度升高会加剧干旱从而引起植被缺水、枯黄甚至脱落, 从而降低NDVI；而在水分较充足的地区, 温度是主要限制因子, 温度升高促进植被生长进而提高NDVI^{[22, 36—37]}。此外, 最高温度直接影响数植被的光合作用, 影响叶绿素、蛋白质、等有机质的形成；最低温度通过影响呼吸作用影响植被生长, 而最高温度也影响呼吸作用^[38], 因此, TMX对植被的影响高于最低温度和平均温度。

植被NDVI的变化是气候因子和人类活动因子共同作用的结果, 然而气候因子的H指数(图 5, H > 0.5)表明大兴安岭林草交错区的人类活动影响较大, 实证分析结果进一步表明开垦面积呈波动上升趋势, 国有牧场和个人农牧场的放牧密度呈增加趋势, 国有农场呈波动上升趋势, 而个人农牧场呈稳定上升趋势。开垦面积和家畜密度波动上升的可能是由于天然草地变为人工草地和农田的经济成本低, 综合政策对农业和畜牧业的调控共同影响的结果^[39—40]。此外, 开垦面积和家畜密度与植被NDVI呈显著负相关, 该结论与杜加强等^[41]研究人类活动和气候变化对植被NDVI的结果一致, 主要原因是过度放牧导致家畜践踏和采食影响地表生态系统物质循环, 降低植被的生产力^[42—43]。开垦面积和家畜密度重要性排序VIP与相关分析结果发现家畜密度对植被NDVI的解释率较高, 而开垦面积的重要性

4 结论

利用遥感区域分析和地面实证分析对大兴安岭林草交错区的植被变化趋势及影响因子进行综合分析, 通过重要性指标定量阐述人类活动因子(开垦面积和家畜密度)与气候因子(温度和降水)的相对影响, 为改善林草交错区生态系统提供理论基础。1982—2015大兴安岭林草交错区的植被呈退化趋势(-0.02/10a), 气候呈变暖变干旱趋势, 总体上降水对大兴安岭林草交错区生态脆弱区植被的影响显著高于温度。人类活动因子对植被退化的影响作用显著高于气候因子(VIP_{人类活动因子} > VIP_降水 > VIP_温度)。因此, 气候变暖背景下, 人类活动因子对植被的影响作用不容忽视, 合理调控农业和畜牧业是改善林草交错区植被生态系统稳定可持续的重要途径。

致谢: 东北师范大学生命科学学院高英志教授提供部分数据资料, 特此致谢。