煤炭资源粗放式、大规模、高强度的开采, 改变了矿区原有的地形地貌, 造成了矿区地表下陷、植被土壤退化、水土流失, 引发了矿区生态环境进一步恶化的风险, 增加了矿区生态环境治理的难度^[1-2]。在矿区生态环境治理中, 植被覆盖度在一定程度上可以量化地表植被状况反映生态恢复程度, 从而为矿区复垦和生态环境重建提供参考依据^[3]。

目前, 植被覆盖度的获取方法有实测法的目估法、采样法、仪器法和模型法；遥感监测方法的回归模型法^[4-5]、植被指数法^[6-7]和像元分解模型法^[8-9]。地表实测法精度高、但范围小且费时费力；遥感监测法精度不及实测法, 但监测范围广, 时间序列长, 能够更好地揭示出植被覆盖度的时空变化特征。

已有研究表明植被覆盖度的空间分布和随时间变化的趋势受气候、地理、人为及因素组合效应的影响^[10-12]。气候因素中的降雨量和温度是对植被覆盖度影响最复杂的2个主导因素。如徐占军等^[13]研究发现在采煤区域, 气候变化对矿区植被覆盖度的影响具有主导作用；穆少杰等^[14]研究表明在年际水平上, 植被覆盖度与降雨量的关系比与温度的关系更为密切；在月际水平上降雨量和温度对植被覆盖度的影响作用相当。随着研究区域和植被类型的改变, 植被覆盖度与降雨量、温度之间统计关系也在发生变化, 比如, 出现正向协同变化、无统计关系、存在互作效应、或负向协同变化^[15-17]。植被覆盖度与降雨量及温度之间出现复杂的统计变化关系, 一方面可能与植被类型没有分类有关, 另一方面与植被所在的高程、坡向和坡度等地形因素的影响有关^[18-19]。如张学玲等^[20]研究发现武功山山地草甸植被覆盖度分布阳坡大于阴坡, 随着坡度的上升先增加后下降, 随高程的增加出现波浪式下降；陈建军等^[21]研究表明疏勒河源区高寒草地坡度、坡向是限制植被分布的主要因子, 坡度越小, 平均植被覆盖度越大, 随坡向由无坡向、阴坡、半阴(阳)坡到阳坡平均植被覆盖度不断减少；贾铎等^[22]研究表明植被覆盖度变化与气候和地形因素具有相关性, 长期植被覆盖度的变化受气候因素的影响, 短期植被覆盖度的变化主要受人为因素的影响。植被覆盖度除了受自然气候和地理因素的影响外, 同时还不同程度地受人为因素的影响, 如煤炭等资源开采对植被破坏及复垦、土地利用方式的改变等, 而且自然因素和人为因素常常交错在一起共同对植被覆盖度发挥作用。Laura J. Sonter等^[23]对矿区土地利用变化过程进行了研究, 揭示采矿活动直接和间接影响着矿区植被变化和生态环境质量；张世文等^[24]分析了露天矿区随着时间变化采矿和生态恢复活动与植被覆盖度间具有一定的响应关系；黄翌等^[25]研究发现大同矿区植被变化的内在因素是煤炭开发的扰动作用。

综上所述, 植被覆盖度时空变化受诸多因素影响, 已有研究得出的结论大多是总的植被覆盖度的平均效应, 具有较大的不确定性, 且研究结果各异。如何解释这些研究结果的多样性, 且进一步认识植被覆盖度时空变化和不同的植被类型与不同因素及其互作效应的关系是本研究拟解决的主要问题。为此, 本研究选取了3个位于山西省不同气候带、不同地形地貌下的典型井工煤矿区与露天煤矿区, 利用2001—2016年7—9月份的Landsat遥感影像数据获取研究区地表 30 m×30 m网格的植被覆盖度数据, 采用线性回归斜率法和相关分析方法分析植被覆盖度的时空变化特征和影响因素, 并进一步讨论自然因素对矿区不同植被类型的植被覆盖度变化的综合效应。

