基于藤Copula多維聯(lián)合分布的CMIP5多模式降雨綜合方法研究

胡義明，羅序義，梁忠民，劉 楊

（河海大學水文水資源學院，南京210098）

0 引 言

全球氣候模式是對氣候進行模擬和預測的重要工具之一［1］，常被用于未來水文情勢的分析預估。耦合模式比較計劃第五階段（CMIP5）新增了一些模式試驗，采用的耦合器技術、通量處理方案及參數(shù)化方案比CMIP3的更為合理，在一定程度上提高了氣候模式的模擬和預估能力［2］。雖然氣候模式在不斷地改進和發(fā)展，但模式的輸出值與真實值之間還是存在較大偏差。Xu 等［3］評估了CMIP5 中18 個GCMs 降水數(shù)據(jù)對中國地區(qū)降水變化的模擬效果，發(fā)現(xiàn)大多數(shù)模式的降水數(shù)據(jù)偏高，且對降水在時間上變化的模擬能力有限。李振朝等［4］評估了CMIP5中8 個GCMs 對青藏高原地區(qū)降水的預估效果，結果表明各模式的年降水預估值均大于實測降水，且預估序列與實測序列間的相關性很小。因此，通過多模式集成綜合以提升GCM的模擬精度，成為常用研究方法。許崇海等［5］評估了22 個氣候模式對東亞地區(qū)氣候的模擬能力，表明多模式集合平均的模擬能力比單一模式更優(yōu)。Santer 等［6］考慮到多模式集合結果通常優(yōu)于單一模式，且不同模式的模擬能力不同，建議使用加權多模式集合（Multi-Model Ensemble）方法對多模式結果進行綜合。王晶［7］利用8個氣候模式的回報資料，采用BMA進行多模式集成，比較發(fā)現(xiàn)BMA 集成后的回報效果優(yōu)于單一模式和多模式集合平均。

然而上述研究沒有考慮不同模式之間的依賴結構，而模式間的相依關系對模式集成效果是有影響的，通過建立聯(lián)合分布可以考慮變量之間的相關性。Copula 函數(shù)是一種靈活的統(tǒng)計工具，允許不同變量有不同的邊緣分布，可用于建立兩個或多個隨機變量之間的聯(lián)合概率分布［8］，其中，藤Copula在描述多變量間復雜相關結構時具有更好的靈活性［9，10］。Pham等［11］基于C藤Copula 函數(shù)建立了隨機蒸散發(fā)模型，得到了與降雨一致的日參考蒸散發(fā)隨機模擬序列。Wang 等［12］基于D 藤Copula 函數(shù)，以與當月徑流相關的歷史徑流和氣候要素（降雨、溫度、蒸發(fā)等）構建條件聯(lián)合分布，較好地模擬了黃河源區(qū)月徑流。

本文采用C 藤Copula 建立GCM 多模式集成方法。應用泰勒圖技術從12 個GCM 模式中挑選6 個較好的模式，構建該6 個模式模擬序列與實測序列的聯(lián)合分布，推求給定模式數(shù)據(jù)時實測值的條件分布，據(jù)此實現(xiàn)多模式結果的集成。采用相關系數(shù)（r）和均方根誤差（RMSE）對集成結果的擬合效果進行評價，并將其與單一模式、多模式集合平均進行對比分析。

1 模型與方法

1.1 GCM模式及優(yōu)選

采用CMIP5 的12 個GCMs 歷史模擬實驗的降雨模擬數(shù)據(jù)，模擬數(shù)據(jù)時段為1961-2005年，12 個GCMs 模式［模式編號］分別為：BCC-CSM1.1［m1］、BCC-CSM-1.1（m）［m2］、CCSM4［m3］、CSIRO-Mk3.6.0［m4］、GFDL-CM3［m5］、GFDL-ESM2G［m6］、IPSL-CM5A-MR［m7］、MIROC5［m8］、MIROC-ESM［m9］、MIROC-ESM-CHEM［m10］、MRI-CGCM3［m11］、NorESM1-M［m12］。

