劉章君，成靜清，溫天福，吳向東，賈 磊

(1.江西省水利科學研究院，江西南昌330029;2.水利部鄱陽湖水資源水生態(tài)環(huán)境研究中心，江西南昌330029;3.江西省鄱陽湖水資源與環(huán)境重點實驗室，江西南昌330029)

0 引 言

設(shè)計洪水的地區(qū)組成是大中型水利水電工程洪水計算的主要研究內(nèi)容之一，它是確定流域梯級開發(fā)和工程設(shè)計規(guī)模的重要依據(jù)[1]。擬定設(shè)計洪水的地區(qū)組成時，需要區(qū)間各種防洪標準的設(shè)計洪量或洪量頻率曲線，并且區(qū)間洪水成果對下游斷面防洪影響較大[2-3]。因此，準確的區(qū)間設(shè)計洪水估算具有重要意義。

區(qū)間流域通常都缺少直接實測的洪水資料，因此只能采用間接資料或用間接方法計算。但在區(qū)間無暴雨資料或暴雨資料不足時，間接資料計算方法的應(yīng)用受到限制。間接方法根據(jù)資料條件的不同主要分為[4]：①當上下游斷面均有實測洪水資料時，可將歷年下游斷面洪水過程減去上游段面洪水過程(包括洪水演進)，得到歷年區(qū)間的洪水過程線。之后，按照傳統(tǒng)的規(guī)范方法進行區(qū)間洪水頻率計算。但這種方法在區(qū)間洪水所占比重較小時，二者相減的誤差將嚴重影響區(qū)間設(shè)計洪水成果的精度。②當區(qū)間某一部分流域面積(如支流)上有水文站時，可參照水文站控制流域面積與區(qū)間總面積以及暴雨、產(chǎn)匯流特性等，將水文站的歷年洪水過程轉(zhuǎn)換成整個區(qū)間的洪水過程，再進行區(qū)間洪水頻率計算。但當區(qū)間有水文站控制的流域面積較小時，使上述轉(zhuǎn)換產(chǎn)生較大誤差。方法①和②一般均只利用實測資料序列進行分析計算，難以考慮歷史洪水資料，應(yīng)用于資料系列較短的站點時，其估計結(jié)果難以滿足工程實際要求。③采用區(qū)間所在的暴雨等值線圖及暴雨徑流查算圖表進行計算，或者由相似流域經(jīng)地區(qū)綜合的經(jīng)驗公式來推求區(qū)間的設(shè)計洪水。這種方法一般精度較低，適用于無法采用上述方法①和②計算區(qū)間設(shè)計洪水的情形。

工程實際中以區(qū)間流域缺乏實測資料，而區(qū)間的控制邊界上下游斷面均有洪水資料及設(shè)計洪水成果的情況最為常見。由于區(qū)間洪量是下游斷面洪量與上游洪量之差，不僅受上、下游洪量的影響，而且上、下游洪量有一定關(guān)系。因此，區(qū)間洪量的頻率計算實質(zhì)上可以看作為上、下游洪量的頻率組合計算。本文從多變量函數(shù)的頻率分析與計算角度，提出一種計算區(qū)間流域設(shè)計洪水的新方法；并以清江流域為例，推求水布埡壩址與隔河巖壩址之間的區(qū)間流域的洪量頻率曲線，對該方法進行了驗證。

1 區(qū)間流域設(shè)計洪水計算的頻率組合模型

如圖1所示，A為下游斷面；B為上游斷面，C為斷面A與斷面B之間的區(qū)間流域。X、Y及Z分別表示斷面A、斷面B和區(qū)間C的天然來水量。

圖1 區(qū)間流域示意

由水量平衡原理得

Z=X-Y

(1)

現(xiàn)欲推求Z的概率分布FZ(z)。Z不僅受X、Y的影響，且X和Y一般不是獨立的，而是具有一定關(guān)系的。Z的頻率計算實質(zhì)上可以看成是X、Y的頻率組合計算。

