考虑人体与水流相互作用的溃坝洪水生命损失评估模型

彭 铭，马福军，沈丹祎，蔡一坚，石振明，周家文，刘西军

(1.同济大学土木工程学院，上海 200092；2.同济大学岩土及地下工程教育部重点实验室，上海 200092；3.浙江大学建筑工程学院，浙江杭州 310058；4.中国电建集团华东勘测设计研究院有限公司，浙江杭州 310000；5.四川大学山区河流保护与治理全国重点实验室，四川成都 610065)

堰塞坝是由地震、降雨、火山喷发等诱发斜坡失稳堵塞河道而形成的天然坝体。据统计，51%的堰塞坝会在形成后的7 d内发生溃坝[1]。堰塞坝溃决通常具有突发、隐秘和难以预控的特点，极易对下游群众生命安全构成威胁。如：1933年，四川叠溪堰塞坝在形成45 d后突然溃决，造成下游2 500多人死亡；2018年，金沙江白格堰塞坝在形成10.6 d后发生溃决，溃决洪水威胁下游500多千米的丽江市，迫使25 000人紧急撤离[2]。因此，定量评估溃坝洪水造成的生命损失，进而指导洪水影响范围内的人员疏散具有重要减灾价值。

溃坝洪水造成人类生命损失的影响因素包括直接因素和间接因素。其中，水流流速和水深是最显著的直接影响因素，对暴露在洪水中人的稳定性至关重要。现有的溃坝生命损失模型包括经验模型、物理模型和混合模型。经验模型主要通过历史数据拟合生命损失率和关键变量间的关系[3–5]，易于应用，但不能揭示实际物理作用机制和变量间的相关关系。物理模型侧重于研究人在水中的行为和稳定性[6–9]，但不能给出生命损失和稳定性间的关系。以上两种模型大多是基于假设几个互相独立的变量建立的定性或半定量研究，没有考虑到变量本身存在的不确定性及变量之间的相关关系。混合模型将经验模型中历史数据、专家经验和物理模型揭示的内在规律统一起来，扬长避短[10–12]。

近年来，很多基于贝叶斯网络建立的混合模型被广泛应用于各行各业的决策、推理和风险评估。贝叶斯网络由节点和弧链接及其（条件）概率组成，可以通过逻辑推理解决不确定性问题。Zhang[13]和Xu[14]等采用贝叶斯网络诊断了堤坝灾害。Peng和Zhang等[11–12]基于贝叶斯网络提出HURAM 1.0模型分析了溃坝生命损失，该模型考虑了水力学参数、地形参数和人员疏散情况等15个变量，采用条件概率考虑变量间的不确定性及相关性。但HURAM 1.0模型没有考虑到人体和洪水间的物理相互作用过程，忽略了水动力因素对人体稳定性的影响，在水流速较大时的预测结果与实际情况偏差较大。

本文基于HURAM 1.0模型，嵌入考虑人体稳定性的物理模型，用基于物理和统计的新模块替代原本假设生命损失是水深函数的经验模型，考虑了水流流速对生命损失的影响，从而得到HURAM 2.0模型；将所提出的模型和HURAM 1.0模型对比并开展参数敏感性分析；同时，将所提出的HURAM 2.0模型应用于唐家山堰塞坝中评估溃坝洪水对北川县造成的生命损失风险。

1 生命损失风险评估模型

1.1 HURAM 2.0模型

HURAM 1.0模型是Peng和Zhang[11]基于贝叶斯网络提出的生命损失风险评估模型，实现了基于物理机理和概率统计的综合定量分析。图1红框部分为建立HURAM 1.0模型的逻辑结构示意图。

图1 HURAM 2.0模型逻辑结构示意图Fig. 1 Logic structure of the HURAM2.0 model

