基于改进遥感生态指数的山西省及煤炭矿区生态环境分析

李 晶 李生财 郭 伟 闫星光 张 瑞

(中国矿业大学(北京)地球科学与测绘工程学院,北京 100083)

随着人类社会的不断发展,资源开发、城市扩张等已经使得诸多地区和全球性生态环境问题日益突出,如资源短缺、大气污染、土壤污染、水土流失等。尤其对于黄河流域生态脆弱区,生态环境保护工作迫切性持续增强[1]。准确监测区域生态环境质量状况及变化趋势,对于加强生态环境治理和支撑区域规划、应对全球环境变化以及实现区域可持续发展目标等具有重要意义。

山西省位于黄河流域黄土高原区,区域内地形起伏大,加之降水比较集中,地面植被覆盖度相对较低,境内侵蚀严重,径流多泥沙,水土流失面积达10.8万km2,为全国水土流失最为严重的区域之一[2]。近年来,随着城市扩张、工矿开发、气候变化等的持续扰动,进一步加剧了山西省的生态环境问题。山西省作为我国重要的煤炭基地,长时间、大范围、高强度的煤炭资源开采对地表最直接的影响是挖损、沉陷或压占,在局部范围内会加速土壤侵蚀,加剧植被破坏、生态系统受损等[3-6]。为推进区域生态保护修复和促进人地协调发展,亟需在山西省及其煤炭国家规划矿区(以下简称“煤炭矿区”)开展生态环境监测和系统性研究。

近年来,遥感技术已在生态环境监测领域得到广泛应用。徐涵秋[7]提出一种完全基于遥感信息的生态指数(Remote Sensing Ecology Index,RSEI),并对福州主城区进行了生态环境质量评价。遥感生态指数RSEI的优势之一在于不需要人为确定各生态指标的权重,通过主成分分析(Principal Component Analysis,PCA),自动确定各个指标对各主成分的贡献度,计算结果具有很强的客观性和稳定性[8]。目前,针对煤炭矿区,一些学者应用该指数对特定矿区进行了生态环境质量分析[9-13],在矿区生态环境监测应用方面取得了一定的进展,但少有考虑特定研究区的具体地质环境和人为干扰特征等因素。李粉玲等[14]采用遥感生态指数RSEI对黄土高原沟壑区的生态环境质量进行了评价,类似研究未充分顾及不同区域生态环境的地域特点。对此,一些学者提出了多种改进的遥感生态指数。王杰等[15]针对特定的干旱区,基于遥感生态指数RSEI,将盐度和土地退化度指标共同参与进行主成分分析,构建了干旱遥感生态指数ARSEI。宋美杰等[16]提出利用贡献度对主成分分量进行相加求和改进RSEI,分析了锡林郭勒草原的生态状况。杨羽佳等[17]基于遥感生态指数RSEI,加入人类活动强度指标(IPOI),结合AHP层次分析,分析了苏州市的生态环境状况。朱冬雨等[18]考虑矿区不同损毁特征用地(采掘场、沉陷区、排土场等)具有不同的生态影响范围,提出了基于遥感生态指数和移动窗口评价单元的改进生态指数MW-RSEI,并将其应用于矿区生态环境评价。

黄土高原特有的地理环境决定了其易于水土流失的现状,矿山开采活动[19]有可能加重这种趋势。本研究在绿度、湿度、干度、热度因子的基础上,提出了一种结合土壤侵蚀因子的改进型遥感生态指数(Modified Remote Sensing Ecology Index,MRSEI)。基于GEE云计算平台[20],采用Landsat TM/OLI数据,对山西省煤炭矿区进行生态环境分析,为实现区域可持续发展和制定资源环境协调开采政策及相关治理规划提供科学依据,也为矿山生态智能开采等提供区域数据支撑。

1 研究区概况

黄河流域是我国重要的生态屏障和能源走廊,据统计,我国14个大型煤炭基地有9个分布在黄河流域,煤炭开采活动与黄河流域矿区周边生态变化息息相关[21]。山西省是黄河流经的重要区域,地处黄河中游,位于34°34′N～40°43′N和110°14′E～114°33′E区域,在太行山与黄河北干流域峡谷之间,地处华北西部的黄土高原东翼,总面积为15.68万km2。同时,山西省是我国重要的煤炭生产基地,这一区域原生生态脆弱,而煤炭开采在一定程度上不可避免地加剧了所在区域的生态环境恶化过程。《土地复垦规定》颁布实施后,尤其是党的十八大以来,矿山企业积极履行土地复垦责任,但仍然存在大量历史遗留矿山废弃地和部分新增矿山损毁土地尚未治理。

