基于P‑M谱和JONSWAP谱的实测风浪情况筛选和理论频谱拟合方法与流程

本发明涉及一种实测海浪情况筛选及针对实测海浪频谱的拟合方法，尤其涉及一种基于P-M谱和JONSWAP谱的实测风浪情况筛选及其频谱拟合方法。涉及专利分类号G06计算；推算；计数G06F电数字数据处理G06F19/00专门适用于特定应用的数字计算或数据处理的设备或方法。

背景技术：

在海浪能量的预测及应用中，研究理论频谱模型是一种重要的分析手段。理论频谱模型能够根据实测或假设得到的海浪要素，确定对应的海浪能量在各频率位置的分布，以此提供峰值能量的频率位置、各频率处的能量大小等信息。目前，通常使用的海浪理论频谱包括：Bretschneider(布氏)谱、Pierson-Moskowitz(P-M)谱、JONSWAP谱等。这些频谱通常都是根据特定海域的实测风浪数据拟合得到，虽然能够较好地表达一些海域的海浪能量分布情况，但是由于其修改余地小，对于部分海域的频谱拟合效果不够理想。同时，这些方法在拟合之前没有给出根据海浪要素的风浪筛选方法。

布氏谱是Bretschneider(1959)根据无因次的波高和波长联合分布，导出的一种适用于已充分成长或正处于成长阶段海浪的频谱。

P-M谱是Pierson和Moskowitz(1964)利用1955-1960年之间的北大西洋观测资料，找出已经充分成长的波浪数据，根据Kitaigorodskii相似定律，通过拟合无因次谱得到适合于充分成长海浪的风浪频谱。该谱即有实测数据的依据，又有一定的理论基础，至今依旧得到广泛使用。

JONSWAP谱是Hasselmann等(1973)通过对北海的观测资料进行分析而提出的。实际上，JONSWAP谱是在P-M谱的基础上多乘了一项谱峰因子γ，因此其谱形较P-M谱而言更集中于谱峰附近，谱峰形状更为尖突。这说明JONSWAP谱的提出者认为，频谱峰值附近应该包含着更多的波浪能量。

这三种谱在一定程度上都能够表达海浪的能量分布情况。但是由于其可修改项过少，导致在部分海况下拟合效果不好，甚至完全失灵。因此在实际海浪频谱的拟合中，需要一种可以根据实测数据，对理论频谱形式进行一定程度修改，以得到合适形式的方法。

技术实现要素：

本发明针对以上问题的提出，而研制的一种基于P-M谱的实测风浪情况筛选方法，包括如下步骤：

—获取目标海区内选定时间段内进行一定频率的海面高程；得到一组海面高程的时间序列x(t)；对所述的时间序列x(t)进行分块，得到多个时间分块；

—通过对每个时间分块进行傅里叶变换DFT，将每个时间分块转化成频谱函数S(f)；通过插值得到每个时间分块中频谱的实测频谱谱峰频率；

—根据所述的海面高程时间序列，求出所述每个时间分块中的三分之一大波高，作为时间分块中的有义波高Hs；

—根据P-M谱和有义波高Hs求得各所述的时间分块对应的P-M谱谱峰频率；

—计算所述实测频谱谱峰频率和对应的P-M谱谱峰频率的差值的绝对值；选取绝对值接近零的数值，作为风浪数据，完成实测风浪的筛选。

作为优选的实施方式，所述的相邻两时间序列分块的起始时间接近，每块内部持续时间长，两块之间允许出现重叠。

作为优选的实施方式，所述的x(n)的N点离散傅立叶变换DFT方法为：

海面高程的时域序列x(t)经过上述离散傅立叶变换，转为表达能量的频谱函数S(f)。

作为优选的实施方式，所述的P-M谱的表达式如下。

一种基于JONSWAP谱的理论风浪拟合方法，包括如下步骤：

—对JONSWAP谱增加可变参数，形成如下理论谱表达式；

—根据权利要求1中得到的风浪谱数据，采用最小二乘法对理论风浪谱进行双参数拟合得到k和γ值；

—根据得到的k和γ值分布散点分布，采用最小二乘法拟合k和γ值，得到使γ为k的二次函数；将k写为k＝f(γ)；

—将双参数的修改JONSWAP谱重新退化为单参数形式：

