02大气风特性与风环境022宋淳宸_W

1、首届中国空气动力学大会2018.08.绵阳基于 MEMD 时频分析的多点非平稳雷暴风速模拟宋淳宸 1，苏延文 2，黄3，刘瑞莉 4，杨成5（1,4,5. 西南交通大学土木工程学院，四川成都 610031；2. 中铁二院工程集团有限责任公司，四川 成都 610031;3. 重庆大学土木工程学院，重庆 400000）1引言强风诸如雷暴风、等产生的风荷载对建筑物和输电线塔等结构造成大量破坏， 近年来已引起广泛的关注。据统计，在美国、澳大利亚、南非等地与气候相关的输电线塔破 坏约 80%是由雷暴风造成的1。为了准确估算结构的动力响应，获得可靠的风速样本至关重要。目前，国内外大部分研究利用演化功率谱2-

2、6来模拟非平稳激励。实测风通常只有 单点记录或者少数几点记录。由于不确定准则而使得演化功率谱估计精度较低，因此难以准确模拟大量风速时程。 本文结合多元经验模态分解(multivariate empirical mode decomposition，MEMD)，调频函数/调幅函数（amplitude modulation/frequency modulation，AM/FM）分解和本征正交分解 （proper orthogonal decomposition，POD）建立了基于时频分析的非平稳多点风速的模拟方法。第一，采用 MEMD 分解多点风速，产生固有模态函数。第二，采用 AM/FM 分解计

3、算各点的固有模态函数的瞬时频率和瞬时幅值。第三，利用 POD 对瞬时频率解耦。最后，将瞬时幅值和解耦后的瞬时频率用于重构多点非平稳风速。实测多点雷暴风的模拟结果表明：相对于已有的方法，MEMD-AM/FM 分解-POD 方法能更好的模拟多点非平稳雷暴风速。 2MEMD 模型2.1 MEMD 模型MEMD7能够同时分解多点风速产生数目相同的 IMFs，因此是分析多点风速的有力工具。为了保持更好的一致性，MEMD 分解中采用低差异性 Hammersley 序列（类似蒙特卡洛 采样）。多点风速 x( t) = x (t), x (t),., x (t)T 通过 MEMD 模型分解得到 IMFs 的算

4、法步骤 12N如下： (1)在(N-1)-球上采样生成 P 个点的 Hammersley 点集； (2)计算信号 X( t ) 沿着所有方向向量vn (n= 1,2,.,P) 的映射向量，用 qn ( t ) 表示； (3)找到映射 q ( t ) P 的极大值、极小值所对应的时间量t (t ) ；P nnn=1in=1i(4)Pn=1对 t , x ( t) 采用三次样条插值求多元信号的包络线e (t )；nnnii(5)(6)求得包络均值 ( t ) ；c ( t ) = x ( t ) - ( t ) ，当c ( t ) 满足 MEMD 的 IMF 的终止条件时，重复步骤(2)-(5)获

5、基金项目：项目名称（基金编号） 首届中国空气动力学大会2018.08.绵阳得更高阶的 IMFs。MEMD 同时解决了 EMD/EEMD 模态混叠和无法应用于多元信号模拟的问题。而且相对于 SWT 等小波变换，MEMD 不需要先验基函数，具有更好的自适应性，因此能够提高多元信号的模拟精度。 2.2 AM/FM 分解尽管传统希尔伯黄变换可以得到 IA 和 IF，但为了获得具有物理含义的 IF，IMFs 不仅需要满足 EMD/MEMD 的终止条件，还必须同时满足著名的 Bedrosian 和 Nuttall8-9定理。因此采用希尔伯黄变换的改进方法即基于 AM/FM 分解来获得 IA 和 IF。AM

