基于正交试验法的梯级水库连溃分析模型参数敏感性分析

李炎隆，佘 磊，，周兴波，王 琳，于 沭

（1.西北旱区生态水利国家重点实验室，西安理工大学，陕西 西安 710048；2.中国水利水电科学研究院 岩土工程研究所，北京 100048；3.水电水利规划设计总院，北京 100120）

1 研究背景

随着我国社会生产力的快速发展和水利水电技术的不断创新，流域梯级水电开发在工农业生产和国民经济的发展中发挥着显著的经济效益。然而，一旦发生由于人为原因或不可抗力的自然灾害导致流域梯级水库连溃事故，会给下游带来毁灭性的后果［1］。因此在水利工程建设中，对梯级水库连续溃坝风险分析问题的研究很有必要［2］。

近些年来，众多学者针对土石坝溃决问题进行了研究［3-4］。陈祖煜等［5］基于溃坝过程的实际物理机制，利用Microsoft Excel开发出DB-IWHR模型，提出溃口侵蚀率的双曲线冲刷模型，引入了边坡稳定计算模拟溃口展宽，同时，在DB-IWHR模型和洪水演进分析的基础上，该团队建立了梯级水库连溃风险分析模型［1］，并应用于2014年云南鲁甸“8·3”地震形成的红石岩堰塞湖梯级风险分析［6］。周兴波等［7］、刘磊等［8］基于该模型对2008年唐家山堰塞湖溃决的实测数据进行了参数反演分析，发现冲刷参数对计算结果影响较大。梯级水库连溃分析模型参数较多，参数反演分析比较困难，计算工作量也大。因此，进行模型参数的敏感性分析，准确、快捷地确定参数取值具有工程实际应用价值。

现行参数敏感性研究方法主要有单因素分析法和多因素分析法两大类。传统的参数敏感性分析采用单因素分析法，即指选定一个指标值，使其中一个参数变化，同时假定其他参数保持不变，通过比较基准指标值随参数变化曲线能直观反映各参数的敏感性大小［9］。但此法适用于参数较少，且参数之间需相对独立，这与工程实际情况不符。因此，本文采用正交试验法，以我国西南某流域上虚拟的“下庄（1）—达里—双屯子”［10］三级梯级水库系统为例，进行梯级水库连溃分析模型参数对下游水库风险程度的敏感性分析，为溃坝参数反演分析参数选取和工程应急处置措施提供技术参考。

2 梯级水库连溃分析模型基本原理

梯级水库连溃分析模型包括3个过程：上游梯级溃决过程、上游溃决洪水向下游演进过程、下游梯级连溃过程。如图1所示［10］，分别位于x1、x2处的上下游梯级M和N发生连溃事件，则可分为以下时段过程进行计算。

（1）上游梯级溃决过程。即M坝在t1时刻发生溃决，可按单一梯级溃决洪水分析模型计算。

溃口断面处的流量可近似采用宽顶堰公式计算［5］，即：

式中：C为流量系数，陈祖煜等［5］在唐家山堰塞湖溃坝反演得到C取值为1.43 m0.5/s；B1为溃口过水断面的宽度；H为库水位高程；z为溃口底部高程。

对于溃口冲刷过程，溃口侵蚀率采用双曲线模型进行计算［11］，其公式为：

图1 梯级水库连溃示意

式中：τ为剪应力；t为时间；a、b为侵蚀率参数；ν为扣除临界剪应力后的剪应力，定义为：

式中：k为单位变换因子，取100［7］；τc为临界剪应力。

对于溃口扩展过程，DB-IWHR模型采用简化的Bishop法进行模拟。通过边坡分析程序STAB自动搜索临界圆弧滑裂面，进行溃口扩展数值计算，相关原理和计算方法见文献［12］。

（2）上游溃决洪水向下游演进过程。自t1时刻M坝溃决洪水开始向下游演进，至t2时刻溃坝洪水演进至N坝水库，可按洪水演进模型计算。

洪水演进模型选用天然河道一维浅水模型，其运动控制方程为［1］：

（3）下游梯级连溃过程。当M水库溃决洪水演进至N水库后（t2时刻），N水库将出现库水位壅高、漫坝和连溃3个阶段（t2～t3），每个阶段计算方法如下［1］。

阶段1：M水库溃决洪水进入N水库，N水库水位上涨，该过程总水量满足能量守恒：

式中：W为水库库容，m3；h为库水位，m；q1、q2分别为入库和出库流量，m3/s。

当水位超过坝顶后，出库流量将包含从坝顶溢流部分。为简化计算，此过程不考虑坝顶冲刷，则水量平衡控制方程为：

