CN102930172A - 一种基于emd的海浪多尺度特征及波动参数的提取方法 - Google Patents

Abstract

本发明公开了一种基于EMD的海浪多尺度特征及波动参数的提取方法，所采用的方法是：读入定点测量的海浪时间序列或仿真的海浪序列；根据数据分析的需要，确定分析海浪尺度数；采用EMD方法，对海浪时间序列进行分解得到多个IMF序列；采用q次关联积分和差分盒维数法计算IMF分量的分形特征，依据矩估计法计算IMF分量的统计特征，根据谱方法计算IMF分量的波动参数。IMF分量的分形特征、统计特征及波动参数可以揭示海浪内部结构和外部表现特征。本发明可用于海浪的建模、分析及预报。

Description

一种基于EMD的海浪多尺度特征及波动参数的提取方法

技术领域

本发明属于物理海洋学技术领域，具体地讲就是涉及一种基于EMD的海浪多尺度特征及波动参数的提取方法。

背景技术

海浪是发生在海表面的一种波动现象，它与人类经济、社会活动密切相关。海浪的研究对于保障海上交通安全、有效地管理海洋环境和可持续利用海洋资源等具有重要意义。

研究表明，海浪具有随机性、非线性、非平稳性及分形的特征。从不同的角度和应用背景，已建立了多种海浪模型和分析方法。如以流体力学方程为基础，建立了深水表面波理论、小振幅波理论及有限振幅波理论，通常采用数值方法求解，计算效率和算法的稳定性是核心问题。双尺度海浪模型认为海浪是由小尺度波(毛细波)骑行在大尺度波(重力波)之上形成的，它在海洋微波遥感领域广泛应用，小尺度波对微波的散射是核心问题。而用三角函数、样条函数等几何曲线的线性组合来模拟海浪也是常用的模型之一，由于三角函数中的频率、幅度及相位与频谱之间的数学关系，采用经验海浪谱模型(如Phillips海浪谱、Pierson-Moskowitz海浪谱、Fung&Lee海浪谱及JONSWAP海浪谱)和随机过程理论相结合可以模拟出具有随机性的海浪，该方法不仅在计算机视觉中的海浪模拟应用广泛，而且在物理海洋学中也广泛应用。然而，该方法未能表现出海浪的非平稳和分形方面的特性。

发明内容

为了克服现有海浪模型和分析方法的不足，本发明提供了一种基于EMD的海浪多尺度特征及波动参数的提取方法。所述的识别方法包括以下步骤：

1.读入定点测量的海浪时间序列或仿真的海浪序列，记为x(m)；其中，m＝1,2,...,M表示测量或仿真海浪的离散间隔，M表示序列的点数。

2.根据数据分析的需要，确定分析海浪尺度数N，通常N取5到10之间的整数。

3.采用EMD方法，对海浪时间序列x(m)进行分解得到N个IMF序列，记为[c₁(m);c₂(m);...;c_N(m)]，包括以下过程：

(1).令n′＝1，c₀(m)＝0，r₀(m)＝x(m),m＝1,2,...,M，m为测量或仿真海浪的间隔；

(2).令k＝1，h₀(m)＝r_n′(m)，从序列r_n′_-1(m)中提取第n′个IMF序列，记为c_n′(m)；搜索序列h_k-1(m)的局部极大值和局部极小值，并用三次样条插值分别拟合局部极大值和局部极小值，得到序列h_k-1(m)的极大值包络线u_k-1(m)和极小值包络线v_k-1(m)；计算瞬时平均序列及序列h_k(m)＝h_k-1(m)-y_k-1(m)；计算序列h_k(m)的极值点个数和零交叉点个数，分别记为j₁和j₂；若序列h_k(m)不满足0≤|j₁-j₂|＜2或|y_k1(m)|＜ε条件，其中，ε为任意小的正数，通常ε取10^-1~10^-3之间的实数；则k＝k+1，对序列h_k(m)重复该步骤，直到序列h_k(m)满足0≤|j₁-j₂|＜2或|y_k-1(m)|＜ε条件；则h_k(m)为一个IMF序列，记为c_n′(m)＝h_k(m)；

(3).计算残余序列r_n′(m)＝r_n′_-1(m)-c_n′(m)，令n′＝n′+1，若n′＜N，转步骤(2)；否则，停止计算，得到了N个IMF序列，记为[c₁(m);c₂(m);...;c_N(m)]。

4.计算每个尺度下IMF的分形特征，记为

对第n个IMF序列c_n(m)，计算分形特征，记为

其中，n＝1,2,...,N，上标T表示转置，

为盒维数，D_q为扩展分形维数，f_α为扩展分形谱；

