基于改进经验包络法的非线性系统参数识别

曾 华， 许福友

(大连理工大学 建设工程学部， 辽宁 大连 116024)

模态参数(包括模态频率、模态阻尼等)是土木(工程结构系统的重要参数，通常通过结构动力试验(包括强迫振动、自由振动和环境激励下的随机振动)获得的振动(位移或加速度)信号来识别。由于大型复杂土木工程结构的强迫振动试验，对激励设备和能量要求较高，且过大的强迫激励会对结构安全与运营造成，难以开展强迫振动试验，因此常采用无需对输入进行测量的自由振动试验或环境激励下的随机振动试验，这种只对输出进行测量的模态识别一般被称作基于输出的模态参数识别(output-only modal analysis, OMA)方法。

OMA方法一般可分为时域分析法、频域分析法和时频分析法三类。其中，时域分析方法是对结构振动时程响应直接分析，如随机减量技术(random decrement technique, RDT)[1]，随机子空间法(stochastic subspace identification, SSI)[2]等；频域分析方法一般通过傅里叶变换到频域再进行模态参数识别，主要有峰值拾取法[3]，频域分解法[4-5]等；时频分析法一般有基于小波变换(wavelet transform, WT)[6]、基于希尔伯特变换(Hilbert transform, HT)[7-8]、基于经验分解(empirical mode decomposition, EMD)[9-10]、基于极点对称模态分解(extreme-point symmetric mode decomposition, ESMD)[11]、基于变分模态分解法(variational mode decomposition, VMD)[12-15]等信号分析方法。其中，WT方法受母小波选择的影响，人为因素较大，没有通用的小波函数。EMD方法是一种自适应的信号分解技术，它与HT方法结合成为一种强大的时频分析技术，通常被称为希尔伯特-黄变换(Hilbert-Huang transform, HHT)[16]，广泛运用于机械土木等领域，并有不错识别的效果。王金良等在HHT基础上提出一种自适应分解的ESMD法，该方法通过极点对称构造内部的局部对称均线，能够降低插值带来的不确定性，但同样存在模态混叠问题。VMD方法是通过迭代搜索变分模型最优解来获取每个分量的中心频率与带宽，克服了EMD方法的模态混叠等缺陷，并已应用在非线性系统模态参数识别及非平稳信号分析等问题中。

对于土木结构的振动系统，其模态参数往往表现为非线性函数。其中瞬时频率和瞬时阻尼比是非线性系统最为关注的模态参数。求解信号瞬时特性的常见方法有：相位差法，能量算子法，小波变换法，Hilbert变换法(HT)以及经验包络法。张贤达[17]提出相位差分法，对线性系统具有无偏估计、零延迟等优点，但对噪音过于敏感。Maragos等[18]利用能量算子将单分量信号分离为调幅和调频分量，并总结了几种有效的瞬时频率估计方法，具有原理简单、计算量小等优点，广泛应用[19]。但能量算子法基于线性相位假设，因此应用于非线性系统时将会产生很大误差。Staszewski[20]运用系统响应首先分解得到时程小波脊和小波骨架，然后求解信号的时程瞬时特性，最后识别非线性系统参数，具有较好的精度。Feldman[21-22]详细论述了HT在非线性系统识别中的应用，该方法会产生明显端点效应、负频率等现象。

郑近德等[23]基于经验调幅调频分解，提出了经验包络法(empirical envelope, EE)，并验证了EE在求解瞬时频率具有计算量小，识别精度高，操作简单等优点。不过，在整个求解过程中存在两点不足：① 调幅调频分解中，为了得到较纯的调频信号，需要多次包络迭代，这样将会导致样条插值的累计误差的增大；② 对信号求导时，由于信号不光滑，可能造成过包络问题，即包络线与信号线相交，导致包络误差增加。因而，EE用于信噪比较小的非线性结构信号时，常出现识别结果不稳定现象。

针对上述两个问题，本文采用滑窗阈值去噪思想与滑动平均技术对其改进，提出了改进的经验包络法(improved empirical envelope, IEE)。通过多种信号参数识别与分析，全面证实了该方法利用自由振动试验识别模态参数的有效性。

1 经验调幅调频分解

任意单频信号y(t)均可分解为调幅函数和调频函数的乘积形式，即：

y(t)=a(t)cos(φ(t))

(1)

式中：a(t)为调幅函数；cos(φ(t))为调频函数；φ(t)为瞬时相位。

经验调幅调频分解具体过程如下

步骤1对y(t)取绝对值|y(t)|，搜寻|y(t)|的所有极值点(ti,|yi|)(i=1,2…,N)，i为极值点编号，N为极值点的个数

步骤2采用三次样条拟合这些极值点(ti,|yi|)，得到|y(t)|的第一轮经验包络函数曲线a1(t)

步骤3对|y(t)|归一化为

y1(t)=|y(t)|/|a1(t)|

(2)

理想条件下，y1(t)应小于等于1。由于y(t)可能存在明显的噪声甚至野点，因此并不能保证在任何时刻a1(t)>|y(t)|，即某些时刻y1(t)有可能大于1。此时，将y1(t)视为y(t)，重复步骤1～2，经n次迭代

y2(t)=y1(t)/|a2(t)|

(3)

……

yn(t)=yn-1(t)/|an(t)|

(4)

直到yn(t)≤1，迭代结束，此时的yn(t)为y(t)调频部分。存在φ(t)，使得：

yn(t)=cos(φ(t))

(5)

y(t)调幅部分可定义为

(6)

至此，原信号y(t)可以分解为调幅调频两部分乘积a(t)cos(φ(t))形式。

以上经验调幅调频分解至少存在两方面不足：

(1) 由于噪声的存在，为了得到纯调频信号，需要多轮包络迭代，这样将会导致样条插值的累计误差增大。

(2) 三次样条插值存在端点效应。

下节将针对这两方面的不足进行改进，提高分解精度。

2 调幅调频分解的改进

2.1 包络迭代的改进

采用滑窗阈值去噪思想筛选剔除异常极值点。具体步骤如下：

步骤1搜寻|y(t)|的所有极值点，并确定其极值点(ti,|yi|)(i=1,2…,N)，i为极值点编号，N为极值点的个数

步骤2对|y(t)|极值点集U={(ti,|yi|)}(i=1,2,…,N)，按“窗体”极值点个数为W1进行滑动，并选取第1个子集Q1={(tj,|yj|)}(j=1,2,…,W1)

步骤3对子集Q1进行最小二乘法拟合

(7)

拟合优度定义为

(8)

步骤4搜索子集Q1内所有异常极值点

η<ηc

(9)

步骤5“窗体”向前滑动一个极值点，选取得到子集第2个Q2={(tj,|yj|)}(j=2,3,…,W1+1)，按同样方法搜寻异常极值点集合ERR2

步骤6“窗体”继续滑动，直至选取最后一个子集QN-W1+1={(tj,|yj|)}(j=N-W1+1,N-W1+2,…,N)，并搜寻异常极值点集合ERRN-W1+1

后续过程与调幅调频分解相同，此处不再赘述。需要说明的是，“窗体”大小W1和阈值因子ηc应该选择适当，过大过小都会影响最后结果。“窗体”(阈值因子)过大将会舍去大量的极值点，出现欠包络现象，从而导致调幅调频信号的失真；“窗体”(阈值因子)过小会导致极值点筛选不充分，出现过包络现象，从而导致调幅调频信号不稳定。一般而言，对于信噪比较小、质量较差的信号，可适当增加W1与ηc的值；对于信噪比较大、质量较好的信号，可适当减小W1与ηc的值。通过对大量理论模拟信号分析可知，“窗体”大小可以在(8～20)范围内取值，阈值因子可在(0.3～0.9)范围内取值。在此范围内，W1与ηc的选取对结果影响较小。

2.2 端点效应的改进

本文采取极值点向外线性延伸手段，来抑制端点效应，具体步骤如下：

步骤1取待拟合的极值点的前两个极值点(t1,|y1|)与(t2,|y2|)；

步骤2令Δ=(|y1|-|y2|)/(t1-t2)，计算Δy=|y1|-t1Δ；

步骤3如果Δy>|y(0)|，则首部向外延伸点为(0，Δy)，否则为(0，|y(0)|)；

步骤4尾部向外延伸点可按类似方法计算，此处不再赘述。

3 经验包络法

获得调幅调频函数后，对调频函数求导可以获得新函数，然后再对该新函数进行经验调幅调频分解，最终获得原始信号y(t)的瞬时频率与瞬时振幅。具体步骤如下：

步骤1对F(t)=cos(φ(t))求导可得F′(t)=-φ′(t)sin(φ(t))，y(t)的瞬时频率f(t)可表示为

(10)

步骤2再对-φ′(t)sin(φ(t))进行调幅调频分解可得

F′(t)=b(t)cos(φ(t))

(11)

式中，b(t)可看作φ′(t)。因此原信号的瞬时频率为

(12)

经验包络法相对小波变换法、Hilbert变换法有计算量小，识别精度高，操作简单等优点[24]。但仍存在问题：对信号求导，如果信号不光滑，可能造成过包络问题，即包络线与信号线相交，导致包络误差增加。

4 经验包络法的改进

本文采用滑动平均算法对经验包络法步骤2中调幅调频分解求导后的信号F′(t)进行滑动平均处理，去除因求导之后信号混有不光滑的成分。具体操作步骤如下：

步骤1将F′(t)看作采样信号S(t)；

步骤2选取固定长度为W2的采样数据队列；

需说明的是，应选取合适的W2，其值越大平滑效果越好，识别结果越稳定，但过大会导致信号失真。一般而言，对于采样频率较高、信噪比较小、质量较差的信号，可适当增加W2的值；采样频率较低、信噪比较大、质量较好的信号，可适当减小W2的值。通过对大量理论模拟信号分析可知，W2可在(2～20)内取值。

为了说明对调幅调频分解与经验包络法改进后的效果，现给出示意图加以说明。如图1所示，给出了异常极值点分别采用调幅调频分解与改进调频调幅分解得到的包络误差影响的示意图。由图1可知，通过对异常极值点的剔除，避免了过包络现象，减小了包络误差，进而减少了包络迭代次数。图2给出了求导后信号不光滑性对包络误差的影响的示意图。由图2可知，通过对求导后信号的滑动平均处理，可以有效的减小包络误差。

(a) 调频调幅包络

(b) 改进调幅调频包络图1 异常极值点对包络影响示意图Fig.1 Schematic diagram of influence of abnormal extreme point on envelope

(a) EE包络

(b) IEE包络图2 求导信号不光滑对包络插值影响示意图Fig.2 Schematic diagram of influence of unsmooth derivative signal on envelope difference

5 模态参数识别

由单分量自由衰减信号的瞬时振幅a(t)可导出阻尼系数为

(13)

式中，m为离散采样极值点。故单分量自由衰减信号的瞬时阻尼比ξ(t)可以表示为

(14)

在低阻尼下，单分量自由衰减信号的瞬时阻尼比可近似计算为

(15)

如图3所示，为绝对误差ξr随阻尼比的变化。其中ξr可按下式给出

图3 阻尼比的误差Fig.3 Relative error of damping ratio

ξr=ξ′-ξ

(16)

式中，ξ为精确结果由式(14)算出，ξ′为近似结果由式(15)算出。

由图3可知，当ξ<0.1时，通过式(15)计算，ξr低于0.000 5，因此单分量自由衰减信号在低阻尼情下，采用式(15)是合理并且方便的。

这里需要说明的是，式(13)～(15)仅适用于单分量自由衰减振动，故而本文提出的非线性系统模态参数识别的方法也仅适用于单分量自由衰减振动。但本文方法可以与VMD、EMD、ESMD、RDT等方法结合，形成一种比较通用的非线性系统识别方法。例如，对于单模态自由振动，可直接利用本文方法识别时变模态参数；对于多模态自由振动，可以现利用带通滤波器(如若存在密集模态问题，该方法欠佳)或者VMD等分解技术将多模态解耦，再利用本文方法识别时变模态参数；对于多模态环境激励[25]，可先通过功率谱分析确定合适的带宽，使用带通滤波器将其分解为多个窄带信号，再对这些窄带信号使用VMD等技术进行解耦，并利用RDT技术将其转化为单分量自由衰减信号，最后利用本文提出的方法识别时变模态参数。下文通过几个算例进一步说明。

6 算 例

本节采用三个算例对比验证EE与IEE识别效果，下文均为端点效应处理后的识别结果。

6.1 算例一

某单自由度非线性系统阻尼比ξ=0.005+0.001A，圆频率为ω=2π(2-0.01A)，其中A(0.5≤A≤5)为振幅，采样频率为500 Hz，初始振幅A0=5 cm，通过Newmark-β法可计算非线性系统的自由衰减振动响应。

对系统响应施加信噪比为80 dB零均值的高斯白噪声n，如图4所示。对于信噪比可表示为

SNR=10lg(Ps/Pn)

(17)

图4 高斯白噪声Fig.4 White Gaussian noise

式中：Ps为信号的平均功率；Pn为噪声的平均功率。添加噪声后的模拟信号如图5所示。

图5 模拟信号Fig.5 Analog signal

对信号A(t)(为保证识别的准确性，截去信号的前三个波峰)分别采用EE与IEE(W1=14,W2=4，ηc=0.5)识别计算，可得瞬时频率与瞬时阻尼比，分别如图6和图7所示。由图6(a)可知，EE识别的瞬时频率波动范围较大，且在小振幅区间偏离真值。这是由于噪声对小振幅区域信号影响更加明显，导致小振幅区间出现大量异常极值点，进而在经验调幅调频分解时，出现过包络现象既增加了三次样条拟合误差，又增加了包络迭代次数，最终导致小振幅处识别结果偏差很大，甚至错误。由图6(b)所示，IEE识别的瞬时频率波动范围较小，识别结果更稳定。

(a) EE识别结果

(b) IEE识别结果图6 幅变频率Fig.6 Variable instantaneous frequency

如图7(a)所示，由于EE识别的瞬时频率在小振幅处有异常波动，直接影响了瞬时阻尼比的识别结果，导致其结果在小振幅处偏离真值。如图7(b)所示，IEE识别的瞬时阻尼比非常稳定，其相对优势显而易见。

(a) EE识别

(b) IEE识别图7 幅变阻尼比Fig.7 Variable instantaneous damping ratio

为了进一步说明噪声强度与噪声均值对识别结果的影响。现对模拟信号施加不同信噪比，瞬时频率识别误差如表1所示。施加不同均值的噪声，瞬时频率识别误差如表2所示。

表1 均值为零不同信噪比瞬时频率识别结果Tab.1 Recognition results of instantaneous frequency with zero mean value and different SNR

表2 SNR为80不同均值瞬时频率识识别结果Tab.2 Recognition results of instantaneous frequency with SNR of 80

瞬时频率识别误差R按下式给出

R=|f′-f|/f×100%

(18)

“窗体”大小W1和阈值因子ηc对信号A(t) 频率识别误差如表3、表4所示。

表3 不同“窗体”大小IEE识别频率误差Tab.3 Frequency error identification by IEE with different sliding window width

表4 不同阈值因子IEE识别频率误差Tab.4 Frequency error identification by IEE with different threshold factors

由表3和4可知，“窗体”大小与阈值因子对识别结果影响较小，因此，本文推荐取值范围是合理的。

6.2 算例二

本信号为单自由度没水圆柱纵摇自由衰减振动角位移时程，如图8所示。

图8 数值模拟信号Fig.8 Numerical simulation signal

对数值模拟信号分别采用EE与IEE(W1=10,W2=8,ηc=0.8)识别计算，可得瞬时频率和瞬时阻尼比。如图9所示，EE识别的瞬时频率出现较大的波动，导致识别结果不可用。这是因为数值误差的影响，对信号求导时，出现过包络现象，即包络线与信号线相交(见图2(b))，导致包络迭代进入死循环。而IEE识别的瞬时频率无任何异常，识别效果很好。

图9 瞬时频率Fig.9 Instantaneous frequency

如图10所示，为两种方法识别幅变瞬时阻尼比。图10(a)为摆角大于1°的识别结果，由于EE识别的瞬时频率波动的影响，导致EE识别阻尼比结果波动较大；IEE识别结果随摆角逐渐增大，结果较为稳定。图10(b)，为摆角小于1°EE的识别结果，由于EE识别的瞬时频率在小摆角区间异常，导致阻尼比识别结果不可用。图10(c)，为摆角小于1°IEE的识别结果，可以看出IEE在小摆角区间识别效果也较好。

6.3 算例三

为说明本文方法能应用于实际结构当中，现以苏通大桥为背景建设的自然风场大比例试验模型[26]为研究对象，试验模型效果图如图11所示。主梁断面布置22个加速度传感器，分别为AL-1～AL-11(主梁左侧)、AR-1～AR-11(主梁右侧)如图12所示。

图11 试验模型效果图Fig.11 Effect diagram of test model

图12 传感器布置示意图Fig.12 Schematic diagram of sensor arrangement

6.3.1 广义单自由度自由衰减振动

为了获取模型一阶对称竖向振动模态参数，在无风环境下，对主梁跨中施加初始位移使其作自由衰减振动，以跨中AL-6加速度(采样频率为400 Hz)响应为例，如图13所示。

图13 AL-6加速度响应Fig.13 AL-6 acceleration response

为保证识别的精度，现对该信号采用低通滤波器(截止频率为2 Hz)滤波处理，并截去前10 s与后15 s响应。如图14所示。

图14 滤波后加速度信号Fig.14 Filtered acceleration signal

现对滤波后信号分别采用EE与IEE(W1=10,W2=10,ηc=0.8)识别计算，可得瞬时频率与瞬时阻尼比。图15(a)为t< 64 s两种方法识别瞬时频率的结果，可以看出两种方法识别结构较为接近。图15(b)为t> 64 s EE识别瞬时频率的结果，可以看出识别结果出现大幅波动。这是因为滤波后残余噪声以及信号非平稳等因素的影响，对信号求导时，出现过包络现象，即包络线与信号线相交(如图2(b)所示)，导致包络迭代进入死循环。图15(c)为t> 64 s IEE识别瞬时频率的结果，可以看出识别结果较为稳定。

图16(a)为a> 0.2 m/s2两种方法识别幅变阻尼比的结果，可以看出两种方法识别结构较为接近。图16(b)为a< 0.2 m/s2EE识别幅变阻尼比的结果，由于EE识别的瞬时频率波动的影响，导致识别结果不可用。图16(c)为a< 0.2 m/s2IEE识别幅变阻尼比的结果，可以看出识别结果较为稳定。

为了进一步说明IEE的有效性，现将识别的模态频率与模态阻尼比代入单自由度动运动微分方程，运用Newmark-β法计算结构位移响应。为了比较的精准性，只截取30 s信号加以对比分析，如图17、图18所示。由图可知，通过IEE识别的模态参数反算结构位移响应更加接近真实的位移响应，验证了IEE的有效性。

图17 EE计算结构响应Fig.17 EE calculation of structural response

图18 IEE计算结构响应Fig.18 IEE calculation of structural response

6.3.2 多自由度自由衰减振动

为了获取模型高阶竖向振动模态参数，在无风环境下，对主梁八分点施加初始位移使其作自由衰减振动，以AL-3加速度(采样频率为400 Hz)响应为例，如图19所示。由图可知，成功激励出高阶竖向振动模态。通过幅值谱分析(如图20所示)，可初步判断该信号存在两阶模态，其频率大约在2.5 Hz与5 Hz左右。

图19 AL-3加速度响应Fig.19 AL-3 acceleration response

图20 加速度响应幅值谱Fig.20 Acceleration response amplitude spectrum

为保证识别的精度，现对该信号采用带通滤波器(通带频率为2 Hz、5.5 Hz)滤波处理，并截去前10 s与后15 s响应。如图21所示。

图21 滤波后加速度信号Fig.21 Filtered acceleration signal

通过VMD法将该信号分解为两个单模态响应(分别记为IMF1、IMF2，结果如图22所示)，具体多模态参数分解过程，本文不重点展开讨论，可参照文献[27]。这里也可以使用带通滤波器将两个模态分解，但受带通滤波器阻带衰减率的影响，当两个模态频率较为接近时，难以用带通滤波器将其分开。故而，VMD等技术在处理模态耦合问题中更加通用。

(a) IMF1模态

(b) IMF2模态图22 单模态响应Fig.22 Single mode response

现采用IEE法分别对两个单模态自由衰减信号进行识别计算，可得瞬时频率与瞬时阻尼比。如图23所示，成功识别出了两阶模态频率。如图24所示，成功识别出了两阶模态阻尼比，但对于IMF2识别的阻尼比波动较大，这是由于利用VMD重构信号时，部分信号振幅失真导致的。

图23 瞬时频率Fig.23 Instantaneous frequency

图24 幅变阻尼比Fig.24 Amplitudevarying damping ratio

7 结 论

针对非线性系统瞬时频率问题，提出了改进经验包络法，其关键技术包括：滑窗阈值去噪与滑动平均。基于该方法对单自由度(一般通过经验模态分解、变分模态分解等方法解耦获取)非线性系统自由衰减振动(一般通过自由振动试验或通过随机减量法从结构随机振动响应中获取)模态参数识别(瞬时频率、瞬时阻尼)。通过几个算例参数识别分析，证实了本文方法在利用自由振动试验进行模态参数识别具有良好的抗噪声性能与较高的识别精度。