基于脆弱性曲线的寿溪河流域村镇建筑洪灾风险评估

摘要

以寿溪河流域的村镇建筑为研究对象，采用MIKE FLOOD耦合水动力模型计算洪水淹没特征，结合结构静力学原理与洪水冲击荷载经验公式，建立村镇建筑的洪灾脆弱性曲线，开展洪灾风险评估。结果表明，寿溪河流域洪水淹没范围呈狭长的条状，洪水流速较大。低重现期下房屋以中低风险为主，高重现期下高风险与极高风险房屋主要分布在洪水发生点、洪水淹没区以及河坝村附近。建立的脆弱性曲线考虑了洪水流速的影响，能较好地刻画村镇建筑的结构损伤程度。

关键词

洪灾风险评估; 村镇建筑; 脆弱性曲线; 寿溪河流域; 情景模拟

洪水是当今世界上发生频率最高、最严重的自然灾害之一^［1］。根据水利部颁布的《中国水旱灾害防御公报》^［2］，1999年―2020年，我国平均每年因洪涝灾害死亡1 110人、倒塌房屋82.63万间，直接经济损失达1 688.28亿元。开展洪水风险管理是应对洪水威胁的有效方法^［3］，其中最关键的措施是进行洪灾风险评估。目前国内外的洪灾风险评估主要针对流域和区域开展，常用方法包括历史灾情法^［4-5］、指标体系法^［6-7］和情景模拟法^［8-9］。在洪涝灾害造成的各类损失中，村镇建筑的损坏对生产和生活的影响最大，也最为严重^［10］。针对房屋建筑的洪灾风险评估在我国并不多见，而在国外有一定研究基础，常用方法包括指标体系法^［11-12］和脆弱性曲线法^［13-14］。指标体系法能够很好地反映区域的致灾特性与建筑的承灾特性，但由于涉及分项指标众多，存在指标选取标准不一、获取难度较大等问题，该方法目前仅在小范围区域内得到应用^［15］。洪灾脆弱性曲线是反映承灾体灾损程度与洪水要素函数关系的曲线（面），具有评估效率高、适用范围广的特点^［16］。国外脆弱性曲线的建立方法主要包括历史灾情法^［17-18］、模型试验法^［19-20］和现场观测法^［21］。其中，最常用的是历史灾情法，即依托历史灾损数据库建立脆弱性曲线（库中记载历史洪灾的淹没范围、淹没水深分布，各建筑的洪水保险赔偿金额或修缮成本等）。由于缺乏建筑结构灾损和洪水流速的历史记录，因此这类洪灾脆弱性曲线常常反映的是建筑内部经济损失与淹没水深的关系。一些研究中虽然考虑了洪水对建筑结构的破坏作用^［22-23］，但是仅给出了建筑完全损毁的临界条件，难以直接用于风险评估。Penning-Rowsell等^［24］指出，流速在对室内经济损失的评估中带来很大的误差，在建筑结构损伤的评估中不可忽视。在我国很多村镇建筑本身质量较低，又面临着山洪的威胁，因此对村镇建筑承载性能开展洪灾风险评估具有重要意义。

以汶川县寿溪河流域村镇建筑为研究对象，开展洪灾风险评估。一方面，寿溪河流域地处山区，村镇在空间分布上离散化程度高，全面获取村镇建筑信息难度较大，难以采用指标体系法；另一方面，流域内洪水流速较大、冲击性较强，现有洪灾脆弱性曲线主要用于室内经济损失评估，而且由于无法获取流速信息，因此不能很好地反映村镇建筑的结构灾损情况。针对以上问题，建立耦合的一、二维水动力模型计算超标溃堤洪水下村镇建筑处的水深与流速，并基于结构静力学原理和洪水冲击荷载经验公式建立村镇建筑的脆弱性曲线，最后根据以上结果对流域内的村镇建筑开展洪灾风险评估。

1 研究区域概况

