基于Copula函數(shù)的高強(qiáng)度降雨設(shè)計(jì)研究

余 璐

(北京市十三陵水庫管理處,北京 102200)

全球氣候變暖和人類活動(dòng)可能引發(fā)極端氣候事件發(fā)生強(qiáng)度和頻率增加，極大地影響了流域的水文特征[1- 2]。強(qiáng)降雨事件作為極端天氣中的一種嚴(yán)重致災(zāi)的因子，容易導(dǎo)致城市洪澇、泥石流和山體滑坡等一系列災(zāi)害。因此，強(qiáng)降雨事件分析對(duì)于城市防災(zāi)減災(zāi)有重要現(xiàn)實(shí)意義[3]。

水文中通常采用降雨的概率評(píng)估水文事件的嚴(yán)重程度。通過降雨量超過一定量級(jí)的頻率來度量。然而，單一降雨要素不能表征降雨的多重屬性。流域內(nèi)一場(chǎng)降雨的降雨量越大，降雨歷時(shí)越小，則控制流域的水庫越難調(diào)節(jié)此場(chǎng)降雨，容易引起洪水災(zāi)害。現(xiàn)有的水庫防洪研究往往只考慮不同頻率的降雨量對(duì)水庫防洪的影響，而沒有考慮到降雨歷時(shí)，這容易忽略降雨量大而降雨歷時(shí)小的高強(qiáng)度降雨。運(yùn)用多變量能更好描述降雨事件，Copula函數(shù)由于可以不受邊緣分布類型拘束，靈活描述多元變量之間復(fù)雜的非線性相關(guān)結(jié)構(gòu)，在多變量分析中廣泛運(yùn)用。曹偉華等通過阿基米德Copula函數(shù)構(gòu)建降水持續(xù)時(shí)間和過程降水量的二元聯(lián)合分布模型，計(jì)算了北京地區(qū)強(qiáng)降水事件條件重現(xiàn)期[4]。何英等基于阿基米德Copula函數(shù)分析了葉爾羌河流域卡群水文站洪峰流量與洪量的聯(lián)合頻率分布特征[5]。張冬冬等應(yīng)用阿基米德Copula函數(shù)探討了洪峰、洪量和歷時(shí)的聯(lián)合概率分布和條件概率分布[6]。

北京市十三陵水庫作為溫榆河重要的防洪建筑，開展十三陵水庫控制流域多變量降雨分析對(duì)北京市防洪減災(zāi)工作具有重要意義。本文以北京市十三陵水庫流域?yàn)槔攸c(diǎn)考慮降雨量大而降雨歷時(shí)小的高強(qiáng)度降雨，采用25種常見Copula函數(shù)全面分析降雨量與降雨歷時(shí)倒數(shù)之間的二元聯(lián)合分布規(guī)律，并以此為基礎(chǔ)設(shè)計(jì)出對(duì)應(yīng)頻率下的最可能降雨場(chǎng)景，以期為北京地區(qū)防洪降險(xiǎn)提供科學(xué)依據(jù)和推求高強(qiáng)度降雨的設(shè)計(jì)降雨提供技術(shù)支持。

1 研究區(qū)概況

十三陵水庫位于溫榆河北支東沙河，流域面積223km2，控制東沙河84%的流域面積，占溫榆河上游山區(qū)流域面積1/4。該流域?yàn)榈湫偷呐瘻貛О霛駶?rùn)大陸性季風(fēng)氣候，夏季炎熱多雨，冬季寒冷干燥，年平均氣溫10～12℃，年平均降雨量550mm，降水季節(jié)分配不均，全年降水的75%集中在夏季，7、8月常有暴雨。水庫總庫容7450萬m3，防洪設(shè)計(jì)洪水標(biāo)準(zhǔn)為100年一遇，校核洪水標(biāo)準(zhǔn)為2000年一遇。水庫上游有德勝口、錐石口、上下口及老君堂四條溝。

十三陵水庫流域設(shè)有老君堂、果莊、上口、碓臼峪、十三陵水庫5個(gè)雨量站。選取5個(gè)雨量站1980—2020年6—9月汛期每場(chǎng)降雨的雨量和歷時(shí)，采用算數(shù)平均法計(jì)算5個(gè)雨量站的平均降雨量作為十三陵水庫流域的降雨量。

