基于概化洪水過程的水庫防洪安全計算

姚瑞虎，覃光華,2，丁 晶,2，曹玲然，李品良

(1.四川大學水利水電學院，四川 成都 610065；2.四川大學 水力學與山區河流開發保護國家重點室，四川 成都 610065)

洪水過程作為防洪安全計算的依據，當前主要采用同頻率法，盡管這樣一種做法已被列入洪水規范[5]，但一直也存在各種質疑，其核心的疑點在于計算結果是否達到水庫防洪安全標準要求。對這一疑點，已先后從各個角度做了一些探討[1]，但迄今仍未取得實質性進展，紛爭依舊存在。本文著眼于概化洪水過程的水庫防洪安全計算，并以此為基礎研究同頻率洪水過程的計算結果是否達標。

水庫工程的防洪安全取決于水庫的防洪最高水位Hm[1]，在一定的泄洪能力和調洪方案下，Hm與入庫的洪水緊密相關。洪水過程由洪峰(x)、洪量(y)和過程線的形狀(s)三要素組成[2]，所以Hm同x、y、s有關。因此，對水庫防洪安全計算和評價，本質上取決于Hm，而實際上取決于x、y、s的匹配。一般而言，對不同情況的水庫常常采用不同的匹配組合方式。對庫容較小的工程，Hm主要取決于x，即x是最為關鍵的因素，對安全起“頂梁柱”的作用；對庫容很大的工程，Hm主要取決于y和s。對多數水庫工程，Hm受x、y、s3個要素的制約。s涉及峰形、主峰位置、漲落速率等，為一個非常復雜的要素。為了便于分析研究，著重討論特殊的單峰情況，即討論一種單峰概化洪水過程[3]。在這樣的前提下深入分析Hm同x、y二者的關系。將復雜的洪水過程線概化成三角形為本文研究的前提。在水文分析計算領域，抓住研究對象特點將復雜問題簡化為單一問題，使分析計算易于進行并獲得實用意義成果，這是行之有效的權宜方法。將洪水過程概化成三角形過程，國內外早有相關研究[2,13-15]。當然，獲得的成果存在一定的局限性。不過單峰為我國多數河流發生的洪水所呈現出的一種形狀[4]，對這種峰形的研究有一定的實用意義。

本文基于概化洪水過程，在一定的泄洪能力和調洪方式下推求出Hm與x、y的關系式。據此，探討x和y如何匹配才能達到水庫工程防洪計算的最終目的，即達到國家規定的防洪安全標準[5]。

1 概化洪水過程和水庫調洪

將曲線構成的洪水過程概化成圖1所示的三角形單峰過程。

圖1 三角形概化洪水過程及其調洪示意Fig.1 Generalized triangle flood hydrograph and its flood regulation

圖中：x為最大洪峰流量；qm為最大泄洪流量；t0為洪量時段；qs為洪水來臨時的基流量；qc為起調水位H2時相應的起調流量；t1為大于qc的洪水歷時；V′為起調后的調節庫容與水庫調洪緊密有關的兩條曲線：

(1)庫容曲線，以下式表示：

V=a1(H-H1)n1

(1)

式中：V為水庫需蓄水庫容；H是與V相應的庫水位；H1是庫底高程；a1和n1為參數，由實測資料估計。

(2)泄洪曲線，以下式表示：

q=a2(H-H2)n2+qc

(2)

式中：q為水庫下泄的流量；H為水庫的水位；a2和n2為參數，由確定的泄流曲線估計。

洪水期，水庫的管理運用按一定的規則運行，一般當來水達到流量qc時開始調節，其相應的水位稱起調水位為H2(qc對應的庫水位)。在調洪時，隨水位H的上升，下泄流量不斷增加，最終到達最大值qm，對應的庫水位達到最高值Hm[6]。洪水導致的Hm是制約水庫安危的最重要變量，是水庫防洪安全最關鍵的因素。Hm在以上論述的條件下，取決于三角形的高—洪峰x和面積—洪量y，即：

Hm=f(x,y)

(3)

