3維干旱變量聯合概率分布研究

晁智龍

(陜西省水文水資源勘測局，陜西 西安710068)

干旱一般是指長期少降水或無降水的現象，一般分為氣象干旱、農業干旱、水文干旱和社會經濟干旱 。由于干旱的特性是隨機的，因此，概率方法成為描述干旱特性的重要方法之一。繼Yevjevich(1967)提出游程理論以來，單變量干旱特性(干旱歷時、干旱烈度、干旱強度和到達間隔)概率分布研究取得了顯著的成果。與其他水文變量一樣，單變量干旱特性一般被認為是獨立的、同分布的。然而，干旱特征一般表現為多個相依變量組成的復雜水文事件，單變量只是描述干旱的一種特征屬性，無法全面反映干旱事件的真實特征。20世紀90年代以來，多變量的聯合分析成為水文分析的研究熱點。但是，多變量聯合分布要求各變量的邊際分布屬于同一類型分布，而實際的水文變量可能服從不同類型的分布;隨著分析變量數的增加，需要較多的數據，數學模型復雜，求解困難。一些學者采用兩變量正態分布、指數分布、Gamma分布和極值分布描述多變量的聯合分布［1-13］。近年來，Copula函數被引入多變量聯合分布分析，克服上述多變量聯合分布構建的缺點，被廣泛地應用于洪水、暴雨等極值水文事件的聯合概率分布分析，成為水文計算領域的一個研究熱點，研究表明，Copula函數能有效地描述水文事件的內在規律和反映特征屬性之間的相互關系［3］。目前，采用3維變量copula函數研究干旱特性較少，多變量的頻率分析大多采用兩變量Copula函數進行推求聯合分布。本文在吸收國外Copula函數研究成果的基礎上，以渭河流域北洛河狀頭站的徑流序列為例，采用游程理論選擇干旱歷時、干旱烈度和烈度峰值為水文干旱特性變量(以下簡稱干旱變量)，應用3維Archimedean copula研究單變量分布的擬合度檢驗、相依變量的度量、Copula函數的擬合度檢驗以及干旱變量的聯合概率分布，通過數學推導，給出了3維條件概率和重現期的計算公式，闡述3維ArchimedeanCopula函數研究干旱特性聯合分布的方法。文中方法同樣適用于其他水文變量的聯合分布計算。

1 干旱定義

根據游程理論，一個干旱事件是指在一定時間段內，徑流小于或等于截取水平值。本文選擇干旱歷時(D)、干旱烈度(S)和烈度峰值(P)來描述干旱事件。截取水平值(Threshold level)選取年徑流或月徑流序列(t=1，2，…，12)的均值。干旱歷時即負游程長度，指徑流連續低于截取水平值的年或月份數。干旱烈度為干旱事件在干旱歷時內低于截取水平值的虧缺水量和，也稱負游程總量。烈度峰值表示干旱事件在干旱歷時內的最大虧缺水量。

2 3維Archimedean copulas

Copula的定義和性質可參見有關文獻，這里不再敘述。對于1個參數的Archimedean Copula函數，Copula可以表示為

式中，φ(t)稱為Archimedean Copula的生成函數(generator function)，它完全決定了Archimedean Copula的形式。

3 應用實例

本文選用渭河流域北洛河狀頭站(1937-2006年)的月徑流序列，截取水平取各月的平均流量值。

3.1 邊際分布和參數計算

干旱變量的經驗分布采用Gringorten position-plotting formula公式計算。干旱歷時假定服從2參數Weibull分布，干旱烈度和烈度峰值服從2參數Gamma分布。采用矩法(MOM)、極大似然法(MLM)，概率權重法 (PWM)和遺傳算法(GA)進行參數估計(見表2)。根據均方根誤差(Root Mean Square Error—RMSE)，Akaike信息準則 (Akaike Ⅰnformation Criteria—AⅠC)和 Bayesian信息準則 (Bayesian Ⅰnformation Criteria—BⅠC)進行最優參數選擇。經計算，GA法對應的AⅠC、BⅠC和RMSE值均為最小，所以，3種干旱變量的邊際分布參數選擇GA法估算的參數。

