半干旱煤炭开采沉陷区植被扰动响应的时间特征

陶文旷，雷少刚

(1.中国矿业大学环境与测绘学院，江苏 徐州 221116；2.中国矿业大学国土资源研究所，江苏 徐州　221116)



半干旱煤炭开采沉陷区植被扰动响应的时间特征

陶文旷1，雷少刚2①

(1.中国矿业大学环境与测绘学院，江苏 徐州 221116；2.中国矿业大学国土资源研究所，江苏 徐州221116)

摘要：利用2005—2012年MODIS影像和TM影像的月均归一化差分植被指数(NDVI)时间序列数据以及神东矿区大柳塔地区气象数据，采用基于Loess算法的时间序列季节和趋势成分分解法(STL)、时滞互相关分析和改进后的残差分析等方法，探讨大柳塔地区2005—2012年植被NDVI随降雨和温度变化而变化的时间规律，以及植被NDVI因采矿等人类活动而变化的时间规律。结果表明：大柳塔地区植被NDVI对降雨和温度变化的响应滞后时间为1个月，与最近2个月降雨和的相关性最高；开采对植被NDVI具有负面影响，植被NDVI响应滞后时间短于1个月。

关键词：大柳塔矿区；气候变化；人类活动；时间特征

煤炭开采在促进区域经济飞速增长的同时也严重破坏了矿区生态环境。矿产资源的开发和利用不可避免地占用和破坏了大量土地和植被,并由此引发一系列生态环境问题[1]。根据相关研究,煤炭开发活动对矿区生态环境的主要影响包括地表沉陷、土地荒漠化、土壤质量下降、地表及地下水污染、植被破坏、生态系统退化、生物多样性丧失、景观破坏以及农作物减产等多个方面[2]。全国矿山占用的土地面积为581万hm2,因煤炭开采导致的大面积塌陷达180余处,累计破坏土地面积为11.5万hm2,并且以200 hm2·a-1的速度逐年递增,占全国土地破坏面积的10%[3];2005年矿井水排放量约为45.4亿t,重复利用率不足50%,2007年矿井水排放量达46亿t,2010年矿井水排放量约为52亿t[4];我国因采矿直接破坏的森林面积累计达106万hm2,破坏草地面积为26.3万hm2[5]。植被是生态完整性损失类的主要指标(植被破坏面积、植被覆盖度、生物丰度、自然生产力等),可表征矿区土地退化状况[6],故笔者将着手于生态环境问题中的植被影响变化方面的研究。目前,国内外有关植被动态变化方面的研究主要集中在植被覆盖度与气候因子的关系[7-10]、年际气候变化和人类活动对大范围归一化差分植被指数(NDVI,INDV)的影响[11-13],而开采活动对植被变化的影响研究主要停留在定性描述上[14],有必要对小范围区域月际尺度下植被对气候因子和开采活动的响应规律进行定量分析。因此,笔者依据多源遥感影像开展矿区植被动态研究,探索矿区气候变化与人为活动影响下的植被变化规律,旨在为生态环境脆弱区煤炭资源开发过程中生态地质环境的保护与修复提供借鉴。

1研究区概况

分矿区和工作面2种尺度进行研究。矿区尺度所选区域为大柳塔矿区,该矿区位于陕西省榆林市神木县大柳塔镇以及乌兰木伦河以西中鸡镇境内,地理坐标为北纬39°27′～110°15′,东经110°05′～110°20′[15],地处毛乌素沙地边缘,北部为风沙滩地,南部和西部为典型盖沙黄土丘陵地貌区,地势北高南低,中间高东西低,海拔为1 000～1 250 m,属中温带半干旱大陆性季风气候区,干旱少雨多风沙。区内土地类型为沙地、乔灌木林地、天然草地和裸岩地等。以沙蒿、沙柳和柠条为代表的沙生植被组合主要生长于半固定沙地、固定沙地和沙地沙丘间低地[16]。近年来,煤炭开采极大地促进了大柳塔地区的经济发展,然而,该地区为半干旱半沙漠的高原地区,气候干燥,地下水贫乏,植被稀少,沙漠化严重,煤炭工业的大规模开发促使本来就脆弱的生态环境更加恶化。

