王 琳,苑鹏飞,钟启明,胡 亮,单熠博,薛一峰
(1. 西安理工大学 省部共建西北旱区生态水利国家重点实验室,陕西 西安 710048; 2. 南京水利科学研究院,江苏 南京 210029;3. 水利部土石坝破坏机理与防控技术重点实验室,江苏 南京 210029; 4. 陕西省水利电力勘测设计研究院, 陕西 西安 710001)
堰塞湖是受地震活动、气候变化、人类活动等影响的一种典型地质灾害[1-2]。据统计[3-8],9%的堰塞湖在1 h内溃决,34%在1 d内溃决,67%在1个月内溃决。其一旦溃决,将形成溃决-致灾灾害链,威胁下游沿岸的工程设施及人民群众生命财产安全。目前,我国堰塞湖溃决灾害频发,例如:2000年4月,易贡堰塞湖[9-11]形成,6月溃决,溃决时峰值流量达到94810 m3/s。2008年最具代表性的唐家山[12-13]和小岗剑堰塞湖[14],溃决洪峰分别达到6500、3950 m3/s,使接近百万人的生命财产安全受到威胁。2018年白格堰塞湖更是2次堵塞金沙江上游干流[15-17],溃决洪峰流量分别达到10000、31000 m3/s[18-19],10月10日形成的堰塞湖更是在3 d内就发生溃决,11月3日堰塞湖溃决导致下游在建的苏洼龙水电站围堰拆除,在建的叶巴滩、拉哇和巴塘等水电站工程延期。堰塞湖应急处置易受到基础资料收集困难、周边环境危险、分析研究时间紧和溃决灾害扩大效应强等因素制约,由此导致应急处置难度极大[20-23]。因此,考虑其致灾灾害链效应,基于极其有限的数据信息建立堰塞湖溃决风险快速、定量评估方法是目前应急处置亟需解决的难题。
目前,堰塞湖溃决风险评估研究主要分为两类:统计分析和定量分析。统计分析是基于堰塞湖统计信息,评估溃决灾害损失,划分不同类型灾害风险等级,代表性的有水利部行业标准SL/T 450—2021《堰塞湖风险等级划分与应急处置技术规范》[24]。该标准规定,可根据堰塞湖坝体规模、物质组成和高度,运用查表法和模糊数学法将其风险划分为极高、高、中和低4级。崔鹏等[25]和XU等[26]基于汶川地震堰塞湖溃决实例,根据堰塞湖坝体高度、库容和物质组成3项指标提出风险等级划分建议;XU等[27]采用模糊数学方法对红石岩堰塞湖的风险等级进行了评估。定量分析模型多基于堰塞湖溃决洪水淹没情况与损失之间的关系构建评估模型。代表性的有石振明等[28]、SHI等[29]构建的综合考虑人员疏散情况、生命损失和溃决参数等因素的定量分析模型。
以上两类评估方法均可实现溃决风险评估。第一类方法主要是基于经验的定性分析,对堰塞湖风险快速评估具有较好的指导作用,但其主要考虑堰塞湖坝体地貌学特征,未考虑坝体的物质组成结构。而目前不同诱因导致堰塞湖坝体级配特征存在极大差异,若未合理考虑其物质组成结构,将无法准确评估堰塞湖风险[30]。第二类方法虽然综合考虑了水力学参数和灾害损失,但其风险评估过程中并未开展坝体稳定性快速评价,无法快速实现稳定性的量化评估;也并未考虑到溃决机理对溃决特征值变化的影响机制,导致其溃决-致灾的灾害链效应尚未明晰。
本文构建了一套堰塞湖溃决-致灾的快速定量评估方法,实现快速定量风险评估。本方法基于极其有限的数据信息构建,考虑了堰塞湖和坝体相关物理力学指标、下游河道和影响区地形数据、人口分布情况,实现了考虑颗粒级配组成的稳定性快速评价、基于冲蚀特性和崩塌过程的溃口流量变化过程分析、洪水演进过程模拟、生命损失评估预警,明晰了溃决-致灾的灾害链效应。该方法包括五部分:三维地理信息模块、稳定性快速评价模块、溃口流量快速定量分析模块、溃决洪水演进模块和生命损失快速评估模块,将该方法应用于拥有完整溃决及演进过程实测资料的唐家山堰塞湖案例,一方面验证了本文模型的可行性,同时也为堰塞湖的风险评估和应急处置提供技术参考。
堰塞湖溃决快速定量风险评估方法主要包括:信息收集—稳定评价—溃决洪水分析—风险评估—损失预警五部分。详细步骤如下:①运用HEC-GeoRAS插件快速获取研究区域的三维地形数据,建立坝体和河道的三维地理信息模块,实现信息快速收集;②基于逻辑回归算法,建立考虑堰塞湖坝体形态特征、物质组成和上游堰塞湖水动力条件的稳定性快速评价模块,实现堰塞湖坝体稳定性快速评价;③基于堰塞湖溃决机理,开发基于物理机制并可在1 h内快速计算溃决洪峰流量、溃口崩塌过程、溃决时间和溃决库容等数据的溃决洪水分析模块,实现溃决参数变化过程快速预测;④将三维地理信息数据和溃决洪水定量分析数据导入HEC-RAS软件,构建溃决洪水演进计算模块,开展风险评估;⑤在地理信息系统(geographic information system, GIS)中构建可视化溃决洪水淹没范围图,实现风险快速评估预警,揭示溃决洪水影响区域内的生命损失分布特征。
基于中科院地理空间数据云(http://www.gscloud.cn)获取堰塞湖区域范围内的高程数据,运用HEC-GeoRAS插件快速获取溃决及演进范围的三维地理信息。主要步骤如下:①基于中科院地理空间数据云平台获取堰塞湖具体的经纬度信息及其等高线数据图;高程数据精度越高,数据量几何倍增加,数据处理效率越低,因此,在实际应用时,应根据具体工程实际选择所需数据精度。对于堰塞湖规模、下游河道信息和周边城镇地区等数据,使用精度30M的高程数据可满足应用需求。对于坝体高度、宽度和坝体体积等数据,可在堰塞湖所在经纬度范围内使用精度为12M的高程数据。②运用HEC-GeoRAS插件在GIS平台将河道数据、高程数据和其它三维地形信息数据直接导入并叠加,实现图像数据和数字数据的转化[31]。迅速获取堰塞湖的几何特征、下游影响区域等关键数据(1 h内)。堰塞湖实际应急处置时间非常紧迫,例如唐家山堰塞湖在1个月内溃决,10月10日形成的白格堰塞湖更是在3 d内发生溃决。应急处置面临着基础资料收集少、分析研究时间紧的困难,在满足精度的要求下快速处理三维地理信息,能够为后续稳定性评价-溃决洪水分析等一系列工作提供数据支持。
本模块采用的稳定性快速评价方法,考虑了堰塞湖坝体级配特征、物质组成以及上游的水动力条件,而现有评价方法多未考虑堰塞湖坝体物质组成。该方法可根据应急处置时获取的级配特征(如块石、砾石、粗细颗粒等),实现精细化和简化稳定性快速评价。具体表达式为
Ls(IVAM)=-0.1981 lg(I)+1.387 lg(Vd)-1.432 lg(Ab)+4.169Mi-8.674
(1)
式中:Ls(IVAM)为堰塞湖坝体稳定性评价指标,当Ls(IVAM)<0时,堰塞湖坝体可判定为不稳定;当Ls(IVAM)>0时,堰塞湖坝体可判定为稳定。I=Hd/W,其中Hd为堰塞湖坝体高度,W为堰塞湖坝体宽度,即顺河向距离;Vd为堰塞湖坝体体积,Ab为堰塞湖流域面积;Mi为颗粒特征参数,Mi的取值方法在文献[32]中具体阐述。
该稳定性快速评价方法的绝对准确率、保守准确率、判错率分别为:88.06%、94.03%、 5.97%,同时与国内外同类方法比较[32](BI法[33]、DBI法[34]、Is法和Ia法[35]),结果表明该方法具有可靠性与优越性。其具体比较过程已在文献[32]中详细说明,文中不再进行赘述。
堰塞湖溃决风险评估方法中的溃决洪水定量分析模块多基于历史上的溃决资料,未考虑溃决机理。虽然参数模型可利用较短计算时间完成溃决洪峰流量的模拟,但无法获取溃决流量过程,且未考虑土水耦合作用机制和冲蚀破坏过程[36]。溃决洪水快速定量分析是风险评估和应急处置方案制定的基础。因此,构建基于溃决机理,考虑坝体材料冲蚀特性和溃口崩塌过程的快速、定量溃决过程计算模型是提高堰塞湖溃决洪水风险评估准确性的必然需求。
本文针对20世纪以来典型的唐家山[37]、易贡[9]、小岗剑[14]、红石岩[38]和白格[39]等堰塞湖的实测溃决数据进行梳理,并在此基础上提出了基于物理机制的溃决洪水分析模型DB-IWHR(https://github.com/ChenZuyuIWHR/DB-IWHR)免费下载使用,如图1所示,H为水库水位(m);t为时间(s);W为库容(m3);B为溃口宽度(m);m为跌落系数;Z为溃口底高程(m);h为水深(m);C为流量系数(m1/2);a、b为双曲线侵蚀模型系数;v、vc分别为流速和初始流速(m/s);γ为容重(kN/m3);c为黏聚力(kPa);φ为摩擦角(deg);m1、m2为双曲线倾角扩展系数;β、β0分别为溃口倾角和初始倾角(°)。此模型采用宽顶堰公式计算溃决流量,基于圆筒型冲蚀试验设备(CETA)的试验数据提出双曲线形式的冲蚀本构模型;在计算溃口崩塌时基于总应力法,采用双曲线图表形式,快速确定和模拟溃口崩塌过程。当溃口不断被冲深时,发生侧向破坏,可通过溃口边坡稳定分析逐级确定临界坡面和滑裂面[40-41]。DB-IWHR溃决分析程序,可保证工程师在1 h内完成目标工况计算,在堰塞湖应急处置时可快速应用。
图1 堰塞湖溃决洪水分析模型(DB-IWHR)Fig.1 Dam breach analysis model for barrier lake (DB-IWHR)
将HEC-GeoRAS插件中构建的河道断面三维模型导入HEC-RAS软件构建河道模型,将其分为左漫滩、中心区域和右漫滩3个区域。确定不同区域的曼宁系数以及演进计算的上下游边界,确定溃决及演进区域的中心线、河岸线、行洪区线、区域面积和糙率等数据。溃决演进计算采用开源的HEC-RAS 4.1.0软件的非恒定流模型计算,实现溃决洪水演进的快速准确模拟。基于非恒定流的方程采用差分法进行求解,通过划分河道断面、水流路径等模拟河道信息,水位-流量关系为边界条件,时间步长设置为30 s,基于四点隐式差分法进行迭代计算。
非恒定流模拟是基于水流的连续方程与动量方程展开的,其连续方程为
(2)
式中:ρ为液体密度;v为流速。
动量方程为
(3)
式中:fi为质量力;p为压力;υ为运动黏滞系数。
堰塞湖溃决产生的洪水会对下游地区造成严重损失,其中以生命损失最为重要[42-43]。溃决洪水造成生命损失的致灾环境复杂,亟需在堰塞湖应急处置时开展预警,计算溃决可能造成的损失,及时撤离群众[44]。目前的生命损失评估方法多存在参数量化过程不统一,影响因素定义模糊,评估精度不高等问题。本文方法基于极限学习机网络模型建立风险人口与生命损失间的函数关系,考虑了主要影响因素与次要影响因素间的作用关系与影响程度;并通过对生命损失影响因素的统一量化,结合前文计算的溃决洪水定量分析数据与地理信息数据,获取溃决洪水深度、淹没面积和水位变化过程,绘制淹没范围图,确定溃决洪水到达城镇时间以及淹没水深和淹没范围。
在堰塞湖灾害所造成的生命损失中,风险人口和人口密度等人口统计学变量常作为关键参数参与生命损失模型的计算,如王志军等[45]的模型。在实际运用中,由于城镇的人口分布通常较为集中,在城镇中心区域的人口密度通常较大,其他区域人口密度较小,对生命损失计算影响较大。在本模型中,通过溃口流量与洪水演进模块所计算出的下游淹没水深,能够快速预估出下游淹没范围,并将该范围作为生命损失模块的计算起始条件。在这个范围内,根据人口密度将受灾区域进一步划分为高密度区域、低密度区域[46],不同密度区域中的修正系数α、β取值不同。在计算过程中,根据式(4)对各区域分别计算求和,最终获得淹没范围内的生命损失L为
L=α×β×10A·lgPAR+B·(lgPAR)2
(4)
式中:L为生命损失;α为主要因素修正系数;β为次要因素修正系数;A为人口密度参数;B取0.06。
2008年5月12日,受汶川8.0级地震影响,四川通口河右岸发生特大滑坡,堵塞河流形成唐家山堰塞湖。其距下游北川县城6 km,绵阳市70 km,下游有将近130万人民,还有宝成铁路、输油管道等重要基础设施,如图2所示[12]。为减少唐家山堰塞湖溃决带来的巨大灾害,应急处置开挖了宽为8 m、深为13 m、总长为475 m的引流槽,开挖后的堰塞湖槽底高程为740.4 m。2008年6月8日,槽内水深和宽度开始逐步增加,此时上游来流量为80 m3/s,堰塞湖水位逐渐抬升。堰塞湖水位至6月10日01∶30达到最大值743.1 m。从01∶30到09∶30,测量的平均流速值在2.3~2.9 m/s之间,相应的水流剪应力在20~30 Pa之间。堰塞湖溃决时颗粒的启动流速大约为2.6 m/s。至2008年6月10日12∶36,溃口达到峰值流量,约为6500 m3/s。溃决过程中槽底的平均侵蚀率约为0.5 mm/s,对应的流速为2.4~4.9 m/s,剪应力为20~80 Pa。6月10日20∶00,溃口流速(约2.4 m/s)及流量基本保持不变,表明堰塞湖溃决基本结束,此时堰塞湖水位为720.3 m,整个过程中唐家山堰塞湖的下泄水量约2.3×108m3[12-13]。溃口最终深度为42 m,顶部宽度约为145~235 m,底部宽度约为80~100 m,溃口洪水最大流速为4.96 m/s[47]。6月10日15∶15,洪水到达绵阳市,未出现较大洪峰,排险成功结束。唐家山堰塞湖拥有珍贵的溃决全过程实测资料,为堰塞湖溃决研究提供了翔实的基础数据。堰塞湖相关参数如表1所示。下游北川(距唐家山坝址7.5 km)、通口(距唐家山坝址31.3 km)、涪江桥(距唐家山坝址78.5 km)等水文站的实测洪水特征值见表2[12]。
表2 唐家山堰塞湖及下游水文站的水文特征值Table 2 Hydrologic characteristics of Tangjiashan barrier lake and downstream gauging stations
根据前述堰塞湖溃决风险快速定量评估方法,将其应用于唐家山堰塞湖风险评估,具体模拟方法如下。
唐家山堰塞湖位于通口河,沿通口河向下游途经北川县城、通口镇、含增镇、青莲镇、龙凤镇、石马镇以及绵阳市[48]等城市和乡镇,约有120万居民。因此,选取的溃决计算及演进计算区域为唐家山堰塞湖所在地沿通口河至绵阳市的矩形区域(104°24′E~104°48′E、31°29′N~31°50′N),基于HEC-GeoRAS插件共确定20个断面,河道中心线长度为84.79 km。
唐家山堰塞湖底高程约为640 m,计算区域高程急剧下降约200 m,如图3所示,另外,图3标注了沿河段的3个典型水文站位置。将基于HEC-GeoRAS插件构建的河道模型导入HEC-RAS软件,构建溃决洪水演进模型。在唐家山至绵阳市河道选择20个断面作为典型断面,其中1号、4号、19号、20号断面分别为唐家山堰塞湖、北川水文站、通口水文站、涪江桥水文站和绵阳市如图4所示。
图3 通口河段剖面图
图4 溃决洪水演进模型构建
唐家山堰塞湖坝体高度为85 m,堰塞湖坝体宽度为802 m、堰塞湖坝体体积为3.26×108m3,堰塞湖流域面积为3550 km2。其颗粒组成以碎石为主,粒径多小于20 cm,粉质黏土达到60%,碎石约为30%~35%,块石约为5%~10%,颗粒特征参数Mi取为0.75。将各参数代入式(1),获得Ls(IVAM)=-10.87<0,堰塞湖坝体稳定性判断结果为“不稳定”,亟需开展溃决洪水模拟分析。
采用DB-IWHR开展溃决洪水过程定量分析,输入参数如表3所示。图5和表4为溃口流量、堰塞湖水位、溃口水流流速、溃口宽度、溃口底高程的实测值和计算值的对比情况。反演分析结果显示,6月10日12∶04,溃口出现峰值流量,为7322.79 m3/s,溃决库容达到2.09×108m3,溃口顶部宽度为131.2 m,溃决水流最大流速为6.36 m/s。若不开挖宽为8 m、深为13 m的引流槽,则在坝顶754 m时堰塞湖将溃决,库容将达到3.14×108m3,峰值将达到9343.35 m3/s,溃口顶部宽度为151.6 m,溃决最大流速为6.66 m/s。开挖引流槽可使溃口峰值流量减少12.6%,洪量减少36.5%。反演分析的误差均在12%以内[26],因此可以认为反演结果能够验证模型的适用性,能够在保证精度需求的前提下进行快速反演。因此,运用基于物理机制的溃决洪水分析程序可快速、定量地评价堰塞湖溃决过程。
表3 溃决过程反演分析DB-IWHR输入参数Table 3 DB-IWHR input parameters of inversion analysis of breach process
图5 唐家山堰塞湖溃决计算结果与实测数据对比Fig.5 Comparison of the measured and calculated data for Tangjiashan barrier lake
表4 唐家山堰塞湖溃决参数计算结果对比Table 4 Comparison of breach output parameters for Tangjiashan barrier lake
溃决洪水演进分析输入参数如表5所示。将计算结果与具有实测资料的下游3个水文站资料进行分析对比,验证洪水演进模块的可行性。唐家山堰塞湖下游3个水文站反演计算与实测数据对比结果如图6及表6所示。
表5 洪水演进输入参数表Table 5 Inputs parameters for flood routing model
图6 唐家山堰塞湖下游3个水文站洪峰流量过程计算与实测数据对比Fig.6 Comparison of measured and calculated data for three gauging stations downstream of Tangjiashan barrier lake
表6 唐家山堰塞湖下游3个水文站实测与计算结果对比Table 6 Measured and calculated results of three gauging stations downstream of Tangjiashan barrier lake
针对北川水文站,实测北川水文站溃决洪峰时间为6月10日13∶00,溃决流量为6540 m3/s,最高水位为620.29 m。反演分析的溃决洪峰时间为13∶22,溃决流量为6414 m3/s,最高水位为641.00 m。溃决洪峰的误差为1.9%,最高水位的误差为3.2%,到达溃决洪峰时间误差为3.12%。反演分析与实测数据间相关系数为0.9496。
针对通口水文站,实测通口水文站溃决洪峰时间为6月10日13∶00,溃决流量为6210.0 m3/s,最高水位为549.73 m。反演分析的溃决洪峰时间为13∶10,溃决流量为6306.5 m3/s,最高水位为572.50 m。溃决洪峰的误差为1.6%,最高水深的误差为4.1%,到达溃决洪峰时间误差为5.5%。反演分析与预测数据间相关系数为0.9770。
针对涪江桥水文站,实测涪江桥水文站溃决洪峰时间为6月10日15∶18,溃决流量为6100.0 m3/s,最高水位为465.28 m。反演分析的溃决洪峰时间为16∶14,溃决流量为6039.6 m3/s,最高水位为480.50 m。溃决洪峰的误差为1.0%,最高水深的误差为3.3%,到达溃决洪峰时间误差为6.1%。反演分析与实测数据间相关系数为0.9365。
唐家山堰塞湖溃决后,到达下游北川、通口、涪江桥水文站的溃决流量、洪峰到达时间、最大水深等溃决参数的计算误差均在10%以内,且反演数据与实测数据的相关系数均到达0.9以上,在洪峰局部与过流全过程均验证了洪水演进模块的可行性。计算误差的来源主要为地形信息误差、河道糙率数据、断面划分密度、理论计算依据等的影响,这是洪水演进研究中不可避免的,且随着模拟时间与洪水演进,这些缺陷将被不可避免地逐渐放大。
基于溃决洪水演进分析数据,在GIS中可视化并构建溃决洪水淹没范围图,获取溃决后影响区域的淹没范围、淹没水深和淹没历时,评估堰塞湖溃决洪水影响区域内的生命损失分布特征[48]。唐家山堰塞湖溃决后最大淹没范围如图7所示。
图7 唐家山堰塞湖溃决后最大淹没范围Fig.7 The maximum inundation area after the breach of Tangjiashan barrier lake
为规避面临的溃决风险,绵阳市大约有27.5万人在大坝溃决前10 d被安全撤离,并没有人员丧生[49]。本文主要研究若没有进行及时疏散,唐家山堰塞湖溃决后将面临的生命损失。采用生命损失模块预测实际溃决过程和相应的预警时间。基于洪水演进计算结果,唐家山堰塞湖下游6个城镇和绵阳市均有可能遭受溃决洪水的威胁,如表7所示。
不同预警条件下的生命损失计算结果,如表8所示。计算结果表明,生命损失与预警时间存在明显的相关性,若预警时间仅为0.25 h时,生命损失率达到0.25%,当预警时间大于2.0 h时,基本上可做到无生命损失。而唐家山堰塞湖提前10 d预警,保证了无人员丧生。
表7 唐家山堰塞湖下游信息Table 7 Information of the subareas downstream of Tangjiashan barrier lake下游城镇距离/km河岸高度/m人口数量/人北川镇8.0615.830000通口镇31.0547.57300含增镇38.0523.410000青莲镇45.9493.816300龙凤镇55.9474.115000石马镇64.0463.020500绵阳市78.0139.01127000表8 不同警报时间下生命损失评估Table 8 Life loss assessment under different alarm time预警时间/h生命损失/人北川镇通口镇含增镇青莲镇龙凤镇石马镇绵阳市0.25512961221791684819290.7539473941381293714841.031559751101033011871.510000032.00000001
1)本文提出了在堰塞湖突发情况下的快速风险评估方法,该方法分为:信息收集—稳定评价—溃决洪水分析—风险评估—损失预警共五部分,可实现坝体稳定性评价、溃决洪水分析、演进过程快速计算和生命损失分布特征预警的溃决快速定量风险评估。方法考虑了堰塞湖坝体的级配特征,实现精细化和简化稳定性快速评价;考虑溃决机理,建立了基于冲蚀特性和崩塌过程的溃口流量变化过程分析物理模型,实现溃决流量过程快速分析(1 h内);基于极限学习机网络模型考虑风险人口与生命损失间的函数关系,实现了生命损失分布特征预警,明晰了溃决-致灾的灾害链效应。
2)堰塞湖的溃决风险与坝体形态特征、物质组成和堰塞湖水动力条件密切相关,基于有限信息,合理评价堰塞湖坝体的稳定性是溃决风险评估的先决条件。构建基于溃决机理,考虑坝体材料冲蚀特性和溃口崩塌过程的快速、定量溃决过程计算模型,获取水位、溃口流量、溃口宽度、流速和溃口底高程等溃决数据的变化过程是堰塞湖风险评估的关键。
3)基于GIS技术,可视化构建溃决洪水淹没范围图,迅速获取溃决后影响区域的淹没范围、淹没水深和淹没历时,是准确开展风险评估的重要保证。根据溃决洪水的淹没情况,充分考虑洪水影响区域内的人口分布特征,开展生命损失评估是风险评估和抢险救灾的重要目标。
4)在唐家山堰塞湖应急处置中,若不开挖引流槽,则在坝顶754 m时堰塞湖将溃决,库容将达到3.14×108m3,峰值将达到9343.35 m3/s,溃口宽度为151.6 m,溃决最大流速为6.66 m/s。开挖引流槽可使洪峰流量减少12.6%,洪量减少36.5%。因此,对于库容较大的堰塞湖,在条件允许的情况下,可通过开挖引流槽降低溃决风险。
5)本文的评估方法基于堰塞湖和坝体的物理力学特性,有效结合三维地理信息,考虑堰塞湖坝体的稳定性、冲蚀特性和溃口崩塌,以及复杂条件下的溃决洪水演进过程,可快速、定量评估堰塞湖溃决洪水风险,为堰塞湖的应急处置提供技术支撑。
仍需考虑经济损失和环境影响等类别的风险后果,需进一步研究考虑致灾后果的堰塞湖溃决风险评估方法。