基于Copula函數(shù)的錦河與連錦河洪水遭遇分析

蔣 楠， 高 成， 夏 歡

(1.河海大學(xué) 水文水資源學(xué)院， 江蘇 南京 210098； 2.河海大學(xué) 農(nóng)業(yè)工程學(xué)院， 江蘇 南京 210098)

1 研究背景

洪水發(fā)生及產(chǎn)生的災(zāi)害是難以避免的自然現(xiàn)象。作為能夠有效抵御洪水的手段之一，防洪工程措施在防治大洪水中起到了不可小覷的作用。在擬定防洪方案過程中，洪水計(jì)算是不可或缺的。因此在進(jìn)行防洪評(píng)價(jià)的洪水計(jì)算部分有必要進(jìn)行洪水遭遇的概率分析。利用Copula函數(shù)進(jìn)行多變量水文概率分析已十分常見。例如：王占海等[1]使用多維Copula函數(shù)研究長(zhǎng)江干流宜昌水文站各設(shè)計(jì)標(biāo)準(zhǔn)洪水空間分布的合理性。侯蕓蕓等[2]在基于三變量Copula函數(shù)算法上建立了峰量以及歷時(shí)的聯(lián)合概率分布模型。陳璐等[3]利用多維Archimedean Copula函數(shù)計(jì)算長(zhǎng)江干支流沿岸水文站洪水發(fā)生時(shí)間及量級(jí)遭遇的概率，李子遠(yuǎn)[4]采用Gumbel Copula函數(shù)模擬兩河洪峰發(fā)生時(shí)間的雙變量聯(lián)合分布，以此計(jì)算洪峰遭遇概率。葉姍姍等[5]運(yùn)用Copula函數(shù)風(fēng)險(xiǎn)概率模型對(duì)宿遷市暴雨雨型風(fēng)險(xiǎn)率進(jìn)行了研究；高玉琴等[6]采用三維G-H Copula函數(shù)對(duì)秦淮河流域的洪峰和洪水位進(jìn)行了風(fēng)險(xiǎn)率的計(jì)算；王亞雄等[7]用最優(yōu)Copula函數(shù)進(jìn)行了飛來峽水庫壩址洪峰和最大7日洪量的聯(lián)合重現(xiàn)期與邊緣分布重現(xiàn)期的比較；Salvadori等[8]運(yùn)用二維的Copula函數(shù)對(duì)降雨歷時(shí)和降雨強(qiáng)度的聯(lián)合概率分布進(jìn)行模擬。Poulin等[9]分析了概率模型的多種邊緣分布函數(shù)特性，并將之應(yīng)用于實(shí)例。Hochrainer-stigler等[10]對(duì)Copula函數(shù)的原理結(jié)構(gòu)進(jìn)行了研究并用于系統(tǒng)整合風(fēng)險(xiǎn)分析。但以上的Copula函數(shù)研究還存在一定局限性，主要為以下幾點(diǎn)：

(1)在水文頻率計(jì)算中，經(jīng)驗(yàn)頻率法是基于實(shí)測(cè)資料進(jìn)行經(jīng)驗(yàn)公式計(jì)算的一種方法，它需要足夠長(zhǎng)的實(shí)測(cè)系列資料才能滿足計(jì)算精度，其次，該方法在外延預(yù)測(cè)方面效果不佳。

(2)大多數(shù)水文變量服從偏態(tài)分布，如果假設(shè)這些變量屬于正態(tài)分布會(huì)使得計(jì)算結(jié)果的精度不高。

(3)部分地區(qū)有高潮出現(xiàn)時(shí)，其降雨類型復(fù)雜，Copula函數(shù)在雨潮組合遭遇的風(fēng)險(xiǎn)分析方面還有待進(jìn)一步研究。

(4)Copula函數(shù)在干支流特大洪水同時(shí)發(fā)生的分布模型還需要探索。

在以上研究基礎(chǔ)上，本文主要研究應(yīng)用Copula函數(shù)計(jì)算錦河和連錦河洪水遭遇概率，收集江西省高安和宜豐兩個(gè)水文站的實(shí)測(cè)洪峰流量資料，采用年最大值法選取年最大流量，并計(jì)算其經(jīng)驗(yàn)頻率，其邊緣分布為P-Ⅲ型曲線。然后選擇4種常用的Archimedean Copula函數(shù)對(duì)錦河和連錦河洪水遭遇的聯(lián)合分布概率進(jìn)行研究，經(jīng)相關(guān)參數(shù)估算后，以AIC[11]和OLS值最小的原則來選出模擬效果最佳的Archimedean Copula函數(shù)，將其應(yīng)用于計(jì)算錦河和連錦河洪水的同現(xiàn)概率及條件概率，對(duì)計(jì)算進(jìn)行分析研究，并可將洪水遭遇概率及河道水面線的計(jì)算成果應(yīng)用于實(shí)際工程項(xiàng)目中。