工作面尺度选取大柳塔矿区22615和22103综采工作面,基本情况分别为:(1)大柳塔矿区22615综采工作面，倾向长度为339.2 m,走向长度为1 457.5 m,煤层倾角为1～3°,上覆基岩厚度在9.2～12.5 m之间,煤层平均厚度为5.15 m,选择2011年9月工作面推进位置正上方100 m×100 m区域为研究区域。(2)大柳塔矿区22103综采工作面，倾向长度为310.6 m,走向长度为1 560.4 m,煤层倾角为1～3°,煤层厚度为3.81～4.21 m,工作面设计采高为3.75 m,选择2011年7月工作面推进位置正上方100 m×100 m区域为研究区域。

2材料与方法

2.1数据来源与预处理

MODIS影像为美国国家航空航天局(NASA)提供的中、低分辨率遥感数据,空间分辨率为250、500和1 000 m,是当前世界上新一代“图谱合一”的光学遥感仪器,其高时间分辨率特性为动态监测提供基础,其产品分5个等级44种。该文所用MODIS NDVI为250 m分辨率的16 d合成数据,产品号为MOD13Q1,为陆地2级标准数据产品,已经过定标和定位,格式为国际标准EOS-HDF格式。选取时段为2005—2012年,平均每年23景。其中,2005年为24景,2008和2010年各为22景,其余年份均为每年23景,共183景数据,获取时间间隔为16 d。用于整个大柳塔地区NDVI时间序列的计算与统计。

TM影像为Landsat卫星搭载的中、高分辨率卫星遥感数据,空间分辨率为30 m,重访周期为16 d,选取云量较少、图像清晰、成像时间相近的数据,经辐射定标、大气校正和几何校正等预处理,根据式(1)计算INDV值:

(1)

式(1)中,RNI为近红外波段反射率;R为红波段反射率;INDV范围为-1～1。但由于云雨等天气因素的干扰很难保证每月1景有效的TM观测数据,采用ESTARFM数据融合算法补全缺失的数据,经验证,补全影像与实际影像相似度最高可达95%[17-19],其分辨率较高，可用于各工作面上方研究区域的INDV时间序列的计算与统计。

气象数据来源于中国气象科学数据共享服务网提供的2005—2012年大柳塔区域逐月气温和降水数据。工作面矢量图则利用矿区已有的井上下对照图得到。

2.2分析方法

2.2.1趋势和周期分析

采用基于Loess算法的时间序列季节和趋势成分分解法(season and trend decomposition using Loess,STL)[20]分析降雨、温度和NDVI的时间序列数据,目的在于分析降雨、温度序列数据和NDVI序列数据的变化趋势以及周期性,观察降雨、温度等气候因子与NDVI变化趋势是否具有一致性。STL分解法的数学模型表达式为

yt=Tt+Pt+t。

(2)

式(2)中,yt为t时刻观测值；Tt为反映该时间序列多年发展趋势的趋势项;Pt为反映时间序列要素多年发展周期性的周期项;t为随机波动项,反映了时间序列要素的随机特征,一般为高斯白噪声,即符合标准正态分布。STL方法是一个递归的过程,每次递归要分别进行3次Loess和滑动平均过程。Loess过程是基于距离越近、相关性越强的假设,赋予不同位置的点不同权重并进行局部加权回归。该过程需要选定局部回归的窗口长度(u)、回归方程阶数(d)和权重函数(W),常采用立方权重函数:

(3)

具体各项的详细检验与提取方法见文献[21-22]。

2.2.2相关性和滞后性分析

生长期内NDVI对气候水、热因素的变化响应呈现出不同程度的滞后现象[23],为进一步探明植被生长对水、热因素响应的滞后特征,采用时滞互相关分析计算生长期内植被生长对月平均降雨量、月积雨、月均温和月积温的滞后期[24]。首先,假设生长期内水、热因素(月平均降雨量、月积雨、月均温和月积温)xt与生长期内月平均NDVI的时间序列yt在时滞为k时彼此相关,互相关系数为

