基于遥感生态指数的龙溪—虹口国家级自然保护区生态环境状况评估

刘立冰,熊康宁,任晓冬①

(1.贵州师范大学喀斯特研究院,贵州 贵阳 550001;2.国家喀斯特石漠化防治工程技术研究中心,贵州 贵阳 550001)

随着社会经济的快速发展以及自然环境的不断改变,生物多样性减少、水土流失、土地荒漠化等环境问题层出不穷,给社会经济可持续发展带来严重挑战。如何快速准确地评估区域生态环境状况,成为当今迫切需要解决的问题。传统的生态环境质量评价方法包括模糊评价法[1]、层次分析法[2]、PSR模型[3]等,这些方法受主观因素影响较大,主成分分析法[4]、人工神经网络法[5]的出现解决了主观权重赋值的问题。2006年原环境保护部发布的《生态环境状况评价技术规范(试行)》为生态环境状况评价提供了标准,生态环境状况指数(EI)得到广泛应用,但其指标权重需要人为设定,且存在评价结果无法空间可视化的问题。伴随着遥感与地理信息技术的迅猛发展,徐涵秋[6]在对福建长汀水土流失区进行生态环境质量评价时,率先提出遥感生态指数(RSEI),由于数据完全通过遥感获取,大大降低了数据收集的难度,还可以对数据进行可视化操作,便于不同时期之间的对比，利用主成分分析法进行权重赋值,提高了评价结果的客观性。近年来,众多专家学者利用该方法对城市[7]、矿区[8-9]、沙漠[10]、黄土流失区[11]等不同区域进行生态环境评价,其中也不乏针对自然保护区的研究[12-13]。

龙溪—虹口国家级自然保护区所在的横断山北段动植物资源极其丰富,是全球34个生物多样性热点区域之一[14]。保护区为岷山山系大熊猫B种群的重要栖息地,是连接2个最大的大熊猫种群(岷山种群和邛崃种群)的通道,是世界自然遗产——四川大熊猫栖息地的重要组成部分。2018年1月,大熊猫国家公园成都管理分局依托龙溪—虹口自然保护区管理局设立大熊猫国家公园试点。保护区位于四川盆地西北部,对于成都平原的水源涵养、水土保持意义重大。保护区距2008年汶川地震的震中映秀镇仅4 km,受地震及滑坡、泥石流等次生灾害的影响,区内生态环境严重退化,大熊猫栖息地严重受损,植被受损面积达100 km2以上,大熊猫生境大量丧失[15]。

目前,对于龙溪—虹口自然保护区的研究主要集中在景观多样性[16]、大熊猫与金丝猴栖息地[17-18]、植被自然演替[19]等领域,鲜见针对保护区进行生态环境状况评估的研究。笔者基于遥感生态指数(RSEI),集成绿度、湿度、热度和干度4项指标,利用主成分分析法建立评估模型,评估龙溪—虹口自然保护区的生态环境状况,综合分析1998—2017年保护区生态环境状况变化,着重分析2008年汶川地震前后保护区的生态环境变化程度,揭示地震10 a后其生态环境恢复状况,并提出针对性的保护对策,以期为保护区生态环境治理提供科学支撑。

1 研究区概况与研究方法

1.1 研究区概况

龙溪—虹口自然保护区位于四川省都江堰市西北部(图1),地处四川盆地向青藏高原的过渡地带,成立于1988年,是以保护大熊猫及其栖息地为主的保护区,1997年升级为国家级自然保护区,地理坐标为31°04′～31°22′ N,103°32′～103°43′ E,保护区总面积310 km2,占都江堰市总面积的1/4,包括核心区203 km2、缓冲区37 km2、实验区70 km2,保护区外围的保护带面积为117 km2。

图1 保护区高程及水系图

保护区属华夏地质构造体系,处于龙门山褶皱带中南部,主要河流为白沙河和龙溪河。保护区最高峰光光山海拔4 582 m,最低海拔1 200 m,相对高差达3 382 m。保护区属亚热带湿润气候区,年降水量约1 600～1 900 mm,且80%集中于3—9月,年平均气温10 ℃左右,1月最低气温可达-10 ℃以下,7月最高气温25 ℃左右[20]。

保护区植被类型多样,海拔2 400 m以下为亚热带常绿阔叶林、常绿落叶阔叶混交林、落叶阔叶针叶混交林,海拔>2 400～3 800 m为亚高山暗针叶林和亚高山灌丛,海拔>3 800～4 582 m为高山草甸和高山流石滩稀疏植被。

