极值Ⅰ型分布求矩适线法在设计潮水位中的应用

许弘 程素珍 尹艳青

(1.兰州大学数学与统计学院，甘肃 兰州 730107；2.山东省水利科学研究院，山东 济南 250013；3.北京市水务工程建设与管理事务中心，北京 100036)

我国是海洋大国，拥有漫长的海岸线，山东省陆域岸线长度3122km，已建大陆海堤1184km，海堤与岸线比值为38%，且现有海堤破损严重，因此山东省已经明确提出要加大海堤建设，提高沿海重点地区的防风暴潮能力，而潮水位是确定海堤规模的关键性参数。《海堤工程设计规范》(SL 435—2008)(以下简称《规范》)提出在海岸地区设计潮水位频率分析的线型可采用极值Ⅰ型分布曲线[1-2]，但《规范》仅列出了频率、数据序列项数等相关系数的表格，通过查表求解相应设计频率的潮水位。查表法计算简单、应用方便，但其方法的局限性和成果的缺陷明显，在大中型工程中的应用受到较大的限制。适线法在水文频率分析中已经普遍和成熟应用，但主要基于皮尔逊Ⅲ型分布的应用[3-9]，有关极限Ⅰ型分布的工程应用很少，现有规范、相关文献也缺少极限Ⅰ型分布适线法的计算程序、适线方法和成果合理性判断标准，使该方法在潮水位频率分析应用中受到较大的限制。本文基于matlab计算软件，结合山东省荣成市石岛湾海堤设计，开发了极值Ⅰ型分布求矩适线法程序，并与规范表格系数查算、理论曲线成果进行比较，验证了该法的合理性、便利性，以期对类似工程提供借鉴。

1 极值Ⅰ型分布原理和适线求解方法

根据《规范》规定，设计潮水位计算应采用频率分析法，在海岸地区频率曲线的线型可采用极值Ⅰ型分布。

潮水位观测站历年最高潮水位观测值为h1，h2，…hn，观测值呈极值Ⅰ型分布，则某一设计潮位hP对应的设计频率值P为

P=P(h≥hP)=1-F(hP)=1-e-e-α(hP-β)

(1)

σ——潮水位序列均方差。

潮水位序列均值和均方差采用矩法进行估算，见式(2)、式(3)。

(2)

(3)

式中hi——序列第i年的年最高潮水位；

n——年最高潮水位序列项数。

对式(1)两边取两次对数，把α、β化简代入得式(4)，得

(4)

令：λPn=-(0.45+0.78×ln(-ln(1-P)))，式(4)转化为：

(5)

根据极值Ⅰ型分布理论对实测经验频率曲线进行适线拟合须借助数学工具。在matlab语句中，有evcdf、evpdf、evinv等函数，分别表示极值Ⅰ型分布的累计函数、密度函数和逆函数，可以通过P=evcdf(-hP,-b,1/a)、hP=-evinv(k,-b,1/a)实现对给定设计频率和相应设计频率的设计水位进行交互计算实现，实现频率适线法。

2 资料分析

2.1 基本资料

图1 设计潮水位极值Ⅰ型分布适线

2.2 表格法计算设计潮水位

根据《规范》，对石岛海阳观测站35年连续的年最高潮水位序列的均值、均方差及年频率为0.1%、0.2%、0.5%、1%、2%、4%、5%、10%、25%、50%、75%、90%、95%、97%、99%、99.9%的最高潮水位分别按式(2)、式(3)和式(5)计算，系数λPn按《规范》表A.0.2查取，不同频率下的设计潮水位见表1，不同频率及相应潮水位点绘在海森概率格纸上并连成曲线(见图1)。

表1 不同频率下设计潮水位计算成果 单位：m

2.3 理论计算设计潮水位

通过理论计算知，曲线形状由均值、均方差两个参数决定。调整均值仅使频率曲线平行上下移动，偏移经验点，通过调整均值进行适线法是不合适的(见图2)。但当验潮零点高程的标准不一致时，可以通过调整序列均值实现。

图2 调整均值的极值Ⅰ型适线

设定理论频率曲线P=42.97%为分界点，该点的设计潮水位等于序列均值，高于该频率的设计潮水位均大于序列均值，小于该频率的设计潮水位均小于序列均值。由于理论曲线以分界点旋转，因此通过调整均方差，能够适配出合理的曲线。

2.4 求矩经验适线法计算设计潮水位