rk(x,y)=ck(x,y)/δxδy+k。

(4)

式(4)中,ck(x,y)为样本的协方差,表示样本的均方差,其计算公式为

(5)

(6)

(7)

式(5)～(7)中,

(8)

(9)

2.2.3残差分析和改进后的残差分析

基于EVANS等[25]和GEERKEN等[26]提出的残差分析法进行人类活动和自然影响的分离,该方法是对每个栅格像元的NDVI与气候指标做回归分析,从而得到每个像元NDVI的预测值,该预测值可以视为气候因子对NDVI的影响,然后利用遥感观测的NDVI真实值减去NDVI预测值,就得到人类活动对NDVI的影响,从而剥离开气候变化和人类活动分别对植被变化的影响[13],其表达式为

ε=INDV，观测-INDV，预测。

(10)

式(10)中,ε为NDVI中人类活动的贡献部分,即NDVI残差值,ε>0表示人类活动产生正面影响,ε<0表示人类活动产生负面影响,ε≈0表示人类活动影响微弱。

其中,气候指标选取采用2.2.2节所述方法计算所得的最佳影响指标:最近2个月降雨和及上月平均温度。由于研究区域较小,可以认为研究区域内部气候条件完全相同,逐个建立工作面上方研究区内像元NDVI平均值与最近2个月降雨和及上月平均温度2个气候指标的多元线性回归方程:

Zt=axt+byt+c。

(11)

式(11)中,Zt为t时间研究区域内NDVI平均值,即NDVI预测值;xt和yt为该区域的气候指标;a、b和c为回归方程系数。

但该方法对观测值的要求较高,观测值数据的质量直接影响了最终结果的质量,如观测当天的薄云等影像采集质量问题有可能导致完全错误的结果,因此对该方法稍作改动,改进后表达式如下。

在开采影响区域:

ε采区=INDV，观测-INDV，预测=ε系统误差+ε人为影响;

(12)

而在非开采区域:

ε非采区=INDV，观测-INDV，预测=ε系统误差。

(13)

从而可以得到人类活动的贡献部分(ε人为影响):

ε人为影响=ε采区-ε非采区。

(14)

为减少生长环境不同所带来的影响,采区与非采区应选取相同或相近的立地条件。立地条件指植被生长的自然条件,包括土壤、地形、降雨、植被覆盖率和乔灌草比等,植被立地条件的好坏直接影响着植被的生长方式。

非采区的选取参照文献[27],主要采用地形、植被覆盖率和乔灌草比3个影响因素,结合增加及衰退演变方式,将荒地演变及稳定性划分为16种立地条件,选取与各工作面相同或相近立地条件的非采区作为对照。

3结果与分析

3.12005—2012年大柳塔地区月平均降雨和温度变化

由图1可见,大柳塔地区2005—2012年月平均降雨波动较明显,每年会出现2次降雨高峰,春季之后,降雨量逐渐增加,第1次高峰在4—5月期间,紧接着在6月中旬该地区迎来干旱期,在8或9月出现另1个降雨高峰,降雨量达到一年中的最大值。

分解结果图(图2)包括3个部分,季节分图反映该时间序列的周期信息,趋势分图为长期趋势,即该时间序列的发展方向,随机分图为噪声部分。图2的3个分图中,最右侧均有高度不等的细长矩形条,该长条是反映各分图中指标值大小的一个可视化度量,它们对应的数值大小是一致的,也就是说,数据分图中很小的长条对应的数值大小幅度与趋势分图中最大坐标差是一致的。可以看出,该地区降雨量在2005—2010年平稳上升,2011年突然降低,2012年回升,但总体变化很小,可以得出该地区近年来降雨量除2011年以外呈小幅上升趋势。

图1　大柳塔地区降雨时序

每个分图中最右侧的细长矩形条是反映各

由图3可见,大柳塔地区2005—2012年月平均温度波动较明显,12或1月达到年温度最小值,温度均匀上升,在每年6、7或8月达到年温度最大值。