2 数据及来源

本研究所用数据包括遥感影像数据、DEM数据、气象数据、土壤数据、土地利用数据等。

2.1 遥感影像数据

Landsat TM/ETM+/OLI数据来源于美国地质调查局(United States Geological Survey,USGS),在GEE平台集成,空间分辨率为30 m,使用的是经过辐射校正和大气校正后的地表反射率(Surface Reflectance,SR)数据,已完成几何精校正[22]。

2.2 土壤侵蚀计算数据

DEM数据为ASTGTM DEM的数字高程模型,来源于地理空间数据云(http:∥www.gscloud.cn/search);气象数据为中国国家级地面气象站基本气象要素日值数据集(V3.0),来源于国家科学气象信息中心(http:∥data.cma.cn/);土壤可蚀性因子(K)栅格数据来源于国家科技资源共享服务平台——国家地球系统科学数据中心黄土高原分中心(http:∥loess.geodata.cn);植被覆盖数据在GEE平台集成,为年度NDVI最大合成数据;土地利用数据为全球地理信息公共产品Globe Land 30数据(http:∥globeland-30.org/)。除了气象数据外,其他4类数据的空间分辨率均为30 m。

3 研究方法

3.1 RUSLE模型

RUSLE是经过修正的通用土壤流失方程,自1997年提出以来被广泛采用,其模型结构相对简单[23]。本研究采用该模型作为研究区土壤侵蚀因子的定量计算模型,公式为

式中,A为年平均土壤侵蚀模数,t/(km2·a);R为降雨侵蚀力因子,MJ·mm/(hm2·h·a);K为土壤可蚀性因子,t/[(MJ·mm/h·hm2·a)·hm];L和S分别为坡长和坡度因子;C为地表植被覆盖和管理因子;P为水土保持措施因子。

根据模型定义可知:土壤侵蚀模数取值越大,所在区域的土壤侵蚀强度就越大;土壤侵蚀模数取值越小,所在区域的土壤侵蚀强度就越小。

3.1.1 坡长坡度因子(LS)

坡度和坡长是地形地貌主要的地理特征因素。通常大尺度区域的坡度坡长因子利用DEM数据进行提取。坡度因子S的计算公式为

式中,S为坡度因子;θ为坡度,(°)。

坡长因子L的计算公式为

式中,L、m分别为坡长因子、坡长指数,m一般在黄土高原取值为0.5[24];l为坡长,m。

3.1.2 降雨侵蚀力因子(R)

R因子是降水对土壤侵蚀影响的反映指标,利用年平均降水量估算区域的降雨侵蚀力[25]。计算公式为

式中,pi为第i月的降雨量,mm;p′为年平均降雨量,mm;n为有侵蚀性降雨的月份数;F为模型修正参数;R为年平均降雨侵蚀力,MJ·mm/(hm2·h·a);α和β为模型参数,根据已有成果[25],取α=0.183 3,β=1.995 7。

根据山西省和周边共109个站点的日降水量数据,使用Arc GIS10.3软件进行克里金插值,得到研究区2000、2010和2020年的R值分布图。

3.1.3 地表植被覆盖和管理因子(C)

水土流失的抑制程度受植被覆盖度或田间管理措施的影响。NDVI指数与地表植被覆盖度密切相关,可以在单位像元内综合反映植被覆盖形式、类型和生长状况。土壤侵蚀与植被覆盖度之间存在显著的负指数关系,其拐点约为78.3%。蔡崇法等[26]提出了一种新的计算方法,通过观察自然降雨和人工模拟降雨条件下试验区的植被覆盖率来估算因子C。计算公式为

式中,C为植被覆盖与作物管理因子;f为植被覆盖度。

植被覆盖度f的计算公式为

式中,NDVI为年度最大合成NDVI值;NDVIsoil为裸土或无植被覆盖区域的NDVI值;NDVIveg为完全被植被所覆盖像元的NDVI值。

3.1.4 水土保持措施因子(P)

水土保持措施是防治土壤侵蚀的有效方式,通常需要实地考察来确定不同土地利用类型的P因子值,一般取值为0～1。P因子可以表示水土保持措施对土壤侵蚀的抑制作用,所以P值越小,说明水土保持措施对土壤侵蚀的抑制作用越突出[27]。目前,国内在土壤侵蚀预报研究中尚未提出对P因子的赋值标准,本研究参考现有成果[24],利用不同的土地利用类型估算了P值,结果见表1。

表1 研究区不同土地利用类型的P值Table 1 P values for different land use types in the study area

