蔣 楠， 高 成， 夏 歡

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

1 研究背景

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

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

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

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

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

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

2 研究區(qū)概況

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

圖1 錦河流域水系圖

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

3 研究方法

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

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

(1)

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

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

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

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

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

3.3 Archimedean Copula函數(shù)

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

估計4種函數(shù)的參數(shù)是一個不可或缺的環(huán)節(jié)。以相關(guān)指數(shù)法來計算上述4種洪峰聯(lián)合分布函數(shù)的參數(shù)。相關(guān)指數(shù)法是依據(jù)Kendall秩相關(guān)系數(shù)τ[13]來計算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)系計算出的參數(shù)θ值見表3。

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

3.5 Archimedean Copula函數(shù)選擇

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

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

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

(6)

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

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

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

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

3.5.2 最優(yōu)Copula函數(shù)選取 為選出對錦河和連錦河洪峰遭遇擬合效果最好的Copula函數(shù)，需采用優(yōu)度評價準(zhǔn)則即AIC準(zhǔn)則及OLS準(zhǔn)則來評價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)計算出4種Copula函數(shù)的AIC值與OLS值如表4所示。

表4 AIC值與OLS值計算成果表

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

4 干支流洪水遭遇分析

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

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

=C(u1,u2)

(10)

錦河和連錦河同時發(fā)生某一設(shè)計標(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ù)量級的洪水概率如公式(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)概率隨著洪水的量級的增大而減小，聯(lián)合概率隨著洪水的量級的增大而增大。在錦河發(fā)生20年一遇的洪水工況下，連錦河發(fā)生同等量級的洪水概率明顯增加。在錦河發(fā)生一特定數(shù)量級的洪水前提下，連錦河發(fā)生洪水的概率隨著洪水量級的增加略有增加。以上規(guī)律經(jīng)驗證符合兩河洪水發(fā)生的概況。表5統(tǒng)計了100年一遇量級以下錦河和連錦河的洪水遭遇的概率。

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

表5 錦河和連錦河洪水遭遇概率計算表

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

5 水面線計算

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

5.1 模型原理

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

連續(xù)方程：

(13)

動量方程：

(14)

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

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

計算時間步長為0.16 h。

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

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

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

5.2 計算結(jié)果分析

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

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

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

6 結(jié)論與展望

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

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

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

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

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