6、/FM 分解通过唯一分解 IMF 为 AM 函数 a ( t ) 和 FM 函数cos ( t ) 来实现，具体步骤如下： (1) 首先获取单元信号 z ( t ) 的绝对值的极值， z ( t ) 为 MEMD 分解产生的单个 IMF；(2) 采用三次样条拟合其绝对值的极值，得到经验包络e ( t ) 。用 z ( t ) 除以e ( t ) 得到 z ( t ) 。 111理论上e1 ( t ) = a ( t ) ，因此满足 z1 ( t ) = cos ( t ) 和 z1 ( t ) 1 。但实际上e1 ( t ) 经常不等于 a ( t ) ， 特别是 a ( t ) 变化很快时，

7、继续迭代至 zn ( t ) 1。此时 FM，AM 分别为： F ( t ) = cos ( t ) = zn ( t )A ( t ) = a ( t ) = z ( t )/ F ( t ) = e1 ( t ) e2 ( t )Len ( t )迭代后，IA 近似等于该 AM 函数。 (1)(3) 对 FM 求导，对导函数 F ( t ) 再次进行 AM/FM 分解。 F ( t ) 可表示为：F ( t ) = A ( t )cos ( t ) = ( t )cos + ( t ) (2)2其中 A ( t ) 和cos ( t ) 分别为导函数 F ( t ) 的 AM 和 FM 部

8、分，因此 IF 为： f ( t ) = A (t) / 2p = ( t ) / 2 (3)综上，通过两次 AM/FM 分解和一次求导得到 IA 和 IF，能够满足 Bedrosian 定理。由于没有涉及到希尔伯特黄变换，能够满足 Nuttall 定理。因此，该算法得到具有物理意义的IF 和 IA，并且没有负频率，能够显著提高风速模拟精度。 2.3 POD 分解POD10通过将多元信号的协方差或互谱对角化，将时间空间相关的多点随机过程解耦为一系列子过程。 Y( t ) = y1 ( t ), y2 ( t ),., yN ( t ) 是均值为零的离散矢量随机过程，现使得： TNY( t )

9、= h( t ) = h (t ) (4)n nn=1其中1, 2 ,K, n ， TI ， h( t ) TY( t ) 。 是 Y( t ) 的协方差矩阵C 的特征向 nY量，即为使得 Y( t ) 为最佳映射的模态函数。将 按特征值降序排列得到新的向量矩阵 = l ,l ,K,l ，仅用 N 阶模态函数重构时程函数： t N1 2tNtY( t ) =l h ( t )(5)n nn=1式中 Nt N ，因较大特征值对应的低阶模态为主模态，故舍去较小特征值对应的高阶模态。POD 通过截断多点风速高阶模态，使能量尽可能的集中在低阶主模态中，剔除噪声首届中国空气动力学大会2018.08.绵阳

10、的影响，从而还原数据特征并显著提高计算效率。3多点非平稳风速模拟3.1 算法步骤首先采用 MEMD 分解多点非平稳雷暴风速 X( t ) = x ( t ), x ( t ),., x ( t ) 产生 IMFs，再 T12N用 AM/FM 分解计算各点的固有模态函数的 IF 和 IA。为了平衡多点模拟方法的准确性和计算复杂性，参照 Wen11中的处理方式，假设实测风速分解得到的 IA 和 IF 能够用于模拟非 平稳过程。则第 k 个风速 xk ( t ) 展开后，并引入随机初始相位角j,k 可重构为： i ( j,k ( t ) d ( t )+j j,k )Mx (t ) = Rea (

11、t )e + AM,k ( t ) (6)kj,kj=1式中 Re 为实数部分，M 为分解的 IMFs 的层数，aj,k ( t ) 为第 j 阶的幅值，w j,k ( t ) 为第 k 个风 速第 j 阶的频率， j,k 为第 k 个风速第 j 阶服从（0,2）均匀分布的初始相位角， AM,k ( t ) 为第 k 个风速时变均值。 MEMD 类似带通滤波器，则可认为其分解得到的不同阶的 IMFs 之间相互独立，同阶的IMFs 之间具有相关性。因此考虑到多点非平稳风速同阶 IMFs 的相关性，采用 POD 对同阶IF 解耦，解耦后的 IF 矩阵为： Ntw ( t ) = h ( t )(7

12、)j,kj,r j,kr=1w j,k ( t ) 为 POD 解耦后第 k 个风速第 j 阶 IF 矩阵；j,r 为重列后第 j 阶第 r 个特征向量，h j,k ( t )为w j,k ( t ) 最佳映射的函数， Nt 为低阶主模态数目（ Nt N ）。解耦后第 k 个风速 xk ( t ) 重构 为： i ( j,k ( t ) d ( t )+j j,k )Mx ( t ) = Rea ( t )e + AM,k ( t ) (8)kj,kj=1上述模拟方法适用于多点非平稳风速模拟，模拟主要分为以下四个步骤：(1)采用 MEMD 分解多点风速X( t ) = x1( t ), x 2

