魏琳，楊學(xué)軍，王旭東

（海河水利委員會水文局，天津 300170）

灤河是海河流域的重要河流之一，研究灤河水系的豐枯規(guī)律，分析灤河流域主要產(chǎn)流區(qū)干支流豐枯遭遇問題，對該地區(qū)防洪抗旱、雨洪資源利用、水資源的合理開發(fā)提供科學(xué)的技術(shù)支撐。因此，本文對灤河干流、潵河支流、青龍河支流3條主要干支流的豐枯遭遇問題開展分析研究，基于1956—2016年近61 a徑流量長系列水文資料，選取干流重要控制斷面潘家口、潵河支流潘大區(qū)間、青龍河支流桃林口為主要控制站（區(qū)間），采用Copula函數(shù)構(gòu)建具有相依關(guān)系的三維聯(lián)合分布模型，研究不同干支流、不同量級來水量的遭遇概率及條件概率分布，實(shí)現(xiàn)準(zhǔn)確定量豐枯遭遇的分析，為灤河流域洪水資源利用及調(diào)度決策提供技術(shù)參考。

目前，傳統(tǒng)豐枯遭遇研究方法主要是通過實(shí)測資料建立變量之間的相關(guān)關(guān)系，或者通過典型年法分析各流域來水量之間的遭遇關(guān)系［1］，難以量化不同條件下的豐枯遭遇；現(xiàn)有的水文多變量概率分析方法主要有多元正態(tài)分布法、特定邊緣分布構(gòu)成的聯(lián)合概率分布法、非參數(shù)法、將多維分布轉(zhuǎn)換成一維方法、經(jīng)驗(yàn)頻率法等，但這些方法不能準(zhǔn)確描述變量內(nèi)部的相依關(guān)系。多維Copula函數(shù)的構(gòu)建既能反映單個隨機(jī)變量分布特征，又能呈現(xiàn)變量之間的相依性，是一種可連接多個變量邊緣分布的聯(lián)合分布函數(shù)，可反映水文事件隨機(jī)性、相依性的特點(diǎn)，是解決多變量聯(lián)合分布問題的值得嘗試的工具。因此，采用Copula函數(shù)對灤河流域干支流豐枯遭遇問題進(jìn)行分析，為灤河流域?yàn)纯h以上水資源合理利用提供指導(dǎo)。

1 基于Copula方法的豐枯遭遇風(fēng)險(xiǎn)分析

1.1 Copula函數(shù)及方法

Copula理論是在1959年由Sklar提出的，可以將任意一個m維聯(lián)合累積分布函數(shù)分解為m個邊緣分布和一個Copula函數(shù)。變量的分布特性采用邊緣分布函數(shù)，變量之間的相關(guān)性采用Copula函數(shù)進(jìn)行銜接[2]，即在變量的邊緣分布函數(shù)已知的前提下存在一個Copula連接函數(shù)，形成多元聯(lián)合分布。

Copula為[0，1]之間的聯(lián)合概率分布函數(shù)。假設(shè)F為一個m維分布函數(shù)，其邊緣分布為F1，F(xiàn)2…Fm，那么存在一個m維Copula函數(shù)C，使得對任意X∈Rm有：

若F1，F(xiàn)2…Fm是連續(xù)的，則該函數(shù)C是惟一的，否則C由RanF1×RanF2×…×RanFm惟一確定；相反，若C是一個m維Copula，F(xiàn)1，F(xiàn)2…Fm為一元分布函數(shù)，則由式（1）定義的函數(shù)F是一個m元分布函數(shù)，其邊緣分布為F1，F(xiàn)2…Fm。

1.2 常用的Copula函數(shù)

Copula函數(shù)總體上可以劃分為3類，即橢圓型、二次型、Archimedean型。目前，水文上常用的Archimedean Copula函數(shù)族應(yīng)用最為廣泛，三維非對稱型Copula函數(shù)[3]主要有以下3種。