汶川县地处四川省阿坝藏族羌族自治州（下文简称阿坝州）东南部（见图1），杂谷脑河、鱼子溪和寿溪河是县内3条主要河流。寿溪河是长江支流岷江上游右岸的一级小支流，发源于汶川县境内，流域面积632 km²，主河道长约56 km，流域内主要水文站点分布如图1所示。由于特殊的自然地理条件，因此汶川县洪水频发，而且洪水造成的次生灾害较为严重^［25］。寿溪河流域在1989年、2003年、2013年、2017年、2018年、2020年均发生过山洪灾害，造成严重的人员伤亡和房屋损毁。本研究以寿溪河干流沿岸的村镇建筑为对象，开展洪灾风险评估。研究范围主要包括三江镇、水磨镇等地区，以村民自建房为主。

图1  研究区域地理位置及水系分布

Fig.1  Position of studied area and water distribution

2 洪水模拟计算

洪水的淹没水深和流速是造成村镇建筑破坏的主要因素，为精确获得影响村镇建筑洪灾风险评估的水深和流速，采用丹麦水力学研究所（Danish Hydraulic Institution， DHI）MIKE FLOOD耦合模型对洪水要素进行模拟。课题组分别于2019年10月与2020年12月赴汶川考察，通过无人机测量、数据资料收集、现场调查等手段，获取了2019年“8·20”和2020年“8·17”洪灾后的建筑受损情况，发现汶川县的洪水灾害大多是因泥石流阻断河道而造成的洪水漫溢或溃堤，受灾建筑主要分布在河道附近。因此，本研究的洪水模拟为溃堤洪水模拟。

2.1　模型原理

MIKE FLOOD模型主要包括一维模型（MIKE 11）与二维模型（MIKE 21）2个模块，通过侧向连接实现河道内与河道外的水量交换。MIKE 11模型利用Abbott-Inoescu六点隐式有限差分格式求解圣维南（Saint-Venant）方程组，进行水动力模拟和计算。圣维南方程组是反映有关物理定律的微分方程组，包括连续方程和动量方程，如下所示：