HURAM 1.0模型的研究对象是所有受洪水影响区域内的人口，包括风险人口和暴露人口。其中，风险人口所处的区域水深大于0，暴露人口是指没有成功从洪水淹没地区撤离的人员[1]。同时，考虑建筑物的损坏和淹没程度，模型将洪水的破坏性分为安全、低破坏性、中破坏性和高破坏性4个等级，并考虑人员在不同洪水强度下的初始分布及预警作用下的人员重新分布，假设房屋建筑均匀分布在所在区域上。安全等级下，建筑物保持结构完整且楼层没有淹没，死亡率假设为0；高破坏性等级下，建筑物结构破坏，且楼层完全被淹没，通过大量数据统计得到死亡率为0.91；低、中破坏性等级下，假设生命损失率FD是自变量为水深d的对数正态分布函数：

通常采用大量案例数据得出中、低破坏性下的生命损失率。其中：低破坏性时， µN=3.376， σN=1.188；中破坏性时， µN=1.649， σN=0.562。

HURAM 1.0模型低破坏性等级下流速较小，主要为淹没风险，故只考虑生命损失率是水深的函数是合理的。但在中破坏性等级下，水流流速相对较大，快速冲击的水流会使人体自身稳定性大幅下降。如在暴雨产生的急性山洪事件中，快速流动的河水经常会造成伤亡事件，这其中的水动力因素作用不可忽视。另外，人在水中有一定的求生能力，应予以考量。

基于此，HURAM 2.0模型在原模型基础上，增加了一条从水流速到生命损失的新连接，引入人体稳定性物理模型，对处于中破坏性等级下的人体基于水深和水流流速进行稳定性判定，考虑了中洪水强度下流速对人体稳定性的作用，进而得到其对生命损失的影响规律（图1）。假定当事人失稳跌倒，即进入溺水状态；溺水状态下，人无法保持稳定站立，需要通过游泳求生，不会游泳的人处于高风险，会游泳的人进入求生状态；求生状态下，人体在水中的游泳能力和水流流速有关，当超过游泳速度的极限后，人体将失去控制，体力不支而处于高风险状态。

1.1.1人体在水中稳定性物理模型

假设人体站在水中仅受到水流的拖曳力D、流水的升力L、人体自身的重力G和水对人体的浮力B（图2（a）），则有：

式（2）～（5）中， ρ 和 ρP分别为水和人体的密度，CD和CL分别为阻力系数和升力系数，v为水流速度，Hp为人体的身高（从肩膀到脚底），w为人体正面的宽度，b为水流速度方向上人体的投影长度。

1）失稳模式

由于个体差异和洪水环境的复杂性，人体的失稳模式多种多样，现有研究大都做了简化处理[6–8]，只考虑滑动和倾倒两种最常见的失稳模式（图2（b）），具体分析如下：

图2 站立在水中的人受力示意图和失稳模式Fig. 2 Forces acting on a submerged human body and instability mode

a.滑动失稳

当脚底与地面之间的摩擦力不足以抵抗正面的动水压力时，人体滑动失稳，其极限平衡状态为：

进一步处理可得：

根据v–d平面（图2）上数据的分布，可以发现人存在最大防洪极限条件，这些极限条件对应于流速较小时相对较大的水深，以及流速较大时相对较小的水深。此时，选取合适的α和 β，利用冲击参数W可以获得水与人体之间物理作用机制和人体稳定性破坏过程。v–d平面上等W值线上的点对特定类别的人体破坏过程的破坏潜力是等价的，通过经验数据拟合得到W=1.00、W=0.35等值线（分别对应成人和儿童），较好地包裹了数据点的上限。以上表明，超过此界限，人体极易失稳，此时YW=1.25m，α=2.0，β=4.0。

3）人体在水中站立稳定性量化评估

当淹没水深处于0～1.5m区间内，基于人所处环境下的水流流速和水深进行失稳判定。通常，水流流速越大，稳定性越差，成人和儿童由于身体素质不同，具有不同的脆弱性（成人稳定性∶儿童稳定性=4∶1）。以成人为例，当W≤1时，人体可保持站立稳定，认为安全；当W> 1时，人体失稳进入溺水状态。每种情况下采用蒙特卡洛模拟100万次，得出人体在流动水中的稳定性量化结果（表1）。

表1 人在水中的站立稳定性Tab.1 Stability for a flooded human body

1.1.2溺水状态下人的生命损失分析

溺水状态的人需要通过游泳逃生，人在水中的自控能力与水流流速有关，当水流速超过游泳速度的极限（以世界冠军游泳速度为基准，一般是普通人的2～3倍），人体将失去控制处于高风险状态。

