基于Copula模型的暴雨綜合指標(biāo)分析
——以萬州為例

馮介玲 祝好

（1 重慶市氣象科學(xué)研究所，重慶 401147；2 重慶市農(nóng)業(yè)氣象與衛(wèi)星遙感工程技術(shù)研究中心，重慶 401147）

0 引言

IPCC第五次評(píng)估報(bào)告提出，由氣候變化引起大多數(shù)地區(qū)強(qiáng)降水事件的發(fā)生頻率呈現(xiàn)增加趨勢(shì)，發(fā)生影響大和范圍廣的暴雨災(zāi)害的可能性也會(huì)大幅度增加，嚴(yán)重地影響了人民的生命和財(cái)產(chǎn)安全以及社會(huì)的可持續(xù)發(fā)展。在全球變化的背景下，暴雨災(zāi)害的發(fā)生過程更加復(fù)雜，暴雨事件通常包含多個(gè)方面的特征屬性，是多重指標(biāo)綜合作用的結(jié)果，各個(gè)要素之間的相關(guān)性也趨于復(fù)雜，因此，通過單一要素的分析很難全面地反映暴雨特征，暴雨多要素綜合分析成為極端降水研究領(lǐng)域的熱點(diǎn)。

Copula函數(shù)具有變量不受邊緣分布形式限制，不基于線性假設(shè)等特點(diǎn)，能夠靈活地應(yīng)用于多維變量的聯(lián)合特征分析中，在金融、保險(xiǎn)行業(yè)應(yīng)用成熟，而近年來也被廣泛地運(yùn)用在氣象和水文等方面的研究中，但是在暴雨方面的研究相對(duì)較少。重現(xiàn)期是一個(gè)表征發(fā)生可能性的指標(biāo)，一般是指在較長(zhǎng)時(shí)期內(nèi)超過某強(qiáng)度閾值的事件發(fā)生的一個(gè)相對(duì)穩(wěn)定的周期，是描述災(zāi)害等事件的一個(gè)重要參數(shù)。重現(xiàn)期可用于描述災(zāi)害出現(xiàn)的嚴(yán)重程度和頻率，氣象或水文領(lǐng)域上常用“50年一遇”“百年一遇”等來表示。研究暴雨重現(xiàn)期對(duì)于城市工程的規(guī)劃設(shè)計(jì)以及防汛抗洪等都具有重要的參考意義。

本文以重慶萬州區(qū)暴雨事件為研究對(duì)象，利用Copula函數(shù)在構(gòu)建多維聯(lián)合分布中的優(yōu)勢(shì)，構(gòu)建二維聯(lián)合分布模型研究歷年暴雨多變量的概率分布和聯(lián)合重現(xiàn)期的變化，揭示多變量聯(lián)合下的暴雨特征，以期利用綜合指標(biāo)分析提高暴雨重現(xiàn)期的評(píng)估精度，為減小暴雨災(zāi)害頻率分析的不確定性，制定合理的防災(zāi)減災(zāi)決策以及提高風(fēng)險(xiǎn)管理水平提供科學(xué)依據(jù)。

1 資料與方法

1.1 資料選取

根據(jù)中國(guó)氣象局《降水量等級(jí)》標(biāo)準(zhǔn)，24 h內(nèi)，降雨量達(dá)到50 mm的降水定為暴雨。研究采用由重慶氣候監(jiān)測(cè)分析平臺(tái)提供的萬州區(qū)1967—2019年的歷年降水資料，參照國(guó)內(nèi)外學(xué)者研究暴雨特征的指標(biāo)選取情況，統(tǒng)計(jì)提取出年暴雨量（

P

）、年暴雨強(qiáng)度（

I

）和年暴雨貢獻(xiàn)率（

R

）作為暴雨要素，從多角度反映暴雨的特征，相應(yīng)暴雨要素的選擇及定義如表1所示。

表1 暴雨要素的選擇及定義Table 1 Selection and definition of rainstorm factors

1.2 研究方法

1.2.1 Copula函數(shù)的基本形式

在實(shí)際應(yīng)用中，單參數(shù)的Copula函數(shù)由于靈活多變，計(jì)算簡(jiǎn)單，應(yīng)用最為廣泛。建立Copula函數(shù)模型時(shí)，需要采極大似然法和分步估計(jì)法對(duì)Copula函數(shù)的參數(shù)估計(jì)，三種常見的Copula函數(shù)的表達(dá)式及參數(shù)范圍見表2。其中

