黃鋒華, 黃本勝，郭 磊，譚 超，胡 培(1.廣東省水利水電科學(xué)研究院，廣州 510610；2.廣東省水動力學(xué)應(yīng)用研究重點實驗室，廣州 510610；3.河口水利技術(shù)國家地方聯(lián)合工程實驗室，廣州 510610)

對于某一特定流域來說，流域防洪排澇設(shè)計不僅與流域集雨范圍內(nèi)降雨形成洪水有關(guān)，還與外江洪水相遭遇的外江洪水位息息相關(guān)；也就是說，某一流域設(shè)計洪水和與之相組合的外江設(shè)計洪水(流量或水位)是分析該流域排澇方案和排澇設(shè)施規(guī)模的依據(jù)。因此，為合理確定流域防洪排澇設(shè)計規(guī)模，必須知道流域暴雨洪水與外江洪水遭遇的概率風(fēng)險規(guī)律，從而選定的流域設(shè)計洪水與外江洪水組合對防洪排澇計算才真正的具有重現(xiàn)期意義。

對于洪水的遭遇分析，常用的方法是根據(jù)研究區(qū)域歷年同步洪水資料進行水文學(xué)分析法或簡單統(tǒng)計[1-3]，這些方法只能對已發(fā)生的洪水進行簡單的統(tǒng)計，對于洪水發(fā)生概率分布、洪水遭遇重現(xiàn)期的風(fēng)險內(nèi)容難以做到定量分析計算。近年，多變量聯(lián)合分布在多變量極端洪水要素(洪峰、洪量、峰現(xiàn)時間)統(tǒng)計分析中成為研究熱點，郭生練[4,5]等總結(jié)了Copula函數(shù)在水文分析計算的應(yīng)用和研究進展，提出了Archimedean Copula 函數(shù)中幾種常用于水文分析計算的函數(shù)及其適用性，其中Frank Copula函數(shù)能夠描述正/負(fù)相關(guān)的隨機變量，且對相關(guān)性的程度沒有限制，適用性廣。Copula函數(shù)因其結(jié)構(gòu)簡單、擬合效果優(yōu)良等特點在多變量水文分析計算中得到了廣泛的應(yīng)用[6-9]，劉曾美等[10-12]將Copula函數(shù)成功應(yīng)用于洪水遭遇分析中?？紤]到望江瀝作為無資料小流域地區(qū)，洪水歷時一般在1天內(nèi)且降雨與洪水是同頻的，本文采用望江瀝(惠陽站為代表站)年最大24 h降雨作為年最大洪水作為望江瀝暴雨洪水計算依據(jù)，采用Von Mises分布擬合博羅站年最大洪水和惠陽站年最大24 h降雨出現(xiàn)時間分布規(guī)律，通過Copula函數(shù)構(gòu)建東江干流與支流望江瀝年最大洪水發(fā)生時間及其量級聯(lián)合分布，探討東江干支流洪水遭遇風(fēng)險特征，為望江瀝流域防洪排澇設(shè)計和風(fēng)險管理提供的理論支撐。

1 研究區(qū)域概況

東江是珠江流域在廣東境內(nèi)的三大水系之一，干流全長562 km，流域總面積35 340 km2，其中廣東省境內(nèi)31 840 km2；博羅水文站是國家基本水文站、東江流域下游干流關(guān)鍵控制站，控制流域面積25 325 km2，占全流域面積的72%，是東江干流匯入珠江三角洲河網(wǎng)的控制性站點，是東江流域的代表性水文站?；蓐栒臼菛|江流域重要水位雨量站，其位于望江瀝河涌匯入東江干流交匯口的下游3 km(見圖1)，可作為望江瀝流域(集雨面積不足20 km2的小流域)代表性雨量站。望江瀝位于東江下游廣東惠州市惠城區(qū)惠州大堤(北堤)的保護范圍內(nèi)，為東江河谷平原區(qū)，集雨面積為17.12 km2，河道平面干流河長8.3 km。望江瀝流域所在區(qū)域位于惠城區(qū)，是惠州市的政府所在地和中心區(qū)，是惠州乃至整個東江流域政治、經(jīng)濟、文化和交通的中心，素有“粵東重鎮(zhèn)”、“嶺南名郡”的美譽。近年來，隨著望江瀝其周邊城市迅速發(fā)展，流域內(nèi)人口密集、經(jīng)濟總量大，對流域內(nèi)防洪排澇設(shè)計提出了更高的要求，特別是惠州市成為第二批全國水生態(tài)文明城市建設(shè)試點后，提出了以“河清、湖美、惠民城”為創(chuàng)建理念和目標(biāo)，將重點完成惠州城區(qū)內(nèi)望江瀝等河涌綜合整治工程，對河涌防洪排澇設(shè)計工作提出更為嚴(yán)格的要求?；葜菔兄行某菂^(qū)有大小河涌10余條，均匯入東江干流，根據(jù)惠州市水生態(tài)文明實施方案，惠州市每年啟動1～2條河涌綜合整治項目，河涌與東江干流洪水遭遇風(fēng)險分析，可為河涌綜合整治防洪排澇方案設(shè)計提供依據(jù)。

