祁连山南麓露天矿区高寒草地植被生态系统稳定性研究
——以青海木里-聚乎更矿区为例

陈 峰,杨利亚,冯海波,3,周建伟*,熊瑞民,朱 怡

(1.中国地质大学(武汉)环境学院,湖北 武汉 430078;2.自然资源部矿山生态效应与系统修复重点实验室,北京 100081;3.中国地质大学(武汉)地理与信息工程学院,湖北 武汉 430078;4.中国地质大学(武汉)地质调查研究院,湖北 武汉 430074;5.中国地质大学(武汉)艺术与传媒学院,湖北 武汉 430074)

青海木里位于祁连山高寒多年冻土区与草甸草原区,是青海省最大的煤矿区,也是西北地区重要的炼焦煤资源产地,其地处黄河重要支流大通河的发源地,是祁连山区水源涵养地和生态安全屏障的重要组成部分,生态地位极为重要。近年来区内大型露天煤矿开采加剧了草地生态系统退化。何芳等[1]基于高分辨率遥感影像监测和野外调查,得出木里矿区地形地貌景观和土地资源被压占与破坏严重,高寒沼泽、草甸退化,冻土与含水层结构遭到严重破坏。高寒草原矿区生态安全问题一直备受关注。因此,开展高寒草原矿区植被生态系统稳定性研究,能够深入揭示露天采矿活动及其生态修复工程对高寒草原矿区植被生态系统的影响,其研究成果可为高寒草原矿区生态保护与修复提供科学依据。

近年来,学者们对于生态系统稳定性的研究大部分集中在评价方面[2],一些学者在对生态系统稳定性概念理解的基础上,通过选取生态系统活力、抵抗力、恢复力、变异性等方面的指标对生态系统稳定性进行了评价[3-5],但由于生态系统的复杂性以及各要素衡量的困难性,目前尚没有建立标准化的计算与量化方法[6-7]。也有一些学者从景观格局动态的角度,通过选取外界干扰因子、景观本身生态学特征因子以及景观格局指数对生态系统稳定性进行评价[8]。但是由于目前景观生态学中固有的景观格局指数更注重景观几何特征、状态的描述,而缺乏动态机制与过程的研究[9-10],使得目前对于生态系统稳定性研究仍以描述性评价为主[11-13]。然而,本文的共同作者冯海波提出了基于归一化植被指数(NDVI)的时空数据序列分析方法,该方法基于植被生态系统稳定性与涨落理论的分析,以表征植被自然属性的NDVI数据为基础,通过对植被生态系统如何由稳定走向失稳的判别,提出能够表征植被生态系统稳定性的3个指数,即NDVI的均值(μ)、均方差(σ)和相对涨落值(δ=σ/μ),并通过NDVI数据的网格化处理、网格内植被稳定性指数计算与时空维度分析,揭示植被生态系统的稳定性特征[14]。该方法已被应用于呼伦贝尔草原陈旗盆地露天矿区植被生态系统稳定性的研究中,由于植被稳定性指数计算简单,可对重点区域、异常涨落区域进行动态精细化分析及稳定性识别,不涉及区位性的干扰,可广泛应用于荒漠化、石漠化、高寒草原矿区、森林以及人造景观等植被生态系统稳定性的研究。

因此,本文以祁连山南麓木里-聚乎更矿区为研究区,以2002—2022年间9期遥感影像为数据基础,采用基于归一化植被指数(normalized difference vegetation index,NDVI)的时空数据序列分析方法,综合运用GIS技术和数理统计方法,对矿业活动影响下的高寒草地植被生态系统开展时空维度动态规律的精细分析,揭示植被生态系统稳定性特征及其在矿山生命周期不同阶段的演化方向与稳定性,这些研究成果可为祁连山高寒草原矿区的资源合理开发与生态建设决策提供理论支撑与科学依据。

1 研究区概况

木里-聚乎更矿区位于青海省东北部,海西州天峻县境内,本文选取木里-聚乎更矿区及其周边区域8 988.54 hm2范围为研究区,见图1。该区域气候寒冷,昼夜温差大,6—8月份为雨季,11月份至来年5月份以降雪为主;年平均气温为-4.2 ℃,年平均降水量为477.1 mm左右,年平均蒸发量为1 049.9 mm,长冬无夏,冻土发育,为典型的高原亚寒带草原半干旱气候区[15]。该地区土壤主要为草甸沼泽土,广泛发育冻土[16];矿区植被类型分为高寒沼泽类和高寒草甸类,以耐寒旱的多年生丛生禾草、披碱草、蒿草、苔草为优势种,但植物群落种类组成贫乏、群落结构简单、植被稀疏,对人类活动的抗干扰能力较弱[17]。

