基于AR-EWT和SSI的桥梁结构模态参数识别

黄天立，詹晨路,2，万熹，赵玉印，任伟新

(1. 中南大学 土木工程学院，湖南 长沙，410075；2. 抚州赣东公路设计院有限公司，江西 抚州，344000；3. 中国中铁二局集团有限公司，四川 成都，610031；4. 深圳大学 土木与交通工程学院，广东 深圳，518060)

随机子空间识别 (stochastic subspace identification，SSI)方法仅需要实测结构动力响应便可识别得到结构的频率、阻尼比和振型等模态参数，在计算过程中不涉及迭代和收敛性问题，具有较好的适用性和较高的识别精度，是目前桥梁结构模态参数识别中应用最广的方法之一[1-2]。然而，SSI 在识别结构模态参数时存在模型定阶、判别并剔除虚假模态等问题[3-6]。采用经验模式分解(empirical mode decomposition，简称EMD)[7]、变分模态分解(variational mode decomposition，简称VMD)[8]等信号分解方法预处理信号后，再利用SSI 识别结构模态参数，为SSI 存在的模型定阶问题、降低模态间相互影响对参数识别的干扰提供了可行的解决方法[9-12]。YU 等[9]提出了结合EMD和SSI的桥梁结构模态参数识别方法，识别了环境激励下西宁北川河桥的前3阶模态参数，研究表明该方法能降低虚假模态和噪声对识别结果的影响。黄天立等[10]针对EMD 方法信号分解能力不足和固有模式函数分量不正交等问题，提出了结合波组信号前处理、正交化经验模式分解和Hilbert 变换的密集模态结构参数识别方法。杨勇等[11]通过解析模式分解将模态密集分布的振动信号分解为一系列单模态子信号，结合随机减量技术提出了大型桥梁结构的模态参数识别方法。孙猛猛等[12]通过VMD 分解将实测振动信号分解为单模态信号，结合自然环境激励技术和直接插值法识别了高层建筑结构自振频率和阻尼比。需指出的是，EMD方法本质上属于经验分解，没有严格的数学理论基础。VMD 分解需要确定分解层数、初始中心频率、限制带宽等参数，参数选取对VMD的分解效果有较大影响。为更好地分解多分量信号，GILLES 等[13]提出了经验小波变换(empirical wavelet transform，简称EWT)自适应多分量信号分解方法。EWT 通过自适应划分频带边界，构建小波滤波器组提取多分量信号中的不同模态分量，由于其计算高效、理论完备，EWT 方法已应用于土木工程结构的模态参数识别和状态评估等领域[14-15]。研究表明[16-17]，EWT 分解高信噪比信号的效果较好，而分解低信噪比信号和非平稳信号的效果较差。这是由于EWT 分解信号时需根据信号的傅里叶谱划分频带，常存在频带划分不准确甚至错误划分频带的情况，导致无法正确分离代表不同频率成分的分量信号。万熹等[18]提出采用自回归(autoregressive，简称AR)功率谱替代傅里叶频谱作为EWT 频带划分的依据。该方法即为AR-EWT 方法，其提高了EWT 自适应分解信号的能力，结合随机减量技术和Hilbert 变换，可识别桥梁结构的自振频率和阻尼比，但是该方法不能识别结构振型。

本文采用信号分解方法预处理信号后再利用SSI识别结构模态参数，提出一种环境激励下基于AR-EWT和SSI的桥梁结构模态参数识别方法。该方法采用AR-EWT 将结构不同测点的振动响应信号分解成一系列仅含单一模态信息的单分量信号，然后，组装包含相同模态信息的单分量信号，利用SSI识别结构的自振频率、阻尼比和振型。采用ASCE Benchmark 模型、简支梁模型和环境激励下西宁北川河桥的模态参数识别结果验证所提方法的有效性。

1 基本原理

1.1 基于自回归功率谱的经验小波变换

