基于蒙特卡罗法的缓变型防洪风险分析

2018-11-20周方明谷明晗

水资源与水工程学报 2018年5期

周方明，王锋，谷明晗，蔡磊

(1.宿州学院资源与土木工程学院，安徽宿州 234000； 2.江西水利职业学院水利工程系，江西南昌 330013； 3.南昌大学建筑工程学院，江西南昌 330031)

1 研究背景

洪涝灾害是世界发生频次最高、损失最大的自然灾害，我国大部分地区均存在着不同程度的洪涝灾害，对人民生命财产、国家经济建设、地域生态环境构成了重大威胁[1-3]。世界坝工委员会[4]于1998年总结了水库防洪安全设计的进展，主要是根据经验和水库的共性粗略分级制定标准列表备查，并未考虑到水库的个性，如水库大小、坝型、库容以及垮坝的危害。外国学者Kwon等[5]采用了Monte-Carlo模拟随机变量分布，用来研究大坝的水文风险概率；国内学者梅亚东等[6]认为水库大坝在长期服役过程中，影响防洪风险的诸多因素(水力、水文、设计)伴随时间历程发生随机不确定性的变化，致使大坝的防洪能力呈现明显的时变特性[7-9]。时变特性具体呈现为不确定性洪水的荷载效应与结构的抗力作用是随机变化的过程。故防洪建筑物的失事因素与时变效应存在着协同作用，导致了诸多破坏形式和失效模式，如洪水漫顶、边坡失稳、渗透破坏等；其中洪水漫顶是最常发生的失效模式。

因此，防洪风险分析显得尤为重要，对水库大坝防洪能力进行评估能最大限度地降低其风险和灾害影响，具有重大意义。故本文针对水库大坝泄流能力不足或遭遇超标准洪水所引起水库漫坝的风险进行研究。一般的风险分析方法并未考虑长期服役下大坝防洪性能的衰退，也未能对不同服役环境下大坝实际性态的

差异给予定量的分析。据此，本文在时变效应理论的基础上，试图构建时变随机变量量化的函数模型，并提出缓变型防洪风险分析模型，同时运用蒙特卡罗法求解风险率，从而定量分析大坝防洪时变安全度。

2 时变风险率计算方法

2.1 蒙特卡罗法的基本原理

蒙特卡罗法通过构建一个概率模型或随机过程，使其参数等于所求问题的解，然后通过随机试验或观察等来计算所求参数的统计特征，由此得出问题的近似解，而解的精度可由估计值的标准差来表示，并且据此可得结果的置信度及置信区间，从而解决实际工程问题[10-11]。首先由多个相互独立的随机变量X1,X2,…,Xn构成功能函数为Z=g(X1,X2,…,Xn)，再分别对其抽样获得各随机变量的分位值x1,x2,…,xn，接着对功能函数式zi=g(x1,x2,…,xn)进行求解，设抽样次数为N，再计算每组对应的功能函数值zi，其中zi<0的次数为M，在足够多的抽样结束后，其失效概率为Pf=M/N。

首先选用一种特定的方式构造随机数，在[0,1]上均匀分布的随机数是最基础的，其他分布的随机数均可在[0,1]范围内的均匀分布基础上经变换获得。混合同余法[12]以速度快、周期长、统计性质好的特点而被广泛运用。

混合同余法的计算公式为：

xi+1=(axi+b)(modm)

(1)

式中：a、b、m都取非负整数，m≠0。

计算时，引入一个参数ki，令：

(2)

式中：符号int表示取整。

则：

xi+1=axi+b-mki

(3)

将xi+1除以m后，即得[0,1]上的均匀随机数ri+1：

(4)

2.2 随机变量的抽样

(5)

随后将抽样值依次代入功能函数式中，求解成果唯有g(xi)>0和g(xi)≤0两种情况，故可定义指标函数：

(6)

依据大数定律计算其失效概率Pf为：

(7)

式中：M为N次模拟计算中g(xi)≤0的总次数。

本文拟采纳95%的置信度来确保蒙特卡洛法求解在容许误差ε内：

(8)

由公式(8)可知，抽样数目N越大，容许误差ε随之减小。若计算结果要满足既定的精度要求，样本数目N需满足一定要求：

(9)

求解风险率的具体形式如下：

(1)依据实测资料定义功能函数中每个随机变量以及其分布规律；

(2)产生0到1之间的均匀随机数ri；

(3)由ri产生各个随机变量的抽样值；

(4)将随机变量的抽样值依次代入功能函数式，求解功能函数值zi；

(5)统计出现zi<0的次数M和总抽样次数N；

3 缓变型防洪风险分析模型

当坝前水深接近坝顶高程时，大坝防洪风险将显著标志，建立其极限状态方程:

Z=H-D=0

(10)

