第39卷第16期2009年8月数学的实践与认识M A TH EM A T I CS I N PRA CT I CE AND TH EO R Y V o l 139 N o 116
A ugu st,2009 唐家山堰塞湖泄洪问题研究
曾向荣, 郝红星, 孙博良
(国防科技大学信息系统与管理学院,湖南长沙 410073)
摘要: 研究了唐家山堰塞湖的泄洪规律,并对其蓄水、溃坝、泛洪等三个过程进行了模型分析.
关键词: 唐家山;堰塞湖;泄洪;蓄水;溃坝;泛洪
1 问题描述和分析
收稿日期:2009204230
基金项目:国防科技大学数学建模创新实践基地(S M 0802)
地震后形成的唐家山堰塞湖对下游人民的生命财产和国家建设构成了巨大的威胁.因此对唐家山堰塞湖水情的关注是非常必要的.分析堰塞湖的各种水情,采取相应的对策,最大程度的减少人民生命财产和国民经济的损失.为了研究唐家山堰塞湖的泄洪规律,需要分析堰塞湖三个位置、三个阶段的水情,如图1所示,三个位置分别是坝前、坝体、坝后,三个阶段则为蓄水阶段、溃坝阶段、泛洪阶段.
图1 堰塞湖三个位置和三个阶段示意图对于蓄水阶段,可以通过数字高程模型
(D E M ),使用连通性分析算法得到蓄水量和水位
高程关系曲线;又由于存在蓄水量和水位高程变
化的历史资料,也可以采用曲线拟合的方法.降水
量对水位的影响则可以使用调洪演算模型得出.
对于溃坝阶段,溃口流量过程可以使用溃口
冲刷扩展速率,溃口动态出流计算公式联合求解
的溃口流量非恒定过程计算模型进行分析.
对于泛洪阶段,基于D E M ,可以采用搜索算
法进行淹没面积的计算.而洪水到达时间等值的计算则需要在经验公式和基于圣维南方程的数值模拟之间进行选择.
2 蓄水阶段的模型建立与求解
2.1 曲线拟合模型
根据地震期间新闻提供的数据[127],以水位高程为自变量,蓄水量为纵坐标得到分布图
2所示.分析图中点的分布特征可以得到蓄水量与水位高程的关系并不是简单的线性关系,而根据实际情况容易得到,由于水位高程的增大会使蓄水面积增大,因此随着水位高程的增大,需水量的增大速度应当逐步增大.考虑到散点的分布特征,本论文引入二次曲线、指数函数和双曲线等几种预测模型,通过已知点用M A TLAB 对上面三种模型进行拟合,可以得到各模型参数以及误差如表1所示.
表1 曲线拟合参数计算结果模型
参数值误差二次曲线模型
(y =ax 2+bx +c (a ≠0))a =
6.
2945b =-8608c =2945900
1203.7
指数函数模型
(y =y 0e Αt +b )y 0=658Α=-3b =17548
31417双曲线模型y =a +b
c +t a =59100
b =-3362400
c =-
6009367从表1可以得出,运用二次曲线进行拟合得到的表达式的拟合误差比较小,因此本文选择水位高程h (m )与蓄水量x (104m 3)的模型为:
h =6.2945x 2-8608x +2.9459×106(x >683)
拟合曲线如图3所示.
图2 水位高程与蓄水量数据点示意图图3 水位高程与蓄水量拟合曲线
从图3可以看出,用二次曲线对水位高程为709～742米时的蓄水量拟合程度比较好,但是二次曲线只能拟合当水位高程大于683米的情况,而且当水位高程很大时,该曲线得到的蓄水量将会和实际差别很大,因为该曲线不能反映实际的地形条件.
2.2 基于D E M 的蓄水量计算[829]
如果考虑更精确地计算唐家山堰塞湖的蓄水量,那么就应该立足于唐家山堰塞湖的具体地形,因此,考虑基于D E M 来计算.
本文将数字地图用规则网格覆盖,然后进行网格抽取,由于规则网格将区域空间切为规则的格网单元,而每个单元网格的高程本文用该网格中心点的高程作为该网格内的平均高程,这样就可以得到该地图的高程矩阵.规则格网的高程矩阵,很容易进行计算机处理,例如可以进行等高线计算,坡度坡向以及自动提取流域地形等.如果数据量不大,可以运用手工83数 学 的 实 践 与 认 识
39卷
填表的方法采集的网格化高程数据,但是假如数据量过大,有必要研究D E M 文件的存储格式,进行程序化处理.
在后面的洪水演进中,
根据网格化高程数据计算地形的连通性是非常重要的,此部分将进行连通性分析.当选定一个边缘网格后,本文按照低于此网格高程数据h 的标准寻找其周围的网格并进行标注.当选定一个网格后,下一步要对其相邻的其他三个网格的高程进行判断,连通性分析如图4所示.
连通算法的表述为:选定一个需要计算联通的网格,其高程值记为h ,将其高程值与周围网格高程值进行比较,按照低于h 的标准判断其可能的方向,并对可能的方向数量进行累加,并将其判断过的网格进行标注,然后取其中一个方向继续判断,把方向数量减1,循环判断.但是在后来的判断过程中可能遇到方向指示已经标注过格网,这样为了不重复标记,对于此方向不进行计数.每次累加方向后判断方向数量是否为零,不为零时继续判断,直到方向数量为零的时候,连通范围也就判断出.
整体算法流程图可以表述如图5.在该流程图中,m 表示该网格周围被判断过的网格的数目,流程图比较重要的是:1)判断该网格与周围的连通关系;2)判断周围网格是否标记;3)继续进行连通性判断.图4 连通性分析算法示意图图5 连通性分析算法流程图
2.3 唐家山堰塞湖蓄水量计算
在提取连通性网格的时候,同时对于网格数量记数为p ,对于提取网格表示为H (i ,j ),这样堰塞湖湖面面积、蓄水量按如下公式进行计算
S y =p ×d x
2V y =d x 2×∑h (i ,j )
其中S y 表示堰塞湖湖面面积,d x 表示正方形格网边长,V y 表示堰塞湖蓄水量,h (i ,j )表示格网高程数值.
对于格网淹没的高度按如下公式计算
H y =h -h (i ,j )
其中H y 表示格网水面高度,h (i ,j )表示格网水位高程数值.D E M 计算得到的蓄水量与水
9316期曾向荣,等:唐家山堰塞湖泄洪问题研究
位高程曲线与拟合曲线对比如图6
.图6 拟合曲线与D E M 计算曲线对比图
前面已经提到,随着水位高程的增加,每增加单位高程的蓄水量的增加逐步变大,
而从上图中的两条曲线分析,曲线的斜率逐渐变大,符合实际情况.另一方面,由于二次曲线拟合是一条抛物线,而抛物线的拐点在水位高程为683米处,由于随着水位高程的降低,蓄水量增加显然不合理,因此在抛物线对称轴左侧不适宜采用拟合公式,而运用调洪演算得到的曲线较好地克服了这种拟合的缺点.
虽然拟合曲线在媒体报道点附近比调洪演算曲线更好地符合了实际数据值(3亿立方米),但是曲线拟合并不能很好地反应拟合点界限以外的数据.而运用调洪演算得到的曲线由于根据出入水量的关系进行演算,能够在较大的范围内很好地反应水位高程与蓄水量的关系.
2.4 调洪演算模型
1)水量平衡方程
在堰塞湖的蓄水过程中遵循水量平衡,即在某一时段内(?t =t n -t n -1),入湖流量与泄流量之差应等于该时段内堰塞湖蓄水量的变化值,由于蒸发水量很少,因此忽略堰塞湖蒸发水量,表示为:
Q n +Q n -1
2?t -q n +q n -1
2?t =W n -W n -1(1)
式中:Q n -1,q n -1分别为时段初入、出湖流量;Q n ,q n 分别为时段末入、出湖流量;W
n -1,W n 分别为时段初、末湖蓄水量;?t 为时段长.
由于不考虑泄洪,因此q n ,q n -1,W n =0,得:
Q n +Q n -12
?t =W n -W n -1(2) 2)数值解析法调洪计算原理
数值解析法采用堰塞湖水量平衡的微分方程.得:
d h d t =Q (t )F (h )
h (t 0)=h 0(起调水位)(3)
上式中,F (h )为堰塞湖水面面积,湖水面的面积由上面讨论的D E M 算法求得;Z (t )为水位高程;Q (t )为入湖流量;h 0,t 0分别为起调水位和起调时间.
04数 学 的 实 践 与 认 识39卷
图7 降水汇入堰塞湖流图
2.5 堰塞湖的蓄水特征
由于要考虑降雨量与水位高程的关系,入湖流量为堰塞湖上游流量、区间入流和直接降到湖面的水量,后两者都与
降雨量有关,如图7所示.
因此,得到
Q =f (R )+M
(4) 式中:Q 为入湖流量,R 为降雨量(mm ),f (R )为关于
R 的函数,M 为堰塞湖上游流量.则第t 日的入湖流量为
Q (t )=f (R (t ))+M
d t
(5)
假设相邻两日的堰塞湖上游流量相等,得
Q (t )=d f (R (t ))d t
(6) 联立(3)和(6)即可推算降雨和水位高程的关系
2.6 由多年平均降雨推求水位高程如果要在5月25日前对接下来17天唐家山堰塞湖水位高程进行预测,那就只能依赖历年的这些天的平均降雨情况了,因为天气预报很难准确预测17天的降雨情况.
根据历史日降雨量系列,分析计算多年平均情况下5月25日至6月12日降水过程,按起调水位723m (5月25日8点实测值),进行唐家山堰塞湖未来堰前水位预估.多年平均来水情况下堰塞湖堰前水位过程如图8.
图8 多年平均来水情况下唐家山堰塞湖堰前水位变化过程
根据上述分析计算结果可知,若唐家山堰塞湖以上流域未来一段时间降雨按多年平均考虑,到6月4日,唐家山堰塞湖蓄水将超过泄洪槽740m ,如果没挖泄洪槽,则6月11日后将会发生漫顶.
2.7 由天气预报降雨推求水位高程
根据天气预报结果,按其预报的上限、中间值和下限整理,得到唐家山堰塞湖以上流域5月25日至6月12日的降雨预报成果,根据多年平均日降水过程,按起调水位723m (5月25日8点实测值),进行唐家山堰塞湖未来堰前水位预估多年平均来水情况下堰塞湖堰前水位过程见图9.
1
416期曾向荣,等:唐家山堰塞湖泄洪问题研究
图9
不同预报来水情况下堰塞湖堰前水位变化过程对比图
根据上述分析结果,若唐家山堰塞湖以上流域未来一段时间降雨按天气预报结果考虑,到6月5日,唐家山堰塞湖蓄水将超过泄洪槽740m ,如果没挖泄洪槽,则6月12日后将会发生漫顶.预测当降雨量为天气预报降雨量中值的50%、80%和150%时,水位的结果如图10:
图10 水位高程预测曲线
从上面图表可以看出,随着降雨量的增加,水位高程增加的越快,天气预报降雨量的50%、80%、100%和150%时,堰塞湖湖水超过泄洪槽的时间依次为:6月9日、6月7日、6月5日和6月5日.当降雨量为天气预报降雨量的150%时,6月11日将发生漫顶,而为天气预报降雨量的50%、80%、100%时,6月12日以后水位高程达到警戒水位.
3 溃坝阶段的模型建立与求解
3.1 溃口流量
溃口流量利用溃口动态出流公式[10]计算
Q b =c v k s [3.1b i (h -h b )1.5+2.45z (h -h b )2.5]
式中,c v 为对行近流速的修正,b i 为溃口瞬时底宽(m ),h 为堰塞湖水位高程(m ),h b 为溃口底高程(m ),z 为溃口边坡(1:z ),k s 为尾水影响出流的淹没修正.
24数 学 的 实 践 与 认 识39卷
当(h l -h b ) (h -h b )>0.67时,k s =1.0-27.8h l -h b h -h b
-0.673,否则,k s =110.
上式中h l 为尾水位高程(m ).
c v 由下式计算
c v =1.0
+0.023Q 2b [B 2d
(h -h b m )2(h -h b )]式中B d 为堰塞湖坝前宽度(m ),h b m 为最终溃口底高程(m ).
如果溃口因管涌引起,那么可以应用孔流方程.Q b 可用以下孔流方程替代:Q b =4.8A Θ(h -h θ)
1 2其中
A Θ=2b i (h p -h b )
式中h p 为管涌中心高程线(m ),h θ=h p ,如果尾水位高程h l >h p ,则h θ=h l .当坝前水位下
降或管涌扩大到满足如下条件
h <3h p -2h b
则溃口流量由孔流转变为宽顶堰流.
3.2 溃口参数选择
早期由于对溃坝机理不甚了解,一般假设溃坝形式为瞬时全溃或顺势部分溃.对混凝土拱坝可以适用,但对土石坝和混凝土重力坝则不太适用.溃口形成如图11所示图11 溃口形成示意图[10]
对逐渐溃来说,主要是确定两个参数:最终溃口底宽b 和形状z (即溃口边坡).最终溃口底宽b 与溃口平均宽度b λ关系如下
b =b λ-z h d
式中h d 为溃口底部以上水高,一般以坝高近似代替.假设溃口底宽从一点开始,然后在溃决时间范围Σ内逐渐以线性或非线性速度增长,直到形成最终的溃口底宽b 和最终溃口底高程h b m .如果溃决时间Σ小于一分钟,那溃口底宽开始于限定的b 值,而不是某一点.与冲刷溃决相比,这更像是崩塌溃决.溃口底高程h b 是时间Σ的函数
h b =h d -(h d -h b m )t b ΣΘ0, 0 最终溃口底高程h b m 通常是水库或泄水河道段底部高程,但不是必须的.t b 为溃口形成开始算起的时间.Θ0为非线性程度参数,1ΦΘ0Φ4.通常假设为线性速度增长,因而Θ0=1.溃口瞬时底宽则由下式确定 b i =b t b ΣΘ0 , 0 溃坝模拟时,当库水位高程h 超过某一指定高程h f 时,溃口开始形成.对于土石坝,溃3 416期曾向荣,等:唐家山堰塞湖泄洪问题研究 口平均宽度b λ为0.5h Φb λΦ8h d ,如果是管涌引起的溃决,模拟管涌初始溃口Θ0Ε2是合适 的. 3.3 溃口冲刷速率分析溃口冲刷速率假定溃口冲刷形成的总时间为T ,总的溃坝深度为Z ,按照下述方法确定,可得到任意时间内的溃口发展数度V d .由以上分析,本文可以假定溃口发展模式符合以下曲线(如图12). 由溃口发展速度基本模式可以得到,当溃口发展时间t ΦT 2时,冲刷速率 V d =2(V m ax T )×t 当溃口发展时间T 2Φt ΦT 时,冲刷速率 V d =-2(V m ax T )×t +2V m ax 当溃口发展时间t > T ,冲刷速率为V d =0,可见V m ax =2Z m ax T . 上述各式中,V m ax 为溃口冲刷过程中的最大冲刷速率,t 为是时间,Z m ax 为给定可能冲刷溃口高度. 3.4 溃口冲刷稳定时间的确定 就目前的研究水平,对溃口流量与冲刷速度的关系仍然不是十分清楚,存在较大的不确定性,很难定量描述.对于如何确定溃口冲刷稳定时间,文章采用与宽顶堰公式和肖克列奇公式进行比较的方式进行确定. 将溃口冲刷速率公式带入溃口动态出流公式,在残留坝底为720m 的条件下,本文可以得到在溃口形成时间变化时,溃口流量的变化曲线,如图13. 图12 溃口形成速率示意图图13 溃口流量随时间变化图 将其于宽顶堰公式和肖克列奇公式得出的结果进行比较. 宽顶堰公式 Q m ax =m b 2g H +Αv 22g 3 2式中,m 为系数;Α为水头损失系数;H 为坝前水深(m );v 为流速(m s );g 为重力加速度 (m s 2). 肖克列奇公式Q m ax =8g 14bH 32 44数 学 的 实 践 与 认 识39卷 式中,B 为坝体宽度;b 为决口宽度(m );H 为坝前水深(m ). 通过计算,可以得到在相应溃口特征尺寸的最大流量分别为15500m 3 s 和13500m 3 s . 与上述曲线结果比较, 曲线中1小时左右溃口达到稳定的结果与这两个公式的计算结果最为接近.因此,可以大致确定在在堰塞湖蓄水漫顶1小时后溃口形态基本达到稳定.同时,这一小时内的最大溃口流量 14000m 3 s 可以作为下游风险评估的基本条件 .3.5 溃坝模型的灵敏度分析 按照上述讨论的非恒定流数学模型,由图3的水位蓄水量关系,假定不同的溃坝水位,不同的残留坝体高程,计算了各种情况下的溃坝最大流量及其溃口洪水过程,其中,各状态下的溃坝最大流量情况如图14～图16. 图14 747水位、不同残留坝高、不同溃口形成时间的溃口最大流量 图15 745水位、不同残留坝高、不同溃口 形成时间的溃口最大流量 据图14、图15、图16各种情况的计算结果可见,如果溃口冲刷至弱风化顶部720m 高程后停止冲刷,以堰顶高程747m 作为起始溃坝水位,得到溃坝最大流量约为16500m 3 s ;以堰顶高程745m 作为起始溃坝水位,得到溃坝最大流量约为14000m 3 s ;以堰顶高程742m 作为起始溃坝水位,得到溃坝最大流量约为11000m 3 s .又当固定水位而改变残留坝高时,溃口最大流量的变化同样非常明显.因此,可以看出由于残留坝高和水位的变化对溃口最大流量的影响还是相当剧烈的. 5416期曾向荣,等:唐家山堰塞湖泄洪问题研究 图16 742水位、不同残留坝高、不同溃口 形成时间的溃口最大流量 4 泄洪阶段的模型建立与求解 4.1 淹没区域计算模型 基于数字高程模型(D E M )[11]求取给定水位条件下的淹没区,应当区分两种情形.第一种:凡是高程值低于给定水位的点,皆计入淹没区;第二种则需考虑“流通”淹没的情形,即洪水只淹没它能流到的地方.例如对于环形山(一种中间低洼、四周环形隆起的地形),第一种淹没计算可能导致环形山内外都生成淹没区,而在第二种淹没中,外来的洪水若未及山顶,只能在山环外形成淹没区. 这两种情形都具有实际意义.第一种相当于整个地区大面积均匀降水的情形,所有低洼处都可能积水成灾;第二种相当于高发洪水向邻域泛滥,例如洪水决堤,或局部暴雨引起的暴涨洪水向四周扩散.通常,将这两种情形分别称为“无源淹没”和“有源淹没”.在有源淹没中,洪流不仅可能被阻于环形山外,也可能为任何高地所阻隔,一条流域的洪水一般不会越山而过在相邻流域形成淹没区.可以看出,唐家山堰塞湖溃坝洪水问题是有源淹没问题.本文采取种子蔓延算法处理唐家山溃坝洪水有源淹没问题. 种子蔓延算法与前面讨论的连通性分析算法相似,是一种基于种子空间特征的扩散探测算法,其核心思想是将给定的种子点作为一个对象,赋予特定的属性,在某一平面区域上沿4个(或8个)方向游动扩散,求取满足给定条件、符合数据采集分析精度、且具有连通关联分布的点的集合.利用种子蔓延算法计算淹没区,就是按给定水位条件,求取满足精度、连通性要求的点的集合,该集合给出的连续平面就是所要求算的淹没区范围.满足水位条件但与种子点不具备连通关联性的其他连续平面,将不能进入集合区内.淹没区计算的准确性在很大程度上依赖D E M 的分辨率.每一个象元都代表着地面上一个区域,且每一个象元都拥有自己的高程值和精度.依据象元的高程值精度,选择恰当的阈值作为判断象元归属的条件,该值称为种子蔓延探测分辨率.设空间数据精度为m ,探测分辨率为w ,经分析和计算实验,其关系式可以表达为: w =1n m +K (K 淹没区域不断扩大. 淹没区面积计算模型如图18, 有源淹没区的淹没范围计算比较复杂,所有像元与种子点必须连通,程序沿4个方向寻找满足条件的像元.有源淹没所获得的区域,必须与地形图上的行政区划层进行叠加计算,才能最后得出淹没区的实地位置.事先应将矢量行政区划层转换成网格数据. 图17 淹没区域计算算法流程图图18 淹没区面积计算模型 由于预先得知坝体发生溃塌造成堰塞湖内1 3的蓄水突然下泻,不考虑下渗等作用的影响,洪水淹没容积是可以计算得到的.由资料可以得到唐家山堰塞湖的最大蓄水容量为3亿立方米,发生溃坝时堰塞湖蓄水量按3亿立方米计算.又由于下流河道以及堰塞湖会储存部分洪水,约为2270万立方米,因此淹没容积为: W ′=30000 3-2270=7730(万立方米) 在种子点的起始位置和水面的水位高程h 的情况下,使用种子蔓延算法可以得到洪水淹没面积、容积和洪水淹没水深分布.因此,可以通过改变水面高程h ,使洪水淹没容积为7730万立方米,从而得到洪水淹没区域. 4.2 淹没面积以及人员撤离方案 计算机实现得到1 3溃坝影响到的地区有北川县城、通口镇、含增镇和青莲镇、绵阳县城(见图19),淹没总面积为1412.64万平方米,其中绵阳市的淹没面积为481.10万平方米,人员撤离18万左右. 4.3 洪水演进计算 要求得洪水到达时间,必须对洪水演进进行推算,由于圣维南方程组的拟线性及水流的不连续性,溃坝问题的分析求解难度非常大,因此采取近似方法和经验公式[12214]. 1)溃坝下游流量计算.坝址下游沿程各河段控制断面溃坝最大流量采用李斯特万公式计算: Q LM =W W Q M +L V K 7416期曾向荣,等:唐家山堰塞湖泄洪问题研究 图19 各地区淹没范围 式中,W 为堰塞湖蓄水量,Q LM 为距溃坝L (m )控制断面溃坝最大流量(m 3 s ),W 为堰塞湖蓄水量(m 3),Q M 为坝址最大流量(m 3 s ),L 为控制断面距堰塞湖坝址的距离(m ),V 为特大洪水的最大流速(山区取3～5m s ,丘陵区取2～ 3m s ,平原区取1～2m s ),K 为经验系数(山区取1.1-1.5,丘陵区取1.0,平原区取0.8～0.9). 考虑唐家山以下河道形态,取V K =5.0. 2)溃坝洪水传播时间和流量过程线.黄河水利委员会水利科学研究所根据实验得出了溃坝洪水传播时间及概化流量过程线. 洪水到达时间计算t 1=K 1L 1.75(10-h 0)1.3W 0.2H 0.350 式中,h 0为溃坝洪水到达前下游计算断面的平均水深,单位是m ,可以根据河道横切面粗略计算;K 1为系数,取为0.70×10-3;t 1是洪水起涨时间,单位是s . 最大流量到达时间计算 84数 学 的 实 践 与 认 识39卷 t 2=K 2L 1.4W 0.2H 0.50h 0.25M 式中,K 2为系数,取值0.8～1.2;h M 是最大流量时的平均水深,首先求出最大流量,再由水位流量关系曲线求得,单位是m ;t 2是最大流量到达时间,单位是s ,计算结果如表2所示. 表2 洪水演进计算结果 城镇名称 所属县市与溃坝距离(km )洪水到达时间(h )最大流量到达时间(h )北川县城北川县2.10.1210.364通口镇 北川县22.52.1254.225含增镇 江油市28.13.0475.613青莲镇 江油市48.75.2416.082绵阳县城绵阳市65.97.8568.510: [1] 唐家山堰塞湖抢险全记录[EB OL ].h ttp : www .new ssc .o rg ,200829228. [2] 地震灾区气象服务专报(第41期)[N ],200825225. [3] 四川省气象台.四川汶川地震灾区堰塞湖气象服务专报(第60期)[N ],20082622. [4] 中央气象台.四川省汶川天气预报[EB OL ].h ttp : www .sina .com .cn ,2008205217. [5] 四川汶川地震灾区气象服务专报(70期)[N ],2008206206. [6] 新华网直播车实时报道.唐家山堰塞湖抢险进展和人员撤离[EB OL ].www .xinhuanet .com .2008205231. [7] 长江委组织专家对唐家山堰塞滑坡体进行应急研究分析[EB OL ].h ttp : www .hw cc .com .cn ,200825221. [8] 陈凯.洪水淹没计算与评估系统[D ].硕士学位,2004.13218. [9] 刘仁义,刘南.基于G IS 的复杂地形洪水淹没区计算方法[J ].地理学报,2001,56(1):126. [10] 彭雪辉.风险分析在我国大坝安全上的应用[D ].硕士学位,2003.52254. [11] 吕杰堂,王治华,周成虎.西藏易贡滑坡堰塞湖的卫星遥感监测方法初探[J ].地球学报,2002,23:3632368. [12] 张振,张光科,刘超,杨中华.小型水库溃坝洪水演算探讨[J ].东北水利水电,2007,25(11):325. [13] 付磊,张洪明,姚激.水库调洪演算的数值解析法[J ].中国水利水电科学研究院《科研简报》特刊,2006,4(4):752 77. [14] 中国水科院.四川汶川抗震救灾[R ].2008,7. The Research of Tangj i ashan Barr ier Lake Sp illway ZEN G X iang 2rong , HAO Hong 2x ing , SU N Bo 2liang (D epartm en t of Info rm ati on System and M anagem en t ,N ati onal U n iversity of D efen se T echno logy ,Changsha 410073,Ch ina ) Abstract : T he paper studied the sp ill w ay law of T angjiashan barrier lake ,and did model analysis of its w ater sto rage ,dam b reak and flooding . Keywords : T angjiashan ;barrier lake ;sp ill w ay ;w ater sto rage ;dam b reak ;flooding 9 416期曾向荣,等:唐家山堰塞湖泄洪问题研究 堰塞湖及其处理 1.什么是堰塞湖 堰塞湖（barrier lake）是由山体滑坡堵塞河道后蓄水而形成的湖泊。堵塞河道的滑坡体, 称“堰塞体”, 实际上是一座天然水坝, 堰塞湖实际上是一座水库。 在堰塞湖形成的当时, 如果它对上下游构成威胁, 就要采取紧急措施进行处理。对于那些对上下游不构成威胁, 或者经处理后, 堰塞体还可加以利用, 那么, 还可将这些堰塞体或堰塞湖永久保留, 并改造成永久性的水坝和水库, 将害转化为利。 造成山体滑坡的原因主要有以下三种： （1）地震 由地震造成的滑坡体形成的堰塞体也称“地震坝”。我国重庆黔江上的小南海，就是在清代（1856年）因地震形成的一座堰塞湖, 坝长约100m, 坝高约30m。由于后来又人为开设了溢洪道，因此它能与人工水坝工程比美, 现已成为旅游胜地, 并享有“深山明珠”等美誉。对其成因，李四光曾有“冰窖”之说。世界上有众多的这类堰塞湖。 （2）火山熔岩流 由火山熔岩流形成的堰塞湖, 特称“熔岩堰塞湖”。 我国有多座著名的熔岩堰塞湖, 如： ①黑龙江省东南部的镜泊湖。是一座经5次火山爆发由熔岩流堵塞河道形成的。 ②黑龙江省的五大连池。由14座火山爆发形成, 在河道上形成了5个一串相邻的堰塞湖, 故名。 ③新疆的天山天池。天山天池距乌鲁木齐110km, 湖面呈半月形, 长3400m, 最宽处1500m, 最深105m, 有“天山明珠”的盛誉。 （3）特殊的地质构造 由特殊的地质构造造成的山体滑坡形成的堰塞体也称“山崩坝”。 我国陕西的天池就是由特殊的地质构造产生山体滑坡堵塞了太乙河形成的堰塞湖。天池也称太乙池。 诚如上述, 滑坡体可以堵塞河流形成水库；然而, 世界上还有完全相反的、也许是绝无仅有的一个史实：滑坡体竟然会使已有的水库消失——这就是意大利的瓦依昂（Vajont）坝库区的一次特大山体滑坡, 滑坡体将该水库填满, 原有水库不复存在, 该坝成了一座不蓄水的被废弃的坝。该滑坡体方量竟达2.4亿m3, 而当时的库容仅1.2亿m3。 还要指出的是, 就堵塞河道形成天然水坝来说, 除了上述的堰塞体外, 还有“冰坝”。由冰坝形成的水库称“冰坝湖”。 现在的研究发现, 在远古时期, 俄罗斯就存在有大型冰坝湖。 我国黄河每值冬春季节, 在多湾多滩的河段, 易发生冰凌插塞堆积现象, 严重时便形成了冰坝, 对上下游形成威胁。 2. 堰塞湖的处理措施 在发生山体滑坡形成堰塞湖的初期, 如果堰塞湖对上下游的人民生命财产构成威胁时, 就要进行紧急处理。处理措施主要有： （1）开挖溢洪道。 在堰塞湖岸开挖溢洪道（或明渠）是处理堰塞湖的常用方法, 如图1所示。溢洪道通常都是将湖水排泄到堰塞体下游的原河道, 也可考虑排泄到相邻河道或相邻河道的水库中（要具有这种地形条件, 同时在所泄湖水对相邻河道或相邻河道的水库不带来危害的情况下才予考虑）。 第39卷第16期2009年8月数学的实践与认识M A TH EM A T I CS I N PRA CT I CE AND TH EO R Y V o l 139 N o 116 A ugust,2009 唐家山堰塞湖泄洪的理论模型及数值模拟 王 莉1, 高毅欣2, 程佩佩1, 魏永生3 (1.西安交通大学理学院,陕西西安 710049) (2.西安交通大学电信学院,陕西西安 710049) (3.北京交通大学理学院,北京 100044) 摘要: 针对2008年全国研究生数学建模竞赛A 题,研究和解决汶川地震中唐家山堰塞湖的泄洪问题.建立了唐家山堰塞湖的蓄水量模型、溃坝模型、洪水演进模型和人员调度模型等理论模型,并给出了这些模型的精确数值模拟.模拟结果显示,提出的模型具有较高的精度,依据该模型提出的调度方案能够合理解决泄洪时的人员撤离问题,具有重要的参考意义. 关键词: 产汇流估计;圣维南方程组;动态模拟;网络流算法 1 问题分析 收稿日期:2009204230 2008年5月12日14:28在我国四川汶川地区发生了8.0级强烈地震,给人民生命财产和国民经济造成了极大的损失.地震引发的次生灾害也相当严重,截至5月22日,震灾区共发现堰塞湖33处,其中以唐家山堰塞湖尤为严重. 唐家山堰塞湖距离北川县城6公里,大坝坝顶高程750.2米,坝高82.8米,顺河长约220米,湖上游集雨面积3550平方公里.地震后每天新增至少500万立方米的水量,随着水位不断上涨,对下游地区造成了严重威胁. 本文工作旨在:1)建立唐家山堰塞湖以水位高程为自变量的蓄水量数学模型;2)在合理的假设下,建立堰塞湖蓄水漫顶后在水流作用下发生溃坝的数学模型,包含缺口宽度、深度、水流速度、水量、水位高程,时间等变量;3)根据数字地图给出坝体发生溃塌、造成堰塞湖内1 3的蓄水突然下泻时的洪水水流速度及淹没区域;4)根据地形和人口数据,提出当上述溃坝情况发生时的有效调度方案,从而极大减少损失. 2 蓄水量模型 设唐家山堰塞湖的水位高程和蓄水量分别为时间t 的函数H (t )和V (t ),并设t 时刻湖面面积为S t .根据流入堰塞湖中的水体积与堰塞湖新增水体积相等,建立微分方程:dV (t )=S t dH (t ) (2.1)根据[4]可知,唐家山堰塞湖的湖面面积S 和时间t 之间是分段线性函数,水位高程H (t )和时间t 之间的关系也是一个分段多项式函数,即: H (t )=a i (t -T i -1)m +H (T i -1),T i -1Φt 汶川地震中唐家山堰塞湖泄洪问题 摘要 本文研究的是唐家山地震次生灾害引发的堰塞湖问题，结合数字高程地图和新闻报道所提供的数据分别建立模型研究了唐家山堰塞湖形成之后湖水高程、蓄水量、溃坝情形、溃坝灾害等一系列问题。 针对第一问，首先对数字高程地图进行等高图像分析求解了堰塞湖不同高程水位（高程间隔为1米）对应的湖区面积，接着本文采用拟合法得到任意高程（710-750）的湖区面积，对高程积分建立了蓄水量体积与堰塞湖水位高程的离散化模型，然后建立了多元线性回归模型研究了北川天气预报3天降雨量与堰塞湖入库流量的关系，继而求解得到不同降雨量下每日堰塞湖水位高程。 针对第二问，首先研究泄洪过程和溃坝过程的区别与联系，从而分别建立正交多项式逼近和仿真模型得到溃坝时的溃口流量随时间变化的关系，继而分析求解得到溃坝时溃口宽度、深度、水位高程和水流速度随时间变化的关系。 针对第三问，综合数字高程地图和行政区域地图，在数字地图中查找地势相对较低区域,进而得到洪水下泄区域及被淹没区域。确定洪水的流速之后，利用数字地图计算了洪水到达各被淹没区域的时间，淹没范围，以及淹没之后的安全区域，并据此制订了初步的撤离方案。 关键字微元积分多元线性回归模型 2008年5月12日14：28在我国四川汶川地区发生了8.0级特大地震，给人民生命财产和国民经济造成了极大的损失。地震引发的次生灾害也相当严重，特别是地震造成的34处高悬于灾区人民头上的堰塞湖，对下游人民的生命财产和国家建设构成巨大威胁。加强对震后次生灾害规律的研究，为国家抗震救灾提供更有力的科学支撑是科技工作者义不容辞的责任。 唐家山堰塞湖是汶川大地震后山体滑坡后阻塞河道形成的最大堰塞湖，位于涧河上游距北川县城6公里处，是北川灾区面积最大、危险最大的堰塞湖，其堰塞体沿河流方向长约803米，横河最大宽约611米，顶部面积约为30万平方米，主要由石头和山坡风化土组成。由于唐家山堰塞湖集雨面积大、水位上涨快、地质结构差，溃坝的可能性极大，从最终的实际情况看，从坝顶溢出而溃坝的可能性比其它原因溃坝的可能性大得多。 经过专家分析，采取有效措施，最终完成了唐家山堰塞湖的成功泄洪。当时的科技工作者记录了大量的珍贵数据，新闻媒体也对唐家山堰塞湖进展情况进行了及时的报道，通过对这些数据的收集（由于数据来源不同，数据有些冲突，以新华社报道的相关数据为准），我们对堰塞湖及其泄洪规律进行了初步研究，完成以下工作： 1.建立唐家山堰塞湖以水位高程为自变量的蓄水量的数学模型，并以该地区天气预报的降雨情况的50%，80%，100%，150%为实际降雨量预计自5月25日起至6月12日堰塞湖水位每日上升的高度(不计及泄洪)。（由于问题的难度和实际情况的复杂性及安全方面的考虑，没有充分追求模型的精度，以下同）； 2.唐家山堰塞湖泄洪时科技人员记录下了大量宝贵的数据。我们在合理的假设下，利用这些数据建立堰塞湖蓄水漫顶后在水流作用下发生溃坝的数学模型，模型中包含缺口宽度、深度、水流速度、水量、水位高程，时间等变量。 3.根据数字地图，给出坝体发生溃塌造成堰塞湖内1/3的蓄水突然下泻时（实际上没有发生）的洪水水流速度及淹没区域（包括洪水到达各地的时间），并在此基础上考虑洪水淹没区域中人口密集区域的人员撤离方案。 4.根据我们所建立的数学模型分析当时所采取对策的正确性和改进的可能性。讨论应对地震后次生山地灾害 (不限堰塞湖) ，科技工作中应该设法解决的关键问题，并提出有关建议。 堰塞湖及其爆破处置方法研究 “堰塞湖”在“5.12”汶川大地震后引起公众的普遍关注。本文论述了采用爆破方式炸开缺口泄流一类堰塞湖应急排险的有效手段, 较详细地分析了堰塞湖滑坡坝治理中常用的工程爆破技术。包括具体方法的使用特点和工程处治典型实例。 新世纪新阶段，我国经济高速发展，针对未来可能发生的突发情况，随着我军的职能使命任务正在发生新的变化，由承担单纯的水电施工逐渐向应急救援领域转变。特别是近年来，由于自然环境恶化，灾害频发，导致重大的人员财产经济损失。作为水电兵，我们在抢险救灾方面拥有得天独厚的专业优势，在抢险救灾中可担任尖刀兵，为抢救人员生命和财产做出重要贡献。 从98特大洪灾到汶川特大地震灾害，我军在抢险救灾方面积累了丰富的经验。然而，在汶川地震中我军面临了一类新的灾害考验——地震堰塞湖。堰塞湖是由于河流被外来物质堵塞而形成的湖泊。常由山崩、地震、滑坡、石流、火山喷发的熔岩流和流动沙丘等造成。长发生在地震活跃，泥石流多发地带。地震堰塞湖的堤坝大多由滑落的松散石块组成，一旦决堤将会对下游民众造成毁灭性的破坏。在治理汶川地震形成的堰塞湖过程中，针对唐家山堰塞湖极难险情，我军联合各领域的专家解决了问题，但同时也反映出在处理这一类灾害上的技能欠缺，经验欠缺，排险措施单一等问题。本文主要是通过对各国处理地震堰塞湖的经验技术进行分析总结，为我军治理堰塞湖提供一定的技术支撑。 在汶川地震后，国内学者专家通过对唐家山等地震堰塞湖的应急处理积累了经验。在治理过程中通过多平台，多时相的卫星遥感影像，利用3S技术掌握了灾区堰塞湖的位置及范围，并利用先进的数学建模技术模拟灾区堰塞湖模型，为治理方案的产生提供了宝贵的参考依据。这些经验对我们应急抢险都有宝贵的借鉴价值，在处理未来的突发灾难时，我们可以借助先进的卫星通信，遥感成像，计算机模拟等技术为我军抢险救援提供技术支撑。 本文以处理堰塞湖方法技术为切入点，通过对爆破处置技术的原理、原则、具 清华大学-汶川地震中唐家山堰塞湖泄洪问题 本页仅作为文档页封面，使用时可以删除 This document is for reference only-rar21year.March 全国第五届研究生数学建模竞赛 题目汶川地震中唐家山堰塞湖泄洪问题研究 摘要： 本文首先根据库容量和水位高程的变化关系建立了一个库容量计算模型，得出库容量和水位高程符合指数关系，再根据降水量和水位高程变化的关系建立了降水量模型，分析了降水量对堰塞湖水位的影响，并给出了50%、80%、100%、150%各种降水情形下的水位变化。 然后利用新闻媒体搜集的数据建立了一个逐渐溃坝模型，该模型包括溃口计算模型，水流量计算模型和库容计算模型，包含了溃口宽度、深度、水流速度、水量、水位高程，时间等变量，并根据该模型计算了唐家山堰塞湖发生漫顶逐渐溃坝时的各种变量的数据。 根据河道内质量守恒定律和能量守恒定律，在假设河道分段逐渐变化的前提下推导得到了溃坝推演模型，并给出了模型的离散形式，在已知河道信息和溃坝处信息情况下的迭代推算方法，并将溃坝推演模型用来推测唐家山发生 1/3溃坝时的水流速度变化和水面高度变化，并提示了可能被洪水淹没的地区。 最后根据本文中的模型分析了当时采取的政策和方案，并提出了我们的建议。 参赛密码 参赛队号清华大学参赛队13 一、问题的提出 2008年5月12日14：28在我国四川汶川地区发生了级特大地震，给人民生命财产和国民经济造成了极大的损失。地震引发的次生灾害也相当严重，特别是地震造成的34处高悬于灾区人民头上的堰塞湖，对下游人民的生命财产和国家建设构成巨大威胁。加强对震后次生灾害规律的研究，为国家抗震救灾提供更有力的科学支撑是科技工作者义不容辞的责任。 唐家山堰塞湖是汶川大地震后山体滑坡后阻塞河道形成的最大堰塞湖，位于涧河上游距北川县城6公里处，是北川灾区面积最大、危险最大的堰塞湖，其堰塞体沿河流方向长约803米，横河最大宽约611米，顶部面积约为30万平方米，主要由石头和山坡风化土组成。由于唐家山堰塞湖集雨面积大、水位上涨快、地质结构差，溃坝的可能性极大，从最终的实际情况看，从坝顶溢出而溃坝的可能性比其它原因溃坝的可能性大得多。 经过专家分析，采取有效措施，最终完成了唐家山堰塞湖的成功泄洪。当时的科技工作者记录了大量的珍贵数据，新闻媒体也对唐家山堰塞湖进展情况进行了及时的报道，通过对这些数据的收集（由于数据来源不同，数据有些冲突，以新华社报道的相关数据为准），我们对堰塞湖及其泄洪规律进行了初步研究，完成以下工作： 1.建立唐家山堰塞湖以水位高程为自变量的蓄水量的数学模型，并以该地区天气预报的降雨情况的50%，80%，100%，150%为实际降雨量预计自5月25日起至6月12日堰塞湖水位每日上升的高度(不计及泄洪)。（由于问题的难度和实际情况的复杂性及安全方面的考虑，没有充分追求模型的精度，以下同）； 2.唐家山堰塞湖泄洪时科技人员记录下了大量宝贵的数据。我们在合理的假设下，利用这些数据建立堰塞湖蓄水漫顶后在水流作用下发生溃坝的数学模型，模型中包含缺口宽度、深度、水流速度、水量、水位高程，时间等变量。 3.根据数字地图，给出坝体发生溃塌造成堰塞湖内1/3的蓄水突然下泻时（实际上没有发生）的洪水水流速度及淹没区域（包括洪水到达各地的时间），并在此基础上考虑洪水淹没区域中人口密集区域的人员撤离方案。 4.根据我们所建立的数学模型分析当时所采取对策的正确性和改进的可能性。讨论应对地震后次生山地灾害 (不限堰塞湖) ，科技工作中应该设法解决的关键问题，并提出有关建议。堰塞湖及其处理
唐家山堰塞湖泄洪的理论模型及数值模拟
汶川地震中唐家山堰塞湖泄洪问题
堰塞湖的处置方法
清华大学-汶川地震中唐家山堰塞湖泄洪问题
相关主题
文本预览