式(3)定量形式的推求為本文研究的重點。

2 水庫最高洪水位Hm和入庫洪峰x，洪量y之間的定量關系式

(4)

當達到qm時，與之相應的庫水位為Hm，與Hm相應的庫容為Vm。由圖1得：

Vm=Vc+V′

(5)

式中：Vc為qc時H2相應的蓄水庫容；V′起調后的調節庫容。由圖1知，V′為1,2,4三點形成的面積減去1,3,4三點形成的面積。

(6)

將式(6)代入式(5)得：

(7)

考慮到Vm和Hm的關系，將式(1)代入式(7)得：

進一步推導得：

(8)

將式(8)的t1表達式代入式(4)后整理得：

由圖1可知，式(9)適用條件為x>qc，當x≤qc，由于來水較小，水庫無需進行調洪。式(9)為Hm與x,y的定量關系式，但難以得到像式(3)那樣的顯式?？傊?，一旦(x,y)給定便有唯一Hm與之對應。

式(9)為Hm和x,y之間的定量關系式，其中參數的確定(選定)的情況如下：①a1，n1和a2，n2由該水庫庫容曲線和泄流曲線定出；②qc為起調流量，由調洪規則定出；③qs為洪水來臨時的基流量，據水文資料或調查信息估計確定。在式(9)條件下，Hm取決于x和y的匹配，Hm由x和y間接推求。下面論述由x,y推求Hm的模擬法。

3 x、y的模擬及Hm的計算

Hm的統計特性取決于x和y的統計特性，即Hm的分布受x和y的二維聯合分布制約,這是通過隨機模擬x和y推求相應Hm的依據。模擬x和y的關鍵在于建立x和y的二維聯合分布F(x,y)以及x和y的一維分布F(x)和F(y)[7-9]。

3.1 F(x)和F(y)的建立

3.2 F(x, y)的建立

由x和y的對應實測資料，推求選定的Copula函數C(u,v)的參數θ，其中：

