基于Copula函数的水库强风暴雨遭遇分析研究

蔡奇霖，徐丹丽，王兆礼，赖成光*

(1.华南理工大学土木与交通学院，广东 广州 510640；2.广东省水利电力勘测设计研究院有限公司，广东 广州 510635)

在全球气候变暖的背景下，不少地区的台风与暴雨频度、规模、时空特征等均发生了较大变化[1-3]，这些变化给地区水库设计、建设、管理等带来诸多挑战。在水库大坝设计过程中，确定合理的设计风速对雍水及风浪爬高计算至关重要，决定了大坝的综合超高以及安全裕度。然而，当前中国大部分技术规范仅考虑单一的风速变量确定设计风速，鲜有定量考虑强风和暴雨遭遇的情况，如SL 274—2020《碾压式土石坝设计规范》和SL 25—2006《砌石坝设计规范》是按大坝的重要性给出相应的风速标准，没有考虑风速与暴雨的遭遇。尽管SL189—2013《小型水利水电工程碾压式土石坝设计规范》要求在沿海地区应该考虑最大风力与暴雨同时出现的条件，但如何执行也没有明确规定。广东省是中国受台风和暴雨影响最严重的沿海地区之一，双台风甚至是三台风同时出现的现象也并非历史罕见。已有研究表明该地区登陆台风呈现出频数减少但强度明显增强的特点[4-6]，而降雨则表现出强降水增加，弱降水减少的特点[7-9]。在降雨和风速均发生改变的背景下，当二者发生遭遇时，按照目前有关规范计算的坝顶超高(特别是沿海地区的水库)是否仍能满足安全要求，这是一个值得深入研究的课题。

在众多的多变量分析方法中，Copula函数因其不受边缘分布限制的优点，已在水文频率分析中得到了广泛应用[10-11]。徐明等[12]利用Copula函数针对上海地区建立了日降雨量极值与极大风速极值的风雨联合分布模型，结果表明联合分布模型得到的风雨同现超越概率能较好地反映风险程度。侯静惟等[13]以海南省计值风速和总降水量为研究样本，建立风雨联合分布模型，并分析了“或”和“且”2种重现期的风雨重现水平。Dong等[14]针对崇明岛地区基于Copula函数建立了风速和雨强的联合分布模型，并通过了实测数据检验。徐宗学等[15]基于Copula函数对深圳市洪潮组合进行风险分析，结果表明降水与潮位的双阈值组合风险率小于单阈值组合风险率，设计时应更加关注降水和潮位的双阈值组合风险率。马川惠等[16]基于Copula函数针对泾河流域建立水沙联合分布模型得到了两变量联合设计值的最可能模式。李天元等[17]以隔河岩水库为例，基于Copula函数构建洪峰、洪量间的二位联合分布，推导了两变量同频率组合和条件期望组合两种联合设计值组合。林娴等[18]以武江流域实测资料为基础，基于Copula函数构建了洪峰、洪量与洪水历时间的联合概率分布，结果表明同频率下联合分布的洪水设计值相对于单变量设计值偏安全。曾明等[19]基于上海徐家汇站的逐日降雨资料，采用Copula函数构建了年最大1 d与3 d降水量的联合分布模型，结果表明考虑最大1 d和3 d降水量的遭遇组合，有利于提升防洪排涝能力。综上可知，Copula函数已被广泛用来构建联合分布连接函数，能够有效地描述水文事件的内在规律和特征属性之间的相互关系。

尽管目前关于风速和降雨的联合分布研究已有较多成果，但关于广东沿海地区大型水库风速与降雨的联合分布研究尚不多见。基于此，本文以广东省台山大隆洞水库为研究案例，基于Archimedean Copula函数和Kendall分布函数分别构建了年最大1 d降雨量、年最大3 d降雨量、年最大7 d降雨量及其对应时段最大风速的二维联合分布函数，推求不同时段最大降雨量以及最大风速在“或”“且”和Kendall重现期的重现期水平，并对比分析3种遭遇情景结果间的差异，最后基于遭遇分析的风速结果与当前技术规范进行对比，并提出相关建议。研究结果以期为沿海地区水库大坝的设计、建设以及管理提供科学的参考依据。

