冰碛湖溃决危险区铁路沿河线方案评价方法

张 聪，姚令侃,2,3，黄艺丹,2,3，邱燕玲，谭 礼

（1.西南交通大学土木工程学院，四川 成都 610031；2.西南交通大学高速铁路线路工程教育部重点实验室，四川 成都 610031；3.西南交通大学陆地交通地质灾害防治技术国家工程实验室，四川 成都 610031；4.中铁二院工程集团有限责任公司，四川 成都 610031）

青藏高原南缘是世界上中低纬度地区最大的现代冰川分布区，发育了世界上罕有的海洋性山谷冰川，并由于冰川的活动形成了大量冰碛湖[1].海洋性冰川具有消融量大、活动性强的特点，使得青藏高原南缘，尤其是喜马拉雅山脉，成为全球冰碛湖溃决灾害事件最密集的地区之一[1-2].冰碛湖位于高海拔地区，一旦溃决，其溃决洪水在重力作用下往往给下游的生命财产和交通设施带来难以估量的损失.例如，1981年喜马拉雅中部波曲河流域的次仁玛错冰碛湖溃决事件，使得中尼公路超过55 km路段损毁，交通中断，并导致数百人死亡[1].因此，预测、评估冰碛湖溃决风险对冰碛湖风险调控具有重大的现实和指导意义.

长期以来，随着国内外学者对冰碛湖溃决灾害研究的不断深入，从溃决的激发因素与溃坝机制出发，对冰碛湖危险性进行评估，取得了一系列成果.徐道明等[3]通过对我国西藏地区溃决冰碛湖的考察，提出了从坝体特征、气候因素判别危险冰碛湖的几类指标.吕儒仁等[4]分析了我国西藏冰碛湖溃决特点，统计了有利于溃决事件发生的7类指标.Huggel等[5]利用遥感影像从冰碛湖、坝、激发条件、下游洪水路径等方面提出冰碛湖潜在危险性评判的三级准则和相应指标.基于上述指标体系可对冰碛湖溃决危险性作出定性评价.

从冰碛湖溃决危险性定量评估角度，吕儒仁等[4]将冰碛湖与母冰川变坡点（或裂隙面）下游冰体的体积比定义为危险性指数，提出了评估由冰崩涌浪造成坝体溃决的定量化指标.陈晓清等[6]从喜马拉雅地区冰碛湖相关因子中选取7个参数作为直接判别的依据，结合冰碛湖溃决危险指数建立起一种定性-半定量的个体冰碛湖危险性评价方法.舒有锋等[7]从冰川、湖盆、终碛堤、气候等特征中提取冰碛湖溃决危险性评价指标，通过粗糙集理论确定各指标权重，建立了冰碛湖溃决危险性评价方法.王欣等[8]将冰碛湖溃决机制的定性描述与溃决概率进行转换，建立了冰碛湖溃决的事故树模型，对我国青藏高原以南冰碛湖溃决危险性进行了分析.

冰碛湖溃决通常是一种因素主导，多种因素共同作用的具有高度不确定性的复杂物理过程，不同诱因组合可能导致不同溃坝模式以及严重程度不同的溃决洪水灾害.上述评价方法对冰碛湖溃决这类高度不确定性、多状态、共因导致的失效问题针对性不强.贝叶斯网络（B-Ns）分析法具有良好的适用性[9-10]，能够很好地弥补以上评价方法的不足，同时可充分利用已有案例作为先验知识，提高预测结果的客观性和准确性，但目前未见其应用于冰碛湖溃决危险性分析.

中巴铁路、中尼铁路等穿越冰碛湖分布区的口岸铁路均列入国际通道的规划中.这些跨境铁路将穿越喜马拉雅山脉，冰碛湖溃决已成为对铁路线路方案起控制作用的重大灾害类型，而冰川覆盖区是我国铁路没有建设经验的地貌单元.冰碛湖一般位于极高海拔区，难以通过工程手段处治，因此以在选线阶段开展风险调控最为主动.鉴此，以冰碛湖溃决危险区铁路方案风险评价为目标，首先对青藏高原南缘74例冰碛湖溃决灾害实例溃决机制进行统计分析；在此基础上针对冰碛湖溃决影响因素众多、关系复杂且具有高度不确定性的特点，利用多态BNs建立了冰碛湖溃决概率预测模型；依据冰碛湖溃决危险性评价结果与沿河线河谷地貌形态特征提出了针对冰碛湖溃决危险性的铁路选线方案评价方法；最后，以中尼铁路跨喜马拉雅段吉隆、樟木局部走向方案比选为例说明评价方法的作业程式.为冰碛湖溃决危险区铁路方案风险评价提供具普适性的技术评价体系.