13、( t ),., x (Nt )T 产生若干数目相同的 IMFs； (2) 采用 AM/FM 分解分别计算各点风速 IMFs 的 IA 和 IF；(3) 对 IF 矩阵进行 POD 分解产生新的 IF 矩阵； (4)采用算式(8)方法模拟多点非平稳风速。首届中国空气动力学大会2018.08.绵阳3.2 非平稳模拟算例本小节将采用基于 MEMD-AM/FM 分解-POD 的方法来模拟 RFD（rear-flank downdraft）。图 1(a)中的三组风速数据为塔四、五、六上高度为 10 m 的实测数据，持续时间为 1800s，采样频率为 1Hz。首先采用 MEMD 同时分解三组风速产生数目

14、相同的 IMFs，图 2 中为 MEMD 分解塔四实测风速得到的 IMFs。图 1(b)给出塔四、五、六经 MEMD 分解后风速的时变均值， 图1(c)中为塔四、五、六上高度为 10 m 风速时程的模拟结果。从图中可见模拟风速与实测风速基本一致，接下来将验证模拟风速与实测风速的统计特征是否保持一致。 404020200塔六0塔六20002000塔五塔五10001000塔四 0塔四 0时间/(s)时间/(s)(a)(b)(c)图 1 (a)塔四、五、六实测风速；(b)MEMD 分解时变均值；(c)基于 MEMD-AM/FM-POD 模拟风速图 2 MEMD 分解塔速产生的IMFs表 1 中列出了

15、实测风速与模拟风速的均值、方差和相关系数，其皆为 1800 个风速数据点的统计指标。r45、r46、r56 分别表示塔四与塔五、塔四与塔塔五与塔六之间的相关系数。表 2 为模拟风速与实测风速的均值、方差和相关系数的误差。塔四模拟风速的均值、方差与实测风速误差为 0.41%和 5.00%的，塔四与塔五相关系数的误差为 0%。可见模拟的风速与实测风速的统计特性基本保持一致。 3.3 与基于 SWT-HT 的风速模拟的对比为了进一步说明 MEMD-AM/FM 分解-POD 模拟非平稳风速的有效性，将与基于SWT-HT 的多点风速模拟进行精度对比。SWT12能够在保持信号的位移不变性的前提下通过平移和

16、伸缩母小波函数来获取不同尺度的子序列。首先通过 SWT 分解以上三组风速数据， 建立 HT 模型计算其 IA 和 IF，再将经 POD 解耦后 IF 采用算式(8)模拟。 表 1 中列出了两种模拟方法得到的均值、方差、相关系数，总的来说基于 MEMD-AM/FM首届中国空气动力学大会2018.08.绵阳分解-POD 的模拟法得到的统计值与原数据更加接近。表 2 中塔四基于 SWT-HT 模拟风速均值和方差的误差分别为 0.65%和 7.51%，塔四与塔五相关系数的误差为 21.88%。可见基于SWT-HT 的模拟法与原数据统计值的误差更大。可以看出 MEMD 与 AM/FM 分解相互结合，是分

