汶川地震中唐家山堰塞湖泄洪问题

汶川地震中唐家山堰塞湖泄洪问题摘要本文研究的是唐家山地震次生灾害引发的堰塞湖问题，结合数字高程地图和新闻报道所提供的数据分别建立模型研究了唐家山堰塞湖形成之后湖水高程、蓄水量、溃坝情形、溃坝灾害等一系列问题。针对第一问，首先对数字高程地图进行等高图像分析求解了堰塞湖不同高程水位（高程间隔为1米）对应的湖区面积，接着本文采用拟合法得到任意高程（710-750）的湖区面积，对高程积分建立了蓄水量体积与堰塞湖水位高程的离散化模型，然后建立了多元线性回归模型研究了北川天气预报3天降雨量与堰塞湖入库流量的关系，继而求解得到不同降雨量下每日堰塞湖水位高程。针对第二问，首先研究泄洪过程和溃坝过程的区别与联系，从而分别建立正交多项式逼近和仿真模型得到溃坝时的溃口流量随时间变化的关系，继而分析求解得到溃坝时溃口宽度、深度、水位高程和水流速度随时间变化的关系。针对第三问，综合数字高程地图和行政区域地图，在数字地图中查找地势相对较低区域,进而得到洪水下泄区域及被淹没区域。确定洪水的流速之后，利用数字地图计算了洪水到达各被淹没区域的时间，淹没范围，以及淹没之后的安全区域，并据此制订了初步的撤离方案。一、 问题重述2008年5月12日14：28在我国四川汶川地区发生了8.0级特大地震，给人民生命财产和国民经济造成了极大的损失。地震引发的次生灾害也相当严重，特别是地震造成的34处高悬于灾区人民头上的堰塞湖，对下游人民的生命财产和国家建设构成巨大威胁。加强对震后次生灾害规律的研究，为国家抗震救灾提供更有力的科学支撑是科技工作者义不容辞的责任。唐家山堰塞湖是汶川大地震后山体滑坡后阻塞河道形成的最大堰塞湖，位于涧河上游距北川县城6公里处，是北川灾区面积最大、危险最大的堰塞湖，其堰塞体沿河流方向长约803米，横河最大宽约611米，顶部面积约为30万平方米，主要由石头和山坡风化土组成。由于唐家山堰塞湖集雨面积大、水位上涨快、地质结构差，溃坝的可能性极大，从最终的实际情况看，从坝顶溢出而溃坝的可能性比其它原因溃坝的可能性大得多。经过专家分析，采取有效措施，最终完成了唐家山堰塞湖的成功泄洪。当时的科技工作者记录了大量的珍贵数据，新闻媒体也对唐家山堰塞湖进展情况进行了及时的报道，通过对这些数据的收集（由于数据来源不同，数据有些冲突，以新华社报道的相关数据为准），我们对堰塞湖及其泄洪规律进行了初步研究，完成以下工作：建立唐家山堰塞湖以水位高程为自变量的蓄水量的数学模型，并以该地区天气预报的降雨情况的50%，80%，100%，150%为实际降雨量预计自5月25日起至6月12日堰塞湖水位每日上升的高度（不计及泄洪）。（由于问题的难度和实际情况的复杂性及安全方面的考虑，没有充分追求模型的精度，以下同）；唐家山堰塞湖泄洪时科技人员记录下了大量宝贵的数据。我们在合理的假设下，利用这些数据建立堰塞湖蓄水漫顶后在水流作用下发生溃坝的数学模型，模型中包含缺口宽度、深度、水流速度、水量、水位高程，时间等变量。根据数字地图，给出坝体发生溃塌造成堰塞湖内1/3的蓄水突然下泻时（实际上没有发生）的洪水水流速度及淹没区域（包括洪水到达各地的时间），并在此基础上考虑洪水淹没区域中人口密集区域的人员撤离方案。根据我们所建立的数学模型分析当时所采取对策的正确性和改进的可能性。讨论应对地震后次生山地灾害（不限堰塞湖），科技工作中应该设法解决的关键问题，并提出有关建议。二、 问题假设1、 降雨量预报数据绝大部分为区间值，而且一天预测多次，我们取其平均值作为一天的降雨量。不计及泄洪，水位上升高度除与降雨量有关，与上游流量亦有关联。考虑用神经网络模型反映前两天的降雨量与当日入湖流量增量的关系。进而利用堰塞湖该地区的降雨量预测湖中水位上升高度；2、 假设单位长度湖面可以看做为矩形，当长度很小的时候鉴于微元的思想是3、 溃坝过程与泄洪过程内在机理是相同的，泄洪相当于一次小型溃坝；4、 溃坝水位高程750米对应的堰塞体构造与泄洪743米对应的堰塞体构造是一样的，因而溃口处沙石泥土等内部结构一样。三、 符号说明和名词解释V：堰塞湖内蓄水量，即总库容，单位：亿立方米h(t)：坝前水位高程，单位：米h:堰塞湖底部高程，常数二667.4米b0R(t):堰塞湖每天的新增水量，单位:亿立方米Q(t):t时刻的单位入湖流量，单位:立方米/秒INQ(t):t时刻的单位泄流量，单位:立方米/秒OUT四、 问题分析研究不同高程水位对应的堰塞湖蓄水量，即研究不规则体的体积问题，可以考虑引入面积对高度的微元进行积分这一思想进行求解。新闻报道中的实际数值可以作为求解时的重要参考和校正点。在研究降雨量和堰塞湖水位的关系之时，认为将这一问题转化为首先研究降雨量和库容变化量之间的关系更为直观，考虑到两者关系的复杂性，并不一定需要建立明确的函数关系式。现实中唐家山并没有发生溃坝，无法获得溃坝的直接数据。这里也十分缺乏利用成型模型求解溃坝问题的大量参数。但是唐家山的泄洪的规模也是非常大的，可以借助研究唐家山的泄洪过程分析这一泄洪和溃坝的区别与联系，继而通过拟合、仿真等手段得出唐家山的溃坝过程。结合数字地图可以求得不同位置的海拔高度，根据水顺势而流的原理就可以推算洪水的流经区域。结合泻出的洪水量就可以结合问题1的思想计算求解淹没范围。