1 多态B-Ns冰碛湖溃决危险性评估方法

B-Ns能够有效地在不确定性环境中实现概率推算与结果预测，利用其进行危险性分析的一般流程是构建B-Ns模型、确定节点变量及相应的概率分布，然后进行推理计算.其中节点变量及网络构建均依赖于已有案例作为先验知识，因此，本节首先对青藏高原南缘冰碛湖溃决灾害案例进行梳理.

1.1 青藏高原南缘冰碛湖溃决案例

结合实地考察、影像资料与历史资料[11-17]，本文共梳理了青藏高原南缘74例冰碛湖溃决灾害事件（图1，附加材料表S1）.可知，该区域冰碛湖溃决事件主要发生在喜马拉雅中、东部以及藏东南地区.表S1中，对冰碛湖溃决机制有记载的案例为62例，统计知其溃决机制可以大致分为5类（表1），其中冰崩涌浪漫顶是最主要溃决机制.此外，有8例为多种溃决机制共同作用导致.

图1 青藏高原南缘冰川、冰湖及冰碛湖溃决灾害分布Fig.1 Distributions of glaciers, glacier lakes and morainedammed lake outburst events in the southern Tibet Plateau

1.2 基于多态B-Ns冰碛湖溃决概率计算方法

1.2.1 构建多态B-Ns模型

构建B-Ns模型的首要工作是确定B-Ns节点变量.将表1中冰碛坝失效的5种主要模式作为一级指标，然后结合冰碛湖溃决灾害形成的基本条件，分别对5个一级指标建立各自的二级分析指标.

表1 青藏高原南缘冰碛湖主要溃决机制统计Tab.1 Main outburst mechanisms of the moraine-dammed lakes in the southern Tibet Plateau

1）冰崩涌浪漫顶（y1）：冰川末端的冰体在气候等激发下，突然脱离母冰川进入下游冰碛湖激发冰碛湖涌浪，同时导致冰碛湖水位骤涨，漫顶冲刷冰碛坝致使坝体溃决.冰崩的断裂面一般位于冰川坡度突变处或者冰舌裂隙发育部位，因此，该部位下游冰体可称为危险冰体.y1的影响因子分析如下：

① 危险冰体坡度（x1）决定着冰体崩落起动的难易程度.由青藏高原南缘的冰碛湖溃决案例统计可知，发生冰崩灾害的冰川末端坡度变化较大，但当冰川末端的危险冰体坡度超过8° 时，容易导致冰崩灾害的发生[6].

② 危险冰体指数（x2）为危险冰体体积与冰碛湖湖水体积的比值.危险冰体滑入冰碛湖后产生涌浪，同时导致冰碛湖水位骤涨.二者综合形成的水头高度越大，坝体溃口处起动泥沙颗粒也越容易起动，相应发生溃决可能性越高[18].因此可以采用危险冰体指数表示其危险性.参照冰碛湖溃决案例，溃决前危险冰体指数一般位于0.11 ～ 0.52[6].

③ 冰舌冰碛湖距离（x3）.冰崩涌浪的首要形成条件是冰川冰碛湖碛处于邻近状态.参照冰碛湖溃决案例，已溃决的冰碛湖x3值大多小于1 km.

④ 冰川裂隙发育程度（x4）.在气候反常年份，冰川容易在裂隙发育面发生断裂，因此冰川裂隙发育程度与冰崩灾害的易发性呈现正相关关系.

⑤ 年际气温变异系数（x5）.气温突增使得冰川底部融水增加，摩阻减小，同时导致冰川内部应力差增加，容易引发冰崩.采用年际气温变异系数代表气温突变程度，系数值越大，气温年际变化越剧烈，也越有利于灾害的形成.由西藏地区冰碛湖溃决区邻近气候站温度统计知，冰碛湖溃决区的气温变异系数值位于0.23 ～ 0.27.