（1）Gumbel Copula函數(shù)，表達(dá)式為：

（2）Clayton Copula函數(shù)，表達(dá)式為：

（3）Frank Copula函數(shù)，表達(dá)式為：

式中：C(u1,u2,u3)為三維的Copula函數(shù)；u1,u2,u3分別為邊緣分布函數(shù)，u1=F1(x1),u2=F2(x2),u3=F3(x3)；θ1,θ2分別為Copula函數(shù)的參數(shù)。

1.3 重現(xiàn)期及遭遇分析

能夠?qū)搅髁控S枯遭遇的概率及條件概率的準(zhǔn)確分析是該方法的明顯優(yōu)勢。假定X1，X2和X3分別為具有相依關(guān)系的變量序列，事件（x1，x2，x3）的聯(lián)合概率分布F可表示為[4]：

相應(yīng)的聯(lián)合重現(xiàn)期可表示為：

式中：F(x1,x2,x3)為聯(lián)合概率分布函數(shù)，x1，x2，x3分別代表X1，X2，X3的取值；P為不超過概率；T(x1,x2,x3)為聯(lián)合重現(xiàn)期。

當(dāng)水文事件相互獨(dú)立時(shí)，聯(lián)合概率分布函數(shù)可表示為：

相應(yīng)的重現(xiàn)期可表示為：

式中：F1(x1),F2(x2),F3(x3)分別為X1、X2和X3的邊緣分布函數(shù)；T為變量相互獨(dú)立時(shí)的重現(xiàn)期；T(x1),T(x2),T(x3)分別為X1，X2和X3的重現(xiàn)期；其余變量含義同上。

1.4 參數(shù)估計(jì)和擬合優(yōu)度檢驗(yàn)

Copula函數(shù)的參數(shù)估計(jì)主要有兩種類型的方法，一種是非參數(shù)法，另一種是參數(shù)法。其中，非參數(shù)法要求有明確的Kendall的τ與Copula參數(shù)θ的表達(dá)式，一般二維Copula函數(shù)參數(shù)估計(jì)用得較多，對資料長度也有一定要求；參數(shù)法相對比較靈活，本文三維Copula函數(shù)的參數(shù)確定采用參數(shù)法中IFM（Inference of Functions for Margins）方法，該方法主要包括兩方面內(nèi)容：①采用極大似然法估計(jì)邊際分布中的參數(shù)；②采用極大似然法估計(jì)Copula函數(shù)中的兩個參數(shù)。

最優(yōu)函數(shù)的選取采用擬合優(yōu)度檢驗(yàn)，先假定不同的Copula函數(shù)，采用不同檢驗(yàn)方法進(jìn)行擇優(yōu)。目前，檢驗(yàn)方法有多種，包括Kolmogorov-Smirnov檢驗(yàn)（K-S檢驗(yàn)）、偏差、均方根誤差、AIC準(zhǔn)則等，本文采用常用的均方根誤差進(jìn)行擬合優(yōu)度檢驗(yàn)。

2 研究區(qū)概況

灤河發(fā)源于河北省豐寧縣巴彥圖古爾山麓，始稱閃電河，流經(jīng)內(nèi)蒙古自治區(qū)正蘭旗至大河口納吐力根河后稱大灤河，至河北省隆化縣郭家屯附近與小灤河匯合后稱灤河，郭家屯以下至潘家口河道蜿蜒曲折，穿行于燕山山脈，過桑園峽口進(jìn)入遷安盆地，至灤縣城關(guān)流出燕山山脈，于樂亭縣兜網(wǎng)鋪入渤海。灤河干流與伊遜河匯合處以上為灤河上游，多為壩上草原，地勢平緩，支流少；灤河三道河子以下至干流灤縣以上為灤河中游，地處燕山腹地，地形起伏較大，山嶺重疊，溝谷縱橫，區(qū)間有6條較大支流匯入，均屬山溪性河流，河道比降大，水流速度快。在干流和支流青龍河建有潘家口、大黑汀和桃林口3座大型水庫。燕山迎風(fēng)坡是暴雨多發(fā)地帶，也是灤河洪水的主要發(fā)生區(qū)和徑流產(chǎn)流區(qū)。因此，選取灤河中游潵河、青龍河及灤河干流為研究對象，開展干支流豐枯遭遇分析。

