基于多期遥感数据的华蓥山中段植被覆盖变化研究

杨 玲，周学铖，熊丽媛，徐争强，胡 祎，吕 品

(1.四川省核工业地质调查院，四川 成都 610052；2.四川川核地质工程有限公司，四川 成都 610052)

0 引言

华蓥山突起于四川盆地底部，以西为川东丘陵，以东为川东平行岭谷，境内森林资源丰富，是川渝重要的生态屏障。近年来，矿业开采、旅游开发及城镇建设等人类工程活动的持续进行对华蓥山生态环境造成了不同程度的破坏[1]。1989年7月10日华蓥市溪口镇发生特大滑坡—溪口滑坡，伤亡221余人，直接经济损失逾600万元[2]。该滑坡发生后，华蓥山地区地质环境开始受到学者的广泛关注，但涉及华蓥山生态环境的相关研究较少。因为，为了填补华蓥山中段生态环境研究空白，减免因人类工程活动对华蓥山生态环境造成破坏，需要对研究区历史生态环境条件进行调查并研究其发展趋势及变化原因。

基于遥感技术进行区域生态环境变化分析及研究已较为广泛。方洪宾等利用20世纪70年代中期MSS、1990年左右TM及2000年左右ETM数据为主，通过对各专项因子调查，分析了青藏高原生态地质环境形势及相关变化趋势[3]；李洪义等利用ETM数据提取了土壤亮度、植被、湿度、热度指数等参数并结合气象和地学资料建立模型，评价了福建省生态环境质量[4]；王俊华等利用GF2，OLI数据，采用人机交互解译的方式收集了西藏多龙矿集区人类工程活动、水文、荒漠化、金属物氧化污染及植被等信息，评价了当地生态环境现状并分析了影响生态环境问题的原因[5]。

前人主要依据单时相遥感数据或宽时相多源遥感数据进行区域生态环境调查及分析研究，而基于窄时相多源遥感数据对生态环境研究相对较少。因此，本文以广安市境内华蓥山中段为研究对象，在充分收集研究区多期遥感影像数据的基础上，通过遥感定量计算和野外调查，分析了各时期研究区植被覆盖情况、变化趋势及变化原因，研究成果不仅可以为华蓥山生态环境规划与防治工作提供科学依据，而且可以为地表覆被率变迁与生态环境承载力间的关系研究提供新思路。

1 研究区概况

研究区主要位于四川省广安市境内华蓥山区，是华蓥山脉中段核心区，地理坐标范围为106°39'56″E-107°03'17″E，北纬30°03'36″N-30°38'52″N之间，面积约776.89km2，包括广安市前锋区东部、华蓥市东部及邻水县西部。

图1 研究区位置图

华蓥山地区属中亚热带湿润季风气候区，气候平和，雨水丰沛，年均降雨量为1259.13mm。区内地表径流发育、沟谷纵横、水网密度较大，主体以华蓥山脉岭脊线为分水岭，以西属长江二级支流渠江流域，以东属长江一级支流御临河流域。天池湖位于研究区北部，为天然间歇性岩溶湖泊，面积2.53km2，是川东海拔最高、面积最大的湖泊。华蓥山平均海拔1200～1400m，山势陡峭，多形成“V”字形峡谷，山谷深而窄，溪沟发育，最高峰为高登山，海拔1704.1m。地貌类型有脊状岩溶中山低山、岩溶槽谷、单面山列峰脊状中低山。

研究区内出露的地层主要有寒武系、奥陶系、志留系、石炭系、二叠系、三叠系、侏罗系及第四系地层。寒武系-三叠系碳酸盐岩出露面积最大，三叠系-侏罗系碎屑岩主要分布于两翼，第四系松散堆积层零星不连续出露。研究区内矿产资源丰富，采矿历史悠久，分布各类矿山共186个，主要为地下开采的煤矿山与露天开采的石灰岩矿山，长期、无序、大规模的采矿活动对华蓥山生态环境造成了巨大威胁。