1.2 数据源及预处理

根据研究需要,从美国地质勘探局网站(https:∥earthexplorer.usgs.gov/)获取1998年7月的Landsat 5 TM影像、地震前2007年9月的Landsat 5 TM影像、地震后2008年8月的Landsat 5 TM影像以及2018年7月的Landsat 8 OLI影像,云量均低于20%,影像质量较好,植被处于生长季,4个时期遥感影像的时相基本一致,可以确保不同时期研究结果的准确性和可比性。统计数据来源于保护区管理局或通过野外实地调查获取。

运用遥感分析软件ENVI 5.3对4个不同时期的影像进行辐射校正(辐射定标和大气校正)和几何校正。首先,利用软件自带的辐射定标工具(Radiometric Calibration)对遥感影像进行辐射定标,其参数可从遥感影像的元文件中自动获取,将像元亮度值(DN)转换为辐射亮度值;其次,利用大气校正工具(FLAASH Atmospheric Correction)减少或消除大气和光照等因素对地物反射的影响;然后,对4个不同时期的遥感影像进行几何校正,利用二次多项式和最邻近像元法配准,为满足研究精度的需要,将均方根误差控制在0.5个像元内;最后,使用龙溪—虹口自然保护区的矢量边界对图形进行裁剪。

1.3 研究方法

选取易获取的归一化植被指数、湿度分量、干度指数和地表温度代表绿度(NDVI,INDV)、湿度(WET,TWE)、干度(NDSI,INDS)和热度(LST,TLS),对各项指标进行标准化处理,利用ArcGIS软件的主成分分析工具将遥感生态指数(RSEI)表示为4项指标合成的函数。

1.3.1绿度指标

NDVI(INDV)作为应用最广泛的植被指数之一,是植被生长及空间分布状态的最佳指示因子[21],计算公式为

INDV=(ρNIR-ρred)/(ρNIR+ρred)。

(1)

式(1)中,ρNIR和ρred分别为近红外波段和红波段的反射率。

1.3.2湿度指标

土壤湿度为陆地和大气能量交换过程中的关键因子,是气候、水文、生态、农业等领域的重要指标。经过缨帽变换后的影像可以较好地区分土壤、植被及作物信息,被广泛应用于生态环境状况评估。缨帽变换的第3分量对土壤湿度很敏感,可以很好地反映植被、水体和土壤湿度状况,因此被称为湿度分量(WET,TWE)。不同传感器的相关参数设置有差异,TM和OLI影像的温度分量计算公式[22-23]分别为

TWE，TM=0.031 5ρblue+0.202 1ρgreen+0.310 2ρred+0.159 4ρNIR-0.680 6ρSWIR1-0.610 9ρSWIR2,

(2)

TWE,OLI=0.151 1ρblue+0.197 2ρgreen+0.328 3ρred+0.340 7ρNIR-0.711 7ρSWIR1-0.455 9ρSWIR2。

(3)

式(2)～(3)中,ρblue、ρgreen、ρred、ρNIR、ρSWIR1、ρSWIR2分别为各影像的蓝波段、绿波段、红波段、近红外波段、短波红外1 波段和短波红外2波段的反射率。

1.3.3干度指标

干度指标可表征土地荒漠化和退化状况,由于研究区域内存在一部分建筑用地和裸地,所以选取裸土指数(SI,IS)和建筑指数(IBI,IIB)来合成干度指标(NDSI,INDS),计算公式[24-25]为

INDS=(IS+IIB)/2,

(4)

IS=(ρSWIR1+ρred-ρblue-ρNIR)/(ρSWIR1+ρred+ρblue+ρNIR),

(5)

IIB=[2ρSWIR1/(ρSWIR1+ρNIR)-ρNIR/(ρNIR+ρred)-ρgreen/(ρgreen+ρSWIR1)]/[2ρSWIR1/(ρSWIR1+ρNIR)+ρNIR/(ρNIR+ρred)+ρgreen/(ρgreen+ρSWIR1)]。

(6)

1.3.4热度指标

热度指标通常通过陆地表面温度计算[26],从遥感卫星数据反演出的地表温度可以很好地反映单个像元的下垫面平均温度,是表征地球表面能量平衡和温室效应的有效指标。运用Landsat TM影像反演地表温度的方法很多,常用的方法主要有大气校正法、单通道算法、分裂窗算法(劈窗算法)和多波段算法等。采用3期影像的热红外波段(Landsat 5影像选择第6波段,Landsat 8选择第10波段)进行辐射定标,得到热红外波段的辐射亮度值L6/10,根据普朗克公式的反函数求得地表真实温度,计算公式为

