基于Copula函数的铁山港雨潮遭遇组合分析

2022-06-25甘富万邵帅杰黄宇明胡秀英

中国农村水利水电 2022年6期

甘富万，邵帅杰，黄宇明，胡秀英，刘旻，陈帅

（1.广西大学土木建筑工程学院，南宁 530004；2.广西防灾减灾与工程安全重点实验室，南宁 530004；3.南京水利科学研究院港口航道泥沙工程交通行业重点实验室，南京 210029；4.南宁学院，南宁 530004）

0 引言

洪、潮组合问题一直是沿海地区研究的热点问题。近年来，Copula 函数广泛应用于水文气象领域，对于具有任意边缘分布、任意相关性结构的多个变量，Copula 函数可以方便地构造出其联合概率分布，并用于分析多变量水文事件。涂新军等［1］基于Copula方法模拟24 h暴雨遭遇日高潮位的联合分布特征，对比雨潮遭遇传统重现期和二次重现期差异，根据同频法和权函数法反推计算雨潮设计组合值，为滨海地区防洪潮设计提供参考。杨星等［2］采用Clayton Copula 构建了深圳市洪潮遭遇组合的风险分析模型，得到了深圳市洪潮组合的风险率，可用于深圳地区河道防洪治理。刘曾美等［3］采用Copula函数构建感潮河段的年最大洪水流量和相应潮位的联合分布以及年最高潮位和相应洪水流量的联合分布，可为感潮河段洪潮遭遇组合的合理选取提供科学依据。

鉴于此，本文基于Copula 函数，构建铁山港海湾暴雨与石头埠站潮位的联合概率分布，引入二次重新期概念，在雨潮组合风险分析的基础上，为铁山港海湾防洪治涝规划提供科学依据。

1 雨潮遭遇组合分析方法

1.1 Copula函数基本理论

Copula 函数是定义域为［0，1］均匀分布的多维联合分布函数，它可以通过连接多个随机变量的边缘分布来构造联合分布。令H为联合分布函数，F和G为其边缘累积分布函数，那么存在唯一的Copula函数C使得对∀x，y∈有：

如果F和G是连续的，则C是唯一的。否则，C在RanF×RanG（Ran 表示值域）唯一确定。相反地，如果C是一个Copula函数，F和G是边缘累积分布函数，则上式定义H 是具有边缘累积分布F和G的一个二维联合累积分布函数。

现阶段水文气象领域的Copula 函数主要有椭圆Copula 函数［4］、Plackett Copula 函数［5］以及阿基米德Copula 函数［6，7］，相较于其他两种，阿基米德Copula 函数具有构造简单且相关参数θ计算较为简便等特点，使得阿基米德Copula 函数应用最为广泛［8，9］。本文选取阿基米德Copula 函数中应用最为广泛的Frank、Clayton、GH、AMH Copula、No.16 Copula 以及No.17 Copula 函数为候选函数，并对其函数模型进行检验，选择最合适的Copula函数模型。

1.2 雨潮遭遇设计组合

在流域防洪规划中，通常需要确定雨潮组合方式，即在给定流域暴雨的条件下，确定河口潮位取值情况。但在实际工程中，当流域暴雨一定时，河口潮位取值并非唯一，并且在实际计算中不同的设计组合值所形成的防雨、潮能力也不尽相同，在现阶段针对多变量水文事件无法确定设计组合的情况，一般是取特殊的设计组合如同频组合、最可能组合等进行分析计算。

雨、潮组合的同频组合是指在计算中令雨、潮的遭遇事件同频发生，即令其边缘分布u=v；另外一种组合就是最可能组合，其表达式如下，

式中：c（u，v）为雨潮联合分布的概率密度函数；f（x）和f（y）则是其边缘分布的概率密度函数。

1.3 雨潮组合风险

为了定量分析不同的雨潮组合发生的概率，本文构建两种不同的风险概率模型，用于评价洪潮遭遇风险：①同现风险率。即降雨和潮位同时超过某一设计标准的概率，该情况是最不利，危险最大的情况。②条件风险率。即当流域发生超过某一标准暴雨时，潮位也发生超过某一标准潮位的概率。设P表示流域暴雨，Z表示潮位，借助Copula 函数，同现风险率和条件风险率的计算公式分别为：

式中：u为P的边缘分布函数；v为Z的边缘分布函数。

2 实例分析