對上述12 個GCMs 模式，采用泰勒圖方法［13］評估各模式對實測系列的模擬能力，并從中優(yōu)選6 個最優(yōu)的模式。泰勒圖方法是將模式模擬系列與實測系列之間的相關系數(shù)、各自標準差及均方根誤差繪在同一張圖上，可以直觀地反映模式模擬值對實測值的模擬效果，在氣候模式的評估研究中應用廣泛［13-15］。本文采用的泰勒圖為標準化泰勒圖，即將模式模擬系列的標準差、均方根誤差除以實測系列的標準差［15］。同時，引入一個綜合評價指標Ms［16］來評價模式的綜合模擬能力。

式中：n為選用的GCMs 個數(shù)（n=12）；k為評估指標個數(shù)（k=3）；ranki為各模式在第i個評估指標下的排名（模擬能力越強排名越靠前，rank值越小）。

顯然，模式的綜合模擬能力越強，其Ms的值越接近于1，或反之［17］。

1.2 藤Copula

Copula 函數(shù)［18，19］是一個定義在［0，1］區(qū)間上均勻分布的多維聯(lián)合分布函數(shù)，通過連接多個隨機變量任意形式的邊緣分布構建聯(lián)合分布。假設n維隨機變量x1，…，xn的聯(lián)合分布函數(shù)為F(x1，x2，…，xn)，則存在一個Copula函數(shù)C使得：

式中：ui=Fi(xi)，i= 1，…，n；Fi(xi)為隨機變量xi的邊緣分布函數(shù)。

假如Fi(xi)，i= 1，…，n都是連續(xù)的，那么存在唯一的C。多維聯(lián)合分布密度函數(shù)f(x1，…，xn)可以表示為：

式中：fi(xi)為隨機變量xi的邊緣密度函數(shù)；c(u1，u2，…，un)為Copula密度函數(shù)。

當Copula函數(shù)絕對連續(xù)時，c(u1，u2，…，un)=?nC(u1，u2，…，un)/?u1?u2…?un。

藤Copula 函數(shù)是一種高維Copula 函數(shù)，能更好地描述高維變量之間的復雜相依關系。對于高維聯(lián)合分布，存在很多可能的Pair-Copula 構造方式［20］。為了便于組織Pair-Copula，Bed?ford 和Cooke［21］引入了正則藤圖形建模工具。一個d維隨機變量的正則藤Copula，有d-1 棵樹，第k棵樹中有d+1-k個節(jié)點和d-k條邊。其中第一棵樹中一個節(jié)點代表一個邊緣分布，一條邊對應一組Pair-Copula，除第一棵樹外，其余樹節(jié)點均為上一棵樹的邊。

C 藤和D 藤是兩種特殊的正則藤，其中C 藤為星型結構，一個中心節(jié)點連接了其他所有節(jié)點；D藤為平行直線型結構，一個節(jié)點最多與兩條邊連接［20］。藤結構邏輯不同，適用的數(shù)據(jù)集類型也有所不同。D 藤適合用于描述變量相互獨立的數(shù)據(jù)集，而當數(shù)據(jù)集含有主導其他變量的關鍵變量時C藤更合適。

根據(jù)Bedford 和Cooke［22］的表述，一個n維C 藤、D 藤Copula的密度函數(shù)表達式分別如下：

式中：n為隨機變量的維數(shù)；fk(xk)為邊緣密度函數(shù)；j為樹的棵數(shù)；i為樹上邊的索引；c?|?(?，?)為藤結構中的二維Copula 的密度函數(shù)；F(?|?)為條件分布函數(shù)。

1.3 條件分布

根據(jù)Joe［23］的表述，式（4）和（5）中的條件分布函數(shù)可以用如下通用公式表示：

式中：v為一個n維的向量；vj為v中的任意一個變量；v-j為v去除vj后剩余向量；Cx，vj|v-j為二維條件Copula 函數(shù)，當v的維數(shù)為1時：

當x和v都服從U[0，1]時，即F(x)=x，F(xiàn)(v)=v，f(x)=f(v)= 1時，可以用h函數(shù)表示F(x|v)：

式中：θ為求x，v之間聯(lián)合分布時Copula函數(shù)的參數(shù)。

Joe［24］綜合介紹了二維單變量的h函數(shù)及其反函數(shù)。根據(jù)式（6），C藤Copula和D 藤Copula的條件分布函數(shù)表達式分別為式（9）和（10）：

