基于RS和GIS的峨眉山风景区生态风险评价

黄诗曼， 胡庆武*， 李海东， 王少华

1.武汉大学遥感信息工程学院， 湖北 武汉 430079

2.生态环境部南京环境科学研究所， 江苏 南京 210042

近年来，为响应乡村振兴战略，各地旅游小镇迎来建设热潮，积极促进乡村产业发展的同时，也为相应风景区生态安全带来了挑战，威胁了生态旅游的可持续发展[1-3]. 进行区域生态风险评价，促进风景区功能分区，为村镇建设和风景区生态系统健康管理提供参考依据，对风景区绿色可持续发展有重要意义[1-6].

生态风险评价经历了由单一风险源到多风险源、单一风险受体到多受体、局地到区域及景观尺度的发展历程，区域生态风险评价是生态风险评价发展的综合阶段. 区域生态风险评价是在区域尺度上，描述和评估自然因素和人为活动对生态系统及其各个组分所产生不利影响的可能性和大小的过程[7-9]. RS (Remote Sensing，遥感)具有大范围、高时效、短周期、综合性、少约束等特性，GIS (Geographic Information System，地理信息系统)具有强大的空间分析和数据管理能力，二者相结合能为生态风险评价提供极大的支持和辅助作用[1,10-11]. 李耀明等[12]基于北京市遥感影像数据，结合GIS手段，选取多指标并对各指标因素进行打分，最终得到北京市生态风险评估图. ZHANG等[13]基于GIS构建了吉林省生态灾害风险评估模型，并绘制吉林省生态灾害风险和生态灾害风险四因素(风险性、暴露度、脆弱性以及应急响应和恢复能力)的区域分布图. Bathrellos等[14]通过AHP (层次分析法)结合GIS生成滑坡、洪水和地震灾害评估地图，合并后得到代表研究区域Xerias流域内城市发展的潜在适宜性地图，实现了生态风险评价在城乡规划、功能区划方面的应用. Agapiou等[15-16]基于多个遥感影像数据，选取多个风险源并在GIS环境中计算，利用AHP确定风险源权重，最后得到塞浦路斯帕福斯区灾害评估图. LIU等[17]利用2000年、2004年、2008年、2013年、2017年获得的多时相遥感影像数据和GIS构建生态风险指数，识别2000—2017年江苏省沿海湿地的生态风险格局的时空变化. 近年来，人类活动对景观格局的影响越来越受到重视，景观分布在遥感中的表现具有直观性和准确性，它可以准确显示出各种生态影响的空间分布和梯度变化特征，使各种空间分析的手段成为可能. 张雅洲等[11]用景观格局指数评价生态系统的地位和受体的易损性，通过计算干旱、洪涝和水体污染的综合风险概率及生态综合损失度得到南四湖的生态风险值. LI等[18]基于RS、GIS将香格里拉县分为10种景观类型，利用景观格局指数构建景观损失指数，综合考虑风险源范围、发生概率和景观损失指数得到了1990年和2006年香格里拉县生态风险评估的空间分布.

经过近20年的研究发展，生态风险评价的模型和方法步骤得到了相应完善，但研究尺度多集中在省市、流域等尺度上，对于风景区的生态风险评价较少. 综上，该研究利用峨眉山风景区土地覆盖数据、DEM数据与同期遥感影像数据，基于RS和GIS技术，综合考虑自然灾害与人类干扰，自然灾害采用传统的因子权重法，人类干扰利用景观指数来衡量，进行灾害概率和生态易损性分析，依据生态风险度量原理和克里金插值法得到峨眉山风景区2015年生态风险空间分布，确定生态风险高值分布区域，以期为峨眉山风景区生态承载力风险预警提供参考，有助于确定风景区内适应村镇建设、产业发展的区域，指导风景区生态风险管理.

1 研究区概况与研究方法

1.1 研究区概况