1 水库简介

大隆洞水库位于广东省台山市东南部，是广东省受台风、暴雨影响最为频繁的大型水库之一。由大隆洞水库大坝基本资料登记表(BHF24000112)和水库管理处有关资料获知，该水库于1959年建成，2003年经加固改造后水库集雨面积达148 km2，水库总库容2.96亿m3，为大(2)型水库，设计重现期为100年一遇。水库主坝坝型为均质土坝，坝顶高程(不含防浪墙)39.2 m，坝顶长度为441 m，坝顶宽度6 m，坝基采用混凝土防渗墙防渗；水库正常蓄水位30.80 m，设计洪水位35.42 m，校核洪水位38.11 m。

降雨和风速数据来源于水库附近的台山站(112°78′E、22°25′N)，其中降雨(逐日)数据序列长度为1960—2018年，最大风速(逐日)数据序列长度为1973—2018年。考虑到规范要求和历史设计资料一致性问题，采用年最大值法(AM)提取每年的最大日降雨量和最大风速作为年最大日降雨量和年最大风速，并重新组成年最大日降雨量和年最大风速序列，二者其变化趋势见图1。此外，分别使用滑动法获得的年最大3 d和7 d降雨量并形成新的暴雨序列，并取对应时间段内日最大风速作为风速序列，以分析二者之间的遭遇规律。

a)年最大日降雨

b)年最大风速

2 研究方法

首先分别采用P-III、Gamma和GEV 3种边缘分布函数对3组降雨和1组风速序列进行拟合，并采用K-S检验和PDE方法选出最优边缘分布函数。然后采用Gumbel、Frank和Clayton Copula 3种Copula函数构建强风暴雨的联合概率模型，并采用OLS和AIC准则选出最优拟合函数。最后分别计算风雨两变量的“或”“且”和Kendall重现期设计值，并基于风速分析结果与现行技术规范进行对比分析。

2.1 边缘分布函数

边缘分布函数均已在风速、暴雨、风暴潮、洪水等序列的极值分布规律分析方面取得了理想的效果[20]。通过综合对比和考虑，最终选择皮尔逊III型(P-III)、伽马分布(Gamma)、广义极值分布(GEV)3种常用函数，各函数见表1。

表1 独立概率分布函数与概率密度函数

边缘分布采用概率分布误差率(Probability Distribution Error，PDE)[21]和Kolmogorov-Smirnov(K-S)检验[22]进行拟合优度评价，以确定最优边缘分布函数。其中概率分布误差率计算公式如下：

(1)

式中Fn(x)——经验分布函数；F(x)——理论分布函数；n——样本的个数。

2.2 Copula函数

水文分析中常用的Copula函数包括经验Copula函数、椭圆Copula函数及Archimedean Copula函数等，其中Archimedean Copula函数因其计算简便、形式多样，已在实际应用中得到了广泛应用，效果良好[23-24]。本研究将选用Gumbel、Frank和Clayton Copula 3种函数(表2)，联合分布采用AIC(Akaike information criterion)信息准则法[25]及OLS(Ordinary least squares)离差平方和最小准则[26]进行拟合优度评价，AIC值与OLS值越小则表明拟合效果越好。

表2 常见Archimedean Copula函数

OLS准则是通过计算理论值和实测值的均方根误差RMSE来定量地评估拟合误差，其计算公式为：

(2)

式中Femp(xi1,xi2…,xin)——经验频率值；C(ui1,ui2…,uin)——理论频率值。

AIC指标用于对Copula函数拟合情况进行评价，其计算公式为：

C(ui1,ui2…,uin)]2}+2k

(3)

式中n——函数维数；k——模型参数个数。

2.3 重现期计算

水文分析中常用的重现期包括“或”重现期和“且”重现期，“或”重现期表示2个变量中至少有1个变量超过特定频率时的重现期，“且”重现期表示2个变量同时超过特定频率时的重现期[27]。“或”重现期和“且”重现期分别用式(4)、(5)表示：

