甘富萬，張華國，覃麗娜，高 揚(yáng)，黃宇明，肖 良

(1.廣西大學(xué)土木建筑工程學(xué)院，廣西 南寧 530004；2.廣西大學(xué)工程防災(zāi)與結(jié)構(gòu)安全教育部重點(diǎn)實(shí)驗(yàn)室，廣西 南寧 530004)

0 引 言

水庫的防洪標(biāo)準(zhǔn)[1]是與大壩的防洪安全密切相關(guān)的。大壩防洪安全最值得關(guān)注的是洪水漫頂[2]，洪水漫頂是影響水利工程安全的一個(gè)重要因素。因此，洪水漫頂風(fēng)險(xiǎn)率的計(jì)算分析對于大壩的防洪安全至關(guān)重要。大壩防洪風(fēng)險(xiǎn)分析的方法主要有直接積分方法、一次二階矩方法和蒙特卡羅模擬方法。直接積分法[3]可求得精確解，但由于涉及多重積分的問題，只適用于低維且調(diào)洪函數(shù)為簡單表達(dá)式的情況。一次二階矩方法要求計(jì)算函數(shù)對變量的偏導(dǎo)數(shù)，只能給出近似結(jié)果，如莫崇勛、董增川等[4]采用一次二階矩法計(jì)算澄碧河水庫漫壩風(fēng)險(xiǎn)。蒙特卡羅模擬方法是通過對實(shí)測資料統(tǒng)計(jì)分析，確定合適的隨機(jī)模型，然后應(yīng)用隨機(jī)模型，生成足夠多的、各種可能組合形式的洪水過程，供水庫調(diào)洪計(jì)算所用，如梅亞東和談廣鳴[5]采用蒙特卡羅模擬法計(jì)算大壩防洪安全綜合風(fēng)險(xiǎn)率；呂彎彎等[6]運(yùn)用蒙特卡羅法計(jì)算土石壩洪水漫頂風(fēng)險(xiǎn)率的模型，分析洪水、風(fēng)浪、水庫運(yùn)行調(diào)度方式和水庫特性對洪水漫頂風(fēng)險(xiǎn)率的影響。

傳統(tǒng)的大壩洪水漫頂風(fēng)險(xiǎn)率計(jì)算方法通常只對洪峰流量進(jìn)行單變量頻率分析經(jīng)調(diào)洪演算計(jì)算漫壩風(fēng)險(xiǎn)率。后來，越來越多學(xué)者認(rèn)為洪水是一個(gè)多變量的隨機(jī)事件，研究了同時(shí)考慮洪峰洪量的兩變量洪水頻率分析的大壩洪水漫頂風(fēng)險(xiǎn)率[7]，表明兩變量洪水頻率分析更合理，以上都是僅基于單一洪水作用下的洪水漫壩風(fēng)險(xiǎn)率研究。然而隨著我國經(jīng)濟(jì)的快速發(fā)展，越來越多沿江城市在穿越河流的支流入?yún)R處修建了大量的水利工程，這些工程同時(shí)起著抵御外江洪水倒灌和保證城市排澇及生態(tài)需水等多重作用。對于這些干支流洪水作用下的水利工程而言，漫頂風(fēng)險(xiǎn)率的計(jì)算需要考慮干支流洪水的共同作用。目前國內(nèi)外研究較多的是入庫洪水、泄流能力、水庫調(diào)度規(guī)則等對漫頂風(fēng)險(xiǎn)率計(jì)算的影響，還未見對于同時(shí)考慮干支流洪水對漫頂風(fēng)險(xiǎn)影響的研究。為此，本文通過引入Copula函數(shù)，開展基于Copula函數(shù)的干支流洪水作用下的水利工程洪水漫頂風(fēng)險(xiǎn)率的研究，分析不同Copula函數(shù)以及不同干支流洪水相關(guān)性對漫頂風(fēng)險(xiǎn)率計(jì)算的影響。

1 干支流洪水作用下的水利工程洪水漫頂風(fēng)險(xiǎn)率計(jì)算模型

1.1 大壩洪水漫頂風(fēng)險(xiǎn)分析模型

洪水漫頂風(fēng)險(xiǎn)率是壩前最高水位Z超過壩頂高程Zd發(fā)生的概率，可表示為

Pf=P(Z>Zd)

