基于阿基米德聯(lián)合函數(shù)的洪水過程隨機(jī)模擬研究

楊立峰，周研來，夏世明，曾令軍

（1.中國水利水電第十四工程局有限公司，云南 昆明 650041；2.武漢大學(xué)水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072）

洪水過程的隨機(jī)性是復(fù)雜水文系統(tǒng)中較難描述的不確定性因素之一。洪水過程隨機(jī)模擬國內(nèi)外已有多種模型[1-3]。近年來，非參數(shù)方法[3]、小波分析理論[4]以及聯(lián)合 （Copula）函數(shù)[5-8]等被先后應(yīng)用到洪水隨機(jī)模擬領(lǐng)域，極大地豐富了洪水隨機(jī)模擬的理論和方法。文獻(xiàn) [7]采用耿貝爾聯(lián)合 （Gumbel-Hougard Copula）函數(shù)對洪峰和洪量建立聯(lián)合隨機(jī)模型，初步探討了Copula函數(shù)在洪水過程隨機(jī)模擬中的應(yīng)用。文獻(xiàn) [8]應(yīng)用Gumbel-Hougard Copula函數(shù)構(gòu)造洪峰與歷時的聯(lián)合分布，以及相鄰截口的聯(lián)合分布，較好地描述了洪峰、歷時的統(tǒng)計特征及相關(guān)關(guān)系，但洪量的統(tǒng)計特征值的模擬精度相對較差。

本文分別采用三種阿基米德聯(lián)合 （Archimedean Copula）函數(shù)對洪峰和洪量建立聯(lián)合隨機(jī)模型，從實(shí)測資料中優(yōu)選峰量比接近的洪水過程進(jìn)行縮放得到模擬的洪水過程線，來探討阿基米德聯(lián)合函數(shù)在洪水過程隨機(jī)模擬中的適用性。實(shí)例表明：Copula模擬模型彌補(bǔ)了常規(guī)模型模擬峰和量的統(tǒng)計特征不能同時與實(shí)測系列保持一致以及生成洪水的頻率分布偏離洪水設(shè)計值的理論頻率分布較大的不足，其模擬系列能夠很好地保持實(shí)測洪水過程線的形狀，截口具有隨機(jī)波動性、無鋸齒現(xiàn)象。

1 洪水過程的隨機(jī)模擬模型

1.1 Copula理論

Sklar在1959年提出了著名的多元分布的Sklar定理[10]：

令 F（·，…，·）為具有邊緣分布 F1（·），F(xiàn)2（·），…，F(xiàn)N（·）的聯(lián)合分布函數(shù)，那么存在一個Copula函數(shù)C（·，…，·）， 滿足

若 F1（·），F(xiàn)2（·），…，F(xiàn)N（·）連續(xù)， 則 C（·，…，·）唯一確定； 反之， 若 F1（·），F(xiàn)2（·），…，F(xiàn)N（·）為一元分布，C（·，…，·）為相應(yīng)的Copula函數(shù)，那么由式（2） 定義的函數(shù) F（·，…，·）為具有邊緣分布 F1（·），F(xiàn)2（·），…，F(xiàn)N（·）的聯(lián)合分布函數(shù)。 Gumbel-Hougard、克萊特 （Clayton）和弗蘭克 （Frank）聯(lián)合函數(shù)是三種常用于隨機(jī)模擬領(lǐng)域的單參數(shù)Archimedean Copula[6，10]。

N元Gumbel-Hougard Copula函數(shù)