通过以上对RUSLE模型中各因子的计算,得出山西省土壤侵蚀所需的各个因子分布图。在此基础上,利用土壤侵蚀模型,计算出土壤侵蚀模数。

3.2 MRSEI模型构建

3.2.1 绿度指标

归一化植被指数(NDVI)是当前应用最为普遍的植被指数,能够有效地反映植物生长状况。NDVI取值区间为[-1,1],值越大,表明植被覆盖度越高[28]。因此,本研究使用归一化植被指数(NDVI)作为绿度指标参与计算,公式为

式中,ρnir为近红外波段的反射率;ρred为红波段的反射率。

3.2.2 湿度指标

湿度指标WET与植被和土壤的含水量紧密相连[28],取值区间为[-1,1],值越大,表明湿度越高。湿度指标可以由K-T变化中的WET分量来代表,Landsat不同传感器影像数据WET分量的计算公式为

式中,ρblue、ρgreen、ρred、ρnir、ρswir1分别为蓝波段、绿波段、红波段、近红外波段和短波红外1的反射率。

3.2.3 热度指标

热度指标LST用地表温度来表示,取值越大,表明地表温度越高。本研究通过对亮度温度进行修正处理来获取地表温度,计算公式为

式中,LTIR为热红外波段的辐射定标值;gain为增益;bias为偏置;T为亮度温度;K1和K2为定标系数;LST为热度指标;λ为热红外波段的中心波长;ε为地表比辐射率。

地表比辐射率采用SOBRINO等[29]提出的NDVI阈值法来计算,公式为

式中,ε为地表比辐射率;f为植被覆盖度。

3.2.4 干度指标

建筑物和裸土会造成土地干化,严重影响流域的生态环境。NDBSI取值区间为[-1,1],值越大,表明干化程度越高。本研究采用建筑物指数(IBI)和裸土指数(SI)的平均值来构建干度指标(NDBSI)。计算公式分别为

式中,ρblue、ρgreen、ρred分别为蓝波段、绿波段、红波段的反射率;ρnir、ρswir1分别为近红外波段、短波红外1的反射率。

3.2.5 MRSEI指数构建

常见的生态因子集成方式有简单相加、相乘、分组求均值后再相加、加权求和等,但各生态因子间的相关性及权重设定的主观性会影响该类整合方法的稳定性[7]。主成分分析(PCA)是一种通过正交线性变换将多个变量压缩出少数主要变量的技术,可以将多维数据的信息集中到少数的变量中,这些少数的变量都具有一定的特征信息。本研究构建MRSEI指数时采用PCA法对5个生态因子进行集成。为避免量纲不统一造成权重失衡,在进行PCA运算之前,将各生态因子归一化为[0,1]区间上的数值,计算公式为

式中,Yi表示归一化后的生态因子值;Ii为生态因子的原始像元值;Imax为生态因子的全局最大像元值;Imin为生态因子的全局最小像元值。

主成分分量的计算公式为

式中,PCi为第i个像元获取的主成分分量;NDVI为归一化植被指数;WET为湿度;NDBSI为建筑和裸土指数;LST为地表温度;Erosion为土壤侵蚀模型的归一化因子。

本研究提取第一主成分分量来构建改进型遥感生态指数MRSEI。在主成分分析结果中,只有当对生态起正面作用的绿度(NDVI)和湿度(WET)指标在PC1的载荷为负值时,才用1减去PC1,公式为

式中,PC1为第一主成分分量。

为便于度量和比较,可进一步将MRSEI归一化到[0,1]区间。当MRSEI值越靠近1时,说明生态环境质量越好;MRSEI值越靠近0时,说明生态环境质量越差。

4 结果分析

4.1 模型检验

山西省3个年份5项指标的主成分PC1分析结果见表2。由表2可知:2000、2010和2020年第一主成分特征值的贡献率分别为84.57%、71.82%和75.26%。这说明PC1集中了5个指标的大部分特征。2000年和2010年,在PC1中,代表绿度的NDVI和代表湿度的WET的主成分载荷为正值,表明它们在生态环境方面发挥了正向作用,而代表热度、干度和土壤侵蚀的LST、NDBSI和Erosion的主成分载荷为负值,表明它们对生态环境产生了负面影响,与现实相符。2020年,主成分对应系数情况与2000年和2010年相反,因此,需要进行“1-PC1”处理后,结果NDVI和WET均为正值,而NDBSI、LST和Erosion均为负值。

表2 主成分分析结果Table 2 Results of principal component analysis results

4.2 MRSEI与RSEI对比分析

