Copula函數在水文多變量分析計算中的問題

宋松柏

（1.西北農林科技大學 水利與建筑工程學院，陜西楊凌712100；2.西北農林科技大學旱區農業水土工程教育部重點實驗室，陜西楊凌712100）

水文事件一般具有多種相依性。選擇一種特征屬性變量，應用單變量頻率分析難以為水利工程規劃、設計和管理提供全面的水文設計值信息。傳統多變量聯合分布要求其變量邊際分布必須為同一類型的分布，而實際水文現象具有高度的非線性和復雜性，同一水文事件的多種特征屬性變量不可能服從同一類型的分布。因而，傳統多變量聯合分布的這種不足限制了其應用［1-5］。1959年以來，隨著 Copula函數理論的發展，基于Copula函數的多變量分析開始受到重視。20世紀 90年代以后，Copula函數被引入洪水［1-14］、暴雨［1-4］、干旱［15-29］、徑流預測［30］、徑流模擬［31］、枯水概率［32］、水 文 預 報 不 確 定性［33-34］、 水文 豐枯 遭 遇 特征［35-36］、水量水質聯合分布［37］、水資源供需風險損失［38］、水沙關系演變［39］、外調水供水補償特性［40］、缺水風險［41］、用水效率［42］、設計潮位［43-45］和大壩洪水漫頂風險率［46］等的聯合概率分析計算。從研究文獻來看，目前在水文學中應用的Copula函數主要有：對稱Archimedean Copula以及非對稱的Archimedean Copula、Plackette Copula、Metaelliptical Copula 和混合Copula 等［2，5］。 Copula 函數參數估計方法有精確極大似然法、邊際函數推斷法和半參數法。其中：邊際函數推斷法也稱兩階段法（two stage method），是目前水文學中廣泛采用的Copula函數參數估計方法。按照邊際函數推斷法，基于Copula函數的水文多變量聯合概率計算步驟主要有［2］：①單變量的邊際分布參數估計和擬合度檢驗；②計算單變量的概率分布值；③變量間的相依性度量指標計算；④選擇Copula函數；⑤Copula函數參數估計和擬合度檢驗；⑥多變量聯合概率、條件概率以及相應重現期計算。武漢大學郭生練教授和熊立華教授在國內率先開展了基于Copula函數的水文分析計算［1-4］。2017年8月，國際統計水文學委員會副主席、河海大學陳元芳教授主持舉辦了河海大學“Copula在水文與環境科學中的應用”國際課程班，培訓期間，國際水文科學協會國際統計水文學委員會主席Salvatore Grimaldi教授來華授課。所有這些工作推動了Copula在水文科學領域的應用。目前，國內代表性的專著有《Copulas函數及其在水文中的應用》《Copula函數理論在多變量水文分析計算中的應用研究》《Copulas and Its Application in Hydrology and Water Resources》以及《非一致性水文概率分布估計理論和方法》等。

基于Copula函數的多變量水文頻率已有20多a的研究和實踐歷史，是水文學研究的活躍領域之一，經過國內外水文科學工作者的不懈努力，豐富了Copula函數的理論體系，拓寬了Copula函數的應用領域。基于Copula函數的多變量水文頻率計算也是一個嶄新的研究領域，涉及高等數學、概率論與數理統計、水文學、數值計算、優化計算等學科的交叉和滲透，面臨一系列亟待解決的科學問題。本文不再重復歸納總結水文頻率計算研究進展，而是重點探討Copula函數進行水文多變量聯合概率計算的幾個問題，以期促進基于Copula函數的多變量水文頻率計算理論的進一步發展。

1 Copula函數的主要類型

Copula函數分類很多。按照參數的多少，可以分為單參數Copula函數和多參數Copula函數。按照變量間相依性特性，可以分為對稱Copula函數和非對稱Copula函數。本節主要介紹水文中應用的幾種Copula函數。

1.1 對稱Archimedean Copula函數

對稱 Archimedean Copula（Symmetric Archimedean Copula）也稱為可交換 Copula（Exchangeable Archimedean Copula，EAC）。 對稱 Archimedean Copula函數主要有20種。這類函數的特點是Copula函數僅含有一個參數，函數構造簡單，相對而言求解容易［2］。其中：Cook-Johnson（Clayton）Copula函數、Gumbel-Hougaard Copula函數、Ali-Mikhail-Haq Copula函數和Frank Archimedean Copula函數是水文多變量頻率計算中應用最多的函數。這4種常見的對稱Archimedean Copula函數的生成函數φ（t）、參數θ取值范圍和Copula分布函數（CDF）表達式見表1（t為生成函數自變量，d為Copula函數維數，uj為第j個變量的邊際分布，R為實數）。

