澜沧江流域重大水电工程扰动灾害风险评价*

李 萍 李丽慧 刘昊磾 祁生文 邓清海

(①中国科学院地质与地球物理研究所，中国科学院页岩气与地质工程重点实验室，北京 100029，中国)(②山东科技大学地球科学与工程学院，青岛 266590，中国)(③中国科学院大学，地球与行星科学学院，北京 100049，中国)(④中国地质大学(北京)工程技术学院，北京 100083，中国)

0 引 言

地质灾害风险评价可定义为对特定影响因子造成暴露于该因子的单体或区域地质灾害发生的概率及对人类社会产生危害的程度、时间或性质进行定量描述的系统过程，是在危险性评价和易损性评价的基础上对风险区遭受地质灾害的可能性和后果进行定量分析与评价，以及采取相应措施来降低风险可能性的一门集自然属性和社会属性并重的交叉学科(向喜琼等，2000，齐信等，2012)。

澜沧江是一条重要的国际河流，在中国境内天然落差4583m，水力资源丰富，为充分利用澜沧江水能资源，规划建设23级梯级水电站，截至目前，已建的水电站有乌弄龙、里底、黄登、大华侨、苗尾、功果桥、小湾、漫湾、大朝山、糯扎渡、景洪水电站，在建的有托巴水电站，未建的水电站有侧格、约龙、卡贡、班达、如美、邦多、古学、曲孜卡、古水、橄榄坝、勐松水电站，各水电站参数如表1所示。

表1 澜沧江流域水电站参数

在库区形成及运行过程中，随着库区水位的大幅涨落，库岸坡体内的渗流场产生巨大变化，导致边坡失稳，地质灾害多发，如塌岸、古滑坡复活、建筑物受损等(图1)，不仅会影响水电站的正常运行，还会对周边居民的生命财产造成威胁。2009年7月20日发生在小湾电站库区的荒田大滑坡，造成14人死亡(失踪)和两艘船沉江失踪(王飞，2018)。水库蓄水后库岸原有自然状态改变，岸坡稳定性降低，导致了灾害多发，已经引起了库区管理者和当地政府重视，为此，开展重大水工程扰动灾害尤其是风险区划工作显得尤为迫切。

图1 研究区内地质灾害

国内一些学者和专业技术人员针对澜沧江水电站工程区做了大量的工作。针对区域性地质灾害，姚鑫等(2020)针对澜沧江流域高山峡谷区活动性滑坡提出了相位共振增强InSAR变形观测技术(PRE-InSAR)，根据解译的活动性滑坡数据，揭示了其空间发育规律和主控因素；周振凯等(2020)针对澜沧江流域梅里雪山段，综合使用PALSAR-1、PALSAR-2和Sentinel-1这3种SAR数据，对工作区进行全面InSAR观测，使用了D-InSAR、IPTA-InSAR和Stacking多种InSAR计算方法，分析了不同SAR数据和InSAR方法在研究区的适用性，解译并分析了重大单体滑坡变形体分布。针对水电站蓄水后工程区边坡稳定性问题，石崇等(2010)、赵永辉等(2016)、毛峰(2020)先后采用FLAC3D、UDEC软件对蓄水作用下争岗滑坡变形破坏特征进行数值分析，分析其稳定性；尹云坤等(2015)、安可君等(2020)根据边坡性质、库岸变形特征分析了小湾、糯扎渡水库蓄水后的稳定性问题；针对地质灾害危险性分区，苏鹏程等(2014)分析了澜沧江流域滑坡泥石流形成的地质、地貌、地震、降雨以及人类活动等影响因素，并进行灾害危险性分区；李浩(2012)分别采用模糊数学综合评判法、层次分析法(AHP)和信息量计算法3种方法对小湾水库库岸边坡进行危险性评价；成旭鹏(2017)分析了小湾库区消落带内崩塌滑坡的分布特征，采用信息量法对库区消落带进行危险性评价；王飞(2018)利用层次分析法对小湾水库库岸崩塌滑坡进行了危险性评价并选取人口密度、经济密度作为承灾体进行风险评估。