经验小波变换(EWT)方法通过自适应设计的小波滤波器组，从信号中提取一系列具有紧凑支撑傅里叶频谱的调幅-调频单分量成分[13]。自适应设计小波滤波器组的关键在于基于傅里叶谱的频带划分，其频带划分如图1所示。基于Shannon准则将信号傅里叶频谱的频率范围标准化到区间[0,π]，根据频带划分方法将其分割为连续的N段，每个区间表示为Λn=[ωn-1,ωn]，其中，Λn为频带范围，ωn、ωn-1分别为频带上、下限(n=1,2,…,N)。令ω0=0，ωN=π，以ωn为中心，定义宽度为2τn的过渡区域，如图1中的阴影部分所示。

图1 傅里叶谱频带划分Fig.1 Frequency band division of Fourier spectrum

根据Little wood-Paley 和Meyer 理论构造小波的思路，按照式(1)和式(2)分别定义EWT中的经验尺度函数ϕn(ω)和经验小波函数ψn(ω)，将每个频带Λn的边界频率代入，构建带通滤波器组。

式中：ω为频率；ωn为第n个边界频率；τn=γωn；γ为确保ωn-1和ωn之间不存在重叠的转换参数，

类似于构造小波变换，将信号分别与经验尺度函数ϕn(ω)和经验小波函数ψn(ω)内积得到近似系数Wsε(0,t)和细节系数Wsε(n,t)，如式(4)和式(5)所示(其中，近似系数代表信号低频成分，细节系数代表信号高频成分)。

式中：F表示傅里叶变换，F-1表示傅里叶逆变换；分别为F[ψ1(ω)] 和F[ψn(ω)]的复共轭。

叠加低频和高频成分的信号分量即可重构原始信号，如式(6)所示。

式中：s0(t) =Wsε(0,t)*ϕ1(t)，sk(t) =Wsε(k,t)*ψk(t)，分别为低频和高频成分的信号分量；“*”表示卷积运算。

应该指出的是，当多分量信号中含有较强噪声时，基于傅里叶频谱难以有效划分信号的频带边界，这将直接影响信号分解的正确性。万熹等[18]指出将信号的AR功率谱代替傅里叶频谱进行频带划分，能更准确地估计频带边界，从而更好地分离信号的各阶模态分量。

假设某线性系统受到均值为0、方差为σ2的高斯白噪声激励w(n)，输出信号序列为x(n)，则其AR 模型应满足方程其中，k=1, 2, …,p，p为预设的模型阶数，ak为AR 模型预测系数。基于Burg 算法求解AR 模型参数ak，代入式(7)得到信号序列x(n)的AR 模型自回归功率谱PAR，具体求解过程见文献[17]。

1.2 协方差驱动随机子空间识别算法

随机子空间算法主要包括协方差驱动随机子空间法(Covariance ‐ driven Stochastic Subspace Identification，SSI-Cov)和数据驱动随机子空间法(Data-driven Stochastic Subspace Identification，SSIData)共2种，其中SSI-Cov根据测量响应数据的协方差矩阵求系统状态矩阵，在一定程度上可提高识别精度，因此，本文采用SSI-Cov算法进行模态参数识别。

1.2.1 数学模型

对于自由度为n1的振动系统，其在物理坐标下的振动方程为

式中：M∊Rn1×n1,C1∊Rn1×n1,K∊Rn1×n1，分别为系统的质量、阻尼和刚度矩阵；分别为系统的加速度、速度和位移列向量；f(t) ∊Rn1×1，为n1阶外部激励列向量；B1∊Rn1×m1，为广义输入矩阵；m1为输入的点数；u(t) ∊Rm1×1，为随时间变化的输入列向量；R为实数集。

定义状态向量式(9)，得到随机子空间离散状态空间模型，如式(10)所示。

