神东矿区植被覆盖度时序变化与驱动因素分析及引导恢复策略

刘 英，雷少刚，陈孝杨，陈 敏，杨英明，李心慧，张旭阳，龙林丽，卞正富

(1.安徽理工大学 地球与环境学院，安徽 淮南 232001；2.中国矿业大学 矿山生态恢复教育部工程研究中心，江苏 徐州 221116；3.深部煤矿采动响应与灾害防控国家重点实验室，安徽 淮南 232001；4.煤炭水资源保护与利用国家重点实验室，北京 102209)

我国重点建设的14个大型煤炭基地主要集中在生态环境脆弱、水土流失严重的黄土高原以及新疆等地区[1]。根据煤炭工业“十三五”发展规划，预计截止到2020年煤炭占我国一次能源消费中的比例仍占近60%，立足于国内是我国能源安全战略要求，煤炭作为能源主体地位和压舱石作用在很长一段时间内不会变化[2]。由于地理位置与气候原因，西部干旱半干旱地区水资源占有量仅占我国的8.3%，水资源十分匮乏，生态环境十分脆弱[3-4]。加之国家煤炭开采战略重心西移及集群化、高强度的开采方式，使得西部矿区的地表生态环境呈现出强烈扰动的态势，使本就脆弱的生态环境日趋退化[5]。西部煤炭资源开发引发的生态环境约束依旧是社会各界关注焦点[6]，大规模的煤炭开采引起地表沉陷，影响塌陷区地形、水文、土壤等植物生境因子，从而导致植物生境改变、影响植被生长[7-8]。

当前我国西部矿区已开展了诸多且富有成效的植被恢复重建研究，形成了以表土改良、地表裂缝治理、植物根际环境改善以及菌根共生等为主的植被恢复与重建技术[9-11]。但是这些技术成本往往过高，大面积种植新进的乔、灌类植物对本地植物物种形成冲击，大量开挖鱼磷坑亦会引起采后生态系统重复扰动，且部分重建植物难以适应恶劣的气候地理条件，管护周期长，成活率低，生长缓慢甚至衰亡，现有西部矿区重建植被可持续性不容乐观，其生态恢复效益以及对生态环境影响还有待商榷，形成了植被恢复资金高投入与恢复低效率的矛盾局面[12]。另一方面，煤炭开采地表塌陷对植被的影响并非全是负面的，部分区域的植被覆盖度、生物多样性呈增长态势[13]。所以当务之急，应严格遵循当地自然生态系统发展规律，避免对采后生态系统的再次扰动，科学配置，优化布局、引导恢复生态系统的自修复能力[14]。弄清煤炭开采扰动区植被覆盖度时序变化，是探索半干旱矿区植被引导型恢复有效方法重要基础与前提。

对于半干旱矿区受损植被引导型恢复来说，需要面临：植被在哪种破坏程度下可以实现自恢复?当需要人工引导干预时，在什么地方干预?怎么干预?干预到何种程度?解决上述基本问题，首先需要对矿区植被覆盖度时序变化及驱动因素进行分析，这是采煤扰动矿区植被受胁迫以及采后植被恢复状况空间感知、采后恢复区植被覆盖度变化趋势和年际波动程度评估的基础，可以为矿区有序、精准进行植被引导恢复工作提供依据。笔者将基于Landsat-NDVI时序数据，以神东中心矿区为研究尺度，分析采区和非采区植被覆盖度的时序演变规律，并依据不同时序植被覆盖度格点回归斜率、相关系数及标准差，对研究区不同采后恢复区植被覆盖度变化趋势以及年际波动程度进行评估，摸清研究区植被生长状况，为矿区有序、有针对地进行植被恢复工作提供基础数据。

1 研究区域概况与数据来源

1.1 研究区域概况

