基于非對(duì)稱極值Copula的設(shè)計(jì)暴雨過(guò)程線分析

陳子燊，趙玲玲,楊興

(1. 中山大學(xué)新華學(xué)院,廣東廣州510520;2.中山大學(xué)地理科學(xué)與規(guī)劃學(xué)院, 廣東廣州510275； 3. 廣州地理研究所，廣東廣州510070;4.安徽理工大學(xué)地球與環(huán)境學(xué)院，安徽淮南232000)

臺(tái)風(fēng)災(zāi)害是全球發(fā)生頻率最高的一種自然災(zāi)害。臺(tái)風(fēng)帶來(lái)的極端降雨過(guò)程可導(dǎo)致山洪暴發(fā)和衍生滑坡等地質(zhì)災(zāi)害，從而造成重大的生命財(cái)產(chǎn)和經(jīng)濟(jì)損失。如何應(yīng)對(duì)極端暴雨是城市和山地防災(zāi)減災(zāi)都需面對(duì)的重大問(wèn)題。一些研究人員從城市應(yīng)對(duì)極端天氣事件與防災(zāi)減災(zāi)的風(fēng)險(xiǎn)管理角度對(duì)雨型作了探索。蔣明[1]指出，雨型是描述降雨過(guò)程和降雨強(qiáng)度在時(shí)間尺度上的分配過(guò)程, 是徑流過(guò)程計(jì)算的基礎(chǔ)。成丹等[2]把設(shè)計(jì)雨型作為制定排水防澇系統(tǒng)設(shè)計(jì)時(shí)的重要因素，應(yīng)用于城市市政排水系統(tǒng)的規(guī)劃和管理及排水分析，為城市流域雨洪調(diào)度計(jì)算提供科學(xué)依據(jù)。葉姍姍等[3]選取宿遷市實(shí)測(cè)的主副型雨峰偏后的暴雨雨型，對(duì)其降雨過(guò)程進(jìn)行同頻率分時(shí)段縮放,采用Copula 函數(shù)的風(fēng)險(xiǎn)聯(lián)合概率模型分析了不同的兩時(shí)段之間出現(xiàn)的暴雨風(fēng)險(xiǎn)。楊星等[4]利用深圳雨量站34 a實(shí)測(cè)逐時(shí)降雨資料，對(duì)比了不同典型暴雨設(shè)計(jì)雨型研究方面的差異，按構(gòu)建的兩變量Copula推求了深圳市不同重現(xiàn)期雨型的風(fēng)險(xiǎn)率和典型暴雨的特征。

山區(qū)中小流域山洪至今仍然是防災(zāi)減災(zāi)的重要研究方向，為此，可借鑒設(shè)計(jì)洪水過(guò)程線的方法，從高維(大于二維)尺度上設(shè)計(jì)典型暴雨過(guò)程，將更有利于山洪風(fēng)險(xiǎn)管理。至今在應(yīng)用Copula函數(shù)分析三變量洪水的聯(lián)合概率分布和設(shè)計(jì)洪水過(guò)程線已有不少研究。侯蕓蕓等[5]和ZHANG等[6]分別應(yīng)用對(duì)稱的單參數(shù)Archimedean Copula函數(shù)分析了洪水三變量的聯(lián)合概率分布和條件概率分布。由于具有不同相關(guān)性的高維隨機(jī)變量，單參數(shù)難以真實(shí)反映其復(fù)雜的不對(duì)稱相關(guān)結(jié)構(gòu)。非對(duì)稱形式的Copula 函數(shù)具有更加靈活的參數(shù)和結(jié)構(gòu)形式，更適合用于擬合高維的隨機(jī)變量[7 ]。為此，Grimaldi等[8]、Ganguli等[9]、陳子燊等[7]分別采用非對(duì)稱的阿基米德Copulas(Asymmetric Archimedean Copulas)構(gòu)建了不對(duì)稱三變量洪水要素聯(lián)合分布模型推算設(shè)計(jì)洪水，以嘗試應(yīng)用于洪水風(fēng)險(xiǎn)規(guī)劃管理。肖義等[10]和李天元等[11]則分別采用兩變量和三變量的Copula函數(shù)建立了聯(lián)合分布的設(shè)計(jì)洪水過(guò)程線的推求方法，為設(shè)計(jì)洪水過(guò)程線提供了一種新思路。本文把非對(duì)稱阿基米德極值Copula用于構(gòu)建山區(qū)中小流域設(shè)計(jì)暴雨過(guò)程線，希望有助于防災(zāi)減災(zāi)的風(fēng)險(xiǎn)管理。