为了进一步验证MRSEI的适用性,选取山西省吕梁市离柳矿区的部分区域进行分析。2020年采用RSEI与MRSEI两种指数在煤炭矿区生态环境的监测效果对比如图1所示。图1局部细节和关键区域差异表明,相较于RSEI,MRSEI显示的空间信息更为丰富,与地面覆盖真实情况较为匹配。A1、A2和A3区域分别是东江煤业、临县胜利煤矿和离石西山晋邦德煤矿所在区域,靠近黄土高原东侧,矿山开采叠加地质环境使得土壤侵蚀较为严重。根据RUSLE模型计算,将土壤侵蚀强度分为微度、轻度、中度、强、极强和剧烈5个等级,A1、A2和A3区域侵蚀“强”及以上的面积占比分别为76.95%、80.15%和67.49%。与RSEI相比,MRSEI可以完整地展示土壤侵蚀的分布形态,土壤侵蚀等级越高的区域,MRSEI值越低,越可以较为精准地表达煤炭矿区地表的真实状态。

图1 2020年区域RSEI与MRSEI的结果对比Fig.1 Comparison of regional RSEI and MRSEI results in 2020

为了进一步分析RSEI与MRSEI应用效果的差异,采用信息熵进一步定量化对比分析局部细节信息。信息熵被用来度量图像的信息量,当熵值增大时,表示图像中复杂程度较高,所含的信息量越大。对选取的区域利用MATLAB软件计算图像的信息熵,RSEI与MRSEI的信息熵分别为1.38、1.563;对整个山西省而言,RSEI与MRSEI的信息熵分别为1.698、1.719。结果表明,MRSEI与RSEI相比,前者图像的复杂程度更高,所包含的信息量更大,更能反映煤炭矿区的真实生态状况。综上分析,MRSEI可以有效集成新加入的Erosion因子,增强了其在生态脆弱区生态环境分析中的适用性。

4.3 模型相关性检验

通过计算各指标间的相关性系数,来检验MRSEI的适用性,相关系数越接近1,表明改进的生态指数的综合代表程度越高,适宜性越好。计算公式为

式中,表示平均相关度;p、q、r、s是进行相关性分析的指标;n为相关性分析指标数量;Xp、Xq、Xr和Xs分别表示各指标之间的相关性系数。

为了进一步分析RSEI的综合代表性,对MRSEI和5个生态指标进行了相关性分析,结果见表3。通过计算各指标间的相关系数,在95%置信水平P值均小于0.05,通过显著性检验。可以看出,平均相关度最大的为MRSEI,3 a的平均值为0.876,在2000年达到最大值0.895,第2位是LST,3个年份的平均值为0.841,NDVI、WET、NDBSI和Erosion的相关度较低,分别为0.804、0.799、0.689和0.362。从单个因子来看,各因子互相之间的平均相关度最高的为LST,在2000年达到最大值0.871,改进的遥感生态指数MRSEI与这5个因子的各年份平均相关度都大于0.860,其最大值比单因子最高LST的0.871高出了2.7%,明显可以看出,改进的生态指数除了综合了各因子的信息外,还比单独因子更具代表性,能更好地代表区域综合生态状况。在改进的生态指数MRSEI与各单一因子的相关分析中,NDVI和WET与MRSEI呈正相关,NDBSI、LST和Erosion与MRSEI呈负相关,进一步说明MRSEI具有较好的适用性。

表3 相关性分析结果Table 3 Results of correlation analysis

4.4 山西省及煤炭矿区生态环境静态比较

采取等差分级的方法,以0.2为间隔,将MRSEI依次分为差、较差、中等、良、优5个等级。山西省及其煤炭矿区各生态级别面积占比见图2。从时序生态级别面积变化来看,2000—2010年山西省生态环境“优良”面积占比减少了1.57%,生态环境为“差”“较差”“中等”的面积比例则相应略有增加;2010—2020年,生态环境级别为“优”“良”的面积占比增加了26.16%,生态环境“差”“较差”“中等”的面积比例相应减少;2000—2020年,生态环境为“优”“良”的面积占比由17.96%增加至42.55%,生态级别为“差”“较差”“中等”的面积占比相应减少。

图2 山西省及煤炭矿区各生态级别面积占比Fig.2 Percentage of area of each ecological level in Shanxi Province and coal mining areas

与山西省相比,2000—2020年煤炭矿区的生态环境质量呈先下降后上升的趋势,2020年的生态环境质量最好。2000年,煤炭矿区生态环境等级以“较差”和“中等”为主,面积占比为89.56%;“优”和“良”面积合计占比为10.44%。2010年,煤炭矿区生态环境等级仍以“中等”和“较差”为主,面积占比为88.06%。2020年,煤炭矿区生态等级以“中等”和“良”为主,面积合计占比为89.73%;“差”和“较差”面积占比为9.98%。与2000年相比,煤炭矿区的生态环境状况正在持续变好。与整个山西省相比,煤炭矿区的生态环境相对较差,但两者差距从2000年、2010年至2020年逐渐缩小。