干旱歷時采用χ2檢驗。干旱烈度和烈度峰值分別采用Kolmogorov-Smirnov(K—S)檢驗 Dn、Cramer-von Mises(C—M)檢驗、Aderson-Darling(A—D)檢驗、修正權重Watson檢驗和 Liao—Shimokawa檢驗 Ln［6，24］。取模擬樣本數NB=3000，5種樣本統計檢驗值均小于相應的檢驗臨界值;所以，本文接受干旱烈度和烈度峰值服從 2參數Gamma分布。

3.2 變量相依性的度量與Copula參數計算

采用Pearson’s古典相關系數 γ，Spearman秩相關系數ρ，Kendall’sτ，Chi-Plots，and K-Plots進行干旱變量相依性的度量［25］。γ、ρ 和 Kendall’sτ均大于 0.5，表明了干旱歷時、干旱烈度和烈度峰值間有較強的相依性。

3維干旱變量的聯合經驗分布仍采用Gringorten position-plotting公式計算。對于3維Archimedean Copula，采用極大似然法進行參數估算。似然函數兩邊取對數，求θ的導數，并令其導數表達式等于0，有:

根據 Gumbel Copula、Frank Copula、Clayton Copula 和 Ali-Mikhail-Haq Copula參數的取值范圍，本文采用二分法求解式(2)所示的非線性方程組，見表1。

表1 三變量Archimedean copulas聯合分布參數和擬合評價計算結果

表1表明 Gumbel Copula的 AⅠC、BⅠC和 RMSE值最小。因此，本文選用3維Gumbel Copula描述干旱變量的聯合分布特性。

3.3 條件概率和重現期計算

干旱變量是按照截取水平獲得的，這實際上是一個部分歷時序列(Partial duration series，PDS)，對于給定的截取水平(Threshold level)，邊際分布為FX(x)的重現期TX應按下式進行計算。

對于2維聯合分布 F(X≤x，Y≤y)=C(u1，u2)和條件分布 F1(X≤x|Y=y)，F2(Y≤y|X≤x)，邊際分布 u1=FX(x)，u2=FY(y)，2維重現期T(X≤x，Y≤y)和條件重現期 T(X≤x|Y=y)、T(X≤x|Y≤y)計算為:

3維重現期T(X≤x，Y≤y，Z≤z)和條件重現期T(X≤x|Y=y，Z=z)、T(X≤x，Y≤y|Z=z)、T(X≤x|Y≤y，Z≤z)有類似計算式，本文不再列舉。

圖1 給定P=p條件，D≤d，V≤v出現的重現期

所以，上式條件重現期關鍵在于如何計算相應的條件概率。由概率論原理，我們可以推導出這些條件概率的計算公式。

圖1分別給出了在P=10、P=15條件下，(D≤d，V≤v)事件發生的重現期。

4 結語

本文以干旱變量聯合分布為例，在吸收當前copula函數研究成果的基礎上，系統地闡述了單變量分布的擬合度檢驗，相依變量度量，copula參數估算，最優 copula選擇以及copula擬合度檢驗等方法和步驟。上述這些問題均為copula函數在多變量水文頻率分析中的幾個關鍵技術問題。最后，通過數學推導，給出了多變量條件概率計算以及相應的重現期計算公式，更正了Zhang L，Singh VP(2007)給出給定兩變量值(X=x，Y=y)條件下，事件Z≤z發生的條件概率計算公式。采用截取水平獲得的水文序列實際上是一個部分歷時序列(PDS)，其重現期應采用Kim(2003)給出的公式進行計算。文中實例應用表明:copula函數能夠推求變量不同組合條件下的概率分布，是推求干旱多變量聯合分布和其他水文多變量聯合分布的有效途徑，文中方法同樣適用其它水文變量的聯合分布計算。

［1］袁超，宋松柏，荊萍.極限水文干旱歷時概率分布解析法研究［J］.西北農林科技大學學報(自然科學版).2008，36(7):212-218.

［2］郭生練，閆寶偉，肖義，等.Copula函數在多變量水文分析計算中的應用及研究進展［J］.水文.2008，38(3):1-7.

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

［4］馮平，毛慧慧，王勇.多變量情況下的水文頻率分析方法及其應用［J］.水利學報.2009，40(1):33-37.

［5］許月萍，李佳，曹飛鳳，等.Copula在水文極限事件分析中的應用［J］.浙江大學學報(工學版).2008.