2）管涌（y2）：坝体厚度和坝体物质组成是管涌的影响因素[1].目前，对已溃决冰碛湖统计发现大部分冰碛坝顶宽（x6）不超过60 m.其次，堰塞坝体的物质组成影响着渗流稳定，特别是当冰碛坝体内部含有埋藏冰时，随着气温升高，埋藏冰融化容易引发冰碛坝的管涌现象.当前并无判别冰碛坝内是否存在埋藏冰的方法，但是徐道明等[3]调查后认为，在小冰期末时期所形成的冰碛堰塞湖的坝体中，埋藏冰是普遍存在的，且受气温变化的影响较为明显.据此，可以采用冰碛坝顶宽、坝体物质状态（x7）与年际气温变异系数对管涌的危险性进行分析.

3）洪水漫顶（y3）：在气候异常年份，冰川融水与降雨导致的超标洪水对冰碛湖补给超过溢流、蒸发的损耗时，冰碛湖容易发生漫流型溃决.采用年际降雨变异系数表征降雨的不均匀度，则冰碛湖漫流溃决影响因子可以采用冰碛湖满溢程度（x8）、年际降雨变差系数（x9）与年际气温变异系数表示.

4）冰碛坝坍塌（y4）：冰碛坝多为松散岩屑或积雪（冰）夹带岩屑，稳定性较差的自然堆积体.坝体的坡度越大，其稳定性越差，在内、外部因素作用下容易坍塌[1].对青藏高原南缘冰碛湖溃决案例调查表明，溃决冰碛堰塞坝背水面的坡度多位于20° ～ 33°.其次，坝顶若存在漫流下切槽道，则坝体容易在溯源侵蚀的作用下坍塌.此外，冰碛坝内埋藏冰消融也可能导致坝体坍塌.因此，可以采用背水坡度（x10）与坝顶漫流槽道分布状态（x11）以及年际气温变异系数分析.

5）外界动力激发（y5）：外界动力主要为地震作用.地震力一方面促使冰碛湖附近岩体、冰雪崩滑激发涌浪，另一方面地震直接激发水体涌浪导致漫顶现象的发生.Ai等[18]利用地震振动台模拟地震与滑坡共同作用下复合涌浪实验发现，地震峰值加速度（PGA，x12）与涌浪波高正相关，当x12> 0.25g（超越概率等于10%）时，涌浪波高可达到冰碛坝颗粒起动临界条件，可能引发冰碛湖溃决灾害.

确定上述5种溃坝机制作为中间节点变量，根据节点变量之间的因果关系建立B-Ns拓扑结构，如图2所示，其中：T为叶节点事件，即冰碛堰塞湖溃决；yj为中间节点事件 （j=1,2,···,m）；xi为根节点事件 （i=1,2,···,n）.依据有利于形成灾害的节点状态分界值，对B-Ns根节点变量进行离散化处理，得到各根节点变量状态及取值范围，如表2所示.表中：状态值0、1分别表示安全、危险状态.

图2 冰碛湖溃决B-Ns模型Fig.2 B-Ns model of moraine-dammedlakes outburstrisk

表2 根节点因素多区间状态划分标准Tab.2 Division of risk factors in multi-interval safety states for root nodes

传统的事故树模型以及二态B-Ns系统仅将节点状态简化为“正常”和“失效”两种，但在实际中中间节点与叶节点往往存在着中间状态，即其节点状态虽然不属于正常状态范畴，但是并未导致大的灾害形成.以冰碛湖溃坝机制中占比最高的冰崩涌浪漫顶为例，漫顶涌浪往往需要达到一定的水头高度才能使得堰塞坝溃口处粗颗粒泥沙连续起动形成溃决连锁反应[18].因此，根据事故严重程度，本文将中间节点y1的状态分为低、中、高危险3类，状态值分别用0、1、2表示；其余溃决机制（y2～y5）由于发生几率较小，为降低贝叶斯网络复杂度，仍将其简化为两个状态（以0、1表示）.对于叶节点T，从堰塞湖溃决物理过程角度，其溃决方式可以分为全溃与部分溃，其中全溃形成的相对洪峰流量最大，危险性也最高[19].根据溃决危险性严重程度，将叶节点T分为低、中、高3种危险状态.

