尤再进

海洋重现期波高统一化计算方法*

尤再进

(大连海事大学港口与航运安全协创中心 辽宁大连 116085)

重现期波高是港口海岸及海洋工程设计中不可回避的一个重要设计参数, 尤其对深水海港、海上平台、海底油气管道、沿海核电站等重大涉海工程设计具有巨大的经济价值和深远的社会效益。但是, 现有重现期波高推算缺乏统一的计算方法, 导致计算结果相差悬殊。研究重现期波高的统一化计算方法, 分析重现期波高计算中存在的各种不确定因素, 提出减少这些不确定因素的新方法, 建立误差小、应用方便、方法统一的重现期波高计算方法。基于澳大利亚悉尼站的长期连续观测波浪数据, 研究发现: 广义帕累托函数(generalized Pareto distribution III, GPD-III)和威布尔(Weibull)是重现期波高计算的最佳候选极值分布函数, 新推导的函数形状参数计算公式较好提高重现期波高的计算精度, 极值波高数据的分析方法和样本大小是影响重现期波高计算精确度的两个重要因素, 短期波浪资料和年极值法可能高估重现期波高值。逐个风暴的极值波高数据分析法及最佳候选极值分布函数GPD-III和Weibull建议应用于涉海工程设计的重现期波高推算。

重现期; 重现期波高; 极值分析; 广义极值函数; 广义帕累托函数

海洋波浪重现期/设计波高是港口海岸和海洋工程设计中不可回避的一个重要设计参数, 尤其对深水海港、沿海核电站、海洋平台、海底油气管线等重大涉海工程设计具有巨大的经济价值和深远的社会效益(尤再进, 2016)。国际学者采纳的设计波高或重现期波高的计算方法可以归纳以下三种: (1) 年个大波法(annual-largest method, ANL; 曹兵等, 2006), 此方法是由年极值法(=1)延伸而来(Sobey, 1995), 从时间序列的波浪数据中每年提出个独立的大浪峰值(≥1), 组成一个极值波高数据样本, 然后应用广义极值函数(generalized extreme value, GEV)推算多年一遇的重现期波高; (2) 阈值法(peaks-over-threshold method, POT; Hosking, 1987), 该方法是设置统一的波高阈值0, 从时间序列波浪数据中提出大于0的独立大浪峰值, 组成一个极值波高数据样本, 然后应用广义帕累托函数(generalized Pareto distribution, GPD)推算极值波高; (3) 个风暴法(storm-by-storm method, SAS; Goda, 1988), 是从时间序列的波浪数据中提取所有风暴波浪峰值, 组成一个风暴波浪峰值数据样本, 然后应用右偏分布的概率函数(如Weibull、Pearson-III等)推算重现期波高。

但是, 这三种重现期波浪计算方法通常会导致不同的计算结果。首先, 这三种极值波浪数据分析方法(ANL, POT, SAS)将会导致显著差异样本的重现期波高数据。Sartini等(2015)采用ANL和POT两种方法来计算百年一遇的极值波高, 发现这两种广义分布函数的计算结果差别很大。POT方法的难点是如何确定一个合适的波高阈值0(Liang, 2019)。You (2012)探讨了如何正确地选择极值波高阈值, 认为SAS是POT的一种特殊情况, 当POT的阈值等于SAS的当地最小台风/风暴波浪峰值时, 这两种方法就变成一种方法。ANL方法每年都提取个独立和随机的大浪峰值, 而不能保证每年的波高阈值相同或者认为波高阈值是时间的变量, 而POT方法保持每年的波高阈值相同, 而每年获取的风暴波浪数不等, 或者认为POT方法每年提取的独立和随机大浪数是一个变量。所以, 这两种不同方法生成不同的极值波高数据样本, 导致重现期波高数据样本的不确定性。但是, 如果ANL方法每年选取的个独立的最大波浪数据点是足够多(如≥每年独立台风浪个数), 这两种方法能够生成几乎相同的重现期波高数据样本(You, 2015)。