表1 常見的4種對稱Archimedean Copula函數形式

因為對稱Archimedean Copula函數只含有一個參數，僅用一個生成函數描述正的相依性，所以要求變量為對稱相依。實際上，這種假定通常是不合理的。如洪水事件可用洪量、洪水歷時和洪峰流量來描述，顯然，它們任意兩兩變量間的相依性是不對等的。表1中二維對稱Archimedean Copula函數參數與Kendall’s相關系數間存在一定的函數關系，可先計算出樣本的Kendall’s相關系數，再推算Copula函數參數。但是，三維以上的對稱Archimedean Copula函數需要利用極大似然法，通過迭代法計算獲得Copula函數參數。

1.2 非對稱Archimedean Copula函數

非對稱Archimedean Copula函數通過若干個生成函數描述變量間的相依性，克服了對稱Archimedean Copula函數的缺陷，是水文多變量分析計算較為合理的選擇函數，這些函數參數的計算需要利用極大似然法。目前，水文分析中最為常用的非對稱Archimedean Copula函數有嵌套 Archimedean（Nested Archimedean Copula，簡稱NAC）構造函數、層次Archimedean Copula（Hierarchical Archimedean Copula，簡稱 HAC）構造函數和配對 Copula（Pair-Copula Copula，簡稱 PCC）構造函數等［2］。

1.2.1 嵌套Archimedean構造函數

這類函數主要有M3函數、M4函數、M5函數、M6函數和 M12函數，見式（1）～式（4）。

（1）M3函數。其表達式為

式中：u1、u2、u3分別為邊際分布值；θ1、θ2分別為 Copula函數參數， θ2≥ θ1∈ ［0，∞） 。

（2）M4函數。其表達式為

式中： θ2≥ θ1∈ ［0，∞） 。

（3）M5函數。其表達式為

式中： θ2≥ θ1∈ ［1，∞） 。

（4）M6函數。其表達式為

式中： θ2≥ θ1∈ ［1，∞） 。

（5）M12函數。其表達式為

式中： θ2≥ θ1∈ ［1，∞） 。

1.2.2 層次Archimedean Copula構造函數

層次 Archimedean Copula（Hierarchical Archimedean Copulas）構造函數又稱為廣義嵌套Archimedean構造函數，其一般形式主要基于嵌套多元Archimedean Copula函數的框架。每一級Archimedean Copula函數由前一級聚集構成，頂層級最后以層次Archimedean Copula函數結束，形成d維標準均勻隨機變量（U1，…，Ud）的聯合分布，聯合分布值可在點 u＝（u1，…，ud）∈［0，1］d估算［2］。

1.2.3 配對Pair-Copula構造函數

配對Pair-Copula構造函數是把多元密度分解為d（d-1）/2個二維Copula密度函數，其中前d-1個為無條件Copula函數，其余為條件Copula函數，主要有Canonical vines和 D-vines兩類結構［2］。 其中：D-vines結構的表達式為

Canonical vines結構的表達式為

式中：fk為第k個變量的邊際概率密度；xk為第k個邊際變量；ci，i+j為第i個邊際變量與第i+j個邊際變量構成的Copula密度函數；xi+j為第i+j個邊際變量；xj為第c14｜23表示 c14｜23（F（x1｜x2，x3），F（x4｜x2，x3）） 。

1.2.4 三維Plackett Copula函數

設三變量（u，v，w），其兩兩交乘比率分別為二維Copula分布 CUV、CVW和 CUW，則三維 Plackett Copula參數 ψUVW定義為［2］

其中

式中：在給定ψUVW（u，v，w）下，ψUV、ψVW和ψUW分別為二維Plackett Copula函數CUV、CVW和CUW的參數。

記z＝ CUVW，則三維Plackett Copulas分布CUVW（u，v，w）為

其中

給定ψUV、ψVW、ψUW和ψUVW，三維Plackett Copula分布 CUVW（ u，v，w）可用式（8）、式（9）表示，但是三維Plackett Copula的分布密度較為復雜，其參數可采用極大似然法進行計算。不難看出，三維Plackett Copula分布有3個參數，屬于不對稱Copula函數。

1.2.5 Metaelliptical Copula函數

參數為 μ （p×1） 和Σ（p×p）（p為邊際變量維數）的d維隨機變量z具有elliptical類分布（elliptical distribution，ED），可定義為［2］