2005—2012年期间,最高温度为23.8 ℃,最低温度为-13.7 ℃。由STL分解图(图4)可以看出,该地区在2005—2012年期间温度均值基本保持稳定,呈缓慢下降趋势,波动不大,且无异常变化年。可以得出,该地区2005—2012年温度变化趋势为缓慢下降,但下降幅度小。

图3　大柳塔地区温度时序

每个分图中最右侧的细长矩形条是反映各

3.22005—2012年大柳塔地区月均NDVI变化

由图5～6可见,大柳塔地区NDVI均值呈缓慢上升趋势,部分年份NDVI值在1和12月下降幅度较大,可以解释为该地区大面积降雪所致,其变化特点与温度和降雨的月均变化总体趋势一致,其峰值和低谷值出现时间也一致。

NDVI均值对每年6月的降雨低谷有所反应,出现图6中所示NDVI上升过程中的“开叉口”,2010和2012年没有出现反应,可以解释为温度的上升使该地区更加适宜植被生长,且此时正处于植物快速生长期,因此NDVI持续上升,并未对短暂的干旱期做出反应。

图5　大柳塔地区NDVI时序

每个分图中最右侧的细长矩形条是反映各

3.3NDVI对气候变化响应的滞后特征

由表1～2可以看出,该地NDVI对水、热响应的滞后期约为1个月,对水、热响应的最敏感因素为上月平均温度和最近2个月降雨和;就对水、热因素响应的敏感度差异而言,NDVI对温度的响应更加敏感,但与对降雨的响应敏感程度相差不大,原因是该地区属于西北半干旱高原气候,同时缺乏水、热2种资源,因此2种气候因素的升高都会导致植物光合作用效率提高,从而导致NDVI升高。

3.4NDVI对人类活动影响的滞后特征

气候变化和人类活动是导致植被变化的2大驱动力,系统评价两者在植被生长变化中的相对作用对于深入理解矿区植被生长变化驱动力作用机理以及开展矿区重建工作具有重要意义[28]。由于该地区较贫瘠,导致人类活动较少,主要的人类活动即为煤炭开采活动。基于改进后的残差分析法进行人类活动和自然影响的分离,剥离开气候变化和人类活动分别对植被NDVI变化的影响。2个工作面上方研究区内像元NDVI平均值与最近2个月降雨和及上月平均温度2个气候指标的多元线性回归系数以及相关系数见表3。可以看出,各回归方程均通过置信度95%的显著性检验。

表1月平均降雨、月平均温度与NDVI的相关系数

Table 1Correlation coefficients of vegetation NDVI with monthly mean rainfall, and monthly mean temperature

滞后时间月平均降雨月平均温度0个月 0.7507 0.80031个月0.76970.91492个月0.52200.79143个月0.15680.46454个月0.18300.00235个月-0.4480-0.46106个月-0.5904-0.8080

滞后时间为0个月对应的是当月平均降雨(或当月平均温度),滞后时间为1个月对应的是上月平均降雨(或上月平均温度),其他滞后时间依此类推。

表2降雨和、温度和与NDVI的相关系数

Table 2The correlation coefficient between NDVI and cumulative rainfall, cumulative temperature

时间降雨和温度和最近2个月0.88060.6958最近3个月0.87780.6785最近4个月0.78910.6172最近5个月0.65170.5080最近6个月0.47400.3570最近7个月0.27570.1741

最近2个月降雨和(或温度和)指当月平均降雨量(或温度)与上月平均降雨量(或温度)之和,最近3、4、5、6和7个月降雨和(或温度和)依此类推。

表32个工作面NDVI平均值与2个最佳影响指标的多元线性回归系数和相关系数

Table 3Regression coefficients, correlation coefficients andFvalue of different mining face

工作面名称系数a系数b系数c复相关系数RF值226150.0012610.0034810.1144100.840.00001020221030.0024020.0036300.1298860.880.00000013