图1 祁连山南麓青海木里-聚乎更矿区地理分布图Fig.1 Geographical distribution map of Qinghai Muli-Juhugeng mining area in the southern Qilian Mountains

祁连山地区气候在暖干化向暖湿化演变的过程中[18],气温与降水是影响高寒草地植被生长的主要自然因素[19-21]。根据2002—2022年区域气象数据,研究区区域降水量与NDVI均值之间具有高度的正相关性,如图2所示。

图2 研究区区域降水量与NDVI之间的相关性Fig.2 Correlation between regional precipitation and NDVI in the study area

木里-聚乎更矿区自20世纪70年代开始就有小规模开采,到2003年先后有8家企业进入矿区从事煤炭资源勘查及开发工作。2015年青海省对祁连山木里-聚乎更矿区煤炭资源开发引起的环境影响开展了渣山回填、渣山刷坡、整形覆土、种草覆绿、格构护坡等生态环境综合整治工作,但由于采坑部分地段坡体松散、草地退化,生态修复治理工程未能达到预期效果。2020年8月,青海省进一步对该矿区进行采坑回填、渣山覆绿、边坡治理、植被恢复、围栏封育等生态修复治理,选取的矿区覆绿植被种类增多[22],并通过对各类型矿区覆绿渣土进行改良[23]、土质优化、提升土壤肥分及坡体稳定性,目前矿区生态修复治理效果良好。

2 数据来源与分析方法

2.1 数据来源

本研究选取USGS官网(usgs.gov)的Landsat遥感影像数据,依据木里-聚乎更矿区开采与整治历史及云量、雪量情况,选取煤矿开采至修复治理期间地表覆盖清晰可见、云量小于10%的9期夏季Landsat遥感影像数据,见表1。

表1 Landsat遥感影像数据来源

气象数据来源于ERA5 monthly averaged data on single levels from 1979 to present数据集,空间分辨率为0.25°×0.25°,时间分辨率为逐月。

本研究使用ENVI 5.3软件对下载获取的各期次Landsat遥感影像数据进行辐射定标、大气校正、裁剪镶嵌等预处理后,进行归一化植被指数(NDVI)计算与合成。

2.2 研究方法

本研究基于前人提出的草原植被生态系统稳定性研究的新方法[9],以NDVI的均值(μ)、均方差(σ)、相对涨落值(δ)3个植被稳定性指数为依据,对研究区高寒草地植被生态系统异常涨落的区域进行识别,揭示其动态过程及稳定性。

基于NDVI的时空数据序列分析方法的流程(图3)如下:①NDVI时间序列数据整理,见表1;②NDVI数据网格化,基于研究区露天煤矿各景观类型土地的面积变化与几何形状,筛选确定300 m×300 m的网格为相对较好的基本单元;③网格内植被稳定性指数的计算,包括能够反映网格内植被天然属性时空序列集中趋势与平均状态的NDVI均值(μ)、能够反映植被波动程度和差异性的NDVI均方差(σ)和反映植被群落时空演替方向和稳定性的NDVI相对涨落值(δ=σ/μ),具体计算公式参见文献[14];④通过3个植被稳定性指数的时空分布特征分析、局部异常涨落区域识别,检测植被生态系统稳定性变化特征。

图3 NDVI的时空数据序列分析方法的流程示意图Fig.3 Process of spatiotemporal data sequence analysis method based on NDVI

3 研究结果与分析

3.1 天然草地植被生态系统稳定性特征

在研究区选取10处受人类活动影响较小的天然草地作为背景区,对其植被稳定性指数的演替规律进行了研究,其结果见图4。

图4 研究区背景区天然草地NDVI均值(μ)、均方差(σ)和相对涨落值(δ)的时间序列动态图Fig.4 Time series of NDVI mean value(μ),mean square deviation(σ) and relative fluctuation(δ) of natural grassland in the background area of the study area