2 数据源及技术路线

2.1 数据源及预处理

为进行华蓥山中段植被覆盖变化状况研究，收集了区内基础地理、地质、矿山、地灾等资料及Landsat-5、8遥感数据。

本文选择1988年至2018年间共7期美国陆地卫星LANDSAT系列遥感数据进行研究区植被覆盖变化变化研究。分别为1988年9月15日、1993年5月24日、2001年7月17日、2006年4月10日、2010年9月21日共五期Landsat-5 TM遥感影像和2014年8月6日、2018年4月4日共两期Landsat-8 OLI遥感影像，影像空间分辨率均为30m。各期影像均无云、雾等影响，影像拍摄日期分布在4月至10月之间，为研究区植被生长较好季节。在具体分析之前，对各期影像分别进行了辐射定标、大气校正及几何校正。影像数据来源于地理空间数据云。

2.2 归一化植被指数

NDVI计算公式为：

NDVI=((NIR-R))/(NIR+R)

(1)

其中：NIR、R分别为近红外波段及红波段。

NDVI可以部分消除或减弱因卫星观测角、太阳高度角、地形起伏、少量云影及大气辐照对成像的影响[6-9]。广泛应用于基于遥感影像的植被研究，是植物生长状态、覆盖度最佳指示因子[6]。

本文使用的TM及OLI数据，近红外波段分别为第4波段、第5波段，红波段分别为第3波段、第4波段，基于ENVI软件Band Math工具可获取研究区不同时期NDVI分布图。

2.3 植被覆盖度

基于像元二分模型进行植被覆盖度计算，即：

VFC=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(2)

其中VFC为像元植被覆盖度，NDVIveg为无植被覆盖区NDVI值， NDVIsoil为100%植被覆盖区NDVI值。

NDVIsoil=(VFCmax·NDVImin-VFCmin·NDVImax)/((VFCmax-VFCmin))

(3)

NDVIveg=((1-VFCmin)·NDVImax-(1-VFCmin)·NDVImin)/((VFCmax-VFCmin))

(4)

式(3)及式(4)中VFCmax、VFCmin、NDVImax、NDVImin分别代表区内植被覆盖最大值、植被覆盖最小值、NDVI最大值、NDVI最小值[6]。

经目视粗解译，研究区多期影像中均存在植被高富集区，且多期影像中均存在无植被覆盖区(天池湖区)，因此认为研究区可近似取VFCmax=100%，VFCmin=0%。

因此公式(2)可改写为：

VFC=(NDVI-NDVImin)/(NDVImax-NDVImin)

(5)

排除不可避免的影像噪声，分别取NDVI统计值累计百分比的5%和90%所对应的值作为置信阈值，认为NDVI值位于0～5%区间的为无植被覆盖区，位于90%～100%区间的为纯植被覆盖区，取对应的值作为NDVImax和NDVImin。

基于ENVI软件对7期NDVI图进行统计，获取相应的NDVImax和NDVImin，按公式(5)通过Band Math工具计算得出各时期研究区植被覆盖度图。

2.4 植被覆盖度分级及植被覆盖贡献值

本文依据《草场资源调查技术规程》、《土地利用现状调查技术规程》、《森林资源规划设计调查主要技术规定》等相关条款及规定[8,10]，结合研究区植被覆盖实际情况，对植被覆盖度分为6级。

植被覆盖贡献值即各级覆盖区域对植被覆盖总面积的贡献值，其中Ⅰ级区为无植被覆盖区，其对植被覆盖贡献值为0；Ⅱ级区至Ⅴ级区以植被覆盖度均值记入；Ⅵ级区为近纯植被区，本文认为Ⅵ级区植被覆盖贡献值为1(表1)。

表1 植被覆盖度分级表

3 结果分析