峨眉山风景区(103°15′E～103°27′E、29°28′N～29°36′N)处于四川省乐山市峨眉山市城郊，位于四川盆地西南边缘，山体自西向东而卧，素有“峨眉天下秀”之美誉，是世界文化与自然双重遗产、国家重点风景名胜区、5A级景区. 峨眉山风景区位于四川省三大暴雨中心之一的盆地西部，有丰富的降水资源，但降水量时空分布极不均匀，主要集中在4—10月，占全年降水量的91.2%，尤以7—8月为多,期间暴雨频发，枯水期为11月—翌年3月，期间降雨量占全年降雨量的不到10%. 峨眉山地表水径流的区域变化与降水的区域变化相似，一般6—9月径流占全年总径流量的一半以上，10月—翌年2月次之，3—5月再次之，汇集的天然河流均具有坡度陡、流程短，易涨易涸的山溪水特点. 景区内居民人数为2 369人，主要为景区工作人员及僧侣. 图1为峨眉山风景区地理位置示意图. 运用RS、GIS技术和生态风险度量理论对峨眉山风景区旅游生态系统进行生态风险分析，对峨眉山风景区的可持续性开发提供理论支持.

图1 峨眉山风景区地理位置示意

1.2 数据来源及预处理

2015年地表覆盖产品源于国家科技基础条件平台——国家地球系统科学数据共享服务平台(http:www.geodata.cn). 数据集原始数据均为广泛使用的遥感影像数据资源，来源可靠，加工后一级类数据精度在70%以上，分辨率为30 m. DEM数据来源于ASTER GDEM第一版本原始数据处理得到的30 m分辨率数字高程产品. DEM数据和2015年10月20日的Landsat8遥感影像数据都在中国科学院计算机网络信息中心地理空间数据云平台(https:www.gscloud.cn)下载得到，对原始遥感影像进行了辐射定标、大气校正、影像裁剪等预处理；同时，由于研究区地形起伏较大，基于DEM数据对遥感影像进行了地形校正.

1.3 研究方法

风险通常定义为在一定区域和时期内有害事件的发生概率及其导致后果的乘积，风险值是区域生态风险的表征[10,19]. 该研究采用文献[7]归纳的区域生态风险评价步骤，即研究区的界定与分析、受体分析、风险源分析、暴露与危害分析以及综合生态风险评价.

为体现受体和风险源在研究区内的空间异质性，首先划分风险小区，该研究按500 m等间距划分风险小区，共790个，其中不同风险小区对应着不同的风险源和风险受体，而每个风险小区内受到不同种类、不同级别风险源的叠加作用而有着一致的综合风险源. 在进行区域生态风险评价时，往往由于生态风险事件可能造成损失难以计算，得不到较为准确的评价结果. 因此，该研究引入生态易损性指标来评估风险受体在不同风险源作用下的损失程度，最终用综合生态风险值表征区域生态风险的相对大小，每一个风险小区对应一个综合生态风险值，即：

Rk=Pk×Tk

(1)

(2)

(3)

式中：Rk、Pk、Tk分别为第k个风险小区的综合生态风险值、综合风险概率与生态易损值，取值区间均为[0,1]；βi为第i类生态风险源的权重；Pki为第k个风险小区内第i类生态风险源的风险概率；N为风险源数；Ak为第k个风险小区景观总面积，km2；Akj为第k个风险小区的第j类景观组分的面积，km2；Fj为第j类景观组分对应的生态易损值；M为景观类型数. 其中，Pki根据风险小区内同一生态风险源不同风险等级面积占比加权求和得到，即：

(4)

式中：l为风险等级；L为风险等级分级数；Akil为第k个风险小区内第i类生态风险源中风险等级为l的区域面积，km2.

对于整个研究区，综合生态风险值根据各风险小区面积占总面积比例进行加权求和得到，即：

(5)

式中：R为研究区的综合生态风险值；A为研究区总面积，km2；K为风险小区个数.

2 结果与讨论

2.1 受体分析

受体即风险承受者，在生态风险评价中指生态系统中可能受到来自风险源不利作用的组成部分. 在生态风险评价中，受体往往是多个生态系统类型，不同的生态系统类型在区域整体的生态功能方面所发挥的作用不同. 根据土地覆盖数据，研究区被分为农地、森林、草地、灌丛、水体、不透水层、裸地，该研究将这7种土地覆盖类型所代表的生态系统作为自然生态风险和人文生态风险二者共同的风险受体.