3 灤河干支流豐枯遭遇分析

3.1 選取干支流特性分析

灤河干流選取潘家口水庫站來水、潵河選用潘大區(qū)間來水、青龍河選取桃林口水庫站來水作為代表系列，包括1956—2016年共61 a數(shù)據(jù)，假定常用水文變量頻率分布為P-Ⅲ型分布，采用現(xiàn)狀參數(shù)估計(jì)較為準(zhǔn)確的線性矩法進(jìn)行初步估計(jì)并進(jìn)行適線，各河流特征值及不同概率下的來水量詳見表1。

表1 灤河干支流代表站特征值及特征來水量（不超過概率）

由表1可以看出，灤河干流多年平均徑流量為15.58億m3，青龍河為6.86億m3，潵河為2.78億m3，其中支流青龍河多年平均徑流量大于潵河。豐水年灤河干流徑流量為19.72億m3，青龍河為9.07億m3，潵河為3.72億m3；枯水年灤河干流徑流量為8.3億m3，青龍河為2.92億m3，潵河為1.24億m3。

3.2 Copula函數(shù)的最優(yōu)選取

3.2.1 邊緣分布函數(shù)選取

日常工作中，常采用P-Ⅲ型分布作為水文變量的分布線型。由于不同的水文氣象特征，各地區(qū)來水量系列也可能服從其他分布，因此我們假定P-Ⅲ型分布、對數(shù)正態(tài)分布、Gumbel分布、極值分布4種常用的分布線型，采用均方根誤差、絕對值誤差、AIC信息準(zhǔn)則、PPCC檢驗(yàn)法4種擬合準(zhǔn)則進(jìn)行線型擬合，選出擬合效果最好的分布線型，其中均方根誤差（RMSE）、絕對值誤差（ABS）及AIC信息準(zhǔn)則越小越好，PPCC越接近于1越優(yōu)，結(jié)果詳見表2。

表2 邊際分布的擬合優(yōu)度檢驗(yàn)

由表2來看，對潵河潘大區(qū)間而言，P-Ⅲ分布具有最小的ABS和RMSE值，表明P-Ⅲ分布為該站最優(yōu)線型；對灤河干流潘家口水庫站而言，P-Ⅲ分布具有最小的ABS和RMSE值，表明P-Ⅲ分布為該站最優(yōu)線型；對青龍河桃林口水庫站而言，P-Ⅲ分布具有最小的ABS、RMSE和AIC值，表明P-Ⅲ分布為該站最優(yōu)線型。經(jīng)綜合考慮，選定P-Ⅲ分布作為研究潘家口、桃林口、潘大區(qū)間來水量數(shù)據(jù)的分布線型。

3.2.2 相關(guān)性分析

不同的Copula函數(shù)有各自的特性，其中Gumbel Copula函數(shù)和Clayton Copula函數(shù)適用于具有正相關(guān)關(guān)系的變量；Frank Copula函數(shù)對相關(guān)關(guān)系為正、負(fù)的變量均可；Ali-Mikhail-Haq Copula函數(shù)適用于具有弱相關(guān)關(guān)系的變量。因此，根據(jù)變量之間的相關(guān)關(guān)系，初選適宜的Copula函數(shù)。本文采用線性相關(guān)系數(shù)及秩相關(guān)系數(shù)，對變量之間的相關(guān)性進(jìn)行分析計(jì)算，結(jié)果詳見表3。

表3 相關(guān)性計(jì)算結(jié)果