用r′×r′的盒子覆盖去第n个IMF序列c_n(m)，根据差分盒计数法，计算盒维数

其中，

为覆盖整个序列所需的盒子数；

扩展分形维数

其中，q为参数，r为尺度参数，G_q(r)为q次关联积分，其计算方法为，

其中，H(·)为阶跃函数，m′为嵌入维数，

r_i′j′为点集中任意一个元素与其它元素间的欧式距离， $r_{i^{\′}{j}^{\′}}=d({\tilde{C}}_n^{i^{\′}},{\tilde{C}}_n^{j^{\′}})=\sqrt{\underset{l=0}{\overset{m^{\′}-1}{\mathrm{\Σ}}}{[c_n(i^{\′}+l{r}^{\′})-c_n(j^{\′}+l{r}^{\′})]}^2},j^{\′}=\mathrm{1,2},...,\tilde{M},j^{\′}\≠i^{\′},$ 点集

是由序列IMF用如下的方法构造的 ${\tilde{C}}_n^{i^{\′}}={[c_n\left(i^{\′}\right),c_n(i^{\′}+r^{\′}),c_n(i^{\′}+2r^{\′}),...,c_n(i^{\′}+(m^{\′}-1)r^{\′})]}^T,i^{\′}=\mathrm{1,2},...,\tilde{M},$

扩展分形谱f_α＝qα-(q-1)D_q，其中，α为参数；则可得到N个IMF的分形特征，记为

5.计算每个尺度下IMF的统计特征，记为

对第n个IMF序列c_n(m)，采用矩估计法计算统计特征，记为

其中，均值 $\mathrm{\μ}=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{c}_n\left(m\right),$ 标准差 $\mathrm{\σ}=\sqrt{\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{[c_n\left(m\right)-\mathrm{\μ}]}^2},$ 偏斜度 $\mathrm{\γ}=\frac{\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{[c_n\left(m\right)-\mathrm{\μ}]}^3}{{\left(\sqrt{\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{[c_n\left(m\right)-\mathrm{\μ})]}^2}\right)}^3},$ 峰度

则可得到N个IMF的统计特征，记为

6.计算每个尺度下IMF的波动参数，记为

对第n个IMF序列c_n(m)，采用谱估计法计算波动参数，记为

其中，ω_peak为谱峰角频率，λ_peak为谱峰波长，H_sw为有效波高，

为平均周期；

首先对IMF序列c_n(m)作付里叶变换，得到频谱C(ω)，然后在频谱C(ω)曲线中搜索其最大值，即谱峰。那么谱峰所对应的角频率为谱峰角频率，记为ω_peak，根据海浪色散关系，谱峰波长

其中，g为重力加速度；有效波高

平均周期

则可得到N个IMF的波动参数，记为

本发明是基于EMD方法，将海浪进行多尺度分解，得到多个IMF分量。采用q次关联积分和差分盒维数法计算IMF分量的分形特征，依据矩估计法计算IMF分量的统计特征，根据谱方法计算IMF分量的波动参数。IMF分量的分形特征、统计特征及波动参数可以揭示海浪内部结构和外部表现特征。本发明可用于海浪的建模、分析及预报。

附图说明

图1为本发明的流程图。

具体实施方式

现结合附图对本发明进行详细的阐述。一种基于EMD的海浪多尺度特征及波动参数的提取方法，包括以下步骤：

1.读入定点测量的海浪时间序列或仿真的海浪序列，记为x(m)；其中，m＝1,2,...,M表示测量或仿真海浪的离散间隔，M表示序列的点数。

2.根据数据分析需求1，确定分析海浪尺度数N，通常N取5到10之间的整数。

3.采用EMD方法2，对海浪时间序列x(m)进行分解得到N个IMF序列，记为[c₁(m);c₂(m);...;c_N(m)]，包括以下过程：

(1).令n′＝1，c₀(m)＝0，r₀(m)＝x(m),m＝1,2,...,M，m为测量或仿真海浪的间隔；

(2).令k＝1，h₀(m)＝r_n′(m)，从序列r_n′_-1(m)中提取第n′个IMF序列，记为c_n′(m)；搜索序列h_k-1(m)的局部极大值和局部极小值，并用三次样条插值分别拟合局部极大值和局部极小值，得到序列h_k-1(m)的极大值包络线u_k-1(m)和极小值包络线v_k-1(m)；计算瞬时平均序列

及序列h_k(m)＝h_k-1(m)-y_k-1(m)；计算序列h_k(m)的极值点个数和零交叉点个数，分别记为j₁和j₂；若序列h_k(m)不满足0≤|j₁-j₂|＜2或|y_k1(m)|＜ε条件，其中，ε为任意小的正数，通常ε取10^-1~10^-3之间的实数；则k＝k+1，对序列h_k(m)重复该步骤，直到序列h_k(m)满足0≤|j₁-j₂|＜2或|y_k-1(m)|＜ε条件；则h_k(m)为一个IMF序列，记为c_n′(m)＝h_k(m)；

