基于G-H copula函數(shù)的秦淮河流域洪水風(fēng)險(xiǎn)分析

高玉琴， 葉 柳， 賴麗娟

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

1 研究背景

隨著洪水風(fēng)險(xiǎn)問題研究越來越深入，從多角度對(duì)洪水進(jìn)行定義和描述的多變量聯(lián)合風(fēng)險(xiǎn)分析方法已經(jīng)越來越受到關(guān)注。Copula函數(shù)因?yàn)槠錁?gòu)造靈活，計(jì)算簡(jiǎn)單，且不限制邊際分布類型等優(yōu)點(diǎn)而備受青睞，已被大量應(yīng)用于降雨、干旱等多變量聯(lián)合概率分析以及洪水遭遇、多站點(diǎn)多水文區(qū)的洪水聯(lián)合風(fēng)險(xiǎn)分析等水文領(lǐng)域。

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

基于次洪變量，如時(shí)段洪量、洪峰、峰現(xiàn)時(shí)間、洪水歷時(shí)等進(jìn)行的聯(lián)合分布研究以及降雨和洪水、降雨和潮位相關(guān)性研究較多，但是目前從流域洪水風(fēng)險(xiǎn)角度將影響風(fēng)險(xiǎn)大小的洪量、洪峰流量和洪水位三者進(jìn)行聯(lián)合分析的研究較少。本文以秦淮河流域?yàn)檠芯繀^(qū)域，選取P-Ⅲ、 GEV分布和兩參數(shù)對(duì)數(shù)正態(tài)分布Ln2描述洪水風(fēng)險(xiǎn)變量洪量、洪峰和洪水位的邊際分布類型，利用G-H copula構(gòu)建二維和三維風(fēng)險(xiǎn)評(píng)價(jià)模型，計(jì)算不同組合洪水事件的聯(lián)合重現(xiàn)期、條件重現(xiàn)期以及二次重現(xiàn)期，為區(qū)域工程設(shè)計(jì)和風(fēng)險(xiǎn)評(píng)估工作提供依據(jù)。

2 數(shù)據(jù)與方法

2.1 研究區(qū)域與數(shù)據(jù)來源

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

2.2 G-H copula函數(shù)

G-H Copula是Archimedean copula函數(shù)的一種，通常用于描述正相關(guān)和上尾相關(guān)性隨機(jī)變量[11]，當(dāng)維數(shù)為2和3時(shí)，G-H Copula表達(dá)式分別為：

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 參數(shù)估計(jì)與擬合優(yōu)度檢驗(yàn)

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

為驗(yàn)證Copula和實(shí)測(cè)樣本經(jīng)驗(yàn)累積概率之間的擬合程度，需對(duì)構(gòu)建的模型進(jìn)行擬合優(yōu)度檢驗(yàn)。常見檢驗(yàn)方法有RMSE法、AIC法和BIC法。本文采用均方根誤差法：

(3)

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

(4)

式中：njkl為同時(shí)滿足Xi1≤xi1、Xi2≤xi2、Xi3≤xi3時(shí)的聯(lián)合觀測(cè)個(gè)數(shù)。

實(shí)施國土資源所巡察“五個(gè)三”行動(dòng)計(jì)劃 整體提升基層國土資源一線工作水平（殷金蘭） ..........................3-12

2.4 重現(xiàn)期

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

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

(5)

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

(6)

對(duì)于三維G-H copula， 可以表示為：

(7)

3 結(jié)果與討論

3.1 相依性度量和邊際分布

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

表1 變量之間的相依性

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

表2 邊際分布參數(shù)估計(jì)和擬合優(yōu)度檢驗(yàn)

圖1 邊際分布擬合圖

3.2 聯(lián)合分布函數(shù)