五、 模型准备5.1.1模型的分析与准备初始数据提取与分析唐家山堰塞湖是由地震活动引起山崩滑坡体堵截河谷或河床后贮水而形成的湖泊，湖泊总蓄水量与湖泊水位高程相关，研究附件提供的数据（见表一）表一唐家山堰塞湖水位（米）与蓄水量（亿立方米）水位蓄水量水位蓄水量水位蓄水量713.410.834717.511738.52.159713.420.835718.061.02738.712.175713.490.837719.481.12738.852.185713.50.838720.251.139738.932.19713.540.8397331.8739.372.222713.570.84733.421.834739.592.237713.640.843735.531.957740.512.328713.680.844736.632.008740.552.307713.790.849736.782.055741.122.351714.010.857737.532.094741.662.389714.130.861737.562.071742.582.457714.510.876737.782.116742.82.469注：数据特点：数据分布不均匀，主要集中在（710□720）和（730□740）两个水位段；数据量较少，整个水位高程对应的数据不完善；随着水位高程增加堰塞湖蓄水量的变化率趋势总体上越来越大，但增长规律不够稳定，有些高程大的水位蓄水量增长反而慢，数据可能有误差。数字高程地图图像处理湖泊总蓄水量由不同水位高程的剖面蓄水量组成，且不同水位高程剖面的蓄水量随着水位高程增加有增大的趋势，图5.1，图5.2为数字高程地图处理后所得730米、750米水位高程对应的横向剖面图：

相对面积的获取可以从附件所给的唐家山数字地图上得到，具体步骤如下：Stepl:利用所给的3DEM软件将tjs.dem文件载入。Step2:利用唐家山的经纬度与高度，从709m开始至750m,以1m为间隔设定某一高度值，例如720m，为高度分界值后，将原数字地图分为两种明显的不同颜色。Step3:保存新得到的数字地图。

Step4:使用Matlab载入新得的数字地图，运用不同颜色的像素点的值不同的原理，提取出此时堰塞湖所在地的地图，并求得湖面的图像面积。Step5:再次设定另一高度值，重复进行Step2—Step4直到高度值到达750m。其流程图如下图所示：图三堰塞湖相对面积处理流程图通过对数字高程地图不同高度进行分色处理结合唐家山的经纬度，利用Matlab获取图像就可获得不同高程下唐家山堰塞湖湖面的图像面积（1个像素为单位1,程序见附录一）如表二所示：表二唐家山不同水位高程（米）对应的横向剖面相对面积水位高程相对面积水位高程相对面积水位高程相对面积水位高程相对面积71111596721133587311529174117239712117667221353673215471742174237131192572313728733156577431759371412093724139197341586174417794715122837251411573516049745180107161245872614323736162327461820171712629727144957371644174718405718128047281468673816625748185917191300972914882739168127491879672013183730150837401702975018992六、 模型建立与求解问题一总蓄水量与坝前水位高程的数学模型和降水模型6.1.1任意高程相对面积的求解本文采用数据拟合的方法，利用mathematica6.0得到相对面积的连续函数:S（h）二0.3145x2-268.52x+141624为了直观的查看拟合效果，本文将散点图和拟合函数图画在一起，如图：