铁山港地区地势较为低洼，容易受到暴雨洪水以及风暴潮的影响。以铁山港工业区为对象的研究区域的汇水面积为19.41 km2，属于小流域汇水区域。现阶段，针对雨、潮组合的水文事件研究选样一般有两种方法：①以雨为主的24 h 暴雨遭遇日高潮位；②以潮为主的日高潮位遭遇日降水。参考年最大值法［1］，选择北海站1967-2008年的年最大24 h 降雨序列作为本次研究的雨量资料，再根据北海站年最大24 h 降雨事件，选择对应时段石头埠站石头埠站1967-2008年的42年日最高潮位序列，作为雨、潮事件进行分析计算。

2.1 雨潮联合分布构建

根据北海站连续42 a 最大24 h 雨资料序列，以及石头埠潮位站对应的日最高潮位的Cs和Cv值，计算雨量资料以及潮位资料各自的四种边缘分布函数：皮尔逊III 型（PE3）、广义极值分布（GEV）、广义Pareto 分布（GPD）和广义Logistic 分布（GLO）的相关参数，并采用RMSE 准则和AIC 准则对雨量、潮位的边缘分布函数进行优选。计算结果表明雨量资料序列的最优边缘分布模型为皮尔逊III 型（PE3），而潮位资料序列的最优边缘分布模型则为广义Pareto 分布（GPD）函数模型，因此分别选取皮尔逊III 型和广义Pareto 分布雨量以及相应潮位的边缘分布模型，并进行下一步计算。

采用Kendall 秩相关系数方法来计算不同Copula 函数模型的相关参数得到北海站历年最大24 h 降雨量与石头埠站相应的潮位资料之间的Kendall秩相关系数τ的值为0.217 2，证明雨量资料与潮位资料之间存在一定的相关性。

计算Frank Copula、Clayton Copula、GH Copula、AMH Copula、No.16 Copula 以及No.17 Copula 的相关参数，计算结果见表2。

表2 不同Copula函数的参数计算结果Tab.2 The result of the parameter calculation of different Copula functions

通过对比不同Copula 函数参数的取值范围，得到上述6 种Copula 函数的参数计算结果均在合理的范围，可以进行下一步的计算。

选择Frank、Clayton、GH、AMH Copula、No.16 Copula 以及No.17 Copula 函数作为候选函数，采用图形评价法，AIC 准则［10，11］，BIC准则［12］和OLS准则［13］进行函数模型评价。

表1 雨、潮边缘分布计算成果Tab.1 Calculated results of rain and tide margin distribution

图形评价法的计算结果如图1。

图1 雨量与潮位联合分布理论频率与经验频率拟合结果Fig.1 Fitting results of theoretical frequency and empirical frequency of joint distribution of rainfall and tide level

由图1可以看出不同的联合分布模型均比较分散地分布在45°线的两侧，同时根据理论频率与经验频率的散点分布情况，并不能直观明显地判断出哪个联合分布模型是最优的雨、潮联合模型，因此需要通过分析其他优度评价结果才可以选择出最合适Copula函数模型。

表3计算结果显示，6 种不同的雨、潮联合分布模型的分别对应的三种函数优度评价计算结果差异不大，均比较接近。通过对比，可以发现GH Copula函数的AIC值，BIC值以及RMSE值均为6 个Copula 函数中的最小值，根据函数优度评价方法——评价准则计算结果越小，函数拟合优度越好的依据，证明了GH Copula 函数是本研究中雨、潮联合分布的最优模型。根据上述结果，雨、潮的联合分布模型为：

表3 Copula函数联合分布模型优度评价结果Tab.3 Results of the evaluation of the superiority of the joint distribution model of Copula function

根据雨、潮联合分布计算结果，采用最优的GH Copula 函数作为雨、潮联合分布模型，并分析雨、潮联合重现期以及同现重现期下的同频组合以及最可能组合等的设计组合值。

分别采用变量X和Y来表示北海站的雨量资料和石头埠站潮位资料序列，以GH Copula 函数构建雨、潮组合的联合分布模型如图2。

图2 雨量与潮位资料联合概率分布模型Fig.2 Joint probability distribution model of rainfall and tide level data

2.2 雨潮组合设计重现期

现阶段对于雨、潮组合的重现期的研究关注较多的是同现重现期以及二次重现期。在单变量条件下，雨量以及潮位超过某一个阈值的计算可以表示如下：

雨、潮同现重现期（“And”重现期）表示为降雨量以及潮位同时超过某一个特定设计阈值的情况。联合Copula 函数，则根据水文重现期概率计算方法，雨、潮组合的同现重现期可以表示如下：