2 研究方法

2.1 邊緣分布函數(shù)

邊緣分布函數(shù)的確定是構(gòu)建降雨歷時(shí)倒數(shù)和降雨量聯(lián)合分布的前提。采用應(yīng)用較為廣泛的16種邊緣分布函數(shù)擬合降雨量和歷時(shí)的倒數(shù)的邊緣分布，并運(yùn)用極大似然法估計(jì)邊緣分布的參數(shù)。16種邊緣分布函數(shù)包括Birnbaum-Saunders分布、指數(shù)分布(exponential)、極值分布(extreme value)、伽瑪分布(Gamma)、廣義極值分布(generalized extreme value)、廣義帕累托分布(generalized Pareto)、逆高斯分布(inverse Gaussian)、邏輯分布(logistic)、對(duì)數(shù)邏輯分布(log-logistic)、對(duì)數(shù)正態(tài)分布(lognormal)、Nakagami 分布、正態(tài)分布(normal)、瑞利分布(Rayleigh)、萊斯分布(Rician)、t位置尺度分布(t location scale)和韋伯分布(Weibull)，模型詳細(xì)情況參考文獻(xiàn)[7- 10]。

為了識(shí)別出最佳的邊緣分布函數(shù)，需要對(duì)擬合的結(jié)果進(jìn)行統(tǒng)計(jì)檢驗(yàn)。本文運(yùn)用赤池信息準(zhǔn)則(AIC)和貝葉斯信息準(zhǔn)則(BIC)對(duì)16種邊緣分布的擬合結(jié)果進(jìn)行優(yōu)選，識(shí)別出降雨量和歷時(shí)的倒數(shù)最適合的邊緣分布函數(shù)。分別以AIC和BIC值最小為最優(yōu)。AIC和BIC公式如下：

AIC=-2In(L)+2k

(1)

BIC=-2In(L)+In(n)×k

(2)

式中，L—極大似然函數(shù)；k—模型的變量的個(gè)數(shù)；n—數(shù)據(jù)數(shù)量。

2.2 Copula函數(shù)

Copula函數(shù)描述的是變量之間的相關(guān)性，是一種將變量聯(lián)合累積分布函數(shù)同變量邊緣累積分布函數(shù)連接起來的函數(shù)。設(shè)H為聯(lián)合分布函數(shù)，F(xiàn)和G為其邊緣分布，那么存在唯一的Copula函數(shù)C，有[11]：

H(x，y)=C[U(X)，V(Y)]

(3)

Copula函數(shù)在水文領(lǐng)域應(yīng)用廣泛。在研究中使用最多的Copula函數(shù)主要有阿基米德Copula族和橢圓Copula族兩大類[11]。為使備選的Copula更全面，本文選用25種具有代表性的Copula函數(shù)構(gòu)建歷時(shí)倒數(shù)和降雨量的聯(lián)合分布函數(shù)，并采用極大似然法估計(jì)參數(shù)θ。25種Copula函數(shù)包括Gaussian、t 2種橢圓Copula函數(shù)，Clayton、Frank、Gumbel、Independence、AMH、Joe、FGM、Plackett、Cuadras-Auge、Raftery、Shih-Louis、Linear-Spearman、Cubic、Burr、Nelsen、Galambos、Marshal-Olkin、Fischer-Hinzmann、Roch-Alegre、Fischer-Kock、BB1、BB5、Tawn 23種Archimedean Copula函數(shù)，Copula函數(shù)、參數(shù)θ取值范圍及聯(lián)合分布函數(shù)見表1。

表1 Copula函數(shù)基本特征

為選擇最優(yōu)的Copula函數(shù)來描述變量間的相關(guān)結(jié)構(gòu)，需要對(duì)各函數(shù)進(jìn)行擬合度評(píng)價(jià)和檢驗(yàn)。選用極大似然法(ML)、赤池信息準(zhǔn)則(AIC)和貝葉斯信息準(zhǔn)則(BIC)檢驗(yàn)Copula函數(shù)。以ML值最大和AIC、BIC值最小作為擬合最優(yōu)的檢驗(yàn)標(biāo)準(zhǔn)。

3 結(jié)果與分析

3.1 降雨量與歷時(shí)倒數(shù)的邊緣分布

