清华大学-汶川地震中唐家山堰塞湖泄洪问题

全全国国第第五五届届研研究究生生数数学学建建模模竞竞赛赛题目汶川地震中唐家山堰塞湖泄洪问题研究摘要：本文首先根据库容量和水位高程的变化关系建立了一个库容量计算模型，得出库容量和水位高程符合指数关系，再根据降水量和水位高程变化的关系建立了降水量模型，分析了降水量对堰塞湖水位的影响，并给出了50%、80%、100%、150%各种降水情形下的水位变化。然后利用新闻媒体搜集的数据建立了一个逐渐溃坝模型，该模型包括溃口计算模型，水流量计算模型和库容计算模型，包含了溃口宽度、深度、水流速度、水量、水位高程，时间等变量，并根据该模型计算了唐家山堰塞湖发生漫顶逐渐溃坝时的各种变量的数据。根据河道内质量守恒定律和能量守恒定律，在假设河道分段逐渐变化的前提下推导得到了溃坝推演模型，并给出了模型的离散形式，在已知河道信息和溃坝处信息情况下的迭代推算方法，并将溃坝推演模型用来推测唐家山发生1/3溃坝时的水流速度变化和水面高度变化，并提示了可能被洪水淹没的地区。最后根据本文中的模型分析了当时采取的政策和方案，并提出了我们的建议。参赛队号清华大学参赛队13参赛密码（由组委会填写）2一、问题的提出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.根据我们所建立的数学模型分析当时所采取对策的正确性和改进的可能性。讨论应对地震后次生山地灾害(不限堰塞湖)，科技工作中应该设法解决的关键问题，并提出有关建议。3二、符号说明W：堰塞湖内蓄水量，即总库容，单位：亿立方米()Ht：坝前水位高程，单位：米0bH：堰塞湖底部高程，常数667.4米()Lht：堰塞湖内水深，单位：米()Rt：堰塞湖每天的新增水量，单位：亿立方米()Jt：第t天的降雨量，单位：毫米()bt：泄流槽的宽度，单位：米()INQt：t时刻的单位入湖流量，单位：立方米/秒()OUTQt：t时刻的单位泄流量，单位：立方米/秒溃口中的水b0HH(t)坝体地基wh(t)b(t)水面坝顶TOPHh(t)4三、模型的建立与求解1．总蓄水量与坝前水位高程的数学模型1.1一般模型一般情况下，在截面积规则的情况下，蓄水量可用水深的二次方或三次方来进行描述，但由于唐家山堰塞湖湖体结构复杂，其蓄水量不能用其水深的二次方或三次方进行简单描述。根据报道中搜集的所有数据，可以基本确定总蓄水量W与水深Lh存在指数关系，假定：nLWah(1.1)其中，a与n为库容特性系数。根据新华社报道“唐家山堰塞湖坝顶高程750.2米，坝高82.8米”，可以假设堰塞湖的底部高程0750.282.8667.4bH米，从而可以得出湖内水深与坝前水位高程H的关系为：0LbhHH（1.2）由方程（1.1）、（1.2）联立可以得出总蓄水量与坝前水位高程的数学模型为：0()nbWaHH（1.3）将方程（1.3）两边取对数，方程变换为：0ln()ln()ln()bWanHH（1.4）根据所给的总蓄水量W与相应坝前水位高程H的实际数据，计算相应的ln()W与0ln()bHH，对方程（1.4）进行线性拟合，可以得出参数a与n的最优解。3.73.83.944.14.24.34.44.5-0.4-0.200.20.40.60.81log(H-667.4)log(W)数据采样图1蓄水量与水位高程采样点数据（取对数处理后）53.73.83.944.14.24.34.44.5-0.6-0.4-0.200.20.40.60.811.21.4log(H-667.4)log(W)数据拟合y=2.1894*x-8.5638data1linear图2采样点数据的数据拟合拟合结果：2.1894nln()8.5638a=41.908910a残差=0.07455故可以得出最终的总蓄水量W与坝前水位高程H的数学模型为：42.18941.908910(667.4)WH（1.5）1.2降雨模型为了预测降雨给堰塞湖水位带来的影响，首先根据实际的降雨量建立降雨-水位模型。经过分析，我们认为，第n天（假设5月25日是第1天）之后堰塞湖每日的新增水体主要来自于以下几个方面：(1)前24小时直接落在堰塞湖面的降雨()aRt；(2)周围山体的降雨经过一定的延时后汇入堰塞湖而形成的水体()bRt；(3)上游地区的积雨经过一定的延时后流入堰塞湖的水体()cRt；(4)河流上游正常流入湖体的水体()dRt；根据附件1提供的数据，我们采用了每天早上8点的堰塞湖水位数据，这样(1)、(2)部分主要考虑的是前一天晚上和前一天白天的降雨带来的影响，而假设(3)部分即上游的积雨要经过一天的时间汇入堰塞湖，这样，(4)部分主要考虑的是前48小时到前24小时的降雨带来的影响。附件一给出的降雨量值，有的分为白天和晚上的降雨量，这里均转换成了一天的降雨量。每部分的模型分别如下：(1)对于()aRt，假设堰塞湖面积为1S，前24小时的降雨量为1(1)Jt，那么可得以下关系式：11()(1)aRtSJt(1.6)6这里，由于堰塞湖处于北川境内，所以1()Jt用北川的降雨量来计算。(2)对于()bRt，假设周围山体集雨面积为2S，前24小时的降雨量为2(1)Jt，最终汇入堰塞湖的雨水占山体总集雨的比例为2，那么：222()(1)bRtSJt(1.7)这里，山体和堰塞湖同处一个地区，所以2()Jt仍用北川的降雨量来计算。而且由下面的卫星照片看出山体的面积远远大于堰塞湖本身的面积，所以这一部分的雨量对于堰塞湖水位上涨的贡献量较大。图3唐家山卫星图片(3)对于()cRt，由下面的行政图可以看出，茂县处于唐家山堰塞湖的上游。设其面积为3S，一天前的降雨量为3(2)Jt（这里假设上游的雨水经过了一天的流程到达唐家山，所以用的是前48至前24小时的降雨量），最终汇入唐家山的雨水占茂县总集雨的比例为3，得到：333()(2)cRtSJt(1.8)图4唐家山地区行政图(4)对于()dRt，这里定义的是上游河流的正常流入量（即不考虑降雨因素），7所以假设它是一个不随时间变化的常数。即()dRtC(1.9)这样，堰塞湖每天的新增水量为：112223331221333()()()()()(-1)(-1)(-2)()(-1)(-2)abcdRtRtRtRtRtSJtSJtSJtCSSJtSJtC(1.10)上式中，1S、2、2S、3、3S、C均为常量，令11222333pSSpSpC(1.11)则式(1.10)变为：11233()(-1)(-2)RtpJtpJtp(1.12)而堰塞湖每天新增水的体积与湖水的水位有关，根据式(1.3)可得：0011233()[()-]-[(-1)-](-1)(-2)nnbbRtaHtHaHtHpJtpJtp(1.13)根据附件一中的数据，可以得到一系列()Ht、1()Jt、2()Jt的实际数值，通过与蓄水量-水位高程曲线拟合，得到参数1p、2p、3p的最佳值分别为0.0025，0.0028，0.0649。经过转换可以得到()Ht与降雨量之间的关系为：0ln[((-1)-)]0()nbsHtHnbHteH(1.14)其中：11233(-1)(-1)pJtpJtpsa123400.0025,0.0028,0.06491.908910,2.1894667.4bpppanH由式(1.14)，通过迭代可以预测出降雨量为实际降雨量的50%、80%、100%、150%时的水位高程变化情况。下面分别是预测出的每日水位高程图：8图5各个降雨量时水位高程预测图由上图可以看出，当降雨量为实际的100%时，至6月12日，水位高程已经达到了坝顶。而当降雨量为实际的150%时，水位高程会超过750.2米的坝顶，这时很容易造成溃坝的发生。下面分别是预测出的每日水位上升图：图6各个降水量时每日水位上升预测图92．蓄水漫顶后发生溃坝的数学模型在这里，我们要根据唐家山堰塞湖泄洪时科技人员记录下的大量宝贵数据，在合理的假设下建立堰塞湖蓄水漫顶后在水流作用下发生溃坝的数学模型，模型包括缺口宽度、深度、水流速度、水量、水位高程，时间变量。溃坝过程可以简单描述如下：在上游不断的有水流和雨水流入堰塞湖的作用下，湖内水位不断上升，当水位超过坝顶时，湖内的水将漫过坝顶不断向下游流去。在水流流过坝顶时，水流也会逐渐的侵蚀大坝，带走部分土石，从而形成一个的缺口，即泄流槽，为简化后面的计算，这里我们假设泄流槽的截面形状为矩形。随着泄流槽内水流的不断加大，泄流槽也会不断地加宽加深，从而进一步促进了泄流流量的增大，形成了一个正反馈的机制，即泄流槽的宽度与深度变化量与泄流流量成正比例关系，但是堰塞湖内的水量会随着泄流量的增大而不断减少，最终限制了泄流量的进一步增加。另外，泄流槽在增大到一定宽度和高度时，会因为大坝本身的结构限制其进一步加宽或加深，由附件二中的数据发现，当大坝宽度达到145米时，宽度不再增加，而是逐步加深。这里我们假设在逐渐溃坝时泄流槽的宽度极限为145米。在单位时间内，堰塞湖中的水量的变化等于入水量与出水量之差，即()()()INOUTWtQtQtt(2.1)其中：()Wt表示总水量；INQ表示单位时间内流入堰塞湖的水量；OUTQ表示单位时间内流出堰塞湖的水量，即泄流量；溃坝时，水流不断冲刷着坝体，并带走坝体的泥沙，造成泄流槽缺口的不断地拉宽、变深，根据带走泥沙的量可以列出以下方程式：0()()()tOUTQtdtcbtht(2.2)其中：()Qt表示t时刻的泄流量；表示泄出水量的平均含砂量；()bt、()ht、c分别表示t时刻泄流槽的宽度、深度和顺河长度，顺河长度为一常数，由实际数据可知803c米另外，泄流量可以用泄流槽的宽度与高及相应的水流速度的乘积来表示：()()()()OUTwQtvtbtht(2.3)其中：()vt表示t时刻的水流速度；()bt表示t时刻的泄流槽宽度；()wht表示t时刻的泄流