2.2 风险源分析

风险源分析是指对区域中可能对生态系统或其组分产生不利作用的干扰进行识别、分析和度量,分为风险识别和风险源描述两部分. 风险识别是指根据评价目的找出具有风险的因素，风险源描述是用区域生态风险事件所发生的概率和强度来描述风险源的危害性. 风险源可分为自然灾害和人类干扰.

2.2.1自然灾害

由于研究区为山岳型风景区，多高山且坡度大，在暴雨作用下，易发生崩塌、滑坡、泥石流等地质灾害，同时区域内降水具有时空分布极不均匀的特征，汛期暴雨频发景区易遭受水力侵蚀，枯水期又给景区带来巨大的干旱压力，因此该研究选择干旱、土壤侵蚀和地质灾害作为主要的自然灾害.

2.2.1.1干旱

单独利用地表温度(Land Surface Temperature, LST)和归一化植被指数(Normalized Difference Vegetation Index，NDVI)进行旱情监测都具有片面性. 基于NDVI-LST二维特征空间进行分析，是干旱监测领域中常见的方法[20-22]. Sandlot等[23]将NDVI-LST二维特征空间简化处理为三角形并提出温度植被干旱指数(Temperature Vegetation Dryness Index，TVDI)，定义如下：

TDVI=(Ts-Tsmin)(Tsmax-Tsmin)

(6)

式中：Ts为给定影像像元s对应的地表温度值，K；Tsmin为影像中与给定影像像元s同等NDVI值对应的最低地表温度，K；Tsmax为影像中给定影像像元s同等NDVI值对应的最高地表温度，K.

利用ENVI依次计算研究区NDVI、LST和TVDI，其中，LST采用辐射传输方程法求取，TVDI利用扩展工具中的TVDI工具计算得到，利用自然断点法将研究区影像像元的TVDI分为5级，分别赋值1、2、3、4、5，值越大表示干旱等级越高，实现了研究区的干旱概率确定.

2.2.1.2土壤侵蚀

土壤侵蚀即水土流失，指在风和水等外力作用下，土壤及其母质被损坏、侵蚀、移动和堆积的过程[24]. 综合考虑影响土壤侵蚀的各种因素，选择地形坡度、土地利用类型、植被覆盖状况构建指标交叉表模型来评价土壤侵蚀敏感性. 参照SL 190—2007《土壤侵蚀分类分级标准》[25]，研究区属于水力侵蚀类型区，水力侵蚀分级标准如表1所示，表中轻度、中度、强度、极强度、剧烈分别赋值为1、2、3、4、5. 其中，坡度数据利用ArcGIS从DEM中提取，植被覆盖度数据在2.2.1.1节计算地表温度时作为中间数据获取.

表1 土壤侵蚀分级

2.2.1.3地质灾害

参考《地质环境评价与地质灾害管理》(2008年出版)，地质灾害的发育与地形地貌、地质构造、地层、水文条件和人类工程活动等有着密切关系[26-29]. 研究[27]指出，地质灾害易发性评价因素主要分为两类，分别为控制因素和影响因素. 峨眉山风景区属于建筑管辖区域，不存在或很少存在矿石开采等工程，且该研究在进行生态风险评价时将人类干扰和自然灾害影响因素分开考虑，故在此不考虑社会属性因素. 此外，由于降雨栅格数据缺失，只获取到国家级气象站点数据，不能利用多个自动站点数据插值得到降雨分布图，因此该研究在计算地质灾害发生概率时将降雨因素忽略.

基于此，根据以往学者研究经验以及研究区特点，在遵循科学性、代表性和易获取性原则上，该研究选取了6个指标作为地质灾害风险性评价指标，各影响因子分级赋值如表2所示，值越大表示风险等级越高. 其中，坡度、坡高、坡向、坡型、河网利用ArcGIS从DEM中提取，NDVI利用遥感影像数据计算得到.