式中：Bo、zo分别为坝顶漫流宽度和坝顶高程。

在时间步长为Δt内的水位增量Δh通过下式求解：

当Δh很小时，可使用下面的近似表达式：

因此，可得到一个直接求解Δh的公式，即：

阶段2：坝顶受到冲刷，形成缺口，溃坝开始，发生梯级连溃。可按过程（1）进行。

阶段3：梯级连溃洪水向下游演进。可按过程（2）进行。

由上述原理可知，梯级水库连溃分析模型输入参数较多。在进行溃坝洪水参数反演分析时，该模型的取值存在参数数量较多、各参数对计算结果影响程度不同和确定参数值困难的问题。因此，梯级水库连续溃决参数反演分析计算的难点问题之一就是该模型的参数取值。

3 正交试验方法

3.1 正交试验设计正交试验设计是基于数理统计学和正交性原理从大量试验点中选择具有代表性的点，利用正交表安排与分析多因素试验的一种设计方法［13］。由于正交表具有均衡分散性和整齐可比性的构造原则，因此，按照此方法设计的试验次数少，并且能反映事物的客观规律，具有较高的计算效率。

正交试验表作为正交试验法基本分析工具的核心，它具有以下两个特点［14］：每一列（因素）的不同水平在试验中出现的次数相同，体现了均衡性；任意两列（因素）的不同水平组合在试验中出现的次数相同，体现了试验点分布的均匀性。正交表用Ln（tc）表示，其中：L为正交表的代号；n为总共试验次数；t为因素的水平数；c为正交表列数，即可以安排的最多的因素个数。以四因素三水平情况为例，其正交试验表如表1所示［15］。

表1 L9（34）正交试验

3.2 正交试验结果分析方法按照正交表设计的试验进行分析，计算各试验的指标值，据此可判断各因素对指标影响的敏感性大小。正交试验的结果分析方法有极差分析和方差分析两种方法，本文采用极差分析方法对正交试验结果进行敏感性分析。

极差分析法中计算统计参数为［16］：

式中：Kij为因素j在i水平下各次试验结果的平均值；n为因素j在i水平下的试验次数；Yk为第k个试验指标值；ȳ为所有试验结果的平均值。

极差值Rj的计算公式为：

极差值Rj可以反映各因素的敏感性大小。极差值Rj越大，表明该因素的改变对试验指标的影响越大，即该因素的敏感性越大；反之亦然。

4 梯级水库连溃算例分析

4.1 工程概况以我国西南某流域上虚拟规划建设的“下庄（2）—达里—双屯子”三级梯级水库为例，其中下庄（2）水库为一级水库，距离达里大坝约85.6 km，高差约318.9 m；达里水库为二级水库，距离双屯子大坝约89.3 km，高差约380.5 m。其中双屯子水库为流域控制性三级水库，平面布置如图2所示［1］。

达里水库不具有调洪能力，故不考虑其泄洪过程。双屯子水库具有防洪功能，其泄洪设施包括溢洪道、竖井泄洪洞、深孔泄洪洞，全闸泄洪流量约8000 m3/s。此外，由于原规划的下庄（2）水库总库容为30.19亿m3，相对于达里水库和双屯子水库较大，一旦其发生溃决，将导致下游的双屯子水库即使在有预警的条件下仍发生漫顶。因此，为了研究“中等库容—小型库容—大型库容”三级串联梯级水库案例，本文将原有的下庄（2）水库的坝高和库容做出相应的折减，其他参数保持不变，得到对应的下庄（1）水库。各梯级水库水库工程特性如表2所示。

图2 梯级水库平面布置（水平距离单位：km；河床高程单位：m）

表2 各水库工程特性

4.2 正交试验设计

4.2.1 选取试验指标 根据参数敏感性分析中试验指标的选取原则，基于梯级水库连溃分析模型基本原理，考虑到梯级水库调洪的实际情况，因此，选择一级下庄（1）水库溃决流量Q、二级达里水库最高水位H1和三级双屯子水库最高水位H2作为参数敏感性分析的主要试验指标。

4.2.2 确定试验因素和因素水平 已有相关文献提到［7，11］，在计算溃坝洪水分析过程中，双曲线侵蚀率模型和洪水调节模型的参数值对计算结果有较大的影响。因此，本文根据梯级水库连溃实际过程选取洪水频率P、侵蚀率参数值b与双屯子水库起调水位H*这3个主要参数进行敏感性分析。本文以虚拟的“下庄（1）—达里—双屯子”三级梯级水库作为敏感性分析的研究对象，以实际工程参数为基础，在参数允许变化范围内设定3个试验水平。则参数敏感性分析的试验因素和因素水平如表3所示。