(4)

(5)

式中P——同现、联合概率；X1、Y1——变量序列；Fx1、Fx2——边缘分布函数；F——联合分布函数；Tor、Tand——“或”“且”重现期。

Kendall重现期考虑了曲线概率等值情况，相比于“或”“且”重现期更为严谨[28]。该重现期计算中引入了Kendall测度KC来计算，其中KC与Copula函数有关，计算公式如下：

(6)

(7)

3 结果分析

3.1 相关性与遭遇统计

分别计算1973—2018年最大风速与最大日降雨量、最大3 d降雨量、最大7 d降雨量序列的相关性，可得其相关系数分别为0.153、0.012和-0.036，表明3组序列相关性极弱。然而，相关性分析难以获悉二者的实际遭遇情况。分别选取年最大风速以及最大日降雨、最大3 d和最大7 d降雨统计大隆洞水库的强风暴雨遭遇情况(表3)。在46年的样本中，当大隆洞水库发生年最大风速时刚好遭遇最大日降雨的次数为6次，占比为13.04%；当发生年最大风速刚好遭遇最大3天降雨的次数为8次，占比为17.39%；当发生年最大风速刚好遭遇最大7 d降雨的次数为7次，占比为15.22%。以上表明，尽管大隆洞水库的风速和暴雨序列相关性极弱，但出现年最大风速和年最大暴雨遭遇的情况并非罕见。

表3 大隆洞水库年最大风速与不同时段最大降雨量遭遇情况

3.2 边缘分布及Copula函数拟合优度评价

以1973—2018年间年最大日降雨、最大3 d降雨、最大7 d降雨以及对应时段内的日最大风速为样本，采用P-III、GEV和Gamma 3种函数进行边缘分布拟合，边缘分布采用K-S、PDE检验进行拟合优度评价(表4)。

表4 各序列边缘分布的参数估计及拟合优度

结果表明，最大日降雨和对应日最大风速的最优函数均为GEV分布，最大3 d降雨和对应日最大风速的最优函数分别为P-III分布和GEV分布，最大7 d降雨和对应日最大风速的最优函数均为GEV分布。选择Gumbel、Frank及Clayton Copula 3种函数分别构建不同序列间的联合分布模型，采用AIC信息准则法和OLS准则进行拟合优度评价，Copula函数拟合优度与OLS、AIC值呈负相关关系，结果见表5。

表5 不同序列Copula联合分布函数拟合优度

结果表明，日最大风速与最大日降雨、最大3 d降雨、最大7 d降雨联合分布的最优函数均为Clayton Copula函数，图2为不同序列间联合概率分布图。随着降雨量上限从最大日降雨的300 mm增加至最大3 d降雨的500 mm时，风速上限有一定增大，但增幅较小；随着降雨量上限继续增加至最大7 d降雨的600 mm时，风速上限反而呈现一定减小，表明风速最大值与降雨量最大值之间并没有显著正相关关系。

a)风速与最大日降雨

3.3 风速和降雨量设计值

3种遭遇情景下，采用等概率法分别计算200、100、50、20、10、5年一遇条件下的最大风速和最大降雨量单变量设计值，“或”“且”与Kendall重现期的结果见表6。结果表明，最大风速与最大日降雨量遭遇情景下，与单变量的边缘分布结果相比，按两变量“或”重现期推算的最大风速与降雨量设计值的相对差异范围分别为25.1%～25.7%、5.8%～12.7%，高于边缘分布设计值。按两变量“且”重现期推算的最大风速与降雨量设计值的相对差异范围分别为26.1%～56.6%、17.6%～24.4%，低于边缘分布设计值。而按两变量Kendall重现期推算的最大风速与降雨量设计值的相对差异范围分别为18.9%～32.1%、10.7%～12.1%，低于边缘分布设计值。可见，相比于单变量风速与降雨量设计值，按风速和降雨量联合分布的“且”、Kendall重现期推算的风速和降雨量低于单变量的边缘分布结果，而采用“或”重现期推算的风速和降雨量则高于单变量的边缘分布结果。最大3、7 d降雨遭遇条件下具有相似的结论。

