基于EVSEWT的下击暴流时变平均风提取*

李春祥， 李 洲

(上海大学土木工程系 上海，200444)

引 言

极端风(热带气旋、下击暴流、龙卷风等)特性研究对于提高结构抗风设计水平、完善风荷载规范等意义重大。与大气边界层(atmospheric boundary layer，简称ABL)风不同，极端风具有明显的非平稳性，即风信号的统计特征值具有时变性。为了便于分析，学者们通常将非平稳风视为时变平均风与零均值随机脉动风的叠加[1]。在非平稳风特性分析中时变平均风提取的准确与否会影响时域中的相关分析或频域中的功率谱分析，从而产生较大误差，甚至使低频谱完全失去真实性。另外，在非平稳脉动风模拟与预测、结构非平稳振动响应分析和结构参数识别等研究中都涉及到信号时变平均项的提取[2]。因此，合理时变平均风速的提取对于非平稳风研究具有重要意义。

在时变均值提取的各类方法中，平均斜率法、差分法、最小二乘拟合等通常需要预先知道信号趋势项的类型，对处理具有复杂变化趋势的信号并不适用[3]。线性拟合、高阶滤波和滑动平均等缺乏统一的定量提取方法或准则，提取效果依赖于使用者的经验以及计算参数的选择，影响了模型的实用[4]。与此同时，由于经验模态分解(empirical mode decomposition，简称EMD)具有自适应性，不需要先验信息，使用方便；离散小波变换(discrete wavelet transform，简称DWT)在处理复杂信号上具有良好的时频分辨率，故这两种方法得到了学者们的普遍使用。陈隽等[4-5]提出使用EMD方法，并引入脉动风平稳度标准和截止频率，来提取台风风速记录中的时变平均风。申建红等[3]基于人工模拟风速数据，提出利用单尺度小波能量的突变确定时变均值分解的准确层次，较EMD分解方法更加准确有效。Su等[6]基于下击暴流实测数据，针对Kernel回归、DWT和集合经验模态分解(ensemble empirical mode decomposition，简称EEMD)方法给出了选择合适的方法和时间窗大小以获得合理的时变平均值的建议，推荐采用高阶的DWT方法用于下击暴流时变非平稳风提取。Wang等[7]将DWT方法用于时变均值提取，用非平稳风速模型对苏通大桥风场数据进行了分析，分析结果与早期结果吻合。但是EMD方法仍存在模态混叠、端部效应和缺乏数学理论等不足，DWT方法对信号分解层次和母小波类型的选择有较严格的条件，相关参数的确定需要依据经验。另一方面，目前的非平稳风研究虽然取得一定进展和突破，但仍存在以下一些挑战：由于下击暴流、龙卷风的发生具有偶然性，相关实测样本仍不丰富；对极端风机理揭示不足，数学分析理论有待创新；传统方法仍需要凭借丰富的经验，新一代精细化时频分析技术急需研究。

本研究的重点是提出一种基于能量的波谷寻找经验小波变换方法(energy valley searching empirical wavelet transform，简称EVSEWT)，并给出了时变平均风的确定准则。首先，简要介绍了时变平均风提取的常用方法，包括EMD，EEMD和DWT方法；其次，针对时变平均风的特征，对经验小波变换(empirical wavelet transform，简称EWT)进行改进，提出EVSEWT，着重于提高对低频信息的分辨率；最后，给出了时变平均风提取的结论和关于EVSEWT的改进和应用方面的展望。为避免以往研究只在信号层面进行效果验证，笔者将提取的时变平均风和脉动风用于结构非平稳响应分析，通过2组下击暴流风速试验，对比EMD，EEMD，DWT(db10，db18)和EVSEWT的效果。

1 数据与方法

1.1 数据说明

2002年6月4日和15日，在图1所示的沿直线均匀分布的7座观测塔上，同时记录两组不同高度的下击暴流风速时间序列。测量了超过40 m/s的阵风风速。一个记录的下击暴流被归类为后翼下沉气流(rear-flank downdrafts,简称RFD)，另一个归类为Derecho。该实测数据具有强烈非平稳特性，其中采样频率为1 Hz，总时长为30 min，具体实测信息可参考文献[8]。笔者选取塔4上高度15m处的RFD数据和Derecho数据，如图2所示。

图1 下击暴流测点布置示意图Fig.1 Diagram of measuring point arrangement of downburstL

图2 两组下击暴流实测数据Fig.2 Two sets of downburst data

1.2 常用方法

在时变均值提取的诸多方法中，本章主要介绍离散小波变换、经验模态分解和经验小波变换方法。

