东湖水库强风暴雨遭遇联合分布及其对坝顶高程设计的影响

蔡奇霖， 涂金良， 王兆礼， 赖成光

(1.华南理工大学 土木与交通学院， 广东 广州 510640； 2.广东省水利水电技术中心， 广东 广州 510635)

1 研究背景

水库在防洪、灌溉、供水、发电等方面对我国国民经济发展起着不可替代的支撑作用，水库大坝的安全问题也是涉及人民生命财产安全的一个重大问题[1-2]。据统计[3]，1954-2006年之间全国有3 498座水库溃坝。溃坝事故主要发生在1990年以前，其中1954-1990年共有3 260座水库溃坝，占总溃坝数的93.2%，年均约88座；1991-2000年，共有227座水库溃坝，年均约23 座；2001-2006年，共有35座水库溃坝，年均约6座。在此期间有2个溃坝高峰，一个是 1960年前后，即1959-1961年间，共计溃坝507座，每年约253座；另一个高峰期在1973年前后，仅1973年就溃坝554座。进入 20 世纪 90年代以来，特别是2001年以后，由于设计、建设以及管理水平的不断提高使得年溃坝数量明显减少。虽然如此，每年仍有溃坝事故发生。毫无疑问，合理的设计是保证大坝安全运行的前提条件，对保障大坝下游人民的生命财产安全意义重大。

对于水面面积大、周边空旷的水库，因其吹程较大容易形成大风浪，导致风壅水位过高以及波浪爬高过大而引起漫顶或跃顶现象，进而威胁大坝的安全。因此，在大坝设计过程中，确定合理的设计风速对雍水及风浪爬高计算至关重要，也决定了大坝的综合超高以及安全裕度。近几十年来，全球气候变暖导致众多地区风速和降雨特征均发生显著变化[4-6]。不少研究发现，我国沿海地区的强风和暴雨事件发生频率增加[7-9]，这也意味着两种事件遭遇的可能性也会进一步提高。与仅考虑单一变量(如降雨或风速)的分析相比，双变量(如降雨和风速)遭遇的联合分布情况更为复杂，这将会为水库设计风速的计算带来极大的挑战。

对于风速和暴雨联合分布方面的研究，国内外学者做了大量工作。武占科等[10]对上海地区的台风风速和降雨量建立Copula联合分布模型，提出了台风条件下风速和雨强联合概率分布函数，结果表明Gumbel Copula函数联合概率分布模型适合描述台风条件下风速和雨强联合概率分布。陈立华等[11]对钦州市台风影响下的风雨联合概率分布模型进行了研究，得到了适用于钦州市的风雨联合概率分布模型，并对两场台风降雨数据进行了检验，结果表明该风雨模型结果可靠，可应用于台风影响下钦州市降雨情势的预测分析。Dong等[12]应用Copula函数对崇明岛进行风速和雨强的联合回归概率分析，建立了G-H Copula联合概率模型，通过实测数据验证了该模型的效果十分理想。王修勇等[13]以台风最大风速和最大1 h时段雨量为样本，选出最优Copula风雨联合分布函数模型，并对琼州海峡大桥风雨极值重现期进行了分析，结果表明边缘分布重现期计算出的风速、雨量值均小于联合分布重现期的计算值。侯静惟等[14]以298场热带气旋最大3 s极值风速与总降水量为研究样本，建立了适用于海南岛的热带气旋风雨联合概率模型并分析了风雨单个致灾因子超过阈值和两个致灾因子超过阈值两种联合重现期，结果表明联合重现期RPand相比于联合重现期RPor与历史热带气旋归一化损失率具有更好的相关性。Um等[15]采用Copula函数分析了风速与台风降雨的相关关系，并将其应用于韩国济州岛的台风降雨的预测。Bushra等[16]应用Gumbel Copula函数分析了孟加拉湾沿岸风速与风暴潮的联合风险，结果表明该模型结果可靠。上述研究均表明，考虑双变量遭遇情况的计算结果比仅考虑单一变量更加接近实际情况；然而，上述研究大部分聚焦于相对较大空间尺度的风速和暴雨遭遇问题，而对于沿海地区水库在风雨遭遇条件下有关设计风速的研究尚不多见。