由于大坝服役过程中内外环境等因素的变化，致使各种随机性因素发生缓慢变异，有学者将该时变效应称之为缓变型时变效应[13-15]。大坝防洪风险的缓变特性随时间逐步发生，或许是结构上的变化，或许是非结构上的。在防洪时变风险分析中，各种不确定性因素将发生缓慢变化，D(n)是坝顶高程一个随机过程，故利用衰减函数模型来表述坝顶高程的沉降缓变历程：

D(n)=a(n,k)D0

(11)

式中：D(n)为服役n年以后随机坝顶高程；D0为大坝的坝顶高程初设值；a(n,k)是衰减函数;k为衰减系数，相应的衰减模型用指数函数来表示：

(12)

式中：T为大坝设计基准期;k需要通过工程经验或试验资料确定。

坝前最高水位变量Hi对极值Ⅰ型分布、对数正态分布和正态分布均不拒绝[16-17]，依据实测资料对各年的坝前最高水位Hi进行排序(H1

(13)

随机变量Hi的均值μH、标准差σH及变异系数δH可应用矩法公式求解：

(14)

(15)

(16)

依照公式(13)求解经验分布F(Hk)，并采用K-S检验法[18]依次检验相应的经验分布F(Hk)，建立统计量：

Dn=

(17)

根据显著性水平(一般取95%)，查K-S检验临界值表得Dn,95%与Dn，通过对比来判辨分布概型的拟合度，从而获得随机变量Hi的最优概率分布。假定H与D是相互独立且服从正态分布的随机变量，D的均值μD和标准差σD由以下公式确定：

μD(n)=μD0·a(n,k)

(18)

σD(n)=σD0·a(n,k)

(19)

随机变量H的均值μH和标准差σH由荷载评估期TE内荷载统计分析方法求解：

(20)

(21)

由此可知，大坝第n年的防洪风险率的计算公式如下：

Pf(n)=Φ(-βn)=

(22)

若随机变量H和D分别服从对数正态分布和正态分布，经公式(4)产生[0,1]上的均匀随机数ri，再运用坐标变换法获得正态分布N(0,1)的两个随机数ui1，ui2：

(23)

(24)

由于随机变量H服从对数正态分布，令H′=lnH，可知H′为正态分布，则有：

(25)

(26)

(27)

具体计算流程见图1。

4 实例分析

4.1 工程概况及计算参数

金溪流域上池潭混凝土坝已正常运行多年，选取该坝的典型坝段进行防洪风险分析，典型坝段剖面示意图见图2。

工程主要建筑物为二级永久建筑物，坝顶高程为384.5 m，岩基建基面高程314.5 m，正常蓄水位382.0 m，设计洪水位(P=1%)为382.0 m，校核洪水位(P=0.1%)为384.0 m，死水位354.0 m，为年调节水库。

图1 计算流程图

图2 典型坝段剖面示意图

选取1990年至2011年间的坝前实测年最高水位以及对应的坝顶实测年沉陷最大值进行分析，如表1所述。

表1 服役期实测资料

对上文所述的3种分布函数进行K-S检验，其检验成果见表2。

查K-S检验临界值表获得允许概率差值D22.95%=0.2807，与表2的计算成果对比分析，3种分布的最大概率差值Dn均小于D22.95%，表明随机变量H和D均服从这3种分布的假定，且正态分布的Dn值小于对数正态分布和极值I型分布的Dn，说明实测序列值与正态分布假定之间差别最小，随机变量均为正态分布时的拟合效果最好。

表2 随机变量分布检验表

本工程中的设计基准期为50 a(TN=50)，结构荷载评估基准期TE可采用类比法式求解：

(28)

随机变量H在继续使用期TM内的统计参数μH、σH可依据公式(20)、(21)计算，其计算结果见表3。

由于缺乏坝顶沉降衰减模型的资料，故选用公式(12)模拟沉降衰减；根据已服役22 a内坝顶沉降统计的的均值近似确定坝顶高程的衰减函数a(n,k)。

已知μD0=384.5 m，μD(22)=384.495 m，σD(22)=7.43×10-4m，则由公式(18)、(19)可得：

μD(n)=μD0·a(n,k)

(29)

σD(n)=σD0·a(n,k)

(30)

4.2 大坝防洪风险率评估

选用表3所得的荷载评估期统计参数，依据公式(20～22)来推求大坝继续使用期内的防洪风险率，并计算现行规范设防标准下的防洪风险率，其求解成果如表4及图3、4所示。

表3 随机变量H在继续使用期内的统计参数表

表4 大坝继续使用期内防洪风险率统计表

图3 大坝继续使用期内的防洪风险率

图4 大坝后续服役期防洪风险率与现行规范设防标准对比图

由图3和4可知，大坝在服役过程中其防洪性能伴随时间不断的衰退，防洪风险率逐步增大，算例中的大坝在后续服役期内的防洪风险率完全符合现行规范设防标准，且存在一定的安全余度。为进一步确保大坝生命周期内的防洪安全，减小洪水漫顶风险率，建议在工程建设前期充分考虑泥沙作用效应，并预设一定的淤积库容；也可在后续服役期内增加坝顶高程，以增大防洪库容来提升大坝防洪能力。