(1)

一般而言，壩頂高程Zd隨機(jī)變動幅度變動不大，可將它視為確定值。壩前最高水位Z則存在著諸多不確定的因素，如入庫洪水過程的不確定性，泄流能力的不確定性等，因此壩前最高水位Z屬于隨機(jī)變量，漫頂率的計(jì)算可轉(zhuǎn)化為壩前最高水位Z的分布研究，一旦累積分布F(Z)確定，漫頂率就可以通過Pf=P(Z>Zd)=1-F(Z)計(jì)算求得。

由于干支流洪水作用下水利工程的設(shè)計(jì)洪水位不僅會受支流洪水的影響還會受到干流洪水頂托影響，利用壩址洪水進(jìn)行調(diào)洪演算反推壩前水位Z時(shí)需考慮干支流洪水流量間的相關(guān)性。因此，本文將采用Copula-Monte Carlo模擬法求解干支流洪水作用下水利工程的洪水漫頂風(fēng)險(xiǎn)率。

1.2 基于Copula的干支流洪水聯(lián)合分布

Copula函數(shù)可用來描述多變量水文要素間的相關(guān)性，可構(gòu)造出邊緣分布服從任意分布的聯(lián)合分布。假設(shè)X干和Y支分別表示干流和支流的洪水流量，則二維Copula函數(shù)為

F(x干，y支)=Cθ(F(x干)，F(xiàn)(y支))=Cθ(u，v)

(2)

式中，F(xiàn)(x干，y支)為干流和支流洪水流量的聯(lián)合分布，u=F(x干)，v=F(y支)分別為邊緣分布函數(shù)。

本文選取5種常用的阿基米德Copula函數(shù)，即GH、Frank 、AMH、Clayton和CClayton Copula函數(shù)來描述干支流洪水間相關(guān)結(jié)構(gòu)，其中，GH和CClayton Copula均具有上尾相關(guān)性，相關(guān)系數(shù)分別為2-2-1/θ和2-1/θ，Clayton Copula具有下尾相關(guān)性，F(xiàn)rank和AMH Copula函數(shù)則是尾部漸近獨(dú)立的。采用K-S檢驗(yàn)法對五種Copula函數(shù)進(jìn)行假設(shè)檢驗(yàn)，利用圖形分析法、AIC信息準(zhǔn)則法、均方根誤差法進(jìn)行優(yōu)選。

1.3 Copula-Monte Carlo方法計(jì)算漫頂率

漫頂風(fēng)險(xiǎn)率的核心內(nèi)容是獲取壩上水位的分布。本文采用的Copula-Monte Carlo方法的模擬思路為：

(1)根據(jù)建立的干支流洪水聯(lián)合分布函數(shù)，可以得到當(dāng)x干為指定值時(shí)Y支的條件分布，數(shù)學(xué)表達(dá)式為

SX干(y支|X干=x干)=?F(x干，y支)/?x干

(3)

(2)產(chǎn)生服從[0，1]均勻分布的兩個(gè)獨(dú)立的隨機(jī)數(shù)n1和n2。

(3)令n1為x干發(fā)生的不超過概率F(x干)，根據(jù)其反函數(shù)得到x干，即

x干=F-1(u)

(4)

(4)令n2為x干發(fā)生時(shí)y支的條件概率值，根據(jù)條件概率的反函數(shù)可得到y(tǒng)支，即

(5)

(5)由步驟(3)和(4)即可抽取一對干支流洪水組合流量(x干，y支)，將該組合洪水組合代入調(diào)洪模型進(jìn)行計(jì)算，可得到該組合對應(yīng)的壩上水位。

(6)重復(fù)上述步驟n次，可以模擬出n個(gè)壩上水位Z值。當(dāng)模擬的樣本容量足夠大，采用數(shù)學(xué)期望公式計(jì)算經(jīng)驗(yàn)頻率，獲得壩上水位的累積分布函數(shù)F(Z)，再利用公式Pf=P(Z>Zd)=1-F(Z)可得洪水漫頂率Pf。

2 實(shí)例應(yīng)用

2.1 工程概況

