基于二次重現(xiàn)期的多變量洪水風(fēng)險(xiǎn)評(píng)估*

黃 強(qiáng)，陳子燊

(中山大學(xué)水資源與環(huán)境系，廣州 510275)

基于二次重現(xiàn)期的多變量洪水風(fēng)險(xiǎn)評(píng)估*

黃 強(qiáng)，陳子燊**

(中山大學(xué)水資源與環(huán)境系，廣州 510275)

由于洪水是一種具有多個(gè)特征屬性的隨機(jī)事件，頻率分析成為洪水風(fēng)險(xiǎn)評(píng)估的一種有效手段，多變量重現(xiàn)期與設(shè)計(jì)值的定義與計(jì)算則是洪水頻率分析中的重點(diǎn)和難點(diǎn).本文通過構(gòu)造洪水歷時(shí)、洪峰與洪量的聯(lián)合分布，介紹了一種新的多變量重現(xiàn)期定義——二次重現(xiàn)期，并探討了“或”重現(xiàn)期、“且”重現(xiàn)期和二次重現(xiàn)期對(duì)安全與危險(xiǎn)域識(shí)別的差異性，以及在洪水風(fēng)險(xiǎn)管理與工程設(shè)計(jì)中的合理性與可靠性.傳統(tǒng)的“或”和“且”多變量重現(xiàn)期對(duì)安全與危險(xiǎn)域的識(shí)別存在局限性，利用Kendall函數(shù)定義的二次重現(xiàn)期則提供了更加合理的安全與風(fēng)險(xiǎn)域識(shí)別，避免了對(duì)安全事件與危險(xiǎn)事件的錯(cuò)誤判定，更有利于指導(dǎo)洪水風(fēng)險(xiǎn)的管理.在給定的二次重現(xiàn)期條件下，依據(jù)出現(xiàn)概率最大原則推算的歷時(shí)、洪峰與洪量設(shè)計(jì)值組合可以滿足工程設(shè)計(jì)以較低成本承受較大風(fēng)險(xiǎn)的追求，相比于單變量設(shè)計(jì)值，考慮了洪水多個(gè)屬性聯(lián)合特征的多變量設(shè)計(jì)值提供了更加全面和可靠的參考信息.

多變量洪水特征；極值分布；安全與危險(xiǎn)域；Kendall函數(shù)；二次重現(xiàn)期；多變量設(shè)計(jì)值

洪水作為一種自然災(zāi)害事件，每年造成的經(jīng)濟(jì)損失數(shù)以億計(jì)，如何控制洪水并減少損失一直是災(zāi)害風(fēng)險(xiǎn)管理的重點(diǎn).由于洪水事件發(fā)生的隨機(jī)性，頻率分析成為洪水風(fēng)險(xiǎn)評(píng)估的一種有效手段，重現(xiàn)期與設(shè)計(jì)值的定義與計(jì)算則是頻率分析中的重點(diǎn).同一個(gè)洪水事件前后兩次發(fā)生的平均間隔時(shí)間即為該事件的重現(xiàn)期，若把該洪水事件造成的破壞和影響(以下稱效應(yīng))作為洪水風(fēng)險(xiǎn)管理與工程設(shè)計(jì)的標(biāo)準(zhǔn)，重現(xiàn)期更大的洪水事件則可視為該標(biāo)準(zhǔn)下的危險(xiǎn)事件.傳統(tǒng)的單變量頻率分析并不能滿足洪水事件多個(gè)特征屬性的特點(diǎn)，隨著Copula函數(shù)在水文分析中的應(yīng)用，多變量的洪水頻率分析已成為近年來的研究熱點(diǎn)[1-4].與單變量重現(xiàn)期相比，多變量重現(xiàn)期的定義則要復(fù)雜得多，“或”和“且”定義法[5]是目前最常用的兩種多變量重現(xiàn)期定義方法，但在安全事件與危險(xiǎn)事件的判定上，兩者都存在局限性[6].Salvadori等[6]在“或”和“且”多變量重現(xiàn)期定義法的基礎(chǔ)上，基于安全與危險(xiǎn)域的概念提出了新的多變量重現(xiàn)期定義——二次重現(xiàn)期，相關(guān)概念與定義已在海岸工程設(shè)計(jì)的研究中有了初步的應(yīng)用[7-8].