圖1 東江干流、望江瀝及相關(guān)站點位置示意圖

2 研究方法

2.1 Copula函數(shù)

Sklar定理：設(shè)X、Y為連續(xù)的隨機變量，其邊緣分布函數(shù)分別為F(x)和F(y)，F(xiàn)(x，y)為變量X和Y的聯(lián)合分布函數(shù)，那么存在唯一的Copula函數(shù)C使得：

F(x,y)=P(X≤x,Y≤y)=Cθ[FX(x),FY(y)],?x,y

(1)

式中：Cθ[FX(x),FY(y)]為Copula函數(shù)；θ為Copula 函數(shù)的參數(shù)。

在水文分析計算中最常用的幾種函數(shù)有Gumbel Copula、Clayton Copula、AMH Copula和Frank Copula。本文選用適用性較強的Frank Copula進行量變量聯(lián)合分布構(gòu)建，其表達(dá)式如下：

(2)

式中：u、v均為單變量邊緣分布；θ為Copula函數(shù)參數(shù)。

采用相關(guān)指標(biāo)法[4]進行參數(shù)估計，即根據(jù)Kendall秩相關(guān)系數(shù)τ與參數(shù)θ之間的關(guān)系間接求出，表達(dá)式如下：

(3)

2.2 年最大洪水發(fā)生時間分布

方彬等[13]首次采用Von Mises分布擬合年最大洪水發(fā)生時間，研究表明Von Mises分布能較好地擬合年最大洪水發(fā)生時間。設(shè)洪水發(fā)生時間隨機變量X服從Von Mises分布，則其密度函為：

(4)

0≤x≤2π, 0≤μ≤2π,k>0

式中：μ為位置參數(shù)(平均角)；k為尺度參數(shù)；I0(k)為第一類0階修正的Bessel函數(shù)。

對于雙峰混合Von Mises分布則是由兩個單一的Von Mises分布按一定比例p相加而得，其密度函數(shù)為：

f(x,μ1,μ2,k1,k2,p)=

p·f1(x,u1,k1)+(1-p)·f2(x,u2,k2)

(5)

式中：f1(x,u1,k1)、f2(x,u2,k2)表示兩個單一的Von Mises分布；p表示兩個單一Von Mises分布在混合分布中的比例。關(guān)于Von Mises分布函數(shù)參數(shù)可采用最小二乘法獲取。

2.3 年最大洪水量級分布

根據(jù)《水利水電工程設(shè)計洪水計算規(guī)范(SL44-2006)》，采用皮爾遜Ⅲ型分布計算設(shè)計洪水，即年最大洪峰流量X和年最大24 h降雨Y服從皮爾遜Ⅲ型分布，其概率密度函數(shù)為：

(6)

α>0,β>0,δ≤y<+∞

式中：α、β和δ分別為皮爾遜Ⅲ型分布的形狀、尺度和位置參數(shù)；Γ(*)為伽馬函數(shù)。皮爾遜Ⅲ型分布參數(shù)估算可采用極大似然估計方法(MLE)[14]。

2.4 聯(lián)合概率分布與重現(xiàn)期

根據(jù)Copula函數(shù)，兩變量聯(lián)合分布為：

F(x,y)=P(X≤x,Y≤y)=C[FX(x),FY(y)]

(7)

設(shè)FX(x)為博羅站年最大洪峰流量X邊緣分布函數(shù)，F(xiàn)Y(y)為惠陽站年最大24 h降雨量Y的分布函數(shù)。

則X和Y的聯(lián)合重現(xiàn)期T0：

(8)

X和Y的同現(xiàn)重現(xiàn)期Tα：

(9)

邊緣分布的重現(xiàn)期：

(10)

則有：

T0(x,y)≤min[T(x),T(y)]≤max[T(x),T(y)]≤Tα(x,y)

(11)

2.5 條件概率分布