式中， un=Fn（xn），n=1，2，…， N 為邊緣分布函數(shù)， α∈[1，∞）為Gumbel-Hougard Copula函數(shù)的參數(shù)。該函數(shù)僅能夠適用于變量存在正相關(guān)的情形，因此，較適用于構(gòu)造如洪峰和洪量、洪量與洪水歷時等變量之間存在的正相關(guān)的聯(lián)合分布。

N元Clayton Copula函數(shù)

式中，β∈（0，∞）為Clayton Copula函數(shù)的參數(shù)。該函數(shù)與Gumbel-Hougard Copula函數(shù)一樣，均僅適合描述正相關(guān)的隨機(jī)變量。

N元Frank Copula函數(shù)

式中， γ∈（-∞，+∞）{0} 為 Frank Copula函數(shù)的參數(shù)。該函數(shù)既能夠描述正相關(guān)的隨機(jī)變量，也能夠描述存在負(fù)相關(guān)的隨機(jī)變量。

通過Sklar定理，可以將一個N維的聯(lián)合分布函數(shù)分解為N個邊緣分布和一個Copula函數(shù)，這個Copula函數(shù)描述了變量間的相關(guān)性 [5，6，10]。 因此， 中參數(shù) C（u1，u2，…，uN；θ）的估計成為了一個非常重要的問題。國內(nèi)外許多學(xué)者對此問題做了廣泛的研究[5，6，9，10]， 主要有非參數(shù)估計方法和參數(shù)估計方法。非參數(shù)估計法包括Kendall秩法、Spearman秩法等。非參數(shù)估計法適用于單參數(shù)的兩變量阿基米德聯(lián)合函數(shù)估計。參數(shù)估計法主要有極大似然估計法 （ML）、 兩步法 （IFM）、 半?yún)?shù)法 （CML）， 其中ML與IFM作參數(shù)估計時，必須對邊緣分布作出假設(shè)和檢驗(yàn)，才能估計出Copula函數(shù)中的參數(shù)，而CML不需對邊緣分布作出假設(shè)[8]。

1.2 Copula聯(lián)結(jié)模型

本次研究采用阿基米德聯(lián)合函數(shù)來描述洪峰和洪量的相關(guān)性結(jié)構(gòu)，在此僅以Clayton Copula函數(shù)為例來表述洪峰和洪量的相關(guān)性。即

式中，β∈（0，∞）為Clayton Copula函數(shù)的參數(shù)，β越大，相關(guān)性越強(qiáng)；當(dāng)β趨近于0時，變量相互獨(dú)立。

根據(jù)洪峰Qm和時段洪量W的邊緣分布FQm（qm）和 FW（w）以及 Copula函數(shù) （式 （5））， 由 Sklar定理（式 （1））， 可以得到聯(lián)合分布 F（qm，w）

式中， FQm（qm）， FW（w）均為 P-III型分布。

聯(lián)合分布的參數(shù)估計先用Kendall秩法初估計，然后用兩步法 （IFM）校正。單參數(shù)的兩變量Clayton Copula函數(shù)的參數(shù)β與Kendall秩τ存在確定的解析關(guān)系[5，6，9，10]

兩步法校正參數(shù)：第1步，根據(jù)Qm和W的觀測值系列，采用單變量的參數(shù)估計邊緣分布FQm（qm）和FW（w）的參數(shù)；第2步，估計Copula函數(shù)的參數(shù)β。最優(yōu)的估計值β應(yīng)使得

式中， Femp（qm，i，wi）為根據(jù)聯(lián)合觀測值 （qm，i，wi）求出的經(jīng)驗(yàn)聯(lián)合分布[5，6]， F＾（qm，i，wi）為所采用的理論分布的估計值。

一般按照峰高量大，主峰靠后的年最大值法選出的洪水樣本，其洪峰和洪量都具有一定的相關(guān)性[7]。采用峰量聯(lián)合分布對洪峰和洪量進(jìn)行聯(lián)合描述，通過聯(lián)合分布的隨機(jī)抽樣方法可以對存在相關(guān)性的峰和量成對取樣，從而達(dá)到隨機(jī)模擬洪峰和洪量的目的。峰量隨機(jī)抽樣方法以及模擬洪水過程線的方法，具體步驟詳見文獻(xiàn)[7]。

2 實(shí)例研究

