楊博石，張貴金

（長沙理工大學水利工程學院，湖南 長沙410114）

1 工程背景

雙牌水庫根據(jù)工程的開發(fā)任務和規(guī)模等，考慮大壩自身和下游防洪安全，采用天然河道年最大洪水系列數(shù)據(jù)，通過經(jīng)驗頻率計算得到設計洪水，通過水庫調洪計算得到水庫防洪和調節(jié)庫容及其特征水位。在設計的各個環(huán)節(jié)都偏于安全保守考慮，導致設計洪水值普遍偏大[1]，另一方面，由于多種原因，多數(shù)水庫在汛期實際調度運行時，為了保證大壩安全及下游防洪安全，又進一步降低了汛限水位，減小了水資源利用效率[2～3]。

洪水漫頂風險分析計算對于水庫大壩安全運行和大壩風險管理具有重要意義，方法有直接積分法、結構可靠度分析法、抽樣模擬分析法等，在工程實踐中，數(shù)據(jù)的可用性、模型的復雜性和結果的精度要求是分析方法選取的重要依據(jù)[4]。對漫頂風險影響最大的是洪水的不確定性，關于洪水分析計算，單變量經(jīng)驗頻率分析是傳統(tǒng)水庫洪水計算的常用計算方法，目前研究較多的是洪峰、洪量、洪水過程線形狀等多變量分析方法[4～7]。

本文采用Copula-Monte Carlo法構建入庫洪峰和洪量的聯(lián)合分布模型及漫頂風險計算模型，并量化了不同起調水位下的漫頂風險，探索了提高水庫水資源利用效率及優(yōu)化調度運行的方法，為相關研究提供了應用案例。

2 基于Monte Carlo的大壩洪水漫頂風險計算模型

2.1 大壩漫頂風險分析模型

洪水漫頂風險率pFO是指壩前最高洪水位Hmax超過壩頂高程hc的概率，計算公式如下[8]：

式中FHmax（hc）——Hmax的概率分布函數(shù)。

可以看出，得到FHmax（hmax）是計算洪水漫頂風險率的關鍵問題。本文將采用Copula-Monte Carlo法來模擬洪水峰量聯(lián)合分布和求解pFO。

2.2 洪峰和洪量的聯(lián)合分布

假設設計洪水洪峰和洪量的邊緣分布分別為u=FQ（q）、v=FW（w）。根據(jù)Sklar定理，參考劉章君等人的計算方法[4]，用一個二維Copula函數(shù)C建立Q和W的聯(lián)合分布函數(shù)：

式中Q、W——分別為洪峰和洪量；

q、w——分別為Q、W對應的取值；

F（q，w）——Q和W的聯(lián)合分布函數(shù)；

θ——Copula函數(shù)的參數(shù)。

選 用Gumbel、Clayton、Frank三 種 兩 變 量Archimedean Copula函數(shù)作為備選函數(shù)，其基本形式、參數(shù)θ及θ與Kendall’sτ之間的關系如表1所示。采用Kolmogorov-Smirnov（K-S）檢驗進行邊緣分布模型的檢驗；對Copula函數(shù)擬合優(yōu)劣進行評價的方法有均方根差（RMSE）、K-S法、AIC準則法和Bias值等；K-S值通過檢驗，RMSE、AIC和Bias值越小，表示Copula函數(shù)的擬合程度越好。

當洪峰Q為指定值q時，通過建立的洪峰洪量聯(lián)合分布Cθ（u，v），可以得到洪量W的條件分布，其數(shù)學表達式為[4]：

表1 三種備選兩變量Copula函數(shù)基本形式、參數(shù)θ及θ與Kendall’sτ的關系[4,9～10]

2.3 漫頂風險計算

采用Copula-Monte Carlo法計算洪水漫頂風險率pFO，具體步驟如下：

1）生成兩個服從[0，1]均勻分布的獨立隨機數(shù)r1和r2。

2）令u=r1，SU（v|U=u）=r2。

3）求解SU（v|U=u）=r2，得到：v=SU-1（r2|U=u）。

4）計算q=FQ-1（u），w=FW-1（v）得到一對隨機洪峰洪量組合（q，w）。

5）選取典型洪水過程線，采用變倍比放大方法[5]來計算設計洪水過程線：

式中QD（t）——設計過程在t時刻的流量；

QT（t）——典型洪水過程在t時刻的流量；

WT——典型洪水歷時T內的洪量；

QDm——典型洪水的洪峰流量；

q、w——分別為步驟4）中所求得的洪峰洪量。

6）結合步驟5）中所求得的設計洪水過程線QD（t）、水庫泄流曲線和庫容曲線，按汛期調度運行方式進行調洪演算，得到壩前最高水位hmax。

