3維Copula函數(shù)在降雨特征多變量頻率分析中的應(yīng)用

劉成林，周玉文，隋 軍，高 琳

(1.北京工業(yè)大學(xué)建筑工程學(xué)院，100124北京;2.廣州市市政工程設(shè)計研究院，510060廣州;3.北京市水質(zhì)科學(xué)與水環(huán)境恢復(fù)重點(diǎn)實(shí)驗(yàn)室，100124北京)

我國大部分城市的洪水主要由暴雨形成，與通過流量資料估算設(shè)計洪水相比，根據(jù)暴雨資料推求設(shè)計洪水是一種間接方法［1］，但在實(shí)際工作中經(jīng)常遇到所在地點(diǎn)流量系列太短無法直接推求設(shè)計洪水的現(xiàn)象.降雨監(jiān)測較為系統(tǒng)，數(shù)據(jù)序列長且完整，因此，通?；诮涤觊_展水文分析，進(jìn)而確定城市洪澇系統(tǒng)設(shè)施規(guī)模.一般地，一場降雨通?？捎糜炅?、歷時和雨強(qiáng)等特征量以及雨量在時空上的變化來反映，設(shè)計降雨是地表徑流計算的重要參數(shù)，主要包括以下特征要素:頻率(重現(xiàn)期)、雨量(設(shè)計降雨量、次降雨量、降雨峰值)、歷時(設(shè)計歷時、總歷時、峰值位置).目前，暴雨強(qiáng)度公式是普遍采用的設(shè)計降雨計算方法，是基于降雨量的單變量分析，不具備計算影響洪澇設(shè)施過流和調(diào)蓄能力的“降雨峰值”和“總降雨量”等特征量的功能.近幾十年來，人們已經(jīng)認(rèn)識到采用單變量極值分布進(jìn)行頻率分析研究暴雨事件存在一定的局限性，開始探討利用兩變量聯(lián)合分布描述暴雨不同特征值之間的相互關(guān)系，以更加全面地描述整個暴雨事件［2］，但用3變量及以上的聯(lián)合分布描述暴雨事件的研究還鮮有報道.

Copula［3-4］是求解多變量概率問題優(yōu)良的數(shù)學(xué)工具，該函數(shù)可將多個隨機(jī)變量的邊際分布連接起來得到其聯(lián)合分布，是研究與變量尺度無關(guān)的相關(guān)性度量的一種途徑，同時可以基于給定的邊際分布構(gòu)造聯(lián)合分布.近幾年，Copula函數(shù)理論廣泛應(yīng)用于降雨［2，5-7］、洪水［8-9］、干旱［10］和枯水［11］的多特征屬性頻率分析等問題.為此，以降雨特征量(設(shè)計降雨量、次降雨量和降雨峰值)頻率研究為背景，引入3維Copula函數(shù)構(gòu)建降雨特征變量的聯(lián)合概率分布模型，研究不同量級降雨特征變量的遭遇概率及條件概率，以期為城市防洪排澇提供科學(xué)參考.

1 研究方法

1.1 Copula 函數(shù)簡介

Copula函數(shù)是用來描述變量間相依關(guān)系的函數(shù)，不限定變量的邊際分布.通過Copula模型可以將k個任意形式的邊際分布連接起來，生成一個多變量聯(lián)合概率分布模型.Copula函數(shù)的原理可參考文獻(xiàn)［3-4，8］.Copula函數(shù)主要有3種類型，即橢圓型、阿基米德型和二次型.目前，在水文領(lǐng)域應(yīng)用最廣泛的是Archimedena Copula函數(shù)，其常見的3維Copula函數(shù)主要有 Clayton copula、Gumbel-Hougaard(GH)copula、Ali-Mikhail-Haq(AMH)copula和Frank copula(見表1).

表1 3維Copula函數(shù)

1.2 參數(shù)估計及擬合檢驗(yàn)

Copula函數(shù)不限定變量的邊際分布，確定適宜的邊際分布是應(yīng)用Copula函數(shù)的第一步.選用3種在水文頻率分析中常用的分布線型(P-Ⅲ型分布曲線(P3)、對數(shù)正態(tài)分布曲線(LN2)和廣義極值分布曲線(GEV))分別對“設(shè)計降雨量”、“降雨峰值”和“次降雨量”進(jìn)行曲線擬合，采用概率點(diǎn)據(jù)相關(guān)系數(shù)檢驗(yàn)法(PPCC)、擬優(yōu)平方和準(zhǔn)則法(RMSE)和擬優(yōu)絕對值準(zhǔn)則法(MAE)3種擬合優(yōu)度檢驗(yàn)方法確定出與各變量數(shù)據(jù)系列擬合效果最好的邊際分布線型［12］.