本文以東江下游博羅站點(diǎn)為例，介紹多變量重現(xiàn)期新的定義與計(jì)算方法，通過分析傳統(tǒng)重現(xiàn)期與二次重現(xiàn)期在安全與危險(xiǎn)域識(shí)別的差異性，探討二次重現(xiàn)期在洪水風(fēng)險(xiǎn)應(yīng)用的合理性與可靠性.本文的目的在于理論方法的介紹與探討，為洪水風(fēng)險(xiǎn)管理與工程設(shè)計(jì)的研究提供一種新的思路和方法.

1 洪水特征屬性

洪水事件具有多個(gè)特征屬性，其中最重要的是歷時(shí)、洪峰和洪量.一場洪水的選取關(guān)鍵在于辨別洪水開始和結(jié)束的時(shí)間，當(dāng)流量開始增加并超過了基流量時(shí)即可認(rèn)為是洪水開始發(fā)生，當(dāng)流量又退落到基流量時(shí)即意味著洪水的結(jié)束.然而，由于每一場洪水的基流量都不一樣，實(shí)際上要準(zhǔn)確地辨認(rèn)出基流量是十分困難的，一般可以將流量突然并持續(xù)增加的時(shí)刻看作是洪水開始的時(shí)間，將流量突然下降并退落至起漲流量或穩(wěn)定流量的時(shí)刻看作是洪水結(jié)束的時(shí)間[9].

洪水的特征屬性可以按照流量過程線來提取，如圖1所示，流量從突然增加到退落至穩(wěn)定所經(jīng)歷的時(shí)間為洪水歷時(shí)D(d)，其間最大的流量即為洪峰流量Q(m3/s)，所增加的水量為洪水總量(洪量)V(斜線部分面積，日m3/s).歷時(shí)和洪量可分別根據(jù)下式計(jì)算[9]：

D=te-ts

(1)

(2)

式中，ts和te分別為洪水開始和結(jié)束的時(shí)間，q為日流量觀測序列.

圖1 洪水過程及相應(yīng)的特征屬性Fig.1 The flood hydrograph and characteristics

2 邊緣分布與聯(lián)合分布

洪水頻率分析中的樣本序列主要通過年最大抽樣的方法從每年發(fā)生的洪水事件中選取，選取的年最大洪水事件可認(rèn)為是極值事件，根據(jù)極值理論，極值事件應(yīng)服從極值分布[10].極值理論中的單元極值分布為廣義極值分布(GEV)，其分布函數(shù)為：

(3)

式中，μ、σ和k分別為位置、尺度和形狀參數(shù).

與傳統(tǒng)的多變量分布函數(shù)相比，Copula函數(shù)具有簡單靈活的結(jié)構(gòu)特點(diǎn)，并且不受邊緣分布的限制，是一種可以適用于多個(gè)變量間不同相關(guān)結(jié)構(gòu)的多變量分布函數(shù)[11].Archimedean copula是水文分析中最常用的Copula函數(shù)族，其基本形式為：

(4)

式中，uj∈[0， 1](j>1)為邊緣分布，φθ為Archimedean copula生成元，θ為參數(shù).當(dāng)C滿足下式時(shí)，則C是一個(gè)多變量極值分布[13].

(5)

可以證明Gumbel-Hougaard copula滿足上式，并且是Archimedean copula中唯一的多變量極值copula[12-13]，其結(jié)構(gòu)形式為：

(6)

僅當(dāng)邊緣分布和Copula函數(shù)均具有極值性質(zhì)時(shí)，構(gòu)造的聯(lián)合分布才是極值分布[13]，可用于極值事件的頻率分析中.

3 多變量重現(xiàn)期與分位值

3.1 傳統(tǒng)多變量重現(xiàn)期

在傳統(tǒng)的多變量重現(xiàn)期定義中，E∨對(duì)應(yīng)的重現(xiàn)期為“或”重現(xiàn)期T∨，E∧對(duì)應(yīng)的重現(xiàn)期則為“且”重現(xiàn)期T∧.年最大洪水事件的3變量“或”重現(xiàn)期和“且”重現(xiàn)期可分別由以下兩式計(jì)算：