為了確定降雨量和歷時(shí)倒數(shù)的邊緣分布，采用極大似然法估計(jì)降雨量與歷時(shí)倒數(shù)的邊緣分布參數(shù)，并利用AIC和BIC對(duì)十三陵流域降雨量與歷時(shí)倒數(shù)的16種邊緣分布進(jìn)行擬合度檢驗(yàn)。

十三陵流域的歷時(shí)倒數(shù)和降雨量16種邊緣分布函數(shù)的AIC和BIC檢驗(yàn)結(jié)果見表2—3。由表2可得，歷時(shí)倒數(shù)邊緣分布函數(shù)檢驗(yàn)中g(shù)eneralized pareto分布的AIC值和BIC值最小，由表3可得，降雨量邊緣分布函數(shù)檢驗(yàn)中l(wèi)ognormal分布的AIC值和BIC值最小，即generalized pareto分布和lognormal分布分別作為十三陵流域歷時(shí)倒數(shù)和降雨量最優(yōu)的邊緣分布。繪制優(yōu)選的擬合分布與經(jīng)驗(yàn)分布的匹配程度的比較圖(圖1a、2a)和變量邊緣分布經(jīng)驗(yàn)頻率與理論頻率關(guān)系圖(圖1b、2b)，如圖1a、2a所示，點(diǎn)據(jù)均分布在擬合線附近，同時(shí)，由圖1b、2b可得，樣本分位數(shù)點(diǎn)和相應(yīng)的理論上擬合的分位數(shù)接近。兩者都直觀得證明擬合結(jié)果較好。

表2 歷時(shí)倒數(shù)的邊緣分布AIC和BIC檢驗(yàn)結(jié)果

表3 降雨量的邊緣分布AIC和BIC檢驗(yàn)結(jié)果

圖1 歷時(shí)倒數(shù)的邊緣分布擬合圖

3.2 二元Copula分布選擇

在優(yōu)選出降雨量和歷時(shí)倒數(shù)的最優(yōu)邊緣分布后，運(yùn)用Copula函數(shù)構(gòu)建降雨量和歷時(shí)倒數(shù)的聯(lián)合分布。本文運(yùn)用極大似然法對(duì)25種Copula函數(shù)進(jìn)行參數(shù)估算，通過極大似然估計(jì)方法找出使估計(jì)的聯(lián)合概率與相關(guān)的經(jīng)驗(yàn)聯(lián)合概率最匹配的Copula參數(shù)。然后根據(jù)ML和AIC、BIC比較25種Copula函數(shù)的擬合優(yōu)度，從而確定十三陵流域的Copula函數(shù)類型。

十三陵水庫流域的降雨量和歷時(shí)倒數(shù)的檢驗(yàn)結(jié)果見表3，Roch-Alegre的ML值最大為1150.2，AIC和BIC值最小分別為-2296.4和-2289.0，即ML、AIC、BIC分析法3種評(píng)價(jià)方法均得出一致的結(jié)果，因此，采用Roch-Alegre函數(shù)描述十三陵流域降雨量和歷時(shí)倒數(shù)的聯(lián)合分布。Roch-Alegre函數(shù)參數(shù)值θ1為0.0001，為1.0973，十三陵流域降雨量與歷時(shí)倒數(shù)聯(lián)合分布函數(shù)表達(dá)式如下：H[u，v]=exp{1-[(((1-ln(u))0.0001-1)1.0973+((1-ln(v))0.0001-1)1.0973)1/1.0973+1]1/0.0001}

圖2 降雨量的邊緣分布擬合圖

表3 Copula函數(shù)檢驗(yàn)結(jié)果

3.3 不同設(shè)計(jì)頻率下的降雨量和降雨歷時(shí)分析