溺水状态下人本身有一定的求生欲望和能力，但受客观环境影响，其求生能力有限。随水流速的增大，死亡率会急剧上升，在水流超过2m/s时，随着水流速继续增大，死亡率会不断地逼近上限100%，出现收敛趋势，即流速和死亡率近似呈指数关系（图3）。

1.1.3考虑人体稳定性和求生能力的生命损失量化

图4为人体稳定性和求生能力量化逻辑结构。

图4 人体稳定性和求生能力量化逻辑结构Fig.4 Quantification logical structure of human stability and survival ability

由图4可知，组合各种不同水流流速和水深条件，各采用蒙特卡洛实验重复100万次，可得到在不同的洪水状况组合情况下人体的生命损失量化值，见表2。

表2 基于人体站立稳定性物理模型死亡率量化结果Tab.2 Mortality quantification results based on a physical model of human stability

将表2的条件概率嵌入到修改的贝叶斯网络中，使用Hugin Lite计算可得到新模型各节点的先验概率分布。

1.2 模型对比分析

使用Hugin Lite计算得到的先验概率是基于统计数据和已发表文献中的成果，因此，结果仅代表溃坝洪水造成的全球总的平均生命损失风险。

HURAM 2.0模型计算出的先验生命损失概率为7.58%，比HURAM 1.0模型计算出的生命损失（4.54%）风险大。其原因主要是HURAM 1.0模型忽略了水动力因素对人体稳定性的不利影响，低估了特定情况下的生命损失风险。改进模型在保持其他节点概率分布都不变的条件下，进一步考虑了水深和水流流速对生命损失的影响。

在具体的溃坝案例中应用时，需要输入具体的水力参数、人员分布信息、建筑物特征等信息，基于贝叶斯定理更新先验概率得到后验分布。

图5为HURAM 2.0模型和HURAM 1.0模型预测的死亡率。

图5中所使用的历史数据信息经过校核和正则化，剔除了一些不符合常理的异常值；斜率为1的对角线代表预测值等于观测值。由图5可知：死亡率较高（大于5%）时，HURAM 2.0模型预测性能明显好于HURAM 1.0模型；死亡率较低时，两模型预测性能相当。整体规律为死亡率低于1.44%时，HURAM-2.0模型预测值比HURAM 1.0略小；当死亡率值较大时，HURAM 2.0模型预测值比HURAM 1.0模型略大。

图5 HURAM 2.0、HURAM 1.0模型的预测结果Fig. 5 Predicted fatality rates using HURAM1.0 and HURAM2.0

1.3 HURAM 2.0模型参数敏感性分析

图6为14个指标对先验概率的影响程度，对比HURAM 2.0和HURAM 1.0模型，主要差异有：1）变化不大的，如洪水严重程度、水深度等影响因素。因为HURAM 2.0模型的贝叶斯网络仅修正了中风险程度下的生命损失函数，不影响安全和低死亡率（FD= 0），以及高风险高死亡率（FD=0.9077）下的生命损失，故洪水严重程度影响范围区间变化不大，且水深始终是导致生命损失的主要因素。2）变化范围扩大的，如疏散率、洪水到达时间段、在建筑物中庇护情况、预警时间、洪水上升速率、水流流速等，由于考虑了流水动力因素对人体稳定性的不利作用，修改后的贝叶斯网络对水流流速的变化更敏感，并导致生命损失风险整体上升，上限提高，变化范围扩大。3）变量的区间向上发生位移的，如居民区住宅类型、淹没区距坝址处距离，这两个相对不重要的影响因素因生命损失整体风险上升而产生取值的偏移。

表3为生命损失的参数敏感性分析及对比。