L6/10=gain×ND+bias，

(7)

BLST=[L6/10-L1-τ(1-ε6/10)L2]/τ×ε6/10，

(8)

TLS=K2/ln (K1/BLST+1)-273。

(9)

式(7)～(9)中,L6/10为热红外波段在传感器处的辐射值,W·m-2·sr-1·μm-1;gain和bias分别为热红外波段的增益与偏置值,W·m-2·sr-1·μm-1;ND为灰度值,其值为常数;K1和K2为定标参数,单位分别为W·m-2·sr-1·μm-1和K;L1和L2分别为大气向上、向下辐射亮度,W·m-2·sr-1·μm-1;τ为热红外波段的透过率,其值为常数;BLST为黑体辐射亮度值,W·m-2·sr-1·μm-1;ε6/10为地表比辐射率，其值为常数。L1、L2、τ值可以在美国宇航局(NASA)提供的大气剖面参数网站(http:∥atmcorr.gsfc.nasa.gov/)或元数据中直接获取,由于该网站无法提供1998年数据,该研究用2007、2008和2017年的平均值代替[27]。

目前，惠州港还没有建设专用防台锚地，生产锚地严重不足，要解决这些问题，对锚地进行尽早重新合理规划，强化防台锚地的规划建设。加强与地方各相关部门的沟通联系，将锚地规划纳入“多规合一”统筹规划，加快实施海洋功能区划，合理布局港口功能，强化锚地用海管理，预留足够防台锚地，协调解决船舶锚泊安全与海洋养殖业、临港经济开发之间的矛盾。是解决锚地资源紧张的合理途径。

地表比辐射率表征地表向外辐射电磁波的能力,是热红外遥感的关键参数之一,其值主要取决于地表的物质结构和遥感器的波段区间,与NDVI具有较高的线性相关性。水体的地表比辐射率为0.995。根据SOBRINO等[28]的研究,利用NDVI阈值法计算地表比辐射率ε,当NDVI值<0.2时,为纯裸土像元,其地表比辐射率εs=0.98;当NDVI值>0.5时,为纯植被像元,其地表比辐射率εv=0.99;当0.2≤NDVI值≤0.5时,为混合像元。其计算公式为

ε=εvPv+εs(1-Pv)+dε,

(10)

Pv=(INDV-INDV,s)/(INDV,v-INDV,s),

(11)

dε=(1-εs)(1-Pv)Fεv。

(12)

式(10)～(12)中,Pv为植被覆盖度;dε为地形几何形状,相对平整的情况下可以忽略,混合像元和粗燥表面则必须考虑其影响;F为地形几何形状参数,取经验值0.55;INDV，s为纯裸土像元的NDVI值,取经验值0.2;INDV，v为纯植被像元的NDVI值,取经验值0.5。

1.3.5遥感生态指数(RSEI)评估模型的构建

(1)指标标准化处理

由于各项指标具有单位和季相差异,需要对指标进行标准化处理,计算公式为

IN=(Imax-I)/(Imax-Imin)。

(13)

式(13)中,IN为标准化后的指标值;I为该指标实际值;Imax和Imin分别为该指标的最大值和最小值。

(2)遥感生态指数计算

主成分分析是数据挖掘中主要的降维方法之一[29]。其权重根据各指标对第1主成分(PC1)的贡献度来确定,可避免人为因素的影响,保证权重赋值的客观性。

将标准化处理后的各项指标组合成4个波段的文件,借助ENVI软件对其进行主成分变换,用1减去PC1来构建初始遥感生态指数,计算公式为

IRSE=(I0-I0，min)/(I0，max-I0，min)。

(14)

式(14)中,I0，max、I0，min分别为初始遥感生态指数的最大值和最小值；IRSE为经过标准化处理的遥感生态指数,取值范围为[0,1]。RSEI值越大表示生态环境状况越好,反之则代表生态环境状况越差。将其值划分为5个等级[30],0～0.2、>0.2～0.4、>0.4～0.6、>0.6～0.8和>0.8～1分别对应差(1级)、较差(2级)、中等(3级)、良(4 级)、优(5 级)。

2 结果与分析

2.1 龙溪—虹口自然保护区生态环境状况总体评估