澜沧江流域梯级水库的建立诱发了诸多古滑坡和堆积体复活变形，从目前研究来看，前人对澜沧江的研究多集中于对单体典型灾害点在蓄水扰动下的稳定性分析，或对某一水电站工程区进行稳定性分析或危险性评价，结合多个水电站进行全区域内扰动灾害调查的研究较少。基于此，本文以干流23个梯级水电站扰动灾害为研究对象，采用信息量模型、基于层次分析法的加权信息量模型，开展澜沧江流域重大水电工程区风险评价，对澜沧江流域在建或已建的水电工程区防灾减灾工作提供帮助。

1 研究区概况及数据源

1.1 研究区概况

澜沧江起源于青藏高原中部，青海省扎多县唐古拉山北侧的查家日马峰南坡的扎曲，主要流经青海省、西藏自治区、云南省3省。上游与通天河相邻，西部与怒江分界线为他念他翁山、怒山(南段碧罗雪山)等山脊线，东部与金沙江、李仙江的分水岭分别为云岭、无量山，是西南横断山区的重要河流之一，干流全长2179km，为世界第7长河。澜沧江流域从北到南共纵跨13个纬度带，流域自北向南呈条带状，上下游较宽阔，中游较为狭窄，涵盖了高山峡谷、中低山宽谷、冲积平原等多种地貌类型(苏鹏程等，2014)(图2)。

图2 研究区概况图

由于研究区南北跨度大，因此流域内气候差异很大，随纬度降低，年平均降雨量以及气温呈现增加的趋势。从行政范围上看，澜沧江干流分别以昌都、功果桥为界分为上、中、下游。上游地区地处青藏高原区，在气候上属于高原温带以及高寒气候，并且由于地势原因，气候具有明显的垂直变化，中游在气候带上属亚热带，澜沧江下游属亚热带或热带气候。

地质构造上，澜沧江流域属于特提斯-喜马拉雅构造域，经历了漫长的地质演变过程，由东印支板块、昌都-兰坪-思茅-南邦中间板块、毛淡棉-金边移动板块，以及掸泰马板块连接而成，受板块运动和地壳变化影响，形成了一系列断裂构造。研究区内地层以南北向走向为主，有元古界、古生界和中生界地层的砂岩、泥岩、变质岩；泥盆系、二叠系、三叠系的碳酸盐岩、砂岩、板岩、片岩、片麻岩、大理岩；沿澜沧江构造带分布少量的上侏罗系、白垩系、古近系、新近系的砂岩和泥岩等(姚鑫等，2020)。

1.2 数据源

本文灾害点数据通过野外调查和基于Google Earth解译获取，共计有崩塌73个，滑坡207个，泥石流6个，分布状况见图3。高程、坡度、坡向数据基于地理空间数据云平台GDEMV3,30m的 DEM 数据，工程地质岩组来源于1︰50万地质图资料，断层线密度来源于1︰20万地质图资料，地表植被覆盖信息来源于30m分辨率的 Landsat8 影像，年均降雨量数据来源于全国降雨统计资料，地震动加速度来自于《中国地震动参数区划图》(GB18306-2015)，人口密度来源于Worldpop2020人口数据集，土地覆盖数据来源于中国土地覆盖遥感监测数据库，道路数据来源于全国路网图。

图3 研究区内主要地质灾害分布图

2 地质灾害危险性评价

2.1 研究方法

当前常见的地质灾害危险评价方法大致可以分为主观推断分析方法，如层次分析法；统计分析模型，如信息量模型、证据权模型、确定性系数方法、判别分析、线性回归和逻辑回归等本质上都属于统计分析；机器学习算法，如支持向量机、多层感知器模型、决策树模型等确定性模型等(胡瑞林等，2013；Chen et al.，2016；杨光等，2019)。由于物理意义明确、操作简单，信息量模型在地质灾害危险性评价研究和实践中得到了广泛的应用。但是信息量模型的统计只反映了特定因子组合对地质灾害的影响，忽略了不同因子对灾害发育的差异。层次分析法作为一种基于专家打分的主观赋权的方法，能够反映研究者对地质灾害的直观认识，但是受主观影响较大，所以经常与其他模型相结合，组成加权的统计模型(单博，2014)。