C

(

u

,

v

)為Copula函數(shù)，

u

和

v

分別為變量的邊緣分布函數(shù)，

θ

為Copula函數(shù)的參數(shù)。

通過計(jì)算擬合優(yōu)度評(píng)價(jià)指標(biāo)RMSE，AIC，Bias和OLS值，可以對(duì)擬合的Copula函數(shù)模型進(jìn)行優(yōu)度評(píng)價(jià)。

RMSE的表達(dá)式為：

表2 三種常見Copula函數(shù)Table 2 Definition of the three Copula functions

AIC的計(jì)算式如下：

Bias的表達(dá)式為：上述公式中，

n

是樣本容量，

m

為Copula函數(shù)模型中獨(dú)立參數(shù)的個(gè)數(shù)，

P

，

P

分別為聯(lián)合分布函數(shù)的經(jīng)驗(yàn)概率和理論概率。

根據(jù)RMSE，AIC，Bias和OLS值越小，所對(duì)應(yīng)的Copula函數(shù)擬合優(yōu)度越高的原則，選取最優(yōu)的Copula函數(shù)來進(jìn)一步構(gòu)建聯(lián)合重現(xiàn)期模型。

1.2.2 基于Copula函數(shù)的重現(xiàn)期計(jì)算

如果

N

是暴雨事件的時(shí)間序列長(zhǎng)度，

n

是事件的發(fā)生次數(shù)，

E

(

L

)為暴雨事件發(fā)生的平均間隔時(shí)間，則有

E

(

L

) =

N

/

n

。

對(duì)于暴雨事件而言，傳統(tǒng)的基于單要素的重現(xiàn)期的計(jì)算公式為：

基于Copula函數(shù)的雙要素聯(lián)合重現(xiàn)期可以表示為：

基于Copula函數(shù)的雙要素同現(xiàn)重現(xiàn)期可以表示為：

當(dāng)同時(shí)考慮強(qiáng)暴雨災(zāi)害的2個(gè)要素時(shí)，聯(lián)合重現(xiàn)期即為某一暴雨要素超過某一特定值時(shí)的重現(xiàn)期，而同現(xiàn)重現(xiàn)期是指2個(gè)要素都超過某個(gè)特定值的重現(xiàn)期。

2 萬州區(qū)暴雨聯(lián)合特征分析

2.1 暴雨要素特征檢驗(yàn)

對(duì)暴雨的3種要素進(jìn)行描述性統(tǒng)計(jì)分析，圖1顯示了暴雨量、暴雨強(qiáng)度和暴雨貢獻(xiàn)率的頻率直方圖以及擬合的正態(tài)分布曲線，各要素的描述性統(tǒng)計(jì)數(shù)據(jù)見表3。由圖1可知，暴雨要素的數(shù)據(jù)分布都具有一定的非對(duì)稱性，而由表3進(jìn)一步可知，與正態(tài)分布相比，3個(gè)暴雨要素總體的數(shù)據(jù)分布較為陡峭，存在右偏。根據(jù)要素的頻率直方圖形態(tài)及偏度、峰度值初步判斷數(shù)據(jù)不服從正態(tài)分布，再利用Matlab R2014a軟件分別對(duì)變量進(jìn)行Jarque-Bera正態(tài)性檢驗(yàn)，返回值都為1，則認(rèn)為數(shù)據(jù)都拒絕服從正態(tài)分布的假設(shè)，表示暴雨量、暴雨強(qiáng)度和暴雨貢獻(xiàn)率數(shù)據(jù)都具有非正態(tài)性，考慮用無需假定邊緣分布形式的Copula函數(shù)來進(jìn)行暴雨多要素的聯(lián)合分布建模。

圖1 暴雨量（a）、暴雨強(qiáng)度（b）和暴雨貢獻(xiàn)率（c）的頻率直方圖Fig.1 Frequency histogram of rainstorm factors(a) quantity, (b) intensity, (c) contribution

表3 描述性統(tǒng)計(jì)數(shù)據(jù)表Table 3 Descriptive statistical data table

2.2 邊緣分布選擇估計(jì)

令暴雨量、暴雨強(qiáng)度和暴雨貢獻(xiàn)率的邊緣分布分別為

F

(

p

)、

F

(

i

)和

F