F(x,y)=P(X≤x,Y≤y)的物理意義為博羅站年最大洪峰流量小于x, 同時惠陽站年最大24 h降雨量小于y的概率，其重現(xiàn)期為年。

給定X>x，Y的條件分布函數(shù)為：

(12)

給定Y>y，X的條件分布函數(shù)為：

(13)

3 東江干支流洪水遭遇風(fēng)險分析

本文研究以1953-2013年(系列長度N1=61)博羅水文站年最大洪峰流量及其出現(xiàn)時間和惠陽站1974-2013年(系列長度N2=40)年最大24 h降雨及其出現(xiàn)時間序列為基礎(chǔ)，假定年最大洪水均發(fā)生在汛期，并按廣東汛期選取4月15日至10月15日(總天數(shù)L=184 d)的流量資料分別選取兩站年最大洪水和年最大24 h降雨出現(xiàn)日期作為樣本。本文所采用的水文氣象數(shù)據(jù)來源于廣東省水文局整編資料，數(shù)據(jù)可靠。

3.1 東江干支流洪水發(fā)生時間的遭遇分析

經(jīng)樣本分析，博羅站年最大洪水出現(xiàn)時間系列屬于單峰過程，惠陽站年最大24 h降雨出現(xiàn)時間系列屬于雙峰過程。采用最小二乘法進行參數(shù)尋優(yōu)，得到博羅站年最大洪水出現(xiàn)時間和惠陽站年最大24 h降雨出現(xiàn)時間Von Mises分布參數(shù)估算結(jié)果，兩站擬合的均方根誤差RMSE均較小，如表1所示。

表1 Von Mises分布參數(shù)估算結(jié)果

以時間弧度(即4月15日對應(yīng)弧度為0，10月15日對應(yīng)弧度為2π)為橫坐標(biāo)，圖2給出了博羅站年最大洪水出現(xiàn)時間概率密度分布fb和惠陽站年最大24 h降雨出現(xiàn)時間概率密度分布fh擬合情況，由圖2可知，Von Mises分布能較好擬合博羅站和惠陽站極端水文事件出現(xiàn)時間分布。博羅站年最大洪水出現(xiàn)時間概率密度曲線為單峰過程，在7月4日出現(xiàn)概率最大，約為34.1%；惠陽站年最大24h降雨出現(xiàn)時間概率密度曲線為雙峰過程，分別在6月14日和8月16日出現(xiàn)概率較大，分別為30.7%和25.2%。

圖2 博羅站年最大洪水和惠陽站年最大24 h降雨出現(xiàn)時間概率密度擬合圖

根據(jù)博羅站和惠陽站同步觀測數(shù)據(jù)分析，fb和fh兩同步序列Kendall系數(shù) 為0.112，呈現(xiàn)弱的正相關(guān)。根據(jù)Copula函數(shù)選取原則，選取Frank Copula函數(shù)構(gòu)建兩站的極端水文事件出現(xiàn)時間聯(lián)合分布，參數(shù)估算采取相關(guān)性指標(biāo)法，Copula函數(shù)系數(shù) 為1.02，均方誤差RMSE為0.189，采用Gringorten公式[15]計算聯(lián)合觀測值的經(jīng)驗頻率，并將博羅站年最大洪峰流量出現(xiàn)時間X和惠陽站年最大24 h降雨Y的經(jīng)驗頻率和理論頻率P-P繪于圖3，從圖中可知散點均落在45°直線對角線附近，表明所選Frank Copula函數(shù)擬合效果較好。在5%的顯著水平下，F(xiàn)rank Copula聯(lián)合分布χ2=0.29<χ20.995(31)=15.46，通過了χ2檢驗。

圖3 年最大洪水和年最大24 h降雨出現(xiàn)時間經(jīng)驗值與理論值P-P圖

根據(jù)Copula函數(shù)定義，望江瀝洪峰與東江干流洪峰遭遇是指在不考慮洪峰量級情況下，兩者洪峰在同一天出現(xiàn)，由此可以定義東江干流(博羅站)與望江瀝年最大洪水發(fā)生時間在第i天遭遇的風(fēng)險為：