u=P(X

(10)

由此建立F(x,y)[10]。

3.3 x，y和Hm的隨機模擬

Hm的模擬取決于x、y模擬，且每對x和y有下列要求：①x必須服從F(x)分布；②y必須服從F(y)分布；③x和y必須服從F(x,y)分布。顯然，前兩個要求易于實現。因為u=1-F(x)，v=1-F(y)中的u和v服從(0～1)均勻分布。對于第三個要求，模擬情況較為復雜。F(x,y)的建立基于C(u,v)，因此滿足F(x,y)制約關系的一對x和y，可由滿足C(u,v)的一對u和v代替。如果C(u,v)服從分布類型已知，則很容易進行模擬。但C(u,v)所隸屬的分布未知，因此，要換一種思路進行模擬。

基于C(u,v)可以求出已知其中一個變量發生的條件下另一個變量發生的條件分布。例如，已知u發生條件下，v發生的條件分布為：

(11)

已知v發生條件下，u發生的條件分布為：

(12)

S(v|u)和S(u|v)均服從(0～1)均勻分布[11]。本文正是依據這一重要特性模擬符合上述要求的一對u和v。若是先模擬出u，則可以利用公式(11)，模擬出相應的v。類似地適用于模擬u。這樣，基于式(11)或式(12)模擬出的u和v符合C(u,v)的要求。這是因為S(v|u)和S(u|v)均源自C(u,v)，從而符合S(v|u)和S(u|v)條件的u和v必然符合C(u,v)分布的要求。因此，基于一對u和v，可模擬出一對x和y，所以x，y和Hm的隨機模擬步驟[12]如下。

(1)模擬出符合(0～1)均勻分布的隨機數R1，R2，令u=R1，S(u|v)=R2，將u=R1代入到S(u|v)=R2得S(R1|v)=R2，進而求得v。

(2)檢驗模擬出的u、S和v是否服從(0～1)均勻分布[11]。

(3)由模擬出的u和v，基于式(10)可模擬出xi和yi；

(4)由模擬出的xi和yi，據式(9)計算出相應的(Hm)i；

(5)重復以上步驟K次，可得K對(xi,yi)及其相應的(Hm)i(i=1, 2, 3,…,K)。

3.4 (Hm)i的統計計算[1]

在(Hm)i,i=1,2,3,…,K的基礎上，重新按大小的次序排列，以一般經驗頻率公式計算出對應的經驗頻率P，據Hm和對應P，點繪在頻率格紙上，得K個經驗點據[2]。由于K很大(例如K=2 000)，故可直接聯結圖上的點據得Hm的頻率曲線。

圖2 Hm頻率曲線示意Fig.2 Frequency curve of Hm

4 (Hm)N作為防洪安全計算依據

(Hm)N的意涵還可通過圖3作進一步論述。對應(Hm)N的x和y，據式(9)有很多組，記為(x1，y1)，(x2，y2)，…，(xh，yh)。將h組數據點在圖3上，可連成一條線(abc表示的一條線)稱之為臨界線。在此線以上的部分為風險區(S1)，在線以下的部分為安全區(S2)。

(13)

圖3 洪水風險示意Fig.3 Sketch map of Flood risk

5 Hm下的x和y之匹配

(Hm)N作為安全計算依據，簡明、形象和直觀，但它的推求非常復雜。本文在三角形概化條件下，基于x和y的Copula函數，以隨機模擬法推求出(Hm)N。實際上我國長期洪水計算事件未以(Hm)N作為安全計算的依據，而是以特定洪水過程線(洪峰、洪量和洪水過程形狀)代替(Hm)N作為計算依據，這樣做的基本理念是以特定洪水過程線經調洪獲得特定最高水位，可以當做(Hm)N。即以間接計算的特定最高水位為依據做防洪計算，其結果使工程達到安全標準。這里“特定”和“標準”緊密相連。本文基于過程線形狀概化為三角形，探討符合“特定”要求(即“達標”要求)時，洪水過程的洪峰x和洪量y如何匹配。

式(9)表明：對于一定的Hm，由相應的x和y與之匹配。x和y如何匹配達標的最根本條件是二者代入式(9)所得的Hm值必須等于(Hm)N。在該條件下有許多組x和y，那么實踐中如何確定一組匹配。如圖3所示，x和y的不利匹配造成風險。不利匹配包括多種情況：較大的x和較小的y；較小的x和較大的y。在確定x時，x要有一個上限(類似地在確定y時，y亦應有一個上限)。在當前防洪標準以洪水重現期N表示的情況下，作者認為這個上限宜定為重現期為N的x，即xN(類似地y的上限為yN)。x和y的最大值定為xN和yN是處理這一難題的權宜策略。如何科學而合理地確定尚需進一步探討。

在上述處理策略下，x的最大取值為xN，y的最大取值為yN。因此在(x,y)的多種匹配中，選擇兩個特殊的匹配：(xN，yx)和(xy，yN)。yx表示與xN相應的y，xy表示與yN相應的x。具體推求方法如下：

將(Hm)N和xN一起代入式(9)，解得yx；將(Hm)N和yN一起代入式(9)，解得xy。在式(9)的前提下，探討達標要求下x和y匹配情況，可以從中得到一些有益的啟示。例如在式(9)的前提下，就可能定量探求同頻率洪水作為防洪安全計算依據，其結果超標的問題。同頻率洪水的實質在于選用xN和yN的匹配?；谑?9)由xN和yN得最高水位(Hm)C。若(Hm)C>(Hm)N，則同頻率法超標；若(Hm)C<(Hm)N，則同頻率法未達標。由圖3可知，凡是 和 符合標準的匹配都在abc臨界線上，若以xN和yN匹配點繪在圖3上便會落在abc臨界線的上面。為清楚起見，特給出圖4。該圖顯示xN和yN的匹配點落在危險區S1內，即(Hm)C>(Hm)N表示超標。x，y的匹配既受x和y本身及其相關關系統計特性的影響，又受制于工程庫容和泄流特性等因素，這是非常復雜的問題。但是就工程的防洪安全計算而言，匹配至關重要。這是因為x，y的匹配制約了Hm(制約工程防洪安全的關鍵特征量)。x和y的匹配結果會存在相應的Hm與之對應。因此，站在工程防洪安全的后果角度，x和y的匹配所造成的工程防洪風險程度可合理地以與之相應的Hm的頻率來表征，本文將此頻率表示為x和y匹配事件所導致的后續事件(壩前最高水位Hm)發生的頻率，稱為導致頻率。顯然它和匹配事件發生的頻率有根本的差異，下面實例計算也說明了這一點。