离散小波分析法是用小波基函数逼近信号的时频分析方法。小波基函数是一组由快速衰减、有限长的母小波经缩放和平移而生成的函数序列，常用的母小波函数有db族小波、Harr小波、sym小波、Morlet小波和Meyer小波等。但离散小波分析法存在分解层数和小波母函数选择的困难，对不同特征的数据，母函数的效果也有较大差距，常需要按经验选择。

EMD方法是一种信号分析方法，依据数据的时间尺度特征进行信号分解，无须预先设定基函数。它将一个信号分离为多个本征模态函数(intrinsic mode function,简称IMF)与一个余量R之和。本征模函数视为信号中的振动成分，反映了信号中不同频率的成分。从原始信号中先分离出来的本征模函数较后分离出来的频率高，余项的频率最低，为趋势成分，这样通过选择性地叠加低频分量就可得到原始信号的低频趋势项。EMD方法使用简单，适用于大多数非平稳信号分析，但EMD分解的结果会出现端部效应和模态混叠的问题，同时存在缺乏数学基础理论、分解的子信号物理意义不明确等不足。为了改进EMD的端部效应和模态混叠，Wu等[9]提出集成经验模态分解，通过在原信号中加入合适比例的白噪声，来补充一些缺失的尺度，在机械、土木、航空航天等领域的信号分解中具有良好表现。

在时变均值提取的应用中，Holmes等[10]通过两个准则确定移动平均法中的窗口大小，即时变平均风速需要反映原始风速低频部分的趋势、峰值以及脉动风的均值应趋近为0。Su等[6]根据如下3个准则确定合理的小波基函数和分解层数：①时变平均风速应反映原始风速的变化趋势；②对应脉动风的演化功率谱密度(evolutionary power spectral density,简称EPSD)估计结果应具有物理意义的解释；③EPSD作用下的高层建筑顺风向结构响应是否相对保守。为了简便，本研究中DWT方法按照尺度小波能量的突变确定时变均值的准确层次[3,11]；EMD和EEMD方法取分解结果中小于截止频率fc的子信号为时变均值项，通过间歇检测准则来完成，截止频率确定取决于高通信号的平稳性，参考文献[12]。

1.3 经验小波变换

为了解决小波变换参数选择困难、经验模态分解缺乏数学理论以及模态混叠等问题，Gilles[13]提出了一种经验小波变换。这种方法可以看作是构造1组带通滤波器和1个低通滤波器，用于将原信号分解成多个频率均一、稳定的子信号。具体实现流程如下：

1) 如图3所示，对原信号进行快速傅里叶变换，得到原信号标准化傅里叶谱(横轴标准化至0～π)；

图3 标准化傅里叶谱的划分Fig.3 Classification of standardized Fourier spectra

2) 在傅里叶谱上确定各带通滤波器的边界，并基于这些边界，利用Meyer小波构造滤波器组；

3) 原信号通过低通滤波器和这些带通滤波器，被分解为频率从高到低的若干子信号集。

1.4 EVSEWT介绍

尽管EWT已经给出了几种傅里叶边界划分的策略，如“locmaxmin”，“adaptive”和“scalespace”等，但由于应用的目的和场景不同，难以结合时变平均风的提取。不过，EWT方法有完整的数学理论，笔者基于EWT的理论框架，在傅里叶谱边界划分上更注重低频频段的提取，提出EVSEWT，图4所示为其提取时变趋势项的流程，具体如下：

图4 基于能量的波谷寻找经验小波变换原理图Fig.4 Schematic diagram of energy-based valley searching empirical wavelet transform

1) 对原始信号计算进行快速傅里叶变换，获得标准化傅里叶谱；

2) 在傅里叶谱的0～X之间，对每一处波谷设立划分边界，边界数量自适应确定；

3) 计算相邻两边界与傅里叶谱围成的面积，作为子信号能量指标，并按照能量从大到小对子信号进行排序；

4) 依据步骤2确定的边界，利用Meyer小波构造滤波器，按能量从大到小逐次提取子信号放入时变趋势项中，剩余部分子信号作为脉动分量；

5) 利用Runtest法确定脉动风的平稳性；

6) 重复步骤4和5，直到脉动风为平稳数据且均值趋近于零，若对0～X的信号处理完后，脉动风的平稳性要求未能满足，则只取0～X之间频段作为时变平均风。

