基于年月聯(lián)合分布的典型年徑流過(guò)程設(shè)計(jì)計(jì)算

劉 進(jìn)，方國(guó)華，方功先，顏 敏

（1.河海大學(xué)水利水電學(xué)院，南京 210098；2.江蘇蘇盛工程造價(jià)咨詢有限公司，南京 210000）

典型年的選取是水能、水利規(guī)劃設(shè)計(jì)的基礎(chǔ).在確定典型年的過(guò)程中，要合理考慮年徑流總量及其年內(nèi)分配，使之具有豐水年、平水年和枯水年的代表性.通常的做法是，首先對(duì)實(shí)測(cè)徑流資料進(jìn)行頻率分析，確定某一設(shè)計(jì)頻率下典型年的總徑流量，然后選擇幾個(gè)年徑流量與其最接近的實(shí)際年份產(chǎn)生典型徑流系列，從中選出對(duì)水資源供需平衡最不利的徑流過(guò)程，按同倍比法或者同頻率法縮放后作為典型年徑流過(guò)程.采用上述這種方法設(shè)計(jì)的徑流過(guò)程，可能會(huì)因人而異，主觀經(jīng)驗(yàn)對(duì)結(jié)果的影響比較大.為了減少人為因素的影響，相關(guān)學(xué)者引入系統(tǒng)聚類法，提出采用模糊數(shù)學(xué)方法選取典型年.夏莉敏等提出一種模糊聚類迭代模型，使水文年選取過(guò)程不受主觀權(quán)重的影響，計(jì)算結(jié)果相對(duì)客觀、準(zhǔn)確［1］.胡斌奇等基于云模型思想實(shí)現(xiàn)了對(duì)徑流模糊特征的“軟劃分”，避免了人為確定分類標(biāo)準(zhǔn)的主觀性，能夠更準(zhǔn)確地描述徑流的水文豐枯特性［2］.這些研究和理論雖然統(tǒng)籌考慮了年徑流總量及其年內(nèi)分配，但未涉及年、月徑流之間的豐枯遭遇概率［3］，使得極端月徑流值會(huì)對(duì)典型年徑流過(guò)程產(chǎn)生較大的影響，且對(duì)于生成的徑流過(guò)程，也不能確定其是否為設(shè)計(jì)頻率下最有可能的年、月徑流組合.

Copula函數(shù)是研究變量相依性的工具，可以把年徑流分布與各月徑流分布有效地結(jié)合起來(lái)，能夠刻畫(huà)出年、月徑流之間的依賴關(guān)系，可作為計(jì)算年、月徑流聯(lián)合分布，考慮豐枯遭遇概率的一種方法.近年來(lái)，Copula函數(shù)由于其適用性強(qiáng)、構(gòu)造靈活、計(jì)算簡(jiǎn)單、對(duì)邊緣分布形式?jīng)]有限制等優(yōu)點(diǎn)，在降雨計(jì)算、洪水設(shè)計(jì)計(jì)算、徑流豐枯遭遇分析等水科學(xué)領(lǐng)域中受到了廣泛的應(yīng)用［4-5］.Zhang等根據(jù)美國(guó)路易斯安那州阿米特河流域的降雨數(shù)據(jù)，采用阿基米德Copula函數(shù)構(gòu)建了雙變量降雨頻率分布模型，并對(duì)降雨的聯(lián)合重現(xiàn)期進(jìn)行了計(jì)算，得出了較為理想的結(jié)果［6］.莫淑紅等建立了渭河干、支流的徑流聯(lián)合分布，證實(shí)了用Copula函數(shù)描述干、支流的豐枯遭遇概率問(wèn)題具有一定的可行性［7］.何慶升等根據(jù)長(zhǎng)江支流羅田站月徑流量資料，研究了高維阿基米德Copula函數(shù)在各月徑流量聯(lián)合分布中的應(yīng)用，并利用Copula函數(shù)對(duì)徑流預(yù)測(cè)值進(jìn)行了有效的評(píng)價(jià)［8］.晉恬等基于阿基米德Copula函數(shù)構(gòu)建了洪峰流量和洪量的二維聯(lián)合分布，并以此為基礎(chǔ)比較了不同重現(xiàn)期標(biāo)準(zhǔn)對(duì)于洪水組合設(shè)計(jì)值的影響［9］.這些研究表明，Copula函數(shù)在多變量水文事件中具有很高的應(yīng)用價(jià)值，基于Copula函數(shù)構(gòu)建的聯(lián)合分布具有較高的可行性.

