不同方法对南黄海重现波高计算值的影响

赵嘉静，冯 曦，冯卫兵，李慧超

(河海大学 港口海岸与近海工程学院，南京 210098)

海岸工程的规划和建设对一个国家的经济发展起到带头作用。为了保证海岸基础设施的设计与建设的安全稳定，如码头与港口的设计与施工[1-2]，需要对波浪参数有确切的掌握。波高极值资料对于工程的设计与施工都不可或缺。

波高的长期变化趋势可以通过累积分布函数图和重现期分布来表征，此方法广泛应用于沿海极端波高的研究[3-5]。极值分析(EVA)在海洋、气象、金融、交通等领域的应用十分广泛，极值分析理论应用于极端波高的分析可追溯于Goda[6]、Mathiesen等人[7]和Menendez等人[8]，由于工程的需要，在中国近海海域已有基本的研究[9-12]。广义极值分布函数(GEV)经常被应用于极值的统计分布，最常用的方法是对于长期数据，采用每年取一个极大值的方法进行拟合，缺点是往往会忽略每年发生的其他极端情况。对此有2种解决办法：分块取极值法和超阈值(peak-over-threshold，POT)方法[13-14]。分块取极值法最常用的是采用取月极值波高的方法，即每个月取一个极值，尽可能统计到发生在全年的极端波高情况；POT方法即选取一个特定的阈值，对超过该阈值的值进行统计并分析。极端波高的长期变化趋势一直以来备受关注，Komar等人[15]使用NOAA浮标测量数据对美国西海岸等地也做了类似的研究。

南黄海地区位于中国江苏、浙江省沿岸，由于辐射沙洲的存在，近岸地形复杂(图1-a)，同时受到季风影响，多台风浪和风暴潮，出现极端水位情况较多[16]。现有的测站多位于沿岸且数据有缺失，不能完全表示出整个南黄海地区的长期波浪分布情况，因此不仅需要合适的极值计算模型，还需进一步分析数据长度对于重现波高计算的影响[13]。本文研究目的在于利用3种方法计算南黄海地区的百年重现波高空间分布，并进行比选，得到不同情况下的最优方法。同时，计算3种方法受到不同时间跨度的影响，为工程建设提供指导。

1 模型及极值计算方法介绍

1-a 区域地形 1-b 网格区域

1.1 模型设置

本研究采用SWAN模型进行模拟，此波浪模型已广泛应用于中国近海海域，模拟结果良好[17-19]。网格经度范围119°E～124°E，纬度范围31°N～35°N(图1-b)，共有10 429个节点和20 097个单元，精度范围从0.02°到0.11°(1 322～12 454 m)。利用欧洲天气预报中心(ECMWF)1979～2018年共40 a的风场数据作为初始条件，时间分辨率为6 h，空间分辨率为12 km。物理过程考虑了包括白帽、破波、底摩擦、波-波相互作用在内的波能耗散机制。SWAN模型的控制方程如下[20]

(1)

S=Sin+Swc+Sbrk+Sbot+Snl4+Snl3

(2)

式中：N为波作用量密度，方程左边第1项为动密度谱对时间的偏导，第2～5项分别为能量在x、y、σ、θ空间的传播，Ci(i=x、y、σ、θ)为各个空间的能量传播速度。方程右边S为源项，包括能量产生项(风能输入Sin)，耗散项(白帽耗散Swc，深度诱导引起的波浪破碎Sbrk，底摩阻耗散Sbot)及转化过程(四波相互作用Snl4, 三波相互作用Snl3)，σ为相对频率。

1.2 模型验证

表1 有效波高的验证参数

为了证明模拟结果的可信性，选取大丰(120.81°E，33.285°N)和蛎蚜山(121.568°E，32.147°N)2个浮标测点2013年的实测波高(Hs)数据进行对比。相关系数R，均方根误差RMSE及偏差Bias如表1所示，大丰测站吻合较好，蛎蚜山测站由于浮标放置在狭窄的潮道中，潮流的影响较明显，故在单纯考虑风场的情况下模拟准确性不高。图2为大丰测站2013年部分模拟值与实测值，可以看出除了实测数据有缺失的部分，模拟波高和实测数据拟合较好。图3中给出了大丰测站模拟值与实测值两者的线性比例的概率密度分布图，实测数据与模拟值贴近45°线，说明SWAN模型能够较准确地模拟本地区波高的时间变化。

图2 大丰测站模拟波高和实测数据对比

Fig.2 Comparison of simulatedHsand measuredHsat Dafeng station

图3 大丰测站模拟波高和实测数据线性比例的概率密度分布图

Fig.3 Linearly scaled probability density function(pdf)of observedHsversus modeledHsat Dafeng station

1.3 极值分布函数

1.3.1 GEV模型

