晁智龍
(陜西省水文水資源勘測局,陜西 西安710068)
干旱一般是指長期少降水或無降水的現(xiàn)象,一般分為氣象干旱、農(nóng)業(yè)干旱、水文干旱和社會經(jīng)濟干旱 。由于干旱的特性是隨機的,因此,概率方法成為描述干旱特性的重要方法之一。繼Yevjevich(1967)提出游程理論以來,單變量干旱特性(干旱歷時、干旱烈度、干旱強度和到達間隔)概率分布研究取得了顯著的成果。與其他水文變量一樣,單變量干旱特性一般被認為是獨立的、同分布的。然而,干旱特征一般表現(xiàn)為多個相依變量組成的復(fù)雜水文事件,單變量只是描述干旱的一種特征屬性,無法全面反映干旱事件的真實特征。20世紀90年代以來,多變量的聯(lián)合分析成為水文分析的研究熱點。但是,多變量聯(lián)合分布要求各變量的邊際分布屬于同一類型分布,而實際的水文變量可能服從不同類型的分布;隨著分析變量數(shù)的增加,需要較多的數(shù)據(jù),數(shù)學(xué)模型復(fù)雜,求解困難。一些學(xué)者采用兩變量正態(tài)分布、指數(shù)分布、Gamma分布和極值分布描述多變量的聯(lián)合分布[1-13]。近年來,Copula函數(shù)被引入多變量聯(lián)合分布分析,克服上述多變量聯(lián)合分布構(gòu)建的缺點,被廣泛地應(yīng)用于洪水、暴雨等極值水文事件的聯(lián)合概率分布分析,成為水文計算領(lǐng)域的一個研究熱點,研究表明,Copula函數(shù)能有效地描述水文事件的內(nèi)在規(guī)律和反映特征屬性之間的相互關(guān)系[3]。目前,采用3維變量copula函數(shù)研究干旱特性較少,多變量的頻率分析大多采用兩變量Copula函數(shù)進行推求聯(lián)合分布。本文在吸收國外Copula函數(shù)研究成果的基礎(chǔ)上,以渭河流域北洛河狀頭站的徑流序列為例,采用游程理論選擇干旱歷時、干旱烈度和烈度峰值為水文干旱特性變量(以下簡稱干旱變量),應(yīng)用3維Archimedean copula研究單變量分布的擬合度檢驗、相依變量的度量、Copula函數(shù)的擬合度檢驗以及干旱變量的聯(lián)合概率分布,通過數(shù)學(xué)推導(dǎo),給出了3維條件概率和重現(xiàn)期的計算公式,闡述3維ArchimedeanCopula函數(shù)研究干旱特性聯(lián)合分布的方法。文中方法同樣適用于其他水文變量的聯(lián)合分布計算。
根據(jù)游程理論,一個干旱事件是指在一定時間段內(nèi),徑流小于或等于截取水平值。本文選擇干旱歷時(D)、干旱烈度(S)和烈度峰值(P)來描述干旱事件。截取水平值(Threshold level)選取年徑流或月徑流序列(t=1,2,…,12)的均值。干旱歷時即負游程長度,指徑流連續(xù)低于截取水平值的年或月份數(shù)。干旱烈度為干旱事件在干旱歷時內(nèi)低于截取水平值的虧缺水量和,也稱負游程總量。烈度峰值表示干旱事件在干旱歷時內(nèi)的最大虧缺水量。
Copula的定義和性質(zhì)可參見有關(guān)文獻,這里不再敘述。對于1個參數(shù)的Archimedean Copula函數(shù),Copula可以表示為
式中,φ(t)稱為Archimedean Copula的生成函數(shù)(generator function),它完全決定了Archimedean Copula的形式。
本文選用渭河流域北洛河狀頭站(1937-2006年)的月徑流序列,截取水平取各月的平均流量值。
干旱變量的經(jīng)驗分布采用Gringorten position-plotting formula公式計算。干旱歷時假定服從2參數(shù)Weibull分布,干旱烈度和烈度峰值服從2參數(shù)Gamma分布。采用矩法(MOM)、極大似然法(MLM),概率權(quán)重法 (PWM)和遺傳算法(GA)進行參數(shù)估計(見表2)。根據(jù)均方根誤差(Root Mean Square Error—RMSE),Akaike信息準則 (Akaike Ⅰnformation Criteria—AⅠC)和 Bayesian信息準則 (Bayesian Ⅰnformation Criteria—BⅠC)進行最優(yōu)參數(shù)選擇。經(jīng)計算,GA法對應(yīng)的AⅠC、BⅠC和RMSE值均為最小,所以,3種干旱變量的邊際分布參數(shù)選擇GA法估算的參數(shù)。
干旱歷時采用χ2檢驗。干旱烈度和烈度峰值分別采用Kolmogorov-Smirnov(K—S)檢驗 Dn、Cramer-von Mises(C—M)檢驗、Aderson-Darling(A—D)檢驗、修正權(quán)重Watson檢驗和 Liao—Shimokawa檢驗 Ln[6,24]。取模擬樣本數(shù)NB=3000,5種樣本統(tǒng)計檢驗值均小于相應(yīng)的檢驗臨界值;所以,本文接受干旱烈度和烈度峰值服從 2參數(shù)Gamma分布。
采用Pearson’s古典相關(guān)系數(shù) γ,Spearman秩相關(guān)系數(shù)ρ,Kendall’sτ,Chi-Plots,and K-Plots進行干旱變量相依性的度量[25]。γ、ρ 和 Kendall’sτ均大于 0.5,表明了干旱歷時、干旱烈度和烈度峰值間有較強的相依性。
3維干旱變量的聯(lián)合經(jīng)驗分布仍采用Gringorten position-plotting公式計算。對于3維Archimedean Copula,采用極大似然法進行參數(shù)估算。似然函數(shù)兩邊取對數(shù),求θ的導(dǎo)數(shù),并令其導(dǎo)數(shù)表達式等于0,有:
根據(jù) Gumbel Copula、Frank Copula、Clayton Copula 和 Ali-Mikhail-Haq Copula參數(shù)的取值范圍,本文采用二分法求解式(2)所示的非線性方程組,見表1。
表1 三變量Archimedean copulas聯(lián)合分布參數(shù)和擬合評價計算結(jié)果
表1表明 Gumbel Copula的 AⅠC、BⅠC和 RMSE值最小。因此,本文選用3維Gumbel Copula描述干旱變量的聯(lián)合分布特性。
干旱變量是按照截取水平獲得的,這實際上是一個部分歷時序列(Partial duration series,PDS),對于給定的截取水平(Threshold level),邊際分布為FX(x)的重現(xiàn)期TX應(yīng)按下式進行計算。
對于2維聯(lián)合分布 F(X≤x,Y≤y)=C(u1,u2)和條件分布 F1(X≤x|Y=y),F(xiàn)2(Y≤y|X≤x),邊際分布 u1=FX(x),u2=FY(y),2維重現(xiàn)期T(X≤x,Y≤y)和條件重現(xiàn)期 T(X≤x|Y=y)、T(X≤x|Y≤y)計算為:
3維重現(xiàn)期T(X≤x,Y≤y,Z≤z)和條件重現(xiàn)期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出現(xiàn)的重現(xiàn)期
所以,上式條件重現(xiàn)期關(guān)鍵在于如何計算相應(yīng)的條件概率。由概率論原理,我們可以推導(dǎo)出這些條件概率的計算公式。
圖1分別給出了在P=10、P=15條件下,(D≤d,V≤v)事件發(fā)生的重現(xiàn)期。
本文以干旱變量聯(lián)合分布為例,在吸收當前copula函數(shù)研究成果的基礎(chǔ)上,系統(tǒng)地闡述了單變量分布的擬合度檢驗,相依變量度量,copula參數(shù)估算,最優(yōu) copula選擇以及copula擬合度檢驗等方法和步驟。上述這些問題均為copula函數(shù)在多變量水文頻率分析中的幾個關(guān)鍵技術(shù)問題。最后,通過數(shù)學(xué)推導(dǎo),給出了多變量條件概率計算以及相應(yīng)的重現(xiàn)期計算公式,更正了Zhang L,Singh VP(2007)給出給定兩變量值(X=x,Y=y)條件下,事件Z≤z發(fā)生的條件概率計算公式。采用截取水平獲得的水文序列實際上是一個部分歷時序列(PDS),其重現(xiàn)期應(yīng)采用Kim(2003)給出的公式進行計算。文中實例應(yīng)用表明:copula函數(shù)能夠推求變量不同組合條件下的概率分布,是推求干旱多變量聯(lián)合分布和其他水文多變量聯(lián)合分布的有效途徑,文中方法同樣適用其它水文變量的聯(lián)合分布計算。
[1]袁超,宋松柏,荊萍.極限水文干旱歷時概率分布解析法研究[J].西北農(nóng)林科技大學(xué)學(xué)報(自然科學(xué)版).2008,36(7):212-218.
[2]郭生練,閆寶偉,肖義,等.Copula函數(shù)在多變量水文分析計算中的應(yīng)用及研究進展[J].水文.2008,38(3):1-7.
[3]熊立華,郭生練,肖義,等.Copula聯(lián)結(jié)函數(shù)在多變量水文頻率分析中的應(yīng)用[J].武漢大學(xué)學(xué)報(工學(xué)版).2005,38(6):16-19.
[4]馮平,毛慧慧,王勇.多變量情況下的水文頻率分析方法及其應(yīng)用[J].水利學(xué)報.2009,40(1):33-37.
[5]許月萍,李佳,曹飛鳳,等.Copula在水文極限事件分析中的應(yīng)用[J].浙江大學(xué)學(xué)報(工學(xué)版).2008.