一种基于MEMD与SRM的单样本非平稳风速模拟方法与流程

本发明涉及多点非平稳风速模拟，具体涉及一种基于memd与srm的单样本非平稳风速模拟方法。
背景技术：
：随着科技与经济的高速增长，高压、超高压输电逐渐成为各国供电的主要模式。中国对能源特别是电能有着更加迫切的需求，电力供应已经进入“大电网、大机组、西电东送、南北互济、全国联网”的新时代。输电塔线体系承担着高压电能运输的重要使命，是事关国家稳固发展的重要工程。高压输电线塔具有结构高、导线截面大、线圈多、跨距长、负荷大、柔韧性高的特点。高压输电线塔对风力作用非常敏感，在风力作用下容易对其造成极大的风致动力响应，导致疲劳和失稳破坏，甚至引起输电线塔的倒塌破坏。输电线塔在极端风(雷暴风、龙卷风)产生的风荷载作用下，极易因疲劳和失稳破坏而导致倒塌。研究表明，近年来雷暴风造成的破坏在强风破坏中占了很大的比例。据统计，在美国、澳大利亚、南非等地与气候相关的输电线塔破坏约80％是由雷暴风造成的。因此为了进一步提高输电线塔的安全性，有必要使其在极端风作用下的抗风设计更加精准可靠。为了保证抗风设计的准确性需要获取大量的极端风样本。然而难以通过直接测量获得大量极端风样本，实际上往往只能测得少数几个样本或单样本。因此可靠获取大量风速样本对高耸结构的风振分析至关重要。但由于极端强风具有非平稳性，一般采用非平稳过程来描述。考虑到结构和气动力的非线性，常采用时频分析来描述极端风等非平稳信号。目前，模拟非平稳信号通常有两类方法，第一类是时间序列法。该类模拟方法计算效率高，但模型定阶和时变参数的识别通常十分复杂。第二类方法是基于演化功率谱的谱表示法。与时间序列方法相比，其更直接更准确。并且随着快速傅里叶变换的应用，基于演化功率谱的谱表示法模拟效率得到了显著提高。然而实际上通常智能获取单样本用来评估时变相干函数或者演化功率谱。因此不确定性准则会导致时域和频域的平滑估计，造成模拟结果偏离原始信号。除了上述模拟方法之外，基于单样本(或基于数据)的时频分析的模拟方法也引起了广泛关注。基于单样本的模拟方法利用数据中的瞬时特征，为多点非平稳而模拟提供直接且具有统计意义的算法。wen采用emd和希尔伯特黄变换来模拟多点非平稳信号。该方法假设任意两个不同点之间的初始相位角之差是常值，以便建立互协方差。然而这种假设可能无法真实反映多点风速之间的相关性。此外，emd可能导致模态不一致性，即其分解多点信号产生数量不同的imf，并且对于不同的信号点在相同尺度上的频带差异也很大。为了解决这个问题，wang等人采用平稳小波变换和希尔伯特变换来模拟多点风速。首先采用本征正交分解来解耦瞬时频率矩阵。然后，考虑多点风速的相关性，添加服从高斯分布的随机频率以形成最终瞬时频率矩阵。在上述两种模拟方法中，假设任何两个不同风速点之间的初始相位角之差为常数，以便使得模拟风速的互协方差与目标值一致。但是这种假设不合理。wang和wu等利用自适应离散小波包分解和希尔伯特变换进行模拟。通过瞬时相位差，相关系数和相干函数之间的关系，提出了基于假定相干函数和单点风速的多尺度空间相关性嵌套的模拟方法。然而这些模拟方法存在两个问题：希尔伯特变换可能产生巨大的频率突变，可能产生大量的没有意义的负频率；这些方法假设不同点同一尺度的初始相位角之差为常数。基于此，需要进一步探究能够基于单样本实测数据获取大量风速样本的方法。技术实现要素：本发明针对现有技术存在的问题提供一种能够基于单个样本生成更多样本的一种基于memd与srm的单样本非平稳风速模拟方法。本发明采用的技术方案是：一种基于memd与srm的单样本非平稳风速模拟方法，包括以下步骤：步骤1：采用多元经验模态分解memd方法将风速样本信号分解产生不同尺度的本征模态函数imfs；步骤2：根据am/fm分解步骤1得到各点的imfs的瞬时频率和瞬时幅值；步骤3：根据步骤2中各点的瞬时频率的自谱和互谱，得到瞬时功率谱密度矩阵；步骤4：根据步骤3中的功率谱密度矩阵修正为正定谱矩阵；步骤5：将步骤4得到的正定谱矩阵进行cholesky分解；步骤6：引入随机初始相位角即可获取模拟的随机风速。进一步的，所述步骤3中假设多点信号分解得到的瞬时频率的自谱和互谱中的频率存在的差异极小，忽略其差异；谱密度矩阵中的频率取加权频率，得到统一频率后的瞬时功率谱密度矩阵。进一步的，所述步骤4中修正算法如下：s41：通过特征值计算瞬时功率谱密度矩阵的特征向量矩阵和对角特征值矩阵；s42：将特征向量矩阵中的负值全部替换为正值即可得到所需正定谱矩阵。进一步的，所述模拟的随机风速用于高耸结构风振分析。进一步的，所述步骤1的分解过程如下：p点的多元信号x(t)＝[x1(t)，x2(t)，...，xp(t)]t分解过程如下：s11：在n-1维单位超球面上采样生成p个点的hammersley(哈莫斯利)点集；s12：计算信号x(t)沿着所有方向向量的映射向量，用表示；其中，θn为向量空间；s13：找到映射向量的极大值、极小值所对应的时间量s14：对采用三次样条插值求多元信号的包络线s15：求得包络均值s16：令c(t)＝x(t)-μ(t)，当c(t)满足memd的imf的终止条件时，重复步骤s12～s15获得更高阶的imfs。进一步的，所述步骤2的分解过程如下：假设y(t)为经memd分解产生的单组分信号；s21：获取单元信号y(t)的绝对值的极值；s22：采用三次样条拟合其绝对值的极值，得到经验包络e1(t)，用y(t)除以e1(t)得到y1(t)：y1(t)＝y(t)/e1(t)继续迭代至yn(t)≤1；y2(t)＝y1(t)/e2(t)；…；yn(t)＝yn-1(t)/en(t)此时单组分信号y(t)的fm，am分别为：a(t)＝a(t)＝y(t)/f(t)＝e1(t)e2(t)…en(t)至此原始信号y(t)被唯一分解为：其中，为相位角；迭代后，瞬时幅值近似等于该am函数；s23：对fm求导，对导函数f′(t)再次进行am/fm分解，f′(t)可表示为：通过对比上式，因此瞬时频率为：本发明的有益效果是：(1)本发明采用多元经验模态分解memd和幅值调制和频率调制am/fm分解算法，其中memd可以分解多点信号生成相同数量的imf；此外，还具有更强大的滤波功能，使相同尺度的imf位于几乎处于相同的频率带，并减少模态混叠；am/fm分解算法分解imf产生瞬时幅度和瞬时频率，有效抑制了负频率；能够显著提高风速模拟精度。(2)本发明方法得到的模拟风速能够很好的保存了原始风速的湍流强度特性和相关性。(3)本发明方法保持原始不同点风速的空间相关性，模拟方法更加简单合理，能够有效模拟多点非平稳信号，可用于高耸结构的抗风设计。附图说明图1为本发明流程示意图。图2为本发明实施例中采用的原始风速信号图。图3为本发明实施例中经memd分解得到的时变均值信号图。图4为本发明实施例中经memd分解得到的imf图。图5为本发明实施例中经am/fm分解得到的瞬时频谱。图6为本发明实施例中瞬时频谱相应的瞬时互谱和统一频率后的瞬时互谱。图7为本发明实施例中得到的模拟风速信号图。图8为本发明实施例中原始风速和模拟风速的时变湍流强度对比示意图。图9为本发明实施例中原始风速和模拟风速的psd对比示意图。图10为本发明实施例中原始风速和模拟风速的cnai对比示意图。图11为本发明实施例中模拟风速的瞬时频谱。图12为本发明实施例中模拟风速的瞬时互谱。图13为本发明实施例中模拟风速的k-l散度图。图14为本发明实施例中模拟风速的协方差图。图15为本发明实施例中模拟风速的互协方差图。具体实施方式下面结合附图和具体实施例对本发明做进一步说明。本发明中出现的符号如下：memd为多元经验模态分解，srm为谱表示方法，am/fm为幅值调制和频率调制；imfs为本征模态函数。如图1所示，一种基于memd与srm的单样本非平稳风速模拟方法，包括以下步骤：步骤1：采用多元经验模态分解memd方法将风速样本信号分解产生不同尺度的本征模态函数imfs；假设已知一个潜在非平稳信号x(t)＝[x1(t)，x2(t)，...，xn(t)]t的一个单样本经memd分解后为其中，xj(t)为第j个原始信号，n为原始信号总数，cjp(t)为emd的imf或swt的第p层细节分量，n为分解的总层数，rj(t)为残量或逼近分量，t为时间。步骤2：根据am/fm分解步骤1得到各点的imfs的瞬时频率和瞬时幅值；根据基于memd和改进的am/fm分解方法的时频分析框架，基于任何信号的基于memd的瞬时频谱为：其中：ajp(t)为第p阶imf得到的瞬时幅值，ω为频率，δ(·)为dirac函数，为信号第p阶imf得到的瞬时频率。步骤3：根据步骤2中各点的瞬时频率的自谱和互谱，得到瞬时功率谱密度矩阵；以三点非平稳风速的模拟为例，其第m层谱密度矩阵sm(ωm，t)为：其中，ajp(t)为第p阶imf得到的瞬时幅值，akp(t)为第p阶imf得到的瞬时幅值，为，为第p阶imf得到的瞬时频率，k为第k个信号，j为第j个信号。由于多点风速同阶imfs的瞬时频率不是完全一致的，因此上述矩阵实际上并不存在。但是因为memd的强大的滤波特性使得同阶imfs分解产生的频率带只存在极小的差异，因此我们可以假设多点信号分解得到的瞬时频率谱与互谱中的频率存在极小的差异，为了统一谱密度矩阵中的频率，其中频率取加权频率：统一频率后的谱密度矩阵为：其中：步骤4：根据步骤3中的功率谱密度矩阵修正为正定谱矩阵；通常第m层的谱密度矩阵为正定矩阵，但由于模型的计算误差等因素矩阵存在负的特征值。通过以下算法调整其为正定矩阵。通常统一后瞬时谱矩阵是非正定非负定矩阵。然而由于数值误差，存在一些极小的负特征值。通过以下算法调整其为正定矩阵。首先，通过特征值分析计算矩阵的特征向量矩阵和对角特征值矩阵然后，将中的负值全部替换为小的正值。本发明中将其设置为0.001，调整后的瞬时频谱矩阵为：其中：为经过调整后的特征值对角化后的矩阵。步骤5：将步骤4得到的正定谱矩阵进行cholesky分解；调整后第p层的谱密度矩阵经cholesky分解为：其中：为cholesky分解得到的下三角矩阵。步骤6：引入随机初始相位角即可获取模拟的随机风速。引入随机的初始相位角，信号xj(t)可重构为：其中：φmp为引入的初始相位角。模拟非平稳信号的均值μj(t)为：其中，自协方差为：方差为：互协方差为：因为同层imfs之间存在相关性，不同层间独立，仅当m1＝m2＝m，且n1＝n2＝n时存在互协方差。因此有：互协方差仅取决于分解得到的谱和实测信号中获取的确定性相位角差，从而可以看出随机初始相位角的假设不是必须的。因此基于单个样本的多变量随机过程模拟更加简单和合理。上述方法适用于多点非平稳风速模拟。本发明中涉及的算法模型包括：1、多元经验模态分解memdmemd通常采用低差异hamersley序列(准蒙特卡罗采样)以便保持更好的一致性。分解一个p点的多元信号x(t)＝[x1(t)，x2(t)，...，xp(t)]t的过程如下：s11：在n-1维单位超球面上采样生成p个点的hammersley点集；s12：计算信号x(t)沿着所有方向向量的映射向量，用表示；其中，θn为向量空间；s13：找到映射向量的极大值、极小值所对应的时间量s14：对采用三次样条插值求多元信号的包络线s15：求得包络均值s16：令c(t)＝x(t)-μ(t)，当c(t)满足memd的imf的终止条件时，重复步骤s12～s15获得更高阶的imfs。2、am/fm分解am/fm分解通过唯一分解imf为am函数a(t)和fm函数cosθ(t)来实现。假设y(t)为经memd分解产生的单组分信号，则信号y(t)的瞬时幅值ia和瞬时频率if可以通过改进的am-fm分解计算得到，其具体的计算步骤如下：s21：获取单元信号y(t)的绝对值的极值；s22：采用三次样条拟合其绝对值的极值，得到经验包络e1(t)，用y(t)除以e1(t)得到y1(t)：y1(t)＝y(t)/e1(t)(1)理论上y1(t)的极值为一个定值并且e1(t)等于a(t)，因此同时满足y1(t)＝cosθ(t)和y1(t)≤1。但实际上e1(t)的值一般不等于a(t)，特别是在a(t)快速变化时。因此应继续迭代至yn(t)≤1；y2(t)＝y1(t)/e2(t)；…；yn(t)＝yn-1(t)/en(t)(2)此时单组分信号y(t)的fm，am分别为：a(t)＝a(t)＝y(t)/f(t)＝e1(t)e2(t)…en(t)至此原始信号y(t)被唯一分解为：其中，为相位角；迭代后，瞬时幅值近似等于该am函数；通常，两次或三次迭代能够满足要求；所以本发明中采用线性归一化和三次迭代；s23：对fm求导，对导函数f′(t)再次进行am/fm分解，f′(t)可表示为：通过对比上式，因此瞬时频率为：通过两次am/fm分解和一次求导得到ia和if，能够满足bedrosian定理；由于没有涉及到希尔伯特变换，能够满足nuttall定理；因此该算法得到具有物理意义的if和ia，并且没有负频率，能够显著提高风速模拟精度。目前基于单样本的模拟方法大致分为三种，即wen和gu(2004)，wang等人提出的方法(2013；2014)和wang和wu(2018)；前两种模拟方法直接基于多点样本，而最后一种模拟方法基于单点实测样本和规定的相干函数。本发明中针对的是基于多点的单样本模拟，介绍前两种模拟方法：假设已知一个潜在非平稳信号x(t)＝[x1(t)，x2(t)，...，xn(t)]t的一个单样本通过emd或swt分解实测样本j＝1，2，…，n的产生频率分布由高到低的子序列。式中cjp(t)为emd的imf或swt的第p层细节分量；n是分解的总层数；rj(t)是残量(emd)或逼近分量(swt)。获得单组分信号后，通过希尔伯特变换计算ia和if：其中ajp(t)和ωjp(t)分别是在第j层分解产生的ia和if。希尔伯特变换不仅在生成ia和if时计算冗余，而且还产生大量没有物理意义的负频率，应该减少这些负频率以获得更具物理意义的结果。wenandgu(2004)提出了以下公式来模拟多元随机过程：式中表示模拟得到的随机信号；φj，k是服从[0，2π]均匀分布随机信号的相位角。和互协方差为：式中e{·}表示数学期望；θjp(t)＝∫ωjp(t)dt是瞬时相位。假设和同尺度的imfs初始相位角的差值为常数，如下：φ1p-φjp＝ψjp，j＝2，3，…，n；p＝1，2，…，n(11)其中ψjp是个常数，则算式(11)就化简为：从上式可以看出，模拟随机信号的互协方差取决于同尺度imf产生的幅值以及初始相位和随机相位的之差。初始相位差异反映了实测样本的确定的空间相关性，只有随机相位能够模拟出随机过程的不确定的空间相关性。因此，模拟随机信号容易受到随机相位差的影响，这很难直接从样本中获取。为了将上述方法应用于多点模拟，理论上需要知道协方差函数的目标值，从而才能获取不同点的初始相位角。实际上大多数情况观测样本只有一个，一定的假设和近似是必需的。wen在模拟地震动时假设初始相位角之差为定值，ψjp＝ψkp＝0。为了增大模拟随机信号的的变异性，wang等采用pod将第p层相关的瞬时频率ωp(t)＝[ω1p(t)，ω2p(t)，…，ωnp(t)]t分解为一系列不相关子过程，即：式中φptφp＝i并且是标准化的正交函数，是展开序列。将视作服从高斯分布随机频率，截断其高阶项，第p层随机瞬时频率可表示如下：其中是第j个信号的用替换式(10)中ωjp(t)：此时，和之间的互协方差可表示为：式中当j＝k时，和独立；j≠k时，和不独立，但同尺度的随机相位独立。假设初始相位角关系如下：φ1p-φjp＝ψjp，j＝2，3，…，n；p＝1，2，…，n(17)式(17)可以化简为因此，互协方差取决于同尺度的幅度函数，和以及随机初始相位角。基于memd和am/fm分解方法的时频分析框架如下：对于多点j＝1，2，…，n，经memd分解后为：然后可以采用希尔伯特变换或其他方法来获得单组分信号的ia和if。采用am/fm分解方法计算单组分信号的ia和if。式(20)可以改写为根据基于memd和改进的am/fm分解方法的时频分析框架(huangetal.2016)，任何信号的基于memd的瞬时频谱为：式中δ(·)为dirac函数，；为信号第p阶imf得到的瞬时频率。由于if是根据不同点的信号计算的得到，因此在相同尺度的不同点imf分解得到的的if通常是不同的；所以式(22)不能直接扩展到瞬时互谱，由于memd的强大的滤波特性使得处于同尺度的imf位于几乎相同的频带中；这些if之间的差异很小。所以来自两个不同信号点imf的加权平均if可用于定义瞬时互谱：式中是第p个imfs的加权瞬时频率：需要说明的是采用加权瞬时频率比直接采用频率的平均值更有意义。具体实施例本实施例中采用2002年6月4日在德州理工大学测量的多点雷暴风数据来证明本发明方法的有效性。选取高度6米、10米、15米的实测数据，持续时间为1800s，采样频率为1hz，如图2所示。首先采用memd将多点雷暴风速分解为7个imf(如图4，a为6米，b为10米，c为15米)和1个时变均值(图3)。从图中可以看出，这些时间均值变化趋势十分类似。此外，不同风速点的imf数量相同，同一尺度的imf具有相似的变化趋势，说明memd具有强大的过滤特性。然后通过改进的am/fm分解获得ia和if。原始风速的if估计平均值如表1所示。从中可以看出，在相同尺度下不同风速点平均if几乎相同，进一步反映了memd的优越性。表1.原始风速的if估计均值根据本发明提出的模拟方法步骤2能够得到瞬时频率。6米和15米高度的风速的瞬时频谱(如图5所示，a为6米，b为15米)和相应的瞬时互谱(如图6所示，a)。为了更好的演示，前两阶imf分解得到的if已被忽略，并通过引入具有20s固定宽度的移动平均窗口，可以看出瞬时频谱中同尺度的幅值和频率几乎一致。通过步骤3和步骤4之后，可以得到统一后的瞬时频率矩阵，其中6米和15米的瞬时互谱(如图6b)。图中瞬时互谱表现出于原始风速的相似的特征，这说明步骤3和步骤4处理方法的有效性。在进行步骤5和步骤6后，可以生成模拟风速。模拟风速和相应的时变平均值如图7和图3所示。与图2所示结果相比，可以发现模拟样本的变化趋势和时变均值与原始风速几乎一致。为了进一步证明所提方法的有效性，采用湍流特性、功率谱密度psd、累积正态化阿利亚斯烈度cnai，基于memd的瞬时频谱/瞬时互谱，k-l散度(kullback-leiblerdivergence)和协方差/互协方差来验证模拟结果。湍流特性湍流强度和阵风因子被广泛用于描述风场的湍流特性。雷暴风的湍流强度是随时间变化的，通常定义为：其中为memd求得的时变均值，σ(t)是通过窗宽为30秒的移动窗口获得的时变标准偏差。需要注意的是，通过该方法可以使获得的σ(t)更平滑。原始风速和模拟风速的时变湍流强度的对比情况如图8所示。可以看出，除了由极低的均值引起的峰值，模拟风速很好的保存了原始风速的湍流强度特性。雷暴风的阵风因子有多种定义，这可能导致不同的结果。本发明中将雷暴风的阵风因子定义为：式中：是t＝3s的移动平均风速的最大值，是时变平均风速的最大值。表2为原始风速和模拟风速的阵风因子的对比情况，可以看出它们总体上彼此接近，并且在6米高处差异的最大白封闭是8.1％。表2.阵风因子的对比高度(m)实测值模拟值差异占比61.231.140.098.1％101.181.140.043.1％151.211.160.053.6％功率谱密度除了湍流特性之外，反映频域总能量的功率谱密度也可以用来验证模拟精度。剔除时变均值后的原始风速和模拟风速之间的psd的比较如图9所示。可以看出，模拟风速和原始风速的psd基本保持一致。累积正态化阿利亚斯烈度相对于频域中psd显示的能量分布状况，相应的cnai用于描述的时域中的累积能量，能够用于验证所提出的模拟方法的准确性。原始风速的cnai定义为：其中0≤t≤td，td表示原始风速的持续时间，表示在时刻τ的原始风速。类似的方法也可以获取模拟风速的cnai。图10中比较了原始风速和模拟风速的cnai，可以看出它们几乎重合，验证了所提出的模拟方法的精确性。基于memd的瞬时频谱和瞬时互谱为了检验模拟风速在时域上的能量分布，采用了基于memd的瞬时频谱和瞬时互谱。6米和15米高度的模拟风速瞬时频谱(如图11所示，a为6米，b为15米)和瞬时互谱(如图12所示)。与图5和图6中原始风速的瞬时频谱相比，除了小的频移之外，它们具有相似的变化趋势。这是因为原始风速和模拟风速的分解得到的imf的数量不完全相同。k-l散度k-l散度(也称为散度，信息散度或相对熵)反映两种概率分布之间的差异程度。对于两个离散概率分布和相对于随机变量。它们的k-l散度可以表示为(kullback和leibler，1951)：式中，dkl(p||q)表示p(u)与q(u)之间的k-l散度，u是概率空间；从式(28)可以看出k-l散度值越小，两个分布差异越小。特别的，当k-l散度为0时，两个分布是相同的。k-l散度的另一个重要特性是不对称的。因此不能用于度量距离，则对称k-l散度的定义(johnson和sinanovic，2008)：式中，为对称k-l散度。此时，可以用于度量两个分布之间距离。本发明对称k-l散度用于描述两个随机过程之间的相关性。原始风速和2000个模拟风速样本的k-l散度的对比如图13所示。模拟风速样本的k-l散度几乎对称地围绕原始风速的波动，这说明了模拟风速很好的保存了原始风速的相关性。协方差和互协方差为了验证所提出方法的有效性，有必要将模拟风速的协方差/互协方差与目标进行比较。本发明中获取了2000个模拟风速的集合协方差/互协方差。图14中给出了在6米(a)和15(b)米的高度处模拟风速的协方差与目标值之间的比较。如15中给出了在6米和15米的高度处模拟风速的互协方差与其目标值之间的比较。可以看出，估计的协方差/互协方差与目标值几乎没有差异，证明了所提方法的有效性。本发获取基于memd的瞬时频谱矩阵，由于基于memd的瞬时频谱矩阵已经涉及空间相关性，可以在不假设随机初始相位角之差为常数的前提下，自然的考虑空间相关性；所提方法更加简单和合理；实施例证明了所提方法的有效性，可以看出，模拟样本的湍流特征，psd、cnai，基于memd的瞬时频谱或瞬时互谱和k-l散度与原始风速的湍流特征一致。估计的协方差和互协方差与目标值一致。所提出的方法在保持所测单个样本的空间相关性方面比假设随机初始相位角之差为常数的处理更方便和合理。具有非平稳特性的雷暴风造成了大约70～80％的输电线塔等结构的倒塌；获取大量雷暴风样本是进行雷暴风作用下结构非线性分析的前提。而实地测量只能获得单样本。进行输电线塔、风机等高耸结构的风振分析时，因结构和气动非线性，必须采用大量的风速样本对结构进行纽马克时域动力分析，计算多个响应时程，从而获得用于结构设计的响应方差和极值等统计参数。当前难以从单个非平稳雷暴风实测记录中获取时域动力分析中需要的大量风速样本。本发明结合多元经验模态分解和谱表示方法，提出一种基于单样本的时频模拟方法来可靠获取大量非平稳雷暴风速样本；采用基于memd的瞬时频率矩阵来考虑空间相关性，采用srm生成样本，进而无需更多假设，使得模拟方法更加简单合理。当前第1页12