开采当年NDVI均值变化和人类影响残差随时间变化见图7。对比可以得出,图7(b)和(f)以及(c)和(g)具有很强的相似性,从而验证笔者观点:单纯的观测值与预测值的残差并不能代表人类活动的影响,该残差主要由获取影像的质量所决定,可以归属为系统误差,通过选取相同立地条件的研究区域,计算该研究区域的残差,并计算2种残差之差来减小该系统误差。由图7(d)和(h)可以看出,开采到该地区前1个月该地区人为影响即残差(ε人为影响)开始增加,该地区开采当月,残差上升至最大值,之后便逐渐下降,导致该现象的原因为开采的超前影响,即在开采工作面推进过程中,工作面前方的地表受采空区的影响而下沉,这种现象被称为超前影响。在开采工作面推进到实验区域的前1个月,由于超前影响的存在,研究区域已经开始受到影响,人类影响值开始增加。可以得出:开采对植被NDVI的扰动为导致其下降的负面影响,植被NDVI下降响应的滞后时间短于1个月,即开采当月人类活动对植被的影响达到最大值。

图7　不同工作面上方残差变化

4结论

对2005—2012年降雨和温度等气候因子和NDVI时序数据进行分解分析,并分别在大柳塔矿区尺度和工作面尺度上建立气候因子和NDVI的相关性,对残差分析方法进行改进并用于分离NDVI变化中的人为影响和气候影响。可以得出,2005—2012年该区域降雨和温度因子总体变化较小,其中降雨呈微弱上升趋势,温度呈微弱下降趋势;在自然影响和开采活动的共同影响下,2005—2012年大柳塔地区植被NDVI呈缓慢上升趋势,其中气候因素为正面影响,开采活动为负面影响。月均NDVI与最近2个月降雨和、上月平均温度呈显著正相关,相关系数分别为0.880 6和0.914 9,得出月均NDVI对降雨和温度因素的响应滞后时间约为1个月。工作面尺度上的月际NDVI受开采影响,响应滞后时间短于1个月,即植被NDVI在开采当月受到的人为扰动最强烈。

参考文献：

[1]孙庆先,胡振琪.中国矿业的环境影响及可持续发展[J].中国矿业,2003,12(7):23-26

[2]范英宏,陆兆华,程建龙,等.中国煤矿区主要生态环境问题及生态重建技术[J].生态学报,2003,23(10):2144-2152.

[3]侯晓丽.废弃煤矿土地复垦研究[D].重庆:西南大学,2010.

[4]李泾妮.煤炭开采环境影响研究[D].太原:太原理工大学,2011.

[5]耿殿明,姜福兴.我国煤炭矿区生态环境问题分析[J].中国煤炭,2002,28(7):21-24.

[6]李海东,沈渭寿,司万童,等.中国矿区土地退化因素调查:概念、类型与方法[J].生态与农村环境学报,2015,31(4):445-451.

[7]韦振锋,王德光,张翀,等.近12年陕甘宁黄土高原区植被物候时空变化特征[J].生态与农村环境学报,2014,30(4):423-429.

[8]刘绿柳,肖风劲.黄河流域植被NDVI与温度、降水关系的时空变化[J].生态学杂志,2006,26(5):477-481，502.

[9]杨建平,丁永建,陈仁升.长江黄河源区高寒植被变化的NDVI记录[J].地理学报,2005,60(3):467-478.

[10]ROERINK G J,MENENTI M.Assessment of Climate Impact on Vegetation Dynamics by Using Remote Sensing[J].Physics and Chemistry of the Earth,2003,28(3),103-109.

[11]许端阳,康相武,刘志丽,等.气候变化和人类活动在鄂尔多斯地区沙漠化过程中的相对作用研究[J].中国科学(地球科学),2009,39(4):516-528.

[12]周伟,刚成诚,李建龙,等.1982—2010年中国草地覆盖度的时空动态及其对气候变化的响应[J].地理学报,2014,69(1):15-30.

[13]刘斌,孙艳玲,王中良,等.华北地区植被覆盖变化及其影响因子的相对作用分析[J].自然资源学报,2015,30(1):12-23.

[14]TOOM1K A,LIBIK V.Oil Stale Mining and Processing Impact on Landscapes in Northeast Estonia[J].Landscape and Urban Planning,1998,41(3):285-292.