桂平航運(yùn)樞紐位于廣西桂平市郁江河段，是一個(gè)集航運(yùn)、發(fā)電、灌溉、交通于一體的綜合利用性航運(yùn)樞紐工程，處于郁江與潯江干流(包括黔江、潯江)匯合口附近的郁江段，主要由船閘、溢洪壩及水閘組成。桂平航運(yùn)樞紐同時(shí)受到郁江支流和潯江干流洪水的作用，正常擋水位29.6 m，死水位28.6 m，壩頂高程為47.8 m，工程等級為2級。

2.2 干支流洪水聯(lián)合分布構(gòu)建

潯江干流選擇大湟江口站1951年～2010年最大洪峰流量系列，郁江以大湟江口站發(fā)生最大洪峰流量的時(shí)間為控制，統(tǒng)計(jì)貴港站相應(yīng)場次洪水的相應(yīng)洪峰流量并進(jìn)行頻率計(jì)算及P-III型曲線適線得到邊緣分布，參數(shù)估計(jì)和檢驗(yàn)結(jié)果見表1。

表1 邊緣分布參數(shù)估計(jì)及檢驗(yàn)結(jié)果

計(jì)算兩變量間的Kendall相關(guān)系數(shù)τ為0.158 5，采用Kendall相關(guān)系數(shù)法計(jì)算5種備選Copula函數(shù)的參數(shù)θ，如表2所示。表2還給出了K-S檢驗(yàn)值，均方根誤差值，AIC值，經(jīng)驗(yàn)值與理論值的擬合效果如圖1所示。從表2可看出，5種備選Copula函數(shù)均通過了K-S檢驗(yàn)，但Clayton Copula函數(shù)的均方根誤差和AIC值最小，被選為最優(yōu)的Copula函數(shù)。

表2 五種候選Copula函數(shù)參數(shù)估計(jì)與擬合檢驗(yàn)結(jié)果

圖1 經(jīng)驗(yàn)頻率與理論頻率擬合

2.3 漫頂計(jì)算結(jié)果

根據(jù)1.3節(jié)中Copula-Monte Carlo方法的模擬思路，隨機(jī)生成50萬組干流流量、支流流量組合，最優(yōu)Copula隨機(jī)模擬的潯江、郁江洪水流量點(diǎn)據(jù)與實(shí)測點(diǎn)據(jù)的對比如圖2所示。查詢壩址水位與潯江流量、郁江流量的關(guān)系表，通過插值可得到50萬組洪水組合所對應(yīng)的壩址水位。根據(jù)桂平航運(yùn)樞紐的調(diào)度規(guī)則[9]進(jìn)行調(diào)蓄計(jì)算可得到對應(yīng)的50萬個(gè)壩上水位。最優(yōu)Copula函數(shù)的壩上水位的累積分布曲線見圖3。為了更好地展示結(jié)果，只給出了累積概率值大于0.95，即重現(xiàn)期20年一遇的上尾部分。由圖3可看出，桂平航運(yùn)樞紐壩頂高程47.8 m所對應(yīng)的累積分布概率F(Z)為0.999 862，相應(yīng)的漫頂率Pf=1-F(Z)=1.38×10-4，表明桂平航運(yùn)樞紐有著較高的防洪安全度。

圖2 最優(yōu)Copula蒙特卡羅模擬點(diǎn)與實(shí)測點(diǎn)比較

圖3 最優(yōu)Copula函數(shù)的壩上水位累積分布曲線

2.4 不同Copula函數(shù)對洪水漫頂影響

一般來說，不同的Copula函數(shù)描述干支流洪水兩變量間的相關(guān)關(guān)系是不一樣的，因此有必要對不同Copula函數(shù)進(jìn)行比較，分析不同Copula函數(shù)構(gòu)造聯(lián)合概率分布函數(shù)對洪水漫頂?shù)挠绊懀鐖D4所示。從圖4中可以看出，不同Copula函數(shù)對洪水漫頂?shù)挠绊懖⒉幻黠@。出現(xiàn)這種現(xiàn)象的主要原因有兩方面：①由于本研究區(qū)域的干支流洪水間相關(guān)性為0.158 5相對較小，由Copula函數(shù)理論可知，當(dāng)變量相對獨(dú)立時(shí)，不同Copula函數(shù)收斂于獨(dú)立Copula函數(shù)，此時(shí)不同Copula函數(shù)之間沒有差異。當(dāng)相關(guān)性增加時(shí)，不同Copula函數(shù)差別隨相關(guān)系數(shù)的增加先增大后減小。因此，不同Copula函數(shù)間差異與相關(guān)性有著一定聯(lián)系。也就是說，當(dāng)模擬的干支流洪水間相關(guān)性較小時(shí)，不同Copula函數(shù)對洪水漫頂?shù)挠绊懖粫蟆＂谒坏戎稻€接近與X軸垂直，略微向左傾斜(見圖5)，因此相比支流流量Q貴港而言，干流流量Q大湟對壩上水位的影響更大。然而，從表2中可以看出，大湟江年最大洪峰流量的變異系數(shù)Cv小于貴港相應(yīng)洪峰流量的Cv，從而導(dǎo)致不同Copula函數(shù)對洪水漫頂?shù)挠绊戄^小。

