基于雙參數Gumbel-Hougaard Copula函數的 設計洪水地區組成分析

牟時宇，石 朋,2，趙蘭蘭，崔彥萍，馮 穎，董豐成，陳 琛

(1.河海大學水文水資源學院，南京 210098；2.水文水資源與水利工程科學國家重點實驗室，南京 210098； 3.水利部信息中心，北京 100053；4.江蘇省水文水資源勘測局，南京 210029)

洪水地區組成是研究流域開發方案，上下游、干支流工程聯合調洪運行時需要分析的問題。傳統分析方法有地區組成法、頻率組合法、隨機模擬法[1]。前2者的共同缺點是難以理清各個分區洪水組成的概率問題，主觀假定可能的組成方案而削弱了洪水的空間相關性[2]。而隨機模擬法難以掌握模擬序列統計特征的穩定性。近年來，Copula函數在水文上的應用[3,4]為設計洪水地區組成研究提供了突破口。Copula函數可以靈活地構建任意多個變量之間的聯合分布函數，模擬變量間非線性關系的性質使推求更加合理的地區組成成為可能[5-7]。

然而，目前設計地區洪水組成的研究均選取單參數的Gumbel-Hougaard Copula函數(以下簡稱GH Copula)作為聯合分布函數的載體，其擬合效果受系列長短的影響較大，參數估計結果具有較大的不確定性。而多參數GH Copula函數保留了單參數GH Copula函數的上尾相關性，且因參數個數的增加，使函數結構更加靈活，能夠更好地反映極值事件的相關性。本文基于多參數GH Copula函數分別運用條件期望組合法、條件最可能組合法、最可能組合法、同頻率組成法擬定淮河流域中渡斷面的設計洪水地區組成方案?；春又猩嫌未笮退畮?、閘壩眾多，分析該區間段設計洪水地區組成對淮河流域綜合治理有重要意義。

1 研究方法

1.1 Copula函數模型

1.1.1 GH Copula函數

水文上常用的GH Copula為單參數Copula函數，其數學表達式如下：

