采煤胁迫下干旱区植被生长周期红边指数动态变化

李晓婷1,杨丽帆2,邹友峰1,3,马 超1,3

(1.河南理工大学 测绘与国土信息工程学院,河南 焦作 454003; 2.中国地质大学 地理与信息工程学院,湖北 武汉 430074; 3.河南理工大学 自然资源部矿山时空信息与生态修复重点实验室,河南 焦作 454003)

摘 要:红边指数对植被叶绿素浓度、叶表面冠层、叶聚丛和冠层结构敏感。研究采矿胁迫下红边指数的响应规律,可为沉陷区植被损伤诊断、植被自然恢复和人工修复效果评估提供生态学指标。以黄河流域鄂尔多斯市泊江海子矿113101首采工作面为例,结合Sentinel-1A 数据D InSAR技术获得的形变场,同步选取2016年10月至2017年10月间11期带有红边波段的Sentinel-2A多时相遥感数据,利用归一化植被指数(NDVI)和5种典型红边植被指数(NDVIre,NDre,MCARI,TCI,MSRre),同步分析沉陷区植被指数的时空变化规律。结果表明:① 一个生命周期的植被指数呈GAUSS单峰变化,研究区4月进入返青期,8月达到峰值,10月进入衰退期。② 沉陷区6种植被指数线性回归斜率都小于对照区,沉陷区植被的增长速度慢于非沉陷区。③ 受沉陷扰动的影响,沉陷区植被指数的增量都小于非沉陷区,6种植被指数平均降幅达18.8%~32.8%。④ 二元Logistic回归显示,除MCARI,TCI以外的4种植被指数在植被生长期的4,5,6,8月份与采矿活动相关或显著相关(显著性水平P为0.001~0.100)。在生长季,采矿活动对工作面上方的植被指数有显著影响。基于红边波段组合的植被指数能敏感反映采矿胁迫下植被生长状况的细微变化,可以作为黄河流域黄土—风积沙矿区植被损伤诊断及生态环境评价的优势指标。

关键词:开采沉陷;采矿胁迫;红边指数;Sentinel-2A;时空变化;Logistic回归

煤矿开采对塌陷区地表植被影响及恢复的程度、范围和时间,一直是学术界争论的焦点问题之一[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波段中心波长/nm带宽/nm分辨率/m波段特点B1-Coastal aerosol4432060可用于大气气溶胶校正B2-Blue4906510可用于监测海洋及浮游生物B3-Green5603510可用于监测水体B4-Red6653010可用于监测植物生长及地表覆盖B5-Vegetation Red Edge7051520B6-Vegetation Red Edge7401520可用于监测植物色素及健康状况B7-Vegetation Red Edge7832020B8-NIR84211510可用于获取植被的叶面指数B8a-Narrow NIR8652020可用于植被监测B9-Water vapor9452060可用于大气校正(水蒸气吸收)B10-SWIR-Cirrus1 3753060可用于监测薄卷云B11-SWIR11 6109020可用于监测云和雪B12-SWIR22 19012020可用于监测植被和地质

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

序号成像日期相差时间/d传感器数据类型帧号相对轨道号12016-10-190MSILevel-1CN0204R06122016-11-1527MSILevel-1CN0204R01832016-12-1833MSILevel-1CN0204R06142017-01-1427MSILevel-1CN0204R01852017-02-1330MSILevel-1CN0204R01862017-03-1833MSILevel-1CN0205R06172017-04-2740MSILevel-1CN0205R06182017-05-1720MSILevel-1CN0205R06192017-06-1327MSILevel-1CN0205R018102017-08-0520MSILevel-1CN0205R061112017-08-2520MSILevel-1CN0205R061

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

植被指数表达式参考文献NDVIρ800-ρ670ρ800+ρ670[39]NDVIreρ730-ρ670ρ730+ρ670[32]NDreρ790-ρ720ρ790+ρ720[33]MCARIρ700ρ670[(ρ700-ρ670)-0.2(ρ700-ρ550)][34]MSRreρ750ρ705-1() ρ750ρ705+1()[35]TCI1.2ρ700-ρ550()-1.5(ρ670-ρ550)× ρ700ρ670[36]

由于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,则为水体;00.6为灌丛及林地[44]

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)