广东省沿海地区属于东亚季风区，降水充沛，台风登陆频繁。广东省沿海地区水库众多，但这些水库无论是在20世纪50、60年代的建设期间还是在2000年前后的除险加固过程中，其设计风速均没有考虑强风暴雨的遭遇影响。在极端暴雨频次与台风强度不断增加的背景下，这些水库大坝的高程是否仍能满足安全要求，需要进一步复核。目前的技术规范如《碾压式土石坝设计规范》(SL 274—2020)和《砌石坝设计规范》(SL 25—2006)是按大坝的重要性给出相应的风速标准，也尚未考虑风速与暴雨的遭遇情况；尽管《小型水利水电工程碾压式土石坝设计规范》(SL 189—2013)要求在沿海地区应该考虑最大风力与暴雨同时出现的条件，但如何执行也没有明确规定。

鉴于此，本文将以位于广东省沿海地区的东湖水库为研究对象，探讨变化环境下强风暴雨遭遇的联合分布特征，应用Copula函数构建强风-暴雨联合概率分布模型，计算设计风速与暴雨的重现水平，最后基于联合分布的设计风速对水库主坝坝顶和防浪墙高程进行安全复核。研究成果可为我国沿海地区水库大坝的设计风速计算和坝顶高程设计提供参考依据。

2 数据来源与研究方法

2.1 数据来源

东湖水库位于广东省阳江市阳东区那龙河上游，于1960年建成，2006年经加固改造后集雨面积达51.15 km2，水库总库容为1.27×108m3，为大(Ⅱ)型水库，设计重现期为100年一遇。水库主坝坝型为均质土坝，坝顶高程为36.7 m，坝顶长度为130 m；水库正常蓄水位为31.00 m，设计洪水位为32.95 m，校核洪水位为33.73 m。

所采用的降雨和风速数据来源于水库附近的阳江气象站(东经111°58′、北纬21°51′)，其中降雨(逐日)数据序列长度为1960-2018年，最大风速(逐日)数据序列长度为1971-2018年。提取每年最大日降雨和最大风速作为年最大日降雨和年最大风速，并重新组成年最大日降雨和年最大风速序列；使用滑动法分别获得年最大3 d和7 d降雨序列，并分别提取二者对应时间段内的日最大风速作为年最大风速序列。上述序列均经过平稳性分析。

2.2 边缘分布函数

目前常用的边缘分布函数包括皮尔逊Ⅲ型(P-Ⅲ)、伽马分布(Gamma)、广义极值分布(generalized extreme value distribution, GEV)、极值Ⅲ型(Weibull)、极值Ⅰ型(Gumbel)等，均取得理想的效果[17]。通过综合对比和考虑，本研究最终选择P-Ⅲ、Gamma、GEV 3种常用函数，各函数表达式如表1 所示。

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

(1)

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

2.3 Copula函数

Copula 函数可描述变量间的相关性，常用的Copula函数包括经验Copula函数、椭圆Copula函数以及Archimedean Copula函数等。Archimedean Copula函数因其计算简便、形式多样，在实际工程中得到了广泛应用[20-23]。本研究将选用Gumbel Copula、Frank Copula 和Clayton Copula 3种常见二元函数，各Copula函数表达式见表2。应用AIC(Akaike information criterion)信息准则法[24-25]及OLS(ordinary least squares)离差平方和最小准则[26]进行拟合优度评价，AIC值与OLS值越小，则拟合效果越好。

注：表达式中的α、β、γ分别为形状、尺度、位置参数。

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

OLS=

(2)

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

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

(3)

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

Copula函数拟合优度与OLS、AIC值呈负相关关系。