其次, GEV和GPD是由不同的极值分析法而推导出的两种广义分布函数。理论上GEV只能由ANL方法生成的极值波高数据样本计算重现期波高, 同样GPD只能由POT方法生成的极值波浪样本来推算重现期波高。ANL或POT方法生成的极值波高数据样本和右偏分布的普通概率函数(如Weibull, Pearson-III)也应用于重现期波高计算(尹宝树等, 2002)。Isaacson等(1981)指出, 没有一个极值分布函数是最佳的, 只有与具体数据比较以后, 才能够判断出哪个函数是最佳函数。

再者, 国内外学者们应用多种不同函数参数估算法计算极值函数参数, 如矩法(method of moment, MOM), 最小二乘(least-square method, LS), 极大似然(maximum likelihood method, ML), 概率权重矩(probability weighted moments, PWM)等参数估算方法, 其中MM, ML和LS是重现期波高计算中常用的三种主要参数估算法, 这些方法的物理含义比较清晰(陈子燊, 2011; You, 2015)。MM法应用比较简单, 从极值波高数据样本中分别估算出样本的均值、方差、偏度等特征参数。这就是概率论中-阶矩的定义, 也就是波高变量的次方对极值概率密度的积分。所以, MM方法通常适用于一些简单的两参数概率函数(如极值I型分布, FT-I), 而且极值波高数据样本量需要很大, 使得计算的极值波浪样本特征值(均值、方差、偏度)趋于稳定。ML方法适用于大多数极值函数, 但是在估算三参数或多参数函数的参数时(如Weibull, GEV, GPD), 其收敛速度比较慢, 甚至不收敛(Mazas, 2011)。You (2011)推广了LS方法, 高精度计算Weibull、GEV、GPD的三函数参数。基于澳大利亚新南威士州海岸采集的长期和高质量波浪数据, You (2012)定量讨论了ML和LS两种方法的关系, 研究发现LS方法是ML方法的一种特殊形式。很多学者认为, LS方法需要计算极值波高的经验频率分布, 而ML方法却不需要, 能够直接基于极值波浪数据样本确定极值函数参数。但是, 事实上所有极值波浪计算方法均需要重现期波高数据, 因为它们的计算结果需要与重现期波高数据进行比较。现有的经验频率计算公式大多是从序列函数推导出来的, 表达形式随着函数不同而改变, 导致生成不同的重现期波高数据样本(Goda, 1988)。基于波浪重现期的定义, You等(2015)推导出唯一的极值波浪经验频率的表达式, 消除重现期波高数据样本的不确定性。

最后, 国内外学者建立了不同的评估标准, 选择最佳候选极值函数推算重现期波高。现今应用最多的评估标准是均方误差(mean squared error, MSE): MSE=Σ(观测的重现期波高–计算的重现期波高)2/, 其中,是重现期波高数据样本量。但是, 也有少部分学者采用观测和计算的极值波高经验频率分布的均方误差MSE作为选择最佳极值函数的评估标准(Liu, 2021), 而极值波浪经验频率分布不是直接观测的极值波浪数据, 而且计算方法也具有不确定性。You等(2015)采用延伸的LS函数参数估算法, 建立了最佳候选极值函数评估标准, 其方法收敛快且精度高。

该论文重点研究重现期波高计算的统一化方法, 分析重现期波高计算中存在影响计算结果的不确定因素, 提出减少这些不确定因素的新方法, 建立一种误差小、应用方便、方法统一的重现期波高计算方法。

1 重现期波高计算控制方程

其中,=/是重现期值定义(CEM, 2002),是年观测波浪资料中共有的极值波高数据大于或等于重现期波高T,是年观测波浪资料的极值波浪总数,=/是年极值波浪的年均个数。将方程(2)代入方程(1), 重现期波高计算控制方程推导为

其中,H是年一遇的重现期波高,(,)是函数变量。方程(3)是重现期波高计算的控制方程, 仅是重现期的函数。

控制方程(3)中的函数参数(,,)能够通过线性回归重现期波高数据(,H)来确定。由于函数变量是形状参数的隐性函数, 函数参数(,)随着的变化而变化。只有当值满足和方差SSE最小要求时,

并将其代入方程(5), 其第一项和第三项相加等于零, 方程(5)最后简化为