根據(jù)定義有

FZ(z)=P(Z≤z)=P(X-Y≤z)

(2)

式中，f(x，y)為二元隨機變量(X，Y)的聯(lián)合概率密度函數(shù)。

對式(2)求導，可得

(3)

若已知聯(lián)合概率密度函數(shù)f(x，y)，在此基礎(chǔ)上再積分計算，就可得到Z的概率分布。我國一般假定X和Y均服從P-III型分布[1]，且它們之間還存在一定的相關(guān)性。若用傳統(tǒng)辦法來確定f(x，y)的解析表達式是非常困難的，因此直接應(yīng)用式(3)積分求解Z的概率分布尚存在一定的困難，實際應(yīng)用中只能采用一些近似的計算方法。

2 多變量函數(shù)頻率分析與計算方法

對于多變量函數(shù)的頻率分析與計算，常用的方法主要有統(tǒng)計參數(shù)法[6]和概率組合法[7]。隨著聯(lián)合分布模擬技術(shù)的不斷發(fā)展，尤其是近年來Copula函數(shù)在水文領(lǐng)域的成功應(yīng)用，學者們提出了基于Copula函數(shù)的蒙特卡羅法[8-10](以下簡稱“Copula-MC法”)和基于Copula函數(shù)的概率組合法[11](以下簡稱“Copula-PC法”)。本文將采用以上4種方法來計算Z的分布函數(shù)FZ(z)，并對各方法進行比較分析。

2.1 Copula函數(shù)、聯(lián)合分布和條件分布

(1)Copula函數(shù)及聯(lián)合分布。Copula函數(shù)是定義域為[0，1]均勻分布的多維聯(lián)合分布函數(shù)，根據(jù)Sklar定理，它可以將多個隨機變量的邊緣分布連接起來構(gòu)造聯(lián)合分布。選用在構(gòu)造上下游斷面洪水聯(lián)合分布時應(yīng)用最為廣泛的Gumbel-Hougaard Copula函數(shù)作為下游斷面A和上游斷面B洪量的聯(lián)合分布函數(shù)，數(shù)學表達式為[12]

F(x，y)=C(u，v)=exp{-[(-lnu)θ+
(-lnv)θ]1/θ}

(5)

式中，u=FX(x)和v=FY(y)分別代表邊緣分布函數(shù)；θ為Copula函數(shù)的參數(shù)，采用可以考慮歷史洪水信息的改進極大似然法估計。

(2)基于Copula函數(shù)的條件概率。當隨機變量X、Y的聯(lián)合分布F(x，y)確定后，給定X=x時，Y≤y的條件分布函數(shù)為

(6)

借助Copula函數(shù)，上式可表示為[13]

(7)

對于式(5)中Gumbel-Hougaard Copula函數(shù)，式(7)可具體表示為[14-15]

FY|X(y|x)=exp{-[(-lnu)θ+(-lnv)θ]1/θ}×

[(-lnu)θ+(-lnv)θ]1/θ-1(-lnu)θ-1/u

(8)

則依據(jù)概率分布的性質(zhì)可知，當給定X=x時，Y≥y的條件概率

Py|x=P(Y≥y|X=x)=1-FY|X(y|x)

(9)

2.2 統(tǒng)計參數(shù)法

統(tǒng)計參數(shù)法就是根據(jù)式(1)計算Z的統(tǒng)計參數(shù)(均值、變差系數(shù)和偏態(tài)系數(shù))；再假定Z的頻率曲線分布線型，即可得到Z的概率分布曲線。因此，該方法的關(guān)鍵是組合變量Z統(tǒng)計參數(shù)的計算。設(shè)隨機變量X、Y均有n年的資料序列，則對應(yīng)的Z也有n年的資料序列。

(1)均值

(10)

(2)變差系數(shù)。Z的方差為

(11)

(12)

(3)偏態(tài)系數(shù)