—针对γ进行单参数拟合，拟合方法为最小二乘法，求出所述每个时间分块对应的γ值；

—对各时间分块的γ值取平均数作为最终的γ值；根据拟合得到S(ω)表达式和所述最终的γ值得到最终的理论风浪谱表达式。

由于采用了上述技术方案，本发明公开的基于P-M谱和JONSWAP谱的实测风浪情况筛选和理论频谱拟合方法，具有如下优点：相对与现实海浪情况中只有部分为充分成长的风浪情况，本发明提出的实测频谱的谱峰频率和P-M谱的谱峰频率之间的关系可以很好地筛选风浪情况，为风浪谱拟合工作提供准确的实测频谱资料。通过为JONSWAP谱增加可修改参数，重新进行拟合处理，能够提高该理论谱的适用范围，实现对风浪情况下频谱的准确表达。

附图说明

为了更清楚的说明本发明的实施例或现有技术的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做一简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是根据P-M谱与实测谱的谱峰周期关系筛选风浪数据的计算方法的流程图。

图2为P-M谱谱峰频率与实测谱谱峰频率差值的绝对值分布散点图。

图3为参数k和γ分布散点及拟合线图。

图4为应用实例的典型拟合结果。

图5为基于JONSWAP谱的理论风浪拟合方法流程图

具体实施方式

为使本发明的实施例的目的、技术方案和优点更加清楚，下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚完整的描述：

如图1-5所示：

结合附图，对本发明的技术方案进行进一步详细的说明：

本风浪情况筛选及频谱拟合方法以海域定点实际测量为依据，根据实测资料中的海面高程时间域序列实现。

图1是P-M谱与实测谱的谱峰周期关系筛选风浪数据的计算方法的流程图。本实施例以已知海面高程时间序列来举例说明该方法可以包括以下几个步骤。

在步骤101中，将实测资料中的海面高程时间序列分块。其分块原则为：相邻两块的起始时间接近，每块内部持续时间长，两块之间允许出现重叠。例如：两块起始时间间隔可以设为2分钟，每块内包含30分钟的数据资料。

在步骤102中，对每个数据分块进行离散傅立叶变换(DFT)，以使资料中的时间序列转化成频域系列。x(n)的N点离散傅立叶变换DFT方法为：

海面高程的时域序列x(t)经过上述离散傅立叶变换可以转为表达能量的频谱函数S(f)。

在步骤103中，根据102中得出的频域形式，通过插值求出每个时间块中频谱的谱峰频率。

在步骤104中，根据海面高程的时间序列x(t)，求出每时间块中的三分之一大波高，记为有义波高Hs。

在步骤105中，根据P-M谱的表示：

和有义波高Hs求得每个时间块对应的P-M谱谱峰频率。

在步骤106中，绘制实测频谱谱峰频率和P-M谱谱峰频率差的绝对值散点图，如图2。能够发现在0附近有一组呈直线趋势的散点，将这组散点筛选出来，组成风浪数据。

根据挑选出的风浪资料，进一步完成理论风浪频谱的拟合，步骤如下：

通常的，在用未增加参数的理论表达式时发现，理论谱在高频处较实测谱低，需要抬高理论谱，通过增加可变参数可以更好的完成拟合，故作为优选的实施方式，对JONSWAP谱增加可变参数k，形成如下理论谱表达式：

根据风浪谱数据，对上述理论谱进行双参数拟合，拟合方法使用最小二乘法。

根据拟合结果，得到k和γ的值并画出它们的分布散点图，如图3。

根据k和γ的分布散点图拟合两者函数关系，拟合方法根据特定函数形式进行最小二乘拟合，通常可以使γ为k的二次函数，于是可以将k写为k＝f(λ)。在本实例中：该函数为拟合后二次函数的反函数形式。

将双参数的修改JONSWAP谱重新退化为单参数形式，如下：

然后针对γ进行单参数拟合，拟合方法为最小二乘法，求出每个时间块对应γ值。

最后对各时间分块的γ值取平均数作为最终的γ值。

根据拟合得到S(ω)表达式和γ值得到最终的理论风浪谱表达式。

应用实例得到的典型拟合结果图如图4，其中蓝线为该时间段内的实测频谱情况，红线为对应的理论频谱形状。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。