C(u,v;θ)=exp {-[(-logu)θ+(-logv)θ1/θ}

(1)

式中：θ為參數并且滿足θ≥1。

GH Copula的上尾相關系數為2-21/θ，是Archimedean Copula中唯一可以用于多變量極值分析的Copula函數，適合刻畫具有上尾相關性的水文變量。

通過對生成元φ0(t)=(1/t)-1增加內外冪的方式可以構造雙參數GH Copula函數[8]：

C(u,v;β1,β2)={[(u-β2-1)β1+(v-β2-1)β1]1/β2+1}-1/β2

(2)

式中：β1、β2為參數，并且滿足β1≥1，β2>0，其中β2是與下尾相關性有關的參數。

雙參數GH Copula函數更具靈活性和外延性，對于相關性一般的變量而言更具優勢。

在單參數GH Copula函數的基礎上引入了Marshall-Olkin Copula函數可構造三參數GH Copula函數[9]：

C(u,v;θ,π2,π3)=exp [-A(-logu,-logv;θ,π2,π3)]

(3)

A(x,y;θ,π2,π3)=

[(π2x)θ+(π3y)θ]1/θ+(1-π2)x+(1-π3)y

(4)

式中：θ≥1，0≤π2,π3≤1。

其尾部相關系數與單參數GH Copula函數相同。三參數GH Copula函數是非對稱Copula函數，用于擬合概率分布具有非對稱性的二維變量。

1.1.2 Copula函數的擬合優度評價

Copula函數的擬合優度評價分為3類方法：

(1)利用離差平方和(OLS)最小準則、AIC準則從整體上評價Copula函數擬合效果：

(5)

(6)

AIC=nlnMSE+2m

(7)

式中：Femp(xi,yi)為經驗頻率；C(ui,vi)為理論頻率；m為模型參數的個數。

(2)通過繪制圖形的方法直觀比較Copula函數擬合效果：包括Genest和Rivest[10]提出的圖形法，即繪制 Copula函數的理論估計與經驗估計的Q～Q圖；Salvadori[11]提出的Pickland's依賴函數法，對比經驗與參數Copula的Pickland's依賴函數曲線，該方法也被用于多參數Copula函數的參數估計。

(3)除了從整體上評價Copula函數的擬合效果外，還應考慮選取的Copula函數能否反映極值事件的相關性。上尾部相關系數是衡量極大值事件發生時變量間的相關程度，CFG[12]為其常用的估計算法：

(8)

式中：(Ui,Vi)是根據實測點據獲得的第i組隨機樣本。

1.2 設計洪水地區組成計算方法

當下游斷面發生設計洪水zp時，上游斷面洪量X與下游區間洪量Y有無窮多種組合情況，下面給出3種具有代表性[5,7]的地區組成計算方法。

1.2.1 條件期望組合

將上游斷面洪量X取xp，將下游區間洪量Y取期望值E(y|xp)，所形成的組合(xp,E(y|xp))記為條件期望組合[5]。E(y|xp)通過下式進行計算：

(9)

f(Y≤y|X=xp)=c(u,v)fY(y)

1.2.2 條件最可能組合

當上游斷面洪量X取設計洪量xp時，下游區間洪量Y取使f(Y≤y|X=xp)達到最大值時對應的yM值，將(xp,yM)記為條件最可能組合。f(Y≤y|X=xp)對y求導得：

(10)

c2=?c/?v

令df(Y≤y|X=xp)/dy=0，整理得：

(11)

一是評先評優有傾斜。對回鄉創業貢獻突出的能人，優先提名為黨代表、人大代表、政協委員候選人，優先安排參加勞動模范、優秀企業家等各類榮譽評選。二是授予特殊榮譽。對不擔任鎮村職務的回鄉能人，可以探索設立“新鄉賢”“榮譽村民”“先鋒黨員”等稱號，充分尊重其社會表現，給于其榮譽地位，引導社會尊重他們。三是給予適當政治待遇。對表現較好的未擔任鎮村職務的回鄉能人，可以安排其作為特別會議代表，列席一定范圍的鎮村黨務政務會議，可以在重要節慶活動時由組織進行走訪慰問。

當y的邊緣分布函數為P-Ⅲ型分布時，fY(y)的表達式如下：

(12)

式中：α、β、γ分別為位置參數、尺度參數、形狀參數。

將式(12)代入式(11)中，可推導邊緣分布為P-Ⅲ型分布的二維條件最可能組合法計算通式：

(13)

通過牛頓迭代法對非線性方程式(13)求解，即可得到條件最可能組合(xp,yM)。

上游斷面洪量X取值為xp時，分別令F(Y≤y|X=xp)等于α/2和1-α/2，求得對應的下游區間洪量Y的值分別為yL和yU，稱[yL,yU]為xp對應的Y的置信區間。洪量置信區間可以定量評價設計洪水地區組成估計值的不確定性，為決策人員明晰決策風險提供了依據[13]。條件最可能組合與條件期望組合并稱為條件地區組成法。

1.2.3 最可能組合

當下游斷面發生設計洪水zp時，上游斷面洪量X與下游區間洪量Y存在一個最可能組合，即求解在X+Y=zp等式約束條件下X與Y聯合分布概率密度函數f(x,y)的最大值。利用牛頓迭代法對非線性方程式(14)求解，即可得到最可能組合(x,zp-x)。

(14)

c1=?c/?u，c2=?c/?v

式中：αx、βx、γx，αy、βy、γy分別為FX(x)、FT(y)對應的參數。

2 應用實例

選擇淮河干流中渡斷面以上流域作為研究區域(見圖1)。中渡站位于洪澤湖出口處，是淮河中游總控制站，控制流域面積15.82 萬km2?；春又杏伟春拥捞匦院偷匦纹鸱譃轸斉_子以上和以下2個區間段，魯臺子以上流域面積8.86 萬km2；魯臺子至中渡段河長335 km，區間控制面積6.96 萬km2。本文以魯臺子斷面為界，研究其上和其下2個區間段的洪水對中渡斷面的設計洪水組成的貢獻率。

圖1 中渡斷面控制水系示意圖Fig.1 Sketch map of drainage of Zhongdu site

選用1951-2010年魯臺子站和中渡站洪水的同步觀測資料及1931年實測大洪水資料，將流域內的大型水庫、湖洼、行蓄洪區的攔蓄水量、灌溉引用水量等演算至魯臺子、中渡斷面，并與斷面實測流量過程相加得到該斷面天然流量過程，根據淮河流域的洪水特點及水利工程的調節特性，選擇設計洪量控制時段為30 d[1]。繪制歷年魯臺子、中渡站天然年最大30 d洪量W30時歷曲線及模比系數累積平均值曲線[14](見圖2)，系列中較好的包含了豐、平、枯的完整過程且在45 a以上具有較強穩定性。魯臺子和中渡(以下簡稱魯中)區間洪水資料采用中渡控制斷面的天然洪水過程扣除魯臺子控制斷面錯開傳播時間演算的天然洪水過程，求得區間洪水過程后再統計各時段洪量。

圖2 魯臺子、中渡斷面W30時歷曲線及模比系數累積平均值曲線Fig.2 Time history and cumulative average modulus coefficient curve of W30 at Lutaizi and Zhongdu site

2.1 邊緣分布的確定

各區W30服從P-Ⅲ型分布，將1931年、1954年作為大洪水年，其洪量按1915年以來排位處理，采用統一處理法計算不連續樣本頻率，運用計算機優化適線法估計E(X)、Cv、Cs，并將其轉化為P-Ⅲ型分布的位置參數α、尺度參數β、形狀參數γ，其轉換關系如式(15)。最終參數估計結果如表1所示。

(15)

表1 各分區P-Ⅲ分布參數估計結果Tab.1 Parameter estimation results for P-Ⅲ distribution of each region

2.2 聯合分布函數的確定

基于單參數、雙參數、三參數GH Copula函數構建魯臺子W30與魯中區間W30的聯合分布函數。采用極大似然法估計3種GH Copula函數的參數，利用OLS、AIC準則來優選GH Copula函數，參數估計結果及對應擬合優度指標列于表2。分別繪制3種GH Copula函數的Q～Q關系圖(見圖3)；計算經驗Pickland's依賴函數與運用3種GH Copula函數分別生成1 000個隨機模擬點估算的Pickland's依賴函數，并畫圖對比(見圖4)。此外，用CFG法估計實測點據的上尾相關系數為0.636。

表2 3種GH Copula函數參數估計結果及擬合優度指標Tab.2 The estimated parameters and goodness of fit of three GH Copula Functions

圖3 魯臺子W30與魯中區間W30聯合分布的理論估計與經驗估計關系Fig.3 The theoretical and empirical distribution for joint distribution of W30 at Lutaizi and Luzhong interval

圖4 基于魯臺子W30與魯中區間W30的 經驗與擬合Pickland's依賴函數關系Fig.4 Plots of empirical and fitted Pickland's dependence functions of W30 at Lutaizi and Luzhong interval

從表2可以看出，雙參數GH Copula函數對應的OLS、AIC值最小，上尾相關系數估計值與樣本估計值最為接近，擬合效果最優。圖3中3種Copula函數的點據均落在45°線附近；圖4中雙參數GH Copula擬合的Pickland's依賴函數的曲線與經驗值最為接近，近乎重合。綜上，選取雙參數GH Copula作為構造魯臺子W30與魯中區間W30聯合分布的連接函數。

3 結果與分析

3.1 條件地區洪水組成

條件地區洪量組成考慮的是上游斷面與區間洪量的空間相關性。其中，條件期望組成強調的是在所有可能出現的洪水地區組合方案中的平均水平，可作為洪水分配方案參考；條件最可能組成則強調實測組合更具備出現的可能性，是一種條件最優分配方案。本實例計算了分別以魯臺子W30、魯中區間W30為條件，另一分區相應W30的條件地區洪量組成，并與工程上常用的同頻率組合法進行對比分析(見圖5)。

圖5 中渡站W30 3種地區組成關系Fig.5 Comparison for three regional flood composition

圖5(a)給出了以魯臺子W30為條件，魯中區間相應的W30估計值(條件期望組成I、條件最可能組成I、同頻率組成I)及其90%的置信區間。4種設計洪水地區組成方案估計值均在90%置信區間內。以條件期望值I作為比較標準，條件最可能組成I以上游設計頻率為40%為界，當頻率高于40%時，其估計值高于條件期望值；反之，其估計值低于條件期望值。而對于同頻率組成I，當設計頻率高于10%時，其估計值隨設計標準的提高逐漸向下偏離條件期望值?；春又杏螀^5～20 a一遇的中小洪水事件發生頻繁，是防范洪災風險中更應該關注的問題。根據圖5(a)分析結果，對于該頻率段洪水，同頻率組成I方案有可能使下游斷面達不到預期的防洪標準；與之相反，條件最可能組成Ⅰ則提供了一種“對防洪不利”的偏安全的組成方案。

圖5(b)給出了以魯中區間W30為條件，魯臺子相應的W30估計值(條件期望組成Ⅱ、條件最可能組成Ⅱ、同頻率組成Ⅱ)及其90%的置信區間。4種設計洪水地區組成方案估計值均在90%置信區間內。以條件期望值Ⅱ作為比較標準，條件最可能組成Ⅱ以區間W30設計頻率為42%為界，當頻率高于42%時，其估計值高于條件期望值；反之，其估計值低于條件期望值。同頻率組成值II與條件期望值接近。對于5～20 a一遇的中小洪水，3種方案計算的結果相近。但當發生設計頻率低于3%的極值洪水事件時，同頻率地區組成計算結果處于一個較低的保證率水平，難以保證下游中渡斷面的防洪安全。

3.2 最可能地區洪量組成

最可能地區組成法是從自然規律的角度出發，按照下游防洪標準設定的一種最優分配方案。表3列出了最可能組成法分配結果，同時對比條件洪水組成法、同頻率組成法估計結果。最可能組成結果位于條件期望組成Ⅰ、Ⅱ，條件最可能組成Ⅰ、Ⅱ，同頻率組成Ⅰ、Ⅱ之間。當中渡斷面發生5～20 a一遇的設計洪水時，若按最可能組成法進行分配，上游魯臺子斷面及魯中區間均發生低于該頻率段的設計洪水，2分區發生了洪水遭遇是造成中渡斷面高頻洪水的主要原因。

表3 中渡斷面設計洪水地區組成計算結果對比 億m3

3.3 設計洪水地區組成方法計算貢獻率對比

魯臺子站控制面積占中渡站控制面積的56%，而根據實測資料統計其洪水來量占中渡站洪水總量的60%以上，是淮河中上游主要的匯水區。計算綜上所有洪水分配方案在不同量級情景下上游設計洪量對中渡設計洪量的貢獻率，相應區間貢獻率也可經簡單計算得到，見圖6。從圖6可以看出在設計頻率高于20%時，條件期望組成II、條件最可能組成II、同頻率組成II計算的貢獻率存在小于56%的情況，與實際水量分配不符，不宜選擇此3種方案。同頻率組成I計算的貢獻值隨設計頻率變化的幅度很大，且在高頻段計算的貢獻值偏離其他分配方案結果，容易低估區間洪水風險，具有明顯的不合理性。設計洪水地區組成的任務是在下游斷面發生設計洪水條件下，擬定不同種洪水組成方案，使得按其確定的水庫防洪設計指標及調洪規劃在水利工程運行中能夠使下游斷面達到預期的防洪標準。條件最可能組成I隨頻率變化幅度較小，計算結果穩定，可作為確定水庫防洪設計指標的依據。條件期望組成I與最可能組成計算的貢獻率隨設計頻率的增大均存在上升趨勢，且變化規律一致。最可能組成計算的區間洪量的貢獻率最大，為中渡斷面W30設計洪水組成提供了對防洪最為不利的組成方案，可供水庫常規調度參考。

圖6 不同量級情景下魯臺子設計洪水的貢獻率Fig.6 Contribution rate of Lutaizi design flood under different magnitude scenarios

4 結 論

以淮河中游控制站中渡站為例，基于雙參數GH Copula函數，構造魯臺子斷面、魯中區間最大30 d洪量的聯合分布函數并對比分析條件期望組成、條件最可能組成、最可能組成及同頻率組成4種設計洪水地區組成分配方案的適用性，得到以下結論。

(1)雙參數GH Copula函數與單參數GH Copula函數相比結構更加靈活，能夠更好地考慮各區洪水的相依性結構。

(2)條件期望組成代表各種組成方案的平均水平，具有統計基礎和參考意義。條件最可能組成是一種條件最優分配方案；最可能組成是一種按照下游防洪標準設定的最優分配方案，可為防洪規劃提供合理的邊界參考。同頻率組成法人為地將存在相關性的上下游洪量分隔開來，再硬性組合，計算結果處于較低保證率水平。

(3)上述4種地區組成的方法都是以單項特定的地區組成情況去概括所有可能出現的多種多樣的組合，計算的洪水貢獻率存在較大差別。對于中渡斷面而言，從洪澤湖的安全及設計合理性的角度出發應選擇最可能組成分配方案作為設計參考。

□