孫乃波, 劉克琳， 程亮， 杜慧華

(1.威海市水文中心，山東 威海 264209； 2.南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029)

隨著對(duì)江河的治理與開(kāi)發(fā)，水庫(kù)群的調(diào)節(jié)對(duì)下游防洪控制點(diǎn)洪水過(guò)程的影響愈來(lái)愈突出[1]。水庫(kù)群的削峰和滯洪作用，一般會(huì)使得控制點(diǎn)洪水洪峰減小，峰現(xiàn)時(shí)間后移，與天然條件下的洪水過(guò)程有較大差異。洪水過(guò)程變化導(dǎo)致控制點(diǎn)防洪風(fēng)險(xiǎn)率變化，科學(xué)評(píng)估水庫(kù)群調(diào)節(jié)作用下控制點(diǎn)的防洪風(fēng)險(xiǎn)率，對(duì)復(fù)核流域防洪標(biāo)準(zhǔn)、優(yōu)化水庫(kù)群防洪調(diào)度規(guī)則具有重要意義[2]。

當(dāng)只考慮洪水不確定性時(shí)，控制點(diǎn)防洪風(fēng)險(xiǎn)率可用一定時(shí)間內(nèi)控制點(diǎn)洪水流量超過(guò)河道安全泄量的概率來(lái)表示。因此，如何用概率描述水庫(kù)群調(diào)節(jié)作用下控制點(diǎn)的洪水過(guò)程，是評(píng)估防洪風(fēng)險(xiǎn)率的核心。設(shè)計(jì)洪水在概念上具有概率屬性，是洪水概率描述的一種常用方式。根據(jù)《水利水電工程設(shè)計(jì)洪水計(jì)算規(guī)范》(SL 44—2006)，探尋水庫(kù)群調(diào)節(jié)影響下防洪控制點(diǎn)設(shè)計(jì)洪水的方法，目前主要有地區(qū)組成法、頻率組合法[3]和隨機(jī)模擬法[4]。其中地區(qū)組成方法是規(guī)劃設(shè)計(jì)中最常采用的方法，包括典型洪水組成法和同頻率洪水組成法。近些年，有學(xué)者基于Copula函數(shù)提出條件期望組成和最可能組成[5]，或采用當(dāng)量正態(tài)化(JC)法擬定一種最不利的地區(qū)組成[6-7]，這些方法都是對(duì)地區(qū)組成法的改進(jìn)。

由于洪水在空間分布上存在不確定性，當(dāng)控制點(diǎn)發(fā)生一定頻率洪水時(shí)，其上游水庫(kù)和區(qū)間洪水可能存在多種組成方式，或者說(shuō)每一種地區(qū)組成都是一個(gè)概率事件[8]。地區(qū)組成法采用單一洪水地區(qū)組成，或者從多個(gè)地區(qū)組成中選擇一個(gè)最不利的地區(qū)組成，這對(duì)洪水地區(qū)組成的不確定性考慮不足。洪水地區(qū)組成法盡管在主觀上考慮對(duì)防洪不利的原則，但客觀上卻使防洪斷面不滿足預(yù)定的防洪標(biāo)準(zhǔn)[9]。頻率組合法和隨機(jī)模擬法雖然考慮了水庫(kù)與區(qū)間洪水的所有可能組成方式，但計(jì)算過(guò)程中沒(méi)有采用控制點(diǎn)天然洪水的概率分布，從而失去了與天然設(shè)計(jì)洪水以及地區(qū)組成法對(duì)比的基礎(chǔ)。

為此，本文結(jié)合地區(qū)組成法和頻率組合法的思路，以單庫(kù)防洪系統(tǒng)為對(duì)象，提出計(jì)算水庫(kù)下游控制點(diǎn)防洪風(fēng)險(xiǎn)率的離散求和-均值一次二階矩法(discretization summation-mean first order second moment method, DS-MFOSM)。該方法根據(jù)水庫(kù)和區(qū)間洪量的二維聯(lián)合分布函數(shù)推導(dǎo)某一設(shè)計(jì)頻率下洪水地區(qū)組成的概率密度函數(shù)，用以描述洪水地區(qū)組成的不確定性；通過(guò)引入水庫(kù)調(diào)洪函數(shù)和洪水峰量關(guān)系函數(shù)，采用均值一次二階矩法計(jì)算控制點(diǎn)組合洪峰流量在某一設(shè)計(jì)頻率下超過(guò)安全泄量的概率；最后，采用離散求和得到控制點(diǎn)的防洪風(fēng)險(xiǎn)率。該方法不僅考慮了洪水地區(qū)組成的不確定性，還能描述該不確定性在防洪系統(tǒng)中的傳遞特征。最后，以沭河流域的沙溝-青峰嶺防洪系統(tǒng)為實(shí)例，將該方法進(jìn)行了應(yīng)用。