雨、潮联合重现期（“OR”重现期）表示为降雨量或潮位其中一个超过某一个特定设计阈值的情况。联合Copula 函数，则根据水文重现期概率计算方法，雨、潮组合的联合重现期可以表示如下：

近年来一些专家学者认为传统重现期并不是很好地描述多变量危险事件，于是提出了二次重现期［14-17］。生存Kendall重现期［14］是在同现重现期的基础上，以雨、潮联合分布概率相等作为的一系列事件作为临界事件，并以此来划分雨、潮联合分布的危险域和安全域。一般情况下，采用基于C（u，v）的Kendall分布函数KC作为生存Kendall 重现期对应的超阈值的联合分布概率函数［15］，因此可以表示为：

式中，q∈(0，1)。KC的分布函数如下：

则对于雨、潮组合来说，其生存Kendall重现期可以表示为：

式中：NE是暴雨事件年平均发生次数，在本文研究中，取历年最大24 h 降雨量及其相对时刻的潮位作为雨、潮序列组合，因此NE取值为1。

同时，针对联合（“OR”）重现期在描述多变量水文事件时存在局限，Salvadori［16，17］等在联合重现期的基础上提出采用Kendall重现期来计算多变量水文研究。即：

式中：φ(p)是阿基米德Copula 函数的生成元；φ′(p)是生成元函数的导数。

Kendall重现期可以定义如下，

2.3 雨潮遭遇设计组合分析

对传统重现期和二次重现期的雨、潮同频组合以及最可能组合进行计算，从表4的计算结果可以看出，同频组合条件下或者最可能组合的二次重现期对应的100、50、20、10、5以及2 a重现期水平下的雨、潮组合均明显大于传统重现期所对应的各重现期下的组合。同现重现期条件下，最可能组合的各设计重现期水平年对应的潮位的重现期范围为91.97～1.60 a，而同频组合下的相应的重现期范围是30.03～1.47 a，且最可能组合对应的各重现期条件下的潮位值均大于同频组合的设计潮位值；同时，可以看到最可能组合条件的各重现期水平年对应的雨量设计的重现期范围为1.85～1.34 a，均明显小于同频组合下雨量设计重现期。生存Kendall重现期条件下，最可能组合的各设计重现期水平年对应的潮位的重现期范围为316.26～3.06 a，而同频组合下的相应的重现期范围是95.24～2.27 a，生存Kendall重现期条件下不同设计水平年的最可能组合的潮位设计值均明显大于同频组合条件的潮位各水平年设计值。由此可以推断在相同的雨量、潮位组合对应下的同现重现期大于生存Kendall重现期，也就是说在相同组合条件下，生存Kendall 重现期要比同现重现期相对来说更加安全。

表4 雨、潮遭遇组合不同重现期设计值Tab.4 Different reproducibility design values for rain and tide encounter combinations

联合重现期同频组合的计算重现期为172.41～3.02 a，相同条件下Kendall 重现期的变化范围为75.19～1.67 a。相同的设计重现期下，联合重现期的各设计值均比Kendall 重现期的雨、潮设计组合要大。即相同的雨、潮设计组合，采用Kendall 重现期的设计标准要比联合重现期高，说明Kendall重现期比联合重现期更加安全。

进一步计算给定设计雨量条件下的各设计水平年条件下的条件风险概率以及同现风险概率，如表5所示。

表5 雨、潮组合同现风险以及条件风险概率Tab.5 Rain，tide combination co occurrence risk and conditional risk probability

表5表明，雨、潮不同设计重现期水平年组合的同现风险概率均比较小，且雨量（或潮位）遭遇重现期越大，遭遇潮位（或雨量）越大的同现风险概率越小；而当雨、潮遭遇均较小的时，雨潮同现风险概率则相对较大，这说明雨、潮遭遇组合在重现期较大的水平下，其相关性不高，但当雨、潮遭遇水平较小的情况下，其相关性相对较大。从条件风险概率来分析，也可以看出当雨量设计重现期较大时，遭遇较大重现期水平的潮位也较大时的条件概率较小，其中当雨量设计重现期为100 a 时，潮位遭遇也为100 a的条件风险概率为0.286；在相同条件下，潮位遭遇为2 a时的条件风险概率则为0.879。综上，可以得到雨、潮组合各自遭遇重现期越小，其同现风险概率越高；而当雨量遭遇重现期越大的条件下，其遭遇潮位重现期越小的事件的条件概率风险越大。