2.4 联合与同现重现期

对水文研究领域的遭遇问题，联合风险率和同现风险率是关注的重点。联合风险率表示至少有1个出现非期望事件的概率，同现风险率表示两个非期望事件同时出现的概率。同现联合概率可用公式(4)表示，联合重现期和同现重现期分别用公式(5)、(6)表示。

P[X1>x1,X2>x2]=1-Fx1(x1)-Fx2(x2)+

F(x1,x2)

(4)

(5)

(6)

式中：P为同现、联合概率；X1和X2为变量序列；Fx1和Fx2为边缘分布函数；F为联合分布函数；Tor和Tand分别为同现、联合重现期。

2.5 坝顶超高计算

选用《碾压式土石坝设计规范》(SL 274—2020)[27]作为坝顶超高计算依据，其规定坝顶高程等于水库静水位加坝顶超高，应分别按4种运用条件计算并取其最大值，坝顶超高由下式确定：

y=R+e+A

(7)

式中：y为坝顶超高，m；R为波浪在坝坡上的爬高，m；e为风雍水面高度，m；A为安全超高，m。

设计波浪爬高值应根据工程等级确定，1、2级坝采用累计频率为1%的波浪爬高值R1%。由该规范附录表A.1.13查得R1%/Rm=2.23。此外，规范的条件1中建议“采用多年平均最大风速的1.5～2.0倍”，本次分别采用上限2.0倍和下限1.5倍计算，以探讨二者计算结果的差异。

以阳江站1971-2018年间的年最大日、最大3 d、最大7 d 3组降雨序列及其对应年最大风速序列作为计算样本，应用3种分布函数和K-S检验进行边缘分布拟合，随后采用Copula函数、AIC、OLS准则进行联合概率分布拟合，最后计算联合与同现重现期下的风速和降雨量。分别采用遭遇风速与规范推荐的设计风速计算坝顶高程，并与实际坝顶高程进行对比，综合复核东湖水库主坝坝顶高程。

3 结果与分析

3.1 年最大日降雨量与最大风速变化趋势

分别对3组降雨及风速序列进行趋势分析，以年最大日降雨量与年最大风速序列为典型进行分析，二者变化趋势如图1、2所示。由图1、2可知，1960-2018年间，年最大日降雨量呈显著增大趋势(P<0.05)；年最大风速总体呈现不显著增大趋势(P<0.05)，其中1971-2000年呈减小趋势而 2001-2018年呈明显增大趋势。

图1 1960-2018年研究区域年最大日降雨量变化趋势

3.2 相关性与遭遇统计

分别计算年最大日、最大3 d与最大7 d降雨量及其对应年最大风速序列的相关性，可得其相关系数分别为0.007、0.128和0.156，表明3组序列相关性极弱。分别选取年最大风速以及最大日降雨量、最大3 d和最大7 d降雨量统计东湖水库的风雨遭遇情况(表3)，当东湖水库发生年最大风速时刚好遭遇最大日降雨的次数为2次，占比为4.17%；当发生年最大风速时恰好遭遇最大3 d降雨的次数为5次，占比为10.42%；当发生年最大风速时刚好遭遇最大7 d降雨的次数为10次，占比为20.83%。以上表明尽管最大风速与最大降雨序列相关性极弱，但出现了年最大风速和年最大暴雨遭遇的情况。

表3 东湖水库年最大风速与不同时段最大降雨量遭遇情况统计

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

以1971-2018年间年最大日、最大3 d、最大7 d降雨量及其对应年最大风速为样本，采用P-Ⅲ、GEV和Gamma 3种平稳分布函数进行边缘分布拟合，应用K-S、PDE检验进行拟合优度评价，各序列的边缘分布拟合优度如表4所示。

表4中的拟合优度结果表明，最大日降雨和对应日最大风速的最优函数均为P-Ⅲ分布，最大3 d降雨和对应日最大风速的最优函数分别为GEV分布和P-Ⅲ分布，最大7 d降雨的最优函数为P-Ⅲ分布但对应日最大风速为Gamma分布。