神东煤田中心矿区位于鄂尔多斯高原东南部及陕北黄土高原北缘和毛乌素沙漠的东南边缘(图1)，地理坐标为东经109.83°～110.34°，北纬39.56°～39.19°，面积约900 km2。该区年平均气温6.2 ℃，极端最低气温-31.4 ℃，极端最高气温36.6 ℃。年降水量为300～400 mm，主要集中在夏季，蒸发量大于降雨量4倍以上，属于典型的干旱、半干旱的高原大陆性气候。地处干草原与森林草原的过渡地带，主要植被类型为干草原、落叶阔叶灌丛和沙生类型植被，主要植物品种有沙柳、杨树、长芒草、芨芨草、油松、柠条、沙棘、紫花苜蓿、沙竹、油蒿等。地貌地形西北高，东南低，平均海拔在+1 200 m左右。矿区东部及东北部为黄土丘陵山区，区内沟壑纵横，形成梁峁、沟壑和土塬3种地貌，侵蚀活跃，水土流失严重；西部及西南部以风沙地为主的流动、半固定及固定沙地上，分布着沙地植被，主要是沙地先锋植物群落和油蒿群落；在洼地、滩地和湖泊周围分布有湿地植被。

图1 研究区域位置Fig.1 Study area geographic location

1.2 数据来源