1.2.2 节点概率分布的获取

构建B-Ns结构后，需要确定根节点各状态先验概率分布及其他节点（中间节点和叶节点）条件概率分布.节点概率的获取方法主要有资料统计法和专家知识法两种.B-Ns方法优势之一在于能够有效地融合多源信息，既可以利用现有资料统计获取客观概率值，也能够利用专家知识对节点概率进行补充.因此，本文在确定B-Ns节点概率分布时结合资料统计法与专家知识法互相验证.

根节点先验概率可依据资料统计法获得.根据所研究区域的冰川、冰碛湖、遥感影像与数字高程模型（DEM）地形对根节点因子落入各状态区间频数进行统计，从而确定根节点各状态的先验概率.本文所使用统计资料中，冰川、冰湖数据来自国际山地综合发展中心（ICIMOD，2015年）与中科院寒旱所（2014年）提供的第二次冰川编目；冰碛坝等几何与结构参数源自对遥感影像（2014年—2019年）与DEM地形（2009年）的解译；气候数据来自中国气象数据网提供的研究区气象站近30年气象资料.冰碛湖溃决案例为本文的附加材料表S1.

通过以上方法获取根节点先验概率之后，仍然需要确定其他节点的条件概率表（CPT）.节点处CPT可利用模糊数学法将专家的定性描述转化为定量的发生概率[8].针对专家知识法获得的CPT存在主观性的问题，可以将获取的条件概率与实际案例互相验证，以降低主观性，操作程式见1.3节案例.

1.2.3 冰碛湖溃决概率计算

根据B-Ns的运算规则，可以由根节点的先验概率分布出发正向推算叶节点事件的发生概率.

在多态B-Ns中，根节点为xi；中间节点为yj；叶节点为T.、和Tq分别代表根节点、中间节点、叶节点的风险状态.其中：ai=0,1,···,ui−1；bj=0,1,···,vj−1；q=0,1,···,r−1，ui、vj、r分别为根节点、中间节点、叶节点的状态数.则当根节点各状态的概率分别为P,P,···,P时，叶节点T处于风险状态Tq的概率计算式为

式中：π(yj)为中间节点yj的父节点集合；π(T)为叶节点T的父节点集合.

依据待评估冰碛湖的证据信息，利用B-Ns的正向推理可以计算叶节点T各状态的发生概率.确定各风险状态的隶属度函数为

式中：p为叶节点各状态发生概率值.

根据最大隶属度原则确定叶节点隶属的状态.评估程式如图3所示.

图3 基于B-Ns的冰湖溃决危险性评估程式Fig.3 Risk assessment procedure for moraine-dammed lake outburst based on B-Ns method

1.3 B-Ns模型的检验

评价模型的精度直接关系到评价结果的可靠性，因此使用前有必要对该模型进行检验.本节利用1981年7月11日波曲河流域次仁玛错冰碛湖案例对B-Ns模型的精度与作业流程进行验证和说明.

根据吕儒仁等[4]的调查，溃决当年喜马拉雅南坡的气候背景为暖干气候，降水因为量小而作用很小，突然升温导致冰崩，叠加管涌效应导致冰碛湖溃决.该冰碛湖的海拔高度为4660 m，溃决前夕冰舌到达海拔4700 m位置，距离冰碛湖非常近.冰川前缘地形为陡崖状，且在陡崖处呈不连续状态，裂纹发育.溃决之前，冰碛湖处于满溢状态，坝顶槽道发育，已有湖水从坝体中连续渗出.溃决危险性评价流程如下：

1）确定节点概率分布

根节点各风险状态区间的先验概率可依据资料统计法获得，如表3所示.

表3 冰碛湖溃决根节点的先验概率Tab.3 Prior probabilities of root nodes of moraine-dammed lake outburst

依据各风险因素间的因果关系，结合案例资料与专家知识可构造中间节点y1～y5及叶节点T各状态区间的条件概率分布.本文仅选取中间节点y1部分条件概率分布作为示例，如表4所示.表中，第一行数字的意义为：当根节点x1～x5状态量均为0时，中间节点y1状态量为0的概率等于1.0，为1、2的概率均为0.其余各行数字含义可类似推出.完整表格详见附加材料表S2.