Copula函數(shù)的參數(shù)估計主要有非參數(shù)法和參數(shù)法，其中非參數(shù)法要求有明確的Kendall'sτ與Copula參數(shù)θ的表達(dá)式，適用于單參數(shù)2維Copula函數(shù)的參數(shù)確定.而參數(shù)法相對比較靈活，本文采用參數(shù)法中的Inference of Functions for Margins(IFM)［13］方法確定 3維 Copula函數(shù)的參數(shù).

IFM方法是將邊際分布的參數(shù)與Copula函數(shù)的參數(shù)分別進(jìn)行估計，該過程由以下步驟完成.

采用極大似然法(ML方法)［14］估計邊際分布中的參數(shù)，即

采用ML方法估計Copula中的參數(shù)

Copula函數(shù)在實(shí)際應(yīng)用中的一個主要問題就是函數(shù)形式的選擇.Embrechts等［15］對不同的Copula函數(shù)模型進(jìn)行了比較，發(fā)現(xiàn)采用不同形式的Copula函數(shù)可能導(dǎo)致不同的分析結(jié)果，因此，選擇合適的Copula函數(shù)尤為重要.檢驗(yàn)與評價Copula函數(shù)的方法較多，本文采用AIC信息準(zhǔn)則法(AIC)及離差平方和最小(OLS)準(zhǔn)則法［9］對Copula函數(shù)的擬合優(yōu)度進(jìn)行評價.

1.3 風(fēng)險概率

給定不同的變量條件可以得到不同的條件概率分布.主要考慮以下兩種條件風(fēng)險概率.

1)給定X3=x3，條件概率分布函數(shù)可表示為

式中fX3(x3)為x3的概率密度函數(shù)，其他符號意義同前.

2)給定X3≤x3，條件概率分布函數(shù)可表示為

2 研究區(qū)域概況

廣州地處廣東省東南部，珠江三角洲北緣，范圍在東經(jīng) 112°57'～114°3'，北緯 22°26'～23°56'，瀕臨南海，毗鄰港澳.廣州地處亞熱帶沿海，屬海洋性亞熱帶季風(fēng)氣候，全年平均氣溫21.97℃，平均相對濕度68%，無霜期300～341 d，日照時數(shù)為1 571～2 053 h，市區(qū)年降雨量在1 700 mm以上.全年中，4～6月為雨季，7～9月天氣炎熱，多臺風(fēng).

以廣州市天河區(qū)五山站1961～2012年逐分鐘自記數(shù)據(jù)為例，采用降雨強(qiáng)度法［16］將連續(xù)的降雨時間序列分割成降雨事件，按照“降雨強(qiáng)度小于0.1 mm/h，且持續(xù)時間超過10 h”的標(biāo)準(zhǔn)進(jìn)行降雨數(shù)據(jù)分割，將1961～2012年的降雨時間序列分割為7 833場降雨事件，然后采用年最大值法進(jìn)行取樣.以設(shè)計歷時720 min(即12 h)為例，取樣樣本52個，進(jìn)而統(tǒng)計“設(shè)計降雨量”、“次降雨量”和“降雨峰值”3個表征雨型的特征變量，計算3者的聯(lián)合概率分布和各條件下的遭遇概率分布.

3 結(jié)果分析

3.1 邊際分布確定

圖1為“設(shè)計降雨量”、“降雨峰值”和“次降雨量”邊際分布的擬合結(jié)果.可以看出，所采用的3種分布線型與實(shí)測點(diǎn)據(jù)吻合較好，不同線型稍有差距，但總體上看均具有較好的擬合效果.

為了定量評價各分布線型的擬合效果，分別計算了各分布線型對于不同觀測數(shù)據(jù)的RMSE、MAE及PPCC，結(jié)果見表2.對設(shè)計降雨而言，LN2分布具有最小的RMSE值，而GEV分布具有最小的MAE值和最大的 PPCC值.綜合考慮，選定GEV分布作為數(shù)據(jù)的分布線型;對降雨峰值來說，GEV分布具有最小的MAE值，而LN2分布具有最小的RMSE值和最大的PPCC值，因此，選定LN2分布作為降雨峰值數(shù)據(jù)的分布線型;GEV分布對次降雨量具有最小的MAE、RMSE值和最大的PPCC值，所以，選定GEV分布作為次降雨數(shù)據(jù)的分布線型.

圖1 設(shè)計降雨量、降雨峰值和次降雨量邊際分布擬合

表2 邊際分布擬合優(yōu)度檢驗(yàn)

3.2 Copula函數(shù)參數(shù)估計及擬合檢驗(yàn)分析