表4 各序列边缘分布拟合优度

采用Gumbel Copula、Frank Copula及Clayton Copula 3种函数分别构建不同序列间的联合分布模型，采用AIC信息准则法和OLS准则进行拟合优度评价，结果如表5所示。

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

注：加粗字体表示该边缘分布函数拟合优度最佳。

注：加粗字体表示该Copula函数拟合优度最佳。

由表5可看出，日最大风速与最大日降雨量联合分布的最优函数均为Frank Copula函数，日最大风速与最大3 d降雨量联合分布的最优函数均为Clayton Copula函数，日最大风速与最大7 d降雨量联合分布的最优函数均为Frank Copula函数，图3为日最大风速与最大日降雨量、最大3 d降雨量、最大7 d降雨量的联合概率分布图。

图3 不同序列间联合概率分布图表6 不同重现期和计算工况下坝顶及防浪墙高程计算结果

3.4 风速与暴雨重现水平

经过多种组合计算，发现强风和暴雨在低重现期遭遇组合条件下对水库大坝安全性影响有限，因此仅选取二者均为100年一遇(等概率法)遭遇的极端情况的计算结果进行分析。同现重现期、联合重现期以及单变量重现期不同时段下，最大风速和最大降雨量相应坝顶及防浪墙高程计算结果如表6所示。

由表6可知，基于不同联合序列所得出的风速以及降雨量均存在较大差异，其中联合重现期计算值最大，同现重现期计算值最小，而单变量重现期位于二者之间。对于最大风速，联合重现期比同现重现期大1.5～1.9倍，但联合重现期最大风速仅比考虑单变量的最大风速大10%左右；而对于最大降雨量，联合重现期则比同现重现期大1.7～1.9倍。

取样条件同现重现期最大风速/(m·s-1)最大降雨量/mm坝顶高程/m联合重现期最大风速/(m·s-1)最大降雨量/mm坝顶高程/m单变量重现期最大风速/(m·s-1)最大降雨量/mm坝顶高程/m2.0倍多年平均最大风速最大风速/(m·s-1)坝顶高程/m1.5倍多年平均最大风速最大风速/(m·s-1)坝顶高程/m实际坝顶高程/m实际防浪墙顶高程/m最大日降雨量14.97393.3435.3823.68681.0836.3021.80616.8636.10最大3 d降雨量16.69557.1135.5627.151106.5036.6724.83960.4736.4235.6737.5926.7536.6336.7037.50最大7 d降雨量16.94653.3635.5833.011221.6237.3129.721106.7936.95

3.5 坝顶高程复核

根据阳江气象站实测1971-2018年最大风速资料可求得多年平均最大风速为17.83 m/s。对于遭遇情景，选取同现重现期、联合重现期以及单变量重现期的风速最大值并根据规范推求坝顶高程。在3种取样条件下不同重现期坝顶高程最大值、规范计算值以及实际坝顶高程如表6所示。通过调查可知，东湖水库实际坝顶高程为36.70 m，防浪墙顶高程为37.50 m，考虑到东湖水库大坝防浪墙防御能力较强，因此以防浪墙的高程作为依据评判是否安全。通过计算发现，设计洪水位加正常运行条件下的坝顶超高值最大，因此选取该条件作为坝顶高程计算条件，且分别取规范的上限2.0倍和下限1.5倍多年平均最大风速进行分析。

结果表明，当设计风速取1.5倍多年平均最大风速时，所计算的坝顶高程为36.63 m，当前坝顶和防浪墙实际高程均能满足要求；而当设计风速取2.0倍多年平均最大风速时，坝顶高程为37.59 m，比防浪墙实际高程略高。对于遭遇中最危险工况，即联合重现期(最大7 d降雨取样)为100年一遇时所计算的坝顶高程为37.31 m，当前防浪墙高程仍能满足要求。

