陳子燊,趙玲玲,楊興
(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ì)稱(chēng)的單參數(shù)Archimedean Copula函數(shù)分析了洪水三變量的聯(lián)合概率分布和條件概率分布。由于具有不同相關(guān)性的高維隨機(jī)變量,單參數(shù)難以真實(shí)反映其復(fù)雜的不對(duì)稱(chēng)相關(guān)結(jié)構(gòu)。非對(duì)稱(chēng)形式的Copula 函數(shù)具有更加靈活的參數(shù)和結(jié)構(gòu)形式,更適合用于擬合高維的隨機(jī)變量[7 ]。為此,Grimaldi等[8]、Ganguli等[9]、陳子燊等[7]分別采用非對(duì)稱(chēng)的阿基米德Copulas(Asymmetric Archimedean Copulas)構(gòu)建了不對(duì)稱(chēng)三變量洪水要素聯(lián)合分布模型推算設(shè)計(jì)洪水,以嘗試應(yīng)用于洪水風(fēng)險(xiǎn)規(guī)劃管理。肖義等[10]和李天元等[11]則分別采用兩變量和三變量的Copula函數(shù)建立了聯(lián)合分布的設(shè)計(jì)洪水過(guò)程線的推求方法,為設(shè)計(jì)洪水過(guò)程線提供了一種新思路。本文把非對(duì)稱(chēng)阿基米德極值Copula用于構(gòu)建山區(qū)中小流域設(shè)計(jì)暴雨過(guò)程線,希望有助于防災(zāi)減災(zāi)的風(fēng)險(xiǎn)管理。
設(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]。
三變量對(duì)稱(chēng)的Archimedean Copula單參數(shù)形式[13]為:
(2)
式中uj∈[0,1](j>1)——邊緣分布;φθ——Archimedean Copula生成元,θ——參數(shù)。
φθ(u3)]
(3)
三變量非對(duì)稱(chēng)Archimedean Copula的形式為[14]:
C(u1,u2,u3)=C1(u3,C2(u1,u2))=
(4)
式中 符號(hào)“0”——函數(shù)組合。
常用的三維非對(duì)稱(chēng)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)
以算符“∨”定義“或”三維極端事件中至少有一個(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)
選取廣東曹江流域?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ò)程線。
分別采用皮爾遜三型分布(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ì)稱(chēng)Archimedean Copula[6]。為了對(duì)比,分別將對(duì)應(yīng)的5種對(duì)稱(chēng)三變量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ì)稱(chēng)形式的M6 Copula 的OLS和AIC值最小,擬合度最高,各點(diǎn)均勻地分布在45°線左右的非對(duì)稱(chēng)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ì)稱(chēng)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ì)稱(chēng)和對(duì)稱(chēng)三維Copula 參數(shù)估計(jì)及擬合優(yōu)度對(duì)比評(píng)價(jià) a) M2 b) M4 c) M5 d) M6 e)M12圖1 5種非對(duì)稱(chēng)Archimedean Copula的概率分布擬合 典型暴雨的特征,包括降雨集中程度、雨峰位置和雨峰大小等。典型暴雨過(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)期 有關(guān)研究指出[11],由于對(duì)任一給定的三變量重現(xiàn)期Tu1,u2,u3,理論上存在無(wú)數(shù)種u1、u2、u3的組合滿(mǎn)足式(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ì)稱(chēng)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ò)程線比較 按式(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)的最主要類(lèi)型。 本文將不同歷時(shí)雨量之間具有相關(guān)關(guān)系的暴雨過(guò)程簡(jiǎn)化為雨峰、6 h雨量和24 h雨量三變量聯(lián)合分布,采用非對(duì)稱(chēng)極值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)更合適。2.3 典型暴雨過(guò)程線選擇
2.4 同頻率法推求設(shè)計(jì)暴雨過(guò)程線
2.5 設(shè)計(jì)暴雨過(guò)程線的條件重現(xiàn)期及危險(xiǎn)率
3 結(jié)論