17、析非平稳数据强有力的时频工具，提高了多点雷暴风速的模拟精度。 表 1 实测风速与模拟风速分析结果 实测MEMDSWT实测MEMDSWT实测相MEMDSWT 相指标塔均值均值均值方差方差方差关系数相关系数关系数r45r46r56塔四12.1712.0812.1212.0812.2512.1136.4631.6034.6431.8833.7227.850.640.610.640.620.780.70塔五塔六12.6312.7112.6525.9525.2724.900.810.730.80表 2实测风速与模拟风速分析结果误差MEMD 均值误差/%SWT 均值误差/%MEM 方差误差/%SWT 方差

18、误差/%MEMD 相关系数误差/%SWT 相关 系数误差/%塔指标r45r46r560.4100.630.650.240.155.000.892.627.5111.874.0401.649.8821.8814.751.23塔四塔五 本文结合 MEMD 和 AM/FM 两种分解方法的优势，并对 IF 进行了 POD 解耦，提出了 一种基于时频域分析的多点非平稳风速模拟方法，得到以下结论： （1）采用 MEMD 能够一次分解多点风速，有效抑制传统 EMD/EEMD 模型模态混叠并产生相同数目的 IMFs；进一步采用 AM/FM 分解得到 IA 和 IF，克服了传统 HT 产生大量没有物理含义负频率

19、的缺陷；引入的 POD 解耦 IF，进一步剔除噪声的影响，提高计算效率。 通过将三种方法的优势相结合，有效保证多点非平稳风速模拟精度。 （2）塔四基于 MEMD-AM/FM-POD 模拟风速的均值、方差与实测风速误差为 0.41%和5.00%，塔四与塔五相关系数的误差为 0%。可见本方法与原始风速的统计特性基本保持一致，模拟精度高。 （3） 分析对比实测风速数据与模拟风速数据的尺度图和互相关尺度图可知， MEMD-AM/FM-POD 模拟法能够保持其多点风速的能量与时变相关性。 （4） 塔四基于 SWT-HT 模拟风速均值和方差的误差分别为 0.65%和 7.51%，塔四与塔五相关系数的误差为

20、 21.88%。可见基于 MEMD-AM/FM-POD 模拟精度在均值、方差、相关性等方面优于基于 SWT-HT 风速模拟。 参考文献1瞿伟廉, 梁政平, 王力争,等. 下击暴流的特征及其对输电线塔风致倒塌的影响J. 地震工程与工程振动, 2010, 30(6):120-126.Huang G. Application of Proper Orthogonal Decomposition in Fast Fourier Transform Assisted Multivariate Nonstationary Process SimulationJ. Journal of Engineerin

21、g Mechanics, 2015, 141(7):04015015.Peng L, Huang G, Chen X, et al. Simulation of Multivariate Nonstationary Random Processes: HybridStochastic Wave and Proper Orthogonal Decomposition ApproachJ. Journal of Engineering Mechanics,23首届中国空气动力学大会2018.08.绵阳2017, 143(9):04017064.Zhao N, Huang G. Fast simul

22、ation of multivariate nonstationary process and its application to extreme windsJ. Journal of Wind Engineering & Industrial Aerodynamics, 2017, 170:118-127.45苏延文, 黄, 彭留留. 与反应谱相容的多点完全非平稳地震动随机过程的快速模拟J. 工程力学,2015(8):141-148.Gao Y, Yongxin W U, Dayong L I, et al. An improved approximation for the spectr

23、al representation method in the simulation of spatially varying ground motionsJ. Probabilistic Engineering Mechanics, 2012, 29:7-15.Rehman N, Mandic D P. Multivariate empirical mode decompositionJ. Proceedings of the Royal Society A, 2010b:1291-1302.Bedrosian, E. A product theorem for Hilbert transformJ. Proceedings o