改进VMD及其在时变结构地震响应瞬时频率识别中的应用.pdf

第 37 卷第 3 期2024 年 3 月振 动 工 程 学 报Journal of Vibration EngineeringVol.37 No.3Mar.2024改进 VMD及其在时变结构地震响应瞬时频率识别中的应用王焱楠，尚旭强，黄天立（中南大学土木工程学院，湖南 长沙 410075）摘要:强地震动作用下，结构动力特性因损伤将随时间变化，变分模态分解（Variational Mode Decomposition，VMD）可用于分析结构地震响应的瞬时频率变化规律，揭示地震过程中结构的损伤状态。针对人为预设分解模态数 K 和二次惩罚因子参数不准确导致 VMD 非平稳响应出现模态混叠的问题，提出了一种改进的 VMD（Improved VMD，IVMD）算法，结合 Hilbert变换（HT）可准确识别非平稳地震激励下时变结构的瞬时频率。采用多重信号分类（Multiple Signal Classification，MUSIC）算法确定分解模态数 K，基于整体正交系数和能量比系数构建综合目标函数，并采用樽海鞘算法（Slap Swarm Algorithm，SSA）优化确定二次惩罚因子。在优化确定 K 和参数的基础上，结合 IVMD 和 HT 识别时变结构地震响应瞬时频率。典型模拟信号和地震激励下 4层时变框架结构数值模拟响应表明，相比 VMD 算法，IVMD 算法识别得到的瞬时频率精度更高；12层钢筋混凝土框架结构模型地震模拟振动台试验数据验证了所提方法的实用性。关键词:地震激励；时变框架结构；振动台试验；瞬时频率识别；变分模态分解中图分类号:TU311.3；TU323.5 文献标志码:A 文章编号:1004-4523（2024）03-0402-12 DOI：10.16385/ki.issn.1004-4523.2024.03.005引 言地震是一种难以预测且破坏性极强的自然灾害，常造成土木工程结构不同程度的破坏。结构一旦发生破坏，其自振频率、阻尼比等动力特性随之变化。利用先进的非平稳信号分析方法，分析非平稳地震激励下结构的地震动响应，准确提取结构的瞬时频率等时变动力特性，对于结构的损伤诊断和安全状态评估具有重要的价值1。Loh等2采用子空间识别算法，分析了一栋 7层钢筋混凝土框架结构、一栋 9 层层间隔震建筑和一座隔震桥梁的输入输出地震响应数据，识别了结构的自振频率和阻尼比等动力特性，并将其应用于结构地震安全评估。Ni等3提出了一种基于结构地震输入输出响应数据的快速贝叶斯频域模态参数识别方法，基于某 12层框架结构模型振动台模拟地震试验数据和 7 层 Van Nuys 酒店建筑地震响应数据识别得到的结构模态参数，验证了方法的有效性和适用性。Ghahari等4指出，地震动输入数据常常难以获得，且考虑土结构相互作用时，记录的自由场地震动响应或基础响应数据都不能作为输入用于基于输入输出数据的模态参数识别方法，因此，盲源分离5、递归随机子空间识别6和频域分解7等仅基于输出响应的模态参数识别方法更适用于识别地震激励下时变结构的模态参数。时频分析方法作为处理非平稳信号的有力工具，通过分析结构地震响应数据的时频谱图，提取结构瞬时频率等特征参数，可用于揭示结构在地震过程中的损伤状态8。樊海涛等9采用经验模态分解（EMD）从结构地震响应信号中提取结构模态响应，利用 Hilbert变换识别了结构瞬时频率等时变特征，掌握了地震过程中结构损伤的发展规律。石春香等10采用 EMD+HT 方法，分析了深圳罗湖商务大厦模型模拟地震振动台响应数据，利用 Hilbert阻尼谱研究了结构在地震过程中的破坏规律。陈隽等11采用 EMD+HT 方法识别了线性时变系统的瞬时频率和瞬时阻尼比特征，并将其应用于地震激励下12层短肢剪力墙结构振动台模拟地震试验数据，识别了模型结构在试验过程中动力性能的时变特性。应该指出，由于 EMD 为经验算法，没有严格的数学理论基础，分解过程中常存在模态混叠等问收稿日期:2022-10-31；修订日期:2022-12-09基金项目:国家自然科学基金资助项目（52078486，U1734208）；中南大学中央高校基本科研业务费专项资金资助项目（2023ZZTS0410）。第 3 期王焱楠，等：改进 VMD及其在时变结构地震响应瞬时频率识别中的应用题12。为更好地识别多分量非平稳信号的时变频率，提出了一种变分模态分解（VMD）方法13，该方法将信号同时分解为不同中心频率的窄带信号，具有较高的时频分辨率和噪声鲁棒性，已广泛应用于机械故障诊断14、模态参数识别1516和风速预测17等非平稳信号分析领域，但 VMD 应用于非平稳地震激励下结构响应分析和损伤诊断方面的研究工作还较少。VMD 需人为预设分解模态数 K 和二次惩罚因子等参数，常采用多次试算方法确定参数取值，由此导致 VMD 在分解非平稳响应时常出现模态混叠问题。为此，许多学者提出了一些改进措施，以实现VMD 参数优化选取。Huang等18和 Ma等19应用尺度空间表示和尺度空间谱分割来估计分解模态数K。Lian等20采用能量损失、置换熵等指标反复搜索分解模态数 K。唐贵基等21以包络熵为优化目标，采用粒子群算法同时搜寻 K 和的全局最优解。何勇等22以峭度和包络熵作为遗传算法的目标函数，对 K和同步联合优化，并应用于轴承故障诊断。应该指出，同时对 K 和两个参数进行优化，计算效率低，且参数间相互影响导致分解模态数 K 出现偏差，即使参数取得最优值，也可能造成分解存在分量丢失或分量混叠等问题。因此，本文提出分别优化参数 K和的思路，保证二者取值均为最优。此外，参数优化取值时常采用单一的局部目标函数23，这对于机械故障诊断是适用的，但不适用于非平稳地震激励下时变结构瞬时频率的识别问题。针对上述 VMD 算法存在的不足，本文提出一种改进 VMD（IVMD）算法，并应用于非平稳地震激励下时变结构响应瞬时频率识别。首先，考虑到分解模态数 K 对 VMD 结果的决定性作用，基于 MUSIC 功率谱确定K；然后，考虑到仅对二次惩罚因子进行优化可提高计算效率，采用整体正交系数和能量比系数构造综合目标函数，并基于樽海鞘算法（SSA）优化确定；最后，基于优化确定的参数 K 和，采用 IVMD 算法分解地震响应得到固有模式函数（IMF），并结合 Hilbert 变换识别结构瞬时频率。采用模拟信号、地震激励下 4 层时变框架结构的数值模拟响应和 12 层框架结构模型振动台模拟地震试验数据验证了所提方法的有效性和实用性。1VMD算法VMD 算法采用一种非递归的变分模型，预设分解模态数 K，将多分量信号s(t)同时分解为K个中心频率为k的单分量调幅调频信号uk(t)，即固有模式函数（IMF）分量。基于每个 IMF 分量围绕其中心频率有最紧密带宽的基本假定，构造uk(t)的带宽方程：min()uk，k k=1K t ()()t+jtuk()te-jkt22 ，s.t.k=1Kuk()t=s(t)（1）式中 j为虚数单位；表示卷积计算；uk(t)为分解得到的第 k个 IMF 分量；k为uk(t)的中心频率；t()为对时间 t求偏导数；(t)为脉冲函数。将uk的解析信号与e-jkt相乘并将其频谱移至基带，计算解调信号梯度L2范数的平方，构成约束变分模型。引入拉格朗日乘子(t)对约束进行控制，二次惩罚因子提高变分问题的收敛性，从而转化为无约束变分问题：L(uk，k，)=k=1K t ()()t+jtuk()te-jkt22+s()t-k=1Kuk()t22+()t，s()t-k=1Kuk()t（2）式中 L表示增广拉格朗日量；表示求内积。通 过 乘 法 器 交 替 方 向 法（ADMM）交 替 求 解un+1k，n+1k和n+1，寻找增广 Lagrange 函数的鞍点，得到最优解：un+1k()=s()-i kuni()+n()21+2()-nk2（3）n+1k=0|un+1k2d0|un+1k2d（4）式中 s()，ui()和()分别为s(t)，ui(t)和(t)的傅里叶变换；n为迭代次数。方程（4）计算模态分量uk()的过程，相当于对残 余 信 号 分 量s()-ui()进 行 维 纳 滤 波。VMD 在 频 域 对 每 个 IMF 分 量 进 行 更 新，对Re uk(t)进行傅里叶逆变换得到 IMF分量的时域表示。对于所有 0，拉格朗日乘子n+1()由下式更新：n+1()=n()+s()-kKun+1k()（5）式中 为递归步长，当无需精确重构原始信号时，=0。不断按式（3），（4）和（5）迭代更新，直到满足以下收敛准则：403振 动 工 程 学 报第 37 卷kKun+1k-unk22/un+1k22（6）式中 为设定的判别精度。2基于改进 VMD 算法识别时变结构瞬时频率VMD 算法分解信号时，采用的分解模态数 K和二次惩罚因子对分解结果影响较大，其中 K 直接决定分解结果的准确性，K 取值过小或过大将分别导致模态混叠或冗余噪声模态；二次惩罚因子控 制 每 个 IMF 的 频 带 宽 度 和 噪 声 含 量。目 前，VMD 算法分解信号时采用的参数 K 和具有随机性，且缺乏选择标准，从而影响了 IMF 分量用于识别瞬时频率的准确性。因此，确定优化的参数 K 和对基于 VMD 识别时变结构地震响应的瞬时频率至关重要。2.1确定及验证分解模态数 K傅里叶频谱常用来确定 VMD 的分解模态数K15，然而，傅里叶变换处理非平稳信号时，其频谱中存在大量虚假峰值，此外噪声干扰引起频谱中出现较多波动24，因此，傅里叶频谱常无法准确确定分解模态数 K。相比傅里叶频谱，MUSIC 功率谱更平滑，分辨率更高25，能很好地处理非平稳信号并减小噪声干扰，基于 MUSIC 功率谱能更准确地确定分解信号的分解模态数 K。MUSIC 功率谱的基本思想是对信号的(p+1)阶相关矩阵作特征分解，并把特征向量张成信号空间V1，Vp+1和噪声空间VM+1，Vp+1，然后根据信号和噪声空间的正交性进行功率谱估计：PMUSIC()=1AH()(l=M+1p+1VlVHl)A()（7）式中 A()=1 ej ejM；VHl为Vl的共轭转置向量；M为信号个数。PMUSIC()在=i处呈现一个尖峰，因此，可根据 MUSIC 功率谱的峰值数确定分解模态数 K。为验证 MUSIC 功率谱确定分解模态数 K 的准确性，构造一个具有时变频率的非平稳信号x(t)，该信号由 3个频率按正弦变化的单分量信号组成：x(t)=x1(t)+x2(t)+x3(t)+(t)x1(t)=2sin(0.4cos(2t)+4t)x2(t)=5sin(0.3cos(t)+8t)x3(t)=sin(0.2cos(2t)+14t)（8）式中(t)为白噪声，其功率为信号功率的 10%，信号x(t)的信噪比为 14 dB，信号时长为 10 s，采样频率为200 Hz。图 1所示为信号x(t)的 MUSIC功率谱、傅里叶频谱和短时傅里叶谱。由图 1（a）可见，MUSIC 功率谱有 3 个明显的峰值，不存在虚假峰值且几乎无波动，由此确定分解模态数 K=3，真实反映了信号x(t)的分量组成。由图 1（b）可见，傅里叶频谱中有7 个明显的峰值（其中峰值，和为虚假峰值）且波动较大，由此确定分解模态数 K=7，这与信号x(t)的分量组成不符。由图 1（c）可见，短时傅里叶谱模糊不清，无法得出分解模态数 K。由此可知，图 1 信号x(t)的 MUSIC 功率谱、傅里叶幅值谱、短时傅里叶变换谱Fig.1 MUSIC power spectrum，Fourier amplitude spectrum and short-time Fourier transform spectrum of x(t)404第 3 期王焱楠，等：改进 VMD及其在时变结构地震响应瞬时频率识别中的应用MUSIC 功率谱有抑制噪声、分解非平稳信号的能力，可用于准确确定非平稳信号分解模态数 K。2.2优化确定二次惩罚因子确定 K 值后，VMD 算法分解结果的准确性主要取决于二次惩罚因子。取值偏小，则各阶 IMF受噪声干扰大；取值偏大，则信号过度分解。为得到仅包含单一频率成分且受噪声干扰小的 IMF，本节基于 VMD 算法分解的后验信息构建目标函数，采用优化算法实现的确定。2.2.1构造目标函数包络熵、峭度、相关性等目标函数常用于机械故障诊断领域，然而，这些目标函数仅针对单一 IMF。基于分解得到的全部 IMF，构造综合目标函数评估所有 IMF 的性能，可提高识别时变结构各阶瞬时频率的准确性。（1）整体正交系数 P整体正交系数 P 表示 IMF 之间的正交性，其表达式如下：P=j=1Ki=1i jKm=1Nfj(m)fi(m)m=1Nf2(m)（9）式中 fi(m)，fj(m)和f(m)分别为第 i，j 个 IMF 和信号x(t)的频域序列；m 表示离散变量；N为信号长度。信号分解完全时，各阶 IMF分量正交，整体正交系数P值最小。应该指出，由于实际信号中不可避免地存在噪声，极易出现整体正交系数 P最小时，相应的取值过小，由此导致各阶频率的能量聚集性差。（2）能量比系数 E为解决单一整体正交系数P易受噪声干扰的问题，引入能量比系数 E，其定义为各阶 IMF能量与原始信号能量之比，表达式如下：E=k=1Ki=1Nu2k(i)Ni=1Nx2(i)N=k=1Km=1Nf2k(m)2Nm=1Nf2(m)2N（10）根据帕塞瓦等式，式（10）有频域和时域两种表达式。当信号分解不完全时，能量比系数 E值偏小。反之，信号过度分解时，能量比系数 E值偏大。（3）综合目标函数 F当整体正交系数与能量比系数的取值合适时，既可保证模态正交性，也可避免能量损失。因此构造综合目标函数 F，表达式如下：F=Pnorm+(1E)norm=Pnorm+Rnorm（11）式中 Pnorm和Rnorm分别为归一化整体正交系数和能量比系数的倒数。F表示各 IMF之间正交系数和分解过程中能量损失之和，当 F取最小值时，减轻了噪声影响且各阶 IMF 分量正交，从而避免了模态混叠和过度分解。通过寻找综合目标函数最小值 Fmin，即可确定，Fmin对应的即为最优值，表达式如下：Fmin=arg min(Pnorm+Rnorm)=arg min(F)（12）2.2.2验证目标函数为验证构造的综合目标函数 F 的有效性，图 2给出了信号x(t)的整体正交系数 P、能量比系数倒数 R 和综合目标函数 F 与取值的关系曲线，取值范围为 500，2.5105。由图 2 可知，整体正交系数 P 在 B 点取最小值，但此时取值偏小，VMD 算法分解得到的模态分量受噪声干扰大，单一 P 值无法确定的最优值。能量比系数倒数 R 在 C 点取最小值，但此时取值偏大，模态混叠严重，单一 R 也无法确定的最优值。本文构造的综合目标函数F，在介于 B 和 C 两点之间的 A 点取最小值，相应的=34000为最优值，此时，VMD算法分解结果噪声干扰小且模态混叠程度轻。由图 2 可知，当综合目标函数 F 在突变区间取最大值时，对应的值用于 VMD 时会造成过度分解，由此表明的取值偏大。尽管图 2中 D点处综合目标函数 F 取值较小且变化平稳，但对应的值过大导致 VMD严重过度分解，从此处选取的值不合适。因此，本文进一步将的搜寻取值范围限定在 500，40000 范围内，以保证取值为局部最优而非全局最优，并提高了搜索效率。2.2.3SSA优化搜索二次惩罚因子樽海鞘算法（SSA）26是一种群体智能优化算法，具有搜索速度快、有效避免前期搜索不充分而陷入局部最优的特点。该算法将种群个体分为领导者图 2 目标函数与取值的关系曲线Fig.2 The relationship curves between objective functions and value405振 动 工 程 学 报第 37 卷和追随者，前面的个体领导后面的个体相继排列成“链”状，迭代更新个体位置以搜寻最优解。本文采用 SSA算法优化搜索二次惩罚因子的最优值。针对信号x(t)，2.1节中已确定分解模态数 K=3，本节采用综合目标函数 F 和 SSA 算法优化搜索二 次 惩 罚 因 子的 最 优 值，的 搜 寻 范 围 为 500，40000。图 3 所示为综合目标函数 F 的迭代曲线。由图 3 可知，SSA 算法迭代 18 次后收敛，此时 Fmin为 0.9952，对应的最优值为 34000，与网格优化结果一致。SSA 算法迭代收敛用时 112.6 s，网格优化（搜寻步长为 1）算法迭代收敛用时 3424.9 s，SSA算法节约计算时间约 96%，搜寻效率高。综上所述，采用 MUSIC 算法确定分解模态数K，构造综合目标函数 F，采用 SSA 算法优化搜索二次惩罚因子，实现了 VMD 参数的优化确定，本文称其为改进变分模态分解（IVMD）算法。2.3Hilbert变换识别时变结构瞬时频率基于 IVMD 算法分解结构响应得到 IMF，基于Hilbert变换即可计算得到结构的各阶瞬时频率，计算公式如下：ki(t)=uk(t)uk(t)t-uk(t)uk(t)tu2k(t)+u2k(t)（13）式中 ki(t)为第k个模态分量的瞬时频率；uk(t)为uk(t)的 Hilbert变换。针对信号x(t)，基于 2.1 节和 2.2 节中确定的VMD 算法分解最优参数，即 K=3，=34000，利用IVMD对其分解，并将分解后的IMF进行HT识别得到各阶瞬时频率，如图 4 所示。由图 4 可知，识别的信号x(t)前2阶瞬时频率与理论值完全吻合；由于噪声干扰，识别的第 3阶瞬时频率存在一定的波动，但其均值与理论值基本一致。此外，由于 HT的端点效应，信号端部识别得到的各阶瞬时频率误差较大。2.4IVMDHT算法识别时变结构瞬时频率流程图 5 给出了基于 IVMDHT 算法识别非平稳信号瞬时频率的流程图，其基本步骤包括：（a）根据图 3 综合目标函数 F 迭代收敛曲线（地震激励下时变框架结构）Fig.3 The iterative convergence curve of composite objective function F（a time varying frame structure under seismic excitation）图 5 基于 IVMD-HT算法识别瞬时频率流程图Fig.5 Flowchart of instantaneous frequency identification based on IVMD-HT algorithm图 4 IVMD 算法结合 HT 识别的瞬时频率与理论瞬时频率（K=3，=34000）Fig.4 Instantaneous frequencies identified by IVMD algorithm combined with HT and theoretical instantaneous frequencies（K=3，=34000）406第 3 期王焱楠，等：改进 VMD及其在时变结构地震响应瞬时频率识别中的应用式（7）的 MUSIC 功率谱估计，得到光滑、分辨率高的 MUSIC 功率谱，根据谱峰确定信号的分解模态数 K；（b）设置二次惩罚因子的搜寻范围，基于综合目标函数 F并采用 SSA 算法搜寻的最优值，Fmin对应的值即为最优；（c）利用最优 K 和，对信号进行 IVMD 得到一系列 IMF，并采用 HT 识别信号各阶 IMF 的瞬时频率。其中，步骤（a）和（b）两部分即为 IVMD算法。2.5分析与讨论为验证本文算法获取优化的参数 K 和后识别瞬时频率的准确性，人为选取了不同的参数 K 和，图 6 给出了 VMD 算法识别的信号x(t)的瞬时频率。由图 6 可知，当 K 随机取值，=34000 时，尽管取最优值，但 K 取值偏小或偏大分别出现缺少高阶频率或存在冗余噪声频率的问题，K 值不准确直接导致结果偏差。当 K=3，随机取值，偏小时，各阶瞬时频率的识别结果受噪声干扰大；偏大时，尽管识别了第 2 阶频率，但由于模态混叠，第 3 阶和第 2阶频率完全重叠，第 3阶频率无法被识别。由此表明，人为随机对 K 和取值对瞬时频率的识别结果影响很大，会导致不能获取清晰准确的频率信息；采用优化取值的参数 K 和，可更准确地识别信号的瞬时频率。为定量评价信号瞬时频率的识别精度，定义均方根误差 RMSE：RMSE=i=1N(fq-fq)2N（14）式中 fq为瞬时频率识别值；fq为瞬时频率理论值；N为响应信号的长度。当 K 取值准确（即 K=3）时，表 1 对比了不同取值下 VMD 识别得到的信号x(t)各阶瞬时频率的RMSE 值。由表 1 可知，当采用优化确定的参数（即=34000）时，识 别 得 到 的 各 阶 瞬 时 频 率 的RMSE 值最小，结果最准确。由此表明，采用优化取值参数的 IVMD 算法，可解决 VMD 算法参数随机取值导致瞬时频率识别不准、精度不高的问题。3数值模拟为验证 IVMDHT 算法识别地震激励下时变结图 6 参数 K和取不同值时，VMD算法识别的瞬时频率与理论瞬时频率Fig.6 Instantaneous identified frequencies by VMD algorithm with different K，and theoretical instantaneous frequencies407振 动 工 程 学 报第 37 卷构瞬时频率的有效性和准确性，构建了一个具有时变刚度的 4层剪切型框架结构，如图 7所示。结构参数取值为：各层质量m1=m2=m3=m4=4 kg；各层阻尼系数c1=c2=c3=c4=0.5 N s/m；第 3，4层刚 度 系 数 为 常 值，k3=1.8 104 N/m，k4=1.5 104 N/m，第 1，2 层刚度系数随时间变化，k1(t)=k2(t)=1-0.1sin(2t)104 N/m。地震激励 Q 采用图 8 所示的 El Centro 地震波，激励时长为 100 s，采样频率为 250 Hz。采用四阶龙格库塔（RungeKutta）进行结构地震响应分析，并添加 10%的白噪声。图 9所示为 El Centro地震激励下时变框架结构顶层的加速度响应及其傅里叶频谱和 MUSIC 功率谱。响应时长为 100 s，其中 520 s信号能量最大。图 9（b）所示的傅里叶频谱中存在虚假峰值和波动干扰，峰值多而杂乱，难以准确确定分解层数 K；图 9（c）所示的 MUSIC 功率谱中可清楚地观察到 4个峰值，由此确定分解模态数 K=4。确定分解模态数 K 后，设置 SSA 算法的种群个数 为 10，最 大 迭 代 次 数 为 100，搜 索 范 围 为表 1 不同取值下瞬时频率识别值的 RMSE（K=3）Tab.1 RMSE of identified instantaneous frequencies with different value（K=3）1阶2阶3阶总和RMSE=50000.350.851.212.41=340000.300.570.921.79=345000.330.47注：“”表示该方法未能识别相应瞬时频率。图 9 地震激励下结构顶层加速度响应及其傅里叶频谱和 MUSIC功率谱Fig.9 Acceleration timehistory，Fourier amplitude spectrum and MUSIC power spectrum of top layer of structure under seismic excitation图 7 4层剪切型时变刚度框架结构Fig.7 A four-story shear frame structure with timevarying stiffness图 8 El Centro地震波Fig.8 El Centro seismic waves408第 3 期王焱楠，等：改进 VMD及其在时变结构地震响应瞬时频率识别中的应用 500，40000。基于综合目标函数 F 搜寻最优值，其迭代收敛曲线如图 10 所示。由图 10 可知，迭代23 次时收敛，最优值Fmin=1.7506，对应的值为10142.9。采用确定的优化参数 K=4，=10142.9，基于IVMDHT 算法识别得到地震激励下时变框架结构的瞬时频率，如图 11 所示。此外，图 11 中还给出了根据“冻结法”确定的瞬时频率理论值，即假设结构物理参数在每个时间间隔内保持不变，得到每个时间间隔内的瞬时频率。由图 11可见，识别得到的框架结构瞬时频率与理论瞬时频率基本吻合，无模态混叠。其中，第 2 阶瞬时频率识别值与理论值吻合最好；第 3，4阶瞬时频率在部分时刻存在突变，这是由于 HT 计算瞬时频率的公式误差造成的11。采用移动平均（Moving Average）法，即选定窗口长度为80，将窗口内的瞬时频率数据进行算数平均处理后，第 3，4阶瞬时频率的移动均值与理论值基本一致。当 K=4 时，表 2 给出了不同取值，VMD 算法识别得到的地震激励下时变框架结构瞬时频率的RMSE 值。由表 2 可知，当采用优化确定的参数（即=10142.9）时，识别得到的结构瞬时频率的RMSE 值最小，结果最准确；当随机取值偏小（即=2000）时，识别的瞬时频率精度较低；当随机取值偏大（即=25000）时，模态混叠导致第 4 阶瞬时频率无法被识别。由此表明，采用优化取值参数的IVMDHT 算法，可用于识别地震激励下时变框架结构的瞬时频率，且结果最准确，精度最高。4试验验证为试验验证 IVMD 算法识别地震激励下时变结构瞬时频率的有效性和准确性，采用同济大学提供的 12 层钢筋混凝土框架模型模拟地震振动台试验数据27。模型缩尺比例为 1/10，总高度 3.6 m，每层高 0.3 m，平面结构尺寸为 0.6 m0.6 m，如图 12所示。模型试验采用了 4 种地震波，测试了 61 个工况，分别在模型底层、2层、4层、6层、8层、10层和 12层的 X，Y，Z 方向布置了 23个加速度传感器。模型的详细几何尺寸、物理参数和模拟地震振动台试验情况参见文献 27。本研究仅分析 Q1 和 Q2 两种工况下模型在 X 向 El Centro 地震激励下结构顶层的加速度响应，Q1 和 Q2 工况地震激励加速度峰值（PGA）分别为 0.646g 和 0.775g，信号采样频率为255.102 Hz。为更精确地识别地震激励下模型的瞬时频率，研究仅考虑模型的前 3 阶模态，采用截止频率为15 Hz 的 6 阶切比雪夫型低通滤波器对模型地震响应信号进行重采样，重采样频率为 100 Hz。图 13所示为 Q1 和 Q2 两种工况下的模型顶层加速度响应及其傅里叶频谱和 MUSIC 功率谱。图 13（b），图 10 综合目标函数 F迭代收敛曲线（地震激励下时变框架结构）Fig.10 The iterative convergence curve of composite objective function F（a time varying frame structure under seismic excitation）图 11 IVMD-HT 算 法 识 别 的 地 震 激 励 下 时 变 框 架 结 构瞬时频率Fig.11 Instantaneous identified frequencies of a time varying frame structure identified by IVMD-HT algorithm under seismic excitation表 2 不 同取 值 下 框 架 结 构 瞬 时 频 率 识 别 值 的 RMSE（K=4）Tab.2 RMSE of identified instantaneous frequencies of a frame structure with different （K=4）1阶2阶3阶4阶总和RMSE=20000.110.060.220.411.09=10142.90.100.060.220.370.75=250000.110.070.50-注：“-”表示该方法未能识别相应瞬时频率。409振 动 工 程 学 报第 37 卷（e）所示的傅里叶频谱中存在虚假峰值和波动干扰，峰值多而杂乱，难以准确确定分解模态数 K；图 13（c），（f）所示的 MUSIC 功率谱中可清楚地观察到 3个峰值，由此确定分解模态数 K=3。确定分解模态数 K 后，设置 SSA 算法的种群个数 为 10，最 大 迭 代 次 数 为 100，搜 索 范 围 为 500，40000。基于综合目标函数 F 搜寻最优值，其迭代收敛曲线如图 14 所示。由图 14 可知，Q1 和Q2 工况下分别迭代 57 次和 58 次后收敛，最优值分别为 4246.7和 6125.7。Q1 和 Q2 工况下，分别采用确定的优化参数组合K=3，=4246.7和K=3，=6125.7，基 于IVMDHT 算法识别得到 12层框架模型的前 3阶瞬时频率，如图 15 所示。由图 15 可见，Q1 和 Q2 工况下框架模型前 3阶瞬时频率成分被准确地识别。由文献 27 可知，模型结构在模拟地震振动台试验的过程中出现了损伤，其频率具有时变特性。采用移动平均法（窗口长度为 100）处理瞬时频率曲线，消除 HT 计算瞬时频率时出现的负频率问题，得到的瞬时频率移动均值（见图 16）真实地反映了地震激励下模型结构频率的时变规律。由图 16可见，Q2 工况下识别得到的模型瞬时频率值低于 Q1 工况，且波动幅度较大。原因在于，相比 Q1 工况，Q2工况下模型损伤更严重，结构刚度下降更多，这与试验 现 象 基 本 一 致。由 此 表 明，采 用 优 化 参 数 的IVMD 算法，可识别地震激励下框架结构损伤引起的频率时变特征，对比不同损伤工况下结构瞬时频率的变化趋势即可得到结构的损伤程度。图 12 框架结构振动台模型Fig.12 Shaking table model of frame structure图 13 Q1和 Q2工况下模型顶层加速度响应及其傅里叶幅值谱和 MUSIC功率谱Fig.13 Acceleration timehistory，Fourier amplitude spectrum and MUSIC power spectrum of top layer of model under Q1 and Q2 cases410第 3 期王焱楠，等：改进 VMD及其在时变结构地震响应瞬时频率识别中的应用5结 论提出了一种基于 IVMDHT 算法识别时变结构地震响应瞬时频率的方法。采用模拟信号、地震激励下 4层时变框架结构数值模拟响应和 12层框架结构模型地震模拟振动台的试验数据验证了所提方法的有效性和实用性。主要结论如下：（1）相比傅里叶频谱和短时傅里叶谱，基于MUSIC功率谱确定的 VMD模态数 K更准确；基于整体正交系数和能量比系数构建的综合目标函数，结合SSA算法可更有效地优化确定二次惩罚因子。（2）基于优化确定 K 和参数的 IVMD 算法，解决了人为预设参数时 VMD 存在的模态混叠和噪声干扰等问题，提高了计算效率；提高了结合 HT 识别时变结构地震响应瞬时频率的精度。（3）应该指出，二次惩罚因子的优化搜索范围仍需通过试算人为确定。参考文献:1Xin Y，Li J，Hao H.Enhanced vibration decomposition method based on multisynchrosqueezing transform and analytical mode decomposition J.Structural Control and Health Monitoring，2021，28（6）：e2730.2 Loh C H，Chao S H，Weng J H，et al.Application of subspace identification technique to longterm seismic response monitoring of structures J.Earthquake Engineering&Structural Dynamics，2015，44（3）：385-402.3Ni Y C，Zhang F L.Fast Bayesian frequency domain modal identification from seismic response data J.Computers&Structures，2019，212：225-235.4Ghahari S F，Abazarsa F，Ghannad M A，et al.Response-only modal identification of structures using 图 14 综合目标函数 F迭代收敛曲线Fig.14 The iterative convergence curves of composite objective function F图 15 IVMD-HT算法识别的12层框架模型前3阶瞬时频率Fig.15 The first three order instantaneous frequencies of a 12-story frame model identified by IVMD-HT algorithm图 16 Q1和 Q2工况下前 3阶瞬时频率移动均值Fig.16 Moving mean of the first three order instantaneous frequencies under Q1 and Q2 cases411振 动 工 程 学 报第 37 卷strong motion data J.Earthquake Engineering&Structural Dynamics，2013，42（8）：1221-1242.5Guo Y L，Kareem A.System identification through nonstationary data using timefrequency blind source separationJ.Journal of Sound and Vibration，2016，371：110-131.6Loh C H，Chen J D.Tracking modal parameters from building seismic response data using recursive subspace identification algorithmJ.Earthquake Engineering&Structural Dynamics，2017，46（13）：2163-2183.7Pioldi F，Rizzi E.Earthquakeinduced structural response outputonly identification by two different Operational Modal Analysis techniques J.Earthquake Engineering&Structural Dynamics，2018，47（1）：257-264.8Chen W H，Hseuh W，Loh K J，et al.Damage evaluation of seismic response of structure through time-frequency analysis technique J.Structural Monitoring and Maintenance，2022，9（2）：107-127.9樊海涛，何益斌，周绪红.基于 Hilbert-Huang变换的结构损伤诊断方法研究J.建筑结构学报，2006，27（6）：114-1