由表3可以看出，三站徑流量數(shù)據(jù)之間都屬于正相關(guān)關(guān)系，從Copula函數(shù)對變量相關(guān)關(guān)系的要求來看，Gumbel Copula、Clayton Copula及Frank Copula函數(shù)適用于該研究，因此對這3種Copula函數(shù)開展擬合優(yōu)度檢驗(yàn)。

3.2.3 Copula函數(shù)參數(shù)確定及擬合優(yōu)度檢驗(yàn)

為檢驗(yàn)Copula函數(shù)的擬合精度，需分析比較各個數(shù)據(jù)點(diǎn)對（xi,yi,zi）的經(jīng)驗(yàn)累積概率和理論累積概率的一致性[3]。設(shè)N為水文系列觀測樣本對數(shù)，按照x升序排列所得數(shù)據(jù)系列為(x1,y1,z1),(x2,y2,z2)…(xi,yi,zi)…(xN,yN,zN)，三維經(jīng)驗(yàn)聯(lián)合概率分布采用下式計(jì)算：

根據(jù)極大似然法計(jì)算各Copula函數(shù)的參數(shù)[5]，結(jié)果詳見表4。根據(jù)式（9）計(jì)算灤河干流、潵河、青龍河的經(jīng)驗(yàn)聯(lián)合頻率值，通過Copula函數(shù)計(jì)算三變量的理論聯(lián)合頻率，將所得的經(jīng)驗(yàn)聯(lián)合頻率與理論聯(lián)合頻率點(diǎn)繪在圖1中。由圖1可以看出，經(jīng)驗(yàn)頻率點(diǎn)和理論頻率點(diǎn)均勻分布在45°線附近，表明三變量的經(jīng)驗(yàn)頻率點(diǎn)和理論頻率點(diǎn)擬合較好，說明選擇的聯(lián)合分布函數(shù)均是合理的。由于Gumbel、Frank、Clayton Copula函數(shù)中采用Frank Copula函數(shù)擬合最好，選取Frank Copula為最優(yōu)函數(shù)。

圖1 灤河干流、潵河、青龍河年來水量經(jīng)驗(yàn)頻率與不同理論頻率的擬合

為確定灤河干流、潵河、青龍河來水量擬合最好的Copula函數(shù)，需要進(jìn)行擬合優(yōu)度檢驗(yàn)，本文采用均方根誤差（RMSE）檢驗(yàn)方法，結(jié)果詳見表4。根據(jù)結(jié)果可以看出，初選的3種Copula函數(shù)中，F(xiàn)rank Copula具有最小的RMSE值，表明Frank Copula具有最好的擬合優(yōu)度，對本文選用的灤河干流、潵河、青龍河年來水量數(shù)據(jù)而言，F(xiàn)rank Copula即為選定的最優(yōu)Copula函數(shù)。

表4 三維Copula函數(shù)擬合優(yōu)度檢驗(yàn)結(jié)果

3.3 干支流豐枯遭遇分析

從分析得到的灤河干支流徑流量的邊緣分布和最優(yōu)Copula函數(shù)，可得到各種聯(lián)合分布和條件概率分布，從而分析不同情況下干支流豐枯遭遇情況，如灤河干流、潵河支流、青龍河支流的年來水量保證率都能達(dá)到75%的概率；或者當(dāng)灤河干流潘家口來水量不足時(shí)，遭遇支流潵河及青龍河桃林口水庫豐水的概率；或者當(dāng)干流潘家口為豐水年，遭遇潵河及青龍河枯水年份的概率。通過三維Copula函數(shù)的構(gòu)建，這種問題將迎刃而解[2]。

（1）當(dāng)灤河干流、潵河支流、青龍河支流的來水量達(dá)到某一概率水平的聯(lián)合概率分布表示如下：