表2 地质灾害影响因子赋值表

该研究利用Yahhp软件中的AHP确定各影响因子权重. 在构建梯级层次结构模型的基础上建立判断矩阵，并通过一致性检验，坡度、坡高、坡向、坡型、NDVI和河网对应权重分别为0.45、0.22、0.06、0.12、0.03、0.12. 利用ArcGIS将研究区各评价因子按对应权重进行信息叠加计算，得到地质灾害风险概率值，采用自然间断点法将其分为五级，等级赋值1、2、3、4、5，值越大表示风险等级越高.

2.2.1.4自然灾害综合概率

各风险源对风险受体的作用强度不同，对形成区域生态风险的作用大小也有差异，根据相关研究[11,18]并结合研究区状况，将干旱、土壤侵蚀、地质灾害对应权重分别设为0.35、0.15、0.50. 为使量纲不同数据具有可比性，将干旱、土壤侵蚀、地质灾害风险等级归一化，再加权求取研究区的自然灾害综合概率.

2.2.2人类干扰

自20世纪70年代以来，峨眉山风景区游客量呈爆发式增长，2015年接待游客量突破330万. 与此同时，旅游开发建设、旅游服务业发展带来的废气、废水、噪声污染等问题给研究区生态系统带来巨大的压力. 该研究以景观生态学理论为基础，选用景观干扰度作为研究区生态系统受到人为干扰程度的定量表征. 根据相关研究[11,13,17,19]，景观干扰度由景观格局指数破碎度、分离度、周长面积分维数倒数加权计算求得. 其中，破碎度指数用于度量土地利用景观斑块空间分布的均匀程度，常用于描述空间异质性和人类干扰程度，值越大表示人类干扰程度越大；分离度指数反映景观类型在地域上的分布特征以及人类活动对景观结构的影响；周长面积分维数表示斑块周长与面积的关系，是反映形状复杂程度的指数，取值范围为1～2，值越靠近2表明斑块形状越繁杂，说明人类干扰程度越小、自然度越高，该研究选用周长面积分维数倒数作为指标之一，分维数倒数值越大表示景观类型受到的人类干扰越大.

结合已有研究[13,17,30]，破碎度、分离度、周长面积分维数倒数权重分别设为0.5、0.3、0.2，利用Fragstats软件求得3个景观格局指数，加权后归一化得到景观干扰度，值越大表示该景观类型受到的人为干扰程度越大. 经计算，研究区农地、森林、草地、灌丛、水体、不透水层、裸地对应的景观干扰度分别为0.801、0.000、0.988、0.520、1.000、0.156、0.495.

2.3 暴露与危害分析

暴露分析研究各风险源在区域中与风险受体的接触暴露关系. 与之相关联的危害分析则要确定风险源对生态系统及其风险受体的损害程度. 不同景观类型代表的不同生态系统在区域内所处的位置不同，在维护生物多样性、完善生态系统的结构和功能、促进景观结构自然演替等方面的作用有一定差异；面对相同的外来干扰时，不同生态系统的损失度不同. 该研究用生态易损性指数表示不同景观类型所代表的生态系统对外部干扰的抵抗能力.

在确定生态易损性指数时，常将其与土地覆盖类型相联系，常用的方法是直接赋值归一化处理和AHP. 直接赋值归一化处理较AHP有较大的主观性，因此该研究利用AHP确定各景观类型的生态易损性指数. 经计算，农地、森林、草地、灌丛、水体、不透水层、裸地生态易损性指数分别为0.201、0.044、0.065、0.093、0.246、0.027、0.324.

2.4 综合生态风险评价

综合生态风险评价要结合受体分析、风险源分析、暴露和危害分析，得到区域范围生态风险值，从而为区域生态风险管理提供理论依据.

根据生态风险值计算原理分别计算研究区内各风险小区的自然灾害风险值、人类干扰风险值，二者之和即为风险小区对应的综合生态风险值，将其作为风险小区中心点属性值，利用克里金插值法求得研究区生态风险空间分布，利用自然断点法分为五类，结果如图2所示. 表3为生态风险综合评价表，根据式(3)～(5)对各风险小区的风险值加权求和计算得到.