采用阿基米德型Copula家族中應(yīng)用較廣泛的Frank Copula和AMH Copula函數(shù)分別構(gòu)建降雨特征變量的聯(lián)合分布.一般地，不同形式的Copula函數(shù)對于變量的相關(guān)性有不同的要求，F(xiàn)rank Copula函數(shù)對正、負(fù)相關(guān)關(guān)系的隨機(jī)變量均適用;AMH Copula函數(shù)適用于弱相關(guān)關(guān)系的隨機(jī)變量［13］.

為檢驗(yàn)Copula函數(shù)的擬合精度，首先比較各個觀測點(diǎn)對(xi，yi，zi)的經(jīng)驗(yàn)累積概率和理論累積概率的一致性［9，13］.圖 2 為經(jīng)驗(yàn)頻率與不同Copula函數(shù)理論頻率的擬合圖，可以看出，經(jīng)驗(yàn)頻率點(diǎn)和理論頻率點(diǎn)均分布在45°線附近，總體上看，經(jīng)驗(yàn)頻率點(diǎn)和理論頻率點(diǎn)擬合較好.Frank Copula函數(shù)的擬合精度(圖 2(b))優(yōu)于 AMH Copula函數(shù)(圖2(a)).

圖2 經(jīng)驗(yàn)累積概率與理論累積概率比較

為了進(jìn)一步確定“設(shè)計降雨量”、“降雨峰值”和“次降雨量”擬合最好的 Copula函數(shù)，利用MSE、AIC和OLS進(jìn)行擬合優(yōu)度評價.從表3可以看出，F(xiàn)rank Copula函數(shù)的3個評價指標(biāo)值均小于AMH Copula函數(shù)，表明Frank Copula函數(shù)對3變量的擬合效果最佳.所以，選取Frank Copula函數(shù)作為“設(shè)計降雨量”、“降雨峰值”和“次降雨量”的聯(lián)合分布函數(shù).

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

3.3 風(fēng)險概率分析

運(yùn)用Frank Copula函數(shù)計算“設(shè)計降雨量”、“降雨峰值”和“次降雨量”的聯(lián)合分布函數(shù)，根據(jù)1.3節(jié)公式計算條件風(fēng)險概率.分別以100 a重現(xiàn)期的“設(shè)計降雨量”、“降雨峰值”和“次降雨量”為條件，計算其他兩兩變量組合的風(fēng)險概率.從圖3可以看出，與單變量等于某一設(shè)計值條件下的遭遇概率分布(圖3(a)、(c)、(e))相比，對應(yīng)條件下單變量小于等于某一設(shè)計值下的遭遇概率分布圖(圖3(b)、(d)、(f))顯得“高”、“胖”，說明小于等于某一設(shè)計值條件下的風(fēng)險概率出現(xiàn)的可能性更大.

以設(shè)計降雨條件為例，給出了幾種重現(xiàn)期下降雨峰值和次降雨量遭遇的風(fēng)險概率.表4為設(shè)計降雨量等于某一設(shè)計值(100，50，20和10 a)的條件下，降雨峰值和次降雨量同時小于等于某一重現(xiàn)期設(shè)計值的概率.可以看出，在設(shè)計降雨不變的情況下，降雨峰值和次降雨量遭遇的風(fēng)險概率隨著重現(xiàn)期的減小而降低，當(dāng)降雨峰值和次降雨量設(shè)計值不變時，其遭遇的風(fēng)險概率隨著設(shè)計降雨重現(xiàn)期的減小而上升.具體來說，當(dāng)設(shè)計降雨等于100一遇時，小于等于100 a一遇降雨峰值和次降雨量同時遭遇的概率為92.18%，而小于等于10 a一遇同時遭遇的概率僅為48.44%.表5為設(shè)計降雨小于等于某一設(shè)計值(100，50，20和10 a)的條件下，降雨峰值和次降雨量同時小于等于某一重現(xiàn)期設(shè)計值的概率.對比表4和表5可知，兩種條件下的風(fēng)險概率變化具有相似性，但設(shè)計降雨小于等于某一設(shè)計值條件下的風(fēng)險概率出現(xiàn)的可能性更大.

從以上分析可以看出，設(shè)計降雨量、降雨峰值和次降雨量3變量之間具有較強(qiáng)的關(guān)聯(lián)性，用多變量聯(lián)合分布描述暴雨不同特征值之間的相互關(guān)系可以更全面地描述整個暴雨事件，細(xì)化水文分析前端輸入，提高分析精度.具體地說，如進(jìn)行50 a一遇水文分析時，根據(jù)表4，除給定設(shè)計降雨量值(218.9 mm)外，可根據(jù)城市等級及對安全度的需要，結(jié)合自身經(jīng)濟(jì)條件，同時選取某重現(xiàn)期的降雨峰值以及次降雨量作為輸入，如50 a一遇的降雨峰值4.6 mm/min或 50 a一遇次降雨量271.6 mm來確定洪澇設(shè)施規(guī)模.