式中：xk=x(kΔt)，为离散的时间状态向量；Δt为采样周期；k为采样点个数；A和C分别为离散系统状态矩阵和输出矩阵；ωk和vk分别为处理过程、建模误差引起的噪声和测量向量，一般假定为均值为零且互不相关的白噪声。

1.2.2 系统定阶

SSI-Cov通过构建Toeplite矩阵求解系统矩阵A和C，其中，Toeplite 矩阵的秩即为系统阶次。确定系统阶次的方法主要有2种。

一种为基于稳定图的定阶法。首先，假定模型阶次为N(N大于模型真实阶次)，不同阶次对应不同的状态空间模型，从而能得出不同的模态参数识别结果(模态频率、阻尼比和振型)；然后，将所有不同的结果绘制在稳定图上，对比高一阶次与低一阶次的模态参数识别结果，若两者差异不超过预设限值，则将高一阶次对应的参数识别结果设置为稳定点；最终，所有稳定点形成一系列稳定轴，1条稳定轴代表系统的一阶模态，结构的真实阶次可取为稳定轴对应的系统阶次。研究表明[6]，该方法并不能完全排除噪声或非平稳激励的影响。

另一种为基于奇异值跳跃的定阶法。首先，对Toeplite矩阵进行奇异值分解，并将得到的奇异值绘制成图；然后，根据奇异值图人工确定奇异值跳跃点，确定系统阶次为跳跃点前所有奇异值数之和的一半。研究表明[6]，该方法不适用于复杂工程结构。

1.2.3 模态参数识别

将系统矩阵A、C特征值进行分解即可获得系统模态参数。对离散系统状态矩阵A特征值进行分解：

式中：Ω=diag[μi] ∊Cn2×n2，为由离散时间复特征值μi组成的对角阵；n2=2n1，为模型阶数；Ψ∊Cn2×n2，为以特征向量为列向量的矩阵。

对连续系统状态矩阵Ac特征值进行分解：

式中：Ωc=diag[ρi] ∊Cn2×n2，为由离散时间复特征值ρi组成的对角阵；Ψc∊Cn2×n2，为以特征向量为列向量的矩阵。

离散系统状态矩阵A和连续系统状态矩阵Ac存在如下关系：

系统频率ωi和阻尼比ζi与系统特征值ρi存在的关系见式(14)，由此可求出系统的第i阶频率ωi、阻尼比ζi和振型Vi，如式(15)所示。

式中：ρi和ρi*为互异复特征值。

1.3 基于AR-EWT 和SSI-Cov 的桥梁结构模态参数识别流程

图2 所示为基于AR-EWT 和SSI-Cov 的桥梁结构模态参数识别流程，其基本步骤如下。

图2 基于AR-EWT和SSI-Cov的桥梁结构模态参数识别流程Fig.2 Flowchart of modal identification for bridge structures using AR-EWT and SSI-Cov technique

1) 利用AR-EWT 分解振动响应信号。EWT 提供了Scalespace、Adaptivereg、Adaptive、Locmax、Locmaxmin 和Locmaxminf 共6 种频带划分算法[19]，其中，只有“Scalespace”方法不需给定初始条件即可完成频谱分割，但该方法计算效率低，计算时间长，不适合应用于数据量大的工程；其余5种方法中，“Adaptivereg”和“Adaptive”2种自适应法方法需人为给定初始边界，干扰因素多；“Locmax”、“Locmaxmin”和“Locmaxminf”这3种极值方法仅需人为给定分量数量即可，且这3种方法的原理和识别结果基本接近。因此，本文选取“Locmaxmin”算法划分频带边界，即在AR 谱上计算局部最大值，然后将边界设置为连续最大值之间的最小值；根据划分的频带边界值构造滤波器组，将信噪比低、模态密集的桥梁振动响应信号分解成一系列仅含单一模态信息的单分量信号。

2) 为了进一步降低噪声对后续模态参数识别的影响，对AR-EWT 分解后的单一模态分量信号进行奇异值(SVD)去噪处理[19]。