求矩经验适线法是在一定的适线准则下，求解与经验点据拟合最优的频率曲线的统计方法，可避免因观测序列短而产生计算结果可靠性差等问题。采用矩法对观测数据序列进行统计，求出一组参数作为初值计算理论曲线，根据理论频率曲线与经验频率点据的配合情况，通过经验判断调整参数，选定一条与经验点据拟合良好的频率曲线。经验频率曲线选择时，尽可能选择拟合点据好、结果相对保守的曲线作为经验适线。适线要尽可能照顾点群的趋势，使频率曲线通过点群中心，点线配合不佳时，可重点考虑中上部点据，并应分析经验点据的精度，使曲线尽可能接近或穿过比较可靠的点据。

通过初估、适线和综合对比分析，可以得到比较合理的、能够满足海堤工程设计要求的潮水位经验频率曲线。通过选定的经验适线，根据极值Ⅰ型理论计算相应设计频率的潮水位。借助经验频率曲线予以延长，推求相应稀遇频率的设计值。

适线法计算通过matlab软件中的evinv、evcdf函数实现，在软件环境下编制程序，能够实现数据的统计、分析、绘图和曲线拟合同时完成，采用norminv函数生成海森概率格纸。根据理论曲线和经验点据的适配情况，调整均方差值为0.9σ、1.1σ、1.2σ、1.4σ倍的初值，进行适配，曲线见图1，根据每条适配线计算年频率为0.1%、0.2%、0.5%、1%、2%、4%、5%、10%、25%、50%、75%、90%、95%、97%、99%、99.9%的最高潮水位，结果见表1。

3 计算成果和方法的合理性选定

表格法是《规范》推荐采用的设计潮水位计算方法，根据极值Ⅰ型分布，编写了P-n-λPn参数表，计算简单，通过手算即可完成。极值Ⅰ型分布为全样本理论曲线，为避免观测系列较短对拟合计算结果的影响，采用观测值序列项数多少调整λPn值，减小因系列短而产生的误差。由图1知，表格法计算成果和1.15倍均方差初值拟合曲线一致，说明表格法由通过调整均方差初值计算生产，是适线中的一条经验频率曲线成果。

表格法把复杂的计算过程简化为查表计算，并把样本序列长度的影响考虑在内，其计算成果能够满足小型工程的需要，但规范表格是在特定序列项数n、频率P下查算相应λPn值，对于表格中没有列出的n、P对应的λPn值，只能通过内插求解，存在一定的误差，同时由于成果是根据固定适配线确定，其合理性、与经验点据的适配性存在较大的局限。

表格法计算的设计潮水位受序列长度的影响较大，观测序列的长度不同，其λPn值不同，时间越短，其系数越大。为验证计算成果受时间序列长短的影响程度，选取该工程20年、25年、30年、35年最大潮水位系列资料，同时计算20年一遇设计潮水位，其计算成果见表2，根据计算成果绘制的序列项数n与潮水位关系曲线见图3。从图3可知，序列长度为25年的计算设计潮水位值最大，说明表格法并不能保证序列长度越小选定的成果越保守，这是该方法的局限性。

表2 表格法计算20年一遇潮水位成果 单位：m

图3 序列项数n与潮水位关系曲线

适线法根据极值Ⅰ型分布和观测资料统计参数进行理论分析和经验判断，通过调整参数，适配出与经验频率点据拟合最好的曲线，得到满足工程设计要求的不同设计频率下的潮水位设计值，计算成果包含了表格法。它是基于matlab软件环境编制的程序，能够完成数据的统计分析，多条适配曲线拟合，图表自动生成，计算速度快，功能强大，能够把经验点据、规范表格计算值、适配曲线同时绘制在同一张海森概率格纸上，并把成果汇总于同一表格中，便于设计人员比较和分析，确定适合设计要求的设计潮水位计算值。

4 结 论

本文采用极值Ⅰ型分布，结合工程设计潮水位实例计算，得出以下结论：

a.规范表格法简单依靠参数表查算完成，过程机械，结果受资料序列长度影响，且不能确定与经验频率点据的拟合程度，具有较大的局限性。

b.表格法计算成果与特定倍数均方差适线成果一致，能够通过适线法计算完成。

c.适线法基于matlab软件的极值Ⅰ型分布函数完成，计算理论符合《规范》要求，计算功能强大，能够把经验点据、《规范》表格计算值、适配曲线同时绘制在同一张海森概率格纸上，并把成果汇总于同一个表格中，便于设计人员比较和分析，确定适合设计要求的设计潮水位计算值。