本文從二元Copula族中選擇最佳的函數(shù)Roch-Alegre來建立降雨量與降雨歷時(shí)倒數(shù)的依賴結(jié)構(gòu)，將只考慮降雨量的單變量的設(shè)計(jì)頻率概念擴(kuò)展到同時(shí)考慮降雨量與降雨歷時(shí)的多變量。這種方法對(duì)于設(shè)計(jì)降雨量大而降雨歷時(shí)短這種極其危險(xiǎn)的降雨具有重要的現(xiàn)實(shí)意義。圖3為Roch-Alegre函數(shù)的聯(lián)合分布圖，圖中的等高線對(duì)應(yīng)著相應(yīng)設(shè)計(jì)頻率的降雨，同一個(gè)等高線上的點(diǎn)構(gòu)成降雨量大而降雨歷時(shí)短二元風(fēng)險(xiǎn)場(chǎng)景。比如，聯(lián)合分布Copula值為0.9的等高線對(duì)應(yīng)著無數(shù)個(gè)設(shè)計(jì)頻率為0.1的降雨量與降雨歷時(shí)組合。降雨量與降雨歷時(shí)的無限組合會(huì)在統(tǒng)計(jì)上是相似的，而它們的影響可能截然不同。如圖3中所示的十三陵水庫的降雨，同一設(shè)計(jì)頻率降雨下的降雨量與降雨歷時(shí)具有無限種組合，在曲線的右側(cè)對(duì)應(yīng)著降雨歷時(shí)短，降雨量小的場(chǎng)景，而在曲線的上側(cè)對(duì)應(yīng)著降雨歷時(shí)短，降雨量大的場(chǎng)景。選擇哪一種組合是亟需解決的問題。從無限的降雨量與降雨歷時(shí)組合中進(jìn)行選擇的一種直觀方法是根據(jù)它們相關(guān)的Copula密度值為它們分配權(quán)重。Roch-Alegre函數(shù)的概率密度值可以衡量這無限種組合的權(quán)重，因此本文采用密度函數(shù)極大值處的組合為設(shè)計(jì)頻率下的降雨量與降雨歷時(shí)組合，即選擇一種最可能的事件。

圖3 十三陵水庫流域降雨量和降雨歷時(shí)倒數(shù)的Roch-Alegre函數(shù)聯(lián)合分布圖

通過對(duì)Roch-Alegre密度函數(shù)求導(dǎo)找到極值點(diǎn)，即設(shè)計(jì)頻率下最可能降雨量與降雨歷時(shí)組合場(chǎng)景。設(shè)計(jì)頻率0.1條件下的最佳組合如圖4所示，降雨量為57mm，降雨歷時(shí)倒數(shù)為1.3(1/h)。而不考慮降雨量與降雨歷時(shí)依賴關(guān)系各自得到設(shè)計(jì)頻率為0.1條件下的降雨量為29mm，降雨歷時(shí)倒數(shù)為1.2mm。可以看出十三陵水庫控制流域容易出現(xiàn)短歷時(shí)大降雨量，即降雨強(qiáng)度大的情況。而這種情況會(huì)給防洪帶來很大的風(fēng)險(xiǎn)，所以在設(shè)計(jì)降雨時(shí)就需要提前考慮這種風(fēng)險(xiǎn)疊加。

圖4 設(shè)計(jì)頻率0.1條件下降雨量和歷時(shí)倒數(shù)的最佳組合

綜上所述，運(yùn)用本文提出的考慮降雨量與歷時(shí)依賴關(guān)系的方法可以選擇具有聯(lián)合頻率的危險(xiǎn)組合事件，可以用于危險(xiǎn)降雨的設(shè)計(jì)與評(píng)估，在十三陵的防洪中具有重要的現(xiàn)實(shí)意義。

4 結(jié)論

(1)通過擬合優(yōu)度評(píng)價(jià)，得出北京市十三陵水庫流域降雨量采用lognormal分布擬合效果較好，而降雨歷時(shí)的倒數(shù)采用generalized pareto分布擬合效果較好，建議選擇Roch-Alegre函數(shù)描述降雨量和歷時(shí)倒數(shù)的聯(lián)合分布。

(2)采用雙變量設(shè)計(jì)降雨模型與分開考慮降雨量與歷時(shí)的傳統(tǒng)方法相比，得到的設(shè)計(jì)降雨量偏大，降雨歷時(shí)偏小。

(3)本文所提的考慮降雨量與歷時(shí)依賴關(guān)系設(shè)計(jì)降雨方法，有助于應(yīng)對(duì)降雨量大而降雨歷時(shí)短的二元風(fēng)險(xiǎn)，在水庫防洪中具有現(xiàn)實(shí)意義，以期為推求高強(qiáng)度降雨的設(shè)計(jì)降雨提供技術(shù)支持。本文未探討短歷時(shí)強(qiáng)降雨所形成的洪水過程，有待繼續(xù)研究。