因此，本文采用阿基米德Copula函數(shù)構(gòu)建年、月徑流的二維聯(lián)合分布（全年與1月徑流組合、全年與2月徑流組合等），引入條件概率，考慮年、月徑流的最有可能組合，研究并驗(yàn)證聯(lián)合分布在典型年徑流過(guò)程設(shè)計(jì)計(jì)算中的可行性.

1 徑流邊緣分布

P-Ⅲ型曲線是我國(guó)洪水設(shè)計(jì)規(guī)范中推薦采用的徑流分布曲線，少數(shù)情況下，頻率曲線可以采用其他線型，但要分析論證其可行性［10］.在國(guó)外，極值Ⅰ型（Gumbel）曲線和對(duì)數(shù)正態(tài)分布曲線在水文頻率計(jì)算中也受到廣泛應(yīng)用.選取合適的邊緣分布對(duì)聯(lián)合分布的生成至關(guān)重要.本文選擇P-Ⅲ型曲線、Gumbel型曲線和對(duì)數(shù)正態(tài)分布曲線這3種常用的曲線，對(duì)年、月徑流的邊緣分布進(jìn)行擬合.采用矩法估計(jì)分布參數(shù)，根據(jù)K-S檢驗(yàn)統(tǒng)計(jì)量（Dn）、均方根誤差（RMSE）和AIC值等統(tǒng)計(jì)指標(biāo)，選取擬合效果最好的邊緣分布線型.

2 年、月徑流聯(lián)合分布

徑流的年際、年內(nèi)分布具有不確定性，豐水年可能遭遇枯水月，枯水年可能遭遇豐水月［11-13］.在設(shè)計(jì)典型年徑流過(guò)程時(shí)，不僅要考慮到年徑流總量，還應(yīng)關(guān)注徑流的年、月相依性關(guān)系，考慮徑流的年內(nèi)分配.Copula函數(shù)作為研究變量相關(guān)性測(cè)度的一種有效的工具，可以將邊緣分布和聯(lián)合分布很好地聯(lián)系起來(lái)［14］.本文基于Copula函數(shù)原理，選擇年徑流和月徑流兩個(gè)變量，構(gòu)建年、月徑流的聯(lián)合分布.

2.1 常用的Copula函數(shù)

總的來(lái)說(shuō)，Copula函數(shù)分為t-Copula函數(shù)、正態(tài)Copula函數(shù)和阿基米德Copula函數(shù)（Archimedean Copulas）.其中，阿基米德Copula函數(shù)由于其形式多樣、參數(shù)少、易于構(gòu)造等優(yōu)點(diǎn)受到廣泛的應(yīng)用，其數(shù)學(xué)表達(dá)式取決于它的構(gòu)造函數(shù)，對(duì)于二維的阿基米德Copula函數(shù)，其定義如式（1）所示：

式中：C(u,v)為生成的Copula函數(shù)；u、v表示變量的邊緣分布；φ為構(gòu)造函數(shù)，又稱為生成元算子，對(duì)于?t∈[ 0,1)，滿足φ′(t)＜0，φ″(t)＞0且有φ(1)=0；φ-1為φ的反函數(shù).可以看出，不同的生成元產(chǎn)生不同的阿基米德Copula函數(shù)，生成元的形式?jīng)Q定了Copula函數(shù)的最終形式.二維阿基米德Copula函數(shù)中常見(jiàn)的主要有Gumbel Copula函數(shù)、Clayton Copula函數(shù)和Frank Copula函數(shù)［15-16］.

1）Gumbel Copula函數(shù)

當(dāng)生成元算子φ=( -lnt)θ，θ∈[1 ,∞)時(shí)，對(duì)應(yīng)的Copula函數(shù)稱為Gumbel Copula函數(shù)，其數(shù)學(xué)表達(dá)式為：

2）Clayton Copula函數(shù)

當(dāng)生成元算子φ=t-θ-1，θ∈[ )

0,∞時(shí)，對(duì)應(yīng)的Copula函數(shù)稱為Clayton Copula函數(shù)，其數(shù)學(xué)表達(dá)式為：

3）Frank Copula函數(shù)

式（2）、（3）、（4）中，θ為Copula函數(shù)的參數(shù)；t為生成元算子的自變量，通常取聯(lián)合分布各組成對(duì)象的邊緣分布；e表示自然對(duì)數(shù)的底.

2.2 Copula函數(shù)的參數(shù)估計(jì)

Kendall秩相關(guān)系數(shù)作為描述變量相關(guān)性的一個(gè)常見(jiàn)指標(biāo)，可以用來(lái)估計(jì)Copula函數(shù)的參數(shù).當(dāng)( )