本文采用信息量模型、加权信息量模型分别进行计算。其中加权信息量的权重值是基于层次分析法获得的。对同一判断矩阵进行计算时，采用不同方法得到的权重可能具有很大差异，不可靠的计算方式会导致最终赋权不准确。为得到最终合理可靠的权重值，本文放弃常规思路中单一的计算方法，改为选取和积法、方根法、特征向量法3种权重计算方法的平均值作为最终权重，再赋值到加权信息量模型中。

2.1.1 信息量模型

信息量法源于信息理论，是一种对数化随机分布概率以解决度量问题的定量统计预测方法。由美国数学家Shannon应用概率论知识和逻辑方法推导出的计算公式，20世纪80年代，晏同珍教授首次将信息论引入到滑坡灾害预测研究中，后来被广泛应用到灾害评价领域(晏同珍等，1989)。其中心思想是基于各影响因子对地质灾害所贡献的信息量值，将信息量值反应在地质灾害发生的难易程度。可以定量描述地质灾害的危险性，一般来说，信息量值越大，地质灾害越容易发生。具体计算公式如下：

(1)

式中：IAj→B为对应因素A、j状态(或区间)下地质灾害B发生的信息量；Nj为对应因素A、j状态(或区间)下地质灾害分布的单元数；N为调查区已知有地质灾害分布的单元总数；Sj为对应因素A、j状态(或区间)分布的单元数；S为调查区单元总数。

2.1.2 加权信息量模型

常规信息量模型可以理解为，所有因素在地质灾害发育过程中贡献值是一样的，即将所有元素的权重赋为1，而加权信息量模型则是根据客观事实，基于专家经验，综合考虑不同影响因子对地质灾害发生影响程度的差异，采用层次分析法对信息量值进行赋权，其表达式如下：

(2)

式中：LHI为危险性指数；Weighti为第i个因子的权重值；Ii为第i个因子的信息量值。

层次分析法由美国学者Saaty et al.(1983)于20世纪70年代提出，其基本原理是基于专家经验比较两个因子间的重要程度，用数字 1～9 及其倒数作为标度，如表2所示，构建判断矩阵，在此基础上进行定性和定量分析，从而确定不同评价因子的权重。

表2 判断矩阵 1～9 级标度取值含义

一般情况下，可采用的权重计算方法包括：和积法、方根法以及特征向量法。

(1)和积法

a.对判断矩阵A的列向量进行归一化处理，得：

(3)

(4)

(5)

d.计算最大特征值

(6)

(2)方根法

a.对判断矩阵A的行向量进行相乘，得：

(7)

b.对向量，Wi进行开方处理，得：

(8)

(9)

d.计算最大特征值

(10)

(3)特征向量法

Aw=λmaxw

(11)

w为判断矩阵A的最大特征根λmax所对应的特征向量。

为避免其他因素对结果产生干扰，保证判断矩阵的准确度，通常情况下，用一致性指标CI的值对结果进行检验，当n>3时，引入一致性指数RI，RI指标取值如表3所示，一般要求CR<0.1，才能表示判断矩阵通过一致性检验(刘丽娜等，2014)。

表3 平均随机一致性指标RI取值表

(12)

(13)

式中：λmax为判断矩阵的最大特征值；n为判断矩阵的阶数。

2.2 评价因子的选取及分级

地质灾害发育受自身地质条件和外界环境因子共同控制和影响。前人研究表明，滑坡灾害的发育主要受地形地貌、基础地质、水文环境、植被覆盖、降雨、地震和人类工程活动等多种因素的影响(Guo et al.，2015；王世宝等，2022)。根据获取的资料，参考当前被广泛应用于地质灾害危险性评价中的因子，选取高程、坡度、坡向、植被归一化指数(NDVI)、距库区的距离、断裂带密度、工程地质岩组、年均降雨量和地震动加速度9个影响因子构建地质灾害危险评价体系。结合研究区地质灾害的发育特征，对选取的地质灾害评价因子进行了子类划分，各评价因子的子类划分结果见图4。

2.3 信息量值计算及权重的确定

根据地质灾害的解译结果和野外调查数据为基础数据，对各评价因子子类的信息量值进行了计算，其计算结果见表4。

表4 影响因子分级及信息量值