(13)

2.3 離散求和法

離散求和法是最常用的概率組合方法，將X與Y的頻率曲線離散化，概化成有限個“狀態(tài)”，即階梯狀，得到Z的所有組合狀態(tài)，再推求所有組合狀態(tài)的概率。該法一般假定X與Y相互獨立。

設(shè)X、Y的每一狀態(tài)對應(yīng)一概率區(qū)間，即X取狀態(tài)xi的概率區(qū)間為ΔPx，i(i=1，2，…，nx；nx為X的狀態(tài)數(shù))，Y取狀態(tài)yj的概率區(qū)間為ΔPy，j(j=1，2，…，ny；ny為Y的狀態(tài)數(shù))，則Z的狀態(tài)zij=xi-yj(zij的狀態(tài)數(shù)nz=nx×ny)，對應(yīng)的概率區(qū)間為ΔPz，ij，則

P(Z=zij)=ΔPz，ij=P(X=xi)P(Y=yj)=ΔPx.iΔPy，j

(14)

于是有

(15)

當X與Y不獨立時，一般應(yīng)進行獨立性處理，以新變量代換X或Y中的一個。在實際應(yīng)用中，一般對均值較小的組合變量作變量代換，效果較好。不失一般性，假定Y的均值較小，可作如下變量變換

E=Y-k·X

(16)

根據(jù)正交性原理，當k為最小二乘估計時，即

(17)

則經(jīng)獨立性處理后，E與X滿足正交條件，從而認為二者不相關(guān)；實際應(yīng)用中一般可視為是相互獨立的，可按獨立隨機變量一樣進行頻率組合計算。

2.4 Copula-PC法

當Y倚X的條件概率確定后，對于X的每一個取值狀態(tài)，都對應(yīng)一條Y倚X的條件概率曲線Py|x。將X的頻率曲線離散化，即概化成階梯狀，離散后X只能取有限個狀態(tài)值。對X的每一種取值狀態(tài)，將其對應(yīng)的Y倚X的條件概率曲線Py|x也離散化。設(shè)X取nx個狀態(tài)，取ny個狀態(tài)，則組合變量Z的狀態(tài)nz=nx·ny。

X、Py|x的每一狀態(tài)都對應(yīng)一概率區(qū)間，設(shè)X取狀態(tài)xi的概率區(qū)間為ΔPx，i，Y取狀態(tài)yj的條件概率區(qū)間為ΔPyj|xi(i=1， 2， …，nx;j=1， 2， …，ny)。Z相應(yīng)狀態(tài)對應(yīng)的概率區(qū)間為ΔPz，ij，則

P(Z=zij)=ΔPz，ij=P(X=xi)·P(Y=yj|X=xi)

=ΔPx，i·ΔPyj|xi

(18)

式中，zij=xi-yj；而

ΔPyj|xi=[1-FY|X(yj+1|xi)]-[1-FY|X(yj|xi)]

=FY|X(yj|xi)-FY|X(yj+1|xi)

(19)

于是式(19)可表示為

P(Z=zij)=ΔPz，ij=[FY|X(yj|xi)-

FY|X(yj+1|xi)]·ΔPx，i

(20)

則

(21)

可以看出，借助Copula函數(shù)可直接對Y倚X的條件概率曲線Py|x進行離散化，不需要對變量Y作獨立性處理，保持了變量間的相關(guān)性特征，避免了離散求和法在進行獨立性處理時可能出現(xiàn)的數(shù)據(jù)失真問題。

2.5 Copula-MC法

通過聯(lián)合分布函數(shù)可以對存在相關(guān)性的隨機變量X和Y成對地進行隨機抽樣，從而達到對二者同時隨機模擬的目的。

步驟1：根據(jù)建立的兩變量聯(lián)合分布，可以得到當X為指定值時Y的條件分布，數(shù)學表達式為

SX(y|X=x)=P(Y≤y|X=x)=?F(x，y)/?x