基于ENVI计算及处理，获取了研究区1988年至2018年共7期植被覆盖度图(图2)，并基于表2统计了各年度植被覆盖情况(表2)，绘制了相应时序曲线(图3)。

3.1 植被覆盖变化情况

依据表2及图3，1988年9月至1993年5月，研究区植被覆盖总体面积减少了125.93km2，总体植被覆盖率降低16.21%，其中Ⅵ级密林面积减少143.96km2，减幅最大，Ⅱ级区面积增加了78.79km2，增幅最大；

1993年5月至2001年7月，研究区植被覆盖总体面积增加了82.72km2，总体植被覆盖率增加10.64%，其中Ⅵ级密林面积增加了91.99km2，增幅最大，Ⅱ级区面积减少了60.07km2，减幅最大；

2001年7月至2006年4月，研究区植被覆盖总体面积增加了10.74km2，总体植被覆盖率增加1.39%，其中Ⅵ级密林面积增加11.73km2，增幅最大，Ⅰ级区面积减少了6.2 km2，减幅最大；

2006年4月至2010年9月，研究区植被覆盖总体面积增加了28.96km2，总体植被覆盖率增加3.72%，其中Ⅵ级密林面积增加36.44km2，增幅最大，Ⅱ级区面积减少了18.35km2，减幅最大；

2010年9月至2014年8月，研究区植被覆盖总体面积增加了27.61km2，总体植被覆盖率增加3.56%，其中Ⅵ级密林面积增加59.6km2，增幅最大，Ⅴ级区面积减少了33.59 km2，减幅最大；

2014年8月至2018年4月，研究区植被覆盖总体面积减少了18.19km2，总体植被覆盖率增加2.34%，其中Ⅵ级密林面积增加43.17km2，增幅最大，Ⅴ级区面积减少了31.94km2，减幅最大；

图2 研究区各期植被覆盖度图注：A：1988年9月15日植被覆盖图；B：1993年5月24日植被覆盖图；C：2001年7月17日植被覆盖图；D：2006年4月10日植被覆盖图；E：2010年9月21日植被覆盖图；F：2014年8月6日植被覆盖图；G：2018年4月4日植被覆盖图

表2 植被覆盖面积统计表

图3 1988年9月至2018年8月研究区各级植被覆盖面积时序曲线图注：图B为图A中Ⅰ、Ⅱ、Ⅲ、Ⅳ曲线图放大显示。Note：Figure B shows the enlarged graphs of I，II，III and IV in Figure A.

3.2 植被覆盖变化原因分析

自1988年9月至2018年4月，研究区植被覆盖率历经两个阶段。

(1)林地大面积破坏—反弹性恢复阶段(1988年至2001年)

本文分析显示，1988年至1993年研究区植被覆盖面积大幅较少，至2001年研究区植被面积较1993年有明显恢复(图4)。据此实地走访并取证，20世纪80年代至90年代华蓥山小煤窑无序开采活动频繁，小煤窑数量一度达500多个(数据来源于广安市国土局)，至2017年研究区尚存93处煤矿，但均已闭坑。煤矿硐采需大量坑木及矿柱，加之大修矿山工业广场、大量煤矸石无序堆积，导致林地林木被毁，森林生态环境遭到极大的破坏。

为扭转森林覆盖面积下降趋势，改善森林生态环境，1998年10月四川省广安市开启了《天然林资源保护》工程，并于次年10月再次开启了《退耕还林》工程[11]，从遥感影像分析成果来看，两大工程对研究区植被覆盖恢复起到了扭转性作用。

(2)全区植被覆盖率缓慢增长阶段(2001年至2018年)

自2001年至2018年，研究区总体植被覆盖率增加量分别为2001年7月至2006年4月增加1.39%，2006年4月至2010年9月增加3.72%，2010年9月至2014年8月增加3.56%，2014年8月至2018年4月增加2.34%。