(7)

(8)

假設(shè)E的安全率為Psafe(危險(xiǎn)率則為1-Psafe)，在Psafe=p的條件下，以兩變量的情形進(jìn)行闡述，“或”重現(xiàn)期和“且”重現(xiàn)期對(duì)安全與危險(xiǎn)域的識(shí)別如圖2所示.在相同安全率的情況下，“或”重現(xiàn)期和“且”重現(xiàn)期標(biāo)準(zhǔn)下不同的分位值組合對(duì)應(yīng)的安全與危險(xiǎn)域都不一樣.

圖2 “或”重現(xiàn)期(a)和“且”重現(xiàn)期(b)的安全與危險(xiǎn)域識(shí)別(灰色部分為危險(xiǎn)域)Fig.2 Identification of safe and dangerous regions for “OR” return period(a) and “AND” return period(b) (the grey part represents the dangerous one)

一般地，設(shè)計(jì)的重現(xiàn)期越大，相應(yīng)的危險(xiǎn)率就越小，設(shè)置的災(zāi)害防御標(biāo)準(zhǔn)遭遇的風(fēng)險(xiǎn)就越低.就此而言，“或”重現(xiàn)期和“且”重現(xiàn)期識(shí)別的安全與危險(xiǎn)域顯然具有局限性[6].如圖3a所示，對(duì)于“或”重現(xiàn)期，A點(diǎn)處的重現(xiàn)期大于B點(diǎn)(p1>p2)，而A點(diǎn)識(shí)別的危險(xiǎn)域卻有一部分在B點(diǎn)識(shí)別的危險(xiǎn)域之外(深灰色部分)，即當(dāng)事件E*發(fā)生時(shí)，若以B點(diǎn)的組合值作為災(zāi)害防御的標(biāo)準(zhǔn)則可以認(rèn)為E*是安全的，E*造成的效應(yīng)在防御