式中：r≥0，為隨機變量；u為Rd上的均勻分布變量，且獨立于r；A為d×d的常數矩陣，且滿足A AT＝Σ；符號“＝d”的含義是式（12）兩邊具有相同的分布。

當r有密度函數時，z的密度函數可以表示為

式中g（·）為一個尺度函數（Scale Function）。

常見的d維對稱elliptical類分布見表2。

表2 常見的d維對稱elliptical類分布

從表2可以看出，elliptical類分布函數非常復雜，也屬于不對稱Copula函數，只能采用極大似然法進行參數求解。

2 Copula函數在水文多變量分析中面臨的幾個問題

根據文獻報道和應用實踐，Copula函數仍然處在不斷發展和完善的階段，基于Copula函數的水文多變量分析計算中面臨的主要問題有以下幾方面。

2.1 變量邊際分布值計算

按照兩階段法估計參數，Copula函數計算首先需要進行邊際變量的分布函數值計算。由于水文事件的概率（頻率）分布函數未知，氣候變化和高強度人類活動影響使水文數據難以滿足獨立和同分布條件，以及水文序列長度有限等影響，邊際變量的分布參數值計算受到許多挑戰［5-6］，因此邊際變量函數值計算是影響Copula函數選擇和參數計算的關鍵步驟。

2.2 Copula函數參數計算

選擇Pearson古典相關系數rn、Spearman秩相關系數ρn、Kendallτ系數、Chi圖和K圖方法進行變量間的相依性度量。若相依性存在，則可計算Copula函數參數，否則按變量獨立進行多變量聯合概率計算。除表2常見的4種二維對稱 Archimedean Copula外，Copula函數參數需要應用極大似然法進行估計。一般來說，對稱 Archimedean Copulas、非對稱 Archimedean copulas和Plackette Copula分布函數已知，但是其密度函數的推求過程比較復雜。因此，正確的樣本對數極大似然函數參數的偏導數方程組是保證極大似然法正確求解的關鍵。Metaelliptical Copula函數的密度函數見表2，以下列出代表性Copula函數的密度函數表達式［2］。

2.2.1 常見的對稱三維Copula函數的密度函數

（1）Gumbel-Hougaard Copula。 其密度函數表達式為

式中：w ＝ （-ln u1）θ+ （-ln u2）θ+ （-ln u3）θ。

（2） Clayton （Cook-Johnson） Copula。 其密度函數表達式為

（3）Frank Copula。其密度函數表達式為

2.2.2 常見的非對稱Archimedean Copula函數的密度函數

（1）M3 Copula。其密度函數表達式為

式中： w ＝ 1- （1-e-θ2）-1（1-e-θ2u1）（1-e-θ2u2） ； G ＝1-e-θ1。

（2）M4 Copula。其密度函數表達式為

（3）M5 Copula。其密度函數表達式為

（4）M6 Copula。其密度函數表達式為

（5）M12 Copula。 其密度函數表達式為

2.2.3 PlacketteCopula密度函數

對于二維PlacketteCopula函數，其密度函數為

式中：ψUV為二維Plackette Copula參數。

對于三維Plackette Copula函數，其密度函數推導極其復雜，函數表達式較長，本文不再列出，具體表達式見文獻［2］。

2.3 Copula函數值計算

Copula函數參數求解后，對稱 Archimedean Copula、PlacketteCopula和完全嵌套的非對稱Archimedean Copula的分布函數具有顯函數的表達式，可直接進行Copula函數值計算。但是，非對稱Archimedean Copula函數中的三維 Pair Copula函數以及Metaelliptical Copula函數需要數值積分才能獲得Copula 函數值［2］。

2.3.1 三維Pair Copula函數值計算

三維Pair Copula函數構造靈活，但是其分布函數值計算起來比較困難，可采用數值積分計算：

通過積分變量代換，有

顯然式（23）為一個二維條件Copula分布的一維積分，可應用高斯數值積分求解出三維概率分布。

2.3.2 Metaelliptical Copula函數值計算

這類Copula函數的密度函數含有特殊類函數，無法表示為顯函數形式，只能通過數值積分進行求解計算。Kotz和 Nadarajah（2001）、Nadarajah和 Kotz（2007）推導了二維對稱Kotz type分布密度和分布的超幾何級數，二維Pearson typeⅡ、Ⅶ類不完全Beta函數的邊際分布表達式，詳細計算步驟見文獻［2］。

2.4 Copula函數的選擇