表3中，Imin和Imax[11]分别反映各变量对生命损失值的正面（从最小值开始增大）和负面（从最大值开始减小）影响的能力大小，是一个梯度变量，值越大代表影响越大。由表3可知：1）HURAM 2.0模型中水流流速的影响权重（图6水流速对应的条件概率取值范围扩大）和敏感性均有所提高，水深敏感性相对下降，原因可能是考虑了深度相对较浅、流速相对较大的人体失稳情况（HURAM 1.0模型未考虑），此时控制生命损失的主导因素是较大的水流流速；2）居民区住宅层数、在建筑物中庇护情况和溃坝时长的影响程度更大，这是因为这3个变量均能直接影响到人员的疏散和庇护情况，既没有撤离又无成功庇护的暴露人口将直接在洪水中受到水流的冲击，从而变得更危险。

图6 先验的HURAM 1.0和HURAM 2.0中的14个影响指标的分析Fig. 6 Influences of 14 factors on the loss of life in both prior HURAM1.0 and HURAM2.0

表3 生命损失的参数敏感性分析及对比Tab.3 Parameter sensitivity analysis and comparison of life loss

综上所述，HURAM 2.0相比HURAM 1.0模型，很多变量的敏感性都有所上升，这意味着贝叶斯反演时更具优势。此外，由于改进后的模型引入了基于实验数据的物理模型，构建更加精细，更加详尽地考察了水动力作用因素，增强了模型的解释性，故改进模型的贝叶斯网络反演更具有说服力。

2 案例分析

2.1 唐家山堰塞坝介绍

2008年5月12日，中国四川发生了里氏8.0级地震，导致唐家山山体滑坡堵塞湔江上游，形成了一个巨型堰塞坝，如图7所示。堰塞坝坝高82m，坝宽（顺河向）802m，坝长（横河向）611m，坝体体积2.04×107m3，库容3.16×108m3，回水长度达20 km，集水面积近3 550 km2；坝体距离北川镇上游3.5 km。唐家山堰塞坝库容量巨大，汇水面积大，水位上涨快，坝体结构松散，溃坝风险高，对下游120万人的生命财产安全构成巨大威胁。

图7 唐家山堰塞坝卫星图像Fig. 7 Tangjiashan landslide dam satellite image

唐家山堰塞坝坝体为3层结构，最上层是厚度为5～15m的碎石土；中间是10～15m强风化碎裂岩；最下层为50～80m厚的弱风化碎裂岩石[15–17]；坝体抗侵蚀性从顶层到底层由弱变强[18]。2008年6月1日，开掘长475m、宽25m、深12m的泄流槽，开挖土方量1.36×105m3。泄流槽将原始的坝顶高程从752.2m降至740.4m，净容量降至2.47×108m3[17]。

2008年6月11日上午6：00，唐家山堰塞坝开始溃决；12：30达到峰值流量6500m3/s，坝顶高程下降至735.8m；14：00湖水位高程下降至714.1m，相应库容降至8.6×107m3。最终溃口深度42 m，顶宽145～235 m，底宽80～100m。

2.2 溃坝洪水分析

2.2.1工况设置

根据对堰塞坝信息的掌握程度由浅到深和工程减灾措施的从无到有，即堰塞坝灾害的发展阶段设置两种工况研究溃坝洪水生命损失。

工况1：滑坡形成的堰塞坝结构松散，本身极不稳定，很快会发生溃决，在最初几个小时很难获取坝体结构和材料参数，需要对堰塞坝的安全状况做快速评估[19–21]。分析中，按照最坏情况考虑，即设定工况1为高侵蚀坝体材料下发生漫顶溃决的工况[21]。

工况2：通过现场勘测获得坝体基本信息，包括坝体材料、结构参数等，得知坝体由3层不同风化程度的碎裂岩构成，堰塞坝的竖向非均质性将影响到坝体的溃决过程[22]。为降低风险，现场通过开挖泄流槽降低了坝顶高程，导致坝体提前发生漫顶溃决。设定工况2为现场勘测和开挖泄流槽工况。

工况1坝顶的高程为752 m，积蓄的最大库容3.16×108m3；工况2开挖泄流槽后，坝体高程降至742 m，湖水最大库容下降至2.47×108m3。

2.2.2洪水演进模拟

首先，获取河道3维地形信息，主要包括：通过Google Earth获取研究区域的DEM高程信息图片和卫星影像图；采用地图分析软件G lobalM apper获取典型断面的3维数据。该方法可以在短时间内获取全球范围内任何堰塞湖所在流域河道的3维地形信息。在此基础上，开展溃坝模拟和洪水演进分析：