方程(8)是函数形状参数计算的控制方程, 乘号的前项是方程(6)的参数或斜率, 乘号的后项相当于函数参数的倒数(1/)。Newton-Raphson迭代法可以应用于方程(8)的求解, 迭代求解函数形状参数。虽然方程(8)与曹兵等(2006, 2007)和You等(2006)推导结果相同, 但是本文在推导过程中有显著差异, 尤其在方程(8)的推导过程中没有给出已知极值函数的具体表达式, 该方法适用于具有极值函数普通表达式=[(–)/,]的显函数形状参数估算。

2 极值波高数据分析的逐个风暴法

重现期波高数据(H)是估算方程(3)函数参数(,,)的基础数据, 是由极值波高数据(,)并通过方程(2)转换而成, 这两组数据(,)和(H)其实是等同的。在现有的三种常用极值波高数据分析方法中, 逐个风暴法(图1)是从长期时间序列的波高数据中逐个提取风暴/台风过程中的风暴波高峰值, 一个风暴过程仅提供一个风暴波高峰值的数据点, 提取的所有风暴波高峰值数据组成了独立和随机的极值波高数据样本。图1中的时间序列波高数据来源于澳大利亚新南威士州的悉尼波浪骑士站(You, 2008), 每个风暴过程中的风暴有效波高峰值均大于风暴波高阈值3 m。逐个风暴法与波高阈值法存在一些明显区别, 尤其逐个风暴法中的风暴阈值是唯一的, 可由大气压场确定风暴刚刚开始的初时海况(You, 2012), 而POT阈值法中的波高阈值至今缺乏取值的统一标准, 人为影响因素比较大。但是, 这两种方法也是相关的, 逐个风暴法可以看作是POT阈值法的一种特殊形式: 当逐个风暴法中的风暴波高阈值等于阈值法中的波高阈值时, 这两种方法是等同的。

图1 逐个风暴法应用于提取风暴过程中的波高峰值生成极值波高数据样本

本研究采用的重现期波高数据来源于澳大利亚新南威士州悉尼观测站的深水波浪数据(You, 2008)。该站的平均水深约80 m, 波浪数据从1976年开始采集, 一直持续至今。1986年前, 波浪数据纪录在图表纸上, 波峰值从图纸上测量。从1986起, 采用波浪骑士每0.5 s自动测量一个波高数据, 连续观测34 min, 每小时重复一次, 数据采集的成功率大于95%。基于悉尼站采集的25 a波浪数据(1987～2011年), 图2给出了逐个风暴法分析的风暴波高峰值数据, 其中小圆圈数据点代表每个风暴过程中的最大有效波高, 方框中的数字代表年最大有效波高(年极值法), 大圆圈中的数字代表年风暴发生数次。该论文采用的悉尼站极值波高数据样本是由图2中的从大到小排列的所有风暴波高峰值数据组成, 并应用方程(2)将极值波高数据(,)转化成重现期波高数据(,H)。

图2 逐个风暴法分析获取的澳大利亚悉尼波浪观测站的极值波高数据

注: 蓝色小圆圈数据点代表每个风暴波高峰值, 黑色方框中的数字代表年最大波高, 紫色大圆圈中的数字代表年风暴发生数次

3 候选极值分布函数的选取标准

重现期波高计算的候选极值函数众多, 大体可以归纳为3大类: 广义极值分布函数GEV、广义帕累托分布函数GPD、普通概率函数(如Weibull、Pearson-III等), 它们的共性是右偏分布的长尾巴函数, 不同点是GEV应用于年大波法的极值波高数据, GPD针对大于某个波高阈值的极值波高数据, 普通概率函数仅要求随机和独立的极值波高数据。因为GEV和GPD可以简化为等同的普通概率函数 (You, 2015), 所以三类候选极值分布函数可以统一归类为右偏分布的普通概率函数。

3.1 广义极值函数

如果暂先不考虑广义分布函数GEV推导的假设条件, GEV的复杂数学表达式可以直接简化为等同的指数函数表达式:

在方程(3)中, 给定的不同函数参数, 其他两个函数参数(,)由最小二乘法唯一确定。基于悉尼站的重现期波高数据, 图3分别给出了方程(4)计算的均方误差MSE()和重现期波高100随简化的FT-II和FT-III形状参数变化的分布规律。由图3可见, FT-II无最大值, 最小值收敛于FT-I; FT-II计算的MSE()存在最小值, 且当>100时逐渐收敛于FT-I计算的MSE(); FT-III无最小值, 最大值收敛于FT-I; 计算的MSE(<100)不存在最小值, 只有当>100时逐渐收敛于FT-I计算的MSE()。所以, GEV的三种极值分布函数只有FT-II适用于重现期波高计算的候选函数, 三参数的FT-II包括二参数的FT-I函数。

图3 简化的广义极值分布函数(generalized extreme value, GEV)计算的均方误差MSE和重现期波高H100随简化的FT-II和FT-III形状参数a变化的分布规律

简化广义极值函数GEV有两个主要目的, 一是将复杂GEV表达式等同地简化为指数函数表达式, 并归类于普通的概率函数类; 二是能够极大简化方程(10)～(11)的函数变量求导, 应用方程(8)精确计算函数形状参数。

3.2 广义帕累托函数

参考广义极值分布GEV的简化目的和方法, 广义帕累托分布GPD的复杂数学表达式也可以改写成简单的幂函数表达式:

由方程(11)计算, 其中是简化后的GPD函数变量。当a确定后, 函数尺度参数A由方程(6)计算, 函数位置参数为。

3.3 概率分布函数

几种常采用的右偏分布概率函数(如Weibull和Pearson-III)也应用于重现期波高计算。虽然简化后的GEV方程(9)的FT-III(>0)与Weibull函数相似, 但不适用于重现期波高的计算(图3), 且与Weibull(>0)分布函数完全不同

基于相同的澳洲悉尼站重现期波浪数据, 图5分别给出了均方误差MSE随Weibull和Pearson-III形状参数变化的分布关系, 以及百年一遇的重现期波高100与形状参数的变化规律。由图5可以看出, Weibull和Pearson-III函数计算的MSE()均存在最小值, 但Weibull计算的MSE()要比Pearson-III收敛要快。

4 候选极值函数的最佳拟合标准

4.1 最佳拟合度评估法

候选极值函数与重现期波高数据(,H)的拟合程度可由方程(4)计算的均方误差MSE或者和方差SSE值来评价。当MSE()最小或方程(8)()=0, 候选极值函数与重现期波高数据的拟合程度达到最佳。基于悉尼波浪站的重现期波高数据, 两种简化的广义极值函数FT-II方程(9)和GPD-III方程(12)与重现期波高数据(,H)的线性拟合度见图6, 其中是极值函数变量, 函数形状参数由方程(8)计算, 方程(1)中的函数尺度和位置参数(,)由最小二乘法确定。

图5 威尔布(Weibull)和皮尔逊3型(Pearson-III)函数计算的均方误差MSE和重现期波高H100随函数形状参数a变化的分布规律

图6中的两种重现期波高数据(,H)和(, H)是等同的, 并具有相同的分布形态。这是因为当形状参数确定后,只是重现期的唯一函数。当函数参数确定后, 重现期是方程(1)中的重现期波高计算的唯一函数变量, 给定重现期的波高H能够从方程(1)获解。图6也给出了计算的和观测的重现期波高的定量比较, FT-II计算的重现期波高曲线具有上凸特点, 计算的重现期波高随增加而较快增加, 而且计算红线均高于重现期大于4 a的波高, 高估了重现期波高。相比较而言, GPD-III的重现期波高曲线具有上凹特点, 计算随增加而较缓慢增加, 函数与数据的拟合度2要比FT-II的高。

图7给出了Weibull和Pearson-III概率函数与重现期波高数据(,H)的拟合度, 其中形状参数先由方程(8)计算, 函数尺度和位置参数(,)再由最小二乘法计算。这两种概率函数与重现期波高的拟合度比较接近, 线性拟合系数2几乎相等。图7也给出了计算的重现期波高与重现期波高数据的定量比较, 计算的重现期波高曲线均具有上凹特点, Pearson-III计算的重现期波高稍微比Weibull计算的大些。这可能是因为Weibull函数形状参数是由方程(8)精确计算的, 而Pearson-III的表达式是一个半隐性式, 函数形状参数只能通过试算方法来确定, 计算精度比较低。再者, 由图5可见, 当MSE()靠近其最小值时, MSE()随形状参数变化缓慢, 收敛性较差, 而Weibull计算的MSE()收敛性要比Pearson-III的好些。