由图4可看出:对于原生背景区,中、高覆盖草地(μ≥0.5)区域(栅格2、3、4、5)植被NDVI的均方差(σ)和相对涨落值(δ)较低覆盖草地(0.5>μ≥0.2)区域(栅格1、6、7、8、9、10)植被NDVI的均方差(σ)和相对涨落值(δ)低,表明其植被生态系统稳定性更强。天然草地植被的NDVI均值(μ)为0.08～0.82,NDVI均方差(σ)为0.01～0.18,NDVI相对涨落值(δ)为0.02～0.45,NDVI相对涨落背景阈值为0.45。在时间尺度上,枯水年(2006年)天然草地植被分布更不均匀,相对涨落值更高,说明天然草地植被生态系统内部的振荡剧烈,系统的稳定性降低。

3.2 矿区各地段植被生态系统稳定性特征

研究区矿山NDVI均值、均方差、相对涨落时间序列动态图见图5,其中矿区渣山、矿区矿坑斜坡、矿区坑底平台3个植被稳定性指数的时间序列动态图,见图6至图8。

图5 研究区矿区NDVI均值(μ)、均方差(σ)和相对涨落值(δ)时间序列动态图Fig.5 Time series of NDVI mean value (μ),mean square deviation(σ) and relative fluctuation(δ) in the mine area of the study area

图6 研究区矿区渣山NDVI均值(μ)、均方差(σ)和相对涨落值(δ)时间序列动态图Fig.6 Time series of NDVI mean value(μ),mean square deviation(σ) and relative fluctuation(δ)of slag mountain area of the study area

3.2.1 矿区渣山

由图5和图6可以看出:

1) 2002年开采前,矿区渣山相对涨落值(δ)小于0.30,说明植被生态系统处于稳定状态。自2003年开采后矿区不断扩大,植被生态系统的稳定性受到干扰,矿山开采对植被生态系统的影响具有一定的滞后性。至2012年,低、中、高覆盖草地均变为裸地,区域植被波动明显,植被NDVI相对涨落值比背景区阈值高0.94,表明植被生态系统稳定性大幅降低,湿地与草地功能退化,植被生态系统被完全破坏。

2) 2015—2017年第一次治理后,矿区渣山Ⅰ.1南部原生低覆盖区(栅格93、94、109)的覆绿效果最为显著,植被分布波动小;而接近矿坑的渣山Ⅰ.1、渣山Ⅰ.2内(栅格56、70、71、95、110、111)覆绿效果并未达到预期;2019年区域植被NDVI相对涨落值较高,在0.28～1.10之间,说明初期采取的治理措施在高寒矿区适应性较差。

3) 2020—2022年第二次整治后,矿区渣山Ⅰ.1内植被NDVI均值恢复程度大于第一次治理,现观测覆绿区(栅格93、94、95)恢复到高NDVI均值低涨落状态,经土地整治但未覆绿处理的渣山南坡(栅格109、110)处于缓慢恢复状态,植被生态系统稳定性不佳,NDVI相对涨落值达0.80,而矿区渣山Ⅰ.2内各栅格区域相对涨落值分布均较低,接近背景区的天然状态。矿区渣山Ⅲ.3堆积时间较短,植被覆盖持续提高,2022年植被NDVI相对涨落值恢复至小于0.40的稳定状态。

4) 3处渣山由原生低、中、高覆盖草地均演变为裸地后演变为人工草地,标志植被生态系统稳定性的NDVI相对涨落值表现出先升高后降低的趋势。

3.2.2 矿区矿坑斜坡

由图5和图7可以看出:

图7 研究区矿区矿坑斜坡NDVI均值(μ)、均方差(σ)和相对涨落值(δ)时间序列动态图Fig.7 Time series of NDVI mean value(μ),mean square deviation(σ) and relative fluctuation(δ)of pit slope area of the study area

1) 在2003年斜坡未形成前,网格内植被NDVI相对涨落值小于0.62,处于自然演替状态;在2012年开采末期,阴坡斜坡Ⅱ.1(栅格51、67、68、83、84)由草地破坏演变为裸地,区域植被NDVI相对涨落值达1.39,比背景区阈值高0.94,说明植被生态系统被完全破坏,见图7。

2) 两处矿坑边坡在停采后维持现状,阳坡斜坡Ⅱ.2的植被NDVI相对涨落值与背景区阈值的差值为0.97,2015年治理后植被生态系统未能稳定,表明单纯地用覆土+人工建植方式对高寒矿区斜坡植被生态系统恢复效果较差。