(

r

)，利用Easyfit5.5軟件對(duì)各暴雨要素進(jìn)行多種邊緣分布類型的擬合以及參數(shù)估計(jì)，再根據(jù)Anderson-Darling檢驗(yàn)的結(jié)果，針對(duì)各個(gè)暴雨要素分別選取最優(yōu)的邊緣分布函數(shù)。由擬合結(jié)果可知，Gen Extreme Value（廣義極值分布）對(duì)于各暴雨要素的擬合效果均較好，通過0.01水平下的假設(shè)檢驗(yàn)的各要素的最優(yōu)邊緣分布類型以及參數(shù)結(jié)果見表4。

2.3 Copula模型構(gòu)建

利用極大似然法，分別選取3種常見的Copula函數(shù)對(duì)暴雨量、暴雨強(qiáng)度和暴雨貢獻(xiàn)率的兩兩組合進(jìn)行二維Copula模型的擬合，并計(jì)算擬合優(yōu)度評(píng)價(jià)指標(biāo)，進(jìn)行優(yōu)度檢驗(yàn)，相關(guān)的參數(shù)和擬合優(yōu)度檢驗(yàn)值見表5。比較表5中各Copula函數(shù)的RMSE，AIC，Bias以及OLS值，根據(jù)擬合優(yōu)度指標(biāo)越小模型擬合度越好的擇優(yōu)原則，選取Frank Copula為暴雨強(qiáng)度和暴雨貢獻(xiàn)率的二維聯(lián)合的最優(yōu)Copula分布函數(shù)，選取Clayton Copula函數(shù)作為暴雨量和暴雨量貢獻(xiàn)率、暴雨強(qiáng)度和暴雨量貢獻(xiàn)率的二維聯(lián)合Copula分布函數(shù)。擬合優(yōu)度檢驗(yàn)的結(jié)果也表明，Copula函數(shù)對(duì)暴雨量和暴雨強(qiáng)度的擬合精度更高。

表4 邊緣分布估計(jì)結(jié)果Table 4 Distribution functions of rainstorm factors

Kendall秩相關(guān)系數(shù)可以用來度量變量變化總體的一致性程度，從而反映變量間的相關(guān)性，分別計(jì)算基于三種暴雨要素組合最優(yōu)二維聯(lián)合分布的Kendall秩相關(guān)系數(shù)，結(jié)果見表6。由表可知，暴雨強(qiáng)度和暴雨貢獻(xiàn)率的同步性較大，各暴雨要素間存在一定的正相關(guān)關(guān)系。

構(gòu)建3種要素組合的最優(yōu)聯(lián)合分布模型，如圖2所示。從圖2中可以看出，各擬合的Copula函數(shù)模型在上尾和下尾的頻率都相對(duì)較大，符合暴雨要素有較多極端值的特征。由Frank Copula模型反映出的要素上下尾頻率都較大的特點(diǎn)，表示在發(fā)生暴雨災(zāi)害的情況下，同時(shí)出現(xiàn)暴雨量小和暴雨強(qiáng)度低、暴雨強(qiáng)度低和暴雨貢獻(xiàn)率小、暴雨量大和暴雨強(qiáng)度高或者暴雨強(qiáng)度高和暴雨貢獻(xiàn)率大的概率更高，而由Clayton Copula模型反映出暴雨量和暴雨貢獻(xiàn)率兩個(gè)暴雨要素在上下尾部表現(xiàn)出一定的非對(duì)稱結(jié)構(gòu)，尤其在極端下尾呈現(xiàn)出較強(qiáng)的正相關(guān)關(guān)系。分析可知，各暴雨要素之間的尾部相關(guān)性能被不同類型的Copula函數(shù)清晰地定性刻畫出來，為了更加深入地認(rèn)識(shí)Copula函數(shù)模型在暴雨多要素綜合分析中的應(yīng)用，需要進(jìn)一步進(jìn)行基于最優(yōu)Copula函數(shù)的重現(xiàn)期評(píng)估。

表5 Copula函數(shù)的參數(shù)估計(jì)和擬合優(yōu)度檢驗(yàn)值Table 5 Parameters and goodness-of-fit measures for the fitted Copula function

表6 基于最優(yōu)Copula函數(shù)的暴雨要素組合Kendall秩相關(guān)系數(shù)值Table 6 The correlation coefficients of the combination of rainstorm factors

3 萬州區(qū)暴雨聯(lián)合重現(xiàn)期分析

3.1 暴雨重現(xiàn)期計(jì)算

根據(jù)重現(xiàn)期的計(jì)算公式（6）—（8），分別計(jì)算出基于最優(yōu)Copula函數(shù)的聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期，結(jié)果如圖3所示。