7）重復步驟1～6共n次，可以模擬出n個hmax，通過數(shù)學期望公式計算經(jīng)驗頻率曲線，得到Hmax的概率分布函數(shù)FHmax（hmax）。

8）計算洪水漫頂風險率pFO=1-FHmax（hc）。

3 雙牌水庫漫頂風險分析

3.1 洪峰洪量分布的確定

參照《水利水電工程設計規(guī)范》（SL-2006），選用雙牌水庫壩址處水文站1959—2019年實測2 h流量數(shù)據(jù)，提取雙牌水庫的年最大洪峰流量Q和3日洪量W系列數(shù)據(jù)，進行經(jīng)驗頻率計算。采用P-Ⅲ分布作為邊緣分布建立洪峰與洪量的聯(lián)合分布，分析洪水漫頂風險。

采用K-S檢驗法對P-Ⅲ函數(shù)的矩法估計參數(shù)進行假設檢驗，由表2可知洪峰和洪量系列的K-S值均低于臨界值為0.246，通過K-S檢驗。

為了選定擬合效果最優(yōu)的Copula函數(shù)，采用表1中的三種常用Copula函數(shù)建立雙牌水庫Q和W之間的二維聯(lián)合分布，利用Kendall相關系數(shù)τ估計其參數(shù)，采用2.2中所提及的方法對Copula函數(shù)模型進行檢驗和擬合優(yōu)度評價，結果見表3。結果表明三種備選Copula函數(shù)K-S值均小于臨界值0.246，通過K-S檢 驗；其 中，Gumbel Copula函 數(shù) 的RMSE、AIC和Bias值最小，為擬合度最優(yōu)的Copula函數(shù)。由圖2可以看出經(jīng)驗頻率值與理論頻率值的擬合情況較好，表明采用Gumbel Copula函數(shù)所建立的聯(lián)合分布是合理可行的。

表2 涔天河水庫洪水統(tǒng)計特征值和P-Ⅲ型分布參數(shù)估計結果

表3 Copula函數(shù)的參數(shù)估計和擬合優(yōu)度檢驗值結果

3.2 洪水漫頂風險計算結果

按2.3節(jié)中的方法步驟，利用3.1中擬合度最優(yōu)的Coupla函數(shù)（Gumbel Copula函數(shù)）建立的條件分布，對比隨機生成10萬組洪峰、洪量組合與實測值的情況如圖3所示，表明Gumbel Copula函數(shù)能夠較好地捕捉洪峰、洪量的相關性結構。

圖2 聯(lián)合觀測值的經(jīng)驗分布和理論分布比較

圖3 隨機模擬和實測的洪峰、洪量點據(jù)比較

圖4 不同運行方式的漫頂風險

3.3 調度運行方式選取

分析雙牌水庫在不同運行方式的漫頂風險率（圖4），對比大壩社會可接受風險10-5和可容忍風險10-4的量級[11～12]，按保守取值，汛限水位確定為167.351 m，比現(xiàn)行設計汛限水位高出1.0 m，發(fā)生漫頂風險的概率在10-5量級，風險可以接受。另一方面，雖然起調水位較高時，漫頂風險在10-4量級，但是水庫在水位高時，下泄流量大，增加了下游的防洪負擔，防洪作用發(fā)揮得小，所以汛限水位可確定為167.351 m，可以保證漫頂風險小相對下泄流量也較?。浑S著洪水預報精準度的提高，可以加強水庫上下游的洪水預報，汛限水位可以動態(tài)控制在167.351～168.351 m之間，隨著起調水位的變化，漫頂風險率變化相對較小，大壩安全狀態(tài)較平穩(wěn)。

若汛限水位確定為167.351 m，比現(xiàn)行設計汛限水位高出1.0 m，經(jīng)測算主汛期期間可多蓄水2 700萬m3，可多發(fā)電1 490萬kW·h。

4 結論

本文基于水庫大壩漫頂風險分析，對雙牌水庫在汛期優(yōu)化調度運行進行了研究。

1）基于Copula函數(shù)的兩變量統(tǒng)計方法建立了洪峰洪量模擬模型，其中擬合度最好的是Gumbel Copula函數(shù)，并基于Copula-Monte Carlo法構建了大壩漫頂風險模型。

2）計算了雙牌水庫不同起調水位的漫頂風險率，結果表明，起調水位越高，水庫大壩漫頂風險率越大，如果雙牌水庫設計汛限水位提高1.0 m，風險可以接受，主汛期期間可多蓄水2 700萬m3，可多發(fā)電量1 490萬kW·h。

隨著社會經(jīng)濟的發(fā)展，水庫管理單位越來越重視水庫資源綜合利用開發(fā)，提高汛限水位可以提高水資源利用效率，獲得巨大的效益，但還須進行更深入更系統(tǒng)的研究，多方論證，獲得更廣泛的認同后方可實施。