1 研究区概况

山西省位于华北平原西部的黄土高原, 海拔多在1000 m以上, 山地和丘陵面积占总面积的80%以上。山西属温带大陆性季风气候, 年平均气温介于4—14 ℃之间, 年降水量在400—600 mm之间, 大部分地区为半干旱地区, 东南部的部分地区为半湿润地区。

综合考虑山西省的气候、地形地貌、土壤类型、植被类型、矿区开采模式等因素, 选取了3个位于不同气候带、不同地形地貌、不同采矿方式的典型煤矿区(图 1)。其中, 长河矿区位于山西省南部泽州县内, 属暖温带半湿润大陆性季风气候, 为丘陵地貌, 严重缺水, 土壤以褐土为主, 煤矿分布集中, 采煤量大。离柳矿区位于山西省中西部柳林县内, 暖温带半干旱大陆性季风气候区, 西北黄土丘陵沟壑地貌, 大部面积为贫水区和极贫水区, 土壤以褐土为主, 煤炭储量多, 分布广。平朔矿区位于山西省西北部朔州市内, 属温带半干旱大陆性季风气候区, 为黄土丘陵地貌, 地表水系比较丰富, 以栗钙土为主, 植被主要呈农业耕作景观, 煤炭储量多。

2 研究方法 2.1 数据源与预处理