在结构动力学分析中，一般认为平均风速最高频率为结构基频的1/5～1/10时，可以忽略其动力放大效应，按拟静力分析[6]。因下面分析的结构基频为0.23 Hz，故平均风速最高频率为0.023 Hz。对采样频率为1 Hz的信号来说，0.023 Hz在标准化傅里叶谱上对应为0.144 Hz。笔者确定X为0.144 Hz，时变趋势项的频率控制在0～X范围中，且X可以视使用情况调整。总体来讲，EVSEWT注重低频部分的分解，即0～1/10结构基频处，这样提取的时变均值项满足按拟静力分析的要求，也能较好地刻画实测数据的趋势。

2 非平稳随机响应分析

2.1 理论介绍

将实测风速分为时变平均风和脉动风

(1)

(2)

其中：M1为结构第1阶振型的广义质量；ω1为结构第1阶振型的角频率。

采用精确且高效的虚拟激励法对结构的随机动力响应进行分析，结构的虚拟激励力可以表示为

(5)

(6)

虚拟位移y(ω,t)可以由Newmark-β逐步精细积分法计算，参见文献[15]，则广义位移脉动分量的EPSD矩阵Sq′(ω,t)=y(ω,t)y(ω,t)T。当非平稳风荷载的EPSD时间变化率大于结构自然周期时，则可以采用拟平稳假设计算Sq′(ω,t)，结果更偏保守，计算快速

Sq′(ω,t)=|H(ω)|2SQ′(ω,t)

(7)

其中：H(ω)为基本振型下的频率响应函数[16]，H(ω)=1/[M1(-ω2+2iξ1ω1ω+ω2)]。

(10)

其中：ρ为空气密度，取1.225 kg/m3；CD为阻力系数，取1.2；B为建筑物长度；χD(ω)为气动导纳函数；g(z,ω,t)为调制函数；dΘ(z,ω)表示复值零均值正交增量随机过程。

笔者选取Davenport气动导纳函数

|χD(ω)|2=2/λy(1-1/λy+1/λye-λy)

其中：λy=kyfD/Umax；ky=8为衰减因子；f为频率；Umax为时变平均风速的最大值。

当只考虑结构的第1阶振型时，广义力的时变均值和脉动分量计算如下

(11)

(12)

其中：φ(z)=(z/H)β为第1阶振型的振型函数；β为振型指数，当振型为线性时取为1。

将式(9)代入式(11)可得

(13)

广义力的脉动分量Q′(t)的EPSD可表示为

当脉动风的EPSD假定为沿着建筑高度不变时，即Su′(z1,ω,t)=Su′(ω,t)，式(14)可表示为

(15)

若时变平均风剖面进一步简化为

(16)

此时，时变平均风速引起的广义力可表示为

(17)

相应的，式(15)可以简化为

(18)

其中：|Jz(ω)|2为联合接收函数。

(19)

由于实测样本少以及在数学上处理的困难，在大多数研究中会采用时不变指数函数或Davenport提出的相干函数。笔者假设z1和z2高度处的相干函数是时不变的

(20)

其中：kz为衰减因子，取值为8。

本研究中，下击暴流风速的竖向风剖面选取Wood[17]提出的时不变经验模型

(21)

以上数值风剖面经验公式已经通过实测风速数据验证[18]，在其他研究中也被广泛采用[6,16]。由于目前尚缺乏深入的下击暴流风速知识以及足够数目的观测样本，虽然上述假定在研究中会产生不确定性，但是这些简化处理不会影响一般性的规律与结论。

除此之外，采用基于穿越率的极值分布理论可以得到响应的平均极值和标准差[19-20]。设R(t)为时间周期T内的非平稳响应，根据穿越率的泊松近似关系，在时刻t和水平r的平均穿越率可由下式计算

(22)

由文献[21]可知，非平稳响应R(t)在某时间周期T内的极值的累积分布函数如下

(23)

(24)

2.2 结果分析

为了验证EVSEWT在时变均值项提取中的有效性，将EVSEWT用于RFD和Derecho两组实测下击暴流风速的时变均值提取试验，并选取EMD，EEMD和DWT做对比分析。根据文献[6]的建议，母小波函数选择db10和db18，最后将提取的结果用于某幢高层钢结构建筑的随机动力响应分析。

已有的研究表明[20-21]，下击暴流风速会随地面高度增加而增加，在50～200 m高度处达到峰值，再往上风速下降缓慢。因此，下击暴流可能对悬臂结构产生显著的风荷载影响[18]。下面以一座高H=200 m，长B和宽D均为40 m的结构为例。结构密度设为192 kg/m3，结构的基频f1=46/H=0.23 Hz，振型阻尼比ξ1=1%。对于高层建筑，风振响应一般以基本振型为主，高层振型的贡献可以忽略，笔者只考虑第1振型的影响，且振型形状假定为线性函数。具体地，以结构顶端均值位移、基于拟平稳和非平稳的RMS位移等结构响应值以及结构非平稳响应的平均极值作为评价指标。