由表4可知，在<4500m的海拔区域范围内，人类活动密集，受人类对自然环境改造的影响，地质灾害发育，>4500m以上区域主要为高原夷平面，冰川、冰碛物分布区，不利于地质灾害的发育。坡度在20°～40°时，灾害发育面积较大，在这个范围之外，要么坡度太陡发育冰川，要么地势平缓不具备灾害发育的势能。E/W和NE/SW坡向滑坡最发育，与构造行迹和澜沧江的总体流向一致，SE和S坡由于日照充分和土壤水分低，地质灾害不发育。植被覆盖度低的区域多为高海拔地区，长年冻土，植被覆盖度高的地区为下游河谷平原区，均不利于地质灾害发育。地质灾害主要发育在较软弱岩组和坚硬岩组，较软弱岩组主要是侏罗-白垩的“红层”泥岩、砂岩，有利于滑坡的发育。灾害的分布与降水量无直接的正相关，在雨量最少的干热河谷地带反而是滑坡分布最密集的区域。断裂带密度大的区域，坚硬岩组较为破碎，易发生崩塌。PGA值在0.20g区域内，地质灾害发育程度更高。地质灾害的发育具有显著的“库区距离效应”，与库区距离呈负相关，受库区距离的影响最大，因此库区距离因子所占的权重也最大。

根据通用的1～9标度方法，对各因子进行两两比较建立判断矩阵，进行层次排序，分别用和积法、方根法以及特征向量法确定各因子的权重值，并求三者的平均值，权重计算结果如表5所示：

表5 不同方法所得权重值

根据平均随机一致性指标取值表，当n=9时，RI值为1.460，根据式(12)、式(13)对判断矩阵进行一致性检验，如表6所示，可以看出，4种计算方法所得的CR值均小于 0.1，表明该判别矩阵具有满意的一致性。

表6 一致性检验结果

2.4 评价结果及分析

将评价因子子类的信息量值进行叠加，得出单一信息量模型的地质灾害危险性指数；将子类信息量值与评价因子权重值代入式(2)，得到加权信息量模型的危险性指数。引入ROC(Receiver Operating characteristic Curve)对不同方法得到的危险性指数图进行准确性检验，ROC曲线下的面积即为AUC(Area Under Curve)(Fawcett，2006)，其值在0～1之间，数值越大，代表预测准确率越高，当AUC值大于0.9时，则模型具有极高的预测精度(高泽民等，2021；文洪等，2021)。

图5为不同方法得到危险性指数的ROC曲线图，横坐标为研究区危险性等级从高到低的区域面积累加百分比，纵坐标为各危险性等级中对应的真实历史灾害面积的累加百分比。经计算单一信息量模型的AUC值为0.860，加权信息量模型中和积法的AUC值为0.888，方根法的AUC值为0.893，特征向量法的AUC值为0.893，平均值法的AUC值为0.901。5种模型对地质灾害危险性评价均具有较好的结果，但相比于传统的信息量模型，加权后的信息量模型评价结果更加准确，同时，采用3种权重计算方法平均值作为最终权重的方式，进一步提升了预测的可信度。

图5 不同评价模型下的ROC曲线

根据加权信息量模型得到的危险性评价指数图，对结果进行归一化，划分为极高危险区(0.75，1]、高危险区(0.5，0.75]、中危险区(0.25，0.5]和低危险区(0，0.5]4类，从图6可以看出，极高危险区和高危险区沿澜沧江干流呈条带状分布，其中极高危险区主要分布在中游高山峡谷区，受断裂带的影响，小湾和糯扎渡库区的危险性也极高；高危险区主要分布在澜沧江下游以及大华桥—功果桥水电站段；中危险区分布在昌都—黄登段；低风险区分布在黄登—橄榄坝电站远离库区的区域。

图6 地质灾害危险性分区图

3 地质灾害易损性评价

易损性是指承灾体可能遭受地质灾害破坏的严重程度。易损性评价是地质灾害风险评价的重要内容，目前用于地质灾害易损性评价的模型非常多，主要包含定性与定量评价两种。本次评价在调查区的地质环境特征、受到工程扰动的地质灾害发育特征的基础上，结合实际情况，选取人口密度(人·km-2)、道路、水电站(距水电站的距离/km)、土地覆盖类型、单位面积GDP(万·km-2)5个承灾体按表5进行分级赋值，分级标准及结果见表7和图7。

