植被是地球系统中的活跃成员, 对生境变化较为敏感。植被在类型、数量或质量方面必然对陆地生态系统的变化有所响应, 被认为是监测生态环境变化的综合指示器^[1]。植被变化长期来都是诸多学者关注和研究的重点问题^[2-5]。而植被覆盖度(Fractional Vegetation Cover, FVC)不仅是描述地表植被生长状况的重要指标, 也是反映生态系统变化的重要变量^[6]。FVC是以归一化植被指数(Normalized Difference Vegetation Index, NDVI)为基础改进而提出的^[7], 指植被冠层在单位面积内垂直投影的占比^[8], 反映了植被的茂密程度, 对于NDVI对高覆盖度的植被易饱和、而对低覆盖度的植被难区分的不足有所改善^[9], 可以更好地反映植被覆盖动态变化特征。

近年来, 植被覆盖度被广泛应用于区域植被变化研究^[10-15]。其中, 基于NDVI数据, 采用像元二分模型估算植被覆盖度的方法效果良好, 已有丰富的成果。如刘春静等^[16]以中国新疆及中亚五国为研究区, 探讨了大陆干旱性气候区近14年的植被覆盖度时空演变特征。张文强等^[5]分析了黄土高原15年植被覆盖变化情况, 黄土高原作为温带大陆性半干旱半湿润季风气候区的典型, 植被覆盖度呈增加趋势。张学玲等^[17]以江西省武功山为例研究了中亚热带湿润气候区植被覆盖度的时空变化情况, 发现植被覆盖总体呈上升趋势, 且与地形因子有密切关系。赵婷等^[18]将秦岭作为研究区, 探究了暖温带半湿润气候与亚热带湿润气候过渡带的植被变化情况, 发现该区近17年植被覆盖度良好, 整体呈上升趋势, 植被覆盖变化具有明显的地形分异效应。郑朝菊等^[19]研究西南地区, 包括亚热带季风气候、热带季雨林气候及高原气候近15年的植被覆盖度时空变化及与气候、人类活动的关系。

上述已有研究对认识各气候区的植被变化和驱动因素具有重要意义, 但对温带大陆性半干旱季风气候区植被变化的研究还不多见。陈效逑, 王子玉等^[20-21]以内蒙古为例研究了气候变化和人类活动对植被变化的影响。但已有研究表明植被变化不仅与气候, 人类活动显著相关, 且对地形因子高度响应^[22-28]。影响植被分布的生境因子众多, 其中地形是最基本的影响因子。地形因子(坡度、坡向和海拔等)不仅会影响气候条件的空间差异, 也会影响人类的生产与生活, 进而影响植被的空间分布格局, 因此研究地形因子与植被覆盖度变化的关系可以很好地揭示植被的空间分布规律^[29]。

内蒙古大青山位于温带大陆性半干旱季风气候区, 是对全球气候变化最为敏感的地区之一, 2008年成为国家级自然保护区, 对于衡量气候变化大背景和重大人为活动影响下区域生态环境变化具有重要的指示作用。因此, 本文以内蒙古大青山为例, 选用4期Landsat影像, 计算植被覆盖度, 分析植被覆盖度时空变化以及地形对植被覆盖分布格局的影响, 以期为温带大陆性半干旱季风气候区植被监测与保护提供参考。

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

大青山位于109°47′— 112°17′E, 北纬40°34′—41°14′N之间, 面积为3918.9km²(图 1)。地处阴山山脉中段, 是蒙古高原草原区与黄土高原草原区的分水岭。大青山海拔在1030—2338m之间, 属于中山山地。山体呈东西走向, 北部较平缓, 地貌南北差异大。气候条件为温带大陆性半干旱季风气候, 冬季寒冷漫长, 夏季温和, 降水集中在夏季, 日照充足, 历年平均气温6.7℃, 年降雨量350—400mm。随地势的升高, 土壤类型呈现出明显的垂直分布特征, 植被类型也呈现出灌丛草原-森林灌丛草原-森林草原-山地草原的分布规律。大青山因其重要的生态服务功能, 受到各级政府的极度重视。2008年, 国家批准建立内蒙古大青山国家级自然保护区。2013年, 内蒙古自治区通过了《内蒙古自治区大青山国家级自然保护区条例》, 自此大青山的生态保护和建设工作进入了一个全新的时期。