1 三變量Copula函數(shù)

1.1 Copula函數(shù)的定義

設(shè)隨機(jī)變量Xi(i=1,2,…,n)的邊緣分布函數(shù)分別為FXi(xi)=P(Xi≤xi)，其中n為隨機(jī)變量的個(gè)數(shù)，xi為隨機(jī)變量Xi的值。依Sklar理論，對(duì)于連續(xù)分布函數(shù)FXi(xi)，存在唯一的聯(lián)合分布函數(shù)[6]：

H(x1,x2,…,xn)=C(FX1(x1),FX2(x2),…,FXn(xn))=C(u1,u2,…,un)

(1)

利用Copula函數(shù)構(gòu)造聯(lián)合概率分布，使得變量的所有信息都存在于邊緣分布函數(shù)里，不會(huì)在轉(zhuǎn)換過(guò)程中產(chǎn)生信息失真。因此，Copula函數(shù)理論是構(gòu)建多變量水文聯(lián)合概率分布的很好的工具[12]。

1.2 三變量Archimedean Copula

三變量對(duì)稱的Archimedean Copula單參數(shù)形式[13]為：

(2)

式中uj∈[0,1](j>1)——邊緣分布;φθ——Archimedean Copula生成元，θ——參數(shù)。

φθ(u3)]

(3)

三變量非對(duì)稱Archimedean Copula的形式為[14]：

C(u1,u2,u3)=C1(u3,C2(u1,u2))=

(4)

式中 符號(hào)“0”——函數(shù)組合。

常用的三維非對(duì)稱Archimedean族Copula函數(shù)如下。

a) M3 (Frank) Copula

(1-e-θ1u3))},θ2>θ1∈[0,∞)

(5)

b) M4 (Clayton) Copula

(6)

c) M5 (Arch2) Copula

(7)

d) M6(Gumbel-Hougaard) Copula

C(u1,u2,u3)=exp{-([(-lnu1)θ2+(-lnu2)θ2]θ1/θ2+(-lnu3)θ1)1/θ1},θ2>θ1∈[1,∞)

(8)

e) M12 (Arch12) Copula

(9)

1.3 三變量聯(lián)合重現(xiàn)期和條件重現(xiàn)期

以算符“∨”定義“或”三維極端事件中至少有一個(gè)被超過(guò)情況下的“或”聯(lián)合重現(xiàn)期為：

(10)

以“∧”定義“且”三維極端事件同時(shí)被超過(guò)情況下的“且”聯(lián)合重現(xiàn)期為:

(11)

2個(gè)不超過(guò)事件發(fā)生下的條件概率為[5]：

(12)

則事件(X1>x1|X2≤x2,X3≤x3)下的條件重現(xiàn)期為:

T(x1|X2≤x2,X3≤x3)=

(13)

2個(gè)等量事件發(fā)生下的條件概率為：

F(x1|X2=x2,X3=x3)=

(14)

則事件(X1>x1|X2=x2,X3=x3)下的條件重現(xiàn)期為:

T(x1|X2=x2,X3=x3)=

(15)

一個(gè)等量事件發(fā)生下的條件概率為:

(16)

則事件至少有一個(gè)為超過(guò)事件下的條件重現(xiàn)期為:

(17)

一個(gè)不超過(guò)事件發(fā)生下的條件概率為:

(18)

則事件至少有一個(gè)為超過(guò)事件下的條件重現(xiàn)期為:

(19)

2 實(shí)例研究

2.1 流域水文氣象背景與基本數(shù)據(jù)