表6 不同取样条件下不同重现期设计值

3.4 与当前技术规范的差异对比

前已述及，当前大多数技术规范没有定量考虑强风和暴雨的遭遇问题，可能会对水库的安全性产生一定的威胁。例如，SL 274—2020《碾压式土石坝设计规范》中坝顶超高的计算所采用的设计风速是多年平均最大风速的1.5～2.0倍。根据风速资料计算可知，大隆洞水库在1973—2018年的多年平均最大风速为13.91 m/s，其2.0倍则为27.83 m/s。然而，根据表6可知，最大风速与最大日降雨量遭遇情景下，等概率发生条件下100年一遇的“或”联合分布最大风速值为36.1 m/s，显然比2.0倍多年平均最大风速还要大29.7%。根据SL 274—2020《碾压式土石坝设计规范》有关坝顶高程的计算，当取2.0倍多年平均最大风速时，大隆洞水库主坝高程(坝顶+防浪墙)为39.85 m，而该水库主坝实际高程则为40.20 m，表明按照当前的规范计算能满足安全要求；然而，当采用100年一遇的“或”联合分布最大风速36.1 m/s时，其主坝高程则为40.71 m，显然超出了实际高程(表7)。

表7 计算坝顶高程与实际高程对比

4 讨论

由上述分析可知，当前大隆洞水库大坝高程是满足现有规范要求的。根据调查得知，大隆洞水库在最近几十年汛期和台风期间未曾出现水位过高导致的险情，这是由于除了高程满足要求外，还得益于近年来天气预报准确性不断提高以及水库的科学预警和调度管理，使水库在台风来临前降至安全水位，大大降低了波浪爬高和水位雍高。然而，由于大隆洞水库位于广东沿海地区，由前述可知该地区的登陆台风呈现频数减少但强度明显增强的特点，而降雨则表现出强降水增加，弱降水减少的特点。图1也表明该水库年最大降雨极值呈显著上升趋势(p<0.05)，而尽管最大风速整体呈下降趋势，但在2000年以后也呈上升趋势。在未来气候变化背景下，强风和暴雨的强度可能会进一步增加，二者遭遇的概率有可能会提高。当前的技术规范仅考虑风速这一单变量，在没有考虑二者遭遇情况下，仅采取2.0倍多年平均最大风速所计算出来的坝顶高程有可能难以满足安全要求。因此，考虑到未来气候变化，建议有关技术规范可考虑增加强风和暴雨遭遇的情况，适当考虑“或”联合分布的分析结果；或者设计人员在应用当前相关规范的基础上适当增加坝顶高程的安全余度，尽量满足二者遭遇时的高程要求，确保水库大坝的安全。

5 结论

以大隆洞水库为研究实例，采用Copula函数构建了强风-暴雨联合概率分布模型，推算了不同遭遇条件不同重现期的设计风速与降雨量设计值，并基于遭遇分析的风速结果与当前技术规范进行对比，得出如下主要结论。

a)大隆洞水库年最大日降雨呈显著上升趋势(p<0.05)，年最大风速则呈总体下降趋势，但2000年以后呈上升趋势；尽管水库的风速和暴雨序列的相关性极弱，但曾出现年最大风速和年最大暴雨遭遇的情况。

b)相比于单变量风速与降雨量设计值，基于Copula函数推算的结果中，联合分布的“且”和Kendall重现期推算的风速和降雨量低于单变量边缘分布的结果，而采用“或”重现期推算的风速和降雨量则高于单变量边缘分布的结果。

c)按照当前规范要求，当取2.0倍多年平均最大风速时所计算出来的高程低于大隆洞水库主坝实际高程，但当采用100年一遇的“或”联合分布的最大风速所计算出来的主坝高程比实际高程要高；考虑到未来气候变化，有关技术规范可考虑增加强风和暴雨遭遇情况并适当考虑“或”联合分布的分析结果，或基于当前规范适当增加坝顶高程的安全余度，尽量满足二者遭遇时的高程要求。