黃河中上游某水庫的流域控制面積為8 847 km2，有20年實(shí)測壩址洪水資料 （1959年至1978年，每年7月1日到10月31日，次洪歷時大多為24 h）。在實(shí)測洪水資料中，按照 “峰高量大、主峰靠后”的原則每年選擇一場歷時為34 h的最大洪水過程，構(gòu)成樣本容量為20的34 h（截口間距為20 min）洪水過程系列。從每一年洪水過程中， 取樣得到一對洪峰Qm和1日洪量W1d，構(gòu)成20個Qm和W1d的聯(lián)合觀測值系列，并據(jù)此建立Qm和W1d的兩變量聯(lián)合分布。Qm和W1d的邊緣分布采用P-III型分布，分別為 FQm（qm）和 FW1d（w1d）。 洪峰和洪量邊緣分布函數(shù)的擬合分別見圖1和圖2。

圖1 洪峰邊緣分布函數(shù)的擬合曲線

圖2 洪量邊緣分布函數(shù)的擬合曲線

采用矩法估計洪峰和洪量實(shí)測系列的統(tǒng)計特征值 （見表1），并假定估計的參數(shù)為該水庫區(qū)間洪水洪峰和洪量總體分布的參數(shù)。對Qm和W1d的聯(lián)合觀測值計算Kendall秩τ等于0.7516，據(jù)Kendall秩與Gumbel-Hougard、Clayton和Frank Copula函數(shù)的解析關(guān)系式[10]初估Copula函數(shù)的參數(shù)分別為α=4.02、β=6.04及γ=14.24。為進(jìn)一步校核Copula函數(shù)的參數(shù)，點(diǎn)繪20年的Qm和W1d的經(jīng)驗(yàn)聯(lián)合分布值與理論聯(lián)合分布值 （在此僅以Clayton Copula模型為例）擬合最佳的曲線如圖3所示。將調(diào)整后的參數(shù)α、β及γ作為Copula函數(shù)的估計值。據(jù)式 （2）～（4），可以得到 Qm和W1d的 Gumbel-Hougard、Clayton和Frank聯(lián)合分布函數(shù)分別為

依據(jù)峰量隨機(jī)抽樣方法以及模擬洪水過程線的方法，可隨機(jī)模擬出M （如10萬）條洪水過程線，用矩法估計模擬系列的統(tǒng)計特征值見表1。

表1 實(shí)測系列和模擬系列洪峰和洪量的統(tǒng)計特征值

圖3 基于Clayton Copula的聯(lián)合觀測點(diǎn)的經(jīng)驗(yàn)分布和理論分布的相關(guān)關(guān)系

由隨機(jī)抽樣得到洪峰和洪量值，可求得對應(yīng)的峰量比值，并依據(jù)該比值優(yōu)選峰量比值最接近的實(shí)測洪水過程，采用變倍比方法進(jìn)行縮放，圖4給出了以1961年實(shí)測洪水過程為例進(jìn)行縮放得到的若干洪水過程線 （在此僅以Gumbel-Hougard Copula模型為例），圖5給出了SAR（1）季節(jié)性階自回歸模型模擬的若干洪水過程線。

模擬生成的洪水頻率分布是否接近理論頻率分布是評價隨機(jī)模型的適用性重要指標(biāo)之一。由上述4個隨機(jī)模型模擬生成10個樣本容量為10萬的洪水過程，結(jié)合綜合判別法來分析洪水的量級 （用重現(xiàn)期表示），即在洪水調(diào)度中按庫水位和入庫流量中哪一項先滿足各自的最大值來判別洪水的級別，便可計算出10個模擬樣本的大于不同重現(xiàn)期設(shè)計值的平均頻次 （四舍五入取整）見表2。

圖4 以1961年實(shí)測洪水過程為例進(jìn)行縮放得到的若干洪水過程線

圖5 SAR （1）模型模擬的若干洪水過程線

表2 大于不同重現(xiàn)期設(shè)計值的頻次 （基數(shù)10萬次）

3 方法的比較分析

本文以SAR（1）模型為參照，從洪水統(tǒng)計特征參數(shù)、洪水過程線及洪水頻率分布三個方面來評價Copula模型的在洪水過程隨機(jī)模擬領(lǐng)域的適用性。

