陳子燊，劉占明，趙青

(1.中山大學地理學院，廣東 廣州 510275；2.佛山科學技術(shù)學院，廣東 佛山 528000；3.中山市南區(qū)水利所，廣東 中山 528455)

政府間氣候變化專門委員會(IPCC)第5次評估報告預測氣候變化會增大包括洪澇在內(nèi)的各種水文極端事件發(fā)生的概率，此對江河堤圍和水庫大壩等水利工程設(shè)施將造成嚴重威脅。如何評估變化環(huán)境下水文極端事件的風險并提出相應的設(shè)計標準是應對洪澇災害的重要科學與工程應用問題。洪水過程是由洪峰流量、洪水總量和洪水歷時等多個特征要素有機組成的一個整體，至今對多變量洪水頻率分析已經(jīng)有較多的研究[1-5]，文獻[6]對此作了述評。李天元等[7]認為目前大多數(shù)研究還限于在構(gòu)建多變量洪水聯(lián)合分布的基礎(chǔ)上進行聯(lián)合重現(xiàn)期和條件概率的分析。郭生練等[8]進一步指出，對于給定的重現(xiàn)期水平，存在無窮多種滿足防洪標準的洪峰、洪水總量組合，如何合理地選擇設(shè)計值是關(guān)鍵問題，如何在多變量框架下進行洪水聯(lián)合設(shè)計值估算和風險評估仍然是一個有爭議的問題，應加強這方面的研究。

至今，“或”和“且”重現(xiàn)期是最常用的兩種多變量重現(xiàn)期定義方法，但在安全與危險事件的判定上兩者都存在著較大局限性[9]。Salvadori等[10]針對“或”重現(xiàn)期的不足提出了劃分安全與危險臨界域的新多變量重現(xiàn)期——Kendall重現(xiàn)期(Kendall return periods，KRP。又稱二次重現(xiàn)期，secondary return periods)，其相關(guān)原理已在海岸工程設(shè)計研究中得到了初步應用[11-12]。隨后，進一步針對“且”重現(xiàn)期的不足，基于生存Kendall分布函數(shù)，提出了劃分安全與危險的臨界域的新多變量重現(xiàn)期—生存kendall重現(xiàn)期(survival Kendall return periods，SKRP)[13]。本文擬通過實例對比分析 “或”首次重現(xiàn)期和二次重現(xiàn)期的聯(lián)合設(shè)計水平、“且”首次重現(xiàn)期和生存Kendall重現(xiàn)期的聯(lián)合設(shè)計水平之間的差異，希翼有助于深化對洪水風險的認識，為防洪設(shè)計標準提供新的參考。

1　理論與方法

1.1　Copula函數(shù)與首次重現(xiàn)期

根據(jù)Sklar定理，若F(·)是一個二維隨機變量(X,Y)的累積分布函數(shù)，其邊緣分布函數(shù)是連續(xù)函數(shù)u=FX(x),v=FY(y)，則有唯一的Copula函數(shù)C使得：

F(x,y)=P(X≤x,Y≤y)=

C(FX(x),FY(y))=C(u,v)

(1)

(2)

(3)

“或”重現(xiàn)期和“且”重現(xiàn)期也統(tǒng)稱為首次重現(xiàn)期(primary return periods)。

1.2　Kendall分布函數(shù)與Kendall重現(xiàn)期

從首次重現(xiàn)期可知，不同的u、v組合只要其出現(xiàn)累積概率(記為t)相同都可產(chǎn)生相同的重現(xiàn)期。針對“或”首次重現(xiàn)期定義的安全域/危險域存在誤判的問題，Salvadori等[10]利用Nelsen[14]定義的Kendall分布函數(shù)劃分出亞臨界(安全域)、臨界(警戒事件)和超臨界(危險域)3種情景。通過判定累積概率是否小于或等于某臨界概率t，將多維的極值事件投射為一維分布。則基于Copula函數(shù)累積概率為t的(u,v)組合值，Kendall分布函數(shù)KC為[15]：

KC(t)=t-φ(t)/φ′(t),0

(4)

式中，φ′(t)為φ(t)的右導數(shù)。由Kendall分布函數(shù)確定的重現(xiàn)期即Kendall重現(xiàn)期為：

(5)

1.3　生存Kendall函數(shù)與生存Kendall重現(xiàn)期

(6)

定義TSK為生存Kendall重現(xiàn)期(survival Kendall return period，SKRP)：

(7)

1.3　聯(lián)合分布設(shè)計值