由表3可见，峨眉山风景区自然灾害风险值、人类干扰风险值分别为 0.025 8、0.002 8，综合生态风险值为0.028 6. 风景区主要受自然灾害影响，其贡献率高达89.9%，但人类干扰给其带来的影响仍不容忽视. 自然灾害中，地质灾害是影响最大的自然灾害，其次是干旱，最后是土壤侵蚀，贡献率依次为49.3%、40.9%、9.7%. 峨眉山风景区降雨集中，全年降水主要集中在6—9月，该研究计算TVDI所用遥感影像为2015年10月获取，这可能是导致干旱贡献率较大的原因之一.

从空间分布上看，峨眉山风景区综合生态风险呈条带状分布，越往东风险值越大，高值聚集在报国寺景区附近(见图2). 与风景区规划总图对比发现，低风险区全部位于峨眉山风景区的一级(A)保护区内，该区域仅可开展科考和资源管理活动. 高风险区主要位于风景区内二、三级保护区，其中，二级保护区包括北部居民点分布较为集中的区域以及金顶至万佛顶一线的游览设施集中的区域，三级保护区包括黄湾镇张坝村平坝地区、报国村旅游服务设施集中的区域、西部龙洞村等. 由图2可见：风景区的自然灾害风险空间分布与综合生态风险空间分布基本一致，说明自然灾害是导致风景区综合生态风险的主要因素. 研究区人文风险整体处于低和较低风险区，仅东部地区有少量较高和高风险区，说明风景区受人类干扰较小，与表3中结果一致；但人类干扰带来的风险不容忽视，黄湾乡旅游镇附近自然灾害风险值低于清音阁景区，但其综合生态风险值却高于清音阁景区，主要原因是该区域人类干扰风险值较大. 从风景区绿色可持续发展角度来看，研究区内行政村多集中在神水阁、清音阁、万年寺景区，且植被较稀疏，属于综合生态风险高风险区，因此需要严格控制该区域内村庄人口规模和建设规模，依规划进行环境整治，加强村庄绿化，条件成熟时可搬迁. 其中，神水阁、万年寺景区干旱风险概率较高，清音阁景区土壤侵蚀和地质灾害的风险概率较高，故峨眉山自然灾害高风险区多集中在该区域；四季坪景区东部干旱风险概率较大导致该区域内存在少量自然灾害高风险区. 针对高干旱风险区域，应增加水源地的涵养能力，加强林草保护；针对清音阁景区，可加大生态林和经济林建设，可达到保持水土、预防地质灾害和降低生态风险的目的. 而峨眉山东部临近省道，是黄湾乡旅游镇所在地，景区停车场多集中于此，外部为城市发展区，受人类干扰较大，但自然灾害风险值较小使得综合生态风险值较小，因此作为旅游服务基地区域，需要依据相关规划进行游览设施建设和村庄建设，应统筹规划用地，优化建设布局.

图2 峨眉山风景区生态风险评价分级

表3 峨眉山风景区生态风险综合评价结果

3 结论

a) 峨眉山风景区综合生态风险值较小，主要受自然灾害影响，其贡献率为89.9%，人类干扰的贡献率为10.1%，但人类活动给风景区带来的影响同样不可忽视. 自然灾害中，风景区主要受到地质灾害、干旱的影响，贡献率分别为49.3%、40.9%. 峨眉山风景区综合生态风险呈条带状分布，越往东综合生态风险值越大，自然灾害风险空间分布与其基本一致.

b) 为降低风景区生态风险、促进绿色可持续发展，对于行政村所在区域，需要严格控制区域内村庄人口规模和建设规模，依规划进行环境整治，加强村庄绿化. 对于高干旱风险区域，如神水阁、万年寺景区，应增加水源地的涵养能力，加强林草保护. 对于清音阁景区，可加大生态林和经济林建设. 对于景区东部黄湾乡，在进行游览设施建设和村庄建设时需要严格依据相关规划.