图6 极值函数形状参数估算法方程(8)应用于简化的GEV-II和GPD-III参数估算和澳洲悉尼波浪站的重现期波高推算

注: 图中表达式为拟合直线,为相关系数, 图7同

图7 极值函数形状参数估算法方程(8)应用于概率函数Weibull和Pearson-III参数估算和澳洲悉尼波浪站的重现期波高推算

4.2 数据量对最佳拟合度影响

基于重现期波高计算控制方程(1)的推导假设条件, 仅当MSE()值最小时, 候选极值函数与重现期数据的拟合度才能达到最佳。但是, 应用于计算MSE()的重现期波高数据量可能影响MSE()的收敛性以及该计算方法的适用性。同样基于悉尼站的重现期波高样本数据量, 不同样本的重现期波高数据(1>2>3, ….,>H)从悉尼站的重现期波高数据中(=485)提取前个重现期波高数据组成, 共计生产20个不同样本的重现期波高数据(=25, 50, 75, 100, …., 450, 475)。其实, 这里的极值波高数据样本量是由波高阈值的取值大小决定的。

图8给出了FT-II计算的MSE()收敛性以及重现期波高100随重现期波高数据量增加的变化规律。由该图可以发现, MSE()均存在最小值, 但是当<150时, 其最小值位置和收敛值随增大而分别向右方移动并明显减小, MSE最小值对应的重现期波高100随着增加而显著减小; 当≥150时, MSE()最小值位置和收敛值随增大而趋近稳定, 对应MSE最小值的重现期波高100随着增大而变化不显著。综上所述, 随着重现期波高数据增加, FT-II计算的MSE()始终存在最小值, 但仅当≥150时, 计算的重现期波高100才趋近稳定。这就要求采集的波高数据的时间长度(年)应该足够长或者数据点足够多, 满足生成的极值波高数据量大于150。

图8 GEV-II函数计算的MSE(a)收敛性和重现期波高H100随重现期波高数据样本量m的变化规律

图9 GPD-III计算的MSE(a)收敛性与重现期波高H100数据随样本大小m的变化规律

图10 Weibull函数计算的MSE(a)收敛性与重现期波高H100数据随样本大小m的变化规律

图11 Pearson-III和Gumbel (FT-I)计算的重现期波高与年极大值法分析的重现期波高数据的比较

注:表示年平均极值波高的数量

5 结论

重现期波高是港口海岸及海洋工程设计中不可回避的一个重要设计参数, 但其计算方法至今缺乏统一性, 计算结果存在显著不确定性。该论文建立了一种误差小、应用方便、方法统一的重现期波高计算方法。本文研究发现: GPD-III和Weibull是重现期波高计算的最佳候选推算函数, 新推导的它们形状参数计算公式较好提高重现期波高的计算精度, 极值波高数据的样本量和分析方法是影响重现期波高计算精度的2个重要因素。基于上述的研究发现, 我国《港口与航道水文规范》(JTS 145—2015)建议采用的年极值法和Pearson-III极值分布函数有可能高估推算的重现期波高, 但需要中国沿海长期且高质量现场波浪数据进一步验证。逐个风暴的极值波高数据分析法与广义帕累托函数GPG-II和Weibull建议应用于涉海工程的重现期波高计算。

尤再进, 2016. 中国海岸带淹没和侵蚀重大灾害及减灾策略[J]. 中国科学院院刊, 31(10): 1190-1196.

尹宝树, 何宜军, 侯一筠, 等, 2002. 推算波浪多年一遇波高的新方法[J]. 海洋与湖沼, 33(1): 30-35.

陈子燊, 2011. 波高与风速联合概率分布研究[J]. 海洋通报, 30(2): 158-163.

曹兵, 王义刚, YOU Z J, 2006. 三种设计波高计算方法比较[J]. 海洋工程, 24(4): 75-80.

曹兵, 王义刚, YOU Z J, 2007. 设计波高分布函数比较[J]. 海洋湖沼通报(4): 1-9.