4.2.3 选择正交表设计试验 假设模型中各参数之间无交互作用，根据试验因素个数和因素水平数，选择L9（34）正交表安排试验，设最后一列为空列，将试验因素随机分配到正交表中前3列，将正交表中原有的每个元素按其对应的因素和水平替换成相应的设计参数值，即得到该模型参数敏感性分析的正交试验表，表中每一行对应因素水平组合即为一个试验方案［9］。

4.3 梯级土石坝数值计算采用梯级水库连溃分析模型，按照本文设计的试验方案，进行梯级土石坝连溃数值计算。限于篇幅，本文只说明试验方案为P=0（天然入流）、b=0.83 mm/s和H*=2485 m（汛期防洪限制水位）的土石坝溃决数值计算具体过程［1］，其他试验方案计算过程相同。输入计算参数如表4所示。

下庄（1）水库溃决流量过程如图3所示，从图3可知，溃决峰值流量为 .57 m3/s，洪峰时间为12.59 h。在下庄（1）水电站遭遇溃决时，达里水库立即全部打开泄洪设施，达里水库的起调水位设置为2686 m（正常蓄水位），则达里水库出流过程和库水位变化过程分别如图4与图5所示。

从图5可以看到，在遭遇下庄（1）大坝溃决时达里水库库水位为2700.07 m，而达里水库坝高为2690 m，则达里水库已经发生漫顶，漫顶时间为7.3 h。由于下庄（1）水库库容相对较小，则双屯子水库不提前打开泄洪设施。双屯子水库在遭遇下庄（1）和达里水库溃决洪水时再全部打开泄洪设施，双屯子水库的起调水位设置为2485 m（汛期防洪限制水位），则双屯子水库出流过程和库水位变化过程分别如图6和图7所示。

由图7可以看到，在遭遇下庄（1）和达里水库溃决洪水时双屯子水库最高水位达2496.15 m，而双屯子水库坝高为2510 m，则双屯子水库不会发生漫顶。

4.4 敏感性分析将4.3节每种试验方案计算得到的下庄（1）水库溃决流量Q、达里水库最高水位H1和双屯子水库最高水位H2这3个试验指标值分别列在正交试验表的最后三列，如表5所示。然后根据各试验方案计算结果采用极差分析法分析各参数对试验指标的敏感性大小。

对试验指标Q的影响因素进行极差分析，分析结果如表6所示。结果显示，各因素对试验指标Q的敏感性由大到小为：b＞P＞H*。

对试验指标H1的影响因素进行极差分析，分析结果如表7所示。结果显示，各因素对试验指标H1的敏感性由大到小为：b＞P＞H*。

表4 参数输入

图3 下庄（1）溃决流量过程

图4 达里水库出流过程

图5 达里库水位变化过程

图6 双屯子水库出流过程

图7 双屯子库水位变化过程

表5 正交试验设计方案及计算结果

表6 试验指标Q影响因素极差分析结果

表7 试验指标H1影响因素极差分析结果

对试验指标H2的影响因素进行极差分析，分析结果如表8所示。结果显示，各因素对试验指标H2的敏感性由大到小为：H*＞P＞b。

对各试验指标Q、H1和H2的极差分析结果进行整理，按照各因素对各试验指标的极差值大小绘制极差值柱状图，如图8所示。由图8可知：梯级水库连溃分析模型中参数P、b对下庄（1）水库溃决流量的影响显著，对达里水库最高水位的影响较大，其中b的敏感性最大，其次是P；参数H*对下庄（1）水库溃决流量和达里水库最高水位没有影响；参数H*、P对双屯子水库最高水位的影响大，其中H*的敏感性最大，其次是P。

表8 试验指标H2影响因素极差分析结果

图8 各试验指标参数敏感性对比结果

5 结论

本文基于正交试验法，进行了梯级土石坝连溃分析模型参数敏感性计算，研究了梯级水库连溃分析模型中各计算参数对下庄（1）水库溃决流量Q、达里水库最高水位H1和双屯子水库最高水位H2的敏感性，研究结果表明：梯级水库连溃分析模型中参数b、H*对计算结果影响显著，参数敏感性高；而参数P对计算结果影响较小，参数敏感性较低。

因此，在实际梯级土石坝溃坝参数反演分析中，可以重点考虑这些敏感性高的参数取值，不仅能够大大减少反演计算量，还能提高反分析的效率，保证反演结果的正确性。本文的研究方法可以为梯级水库连溃分析模型参数取值和工程应急处置措施提供一种思路。