以一個簡單的4 維C 藤Copula 條件分布函數(shù)為例，根據(jù)式（6）～（9），F(xiàn)(x4|x1，x2，x3)可以表示如下：

1.4 隨機模擬

根據(jù)藤Copula 模型建立的實測序列與模式模擬序列之間的聯(lián)合分布，可以得到給定模式模擬值時實測值的條件分布函數(shù)，并可取條件分布的均值作為綜合后的模擬值。假設第m年實測值為Pm，同期d個模式的模擬值分別為Pm1，Pm2，…，Pmd，則實測值可以用條件分布的反函數(shù)和隨機數(shù)得到：

式中：wm為隨機數(shù)，wm～U[0，1]；Fm-1(?|?)為第m年給定模式模擬值時實測值的條件分布函數(shù)的反函數(shù)。

1.5 模擬效果評價

采用Pearson 線性相關系數(shù)（r）、均方根誤差(RMSE）評價單一模式、多模式集合平均及藤Copula 多模式綜合對實測系列的模擬效果。假設xi，i= 1，2，…，n為實測序列，對應的模擬序列為yi，i= 1，2，…，n，則上述兩個評價指標的表達式分別為：

2 實例研究

梧州水文站是珠江流域西江干流的主要控制站，位于111°20′E，23°28′N，控制廣西境內85%的集水面積和年徑流總量［25］。西江流域位于我國華南沿海，屬于熱帶、亞熱帶季風氣候，降水量的年際變化較大，雨季長且年內分配極不均勻，降水主要集中在4-9月［26］。本次研究以西江梧州站以上流域1961-2005年7月面平均雨量作為研究對象。

2.1 模式優(yōu)選

構建藤Copula 模型時，隨著維數(shù)的增加，可能的藤結構個數(shù)迅速增加，一個d維隨機變量構成的藤Copula，有種可能的C（D）藤，種可能的R 藤［27］。為了避免由于維數(shù)過高導致計算效率不高，同時又盡可能多地考慮多個不同模式的信息，通過泰勒圖指標綜合選用了6個模式。評價原則是：模式模擬值與實測值相關系數(shù)越大、均方根誤差與實測標準差之比越小、標準差與實測標準差之比越接近1，模擬效果越好［28］。圖1給出了12個模式模擬降雨相對于實測值的標準化泰勒圖。

圖1 12個GCMs模擬降雨相對于實測值的泰勒圖Fig.1 Taylor diagram of simulated precipitation of 12 GCMs relative to observed values

基于泰勒圖結果，計算了各模式的綜合評價指標Ms。根據(jù)表1中的計算結果最終選擇了編號為m1，m2，m3，m4，m9，m10的6個GCM模式。

表1 12個GCMs的綜合評價指標Tab.1 Comprehensive evaluation index of 12 GCMs

圖2給出了1961-2005年實測、6個模式原始模擬系列和多模式平均模擬系列。從圖2中可以看出：不同模式模擬值序列差異較大，單個模式的模擬效果不甚理想，而將多個模式進行集合平均處理，又使得模擬序列過于平滑，顯著削弱了對極值的估計，不能反映真實月降雨的年際變化情況。

圖2 實測降雨和模式模擬降雨時間序列Fig.2 Observed precipitation and model simulated precipitation time series

2.2 最優(yōu)邊緣分布擬合

通過藤Copula 建立聯(lián)合分布的第一步是選擇擬合最優(yōu)的邊緣分布。本次研究的是月降雨數(shù)據(jù)，選取了常用的GEV、Weibull、Lognormal、Normal和PE3分布函數(shù)作為備選邊緣分布。采用K-S 檢驗法在P＝0.05 顯著性水平下對各個分布函數(shù)的擬合效果進行了分析。以K-S 檢驗結果中的P值優(yōu)選分布，當D值比0.05 顯著性水平下D的臨界值小時，表明通過顯著性檢驗。各模式最終確定的邊緣分布及K-S檢驗結果見表2。

表2 選用邊緣分布及K-S檢驗結果Tab.2 Selection of edge distribution and K-S test results

2.3 建立藤Copula模型

通過藤Copula 建立實測序列與模式模擬序列之間的聯(lián)合分布，采用在水文領域常用的Archimedean Copula 函數(shù)［29］中的Clayton、Gumbel、Frank 型，以及它們的旋轉結構作為二維Copu?la的選用類型。以AIC準則選出組成藤Copula的最優(yōu)二維Cop?ula 類型，采用極大似然法估計相應參數(shù)。以45年模式模擬值和實測值的邊緣分布作為輸入，分別建立C 藤Copula 和D 藤Copula 模型。依據(jù)AIC 準則、BIC 準則和對數(shù)似然值（LL）［30］，確定最優(yōu)擬合的模型，見表3。采用留一交叉驗證方法進行效果分析，具體為：對于總共45年的數(shù)據(jù)資料，每次采用44 a數(shù)據(jù)構建模型，在構建的模型基礎上對剩余1 a 的實測值進行模擬，重復進行45 次，最后將得到的模擬系列與實測系列進行對比分析，評價模擬效果。對于AIC 和BIC 而言，C 藤Copula 的指標結果更小，并且C 藤Copula 的對數(shù)似然值更大。為此，本研究采用擬合更優(yōu)的C藤Copula模型。

表3 藤Copula擬合優(yōu)度檢驗結果Tab.3 Goodness of fit test results of Vine Copula

2.4 隨機模擬

根據(jù)C 藤Copula 模型構建的實測與模式模擬之間的聯(lián)合分布，在給定GCMs模擬值的情況下，可以獲得實測值的條件分布。由于難以獲得條件分布的解析解，為此采用隨機模擬的方式進行求解。具體為針對每年的條件分布函數(shù)采用馬爾科夫鏈蒙特卡洛方法從條件分布函數(shù)中隨機抽取1×104個樣本，計算1×104個樣本的均值并將其作為對應年份的月降雨模擬值，隨機模擬系列及其均值模擬序列如圖3所示。

圖3 隨機模擬值及其均值與實測值時間序列對比圖Fig.3 Time series comparison diagram of random simulated values and its mean value with observed values

采用相關系數(shù)和均方根誤差指標對均值模擬序列的模擬效果進行評價，結果如圖4所示。

從圖4中可以看出，通過采用C 藤Copula 集成多模式模擬結果得到的均值模擬系列與實測系列間的相關系數(shù)比6個模式的都大，從單一模式最大的0.286 提高至0.363；對于均方根誤差，由93.487 減少至72.695，均優(yōu)于多模式集合平均的兩個指標評價結果（0.050 和77.903）。由此表明，基于C 藤Copula 的多模式綜合方法得到的模擬值比任意單一模式表現(xiàn)更好，且優(yōu)于多模式集合平均，起到了改善模擬效果的作用。雖然改善的幅度受限于選用的單一模式，但是該方法可以在多模式的基礎上提高對月降雨數(shù)據(jù)的擬合效果。

圖4 單一模式、多模式集合平均及多模式綜合的指標比較Fig.4 Index comparison of single-model，multi-model ensemble mean and multi-model synthesis

3 結 論

本文基于12 個GCMs 模式提供的模擬降雨，采用泰勒圖指標優(yōu)選了6 個模擬能力較好的GCMs，通過藤Copula 構建了CMIP5月降水量多模式綜合校正模型，并以西江梧州站以上流域的面雨量為對象進行了示例研究。

（1）基于AIC、BIC 和LL 指標，評估了C 藤Copula 和D 藤Copula 模型的應用效果，結果表明C 藤Copula 的校正效果要優(yōu)于D藤Copula。

（2）基于相關系數(shù)和均方根誤差指標，分析了C 藤Copula校正結果與任一原始模式結果的差異。結果表明校正結果優(yōu)于任一模式的原始結果，也優(yōu)于多模式的平均結果。此外，通過多模式綜合校正后的系列，比多模式平均的系列更能反映真實月降雨的年際變化趨勢。

（3）藤Copula 在應用過程中，其結構數(shù)目會隨著變量維數(shù)的增加而急劇增加，會降低計算效率。但若采用的GCMs 模式數(shù)目太少，又可能導致數(shù)據(jù)信息利用不足。為此，在保證校正效果前提下，如何平衡GCM數(shù)目選取和藤Copula高效構建有待進一步深入研究。 □