将表3中根节点x1～x5的先验概率与表4中y1的条件概率代入式（1），可求得节点y1各风险状态的概率为同理可得，P(y1=1)=0.2655，P(y1=2)=0.0526.

表4 中间节点y1的条件概率Tab.4 Conditional probabilities of node y1

上述计算结果表明，在没有相关证据信息输入B-Ns模型条件下，中间节点y1的状态为低危险的概率很大，中、高危险概率均较小，其中为高危险的概率仅为0.0526.据统计，波曲河流域冰湖数量为121个，参照附加材料表S1，该流域发生严重冰崩涌浪漫顶溃决灾害的案例为4次，概率值计算为0.0331，与本文计算的P(y1=2)=0.0526 较为接近.可初步验证B-Ns模型结构及参数.根据类似方法，可构造其余节点的条件概率.

2）依据证据信息对待评估冰碛湖溃决危险性进行评价.

根据溃决前次仁玛错冰碛湖状态，更新贝叶斯网络中各根节点状态量，结合中间节点的条件概率进行推理计算，可得P(T=0)=0.0670，P(T=1)=0.2381，P(T=2)=0.6949.根据最大隶属度原则，该冰碛湖的溃决危险性评价结果为高危险，与实际情况相一致，从而对本文B-Ns模型准确性进行了进一步验证.

2 冰碛湖溃决危险区沿河线铁路方案比选评价指标体系

2.1 各类风险区线路长度计算方法

各类风险区线路长度应是溃决危险区铁路选线评价时最重要的指标，其计算需综合考虑冰碛湖的危险性与线路工程所处河谷段地形.

青藏高原地区，冰碛堰塞湖均分布于主河两侧的支沟内，而线路工程多沿主河布设.冰碛湖一旦溃决，溃决洪水先沿着支沟演进，进入主河后若形成超标洪水才会对线路工程造成灾害.能否形成超标洪水的必要条件是线路工程位于溃决洪水的淹没范围内，而每个冰碛湖的淹没范围与其体积（V）有关.现借用Washakh等[20]提出的冰碛湖溃决洪水淹没长度（L）计算式 （L=0.52V+36.13），若定义溃决洪水到线路工程的路径长度为Ld，则筛选出需进行风险评估的冰碛湖筛选判据为

能够形成超标洪水的另一必要条件是沿河线所在河谷地形，显然同样流量的条件下越是地形狭窄的河道形成超标洪水的情况越严重.

综上，各类风险区线路方案长度指标计算不仅需要考虑冰碛湖溃决危险性，还需要考虑冰碛湖溃决影响范围、沿河线河谷形态等因素，具体计算流程如图4所示.

图4 各类风险区线路长度计算流程Fig.4 Calculation process of the railway line length in different risk areas

2.2 铁路选线方案冰碛湖溃决危险性评价指标体系

铁路方案冰碛湖溃决危险性评价不仅需要考虑流域内冰碛湖的自然状态（数量、总面积、历史溃决次数），也要考虑危险性冰碛湖对铁路线路的影响.结合这两类考虑因素，提出铁路方案冰碛湖溃决危险性评价指标体系，如表5所示.以中尼铁路为例进行作业程式说明.

表5 铁路方案冰碛湖溃决危险性评价指标Tab.5 Evaluation indexes of railway schemes considering moraine-dammed lake outburst risk

3 中尼铁路局部走向方案评价案例分析

中尼铁路是拟建穿越喜马拉雅山脉的铁路，线路从日喀则出发先在高原面上行走，而后下高原面至加德满都，这时面临经吉隆廊道与经樟木廊道这两个局部走向方案的比选问题（图5）.吉隆、樟木廊道均位于喜马拉雅中部希夏邦马峰冰川群分布区，冰碛湖数量众多，冰碛湖溃决洪水已成为对线路方案起控制作用的重大灾害类型.

图5 吉隆、波曲流域危险冰碛湖与线路方案示意Fig.5 Railway schemes and distribution of dangerous moraine-dammed lakes in the Gyirong and Poiqu river

1）研究区自然条件