1）工况1中，对堰塞坝的信息掌握不充分，只有通过航空照片估算坝体几何参数，采用Peng等[1]根据52例堰塞坝破坏信息数据库提出的多参数经验模型估算坝体溃决的参数。使用基于稳定流能量守恒和非稳态流的质量和动量守恒开发的1维水力演进计算模型HEC–RAS进行洪水演进计算，模拟得到上、下游（北川）水位高程和溃坝流量如图8所示。溃坝持续时间为8.9 h，坝址处的峰值流量达到31 749m3/s；流经北川时峰值流量为31 670 m3/s；最大淹没水深20 m，流速达到1.9m/s。

图8 工况1和工况2上下游水位高程变化和北川流量过程曲线Fig. 8 Upstream and downstream water level elevation and Beichuan flow process curve in Case 1 and 2

2）工况2中，通过现场勘测得到坝体的基本参数（坝体材料、结构参数），在开挖泄流槽后，使用物理模型DABA[23]模拟坝体溃决，该模型基于侵蚀理论和浅水流模型模拟水土之间的相互作用。通过侵蚀系数和临界侵蚀剪应力描述土壤侵蚀性。通过求解河流的连续性方程及宽顶堰的溃口水流模拟流体动力学过程。输入大坝几何信息、土体参数及上游的日平均入流量，即可输出每个时步的溃决参数。将坝址处溃决流量输入HEC–RAS模型，得到下游各位置流量过程线。北川的峰值流量为6 567m3/s，相比工况1大大减小；溃坝持续时间为14 h；北川最大淹没水深为6.1 m，平均水流速为1.1m/s。

2.3 唐家山堰塞坝溃决生命损失风险评估

根据疏散距离和水深将断面信息离散成输入贝叶斯网络的i个区段。假设人员均匀分布在居民区内，那么在不同工况Sj（本文工况1、2）下，人类生命损失风险可计算为：

式中，PPARi为第i个区段上的风险人口，PiL|S j为工况S j下区段i上的死亡率。

以工况2为例，将典型断面分割成6个区段（图9），各区段输入改进的贝叶斯网络的变量取值见表4，计算结果见表5。同样地，计算出工况1、2的生命损失，结果见表6。

图9 工况2风险人群划分Fig. 9 Population at risk in Beichuan in Case 2

表4 工况2输入改进的贝叶斯网络的参数Tab.4 Inputs for improved Bayesian network in Case 2

表5 HURAM 1.0和HURAM 2.0计算得到的北川地区各区间生命损失计算结果对比（工况2）Tab.5 Comparison of life loss in Beichuan area by HURAM1.0 and HURAM2.0 models(Case2)

表6 溃坝洪水流量及其造成的北川人员生命损失Tab.6 Dam-break floods and caused loss of life in Beichuan Town

2.3.1工况1分析

工况1情况下，溃坝时，改进模型计算出的北川人类生命风险分布如图10所示，IR为个体生命损失风险，下同。

北川几乎全部区域被洪水所淹，淹没最深处达20m，淹没区域风险人口达30 000人。即使在发出疏散警报（贝叶斯网络预测）的情况下，生命损失风险也达到4 000左右，属于极大的灾害事故。这主要是由于在工况1中，假设坝体具高侵蚀性，溃决时间很短，最终溃口深度大，漫顶溃坝时积蓄的库容极大，短时间内北川遭遇严重洪水，所有未撤离的人群都将遭遇致命的洪水侵害，风险极大。

图11为HURAM 1.0和HURAM 2.0模型计算的生命损失率。工况1水深20m，水流速为1.94m/s时，属于中破坏性洪水的范畴。在HURAM 1.0模型中（图11虚线），暴露人口在中等洪水下，仅采用有限的几个历史数据拟合，当流速超过9 m/s时，死亡率即高达0.892 1；对不同的水流速状况未做出进一步区分，即不管水流速有多大，只要水深确定，其预测的结果都将保持不变，与事实不符。HURAM 2.0模型采用人体稳定性物理模型，暴露人口在中等洪水下（图11实线和短划线），当流速超过9.00 m/s（大于1.5m的参考深度），考虑到溺水人群在1.00～2.00m/s的流速中具有一定的求生能力，死亡率仅0.819 2，小于HURAM 1.0模型的预测结果。同时，当水流速降低时，死亡率会降低；当水流速超过人体游泳能力的极限2m/s时，死亡率急剧上升，改进后的模型能够区分不同水流条件下的生命损失计算结果。