研究中涉及到Landsat影像数据和DEM数据, 均来源于中国科学院遥感与数字地球研究所及地理空间数据云数据共享平台(http://www.gscloud.cn)。本研究3个矿区涉及到的遥感数据如表 1所示。每期遥感影像的研究区范围内云覆盖率均为0, 空间分辨率均为30 m, 空间参考坐标为WGS 1984 UTM Zone 49N。影像质量较好, 且时相选择都处于植被生长期的7月、8月、9月。对获取的数据分别进行了去条带处理、辐射定标得到表观反射率、按照研究区域进行裁剪处理。

DEM数据的空间参考为西安80坐标系, 转换为WGS84, 对转换后的DEM数据基于影像数据进行几何纠正, 并保证纠正误差在0.5像元之内。通过ArcGIS 10.0基于DEM数据求得坡度、坡向数据。2001—2016年的温度、降雨量数据通过国家气象科学数据共享服务平台(http://data.cma.cn/)下载获得；植被类型数据来自项目组前期的研究成果；煤矿开采生产信息来自山西省国土资源厅发表的山西省煤炭资源相关数据。

2.2 研究方法 2.2.1 地形调节植被指数

本研究采用江洪等人^[26]提出的植被指数计算方法, 用红光波段数据构建阴影植被指数(SAV)与优选比值植被指数(RVI)形成线性组合, 构建基于波段比模型的地形调节植被指数TAVI, 见公式(1), 以此来消除地形对植被指数计算的影响。

(1)

式中, B_nir表示遥感影像近红外波段数据；B_r表示遥感影像红光波段数据。fΔ为地形条件因子, fΔ的确定采用“极值优化”算法^[27], 通过程序实现。首先选取具有强烈地形影响的样区, 通过监督分类分为阴坡和阳坡2类, 然后设计循环程序, 令fΔ从0开始, 以0. 000 1为间隔, 依次递增, 同时考察阴坡部分TAVI的最大值与阳坡部分TAVI的最大值, 当两者相差小于0.01时, 退出循环, 得到fΔ优化结果。

对TAVI数据进行极差标准化处理

(2)

式中, TAVI′是标准化后的植被指数值, 取值0—1；TAVI是标准化前植被指数值；TAVI_max是标准化前植被指数最大值；TAVI_min是标准化前植被指数最小值。

2.2.2 植被覆盖度计算

在地形调节植被指数的基础上, 采用像元二分模型获得研究区植被覆盖度。

(3)

式中, FVC表示植被覆盖度；TAVI_soil为裸土或无植被覆盖区的TAVI值；TAVI_veg表示纯植被像元的TAVI值。根据经验, TAVI_soil与TAVI_veg分别取累积频率5%和95%处的值。

2.2.3 植被覆盖度的时空变化特征

采用线性回归斜率法^[28], 研究像元尺度下植被覆盖度随时间的变化趋势。对每个像元, 在选定的n个年份下分别计算其植被覆盖度, 计算植被覆盖度对年份的回归系数(公式4)。如果回归系数大于一定的临界值, 表明该像元对应的实际矿区位置上植被覆盖度逐年增加；介于一定的范围内, 表明植被覆盖度变化不明显；小于一定的临界值, 表明植被覆盖度减少^[29]。

(4)

式中, b表示植被覆盖度对年份的回归系数；n为监测年数；本文n为7；n_i表示年份；f_i表示第i年的植被覆盖度；n表示总年份数的平均值；f表示监测时段植被覆盖度的平均值。

2.2.4 植被覆盖度稳定性特征

变异系数(coefficient of variation)可以量化一组数据的离散程度^[30], 并且可以消除测量尺度和量纲的影响。本研究采用变异系数对研究区植被覆盖度的稳定性进行研究。

(5)

式中, CV表示植被覆盖度的变异系数；S表示植被覆盖度的标准差；x表示植被覆盖度的平均值；x_i表示第i年的植被覆盖度；n表示年数。

2.2.5 相关性分析

相关分析可以研究现象之间是否存在某种依存关系, 并对具有依存关系的现象探讨其相关方向以及相关程度。本论文采用Pearson相关系数(公式6)计算植被覆盖度变化与各因素的依存关系。

(6)

式中, r表示相关系数；n表示样本数；x_i表示样本中第i个因素的值；x表示样本中某个因素值的平均值；y_i表示样本中第i个样本植被覆盖度的值；y表示样本植被覆盖度的平均值。

3 结果与分析 3.1 植被覆盖度空间分布

采用基于地形调节植被指数的像元二分模型分别计算了3个矿区2001—2016年各像元植被覆盖度的平均值, 为了反映出不同研究区16 a间各像元植被覆盖度之间存在的差异, 依据李恒凯等^[31]对植被覆盖度的分级标准, 将研究区植被覆盖度分为5个等级(图 2)。

长河矿区16 a间植被覆盖度总体呈现中间低, 东西两侧高的格局, 多年均值分布由低到高分别占整个研究区面积比例的21.05%、21.97%、22.45%、17.72%, 16.81%, 中等以下等级占到65.47%。其中, 低植被覆盖度区域主要分布在长河及其两岸高程比较小的区域。离柳矿区16 a间植被覆盖度总体呈现北部低, 南部高的格局, 多年均值分布由低到高分别占整个研究区面积比例的13.01%、24.00%、31.13%、20.08%、11.78%, 中等以下等级占到68.14%。平朔露天煤矿区16 a间的植被覆盖度总体呈现中西部低, 北部和东南部高的格局(图 2), 多年均值分布由低到高分别占整个研究区面积比例的16.77%、18.69%、28.00%、22.49%、14.05%, 中等以下等级占到63.46%。3个矿区的植被覆盖度总体上偏低, 植被覆盖度大于0.6的面积比例仅占到35%左右。

提取每个研究区各像元的高程值和16 a的植被覆盖度均值, 通过对各研究区高程数据建立频次分布表进行分组, 求出每组内像元的植被覆盖度均值和对应的高程均值, 并建立二者的关系(图 3), 可以看出：长河矿区植被覆盖度随着高程的增加呈显著增加趋势(R²=0.78, P < 0.01), 离柳矿区植被覆盖度随着高程的增加也呈显著增加趋势(R²=0.89, P < 0.01), 平朔矿区植被覆盖度随着高程的增加趋势不显著(P>0.05)。

植被覆盖度的时间变化特征通过各像元16 a间对应年份的植被覆盖度与年份进行回归斜率分析, 得到不同像元上植被覆盖度16 a间植被覆盖度增减的空间趋势变化特征, 并按各像元回归得到的斜率值, 将植被覆盖度变化趋势分为严重退化、轻微退化、变化不明显、轻微改善及明显改善5类(图 4、表 2)。

3个典型矿区不同像元植被覆盖度16 a间的趋势变化特征如图 4所示。3个矿区中植被覆盖度均呈现出变化不明显的区域, 面积分别占到66.63%、59.90%和62.25%, 离柳矿区植被覆盖度变化不明显的像元所占比例较小；3个矿区均有植被覆盖度改善的区域出现, 所占比例分别为1.76%、2.00%、1.91%, 轻微改善所占比例分别为26.38%、30.55%、25.90%, 离柳矿区植被覆盖度改善的像元所占比重相对较大；同样3个矿区均有植被覆盖度退化的区域出现, 严重退化所占比例分别为0.11%、0.67%、1.54%, 轻微退化所占比例分别为5.12%、6.88%、8.40%, 其中平朔露天矿区植被覆盖度退化的像元所占比例最大。

3.2 植被覆盖度波动性

2001—2016年3个矿区植被覆盖度总体上均呈增加趋势(图 5), 长河矿区的增长速率为0.09%/10 a(R²=0.139, P=0.22)；离柳矿区的增长速率为0.10%/10 a(R²=0.296, P=0.12)；平朔矿区的增长速率为0.08%/10 a(R²=0.432, P=0.06)；3个矿区均没有通过0.05置信度的显著性检验。

把各像元计算得到变异系数分为低(CV≤0.3)、较低(0.3＜CV≤0.6)、中(0.6＜CV≤1)、较高(1＜CV≤1.4)、高(CV≥1.4)五个等级, 来反映各矿区变异程度的空间特征和稳定性程度(图 6)。3个矿区植被覆盖度变异程度的空间分布存在较大差异, 3个矿区不同程度有变异程度大的区域, 平朔露天矿区覆盖度变异程度较高的区域分布集中, 长河矿区与离柳矿区覆盖度变异高值区分布比较分散。3个矿区中变异程度级别处于中等所占比例分别为以下, 长河矿区占92%、离柳矿区占94%、平朔矿区占93%。

3.3 影响因素分析

随气候、地形特征的不同, 不同植被的响应会存在较大的差异^[32]。为此对3个矿区植被类型进行划分, 长河矿区和离柳矿区的植被类型主要分为草地、灌木林、旱地、疏林地、有林地5类, 平朔矿区的植被类型主要分为草地、旱地、灌木林、有林地4类。潜在的影响植被覆盖度变化的因素主要有降雨量、温度、高程、坡度、坡向及不同因素的互作效应等, 本研究从总的植被和区分不同植被类型两个角度分别讨论植被覆盖度变化的影响因素。考虑到实际取样中, 降雨量和温度数据在同一区域差别不大, 为此采用“以时间换空间”的方法, 得到不同降雨量和温度组合下的植被覆盖类别数据, 进而得到不同因素组合下的植被覆盖度数据。

本研究的气温和降雨量数据是基于影像获取时间的前3个月的积温和累积降雨量(图 7), 离柳矿区温度呈显著下降趋势(R²=0.53, P < 0.1), 长河和平朔矿区温度均呈下降趋势但不显著, 3个矿区的降雨量均呈不显著增加趋势。平朔矿区的积温较低, 长河矿区和离柳矿区的积温存在差异, 但差异不大, 但是均高于平朔矿区的积温；3个矿区累积降雨量有较大差异。采用相关分析分析了积温、累积降雨量、高程、坡度、坡向以及各因素间的互作效应对3个矿区16 a间植被覆盖度变化的影响因素(表 3)及不同植被覆盖类型植被覆盖度变化的影响因素, 见表 4。

从表 3可以看出, 长河矿区总的植被覆盖度与高程、坡度的相关系数分别为0.396、0.315(P < 0.01), 植被覆盖度与温度和坡度、高程的交互作用的相关系数分别为0.313、0.356(P < 0.01), 与降雨量和坡度、高程的交互作用的相关系数分别为0.253、0.098(P < 0.01), 与单因素的温度或降雨量没有表现出显著的相关性；离柳矿区总的植被覆盖度与温度、降雨量、坡度、高程的相关系数分别为-0.074、0.070、-0.190、0.138(P < 0.01), 与温度和降雨量的交互作用的相关系数为0.068, 植被覆盖度与温度和坡度、高程的交互作用的相关系数分别为0.186、0.108(P < 0.01), 植被覆盖度与降雨量和坡度、高程的交互作用的相关系数分别为0.160、0.088(P < 0.01)；平朔矿区总的植被覆盖度与高程的相关系数为-0.054、与高程和温度的互作效应的相关系数为-0.049(P < 0.01)。

对各个矿区的植被区分植被类型后, 对其进行影响因素相关性分析, 结果表明不同的植被类型对影响因素表现出不同的相关性, 具体相关性见表 4。对于草地, 在3个矿区中都表现出与降雨量、高程、温度与降雨量的互作效应、高程与温度的互作效应的相关性(P < 0.01), 相关性表现为长河矿区>离柳矿区>平朔矿区, 在长河和离柳矿区还表现出与坡度、坡度与温度、坡度与降雨量、高程与降雨量的相关性(P < 0.01), 相关性表现为长河矿区>离柳矿区, 在离柳矿区同时表现为与温度、坡向与降雨量的互作效应的相关性；对于灌木林, 在3个矿区均表现出与坡向、坡向与温度、坡度与降雨量、高程与降雨量的互作效应的相关性(P < 0.05), 在长河和平朔矿区表现出了与降雨量、温度与降雨量的互作效应较好的相关性, 在长河和离柳矿区表现出与高程、高程与温度的互作效应的良好相关性, 且长河矿区>离柳矿区, 在离柳和平朔矿区表现出与坡度、坡度与温度的互作效应的相关性；对于旱地, 在3个矿区表现出明显的差异, 在长河矿区表现出与高程、高程与温度的互作效应相关(P < 0.05), 在平朔矿区表现出与坡向、坡向与温度、坡向与降雨量互作效应相关。对于疏林地, 在长河矿区主要受温度、降雨量及其互作效应的影响, 而离柳矿区主要是受坡度和高程的影响。

4 讨论

煤矿开采区生态退化严峻, 为此从2006年起山西省开始矿区生态修复。本研究中3个矿区在近16 a间植被覆盖度均呈增加趋势, 主要是由于随着采矿的发展, 人类活动加强, 露天采煤区及新的排土场扩大, 导致植被覆盖度显著下降, 与此同时, 在生态修复的干预下, 植被恢复范围也逐渐扩大, 导致植被覆盖度明显改善, 表明生态修复初见成效, 这与前人的研究结果一致^[33-34]。对于植被覆盖度的影响因素人类活动的贡献显而易见, 对于自然因素的影响也不容忽视, 但多集中于单因子对总植被覆盖度的影响^[35-36], 实际是诸因子综合作用的结果, 本研究侧重从总的植被覆盖度角度和区分植被类型角度研究自然因素的影响程度, 研究结果表明, 对于井工矿区的长河和离柳区域受自然因素影响要多于露天开采的平朔矿区, 不区分植被类型时, 长河矿区和离柳矿区的植被覆盖度均与高程、坡度呈显著的相关性, 这与前人的研究结果一致^[37], 区分植被类型时, 不同的植被类型对不同的因子做出不同的响应。

长河矿区高程范围为710—1075 m, 平均为852.17 m, 当高程达到990 m以上植被覆盖度的增加趋势减小逐渐平稳, 离柳矿区的高程范围为639—1268 m, 平均为927.81 m, 高程达到1135 m植被覆盖度达到最大值, 大于1135 m植被覆盖度呈较少趋势, 总体植被覆盖度随海拔高度的增加呈先增加后降低的趋势, 这与张诗羽等^[18]的研究结果一致, 这主要是由于高海拔且坡度较大的区域人烟稀少、交通不便, 人类活动少, 植被类型相对稳定, 植被覆盖度相对较高的原因；平朔矿区的高程范围为1106—1564 m, 平均高程为1351.31 m, 但随着高程的增加植被覆盖度变化趋势不明显, 主要是由于相对于自然因素对植被覆盖度影响的缓慢过程, 受人类活动影响更直接。

长河矿区和离柳矿区主要是井工开采方式, 除了部分区域由于采煤引起地表裂缝与塌陷等严重破坏的情况之外, 多数耕地仍然可以继续从事农业生产, 整体上这些区域的植被覆盖度变化不明显。对于长河矿区, 植被覆盖度呈显著退化的区域主要位于海拔较高的成庄矿、天安普瑞矿区及周围的疏林地、部分高覆盖草地以及位于长河矿区东南部的疏林地, 本研究表明长河矿区疏林地与温度呈正相关性(P < 0.01)、与降雨量呈显著负相关(P < 0.05), 与高程与降雨量的交互作用呈显著负相关(P < 0.05), 长河矿区近16 a的温度呈下降趋势, 降雨量呈增加趋势, 这与本研究的疏林地影响因素分析结果一致, 植被覆盖度呈显著增加的区域主要分布于研究区的东北部的耕地, 属于泽州县大东沟乡的中东部和下村镇的南部区域, 这与该乡镇的发展高效生态农业有密切关系；对于离柳矿区, 植被覆盖度呈显著退化的区域主要位于平均高程为669 m的低海拔地形上, 分析可知低海拔区域主要位于离柳矿区北部的临县啧口镇与中部的柳林县穆村镇周边及其河谷地带、城镇周边居民地、建设用地, 农业用地相对强度要大, 受人为干扰是植被覆盖度显著降低的直接原因, 植被覆盖度改善区主要位于柳林县北部的孟门镇、成家庄镇和王家沟乡, 这与政府重视生态环境建设密切相关；对于平朔矿区植被覆盖度显著退化的区域主要位于安家岭露天采矿区、森泰煤矿区和东露天采矿区；植被覆盖度显著增加的区域主要位于安太堡露天采矿生态恢复区。

研究结果表明, 研究植被覆盖度变化不仅要从总的植被覆盖度变化角度研究, 而且要从不同植被类型的植被覆盖度变化角度进行研究。同时, 不同类型植被的植被覆盖度变化, 如退化、稳定和改善趋势对自然因素的响应是急需解决的科学问题。另外, 本研究只关注了自然因素对植被覆盖度变化的影响, 人为因素作为矿区一个重要的影响因素进行定量化研究值得进一步研究。如何定量区分气候变化与人类活动对矿区不同植被变化的贡献率仍将是植被变化对气候响应研究的重要内容。

5 结论

从时间和空间趋势分析, 16 a间3个矿区的植被覆盖度总体呈缓慢增加趋势(P < 0.05)；在重视生态环境建设的矿区复垦区域植被覆盖度得到改善, 煤矿开采用地的植被覆盖度明显降低, 露天开采的平朔矿区植被覆盖度退化严重, 井工开采的离柳和长河矿区的植被覆盖度变化相对平缓。从影响因素分析, 不区分植被覆盖度, 长河矿区表现出与坡度、高程、坡度与温度、高程与温度互作效应显著相关性(P < 0.01)；离柳矿区主要表现出与温度、降雨量、坡度、坡向及其他们之间的互作效应的影响(P < 0.01)；平朔矿区主要表现出与高程、高程与温度互作效应的相关性(P < 0.01), 长河和离柳矿区较平朔矿区受自然因素影响突出；区分植被类型, 同一植被类型在不同矿区与影响因素表现出不同的相关性, 同一矿区不同的植被类型对影响因素也表现出不同的相关性。