简洁起见，只给出了RFD的情况作图示说明。图5是由EMD，EEMD，DWT(db10)，DWT(db18)和EVSEWT方法，依据各自的准则提取的时变平均风速和脉动风速。图6是脉动风速的EPSD图，EPSD图刻画了信号在不同时刻、不同频率上的能量分布，一般认为大多数嵌入在脉动风中的能量分布在0～0.1 Hz的低频范围内。同时，在一定时间内，随着频率的增加，大多数EPSD都呈现出衰减的趋势。脉动风EPSD也是后续结构响应分析中需要用到的物理量之一。从图2中RFD的实测数据可以看出，信号的能量主要集中在150～700s之间，这与图6结果相吻合，且能量集中在0～0.1 Hz频段，与经验相符。另外，从图6的结果来看，EMD和EEMD的低频部分能量最高，可能导致更大的结构动力响应，EVSEWT次之，DWT的最小。对应图5中提取的时变平均风，DWT方法的波动比较明显，EMD和EEMD的波动平缓，EVSEWT介于两者之间。

图6 不同分解方法对应脉动风的演化功率谱密度Fig.6 Evolutionary power spectral density of fluctuating wind speeds extracted by different methods

图5 不同分解方法获得的平均风速和脉动风速Fig.5 Time-varying mean and fluctuating wind speeds obtained by different decomposition methods

与ABL风相似，非平稳风激发的结构响应一般分为静风响应和动力响应两部分。前者由时变平均风速引起，通过拟静力方法简单计算；后者由平稳/非平稳脉动风引起，通过动态分析进行计算[22]。图7是不同信号分解下，RFD风速在结构顶端产生的均值位移(Mean)、基于拟平稳(Pseudo-stationary)和非平稳的(RMS)位移。可以发现，均值位移与提取的时变平均风的趋势类似，DWT方法的波动较大，EMD和EEMD方法的均值位移相对平缓。需要指出，基于非平稳假设计算的RMS位移的最大值略小于基于拟平稳假设的方法，且在时间上，非平稳响应最大值要滞后于拟平稳响应和均值位移的最大值，这与文献[6,15]的结论一致。

表1和表2分别是RFD和Derecho两组下击暴流风速在不同方法下得到的结构200 m高度处的均值位移最大值、RMS位移最大值和平均极值位移。综合图7、表1和表2，可以得出以下结论：①在RFD风速中，最大相差发生在RMS位移，为11.24%(EMD法和EVSEWT法之间)，均值位移相差最大达9%(EMD法和EVSEWT法之间)，相对来看，Derecho风速中各方法相差较小，为RMS位移的7.81%(EMD法和EVSEWT法之间)，说明EVSEWT方法对结构响应分析相对保守；②从3项位移指标来看，2组试验均表现出EVSEWT方法最佳、DWT方法次之、EEMD和EMD略差的结果，这可能是由于EVSEWT方法是基于信号傅里叶谱进行的分解，在信号低频部分具有更好的分辨率，配合相应的准则能比较合理地提取时变平均风速。

图7 顶端位移均值和RMS位移Fig.7 Mean and RMS responses of tip displacement

表1 试验1(RFD)结构顶端位移的最大均值、最大RMS和平均极值位移

表2 试验2(Derecho)结构顶端位移的最大均值、最大RMS和平均极值位移

3 结 论

1) EWT提取的时变均值能兼顾实测信号趋势和高频波动，对低频信息的分辨率更优秀。

2) 在时变平均风极大值和脉动风的能量都不突出的情况下，能计算得到偏于安全的结构响应位移，说明EVSEWT的确定准则更加合理。

3) 3个响应位移之间相互关系符合经验[19,23]，即非平稳RMS位移极大值略小于拟平稳RMS位移极大值，且滞后于拟平稳RMS位移极大值和均值位移最大值。

4) 总体来说，EVSEWT具有非常好的信号低频分辨率，对2组试验数据的时变均值提取试验效果较好。后续的研究可以考虑用于桥梁和高层建筑结构健康监测信号时频分析[24]、机械故障诊断和生物医学信号处理等领域。另一方面，受到多元经验模态分解的启发，将EVSEWT拓展到空间信号处理问题，对于注重空间相关性的风场分析也很有意义。

致谢本研究风速数据来源于德州理工大学风工程研究中心，在此表示衷心感谢。