(3).计算残余序列r_n′(m)＝r_n′_-1(m)-c_n′(m)，令n′＝n′+1，若n′＜N，转步骤(2)；否则，停止计算，得到了N个IMF序列，记为[c₁(m);c₂(m);...;c_N(m)]。

4.计算每个尺度下IMF的分形特征，记为

对第n个IMF序列c_n(m)，采用差分盒维数法3和q次关联积分法4计算分形特征，记为

其中，n＝1,2,...,N，上标T表示转置，

为盒维数，D_q为扩展分形维数，f_α为扩展分形谱；

用r′×r′的盒子覆盖去第n个IMF序列c_n(m)，根据差分盒计数法，计算盒维数

其中，

为覆盖整个序列所需的盒子数；

扩展分形维数

其中，q为参数，r为尺度参数，G_q(r)为q次关联积分，其计算方法为，

其中，H(·)为阶跃函数，m′为嵌入维数，

r_i′j′为点集中任意一个元素与其它元素间的欧式距离， $r_{i^{\′}{j}^{\′}}=d({\tilde{C}}_n^{i^{\′}},{\tilde{C}}_n^{j^{\′}})=\sqrt{\underset{l=0}{\overset{m^{\′}-1}{\mathrm{\Σ}}}{[c_n(i^{\′}+l{r}^{\′})-c_n(j^{\′}+l{r}^{\′})]}^2},j^{\′}=\mathrm{1,2},...,\tilde{M},j^{\′}\≠i^{\′},$ 点集

是由序列IMF用如下的方法构造的 ${\tilde{C}}_n^{i^{\′}}={[c_n\left(i^{\′}\right),c_n(i^{\′}+r^{\′}),c_n(i^{\′}+2r^{\′}),...,c_n(i^{\′}+(m^{\′}-1)r^{\′})]}^T,i^{\′}=\mathrm{1,2},...,\tilde{M},$ 扩展分形谱f_α＝qα-(q-1)D_q，其中，α为参数；则可得到N个IMF的分形特征，记为

5.计算每个尺度下IMF的统计特征，记为

对第n个IMF序列c_n(m)，采用矩估计法计算统计特征，记为

其中，均值 $\mathrm{\μ}=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{c}_n\left(m\right),$ 标准差 $\mathrm{\σ}=\sqrt{\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{[c_n\left(m\right)-\mathrm{\μ}]}^2},$ 偏斜度 $\mathrm{\γ}=\frac{\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{[c_n\left(m\right)-\mathrm{\μ}]}^3}{{\left(\sqrt{\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{[c_n\left(m\right)-\mathrm{\μ})]}^2}\right)}^3},$ 峰度

则可得到N个IMF的统计特征，记为

6.计算每个尺度下IMF的波动参数，记为

对第n个IMF序列c_n(m)，采用谱估计法5计算波动参数，记为

其中，ω_peak为谱峰角频率，λ_peak为谱峰波长，H_sw为有效波高，

为平均周期；

首先对IMF序列c_n(m)作付里叶变换，得到频谱C(ω)，然后在频谱C(ω)曲线中搜索其最大值，即谱峰。那么谱峰所对应的角频率为谱峰角频率，记为ω_peak，根据海浪色散关系，谱峰波长

其中，g为重力加速度；有效波高

平均周期

则可得到N个IMF的波动参数，记为

Claims (6)

1.一种基于EMD的海浪多尺度特征及波动参数的提取方法。其特征在于，包含以下步骤：

步骤1：读入定点测量的海浪时间序列或仿真的海浪序列，记为x(m)；

步骤2：根据数据分析的需要，确定分析海浪尺度数N；

步骤3：采用EMD方法，对海浪时间序列x(m)进行分解得到N个IMF序列，记为[c₁(m);c₂(m);...;c_N(m)]；

步骤4：计算每个尺度下IMF的分形特征，记为

步骤5：计算每个尺度下IMF的统计特征，记为

步骤6：计算每个尺度下IMF的波动参数，记为

2.根据权利要求1所述的基于EMD的海浪多尺度特征及波动参数的提取方法，其特征在于，所述步骤1按如下过程进行：

读入定点测量的海浪时间序列或仿真的海浪序列，记为x(m)，其中，m＝1,2,...,M表示测量或仿真海浪的离散间隔，M表示序列的点数。

3.根据权利要求1所述的基于EMD的海浪多尺度特征及波动参数的提取方法，其特征在于，所述步骤3按如下过程进行：