由表1和圖3可知，三種Copula模型雖然在洪峰和洪量的統(tǒng)計特征指標(biāo)上存在差異，但是均能保持實(shí)測洪水系列的統(tǒng)計特征，且三種模擬模型在峰和量的統(tǒng)計特征上均優(yōu)于SAR（1）模擬模型。

由圖4和圖5可知，SAR（1）模擬的洪水過程線型具有多樣性和各個截口有較好隨機(jī)波動性，但存在明顯的鋸齒現(xiàn)象；而Copula模擬的洪水過程線型不但保持了實(shí)測洪水過程線型，而且無鋸齒現(xiàn)象。

由表2可以看出，四種模擬生成的100年一遇至10 000年一遇的大洪水平均頻率分布上均高于理論頻率分布，小洪水平均頻率分布上均略低于洪水設(shè)計值的理論頻率分布，但Copula模型的洪水平均頻率分布較SAR（1）模型更接近洪水設(shè)計值的理論頻率分布。

4 結(jié)論

本文構(gòu)建了洪峰和洪量的聯(lián)合分布，分別采用了三種Archimedean聯(lián)合分布隨機(jī)抽樣方法模擬了洪水過程線，以SAR（1）模型為參照，探討了Archimedean Copula函數(shù)在洪水過程隨機(jī)模擬領(lǐng)域的適用性。雖然三種Archimedean Copula模型生成的洪水過程的洪峰和洪量統(tǒng)計特征有細(xì)微的差別，但均具有保持實(shí)測系列洪峰和洪量統(tǒng)計特征的優(yōu)點(diǎn)，且能夠反映實(shí)測洪水過程的形狀及洪水頻率分布特征。該模型也存在不足，如無法模擬洪峰與歷時的相關(guān)特性，且洪峰和洪量的時程分配形式受實(shí)測資料長度的限制。

［1］ 梅亞東，談廣鳴.大壩防洪安全的風(fēng)險分析［J］.武漢大學(xué)學(xué)報（工學(xué)版）， 2002， 35（6）:11-15.

［2］ 丁晶，鄧育仁，侯玉，等.水庫防洪安全設(shè)計時設(shè)計洪水過程線適用性的探討［J］.水科學(xué)進(jìn)展， 1992， 3（1）:46-52.

［3］ 袁鵬，王文圣，丁晶.洪水隨機(jī)模擬的非參數(shù)擾動最近鄰自展模型［J］.四川大學(xué)學(xué)報 （工程科學(xué)版）， 2000， 32（1）:82-86.

［4］ 王文圣，袁鵬，丁晶.小波分析及其在日流量過程隨機(jī)模擬中的應(yīng)用［J］.水利學(xué)報， 2000， 31（11）:43-48.

［5］ 熊立華，郭生練，肖義，等.Copula聯(lián)結(jié)函數(shù)在多變量水文頻率分析中的應(yīng)用［J］.武漢大學(xué)學(xué)報 （工學(xué)版）， 2005， 38（6）:16-19.

［6］ 郭生練，閆寶偉，方彬，等.Copula函數(shù)在多變量水文分析計算中的應(yīng)用及研究進(jìn)展［J］.水文， 2008， 28（3）:1-7.

［7］ 肖義，郭生練，熊立華，等.一種新的洪水過程隨機(jī)模擬方法研究［J］.四川大學(xué)學(xué)報 （工程科學(xué)版）， 2007， 39（2）:55-60.

［8］ 張濤，趙春偉，雒文生.基于Copula函數(shù)的洪水過程隨機(jī)模擬［J］.武漢大學(xué)學(xué)報 （工學(xué)版）， 2008， 41（4）:1-4.

［9］ 閆寶偉，郭生練，劉攀，等.基于Copula函數(shù)的徑流隨機(jī)模擬［J］.四川大學(xué)學(xué)報 （工程科學(xué)版）， 2010， 42（1）:5-9.