選取廣東曹江流域?yàn)閷?shí)例研究。曹江是廣東獨(dú)流入海鑒江的一級(jí)支流，發(fā)源于高州馬貴鎮(zhèn)山心村的藍(lán)蓬嶺。出口斷面大拜水文站集水面積394 km2，屬于典型的中小流域。曹江流域多年平均年雨量約2 160 mm，最大年雨量可達(dá)3 150 mm，是廣東的臺(tái)風(fēng)暴雨高區(qū)之一。1967年11月7日代號(hào)為6720 的“Emma”超強(qiáng)臺(tái)風(fēng)，在流域西側(cè)的湛江市登陸，最大風(fēng)速65 m/s，中心氣壓912 hPa。2013年8月14日代號(hào)為“Utor”的超強(qiáng)臺(tái)風(fēng)，在流域東側(cè)的陽(yáng)江市登陸，最大風(fēng)速60 m/s，中心氣壓925 hPa。大拜雨量站測(cè)得二者最大24 h暴雨分別為419.9、412.1 mm，均達(dá)到特大暴雨級(jí)別，也是1967—2013年2個(gè)最大的24 h小時(shí)雨量。

本文根據(jù)曹江流域出口斷面大拜雨量站1967—2013年逐時(shí)降水記錄數(shù)據(jù)，首先提取歷年最大24 h雨量(R24)，再分別提取最大1 h雨量(雨峰：R1)和連續(xù)最大6 h雨量(R6)數(shù)據(jù)，由R1、R6和R24作為實(shí)例分析的樣本，分別構(gòu)建1967、2013年這3個(gè)歷時(shí)雨量聯(lián)合分布的2場(chǎng)設(shè)計(jì)暴雨過(guò)程線。

2.2 邊緣分布與聯(lián)合分布

分別采用皮爾遜三型分布(PE3)和廣義極值分布(GEV)對(duì)R1、R6和R24樣本加以擬合。參數(shù)估計(jì)使用線性矩(L-矩)方法。經(jīng)驗(yàn)頻率分布使用Gringorten公式計(jì)算。擬合結(jié)果采用均方根誤差(RMSE)和概率點(diǎn)據(jù)相關(guān)系數(shù)(PPCC)檢驗(yàn)其擬合優(yōu)度。表1擇優(yōu)對(duì)比結(jié)果表明，R1、R6和R24都以GEV分布相對(duì)更優(yōu)。廣義極值(GEV)分布函數(shù)為：

FX(x)=P(X

(20)

式中ξ、β、μ——形態(tài)參數(shù)、尺度參數(shù)和位置參數(shù)。

表1 曹江大拜站三變量暴雨樣本的 邊緣分布參數(shù)與優(yōu)度檢驗(yàn)值

計(jì)算的R1、R6和R24兩兩間的Kendall 相關(guān)系數(shù)表明大拜雨量站不同歷時(shí)暴雨間都存在正相關(guān)性，其中R6和R24最強(qiáng)，τ= 0. 484;R1和R6相關(guān)性次之，τ=0.399;R1和R24相關(guān)性較弱，τ= 0.171。采用Kendall相關(guān)系數(shù)τ與Copula參數(shù)θ的關(guān)系式，構(gòu)建5種非對(duì)稱Archimedean Copula[6]。為了對(duì)比，分別將對(duì)應(yīng)的5種對(duì)稱三變量Archimedean Copula通過(guò)MLM法計(jì)算其參數(shù)θ。采用Akaike 信息準(zhǔn)則(AIC)、最小(OLS)準(zhǔn)則和Genest-Rivest圖形法驗(yàn)證理論聯(lián)合分布函數(shù)與經(jīng)驗(yàn)聯(lián)合分布函數(shù)的擬合程度，結(jié)果見(jiàn)表2和圖1?？梢?jiàn)以二維Gumbel-Hougaard 為基Copula的三維非對(duì)稱形式的M6 Copula 的OLS和AIC值最小，擬合度最高，各點(diǎn)均勻地分布在45°線左右的非對(duì)稱Archimedean M6 Copula具有相對(duì)最優(yōu)的擬合度。Nelson[13]和Salvadori等[14]證明當(dāng)且僅當(dāng)邊緣分布和Copula函數(shù)均為極值分布時(shí)，構(gòu)造的聯(lián)合分布才是極值分布，而Gumbel-Hougaard Copula是Archimedean Copula函數(shù)族中的唯一多變量極值Copula函數(shù)，適用于極端事件的頻率分析?？紤]到R1、R6和R24之間的相關(guān)性存在明顯差別，因此，選用非對(duì)稱Archimedean M6 Copula構(gòu)建大拜雨量站歷年最大24 h暴雨量不同歷時(shí)R1、R6和R24之間的三維聯(lián)合分布:

C(u1,u2,u3)=exp{-([(-lnu1)1.937+

(-lnu2)1.937]1.395/1.937+(-lnu3)1.395)1/1.395}

(21)

表2 非對(duì)稱和對(duì)稱三維Copula 參數(shù)估計(jì)及擬合優(yōu)度對(duì)比評(píng)價(jià)

a) M2 b) M4 c) M5 d) M6 e)M12圖1 5種非對(duì)稱Archimedean Copula的概率分布擬合

2.3 典型暴雨過(guò)程線選擇

典型暴雨的特征，包括降雨集中程度、雨峰位置和雨峰大小等。典型暴雨過(guò)程線的選擇采用以下原則: ①選擇雨峰量大具有一定代表性的實(shí)測(cè)暴雨過(guò)程線；②從防洪安全考慮，對(duì)主峰靠后和主峰靠前的2種雨型的風(fēng)險(xiǎn)概率加以對(duì)比;③設(shè)計(jì)暴雨過(guò)程線采用同頻率放大法，以降水主峰對(duì)流域洪水形成為首要影響因子，選定時(shí)段為1 h的設(shè)計(jì)雨峰為設(shè)計(jì)標(biāo)準(zhǔn)，使得放大的過(guò)程線形狀能與原來(lái)的典型過(guò)程一致。

按照短歷時(shí)強(qiáng)降水強(qiáng)度20 mm/h劃分雨峰，根據(jù)典型年暴雨過(guò)程的雨峰位置，選取1967、2013年作為曹江流域設(shè)計(jì)暴雨的典型年,24 h暴雨過(guò)程線見(jiàn)圖2。圖2顯示，1967年最大24 h暴雨過(guò)程為主副多峰雨型，主峰靠后，2013年最大24 h暴雨過(guò)程為單峰雨型，主峰靠前。這2個(gè)典型暴雨R1、R6、R243個(gè)時(shí)段最大降水量與相應(yīng)的重現(xiàn)期見(jiàn)表3。

a) 1967年

b) 2013年圖2 2個(gè)典型年的最大24 h暴雨過(guò)程

>

表3 2個(gè)典型年R1、R6、R24最大降水量及重現(xiàn)期

典型年R1R1/mm重現(xiàn)期/aR6R6/mm重現(xiàn)期/aR24R24/mm重現(xiàn)期/aR1-R6-R24“或”重現(xiàn)期/a“且”重現(xiàn)期/a1967年40.83.2168.612.3419.9132.83.1111.32013年74.533.5285.0103.8412.1115.728.2331.7

從表3可見(jiàn)，1967、2013年的R1-R6-R24組合雨量的“或”聯(lián)合重現(xiàn)期小于單一時(shí)段雨量重現(xiàn)期，此說(shuō)明考慮多時(shí)段組合條件下某一時(shí)段雨量致災(zāi)的可能性最高，相比較同時(shí)出現(xiàn)三時(shí)段組合雨量的“且”聯(lián)合重現(xiàn)期可能性很小。表4為不同時(shí)段雨量組合的“或”聯(lián)合重現(xiàn)期，可見(jiàn)同頻率下R1-R6-R24三時(shí)段組合雨量的“或”聯(lián)合重現(xiàn)期最小，危險(xiǎn)率最高。因此，如果以三時(shí)段雨量組合的“或”聯(lián)合重現(xiàn)期作為流域的防雨洪標(biāo)準(zhǔn)，由此設(shè)計(jì)的暴雨過(guò)程線對(duì)于應(yīng)對(duì)流域雨洪風(fēng)險(xiǎn)更合適。

表4 不同時(shí)段雨量組合的“或”聯(lián)合重現(xiàn)期