的標(biāo)準(zhǔn)之下；若以重現(xiàn)期更大的A點(diǎn)組合值作為災(zāi)害防御的標(biāo)準(zhǔn)則E*顯然是危險(xiǎn)事件.高風(fēng)險(xiǎn)的防御標(biāo)準(zhǔn)識(shí)別的安全事件卻被更低風(fēng)險(xiǎn)的防御標(biāo)準(zhǔn)識(shí)別為危險(xiǎn)事件，這顯然是矛盾的.對(duì)于“且”重現(xiàn)期(圖3b)，A點(diǎn)處的重現(xiàn)期小于B點(diǎn)(p1

圖3 “或”重現(xiàn)期(a)和“且”重現(xiàn)期(b)對(duì)安全危險(xiǎn)域識(shí)別的局限性Fig.3 Limitations on the identification of safe and dangerous regions for “OR” return period(a) and “AND” return period(b)

3.2 分位值曲面

3.3 二次重現(xiàn)期

(9)

TK稱為二次重現(xiàn)期(Secondary return period)或Kendall重現(xiàn)期，KC為Kendall分布函數(shù)[16]：

(10)

(11)

對(duì)于三維Gumbel-Hougaard copula，Kendall分布函數(shù)可以表達(dá)為：

(12)

由公式(7)、(8)和(9)以及C的非遞減性可知，單變量重現(xiàn)期、“或”重現(xiàn)期、“且”重現(xiàn)期和二次重現(xiàn)期四者的關(guān)系[14，17]為：T∨≤T≤TK≤T∧，當(dāng)且僅當(dāng)uj=0時(shí)，等號(hào)成立.

圖4 二次重現(xiàn)期的安全與危險(xiǎn)域識(shí)別 (灰色部分為危險(xiǎn)域)Fig.4 Identification of safe and dangerous regions for secondary return period (the grey part represents the dangerous one)

3.4 多變量設(shè)計(jì)值

(13)

f(d，q，v)=c(ud，uq，uv)fd(d)fq(q)fv(v)

(14)

式中，c為三維Gumbel-Hougaard copula的概率密度函數(shù)，fd(d)、fq(q)和fv(v)則分別為歷時(shí)、洪峰和洪量的概率密度函數(shù).

綜上所述，二次重現(xiàn)期條件下的洪水多變量分位值估算步驟可以總結(jié)為：

(2) 利用式(4)，計(jì)算所有使C(ud，uq，uv)=p成立的概率組合(ud，uq，uv)；

(3) 由式(15)計(jì)算出使f(d，q，v)達(dá)到最大值的一組(ud，uq，uv)；

(4) 最后分別再根據(jù)邊緣分布的反函數(shù)推求出歷時(shí)、洪峰和洪量的分位值(dml，qml，vml)，dml=F-1(ud)、qml=F-1(uq)、vml=F-1(uv).

4 案例研究

博羅站位于東江下游，是東江流域重要的水文控制站點(diǎn)之一.本文選用博羅站1953-2010年的日流量觀測數(shù)據(jù)作為研究資料，依據(jù)年最大極值理論，首先提取每年的最大日流量作為該年極端洪水事件的洪峰流量，以洪峰流量為基準(zhǔn)，再按流量過程線提取并計(jì)算出該場洪水的歷時(shí)與洪量，總共提取了58場洪水作為博羅站點(diǎn)多變量洪水頻率分析的樣本序列.歷時(shí)與洪峰、洪量樣本兩兩間的Kendall秩相關(guān)系數(shù)表明洪水特征屬性間存在較強(qiáng)的相關(guān)性，其中，洪峰與洪量的相關(guān)性最大(0.665)，歷時(shí)與洪量次之(0.623)，歷時(shí)與洪峰則相對(duì)較小(0.355).

為了檢驗(yàn)GEV和Gumbel-Hougaard copula構(gòu)造洪水事件聯(lián)合分布的合理性，選用水文分析中常用的對(duì)數(shù)正態(tài)分布(GNO)、廣義邏輯斯特分布(GLO)和皮爾遜三型分布(PE3)，Clayton和Frank copula[12]分別對(duì)歷時(shí)、洪峰和洪量樣本進(jìn)行擬合優(yōu)度檢驗(yàn)，并與GEV和Gumbel-Hougaard copula進(jìn)行對(duì)比.選用的分布參數(shù)均采用極大似然法進(jìn)行估計(jì)，擬合優(yōu)度檢驗(yàn)則依據(jù)Akaike信息準(zhǔn)則(AIC)[7]，AIC值越小，代表擬合優(yōu)度越高.選用的備選分布都通過了K-S檢驗(yàn)，表明選用的分布都可用于樣本的擬合.各種分布的擬合優(yōu)度AIC值如表1所示，雖然擬合優(yōu)度檢驗(yàn)的結(jié)果并沒有顯示GEV和Gumbel-Hougaard copula對(duì)所有歷時(shí)、洪峰和洪量樣本以及它們間組合的擬合都是最優(yōu)的(GEV對(duì)歷時(shí)和Gumbel-Hougaard copula對(duì)歷時(shí)與洪量組合的擬合為次優(yōu))，但從總體角度和出于極值分布的考慮，我們選用GEV和Gumbel-Hougaard copula分別作為邊緣分布和多變量分布來構(gòu)造洪水特征屬性間的聯(lián)合分布.圖5和圖6所示的GEV和Gumbel-Hougaard copula對(duì)歷時(shí)、洪峰和洪量樣本的擬合效果亦表明兩者用于構(gòu)造洪水聯(lián)合分布是合理的.

表1 邊緣分布與copula函數(shù)的擬合優(yōu)度AIC值Tab.1 The AIC values for different marginal distributions and copulas

圖5 GEV對(duì)洪水歷時(shí)、洪峰和洪量樣本的擬合優(yōu)度Fig.5 The goodness of fit of GEV for flood duration， peak and volume

圖6 Gumbel-Hougaard copula對(duì)洪水歷時(shí)、洪峰和洪量樣本的擬合優(yōu)度Fig.6 The goodness of fit of Gumbel-Hougaard copula for flood duration， peak and volume

歷時(shí)、洪峰和洪量保證率分別為96%、98%、99%和99.5%時(shí)的單變量重現(xiàn)期、“或”重現(xiàn)期、“且”重現(xiàn)期和二次重現(xiàn)期計(jì)算值如表2所示，在同一保證率水平下(保證率不為0)，二次重現(xiàn)期與“且”重現(xiàn)期相差較小，與“或”重現(xiàn)期則相差較大；保證率越小，二次重現(xiàn)期與“或”重現(xiàn)期和“且”重現(xiàn)期之間的差值則越小，隨著保證率的增大，三者的差值逐漸增大.根據(jù)3.4所述步驟計(jì)算的25、50、100和200年重現(xiàn)期歷時(shí)、洪峰和洪量的多變量設(shè)計(jì)值以及單獨(dú)推算的單變量設(shè)計(jì)值如表3所示，單獨(dú)推算的歷時(shí)與洪量設(shè)計(jì)值均比多變量推算的大，而多變量洪峰設(shè)計(jì)值則要大于單變量設(shè)計(jì)值，并且隨著重現(xiàn)期的增大，兩者的差值逐漸減小.一場洪水造成的效應(yīng)并不是由單個(gè)特征屬性決定的，而是多個(gè)特征屬性共同作用的結(jié)果，如高洪峰伴隨大洪量的洪水造成的影響和破壞要比高洪峰伴隨小洪量的洪水大，長歷時(shí)高洪峰的洪水造成的影響和破壞要比短歷時(shí)高洪峰的洪水大.再者，歷時(shí)、洪峰與洪量之間的組合以及出現(xiàn)的概率又取決于它們間的相關(guān)結(jié)構(gòu)和聯(lián)合分布.因此，相比于單變量設(shè)計(jì)值，考慮了洪水多個(gè)屬性聯(lián)合特征的多變量設(shè)計(jì)值為防洪工程設(shè)計(jì)提供了更加全面和可靠的參考信息.

表2 同一保證率水平下的單變量、“或”、“且”和二次重現(xiàn)期計(jì)算值Tab.2 The calculation values of univariate， “OR， “AND” and secondary return periods in the same probability level

表3 25、50、100和200年重現(xiàn)期歷時(shí)、洪峰和洪量的單變量與多變量設(shè)計(jì)值Tab.3 The design values for flood duration， peak and volume in a multivariate framework under the secondary return periods of 25， 50， 100 and 200 years， respectively

5 結(jié)論

1) 傳統(tǒng)的“或”和“且”多變量重現(xiàn)期對(duì)安全與危險(xiǎn)域的識(shí)別具有局限性，重現(xiàn)期與危險(xiǎn)域范圍大小的矛盾會(huì)造成對(duì)安全事件與危險(xiǎn)事件的錯(cuò)誤識(shí)別.

