基于G-H copula函數的秦淮河流域洪水風險分析

高玉琴， 葉 柳， 賴麗娟

(河海大學 水利水電學院， 江蘇 南京 210098)

1 研究背景

隨著洪水風險問題研究越來越深入，從多角度對洪水進行定義和描述的多變量聯合風險分析方法已經越來越受到關注。Copula函數因為其構造靈活，計算簡單，且不限制邊際分布類型等優點而備受青睞，已被大量應用于降雨、干旱等多變量聯合概率分析以及洪水遭遇、多站點多水文區的洪水聯合風險分析等水文領域。

二維copula研究理論及應用已經相對成熟，三維及更高維數copula函數由于維數的拓展和各變量間相依性不對等，理論更為復雜。涂新軍等[1]應用Archimedean copula進行了濱海城市降雨和潮位聯合分布的模擬和設計，并對比基于二次重現期的設計值和傳統聯合重現期設計值，證明二次重現期設計更為安全。黃錦林等[2]基于最大熵-copula方法分析降雨和潮位關聯性，得出多種雨潮設計值遭遇組合風險概率。閆寶偉等[3]采用clayton copula構建長江和清江峰現時間與量級二維分布，進而探索洪水遭遇風險特性。Chowdhary等[4]詳細論述了copula函數類型優選方法并介紹其在二元洪水頻率計算中的應用。Zhang等[5]利用Archimedean copula函數研究了降雨強度、徑流深、降雨歷時兩兩變量間的聯合分布，并分析相依性和邊緣分布對結果的影響；Reddy等[6]具體闡述了二維洪水變量頻率分析的步驟流程以及用Monte Carlo數值模擬法判斷copula擬合程度的方法。在三維copula方面，陳子燊等[7]分析洪量、歷時、洪峰三變量聯合分布及各種重現期的設計分位數；侯蕓蕓等[8]考慮洪峰、洪量和歷時之間的相依性，構造了三維對稱性Archimedean copula函數，探討洪水聯合概率分布和條件概率分布；Ganguli等[9]利用多種聯結函數類型描述洪峰流量、洪量、洪水歷時的兩變量、三變量聯合概率分布，并討論了兩種重現期特性和尾部相關性。

基于次洪變量，如時段洪量、洪峰、峰現時間、洪水歷時等進行的聯合分布研究以及降雨和洪水、降雨和潮位相關性研究較多，但是目前從流域洪水風險角度將影響風險大小的洪量、洪峰流量和洪水位三者進行聯合分析的研究較少。本文以秦淮河流域為研究區域，選取P-Ⅲ、 GEV分布和兩參數對數正態分布Ln2描述洪水風險變量洪量、洪峰和洪水位的邊際分布類型，利用G-H copula構建二維和三維風險評價模型，計算不同組合洪水事件的聯合重現期、條件重現期以及二次重現期，為區域工程設計和風險評估工作提供依據。

2 數據與方法

2.1 研究區域與數據來源

秦淮河流域地處長江下游地區，面積為2 631 km2，干流長度為34 km，屬于亞熱帶濕潤、半濕潤季風氣候區，汛期雨量充沛，洪水災害頻繁。根據1986-2006年流域出口處秦淮新河和武定門閘水文站的日流量資料和1986-2006年東山站站前水位資料，采用年最大值選樣法(AM)，基于袁玉等[10]根據秦淮河流域水文、氣象、降雨等資料建立的HEC-HMS水文模型及HEC-RAS水力模型，模擬年最大洪水事件演進到東山站前的水位。洪量、洪峰均為秦淮新河閘與武定閘兩處疊加值，以全流域徑流深的形式表示洪量，洪水位為東山站點前水位。

2.2 G-H copula函數

G-H Copula是Archimedean copula函數的一種，通常用于描述正相關和上尾相關性隨機變量[11]，當維數為2和3時，G-H Copula表達式分別為：

C(u1,u2)=exp{-[(-ln(u1))θ+

(-ln(u2))θ]1/θ}θ∈(1,∞)

(1)

C(u1,u2,u3)=exp{-[(-ln(u1))θ+(-ln(u2))θ+

(-ln(u3))θ]1/θ}θ∈(1,∞)

(2)

2.3 參數估計與擬合優度檢驗

核密度估計、相關性指標法、IFM法、EML法等都可進行copula參數估計。對于二維ArchimedeanCopula常用相關性指標法，根據τ和θ的關系由τ來估計。但三維copula顯然不能使用相關性指標法估算參數，本研究采用IFM法。

為驗證Copula和實測樣本經驗累積概率之間的擬合程度，需對構建的模型進行擬合優度檢驗。常見檢驗方法有RMSE法、AIC法和BIC法。本文采用均方根誤差法：

(3)

式中：Femp為經驗概率值；C為理論概率值；n為系列長度。RMSE值越小，擬合情況越好。三維經驗頻率的計算公式為：

(4)

式中：njkl為同時滿足Xi1≤xi1、Xi2≤xi2、Xi3≤xi3時的聯合觀測個數。