1.2 数据来源

研究区6月至8月降水集中, 为减小云量对研究分析的影响, 特选取9月晴朗无云且日期相近的四期数据。本文采用的遥感数据为2000、2006、2011年landsat-5 TM影像和2017年landsat-8 OLI影像, 轨道编号为126/031、126/032、127/032。四期数据云量均小于10%, 质量良好。此外, 研究还用到DEM数据(30m分辨率), 用于提取地形信息。所有数据都来源于地理空间数据云平台(http://www.gscloud.cn/), 投影坐标系为UTM-WGS 84。

1.3 数据处理

在ENVI 5.3中对TM和OLI数据做几何配准, 校正误差均严格限制在0.5个像元以内, 重采样分辨率为30m。然后进行辐射定标和大气校正, 先将影像灰度值转化为辐射亮度值, 再进行大气校正以去除大气对地物辐射的影响, 呈现地物较为真实的反射光谱特征。标注大气校正后数据中的负值和大于100的高值, 并改正这些异常值。对每一期的多幅影像, 拼接后用大青山边界矢量裁剪, 得到研究区四期遥感数据。提取大青山的DEM, 并在ArcGIS 10.6中使用3D分析模块, 生成研究区的高程图、坡度图及坡向图。通过重新划分高程、坡度、坡向等级, 再与研究区植被覆盖度分布图叠加, 探讨大青山地形因子对植被覆盖度分布格局的影响。

1.4 研究方法

已有的诸多研究表明将NDVI代入像元二分模型计算得到的植被覆盖度, 弥补了NDVI研究高植被覆盖和低植被覆盖效果欠佳的缺点, 监测植被覆盖变化效果良好, 可以更好地表征区域植被的生长和变化情况^[30-32]。因此, 采用像元二分模型计算研究区植被覆盖度, 并在ArcGIS10.6中对4期植被覆盖度图进行差值运算, 分析研究区植被覆盖度时空动态变化, 进而对地形因子(高程, 坡度, 坡向)进行重分类, 探测植被覆盖度的地形分异规律。

1.4.1 植被覆盖度计算和定级方法

采用像元二分模型, 基于NDVI数据来计算植被覆盖度。NDVI一般定义为近红外波段(0.7—1.1 μm) 与可见光红光波段(0.4—0.7μm) 反射率之差与反射率之和的比值。NDVI值的范围为[-1, 1], 植被覆盖越好, 值越大, 值为1表示该像元为纯植被覆盖, 值为0表示该像元为纯裸土覆盖。

而像元二分模型假定一个像元的信息由植被覆盖和裸土覆盖提供。遥感影像上任一像元的遥感信息记为S, 若像元信息完全由植被提供的纯像元遥感信息记为S_veg, 像元信息完全由裸土提供的纯像元遥感信息记为S_soil, 则植被覆盖度的计算公式为：

(1)

将NDVI带入像元二分模型, 得到基于NDVI的植被覆盖度计算公式：

(2)

式中, NDVI为每个像元的归一化植被指数值, NDVI_veg为纯植被覆盖像元的归一化植被指数值, NDVI_soil为纯裸土覆盖像元的归一化植被指数值。

实际计算中, 受数据预处理影响, 部分像元的NDVI值在[-1, 1]之外, 这部分NDVI值为异常值, 为便于后期计算, 本文将异常值转化为背景值, 公式如下：

(3)

式中, NDVI_g表示原始NDVI值, NDVI₁为小于-1的原始NDVI值, NDVI₂为大于1的原始NDVI值, NDVI₃为[-1, 1]的原始NDVI值, NDVI表示去除异常值后的NDVI值。

受年份、大气等条件影响, NDVI_soil和NDVI_veg随时间和空间的不同会产生变化, 因此取累计直方图的5%和95%为置信度区间, 确定研究区的NDVI_soil和NDVI_veg值^[33-34](表 1)。

FVC值的范围为[0, 1], 为分析植被覆盖度动态变化情况, 结合研究区实际情况, 将FVC等间距划分为3个等级^[17]。1等：FVC < 0.3, 低植被覆盖度；2等：0.3≤FVC < 0.6中植被覆盖度；3等：FVC≥0.6, 高植被覆盖度。生成植被覆盖度分级图并统计面积。

1.4.2 植被覆盖度动态分析方法

采用差值法分析植被覆盖度动态变化特征, 公式如下^[35]：

(4)

ΔFVC表示植被覆盖度变化等级, FVC₁和FVC₂分别表示前后两个时期的植被覆盖度等级。ΔFVC值表征的含义如下：ΔFVC=2时, 记为重度改善；ΔFVC=1时, 记为轻度改善；ΔFVC=0时, 记为未变化；ΔFVC=-1时, 记为轻度退化；ΔFVC=-2时, 记为重度退化。

1.4.3 植被覆盖度对地形因子响应分析方法

地形特征是一个多维变量, 不同高程、坡度和坡向具有的不同外部形态造成不同的水热条件分布和养分移动堆积, 在某一特征尺度上不同的地形特征对植被覆盖度的影响强度会存在空间上的差异^[17]。为定量研究植被覆盖度对地形因子的响应, 对大青山高程、坡度、坡向数据重分类(图 2), 再对植被覆盖度分区统计, 分析高程、坡度、坡向各区间内的植被覆盖度分布规律。

采用等间距法将研究区高程重分类为10个等级：< 1200m；1200—1300m；1300—1400m；1400—1500m；1500—1600m；1600—1700m；1700—1800m；1800—1900m；1900—2000m；>2000m, 然后统计各高程带的植被覆盖度。

根据《土壤侵蚀分类分级标准》(SL 190—2007)中的临界坡度分级法, 将研究区的坡度重分类为6级：0°—5°；5°—10°；10°—15°；15°—25°；25°—45°；45°—90°。按照划定的坡度等级统计大青山的植被覆盖度分布规律。

根据国土资源部颁布的《第二次全国土地调查技术规程》, 将研究区的坡向分为平地(0°)、阴坡(315°—45°)、半阴坡(45°—135°)、阳坡(135°—225°)、半阳坡(225°—315°) 五类。统计研究区各坡向的植被覆盖度。

2 结果与分析 2.1 大青山植被覆盖度空间分布格局

基于NDVI的像元二分模型, 计算得到2000—2017年大青山的平均植被覆盖度(图 3)。由图 3知, 大青山植被覆盖度空间分布特征为“东部高, 西部低；南部高, 北部低”。2000—2017年, 大青山的植被覆盖类型以中高植被覆盖度(FVC≥0.3)为主, 平均占64.19%的研究区域。

为详细研究各时期植被覆盖度的分布情况, 依据1.4.3的等级标准划分大青山4个时期的植被覆盖度, 得到如图 4和表 2的2000、2006、2011、2017年大青山植被覆盖度等级图及分布特征。

从图 4来看, 植被覆盖度呈现东北高西南低的空间分布特征。2000年, 中高植被覆盖主要在乌兰察布段和呼和浩特段, 低植被覆盖以包头段和呼和浩特段南部为主。2006年和2011年, 高植被覆盖面积较小, 中植被覆盖主要分布在乌兰察布段和呼和浩特段, 低植被覆盖以包头段为主。2017年, 高植被覆盖聚集在东北部；植被覆盖度等级自东北向西南依次从高、中到低植被覆盖过渡。

由表 2知：虽然不同等级的植被覆盖面积占比都随时间呈动态变化, 但总体来看中植被覆盖区在研究区总面积中占比最大, 其次是低植被覆盖区的比例, 而高植被覆盖区占比最小。此外, 中低植被覆盖类型的面积比例在四个时期分别占到了77.70%、92.94%、99.67%和80.12%, 是研究区主要的植被覆盖度类型。这与研究区地理区位密切相关, 大青山地处大陆性半干旱季风气候区, 植被类型以北方草原为主, 所以中低植被覆盖区域广。而高植被覆盖主要由山地森林贡献, 主要受人为因素影响。

18年间, 低植被覆盖区的面积先减后增再减, 总体呈减少趋势。中植被覆盖区面积变化最大, 共增加24.47%。高植被覆盖区总体变化较平稳。说明大青山地区十八年间植被向着恢复的方向发展。这是因为大青山地区自2013年施行《内蒙古自治区大青山国家级自然保护区条例》以来, 依据完善的保护条例, 加强保护区内各类开发建设活动的监督管理, 清理整顿工矿企业, 加大对私自放牧的管理力度, 加强森林草原防灭火工作, 加快了大青山植被的恢复。

2.2 大青山植被覆盖度动态变化

为研究大青山地区18年间植被覆盖在空间上的变化情况, 在ArcGIS 10.6中对四期植被覆盖度等级图进行叠加运算, 得到如图 5和表 3、表 4、表 5、表 6所示的大青山植被覆盖空间分布格局及其演变特征。

由图 5和表 3知：2000至2006年大青山植被趋于退化。高植被覆盖区转出693.8km², 仅转入100.34km², 损失面积大。这一时期大青山保护范围主要在呼和浩特段, 其余区段保护力度小, 开发力度大。研究区内矿山企业挖掘地表类的开采直接破坏地表植被, 山区农牧民开垦放牧。还有, 用于山林防火、生态保护的经费不足。这些因素综合作用导致了大青山植被覆盖度下降。

由图 5和表 4知：2006至2011年大青山植被趋于退化。中高植被覆盖面积减小, 转化为低植被覆盖度。原因有二：一方面, 降水及蒸散发量制约当年植被覆盖, 统计年鉴显示2011年降水量较往年有所下降；另一方面, 人类活动在短期内对区域植被生长影响显著。大青山中南部紧邻呼和浩特市的人口聚居区, 受人类生产生活影响, 该区域植被退化严重。而西部受人为干扰较小, 植被变化较稳定。

由图 5和表 5知：2011至2017年大青山植被趋于改善。低植被覆盖大量转化为中高植被覆盖, 中植被覆盖度主要转出为高植被覆盖。6年间植被覆盖度增加, 生态向好发展。2013年《内蒙古自治区大青山国家级自然保护区条例》的颁布施行, 使大青山地区保护工作有例可循。且大青山管理部门响应国家的生态建设, 紧随其后制定了针对矿山企业, 私自放牧以及生态建设等的数个方案, 加强执行力度, 使大青山生态恢复迅速, 中高植被覆盖面积大幅增加。

由图 5和表 6知：18年间大青山植被覆盖度总体呈上升趋势。中植被覆盖度的变化量最显著, 且为正值。低覆盖度呈负向变化, 高覆盖度变化量较小。低植被覆盖度主要的转出类型为中植被覆盖度, 中植被覆盖度主要的转出类型为高植被覆盖度, 说明18年间, 植被覆盖度由低、中覆盖度类型向中、高覆盖度类型转变, 变化趋势向好。

2.3 大青山植被覆盖度变化的空间差异性分析

图 6为大青山2000年至2017年植被覆盖度等级变化情况。18年间, 大青山西部地区以轻度改善为主。中部地区变化小, 以轻度改善为主, 有小面积的轻度退化。东部以轻度退化为主, 小范围内存在重度改善和重度退化现象。

植被覆盖等级变化以轻度改善面积最大, 达到1152.45km², 分布最广, 自东北向西南均有轻度改善趋势。轻度退化面积次之, 有395.29km², 主要分布在中部局地和东部。植被重度改善与重度退化区面积较小, 聚集性强。植被重度退化区域最小, 说明研究区植被无严重的恶化情况。总体上, 改善类型占32.46%, 退化类型占12.92%, 表明大青山18年间植被覆盖度总体改善。

2.4 大青山植被覆盖度的地形因子分异研究 2.4.1 大青山植被覆盖度海拔梯度差异

图 7表明, 大青山植被覆盖度随着海拔的升高持续增长, 在2000—2296m最高, 936—1200m最低。且在低于1500m的高程带, 植被覆盖度增长缓慢；1500m以上的区域增长速度较快。海拔不同, 区域的水热条件也会有相应改变, 从而影响植被生长。低海拔区水热条件较好, 但受到人类活动的干扰, 植被覆盖度较小。海拔高于1500m, 人类活动受到相当程度的限制, 植被覆盖好转。且大青山最高海拔2296m, 属于中山山地, 并不会使限制植物生长的热量条件受到影响。因此, 大青山植被覆盖随海拔升高而增大。

2.4.2 大青山植被覆盖度坡度差异

图 8表明, 随坡度增加, 大青山的植被覆盖度逐渐增高, 在坡度0°—5°最低, 45°—90°最高。这是由于坡度较缓的地区方便人类活动, 而坡度越高, 人类活动受到地形限制逐渐减少, 对植被生长的干扰也就越小, 因此植被覆盖度随坡度增加而呈上升趋势。

2.4.3 大青山植被覆盖度坡向差异

图 9表明, 阴坡植被覆盖度最高, 其次是半阴坡和半阳坡; 植被覆盖度最低的是阳坡。主要是因为阳坡接收的太阳辐射较其他坡向更强, 而研究区属于半干旱季风气候, 昼夜温差较大, 水分蒸发量较多, 植被分布以山地草甸为主。山地森林、灌丛主要分布在阴坡, 其次分布在半阴坡和半阳坡, 所以植被覆盖度的分布规律为阴坡＞平地＞阳坡。

3 结论与讨论

使用基于NDVI的像元二分模型计算植被覆盖度, 并对4期植被覆盖度进行分级、计算平均覆盖度、差值运算, 以及对地形因子(高程, 坡度, 坡向)重分类, 据此获得研究区的植被覆盖度时空演变过程及其与地形因子的响应情况。2000—2017年, 大青山植被总体覆盖状况良好且呈上升趋势, 植被覆盖空间分布、变化趋势及其与地形因子的相关性如下：

(1) 大青山以中高植被覆盖度为主, 64.19%的区域植被覆盖类型为中高植被覆盖度。不同时期的植被覆盖度有所波动, 18年来大青山植被覆盖度先降低后升高, 总体呈上升趋势。前期植被覆盖度降低成因是大青山地区的诸多经济开发活动, 如旅游业的发展、工矿企业的开设、私自放牧等。而自大青山国家自然保护区成立, 尤其2013年保护区条例通过后, 大青山保护不断完善, 生态逐渐恢复, 植被覆盖度升高。

(2) 植被覆盖度空间格局特征明显, 总体上呈“东部高, 西部低；南部高, 北部低”的分布特征。中高植被覆盖度集中在大青山呼和浩特段南部和乌兰察布段, 而低植被覆盖度主要分布在西段山体。

(3) 空间上, 植被轻度改善分布最广, 自东北向西南均有轻度改善趋势。轻度退化主要分布在中部局地和东部。重度改善与重度退化区面积较小, 聚集在东部小范围内。数量上, 研究区植被改善类型占32.46%, 退化类型占12.92%, 其中以轻度退化为主。说明研究区18年间植被覆盖度为改善趋势。

(4) 植被覆盖度的地形分异变化明显：植被覆盖度随海拔升高呈增加趋势, 在2000—2296m高程带最高。植被覆盖度与坡度正相关, 坡度越高, 植被覆盖度越大。植被覆盖度在不同坡向上差异明显, 呈现阴坡＞半阳坡＞半阴坡＞平地＞阳坡的分布规律。

当然, 本文还存在一些不足之处。首先, 在数据源上, 本研究采用4期TM/OLI影像图提取NDVI, 未使用多源NDVI数据集进行相互验证。另外, 本研究中对植被覆盖度的地形分异规律进行了分析, 但地形对植被覆盖的影响还可以引入植被类型、土地利用、地质数据等来进一步展开。且限于时间与篇幅, 本文未能进行地形面积差异修正以评价地形对植被变化的影响效应, 这将成为下一阶段的研究工作。