图10 工况1北川人类生命风险分布Fig. 10 Human life risk distribution for Beichuan in Case 1

图1 1 HURAM 1.0和HURAM 2.0模型计算的生命损失率Fig. 11 Life of loss calculated by HURAM1.0 and HURAM2.0 models

2.3.2工况2分析

工况2中，北川人类生命损失风险的强度和分布范围都大幅度降低，如图12所示。

开挖泄流槽后，坝顶高程降低，峰值流量大大减小，洪水强度大幅度下降，减轻了堰塞坝溃决对下游的影响，最大淹没深度为6.1m，北川只有一部分被水淹没，风险人口为18 300，生命损失风险下降至20人左右。此时，多数未撤离人员都可以在3楼寻求庇护，只有少部分来不及采取庇护措施的人遭受一定的洪水风险。HURAM 2.0模型预测的人类生命损失风险为21.96，比HURAM 1.0模型预测结果（15.31）略大，主要是考虑了水流流速对人体稳定性的不利影响。

图12 工况2北川人类生命风险分布Fig. 12 Human life risk distribution for Beichuan in Case 2

由以上可知，及时到现场勘探查明堰塞坝的资料，尽早采取工程减灾措施，可以在很大程度上降低堰塞坝灾害造成生命损失的风险；且在工况1中，HURAM 2.0能区分不同水流速度下的生命损失。

3 模型局限性

HURAM 2.0模型是一个宏观模型，在统计学平均意义上预测洪水水流对生命损失的影响，受限于当前的基础和水平，做了一些必要的假设和简化，可能会降低模型的精度和适用性。

1）溃坝洪水中含有大量泥沙，一般距离坝址处越近，洪水泥沙含量的不利影响就越大；距离坝址越远，泥沙的影响越小。本文暂未考虑水沙混合流对人体稳定和求生能力的影响。在后续研究中将进一步考虑水沙混合流与人体之间的相互作用。

2）只考虑了水流直接作用导致的滑动和倾倒失稳模式，区分了儿童和成人两类群体的差异，在区域统计学的平均意义上预测效果要好于对局部的单独个体的预测效果。后续研究将进一步关注和考虑复杂环境（如被树枝和大石块绊倒、被水沙流体大石块撞击失稳）对特定类别个体的影响。

3）HURAM 2.0模型在计算溃坝生命损失方面具有较大优势，但由于所建立的贝叶斯网络只能接受离散数据信息条件的输入，离散化的过程会不可避免地降低数据的分辨率和敏感性。未来将进一步考虑建立连续变量的贝叶斯网络。

4 结 论

1）在溃坝风险评估模型HURAM 1.0的基础上，进一步考虑水流作用下人体稳定性，提出HURAM 2.0模型。该模型考虑了多因素相互作用及不确定性的影响，实现了在较强洪水强度条件下的生命损失更准确的预测，增强了模型的解释性。

2）在生命损失易损性函数中引入水流速度参数，更精确地刻画了人体在水流中的稳定性和求生能力。相比HURAM 1.0模型仅适合低、高洪水强度生命损失预测，HURAM 2.0模型可以实现低、中、高3类洪水强度下生命损失的准确预测，并且能在中等破坏性洪水下区分出水流流速的影响，更符合实际情况。

3）将本文模型应用到唐家山堰塞坝溃坝洪水风险评估中，结果表明：在未查明堰塞坝基本资料也未开挖泄流槽的情况下，坝体溃决导致的生命损失风险较大；在获得坝体3层结构并开挖泄流槽降低坝顶高程后，生命损失风险大幅度降低。这说明及早查明堰塞坝的情况并采取相应的工程措施（如开挖泄流槽）能在很大程度上降低堰塞坝灾害的风险。