(22)

步驟2：產(chǎn)生服從[0，1]均勻分布的兩個獨立隨機數(shù)r1和r2。

步驟4：令r2等于X=x時Y的條件分布值，即r2=SX(y|X=x)，為求解方便，將式(7)表達為

SU(v|U=u)=P(V≤v|U=u)=?C(u，v)/?u

(23)

步驟5：通過z=x-y計算得到z。

步驟6：重復(fù)步驟2～5共n次，可以模擬出n個z值，再由期望公式計算經(jīng)驗頻率，得到Z的頻率曲線。P(z>zi)=mi/(n+1)，mi為滿足z

3 實例研究

3.1 研究區(qū)域概況

清江是湖北省境內(nèi)最大的一條長江支流，發(fā)源于湖北省利川市齊岳山，流域面積1.7萬km2，中下游建有水布埡-隔河巖梯級水庫。如圖1，A為隔河巖壩址斷面，B為水布埡壩址斷面，C為水布埡、隔河巖之間的區(qū)間流域(簡稱水-隔區(qū)間)。隔河巖壩址控制面積14 430 km2，其洪量由上游水布埡壩址(10 860 km2)和水-隔區(qū)間(3 570 km2)兩個分區(qū)組合形成。水-隔區(qū)間左岸有泗渡河，右岸有泗揚河、天池河等較大支流，盡管泗渡河上設(shè)有招徠河水文站(1960年～1983年)，集水面積792 km2，但只占水-隔區(qū)間面積的22%，資料系列也較短，因此間接方法②不宜采用?？紤]到水布埡水電站壩址和隔河巖水電站壩址均有1965年～2010年的實測流量資料，且水-隔區(qū)間面積占隔河巖壩址控制面積的比例較大(24.7%)，可采用間接方法①計算。本文將采用4種基于頻率組合的方法計算水-隔區(qū)間最大3 d洪量的頻率曲線，與采用傳統(tǒng)間接方法①結(jié)果進行比較。

3.2 邊緣分布、聯(lián)合分布和條件分布

本文基本資料為1965年～2010年(共46 a)3 h時段的水布埡水電站壩址和隔河巖水電站壩址流量。計算最大3 d洪量頻率曲線時隔河巖壩址1935作為特大洪水，水布埡壩址1883年和1997年作特大值處理。根據(jù)《水利水電工程設(shè)計洪水計算規(guī)范》，假設(shè)水布埡壩址和隔河巖壩址最大3 d洪量均服從P-Ⅲ分布，參數(shù)采用線性矩法初估，經(jīng)適線法微調(diào)得到，估計結(jié)果見表1。

表1 清江流域水布埡和隔河巖壩址最大3 d洪量統(tǒng)計參數(shù)

分別采用1965年～2010年水布埡壩址和隔河巖壩址最大3 d洪量資料，利用可考慮歷史洪水信息(實際上本實例無同步的歷史洪水，但如果有的話可以考慮)的改進極大似然法計算得到聯(lián)合分布的參數(shù)為6.05。將46 a的經(jīng)驗聯(lián)合分布值與理論聯(lián)合分布值點繪見圖2，其確定性系數(shù)為0.993 3，經(jīng)驗頻率與理論頻率值的擬合情況很好。圖3是另外一種形式對聯(lián)合觀測變量的經(jīng)驗聯(lián)合分布值與理論聯(lián)合分布值進行了對比。圖2和圖3表明所建立的聯(lián)合分布是合理可行的。有了聯(lián)合分布參數(shù)，分別代入式(5)和式(9)就可以得到聯(lián)合分布函數(shù)和條件分布函數(shù)。

圖2 聯(lián)合觀測值的經(jīng)驗分布和理論分布相關(guān)關(guān)系

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

3.3 水-隔區(qū)間最大3 d洪量頻率曲線

分別采用統(tǒng)計參數(shù)法、離散求和法、Copula-PC法和Copula-MC法，分別推求水-隔區(qū)間最大3 d洪量的頻率曲線及設(shè)計值。