pi=P(ti≤X

F(ti,ti)+F(ti+1,ti+1)-F(ti,ti+1)-F(ti+1,ti)

(14)

式中：pi為東江干流(博羅站)與其支流望江瀝年最大洪水發(fā)生時間在第i天遭遇的概率；F(ti,tj)為東江干流(博羅站)和其支流望江瀝洪水分布發(fā)生在第i天和第j天前的聯(lián)合概率。

圖4給出了東江干流(博羅站)與其支流望江瀝洪峰日遭遇的概率分布曲線。從圖中可以看出，東江干流與其支流望江瀝遭遇最大可能發(fā)生在6月17日，遭遇風(fēng)險約為0.011%，6月2日至7月9日是東江干支流出現(xiàn)洪峰遭遇的集中時期，7月10日至8月10有一個相對平穩(wěn)的遭遇風(fēng)險，約為0.006%，東江干支流在5月9日之前和10月6日之后發(fā)生遭遇的概率很小，幾乎為0。

圖4 東江干流與其支流望江瀝洪峰日遭遇的風(fēng)險分布

3.2 東江干支流洪水量級的遭遇分析

博羅站年最大洪水和惠陽站年最大24 h降雨邊緣分布及聯(lián)合分布參數(shù)估算成果見表2，由表2可知，兩站皮爾遜Ⅲ型曲線對經(jīng)驗樣本擬合的均方根誤差RMSE較小。統(tǒng)計表明，東江干流博羅站年最大洪峰流量和惠陽站年最大24 h降雨量成弱的正相關(guān)，F(xiàn)rank Copula函數(shù)構(gòu)建了兩者的聯(lián)合概率分布，其參數(shù)θ依據(jù)Kendall系數(shù)求得為2.69，均方誤差RMSE為0.296，在5%的顯著水平下，F(xiàn)rank Copula聯(lián)合分布χ2=0.28<χ20.995(32)=15.13，通過了χ2檢驗，其經(jīng)驗值和理論值的P-P圖集中分布在45°直線對角線附近，為省去篇幅不再列出圖。

表2 博羅站年最大洪水和惠陽站年最大24 h降雨邊緣分布及聯(lián)合分布參數(shù)估算表

根據(jù)Copula函數(shù)定義，東江博羅站出現(xiàn)年最大洪峰流量大于某一量級洪水，遭遇東江支流望江瀝(惠陽站)發(fā)生年最大24 h降雨量大于某一量級暴雨的風(fēng)險或重現(xiàn)期，可由式(9)表示。圖5給出了博羅站年最大洪峰流量X和惠陽站年最大24 h降雨量Y聯(lián)合分布三維圖、聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期的等值線圖，從圖5中可以更加清楚地了解聯(lián)合分布函數(shù)。從表3可知，東江干流(博羅站)發(fā)生年最大洪峰流量超過6 579 m3/s的重現(xiàn)期為5年，東江支流望江瀝(惠陽站)發(fā)生年最大24 h降雨量超過210 mm的重現(xiàn)期為20年，則其聯(lián)合重現(xiàn)期T0為4年，即東江干流博羅站發(fā)生超5年一遇洪水或者東江支流望江瀝(惠陽站)發(fā)生超20年一遇暴雨洪水的重現(xiàn)期為4年；而東江支流望江瀝(惠陽站)發(fā)生超20年一遇暴雨同時遭遇東江干流博羅站發(fā)生超5年一遇洪水的重現(xiàn)期為47年，即同現(xiàn)重現(xiàn)期Ta為47年。

圖5 年最大洪水和年最大24 h降雨聯(lián)合分布三維圖、聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期等值線圖

表3 博羅站年最大洪水和惠陽站年最大24 h降雨聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期

3.3 東江干支流洪水條件分布

當(dāng)東江支流望江瀝(惠陽站)發(fā)生某一量級以上的暴雨洪水的條件下，東江干流(博羅站)的概率分布，對研究望江瀝設(shè)計洪水與東江干流洪水遭遇具有重要意義。圖6給出了東江干流(博羅站)發(fā)生特定量級以上洪水條件下，望江瀝(惠陽站)發(fā)生某一量級以上暴雨洪水的概率分布曲線；圖7給出了望江瀝(惠陽站)發(fā)生特定量級以上暴雨洪水條件下，東江干流(博羅站)發(fā)生某一量級以上洪水的概率分布曲線。表4給出了望江瀝(惠陽站)發(fā)生P=5%暴雨洪水條件下，東江干流(博羅站)年最大洪水洪峰量級的概率分布。例如，望江瀝(惠陽站)發(fā)生20年一遇暴雨洪水的條件下，東江干流(博羅站)發(fā)生洪峰量級超過6 579 m3/s(5年一遇，見表3)洪水重現(xiàn)期為4年，發(fā)生洪峰量級超過10 996 m3/s(100年一遇，見表3)洪水重現(xiàn)期僅為38年；表明望江瀝發(fā)生暴雨洪水的情況下，東江干流發(fā)生某一量級洪水的可能性提高。

圖6 Y的條件頻率P(X>y|X

圖7 X的條件頻率P(X>x|Y>y)曲線

表4 惠陽站發(fā)生大于P=5%暴雨條件下博羅站洪水量級概率分布

4 結(jié) 論

以東江下游博羅站年最大洪水和望江瀝臨近的惠陽站年最大24 h降雨資料為研究基礎(chǔ)，采用Von Mises分布擬合暴雨洪水發(fā)生時間的概率分布，采用皮爾遜Ⅲ型分布擬合年最大洪水量級的概率分布，采用Frank Copula 函數(shù)構(gòu)建了東江干流與支流望江瀝洪水遭遇的聯(lián)合分布，各分布函數(shù)的擬合效果較好，在0.5%顯著水平下均通過了χ2檢驗。本文從暴雨洪水出現(xiàn)時間和量級分析了東江干流與支流望江瀝聯(lián)合概率分布特征以及遭遇分析分析，研究表明：東江干流(博羅站)和望江瀝暴雨洪水出現(xiàn)時間和量級之間存在弱的正相關(guān)；東江干流(博羅站)與望江瀝暴雨洪水在同一天遭遇最大可能發(fā)生在6月17日，遭遇風(fēng)險約為0.011%，在5月9日之前和10月6日之后發(fā)生遭遇的可能性極小，幾乎不會發(fā)生；望江瀝發(fā)生超20年一遇暴雨同時遭遇東江干流超5年一遇洪水的重現(xiàn)期為47年；條件概率分布上看，在望江瀝發(fā)生暴雨洪水條件下，東江干流(博羅站)某一量級的洪水發(fā)生的可能性相對提高。本研究較好地揭示了東江支流望江瀝河涌與外江洪水遭遇的概率風(fēng)險規(guī)律，可為望江瀝河涌排澇布置方案和排澇設(shè)施規(guī)模設(shè)計提供參考。

□

[1] 肖天國. 金沙江、岷江洪水遭遇分析[J]. 人民長江,2001,32(1):30-32.

[2] 李孝兵. 連江中上游干支流洪水遭遇問題分析[J]. 廣東水利水電,2008,(5):42-48.

[3] 戴明龍, 沈燕舟.長江上游與洞庭湖洪水遭遇規(guī)律研究[J]. 水資源與水工程學(xué)報,2010,21(6):117-120.

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

[5] 閆寶偉，郭生練，郭 靖，等. 多變量水文分析計算方法的比較[J]. 武漢大學(xué)學(xué)報(工學(xué)版),2009,42(1):10-15.

[6] 閆寶偉，郭生練，陳 璐，等.Copula函數(shù)在水文計算中的適用性分析[J]. 數(shù)學(xué)的實踐與認(rèn)識,2012,42(3):85-93.

[7] 陸桂華,閆桂霞,吳志勇,等．基于Copula函數(shù)的區(qū)域干旱分析方法[J]..水科學(xué)進展，2010，21(2):188-93.

[8] 陳子燊, 劉曾美,劍飛洪. 廣東西江北江洪水聯(lián)合概率分布研究[J].中山大學(xué)學(xué)報(自然科學(xué)版)，2011,50(2):110-115.

[9] 肖 義, 郭生練, 劉 攀，等. 基于兩變量分布的峰量聯(lián)合分析[J].長江科學(xué)院院報, 2007,24(2):13-16.

[10] 劉曾美, 陳子燊. 區(qū)間暴雨和外江洪水位遭遇組合的風(fēng)險[J].水科學(xué)進展,2009,20(5):619-625.

[11] 閆寶偉, 郭生練,陳 璐,等. 長江和清江洪水遭遇風(fēng)險分析[J].水利學(xué)報,2010,41(5):553-559.

[12] 閆寶偉, 郭生練,余 維. 長江和清江洪水過程遭遇風(fēng)險分析[J].水力發(fā)電學(xué)報,2013,32(1):50-53.

[13] 方 彬，郭生練，肖 義,等．年最大洪水兩變量聯(lián)合分布研究[J].水科學(xué)進展，2008,19(4):505-511.

[14] 程銀才,范世香,李明華.一種新的水文頻率計算方法[J].水文,2008,28(1):59-69.

[15] GUO S L.A discussion on unbiased plotting positions for the general extreme value distributeon[J].Journal of Hydrology,1990,121(1-4):33-44.