實施國土資源所巡察“五個三”行動計劃 整體提升基層國土資源一線工作水平（殷金蘭） ..........................3-12

2.4 重現期

洪水特征變量的聯合概率分布、條件概率分布以及相應的重現期計算是水文頻率分析的重要成果，超過某一重現期設計閾值的洪水事件被定義為危險事件。傳統同現重現期和聯合重現期對洪水事件安全和危險域識別存在局限性，有可能發生識別錯誤的情況，造成風險管理和工程設計的成本增加或風險抵抗能力不足，而Kendall重現期能在同一臨界水平下對任意洪水組合的安全和危險域識別保持一致性[12-13]。本文將Kendall重現期應用到多因素風險分析中，并與傳統重現期進行對比。傳統重現期定義見參考文獻[7]，在此不再贅述。

Kendall重現期，又稱二次重現期，根據Genest和Rivest定義的Kendall分布函數，對于給定的t，通過求解特征量聯合分布C小于或等于t的累積概率將多維變量信息投射為一維。Kendall重現期定義為：

(5)

式中：KC為Kendall分布函數，KC(t)=P(C(u1,u2)≤t))。對于Archimedean copula函數族中的二維G-H copula，KC可以表示為：

(6)

對于三維G-H copula， 可以表示為：

(7)

3 結果與討論

3.1 相依性度量和邊際分布

采用Pearson′γ，Spearman′ρ和Kendall′τ進行洪量和洪峰、洪量和洪水位，洪峰和洪水位兩兩變量間的相依性度量，并在括號中給出顯著性水平為0.05時對應的P值，計算結果如表1。由表1可以看出，變量之間的相關性均較大，表明洪量(W)、洪峰(Q)、洪水位(Z)之間存在較強相關性，因而可進行洪水變量間的聯合概率特性分析。

表1 變量之間的相依性

選取P-Ⅲ、GEV和ln2描述變量的邊際分布，邊際分布概率密度表達式見文獻[14]。采用最大似然法估計參數并通過KS檢驗判斷擬合優度，KS統計量值越小，擬合情況越好；P值越大，則有越大可能接受該分布。表2為洪水特征變量邊際分布參數和擬合優度檢驗結果。由表2可知，3種邊際分布均能被接受，基于最小KS統計量分別選擇GEV，ln2和P-Ⅲ分布為洪量、洪峰、洪水位的邊際分布。繪制變量邊際分布經驗頻率與理論頻率關系圖，如圖1所示，由圖1可以看出，點據分布在1∶1斜線附近，證明擬合程度較好。

表2 邊際分布參數估計和擬合優度檢驗

圖1 邊際分布擬合圖

3.2 聯合分布函數

采用G-HCopula描述二維變量間相依結構，用IFM法估計參數，計算得洪量洪峰、洪量洪水位和洪峰洪水位二維聯合分布函數的參數分別為3.1231、2.8655、6.7281，對應的均方根誤差RMSE分別為0.0090、0.0106、0.0149。以相關性最強的洪峰和洪水位為主進行二維變量洪水風險分析，探討其聯合重現期、同現重現期和Kendall重現期性質，并分析洪峰不超過一定重現期時洪水位風險率。

分別計算洪峰和洪水位在5～200a重現期內的聯合T∪，同現T∩和二次重現期TK值，并給出單變量條件下的設計值，見表3。Q-Z二維分布條件下，相較于單因素重現期，T∪更小，T∩更大，而Kendall重現期在T∪和T∩之間且與單因素重現期誤差最小，大小順序為T∪

繪制洪峰和洪水位的概率分布圖和相應的重現期等值線圖，并加入實測數據進行對比，如圖2和圖3。由圖2、3可知，Q-Z最大聯合重現期約為15 a，最大同現重現期則超過20 a?？捎嬎闳我釷-Z組合事件的重現期，如發生50年一遇洪峰和20年一遇洪水位的聯合重現期和同現重現期分別為19.99 a和50.03 a。繪制W-Q和W-Z的同現重現期等值線圖，如圖4。Q-Z的同現重現期較W-Q和W-Z的同現重現期小，說明秦淮河流域洪峰和洪水位同現風險率最大。

繪制洪峰不超過10、20、50年一遇時洪水位的超越概率FZ|Q=P(Z>z|Q≤q)，如圖5。洪峰不超過10、20、50年一遇設計值時，洪水位超過5年一遇 (8.29m)的概率分別為0.1113、0.1579、0.1837。給定某一洪水位閾值，洪峰量級越小，洪水位風險率越小，即洪峰量級小，洪水位也越不容易超過該閾值。

表3 單變量設計值及各種重現期

圖2 Q-Z聯合概率分布圖及聯合重現期等值線圖

圖3 Q-Z同現概率分布圖及同現重現期等值線圖

圖4 W-Q、W-Z同現重現期等值線圖

圖5 洪峰不超過某一重現期時洪水位風險率