Copula函數實際上是把幾個相依邊際變量的分布函數值［0，1］通過某一函數連接起來的函數，本身不具有嚴格的水文物理基礎，也就是說Copula函數同樣與單變量分布函數一樣，它們都不能從物理意義上解釋幾個相依邊際變量的聯合發生概率。因此，Copula函數只能根據計算概率與實測數據變量的聯合經驗概率擬合效果進行定性評估，在此基礎上應用分布擬合度法進行假設檢驗。聯合經驗頻率與相依邊際變量的聯合分布有關，實際應用中，聯合經驗概率可采用Gringorten公式。顯然，這是一種聯合經驗頻率的近似計算。

Copula擬合度檢驗通用的方法為CPI Rosenblatt轉換法。 設邊際分布 FX1（x1），…，FXd（xd） 有聯合分布 Copula 函數C（FX1（x1），…，FXd（xd））＝F（x1，x2，…，xd） ，按照 Copula函數模擬， Φ-1（Zi） （其中 i＝ 1，2，…，d）服從標準正態分布 N（0，1） ， Φ 為標準正態分布； Φ-1（Zi） 為標準正態分布 N（0，1） 的逆函數，則Rosenblatt轉換檢驗法的步驟：①提出原假設 H0，即（X1，X2，…，Xd） 具有 C（FX1（x1），…，FXd（xd））＝ C（u1，u2，…，ud）；②選擇檢驗方法，如 Kolmogorov檢驗、Cramér von Mises檢驗和Anderson Darling（AD）檢驗；③根據顯著水平α，確定相應的臨界值；④根據樣本計算統計量的觀測值；⑤比較統計量的觀測值與臨界值，對原假設H0進行判斷。

在單變量分布擬合度檢驗中，當實測數據序列長度較大時，樣本統計量分布顯著水平α相應的臨界值需要提取模擬序列第（1-α）%的分位數來確定。同樣，Copula擬合度檢驗沒有給出顯著水平α相應的臨界值表，也需要進行模擬試驗進行確定。

3 基于Copula函數的水文多變量分析計算的建議

根據現有基于Copula函數的水文多變量分析面臨的問題，建議深入開展以下研究工作，以期為我國水文頻率計算提供參考和支撐。

3.1 非平穩多變量聯合概率分布與設計值計算的實用化

基于Copula函數的水文多變量聯合概率計算仍然假定邊際變量滿足獨立、同分布條件。若邊際變量不滿足平穩性，則現有Copula函數的計算方法體系不能用于計算非平穩變量的聯合概率與設計值。雖然一些學者提出了時變Copula函數的計算方法，但是其模型求解難度增大。因此，這些方法仍需在實用化方面進一步研究。另外，單變量情況下設計標準（設計重現期）的設計值可以進行清晰的確定并廣泛應用于工程實踐，但是多變量設計值則無法明確確定。Salvadori G等提出了基于聯合概率值定義危險區域的Kendall重現期計算方法［47-50］，他們認為該方法比傳統的重現期計算更為合理，其基本原理是將事件發生區域劃分為超臨界區域、臨界面和亞臨界區域3類，事件發生在超臨界區域上的平均時間間隔長度為多變量事件重現期。這種計算非常復雜，目前仍然停留在研究層面，需要在實用化方面進一步研究。

3.2 提高Copula函數的參數估計精度

如上所述，Copula函數的參數估計方法大多采用極大似然法。按照極大似然法原理，一般需要求解樣本似然函數對Copula函數參數偏導數的非線性方程組。若Copula函數參數較少、不含有特殊函數時，則非線性方程組求解相對容易。當Copula函數參數較多、非線性方程組含有特殊函數和邊際變量非顯式分布函數時，則求解非常復雜。另外，非線性方程組的初值選擇有時會影響非線性方程組的求解。因此，當非線性方程組求解困難時，以樣本似然函數為最大目標函數，考慮Copula函數參數取值范圍，采用遺傳算法、粒子群算法、螞蟻群算法、混合蛙跳算法、頭腦風暴算法、蜻蜓算法和水循環算法等智能優化求解方法，可能不失為一種求解途徑。

3.3 提高隨機變量和、差、積、商分布計算的精度

流域設計斷面各部分洪水的地區組成、干旱分析中特征變量的復合運算、水資源評價和配置中研究斷面來水與區間耗水量的組合分布分析等都可以歸結為水文變量和、差、積、商的分布計算。傳統多變量理論上可以通過數學變換和高維數值積分來獲得，但是其求解非常復雜，難以保證積分值的可靠性。因此，應用Copula函數原理，研究提高水文變量和、差、積、商的分布計算精度是目前急需解決的科學問題。