基于EMD和修正高斯函数的时变平均风提取方法

技术领域

本发明涉及风场特性提取技术领域，尤其涉及一种基于EMD和修正高斯函数的时变平均风提取方法。

背景技术

时变平均反映了风速时程的非平稳特征，因而基于风速时程提取时变平均风成为构建非平稳风速模型的重要环节。实践中一般将风速的低频分量视作时变平均风，可采用多项式曲线拟合,滑动平均,高阶滤波和小波变换等方式从实测风速记录中提取，这些处理方式往往缺乏统一的定量提取准则，提取效果依赖于使用者的经验以及计算参数的选择，例如多项式的阶数，滤波器的类型和截止频率。近年来，Chen and xu(2004)将EMD技术应用于提取实测风速的时变平均风，将分解得到的余量和最后几阶本征模态函数(IMF)之和作为时变平均风。Chen er al.(2007)用EMD将传统平稳和非平稳风速模型进行对比研究.Chen and Wu(2012)结合EMD算法提出了用脉动风平稳度指标来确定时变平均风的方法。Su et al(2015)运用集成经验模态分解获得了两组全尺度非平稳下击暴流风速数据的时变平均风。然而需要说明的是，基于EMD提取时变平均风的效果同样取决于使用者的工程经验，并且EMD分解所得的余量并不一定就是时变趋势,仍然需采用相应判定条件进行识别。对于具有复杂变化趋势或者具有随机变化趋势的非平稳风速时程来说，如何确定相应判定条件并非易事。因此，研究有效的风速时变平均风提取方法成为建立非平稳风速模型的一项迫切的任务。

发明内容

为了解决上述现有技术中缺乏有效的风速时变平均风提取方法的缺陷，本发明提出了一种基于EMD和修正高斯函数的时变平均风提取方法。

本发明采用以下技术方案：

一种基于EMD和修正高斯函数的时变平均风提取方法，包括以下步骤：

A1、利用EMD将实测非平稳风速信号u(t)分解为若干个固有模态分量IMF和余量r

A2、将最后S-1个固有模态分量加上余量作为u(t)的时变平均风TVM