2 研究區(qū)概況

錦河流域位于江西省高安市西部。本文研究?jī)?nèi)容為流域內(nèi)的錦河與連錦河的洪水遭遇情況。錦河起源于湘贛交界的幕阜山脈海拔高程為628.6 m的東麓坪子嶺。全長(zhǎng)304.55 km，流域面積7 884 km2，是贛江下游西岸的一條重要支流，也是流經(jīng)高安的第一大河流。錦河干流城區(qū)段河長(zhǎng)約13 km，寬約150～300 m，河床底高程20～24.5 m，較河岸低8.0～12.0 m，河道縱坡約為(0.6～2.0)/10000。其流域水系圖見圖1。

圖1 錦河流域水系圖

連錦河，又名港溪水，位于錦河干流左岸，為高安城區(qū)主要支流之一，其河道左岸地勢(shì)較高，右岸城區(qū)段的地勢(shì)低洼。該河流發(fā)源于高安市奉新縣交界的米嶺南麓，自北向南流經(jīng)倉下村、汪家圩、港溪、水口庵，至連錦橋進(jìn)入城區(qū)，于五口橋匯入錦河。流域面積89.85 km2，主河道全長(zhǎng)23.01km，主河道縱比降1.91‰，高安城市規(guī)劃內(nèi)連錦河的河段長(zhǎng)約5.5 km，一般河寬50～80 m。

3 研究方法

3.1 邊緣分布函數(shù)計(jì)算

Copula函數(shù)計(jì)算概率之前首先要確定水文變量的邊緣分布函數(shù)。P-Ⅲ型曲線能夠較好地?cái)M合暴雨洪水等變量的概率分布，是國(guó)內(nèi)一種應(yīng)用較為廣泛的邊緣分布函數(shù)。本文給定高安站及宜豐站的洪峰流量服從P-Ⅲ分布函數(shù)[12]，其概率密度函數(shù)如公式(1)所示。

(1)

錦河流域內(nèi)有重要的水文測(cè)站高安站，其實(shí)測(cè)洪水能較好地代表錦河洪水過程。該水文站位于高安市城區(qū)高安大橋上游200 m處，控制集水面積為6 215 km2，有1964、1991-2016年歷年實(shí)測(cè)洪峰流量資料，經(jīng)附近賈村站插補(bǔ)延展有1964年至2016年洪峰流量資料。宜豐水文站位于錦河左岸支流連錦河，控制流域面積為519 km2，河道坡降為1.91‰，現(xiàn)有1957-2016年實(shí)測(cè)洪峰流量資料。資料條件系列長(zhǎng)，數(shù)據(jù)精度高，可作為連錦河水文站，由于連錦河上游有一碧山水庫，其上游集水面積為31.7 km2，壩址以上河長(zhǎng)8.70 km，平均坡降為5.58‰，水庫正常蓄水位58.01 m，總庫容1 979×104m3，考慮其調(diào)蓄作用，進(jìn)行了調(diào)洪計(jì)算，在邊緣分布計(jì)算中洪峰流量將減去其調(diào)蓄流量。對(duì)宜豐站、高安站流量進(jìn)行排頻分析，并將高安站1913年調(diào)查最大流量作特大值處理，采用P-Ⅲ曲線適線，可進(jìn)行洪水頻率計(jì)算并得到錦河和連錦河P-Ⅲ分布函數(shù)的統(tǒng)計(jì)參數(shù)見表1。

表1 錦河和連錦河P-Ⅲ分布函數(shù)統(tǒng)計(jì)參數(shù)值表

3.2 邊緣分布函數(shù)的檢驗(yàn)

水文變量的理論分布是否能夠準(zhǔn)確代表變量總體的分布，需要進(jìn)行假設(shè)檢驗(yàn)。為保證檢驗(yàn)的準(zhǔn)確性，本文采用Kolmogorov-Smirnov檢驗(yàn)、Anderson-darling檢驗(yàn)以及均方根誤差(RMSE)來分別衡量高安站和宜豐站洪峰流量的邊緣分布函數(shù)。