多變量聯(lián)合分布同一個重現(xiàn)期可以有不同的分位數(shù)組合與之對應，這些分位數(shù)組合無法通過概率分布的反函數(shù)直接計算。具有相同重現(xiàn)期Tp(C(u,v)=p)的分位值組合構(gòu)成了一個二維點集(等曲線)，設(shè)計分位數(shù)組合必然存在一個使聯(lián)合概率密度f(u,v)達到最大值的組合(um,vm)，即該組合出現(xiàn)的可能性最大。因此，在設(shè)定重現(xiàn)期條件下出現(xiàn)最大可能組合的設(shè)計值可作為工程設(shè)計與風險管控標準的合理選擇。

(8)

f(u,v)=c(u,v)f(u)f(v)

(9)

式中，c為二維Copula的概率密度函數(shù)。

2　示例研究

2.1　基本數(shù)據(jù)

位于珠江三角洲河口區(qū)的增江是珠江流域三大干流之一東江的一級支流，河流全長203 km, 流域面積約3 160 km2。1954年設(shè)立的麒麟咀水文站是增江流量控制站和國家重要水文站，距河口距離33 km，集水面積2 866 km2。本研究采用麒麟咀水文站1954-2011年共58 a的逐日流量數(shù)據(jù)，按年最大值提取此水文站歷年洪水的洪峰流量Q與同場洪水總量W作為研究樣本。

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

采用水文頻率分析中常用的3種三參數(shù)概率分布：皮爾遜三型分布(P-III)、廣義極值分布(GEV)、對數(shù)正態(tài)分布(GNO)分別擬合洪峰流量、洪水總量和洪峰水位樣本。參數(shù)估計使用線性矩(L-矩)方法。經(jīng)驗頻率分布使用Gringorten公式。擬合結(jié)果采用均方根誤差(RMSE)、和概率點據(jù)相關(guān)系數(shù)(PPCC)檢驗其擬合優(yōu)度。如表1所示，洪峰流量序列選用GEV分布, 洪量序列選用P-III分布。

表1　洪峰和洪量的概率分布參數(shù)與擬合優(yōu)度檢驗值Table 1　The parameters of the marginal distribution and the values of goodness of fit test

洪峰Q和洪量W的Kendall相關(guān)系數(shù)為0.565。采用相關(guān)性指標法計算Q-W聯(lián)合分布的Gumbel-Hougaard copula、Frank copula、A-M-H copula和Clayton copula的參數(shù)θ及相應的AIC和OLS值，計算結(jié)果見表2。以O(shè)LS和AIC值最小，擬合度最高的二維Gumbel -Hougaard copula構(gòu)建的洪水峰量聯(lián)合分布的Copula模式如下：

C(FQ(q),FW(w))=

exp{-[(-lnFQ(q))2.299+

(-lnFW(w))2.299]1/2.299}

(10)

表2洪水峰量聯(lián)合分布的參數(shù)估計與擬合優(yōu)度評價結(jié)果
Table 2Results of parameter estimation and goodness-of-fit tests for four Copulas

ArchimedeanCopulaθOLSAICClayton2.5980.031-173A-M-H0.9900.047-152Gumbel-Hougaard2.2990.025-183Frank7.0550.028-178

2.3　條件概率分布

分析某特定洪峰流量條件下，出現(xiàn)洪水總量的概率分布。當給定Q≥q時，W≥w的條件概率為：

(11)

條件重現(xiàn)期為：

(12)

從表3可看出：① 當麒麟咀洪峰流量大于等于某一特定頻率設(shè)計值時，洪水總量出現(xiàn)大于等于該頻率設(shè)計值的條件概率隨著超值概率的減小而減小，以Q出現(xiàn)大于等于概率P1%時的設(shè)計洪峰流量Q=3 635 m3/s為例，洪水總量出現(xiàn)大于等于P10%、P5%、P2%、P1%、P0.5%和P0.2%的條件概率分別為0.981、0.950、0.836、0.651、0.417和0.189，其相應的條件重現(xiàn)期為：1.0、1.1、1.2、1.5、2.4和5.3 a。② 反之，洪量大于等于某一特定頻率設(shè)計值時，洪峰流量出現(xiàn)大于等于該頻率設(shè)計值的條件概率隨著超值概率的減小而增大。如，百年一遇的設(shè)計洪量W=1 702×106m3，其與洪峰流量的遭遇隨超值概率的減小而增大，洪峰流量出現(xiàn)大于等于P10%、P5%、P2%、P1%、P0.5%和P0.2%的條件概率分別為0.098、0.190、0.418、0.651、0.834、0.947，其相應的條件重現(xiàn)期為：10.2、5.3、2.4、1.5、1.2和1.1 a。此外，洪水峰量同頻率遭遇的概率超過64.9%，表明二者遭遇的概率非常高。上述分析表明，可能存在著滿足多種防洪標準的洪水峰量組合，需要合理地推算其設(shè)計值。