圖4 xN和yN匹配顯示Fig.4 Display of matching xN with yN

6 實例分析

本文以四川官帽舟水庫工程的基本資料以及其下游馬邊水文站的洪水資料為依據，對概化洪水過程進行推求計算。官帽舟電站工程壩址位于距馬邊縣20 km的蘇壩鎮上游約1 km處，集水面積1 449 km2,設計防洪標準P=1%。馬邊水文站位于官帽舟電站壩址下游21 km處，壩址以上集水面積占馬邊水文站以上面積(1 831 km2)的79.1%。所以，以馬邊河水文站作為官帽舟電站水文分析計算的代表站，馬邊站的洪峰x和洪量y分別以面積比的2/3次方和1次方換算到壩址處。具體步驟如下。

6.1 洪水統計模型及其參數估計

6.1.1 峰和量頻率分布和參數估計

基于馬邊站1957-2013年共57 a的年最大洪峰流量x及其相應的24 h洪量y洪水資料，選用P-Ⅲ型分布函數，以洪水規范規定的方法估計參數得F(x)和F(y)，如圖5，圖6所示。

圖5 洪峰頻率曲線Fig.5 Frequency curve of flood peak discharge

圖6 洪量頻率曲線Fig.6 Frequency curve of flood volume

6.1.2 峰和量二維聯合分布選用和參數估計

當前，表征x和y二維聯合分布函數常選用Copula函數作為其連接函數。該函數有多種類型[16]。因此，對x和y的相關特性進行分析，最終選擇與x和y相關結構十分相似的Clayton Copula函數，作為x和y的二維分布函數，由式(10)可知，連接函數形式為：

C(u,v)=(u-θ-v-θ-1]

(14)

那么，以u為條件，滿足Clayton Copula函數的條件v的條件分布為：

S(v/u)=u-θ-1(u-θ+v-θ-1)

(15)

式中:u和v意義同上；θ為參數，與Kendall秩相關系數τ的關系為：

(16)

基于x和y的二維分布函數統計特性可得，相應的二維聯合分布函數F(X>x,Y>y)，即：

F(x,y)=1-u-v+C(u,v)

(17)

將式(10)和式(14)帶入式(17)得：

F(x,y)=F(x)+F(y)-1-

[(1-F(x))-θ+(1-F(y))-θ]

(18)

式(18)中的F(x)和F(y)為x和y的頻率分布函數。式中參數θ可通過τ計算得出，τ的計算式為：

(19)

式中：n是變量x和y的樣本長度；xi、xj、yi、yj是樣本取值；sign(·)為符號函數。當(xi-xj)(yi-yj)>0時，sign(·)=1；當(xi-xj)(yi-yj)<0時，sign(·)=-1，否則sign(·)=0。基于馬邊站57 a年x和y樣本資料，估計得τ=0.59(可知θ=2.88)，這表明x和y是正相關關系并且關系較密切，這符合馬邊水文站洪水資料特征。

6.2 調洪基本參數的確定

由該工程的《蓄水安全鑒定設計自檢報告》獲得相關水庫調洪調度的相關參數值。公式(9)中：基流量qs=40 m3/s，起調流量qc為起調水位時所對應的下泄流量為1 004 m3/s，洪量時段t0=24 h，庫底高程H1=580 m，起調水位H2=665 m。根據該工程的《蓄水安全鑒定設計自檢報告》中所列的水庫庫容和泄流資料，通過適線法得出庫容曲線和泄流曲線方程：庫容曲線V=0.193 3×(H-580)2.36萬m3；泄流曲線q=137.75(H-665)1.24m3/s，因此，a1=0.199 3,，n1=2.36，a2=137.75，n2= 1.24。該工程最高庫水位Hm和入庫的洪峰x和洪量y的關系式，為：