1 基本原理與方法

單庫(kù)防洪系統(tǒng)如圖1所示，水庫(kù)A、區(qū)間和防洪控制點(diǎn)的天然洪量分別用隨機(jī)變量X、Y和Z表示,水庫(kù)下泄流量用qA表示。防洪控制點(diǎn)洪峰流量頻率曲線如圖2所示，圖2中線1為天然洪峰流量的頻率曲線，線2為在某一洪水地區(qū)組成下水庫(kù)調(diào)洪影響后洪峰流量的頻率曲線。傳統(tǒng)的地區(qū)組成法、頻率組合法和隨機(jī)模擬法都能得到一條確定的線2，進(jìn)而根據(jù)線2和安全泄量得到控制點(diǎn)的防洪風(fēng)險(xiǎn)率。由于洪水地區(qū)組成具有不確定性，對(duì)于某一設(shè)計(jì)頻率pi，水庫(kù)調(diào)洪影響后防洪控制點(diǎn)的組合洪峰流量不是一個(gè)確定值，而是一個(gè)隨機(jī)變量，這是本文與傳統(tǒng)方法的主要區(qū)別所在。

圖1 單庫(kù)防洪系統(tǒng)示意圖

圖2 防洪控制點(diǎn)洪峰流量-頻率曲線圖

假定隨機(jī)變量X、Y和Z的取值分別為x、y和z。防洪控制點(diǎn)設(shè)計(jì)頻率pi的天然洪量設(shè)計(jì)值為zpi，根據(jù)洪水地區(qū)組成的原理，存在以下關(guān)系式：

x+y=zpi，x≥0，y≥0。

(1)

隨機(jī)變量X、Y一般服從P-Ⅲ型分布，其概率密度函數(shù)分別為：

(2)

(3)

式中αx、βx、γx和αy、βy、γy分別為隨機(jī)變量X和Y的P-Ⅲ型分布的參數(shù)。

X、Y的概率密度分布函數(shù)分別表示為FX(x)、FY(y)，其聯(lián)合分布函數(shù)可用二維Copula函數(shù)表示：

FX,Y(x,y)=Cθ(u,v)=Cθ(FX(x),FY(y))。

(4)

對(duì)應(yīng)的二維概率密度函數(shù)表示為：

fX,Y(x,y)=cθ(FX(x),FY(y))fX(x)fY(y)。

(5)

根據(jù)式(1)和式(5)，洪水地區(qū)組成(x，zpi-x)的概率分布可用歸一化的條件概率密度函數(shù)g(x)進(jìn)行描述，其表達(dá)式為[8]：

(6)

引入水庫(kù)調(diào)洪函數(shù)h(x)，依據(jù)水庫(kù)調(diào)洪規(guī)則，通過(guò)擬合水庫(kù)最大下泄流量h(x)與入庫(kù)洪量x的關(guān)系獲得調(diào)洪函數(shù)。對(duì)于下游有防洪任務(wù)的水庫(kù)而言，調(diào)洪函數(shù)一般為分段函數(shù)。此外，建立區(qū)間洪峰流量與洪量關(guān)系函數(shù)u(y)，可通過(guò)擬合區(qū)間洪峰流量u(y)與區(qū)間洪量y的散點(diǎn)系列得到峰量關(guān)系函數(shù)。

設(shè)防洪控制點(diǎn)的安全流量為qan，假定水庫(kù)下泄洪峰與區(qū)間洪峰在控制點(diǎn)相遇，則設(shè)計(jì)頻率為pi時(shí)的防洪控制點(diǎn)洪峰流量超過(guò)安全泄量的概率Pfi為：

Pfi=P(h(x)+u(zpi-x)>qan)。

(7)

采用均值一次二階矩法求解式(7)，并建立功能函數(shù)w(x)：

w(x)=qan-h(x)-u(zpi-x)。

(8)

將功能函數(shù)w(x)在均值點(diǎn)μx處展開(kāi)為泰勒級(jí)數(shù)并保留至一次項(xiàng)，得到：

w(x)=w(μx)+w′(μx)(x-μx)。

(9)

式中w′(x)為w(x)對(duì)x的導(dǎo)數(shù)。

根據(jù)式(8)和式(9)，統(tǒng)計(jì)w(x)的均值μw和均方差σw，計(jì)算式分別如下：