吉隆流域位于喜马拉雅中部.东西走向的喜马拉雅山脉将吉隆流域分割成南北两部分，北部海拔较高、地势平缓；南部海拔低，但是山体陡峻，高差较大.远古时期吉隆地区分布有吉隆、沃玛等湖泊，由于青藏高原的不均匀抬升，导致吉隆藏布溯源侵蚀切穿原有湖泊，湖水下泄，这样原古湖盆与现代河流形成了吉隆藏布上游段宽窄相间的河谷地貌形态.此外，由于海洋性冰川塑造地貌效应，中游段也形成了一些宽谷地形.宽谷和峡谷段分布如图5所示.

樟木廊道位于波曲流域，波曲河发源于西藏聂拉木县波绒乡，经聂拉木县城到樟木镇.流域内地形高差较大，最高峰为海拔超过8 km的希夏邦马峰.波曲的河谷地貌以聂拉木县城为分界线，其北基本位于高原面上，主要以宽谷地貌为主，其南属于溯源侵蚀河段，以峡谷地貌为主（图5）.

2）研究区冰碛湖溃决危险性评价

首先利用式（2）筛选出需要进行危险性评价的冰碛湖，然后利用本文建立的B-Ns评价模型，对吉隆、樟木廊道流域内冰碛湖溃决危险性进行评价（如图5）.可知：吉隆流域危险性较大的冰湖总数为7个（其中高危险冰碛湖4个，中危险冰碛湖3个）；樟木流域危险性较大的冰碛湖总数为13个（其中高危险冰碛湖6个，中危险冰碛湖7个）.

3）中尼铁路樟木、吉隆局部走向方案比选

中尼铁路跨喜马拉雅段的樟木、吉隆局部走向方案如图5所示.

依据本文的铁路方案冰碛湖溃决危险性评价指标对两方案进行比选，各项指标值如表6所示.由表可知：虽然樟木方案在冰碛湖分布区的线路长度小于吉隆方案，从选线角度属于线路短直方案，但是樟木方案的高、次高、中风险区段线路长度都远超过吉隆方案.因此，从冰碛湖溃决风险角度，吉隆方案优于樟木方案.

表6 吉隆、樟木线路方案比选表Tab.6 Comparison of the railway schemes in the Gyirong and Poiqu river basins

4 结 论

本文以冰碛湖溃决危险区铁路方案风险评价为目标，建立了基于多态B-Ns的冰碛湖溃决危险性分析方法及冰碛湖溃决风险区新建铁路选线方案风险评价指标体系.主要结论有：

1）通过对青藏高原南缘74个冰碛湖溃决灾害实例统计分析，在此基础上针对冰碛湖溃决影响因素众多、关系复杂且具有高度不确定性的特点，建立了多态B-Ns冰碛湖溃决灾害发生的概率预测模型；以冰碛湖分布区沿河线路工程为承灾体，依据冰碛湖溃决危险性评价结果与沿河线河谷地貌形态特征，提出了铁路选线方案冰碛湖溃决危险性评价指标体系，从而为冰碛湖溃决危险区铁路方案风险评价提供了具普适性的技术评价体系.

2）本文提出的铁路方案冰碛湖溃决危险性评价指标，仅是从冰碛湖角度对现行评价体系进行增补.在实际工程中，仍需在确定常规方案技术经济指标的基础上，结合冰碛湖溃决危险性评价指标进行综合分析，评选最合理方案.

3）铁路选线一般需经过从线路走向的拟定到空间定线等从粗到细逐步优化的阶段，直至确定合理的线路位置.本文提出的铁路选线方案冰碛湖溃决危险性评价指标体系主要用于铁路行经地区局部走向方案的选择；此外针对空间定线作业阶段，提出在冰碛湖溃决高危险区的峡谷段，宜采取抬高程、尽量不跨河、酌情预留采取提高限坡措施条件等定线要点；二者共同构成冰碛堰塞湖溃决减灾选线风险调控技术体系.

备注：附加材料在中国知网本文的详情页中获取.

致谢：本文研究得到了中国科学院水利部成都山地灾害与环境研究所苏立君研究员，中铁第一勘察设计院集团有限公司苗晓岐教授级高工，梁栋、孙先锋高级工程师等热情的帮助和指教，在此一并深致谢意；感谢中国科学院国际伙伴计划“一带一路”科技合作项目“中尼交通廊道灾害风险及其应对策略研究”（131551KYSB20180042）、中铁二院工程集团有限责任公司科研项目（KYY2020054（20-22）对本文研究的资助.