X,Y連續(xù)時(shí)，Kendall秩相關(guān)系數(shù)有如下定義：

設(shè)(X,Y)的Copula函數(shù)為C(u,v)，則有：

聯(lián)立式（5）、式（6），得：

式（7）建立起了Copula函數(shù)與樣本的Kendall秩相關(guān)系數(shù)之間的數(shù)量關(guān)系，通過(guò)這種對(duì)應(yīng)關(guān)系來(lái)估算Copula函數(shù)的參數(shù)，是阿基米德Copula函數(shù)參數(shù)估計(jì)的常用方法［17-19］.對(duì)于常見(jiàn)的3種二維阿基米德Copula函數(shù)，其參數(shù)與Kendall秩相關(guān)系數(shù)的對(duì)應(yīng)關(guān)系如表1所示.

表1 常見(jiàn)的二維阿基米德Copula函數(shù)的參數(shù)估計(jì)方法Tab.1 Parameter estimation methods of common two-dimensional Archimedean Copulas

3 條件最有可能徑流組合

考慮年、月徑流的豐枯遭遇概率，某一特定年徑流x下，月徑流量yi服從于條件分布FY|X(yi|x).假設(shè)條件分布的概率密度函數(shù)為fY|X(yi,x)，則fY|X(yi,x)的值越大，表示月徑流yi發(fā)生的可能性越大，當(dāng)fY|X(yi,x)的值取到最大值時(shí)，認(rèn)為該徑流yi發(fā)生的概率最大，對(duì)應(yīng)的(x,yi)為一組最有可能的年、月徑流組合［20］.

條件分布概率密度函數(shù)如式（8）所示.

4 實(shí)例計(jì)算

黃河是我國(guó)境內(nèi)的第二長(zhǎng)河，發(fā)源于青藏高原，全長(zhǎng)5464 km，流域面積約752 443 km2，水面落差4480 m，蘊(yùn)含著豐富的水能資源.唐乃亥水文站位于黃河上游，站前河長(zhǎng)共1553 km，擁有200.4億m3的多年平均徑流量，占流域天然徑流總量的34.2%，是黃河上游的重要水文站.該水文站為黃河的防洪抗旱、水量調(diào)度以及水能資源的開(kāi)發(fā)利用提供了基本依據(jù)，為下游龍羊峽水調(diào)等重點(diǎn)水利工程單位提供了重要的水情信息.本文選取唐乃亥水文站1956—1997年（共42年）的徑流資料，對(duì)黃河上游的徑流進(jìn)行頻率分析，以探索年、月徑流聯(lián)合分布在典型年徑流過(guò)程設(shè)計(jì)計(jì)算中的作用.

4.1 邊緣分布

首先，建立年徑流和各月徑流的邊緣分布.選取水文頻率計(jì)算中常用的P-Ⅲ型曲線、Gumbel型曲線和對(duì)數(shù)正態(tài)分布曲線擬合樣本數(shù)據(jù)，采用矩法估計(jì)分布參數(shù).運(yùn)用K-S檢驗(yàn)方法對(duì)邊緣分布的擬合優(yōu)度進(jìn)行檢驗(yàn)，取顯著水平α=0.05，則檢驗(yàn)統(tǒng)計(jì)量的臨界值Dα=0.210，若檢驗(yàn)統(tǒng)計(jì)量Dn＜Dα，則認(rèn)為可以接受該分布曲線；若有多個(gè)分布線型通過(guò)了K-S檢驗(yàn)，則按照均方根誤差（RMSE）和AIC值最小原則選擇最優(yōu)分布函數(shù).計(jì)算結(jié)果見(jiàn)表2.

表2 徑流邊緣分布擬合優(yōu)度檢驗(yàn)Tab.2 Goodness of fit test of runoff marginal distribution

由表2可知，所有曲線分布的檢驗(yàn)統(tǒng)計(jì)量均小于臨界值0.210，都通過(guò)了K-S檢驗(yàn)，可按照RMSE和AIC值確定最終的邊緣分布.根據(jù)表2結(jié)果可知，1956—1997年這42年中，1月、8月、9月、11月的徑流量服從對(duì)數(shù)正態(tài)分布，4月、5月和10月的徑流量服從Gumbel分布，其他月份以及全年的徑流量服從P-Ⅲ型分布.

4.2 聯(lián)合分布