將兩個(gè)水文站的洪峰流量系列導(dǎo)入SPSS中進(jìn)行K-S非參數(shù)檢驗(yàn)，得出兩個(gè)測(cè)站的P值均為0.2，大于0.05，表明該邊緣分布通過K-S檢驗(yàn)。Anderson-darling檢驗(yàn)中計(jì)算出的高安站與宜豐站的P值分別為0.26和0.29，大于0.05，該分布通過了Anderson-darling檢驗(yàn)。計(jì)算得出兩個(gè)測(cè)站的RMSE分別為0.025和0.022，表明二者的均方根誤差均較小。根據(jù)上述檢驗(yàn)，可以確定兩站的P-Ⅲ分布可作為邊緣分布函數(shù)。

3.3 Archimedean Copula函數(shù)

Archimedean Copula函數(shù)可用于構(gòu)建水文界常見的聯(lián)合概率分布模型，它的結(jié)構(gòu)形式多種多樣，不同結(jié)構(gòu)對(duì)不同類型的洪水?dāng)M合效果有較大差異，本文選用以下4種對(duì)稱式的結(jié)構(gòu)形式。

(1)Gumbel-Hougaard Copula[7]：

C(u1,u2)=exp{-[(-lnu1)θ+

(-lnu2)θ](1/θ)}θ∈[1,+∞)

(2)

(2)Clayton Copula[7]：

C(u1,u2)=(u1-θ+u2-θ-1)-1/θ

(3)

θ∈(0,+∞)

(3)Ali-Mikhail-Haq Copula[7]：

C(u1,u2)=u1u2/[1-θ(1-u1)(1-u2)]

(4)

θ∈[1,+∞)

(4)Frank Copula[7]：

(5)

θ∈R

式中：C為Copula函數(shù)；u1和u2為邊緣分布函數(shù)，u1=F1(x1)，u2=F2(x2)；θ為相關(guān)參數(shù)。

3.4 Archimedean Copula函數(shù)參數(shù)估計(jì)

估計(jì)4種函數(shù)的參數(shù)是一個(gè)不可或缺的環(huán)節(jié)。以相關(guān)指數(shù)法來計(jì)算上述4種洪峰聯(lián)合分布函數(shù)的參數(shù)。相關(guān)指數(shù)法是依據(jù)Kendall秩相關(guān)系數(shù)τ[13]來計(jì)算Copula函數(shù)的參數(shù)θ，兩參數(shù)的函數(shù)關(guān)系式見表2。

表2 Copula函數(shù)的τ與θ函數(shù)關(guān)系[14]表

根據(jù)SPSS軟件分析得出兩站洪峰流量Kendall秩相關(guān)系數(shù)τ為0.174，根據(jù)函數(shù)關(guān)系計(jì)算出的參數(shù)θ值見表3。

表3 Copula函數(shù)的θ計(jì)算成果表

3.5 Archimedean Copula函數(shù)選擇

由上節(jié)計(jì)算出函數(shù)的θ，可建立兩水文站洪峰流量的4種Archimedean Copula函數(shù)的聯(lián)合概率分布，以上函數(shù)是否能夠較為準(zhǔn)確地?cái)M合兩變量之間的關(guān)系，需要進(jìn)行函數(shù)擬合檢驗(yàn)，然后通過優(yōu)度評(píng)價(jià)準(zhǔn)則選擇出一種擬合效果最好的Copula函數(shù)。

3.5.1 圖形評(píng)價(jià)分析法 本文主要利用Q-Q圖的方法直觀表示擬合效果的優(yōu)劣。將經(jīng)驗(yàn)和理論聯(lián)合概率值點(diǎn)據(jù)呈現(xiàn)在一張圖上，若各點(diǎn)大致位于y=x線左右，則說明構(gòu)建的模型有一定可靠度。二維經(jīng)驗(yàn)聯(lián)合概率計(jì)算見公式(6)：

Femp(xi1,xi2)=P(X1≤xi1,X2≤xi2)

(6)

式中：Femp為經(jīng)驗(yàn)聯(lián)合概率；ng,k為滿足X1≤xi1,X2≤xi2的聯(lián)合觀測(cè)數(shù)量；n為系列長(zhǎng)度。

將4種Copula函數(shù)的理論與經(jīng)驗(yàn)聯(lián)合概率值點(diǎn)繪為散點(diǎn)圖，見圖2。

圖2 4種Copula函數(shù)的理論與經(jīng)驗(yàn)值Q-Q圖

根據(jù)Q-Q圖可見以上4種Copula函數(shù)的理論與經(jīng)驗(yàn)聯(lián)合概率點(diǎn)據(jù)均位于45°線上下，則可說明這4種Copula函數(shù)構(gòu)建的聯(lián)合概率分布模型的模擬效果較為滿意。