CEM, 2002. Hydrodynamic Analysis and Design Conditions [M]. U.S Army Corps of Engineers, Coastal and Hydraulics Laboratory.

GODA Y, 1988. On the methodology of selecting design wave height [C] // Proceedings of the 21st International Conference on Coastal Engineering. Costa del Sol-Malaga: ASCE: 899-913.

HOSKING J R M, WALLIS J R, 1987. Parameter and quantile estimation for the generalized Pareto distribution [J]. Technometrics, 29(3): 339-349.

ISAACSON M S Q, MACKENZIE N G, 1981. Long-term distributions of ocean waves: a review [J]. Journal of the Waterway, Port, Coastal and Ocean Division, 107(2): 93-109.

LIANG B C, SHAO Z X, LI H J,,2019. An automated threshold selection method based on the characteristic of extrapolated significant wave heights [J].Coastal Engineering, 144: 22-32.

LIU G L, CUI K, JIANG S,, 2021. A new empirical distribution for the design wave heights under the impact of typhoons [J]. Applied Ocean Research, 111: 102679.

MAZAS F, HAMM L, 2011. A multi-distribution approach to POT methods for determining extreme wave heights [J]. Coastal Engineering, 58(5): 385-394.

SARTINI L, MENTASCHI L, BESIO G, 2015. Comparing different extreme wave analysis models for wave climate assessment along the Italian coast [J]. Coastal Engineering, 100: 37-47.

SOBEY R J, ORLOFF L S, 1995. Triple annual maximum series in wave climate analyses [J]. Coastal Engineering, 26(3/4): 135-151.

YOU ZJ, 2011. Extrapolation of historical coastal storm wave data with best-fit distribution function [J]. Australian Journal of Civil Engineering, 9(1): 73-82.

YOU Z J, 2012. Discussion of “a multi-distribution approach to POT methods for determining extreme wave heights” by Mazas and Hamm, [coastal engineering, 58: 385-394] [J]. Coastal Engineering, 61: 49-52.

YOU Z J, LORD D, 2008. Influence of the El Niño–southern oscillation on NSW coastal storm severity [J]. Journal of Coastal Research, 24(sp2): 203-207.

YOU Z J, YIN B S, 2006. Estimation of extreme coastal wave heights from time series of wave data [J]. China Ocean Engineering, 20(2): 225-241.

YOU Z J, YIN B S, JI Z Z,,2015. Minimisation of the uncertainty in estimation of extremecoastal waveheights [J]. Journal of Coastal Research, 75(sp1): 1277-1281.

UNIFIED APPROACH FOR ESTIMATION OF RETURN OCEAN WAVE HEIGHT

YOU Zai-Jin

(Centre for Ports and Maritime Safety, Dalian Maritime University, Dalian 116085, China)

Return wave height is an important parameter in the design of coastal and ocean engineering. Accurate calculation of design or return wave height is of enormous economic value and social value especially for deep-water harbors, ocean platforms, subsea gas and oil pipelines, and coastal nuclear power stations. However, there is no unique approach for calculation of return wave height, and results from different methods are significantly different. The present study is undertaken to analyze uncertainty in extrapolation of return wave height, minimize errors, and develop a unified methodology for calculation of return wave height. Based on long-term wave data continuously collected at Sydney permanent wave station on the coast of New South Wales in Australia, it is found that GPD-III and Weibull are the most suitable candidate distribution functions for extrapolation of return wave heights, newly derived formulation Eq.(8) will enable us to accurately estimate the shape parameters of the two distribution functions. The sample size and analysis method of extreme wave data are two important factors affecting the accuracy of return wave height extrapolated. Short wave records and the use of annual maximum method for analysis of extreme wave data could underestimate the return wave heights. The storm-by-storm method for analysis of extreme wave data and GPD-III or Weibull for extrapolation of return wave heights are highly recommended for coastal and ocean engineering design purposes.

return period; return wave height; extreme value analysis; GEV (generalized extreme value); GPD (generalized Pareto distribution)

* NSFC-山东联合重点基金, U1806227号。尤再进, 博士生导师, 教授, E-mail: b.you@dlmu.edu.cn

2022-03-08,

2022-04-05

TV139.2

10.11693/hyhz20220300051