在极值理论中，已经证明，对于足够长的独立和相同分布的随机变量序列，大小为n的样本的最大值可以拟合到广义极值(GEV)分布中[21]，该分布具有以下累积分布函数

(3)

μ、σ和ξ是位置、尺度和形状参数，与重现期T相对应的重现波高的求解可以使用以下方程

(4)

GEV经常被应用于极值的统计分布，最常用的方法是对于长期数据，采用每年取一个极大值的方法进行拟合，缺点是往往会忽略每年发生的其他极端情况。

1.3.2 基于POT方法的广义帕累托分布模型(GP模型)

对于超过一定阈值的数据，满足如下广义帕累托分布函数

(5)

与重现期T相对应的重现波高的求解如下

(6)

POT方法比每月取一个极值的方法更能反映出极端天气的发生情况，但是需要慎重给定取样时间间隔和阈值。时间间隔的选取应结合不同地区的极端天气发生情况，要保证尽可能多地考虑到极端天气发生次数，且不至于在同一次极端天气中选取多个极值。阈值的选取会直接决定取值的数量，过低的阈值会造成数据取样数量大进而导致模型预测结果的偏低，过高的阈值会造成分布参数计算结果的剧烈变动[22]。

对于时间间隔的选取，已经有学者研究发现，极端天气事件的发生间隔是随着不同地区、不同年限而相互独立的[22-23]，因此只能确定一个常用的合理值，常用3 d[13-14]。

1.3.3 极值分布模型的适用性

本文以南黄海地区划分的20个区域的40 a风浪推算数据为基础，采用基于取年极值与月极值的GEV分布模型和基于POT的GP分布模型，分析两种模型的适用性，并推算对应100 a重现期的重现波高。限于篇幅，本文仅选取了20个划分区域中具有代表性5个区块的进行详细说明：辐射沙洲北部区块N2、辐射沙洲中部区块R2、辐射沙洲边缘区块P4共3个区块(图1-b)。

通常采用累积分布函数(CDF)图和Q-Q图来判断数据与分布函数拟合的优劣。GEV分布的CDF图如图4所示，当各区域取年极值波高(图4-a～图4-c)和月极值波高(图4-d～图4-f)时，模拟值(虚线)与理论值(实线)均符合分布趋势，说明GEV分布适用于整个南黄海区域的波高拟合。在Q-Q图(图5)上可以看出年极值和月极值取样基本符合分布趋势，基本分布在45°线(虚线)附近，说明两者吻合良好。但在曲线尾部会有出入，预测曲线比实际偏低，即预测波高会小于实际值，这种现象在取月极值(图5-d～图5-f)时更明显。

4-a N2区域取年极值 4-b R2区域取年极值 4-c P4区域取年极值

4-d N2区域取月极值 4-e R2区域取月极值 4-f P4区域取月极值

图4 代表区块GEV分布的CDF图

Fig.4 The CDF plot of theGEVdistribution in the representative block

5-a N2区域取年极值 5-b R2区域取年极值 5-c P4区域取年极值

5-d N2区域取月极值 5-e R2区域取月极值 5-f P4区域取月极值

图5 代表区块GEV分布的Q-Q图

Fig.5 The Quantile-Quantile(Q-Q)plot derived fromGEVdistribution in the representative block

基于POT方法的GP模型中时间间隔选为3 d，对于模型中阈值的选取，则根据超额均值函数确定。图6所示为3个代表性区块的超额均值函数图。图7所示为3个代表性区块的平均剩余值图，图中实线代表形状参数值Xi，上下两条虚线为形状参数估计值的95%置信区间的上下限。在超额均值函数和形状参数满足稳定线性状态后来看，3个区域阈值结果分别为1.9 m、1.85 m和2.45 m。阈值确定之后，绘制满足超过阈值GP分布的拟合图，如图8所示，真实值由点组成，而黑色线表示分布趋势，纵轴表示GP分布函数。可见，选择最佳阈值后，在3个代表区块中，超出相应阈值的样本均符合GP分布，除了中间部分稍微有所出入，头尾两部分预测状况较好。

图6 代表区块超额均值函数图(时间间隔为3 d)

Fig.6 The Mean Excess Function(MEF)in the representative block(Time span is △t=3 d)

图7 代表区块95%置信区间的平均剩余值图(时间间隔为3 d)

Fig.7 The mean residual life plot with 95% confidence intervals in the representative block(Time span is △t=3 d)

图8 代表区块满足超过阈值的GP分布(时间间隔为3 d)

Fig.8 TheGPdistribution function plot meeting the peak over threshold in the representative block(Time span is △t=3 d)

图9 20个区块取年极值、月极值方法与POT方法计算的百年重现波高对比

综合之前采用取年、月极值的GEV分布和采用POT方法的GP分布来看，2种分布拟合均符合南黄海地区的波高分布趋势。在所研究的分布均适用的情况下，分别用取年极值、月极值、POT方法获得的数据，可以进一步计算得到对应分布函数下的重现波高。