采用G-HCopula描述二維變量間相依結(jié)構(gòu)，用IFM法估計(jì)參數(shù)，計(jì)算得洪量洪峰、洪量洪水位和洪峰洪水位二維聯(lián)合分布函數(shù)的參數(shù)分別為3.1231、2.8655、6.7281，對(duì)應(yīng)的均方根誤差RMSE分別為0.0090、0.0106、0.0149。以相關(guān)性最強(qiáng)的洪峰和洪水位為主進(jìn)行二維變量洪水風(fēng)險(xiǎn)分析，探討其聯(lián)合重現(xiàn)期、同現(xiàn)重現(xiàn)期和Kendall重現(xiàn)期性質(zhì)，并分析洪峰不超過一定重現(xiàn)期時(shí)洪水位風(fēng)險(xiǎn)率。

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

繪制洪峰和洪水位的概率分布圖和相應(yīng)的重現(xiàn)期等值線圖，并加入實(shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比，如圖2和圖3。由圖2、3可知，Q-Z最大聯(lián)合重現(xiàn)期約為15 a，最大同現(xiàn)重現(xiàn)期則超過20 a?？捎?jì)算任意Q-Z組合事件的重現(xiàn)期，如發(fā)生50年一遇洪峰和20年一遇洪水位的聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期分別為19.99 a和50.03 a。繪制W-Q和W-Z的同現(xiàn)重現(xiàn)期等值線圖，如圖4。Q-Z的同現(xiàn)重現(xiàn)期較W-Q和W-Z的同現(xiàn)重現(xiàn)期小，說明秦淮河流域洪峰和洪水位同現(xiàn)風(fēng)險(xiǎn)率最大。

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

表3 單變量設(shè)計(jì)值及各種重現(xiàn)期

圖2 Q-Z聯(lián)合概率分布圖及聯(lián)合重現(xiàn)期等值線圖

圖3 Q-Z同現(xiàn)概率分布圖及同現(xiàn)重現(xiàn)期等值線圖

圖4 W-Q、W-Z同現(xiàn)重現(xiàn)期等值線圖

圖5 洪峰不超過某一重現(xiàn)期時(shí)洪水位風(fēng)險(xiǎn)率

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

因子后，傳統(tǒng)重現(xiàn)期、二次重現(xiàn)期與單因素重現(xiàn)期之間的差值明顯增大，風(fēng)險(xiǎn)率變化幅度增大。

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

表4 3變量同頻率設(shè)計(jì)值及各種重現(xiàn)期

圖6 Z≤10.8 m時(shí)W-Q條件概率分布及條件重現(xiàn)期等值線

4 結(jié) 論

以秦淮河流域?yàn)槔谶呺H分布擬合，選擇廣義極值分布GEV，對(duì)數(shù)正態(tài)分布ln2和P-Ⅲ分布作為洪量(W)、洪峰(Q)、洪水位(Z)的邊際分布類型，構(gòu)建二維和三維G-Hcopula洪水風(fēng)險(xiǎn)評(píng)價(jià)模型，并計(jì)算多變量聯(lián)合概率、條件概率以及傳統(tǒng)重現(xiàn)期和Kendall重現(xiàn)期，得出了以下結(jié)論：

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

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

(3)計(jì)算3變量聯(lián)合分布和限制洪水位條件下W-Q聯(lián)合概率分布，3變量情況下的同頻率設(shè)計(jì)值大于相應(yīng)單變量設(shè)計(jì)值，3變量傳統(tǒng)重現(xiàn)期、Kendall重現(xiàn)期與單因素重現(xiàn)期之間的差值明顯增大。

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

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

[3] 閆寶偉，郭生練，陳 璐，等. 長江和清江洪水遭遇風(fēng)險(xiǎn)分析[J]. 水利學(xué)報(bào)，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] 陳子燊，黃 強(qiáng)，劉曾美. 基于非對(duì)稱ArchimedeanCopula的三變量洪水風(fēng)險(xiǎn)評(píng)估[J]. 水科學(xué)進(jìn)展，2016,27(5)：763-771.

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

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

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

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

[12] 黃 強(qiáng)，陳子燊. 基于二次重現(xiàn)期的多變量洪水風(fēng)險(xiǎn)評(píng)估[J]. 湖泊科學(xué)，2015,27(2)：352-360.

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

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