選擇3種常見(jiàn)的阿基米德Copula函數(shù)模擬年、月徑流的聯(lián)合分布，根據(jù)樣本的秩相關(guān)系數(shù)估算生成元參數(shù)θ.對(duì)生成的聯(lián)合分布，采用K-S檢驗(yàn)方法進(jìn)行檢驗(yàn)（置信水平為95%），并根據(jù)RMSE和AIC值選擇最優(yōu)的Copula函數(shù)，計(jì)算結(jié)果見(jiàn)表3.

表3 徑流聯(lián)合分布擬合優(yōu)度檢驗(yàn)Tab.3 Goodness of fit test of runoff joint distribution

結(jié)果表明，在95%置信水平下，初選的3種阿基米德Copula函數(shù)均通過(guò)了K-S檢驗(yàn).在全年與6月徑流聯(lián)合分布擬合結(jié)果中，Gumbel Copula函數(shù)具有最小的RMSE和AIC值，則選擇Gumbel Copula為最優(yōu)Copula函數(shù)；類似的，對(duì)于其他組合的徑流聯(lián)合分布，F(xiàn)rank Copula的RMSE和AIC值最小，則選擇Frank Copula為最優(yōu)Copula函數(shù).圖1以全年與1月的徑流組合、全年與6月的徑流組合為代表給出了經(jīng)驗(yàn)聯(lián)合概率和Copula理論概率的對(duì)比結(jié)果.從圖1可以看出，數(shù)據(jù)點(diǎn)均勻、緊密地分布在45°線的左右兩側(cè)，表明選定的Copula函數(shù)具有比較好的擬合效果.

圖1 經(jīng)驗(yàn)聯(lián)合概率與二維Copula理論概率對(duì)比圖Fig.1 Comparison of empirical joint probability and two-dimensional Copulas theoretical probability

4.3 典型年徑流過(guò)程計(jì)算

根據(jù)年、月徑流的邊緣分布和阿基米德Copula函數(shù)，計(jì)算得到徑流的聯(lián)合分布以及條件分布.取豐水年、平水年、枯水年的設(shè)計(jì)頻率（超越概率）分別為25%、50%、75%，按照條件分布概率密度函數(shù)值最大原則，分別計(jì)算年徑流在25%、50%和75%設(shè)計(jì)頻率下，各月徑流的最有可能值，從而模擬徑流的年內(nèi)分配，結(jié)果見(jiàn)表4.

表4 不同設(shè)計(jì)頻率下年、月徑流最有可能組合Tab.4 The most likely combination of yearly and monthly runoff under different design frequencies

取上述過(guò)程模擬的結(jié)果作為一組初設(shè)的典型年徑流過(guò)程.其中，豐水年徑流總量Q25%,d=225.70億m3，而根據(jù)年徑流量排頻計(jì)算的結(jié)果Q25%,p=234.62億m3，按照同倍比法對(duì)該典型年徑流過(guò)程進(jìn)行縮放，倍比系數(shù)K25%=Q25%,p/Q25%,d=1.040，即得豐水年徑流過(guò)程；同樣的，可以獲得平水年、枯水年的徑流過(guò)程，結(jié)果見(jiàn)表5.

表5 黃河上游典型年徑流過(guò)程計(jì)算結(jié)果（考慮年、月徑流聯(lián)合分布）Tab.5 Calculation results of runoff process in the upper reaches of the Yellow River in deferent typical years（considering the combined distribution of yearly and monthly runoff）

4.4 考慮聯(lián)合分布與未考慮聯(lián)合分布的徑流過(guò)程比較

關(guān)于典型年徑流過(guò)程設(shè)計(jì)計(jì)算，傳統(tǒng)的方法是，首先對(duì)長(zhǎng)期實(shí)測(cè)徑流資料進(jìn)行頻率分析，確定設(shè)計(jì)年徑流量，然后選擇年徑流量與其最接近的實(shí)際年份為典型年，對(duì)典型年徑流過(guò)程按同倍比法或者同頻率法縮放得到最終的徑流過(guò)程.圖2 a給出了按照傳統(tǒng)方法求得的豐水年、平水年、枯水年的徑流過(guò)程.其中，4月—6月的平水年徑流量比豐水年高；10月—11月的枯水年徑流量比平水年高，12月的枯水年徑流量比豐水年高.這是由于傳統(tǒng)方法在選定設(shè)計(jì)典型年徑流過(guò)程時(shí)，并未完全考慮到徑流的年內(nèi)分布，在實(shí)測(cè)歷史徑流資料中選取代表年時(shí)，平水年遭遇了豐水月，豐水年遭遇了枯水月，這樣求得的徑流過(guò)程，容易受極端月徑流值的影響，使結(jié)果不具有代表性.