其中:SSR为回归平方和;SST为总平方和;yi(i=1,2,…,n)为待拟合数值,其均值为拟合值为

自变量个数越多,决定系数R2越大,当增加自变量没有意义时,决定系数R2还是会越大。所以在比较回归方程时,还需要看自变量个数,考虑方程中自变量个数的影响。因此需对其进行校正,从而形成了校正决定系数R2不同的是,当回归模型中增加的自变量没有意义时,不会增大,反而会减小,因此是衡量构建的回归模型优劣的一个重要指标,校正R2越大,模型拟合效果越优。但当p/n 很小时,比如(p/n)< 0.05时,校正作用减弱甚至趋向消失[47]

(7)

其中,n为样本数量;p为特征数量。即样本为n个(x1x2x3,…,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月开始返青,8月达到峰值,对4—8月的植被指数进行回归分析,这期间植被指数线性比较明显,回归方程为线性方程,则回归系数(方程斜率)可反映植被指数在该时间段内的倾向。对比发现,3个剖面6种植被指数共获得18个回归方程中,12条对比区剖面中有10条剖面斜率大于沉陷区斜率,约占83%(表4),沉陷区大部分小于对比区,表明沉陷区植被的增长速度慢于对比区(即非沉陷区)。

表4 4—8月间6种植被指数的线性回归
Table 4 Linear regression of six Vegetation Indices from April to August

植被指数剖面线回归方程R2Adj植被指数剖面线回归方程R2AdjCK1y=0.002 03x-4 978.230.943 3CK1y=0.000 35x-851.1760.998 0NDVIS1y=0.001 88x-4 631.880.937 1MCARIS1y=0.000 28x-705.6480.988 8CK2y=0.001 87x-4 603.030.952 7CK2y=0.000 28x-701.0160.983 5CK1y=0.002 12x-5 204.190.974 4CK1y=0.000 36x-886.3010.998 1NDVIreS1y=0.002 15x-5 292.530.980 6TCIS1y=0.000 30x-742.7320.989 5CK2y=0.002 30x-5 657.340.977 7CK2y=0.000 31x-751.1500.970 4CK1y=0.001 22x-2 911.460.943 5CK1y=0.002 26x-5 562.000.978 7NDreS1y=0.000 83x-2 053.580.975 3MSRreS1y=0.002 02x-4 958.240.978 5CK2y=0.001 09x-2 680.870.965 7CK2y=0.002 17x-5 342.940.967 5

均值回归结果表明,开采沉陷造成沉陷区植被的增长速度低于对比区,但绝对影响量有多少,有待开展增量研究。

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

日期起点像元位置终点像元位置涉及像元数量/个日期起点像元位置终点像元位置涉及像元数量/个2016-12-18238252152017-05-17168234672017-01-14213262502017-06-13140234952017-02-13213254422017-08-051272341082017-03-18172255842017-08-251202341152017-04-2717823962

注:剖面栅格编号由北向南,工作面开采影响的栅格编号与此相同。

用每一期沉陷区和非沉陷区像元对应植被指数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

数据获取日期总样本数未受采动影响样本数受采动影响样本数显著性水平(P值)NDVINDVIreNDreMCARITCIMSRre植被物候期2016-12-18334317170.6370.5910.5040.7580.8980.461非生长季2017-01-14334284500.6220.2550.8170.7690.8920.979非生长季2017-02-13334292420.6160.1390.2800.1110.0600.912非生长季2017-03-18334250840∗∗0.3680.014∗0.9990.6410.505非生长季2017-04-27334272620.014∗0.0750.009∗∗0.3890.1710.010∗∗生长季2017-05-17334267670.001∗∗0.032∗0.010∗∗0.2600.1070.011∗生长季2017-06-13334249950.001∗∗0.006∗∗0.001∗∗0.1910.0850.001∗∗生长季2017-08-053342261080.004∗∗0.009∗∗0.018∗0.004∗∗0.001∗∗0.026∗生长季2017-08-253342191150∗∗0∗∗0.005∗∗0∗∗0∗∗0.005∗∗生长季

注:*为显著(P<0.05);**为非常显著(P<0.01)。

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种植被指数与采矿活动相关或显著相关,采矿活动对植被指数有显著影响。

黄河流域中上游的黄土-风积沙矿区,气候干旱,生态脆弱,植被红边指数对采矿胁迫极为敏感。这为研究黄河流域中上游煤矿生态修复的理论与方法,实现黄河流域高质量发展,提供了敏感性极佳的生态指标。本文的研究结果,回答了开采沉陷对地表植被有没有影响?且部分回答了有多少影响?的问题。但对于影响多久?和长期影响规律如何?等问题的解答,仍有待长期实践和深入研究。

致谢 感谢中国-意大利合作培养博士、边山大学访问学者刘培副教授对英文摘要的润色。

参考文献(References):

[1] 胡振琪,龙精华,王新静.论煤矿区生态环境的自修复、自然修复和人工修复[J].煤炭学报,2014,39(8):1751-1757.

HU Zhenqi,LONG Jinghua,WANG Xinjing.Self-healing,natural restoration and artificial restoration of ecological environment for coal mining[J].Journal of China Coal Society,2014,39(8):1751-1757.

[2] 卞正富,雷少刚,金丹,等.矿区土地修复的几个基本问题[J].煤炭学报,2018,43(1):190-197.

BIAN Zhengfu,LEI Shaogang,JIN Dan,et al.Several basic scientific issues related to mined land remediation[J].Journal of China Coal Society,2018,43(1):190-197.

[3] SOBRINO J A,JULIEN Y.Global trends in NDVI-derived parameters obtained from GIMMS data[J].International Journal of Remote Sensing,2011,32(15):4267-4279.

[4] 侯美亭,赵海燕,王筝,等.基于卫星遥感的植被NDVI对气候变化响应的研究进展[J].气候与环境研究,2013,18(3):353-364.

HOU Meiting,ZHAO Haiyan,WANG Zheng,et al.Vegetation responses to climate change by using the satellite-derived normalized difference vegetation index:A review[J].Climatic and Environmental Research,2013,18(3):353-364.

[5] MEZNED N,MECHERGUI M,ALAYET F,et al.Hyperspectral spectroscopy based vegetation stress detection around Pb-Zn-Cd abandoned mining site in Northern Tunisia[A].Geoscience & Remote Sensing Symposium[C].Waikdaa:IEEE,2016.

[6] LI Xuqing,LI Long,LIU Xiangnan.Collaborative inversion he-avy metal stress in rice by using two-dimensional spectral feature space based on HJ-1 A HSI and radarsat-2 SAR remote sensing data[J].International Journal of Applied Earth Observation & Geoinformation,2019,78:39-52.

[7] FU Pingjie,ZHANG Wei,YANG Keming,et al.A novel spectral analysis method for distinguishing heavy metal stress of maize due to copper and lead:RDA and EMD-PSD[J].Ecotoxicology and Environmental Safety,2020,206:111211.

[8] WU Ling,LIU Xiangnan,WANG Ping,et al.The assimilation of spectral sensing and the WOFOST model for the dynamic simulation of cadmium accumulation in rice tissues[J].International Journal of Applied Earth Observation and Geoinformation,2013,25:66-75.

[9] FENG Xuping,CHEN Houming,CHEN Yue,et al.Rapid detection of cadmium and its distribution in miscanthus sacchariflorus based on visible and near-infrared hyperspectral imaging[J].Science of the Total Environment,2019,659:1021-1031.

[10] TESNEROV C,ZADINOV R,PIKL M,et al.Predicting the toxicity of post-mining substrates,a case study basedon laboratory tests,substrate chemistry,geographic information systems and remote sensing[J].Ecological Engineering 2017,100:56-62.

[11] WANG Fenghe,GAO Jay,ZHA Yong.Hyperspectral sensing of heavy metals in soil and vegetation:Feasibility and challenges[J].ISPRS Journal of Photogrammetry and Remote Sensing,2018,136:73-84.

[12] STAGAKIS S,GONZLEZ-DUGO V,CID P,et al.Monitoring water stress and fruit quality in an orange orchard under regulated deficit irrigation using narrow-band structural and physiological remote sensing indices[J].ISPRS Journal of Photogrammetry & Remote Sensing,2012,71:47-61.

[13] FANG Wei,HUANG Shengzhi,HUANG Qiang,et al.Probabilistic assessment of remote sensing-based terrestrial vegetation vulnerability to drought stress of the Loess Plateau in China[J].Remote Sensing of Environment,2019,232:111290.

[14] LEES K J,ARTZ R R E,CHANDLER D,et al.Using remote sensing to assess peatland resilience by estimating soil surface moisture and drought recovery[J].Science of the Total Environment,2020,12:143312.

[15] JAVED T,LI Y,RASHID S,et al.Performance and relationship of four different agricultural drought indices for drought monitoring in China’s mainland using remote sensing data[J].Science of the Total Environment,2021,759:143530.

[16] CRLAN I,MIHAI B A,NISTOR C,et al.Identifying urban vegetation stress factors based on open access remote sensing imagery and field observations[J].Ecological Informatics,2020,55:101032.

[17] WEI Chuanwen,HUANG Jingfeng,WANG Xiuzhen,et al.Hyperspectral characterization of freezing injury and its biochemical impacts in oilseed rape leaves[J].Remote Sensing of Environment,2017,195:56-66.

[18] DASH J P,WATT M S,PEARSE G D,et al.Assessing very high resolution UAV imagery for monitoring forest health during a simulated disease outbreak[J].ISPRS Journal of Photogrammetry and Remote Sensing,2017,131:1-14.

[19] AVTAR R,KOMOLAFE A A,KOUSER A,et al.Assessing sustainable development prospects through remote sensing:A review[J].Remote Sensing Applications:Society and Environment,2020,20:100402.

[20] TIAN Fei,HOU Mengjie,QIU Yuan,et al.Salinity stress effects on transpiration and plant growth under different salinity soil levels based on thermal infrared remote (TIR) technique[J].Geoderma,2020,357:113961.

[21] YUAN Lin,BAO Zhiyan,ZHANG Haibo,et al.Habitat monitoring to evaluate crop disease and pest distributions based on multi-source satellite remote sensing imagery[J].Optik International Journal for Light & Electron Optics,2017,145:66-73.

[22] ZHANG Jingcheng,HUANG Yanbo,PU Ruiliang,et al.Monitoring plant diseases and pests through remote sensing technology:A review[J].Computers and Electronics in Agriculture,2019,165:104943.

[23] ACHARYA B S,KHAREL G.Acid mine drainage from coal mining in the United States:An overview[J].Journal of Hydrology,2020,588:125061.

[24] 彭苏萍,毕银丽.黄河流域煤矿区生态环境修复关键技术与战略思考[J].煤炭学报,2020,45(4):1211-1221.

PENG Suping,BI Yinli.Strategic consideration and core technology about environmental ecological restoration in coal mine areas in the Yellow River basin of China[J].Journal of China Coal Society,2020,45(4):1211-1221.

[25] 马超,田淑静,邹友峰,等.神东矿区AVHRR/NDVI的时空、开采强度和气候效应[J].中国环境科学,2016,36(9):2749-2756.

MA Chao,TIAN Shujing,ZOU Youfeng,et al.Dynamic responses of the coalfield ecosystem to mining intensitu,spatio-temporal,and climate change derived from AVHRR/NDVI in Shendong coalfield[J].China Environmental Science,2016,36(9):2749-2756.

[26] 徐占军,侯湖平,张绍良,等.采矿活动和气候变化对煤矿区生态环境损失的影响[J].农业工程学报,2012,28(5):232-240.

XU Zhanjun,HOU Huping,ZHANG Shaoliang,et al.Effects of mining activity and climatic change on ecological losses in coal mining areas[J].Transactions of the Chinese Society of Agricultural Engineering,2012,28(5):232-240.

[27] 吴立新,马保东,刘善军.基于SPOT卫星NDVI数据的神东矿区植被覆盖动态变化分析[J].煤炭学报,2009,34(9):1217-1222.

WU Lixin,MA Baodong,LIU ShanJun.Analysis to vegetation coverage change in Shendong mining area with SPOT NDVI data[J].Journal of China Coal Society,2009,34(9):1217-1222.

[28] 翟孟源,徐新良,江东,等.1979~2010年乌海市煤矿开采对生态环境影响的遥感监测[J].遥感技术与应用,2012,27(6):933-940.

ZHAI Mengyuan,XU Xinliang,JIANG Dong,et al.Remote sensing monitoring of the ecological environment in Wuhai mining area since 1979[J].Remote Sensing Technology and Application,2012,27(6):933-940.

[29] 雷少刚.荒漠矿区关键环境要素的监测与采动影响规律研究[D].徐州:中国矿业大学,2009:45-52.

LEI Shaogang.Monitoring and analyzing the mining impacts on key environmental elements in desert area[D].Xuzhou:China University of Mining Science and Technology,2009:45-52.

[30] 马雯思,马超,刘玮玮.陕西彬长矿区NDVI3g(1982—2013)变化趋势及气候响应[J].煤炭学报,2019,44(4):1197-1206.

MA Wensi,MA Chao,LIU Weiwei.Variation and climate response of NDVI3g (1982—2013) in Binchang mining area of Shaanxi Province[J].Journal of China Coal Society,2019,44(4):1197-1206.

[31] MA Chao,GUO Zengzhang,ZHANG Xiaoke,et al.Annual integral changes of time serial NDVI in mining subsidence area[J].Transactions of Nonferrous Metals Society of China,2011,21:583-588.

[32] ELSAYED S,RISCHBECK P,SCHMIDHALTER U.Comparing the performance of active and passive reflectance sensors to assess the normalized relative canopy temperature and grain yield of drought-stressed barley cultivars[J].Field Crops Research,2015,177:148-160.

[33] BARNES E M,CLARKE T R,Richards S E,et al.Coincident detection of crop water stress,nitrogen status and canopy density using ground-based multi-spectral data[A].International Conference on Precision Agriculture and Other Resource Management[C].Bloomington,2000:16-19.

[34] DAUGHTRY C S T,WALTHALL C L,KIM M S,et al.Estimating corn leaf chlorophyll concentration from leaf and canopy reflectance[J].Remote Sensing of Environment,2000,74(2):229-239.

[35] WU Chaoyang,NIU Zheng,TANG Quan,et al.Estimating chlorophyll content from hyper-spectral vegetation indices:Modeling and validation[J].Agricultural and Forest Meteorology,2008,148(8):1230-1241.

[36] HABOUDANE D,TREMBLAY N,MILLER J R,et al.Remote estimation of crop chlorophyll content using spectral indices derived from hyperspectral data[J].IEEE Transactions on Geoscience and Remote Sensing,2008,46(2):423-437.

[37] VERRELST J,MUOZ J,ALONSO L,et al.Machine learning regression algorithms for biophysical parameter retrieval:Opportunities for Sentinel-2 and-3[J].Remote Sensing of Environment,2012,118:127-139.

[38] CLEVERS J G P W,GITELSON A A.Remote estimation of crop and grass chlorophyll and nitrogen content using red-edge bands on Sentinel-2 and-3[J].International Journal of Applied Earth Observations and Geoinformation,2013,23:344-351.

[39] 郑阳,吴炳方,张淼.Sentinel-2数据的冬小麦地上干生物量估算及评价[J].遥感学报,2017,21(2):318-328.

ZHENG Yang,WU Bingfang,ZHANG Miao.Estimating the above ground biomass of winter wheat using the Sentinel-2 data[J].Journal of Remote Sensing,2017,21(2):318-328.

[40] 赵玉倩.泊江海子煤矿井工开采引起地表变形的问题研究[D].北京:中国地质大学,2015.

ZHAO Yuqian.Research on surface deformation caused by underground mining in BOJIANG HAIZI coal mine[D].Beijing:China University of Geosciences,2015.

[41] 田颖,陈卓奇,惠凤鸣,等.欧空局哨兵卫星Sentinel-2 A/B数据特征及应用前景分析[J].北京师范大学学报(自然科学版),2019,55(1):57-65.

TIAN Ying,CHEN Zhuoqi,HUI Fengming,et al.ESA Sentinel-2 A/B satellite:Characteristics and applications[J].Journal of Beijing Normal University(Natural Science),2019,55(1):57-65.

[42] 苏伟,张明政,蒋坤萍,等.Sentinel-2卫星影像的大气校正方法[J].光学学报,2018,38(1):322-331.

SU Wei,ZHANG Mingzheng,JIANG Kunping,et al.Atmospheric correction method for Sentinel-2 satellite imagery[J].Acta Optica Sinica,2018,38(1):322-331.

[43] 方灿莹,王琳,徐涵秋.不同植被红边指数在城市草地健康判别中的对比研究[J].地球信息科学学报,2017,19(10):1382-1392.

FANG Canying,WANG Lin,XU Hanqiu.A comparative study of different red edge indices for remote sensing detection of urban grass land health status[J].Journal of Geo-Information Science,2017,19(10):1382-1392.

[44] 马超,蔡盼丽.徂徕山-莲花山地区环境生态指数时空变化与驱动因素分析[J].国土资源遥感,2019,31(4):199-208.

MA Chao,CAI Panli.Spatio-temporal changes and driving factors of environmental and ecological index in Culai-Lianhua area[J].Remote Sensing for Land and Resources,2019,31(4):199-208.

[45] 王杰,刘实,兰玉彬,等.基于 HJ-1A/B CCD 数据的玉米倒伏遥感监测方法[J].中国农业气象,2020,41(2):121-128.

WANG Jie,LIU Shi,LAN Yubin,et al.Method of maize lodging recognition based on HJ-1A/B CCD Data[J].Chinese Journal of Agrometeorology,2020,41(2):121-128.

[46] 李敏,盛毅.高斯拟合算法在光谱建模中的应用研究[J].光谱学与光谱分析,2008,28(10):2352-2355.

LI Min,SHENG Yi.Study on application of Gaussian fitting algorithm to building model of spectral analysis[J].Spectroscopy and Spectral Analysis,2008,28(10):2352-2355.

[47] 张文彤,董伟.SPSS统计分析高级教程.第2版[M].北京:高等教育出版社,2013:117-121.

Dynamic change of red edge vegetation index within a growth cycle in arid area under coal mining stress

LI Xiaoting1,YANG Lifan2,ZOU Youfeng1,3,MA Chao1,3

(1.College of Surveying & Land Information Engineering,Henan Polytechnic University,Jiaozuo 454003,China; 2.School of Geography and Information Engineering,China University of Geosciences,Wuhan 430074,China; 3.Key Laboratory of Spatio-temporal Information and Ecological Restoration of Mines (Ministry of Natural Resources,MNR),Henan Polytechnic University,Jiaozuo 454003,China)

Abstract:Red Edge Vegetation Index (VIre) is sensitive to chlorophyll concentration,leaf surface canopy,leaf cluster and canopy structure.The research on the response characteristics of VIre under mining stress can provide quantitative ecological indicators for vegetation damage diagnosis,natural vegetation restoration and artificial reclamation effectiveness evaluation in mining subsidence area.Taking the first mining face (No.113101) of Bojianghaizi Mine in Ordos City along the upper and middle Yellow River Basin as research instance,the mining-induced deformation field was obtained by DInSAR technology using Sentinel-1A SAR remotely sensed data,meanwhile,based on the traditional vegetation index,the red edge band was introduced to improve the existing vegetation index,the spatiotemporal variation of the vegetation indices in subsidence was analyzed with the help of normalized difference vegetation index (NDVI) and 5 typical VIre:NDVIre (Red Edge Normalized Difference Vegetation Index),NDre(Normalized Difference Red Edge),MCARI (Modified Chlorophyll Absorption in Reflectance Index),TCI (Triangular Chlorophyll Index),MSRre (Modified Red Edge Simple Ratio)-inversed from 11 periods of Sentinel-2A multi-temporal remote sensing data from October 2016 to October 2017.The results showed that ① the vegetation indices within a life cycle exhibit a single peak change of GAUSS.In this research area,the starting date of growing season (SOS) is in April,the peak is in August,and the ending date of the season (EOS) is in October.② The linear regression of 6 vegetation indices in subsidence area is smaller than that in check area (CK),which indicates that the growth rate of vegetation in subsidence area is slower than that in non-subsidence area.③ By studying the decline of subsidence area relative to non-subsidence area,it is found that affected by the subsidence disturbance,the increment of vegetation indices (δVIs) in the subsidence area is smaller than that in the non-subsidence area,and the average decrease of 6 vegetation indices is 18.8%-32.8%.The fitting effect of TCI is the best (R2=0.982 1).It further shows that the vegetation growth speed above the subsidence area is slower due to mining,and the red edge index is more sensitive to the vegetation growth in the mining subsidence area.④ Binary logistic regression shows that 4 vegetation indices (except for MCARI and TCI) are significantly correlated with coal mining activities in April,May,June and August (Significant level 0.001-0.100).However,in August,the 6 vegetation indices are significantly correlated (Significant level 0.001-0.026),and the significance level of TCI model is the highest.In the peak season of vegetation growth,coal mining activities have a significant impact on vegetation index above the working face.The vegetation indices based on red edge band combination can sensitively reflect the subtle changes of vegetation growth under mining stress,and can be used as the dominant indicators for the vegetation damage diagnosis and eco-environmental assessment of loess-aeolian coal mine in the Yellow River Basin.

Key words:mining subsidence;coal mining stress;Red Edge Vegetation Indices (VIre);Sentinel-2A;spatiotemporal change;logistic regression

中图分类号:TD88;Q948.1

文献标志码:A

文章编号:0253-9993(2021)05-1508-13

移动阅读

收稿日期:20210127

修回日期:20210428

责任编辑:黄小雨

DOI:10.13225/j.cnki.jccs.ST21.0186

基金项目:国家自然科学基金项目重点资助项目(U1261206);国家自然基金委与英国皇家学会合作资助项目(4191101440)

作者简介:李晓婷(1998—),女,河南郑州人,硕士研究生。E-mail:212004010013@home.hpu.edu.cn

通讯作者:马 超(1967—),男,内蒙古克什克腾旗人,教授,博士生导师。E-mail:mac@hpu.edu.cn

引用格式:李晓婷,杨丽帆,邹友峰,等.采煤胁迫下干旱区植被生长周期红边指数动态变化[J].煤炭学报,2021,46(5):1508-1520.

LI Xiaoting,YANG Lifan,ZOU Youfeng,et al.Dynamic change of red edge vegetation index within a growth cycle in arid area under coal mining stress[J].Journal of China Coal Society,2021,46(5):1508-1520.