2) 根據(jù)C(ud，uq，uv)=p臨界條件并利用Kendall函數(shù)定義的二次重現(xiàn)期，具有同一臨界水平下任意歷時(shí)、洪峰與洪量組合對(duì)安全與危險(xiǎn)域識(shí)別的一致性，重現(xiàn)期越大，對(duì)應(yīng)的危險(xiǎn)域則越小，避免了對(duì)安全事件與危險(xiǎn)事件的錯(cuò)誤識(shí)別，更有利于指導(dǎo)洪水風(fēng)險(xiǎn)的管理.

3) 給定二次重現(xiàn)期條件下的洪水事件可以具有不同的歷時(shí)、洪峰與洪量組合，在工程設(shè)計(jì)中，與每個(gè)特征屬性同時(shí)取較大值相比，依據(jù)出現(xiàn)概率最大的原則推算的設(shè)計(jì)值組合可以滿足以較低成本承受較大風(fēng)險(xiǎn)的追求.相比于單變量設(shè)計(jì)值，考慮了洪水多個(gè)屬性聯(lián)合特征的多變量設(shè)計(jì)值為防洪工程設(shè)計(jì)提供了更加全面和可靠的參考信息.

[1] 郭生練，閆寶偉，肖 義等.Copula函數(shù)在多變量水文分析計(jì)算中的應(yīng)用及研究進(jìn)展.水文，2008，28(3)：1-6.

[2] Zhang L， Singh VP. Bivariate flood frequency analysis using the copula method.JournalofHydrologicEngineering， 2006， 11(2)：150-164.