圖3 降雨特征變量組合的條件風(fēng)險概率

表4 設(shè)計降雨量等于某一設(shè)計值的條件風(fēng)險概率

表5 設(shè)計降雨量小于等于某一設(shè)計值的條件風(fēng)險概率

4 結(jié)論

1)基于3維Copula函數(shù)的降雨多變量分析方法，根據(jù)歷史降雨數(shù)據(jù)，利用Copula聯(lián)結(jié)函數(shù)，實(shí)現(xiàn)了3變量的聯(lián)合分析，有效解決了單變量分析方法的不足，實(shí)現(xiàn)了對降雨事件特征真實(shí)全面的反映.

2)實(shí)際應(yīng)用中，可基于設(shè)計降雨量確定設(shè)施標(biāo)準(zhǔn)，利用三者的組合風(fēng)險概率復(fù)核次降雨量和降雨峰值的失效概率，為更科學(xué)合理地分析、評估和確定城市防洪排澇設(shè)施能力提供有效支撐.

3)如何基于地區(qū)特征、安全性和經(jīng)濟(jì)性的定量分析，優(yōu)化降雨變量組合，科學(xué)合理地確定防洪排澇設(shè)施標(biāo)準(zhǔn)，是該方法需進(jìn)一步研究的重點(diǎn).

［1］ZHANG N，GUO S L，XIAO Y，et al.Design storm method based on bivariate joint distribution［J］.Water Power，2008，34(1):18-21.

［2］FONTANAZZA C M，F(xiàn)RENI G，LA LOGGIA G，et al.Uncertainty evaluation of design rainfall for urban flood risk analysis［J］.Water Science ＆ Technology，2011，63(11):2641-2650.

［3］SALVADORI G，DE MICHELE C.Frequency analysis via copulas:theoretical aspects and applications to hydrological events ［J］.Water Resources Research，2004，40(12):W12511.

［4］NELSEN R.An introduction to copulas［M］.New York:Springer Verlag，2006.

［5］ZHANG L，SINGH V P.Gumbel-Hougaard copula for trivariate rainfall frequency analysis［J］.Journal of Hydrologic Engineering，2007，12(4):409-419.

［6］武傳號，黃國如，吳思遠(yuǎn).基于Copula函數(shù)的廣州市短歷時暴雨與潮位組合風(fēng)險分析［J］.水力發(fā)電學(xué)報，2014，33(2):33-41.

［7］ZHANG Q，LI J，SINGH V P，et al.Copula-based spatio-temporal patterns of precipitation extremes in China［J］.International Journal of Climatology，2013，33(5):1140-1152.

［8］ ZHANG L，SINGH V P.Bivariate flood frequency analysis using the copula method ［J］.Journal of Hydrologic Engineering，2006，11(2):150-164.

［9］侯蕓蕓，宋松柏，趙麗娜，等.基于Copula函數(shù)的3變量洪水頻率研究［J］.西北農(nóng)林科技大學(xué)學(xué)報:自然科學(xué)版，2010，38(2):219-228.

［10］SHIAU J T.Fitting drought duration and severity with two-dimensional copulas ［J］. Water Resources Management，2006，20(5):795-815.

［11］YU K X，XIONG L，GOTTSCHALK L.Derivation of low flow distribution functions using copulas ［J］.Journal of Hydrology，2014，508:273-288.

［12］李興凱，陳元芳.暴雨頻率分布線型優(yōu)選方法的研究［J］.水文，2010，30(2):50-53.

［13］謝華，羅強(qiáng)，黃介生.基于三維copula函數(shù)的多水文區(qū)豐枯遭遇分析［J］.水科學(xué)進(jìn)展，2012，23(2):186-193.

［14］MYUNG I J.Tutorial on maximum likelihood estimation［J］.Journal of Mathematical Psychology，2003，47:90-100.

［15］EMBRECHTS P，H?ING A，JURI A.Using copulae to bound the value-at-risk for functions of dependent risks［J］.Finance and Stochastics，2003，7(2):145-167.

［16］POWELL D N，AZIZ N M，KHAN A A.Impact of new rainfall patterns on the design of hydraulic structures［C］// Design of Hydraulic Structures. World Environmental and Water Resource Congress.Omaha:［s.n.］，2006:1-10.