2.4 同頻率法推求設(shè)計(jì)暴雨過(guò)程線

有關(guān)研究指出[11]，由于對(duì)任一給定的三變量重現(xiàn)期Tu1,u2,u3，理論上存在無(wú)數(shù)種u1、u2、u3的組合滿足式(10)，如果按照同頻率放大法的思路，假定R1、R6、R243個(gè)時(shí)段雨量同頻率，即令u1=u2=u3，可得到基于某一聯(lián)合重現(xiàn)期Tu1,u2,u3的頻率組合(u1,u2,u3)。根據(jù)此組合，按照各變量的邊緣分布函數(shù)反推可得到3個(gè)不同時(shí)段雨量的聯(lián)合設(shè)計(jì)值組合(r1、r6、r24)，進(jìn)而以此設(shè)計(jì)值組合放大典型暴雨過(guò)程，即得到基于三變量聯(lián)合分布的設(shè)計(jì)暴雨過(guò)程線。采用非對(duì)稱M6函數(shù)推算R1、R6、R243個(gè)時(shí)段雨量同頻率分布聯(lián)合設(shè)計(jì)值公式如下：

u1=u2=u3=[1-(1/Tu1,u2,u3)]α;

(22)

按相同原理，可分別推算兩變量u1、u2的重現(xiàn)期Tu1、u2，u1、u3的重現(xiàn)期Tu1、u3和u2、u3的重現(xiàn)期Tu2,u3的同頻率分布聯(lián)合設(shè)計(jì)值。

從表5多變量同頻率設(shè)計(jì)值計(jì)算結(jié)果可見(jiàn)，R1-R6-R24組合同頻率設(shè)計(jì)暴雨設(shè)計(jì)值明顯大于其它同一重現(xiàn)水平組合和單一時(shí)段暴雨的設(shè)計(jì)值。由于多變量方法是基于多個(gè)時(shí)段組合的聯(lián)合重現(xiàn)期，考慮了變量之間的相關(guān)性，設(shè)計(jì)值會(huì)大于單變量同頻率設(shè)計(jì)值。有關(guān)研究結(jié)果顯示[7,15]，三變量同頻率設(shè)計(jì)值十分接近于按聯(lián)合概率密度最大值推算的三變量“或”重現(xiàn)期設(shè)計(jì)值。作為工程設(shè)計(jì)與風(fēng)險(xiǎn)管理，盡管存在偏向安全問(wèn)題，但采用R1-R6-R24組合同頻率設(shè)計(jì)暴雨值為更高安全標(biāo)準(zhǔn)的防雨洪工程設(shè)計(jì)或風(fēng)險(xiǎn)預(yù)警提供了科學(xué)依據(jù)。

表5 樣本設(shè)計(jì)值

選取1967、2013年的受臺(tái)風(fēng)影響的2場(chǎng)典型暴雨過(guò)程進(jìn)行同頻率分時(shí)段縮放。放大系數(shù)公式:

K=X設(shè)計(jì)/X典型

(23)

式中X設(shè)計(jì)——不同重現(xiàn)期的設(shè)計(jì)降雨量；X典型——典型暴雨降雨量。

以雨峰同頻率放大法求重現(xiàn)期200 a(P=0.05%)R1-R6-R24三時(shí)段雨量聯(lián)合分布的設(shè)計(jì)暴雨過(guò)程線。為了比較，另推求了R1-R6和R1-R24兩變量聯(lián)合分布以及以雨峰同頻率放大的設(shè)計(jì)暴雨過(guò)程線。由圖3可見(jiàn)，多變量方法與單變量方法所推求的200年一遇設(shè)計(jì)暴雨過(guò)程線的比較顯示，采用R1-R6-R24組合方法推求的3個(gè)時(shí)段雨量的設(shè)計(jì)值均大于相應(yīng)單一時(shí)段樣本推算的設(shè)計(jì)值，也大于采用2個(gè)時(shí)段雨量組合的設(shè)計(jì)值。可見(jiàn)，采用R1-R6-R24組合方法放大的過(guò)程線對(duì)流域防雨洪設(shè)計(jì)工程更安全，采用R1-R6-R24組合的設(shè)計(jì)暴雨過(guò)程線也更加符合流域水文現(xiàn)象的內(nèi)在規(guī)律和防洪工程實(shí)際的要求。