[3] Zhang L， Singh VP. Trivariate flood frequency analysis using the Gumbel-Hougaard copula.JournalofHydrologicEngineering， 2007， 12(4)：431-439.

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

[5] Shiau JT. Fitting drought duration and severity with two-dimensional copulas.WaterResourcesManagement， 2006， 20：795-815.

[6] Salvadori G， De Michele C， Durante F. On the return period and design in a multivariate framework.HydrologyandEarthSystemSciences， 2011， 15：3293-3305.

[7] Corbella S， Stretch DD. Multivariate return periods of sea storms for coastal erosion risk assessment.NaturalHazardsandEarthSystemSciences， 2012， 12：2699-2708.

[8] Salvadori G， Tomasicchio GR， Alessandro FD. Multivariate approach to design coastal and off-shore structures.JournalofCoastalResearch， 2013， 65：386-391.

[9] Yue S， Ouarda TBMJ， Bobée Betal. The Gumbel mixed model for flood frequency analysis.JournalofHydrology， 1999， 226(1)：88-100.

[10] Salvadori G， de Michele C， Kottegoda Netal. Extremes in Nature： An approach using copulas. Dordrecht： Springer， 2007：15-30.

[11] Favre AC， Adlouni SE， Perreault Letal. Multivariate hydrological frequency analysis using copulas.WaterResourcesResearch， 2004， 40： W01101. doi：10.1029/2003WR002456.

[12] Nelson RB. An introduction to copulas： Second edition. New York： Springer， 2006：115-143.

[13] Salvadori G， de Michele C. Multivariate multiparameter extreme value models and return periods： A copula approach.WaterResourceResearch， 2010， 46： W10501. doi：10.1029/2009WR009040.

[14] Salvadori G， de Michele C. Frequency analysis via copulas： Theoretical aspects and applications to hydrological events.WaterResourcesResearch， 2004， 40： W12511. doi：10.1029/2004WR003133.

[15] Chebana F， Ouarda TBMJ. Multivariate quantiles in hydrological frequency analysis.Environmetrics， 2011， 22：63-78.

[16] Genest C， Quessy JF， Rémillard B. Goodness-of-fit procedures for copula model based on the probability integral transformation.ScandinavianJournalofStatistics， 2006， 33(2)：337-366.

[17] Graler B， van den Berg MJ， Vandenberghe Setal. Multivariate return periods in hydrology： a critical and practical review focusing on synthetic design hydrograph estimation.HydrologyandEarthSystemSciences， 2013， 17： 1281-1296.

Multivariate flood risk assessment based on the secondary return period

HUANG Qiang & CHEN Zishen

(DepartmentofWaterResourcesandEnvironment，SunYat-senUniversity，Guangzhou510275，P.R.China)

Frequency analysis is a useful tool for flood risk assessment， but the definition and calculation of return periods and design values in a multivariate framework are difficult tasks. In this paper， by constructing the joint distribution of flood duration， peak discharge and volume， we introduced the definitions of multivariate “OR”， “AND” and secondary return periods and discuss their differences in identifying the safe and dangerous regions in a critical level. The rationality and reliability of information provided from them in flood risk management and engineering design were then analyzed， respectively. The traditional “OR-AND” approach in multivariate return periods definition is limited in the identification of safe and dangerous regions， while the secondary return period based on the Kendall measure would be more rational. Consequently， mistakes in the recognition of safe and dangerous event can be avoided by the use of second return period， which， obviously， is better for the flood risk management. For a given second return period， the combination of flood duration， peak discharge and volume design values calculated from the most-likely approach， may have an advantage to meet the requirement that with low cost to bear greater risks in engineering design. Furthermore， the information provided from multivariate design values would be more considerate and reliable than the univariate ones.

Multivariate flood characteristics； extreme value distribution； safe and dangerous regions； Kendall function； secondary return period； multivariate design values

*國家自然科學(xué)基金項(xiàng)目(41371498)和廣東水利創(chuàng)新基金項(xiàng)目(2009-41)聯(lián)合資助.2014-02-27收稿；2014-08-11收修改稿.黃強(qiáng)(1989～)，男，博士研究生；E-mail：huangq52@mail2.sysu.edu.cn.

**通信作者；E-mail：eesczs@mail.sysu.edu.cn.