表3Q-W洪水組合條件概率
Table 3Conditional probabilities ofQ-W

P(W≥w|Q≥q)或P(Q≥q|W≥w)W/Q2396278532763635398544369600.6720.8510.9540.9810.9920.99811890.4250.6600.8800.9500.9800.99414830.1910.3520.6530.8360.9310.97917020.0980.1900.4180.6510.8340.94719180.0500.0980.2330.4170.6490.87322000.0200.0400.0980.1890.3490.649

圖1　P(W≥w|Q≥q)條件概率分布圖Fig.1　Diagram of conditional probability distributions of joint distribution between flood discharges and volumes

2.4　聯(lián)合分布重現(xiàn)期和危險率

2.5　洪水設(shè)計值

按設(shè)定重現(xiàn)期(500、200、100、50、20和10 a)推算麒麟咀測站的洪峰流量和洪水總量單變量設(shè)計值，以出現(xiàn)概率最大的原理推算洪水峰量聯(lián)合分布的“或”重現(xiàn)期、“且”重現(xiàn)期、Kendall 重現(xiàn)期和生存Kendall 重現(xiàn)期設(shè)計值列于表5。結(jié)果顯示，對于設(shè)定的10 ～ 500 a重現(xiàn)期，按Kendall重現(xiàn)期推算的設(shè)計洪峰流量和設(shè)計洪水總量分別小于“或” 重現(xiàn)期設(shè)計值和邊緣分布設(shè)計值，按Kendall重現(xiàn)期推算的洪峰流量和設(shè)計洪水總量設(shè)計值和相應邊緣分布設(shè)計值的相對誤差分別為-2.8% ～ -5.4%和-3.8% ～ -8.1%；推算的生存Kendall重現(xiàn)期設(shè)計洪峰流量和設(shè)計洪水總量分別大于“且”重現(xiàn)期設(shè)計洪水值和邊緣分布設(shè)計值，生存Kendall重現(xiàn)期設(shè)計洪水值和相應的邊緣分布設(shè)計值二者的相對誤差分別為2.4% ～ 3.8%和0.3% ～ 4.8%。這一結(jié)果也表明，按目前有關(guān)規(guī)范設(shè)計要求的單變量洪水要素設(shè)計值已達到設(shè)計洪水安全標準，按洪水峰量聯(lián)合分布的“或”重現(xiàn)期推算的洪水設(shè)計值存在高估問題，導致投入費用偏大，按洪水峰量聯(lián)合分布“且”重現(xiàn)期推算的洪水設(shè)計值存在低估問題，可能增大防洪工程損毀風險。按洪水峰量聯(lián)合分布的Kendall重現(xiàn)期或生存Kendall重現(xiàn)期推算的洪水設(shè)計值可為防洪工程安全與設(shè)計標準提供更堅實的理論依據(jù)。為了比較，進一步推算洪峰—洪量同頻率分布設(shè)計值[16]：

表5　不同重現(xiàn)期下洪峰流量與洪水總量的設(shè)計值Table 5　Design values of flood discharges and flood volumes at different return periods

u1=u2=[1-(1/Tu1,u2)]α;

(13)

3　結(jié)　論

本文對比分析了增江流域麒麟咀斷面洪峰洪量之間的二維聯(lián)合分布及其4種重現(xiàn)水平，有以下結(jié)論：

1)對于預設(shè)的重現(xiàn)期，相對于“或”聯(lián)合重現(xiàn)期，Kendall重現(xiàn)期可更準確地識別單個超臨界區(qū)域的風險率；相對于“且”聯(lián)合重現(xiàn)期，生存Kendall重現(xiàn)期可更好地描述洪峰洪量同時超閾值情況下的風險率。

2)按目前有關(guān)規(guī)范設(shè)計要求的單變量洪水要素設(shè)計值基本達到設(shè)計標準，按洪峰洪量“或”重現(xiàn)期和同頻率推算的設(shè)計洪水分位數(shù)偏高，工程費用偏大，采用Kendall重現(xiàn)期和生存Kendall重現(xiàn)期推算的設(shè)計洪水值可為防洪工程風險管理與設(shè)計提供新的選擇與參考。