3.5.2 最優(yōu)Copula函數(shù)選取 為選出對(duì)錦河和連錦河洪峰遭遇擬合效果最好的Copula函數(shù)，需采用優(yōu)度評(píng)價(jià)準(zhǔn)則即AIC準(zhǔn)則及OLS準(zhǔn)則來評(píng)價(jià)4種Copula函數(shù)。

AIC準(zhǔn)則結(jié)構(gòu)如下所示：

C(ui1,ui2,…,uim))2

(7)

AIC=nln(MSE)+2k

(8)

式中:m為函數(shù)的維數(shù);k為模型參數(shù)的數(shù)量。AIC值越小，表明函數(shù)的聯(lián)合概率分布擬合效果越佳。

OLS準(zhǔn)則結(jié)構(gòu)如下：

OLS=

(9)

根據(jù)公式(8)與(9)計(jì)算出4種Copula函數(shù)的AIC值與OLS值如表4所示。

表4 AIC值與OLS值計(jì)算成果表

由表4可知AIC值與OLS值最小的均為Frank Copula函數(shù)，則可說明Frank Copula函數(shù)是模擬該流域洪水洪峰遭遇概率效果最佳的函數(shù)，因此，本文選擇Frank Copula函數(shù)進(jìn)行錦河與連錦河洪水遭遇概率分析。

4 干支流洪水遭遇分析

洪峰流量(x1,x2)的聯(lián)合分布為F(x1,x2)，則：

F(x1,x2)=P(X1≤x1,X2≤x2)

=C(u1,u2)

(10)

錦河和連錦河同時(shí)發(fā)生某一設(shè)計(jì)標(biāo)準(zhǔn)的洪水的概率為重現(xiàn)概率，以H0表示，則:

H0(x1,x2)=P(X1>x1,X2>x2)=1-u1-u2+C(u1,u2)

(11)

在錦河發(fā)生一種重現(xiàn)期的洪水前提下，連錦河發(fā)生某一數(shù)量級(jí)的洪水概率如公式(12)：

(12)

錦河和連錦河的洪峰聯(lián)合概率分布、同現(xiàn)概率與條件概率如圖3、4、5所示。

圖中Qxi為錦河的洪峰流量(m3/s)，Qyi為連錦河的洪峰流量(m3/s)。由圖4可知，錦河的洪峰流量大于5 000 m3/s(約100年一遇)和連錦河的洪峰流量大于1 200 m3/s(約100年一遇)的洪水發(fā)生的概率極小。錦河和連錦河的洪水同現(xiàn)概率隨著洪水的量級(jí)的增大而減小，聯(lián)合概率隨著洪水的量級(jí)的增大而增大。在錦河發(fā)生20年一遇的洪水工況下，連錦河發(fā)生同等量級(jí)的洪水概率明顯增加。在錦河發(fā)生一特定數(shù)量級(jí)的洪水前提下，連錦河發(fā)生洪水的概率隨著洪水量級(jí)的增加略有增加。以上規(guī)律經(jīng)驗(yàn)證符合兩河洪水發(fā)生的概況。表5統(tǒng)計(jì)了100年一遇量級(jí)以下錦河和連錦河的洪水遭遇的概率。

圖3錦河與連錦河聯(lián)合概率分布圖 圖4錦河與連錦河同現(xiàn)概率分布圖 圖5錦河與連錦河條件概率分布圖

表5 錦河和連錦河洪水遭遇概率計(jì)算表

由表5可知，錦河和連錦河同量級(jí)洪水同時(shí)出現(xiàn)的概率較小，在錦河出現(xiàn)100年一遇工況的洪水時(shí)，連錦河發(fā)生同一設(shè)計(jì)標(biāo)準(zhǔn)洪水的幾率僅有2.662%。兩河發(fā)生20年一遇量級(jí)洪水的概率較大，在錦河發(fā)生10年一遇以上洪水的情況下，連錦河同時(shí)發(fā)生20年一遇以上洪水的概率高達(dá)23.01%。由此可見錦河和連錦河大量級(jí)洪水遭遇的可能性較小，而小量級(jí)洪水遭遇的機(jī)率較大，該結(jié)論與《贛江流域水文年鑒》中近幾十年的洪水統(tǒng)計(jì)資料相符合。

5 水面線計(jì)算

根據(jù)水文年鑒上的洪水實(shí)測(cè)資料，估算連錦河洪水傳播時(shí)間約為2.8 h。錦河上高水文站到高安站距離為79 km，洪水傳播時(shí)間約為18 h。因此，錦河干流與連錦河支流洪峰在時(shí)間上是錯(cuò)開的，同時(shí)遭遇的幾率甚小，所以不按同頻遭遇考慮。工程研究的河段對(duì)于錦河而言是從徐王村到筠州大橋，樁號(hào)為-K5+910～K1+523，連錦河的研究河段為平安大道到連錦河河口，樁號(hào)為-K2+966～K2+027。

