3維Copula函數在降雨特征多變量頻率分析中的應用

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

(1.北京工業大學建筑工程學院，100124北京;2.廣州市市政工程設計研究院，510060廣州;3.北京市水質科學與水環境恢復重點實驗室，100124北京)

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

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

1 研究方法

1.1 Copula 函數簡介

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

表1 3維Copula函數

1.2 參數估計及擬合檢驗

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

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

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

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

采用ML方法估計Copula中的參數

Copula函數在實際應用中的一個主要問題就是函數形式的選擇.Embrechts等［15］對不同的Copula函數模型進行了比較，發現采用不同形式的Copula函數可能導致不同的分析結果，因此，選擇合適的Copula函數尤為重要.檢驗與評價Copula函數的方法較多，本文采用AIC信息準則法(AIC)及離差平方和最小(OLS)準則法［9］對Copula函數的擬合優度進行評價.

1.3 風險概率

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

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

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

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

2 研究區域概況

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

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

3 結果分析

3.1 邊際分布確定

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

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

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

表2 邊際分布擬合優度檢驗

3.2 Copula函數參數估計及擬合檢驗分析

采用阿基米德型Copula家族中應用較廣泛的Frank Copula和AMH Copula函數分別構建降雨特征變量的聯合分布.一般地，不同形式的Copula函數對于變量的相關性有不同的要求，Frank Copula函數對正、負相關關系的隨機變量均適用;AMH Copula函數適用于弱相關關系的隨機變量［13］.

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

圖2 經驗累積概率與理論累積概率比較

為了進一步確定“設計降雨量”、“降雨峰值”和“次降雨量”擬合最好的 Copula函數，利用MSE、AIC和OLS進行擬合優度評價.從表3可以看出，Frank Copula函數的3個評價指標值均小于AMH Copula函數，表明Frank Copula函數對3變量的擬合效果最佳.所以，選取Frank Copula函數作為“設計降雨量”、“降雨峰值”和“次降雨量”的聯合分布函數.

表3 3維Copula函數擬合檢驗結果

3.3 風險概率分析

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

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

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

圖3 降雨特征變量組合的條件風險概率

表4 設計降雨量等于某一設計值的條件風險概率

表5 設計降雨量小于等于某一設計值的條件風險概率

4 結論

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

2)實際應用中，可基于設計降雨量確定設施標準，利用三者的組合風險概率復核次降雨量和降雨峰值的失效概率，為更科學合理地分析、評估和確定城市防洪排澇設施能力提供有效支撐.

3)如何基于地區特征、安全性和經濟性的定量分析，優化降雨變量組合，科學合理地確定防洪排澇設施標準，是該方法需進一步研究的重點.

［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，FRENI 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］武傳號，黃國如，吳思遠.基于Copula函數的廣州市短歷時暴雨與潮位組合風險分析［J］.水力發電學報，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函數的3變量洪水頻率研究［J］.西北農林科技大學學報:自然科學版，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］李興凱，陳元芳.暴雨頻率分布線型優選方法的研究［J］.水文，2010，30(2):50-53.

［13］謝華，羅強，黃介生.基于三維copula函數的多水文區豐枯遭遇分析［J］.水科學進展，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.