[15]赵红梅,张发旺,宋亚新,等.大柳塔采煤塌陷区土壤含水量的空间变异特征分析[J].地球信息科学学报,2010,12(6):753-760.

[16]徐友宁,李智佩,陈社斌,等.大柳塔煤矿采煤塌陷对土地沙漠化进程的影响[J].中国地质,2008,35(1):157-162.

[17]赵永光,黄波,汪超亮.基于超分辨率重建的多时相MODIS与Landsat反射率融合方法[J].遥感学报,2013,17(3):590-608.

[18]王昆,张丽,王志勇,等.基于半方差函数的STARFM改进模型[J].测绘科学,2013,24(3):140-142,92.

[19]张开,周红敏,王锦地,等.融合Landsat ETM+和MODIS数据估算高时空分辨率地表短波反照率[J].遥感学报,2014,18(3):497-517.

[20]中国人民银行调查统计局.时间序列X-12-ARIMA季节调整:原理与方法[M].北京:中国金融出版社,2006:133-149.

[21]黄显峰,邵东国,阳书敏.降雨时间序列分解预测模型及应用[J].中国农村水利水电,2007,28(9):6-8.

[22]王妍.中长期降雨预测方法研究[D].合肥:合肥工业大学,2003.

[23]FABRICANTE I,OESTERHELD M,PARUELO J M.Annual and Seasonal Variation of NDVI Explained by Current and Previous Precipitation Across Northern Pata-Gonia[J].Journal of Arid Environments,2009,73(8):745-753.

[24]李学梅,任志远,张翀,等.重庆市NDVI对水热条件变化的响应及其空间特征[J].水土保持通报,2013,33(4):166-169,175,317.

[25]EVANS J,GEERKEN R.Discrimination Between Climate and Human-Induced Dryland Degradation[J].Journal of Arid Environments,2004,57(3):535-554.

[26]GEERKEN R,IIAIWI M.Assessment of Rangeland Degradation and Development of a Strategy for Rehabilitation[J].Remote Sensing of Environment,2004,90(3):490-504.

[27]申艳琴.半干旱区煤炭开采对植被扰动规律的研究[D].徐州:中国矿业大学,2014.

[28]许端阳,李春蕾,庄大方,等.气候变化和人类活动在沙漠化过程中相对作用评价综述[J].地理学报,2011,66(1):68-76.

(责任编辑： 李祥敏)

Temporal Characteristics of Response of Vegetation to Disturbance of Mining and Subsidence in Semi-Arid Coal Mining Area.

TAOWen-kuang1,LEIShao-gang2

(1.School of Environment Science and Spatial Informatics, China University of Mining and Technology, Xuzhou 221116, China;2.Institute of Land Resources, China University of Mining and Technology, Xuzhou 221116, China)

Abstract:STL(season and trend decomposition using Loess) time sequence decomposition(locally weighted regression scatter smoothing)method, time lag correlation analysis and modified residual analysis were performed of the MODIS images, the mean monthly NDVI(normalized difference vegetation index) time sequence data of TM images and the meteorological data of Daliuta of the Shendong Mining Area of the time period of 2005-2012, to explore laws of the temporal variation of vegetation NDVI in Daliuta with changes in precipitation and temperature during the period from 2005 to 2012 as well as laws of temporal variation of vegetation NDVI in response to human activities, especially mining activities. Results show that the vegetation NDVI in Daliuta responded one month behind the changes in precipitation and temperature, and was the most closely related to rainfall of the past two months; coal mining did have some negative impacts on vegetation NDVI, which responded a bit less than one month later.

Key words:Daliuta mining area;climate change;human activity;temporal characteristics

作者简介：陶文旷(1993—),男,安徽马鞍山人,硕士生,主要从事环境遥感方面的研究。E-mail: tao129tao@163.com

DOI：10.11934/j.issn.1673-4831.2016.02.005

中图分类号：X87

文献标志码：A

文章编号：1673-4831(2016)02-0200-07

通信作者①E-mail: lsgang@126.com

基金项目：教育部新世纪优秀人才支持计划(NCET-12-0956);国家科技基础性工作专项(2014FY110800)

收稿日期：2015-11-09