表7 地质灾害易损性评价因子分级标准

图7 承灾体易损性分级结果

基于层次分析法对选取的承灾体构造判断矩阵，确定权重，定量分析承灾体的抗灾能力，权重值见表8，对结果进行一致性检验，经计算λmax=5.12，CR=0.026<0.1，结果满足一致性检验要求。

表8 易损性评价因子判断矩阵及权重

对不同类型承灾体按式(14)进行叠加，根据自然断点法将结果划分为极高易损区，高易损区，中易损区和低易损区4类，划分结果如图8所示。

图8 地质灾害易损性分区图

(14)

式中：LVI为易损性指数；Weighti为承灾体的权重值；Xi为各分级的赋值。

受地形影响，中游高山峡谷段人口分布稀少，经济条件相对落后，因此地质灾害发生时，所造成的人民生命财产和经济损失较小，除水电站附近，多为中、低易损区；下游地势平坦，交通、经济发达，人口分布密集，故高、极高易损区分布面积大。

4 地质灾害风险性评价

地质灾害风险是指地质灾害产生不良后果的可能性，包括地质灾害发生破坏的可能性及其产生的后果(损失)两个方面(Varnes，1984)。本文采用如下的通用风险评估模型开展澜沧江流域地质灾害风险评价(Varnes，1984；Fell, 1994；兰恒星等，2021)。

(15)

式中：LRI为地质灾害风险度；LHI为地质灾害危险度；LVI为地质灾害易损度。

采用风险分级矩阵(表9)来描述地质灾害危险性和承灾体易损性在不同组合状态下的风险性。利用ArcGIS栅格计算器功能将地质灾害危险性分区图层和承灾体易损性分区图层按式(15)进行信息融合，最终完成地质灾害风险性评价图，划分结果如图9所示。

表9 地质灾害风险性分级矩阵

图9 地质灾害风险性分区图

对风险性分区面积进行统计，统计结果如表10所示，低风险区面积为40746.02km2，占全区总面积55.09%；中风险区面积18306.42km2，占全区总面积的24.75%；高风险区面积12524.55km2，占全区总面积的16.93%；极高风险区面积2388.34km2；占全区总面积的3.23%。

表10 地质灾害风险性分区说明表

5 结 论

本文针对澜沧江流域重大水电工程区地质灾害扰动问题，对信息量模型采用不同的加权方式，构建地质灾害危险性评价模型，考虑地质灾害发生时造成人员伤亡及经济损失，构建研究区地质灾害易损性评价模型，在此基础上，进行地质灾害风险划分，本文的具体研究成果如下：

(1)选用层次分析法3种权重计算方法的平均值作为最终权重对信息量模型进行加权，提高了模型的准确度；从危险性评价图可以看出，极高危险区和高危险区沿澜沧江干流呈条带状分布，其中极高危险区主要分布在黄登水电站上游水电站库区，受断裂带的影响，小湾和糯扎渡库区的危险性也极高；高危险区主要分布在澜沧江下游以及大华桥、苗尾、功果桥水电站段；中危险区分布在托巴上游远离库区的区域；低风险区分布在大华桥、苗尾、功果桥、景洪、橄榄坝、勐松水电远离库区的区域。

(2)选取人口密度、道路、水电站、土地覆盖类型、GDP 5个因子基于层次分析法，进行地质灾害易损性评价，评价结果显示，极高和高易损区主要分布在水电站附近以及经济发达，人口分布密集区域。

(3)基于地质灾害危险性和易损性评价分区结果，采用风险分级矩阵对研究区风险性进行划分，结果表示：水电站坝址附近以及约龙、里底、托巴、黄登、大华桥、苗尾、功果桥、小湾、漫湾、大朝山糯扎渡、景洪、橄榄坝、勐松库区均为高风险区或极高风险区；中风险区主要分布在侧格、卡贡、班达、如美、邦多、古学、曲孜卡、古水和乌弄龙电站库区以及里底、托巴电站远离库区的区域；低风险区分布在远离库区的其他区域。