圖2 黃河上游典型年的徑流過(guò)程Fig.2 Runoff process in the upper reaches of the Yellow River in deferent typical years

引入年、月徑流聯(lián)合分布后，基于條件分布考慮年、月徑流的豐枯遭遇概率，以條件最有可能徑流值代替實(shí)測(cè)值，能夠有效地避免極端月徑流值的影響，使徑流過(guò)程具有更典型的年內(nèi)分布特征（如圖2 b）.以枯水年為例對(duì)比兩種徑流過(guò)程，如表6所示，其中條件概率密度表示月徑流xi相對(duì)于年徑流（163.21億m3）的條件概率密度值.由表6可以看出，考慮聯(lián)合分布的各月徑流條件概率密度比不考慮聯(lián)合分布的各月徑流條件概率密度大得多，從概率理論上來(lái)說(shuō)，這代表著前者發(fā)生的可能性更大，是更容易出現(xiàn)的年、月徑流組合.

表6 考慮聯(lián)合分布與不考慮聯(lián)合分布所得徑流過(guò)程的對(duì)比（以枯水年為例）Tab.6 Comparison of the runoff process with and without considering the joint distribution（taking dry years for example）

另一方面，對(duì)兩種方法生成的徑流過(guò)程進(jìn)行頻率計(jì)算，得到豐水年、平水年、枯水年的年、月徑流頻率對(duì)比情況，結(jié)果見(jiàn)圖3.不考慮年月徑流聯(lián)合分布時(shí)，各月徑流頻率與年徑流設(shè)計(jì)頻率相差較大，容易出現(xiàn)豐水年遭遇枯水月（圖3 a中1月—6月）、枯水年遭遇豐水月（圖3 c中10月—12月）的極端情況；考慮聯(lián)合分布之后，豐水年1月—4月設(shè)計(jì)頻率高于25%，5月—12月設(shè)計(jì)頻率基本在25%上下（圖3 a）；平水年5月—10月設(shè)計(jì)頻率基本在50%上下，其他月份的徑流頻率低于50%（圖3 b）；枯水年6月—10月設(shè)計(jì)頻率基本在75%上下，其他月份的徑流頻率低于75%（圖3 c）.總的來(lái)說(shuō)，考慮年月徑流聯(lián)合分布后，各典型年下的年、月徑流豐枯情況存在一定的“交叉”現(xiàn)象，但避免了豐水年遭遇枯水月、枯水年遭遇豐水月的極端情況；從年內(nèi)過(guò)程看，各典型年下汛期的年、月徑流具有較高的同頻性和一致性，這是流域的季風(fēng)性氣候所致，降雨和徑流主要發(fā)生在汛期，枯期則無(wú)明顯規(guī)律.

圖3 兩種徑流過(guò)程的年月同頻情況對(duì)比Fig.3 Comparison of co-frequency between the two runoff processes

5 結(jié)論

本文針對(duì)水文典型年徑流過(guò)程設(shè)計(jì)計(jì)算中的年、月徑流豐枯遭遇概率問(wèn)題，基于二維阿基米德Copula函數(shù)構(gòu)建徑流的聯(lián)合分布，用條件分布概率密度函數(shù)確定了年、月徑流的最有可能組合，并得到了典型年徑流過(guò)程.同時(shí)，以黃河唐乃亥水文站的實(shí)測(cè)徑流為例，比較了基于年、月徑流聯(lián)合分布的徑流過(guò)程與按傳統(tǒng)方法求得的徑流過(guò)程，得出結(jié)論如下：

1）考慮聯(lián)合分布的徑流過(guò)程有效地減小了極端月徑流值對(duì)典型年徑流過(guò)程的影響，使徑流的年內(nèi)分布具有更加典型的豐枯特征，且三種代表年的徑流過(guò)程一致性較強(qiáng).

2）從年、月豐枯遭遇概率看，考慮了聯(lián)合分布后，汛期徑流的年、月同頻性有了較大的改善，年內(nèi)分布過(guò)程更具有代表性.

3）總體而言，本文基于二維阿基米德Copula函數(shù)構(gòu)建的年、月徑流聯(lián)合分布在水文典型年徑流過(guò)程的設(shè)計(jì)計(jì)算中能夠產(chǎn)生較好的結(jié)果，該方法可以應(yīng)用于其他地區(qū)的水文頻率計(jì)算，為其他地區(qū)的典型年徑流過(guò)程設(shè)計(jì)計(jì)算提供借鑒和參考.