μw=w(μx)=qan-h(μx)-u(zpi-μx)，

(10)

σw=w′(μx)σx=(1-h′(μx)+u′(zpi-μx))σx。

(11)

式中：h′(x)為h(x)對(duì)x的導(dǎo)數(shù)；u′(x)為u(x)對(duì)x的導(dǎo)數(shù)。

式(11)描述了洪水地區(qū)組成不確定性通過(guò)方差在防洪系統(tǒng)中的傳遞過(guò)程。Pfi可由下式估算：

(12)

式中Φ(x)為標(biāo)準(zhǔn)正態(tài)分布函數(shù)；μx、σx分別根據(jù)式(6)提供的條件概率密度函數(shù)g(x)得到：

(13)

(14)

由式(13)和(14)一般無(wú)法獲取解析解，可采用數(shù)值積分法得到數(shù)值解。

根據(jù)離散求和法的思想，在頻率曲線的橫坐標(biāo)上選擇n個(gè)典型頻率值，從小到大排序分別為pi(i=1、2、…、n)，以此為中心把整個(gè)頻率的0～1區(qū)間劃分成n個(gè)小區(qū)間，每個(gè)小區(qū)間的概率為：

(15)

水庫(kù)下游控制點(diǎn)的防洪風(fēng)險(xiǎn)率Pf估算式為：

(16)

2 計(jì)算流程

離散求和-均值一次二階矩法(DS-MFOSM)計(jì)算水庫(kù)下游控制點(diǎn)防洪風(fēng)險(xiǎn)率的流程如圖3所示，具體如下：

圖3 離散求和-均值一次二階矩法計(jì)算流程圖

1)根據(jù)水庫(kù)和區(qū)間洪量系列計(jì)算邊緣分布函數(shù)FX(x)、FY(y)和概率密度函數(shù)fX(x)、fY(y)；

2)利用二維Copula函數(shù)構(gòu)造水庫(kù)和區(qū)間洪量的聯(lián)合分布函數(shù)FX,Y(x,y)和概率密度函數(shù)fX,Y(x,y)；

3)將控制點(diǎn)天然洪量在頻率曲線圖上進(jìn)行離散，得到pi，利用式(15)計(jì)算Δpi，并在頻率曲線上查找zpi；

4)對(duì)于任一個(gè)zpi，采用數(shù)值計(jì)算的方法，根據(jù)式(6)計(jì)算地區(qū)組成的概率密度函數(shù)g(x)，再根據(jù)式(13)、(14)計(jì)算其均值μx和均方差σx；

5)構(gòu)建水庫(kù)調(diào)洪函數(shù)h(x)和區(qū)間洪水峰量關(guān)系函數(shù)u(y)，分別推導(dǎo)其導(dǎo)數(shù)函數(shù)h′(x)、u′(y)；

6)確定控制點(diǎn)安全流量qan，根據(jù)式(12)計(jì)算設(shè)計(jì)頻率為pi的風(fēng)險(xiǎn)率Pfi；

7)根據(jù)式(16)計(jì)算水庫(kù)下游控制點(diǎn)的防洪風(fēng)險(xiǎn)率Pf。

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

以沭河流域的沙溝-青峰嶺防洪系統(tǒng)作為應(yīng)用實(shí)例。沙溝水庫(kù)位于沂沭河水系沭河干流的上游，控制流域面積163 km2，總庫(kù)容1.02億m2，汛限水位234 m，百年一遇設(shè)計(jì)洪水位237.32 m，千年一遇校核洪水位239.30 m。沙溝水庫(kù)防洪調(diào)度規(guī)則是20年一遇洪水控泄500 m3/s，超過(guò)20年一遇洪水水庫(kù)敞泄，水庫(kù)下游38 km處的青峰嶺水庫(kù)為沙溝水庫(kù)的下游防洪控制點(diǎn)。

防洪系統(tǒng)各分區(qū)設(shè)計(jì)洪水歷時(shí)均為3 d，因此選擇最大3 d洪量作為洪水地區(qū)組成的特征量，其天然設(shè)計(jì)成果見(jiàn)表1,表中Cv、Cs分別為變差系數(shù)和偏態(tài)系數(shù)。

表1 各分區(qū)最大3 d洪量天然設(shè)計(jì)成果

利用1964—2011年沙溝水庫(kù)和沙溝-青峰嶺區(qū)間(沙青區(qū)間)最大3 d洪量統(tǒng)計(jì)系列，采用Gumbel-Hougaard Copula函數(shù)構(gòu)造二維聯(lián)合分布函數(shù)，參數(shù)θ估計(jì)值為3.125，分布函數(shù)如圖4所示。