通过图四，本文可知拟合效果很好，可以以此拟合函数代表s(H)通过图四，本文可知拟合效果很好，可以以此拟合函数代表s(H)。6.1.2一般蓄水量和比值K的求解V(H)=kfS(H)DH,K为实际大小与相对大小的比值h0根据表5.1可知部分AV，接着本文利用微元积分思想得到对应高程段的AV'，具体如下:TOC\o"1-5"\h\zAV AVk= =—AV' hfS(h)dhh0V(H)二V(H)二5.08x10-3h0按上述步骤求解，得到实际大小与相对大小的比例系数K值，结果见表三。表三比例系数K值0.055470.003390.00057442.16E-00.00035241550.00260.0006856270.002260.003467.5742E-04450.00100.003980.0006585880.055450.00036120.000140.004750.0008460.0022953510.0004050830.003460.000910.00027000.000210.000420.001576567460.0012478510.001530.00017380.010710.00054avrage(k)=0.00508770.006167560.0018561综上得到唐家山堰塞湖以水位高程为自变量的蓄水量的数学模型:堰塞湖水位上升模型6.2.1堰塞湖水位上升模型分析堰塞湖每日上升的高度，与湖水上游的降雨量、湖水上游河流流入湖水的流量、湖水地区的降雨量、湖水的泄流量、湖水的蒸发等条件有关，为了降低问题的复杂性，简化运算，可忽略湖水的蒸发等对湖水高度影响较小的条件；题目中要求不考虑湖水泄流，故湖水泄流对湖面上升高度的影响也可不考虑；湖水上游的降雨量、湖水上游河流流入湖水的流量可统一作为上游对湖水的影响考虑。所以最终考虑其中的两个主要因素：湖水地区降雨量、上游流水的影响。数据来源及处理唐家山堰塞湖位于四川省北川县，收集北川地区的天气情况(降雨量)即可了解堰塞湖区域的天气情况(降雨量)，根据研究人员提供的资料[2][3][4][5][6][7][8]分析得到5月25日到6月12日的天气情况，如表四所示：表四北川降水情况统计表日期5月255月265月275月285月295月305月31降水量(mm)8-1810-2501-22-500期望值(mm)1317.501.53.500日期6月16月26月36月46月56月66月7降水量(mm)00-10-1012-285-111-3期望值(mm)00.50.502082日期6月86月96月106月116月12降水量(mm)0000-10-1期望值0000.50.5(mm)其中天气预报的降雨量5为如表2所示的范围值，我们假设实际降水量满足均匀分布的规律，则预测值的期望值为最低降雨量与最高降雨量的平均值产，实际的降雨量5实际mx5(m分别取50%、80%、100%、150%)。根据研究人员提供的资料知道5J实际4日唐家山堰塞湖湖水高程为h=721.07m〔i]o5月25日至6月12日期间唐家山堰塞湖上流日均入库流量0导致湖水日均上升的高度Ah二0.9m。6.20.3模型建立及求解设唐家山堰塞湖湖水高程用h(t)表示，则后一天的湖水高程用h(t+1)表示；当天湖水高程变化量为Ah；设湖水集水面积为S(t)，单位：m2；湖水地区日降雨量5实际，单位：mm；湖水上流日均入库高度：Ah，单位：m； 5月24日唐家山堰塞湖湖水高程为h=721.07m。由于在水位高程增加量较小的情况下，湖水蓄水量的增加量约等于此时湖水横截面积乘以水位高程的增加量，故可列下列方程：5=h(t+1)S(h)dh实际h(t)