0.193 3(Hm-580)2.36-6 913=

(20)

式(20)適用條件為x>1 004 m3/s，y>4 510 萬m3。

6.3 (Hm)N的計算

按照上述模擬法的步驟得到變量ui、Si和vi(i=1,2,…,K)，本次取K=2 000，進一步檢驗可知，ui和Si均服從(0～1)均勻分布，表明模擬過程正確可行。然后依據公式(10)和式(20)模擬出大量的(Hm)i(i=1,2,…,K)，從而獲得的Hm頻率點據分布如圖7所示。

圖7 Hm頻率曲線Fig.7 Frequency curve of Hm

從中挑選幾個較為典型的概率P(重現期N)，以及所對應的Hm列于表1。

表1 典型概率P所對應的Hm值Tab.1 The value of Hm corresponding to typical probability

由圖7和表1可知，當防洪標準P=1%，即重現期N=100 a時， (Hm)N=668.96 m。

6.4 x和y的匹配推求

基于本文論述的兩種匹配方式：(xN，yx)以及在(xy，yN)和(Hm)N=668.96 m的條件下，利用公式(20)得出：xN=2 320 m3/s，yx=7 720 萬m3；xy=1 920 m3/s，yN=13900 萬m3。由xy=1 920 m3/s，可得F(xy)=4.18%，由yx=7 720 萬m3，可得F(yx)=15.7%。此外，將xN=2 320 m3/s和yN=13 900 萬m3代入到公式(20)得，(Hm)C=670.20 m，根據圖7，查得相應的Pc=0.223%(重現期為N=450 a)。結果表明：該水庫工程達到風險為1%的防洪標準時，須選用100 a一遇的洪峰和約6.4 a一遇的24 h洪量的匹配，或者選用約24 a一遇的洪峰和100 a一遇的24 h洪量的匹配作為管運值，若洪峰和洪量都采用100 a一遇的標準，則匹配結果將達到0.223%(重現期N=450 a)，其防洪安全度超過標準。必須指出，xN和yN匹配經式(20)計算得到的導致頻率Pc和事件(X>xN∩Y>yN)發生的概率Ps(X>xN∩Y>yN)有本質的區別。前者和工程防洪特性和工程調度規則等有關，是指工程以xN和yN匹配作為防洪計算的依據時，工程潛在的洪水風險。后者與工程毫無關系，不涉及工程的特性和工程的調洪調度規則等，是指事件(X>xN∩Y>yN)出現的可能性大小(頻率)。在本例中，由二維聯合分布F(x,y)計算得到，Ps(X>2 320∩Y>1 3900)= 0.038%，顯然Pc明顯大于Ps。

6.5 工程防洪特性對Hm的影響

當水庫調度方式、洪峰x和洪量y以及洪水過程形狀確定以后，影響水庫防洪最高水位Hm的主要因素為工程的防洪特性——泄流特性和庫容特性。本文通過合理的調整參數n1=2.26和n2=1.14，探究工程防洪特性對Hm的影響。計算結果表明調整參數后計算得到的Hm值較未調整之前均有所增大。防洪標準P=1%所對應的(H′m)N=669.85 m較(Hm)N=668.96 m增大0.89 m。合理性分析可知，當n1減小時，水庫的庫容變小，在相同來水情況下，較未調整n1之前水庫庫水位H將抬高。由圖1可知，當庫容變小后，最大泄流量qm變大，由于n2的減小使得水庫泄流能力減小，因此需要更高的水頭來滿足泄流量增加帶來的影響。因此，水位將產生較大幅度的增長，由此可知，上述對Hm的計算結果合理。

綜上所述，通過改變參數n1和n2，進一步計算可知，工程防洪特性改變對水庫工程的防洪有較大的影響。

7 結 語

