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

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

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

1 研究區(qū)域概況

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

圖1 東江干流、望江瀝及相關站點位置示意圖

2 研究方法

2.1 Copula函數

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

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

(1)

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

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

(2)

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

采用相關指標法[4]進行參數估計，即根據Kendall秩相關系數τ與參數θ之間的關系間接求出，表達式如下：

(3)

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

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

(4)

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

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

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

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分布在混合分布中的比例。關于Von Mises分布函數參數可采用最小二乘法獲取。

2.3 年最大洪水量級分布

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

(6)

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

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

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

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

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

(7)

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

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

(8)

X和Y的同現重現期Tα：

(9)

邊緣分布的重現期：

(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的概率，其重現期為年。

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

(12)

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

(13)

3 東江干支流洪水遭遇風險分析

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

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

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

表1 Von Mises分布參數估算結果

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

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

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

圖3 年最大洪水和年最大24 h降雨出現時間經驗值與理論值P-P圖

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

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日，遭遇風險約為0.011%，6月2日至7月9日是東江干支流出現洪峰遭遇的集中時期，7月10日至8月10有一個相對平穩(wěn)的遭遇風險，約為0.006%，東江干支流在5月9日之前和10月6日之后發(fā)生遭遇的概率很小，幾乎為0。

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

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

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

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

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

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

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

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

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

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

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

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

4 結 論

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

□

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

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

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

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

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

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

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

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

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

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

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

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

[13] 方 彬，郭生練，肖 義,等．年最大洪水兩變量聯(lián)合分布研究[J].水科學進展，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.