基于阿基米德聯合函數的洪水過程隨機模擬研究

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

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

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

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

1 洪水過程的隨機模擬模型

1.1 Copula理論

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

令 F（·，…，·）為具有邊緣分布 F1（·），F2（·），…，FN（·）的聯合分布函數，那么存在一個Copula函數C（·，…，·）， 滿足

若 F1（·），F2（·），…，FN（·）連續， 則 C（·，…，·）唯一確定； 反之， 若 F1（·），F2（·），…，FN（·）為一元分布，C（·，…，·）為相應的Copula函數，那么由式（2） 定義的函數 F（·，…，·）為具有邊緣分布 F1（·），F2（·），…，FN（·）的聯合分布函數。 Gumbel-Hougard、克萊特 （Clayton）和弗蘭克 （Frank）聯合函數是三種常用于隨機模擬領域的單參數Archimedean Copula[6，10]。

N元Gumbel-Hougard Copula函數

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

N元Clayton Copula函數

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

N元Frank Copula函數

式中， γ∈（-∞，+∞）{0} 為 Frank Copula函數的參數。該函數既能夠描述正相關的隨機變量，也能夠描述存在負相關的隨機變量。

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

1.2 Copula聯結模型

本次研究采用阿基米德聯合函數來描述洪峰和洪量的相關性結構，在此僅以Clayton Copula函數為例來表述洪峰和洪量的相關性。即

式中，β∈（0，∞）為Clayton Copula函數的參數，β越大，相關性越強；當β趨近于0時，變量相互獨立。

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

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

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

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

式中， Femp（qm，i，wi）為根據聯合觀測值 （qm，i，wi）求出的經驗聯合分布[5，6]， F＾（qm，i，wi）為所采用的理論分布的估計值。

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

2 實例研究

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

圖1 洪峰邊緣分布函數的擬合曲線

圖2 洪量邊緣分布函數的擬合曲線

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

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

表1 實測系列和模擬系列洪峰和洪量的統計特征值

圖3 基于Clayton Copula的聯合觀測點的經驗分布和理論分布的相關關系

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

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

圖4 以1961年實測洪水過程為例進行縮放得到的若干洪水過程線

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

表2 大于不同重現期設計值的頻次 （基數10萬次）

3 方法的比較分析

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

由表1和圖3可知，三種Copula模型雖然在洪峰和洪量的統計特征指標上存在差異，但是均能保持實測洪水系列的統計特征，且三種模擬模型在峰和量的統計特征上均優于SAR（1）模擬模型。

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

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

4 結論

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

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

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

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

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

［5］ 熊立華，郭生練，肖義，等.Copula聯結函數在多變量水文頻率分析中的應用［J］.武漢大學學報 （工學版）， 2005， 38（6）:16-19.

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

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

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

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