A3、将修正标准差σ＝k·α引入到标准高斯分布函数中，获得修正高斯函数 g[m

A4、获得用于量化实测脉动风速m

A5、获得多个不同的修正阶数k，针对每一个修正阶数k循环执行步骤A2 到A4，计算各修正阶数k在S取不同值时对应的修正高斯偏差系数p

A6、获取各修正阶数k对应的修正高斯偏差系数p

A7、获取极值点数量阈值λ，从各修正阶数k对应的备选项中筛选包含的目标极值点个数小于或等于λ的备选项作为最优时变平均风；若有多个备选项包含的目标极值点个数小于或等于λ时，从中选取目标极值点个数最多的备选项作为最优时变平均风；目标极值点为极大值或者极小值。

优选的，步骤A4中，修正高斯偏差系数的计算公式为：

优选的，极值点数量阈值λ的计算方式为：λ＝T

优选的，风速荷载周期T根据检测到的风速频率w计算获得。

优选的，步骤A1中，

优选的，步骤A2中，时变平均风的计算公式为：

优选的，步骤A2中，脉动风速的计算公式为：

本发明的优点在于：

(1)本发明中新定义了修正高斯函数，以修正高斯偏差系数来表征脉动的概率分布与修正高斯函数的吻合度，结合最小修正高斯偏差系数存在的唯一性，解决了趋势项最优判定的问题，从而为最优态的时变平均风的提取提供了判别条件，保证了最终获得的时变平均风的精确。

(2)本发明中，进一步结合极值点数量阈值λ对分析时距内非平稳风速的备选中的极大值点或极小值点数量进行约束，给出了提取最优的时变平均风的方法。

附图说明

图1为一种基于EMD和修正高斯函数的时变平均风提取方法流程图；

图2为数值算例及其理论趋趋势项；

图3为EMD分解结果；

图4a为新方法得到的趋势项和理论值；

图4b为指数趋势

图4c为二项式趋势

图5为实测120米高度处顺风向风速时程；

图6为120米高度处顺风向风速EMD分解imf分量；

图7a为residual作为趋势项分布图；

图7b为图7a对应的脉动概率分布图；

图8a为residual+imf10作为趋势项分布图

图8b为图8a对应的脉动概率分布图；

图9a为residual+imf10～8作为趋势项分布图

图9b为图9a对应的脉动概率分布图；

图10a为residual+imf10～5作为趋势项分布图

图10b为图10a对应的脉动概率分布图；

图11为47米高度处偏差系数随不同趋势项变化图；

图12为120米高度处偏差系数随不同趋势项变化图；

图13为280米高度处偏差系数随不同趋势项变化图；

图14为47米高实测风速数据随k值变换下的备选项图；

图15为120米高实测风速数据随k值变换下的备选项图；

图16为280米高实测风速数据随k值变换下的备选项图。

具体实施方式

本实施方式提出的一种基于EMD和修正高斯函数的时变平均风提取方法，包括以下步骤。

A1、利用EMD将实测非平稳风速信号u(t)分解为若干个固有模态分量IMF和余量r

具体的，本步骤中，

imf

A2、将最后S-1个固有模态分量加上余量作为u(t)的时变平均风TVM

本步骤A2中，时变平均风的计算公式为：

脉动风速m

以上公式(1)、(2)和(3)均为现有技术的直接应用。

A3、将修正标准差σ＝k·α引入到标准高斯分布函数中，获得修正高斯函数 g[m

其中，k为修正阶数，当k＝1时，修正高斯函数为高斯分布的概率密度函数；α是m

A4、获得用于量化实测脉动风速m

本步骤A4中，修正高斯偏差系数的计算公式为：

其中：f[m

A5、获得多个不同的修正阶数k，针对每一个修正阶数k循环执行步骤A2 到A4，计算各修正阶数k在S取不同值时对应的修正高斯偏差系数p

A6、获取各修正阶数k对应的修正高斯偏差系数p

即，针对任一个修正阶数k，获取与各S值对应的修正高斯偏差系数p

A7、获取极值点数量阈值λ，从各修正阶数k对应的备选项中筛选包含的目标极值点个数小于或等于λ的备选项作为最优时变平均风；若有多个备选项包含的目标极值点个数小于或等于λ时，从中选取目标极值点个数最多的备选项作为最优时变平均风；目标极值点为极大值或者极小值。

即，本步骤中，当仅有一个备选项包含的目标极值点个数小于或等于λ时，则选择该备选项作为最优时变平均风；当存在多个备选项包含的目标极值点个数均小于或等于λ时，则从包含的目标极值点个数小于或等于λ的备选项中获取包含的目标极值点最多的备选项作为最优时变平均风。

值得说明的时，本实施方式中，目标极值点取极大值时和取极小值时，时变平均风的提取结果相同。

本实施方式中，极值点数量阈值λ的计算方式为：λ＝T

以下，结合结构在简谐荷载作用下的动力特征，根据结构的动力放大因子D 给出风速荷载的周期T的具体计算方法。

D＝[(1-β

其中，ξ是阻尼比；β是动态荷载频率ω和结构自振频率ω

在结构风工程领域，超高层建筑及大跨结构风效应问题是其研究的重点。而这两类结构的基本自振频率往往很小，通常略高于0.1Hz。当结构阻尼比为 2％，并假定非平稳风速时变平均风的变化频率在0.01Hz以下时，基于公式可得到时变平均风作用下超高或大跨结构的动力放大因子D为：

由式(7)可知，当时变平均风的频率上限低于0.01Hz时，对应的动力放大因子接近于1，因此对应的时变平均风基本不会使结构产生动力响应。由风速的频率ω＝1/T＜0.01,可以得到相应的周期T≥100s。因此，当统计分析时距为T

如果非平稳风速时变平均风的局部目标极值点个数不超过数值λ，其通过式 (7)计算的D值约等于1，也就是说对应的时变平均风既能够体现非平稳风场的时变特征，又不会对结构产生动力作用。基于这一推论，对得到的系列时变平均风分量进行筛选获取最优时变平均风。然而需要说明的是，当结构阻尼比发生变化时，其对应的目标极值点个数限值将需要重新计算。

以下结合具体的实施例对本发明进行说明。

本实施例中，利用Matlab软件模拟生成了一均值为零的随机信号，并采用常幅值正弦函数模拟生成时变平均风分量，将两种信号叠加在一起形成非平稳信号：

y(k)＝1.5sin(0.025t)+2sin(0.05t)+randn (8)

式中y(k)为模拟得到的非平稳信号；randn表示均值为零的平稳随机序列。

上式(8)经离散生成的时域信号如(图2)，图中同时给出了时变平均风的变化趋势。对该非平稳信号进行EMD分解得到各IMF分量及余量如(图3)所示。对比两图可知，经EMD分解得到的余量与非平稳信号的理论趋势项差别较大，显然不能代表该信号的时变平均风。

采用本发明的方法对该非平稳信号进行处理，在分析时修正阶数k分别取1 至24的所有整数，并假定分析时距为10min，即T1＝10min＝600s。T≥100s，λ＝T

分析该组数据时，发现当k的取值可分为四个区间段：1～6、7～9、10～20、 21～24，同一区间段上的k值对应的时变平均风TVM

本实施例中，以k＝4、8、16、24为例。当k取4时，统计得到的备选项的极大值点数和极小值点数均超过了极值点数量阈值6，因此k取4时的备选项不满足要求。当k＝8、16和24时，对应的备选项的极值点数都满足条件，但k＝8 时对应的备选项的目标极值点数最多。即，k＝8、16和24时，对应的备选项中的极大值数量和极小值数量均小于或极值点数量阈值6。本实施例中，以极大值作为目标极值点，k＝8时对应的备选项包含的极大值多于k＝16、24时对应备选项包含的极大值。即确定当k＝8时对应的备选项为该实施例计算获得的最优时变平均风。且该最优时变平均风与时变平均风的理论趋势基本一致，验证了本发明的有效性。

为进一步验证本发明的有效性，利用下式构造不同类型趋势项的离散时域信号，应用本发明所提方法分解该信号，并将分解得到的时变平均风和预设的趋势项进行比较。

y(k)＝x(k)+randn k＝1,2,…,N

其中：y(k)为所构造的离散时域信号；randn为零均值的平稳随机序列；x(k)为预设的时变平均风。这里讨论当x(k)为线性形式、二次多项式形式和指数形式时本文提取方法的有效性。三种时变平均风的模型定义如下，

线性趋势，

x(k)＝a

二次多项式，

x(k)＝a

指数趋势，

x(k)＝ae

其中，线性趋势所取计算参数为a

从图4a/b/c的比较结果可见，对于线性、二项式和指数这三种设定的时变平均风类型来说，应用本发明的方法提取结果和理论值在两端处有一定的偏差，这是EMD分解本身存在端点效应引起的。除此之外，时变平均风的变化与理论值非常接近，说明此方法的有效性。

下面结合实例对该发明做进一步详细的说明，但是此说明不会构成对本发明的限制

验证实例：基于某气象塔的实测数据趋势项的提取

本实施例中，在325m的高塔上设置两组仪表。第一组包含三十个高分辨率 (0.1m/s)和0.05Hz采样频率的叶片风速仪(即机械风速仪)。它们分别安装在塔高8m，15m，32m，47m，63m，80m，100m，120m，140m，160m，180m，200m， 240m，280m和320m处，以测量平均风速数据，为了减少塔架对风速计记录的干扰，将每层的两个风速计安装在沿西北和东南方向长4米的钢悬臂梁上。另一组由三个超声风速仪组成，分别安装在塔高47m，120m和280m的铁架上，以10Hz 的采样频率记录三维均值和脉动风速。超声波风速仪安装在东北方向的钢悬臂梁，且传感器架设方向X轴朝东，Y轴朝北，Z轴朝上。即，本实施方式中，在塔高47m，120m和280m处，均同时安装有叶片风速仪和超声波风速仪。

第一步：利用EMD将实测非平稳风速信号分解为若干固有模态分量与余量。

实测风速数据来自在某次沙尘暴期间，架设在铁塔47米、120米、280米处的超声风速仪进行了全天24小时连续观测的数据。得到的风向以X轴朝东为正，Y轴朝北为正，Z轴朝上为正。

通过轮次检验法(Rotation testing method)分别挑选出47m,120m和280m 高度处的非平稳实测风速数据。然后用EMD理论对120米高度处顺风向十分钟实测风速进行分解，分解结果如图5所示。原信号被分解为10个imf分量和余量之和的形式，各imf分量如图6所示。

第二步：根据不同的高斯修正系数k，来计算偏差系数p

计算公式为：

计算该风速数据可能的趋势项。

当k取定值时，上式中，随着叠加的imf分量的增加，趋势项越来越接近原信号，越来越能表征信号的局部趋势特性。假设k去定值，n＝10；则图7所示为S＝1时对应的趋势项；图8所示为S＝2时对应的趋势项；图9所示为S＝4时对应的趋势项；图10所示为S＝7时对应的趋势项。

可见，通过叠加IMF的数量增大，对应脉动项的风速范围也会逐渐减小，而概率密度函数与X轴所围成的面积为1。结合图10、图9、图8、和图7可知，随着随着叠加的imf分量的增加，概率密度函数分布将变得越来越窄和尖。这从数学角度进行直观展示了计算不同趋势项下脉动的概率分布的均值和方差，得到如下表1。可以看到，随着IMF的叠加，对应脉动的均值和方差都会增加。

表1脉动分量的统计特征值变化规律

利用上面确定大致的高斯修正系数k由式

计算对应的偏差系数并保存偏差系数最小值对应的趋势项。

其中：f[m

本实施例中，取十分钟风速数据为一个计算单元，分别取6组塔高47m、120m 和280m高处的非平稳实测数据，编号为S1,S2,S3,S4,S5和S6.对每组数据进行EMD分解，计算每组风速数据在k＝4时的修正高斯偏差系数，分别参照图11到图13。

第三步：重复上一步，保存不同k值提取得到的趋势项。

本实施例中，分别取k＝4、8、12、16、20和24计算上述风速数据的备选项。同时以同样的方法计算其他高度实测风速数据在不同k值下的备选项，如图14～图16，由图可知，随着S值的递增即时变平均风TVM

第四步：根据条件来选择备选项。

在结构风工程领域，受风荷载影响较大的两类大型建筑，一类是大跨桥梁结构，一类是超高层建筑结构，这两种结构的一阶主频w

可以看到，当所取趋势项频率w

以上仅为本发明创造的较佳实施例而已，并不用以限制本发明创造，凡在本发明创造的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明创造的保护范围之内。