3) 对经去噪后的单分量信号进行Fourier变换，确认其所包含的模态信息，然后将包含同一模态信息的单分量信号组装成新的数据矩阵，利用SSICov对重组的各数据矩阵进行处理，识别得到结构的各阶模态频率、阻尼比和振型。

应该指出的是，基于AR-EWT和SSI-Cov的桥梁结构模态参数识别方法由于首先通过AR-EWT将结构响应信号进行了分解，再对分量信号进行重组，采用SSI识别的各组单分量信号数据时，信号中仅包含一阶模态信息，因此，在应用SSI算法时，可不考虑系统定阶问题，直接将系统阶次设置为2，降低了SSI 计算的复杂度，避免了人为确定系统阶次时的不可靠性，同时也避免了各阶模态之间的相互影响。

2 数值算例

2.1 ASCE Benchmark模型

为了验证结构健康监测领域中模态参数识别算法和损伤识别算法的有效性，美国土木工程师协会(ASCE)提供了一个基准模型，即ASCE Benchmark模型，如图3所示。该模型为1个4层2跨×2 跨的钢框架，层高为0.90 m，跨径为1.25 m。本文采用该钢框架的12 个自由度简化有限元模型[20]，该模型采用剪切梁假设，每层仅包括3个自由度，即X和Y这2个方向的平动自由度和扭转自由度。本文仅考虑模型Y方向的平动自由度，选取模型柱6(Column6)第1层至第4层的节点，每一节点仅考虑Y方向平动自由度，共4个自由度。

图3 ACSE Benchmark模型Fig.3 ACSE Benchmark model

在Y方向对模型各层施加高斯白噪声激励，计算得到模型各层的加速度响应，在响应中添加均方根为10%的白噪声。模态参数识别时采用模型柱6(Column6)各层Y方向包含噪声的加速度响应信号，信号采样频率为1 000 Hz，时长为40 s。采用AR-EWT 分解加速度时程响应， 选择“Locmaxmin”算法检测并划分频率边界，其中分解层数设置为6。图4 和图5 所示分别为模型柱6(Column6)第4 层Y方向的含噪加速度响应信号及其自回归功率谱和频带划分边界。从图5 可以看出，结构的前4 阶模态响应被分别划到第2、3、4和5 阶频带中，第1 和第6 阶频带用于剔除信号噪声。

图4 ACSE Benchmark模型柱6含噪加速度时程响应(第4层Y方向)Fig.4 Noisy acceleration response of ACSE Benchmark model column 6 (4th floor, Y direction)

图5 自回归功率谱及频带边界Fig.5 Autoregressive power spectrum and band boundary

基于划分的频带构建小波滤波器组，将信号分解为4个单分量信号，分别代表模型的4阶模态响应。图6 所示为基于AR-EWT 分解模型柱6(Column6)第4 层Y方向含噪加速度响应信号得到的4阶模态分量信号及其傅里叶幅值谱。同理，分别将模型柱6(Column6)第1 至第3 层Y方向含噪加速度响应信号进行分解，并采用SVD 对所有分量信号进行去噪，组装模型第1至第4层数据中包含相同频率成分的模态分量信号，分别采用SSI-Cov识别得到模型结构的4阶模态参数。由于各组组装数据中仅包含1 阶模态信息，因此，在应用SSICov算法时，无需估计模型阶次，直接将模型阶次设置为2，由此避免了SSI算法中的模型定阶问题，提高了计算效率，同时也避免了各阶模态之间的相互影响。

图6 ACSE Benchmark模型柱6含噪加速度时程响应信号分解结果(第4层Y方向)及其傅里叶幅值谱Fig.6 Decomposition results of noisy acceleration responses of ACSE Benchmark model column 6(4th floor, Y direction) and corresponding Fourier spectrum