选取2000-01—2016-12月的共337景Landsat 5 TM和8 OLI遥感影像数据，时间分辨率为16 d，空间分辨率为30 m，影像下载列/行号为127/33。所有Landsat影像均从美国地质调查局(USGS)网站(https://glovis.usgs.gov/)下载。由于遥感数据不可避免的而受到云层覆盖的影响，数据选取时剔除云量大于5%的Landsat影像。为了提高精度和减少地形和大气的影响，利用软件ENVI 5.1(Exelis Visual Information Solutions USA)对遥感影像进行预处理。所有获取的Landsat TM和OLI遥感影像数据经辐射定标、大气校正、配准、裁剪等预处理，并根据研究区域的矢量边界进行裁剪。所有影像统一采用UTM投影，WGS84坐标系。所有影像以2000年1月份的影像为基准配准，以保证影像之间误差在1个像元以内。

1.3 矿区开采工艺及开采情况

该矿区煤炭资源赋存条件好，目前开采的矿井主要包括大柳塔矿、补连塔矿、哈拉沟矿、上湾矿等，以大采高一次采全高综合机械化采煤为主要采煤方法，300 m宽工作面较为普遍，走向长度1 500～5 220 m不等，采高通常在2.0～7.0 m，由于综合机械化程度高，平均采速每日可推进10.0 m以上，有的甚至达到20 m/d，地表最大下沉值在2 320～5 887 mm，最大水平位移可达3 284 mm，最大下沉速度可达837 mm/d，地表综合边界角、超前影响角与综合边界角一般分布在50°～70°，综合移动角一般为60°～70°，裂缝角一般为65°～85°，研究区部分工作面属于初采，大量工作面属于重复采动。煤炭资源大规模、高强度开采引起的地表沉陷具有下沉速度快、下沉量集中，造成地表破损严重，地裂缝、塌陷坑等非连续形变发育等特点。

2 研究方法

2.1 植被指数与覆盖度的计算方法

归一化植被指数(Normalized Difference Vegetation Index,NDVI)可描述植被的生长状况，是定量评价植被盖度常用指标之一[15]。其计算时由于植物叶片组织对绿光，尤其是近红外光有很强的反射，对蓝光(470 nm)和红光(650 nm)有很强的吸收。因此，在Landsat影像中提取植被信息时多采用红外与红波段图像反射值比值，这样可消除部分辐射误差。其计算公式为

(1)

式中，INDV为归一化植被指数；NNIR,NR分别为近红外和红光带的实际表面反射率。

植被覆盖度(Vegetation Coverage Fraction，VCF)指某一地域植物垂直投影面积与该地域面积的比值。目前应用最广泛的是利用归一化植被指数(如INDV)，根据像元二分模型计算得到植被覆盖度指数[16]，其计算公式为

(2)

(3)

(4)

2.2 植被覆盖度时序变化趋势计算方法

依据神东矿区2000—2016年时序植被覆盖度数据，以格网中包含像素植被覆盖度的平均值，通过对每一个像元进行一元线性回归方程拟合，得到每个像元回归直线的斜率θslope，当θslope<0时，表示植被覆盖度呈减少趋势；当θslope>0，表示植被覆盖度呈升高趋势。采用拟合相关系数R来判定变化趋势显著性，当R<0时，表示植被覆盖度呈减少的趋势；当R>0，表示植被覆盖度呈升高的趋势[18]。用t分布检验2个变量是否真正相关及其相关性的显著水平[19]。计算公式为

(5)

(6)

其中，θslope为一元线性回归方程斜率；ti为年份序号，i=1，2，…，15；xi为第i年生长季植被覆盖度均值；R为相关系数；n为总年份数，在本研究中n=15。依据式(5)计算得到每个像元回归直线的斜率θslope，根据t检验结果对应的拟合相关系数R，将变化趋势及显著性分成10个等级，即：①θslope<0，0.010，0.10，0.050，0.0250，0.010，0.01

2.3 植被覆盖度时序波动程度计算方法

利用时序变异系数(Standard Deviation，SD)来表示植被时序波动程度，时序变异系数值越大，说明植被覆盖度年际波动越大，受外界干扰程度越严重；反之，则植被覆盖度年际波动越小，受干扰程度低。其计算公式为

(7)

其中，DS为植被时序变异系数。结合式(5)得到的θslope，将植被覆盖度波动程度按正负分成10个等级，即：①θslope<0，DS≥0.4，划分为“极高幅度降低”；②θslope<0，DS∈[0.30,0.40)，划分为“高幅度降低”；③θslope<0，DS∈[0.15,0.30)，划分为“中幅度降低”；④θslope<0，DS∈[0.10,0.15)，划分为“低幅度降低”；⑤θslope<0，DS∈[0,0.10)，划分为“极低幅度降低”；⑥θslope>0，DS∈[0,0.10)，划分为“极低幅度升高”；⑦θslope>0，DS∈[0.10,0.15)，划分为“低幅度升高”；⑧θslope>0，DS∈[0.15,0.30)，划分为“中幅度升高”；⑨θslope>0，DS∈[0.30,0.40)，划分为“高幅度升高”；⑩θslope>0，DS≥0.4，划分为“极高幅度升高”。

3 矿区植被覆盖度变化及驱动因素分析

3.1 基于Landsat的神东中心矿区植被NDVI时序变化分析

由于各工作面开采时间不同，主要选取各矿井2014,2011,2009,2007,2004和2001年开采的典型工作面区域为采区，以神东中心矿区为对照区，对2000—2016年间植被NDVI时序变化进行统计(图2)，图中红色箭头为不同煤炭开采年份。从整体时序变化来看，2000—2016年神东中心矿区植被INDV整体呈现出明显物候性周期变化。植被INDV的最大值出现在每年7月下旬至8月中旬，最小值出现在每年12月份到1月份。从植被INDV年际变化来看，2000—2005年植被INDV变化幅度变化较小，总体在0.28左右；2005—2008年植被INDV小幅上升，然后2008—2010植被INDV呈下降趋势，相比于2008年峰值降低了近19.15%；2010—2014年植被INDV出现明显上升，相较于2010年峰值增加了近31%，且以2014年的峰值最大；2014—2016年植被INDV有小幅降低，整体在0.40左右波动。但是从整个时序区间来看，神东中心矿区植被INDV呈明显的增大趋势，这说明植被覆盖情况总体上有所好转。原因在于：地方政府和神华集团进行了大面积矿区植被重建，完成水土保持以及绿化、大面积防风固沙和多项小流域水土保持综合治理措施，对整个矿区植被覆盖度有明显提升，由图2可以看出，截止到2016年，神东中心矿区植被INDV从2000年0.20左右提高到0.45～0.60。

图2 2000—2016年神东中心矿区与不同开采年份工作面INDV(30 m)时序变化Fig.2 Time-series changes of INDV (30 m) of Shendong central mining area and working face in different mining years from 2000 to 2016

由于研究区不同开采年限不同，由此引起的植被INDV变化具有差异性，为此将不同开采年限工作面时序INDV求平均值，减去中心矿区时序INDV，得到不同开采年限工作面INDV差异值时序变化情况。统计发现，INDV差异值时序变化与该工作面煤炭开采年份密切相关。煤炭开采前工作面地表INDV约等于甚至大于中心矿区INDV，INDV差异值均在0附近上下波动，变化幅度为±0.03；时序INDV差异值大于0，说明该工作面采前INDV高于作为对照的中心矿区INDV，这主要由于这些被统计的工作面主要分布在大柳塔矿井硬梁地覆盖区，特别是5-2煤层开采区，植被长势相对较好，分布着大面积的杨树、柠条、沙柳、油蒿、长芒草、百里香、沙棘等，而中心矿区还包括大面积风沙覆盖区，这些地区植被多以低矮沙生草本、灌木为主，其植被覆盖度相对较低。随着煤炭开采活动的推进，INDV差异值逐渐增大，且小于0，说明该煤炭开采工作面INDV小于中心矿区INDV，从工作面采前到采后INDV差异值从“大于0”或者“0附近波动”转变“小于0”，对应的采区非采区INDV表现为“采区INDV不小于等于非采区”转变为“采区INDV小于非采区”，这体现了随着采煤扰动影响地表植被INDV相对于对照区降低的变化趋势。通过对采后INDV差异值时序变化趋势进行拟合发现：采后5 a，工作面地表INDV差异值呈增大趋势；采后7 a，工作面地表INDV差异值开始降低；采后12 a工作面地表INDV差异值逐渐趋向于0。说明采后5 a内，相对于对照区，工作面采区植被INDV的变化表现为持续降低的过程；采后7 a，受扰动区植被开始恢复，工作面地表INDV与对照区INDV差异值开始降低；至采后12 a，工作面采区植被INDV基本能够恢复至对照区水平。

3.2 神东中心矿区植被覆盖度时序变化趋势分析

采用式(5)，(6)计算得到不同格点植被覆盖度一元线性回归斜率θslope(图3(a))和对应的相关系数R(图3(b))，利用ArcGIS叠置分析功能对计算结果进行叠加得到神东矿区植被覆盖度空间变化趋势图(图3(c))。

对图3结果进行统计得到：随着时间的推移，神东中心矿区植被覆盖度均值整体呈现小幅增加趋势，植被覆被度呈升高与降低的区域面积分别占中心矿区总面积的72.35%和27.65%，其中θslope>0.01的面积占研究区总面积的52.03%，研究区域植被覆盖度以正向变化为主，由此可见，煤炭开采引起的地表塌陷对植被覆盖度的影响并不全是负面的。θslope<0，且R∈[-1,-0.89)内，通过了0.01显著性检验，植被受损最为严重，植被覆盖度极显著性降低，面积约为4.81 km2，占研究区域面积的0.62%。θslope<0，且R在[-0.89,-0.71)内，通过了0.025显著性检验，植被受损较为严重，植被覆盖度较极显著降低，其面积约为9.39 km2，占整个研究区域面积的1.21%。θslope<0，且R在[-0.71,-0.58)内，通过了0.025显著性检验，植被覆盖度显著降低，其面积约为20.89 km2，占整个研究区域面积的2.69%。θslope<0，且R在[-0.58,-0.32)内，通过了0.05显著性检验，植被覆盖度较显著降低，其面积约为108.26 km2，占整个研究区域面积的13.92%。值得注意的是，还有大约占整个中心矿区面积的20.78%，161.69 km2的区域地表植被覆盖度变化θslope在[-0.01,0.01]内，尽管整体上神东中心矿区植被覆盖度呈上升趋势，但是在煤炭开采活动最为集中、扰动最为频繁的中心矿区植被覆盖度升高幅度较低，甚至有大面积的植被退化区域。上述植被覆盖度降低的区域主要在神东中心矿区各矿井工作面上覆地表小区域成片分布。因此，在矿山植被修复重建时，有针对性的对植被退化区域进行识别，然后对识别出的退化区域植被进行引导修复，可以降低矿山植被修复的成本，提高植被重建的效率，这也是引导型恢复提出的立足点。

图3 神东中心矿区植被覆盖度时序变化趋势Fig.3 Time-series variation trend of vegetation coverage in Shendong central mining area

3.3 神东中心矿区植被覆盖度时序波动程度分析

采用式(7)计算得到不同格点植被覆盖度时序波动程度(图4(a))，并利用ArcGIS叠置分析功能对计算结果与一元线性回归斜率θslope的正负情况(图3(a))进行叠加，得到神东矿区植被覆盖度正向、负向波动情况(图4(b))。

中心矿区年际间植被覆盖度变异系数(DS)处于0.04～0.46 (图4(a))，根据划分的植被覆盖度波动等级，极低幅度、低幅度、中幅度、高幅度和极高幅度波动区的面积分别占中心矿区总面积的1.11%，34.71%，56.71%，6.09%和1.38%，以中、低幅度变化为主。结合θslope的正负，研究区域东北部、中部、南部大面积地区植被波动等级为低幅度降低，该地区为荒漠化严重地区，植被覆盖度本来就很低，加之煤炭井工开采地下水会被抽空，影响地下水对地表水的补给，植被生长的立地条件会受到极大程度的破坏，进而影响地表受扰动区植物的生长与恢复，导致该区域植被生长与恢复状况较差。中心矿区东南部、西南部、西北部年际间植被覆盖度变化主要以中幅度变化为主，结合区域θslope的正负，植被覆盖度变化表现为中幅度升高，这与神东矿区大力开展植被恢复活动，先后建立完善了专门水保绿化机构，制定了矿区防风固沙与水土保持规划，以林草植被建设为主，完成矿区植被重建以及大面积常规防风固沙和多相小流域水土保持综合治理工程，致力于提高该地区植被覆盖度与多样性有关[3]，如，在大柳塔、补连塔、哈拉沟等采煤塌陷区建立了植被修复重建应用的国家级示范区。此外，居民区、城市道路与河道两边以及煤矿周边辅助设施建设区域，植被受人为活动影响较大，植被覆盖度变化幅度较大。

图4 神东中心矿区植被覆盖度时序波动程度Fig.4 Time series fluctuation degree of vegetation coverage in Shendong central mining area,coefficient of variation

3.4 植被覆盖度变化驱动因素分析

3.4.1植物覆盖度与土壤理化性质关系分析

以大柳塔矿5-2和2-2部分工作面作为实验样地，选取4个采煤工作面，共建立18个植被调查样方，在每个样方内随机建立10 m×10 m植被群落学调查大样方，并利用5点取样法在样地区域内设置5个2 m×2 m小样方，并在每个小样方利用对角线法设立1个1 m×1 m基本调查单元，同时也作为土壤理化性质数据取样的网格单元，然后利用Canoco 4.5软件大柳塔矿井18个样方植被与土壤理化性质典范对应关系进行分析。主要考虑各样方中与植被生长密切相关的土壤含水量(SWC)、有机质(OM)、容重(BD)、pH值、全氮(TN)、总磷(TP)、全钾(TK)、总有机碳量(TOC含量)、土壤电导率(EC)及各样方中24种主要建群种植物、个体数量与盖度，分析结果如图5所示，图5中红色、蓝色、青色分别为土壤理化性质、植物名称和样方编号，红色剪头线段的长短代表各土壤理化性质因子对不同植物及样方的影响大小，红色与蓝色箭头线段之间的夹角代表相关性大小，锐角表示呈正相关，钝角则表示呈负相关，夹角越小相关性越大。

图5 植被-土壤理化性质典范对应关系Fig.5 Canonical correspondence between vegetation and soil properties

从植被与土壤理化性质典范对应关系分析结果来看，第3象限上的因子可以解释67.2%的植被分布状况。土壤含水量SWC与第1排序轴(横轴)第3象限相关系数达到0.727 2，与第2排序轴(横轴)第3象限相关系数达到0.746 7，说明SWC为影响植被生长的第1重要土壤理化性质因子。从植被样方在排序图上的分布，可以看出，有12个样方主要分布在第3象限，并且24种主要植物建群种有13种分布在第3象限，如：紫花苜蓿、杨树、柠条、华北白前、沙棘、斜茎黄耆、狗娃花、胡枝子、地角儿苗、芨芨草等。此外，土壤有机质(OM)与第1排序轴相关系数为0.649 4，表明其也是影响植被分布较为重要的因子，其他重要影响因子还包括容重(BD)、全氮(TN)等。诸多研究表明采煤沉陷裂缝破坏了土体结构，增加了土壤水分的蒸发面，是裂缝区土壤水分降低的主要原因[12,20-21]。

3.4.2植物覆盖度与地下水位埋深关系分析

植被生长过程中，水资源扮演着举足轻重的角色，而对于生长在半干旱区的植被而言，区域内降水量低，因此地下水位埋深对其生长状态起着至关重要的作用，是控制半干旱矿区地表植被生长状况的又一关键决定因素[7]。有诸多学者对我国半干旱地区植物覆盖度与地下水位埋深关系进行了研究。范立民团队较早关注到神东矿区煤炭资源大规模开采与地下水资源破坏这一矛盾，指出该矿区煤层浅埋区采煤会引起大面积地下水漏失、减少河流基流量，从而破坏植被生态系统[22]。徐海量等[23]对我国西部干旱地区不同地下水位对植被的影响研究后认为3.5 m为该地区草本植被生态恢复的最低水位阈值。李丹[24]利用Gamma分布函数对位于鄂尔多斯高原中部的海流兔河流域生长季多年平均INDV随地下水埋深的变化，发现地下水埋深在小于10.0 m时，INDV随地下水埋深增大而减小；地下水埋深大于10.0 m时，INDV与地下水埋深基本没有关系。金晓媚等[25]对海流兔河流域植被分布进行研究时认为，地下水位埋深对植被的影响范围是1.0～5.0 m，当地下水位大于5.0 m 时，决定植被的主要因素就变成了降雨量和土壤含水量。由此可见，不同研究区地下水埋深对植被的影响阈值不同，但是地下水埋深对植被的影响阈值总体处于0～10.0 m。

对于神东中心矿区，根据卞正富[12]、雷少刚[3]等对神东矿区地下水埋深与植物群落关系长期现场调查结果，认为地下水埋深3.2 m和8.0 m为影响植被INDV的2个重要阈值(图6)。

图6 地下水位埋深与植被群落类型的关系 (依文献[26]补充)Fig.6 Relationship between groundwater depth and vegetation community types (supplemented by reference[26])

地下水埋深在0～3.2 m内，对根系较浅的湿生植被分布影响最大，随着地下水位埋深的增加，与湿生植被INDV相关性降低，此时土壤水与地下水之间作用强烈，土壤含水量饱和与非饱和状态交替频繁，湿生植被的生长与地下水的变化密切相关；当地下水埋深在3.2～8.0 m时，这些根系较浅的湿生植被与地下水位埋深之间不再有显著的联系，取而代之的是一些根系较长的旱生植被，此时地下水和土壤水之间双向联系，地下水仍可以通过蒸发与毛细管作用对土壤水进行补给；当地下水埋深大于8.0 m时，植被生长主要依靠的是少量存储在土壤中的大气降水和大气凝结水，而与地下水埋深之间的关系不再显著，地表植被覆盖格局将进一步由旱生植被将演替为沙生植被[3,26]。

3.4.3植被NDVI时序变化与降雨量、温度关系分析

气候条件(降雨量、温度)是决定地表植被生长和分布的重要因素[27]。利用位于矿区的大柳塔气象站收集到的研究区域2001—2016年的气象数据，分析月平均温度和降雨量与神东中心矿区月植被INDV的相互关系(图7)。神东中心矿区植被INDV与月平均气温和降雨量具有相似的年周期变化，如图7所示，这也使得植被INDV与月平均气温(R2=0.640 1，P<0.05)和降雨量(R2=0.651 1，P<0.05)具有较高的相关性。从2000—2016年，神东中心矿区降雨量和INDV均呈增加趋势(线性拟合斜率分别为0.004 8,0.000 38)，而温度相对较稳定(线性拟合斜率为-0.009 2)，且INDV和平均降雨量在每年6—8月份会有较大幅度的升高，温度在此期间平均气温却无明显差异，根据图7，2007—2016年间中心矿区植被INDV的升高、INDV差异值的变化与降雨量的增加也极具同步性，说明降雨量是植被INDV波峰值的最敏感气候因子，对中心矿区植被覆盖度的变化有极其重要的影响。

图7 神东中心矿区植被INDV与月平均温度、降雨量的变化Fig.7 Changes of INDV and monthly mean temperature and rainfall in Shendong mining area

3.4.4采矿与土地损毁植被覆盖的影响传递过程分析

(1)采空区上方岩层破断与运动过程。井工煤炭开采在地底形成大面积采空区，破坏岩层应力的平衡状态，引起岩体内部应力重新分布，导致岩层破断和运动，但是不同类型岩层的破坏变形又具有差异性、非同步的特点。随着工作面推进软岩层破断和运动形式表现为同步断裂，而对于硬岩层，其岩层破断和运动形式表现为间断式断裂。硬岩层破断移动后形成的块体重新分布，势必会改变采空区整个上覆岩层初始状态，因此，硬岩层(关键层)破断形成的块体主导着岩层采空区上覆岩层运动。随着采空区硬岩层破断与关键层破断为岩块，则上覆岩层运动最终形成如图8所示的关键层主导的采动覆岩运动概貌，通常在采空区上方形成垮落带、断裂带和弯曲下沉带。当关键层岩块破断后岩层运动重新趋于平稳时，岩层裂隙将重新闭合，形成重新压实区；对于岩层发生剧烈运动，矿压、裂隙与塌陷剧烈变化区域称为裂缝发育区[28]。

图8 采空区上覆岩层运动分区[28]Fig.8 Movement zoning of overlying strata in goaf[28]

(2)岩层破断与运动到土体沉降变形过程。随着岩层破断和运动过程之后是地表土体沉降变形过程，通常在采空区上方地表形成大范围的塌陷盆地。根据开采沉陷学理论，塌陷盆地中不同区域的下沉程度不同，一般将下沉盆地划分为拉伸变形区、压缩变形区和中性区[29]。“三区”的形成原因以及对土体结构的影响有所差异，对于下沉盆地中间区被称为中性区，该区地表下沉值最大，一般先会出现临时性裂缝，煤炭开采过程中随工作面的推进同时发育，当工作面推过裂缝后，该区破断后岩层运动重新趋于平稳，大部分裂缝将逐步闭合，其相对移动和变形值较小；塌陷盆地内边缘区被称为压缩变形区，位于采空区边界附近到最大沉降点之间，地表土体移动向盆地中心倾斜，使得土体产生压缩变形而形成挤压型裂缝；塌陷盆地外缘区被称为拉伸变形区，位于采空区边界到盆地边界之间，地表沉降程度不均匀使得地表土体移动向盆地中心方向倾斜，土体产生拉伸变形，当拉伸变形超过极限抗拉强度后，地表将产生拉伸型裂缝，其特点为发育宽度、深度大、难以自恢复，破坏原始土地结构[5]，对表土水分的保持以及植物的生长的影响较大，是塌陷区植被恢复需要重点采取工程措施加以整治的区域。

(3)土体沉降变形到植物生长立地条件变化过程。随着采空区土体塌陷变形过程后，地表植物生境条件随之发生改变，主要表现为地形地貌条件、地下水位埋深、SWC、土壤理化性质等的改变[5,12]。地下煤炭开采对SWC的影响在不同塌陷区位存在空间异质性。中性区位于地表下沉盆地的中心，沿开采工作面倾向上其地形往往呈“两边高，中间低”状，该地形有利于地表径流的汇集，导致中性区SWC高于压缩区、拉伸区。压缩区土体被挤压，土壤孔隙度降低、硬度升高、入渗能力降低，不利于土壤水分保持，所以相较于中性区，压缩区的SWC通常更低。拉伸区地表存在永久性裂缝，在一定程度上改变地表径流方向和汇水条件，当地表裂缝与采空区直接连通时，地表径流将直接流进采空区(图9)。此外，永久性裂缝还会导致土壤水分的蒸发面增加，土壤水的蒸发散失速度加快，而裂缝周边的土体受拉伸作用导致结构破坏，也会增大土壤的孔隙比，不利于土壤水分保持，地表水分流失进一步加重[20]。

图9 采煤塌陷对植被生长立地条件变化过程示意Fig.9 Change process diagram of plant growth site conditions caused by mining

通过沉陷过程不同变形区位土壤含水量进行时序监测后发现：土体塌陷变形对“三区”土壤含水量的影响总体表现为，中性区SWC>压缩区SWC>拉伸区SWC(图10)。采后1个月中心区、压缩区、拉伸区土壤水开始出现了较为明显的降低，尤其是永久性裂缝大量发育的拉伸区土壤水下降最为明显，从采前的11.8%最低降至5.5%；压缩区土壤含水量从采前的11.9%最低将至8.4%；中性区随着临时性裂缝产生与自恢复，土壤含水量呈波动状，在8.9%～12.1%波动。SWC在煤炭开采地表不同沉陷区位的空间异质性也是沉陷区植被覆盖度时空变化特征不一致的主要原因。

图10 开采前后不同塌陷区位土壤水分时序变化Fig.10 Temporal variation of SWC at different subsidence locations before and after mining

4 主要结论与建议

4.1 主要结论

(1)从2000—2016年神东中心矿区植被NDVI时序上整体呈物候性周期变化，2005年后中心矿区植被INDV明显提升，在2014年的INDV峰值最大。INDV差异值时序变化与该工作面煤炭开采年份密切相关。

(2)植被覆盖度呈升高与降低的区域面积分别占中心矿区总面积的72.35%和27.65%。中心矿区年际间植被覆盖度标准差处于0.04～0.46，植被覆盖度波动以中、低幅度变化为主。

(3)土壤含水量是影响植被生长的第1重要土壤理化性质因子，影响植被生长较为重要的因子包括有机质含量、容重和全氮含量等。地下水埋深3.2 m和8.0 m为影响神东矿区植被INDV的2个重要阈值，地下水埋深通过影响土壤含水量影响地表植被类型，当地下水埋深大于3.2 m后，根系较浅的湿生植被演替为旱生植被；当地下水埋深大于8.0 m后，旱生植被演替为沙生植被。降雨量是中心矿区植被覆盖度最敏感气候因子。

4.2 神东矿区植被引导型恢复建议

半干旱矿区植被引导型恢复首要任务需要对工作面不同沉陷区位进行划分，在此基础上重点对拉伸区地表永久裂缝进行识别和治理，这是采后沉陷区土壤水分条件恢复重要基础与前提。需要加强对沉陷区植物生长与土壤水分之间的耦合关系研究，分析植物生长受到干旱胁迫及植物死亡的相关土壤含水量阈值，可以为需要引导恢复时，土壤含水量提升到何种程度提供依据。

对于地下水位埋深在3.2 m以内的区域，植被恢复引导主要以降低植被对潜水的依赖为主，比如：种植抗旱性较强的本地植物、改善植被空间结构，但是这种工程类尝试还比较少，其可行性还有待进一步研究论证。对于地下水位大于3.2 m的区域，笔者建议可以采用高空间分辨率的遥感影像对研究区地表植被进行精细分类，结合本文得到的植被INDV时序变化趋势及波动程度，兼顾对植被生长立地条件的考虑(包括坡度、坡向、高程、土壤类型及其他理化性质)，以植被INDV呈正向波动区植被组成情况作为植被引导恢复时的植物配置、覆盖度控制的参考依据，且植物配置时尽量选择具有耐土壤贫瘠、抗旱性强、根系发达的本地植物物种。此外，应该加强对矿区地下水资源分布的探测，综合考虑地下水的采前埋深和采后埋深，对采矿引起地下水变化导致的植被影响评价，加强对采空区积水的循环利用，有计划地优先开采无、少地下水的区域，从源头减少或延缓对地下水的破坏。