圖4 水庫(kù)與區(qū)間3 d洪量聯(lián)合分布函數(shù)圖

根據(jù)沙溝水庫(kù)調(diào)洪規(guī)則和設(shè)計(jì)洪水成果，構(gòu)造沙溝水庫(kù)調(diào)洪函數(shù)：

(17)

根據(jù)1964—2011年沙青區(qū)間的洪峰流量與年最大3 d洪量系列，構(gòu)造洪水峰量與洪量關(guān)系函數(shù)：

u(y)=22.779y-61.757。

(18)

選擇沙溝水庫(kù)防洪標(biāo)準(zhǔn)對(duì)應(yīng)的設(shè)計(jì)頻率5%，根據(jù)表1，該頻率下青峰嶺天然設(shè)計(jì)3 d洪量z5%=115.20×106m3，青峰嶺安全泄量為2 500 m3/s。根據(jù)式(6)計(jì)算出設(shè)計(jì)頻率為5%的地區(qū)組成概率密度函數(shù)g(x)，如圖5所示，根據(jù)數(shù)值積分法計(jì)算出該概率密度函數(shù)的均值為34.31×106m3，均方差為8.68×106m3。

圖5 設(shè)計(jì)頻率為5%的地區(qū)組成密度函數(shù)圖

將0～1區(qū)間頻率進(jìn)行離散，在防洪標(biāo)準(zhǔn)頻率附近適當(dāng)加密，計(jì)算青峰嶺斷面各設(shè)計(jì)頻率下的防洪風(fēng)險(xiǎn)率及總的防洪風(fēng)險(xiǎn)率，結(jié)果見(jiàn)表2。

表2 青峰嶺斷面防洪風(fēng)險(xiǎn)率計(jì)算結(jié)果

將表2中最后一行數(shù)據(jù)相加可得，青峰嶺斷面的防洪風(fēng)險(xiǎn)率為4.02%，小于沙溝水庫(kù)防洪標(biāo)準(zhǔn)頻率5%，表明在現(xiàn)行的沙溝水庫(kù)防洪調(diào)度規(guī)則下，沙溝-青峰嶺防洪系統(tǒng)滿足防洪標(biāo)準(zhǔn)要求。

圖6為青峰嶺斷面設(shè)計(jì)頻率pi與風(fēng)險(xiǎn)率Pfi相關(guān)圖。從圖6中可知，隨著設(shè)計(jì)頻率pi的增大，設(shè)計(jì)洪水變小，設(shè)計(jì)頻率下pi的風(fēng)險(xiǎn)率Pfi隨之減小。

圖6 青峰嶺斷面設(shè)計(jì)頻率pi與風(fēng)險(xiǎn)率Pfi相關(guān)圖

4 結(jié)論

以單庫(kù)防洪系統(tǒng)為研究對(duì)象，提出了計(jì)算水庫(kù)下游控制點(diǎn)防洪風(fēng)險(xiǎn)率的離散求和-均值一次二階矩法，該方法以某一設(shè)計(jì)頻率下的地區(qū)組成概率密度函數(shù)描述洪水地區(qū)組成的不確定性，結(jié)合水庫(kù)調(diào)洪函數(shù)和洪水峰量與3 d洪量關(guān)系函數(shù)，采用均值一次二階矩法計(jì)算某一設(shè)計(jì)頻率下的控制點(diǎn)洪峰流量超過(guò)安全流量的概率，進(jìn)而采用離散求和法計(jì)算控制點(diǎn)總的防洪風(fēng)險(xiǎn)率。最后將該方法應(yīng)用于沭河的沙溝水庫(kù)-青峰嶺斷面防洪系統(tǒng)。主要結(jié)論有：

1)離散求和-均值一次二階矩法既采用了地區(qū)組成法中的先分配、再組合的思路，又考慮了設(shè)計(jì)洪水地區(qū)組成的不確定性，并且能夠描述該不確定性在水庫(kù)防洪系統(tǒng)中的傳遞特征。

2)實(shí)例應(yīng)用表明，在沙溝水庫(kù)防洪調(diào)度規(guī)則下，沙溝-青峰嶺防洪系統(tǒng)滿足防洪標(biāo)準(zhǔn)要求。

3)離散求和-均值一次二階矩法給出了計(jì)算防洪風(fēng)險(xiǎn)率的解析式，但受限于聯(lián)合分布函數(shù)的維數(shù)，該方法僅適用于包含區(qū)間洪水的單庫(kù)防洪系統(tǒng)。