1)令n′＝1，c₀(m)＝0，r₀(m)＝x(m)；

2)令k＝1，h₀(m)＝r_n′(m)，从序列r_n′_-1(m)中提取第n′个IMF序列，记为c_n′(m)；首先，搜索序列h_k-1(m)的局部极大值和局部极小值，并用三次样条插值分别拟合局部极大值和局部极小值，得到序列h_k-1(m)的极大值包络线u_k-1(m)和极小值包络线v_k-1(m)；计算瞬时平均序列

及序列h_k(m)＝h_k-1(m)-y_k1(m)；计算序列h_k(m)的极值点个数和零交叉点个数，分别记为j₁和j₂；若序列h_k(m)不满足0≤|j₁-j₂|＜2或|y_k-1(m)|＜ε条件，其中，ε为任意小的正数，则k＝k+1，对序列h_k(m)重复该步骤，直到序列h_k(m)满足0≤|j₁-j₂|＜2或|y_k-1(m)|＜ε条件；则h_k(m)为一个IMF序列，记为c_n′(m)＝h_k(m)；

3)计算残余序列r_n′(m)＝r_n′_-1(m)-c_n′(m)，令n′＝n′+1，若n′＜N，转步骤(2)；否则，停止计算，得到了N个IMF序列，记为[c₁(m);c₂(m);...;c_N(m)]。

4.根据权利要求1所述的基于EMD的海浪多尺度特征及波动参数的提取方法，其特征在于，所述步骤4按如下过程进行：

对第n个IMF序列c_n(m)，计算分形特征，记为

其中，n＝1,2,...,N，上标T表示转置，

为盒维数，D_q为扩展分形维数，f_α为扩展分形谱；

用r′×r′的盒子覆盖去第n个IMF序列c_n(m)，根据差分盒计数法，计算盒维数

其中，

为覆盖整个序列所需的盒子数；

扩展分形维数

其中，q为参数，r为尺度参数，G_q(r)为q次关联积分，其计算方法为：

其中，H(·)为阶跃函数，m′为嵌入维数，

r_i′j′为点集

中任意一个元素与其它元素间的欧式距离， $r_{i^{\′}{j}^{\′}}=d({\tilde{C}}_n^{i^{\′}},{\tilde{C}}_n^{j^{\′}})=\sqrt{\underset{l=0}{\overset{m^{\′}-1}{\mathrm{\Σ}}}{[c_n(i^{\′}+l{r}^{\′})-c_n(j^{\′}+l{r}^{\′})]}^2},j^{\′}=\mathrm{1,2},...,\tilde{M},j^{\′}\≠i^{\′},$ 点集

是由序列IMF用如下的方法构造的 ${\tilde{C}}_n^{i^{\′}}={[c_n\left(i^{\′}\right),c_n(i^{\′}+r^{\′}),c_n(i^{\′}+2r^{\′}),...,c_n(i^{\′}+(m^{\′}-1)r^{\′})]}^T,i^{\′}=\mathrm{1,2},...,\tilde{M},$

扩展分形谱f_α＝qα-(q-1)D_q，其中，α为参数；则可得到N个IMF的分形特征，记为

5.根据权利要求1所述的基于EMD的海浪多尺度特征及波动参数的提取方法，其特征在于，所述步骤5按如下过程进行：

对第n个IMF序列c_n(m)，采用矩估计法计算统计特征，记为

其中，均值 $\mathrm{\μ}=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{c}_n\left(m\right),$ 标准差 $\mathrm{\σ}=\sqrt{\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{[c_n\left(m\right)-\mathrm{\μ}]}^2},$ 偏斜度 $\mathrm{\γ}=\frac{\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{[c_n\left(m\right)-\mathrm{\μ}]}^3}{{\left(\sqrt{\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{[c_n\left(m\right)-\mathrm{\μ})]}^2}\right)}^3},$ 峰度

则可得到N个IMF的统计特征，记为

6.根据权利要求1所述的基于EMD的海浪多尺度特征及波动参数的提取方法，其特征在于，所述步骤6按如下过程进行：

对第n个IMF序列c_n(m)，采用谱估计法计算波动参数，记为其中，ω_peak为谱峰角频率，λ_peak为谱峰波长，H_sw为有效波高，

为平均周期；

首先对IMF序列c_n(m)作付里叶变换，得到频谱C(ω)，然后在频谱C(ω)曲线中搜索其最大值，即谱峰。那么谱峰所对应的角频率为谱峰角频率，记为ω_peak，根据海浪色散关系，谱峰波长

其中，g为重力加速度；有效波高

平均周期

则可得到N个IMF的波动参数，记为