5.1 模型原理

MIKE11一維河網(wǎng)匯流模型采用非恒定流圣維南(Saint-Venant)方程組構(gòu)建[19]。

連續(xù)方程：

(13)

動(dòng)量方程：

(14)

式中：A為河道過水面積，m2；Q為流量，m3/s；t為時(shí)間，h；x為與水流方向一致的橫向坐標(biāo)，m；q為河道旁側(cè)入流，m3/s；g為重力加速度，m/s2；h為水位，m。

為了推求高安城區(qū)錦河和連錦河的水面線，將對(duì)河道進(jìn)行概化[20]，利用MIKE11一維水動(dòng)力模型對(duì)錦河10年一遇遭遇連錦河20年一遇洪水情況進(jìn)行模擬，并計(jì)算得出該工況下的洪水下水面線。將河網(wǎng)概化，生成河網(wǎng)文件，編輯河道參數(shù)(長(zhǎng)度，斷面)，輸入初始條件和邊界條件。

計(jì)算時(shí)間步長(zhǎng)為0.16 h。

邊界條件：錦河上游給定10年一遇的流量時(shí)間序列，連錦河上游給定20年一遇的流量時(shí)間序列，錦河下游河口位置給定水位流量關(guān)系曲線。

根據(jù)《水力計(jì)算手冊(cè)(第二版)》對(duì)糙率取值要求，結(jié)合計(jì)算區(qū)域土地利用等情況，需要對(duì)錦河、連錦河一維河道模型糙率進(jìn)行選取和率定。根據(jù)歷史相關(guān)資料和經(jīng)驗(yàn)，各河段糙率選取0.033。

模型斷面迭代的初始水深為1.84 m。

5.2 計(jì)算結(jié)果分析

由第4節(jié)可知，錦河發(fā)生10年一遇以上洪水遭遇連錦河發(fā)生20年一遇以上洪水的概率較大，屬于河道工程常見工況，水利模型構(gòu)建必不可少。因此用MIKE11模型計(jì)算錦河與連錦河這種工況的洪水水面線。如表6所示。

由表6可知，當(dāng)錦河出現(xiàn)10年一遇的洪水時(shí)，連錦河出現(xiàn)20年一遇情況的水位較高，這種工況的洪水同期遭遇概率相對(duì)較大，結(jié)合實(shí)際詢問調(diào)查，得出連錦河河口段高水位往往是受錦河洪水頂托影響的結(jié)論，該工況下的水面線計(jì)算結(jié)果對(duì)高安市的防洪工程具有重大意義。

表6 錦河10年一遇遭遇連錦河20年一遇洪水水面線表

6 結(jié)論與展望

本文以高安市錦河與連錦河的洪峰流量資料為研究對(duì)象，采用Copula 函數(shù)建立了兩水文站洪峰的同現(xiàn)概率分布與條件概率分布來研究干支流的洪水遭遇問題，對(duì)高安市的防洪評(píng)價(jià)工程有重要意義。得出以下4點(diǎn)結(jié)論：

(1)通過邊緣分布函數(shù)的K-S和A-D檢驗(yàn)，錦河和連錦河的洪峰分布服從P-Ⅲ型分布。

(2)Archimedean Copula函數(shù)的4種形式對(duì)洪峰的聯(lián)合概率分布均能較好地?cái)M合，通過優(yōu)度評(píng)價(jià)準(zhǔn)則評(píng)價(jià)后，得出擬合效果最好的函數(shù)是Frank Copula函數(shù)的結(jié)論。

(3)錦河和連錦河的洪水同現(xiàn)概率隨著洪水的量級(jí)的增大而減小，聯(lián)合概率隨著洪水的量級(jí)的增大而增大。在錦河發(fā)生某一重現(xiàn)期的洪水情況下，連錦河發(fā)生洪水的概率隨著洪水量級(jí)的增加略有上升。

(4)錦河和連錦河大量級(jí)洪水遭遇的可能性較小，而小量級(jí)洪水遭遇的幾率較大，錦河和連錦河同時(shí)發(fā)生20年一遇的洪水概率為12.14%，在錦河發(fā)生10年一遇以上洪水的情況下，連錦河發(fā)生20年一遇以上洪水的概率為23.01%，并對(duì)該工況下的洪水水面線進(jìn)行了計(jì)算，這對(duì)錦河流域的防洪評(píng)價(jià)等工程具有重要的實(shí)踐價(jià)值。