a) 1967年

b) 2013年圖3 200年一遇設(shè)計(jì)暴雨過(guò)程線比較

2.5 設(shè)計(jì)暴雨過(guò)程線的條件重現(xiàn)期及危險(xiǎn)率

按式(12)—(19)分別推算了1967、2013年典型暴雨過(guò)程的條件重現(xiàn)期及相應(yīng)的條件概率，結(jié)果見(jiàn)表6。結(jié)果顯示：①等量事件發(fā)生條件下的條件重現(xiàn)期小于不超過(guò)事件發(fā)生下的條件重現(xiàn)期，出現(xiàn)的危險(xiǎn)率P(超值條件概率)則大之；②1967年典型暴雨出現(xiàn)的4個(gè)條件重現(xiàn)期都分別小于2013年主峰靠前的單峰雨型的對(duì)應(yīng)條件重現(xiàn)期，危險(xiǎn)率則反之。其中，等于24 h雨量(419.9 mm) 一個(gè)等量事件條件下出現(xiàn)的重現(xiàn)期最小，危險(xiǎn)率最大。這表明主副多峰雨型且主峰靠后1967年暴雨過(guò)程對(duì)于流域防洪安全具有更大的威脅。

表6 2個(gè)典型設(shè)計(jì)暴雨過(guò)程的條件重現(xiàn)期與危險(xiǎn)率

因此，對(duì)于主副多峰雨型且主峰靠后的暴雨過(guò)程，由于前期降雨首先使得流域下墊面土壤水趨于飽和產(chǎn)生超滲產(chǎn)流，疊加在后期雨峰形成的坡面流將匯集形成更強(qiáng)的洪水過(guò)程，流域出現(xiàn)洪水的風(fēng)險(xiǎn)更大，是此流域防范雨洪風(fēng)險(xiǎn)的最主要類型。

3 結(jié)論

本文將不同歷時(shí)雨量之間具有相關(guān)關(guān)系的暴雨過(guò)程簡(jiǎn)化為雨峰、6 h雨量和24 h雨量三變量聯(lián)合分布，采用非對(duì)稱極值Copula構(gòu)建曹江流域典型暴雨過(guò)程線，并與由2個(gè)時(shí)段和由單一雨峰的同頻率設(shè)計(jì)暴雨過(guò)程線方法進(jìn)行了比較。研究結(jié)果有以下結(jié)論。

a) 采用3個(gè)歷時(shí)雨量推求的曹江流域設(shè)計(jì)暴雨值大于2個(gè)時(shí)段和單一時(shí)段設(shè)計(jì)暴雨值，由此放大的設(shè)計(jì)暴雨過(guò)程線，整體效果相對(duì)最優(yōu)，對(duì)設(shè)計(jì)雨型的研究方法提供了新思路。由得到的典型設(shè)計(jì)暴雨過(guò)程線推算的流域洪水過(guò)程更符合流域水文現(xiàn)象的內(nèi)在規(guī)律和防洪工程的實(shí)際要求。

b) 按24 h最大雨量選取的2013年的單峰雨型與1967年的主副多峰雨型都具有較高代表性。但與2013年主峰靠前雨型比較，主峰靠后的1967年的暴雨過(guò)程危險(xiǎn)率更大，對(duì)流域防洪安全具有更大的威脅，構(gòu)建的典型設(shè)計(jì)暴雨過(guò)程線更具代表性。

c) 1967、2013年2個(gè)典型年的R1-R6-R243個(gè)時(shí)段雨量聯(lián)合分布的“或”聯(lián)合重現(xiàn)期都小于單一時(shí)段雨量重現(xiàn)期，危險(xiǎn)率最大，以多時(shí)段雨量組合的“或”聯(lián)合重現(xiàn)期作為流域的設(shè)計(jì)標(biāo)準(zhǔn)，對(duì)于應(yīng)對(duì)此流域雨洪風(fēng)險(xiǎn)更合適。