楊博石,張貴金
(長沙理工大學(xué)水利工程學(xué)院,湖南 長沙410114)
雙牌水庫根據(jù)工程的開發(fā)任務(wù)和規(guī)模等,考慮大壩自身和下游防洪安全,采用天然河道年最大洪水系列數(shù)據(jù),通過經(jīng)驗頻率計算得到設(shè)計洪水,通過水庫調(diào)洪計算得到水庫防洪和調(diào)節(jié)庫容及其特征水位。在設(shè)計的各個環(huán)節(jié)都偏于安全保守考慮,導(dǎo)致設(shè)計洪水值普遍偏大[1],另一方面,由于多種原因,多數(shù)水庫在汛期實際調(diào)度運行時,為了保證大壩安全及下游防洪安全,又進一步降低了汛限水位,減小了水資源利用效率[2~3]。
洪水漫頂風(fēng)險分析計算對于水庫大壩安全運行和大壩風(fēng)險管理具有重要意義,方法有直接積分法、結(jié)構(gòu)可靠度分析法、抽樣模擬分析法等,在工程實踐中,數(shù)據(jù)的可用性、模型的復(fù)雜性和結(jié)果的精度要求是分析方法選取的重要依據(jù)[4]。對漫頂風(fēng)險影響最大的是洪水的不確定性,關(guān)于洪水分析計算,單變量經(jīng)驗頻率分析是傳統(tǒng)水庫洪水計算的常用計算方法,目前研究較多的是洪峰、洪量、洪水過程線形狀等多變量分析方法[4~7]。
本文采用Copula-Monte Carlo法構(gòu)建入庫洪峰和洪量的聯(lián)合分布模型及漫頂風(fēng)險計算模型,并量化了不同起調(diào)水位下的漫頂風(fēng)險,探索了提高水庫水資源利用效率及優(yōu)化調(diào)度運行的方法,為相關(guān)研究提供了應(yīng)用案例。
洪水漫頂風(fēng)險率pFO是指壩前最高洪水位Hmax超過壩頂高程hc的概率,計算公式如下[8]:
式中FHmax(hc)——Hmax的概率分布函數(shù)。
可以看出,得到FHmax(hmax)是計算洪水漫頂風(fēng)險率的關(guān)鍵問題。本文將采用Copula-Monte Carlo法來模擬洪水峰量聯(lián)合分布和求解pFO。
假設(shè)設(shè)計洪水洪峰和洪量的邊緣分布分別為u=FQ(q)、v=FW(w)。根據(jù)Sklar定理,參考劉章君等人的計算方法[4],用一個二維Copula函數(shù)C建立Q和W的聯(lián)合分布函數(shù):
式中Q、W——分別為洪峰和洪量;
q、w——分別為Q、W對應(yīng)的取值;
F(q,w)——Q和W的聯(lián)合分布函數(shù);
θ——Copula函數(shù)的參數(shù)。
選 用Gumbel、Clayton、Frank三 種 兩 變 量Archimedean Copula函數(shù)作為備選函數(shù),其基本形式、參數(shù)θ及θ與Kendall’sτ之間的關(guān)系如表1所示。采用Kolmogorov-Smirnov(K-S)檢驗進行邊緣分布模型的檢驗;對Copula函數(shù)擬合優(yōu)劣進行評價的方法有均方根差(RMSE)、K-S法、AIC準(zhǔn)則法和Bias值等;K-S值通過檢驗,RMSE、AIC和Bias值越小,表示Copula函數(shù)的擬合程度越好。
當(dāng)洪峰Q為指定值q時,通過建立的洪峰洪量聯(lián)合分布Cθ(u,v),可以得到洪量W的條件分布,其數(shù)學(xué)表達(dá)式為[4]:
表1 三種備選兩變量Copula函數(shù)基本形式、參數(shù)θ及θ與Kendall’sτ的關(guān)系[4,9~10]
采用Copula-Monte Carlo法計算洪水漫頂風(fēng)險率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]來計算設(shè)計洪水過程線:
式中QD(t)——設(shè)計過程在t時刻的流量;
QT(t)——典型洪水過程在t時刻的流量;
WT——典型洪水歷時T內(nèi)的洪量;
QDm——典型洪水的洪峰流量;
q、w——分別為步驟4)中所求得的洪峰洪量。
6)結(jié)合步驟5)中所求得的設(shè)計洪水過程線QD(t)、水庫泄流曲線和庫容曲線,按汛期調(diào)度運行方式進行調(diào)洪演算,得到壩前最高水位hmax。
7)重復(fù)步驟1~6共n次,可以模擬出n個hmax,通過數(shù)學(xué)期望公式計算經(jīng)驗頻率曲線,得到Hmax的概率分布函數(shù)FHmax(hmax)。
8)計算洪水漫頂風(fēng)險率pFO=1-FHmax(hc)。
參照《水利水電工程設(shè)計規(guī)范》(SL-2006),選用雙牌水庫壩址處水文站1959—2019年實測2 h流量數(shù)據(jù),提取雙牌水庫的年最大洪峰流量Q和3日洪量W系列數(shù)據(jù),進行經(jīng)驗頻率計算。采用P-Ⅲ分布作為邊緣分布建立洪峰與洪量的聯(lián)合分布,分析洪水漫頂風(fēng)險。
采用K-S檢驗法對P-Ⅲ函數(shù)的矩法估計參數(shù)進行假設(shè)檢驗,由表2可知洪峰和洪量系列的K-S值均低于臨界值為0.246,通過K-S檢驗。
為了選定擬合效果最優(yōu)的Copula函數(shù),采用表1中的三種常用Copula函數(shù)建立雙牌水庫Q和W之間的二維聯(lián)合分布,利用Kendall相關(guān)系數(shù)τ估計其參數(shù),采用2.2中所提及的方法對Copula函數(shù)模型進行檢驗和擬合優(yōu)度評價,結(jié)果見表3。結(jié)果表明三種備選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ù)估計結(jié)果
表3 Copula函數(shù)的參數(shù)估計和擬合優(yōu)度檢驗值結(jié)果
按2.3節(jié)中的方法步驟,利用3.1中擬合度最優(yōu)的Coupla函數(shù)(Gumbel Copula函數(shù))建立的條件分布,對比隨機生成10萬組洪峰、洪量組合與實測值的情況如圖3所示,表明Gumbel Copula函數(shù)能夠較好地捕捉洪峰、洪量的相關(guān)性結(jié)構(gòu)。
圖2 聯(lián)合觀測值的經(jīng)驗分布和理論分布比較
圖3 隨機模擬和實測的洪峰、洪量點據(jù)比較
圖4 不同運行方式的漫頂風(fēng)險
分析雙牌水庫在不同運行方式的漫頂風(fēng)險率(圖4),對比大壩社會可接受風(fēng)險10-5和可容忍風(fēng)險10-4的量級[11~12],按保守取值,汛限水位確定為167.351 m,比現(xiàn)行設(shè)計汛限水位高出1.0 m,發(fā)生漫頂風(fēng)險的概率在10-5量級,風(fēng)險可以接受。另一方面,雖然起調(diào)水位較高時,漫頂風(fēng)險在10-4量級,但是水庫在水位高時,下泄流量大,增加了下游的防洪負(fù)擔(dān),防洪作用發(fā)揮得小,所以汛限水位可確定為167.351 m,可以保證漫頂風(fēng)險小相對下泄流量也較??;隨著洪水預(yù)報精準(zhǔn)度的提高,可以加強水庫上下游的洪水預(yù)報,汛限水位可以動態(tài)控制在167.351~168.351 m之間,隨著起調(diào)水位的變化,漫頂風(fēng)險率變化相對較小,大壩安全狀態(tài)較平穩(wěn)。
若汛限水位確定為167.351 m,比現(xiàn)行設(shè)計汛限水位高出1.0 m,經(jīng)測算主汛期期間可多蓄水2 700萬m3,可多發(fā)電1 490萬kW·h。
本文基于水庫大壩漫頂風(fēng)險分析,對雙牌水庫在汛期優(yōu)化調(diào)度運行進行了研究。
1)基于Copula函數(shù)的兩變量統(tǒng)計方法建立了洪峰洪量模擬模型,其中擬合度最好的是Gumbel Copula函數(shù),并基于Copula-Monte Carlo法構(gòu)建了大壩漫頂風(fēng)險模型。
2)計算了雙牌水庫不同起調(diào)水位的漫頂風(fēng)險率,結(jié)果表明,起調(diào)水位越高,水庫大壩漫頂風(fēng)險率越大,如果雙牌水庫設(shè)計汛限水位提高1.0 m,風(fēng)險可以接受,主汛期期間可多蓄水2 700萬m3,可多發(fā)電量1 490萬kW·h。
隨著社會經(jīng)濟的發(fā)展,水庫管理單位越來越重視水庫資源綜合利用開發(fā),提高汛限水位可以提高水資源利用效率,獲得巨大的效益,但還須進行更深入更系統(tǒng)的研究,多方論證,獲得更廣泛的認(rèn)同后方可實施。