h（0）二h,Ah（t+1）=h（t+1）—h（t）,t代表第几天0后面一天湖水的水位高程是前一天湖水水位高程加上由于降雨造成的水位高程的增加量和湖水上流入库引起的水位增加量，则：h（t+1）二h（t）+Ah（t+1）+Ah0根据初始条件和方程解m分别为50%、80%、100%、150%时，每天的湖水水位升高高度和每天的湖水水位高程，具体数据如表五所示：表五湖水水位升高高度和湖水水位高程数据表日期5月25日5月26日5月27日5月28日5月29日增加高度m（150%）3.83834.72210.91.2143水位高程（150%）721.07724.91729.63730.53731.74增加高度（100%）2.85893.46920.91.1128水位高程（100%）721.07723.93727.4728.3729.41增加高度（80%）2.46712.96240.91.0713水位高程（80%）721.07723.54726.5727.4728.47增加高度（50%）1.87942.19570.91.0082水位高程（50%）721.07722.95725.15726.05727.05日期5月30日5月31日6月1日6月2日6月3日增加高度m（150%）1.62770.90.90.91.0015水位高程（150%）733.37734.27735.17736.07737.07增加高度（100%）1.39250.90.90.90.96858水位高程（100%）730.8731.7732.6733.5734.47增加高度（80%）1.29660.90.90.90.9552水位高程（80%）729.77730.67731.57732.47733.42增加高度（50%）1.15050.90.90.90.93483水位高程（50%）728.2729.1730730.9731.84日期6月4日6月5日6月6日6月7日6月8日增加高度m（150%）1.0010.94.90492.47381.2908水位高程（150%）738.07738.97743.88746.35747.64增加高度（100%）0.968210.93.60231.96371.164水位高程（100%）735.44736.34739.94741.91743.07增加高度（80%）0.954890.93.07351.75671.1126水位高程（80%）734.38735.28738.35740.11741.22增加高度（50%）0.934630.92.27031.44151.0344水位高程（50%）732.77733.67735.94737.39738.42日期6月9日6月10日6月11日6月12日6月13日增加高度m（150%）0.90.90.90.996960.99683水位高程（150%）748.54749.44750.34751.34752.34增加高度（100%）0.90.90.90.965230.96508水位高程（100%）743.97744.87745.77746.74747.7增加高度（80%）0.90.90.90.952450.95231

水位高程(80%)742.12743.02743.92744.87745.83增加高度(50%)0.90.90.90.933090.93298水位高程(50%)739.32740.22741.12742.05742.99为了更为直观地观察湖水水位升高高度，可拟出m=50%、80%、100%、150%分别情况下水位高度变化的数据走势图，如图五所示：65465432101 4 7 10 13 16 19—增加高度m(150%)—增加高度(100%)4增加高度(80%)—增加高度(50%)图五水位高度变化走势图通过图示可知，5月25、26日，6月5日雨量较大，湖水水位高程增加明显；雨量较小时，水位有较小的增加；不下雨的情况下，水位仍有提高，所以在5月25日至6月12日期间，唐家山堰塞湖形势严峻，有较大的溃坝危险。4．模型结果分析关于降雨量计算的不准确有可能导致模型建立的不正确，为了对模型结果进行分析，需要将此模型的水位高程与材料中的实际值进行对比，选择出最优的k值，作为模型中实际降雨量的比例常数。通过计算，水位高程的走势图如图六所示：图六水位高程走势图将表3与表1在5月25日至6月12日19天内的水位高程值进行对比，当k=80%时，预计水位高程与实际水位高程相对较接近。6月10日时，预计水位高程为743.02，实际水位高程为743.13，唐家山堰塞湖水位已经超过了橙色警戒线，形势相当严峻。6.3问题二的模型建立与算法设计