2000—2020年,MRSEI值的空间分布如图3所示。从3期生态质量图对比可以看出,山西省生态质量整体上一直是东部高于中部和西部,2020年生态级别为“优”和“良”的区域主要分布在阳泉、晋中、长治和晋城市的东部。生态级别为“差”和“较差”的区域主要分布在忻州、大同、朔州和吕梁市的西侧。综合来看,生态环境质量较差的地区主要分布在山西省西部、西北部和中部,中部生态环境质量后期有所改善,西北部生态质量始终相对较差。

图3 2000年、2010年和2020年山西省生态质量变化Fig.3 Ecological quality comparison in Shanxi Province in 2000,2010 and 2020

4.5 生态环境时空变化特征

为了进一步分析山西省及其煤炭矿区生态环境变化的时空分布状况,将MRSEI的5个生态级别,采用差值法对3个年份的生态环境进行了变化监测[22],结果见图4。由图4可知:2000—2010年,生态环境变好的区域占11.65%,分布在临汾、运城市的部分区域;生态环境变差的区域占19.92%,分布在吕梁、晋城、长治、朔州市的大部分区域。2010—2020年,生态环境变好的区域占41.03%,分布在大同、晋中、阳泉、太原、长治、晋城市的大部分区域;生态环境变差的区域占5.68%,分布在忻州、吕梁、运城市的西南、西北局部区域。2000—2020年,生态环境变好的区域占35.67%,主要分布在长治、临汾、晋中和吕梁等中部区域。生态环境变差的区域占7.67%,主要分布吕梁、忻州、大同、朔州、运城的西部地区。从变化检测的结果来看,2000—2020年,总体上山西省的生态环境先变差而后得到改善。

图4 山西省生态质量变化监测Fig.4 Monitoring of ecological quality variation in Shanxi Province

从煤炭矿区来看,2000—2010年,生态环境变好的区域占15.86%,生态环境变差的区域占19.05%。2010—2020年,生态环境变好的区域占43.80%,生态环境变差的区域占5.15%。2000—2020年,生态环境变好的区域占43.17%,主要分布在霍州矿区、河保偏矿区、武夏矿区、潞安矿区和汾西矿区的中部地区。生态环境变差的区域占6.58%,主要分布平朔朔南矿区、轩岗岗县矿区、离柳矿区西部地区,尤其以离柳矿区相对严重。

5 结 论

针对黄河流域生态脆弱和强人类活动复合型水土流失易发区生态环境质量动态监测的现实需求,提出了一种改进的生态指数MRSEI,并在山西省及其煤炭矿区进行了应用。主要取得如下结论:

(1)改进的生态指数MRSEI在山西省等生态脆弱区的生态环境分析中具有较好的适用性。该指数可以有效集成绿度、湿度、干度、热度和土壤侵蚀5个生态因子的综合信息,所引入的土壤侵蚀指标Erosion使得MRSEI可以综合考虑土壤侵蚀对生态环境的影响,与现有的遥感生态指标相比可以更好地兼顾区域环境特点。

(2)2000—2020年,山西省及其煤炭矿区的生态环境质量均呈现先下降后上升的趋势,尽管煤炭矿区上升趋势相对更为明显,但其MRSEI均值始终低于山西省MRSEI均值。另一方面,生态环境质量变化空间差异特征相对明显,生态环境质量较好的区域主要分布在山西省东部和南部,较差的区域主要集中在西部地区。此外,无论是山西省还是其煤炭矿区,生态环境质量较好及以上区域面积占比均增加,而质量较差及以下的区域面积占比均下降,这与气候暖湿化变化特征及生态文明政策措施落地实施等因素密不可分。

(3)本研究提出的一种改进型遥感生态指数MRSEI,尽管可以相对客观地揭示山西省及煤炭矿区生态环境质量变化特征,但研究中也存在不足。例如:土壤侵蚀因子分析时未能获取3期的DEM数据,这是导致主成分分析时土壤侵蚀因子贡献度较低的一个原因。此外,气温、降水、地形等是导致生态环境质量变化最基本、影响最广泛的因素,矿山开发的影响局限在煤炭国家规划矿区内的有限范围,城镇化、工业化等是影响范围相对较广的人类活动,上述自然、人为因素对长时序生态环境质量变化的深刻驱动机制有待进一步研究。