3) 经后续2020年边坡整治、覆土改良、人工建植后植被覆盖度与稳定性发生了明显变化,目前NDVI相对涨落值回降到接近背景区天然草地的水平,区域植被NDVI相对涨落值分布范围为0.18～0.49,且阳坡斜坡Ⅱ.2的生态恢复效果相对阴坡斜坡Ⅱ.1更好。这表明覆土可以有效促进植物的生长,提高植被生态系统的稳定性,是高寒草原煤矿区比较理想的生态修复治理措施。

3.2.3 矿区坑底平台

由图5和图8可以看出:

图8 研究区矿区坑底平台NDVI均值(μ)、均方差(σ)和相对涨落值(δ)时间序列动态图Fig.8 Time series of NDVI mean value(μ),mean square deviation(σ) and relative fluctuation(δ) of pit bottom platform of the study area

1) 矿区5号井坑底平台在2002—2003年开采面积较小,处于自然涨落状态,网格内植被NDVI相对涨落值小于0.6。在2012年开采末期,矿区矿坑中西部(栅格53)早已停采,植被相对于2006年几乎无波动,植被NDVI相对涨落值为0.2,植被生长处于中覆盖草地稳定发展状态,矿坑在此期间向东部扩展,矿坑东部植被NDVI均值降至0.2,NDVI均方差小幅波动,但NDVI相对涨落值急剧增大,区域(栅格69、85、102、103、118)植被处于向裸地发展的不稳定过渡阶段。

2) 2015年采坑回填覆绿,植被覆盖度大幅增加,NDVI相对涨落值涨至1.41,比天然草地NDVI相对涨落值阈值高0.96,在停止维护后植被NDVI均值下降至比采矿末期更低的范围,表明此时的植被生态系统并未持续稳定;在2019年、2020年,草地(栅格68、85、102、119)完全失稳,区域处于趋向裸地发展的近裸地状态;直至2021年、2022年覆绿,植被NDVI均值回升,区域NDVI均方差达到0.30后下降,NDVI相对涨落值的范围为0.31～0.80,小幅波动后下降,向良好态势发展,植被生态系统处于向自然草地演进的近自然草地阶段,植被NDVI均方差与相对涨落值接近天然草地。

以上结果表明:高寒草原开采后矿区坑底平台区域植被生态系统几乎彻底失稳,在2015年退化到近裸地状态;高寒草原矿区2015—2017年边采边治的策略持续加大了对区域草地及冻土的破坏;在彻底停止开采后2017—2019年,坑底平台依然不能有效恢复植被;经2020—2022年的“覆土+施肥+人工建植”生态恢复治理,区域植被NDVI相对涨落值现已降至近天然草地背景区的NDVI相对涨落值范围(0.20～0.59),表明该措施是恢复高寒草原矿区植被生态系统的有效途径。

4 讨 论

对比研究区矿山两次生态修复治理效果发现,土壤改良在高寒草地矿区不同地形条件下均能显著提高植被生态系统稳定性,促进植被生态系统正向演替发展,地形是影响矿区植被生态环境恢复的主要因子,高寒矿山受损植被生态系统的治理可通过优化矿区渣山与阴坡斜坡的生态修复治理措施来进一步提升其生态环境恢复效率。

由于生态环境修复时间较短,矿区植被生态系统现阶段表现出趋于稳定化的态势,还需要后期多年的监测才可以完整定论。

5 结 论

1) 在宏观尺度上,地形与气候条件决定着植被空间分布格局,降水与植被NDVI密切相关。采矿与生态修复等人类活动尚未影响到区域植被生态系统的稳定性,背景区天然草地的NDVI均值(μ)范围为0.08～0.82,NDVI均方差(σ)小于0.18,NDVI相对涨落(δ)阈值为0.45。

2) 与背景区相比,矿山开采后矿区植被NDVI相对涨落值明显增加,植被生态系统失稳,矿区植被NDVI相对涨落值与背景区阈值的差值可达0.97。

3) 2015—2017年第一次生态修复治理后,矿区植被NDVI相对涨落值未能有效降低,2019年该值范围为0.28～1.10,覆土+人工建植方式的生态修复治理效果不明显。

4) 2020—2022年第二次生态修复治理后,经矿区渣土改良、土地增稳、养分提升,植被生态系统覆盖度与稳定性明显改善,治理区植被NDVI相对涨落值降低至天然草地背景区相似水平(0.07～0.59),植被生态系统稳定性排序为矿区坑底平台>矿区矿坑斜坡>矿区渣山。