基于提出的AR-EWT和SSI-Cov模态参数识别方法，识别得到ASCE Benchmark 模型Y方向的4阶频率和阻尼比，如表1 所示，相应的振型如图7所示。此外，表1 中给出了ASCE Benchmark 模型频率和阻尼比的理论值，以及基于SSI-Cov 和文献[21]基于希尔伯特-黄变换(HHT)方法识别得到的频率和阻尼比。为说明SVD去噪对AR-EWT+SSICov模态识别结果的影响，表1还给出了未对单分量信号进行去噪时识别得到的频率和阻尼比。

表1 ACSE Benchmark 模型频率和阻尼比识别结果Table 1 Identified natural frequencies and damping ratios of the ACSE Benchmark model

图7 识别的ACSE Benchmark模型振型Fig.7 Identified mode shapes of ACSE Benchmark model

不同方法识别得到的频率和阻尼比结果与理论计算值之间的相对误差见表2。从表1 和表2 可以看出：HHT 方法识别得到的频率和阻尼比结果误差较大；基于AR-EWT+SSI-Cov 方法进行模态参数识别时，不论是否采用SVD 去噪所有单分量信号，频率识别精度均与SSI-Cov算法识别精度相差不大；SVD 去噪所有单分量信号后，基于AREWT+SSI-Cov方法识别的阻尼比精度要略比另外2种方法的高。这说明本文提出的AR-EWT+SSICov模态参数识别方法能有效识别出结构的频率和阻尼比，且采用SVD 去噪可提高阻尼比的识别精度。因此，本文后续算例中采用提出的AR-EWT+SSI-Cov 方法时，都采用SVD 去噪所有单分量信号。

表2 ACSE Benchmark 模型频率和阻尼比识别相对误差Table 2 Identification error of ACSE Benchmark model frequency and damping ratio%

基于SSI-Cov 算法和本文方法识别的ASCEBenchmark 模型Y方向前4 阶弯曲振型、相应的理论振型及其二者之间的振型相关系数φ如图7 所示。从图7可以看出：本文方法识别的振型与SSICov方法识别的振型相符，与理论振型基本一致。

2.2 简支梁模型

为进一步验证本文提出的模态参数识别方法识别桥梁结构模态参数的有效性，考虑一简支梁模型，其物理参数为：桥长L=10 m，弹性模量E=2.06×1011N/m2，单位长度质量m=37.68 kg/m，采用Matlab建立有限元模型，如图8所示，简支梁划分为10 个单元，11 个节点；结构采用Rayleigh 阻尼系数，阻尼矩阵C=αM+βK，式中，α和β为比例系数，α=0.143，β=0.0001，M和K分别为质量和刚度矩阵。通过Matlab有限元分析得到简支梁前4阶模态频率和阻尼比的理论值，如表3所示。

图8 简支梁有限元模型Fig.8 Finite element model of the simply-supported beam

在简支梁节点8 上竖向施加白噪声激励，由Matlab 程序计算得到简支梁节点2 至节点10 的竖向加速度，采样频率为200 Hz，采样时间为30 s。图9 和图10 所示分别为白噪声激励时程和计算得到的节点5竖向加速度响应时程图。为考虑噪声对模态参数识别结果的影响，在图10 所示计算节点响应中添加了均方根为10%的白噪声。图11 所示为节点5加速度响应信号的自回归功率谱及相应的频带划分边界。从图11 可以看出：简支梁的前4阶模态被清楚地划分开，由此可有效地提取4组含单一模态信息的单分量信号。同理，对节点2～4、节点6～10 的加速度响应进行分解，提取相应的单分量信号；利用SVD 去噪所有单分量信号，重组所有节点中包含相同模态信息的单分量信号，最后，利用SSI-Cov算法得到简支梁模型的频率和阻尼比，如表3所示。表3还给出了基于SSI-COV方法识别的频率、阻尼比及其与理论值的相对误差。从表3可以看出：AR-EWT+SSI-Cov方法能有效识别简支梁模型的频率和阻尼比，且阻尼比识别精度比SSI-Cov方法的高。基于AR-EWT+SSI-Cov方法识别得到的简支梁模型振型及其与理论振型之间的振型相关系数见图12。从图12 可以看出：AR-EWT+SSI-Cov 方法识别振型与理论振型较吻合。