圖4 5種方法計算的水-隔區(qū)間最大3 d洪量頻率分布曲線

zij=0.278xi-ej

(24)

該zij流量發(fā)生的概率為

圖5 水-隔區(qū)間E的頻率曲線

(25)

對X和E的全部取值狀態(tài)作組合計算，可得Z的頻率曲線(見圖4)。

(3)Copula-PC法。以X的設(shè)計頻率為0.1%和1%為例，即當X=x0.1%和X=x1%時，Y≥y的條件概率曲線分別見圖6。將X的頻率曲線和Y倚X的條件概率曲線都離散成40種狀態(tài)，設(shè)變量X第i種狀態(tài)的取值為xi，Y取狀態(tài)yj的條件概率區(qū)間為ΔPyj|xi，則

zij=xi-yj

(26)

圖6 給定X=x0.1%和X=x1%時Y的條件頻率曲線

該zij流量發(fā)生的概率為

ΔPz，ij=ΔPx，i·ΔPyj|xi=ΔPx，i·[FY|X(yj|xi)-

FY|X(yj+1|xi)]

(27)

對X和Y的全部取值狀態(tài)作組合計算，可得Z的頻率曲線(見圖7)。

圖7 水布埡與隔河巖最大3 d洪量隨機生成散點

(4)Copula-MC法。運用上述2.5節(jié)中的Copula-MC法計算步驟，模擬15 000組水布埡壩址和隔河巖壩址最大3 d洪量。將模擬結(jié)果和46 a實測資料結(jié)果點繪于圖7。由圖7可知，實測資料點據(jù)都位于隨機模擬散點圖中，對實測資料的吻合效果較好。Copula-MC法不但實現(xiàn)了水-隔區(qū)間最大3d洪量的隨機生成，同時增加了隔河巖壩址洪水地區(qū)組成隨機性和多樣性。利用上述隨機生成的大量水-隔區(qū)間最大3 d洪量，采用期望公式計算的經(jīng)驗頻率，得到Z的頻率曲線。

《隔河巖水庫汛限水位設(shè)計與運用技術(shù)報告》[16]中采用傳統(tǒng)的間接方法①計算水-隔區(qū)間最大3 d洪量的頻率曲線。將以上4種頻率組合方法和文獻[16]中的結(jié)果分別繪于圖4，同時將5種常用設(shè)計頻率的結(jié)果匯總于表2。

從圖4和表2的計算結(jié)果可以發(fā)現(xiàn)：與兩種基于Copula函數(shù)的Copula-PC法和Copula-MC法相比，傳統(tǒng)的統(tǒng)計參數(shù)法和離散求和法計算的結(jié)果明顯系統(tǒng)偏小，且偏小幅度隨著洪水重現(xiàn)期的增大而增大，會低估設(shè)計斷面的防洪風險。統(tǒng)計參數(shù)法只考慮上、下游斷面的線性相關(guān)關(guān)系，而未真實反映兩者復(fù)雜的非線性相關(guān)結(jié)構(gòu)；同時假定Z服從P-Ⅲ分布理論上也缺乏依據(jù)。離散求和法在獨立性轉(zhuǎn)換的過程中，區(qū)間E系列出現(xiàn)大量負值并人為假定E服從P-Ⅲ分布進行強制適線，以致出現(xiàn)數(shù)據(jù)失真，造成誤差較大。

表2 水-隔區(qū)間最大3 d洪量設(shè)計值對比

Copula-PC法采用Copula函數(shù)求解條件概率曲線，直接對其進行離散組合，不需要變量進行獨立性處理，保持了變量間的相關(guān)性特征。Copula-MC法則通過Copula函數(shù)構(gòu)建X和Y的聯(lián)合分布函數(shù)，采用蒙特卡羅方法對其進行聯(lián)合抽樣，隨機生成大量的水-隔區(qū)間最大3d洪量，從而得到頻率曲線，計算結(jié)果合理，可操作性強。