由圖3可知，隨著暴雨量、暴雨強(qiáng)度和暴雨貢獻(xiàn)率的增加，基于不同最優(yōu)Copula函數(shù)模型的聯(lián)合重現(xiàn)期及同現(xiàn)重現(xiàn)期增大。不同暴雨要素組合的聯(lián)合重現(xiàn)期模型的形態(tài)不同，總體來說聯(lián)合重現(xiàn)期對(duì)暴雨強(qiáng)度的變化更敏感。暴雨量與暴雨貢獻(xiàn)率的聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期主要分布在0～400和0～20000年，暴雨量與暴雨貢獻(xiàn)率的聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期主要分布在0～5和0～800年，暴雨強(qiáng)度和暴雨貢獻(xiàn)率的聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期主要分布在0～5和0～1500年。聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期對(duì)暴雨要素滿足的取值要求不同，不同要素組合的重現(xiàn)期變化趨勢(shì)不同，反映了在分析暴雨的重現(xiàn)期時(shí)，需要考慮多個(gè)要素的影響。

圖2 暴雨要素的Copula函數(shù)分布圖（a）暴雨量和暴雨強(qiáng)度；（b）暴雨量和暴雨貢獻(xiàn)率；（c）暴雨強(qiáng)度和暴雨貢獻(xiàn)率Fig. 2 The joint cumulative distribution of rainstorm factors(a) P&I, (b) P&R, (c) I&R

3.2 暴雨災(zāi)害聯(lián)合重現(xiàn)期評(píng)估結(jié)果比較

表7 顯示了在單變量重現(xiàn)期分別為2、5、10、20、50和100 年的情況下，根據(jù)單變量重現(xiàn)期公式（6）推導(dǎo)出暴雨量、暴雨強(qiáng)度和暴雨貢獻(xiàn)率的對(duì)應(yīng)取值，再根據(jù)公式（7）和（8）計(jì)算出的不同Copula函數(shù)模型對(duì)應(yīng)的暴雨多要素聯(lián)合重現(xiàn)期以及同現(xiàn)重現(xiàn)期取值。

由表7可以看出，聯(lián)合重現(xiàn)期均比單要素重現(xiàn)期小，同現(xiàn)重現(xiàn)期均比單要素重現(xiàn)期大，且隨著單要素重現(xiàn)期取值的增大，多要素重現(xiàn)期與單要素重現(xiàn)期的差別也逐漸增大，同現(xiàn)重現(xiàn)期的增速遠(yuǎn)超聯(lián)合重現(xiàn)期的增速。在單變量重現(xiàn)期小于等于2年的情況下，聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期的差別不明顯，但當(dāng)單要素重現(xiàn)期大于2年時(shí)，基于Copula模型的兩種重現(xiàn)期差別較大。

4 結(jié)論與討論

圖3 暴雨要素的二維Copula重現(xiàn)期Fig. 3 The joint return period of the combination of rainstorm factors

1）廣義極值分布對(duì)暴雨量、暴雨強(qiáng)度和暴雨貢獻(xiàn)率的擬合效果最好。根據(jù)邊緣分布情況以及擬合優(yōu)度評(píng)價(jià)結(jié)果，Clayton Copula函數(shù)適合于構(gòu)建暴雨量和暴雨強(qiáng)度以及暴雨強(qiáng)度和暴雨貢獻(xiàn)率兩個(gè)的聯(lián)合分布模型，F(xiàn)rank Copula函數(shù)適合于構(gòu)建暴雨強(qiáng)度和暴雨貢獻(xiàn)率的聯(lián)合分布模型。

2）各暴雨要素之間具有一定的正相關(guān)，Copula函數(shù)可以直觀清晰地刻畫要素之間的相關(guān)關(guān)系，適用于分析暴雨多要素的聯(lián)合概率。

3）Copula函數(shù)模型能夠從綜合指標(biāo)層面反映暴雨多方面的特征，基于暴雨多要素的聯(lián)合重現(xiàn)期以及同現(xiàn)重現(xiàn)期比基于單要素的重現(xiàn)期能夠提供暴雨的綜合信息，可以反映出暴雨災(zāi)害發(fā)生時(shí)的降水集中程度以及暴雨的強(qiáng)度等級(jí)。

表7 基于暴雨單要素的重現(xiàn)期和最優(yōu)Copula函數(shù)模型的暴雨多要素重現(xiàn)期Table 7 The return period of rainstorm based on single element and optimal Copula function model