本文計算分析論證的結果可以歸納為以下幾點：

(1)工程壩前年最高水位為制約工程防洪安全的關鍵特征量，既受洪水影響又與洪水調度有關，是綜合反映洪水和工程特性的隨機變量。該隨機變量的頻率直接度量工程的洪水風險，可以視為防洪標準。

(2)將洪水過程概化為三角形的基礎上，借助工程的庫容和泄流特性可以建立壩前最高水位Hm和相應洪峰x以及洪量y的定量概化關系式，這大大方便了工程的防洪安全計算分析。

(3)基于概化關系式，建議以隨機模擬法推求符合防洪標準的年最高水位。

(4)基于概化關系式可分析在達標要求下，洪峰和洪量的匹配情況。建議當前選用(xN，yx)和(xy，yN)的匹配作為工程防洪安全計算的依據。同頻率洪峰和洪量為一組特殊的匹配，其防洪計算結果超標。

(5)基于匹配頻率的新概念洪峰和洪量的匹配頻率Pc同事件(X>xN∩Y>yN)發生的概率Ps存在本質區別。前者與工程防洪特性等因素緊密相關，后者則毫無關系。

(6)工程防洪特性的改變，對工程防洪安全有較大的影響。

(7)在概化關系式下獲得的成果存在一定的局限性。不過單峰型為我國多數河流發生大洪水時所呈現出的一種形式，對這種峰形做三角形概化有一定的實用意義。

□

[1] 丁 晶, 鄧育仁, 侯 玉, 等. 水庫防洪安全設計時設計洪水過程線法適用性的探討[J]. 水科學進展, 1992,3(1):45-52.

[2] 詹道江, 徐向陽, 陳元芳. 工程水文學[M]. 北京:中國水利水電出版社, 2010.

[3] 楊榮富, 丁 晶. 單峰型洪水過程線的概化及隨機模擬[J].成都科技大學學報, 1990,(5):67-72.

[4] 崔振才, 翟國靜. 河流洪水過程線的隨機模擬及其在推求調洪庫容中的應用[J]. 重慶交通學院院報, 1990,9(1)：103-114.

[5] 中華人民共和國水利部. 水利水電工程設計洪水計算規范[M].北京: 中國水利水電出版社, 2006.

[6] 葉守澤. 水文水利計算[M]. 北京: 中國水利水電出版社, 2008.

[7] 肖 義, 郭生練, 劉 攀, 等. 基于Copula函數的設計洪水過程線方法[J]. 武漢大學學報(工學版), 2007,40(4):13-17.

[8] 張冬冬, 魯 帆, 嚴登華, 等. 基于Archimedean Copula函數的洪水多要素聯合概率分布研究[J]. 中國農村水利水電, 2015,(1):68-73.

[9] 蘇夏弈, 張 鑫, 王 云, 等. 基于SPI和Copula的湟水流域干旱趨勢研究[J]. 中國農村水利水電, 2016,(12):151-154.

[10] 陳 璐. Copula函數理論在多變量水文分析計算中的應用研究[M].武漢: 武漢大學出版社, 2013.

[11] 韋艷華, 張世英.Copula函數理論及其在金融分析上的應用[M]. 北京: 清華大學出版社, 2008.

[12] 丁 晶, 朱宏江. 對現行設計洪水過程線方法適用性的統計實驗研究[J]. 四川水力發電, 1990,(1):17-24.

[13] 廖 松, 王燕生, 王 路. 工程水文學[M]. 北京, 清華大學出版社, 1991.

[14] G Salvadori, C De Michele F. Durante. On the return period and design in a Multivariate framework [J].Hydrologyand Earth System Sciences, 2012,12:2 699-2 708.

[15] C De Michele. Bivariate Statistical Approach to Check Adequacy of Dam Spillway [J]. Journal of Hydrologic Engineering, 2005,10(1):50-57.

[16] 許月萍, 童楊斌, 富 強,等. 幾種Copula模擬不同歷時降雨量的影響分析[J]. 浙江大學學報(工學版),2009,43(6):1 107- 1 111.