$$\left\{\begin{array}{l}\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=q\\ \frac{\partial Q}{\partial t}+\frac{\partial \left(\alpha \frac{Q^{\mathrm{2}}}{A}\right)}{\partial x}+gA\frac{\partial \eta }{\partial x}+\frac{gQ|Q|}{C^{\mathrm{2}}AR}=\mathrm{0}\end{array}\right.$$

式中：t为时间；x为距离坐标；Q为过流流量；q为侧向入流；A为过水断面面积；η为水位；R为水力半径；C为谢才系数；g为重力加速度；α为动量修正系数。

基于三向不可压缩和雷诺数均匀分布的Navier-Stockes方程，MIKE 21模型控制方程为二维非恒定浅水方程组，该方程组由1个连续性方程和2个动量方程组成，如下所示：

$$\left\{\begin{array}{l}\frac{\partial D}{\partial t}+\frac{\partial (D\overline{\mu })}{\partial x}+\frac{\partial (D\overline{v})}{\partial y}=DS\\ \frac{\partial (D\overline{\mu })}{\partial t}+\frac{\partial (D{\overline{\mu }}^{\mathrm{2}})}{\partial x}+\frac{\partial (D\overline{\mu \nu })}{\partial y}=\overline{\nu }D\omega -gD\frac{\partial \eta }{\partial x}-\frac{D}{{\rho }_{\mathrm{0}}}\frac{\partial p_{\mathrm{a}}}{\partial x}-\frac{g{D}^{\mathrm{2}}}{\mathrm{2}{\rho }_{\mathrm{0}}}\frac{\partial \rho }{\partial x}+\frac{{\tau }_{\mathrm{s}x}-{\tau }_{\mathrm{b}x}}{{\rho }_{\mathrm{0}}}+F_{\mu }+D{\mu }_{\mathrm{s}}S\\ \frac{\partial (D\overline{\nu })}{\partial t}+\frac{\partial (D{\overline{\nu }}^{\mathrm{2}})}{\partial y}+\frac{\partial (D\overline{\mu \nu })}{\partial x}=\overline{\mu }D\omega -gD\frac{\partial \eta }{\partial y}-\frac{D}{{\rho }_{\mathrm{0}}}\frac{\partial p_{\mathrm{a}}}{\partial y}-\frac{g{D}^{\mathrm{2}}}{\mathrm{2}{\rho }_{\mathrm{0}}}\frac{\partial \rho }{\partial y}+\frac{{\tau }_{\mathrm{s}y}-{\tau }_{\mathrm{b}y}}{{\rho }_{\mathrm{0}}}+F_{\nu }+D{\nu }_{\mathrm{s}}S\end{array}\right.$$

式中：x、y为笛卡尔坐标系的坐标；D为总水深，D=d+η，其中d为静水深；μ、ν分别为x、y方向上的流速分量，$\overline{\mu }$、$\overline{\nu }$为平均流速；ω为科氏力系数；ρ为水的密度；ρ₀为水的相对密度；p_a为大气压强；S为点源流量；μ_s、ν_s为源项水流流速；τ_sx、τ_sy分别为x、y方向上的水面风剪应力分量；τ_bx、τ_by分别为x、y方向上的河床切应力分量；F_μ、F_ν分别为x、y方向上的水平涡黏应力分量。

2.2　一维模型的建立

综合考虑寿溪河流域内水系的条件与特点，一维模型概化11条河流。西河、中河、黑石江处为模型上边界，采用流量过程；模型下边界取在郭家坝水文站下游约3 km处，采用水位‒流量关系曲线。根据2013年7月9日20时至7月11日18时洪水期间（下文简称“710”洪水）郭家坝水文站实测水位过程进行模型率定（将模拟开始时刻设为计时零点，下同），西河、中河、黑石江和寿溪河上游河道糙率取值为0.049~0.050，寿溪河中下游河道糙率取值为0.045~0.048。

2.3　二维模型的建立

考虑到山区河道洪水淹没范围有限的特点，二维模型的边界大致按等高线（1 120 m）划定，对村镇建筑所在区域进行网格加密（空间步长5 m）。高程选用2015年的高精度（5 m）数字高程模型（DEM）数据（见图1）；模型采用动边界处理，干水深、淹没水深和湿水深分别取为0.005 m、0.050 m和0.100 m；糙率根据不同的土地利用类型按经验值取定^［26］，村镇用地为0.017、耕地为0.03、河滩为0.04、草地为0.05、林地为0.18。考虑村镇建筑群对水流的作用，对建筑所在处网格进行局部加密。

2.4　耦合模型的建立

耦合模型中设定Side Structure连接模拟河堤溃决。寿溪河流域河道防洪标准普遍较高，可防十年一遇洪水。根据历史洪灾数据和实地调研结果，模拟工况假定溃口发生在三江镇（位于108.34°E、30.92°N，历史溃口）与水磨镇（位于103.42°E、30.93°N，堤防薄弱）附近，当洪水达十年一遇时发生溃堤，分别计算重现期为20年、50年、100年下研究区域内的水力要素，并对村镇建筑开展洪灾风险评估。

2.5　模型验证

选用2017年8月27日23时至8月29日11时洪水期间（下文简称“827洪水”）郭家坝水文站实测水位过程对模型进行验证，将模拟开始时刻定位计时零点，结果如表1所示。模型计算水位与实测值的误差均在0.02 m以下；在较长的模拟时段内（36 h以上）洪峰时刻误差较小（20 min以内），表明该耦合模型能够模拟寿溪河流域洪水溃坝过程，可用于本研究中村镇建筑的洪灾风险评估。

表1  水动力模型率定验证结果

Tab.1  Verification of hydrodynamic model

洪水场次	模拟 时长/h	结果	洪峰 水位/h	洪峰 时刻/h
2013年“710”洪水	46	实测	902.83	10.00
		模型	902.81	10.15
2017年“827”洪水	36	实测	902.41	10.00
		模型	902.42	9.64

2.6　洪水计算结果

图2和表2为各重现期洪水下郭家坝水文站的水位过程和溃口处的体积流量过程（重现期10年时不发生溃堤，图2中省略此时的流量过程）。随着洪水重现期增大，洪水水位升高、洪峰时刻提前，溃口发生时刻提前、溃口处峰值体积流量增大、溃口处体积流量达峰时刻也提前。当洪水重现期为20年时，三江镇、水磨镇溃口分别发生在模拟时间段的17.93 h、7.67 h时刻；重现期提高到50年时，分别提前到15.87 h、6.33 h时刻；重现期为100年时，分别提前到14.93 h、4.57 h时刻。流域内上游河道洪水呈现单峰型，下游河道洪水呈现双峰型，随着洪水重现期提高，双峰中前后峰的相隔时间逐渐增大。水磨溃口发生在前峰，此时洪水流量并不大，而在后峰时显著变大。

图2  验证站点与溃口处水文过程

Fig.2  Hydrologic process of verification station and dam-breaks

表2  一维模型计算结果

Tab.2  Results of one-dimentional model

重现期	郭家坝水文站 洪峰水位/m	郭家坝水文站 洪峰时刻/h	三江镇溃口峰值体积 流量/（m3·s-1）	水磨镇溃口峰值体积 流量/（m3·s-1）
20年	903.30	19.37	61.59	23.04
50年	903.92	17.40	104.88	60.53
100年	904.32	16.33	153.15	90.60

图3为洪水淹没范围，当洪水重现期分别为20 年、50年、100年时，淹没范围分别为31.035 hm²、133.203 hm²、186.525 hm²。寿溪河流域地处山区，地形坡度较大，而河岸两侧地势较高，洪水影响范围十分有限，淹没区域呈狭长的条状。三江镇溃口的洪水在下游约1.0~2.0 km处、水磨镇溃口的洪水在下游约1.5~3.0 km处，均形成范围较大的主淹没区，在其余地方洪水主要影响河岸沿边区域。在模拟时间段内，三江镇发生溃堤，洪水影响范围在五十年一遇时可到达河坝村、在百年一遇时逼近水磨镇。

图3  各重现期下洪水淹没范围

Fig.3  Inundation areas in different return periods

3 脆弱性曲线的建立

洪水流速对建筑结构损伤具有重要影响，在村镇建筑洪灾风险评估中有必要合理考虑水深和流速的共同作用。

本研究的脆弱性曲线基于结构静力学原理与洪水荷载经验公式建立，将建筑结构灾损程度与洪水水深、流速建立函数关系，结合水动力模型成果可以快速开展洪灾风险评估。根据实地考察，洪水对寿溪河流域村镇建筑的破坏主要集中在墙体的冲击上，即墙体的弯曲破坏和剪切破坏。参考Fang等^［20］的研究，假定洪水流速在竖直方向上均匀分布，考虑如图4所示的一次超静定墙体，取单位宽度进行分析。

图4  洪水冲击荷载作用下墙体受力情况

Fig.4  Stress condition of wall under flood impacting

洪水冲击荷载与水深、流速的关系采用孙云飞^［27］提出的砌体结构洪水冲击荷载计算公式，如下所示：

$$P(x,v)=\mathrm{9}\mathrm{637.63}(h-x)+\mathrm{391.07}{v}^{\mathrm{2}}+\mathrm{905.91}v-\mathrm{206.63}=A_{\mathrm{0}}(h-x)+f(v)$$

式中：P为冲击荷载，Pa；A₀为冲击荷载中水深项的系数，取值9 637.63；f（v）为冲击荷载中流速项的作用，f（v）=391.07v²+905.91v-206.63；h为淹没水深。根据图4，洪水冲击荷载引起墙体破坏的因素包括单宽底部剪力F_m1（N·m^-1）、单宽顶部剪力F_m2（N·m^-1）、单宽底部弯矩M_m（N），易知F_m1>F_m2。F_m1和M_m的计算式如下所示：

$$F_{\mathrm{m}\mathrm{1}}=\left|\frac{A}{\mathrm{40}{L}^{\mathrm{3}}}{h}^{\mathrm{5}}+\frac{f(v)-AL}{\mathrm{8}{L}^{\mathrm{3}}}{h}^{\mathrm{4}}-\frac{f(v)}{\mathrm{2}{L}^{\mathrm{2}}}{h}^{\mathrm{3}}+\frac{A}{\mathrm{2}}{h}^{\mathrm{2}}+f(v)h\right|$$

$$M_{\mathrm{m}}=\left|\frac{A}{\mathrm{40}{L}^{\mathrm{2}}}{h}^{\mathrm{5}}+\frac{f(v)-AL}{\mathrm{8}{L}^{\mathrm{2}}}{h}^{\mathrm{4}}-\frac{\mathrm{3}f(v)+AL}{\mathrm{6}L}{h}^{\mathrm{3}}+\frac{A+f(v)}{\mathrm{2}}{h}^{\mathrm{2}}\right|$$

式中：L为单层墙体高度，根据实地考察结果取L=3.5 m。考虑墙体自重，根据《砌体结构设计规范》（GB 50003—2011）中受弯构件的计算要求^［28］，假定墙体砂浆强度为M10，则墙体极限单宽承载剪力F_c（N·m^-1）和极限单宽承载弯矩M_c的计算式如下所示：

$$F_{\mathrm{c}}=(f_{\mathrm{v}}+{\rho }_{\mathrm{w}}gL)\frac{\mathrm{2}}{\mathrm{3}}{t}_{\mathrm{w}}$$

$$M_{\mathrm{c}}=(f_{\mathrm{t}\mathrm{m}}+{\rho }_{\mathrm{w}}gL)\frac{\mathrm{1}}{\mathrm{6}}{t}_{\mathrm{w}}^{\mathrm{2}}$$

式中：ρ_w为墙体密度，取值为1 800 kg·m^-3；t_w为墙体厚度，取值为0.24 m；f_v、f_tm分别为烧结普通砖砌体的极限抗剪强度和弯曲抗拉极限强度，取值分别为0.27 MPa和0.27 MPa。

定义砌体结构建筑洪灾脆弱性为

$$\xi =\mathrm{m}\mathrm{i}\mathrm{n}\left\{\mathrm{m}\mathrm{a}\mathrm{x}\left\{\frac{F_{\mathrm{m}\mathrm{1}}}{F_{\mathrm{c}}},\frac{M_{\mathrm{m}}}{M_{\mathrm{c}}}\right\},\mathrm{1}\right\}$$

图5为ξ与h和v的三维曲面函数关系。考虑到国内外常使用“脆弱性曲线（vulnerability curve）”的概念，将三维曲面看作是二维曲线的推广，因此仍称作脆弱性曲线。ξ是砌体结构所受洪水冲击荷载与其极限荷载的比值，反映砌体结构在洪水作用下的安全冗余度。当 ξ=1.0时，表明结构已经达到极限状态并将发生破坏。图5表明，在流速v=0 m·s^-1下，淹没水深h达到0.81 m时建筑为高风险、达到1.12 m时发生破坏；在流速 v=3 m·s^-1下，h达到0.58 m即为高风险、达到0.74 m即发生破坏。该结果与美国陆军工程师团^［23］（淹没水深达0.9 m时建筑大概率破坏，流速达3 m·s^-1时需采取专门措施）、Xiao等^［29］（均布荷载为23.7 kPa时结构破坏，本研究中ξ=1.0对应的均布荷载为20.8 kPa）的研究成果相近。当h≤2.0 m时，破坏状态为弯矩引起的受拉破坏；当2.0 m<h<2.5 m时，破坏状态转变为弯矩引起的受剪破坏，此时安全系数较高，脆弱性降低；当 h≥2.5 m时，洪水产生的静水压强较大，脆弱性较高。

图5  砌体结构脆弱性曲线

Fig.5  Vulnerability curve of masonry structure

课题组于2020年12月赴汶川实地考察，走访了绵虒镇草坡乡、大禹村、羌锋村，三江镇草坪村、街村和水磨镇马家营村6个村庄。在调研的415间村镇建筑中，钢混结构、砖混结构、砖木结构和土石木结构分别有214间、158间、22间和21间。向当地居民了解后得知，汶川县内各村镇的建筑类型相近、材料相似，其中所有的砖木、砖混结构和绝大多数的钢混结构，其外墙均为砖砌体墙，数量约占调研建筑总数的94.5%。考虑到村镇建筑的破坏以墙体破坏为主（见图6），本研究所建立的砌体墙脆弱性曲线能够较好地代表研究区域内的建筑特性。

图6  不同类型建筑墙体破坏实例

Fig.6  Examples of damaged building walls in different structures

4 村镇建筑洪灾风险评估

村镇建筑的洪水模拟结果与村镇建筑的脆弱性曲线相结合即可进行村镇建筑的风险评估。首先，根据水动力模型的计算结果，按式（8）计算各网格处的脆弱性；然后，对每一幢村镇建筑，取其周围网格脆弱性数值，按网格面积加权计算该建筑的脆弱性，并根据其脆弱性分为低风险、中风险、高风险、极高风险4类。将ξ=1.0作为极高风险的判别标准；当0≤ξ<1.0时，基于等距分割原则划定低、中、高风险，各风险等级对应的水深和流速如图7所示。

图7  风险等级划分

Fig.7  Risk class division

基于卫星遥感影像（2015年数字正射影像）对模型区域内村镇建筑进行概化，共计1 688幢。研究区域内各重现期洪水下村镇建筑的评估结果如图8~10所示。与洪水淹没范围相对应，受洪水影响的村镇建筑主要集中在离河道轴线100 m的范围内，20 年、50年和100年重现期下受影响的建筑分别为19幢、51幢、88幢。洪水重现期为20年时，没有极高风险的建筑，中高风险建筑集中分布在溃口下游的主淹没区内；重现期达50年时，中风险建筑分布在溃口附近和溃口下游约2 km处，高风险和极高风险建筑主要出现在溃口下游的主淹没区内；重现期达100年时，高风险和极高风险建筑主要分布于洪水发生点、溃口下游的主淹没区、溃口附近以及河坝村地区。

图8  村镇建筑洪灾风险评估结果（重现期20年）

Fig.8  Results of flood risk assessment on rural buildings (scenario of 20-year return period)

图9  村镇建筑洪灾风险评估结果（重现期50年）

Fig.9  Results of flood risk assessment on rural buildings (scenario of 50-year return period)

图10  村镇建筑洪灾风险评估结果（重现期100年）

Fig.10  Results of flood risk assessment on rural buildings (scenario of 100-year return period)

根据村镇建筑的洪灾风险评估结果，对寿溪河流域村镇建筑的防洪减灾提出如下建议：

（1）溃口对临河距离产生的不良影响较大，应尽快采取措施探明堤防薄弱段，并加以修复与改建。堤防相对薄弱处下游2 km范围内的建筑应采取必要措施，加强排水能力或降低洪水流速，如修建防洪围墙。

（2）村镇建筑的布局应距离河道中轴线100 m以上，这样既可以给河流充分的调洪空间，又能降低居民受到洪灾威胁的风险。

（3）流速对建筑的损伤有较大影响，应定期对村镇建筑结构安全开展检查。

（4）在临河建筑多、洪水威胁较大的寿溪河流域，应注重事前风险控制，提高“四预”能力，即加强洪水预报、提高预警时效、加强专项预案编制以及拟定转移路线，并针对洪灾抢险和撤离定期开展预演。

5 结论

（1）利用水动力模型模拟了寿溪河流域的洪水过程。研究区域内各重现期下的洪水淹没范围均呈狭长的条状，在溃口下游约2 km处形成面积较大的淹没区域，在其余区域洪水主要影响河道沿岸区域。模型结果合理准确（水位误差均在0.02 m以内），能够用于村镇建筑的洪灾风险评估。

（2）基于结构静力学原理建立了砌体结构的洪灾脆弱性曲线，将砌体建筑的脆弱性表示为水深和流速的函数曲线。该曲线考虑了洪水流速的影响，能够表征砌体建筑在洪水作用下的结构灾损严重程度，结合水动力模型结果即可快速开展大范围的村镇建筑洪灾风险评估。

（3）结合洪水计算结果与脆弱性曲线对寿溪河流域的村镇建筑开展了洪灾风险评估。洪水重现期为20年时（郭家坝水位903.30 m）区域内村镇建筑以中低风险为主；重现期为50年时（郭家坝水位903.92 m）高风险和极高风险建筑集中分布于三江镇和水磨镇溃口下游的主淹没区内；重现期达100年时（郭家坝水位904.32 m），高风险和极高风险建筑出现在主淹没区内、洪水发生点和河坝地区。

（4）寿溪河流域的村镇建筑受溃堤洪水影响较大，在防洪减灾时应首要考虑对堤防薄弱处进行修复与加固，在此基础上提高“四预”能力并采取措施降低沿河村镇建筑周围的洪水流速。

作者贡献声明

刘曙光：提供研究思路，校核公式，撰写论文。

郑伟强：数据整理，公式推导，结果处理与分析，论文撰写。

钟桂辉：数据校核，结果分析，论文修订。

甄亿位：数据整理、校核，论文修订。

韩 征：数据校核，论文修订。

李艳鸽：数据校核，论文修订。