劉章君，郭生練，許新發(fā)，成靜清，溫天福，尹家波

（1.江西省水利科學研究院，江西 南昌 330029；2.武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢 430072）

1 研究背景

洪水重現(xiàn)期與設(shè)計值的分析計算是水庫防洪安全設(shè)計及風險分析的關(guān)鍵問題[1］。單變量情況下，重現(xiàn)期和設(shè)計值存在唯一的對應(yīng)關(guān)系，計算方法清晰簡單。在兩變量情況下，由于涉及到洪峰、洪量的組合，因而重現(xiàn)期的定義和聯(lián)合設(shè)計值的計算都要比單變量情形復雜得多，近年來相關(guān)的研究受到了國內(nèi)外水文學者的高度關(guān)注[2-4］。

根據(jù)兩變量框架下危險區(qū)域的劃分方式不同，目前對兩變量洪水重現(xiàn)期的定義也相應(yīng)有很多種?！盎颉保ā癘R”）和“且”（“AND”）重現(xiàn)期是目前使用最廣泛的兩種定義方法，它們均根據(jù)某一特定洪峰、洪量聯(lián)合設(shè)計值定義危險區(qū)域，大量研究表明，在安全與危險事件的判定上兩者都存在著較大局限性[5-9］。為了解決上述問題，Salvadori等[5］根據(jù)聯(lián)合概率等值線定義危險區(qū)域，基于Kendall函數(shù)提出了兩變量洪水Kendall重現(xiàn)期，認為只有當洪峰、洪量聯(lián)合概率值超過了設(shè)定的閾值時才是危險事件。Kendall重現(xiàn)期使得洪峰、洪量聯(lián)合概率等值線上的任意一點都具有相同的危險區(qū)域，對安全事件與危險事件的判定更加合理[10-12］。此外，考慮到對于洪水這種對于越大越不利型變量，理論上安全區(qū)域應(yīng)該是有界的。針對Kendall重現(xiàn)期安全區(qū)域無界的問題，Salvadori等[13］在此基礎(chǔ)上又提出了兩變量洪水生存Kendall重現(xiàn)期，使得兩變量洪水重現(xiàn)期定義更加嚴謹。

事實上，就水庫防洪安全而言，最重要的因素是壩前最高水位，真正的危險事件應(yīng)該是壩前最高水位超過某一指定高程[14］。前述4種兩變量重現(xiàn)期所判定的危險事件和表征的水庫大壩水文失事情景都可能不符合實際情況。因此，有些學者更加突出強調(diào)水文荷載（如洪峰、洪量）與水工結(jié)構(gòu)的交互作用，提出了基于水工結(jié)構(gòu)安全的兩變量洪水重現(xiàn)期，并依此進行相關(guān)的洪水聯(lián)合設(shè)計和風險評估[15-17］。國內(nèi)外相關(guān)研究仍處于起步階段，尚未形成統(tǒng)一方法構(gòu)架，且大多數(shù)都進行了一定的概化或簡化處理，例如將入庫洪水過程形狀概化為三角形、水庫調(diào)洪演算簡化為線性水庫模型或自由泄流模型等。目前對于具有實際洪水過程線形狀和復雜調(diào)洪規(guī)則的情形在國內(nèi)外尚未見諸文獻報道，亟待深入研究。

本文將針對具有實際洪水過程線形狀和復雜調(diào)洪規(guī)則的情形，以超過壩前最高水位的頻率來度量防洪安全設(shè)計標準，通過壩前最高水位等值線劃分危險區(qū)域，開展兩變量洪水結(jié)構(gòu)荷載重現(xiàn)期（TSL）與最可能聯(lián)合設(shè)計值研究。以清江流域隔河巖水庫為例進行應(yīng)用，并與常用的“OR”重現(xiàn)期（TOR）、“AND”重現(xiàn)期（TAND）、Kendall重現(xiàn)期（TK）和生存Kendall重現(xiàn)期（TSK）等4種兩變量重現(xiàn)期最可能聯(lián)合設(shè)計值進行比較，以期為水庫防洪安全設(shè)計和風險分析提供更加科學合理的依據(jù)。

2 兩變量洪水結(jié)構(gòu)荷載重現(xiàn)期計算方法

2.1 基于Copula函數(shù)的洪水峰量聯(lián)合分布假設(shè)Q、W分別表示年最大洪峰流量和時段洪量，對應(yīng)的設(shè)計值分別為q、w，其邊緣分布分別為由Sklar定理知，Q和W的聯(lián)合分布函數(shù)可以用一個二維 Copula函數(shù) C表示[18-19］：

式中：F（q，w）為Q和W的聯(lián)合分布函數(shù)；θ為Copula函數(shù)的參數(shù)。

一般而言，洪水峰、量之間存在較強的正相關(guān)性，通常采用Gumbel-Hougaard Copula函數(shù)來構(gòu)造二者的聯(lián)合分布，其表達式如下[20-22］：

2.2 兩變量洪水結(jié)構(gòu)荷載重現(xiàn)期基本原理重現(xiàn)期是指給定事件連續(xù)兩次發(fā)生的平均時間間隔長度。對于洪水重現(xiàn)期而言，通常定義為洪水危險事件發(fā)生概率的倒數(shù)[17］。就水庫防洪安全而言，最重要的因素是壩前最高水位，真正的危險事件應(yīng)該是壩前最高水位超過某一指定高程。因此，指定防洪標準T年一遇設(shè)計洪水的實質(zhì)內(nèi)涵是以此為據(jù)而規(guī)劃設(shè)計出的水庫工程，其防洪安全事故的風險率應(yīng)恰好等于指定的設(shè)計標準（1/T）[23］。也即兩變量洪水重現(xiàn)期可以通過壩前最高水位的單變量重現(xiàn)期來間接度量。采用這種方式定義的兩變量洪水重現(xiàn)期稱之為結(jié)構(gòu)荷載重現(xiàn)期（TSL），它考慮了水文荷載（如洪峰、洪量）與水工結(jié)構(gòu)的交互作用，表征的水庫大壩水文失事情景更符合實際。

壩前最高水位Z與兩變量水文荷載洪峰Q、洪量W的交互作用可以通過壩前最高水位函數(shù)z=g(q，w)來表征。根據(jù)壩前最高水位等值線{g(q，w)=z}定義危險區(qū)域，認為只有當洪峰、洪量組合對應(yīng)的壩前最高水位值超過設(shè)定的閾值時才是危險事件。

對應(yīng)的結(jié)構(gòu)荷載重現(xiàn)期計算公式為：

2.3 基于Copula-MC法的結(jié)構(gòu)荷載重現(xiàn)期計算步驟式（5）表明結(jié)構(gòu)荷載重現(xiàn)期及設(shè)計值計算的關(guān)鍵問題是確定聯(lián)合概率密度函數(shù)和壩前最高水位函數(shù)z=g(q ，w )。洪峰、洪量的邊緣分布通常是非正態(tài)的且兩者存在一定的相關(guān)性，另外在具有實際洪水過程線形狀和復雜調(diào)洪規(guī)則的情形下積分區(qū)域邊界函數(shù)通常具有非線性及分段特征，這些復雜因素都導致形如式（5）的多重積分方法在實際應(yīng)用中非常困難。本文采用基于Copula函數(shù)的Monte Carlo（Copula-MC）法進行求解，具體步驟如下：

（1）根據(jù)年極值取樣數(shù)據(jù)建立的洪峰洪量聯(lián)合分布Cθ(u，v)，可以得到當Q為指定值q時W的條件分布，數(shù)學表達式為[20-21］：

（2）產(chǎn)生服從[0，1］均勻分布的兩個獨立隨機數(shù)r1和r2；

（3）令 u=r1， SU(v|U=u)=r2；

（6）選取有代表性的典型洪水過程，采用變倍比放大方法[24］獲得設(shè)計洪水過程線：

式中：DFH（t）、TFH（t）分別為設(shè)計洪水過程和典型洪水過程在t時刻的流量；QD、WD分別為典型洪水的洪峰流量和洪水歷時DT內(nèi)的洪量；q、w表示步驟5中隨機生成的洪峰、洪量。

這種方法可以保證放大后的洪峰洪量都等于指定的設(shè)計值，而且可以較好地保持典型洪水過程線的形狀；

（7）將放大后的設(shè)計洪水過程線DFH（t）輸入水庫按調(diào)洪規(guī)則進行調(diào)洪演算，得到壩前最高水位z值；

（8）重復步驟（2）—（7）共n次，可以模擬出n個 z值，采用數(shù)學期望公式計算經(jīng)驗頻率，獲得壩前最高水位Z的概率分布曲線；

（9）對于任意給定的洪峰洪量組合（q，w），依據(jù)步驟（6）—（7）得到相應(yīng)的壩前最高水位 z值，再根據(jù)壩前最高水位Z的概率分布曲線就可以查詢得到對應(yīng)的兩變量洪水結(jié)構(gòu)荷載重現(xiàn)期TSL；或者給定防洪標準T，令結(jié)構(gòu)荷載重現(xiàn)期TSL等于T，根據(jù)壩前最高水位Z的概率分布曲線查詢得到對應(yīng)的壩前最高水位設(shè)計值ZT。

3 5種兩變量洪水重現(xiàn)期比較及最可能聯(lián)合設(shè)計值

3.1 兩變量洪水重現(xiàn)期比較分析對目前常用的“OR”、“AND”、Kendall和生存Kendall等4種以及結(jié)構(gòu)荷載（共5種）兩變量洪水重現(xiàn)期的概念、危險區(qū)域劃分依據(jù)和同一防洪標準對應(yīng)的危險區(qū)域的唯一性進行總結(jié)歸納[6-9］，結(jié)果列于表1。

表1 5種兩變量洪水重現(xiàn)期的概念、危險區(qū)域劃分依據(jù)和唯一性

表2給出了5種兩變量洪水重現(xiàn)期對應(yīng)的危險事件數(shù)學描述和重現(xiàn)期計算公式[9，16］。其中，洪峰、洪量的單變量生存函數(shù)生存聯(lián)合概率分布可以表示為：

從表1和表2可以看出，“OR”和“AND”重現(xiàn)期均根據(jù)某一特定洪峰、洪量聯(lián)合設(shè)計值定義危險區(qū)域，同一防洪標準對應(yīng)的危險區(qū)域是不唯一的，大量研究表明，在安全與危險事件的判定上兩者都存在著較大局限性[5-9］。與“OR”和“AND”重現(xiàn)期相比，Kendall、生存Kendall和結(jié)構(gòu)荷載重現(xiàn)期的核心思想都是分別通過等值線和將多維信息組合轉(zhuǎn)化為單維信息，將兩變量洪水重現(xiàn)期表示為轉(zhuǎn)化變量的單變量重現(xiàn)期，使得同一防洪標準對應(yīng)的危險區(qū)域是唯一的，對安全事件與危險事件的判定在理論上更加協(xié)調(diào)一致。生存Kendall重現(xiàn)期可以解決Kendall重現(xiàn)期安全區(qū)域無界的問題[13］。結(jié)構(gòu)荷載重現(xiàn)期采用壩前最高水位等值線劃定危險區(qū)域，比Kendall用聯(lián)合概率和生存Kendall用生存聯(lián)合概率等值線更有物理基礎(chǔ)，表征的水庫大壩水文失事情景更符合實際。

3.2 最可能聯(lián)合設(shè)計值當確定了兩變量洪水重現(xiàn)期的定義方式后，一般通過假定兩變量洪水重現(xiàn)期等于指定防洪標準T，根據(jù)一定的準則計算得到兩變量洪水聯(lián)合設(shè)計值。然而，對于給定的洪水兩變量重現(xiàn)期，存在無窮多種滿足設(shè)計標準的洪峰、洪量組合，它們構(gòu)成了一條等值線，如何根據(jù)一定的準則科學合理地選擇聯(lián)合設(shè)計值成為另一個關(guān)鍵問題。目前洪水聯(lián)合設(shè)計值選擇方法主要有最可能組合法、超越概率最大組合法、條件期望組合法等，其中尤以最可能組合法得到廣泛采用[6-7］。最可能聯(lián)合設(shè)計值是指滿足指定防洪標準T條件下，洪峰、洪量聯(lián)合概率密度最大對應(yīng)的聯(lián)合設(shè)計值。

組合（qm，wm）即為推求的兩變量洪峰、洪量最可能聯(lián)合設(shè)計值。

綜上，本文采用的研究方案流程如圖1所示。

4 實例應(yīng)用

如圖2所示，隔河巖水庫位于清江流域下游，壩址以上流域面積14 430 km2，具有年調(diào)節(jié)能力。該斷面單峰歷時一般為3～5 d，復峰可達10 d，在推求設(shè)計洪水過程線或進行水庫防洪安全設(shè)計時，一般選取歷時為7 d的典型洪水過程[25］。因此，本文將分析洪峰與7 d洪量的聯(lián)合分布，推求洪峰與7 d洪量的兩變量洪水重現(xiàn)期與聯(lián)合設(shè)計值。隔河巖水庫汛限水位為192.2 m，防洪庫容5億m3。

圖1 研究方案流程

圖2 清江流域隔河巖水庫示意

4.1 邊緣分布及聯(lián)合分布的確定隔河巖水庫洪峰和7 d洪量的設(shè)計參數(shù)及設(shè)計值見表3。采用Gum?bel-Hougaard Copula建立洪峰與7 d洪量的二維聯(lián)合分布，基于1958—2010年同步資料計算得到Ken?dall相關(guān)系數(shù)為0.663，Copula函數(shù)的參數(shù)估計結(jié)果為θ=2.97。Cramér-von Mises檢驗統(tǒng)計量Sn=0.0232，對應(yīng)的p值為0.8267大于0.05，表明建立的Gumbel-Hougaard Copula函數(shù)通過了統(tǒng)計檢驗。圖3對聯(lián)合觀測變量的經(jīng)驗聯(lián)合分布值與理論聯(lián)合分布值進行了對比，為直觀起見，按理論頻率的升序排列，經(jīng)驗頻率與理論頻率值的擬合情況很好，表明所建立的聯(lián)合分布是合理可行的。

表3 隔河巖水庫設(shè)計洪水統(tǒng)計參數(shù)和設(shè)計值

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

4.2 隔河巖水庫壩前最高水位的概率分布按2.3節(jié)中的方法步驟，利用建立的聯(lián)合分布和條件分布隨機生成10萬組洪峰、洪量組合，隨機模擬的洪峰、洪量點據(jù)與實測值的對比見圖4。由圖4可以發(fā)現(xiàn)，Gumbel-Hougaard Copula函數(shù)能很好地捕捉洪峰、洪量的相關(guān)性結(jié)構(gòu)，可以用來隨機生成洪峰、洪量組合用于水庫防洪安全設(shè)計。選擇峰高量大、主峰靠后的1997年實測洪水為代表性典型洪水過程，放大得到10萬條設(shè)計洪水過程線，并輸入到隔河巖水庫進行調(diào)洪演算，從汛限水位開始起調(diào)，得到10萬個對應(yīng)的壩前最高水位值。基于數(shù)學期望公式獲得的壩前最高水位經(jīng)驗頻率曲線（由于樣本數(shù)量足夠大可以近似認為是概率分布曲線），見圖5。為了更好地展示結(jié)果，只給出了頻率曲線的上尾部（累積概率值大于0.80，即重現(xiàn)期5年）。

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

圖5 壩前最高水位概率分布曲線

4.3 隔河巖水庫結(jié)構(gòu)荷載重現(xiàn)期及設(shè)計值計算本文以防洪標準T=10、100、200和1000年為例，根據(jù)圖5壩前最高水位概率分布曲線就可以查詢得到對應(yīng)的壩前最高水位設(shè)計值ZT分別為192.45、195.86、197.19和202.88 m。為了得到滿足ZT的洪峰、洪量組合，結(jié)合洪峰、洪量邊緣概率分布的特點，按洪峰[10 100，27 800］m3/s以 100 m3/s等間距離散，洪量[2200，7800］×106m3以 50×106m3等間距離散，共得到20 114組洪峰、洪量組合，按1997年典型洪水過程線放大得到設(shè)計洪水過程線，輸入到隔河巖水庫進行調(diào)洪演算得到對應(yīng)的壩前最高水位值。從壩前最高水位系列中分別統(tǒng)計水位192.45、195.86、197.19和202.88 m對應(yīng)的洪峰、洪量組合，圖6是給定壩前最高水位條件下洪量與洪峰的關(guān)系散點圖?？梢园l(fā)現(xiàn)，洪量與洪峰存在顯著的線性或分段線性關(guān)系，對散點進行擬合回歸得到洪量與洪峰的函數(shù)關(guān)系w=h（q），具體數(shù)學表達式如圖6所示。

圖6清晰地揭示了水庫防洪安全設(shè)計中洪水過程線洪峰、洪量存在的“異組同效”現(xiàn)象，即不同洪峰、洪量組合得到的壩前最高水位相同，具有等效性?；貧w直線的斜率值為負數(shù)表明對于給定的壩前最高水位，洪峰、洪量呈現(xiàn)負相關(guān)關(guān)系，洪峰越大，對應(yīng)的洪量必然越小，反之亦然。如果斜率絕對值趨近于0表明調(diào)洪過程完全受洪量控制，洪峰不起作用，斜率絕對值趨近于+∞表明調(diào)洪過程完全受洪峰控制，洪量不起作用。圖6中斜率絕對值均大于0小于+∞，表明調(diào)洪過程同時受洪峰洪量控制，兩者都起作用。進一步分析發(fā)現(xiàn)，如果斜率絕對值小于1表明洪量起主要作用，洪峰起次要作用，反之大于1則表明洪峰起主要作用，洪量起次要作用。對照圖6的結(jié)果，T=10年洪峰起主要作用，洪量起次要作用；T=1000年洪量起主要作用，洪峰起次要作用；而T=100年和T=200年還表現(xiàn)出了分段特性，隨著洪峰流量的增大，由洪量起主要作用過渡到洪峰起主要作用。

結(jié)合洪量與洪峰的函數(shù)關(guān)系w=h（q），得到不同重現(xiàn)期結(jié)構(gòu)荷載重現(xiàn)期的最可能聯(lián)合設(shè)計值，列于表4。

圖6 給定壩前最高水位條件下洪量與洪峰的關(guān)系

表4 隔河巖水庫5種不同重現(xiàn)期洪峰、洪量聯(lián)合設(shè)計值（單位：洪峰：m3/s；洪量：106m3）

4.4 洪水重現(xiàn)期及設(shè)計值對比研究圖7給出了“OR”重現(xiàn)期、“AND”重現(xiàn)期、Kendall重現(xiàn)期（KEN）和生存Kendall重現(xiàn)期（SKEN）共4種常用的兩變量重現(xiàn)期與結(jié)構(gòu)荷載重現(xiàn)期（SL）洪峰洪量等值線圖的對比情況，其中黑色點據(jù)為實測值。可以看出，4種重現(xiàn)期等值線圖形狀與結(jié)構(gòu)荷載重現(xiàn)期均存在較大差異。其中，“OR”重現(xiàn)期等值線均位于結(jié)構(gòu)荷載重現(xiàn)期的上方，而“AND”重現(xiàn)期等值線均位于結(jié)構(gòu)荷載重現(xiàn)期的下方，表明這兩種重現(xiàn)期都存在系統(tǒng)性偏差，理論上存在不足。Kendall重現(xiàn)期和生存Kendall重現(xiàn)期等值線則與結(jié)構(gòu)荷載重現(xiàn)期存在交叉，相對“OR”和“AND”重現(xiàn)期更加合理。

采用3.2節(jié)中方法計算得到“OR”重現(xiàn)期、“AND”重現(xiàn)期、Kendall重現(xiàn)期和生存Kendall重現(xiàn)期的最可能聯(lián)合設(shè)計值，結(jié)果列于表4中。將得到的聯(lián)合設(shè)計值組合，按1997年典型洪水過程線放大得到設(shè)計洪水過程線，輸入到隔河巖水庫進行調(diào)洪演算得到對應(yīng)的壩前最高水位值，結(jié)果見表5。根據(jù)圖5壩前最高水位概率分布曲線查詢得到的不同重現(xiàn)期壩前最高水位對應(yīng)的設(shè)計標準見表6。

從表5和表6可知，由于結(jié)構(gòu)荷載重現(xiàn)期考慮了洪峰、洪量與水庫調(diào)洪規(guī)則的交互作用，其定義是基于設(shè)計洪水的實質(zhì)內(nèi)涵，表征的水庫大壩水文失事情景符合實際，按此進行設(shè)計的水庫工程正好達到指定防洪標準。與結(jié)構(gòu)荷載重現(xiàn)期相比，“OR”重現(xiàn)期設(shè)計標準系統(tǒng)偏高，“AND”重現(xiàn)期系統(tǒng)偏低。Kendall重現(xiàn)期最可能聯(lián)合設(shè)計值在T=10、100、200年設(shè)計標準偏低，而T=1000年偏高。生存Kendall重現(xiàn)期最可能聯(lián)合設(shè)計值在T=10年設(shè)計標準略微偏高，T=100、200年略微偏低，而T=1000年偏高。

圖7 4種常用兩變量重現(xiàn)期與結(jié)構(gòu)荷載重現(xiàn)期洪峰洪量等值線對比

表5 隔河巖水庫5種不同重現(xiàn)期壩前最高水位設(shè)計值 （單位：m）

表6 隔河巖水庫5種不同重現(xiàn)期壩前最高水位對應(yīng)的設(shè)計標準 （單位：年）

5 結(jié)論

本文以超過壩前最高水位的頻率來度量防洪安全設(shè)計標準，開展了具有實際洪水過程線形狀和復雜調(diào)洪規(guī)則情形下的兩變量洪水結(jié)構(gòu)荷載重現(xiàn)期與最可能聯(lián)合設(shè)計值計算方法研究。以隔河巖水庫為例進行應(yīng)用，并與常用的“OR”、“AND”、Kendall和生存Kendall兩變量重現(xiàn)期最可能聯(lián)合設(shè)計值進行比較，主要研究結(jié)論如下：（1）兩變量洪水結(jié)構(gòu)荷載重現(xiàn)期的定義考慮了洪峰、洪量與水庫調(diào)洪規(guī)則的交互作用，符合設(shè)計洪水的實質(zhì)內(nèi)涵，能夠達到指定的防洪標準。以壩前最高水位等值線定義危險區(qū)域，表征的大壩水文失事情景更符合實際，可為水庫防洪安全設(shè)計提供科學依據(jù)。（2）在給定壩前最高水位條件下，隔河巖水庫的洪量與洪峰存在顯著的線性或分段線性關(guān)系，調(diào)洪過程同時受洪峰洪量控制，并表現(xiàn)出明顯的“異組同效”現(xiàn)象，即在水庫防洪安全設(shè)計中不同洪峰、洪量組合得到的壩前最高水位相同。（3）本文隔河巖水庫實例中，“OR”和“AND”重現(xiàn)期最可能聯(lián)合設(shè)計值與結(jié)構(gòu)荷載重現(xiàn)期相比設(shè)計標準系統(tǒng)偏高或偏低，Kendall和生存Kendall重現(xiàn)期最可能聯(lián)合設(shè)計值設(shè)計標準也存在不同程度地偏低或偏高的情況。

然而，本文所提方法在推求洪峰、洪量最可能聯(lián)合設(shè)計值時是先確定洪峰，還是先確定洪量，順序不同結(jié)果可能存在差異，有待進一步研究。另外，本文只選擇單次典型洪水過程進行放大，所得結(jié)果具有一定的局限性，因此不同典型洪水過程對設(shè)計結(jié)果有何影響也需要進一步探討。