圖4 五種不同Copula函數(shù)壩上水位累積分布曲線

圖5 不同設(shè)計(jì)重現(xiàn)期對應(yīng)的壩上水位

由圖4中還可看出，GH、CClayton、AMH、Frank Copula均位于Clayton Copula的上方，會低估洪水風(fēng)險(xiǎn)。具有上尾相關(guān)性GH、CClaytonCopula函數(shù)計(jì)算的洪水漫頂率略高于上下尾部均為0的AMH、Frank Copula，可見尾部相關(guān)性對洪水漫頂率的計(jì)算也產(chǎn)生了一定影響。5種Copula函數(shù)計(jì)算的洪水漫頂率結(jié)果見圖6。由于大壩失事將會帶來十分嚴(yán)重的后果，工程上往往采用更保守的設(shè)計(jì)，因此采用最優(yōu)的Clayton Copula函數(shù)計(jì)算的洪水漫頂率更符合工程設(shè)計(jì)要求。

圖6 不同Copula函數(shù)計(jì)算的洪水漫頂率

2.5 不同干支流洪水相關(guān)性對洪水漫頂影響

一般而言，干支流洪水變量間往往存在著較弱的正相關(guān)性[10-11]，因此，本文主要討論τ=0.1，0.2，0.3，0.4，0.5等5種干支流洪水相關(guān)情況對洪水漫頂?shù)挠绊憽D7分別給出基于Clayton Copula函數(shù)的五種不同干支流洪水相關(guān)性下的壩上水位累積分布曲線，相應(yīng)的洪水漫頂風(fēng)險(xiǎn)率計(jì)算結(jié)果見圖8。由圖7可以看出，當(dāng)τ=0.1時(shí)，壩上水位累積分布曲線位于最上端，當(dāng)τ=0.5時(shí)，累積分布曲線位于最下端，也就是說隨著干支流洪水相關(guān)性的增強(qiáng)，漫頂率也在增加。因此在干支流洪水共同作用下的水利工程設(shè)計(jì)中，若忽略干支流洪水間的相關(guān)性，將會低估洪水漫頂風(fēng)險(xiǎn)，使得工程偏危險(xiǎn)。

圖7 不同干支流洪水相關(guān)性壩上水位累積分布曲線

圖8 不同干支流洪水相關(guān)性的洪水漫頂風(fēng)險(xiǎn)率

3 結(jié) 論

本文以桂平航運(yùn)樞紐為例，采用Copula-Monte Carlo模擬方法構(gòu)建干支流洪水作用下大壩的洪水漫頂率計(jì)算模型，并分析不同Copula函數(shù)、不同干支流洪水相關(guān)性對漫頂風(fēng)險(xiǎn)率的影響。得到以下結(jié)論：

(1)最優(yōu)的Clayton Copula計(jì)算的洪水漫頂率為1.38×10-4。

(2)由于干支流洪水間相關(guān)性較弱，不同Copula函數(shù)對洪水漫頂?shù)挠绊懖⒉缓苊黠@，GH、CClayton、AMH、FrankCopula均會低估洪水漫頂風(fēng)險(xiǎn)，從防洪安全角度考慮，采用最優(yōu)的Clayton Copula計(jì)算的洪水漫頂率更合理。

(3)隨著干支流洪水流量相關(guān)性的增強(qiáng)，洪水漫頂率也會相應(yīng)的增大，若忽略干支流洪水間的相關(guān)性，將會低估洪水漫頂風(fēng)險(xiǎn)，使得工程偏危險(xiǎn)。