图9 白噪声激励时程Fig.9 Excitation of white noise

图10 计算得到的节点5的竖向加速度时程响应Fig.10 Caculated vertical acceleration time history response of Node 5

图11 自回归功率谱及频带边界Fig.11 Autoregressive power spectrum and band boundary

图12 简支梁模型振型识别结果Fig.12 Identified mode shapes of the simply-supported beam model

3 实桥应用：西宁北川河桥模态参数识别

为验证AR-EWT+SSI-Cov 方法识别实际桥梁结构模态参数的适用性，以西宁北川河桥通车前进行的动力试验结果为例，对测试结果进行处理，识别该桥的模态参数。西宁北川河桥是位于西宁至湟源一级公路中祁连路高架桥上的1 座连线桥，如图13 所示，为中承式钢管混凝土系杆拱桥，净跨为90 m。2002 年对该桥进行动力试验，在桥梁每根吊杆两侧以及桥面起始处分别布置32 个测点和1个参考点，测点布置位置如图14所示。

图13 西宁北川河桥Fig.13 Xining Beichuang River Bridge

图14 测点布置Fig.14 Layout of measuring points

试验共分为4 组，每组8 个测点和1 个共用参考点。测试时，用1辆重车在桥梁上来回行驶，从而激励桥梁振动，通过竖向加速度传感器记录桥梁的加速度信号，信号采样频率为80 Hz，采样时间为15 min。根据文献[22]中有限元分析结果，西宁北川河桥的前4阶竖向频率和前5阶扭转频率均在0～8 Hz 之间。为更好地计算桥梁自振频率并降低计算成本，采用截止频率为16 Hz的10阶巴特沃斯低通滤波器对结果进行滤波处理后再进行重采样处理，将数据采样频率降到16 Hz。

对33 组数据分别采用AR-EWT 分解，得到33×9 组含有不同频率成分的单分量信号。采用SVD 去噪后，重组具有相同频率成分的单分量信号，分别采用SSI-Cov 识别得到该桥的前9 阶自振频率、阻尼比和前4 阶桥面竖向振型，分别如表4和图15所示。表4还给出了文献[22]中基于峰值法和SSI算法识别得到的自振频率、阻尼比和有限元模型计算得到的理论频率。从表4可以看出：本文方法与其他2种方法的识别结果以及理论值结果基本一致。从图15 可以看出：本文方法识别得到的振型与文献[22]中的有限元模型计算振型结果较吻合。由此表明，本文提出的AR-EWT+SSI-Cov 方法能准确识别实际桥梁结构的频率、阻尼比和振型。

4 结论

1) 提出了一种结合AR-EWT和SSI-Cov的桥梁结构模态参数识别方法。该方法采用AR-EWT 将结构响应信号分解为多个仅含单一固有模态信息的单分量信号，组装包含同一模态信息的单分量信号后，利用SSI-Cov 算法分别处理各组分量信号，识别得到各单一模态参数。该方法在保证模态参数识别精度的前提下，有效避免了SSI算法的模型定阶问题，降低了不同模态之间相互干扰对识别结果的影响，提高了模态参数识别效率。

2) 基于AR-EWT和SSI-Cov的桥梁结构模态参数识别方法，准确识别了ASCE Benchmark 模型、简支梁数值模型的自振频率、阻尼比和振型，识别结果与理论结果基本一致，相对其他识别方法，本文方法识别的阻尼比精度更高；识别了西宁北川河桥的前9阶自振频率、阻尼比和前4阶竖向振型，识别结果与基于峰值法、SSI方法的识别结果以及有限元计算结果较吻合。