水库调洪过程库水位随机分布特征分析.doc

此文档收集于网络，如有侵权，请联系网站删除水库调洪过程库水位随机分布特征分析摘 要： 水库坝前最高水位概率分布的推求是水库防洪安全风险计算中的一个关键环节。根据调洪过程各时段库水位的随机分布都可看作是调洪最高水位随机分布的研究实例的观点，本文首先分析了调洪过程由于风险因子影响而导致的时段风险传递机理，给出了时段分布风险的计算公式，接着从同一场洪水的调洪演算角度入手，应用随机模拟和最大熵分布的方法，进一步探讨了调洪过程库水位的随机变化特征，以此分析水库坝前最高水位的分布规律。最后，结合安康水库设计洪水的分析，说明了该随机模拟调洪演算模型的实施过程，结果表明，调洪最高水位有时并不具有正态分布的特征，所服从的分布也不具有通用性，实际中应根据各个水库和洪水的特性具体分析。关键词： 调洪演算；漫顶风险；风险传递；最大熵分布；随机模拟The characteristics analysis of stochastic probability distributions of the reservoir levels during the flood routing Xi Qiuyi, Huang Qiang, Wang Yimin, Lin leilei( Insitute of Water Resources and Hydroelectric Engineering, Xian University of Technology, Xian 710048 )Abstract: The derivation of probability distribution of highest flood level is a key step in the calculation of risk for reservoir flood-control safety. According to the viewpoint that other intervals flood levels may be seen the samples of highest flood level excluding the initial reservoir level during flood routing of reservoir, firstly, the mechanism of risk transfer induced by risk factors is analyzed and the formulation of interval distribution risk is given, then to the same flood, with the stochastic simulation and Maximum Entropy Distribution methods, the characteristics of random distribution of reservoir levels is further investigated in the paper. At last, taking the Ankang reservoir as a case study, the implementing procedure of the simulation model is demonstrated. It is shown from the example that the highest reservoir flood level doesnt always conform to the normal distribution, and it hasnt a general distribution form in all reservoirs, the analysis should be taken for various reservoirs and inflow floods in practice. key words: flood routing; overtopping risk; risk transfer; Maximum Entropy distribution; Stochastic simulation 1 引言水库是重要的防洪工程措施，一方面可减轻水资源短缺和水库下游遭受洪灾等风险，另一方面又不可避免地形成了大坝失事的新风险。大坝失事具有损失巨大，后果严重，发生概率低的特点。对于水库大坝的一种极值荷载情况极值洪水作用下的防洪安全（漫坝）风险进行评估，国内外已做过广泛的研究1-4，在对风险因素的认识和描述以及计算方法上都取得了很大进展。根据水库防洪安全风险的定义，计算漫坝风险，需要知道水库调洪最高水位的概率分布和坝顶高程（或防浪墙高程）的概率分布，坝顶高程（或防浪墙高程）的概率分布相对比较容易确定，而水库坝前调洪最高水位由于实际资料较少，且往往不包括极值洪水，代表性不足，理论分布常常难以确定。吴世伟、张思俊5等采用K-S非参数检验法检验了我国82座具有19年运行资料的水库坝前年最高水位，认为一般水库的坝上游年水位峰值的分布规律为正态分布或对数正态分布，少数湖泊型水库的库水位峰值分布规律为极值I型分布或正态分布。丁晶、邓育仁等6统计分析了160座遍及我国大部分地区和各大江河，代表各种不同类型、性能的水库的坝前年最高水位系列，认为坝前年最高水位变量分布呈现负偏，可视作服从P-III型分布和三参数的对数正态分布，推荐用三参数的对数正态分布来表征。熊明7则认为坝前水位序列是受人类活动影响的序列，其能否符合一般的理论分布是一个值得探讨的问题，直接采用曲线拟合方法推求这些序列的分布缺乏必要的科学依据，建议通过模拟的方法，推求其经验分布来代替理论分布。麻荣永8在计算土石坝的漫坝风险时，假定由确定的一场洪水引起的坝前最高水位服从于正态分布。姜树海，范子武等2,4采用随机微分方程计算水库的泄洪风险时，依据的前提是水库蓄洪量的变化符合Wiener过程的定义，Wiener过程同时是正态过程，因此，实质上也认为坝前调洪最高水位服从正态分布。综合上面研究成果，可以看出，多数学者认同直接分析实测资料的统计方法，但该序列是否具有理论分布或服从何种常用理论分布存在分歧，坝前年最高水位的分布在各水库有没有通用性，过去研究也没有予以回答。坝前调洪最高水位取决于入库洪水大小和调洪过程。从调洪演算的过程看，水库调洪最高水位的随机变化，是由调洪最高水位出现前逐时段风险因素的随机影响经过累积和抵消作用形成的，调洪演算其余时段库水位的随机变化类似于最高水位，来自于其各自时段前面时段风险的传递，只不过传递的时段数不同，从而，各时段（初时段除外）库水位的随机分布特征可看作是调洪最高水位随机分布特征的一个实例，因此，研究和描述水库调洪演算全过程库水位的随机分布特性和风险传递机理，可以用来推求坝前调洪最高水位的概率分布，具有理论上与实用上的价值。本文基于这种研究思想，采用随机模拟调洪演算，结合最大熵分布的理论方法建立了调洪演算过程库水位随机分布计算模型，通过求解此模型，获得调洪过程库水位的随机分布特征，以进一步探讨坝前年最高水位的概率分布和水库自身的防洪风险计算问题。2 水库调洪演算过程风险传递分析与计算2.1 水库调洪演算过程的风险传递分析水库调洪演算是对入库洪水在水库滞蓄作用下的出库变形过程和在此期间库水位变化过程的数值模拟，其依据的基本资料有设计洪水过程、泄洪建筑物的泄流能力曲线、水位库容关系和防洪调度规则等资料，通常采用的计算方法有试算法和数值调洪方法（最早由陈守煜教授提出），水库调洪演算广泛应用于水库规划设计和调洪方案的模拟中。传统水库调洪演算过程是一确定过程，不考虑调洪所依据各种资料的不确定性，从而当调洪演算成果用于设计和运行管理时，带有一定程度的风险。综合考虑调洪演算基础资料中水文、水力以及调度规则，初始条件等不确定对调洪演算的影响，用随机的观点审视调洪演算过程，对于改进设计成果，合理评估水库大坝的防洪安全均有良好的指导作用。严格意义上的水库水量平衡方程是微分方程，而该方程很难采用解析方法求解，通常的做法是采用差分的方式求其近似解。这样，水库调洪演算过程在具体实现时是一个离散的逐时段求解过程。随着调洪过程中入库与出库水量的差异和调洪过程中各种风险因素的影响，水库蓄水量W的变化过程，显然是一随机过程，任意调洪时段的蓄水增量W不仅相互独立，且整个变化过程具有Markov过程的性质，即任意时刻蓄水量W(tn)的条件概率，仅与最近的值有关，而不依赖于水库较早时的蓄水量值，可用下式表示： （1）其中，f(Wn,tn|Wn-1,tn-1)为条件转移概率密度函数。根据以上水库蓄水过程的性质和水量平衡公式： （2）用随机的观点可以看出，调洪时段末水库蓄水量的不确定性来自于两部分，一是来自于调洪时段初水库蓄水量的不确定性，另一来自于当调洪时段蓄水增量的不确定性。调洪时段初水库蓄水量的不确定性又来自于前一时段水库蓄水增量的不确定性和前一时段初水库蓄水量的不确定性，依次累推，则某调洪时段末水库蓄水量的随机变化，是由初始调洪时刻水库蓄水量在调洪过程中，逐步传递和抵消该调洪时段前各时段蓄水增量的随机不确定性而引起的。水库水位与蓄水量存在着对应关系，对水库蓄水量随机变化的分析，很容易拓展到水库水位的随机分析。由此可以得出：调洪演算过程中，由于调洪所依据资料的不确定性，致使库水位的演算结果呈现出随机性，库水位的变化过程是一随机过程。这里需要注明的是调洪演算模型的不确定性在本研究中不予考虑。调洪模型的不确定性来自于用公式或方程表达水库洪水调节物理变化过程所做的概化，本文认为在资料准确的前提下，数值调洪演算真实模拟了水库洪水蓄泄过程。2.2 水库调洪演算时段风险计算公式水库数值调洪演算方程（Euler法）如式（3）所示： （3）其中：hn-1,hn分别表示时段n-1的初、末水位；In-1,Qn-1分别表示时段n-1的入库流量和出库流量；A(hn-1)表示库水位hn-1对应的水面面积；t表示调洪时段步长。将式（3）转换成极限状态方程，则有： （4）若G0则表示破坏，表明在当前时段末实际水库水位超过了预期调洪水位，存在着风险。本时段的末水位又是下一时段的初水位，当前时段风险通过时段末水位传递到下一时段，下一时段又传递到再下一时段，依次累推，实际调洪最高水位则存在大于模拟调洪最高水位的风险。当前时段的风险可按下式计算： （5）其中，表示当前时段的风险；表示当前时段的分布风险；表示上时段的传递风险；表示当前时段的风险与当前时段的分布风险和上时段的传递风险之间的函数关系。当前时段库水位变化的概率分布是在时段初水位确定下的概率分布，即以时段初水位为条件的条件概率分布。在入库洪水确定的情况下，导致当前时段的分布风险的风险因子有两个，出库流量Q和水库水面面积A(h)，这两个风险因子可看作是相互独立的，并且均服从于正态分布。当前时段风险的表达式如下： （6）若分别用随机变量X1，X2代替出库流量Q和水面面积A(h)，用a代替入库流量I，b代替，令： （7） （8）那么， （9） （10）则 （11）当前时段下泄流量Q和水面面积取决于时段初库水位，由于当前时段初库水位的风险（即上一时段的传递风险）与当时段的分布风险之间的函数关系无法确定，对于这种风险变量存在复杂影响机制的问题，采用随机模拟的方法，模拟时段初库水位与水位增量变化的组合关系，显然在处理这种问题时有很大的优势。下节将通过随机模拟的思想，进一步研究风险变化问题。3 水库调洪过程库水位的随机分布特性分析对于水库调洪过程库水位的随机分布的分析，文献4提出了水库调洪演算的随机微分方法模型，该模型的基本前提是认为水库调洪演算过程中蓄洪量的变化过程是一Wiener过程。根据维纳过程定义9，各时段蓄洪增量应服从均值为0，方差为2t的正态分布，若调洪时段间隔t不变，那么各调洪时段2相等，显然这不符合水库实际情况水库不同水位对应的泄洪能力和库容随机分布特性事实上并不相同，可以得出，调洪演算随机微分方程中的随机作用强度系数的取值，是对调洪过程蓄水量的随机离散程度作了平均，进行了相当大的简化；另一方面，对该模型的求解采用有限差分的数值方法，这并不能保证模型解各时刻库水位的分布密度一定为非负值，即满足概率密度函数本身的特点，特别是当水库调洪规则包含有复杂的控泄条件和水库时段水位变幅较大，不符合模型所作假设的情况时。针对上述不足，本文提出了模拟调洪演算的思想，并建立了相应模型。3.1 二阶Runge-Kutta调洪演算方程从水库调洪作用的分析来看，水库之所以能够削减洪峰流量，主要是由于水库蓄洪和泄洪建筑物控制了下泄流量。当水库泄洪建筑物的形式、尺寸大小一定时，下泄的流量受控于泄流水头，泄流水头又受到水库蓄水量的控制，而水库蓄水量的大小又紧紧依赖于入库与出库水量的变化。因此，水库调洪演算实质上是对两个方程的求解，水量平衡方程和水库蓄泄方程。传统的调洪演算图解法现已较少使用，试算法因其是隐式求解方法，并且在每一时段计算时，常需反复计算同一函数，在进行调洪演算各时段风险分析时不方便，本文选用精度相对较高的显式数值调洪二阶Runge-Kutta法作为模拟模型。记 （12）则二阶Runge-Kutta调洪演算数值求解公式则为： （13）其中,分别表示时段初、末的入库洪水流量；,分别指时段初水位和第一次时段末计算水位对应的水面面积；,分别表示时段初库水位对应的泄流能力和第一次时段末计算水位对应的泄流能力。二阶Runge-Kutta数值调洪演算在每一时段进行两步运算，第1步按照（12）式计算时段末水位值，第2步是将第1步的水位计算值代入（13）式，从而得到修正后的时段末水位值。本时段计算完成后，本时段计算水位将作为下一时段初的已知条件，重复上述计算，直到洪水演算过程结束。在入库洪水确定的情况下，（12）式和（13）式中除时段步长和入库流量确定外，其余因子均是随机变量，因此时段末水位hn不再是一确定的量，作为下一时段初水位，水位hn的随机性又对下一时段末水位hn+1产生影响，hn+1还受其它随机变量影响，同样呈随机变化特性，hn+1进而又影响到hn+2等等。3.2 调洪过程风险影响因子随机特性分析及模拟调洪演算方程包括入库洪水、泄流能力和库水位-水面面积关系这三个风险因素。其中入库洪水的不确定性是由于水文资料的限制和人们主观选择典型洪水过程的差异，采用不同的计算方法等原因引起的。泄流能力的不确定性是由所采用的水力计算公式和公式参数不准确，施工的误差，水流的冲涮等原因引起的，可认为每一水位的泄洪能力服从于正态分布。库水位-水面面积对应关系受测量误差和泥沙冲淤等影响，分析起来难度很大，根据水位与库容的对应关系，也可假设在某一水位上对应的库水面面积服从正态分布。在本文中，入库洪水过程按确定过程考虑，因而，所分析的概率密度函数是对同一场洪水而言的。模拟的具体流程如下：Step1：按各自的随机分布特性随机生成m条泄流能力曲线和k条水位水面面积关系曲线，剔除其中不符合物理背景的曲线，再用同样的方法生成予以补充。Step2：在m条泄流能力曲线中和k条水位水面面积关系曲线中各随机选取一条与同一场入库洪水过程进行组合，经确定的调洪演算后得到时段末库水位的一个样本。Step3：根据模拟次数要求重复步骤2 ，得到时段末库水位的一组大容量样本。Step4：采用最大熵法分别分析各时段末的库水位样本，得到各时段末库水位的概率密度函数的解析表达式。3.3 最大熵分布采用随机模拟结合调洪演算生成的库水位样本反映了随机变量的基本变化信息，分析这些样本服从的概率分布，信息论中的最大熵法具有很大的优势，因为当熵取极大值时，能够在满足给定信息所提供的约束条件下，得到其最小有偏概率分布【10】。这样不仅有效减小了信息不足的影响，而且能得到具体的分布形式的解析表达。最大熵概率分布的基本模型如下：O.b. （14）S.t. （15） （16）其中：f(x)是样本的概率密度函数，m是所使用的矩数，mi表示第i阶原点矩，它们是由样本数据以数字形式定义的。要调整f(x)使熵达到最大,令为目标函数，利用拉格朗日乘子，有 （17）推导（17）式，可得到概率密度函数的解析表达式： （18）只需确定其中的参数就可以完全确定f(x)。经数学变换，推导出求解拉格朗日乘子的联立方程为： （19） （21）求上式残差平方和 （22）的最小值可得到拉格朗日乘子的解。4 漫坝风险计算水库大坝漫顶的极限状态方程为： （23）其中：Z1表示坝顶高程（或防浪墙顶高程），通常认为服从于正态分布，均值是设计坝顶高程，方差由具体施工与运行情况确定；Z表示坝前调洪最高水位，与坝顶高程变量相互独立，概率分布函数可取由前述方法计算出来的最大熵分布。通过JC法求解式（23），可得到某一场洪水下的年漫坝风险Pi。用随机水文学的方法模拟生成多场入库洪水，每场洪水均计算Pi，再通过累加来计算年漫坝总风险。计算公式如下： （24） （25）式中，f(Qi)是洪峰流量的概率密度函数在Qi处的值；n为模拟的有限洪水场次（洪峰流量不同），若模拟的洪水中有洪峰流量相同的情况，则洪峰流量的概率密度函数值仅统计一次。Pi第i场洪水的漫坝概率，当洪峰流量相同时，Pi取根据这些相同洪峰流量的洪水过程计算出的最大漫顶概率，即取maxPj，Pj是洪峰流量相同的第j场洪水过程的漫顶概率。5 计算实例与结果分析本文以安康水库调洪演算为例，通过分析调洪过程各时段末的库水位随机特性来分析水库调洪演算调洪最高水位所服从的概率分布。其中，入库洪水采用设计洪水过程（千年一遇），泄洪建筑物的泄流能力曲线以设计时的泄流能力曲线为均值在各水位按正态分布进行随机生成，其中变差系数取0.01；水位与水面面积关系曲线按同样的方式生成，变差系数取0.01。对安康水库设计洪水的调洪演算中，最高水位（均值332.71m）出现在第23时段（时段长3小时），根据调洪最高水位的概率密度函数，实际调洪最高水位大于设计洪水位333.10m的风险为0.39%，漫顶风险接近于0。部分计算结果如表1和图15所示。表1库水位概率分布函数参数Tab.1 The PDF parameters of the reservoir levels during flood routing时段序号概率密度函数参数0123-273884350.51685443.838-2592.9930747-46227206.35284454.7968-437.591482912-552424.92913377.315888-5.16189978818-2982776.6518086.70019-27.418127523-4032284.14924238.72303-36.4257262826-4502448.44627237.88803-41.19438232 图1 库水位变化过程 图2 调洪演算时段风险变化过程（相对于各时段期望水位）Fig. 1 The reservoir levels h with time t Fig. 2 The change of risks during flood routing relative to expected water levels图3 库水位（12时段）的概率密度函数 图4 调洪演算过程库水位分布的方差和偏态系数Fig. 3 The PDF of reservoir level at the end of Fig. 4 The variance and coefficient of skewness of reservoirtime interval 12 levels during flood routing 其中，图4是第12时段末库水位的概率密度函数曲线。分析以上各图可以得到如下结论：（1）从图1和图4可以看出，对同一场洪水过程，采用本文所述随机模拟调洪演算方法得到的库水位的均值过程，与常规调洪演算方法得到的结果基本一致，随机调洪模拟演算还能额外给出反映库水位过程离散程度的方差和反映分布外形的偏态系数信息。（2）以常规调洪演算得到的各时段库水位作为预期水位，采用随机模拟调洪演算分析得到各时段末库水位可能大于预期库水位的风险，从图2可以看出，这些风险值接近50%，与前面解析法得到的时段分布风险值相差不大（时段分布风险值末列出），说明上时段的风险通过时段初水位传递到末水位时，不是经过简单的累加运算和乘积运算，存在着抵消作用。（3）图4反映的是调洪演算过程库水位的方差和偏态系数变化情况，可见，库水位的方差并不像库水位的变化过程，而是具有波动的现象，分析其原因，一方面是由于调度规则中包含有多级控泄要求，另一方面是由于入库流量I、出库流量Q存在着相互影响，相互制约关系，风险因素的随机变动在某些水位对库水位的变动相对影响不大。（4）图3所示的是12时段末库水位的概率分布曲线，其它时段的结果限于篇幅没有列出，直观上看，各时段库水位概率密度曲线外形与正态分布相似，但有些时段却呈有偏分布，如图3。从图4中的库水位偏态系数也能看出，一些偏态系数接近于0，一些偏态系数偏离0较大，且多属于负偏分布。除调洪演算初始起调水位外，调洪演算最高水位前各时段库水位同样是经过许多时段逐步求得的，因而其随机分布特性，实际上反映了调洪最高水位的分布特性，可以得出，调洪最高水位的随机分布虽在分布曲线外形上相似于正态分布，但有些情况相对于均值情况不对称，呈有偏分布，推广到其它水库，则水库坝前调洪最高水位的概率分布不具有通用性，要根据水库具体特性具体分析确定。6 结语本文在传统随机模拟调洪演算的基础上，通过对调洪演算过程的风险因子进行随机模拟生成，研究了同一场洪水调洪过程库水位的概率分布特性，得到了库水位服从的最大熵分布解析表达式。这种研究思路体现了同一场洪水可能对应多个调洪最高水位的观点，分析同一场洪水得到的坝前调洪最高水位的理论分布，也能有效避免人们对根据所有发生