1998—2017年4个指标的主成分分析结果见表1。

表1 1998—2017年4个指标的主成分分析结果

由表1可知,1998—2017年绿度和湿度均为正值,两者与生态环境状况呈正相关关系,共同对保护区的生态环境状况起正面的促进作用。干度和热度均为负值,它们与生态环境状况呈负相关关系,两者共同起负面作用。1998、2007、2008和2017年绿度、湿度、干度和热度在第1主成分的贡献率分别为89.61%、95.25%、89.28%和90.24%,均大于85%,表明第1主成分已经集成4个指标的大部分特征。由各指标的第1主成分载荷可知,与该区域生态环境状况呈正相关关系的湿度、绿度中,绿度的贡献率更高,1998、2007、2008和2017年绿度贡献率分别为0.750 8、0.885 1、0.652 5和0.623 4,表明植被生长及分布情况是研究区生态环境状况的直接影响因素。

由表2可知,1998—2017年保护区遥感生态指数从0.748 1变为0.599 1,下降19.92%。2007年生态环境状况最佳,遥感生态指数为0.752 1;2008年最差,遥感生态指数仅为0.411 0。对生态环境状况有正向影响的绿度均值呈升高—降低—升高趋势,湿度均值呈降低—降低—升高趋势;对生态环境状况有负向影响的干度和热度指标均表现为降低—升高—降低趋势。近20 a来,龙溪—虹口自然保护区的生态环境状况总体并未保持稳定,表现出升高—降低—升高的趋势,整体生态环境呈退化趋势。

表2 1998—2017年保护区遥感生态指数(RSEI)

2.2 生态环境状况分级

从表3可知,1998年保护区生态环境状况以优和良为主,合计占比达89.78%;等级为中的区域面积占比为7.64%;较差和差等级区域面积很小,只占2.59%。2007年保护区生态环境状况仍以优和良为主,占比高达91.52%;中、较差和差等级区域面积很小,合计占比为8.48%。2008年总体生态环境状况以较差为主,面积占比为52.91%;其次为中,面积占比为44.89%;优、良和差等级区域面积都较小,合计占比为2.21%。2017年保护区生态环境状况以良和中2个等级为主,合计占比达92.03%;较差和差等级区域面积下降幅度较大,占比仅为6.39%。

表3 1998—2017年保护区不同等级遥感生态指数(RSEI)区域面积及占比

2.3 地震前后保护区生态环境状况变化

从空间分布(图2)来看,1998年生态环境状况为优和良的区域主要位于保护区西部、东部和中偏北部等阔叶林、落叶阔叶混交林、亚高山暗针叶林分布区域,保护区北部为高山流石滩,植被稀疏,覆盖度低,海拔高且气温低,生态环境状况等级为中、较差和差。地震前的2007年,生态环境状况为优和良的区域遍布整个保护区。这是由于1997年保护区升级为国家级自然保护区,2003年完成勘界定桩,通过开展专项武装巡护、在各入山口设立保护区宣传标牌等措施,使绝大部分区域的保护效果较好。

2008年汶川地震发生后,因滑坡、崩塌而出现的裸露土地零散分布于整个保护区内,森林以及灌木林面积大幅减少[31],保护区最北部由于本身环境质量较差,加之地震影响,退化程度加重,海拔2 500 m以下的大部分区域生态环境状况为中和较差,生态环境状况为优的区域主要分布在河道附近。

到2017年,通过近10 a的恢复,原本生态环境状况差和较差的区域也逐渐被中等、良取代,差和较差的区域零散分布于保护区北部、中部及河道附近区域,生态环境状况较2008年有好转趋势,但仍未恢复到地震前2007年的水平。

图2 1998—2017年保护区不同等级遥感生态指数的空间分布

利用ENVI软件自带的变化检测工具,对4个不同时期的遥感生态指数进行变化分析。由表4可知,1998—2007年保护区生态环境动态变化较为明显,以上升1个等级〔变好(+1)〕和下降1个等级〔变差(-1)〕为主。变好区域面积为168.05 km2,占比超过45.777%;生态环境状况不变的区域面积不足1 km2,占比也未超过1%;生态环境状况变差的区域面积为141.9 km2。对地震前后的2007—2008年进行对比,地震后生态环境状况明显变差,面积达296.34 km2,占比达95.594%,以下降2个等级〔变差(-2)〕为主;变好和不变的区域占比很少。对比2007—2017年,生态环境状况总体体表现为变差,变差区域面积为262.99 km2,占比超过84.836%,而变好区域面积仅为45.27 km2。