4 讨 论

不少研究表明[7-9]，最近几十年来台风的强度在不断提高，风速的极值也在不断增大，而当风速序列增加这些极值样本时会大大提高多年平均最大风速值。由于东湖水库在2006年进行了除险加固，当年所使用的风速序列相对较短；由最大风速变化趋势(图2)可知在2000年以后年最大风速呈现显著上升趋势，但后续这些极值样本大部分没有被纳入除险加固的设计风速计算，这可能是导致当前防浪墙高程略小于2.0倍多年平均最大风速时所推算高程的主要原因。而对于遭遇情况，当考虑强风和暴雨联合分布时，同现重现期和联合重现期的最大风速与单变量最大风速存在较大差异，其中联合重现期最大风速比仅考虑单变量的最大风速大10%左右。通过对比分析可知，东湖水库的2.0倍多年平均最大风速仍大于联合重现期最大风速，表明《碾压式土石坝设计规范》(SL 274—2020)在考虑强风暴雨遭遇时仍可满足安全要求。此外，通过调查可知，尽管主坝的防浪墙高程未能满足2.0倍多年平均最大风速的高程要求，但最近几十年在台风登陆前后基本上没有出现水位过高而导致险情的情况，这主要是因为水库从未遭遇到2.0倍多年平均最大风速的情景；此外，日益精准的天气预报以及水库的科学调度也有利于水库提前做好防御措施，例如在获知台风预报后可通过迅速泄洪以降低水位，使波高、风浪雍高大大降低，从而保证水库安全度过险情。当然，本研究仅能说明东湖水库在考虑强风暴雨遭遇时能满足安全要求，其它水库是否满足仍需要具体分析。尽管通过科学的管理和预报手段可以降低风险，但在全球变暖驱动下，极端暴雨与台风强度有可能进一步加强，强风和暴雨联合出现的概率也可能会进一步增大，因此在水库大坝特别是沿海地区的水库大坝设计和建设过程中，建议增加考虑强风和暴雨遭遇的情况，根据计算结果可适当增加坝顶高程的安全裕度，以应对未来发生极端暴雨和台风的情况。

图2 1971-2018年研究区域年最大风速变化趋势

5 结 论

以东湖水库为研究对象，利用广义极值分布(GEV)、伽马分布(Gamma)和皮尔逊Ⅲ型分布(P-Ⅲ)分别对最大风速和降雨序列进行边缘分布拟合，基于Copula函数构造最大风速和降雨序列的联合概率分布，估算100年一遇重现期下最大风速和降雨量设计值，并复核了水库主坝坝顶高程，得出如下结论：

东湖水库1960-2018年间年最大日降雨量呈显著上升趋势(P<0.05)，最大风速则总体呈现不显著上升趋势(P<0.05)，其中1971-2000年呈下降趋势，2001-2018年呈明显上升趋势，且出现年最大风速与年最大暴雨遭遇的情况。遭遇分析结果表明，联合重现期计算值最大，同现重现期计算值最小，而单变量重现期位于二者之间。对于最大风速，联合重现期比同现重现期大1.5～1.9倍，但联合重现期最大风速仅比考虑单变量的最大风速大10%左右。当设计风速取1.5倍多年平均最大风速时，当前坝顶和防浪墙高程均能满足要求；当设计风速取2.0倍多年平均最大风速时，所推算的高程比当前实际防浪墙高程略高。对于遭遇最危险工况(联合重现期)，当前的防浪墙高程能满足要求。东湖水库的2.0倍多年平均最大风速所计算的防浪墙高程仍大于联合重现期最大风速所推算的高程，表明《碾压式土石坝设计规范》(SL 274—2020)在考虑强风暴雨遭遇时仍可满足要求。但考虑到未来气候变化，建议增加考虑强风和暴雨遭遇的情况，可适当增加坝顶高程的安全裕度以应对未来发生极端暴雨和台风的情况。