按照公式(3)，通過各變量邊緣分布構建三維G-Hcopula模型，用IFM法由洪水資料估計參數，計算得出θ=3.4472，用均方根誤差法對G-Hcopula函數進行擬合優度檢驗，RMSE=0.0147，誤差值較小 ，滿足擬合要求。分別計算洪量、洪峰和洪水位在5～200a重現期內的T∪、T∩、TK值和3變量情況下的同頻率設計值，如表4。由表4可以看出，3變量同頻率設計值大于相應單變量設計值，洪量洪峰洪水位三維分布條件下，各重現期大小排列依然是T∪

因子后，傳統重現期、二次重現期與單因素重現期之間的差值明顯增大，風險率變化幅度增大。

由《南京城市防洪規劃報告(2011-2020年)》可知，東山站20年一遇規劃水位為10.8m，繪制洪水位不超過10.8m條件下洪量洪峰的條件概率FWQ|Z=P(W≤w,Q≤q|Z≤10.8)及相應的重現期等值線圖，如圖6。洪水位不超過10.8m條件下，實測資料中大部分W-Q聯合重現期小于或等于5a，說明較高洪水位條件下，洪量和洪峰的任意一個超過風險率較大。

表4 3變量同頻率設計值及各種重現期

圖6 Z≤10.8 m時W-Q條件概率分布及條件重現期等值線

4 結 論

以秦淮河流域為例，基于邊際分布擬合，選擇廣義極值分布GEV，對數正態分布ln2和P-Ⅲ分布作為洪量(W)、洪峰(Q)、洪水位(Z)的邊際分布類型，構建二維和三維G-Hcopula洪水風險評價模型，并計算多變量聯合概率、條件概率以及傳統重現期和Kendall重現期，得出了以下結論：

(1)相較于單因素重現期，Q-Z二維T∪更小，T∩更大，而Kendall重現期在兩者之間且與單因素重現期誤差最小。

(2)Q-Z的T∪小于W-Q和W-Z的T∩，說明秦淮河流域洪峰和洪水位同現風險率最大；Q-Z在二維分布條件下，洪峰量級越小，則洪水位風險率越小。

(3)計算3變量聯合分布和限制洪水位條件下W-Q聯合概率分布，3變量情況下的同頻率設計值大于相應單變量設計值，3變量傳統重現期、Kendall重現期與單因素重現期之間的差值明顯增大。

[1] 涂新軍，杜奕良，陳曉宏，等. 濱海城市雨潮遭遇聯合分布模擬與設計[J]. 水科學進展，2017,28(1)：49-58.

[2] 黃錦林，范嘉煒，唐造造. 基于最大熵-Copula方法的降雨潮位關聯性分析——以廣州為例[J]. 災害學，2017,32(1)：65-71.

[3] 閆寶偉，郭生練，陳 璐，等. 長江和清江洪水遭遇風險分析[J]. 水利學報，2010，41(5)：553-559.

[4]CHOWDHARYH,ESCOBARLA,SINGHVP.Identificationofsuitablecopulasforbivariatefrequencyanalysisoffloodpeakandfloodvolumedata[J].HydrologyResearch，2011, 42(2-3)：193.

[5]ZHANGL,SINGHVP.BivariaterainfallfrequencydistributionsusingArchimedeancopulas[J].JournalofHydrology，2007, 332(1)：93-109.

[6]REDDYMJ,GANGULIP.Bivariatefloodfrequencyanalysisofuppergodavaririverflowsusingarchimedeancopulas[J].WaterResourcesManagement，2012, 26(14)：3995-4018.

[7] 陳子燊，黃 強，劉曾美. 基于非對稱ArchimedeanCopula的三變量洪水風險評估[J]. 水科學進展，2016,27(5)：763-771.

[8] 侯蕓蕓，宋松柏，趙麗娜，等. 基于Copula函數的3變量洪水頻率研究[J]. 西北農林科技大學學報(自然科學版)，2010,38(2)：219-228.

[9]GANGULIP,REDDYMJ.Probabilisticassessmentoffloodrisksusingtrivariatecopulas[J].Theoretical&AppliedClimatology, 2013,111(1-2):341-360.

[10] 袁 玉，高玉琴，吳 錫. 基于HEC_HMS水文模型的秦淮河流域圩垸式防洪模式洪水模擬[J]. 三峽大學學報(自然科學版)，2015,37(5)：34-39.

[11] 潘璀林，陳子燊. 基于GHCopula的韓江水文干旱聯合概率分布研究[J]. 中山大學學報(自然科學版)，2015,54(1)：110-115.

[12] 黃 強，陳子燊. 基于二次重現期的多變量洪水風險評估[J]. 湖泊科學，2015,27(2)：352-360.

[13] 范嘉煒，黃錦林. 基于Kendall重現期的降雨潮位風險分析[J]. 水電能源科學，2017,35(5):21-24+20.

[14] 高玉琴，陳釔西，趙立梅，等. 秦淮河流域不同頻率降雨聯合概率分析[J]. 水電能源科學，2016,34(3):1-5+23.