4 結(jié) 論

本文提出的方法用清江流域水布埡-隔河巖區(qū)間流域?qū)嵗M行驗證的結(jié)果表明：

(1)所提出的基于上下游斷面洪水頻率組合計算區(qū)間流域設(shè)計洪水的思路合理可行，推薦采用基于Copula函數(shù)的Copula-PC法和Copula-MC法。

(2)統(tǒng)計參數(shù)法和離散求和法理論上存在不足，計算的結(jié)果系統(tǒng)偏小，可能會低估下游斷面的防洪風險。

(3)本文方法對資料條件要求較低，一般容易滿足，且在頻率組合過程中可以考慮歷史洪水信息，更加符合工程實際要求。

參考文獻：

[1] 郭生練, 劉章君, 熊立華. 設(shè)計洪水計算方法研究進展與評價[J]. 水利學報, 2016, 47(3)： 302-314．

[2] 周祥林, 燕文明. 混聯(lián)型梯級水庫群的典型洪水地區(qū)組成[J]. 水力發(fā)電, 2016, 42(11)：91-94．

[3] 譚喬鳳, 雷曉輝, 王浩, 等. 考慮梯級水庫庫容補償和設(shè)計洪水不確定性的汛限水位動態(tài)控制域研究[J]. 四川大學學報： 工程科學版, 2017, 49(1)： 60- 69．

[4] 水利部長江水利委員會水文局, 水利部南京水文水資源研究所. 水利水電工程設(shè)計洪水計算手冊[M]. 北京： 中國水利水電出版社, 2001．

[5] SL 44—2006 水利水電工程設(shè)計洪水計算規(guī)范[S]．

[6] 金光炎. 一個計算組合頻率的新方法[J]. 中國水利, 1958(4)： 46- 52．

[7] 王銳琛, 陳源澤, 孫漢賢. 梯級水庫下游洪水概率分布的計算方法[J]. 水文, 1990(1)： 1- 8．

[8] 肖義, 郭生練, 熊立華, 等. 一種新的洪水過程隨機模擬方法研究[J]. 四川大學學報： 工程科學版, 2007, 39(2)： 55- 60．

[9] 張超, 胡志根, 劉全. 上游水電站控泄條件下施工導流風險分析[J]. 水利學報, 2012, 43(11)： 1328- 1333.

[10] 劉章君, 郭生練, 胡瑤, 等. 基于Copula-Monte Carlo法的水庫下游洪水概率分布研究[J]. 水力發(fā)電, 2015, 41(8)： 17- 22.

[11] 李天元, 郭生練, 劉章君, 等. 梯級水庫下游設(shè)計洪水計算方法研究[J]. 水利學報, 2014, 45(6)： 641- 648.

[12] 郭生練, 閆寶偉, 肖義, 等. Copula 函數(shù)在多變量水文分析計算中的應(yīng)用及研究進展[J]. 水文, 2008, 28(3)： 1- 7．

[13] 史黎翔, 宋松柏. 基于Copula函數(shù)的兩變量洪水重現(xiàn)期與設(shè)計值計算研究[J]. 水力發(fā)電學報, 2015, 34(10)： 27- 34．

[14] 馮平, 李新. 基于Copula函數(shù)的非一致性洪水峰量聯(lián)合分析[J]. 水利學報, 2013, 44(10)： 1137- 1147．

[15] 劉和昌, 梁忠民, 姚軼, 等. 基于Copula函數(shù)的水文變量條件組合分析[J]. 水力發(fā)電, 2014, 40(5)： 13- 16．

[16] 武漢大學水資源與水電工程科學國家重點實驗, 長江勘測規(guī)劃設(shè)計研究院, 長江水利委員會水文局. 隔河巖水庫汛限水位設(shè)計與運用[R]. 武漢： 武漢大學, 2010．