采煤胁迫下干旱区植被生长周期红边指数动态变化
2021-06-18李晓婷杨丽帆邹友峰
李晓婷，杨丽帆，邹友峰,3，马 超,3
(1.河南理工大学 测绘与国土信息工程学院，河南 焦作 454003; 2.中国地质大学 地理与信息工程学院，湖北 武汉 430074; 3.河南理工大学 自然资源部矿山时空信息与生态修复重点实验室，河南 焦作 454003)
煤矿开采对塌陷区地表植被影响及恢复的程度、范围和时间，一直是学术界争论的焦点问题之一[1-2]。由于受迫植被生理生态指标的甄别需要大范围、长时间的野外调查和实验室工作，相关研究进展缓慢。自20世纪70年代初以来，遥感已经帮助我们极大地提高了对我们所居住的星球的认知[3]，基于遥感手段的很多研究已经提高了我们对植被响应外在条件变化的理解[4]。近10 a来，基于星载、机载和地面平台的多光谱、高光谱遥感技术的植被胁迫研究，一直是植物生态毒理学、病理学研究的热点，国内外关于重金属胁迫如铜和铅[5-7]、镉[8-9]、锌、砷、锰[10-11]等，水热条件胁迫如降水[12]、干旱[13-15]、高温[16]、冻害[17]等，以及农药胁迫[18]，交通胁迫[19]，盐碱胁迫[20]，病虫害胁迫[21-22]和酸性矿井水胁迫[23]导致的植被叶绿素、胡萝卜素的荧光特征、反射光谱特征、红边响应特征研究屡见报导。
开采沉陷对地表植被的直接损伤主要是根系损伤和土壤物理性质的变化，沉陷造成土壤表面裂缝生成，土壤持水能力减弱，进而形成断根胁迫和干旱胁迫[24];长期而言，更会造成土壤氮、磷、钾等养分元素损失，废矸废水废气排放引起植被健康状况恶化，导致光合作用能力变差，生产力下降[25-26]。开采沉陷造成植物生长层受到损伤，对植被的影响是持久的、长期的，植被作为陆地生态系统的主体，是联结土壤、大气和水分等要素的自然纽带，是矿区生态环境变化最直观的“指示器”。而归一化植被指数NDVI是公认的表征植被变化最有效的生态指标之一，基于连续一致的NDVI时间序列数据，沉陷区植被生长状况及其动态变化研究取得了众多成果[27-31]。但由于植被监测的多光谱遥感大多缺少红边波段，鲜有其他遥感生态指标应用于矿区植被生态胁迫的范畴，单一指标表征矿区生态走势尚显薄弱。进入21世纪，以哨兵2号(Sentinel-2)卫星为代表的一系列红边波段传感器研究取得突破性进展，作为植被生长变化的最敏感的波段，一系列创新性红边指数被提出[32-36]。有研究表明，利用Sentinel-2A红边波段估算植被叶绿素时，其决定系数R2最高达到了0.99[37];利用Sentinel-2A数据估算的大豆红边叶绿素指数(CIre)与其叶绿素浓度具有非常高的相关性(R2= 0.94)[38];使用Sentinel-2A多光谱卫星数据提取的红边叶绿素指数(CIre)、红边比值指数(SRre)和叶绿素敏感指数(MTCI)在冬小麦全生长季内对生物量的变化保持高灵敏性[39]。
目前，基于红边波段的高敏感指数的采矿胁迫遥感研究尚属空白，开采沉陷影响地表植被的机理有待深入研究。笔者以地处黄河流域中上游的鄂尔多斯市泊江海子矿为工程示范，针对113101 首采工作面地表覆被，开展一个生长周期的采矿胁迫实验研究。利用归一化植被指数(NDVI)和改进红边归一化植被指数(Red Edge Normalized Difference Vegetation Index,NDVIre)、归一化差值红边(Normalized Difference Red Edge,NDre)、修改的叶绿素反射率指数(Modified Chlorophyll Absorption in Reflectance Index,MCARI)、修改的红边简单比值指数(Modified Red Edge Simple Ratio,MSRre)和三角叶绿素指数(Triangular Chlorophyll Index,TCI)共5种高敏感红边指数，诊断煤炭地下开采导致生长层生态损伤的因果关系，定量分析黄河流域脆弱生态区受迫植被退化程度与恢复过程，旨在阐明采矿与复垦活动对自然生态的扰动与植被生态系统修复的时空关系，期望丰富煤炭开采生态损伤与生态环境协同保护评价指标体系，为黄河流域脆弱生态环境保护与生态重建提供一定的科学依据。
1 研究区与数据
1.1 研究区
泊江海子矿地处黄河流域中上游的黄土-风积沙区，库布齐沙漠以南，毛乌素沙地以北，资源同属于我国十四大煤炭基地的神东基地(图1)。行政隶属内蒙古鄂尔多斯市东胜区泊江海子镇。地理经度109.246°～109.417°，地理纬度39.838°～39.879°。井田东西长约20.16 km，南北宽约5.94 km，面积115.62 km2。113101工作面呈正南正北倾向长壁布置，于2016-12-01联合试采。工作面地表为荒漠化地区，有少量农田和陈家社、瓦窑社村。研究区气候特点为高原大陆性干旱-半干旱气候，日照时长较长，干燥，风沙大，降水少[40]。
图1 研究区地理地貌Fig.1 Geography and geomorphology in the study area
1.2 数 据
哨兵二号(Sentinel-2)A，B星同处一条轨道上，相差180°。其中Sentinel-2A于2015年6月发射。Sentinel-2A重访周期为10 d。
载荷多光谱成像仪(Multi-Spectral Imager,MSI)的13个波段覆盖了可见光-近红外到短波红外，扫描带宽度290 km，具有3种分辨率(10，20和60 m)，见表1。
表1 Sentinel-2A卫星多光谱波段信息Table 1 Multispectral band information of Sentinel-2A satellite
Sentinel-2A特有的3个红边波段(即B5，B6，B7)具有高空间分辨率、高光谱分辨率、宽幅宽、高重访频率特征，可为监测地表植被的生长状况提供更加有利的数据支持。可用于陆地、海洋、极地以及全球变化监测，还可应对洪涝、森林灾害和地震等突发情况[41]。
Sentinel-2A影像每月1次覆盖研究区域。主要选择天气晴朗且云量少的高质量影像。由于7月和9月云量覆盖较大，影像质量较差，因此未选用这2期数据。11期数据信息见表2。
表2 研究采用Sentinel-2A数据Table 2 Data of Sentinel-2A used in study
2 方 法
2.1 数据预处理
Sentinel-2A的Level-1C产品数据是经过辐射校正和几何精校正的大气表观反射率产品，没有进行大气校正。采用欧空局提供的Sentinel-2A/B大气校正工具Sen2cor(Sentinel-2 Level-2A Atmospheric Correction Processor)实施高精度大气校正[42]。Sen2cor可以将L1C产品处理成12个波段的L2A产品，同时将大气上层表观反射率转化为大气下层地表反射率(图2)。
图2 典型地物大气校正前后的光谱曲线Fig.2 Spectral curves of typical surface features before and after Atmospheric Correction
图2为典型地物(植被、水体、建筑物、裸地)大气校正后效果，蓝波段B2经大气校正后相比于大气上层表观反射率TOA显著下降，红边、近红外波段大幅增加，其光谱曲线与典型地物的光谱曲线的基本吻合。
大气校正完成后，因不同波段的空间分辨率不同，利用Sentinel Application Platform(SNAP)对L2A数据采用最邻近像元法进行重采样，将所有的波段重采样为20 m分辨率，并输出为ENVI格式。最终去除B10波段(1 375 nm)，得到相同空间分辨率的12个波段。Sentinel-2A的1C和2A 级产品的反射率会乘上一个固定的系数。该系数默认为10 000(从头文件的QUANTIFICATION Value语句中查得)[43]。将重采样后的结果导入ENVI进行波段融合(Layer Stacking)，得到一个多波段文件。将其各波段的灰度值DN除以10 000，还原为地表反射率。计算公式为
ρλ=Qcal/10 000
(1)
其中，ρλ为λ波段的地表反射率;Qcal为影像以16 bit量化的DN值。正常的DN值应该在0～1，上式计算后如果存在DN<0,DN>1，可将小于0赋值为0，大于1的赋值为1。
2.2 植被红边指数
植被指数(VI)是一种利用遥感影像来监测地表植被状态的方法。红光和红外区间的反射率会随绿色植被的不同覆盖而变化。组合这2个波段信息可以用来区分土壤和植被。一定条件下，不同的植被指数在可以用来定量说明植被在某一阶段的生长状况。通过建立VI，综合有关的光谱信号，使与植被有关信息增强，而使与植被无关信息减弱。
目前，在科学文献中发布了超过150种植被指数模型，主要筛选出归一化植被指数(NDVI)及5种红边植被指数—改进红边归一化植被指数(Red Edge Normalized Difference Vegetation Index,NDVIre)、归一化差值红边(Normalized Difference Red Edge,NDre)、修改的叶绿素反射率指数(Modified Chlorophyll Absorption in Reflectance Index,MCARI)、修改的红边简单比值指数(Modified Red Edge Simple Ratio,MSRre)和三角叶绿素指数(Triangular Chlorophyll Index,TCI)。通过IDL批量计算11个期次的NDVI,NDVIre,NDre,MCARI,MSRre和TCI指数。NDVI和5种红边植被指数表达式，见表3。
表3 植被指数计算公式及对应波段Table 3 Calculation formula and corresponding bands of Vegetation Index
由于Sentinel-2A的波段的中心波长与所选用植被指数对应的中心波长不完全重合，因此利用与Sentinel-2A相邻近波段进行波段运算见表3。
2.3 剖线与信息提取
113101工作面是泊江海子矿首采工作面，开采时间2016-12-01—2017-08-31，其沉陷区域随开采由南向北移动。参考InSAR结果，沿113101工作面中心位置设计1条从北走向南(N—S走向)的剖面线，并在工作面东西两侧各设计1条平行的剖面线，剖面线长约6.6 km，每1条剖线包含334个像元。从西向东分别命名为对比区1(CK1)、沉陷区(S1)、对比区2(CK2)(图3)，其中，图3(a)中蓝色箭头为开采方向，图3(b)中黄色箭头为栅格排序方向，334个像元从北向南编号为1，2，…，334。
图3 InSAR剖面线与光学影像叠加分析Fig.3 Superposition analysis of InSAR section line and optical image
2.4 趋势分析与检验
根据其研究区剖线对应的植被指数，利用长时序的遥感影像数据，建立沉陷区内外植被指数随时间的变化曲线。通过比较其植被指数在一时间段内的各种变化检测指标的生长期曲线的差异，获取研究区植被的物候信息。主要采用剖面分析法和非对称高斯函数拟合法。
2.4.1剖面分析
剖面分析法主要为了更精确的分析沉陷区内外植被分布情况，是多时相趋势分析的数据基础。根据植被指数NDVI的特点，当旺季植被指数NDVI<0，则为水体;0
2.4.2回归分析
(1)线性回归。在剖面数据的基础上，建立自变量与因变量之间的相关关系，从而建立一个具有良好相关性的回归方程，即函数表达式。主要采用最小二乘法对数据进行一元线性回归，回归系数k主要是用来反映植被指数长期来的变化趋势。
y=a+kxi+εi
(2)
式中，a为常数项;k为回归系数;εi为残差。
(2)二元Logistic回归。二元Logistic回归是因变量为二分类的广义线性回归模型。二分响应变量取值为0或者1，0表示事件未发生，1表示事件发生，正好用来检验自变量(如NDVI)是否受开采沉陷影响。我们设剖面上受采动影响的像元值为1，则未受采动影响的像元值为0，构建如下logistic回归模型[45]:
(3)
式中，Y为二分响应变量;P为二分类变量取特定值时的概率;α0为常数项；Xi(i=1，2，…，m)为剖线上第i个像元的NDVI值;βi为自变量Xi的回归系数。
2.4.3高斯拟合
高斯函数是表示连续随机变量的概率密度函数。它通常用于表达自然物候的周期性演变，可以表示植被枯—荣—枯(或荣—枯—荣)过程[46]。函数(单峰型)的一维表达式为
(4)
式中，y0为基线的偏移量;A为钟形曲线下侧的积分面积(NDVI等效生物量);w大约等于峰高一半处宽度的0.849倍(当峰高为1/2时，钟形曲线的宽度),该值可以表示植被生长繁茂时的生长期的长度;x0为峰值位置。
2.4.4变化幅度
变化幅度指的是沉陷区变化值与非沉陷区变化值相比，增长或者下降的比率。函数表达式为
(5)
式中，V为变化幅度，当V为正值时为增长幅度，为负值时为下降幅度；δ1为沉陷区差值；δ2为非沉陷区差值。
2.4.5决定系数
决定系数R2，又称拟合优度。R2越高，模型的拟合效果越好，自变量对因变量的解释程度就越好，自变量本身所造成的变化百分比也就越高。R2的取值范围(0，1)，没有单位。若观测点直接落在拟合曲线上，则R2=1，此时拟合曲线可以解释因变量的所有变化。
(6)
(7)
其中,n为样本数量;p为特征数量。即样本为n个(x1，x2，x3，…,xp,y)。
3 结果与分析
3.1 时序NDVI剖面
获取研究区2016年10月至2017年10月的植被指数时序变化剖面图(S1，CK1，CK2)有3个作用:① 反映剖面植被类型及植被覆被情况;② 认识剖面植被的物候特征及生物量(和值、均值)波动情况;③ 分析各剖面相关性，形成定量化比较的数据基础。
NDVI是使用最广泛的植被指数，已被证明是监测植被状况的最佳指标之一。以此为例绘制S1，CK1，CK2时间序列剖面，如图4所示，其中，黄色线段部分为开采沉陷影响区间，其余像元为非沉陷部分。
沉陷区剖面线(S1)的NDVI值时序图(图4(b))，该剖线上从12月开始形成沉陷，沉陷区域由南向北移动。其中，S为该期该剖线NDVI总和，M为均值。黄色线段部分为开采沉陷影响区间，其余像元为非沉陷部分。根据生长旺季(2017年8月)NDVI阈值，可知研究区地表覆被主要中等覆被类型的草地和湿草地;每年的4—10月为植被的生长期，S1的NDVI和值较高，但在S1开采前后没有明显变化，应是季节变化掩盖了采动影响;3条剖面NDVI同一季节相关性较高，说明地表覆被类型相近。
图4 NDVI剖线时序Fig.4 Sequence diagram of NDVI section line
3.2 时序植被指数均值
将3条剖线(S1，CK1，CK2)上的植被指数的每1期均值(ΣVI/334)做一个散点图，按照每月1期的时间间隔获取该研究区2016年10月至2017年10月的植被指数时序变化(图5)。可见在一个生命周期内的植被指数呈单峰高斯分布，因此采用高斯函数拟合曲线。
图5 6种植被指数均值的高斯分布曲线Fig.5 Gaussian distribution curves of mean value of 6 vegetation indexes
根据变化趋势可知，研究区所处的荒漠化草地在4月之后进入返青期，7，8月份植被的生长最为旺盛，大约在8月初达到峰值，10月之后进入衰退期。沉陷区与对比区的植被生长状况有所差异，但差异不大。
表4 4—8月间6种植被指数的线性回归Table 4 Linear regression of six Vegetation Indices from April to August
均值回归结果表明，开采沉陷造成沉陷区植被的增长速度低于对比区，但绝对影响量有多少，有待开展增量研究。
3.3 时序植被指数差值
为了定量化表达沉陷前后植被指数的变化，设计了差值回归方法。主要是针对采动影响，分析采动期和非采动期两个时期植被指数的变化。即用采动时期的所有像元减去非采动时期相对应像元的植被指数，求出采动引起的VI增量。
考虑到剖线上的植被分布状况有所差异，仅利用S1剖面线的数据，由InSAR获得的时间序列沉陷区可知，第1次卫星过境时(2016-10-19)未发生沉陷，以这一期数据作为基准数据，结合InSAR技术获得的沉陷区具体位置，剖面线上落入沉陷区的像元位置见表5，其余的为非沉陷像元。
表5 沉陷区剖面线上涉及沉陷的像元位置Table 5 Pixel position of subsidence on the section line of subsidence area
用每一期沉陷区和非沉陷区像元对应植被指数VI的均值，减去2016-10-19那一期对应位置像元植被指数VI的均值，得到植被指数增量δVI，分析沉陷区和非沉陷区增量的时序变化趋势。
以卫星数据接收时间为横轴，绘制植被指数增量的散点图并采用高斯回归(图6)。可见6种植被植数增量δNDVI，δNDVI，re，δND，re,δMCARI,δTCI和δMSR,re在沉陷区增量均小于非沉陷区增量。这表明受沉陷胁迫的影响，沉陷区植被的长势要比非沉陷区差。这与3.2节结论一致。通过研究沉陷区相对于非沉陷区的下降幅度，发现6种植被指数增量δNDVI，δNDVI，re，δND，re,δMCARI,δTCI和δMSR,re平均下降幅度分别为24.1%，18.8%，25.9%，32.8%，28.4%，23.1%，沉陷区的植被的增长速率均低于非沉陷区的增长速率。
图6 6种植被指数差值时序变化分析Fig.6 Time series analysis of difference of 6 Vegetation Indexes
对6种植被指数增量(δVI)时间序列进行了高斯回归，δNDVI，δNDVI，re，δND，re,δMCARI,δTCI和δMSR,re回归方程的决定系数分别为0.954 91，0.961 01，0.945 52，0.974 5，0.982 11和0.963 32。易见，当NDVI中的近红外波段替换为红边波段变为NDVIre后，决定系数R2有明显改善。
3.4 Logistic 回归分析
为了确认各植被指数是否受开采沉陷影响，笔者构建了6种植被指数每一期的Logistic回归模型，然后进行逻辑斯蒂回归分析。检验情况见表6。其中显著性水平Sig又称为P值，当P<0.05时，表明该模型是显著的。
表6显示，采动对各植被指数响应在非生长季没有显著影响，即在非生长季其统计学意义不大;而在生长季会有显著影响。且除MCARI，TCI以外的4种植被指数模型在植被生长较为旺盛的4，5，6，8月份均是显著的，在植被生长最旺盛的8月份，6个植被指数均显著相关。表明在生长季，采矿活动对植被指数有显著影响。
表6 6种植被指数响应采动影响的Logistic回归的显著水平Table 6 Significant level of Logistic regression of six vegetation indexes in response to mining
4 结 论
(1)在植被一个生命周期内，6种植被指数值曲线呈明显单峰变化。研究区地表的荒漠化草地在4月启动返青期，7，8月份进入生长旺季，10月之后进入衰退期。
(2)对沉陷区(S1)、对比区1(CK1)、对比区2(CK2)的4—8月的6种植被指数均值进行线性回归分析，沉陷区的回归方程的斜率83%小于对比区(非沉陷区)，沉陷区植被的增长速度慢于非沉陷区。
(3)沉陷区和非沉陷区1 a内植被指数增量(δVI)的时序分析表明，沉陷区增量均小于非沉陷区增量。表明植被受沉陷胁迫，沉陷区植被的长势要比非沉陷区差。采矿活动使6种植被指数增幅减少18.8%～32.8%。
(4)在植被生长期，二元Logistic回归证实6种植被指数与采矿活动相关或显著相关，采矿活动对植被指数有显著影响。
黄河流域中上游的黄土-风积沙矿区，气候干旱，生态脆弱，植被红边指数对采矿胁迫极为敏感。这为研究黄河流域中上游煤矿生态修复的理论与方法，实现黄河流域高质量发展，提供了敏感性极佳的生态指标。本文的研究结果，回答了开采沉陷对地表植被有没有影响？且部分回答了有多少影响？的问题。但对于影响多久？和长期影响规律如何？等问题的解答，仍有待长期实践和深入研究。
致谢感谢中国-意大利合作培养博士、边山大学访问学者刘培副教授对英文摘要的润色。
