甘富萬，張華國，覃麗娜，高 揚，黃宇明，肖 良

(1.廣西大學土木建筑工程學院，廣西 南寧 530004；2.廣西大學工程防災與結構安全教育部重點實驗室，廣西 南寧 530004)

0 引 言

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

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

1 干支流洪水作用下的水利工程洪水漫頂風險率計算模型

1.1 大壩洪水漫頂風險分析模型

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

Pf=P(Z>Zd)

(1)

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

由于干支流洪水作用下水利工程的設計洪水位不僅會受支流洪水的影響還會受到干流洪水頂托影響，利用壩址洪水進行調洪演算反推壩前水位Z時需考慮干支流洪水流量間的相關性。因此，本文將采用Copula-Monte Carlo模擬法求解干支流洪水作用下水利工程的洪水漫頂風險率。

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

Copula函數(shù)可用來描述多變量水文要素間的相關性，可構造出邊緣分布服從任意分布的聯(lián)合分布。假設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ù)來描述干支流洪水間相關結構，其中，GH和CClayton Copula均具有上尾相關性，相關系數(shù)分別為2-2-1/θ和2-1/θ，Clayton Copula具有下尾相關性，F(xiàn)rank和AMH Copula函數(shù)則是尾部漸近獨立的。采用K-S檢驗法對五種Copula函數(shù)進行假設檢驗，利用圖形分析法、AIC信息準則法、均方根誤差法進行優(yōu)選。

1.3 Copula-Monte Carlo方法計算漫頂率

漫頂風險率的核心內容是獲取壩上水位的分布。本文采用的Copula-Monte Carlo方法的模擬思路為：

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

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

(3)

(2)產生服從[0，1]均勻分布的兩個獨立的隨機數(shù)n1和n2。

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

x干=F-1(u)

(4)

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

(5)

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

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

2 實例應用

2.1 工程概況

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

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

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

表1 邊緣分布參數(shù)估計及檢驗結果

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

表2 五種候選Copula函數(shù)參數(shù)估計與擬合檢驗結果

圖1 經驗頻率與理論頻率擬合

2.3 漫頂計算結果

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

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

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

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

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

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

圖5 不同設計重現(xiàn)期對應的壩上水位

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

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

2.5 不同干支流洪水相關性對洪水漫頂影響

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

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

圖8 不同干支流洪水相關性的洪水漫頂風險率

3 結 論

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

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

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

(3)隨著干支流洪水流量相關性的增強，洪水漫頂率也會相應的增大，若忽略干支流洪水間的相關性，將會低估洪水漫頂風險，使得工程偏危險。