參考文獻：

[1]馮平,毛慧慧,王勇.多變量情況下的水文頻率分析方法及其應用[J].水利學報,2009,40(1):33-37.

FENG Ping, MAO Huihui, WANG Yong. Method for hydrological reoccurrence frequency analysis under the condition of multivariate[J]. Journal of Hydraulic Engineering, 2009,40(1):33-37 .

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

[3]熊立華,郭生練,肖義,等. Copula 聯(lián)結(jié)函數(shù)在多變量水文頻率分析中的應用[J]. 武漢大學學報(工學版), 2005,38(6):16-19.

XIONG Lihua, GUO Shenglian, XIAO Yi, et al. Application of Copulas to multivariate hydrological frequency analysis[J]. Engineering Journal of Wuhan University, 2005, 38(6):16-19.

[3]方彬, 郭生練, 肖義, 等. 年最大洪水兩變量聯(lián)合分布研究[J]. 水科學進展, 2008, 19( 4) : 505-511.

FANG Bing, GUO Shenglian, XIAO Yi, et al. Annual maximum flood occurrence dates and magnitudes frequency analysis based on bivariate joint distribution[J]. Advances in Water Science，2008，19(4) : 505-511.

[5]陳子燊, 黃強, 劉曾美. 基于非對稱Archimedean copula與kendall分布函數(shù)的三變量洪水風險評估[J]. 水科學進展, 2016, 27(5):763-771.

CHEN Zhisheng, HUANG Qiang, LIU Zengmei. Risk assessment of trivariate floods based on asymmetric copula and kendall distribution function[J]. Advances in Water Science，2016, 27(5):763-771.

[6]郭生練,閆寶偉,肖義,等. Copula 函數(shù)在多變量水文分析計算中的應用及研究進展[J]. 水文,2008,28(3):1-7.

GUO Sheng lian, YAN Baowei, XIAO Yi, et al. Multivariate hydrological analysis and estimation[J]. Journal of China Hydrology, 2008,28(3):1-7.

[7]李天元,郭生練,劉章君,等. 基于峰量聯(lián)合分布推求設(shè)計洪水[J].水利學報,2014,45(3):269-276.

LI Tianyuan，GUO Shenglian，LIU Zhangjun，et al. Design flood estimation based on bivariate joint distribution of flood peak and volume[J]. Journal of Hydrologic Engineering,2014,45(3):269-276.

[8]郭生練,劉章君,熊立華.設(shè)計洪水計算方法研究進展與評價[J].水利學報,2016,47(3):302-314.

GUO Shenlian, LIU Zhangjun, XIONG Lihua. Advances and assessment on design flood estimation methods[J]. Journal of Hydraulic Engineering,2016,47(3):302-314.

[9]SALVADORI G, de MICHELE C. Frequency analysis via copulas: theoretical aspects and applications to hydrological events[J]. Water Resources Research, 2004, 40(12):229-244.

[10]SALVADORI G, MICHELE D C, DURANTE F. On the return period and design in a multivariate framework[J]. Hydrology and Earth System Sciences, 2011, 15: 3293-3305.

[11]CORBELLA S, STRETCH D D. Multivariate return periods of sea storms for coastal erosion risk assessment[J]. Natural Hazards and Earth System Sciences, 2012, 12: 2699-2708.

[12]SALVADORI G, TOMASICCHIO G R, ALESSANDRO F D. Multivariate approach to design coastal and off-shore structures[J]. Journal of Coastal Research, 2013, 65: 386-391.

[13]SALVADORI G, DURANTE F, de MICHELE C. Multivariate return period calculation via survival functions[J]. Water Resources Research, 2013, 49: 2308-2311.

[14]NELSEN R B. An introduction to copulas (springer series in statistics) [M]. New York： Springer,2006：216.

[15]GRALER B, VANDENBERG M J, VANDENBERGHE S, et al. Multivariate return periods in hydrology: a critical and practical review focusing on synthetic design hydrograph estimation[J]. Hydrology and Earth System Sciences, 2013, 17(4): 1281-1296.

[16]肖義, 郭生練, 劉攀, 等. 分期設(shè)計洪水頻率與防洪標準關(guān)系研究[J].水科學進展, 2008, 19(1) : 54- 60.

XIAO Yi, GUO Shenlian, LIU Pan, et al. Seasonal flood frequency analysis and flood prevention standard[J]. Advances in Water Science,2008,19(1) :54-60.