式中：θ1=8.392 3,θ2=8.398 6，為Frank Copula函數(shù)的參數(shù)；X，Y，Z依次為桃林口水庫站，潘大區(qū)間，潘家口水庫站來水量數(shù)據(jù)系列；U1=Fx（x），U2=FY（y），U3=Fz（z）分別為各站的邊際分布函數(shù)，u1、u2、u3為其取值。

由式（10）可得灤河干流、潵河、青龍河同時(shí)出現(xiàn)枯水年的概率為F(x,y,z)=C(25%,25%,25%)=P(X≤2.92,Y≤1.24,Z≤8.3)=0.134 2。

灤河干流豐水年遭遇支流潵河、青龍河同時(shí)出現(xiàn)枯水年的概率為F(x,y,z)=C(25%,25%,75%)=P(X≤2.92,Y≤1.24,Z≤19.72)=0.174 4。

因此，灤河干流、潵河、青龍河同時(shí)出現(xiàn)枯水年的概率為0.13，灤河干流豐水年遭遇支流潵河、青龍河同時(shí)出現(xiàn)枯水年的概率為0.17。也就是說，3條河同時(shí)出現(xiàn)枯水的概率不大，灤河干流豐水年兩條支流枯水年出現(xiàn)的概率也不大。

（2）假定某一干支流的來水量水平已知，需要據(jù)此預(yù)測其他干支流的來水量時(shí)條件概率出現(xiàn)的情況，本文亦可以準(zhǔn)確計(jì)算得出。各種水平年遭遇組合的條件概率如下。

已知灤河干流為枯水年時(shí)預(yù)測支流潵河及青龍河同時(shí)發(fā)生枯水年的概率為F(y,x|z)=P(Y≤1.24,

已知灤河干流為特豐水年時(shí)預(yù)測支流潵河及青龍河同時(shí)發(fā)生枯水年的概率為F(y,x|z)=P(Y≤1.24,=0.184 2，即灤河干流出現(xiàn)枯水年份兩條支流同時(shí)發(fā)生枯水年的概率為50%，灤河干流為特豐水年支流發(fā)生枯水年的概率為0.18。

圖2和圖3為灤河干流潘家口水庫分別為枯水年（P=25%）和特豐水年（P=95%）時(shí)遭遇潵河、青龍河不同量級水量的概率分布，由此可以簡單直觀得到相應(yīng)的條件概率，為灤河流域水資源的決策分析提供技術(shù)支撐。

圖2 灤河干流枯水年（P=25%）條件下潵河及青龍河條件概率曲線

圖3 灤河干流特豐年（P=95%）條件下潵河及青龍河條件概率曲線

4 結(jié)論及展望

（1）三維Copula函數(shù)通過任意的邊緣分布和相關(guān)性結(jié)構(gòu)來構(gòu)造多維聯(lián)合概率分布，在水文豐枯遭遇的應(yīng)用上改變了傳統(tǒng)方式對該問題難以準(zhǔn)確量化的缺陷，該函數(shù)方法形式靈活多樣，考慮了變量間的相關(guān)關(guān)系，為研究多維聯(lián)合概率分布及豐枯遭遇提供了新的有效手段，可在生產(chǎn)實(shí)踐中應(yīng)用。

（2）針對灤河干支流豐枯遭遇的分析，采用三維Frank Copula函數(shù)擬合最好，從而可以得到不同量級下的干支流豐枯條件下的概率和重現(xiàn)期，為灤河流域洪水資源利用及科學(xué)管理提供必要的理論依據(jù)參考。

（3）基于Copula函數(shù)的三變量聯(lián)合分布的來水頻率分析方法能更好地描述來水特征量之間的關(guān)系，且對邊緣分布類型沒有限制，在現(xiàn)有實(shí)際來水頻率分析計(jì)算中有較大的優(yōu)越性，但未來仍需對參數(shù)估計(jì)、函數(shù)選擇進(jìn)一步開展研究，使得模型更加完善。