基于该组数据做出如下判断：研究区植被覆盖总量处于缓慢上升阶段，其增长速率呈缓慢(2001年至2006年)—加快(2006年至2010年)—增速轻微放缓(2010年至2014年)—增速放缓(2014年至2018年)态势。

笔者认为，维持该阶段研究区植被覆盖面积缓慢但相对稳定增长有两方面因素：

图4 1988年9月至2001年7月研究区局部植被覆盖图注：A：1988年植被覆盖图；B：1993年植被覆盖图；C：2001年植被覆盖图；D：1988年影像图；E：1993年影像影像图；F：2001年影像影像图。

(1)促进因素。“两大工程”持续发挥的直接促进作用；两大工程实施过程中造成研究区农民收入减少[11]，随后农村人员外流，使得区内农村社会经济规模趋于变小，对植被覆盖面积增加起到了间接促进作用[12]；环境的自我恢复力，研究区存在大量废弃露天开采矿山，废弃时间5-10年不等，其开采平台大部分已自然复绿，开采面也存在草本植物生长现象。

(2)抑制因素。露天矿山开采造成的影响：通过对研究区7期遥感影像进行目视解译，获取了各期影像中露天开采矿山图斑面积，分别为1988年16.98ha，1993年24.9ha，2001年79.16ha，2006年118.32ha，2010年287.81ha，2014年379.39ha，2018年412.37ha。由数据可知，自2001年开始研究区露天矿山开采活动频率加快，对植被覆盖面积造成不利影响；

煤矿地下矿山开采造成的影响：据2017年广安市矿山地质环境详查，研究区煤矿开采造成区内含水层破坏与疏干面积约364.33km2，占研究区总面积的46.9%，大面积的含水层破坏对植被生长造成不利影响[13]。

其他直接性人类工程活动造成的无植被覆盖区面积增加，如近年来研究区开发的华蓥山旅游文化景区、华蓥山森林公园、华蓥山地质公园等工程项目。

总体来看，在不受不可抗力因素的干扰下，人类活动是促使研究区植被覆盖变化的主导因素，即人类活动即可造成林区面积大幅减少或退化，亦可快速恢复受破坏区植被，同时研究区环境本身的自愈力也在时刻发挥着促进植被覆盖增长的重要作用。

3.3 植被覆盖变化趋势分析

依据现有数据推测，假定研究区人与自然和谐共处且其他自然条件均处理想状态，则随着时间无限推移，研究区植被覆盖结构将减少至两类：Ⅵ级区和Ⅰ级区，即稳定的无植被覆盖区和稳定的纯植被覆盖区。

沿着上述时间轴反推，即假定在理想状态下，时间推移无限年后研究区植被覆盖结构仅为Ⅵ级区和Ⅰ级区，则沿该时刻反推，在无限年之前研究区植被覆盖结构亦为Ⅵ级区和Ⅰ级区。

在不彻底破坏生态结构的前提下，外界因素可导致Ⅵ级区和Ⅰ级区可相互转化。

4 结论

通过对华蓥山中段植被覆盖度进行调查研究，获取了自1988年9月至2018年4月共7期植被覆盖情况图，基于影像计算成果进行分析，得出了以下结论：

(1)基于TM、OLI等遥感数据对区域植被覆盖度定量探查与分析是可行的，其计算成果较好反映了当时植被覆盖情况，并通过多期计算成果对比可以反映区域植被覆盖变化总体趋势。

(2)人类活动是导致生态环境变化的重要因素。大量获取林木资源、无序矿山开采，易造成植被覆盖面积持续减少，对生态环境不利；较少人为干扰环境下，林区植被呈自然缓慢恢复态势；在有利政策或工程推动下，植被覆盖面积迅速增长，生态环境持续向好。

(3)植被覆盖等级是可以相互转化的，对林区资源的合理利用，并及时加以修复，即可达到社会经济发展的需要，也可维持生态环境持续稳定，最终可以达成人与自然双赢的形势。