逐步溃坝模型这里采用美国国家气象局编制的溃坝洪水预报DAMBRK模型。该模型由三部分组成：（1） 大坝溃口形态描述。用于确定大坝溃口形态随时间的变化，包括溃口底宽、溃口顶宽、溃口边坡及溃决历时；（2） 水库下泄流量的计算；（3） 溃口下泄流量向下游的演进；溃口是大坝失事时形成的缺口，溃口的形态主要与坝型和筑坝材料有关。目前对于实际溃坝机理仍不是很清楚，因此溃口形态主要通过近似假定来确定。考虑到模型的直观性、通用性和适应性，一般假定溃口底宽从一点开始，在溃决历时内，按线性比率扩大，直至形成最终底宽。若溃决历时小于10分钟，则溃口底部不是从一点开始，而是由冲蚀直接形成最终底宽。溃口形态描述主要由四个参数确定：溃决历时、溃口最终底宽、溃口边坡、溃口底部高程。由第一个参数可以确定大坝溃决是瞬溃还是渐溃。由后面三个参数可以确定溃口断面形态为矩形、三角形或梯形及局部溃或全溃。6.3.1唐家山堰塞湖基本情况根据现场查勘资料，唐家山堰塞湖基本情况如下：唐家山堰塞湖坝体位于北川县城上游6km的通口河上，下距苦竹坝约1km,集雨面积3550km2，堰塞湖总容积约3.15亿m3。坝址区通口河为不对称的V型谷，右岸较陡，坡度60°左右；左岸坡较缓，坡度约为30°。坝区两岸基岩为寒武系下统清平组薄层硅质岩、砂岩、泥灰岩、泥岩组成。岩层软硬相间，产状N60°E/NWZ60。。左岸为逆向坡，右岸为顺向坡。基岩裂隙较发育，岩体较破碎，强风化带厚度5〜10m，其下为弱风化，高程大体为720m。两岸坡分布有残坡积的碎石土，厚度0〜15m。碎石土由粉质壤土、岩屑和块石组成，其中粉质壤土占60%左右，岩屑占30%〜35%，块石占5%〜10%。6.3.2溃口流量公式的推导土石坝、拱坝、混凝土重力坝的溃决形式以及溃决过程是不一样的。拱坝的溃决过程较快，几乎在瞬间发生，故模型按瞬溃考虑；混凝土重力坝则可能有一个或者多个坝段溃决，因此其溃决时间稍长，约为几分钟；土石坝不会全部溃决，也不会瞬间溃决，溃决时间相对最长，约从几分钟到几个小时，溃决的快慢，取决于土石坝的土力学参数。这里由唐家山堰塞湖坝体的成分可知，堰塞湖坝体为土石坝。在溃决过程中和溃决后的溃口的形状不同，其流量和水位的过程也不同，因此有必要分成两个过程分别进行研究。图7为溃口大体形状。Wb・・・・溃决过程中：在发生地震时间T内溃决・・・・溃决过程中：在发生地震时间T内溃决口的形状形成并且达到瓠我们假设坝的高度和溃口宽度随着时间呈线性变化，即h=kt,w=kt，又v=cv'2gh，由此便TOC\o"1-5"\h\zd 2up1 v得到了在发生地震溃决口的形状的时间t内的流量与水位的关系：Q二Axv二kkt2xCJ2gH2 12 vw全部溃决后：根据流体力学中的基本知识，在全部溃决后,其缺口面积基本保持不变，由此建立贝努力流体方程为：p v2 p v2O+4+H=厶+1

Y2g Y2g其中p:水库的液面上的大气压；0p:大坝缺口处的侧面液面大气压强；1v:水库的水流速；0v:大坝缺口处的水流速；1H:水库与大坝缺口处的高度差；因为水库的水可以看作静止的，所以v二0，水库的大气压与大坝缺口处的0大气压相等，p=p01将上述条件代入式(5.12)可以得到洪水流出速度为：v^/2gH由方程(5.13)推出的洪水流出速度v,是在理想情况下得到的。在现实生活1中，由于各种因素影响，我们必须在公式中加入系数修整，令C为行进流速的v修正，通常为0.96〜0.99之间，即v=C。1八根据溃口处洪水流量与洪水流出速度及溃口宽度之间的物理关系，建立微分方程模型：AQ二AWAvb于是得到时刻t的溃决口流量为：H(t) H(t) Q(t)二JW(H)v(H)dH=JW(H)C弋2gHdHvHHbb模型的简化处理：由于土石坝发生溃决其溃口有不同的形状，可能是三角形、梯形，也可能是矩形(我们对溃口进行了简化，限定其为几种基本的几何图形)，以下通过不同的缺口模型计算其流量。为了简化计算，我们用参数K来表示溃口的上底宽与堤口底宽的比值，如上图所示，根据上述流量的物理关系模型及图形的几何关系，可以设计以下理想模型：W(H)二W+(W-W)x(H-H)/Hdownupdwdd令W二kxW,W(H)二W+(W-W)x(H-H)/H，可得流量updowndownupdwddQ(t)(WC\.：2gH+(k-1)WCp2gH(H-H)/H)dHdownv downv b bHb当时k*,缺口模型为一个三角形，如图所示；当1<k5时，缺口模型为一个梯形，如图所示；当k二1时，缺口模型为一个矩形，如图所示。引入具体的流体力学中的参数：K:尾水影响出流的淹没修正；b:溃口的瞬时底宽；Z:溃口边坡比。则si积分(5.17)变为:Q(t)二KC(3.1b(H-H)1.5+2.45z(H-H)2.5)svi b b5.2.2逐步溃坝模型的求解6.4第三问的模型建立与求解假设某段河道的出水流速为Qout，如水流速为Qin，河道宽度为B，河道长度为S，水的密度p，引起的水面高度变化为Ah，根据质量守恒定律，(Qout-Qin)tp二pjS(h)dh假设某质量为m的水从海拔H的地方流到海拔H'，相应水面高度由h变化为h'，速度由v变化为v'。损失的重力势能为:mg(H+h-H'-h)22增加的动能为:0.5m(v'2-v2)损失的摩擦能量为:卩mvs其中，v二0.5(v+v'),卩是摩擦能量损失系数，表示单位质量的水以