2 结果

2.1 基于不同方法的百年重现波高空间分布

使用本文方法利用40 a的数据计算得到20个区块的100 a重现波高统计如图9所示。可见，采用GEV分布取年极值的方法比取月极值方法计算的百年重现波高明显高，这是因为年极值仅考虑了每年最大波高的情况，忽略了其他极端情况。POT方法采用的GP分布拟合出的百年重现波高和采用取年极值的GEV方法相差不大，互有高低，显然也比取月极值的GEV方法计算值高，这是由于POT方法平均每年取的极值数约为3个，并没有像取月极值那样每年取12个之多，忽略相对小的极端波高。

为了直观地表述南黄海地区的100 a重现波高的分布特征，图10给出了3种不同取极值方法所计算的百年重现波高分布图。其共同特征是，重现波高以辐射沙洲为中心向呈辐射状外围递增。不同点是，重现波高的较大值集中分布区域有所不同，取年极值方法的重现波高最大值出现于东北部深水区，POT方法位于东南部区域但范围较广，月极值方法大致分布区域则较为均匀。究其原因，可能是不同区域出现极端波况的频率和极端波高值各有不同。东北部地区更偏向于每年发生最极端情况，其他深水区更多发生次极端情况且发生次数更为频繁。

10-a 年极值 10-b 月极值 10-c POT

图10 三种方法分别计算的百年重现波高分布

Fig.10 The 100-year return wave height distribution calculated by three methods

图11 每个区域重现波高计算方法的最保守选择

在图11中给出了3种方法中，每个区域得到百年重现波高最大值时，即工程上最保守时采用的方法，可见在南北部地区取年极值方法计算的百年重现波高大，其他中部地区和辐射沙洲区域则基本以POT方法为大，说明这些地区在个别年份出现了多次高于正常年份的极端波况。

2.2 不同时间跨度对百年重现波高计算的影响

遵循工程上的取值思路，从最近的10 a数据开始，逐次增加数据长度，使用对应数据长度为10 a、20 a、30 a和40 a的2009～2018年、1999～2018年、1989～2018年和1979～2018年期间的数据。采用基于年、月极值的GEV分布和基于POT方法的GP分布计算了辐射沙洲内部的区域的百年重现波高，同时计算20～40 a相对10 a计算结果的增长率(表2～表4)。

表2 基于年极值的不同时间跨度的百年重现波高及增长率

表3 基于月极值的不同时间跨度的百年重现波高及增长率

表4 基于POT方法的不同时间跨度的百年重现波高及增长率

从不同取值长度计算的百年重现波高增长率看，取月极值的GEV分布计算结果对时间跨度的改变不敏感，而取年极值和POT方法均受制于时间跨度的选取。基于年极值的方法取近10 a的数据计算的百年重现波高最大，且不同数据长度计算出来的值差别很大，这是因为每年取一个极值的情况偶然性大，所记录的每一个极端波况对百年重现波高结果都会产生明显的影响。取年极值的方法计算的重现波高值随着取值年段长度增加而减小，即较长的数据集计算产生较小的重现波高，这与Mazas等[25]的发现相同，一方面可能是短期数据集的偶然性导致，另一方面可能是近10 a的极端天气情况发生较之前更为频繁。

取月极值的方法计算结果最为稳定，但是由于取样数量较多，影响相对小的波高也被考虑到，导致计算结果偏于不安全。POT方法的计算结果的稳定程度介于前两者之间。

取月极值方法计算的百年重现波高随着取值年段长度变化虽然稳定，但是明显偏小，取年极值计算结果最大，而POT方法计算的大小介于前二者之间。如果考虑到工程上偏于安全的情况，即在POT方法和取年极值方法中在不同地区选择各自最安全的重现波高值。

3 结论

本文利用第三代近岸波浪模型SWAN在中国南黄海地区模拟了长达40 a的波高数据，模拟波高和实测数据拟合较好。

基于年极值和月极值的广义极值分布(GEV)和超阈值取值方法(POT)的广义帕累托分布模型(GP)在南黄海地区适用。利用40 a模拟波高计算得到的百年重现波高在辐射沙洲北部地区计算的结果差别最大，采用月极值所得极值偏小，辐射沙洲南北外围地区采用年极值计算的重现波高大，其余地区则用POT方法计算出的结果较大。如果考虑到工程上偏于安全的情况，即在POT方法和取年极值方法中在不同地区选择各自最安全的重现波高值。

不同时间跨度的计算结果表明，取月极值的GEV分布计算结果对时间跨度的改变不敏感，而取年极值受之影响最大，POT方法介于两者之间。取月极值方法计算的百年重现波高随着取值时间跨度变化虽然稳定，但是在不同时间跨度下百年重现波高值明显偏小，取年极值的百年重现波高值最大，而POT方法计算的大小介于前二者之间。