表4 不同时段保护区遥感生态指数等级变化

2.4 生态环境状况变化成因分析

从图3可知,遥感生态指数等级变化的空间分布在不同时期有不同特征。1998—2007年,生态环境状况变好的区域主要分布于保护区内的阔叶林、针叶林等植被茂盛区域,这主要是由于保护区对区内资源的保护得力,特别是在“天保”和“退耕还林”工程实施后,生态环境明显好转。2007—2008年,由于地震造成生态环境恶化,保护区内景观离散度增大,破碎化严重,中山和亚高山(海拔1 400～2 600 m)区域生态系统森林及灌木林受到破坏,土壤保持能力与水源涵养功能减弱,水土流失严重。地震还造成大量熊猫栖息地损失,地震前保护区内有大熊猫栖息地272.99 km2,地震造成栖息地损毁97.06 km2,损毁比例高达35.56%。四川省第四次大熊猫调查(2011—2014年)结果显示,其潜在栖息地几乎全部丧失,并且形成多个堰塞湖,由于水域面积增加,该区域的生态环境质量反而好转。2007—2017年,大部分区域生态环境状况表现为变差,生态环境状况变好的区域集中分布在保护区最北部,地震对北部的高山流石滩稀疏植被影响较小,植被得到一定程度的恢复;变差区域主要分布于中山和亚高山区域,以白沙河及其附近区域最为明显,大量的泥沙和石砾冲刷到河道中,改变了生物廊道内部的结构和功能,进而影响研究区整体的生境、传输通道以及过滤和阻碍等廊道功能。

保护区建设初期,社区居民对生物资源造成了较大程度的破坏。1999—2001年,保护区内的挖药、采笋等行为仍大量存在,后随着旅游业的兴起,进入保护区采药和采笋的人数显著降低。地震后,保护区地貌复杂,坡度大,人类的可进入性低,所以人类干扰并不严重。总体来看,虽然保护区生态环境状况还未完全恢复至地震前的水平,但大熊猫数量也由四川省第三次大熊猫调查(1998—2002年)时的6只增长为第四次调查(2011—2014年)时的9只[32]。

括号内数字为遥感生态指数上升或下降的等级数。

3 结论

自然保护区具有显著的环境正效应,利用4期Landsat TM/OLI遥感影像,运用主成分分析法构建由绿度、湿度、干度、热度这4项指标组成的遥感生态指数评估模型,能够快速定量地对1998—2017年龙溪—虹口自然保护区的生态环境状况进行总体评价,并对其空间分布和变化进行监测,结果表明:

(1)龙溪—虹口自然保护区1998、2007、2008和2017年的遥感生态指数值分别为0.748 1、0.752 1、0.411 0和0.599 1,呈现先升高后降低再升高的的变化过程,整体来看生态环境状况表现为退化趋势。

(2)地震前后保护区生态环境状况变化明显,海拔1 400～2 600 m的中山和亚高山地区生态环境状况呈明显恶化趋势,主要原因是其生态环境较为脆弱,加之地震后大量植被遭受破坏,水域面积减少,同时泥石流、滑坡等地质灾害频繁致使地表破碎化程度增高。

(3)绿度和湿度对生态环境状况起正向作用且影响较大,主要是因为研究区植被覆盖度较高且降水量较为丰富,而干度和热度对生态环境状况起负向作用。

(4)该研究通过分析遥感数据可以快速地对大面积区域进行生态环境状况评估,同时运用主成分分析法避免了人为主观因素的影响,提高了分析结果的客观性,对于其他自然保护区进行相关方面的研究有一定的指导意义。但由于传感器和软件设置等原因,在数据处理过程中可能出现一定的误差,导致部分数据不能完全反映区域生态环境状况,有待进一步研究。

4 建议

提升地表覆盖度和防范因地质灾害产生裸地是保护区生态环境治理的关键。保护区要继续实施“天保”和“退耕还林”工程,不断提高区域植被覆盖率；其次,保护区内滑坡、泥石流等地质灾害频发,要充分运用遥感技术对气象及地质灾害易发区进行监控,辅之以必要的工程措施,减轻灾害对生态系统的破坏程度；再次,建议发展旅游业带动周边猕猴桃、药材种植等其他产业的发展,继续降低对生物资源的破坏程度，保护区应同社区制定自然资源管理计划,共同促进社区自然资源的管理。