CN104964837A - 基于emd的结构刚度损伤监测方法及系统 - Google Patents

Abstract

本发明公开了结构刚度损伤的监测方法及系统，方法包括步骤：监测结构振动响应并获取结构加速度响应信号；对加速度响应信号进行EMD分解，确定各个时刻第一个IMF分量斜率，并计算结构不同位置的监测因子；利用约束条件剔除虚假监测因子，监测因子与结构刚度损伤程度呈线性关系；根据监测因子随时间的变化特征确定结构刚度损伤发生的时刻，比较结构不同位置的监测因子的分布确定损伤发生的位置；根据刚度损伤发生时刻对应的监测因子的幅值确定结构的刚度损伤程度。本发明特别适用于冲击荷载作用下的结构微小刚度损伤监测判定，能准确的识别微小程度的刚度损伤，具有较强的抗噪能力，通过计算监测因子的大小可以判定损伤程度，具有较高的应用价值。

Description

基于EMD的结构刚度损伤监测方法及系统

技术领域

本发明涉及一种基于EMD(经验模态分解算法)的结构刚度损伤监测方法。本发明还涉及一种专用于上述方法的基于EMD的结构刚度损伤监测系统。

背景技术

工程结构在外荷载和环境下的安全服役具有重要的意义。目前国内外正出现越来越多的各类大型高层高耸结构、大跨度桥梁和大跨度空间结构。由于恶劣自然环境和强烈外荷载的作用，工程结构不可避免的发生损伤累积并引起性能退化和损伤破坏。强烈外荷载、地震作用和爆炸冲击荷载等容易引起工程结构的性能退化，并进而导致结构发生突变的刚度损伤，如焊接断裂、构件屈曲、支撑破裂等。因此，各国都非常重视结构的安全监测和性能评估工作。采用合理的监测评估手段确保工程结构的服役安全是各国学者和工程技术人员面临的一个现实问题，具有重要的科学意思和实际工程意义。

目前结构损伤监测识别方法的基本思想是利用监测系统获取结构的动静力响应，通过各类监测方法提取出结构响应中的损伤信号，再通过对损伤信号进行分析评估，建立合理的监测评估方法和系统来评判结构的损伤事件。现有技术中损伤监测识别方法从信号来源角度可以主要分为基于频域信息的方法、基于时域信息的方法和基于时频域信息的方法三大类。

工程结构在服役过程中，容易遭受恶劣自然环境和地震、强风、爆炸等荷载的作用。因此容易引起构件由于损伤累积而发生的屈曲和失稳，这将进一步引起结构发生具有突然性的刚度损伤事故。国外多次发生由于结构失稳引发刚度损伤而发生损伤破坏的事故。目前对于这类由于失稳引起的突发性刚度损伤主要采用时频分析工具如小波变换和希尔伯特黄变换等进行损伤监测和识别。这类方法其基本思想均是将包含突变损伤信息的结构动力响应进行高通滤波，通过观察滤波后信号是否出现峰值来判断结构损伤事件。但目前这类方法在识别突发性的刚度损伤上均具有较为明显的弱点：

(1)通常只能识别较大程度的结构刚度损伤(如10％以上)，对于较小的结构刚度损伤(如5％以下)其识别效果很不稳定，取决于结构本身的特性和外荷载的特性。

(2)只能识别突发性刚度损伤发生的位置和发生的时间，却无法定量的确定突变损伤的程度。因此只能进行定性识别而不能进行定量的识别，这是目前方法的一个主要的缺陷。

(3)目前方法不能有效进行结构在冲击荷载作用下的微小刚度损伤的识别。由于损伤信号过于微弱，因此无法判别损伤发生的时间和位置，更无法判别其损伤程度。

(4)目前方法对噪声的抗干扰能力不强。在结构发生较小刚度损伤后，根本无法有效监测到损伤事件。

发明内容

本发明所要解决的第一个技术问题，就是提供一种基于EMD的结构刚度损伤的监测方法。

本发明所要解决的第二个技术问题，就是提供一种专用于上述方法的基于EMD的结构刚度损伤监测系统。

本发明的系统及方法，具有适用范围广、监测效率高、监测精度高、识别效果稳定和抗噪性好的优点。

解决上述第一个技术问题，本发明采用的技术方案是：

一种基于EMD的结构刚度损伤监测方法，其特征是包括如下步骤：

步骤S1，在结构不同位置上安装多个加速度传感器实时监测所述结构的动力响应，获取结构不同位置的加速度响应信号；

步骤S2，对所述结构不同位置的加速度响应信号进行EMD分解，并获取具有最高频信号成分的第一个IMF分量；

步骤S3，计算各个时刻IMF分量的斜率，并计算结构不同位置的一种用以判断结构刚度损伤的监测因子；

步骤S4，利用约束条件剔除虚假的监测因子，所述监测因子与结构的刚度损伤程度呈线性关系；

步骤S5，根据所述监测因子随时间的变化特征确定结构刚度损伤发生的时刻，通过比较结构不同位置的监测因子的分布确定损伤发生的位置；

步骤S6，根据刚度损伤发生时刻对应的监测因子的幅值确定结构的刚度损伤程度。

所述的步骤S2具体为：通过以下公式对所述加速度响应信号进行EMD分解，并确定各个不同的IMF分量：

$\stackrel{\·\·}{x}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{f}_{\mathrm{EMD}}^i\left(\stackrel{\·\·}{x}\left(t\right)\right)+r_n\left(t\right);$

其中：为加速度响应信号；为加速度信号经过EMD分解后的第i个IMF分量；r_n(t)为加速度信号经过EMD分解后的残余趋势项。

所述的步骤S3计算结构不同位置的监测因子通过以下公式进行：

MI_i＝|(D_i-D_i-1)+(D_i-D_i+1)|＝|2D_i-D_i-1-D_i+1| (i＝2,3,...,t_max-1)

其中，D_i表示i时刻的结构加速度响应第一个IMF分量的变化率；D_i-1和D_i+1表示i-1和i+1时刻的加速度响应第一个IMF分量的变化率；t_max为加速度响应信号的最大时间长度。

所述的步骤S4的监测因子约束条件为：

${\mathrm{MI}}_{i-1}+{\mathrm{MI}}_{i+1}\stackrel{\·}{=}{\mathrm{MI}}_i$

将监测因子与约束条件联合使用剔除虚假的监测因子。

所述的步骤S3计算结构不同位置的监测因子更具体的步骤如下：

S3-1，建立发生和不发生刚度损伤结构的频率之间的相互关系

若等效单自由度结构体系由于杆件屈曲或失稳，其结构的刚度在时刻t_i发生了损伤和减小，则结构刚度由正常值K₀减小了ΔK，变为K_s：

$K=\left\{\begin{array}{cc}{K}_0& (0\&le;t\&le;t_i)\\ K_s& (t_i<t)\end{array}\right.;$

ΔK＝K₀-K_s；

则无刚度损伤结构的频率f₀和有刚度损伤结构的频率f_s分别表示为：

$f_0=\frac{1}{2\mathrm{\&pi;}}\sqrt{\frac{K_0}{M}};$

$f_s=\frac{1}{2\mathrm{\&pi;}}\sqrt{\frac{K_0}{M}};$

式中，M为结构的质量；

结构的刚度变化用无刚度损伤结构和有刚度损伤结构的频率表示：

$\mathrm{\&Delta;K}=K_0-K_s=4{\mathrm{\&pi;}}^{2}M({f_0}^2-{f_s}^2);$

S3-2，建立无失稳损伤的原始结构的动力响应计算方法

没有发生失稳损伤的结构不存在刚度损伤，因此其等效单自由度体系的运动方程表示为：

$\stackrel{\&CenterDot;\&CenterDot;}{x}+4\mathrm{\&pi;\&xi;}{f}_0\stackrel{\&CenterDot;}{x}+4{\mathrm{\&pi;}}^2{f_0}^{2}x=0;$

式中，ξ为结构体系的阻尼比；

结构在外荷载引起的脉冲作用下将发生具有初始速度为v₀的振动，则结构的位移、速度和加速度响应分别计算为：

$x\left(t\right)=\frac{v_0\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)}{2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d}\&CenterDot;e^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}};$

$\stackrel{\&CenterDot;}{x}\left(t\right)=v_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}(\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)-\frac{\mathrm{\&xi;}}{{\mathrm{\&xi;}}_d}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right));$

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=-2\mathrm{\&pi;}{f}_0{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}\&CenterDot;\frac{[\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)(2{\mathrm{\&xi;}}_d^2-1)+2\mathrm{\&xi;}{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)]}{{\mathrm{\&xi;}}_d};$

式中：

${\mathrm{\&xi;}}_d=\sqrt{1-{\mathrm{\&xi;}}^2};$

S3-3，建立发生失稳及刚度损伤的结构的动力响应计算方法

在结构在外荷载作用下发生振动中，假设在t_i时刻发生了构件失稳，由于构件失稳过程具有突然性，因此发生过程很短，这将导致构件的刚度在很短的时间内发生减小；为描述该刚度变化过程，采用一个新的时间坐标轴t₁＝t-t_i来描述发生失稳的结构的振动状况；因此，有刚度损伤结构的运动方程表示为：

${\stackrel{\&CenterDot;\&CenterDot;}{x}}_s+4\mathrm{\&pi;\&xi;}{f}_s{\stackrel{\&CenterDot;}{x}}_s+4{\mathrm{\&pi;}}^2{f_s}^2{x}_s=0,(t>t_i);$

则有损伤结构振动的初始条件由无损结构在时刻t_i的位移和速度响应确定：

$x_s\left(0\right)=x\left(t_i\right)=\frac{\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)}{2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d}\&CenterDot;v_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i};$

${\stackrel{\&CenterDot;}{x}}_s\left(0\right)=\stackrel{\&CenterDot;}{x}\left(t_i\right)=v_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}(\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-\frac{\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)}{{\mathrm{\&xi;}}_d});$

由此计算得到在t₁时刻有损伤结构的加速度响应为：

${\stackrel{\&CenterDot;\&CenterDot;}{x}}_s\left(t_1\right)=-\frac{2\mathrm{\&pi;}{v}_0{f}_s(E_1+E_2+E_3)}{{\mathrm{\&xi;}}_d^2{f}_0}\&CenterDot;e^{-2\mathrm{\&pi;\&xi;}(f_s{t}_1+f_0{t}_i)};$

式中：

E₁＝f_s sin(2πf₀ξ_dt_i)[ξ_d cos(2πf_sξ_dt₁)-ξsin(2πf_sξ_dt₁)]；

$E_2=f_0\sin \left(2\mathrm{\&pi;}{f}_s{\mathrm{\&xi;}}_d{t}_1\right)(1-2{\mathrm{\&xi;}}_d^2)\left(\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right){-\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)\right);$

$E_3=-2f_0\mathrm{\&xi;}\cos \left(2\mathrm{\&pi;}{f}_s{\mathrm{\&xi;}}_d{t}_1\right)[\mathrm{\&xi;}{\mathrm{\&xi;}}_d\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d^2\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)];$

由于工程结构的阻尼比往往很小，因此构件发生失稳破坏而导致的刚度损伤发生的时间很短，则有：

t₁＝t_i+1-t_i＝Δt≈0；

sin(2πf_sξ_dt₁)≈0；

cos(2πf_sξ_dt₁)≈1；

由此得结构在发生刚度损伤后后t_i+1的时刻的加速度响应为：

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t_{i+1}\right)={\stackrel{\&CenterDot;\&CenterDot;}{x}}_s\left(\mathrm{\&Delta;t}\right)=-\frac{2\mathrm{\&pi;}{f}_s{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}}{f_0{\mathrm{\&xi;}}_d}\{f_s\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-2f_0\mathrm{\&xi;}[\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)\left]\right\};$

S3-4，进行结构振动响应的EMD分解

为了建立损伤监测因子，则需对结构的加速度响应进行经验模态分解(Empirical Mode Decomposition，EMD)；

具体处理方法是：首先，确定加速度响应的多个局部极大值和局部极小值：采用多次样条函数将的局部极大值点与局部极小值点分别拟和得到其峰值的上包络曲线与下包络曲线然后计算两包络线的均值m₁(t)

$m_1\left(t\right)=\frac{{\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)}_{\mathrm{en\; velope}}^u+{\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)}_{\mathrm{en\; velope}}^l}{2};$

将原加速度序列减去该平均包络m₁(t)后即得一个去掉低频的新加速度时程序列h₁(t)：

$h_1\left(t\right)=\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)-m_1\left(t\right);$

对得到的h₁(t)重复以上数据过程，重复k次直至所得到的平均包络趋于零为止：

h_1k(t)＝h_1(k-1)(t)-m_1k(t)；

其中：h_1k(t)为第k次处理所得加速度数据；h_1(k-1)(t)为第k-1次处理所得加速度数据；m_1k(t)为h_1(k-1)(t)上下包络线的均值；

由此得到该加速度响应的第一个内敛模函数分量(intrinsic modefunction,IMF)c₁(t)；

c₁(t)＝h_1k(t)；

第一个IMF分量c₁(t)代表了原始加速度信号中的最高频成分；将原始加速度响应减去第一个IMF分量c₁(t)，得去掉高频成分的加速度响应时程r₁(t)；将r₁(t)再作为要分解的信号重复上述过程，直至所剩余信号r₁(t)已是一单调函数时停止此分解过程；此时的参与量r_n(t)代表原始加速度响应的低频趋势项；由此确定加速度响应的一组IMF分量c₁(t),c₂(t)…c_n(t)；原始的加速度响应由全部IMF分量和一个趋势项的叠加来表示：

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i\left(t\right)+r_n\left(t\right);$

将结构振动信号的EMD分解过程采用一个隐函数来表示；

则无损伤结构加速度相应的EMD分解信号表示为：

$\begin{array}{c}\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i\left(t\right)+r_n\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i(\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)+r_n\left(t\right))\\ =\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i(-2\mathrm{\&pi;}{f}_0{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}\&CenterDot;\frac{[\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)(2{\mathrm{\&xi;}}_d^2-1)+2\mathrm{\&xi;}{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)]}{{\mathrm{\&xi;}}_d})+r_n\left(t\right)\end{array}$

其中：第i个IMF分量c_i(t)表示为：

$c_i\left(t\right)=f_{\mathrm{EMD}}^i(-2\mathrm{\&pi;}{f}_0{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}\&CenterDot;\frac{[\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)(2{\mathrm{\&xi;}}_d^2-1)+2\mathrm{\&xi;}{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)]}{{\mathrm{\&xi;}}_d});$

同理，发生刚度损伤后结构的加速度响应表示为：

$\begin{array}{c}{\stackrel{\&CenterDot;\&CenterDot;}{x}}_s\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i^s\left(t\right)+r_n^s\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i\left(\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)\right)+r_n^s\left(t\right)\\ =\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i(-\frac{2\mathrm{\&pi;}{f}_0{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}}{f_0{\mathrm{\&xi;}}_d}\{f_s\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-2f_0\mathrm{\&xi;}[\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)]\left\}\right)+r_n^s\left(t\right)\end{array}$

；其中：第i个IMF分量表示为：

$c_i^s\left(t\right)=f_{\mathrm{EMD}}^i(-\frac{2\mathrm{\&pi;}{f}_s{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}}{f_0{\mathrm{\&xi;}}_d}\{f_s\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-2f_0\mathrm{\&xi;}[\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)]\left\}\right);$

S3-5，确定损伤监测因子

EMD分解基于加速度响应局部特征时间尺度，从原加速度时程中提取固有模态函数，其本质是将加速度信号中不同频率和尺度的波动或趋势逐级分解开来；所分解出的各IMF分量分别包含了原加速度信号的不同时间尺度和频率特征的局部特征信息；

突变损伤时刻的信号不连续具有两个明显的特点：(1)信号的幅值在损伤时刻t_i到时刻t_i+1发生了很大的跳跃；(2)在时刻t_i-1(或其它任何之前时刻)和时刻t_i+1(或其它任何之后时刻)信号的斜率远小于损伤时刻t_i的信号斜率；事实上研究表明，发生刚度突然损伤时的结构加速度响应的第一个IMF分量也具有上述的两个相同的特点。

显然可知，发生瞬时失稳事故时，结构构件刚度突然减小，结构的加速度响应出现了一个突然的跳跃；这个突然的跳跃信号具有明显的高频特征和大振幅特点；由于突变刚度损伤具有高频特性，因此，其加速度响应的突变信号只保留在具有最高频成分的第一个IMF分量中；

损伤前后的结构加速度响应第一个IMF的变化率D_i表示为：

$\begin{array}{c}{D}_i=\frac{\mathrm{\&Delta;}{c}_1}{\mathrm{\&Delta;t}}=\frac{c_{1,i+1}-c_{1,i}}{\mathrm{\&Delta;t}}=f_{\mathrm{EMD}}^1({\stackrel{\&CenterDot;\&CenterDot;}{x}}_{i+1}-{\stackrel{\&CenterDot;\&CenterDot;}{x}}_i)\\ \stackrel{\&CenterDot;}{=}{f}_{\mathrm{EMD}}^1(-\frac{2\mathrm{\&pi;}{\mathrm{\&xi;}}_d{v}_0({f_s}^2-{f_0}^2)}{f_0}\&CenterDot;e^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}\&CenterDot;\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right))\end{array},(i=\mathrm{1,2},...,t_{\max }-1);$

式中：Δt为加速度响应信号的时刻间距，t_max为加速度响应信号的最大时间长度；由于：

${f_s}^2-{f_0}^2=\frac{\mathrm{\&Delta;K}}{4{\mathrm{\&pi;}}^{2}M};$

则损伤前后的结构加速度响应第一个IMF的变化率D_i表示为：

$D_i\stackrel{\&CenterDot;}{=}{f}_{\mathrm{EMD}}^1(-\frac{\mathrm{\&Delta;K}{\mathrm{\&xi;}}_d{v}_0}{2\mathrm{\&pi;}{f}_{0}M}\&CenterDot;e^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}\&CenterDot;\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)),(i=\mathrm{1,2},...,t_{\max }-1);$

由于EMD分解过程是一个线性过程，原始加速度信号表示为所有IMF分量和残余分量的线性叠加，因此存在如下关系：

|D_i|∝|ΔK|；

前述的加速度响应信号不连续的第二个特点在数学表示为：

$\left\{\begin{array}{c}\left|D_i\right|>>\left|D_{i-j}\right|\\ \left|D_i\right|>>\left|D_{i+j}\right|\end{array}\right.,(j=\mathrm{1,2},...,t_{\max });$

由此，得到一种结构由于瞬时失稳事故引起的刚度损伤的监测因子(Monitoring Index)MI_i：

MI_i＝|(D_i-D_i-1)+(D_i-D_i+1)|＝|2D_i-D_i-1-D_i+1| (i＝2,3,...,t_max-1)；

由于存在如下关系：

$\left\{\begin{array}{c}\left|D_i\right|>>\left|D_{i-1}\right|\\ \left|D_i\right|>>\left|D_{i+1}\right|\end{array}\right.,(i=\mathrm{2,3},...,t_{\max }-1);$

因此有：

${\mathrm{MI}}_i=|2D_i-D_{i-1}-D_{i+1}|\stackrel{\&CenterDot;}{=}2\left|D_i\right|;$

由上述推导过程可知，所采用的基于EMD的监测因子与结构刚度损伤的程度呈正比关系，即：

MI_i∝|ΔK|；

基于该损伤指标，对应于时刻t_i-1和t_i+1的损伤指标MI_i-1和MI_i+1表示为：

MI_i-1＝|2D_i-1-D_i-2-D_i|；

MI_i+1＝|2D_i+1-D_i-D_i+2|；

考虑i时刻前后的第一个IMF分量斜率均小于D_i,则时刻t_i-1和t_i+1的损伤指标MI_i-1和MI_i+1之和近似等于损伤时刻t_i,的损伤指标值MI_i，即得出监测因子MI_i的约束条件：

${\mathrm{MI}}_{i-1}+{\mathrm{MI}}_{i+1}\stackrel{\&CenterDot;}{=}{\mathrm{MI}}_i;$

S3-6，确定损伤监测因子与损伤程度的对应关系

建立线性模型：

MI＝α*S+β；

式中：MI为损伤指标的幅值；S为损伤程度大小；α、β均为线性模型中的常数参数；

实际应用过程中首先确定待监测结构的基本信息如：质量、刚度、阻尼比等；在预先设定的荷载作用下，采用数值分析方法或模型试验模拟结构发生多种不同程度的损伤，然后可以建立损伤程度与损伤指标的之间的数值关系，并采用数值回归的方法确定参数α和β的数值；由此可以得到监测因子与损伤程度之间的定量关系。

解决上述第二个技术问题，本发明采用的技术方案是：

一种专用于上述方法的基于EMD的结构刚度损伤监测系统，其特征是：包括依次连接的下面五个模块：

振动监测模块，用于实时监测所述结构的振动状况，获取所述结构各层的加速度响应信号；

基于EMD的监测信号处理模块，用于对所述加速度响应信号进行EMD分解、并将所述加速度响应信号的第一个IMF分量提出出来；利用监测因子和约束条件，确定监测因子的具体数值并剔除虚假的监测因子；

刚度损伤时间判别模块，用于根据监测因子时程曲线，并联合采用约束条件，判别刚度损伤事件发生的时间；

刚度损伤位置判别模块，用于根据结构不同位置的监测因子的空间部分，并联合采用约束条件，判别刚度损伤事件发生的位置；

刚度损伤程度判别模块，用于根据损伤发生时刻对应的监测因子的幅值确定结构的损伤程度。

本发明通过在建筑结构上安装多个加速度传感器实时监测结构的振动响应。首先获取结构不同位置的加速度响应信号，然后对加速度响应信号进行EMD分解，对分解得到的第一个IMF分量进行分析，提出一种用以判断结构刚度损伤的监测因子。

若有结构由于杆件屈曲失稳发生刚度损伤，则可通过对振动信号进行EMD分解，然后计算第一个IMF分量的监测因子，判断刚度损伤发生的时间、位置以及损伤程度。本发明的监测方法解决了传统方法识别精度不高、不能适用于冲击荷载作用条件、不能识别微小刚度损伤、抗噪性差和无法定量判别损伤程度的缺陷。

有益效果：本发明能准确的识别微小程度的刚度损伤，具有较强的抗噪能力，而且本发明的监测因子与刚度损伤程度呈线性关系，通过计算监测因子的大小可以判定损伤程度，克服了传统监测方法不能识别结构损伤程度和难以判别冲击荷载作用下的微小刚度损伤的不足，具有较高的应用价值。

本发明适用于各种不同类型荷载作用下的结构由于构件屈曲失稳引起的刚度损伤监测判别，特别适用于冲击荷载作用下的结构微小刚度损伤监测判定。

附图说明

图1为本发明的基于EMD的结构刚度损伤监测方法的流程示意图；

图2为发生刚度突然损伤时的结构加速度响应的变化模式；

图3为多层建筑架结构的立体示意图；

图4a为本发明案例一中结构第一层的加速度响应曲线图；

图4b为图4a的前三个IMF分量曲线图之一；

图4c为图4a的前三个IMF分量曲线图之二；

图4d为图4a的前三个IMF分量曲线图之三；

图5a为本发明案例一中结构不发生不同程度刚度损伤时，加速度响应的第一个IMF分量曲线之一；

图5b为本发明案例一中结构不发生不同程度刚度损伤时，加速度响应的第一个IMF分量曲线之二；

图5c为本发明案例一中结构不发生不同程度刚度损伤时，加速度响应的第一个IMF分量曲线之三；

图5d为本发明案例一中结构不发生不同程度刚度损伤时，加速度响应的第一个IMF分量曲线之四；

图5e为本发明案例一中结构不发生不同程度刚度损伤时，加速度响应的第一个IMF分量曲线之五；

图5f为本发明案例一中结构不发生不同程度刚度损伤时，加速度响应的第一个IMF分量曲线之六；

图6a为本发明案例一中结构第一楼层的监测因子结果；

图6b为本发明案例一中结构第二楼层的监测因子结果；

图6c为本发明案例一中结构第三楼层的监测因子结果；

图6d为本发明案例一中结构第四楼层的监测因子结果；

图6e为本发明案例一中结构第五楼层的监测因子结果；

图7a为本发明案例一中不同损伤程度下监测因子结果之一；

图7b为本发明案例一中不同损伤程度下监测因子结果之二；

图7c为本发明案例一中不同损伤程度下监测因子结果之三；

图7d为本发明案例一中不同损伤程度下监测因子结果之四；

图7e为本发明案例一中不同损伤程度下监测因子结果之五；

图7f为本发明案例一中不同损伤程度下监测因子结果之六；

图8a为本发明案例一中小损伤程度下联合利用监测因子和约束条件的监测结果之一；

图8b为本发明案例一中小损伤程度下联合利用监测因子和约束条件的监测结果之二；

图9为本发明案例一中损伤程度与监测因子之间的相互关系；

图10a为本发明案例一中不同噪声水平下的监测因子结果之一；

图10b为本发明案例一中不同噪声水平下的监测因子结果之二；

图11a为本发明案例二中在地震激励作用下的结构第一层的加速度响应曲线；

图11b为图11a的前三个IMF分量曲线之一；

图11c为图11a的前三个IMF分量曲线之二；

图11d为图11a的前三个IMF分量曲线之三；

图12a为本发明案例二中结构第一楼层的监测因子结果；

图12b为本发明案例二中结构第二楼层的监测因子结果；

图12c为本发明案例二中结构第三楼层的监测因子结果；

图12d为本发明案例二中结构第四楼层的监测因子结果；

图12e为本发明案例二中结构第五楼层的监测因子结果；

图13a为本发明案例二中不同损伤程度下监测因子结果之一；

图13b为本发明案例二中不同损伤程度下监测因子结果之二；

图13c为本发明案例二中不同损伤程度下监测因子结果之三；

图13d为本发明案例二中不同损伤程度下监测因子结果之四；

图13e为本发明案例二中不同损伤程度下监测因子结果之五；

图13f为本发明案例二中不同损伤程度下监测因子结果之六；

图14为本发明案例二中损伤程度与监测因子之间的相互关系；

图15a为本发明案例二中不同噪声水平下的监测因子结果之一；

图15b为本发明案例二中不同噪声水平下的监测因子结果之二；

图16a为本发明案例三中在地震激励作用下，在小损伤程度下采用db2和db4小波的监测结果之一；

图16b为本发明案例三中在地震激励作用下，在小损伤程度下采用db2和db4小波的监测结果之二；

图16c为本发明案例三中在地震激励作用下，在小损伤程度下采用db2和db4小波的监测结果之三；

图16d为本发明案例三中在地震激励作用下，在小损伤程度下采用db2和db4小波的监测结果之四；

图16e为本发明案例三中在地震激励作用下，在小损伤程度下采用db2和db4小波的监测结果之五；

图16f为本发明案例三中在地震激励作用下，在小损伤程度下采用db2和db4小波的监测结果之六；

图17a为本发明案例三中在冲击荷载作用下，不同损伤程度下采用db4小波的监测结果之一；

图17b为本发明案例三中在冲击荷载作用下，不同损伤程度下采用db4小波的监测结果之二；

图17c为本发明案例三中在冲击荷载作用下，不同损伤程度下采用db4小波的监测结果之三；

图17d为本发明案例三中在冲击荷载作用下，不同损伤程度下采用db4小波的监测结果之四；

图17e为本发明案例三中在冲击荷载作用下，不同损伤程度下采用db4小波的监测结果之五；

图17f为本发明案例三中在冲击荷载作用下，不同损伤程度下采用db4小波的监测结果之六；

图18为本发明中的基于EMD的结构刚度损伤的监测系统的结构示意图。

具体实施方式

下面结合附图和具体实施例对本发明内容作进一步的详细描述。

如图1所示为本发明的基于EMD的结构刚度损伤监测方法实施例，其具体流程如下：

S1在结构不同位置上安装多个加速度传感器实时监测所述结构的动力响应，监测结构的振动响应，并获取架结构不同位置的加速度响应信号；

S2对所述结构不同位置的加速度响应信号进行EMD分解，确定各个时刻第一个IMF分量的斜率；

S3计算各个时刻IMF分量的斜率，并计算结构不同位置的监测因子；

S4利用约束条件剔除虚假的监测因子，所述监测因子与结构的刚度损伤程度呈线性关系；

S5根据所述监测因子随时间的变化特征确定结构刚度损伤发生的时刻，通过比较结构不同位置的监测因子的分布确定损伤发生的位置；

S6根据刚度损伤发生时刻对应的监测因子的幅值确定结构的刚度损伤程度。

在一种具体实施方式中，可通过以下公式对所述加速度响应信号进行EMD分解，并确定各个不同的IMF分量：

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i\left(\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)\right)+r_n\left(t\right);$

其中：为加速度响应信号；加速度信号经过EMD分解后的第i个IMF分量；r_n(t)为加速度信号经过EMD分解后的残余趋势项。

通过上述公式确定了加速度响应的IMF分量后，本实施例还提供了一种监测因子计算方法，具体通过以下公式进行计算：

MI_i＝|(D_i-D_i-1)+(D_i-D_i+1)|＝|2D_i-D_i-1-D_i+1|(i＝2,3,...,t_max-1)；

其中，D_i表示i时刻的结构加速度响应第一个IMF分量的变化率；D_i-1和D_i+1表示i-1和i+1时刻的加速度响应第一个IMF分量的变化率；t_max为加速度响应信号的最大时间长度。

同时还提出了一种监测因子的约束条件：

${\mathrm{MI}}_{i-1}+{\mathrm{MI}}_{i+1}\stackrel{\&CenterDot;}{=}{\mathrm{MI}}_i;$

从以上公式可知：监测因子与结构加速度响应第一个IMF分量的变化率密切相关的。

该监测因子还与结构刚度损伤程度呈对应关系，下面从理论上给出证明。

具体而言通过以下步骤建立结构由于失稳导致的刚度损伤的监测因子：

步骤一：建立发生和不发生刚度损伤结构的频率之间的相互关系

若等效单自由度结构体系由于杆件屈曲或失稳，其结构的刚度在时刻t_i发生了损伤和减小，则结构刚度由正常值K₀减小了ΔK，变为K_s：

$K=\left\{\begin{array}{cc}{K}_0& (0\&le;t\&le;t_i)\\ K_s& (t_i<t)\end{array}\right.$

ΔK＝K₀-K_s

则无刚度损伤结构的频率f₀和有刚度损伤结构的频率f_s则分别表示为：

$f_0=\frac{1}{2\mathrm{\&pi;}}\sqrt{\frac{K_0}{M}}$

$f_s=\frac{1}{2\mathrm{\&pi;}}\sqrt{\frac{K_0}{M}}$

式中，M为结构的质量；结构的刚度变化用无刚度损伤结构和有刚度损伤结构的频率表示：

ΔK＝K₀-K_s＝4π²M(f₀ ²-f_s ²)；

步骤二：建立无失稳损伤的原始结构的动力响应计算方法

没有发生失稳损伤的结构不存在刚度损伤，因此其等效单自由度体系的运动方程可表示为：

$\stackrel{\&CenterDot;\&CenterDot;}{x}+4\mathrm{\&pi;\&xi;}{f}_0\stackrel{\&CenterDot;}{x}+4{\mathrm{\&pi;}}^2{f_0}^{2}x=0$

式中，ξ为结构体系的阻尼比；

结构在外荷载引起的脉冲作用下将发生具有初始速度为v₀的振动，则结构的位移、速度和加速度响应可分别计算为：

$x\left(t\right)=\frac{v_0\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)}{2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d}\&CenterDot;e^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}$

$\stackrel{\&CenterDot;}{x}\left(t\right)=v_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}(\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)-\frac{\mathrm{\&xi;}}{{\mathrm{\&xi;}}_d}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right))$

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=-2\mathrm{\&pi;}{f}_0{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}\&CenterDot;\frac{[\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)(2{\mathrm{\&xi;}}_d^2-1)+2\mathrm{\&xi;}{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)]}{{\mathrm{\&xi;}}_d}$

式中：

${\mathrm{\&xi;}}_d=\sqrt{1-{\mathrm{\&xi;}}^2};$

步骤三：建立发生失稳及刚度损伤的结构的动力响应计算方法

在结构在外荷载作用下发生振动中，假设在t_i时刻发生了构件失稳，由于构件失稳过程具有突然性，因此发生过程很短，这将导致构件的刚度在很短的时间内发生减小；为描述该刚度变化过程，采用一个新的时间坐标轴t₁＝t-t_i来描述发生失稳的结构的振动状况；因此，有刚度损伤结构的运动方程表示为：

${\stackrel{\&CenterDot;\&CenterDot;}{x}}_s+4\mathrm{\&pi;\&xi;}{f}_s{\stackrel{\&CenterDot;}{x}}_s+4{\mathrm{\&pi;}}^2{f_s}^2{x}_s=0,(t>t_i);$

则有损伤结构振动的初始条件由无损结构在时刻t_i的位移和速度响应确定：

$x_s\left(0\right)=x\left(t_i\right)=\frac{\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)}{2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d}\&CenterDot;v_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i};$

${\stackrel{\&CenterDot;}{x}}_s\left(0\right)=\stackrel{\&CenterDot;}{x}\left(t_i\right)=v_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}(\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-\frac{\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)}{{\mathrm{\&xi;}}_d});$

由此计算得到在t₁时刻有损伤结构的加速度响应为：

${\stackrel{\&CenterDot;\&CenterDot;}{x}}_s\left(t_1\right)=-\frac{2\mathrm{\&pi;}{v}_0{f}_s(E_1+E_2+E_3)}{{\mathrm{\&xi;}}_d^2{f}_0}\&CenterDot;e^{-2\mathrm{\&pi;\&xi;}(f_s{t}_1+f_0{t}_i)};$

式中：

E₁＝f_s sin(2πf₀ξ_dt_i)[ξ_d cos(2πf_sξ_dt₁)-ξsin(2πf_sξ_dt₁)]；

$E_2=f_0\sin \left(2\mathrm{\&pi;}{f}_s{\mathrm{\&xi;}}_d{t}_1\right)(1-2{\mathrm{\&xi;}}_d^2)\left(\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right){-\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)\right);$

$E_3=-2f_0\mathrm{\&xi;}\cos \left(2\mathrm{\&pi;}{f}_s{\mathrm{\&xi;}}_d{t}_1\right)[\mathrm{\&xi;}{\mathrm{\&xi;}}_d\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d^2\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)];$

由于工程结构的阻尼比往往很小，因此构件发生失稳破坏而导致的刚度损伤发生的时间很短，则有：

t₁＝t_i+1-t_i＝Δt≈0；

sin(2πf_sξ_dt₁)≈0；

cos(2πf_sξ_dt₁)≈1；

由此得结构在发生刚度损伤后后t_i+1的时刻的加速度响应为：

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t_{i+1}\right)={\stackrel{\&CenterDot;\&CenterDot;}{x}}_s\left(\mathrm{\&Delta;t}\right)=-\frac{2\mathrm{\&pi;}{f}_s{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}}{f_0{\mathrm{\&xi;}}_d}\{f_s\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-2f_0\mathrm{\&xi;}[\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)\left]\right\};$

步骤四：进行结构振动响应的EMD分解

为了建立损伤监测因子，则需对结构的加速度响应进行经验模态分解(Empirical Mode Decomposition，EMD)；

具体处理方法是：首先，确定加速度响应的多个局部极大值和局部极小值：采用多次样条函数将的局部极大值点与局部极小值点分别拟和得到其峰值的上包络曲线与下包络曲线然后计算两包络线的均值m₁(t)

$m_1\left(t\right)=\frac{{\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)}_{\mathrm{en\; velope}}^u+{\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)}_{\mathrm{en\; velope}}^l}{2};$

将原加速度序列减去该平均包络m₁(t)后即得一个去掉低频的新加速度时程序列h₁(t)：

$h_1\left(t\right)=\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)-m_1\left(t\right);$

对得到的h₁(t)重复以上数据过程，重复k次直至所得到的平均包络趋于零为止：

h_1k(t)＝h_1(k-1)(t)-m_1k(t)；

其中：h_1k(t)为第k次处理所得加速度数据；h_1(k-1)(t)为第k-1次处理所得加速度数据；m_1k(t)为h_1(k-1)(t)上下包络线的均值；

由此得到该加速度响应的第一个内敛模函数分量(intrinsic modefunction,IMF)c₁(t)；

c₁(t)＝h_1k(t)；

第一个IMF分量c₁(t)代表了原始加速度信号中的最高频成分；将原始加速度响应减去第一个IMF分量c₁(t)，得去掉高频成分的加速度响应时程r₁(t)；将r₁(t)再作为要分解的信号重复上述过程，直至所剩余信号r₁(t)已是一单调函数时停止此分解过程，此时的参与量r_n(t)代表原始加速度响应的低频趋势项；由此确定加速度响应的一组IMF分量c₁(t),c₂(t)…c_n(t)；原始的加速度响应由全部IMF分量和一个趋势项的叠加来表示：

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i\left(t\right)+r_n\left(t\right);$

将结构振动信号的EMD分解过程采用一个隐函数来表示；

则无损伤结构加速度相应的EMD分解信号表示为：

$\begin{array}{c}\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i\left(t\right)+r_n\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i(\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)+r_n\left(t\right))\\ =\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i(-2\mathrm{\&pi;}{f}_0{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}\&CenterDot;\frac{[\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)(2{\mathrm{\&xi;}}_d^2-1)+2\mathrm{\&xi;}{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)]}{{\mathrm{\&xi;}}_d})+r_n\left(t\right)\end{array}$

其中：第i个IMF分量c_i(t)表示为：

$c_i\left(t\right)=f_{\mathrm{EMD}}^i(-2\mathrm{\&pi;}{f}_0{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}\&CenterDot;\frac{[\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)(2{\mathrm{\&xi;}}_d^2-1)+2\mathrm{\&xi;}{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)]}{{\mathrm{\&xi;}}_d});$

同理，发生刚度损伤后结构的加速度响应表示为：

$\begin{array}{c}{\stackrel{\&CenterDot;\&CenterDot;}{x}}_s\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i^s\left(t\right)+r_n^s\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i\left(\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)\right)+r_n^s\left(t\right)\\ =\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i(-\frac{2\mathrm{\&pi;}{f}_0{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}}{f_0{\mathrm{\&xi;}}_d}\{f_s\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-2f_0\mathrm{\&xi;}[\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)]\left\}\right)+r_n^s\left(t\right)\end{array}$

；其中：第i个IMF分量表示为：

$c_i^s\left(t\right)=f_{\mathrm{EMD}}^i(-\frac{2\mathrm{\&pi;}{f}_s{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}}{f_0{\mathrm{\&xi;}}_d}\{f_s\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-2f_0\mathrm{\&xi;}[\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)]\left\}\right);$

步骤五：确定损伤监测因子

EMD分解基于加速度响应局部特征时间尺度，从原加速度时程中提取固有模态函数，其本质是将加速度信号中不同频率和尺度的波动或趋势逐级分解开来；所分解出的各IMF分量分别包含了原加速度信号的不同时间尺度和频率特征的局部特征信息；

图2给出了发生刚度突然损伤时的结构加速度响应的变化模式。由图2中结果可知，突变损伤时刻的信号不连续具有两个明显的特点：(1)信号的幅值在损伤时刻t_i到时刻t_i+1发生了很大的跳跃；(2)在时刻t_i-1(或其它任何之前时刻)和时刻t_i+1(或其它任何之后时刻)信号的斜率远小于损伤时刻t_i的信号斜率。事实上研究表明，发生刚度突然损伤时的结构加速度响应的第一个IMF分量也具有上述的两个相同的特点。

显然可知，发生瞬时失稳事故时，结构构件刚度突然减小，结构的加速度响应出现了一个突然的跳跃。这个突然的跳跃信号具有明显的高频特征和大振幅特点。由于突变刚度损伤具有高频特性，因此，其加速度响应的突变信号只保留在具有最高频成分的第一个IMF分量中。因此，损伤前后的结构加速度响应第一个IMF的变化率D_i可表示为：

$\begin{array}{c}{D}_i=\frac{\mathrm{\&Delta;}{c}_1}{\mathrm{\&Delta;t}}=\frac{c_{1,i+1}-c_{1,i}}{\mathrm{\&Delta;t}}=f_{\mathrm{EMD}}^1({\stackrel{\&CenterDot;\&CenterDot;}{x}}_{i+1}-{\stackrel{\&CenterDot;\&CenterDot;}{x}}_i)\\ \stackrel{\&CenterDot;}{=}{f}_{\mathrm{EMD}}^1(-\frac{2\mathrm{\&pi;}{\mathrm{\&xi;}}_d{v}_0({f_s}^2-{f_0}^2)}{f_0}\&CenterDot;e^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}\&CenterDot;\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right))\end{array},(i=\mathrm{1,2},...,t_{\max }-1);$

式中：Δt为加速度响应信号的时刻间距，t_max为加速度响应信号的最大时间长度；由于：

${f_s}^2-{f_0}^2=\frac{\mathrm{\&Delta;K}}{4{\mathrm{\&pi;}}^{2}M};$

则损伤前后的结构加速度响应第一个IMF的变化率D_i表示为：

$D_i\stackrel{\&CenterDot;}{=}{f}_{\mathrm{EMD}}^1(-\frac{\mathrm{\&Delta;K}{\mathrm{\&xi;}}_d{v}_0}{2\mathrm{\&pi;}{f}_{0}M}\&CenterDot;e^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}\&CenterDot;\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)),(i=\mathrm{1,2},...,t_{\max }-1);$

由于EMD分解过程是一个线性过程，原始加速度信号表示为所有IMF分量和残余分量的线性叠加，因此存在如下关系：

|D_i|∝|ΔK|；

前述的加速度响应信号不连续的第二个特点在数学表示为：

$\left\{\begin{array}{c}\left|D_i\right|>>\left|D_{i-j}\right|\\ \left|D_i\right|>>\left|D_{i+j}\right|\end{array}\right.,(j=\mathrm{1,2},...,t_{\max });$

由此，得到本实施例中的一种结构由于瞬时失稳事故引起的刚度损伤的监测因子(Monitoring Index)MI_i：

MI_i＝|(D_i-D_i-1)+(D_i-D_i+1)|＝|2D_i-D_i-1-D_i+1| (i＝2,3,...,t_max-1)；

由于存在如下关系

$\left\{\begin{array}{c}\left|D_i\right|>>\left|D_{i-1}\right|\\ \left|D_i\right|>>\left|D_{i+1}\right|\end{array}\right.,(i=\mathrm{2,3},...,t_{\max }-1);$

因此有：

${\mathrm{MI}}_i=|2D_i-D_{i-1}-D_{i+1}|\stackrel{\&CenterDot;}{=}2\left|D_i\right|;$

显然由上述推导过程可知，本实施例中所采用的基于EMD的监测因子与结构刚度损伤的程度呈正比关系，即：

MI_i∝|ΔK|；

基于该损伤指标，对应于时刻t_i-1和t_i+1的损伤指标MI_i-1和MI_i+1表示为：

MI_i-1＝|2D_i-1-D_i-2-D_i|；

MI_i+1＝|2D_i+1-D_i-D_i+2|；

考虑i时刻前后的第一个IMF分量斜率均小于D_i,则时刻t_i-1和t_i+1的损伤指标MI_i-1和MI_i+1之和近似等于损伤时刻t_i,的损伤指标值MI_i，即得出监测因子MI_i的约束条件：

${\mathrm{MI}}_{i-1}+{\mathrm{MI}}_{i+1}\stackrel{\&CenterDot;}{=}{\mathrm{MI}}_i;$

由分析结果可知，一个典型判别瞬时刚度损伤引起的加速度信号不连续的监测因子模式包括一个在损伤时刻t_i的相对较大的MI_i和两个在时刻t_i-1和t_i+1的相对较小的MI_i-1和MI_i+1,它们符合上述方程的近似关系。基于监测因子和内在约束条件的联合运用，可以在识别过程中剔除一些虚假的损伤状况，这在微小损伤的识别过程中非常具有价值。

步骤六：确定损伤监测因子与损伤程度的对应关系

由前述分析和推导可知本实施所提出的监测因子与结构刚度损伤程度间存在对应关系。在实际应用过程中，可建立刚度损伤程度与监测因子的定量线性关系。针对结构不同的刚度损伤变化程度，分析结构在荷载作用下的损伤监测因子，在此基础上可建立线性模型：

MI＝α*S+β；

式中：MI为损伤指标的幅值；S为损伤程度大小；α、β均为线性模型中的常数参数。

实际应用过程中可以首先确定待监测结构的基本信息如：质量、刚度、阻尼比等。在预先设定的荷载作用下，可以采用数值分析方法或模型试验模拟结构发生多种不同程度的损伤，然后可以建立损伤程度与损伤指标的之间的数值关系，并可以采用数值回归的方法确定参数α和β的数值。由此可以得到监测因子与损伤程度之间的定量关系。

在实际应用中可在结构不同部位安装加速度传感器，实时监测结构的加速度响应信号。若结构在服役过程中由于遭受外荷载发生了构件失稳，则其引起的刚度损伤会导致加速度信号发生突变。整个过程的加速度响应将被实现安装于结构上的加速度传感器所监测到。监测系统将自动将所得的加速度信号进行EMD分解并提取出第一个IMF分量。

在此基础上，计算各个时刻的监测因子MI_i的结果。通过对比监测因子在时间上变化，可以判定刚度损伤发生的时刻。通过对比监测因子在结构上的空间分布可以确定刚度损伤发生的位置。最后可将监测因子入上述的损伤程度模型中，即可确定结构发生的刚度损伤程度：

$S=\frac{\mathrm{MI}-\mathrm{\&beta;}}{\mathrm{\&alpha;}}$

通过上述一系列步骤，即可实现结构由于构件失稳所引发的刚度损伤的智能监测和评估。

综上所述，本发明与现有技术相比具有以下优点：

一、本发明提出的监测因子具有物理概念清晰、分析快速有效的的优点。该监测因子和监测系统具有很好的适用性，适用于各种不同荷载作用下的结构由于失稳所引发的刚度损伤识别。

二、传统的刚度损伤监测识别方法在不具有很好的抗噪性，噪声的强度和频谱范围对识别效果的影响很大。在小损伤和强噪声作用下，往往不能监测出刚度损伤事件。而本发明提出的监测因子很好的抗噪性，在微小损伤和较高的噪声水平下仍具有很好的监测效果。

三、现有的结构刚度损伤的监测方法通常只能判别较大程度损伤事件，对于非常小的结构刚度损伤，如2％以下的刚度损伤识别效果不佳。而本发明提出的监测因子对微小损伤敏感，可以准确的监测出微笑刚度损伤事件，克服了目前方法的一个明显不足。

四、现有的结构刚度损伤的监测方法通常不能监测到由于冲击荷载所引起的微小损伤。因为此时加速度响应中的高频损伤信号能量微弱，难以用常规方法进行识别。而本发明提出的监测指标和约束条件联合使用，可以利用前后三个时刻的幅值关系，剔除识别出的虚假信息，实现冲击荷载作用下的微小刚度损伤的准确识别，克服了目前监测方法的一个明显不足。

五、现有的刚度损伤监测方法不能够定量的判别损伤程度。特别是在冲击荷载作用下的小损伤识别时，更是难以确定结构的损伤程度。而本发明采用的监测因子明确了刚度损伤程度与监测因子之间的定量对应关系，确立了两者的数学定量模型，这克服了现有方法的一个不足。

下面以几个具体案例来说明本发明的结构刚度损伤监测方法及系统的有效性。

案例一

本案例一考察在冲击荷载作用下的结构刚度损伤的监测识别效果。图3为一个五层建筑结构的示意图，以该结构为对象，描述刚度损伤监测的实施过程。该建筑结构层高2.8m，每层有10个房间，各层的质量为m＝1.35×10⁶kg，刚度为k＝4.2×10⁶kN/m。该建筑结构第一层在冲击荷载作用下由于构件破坏发生了突发的刚度损伤，导致第一节点层的刚度在0.2秒时刻发生了20％的突然减小，即第一节点层刚度由4.2×10⁶kN/m减小为3.36×10⁶kN/m。

图4a-图4d显示了该结构第一层的加速度时程曲线及其前三个IMF曲线。从图4a-图4d中显示的加速度时程曲线中并不能直接发现在0.2s发生了损伤事件。但是对比不同的IMF曲线可知，在IMF1分量的0.2s信号出现了一个不连续。而IMF2和IMF3分量重均未发现信号不连续。事实上其它IMF分量和参与趋势项中也没有发现由于刚度损伤所引起的加速度信号不连续。而由EMD分解过程可知：第一个IMF1分量包含了原始加速度信号中的最高频成分。由此验证了前述所指出的结论：加速度信号中的突变成分保留在具有最高频信号成分的第一个IMF分量中。

表一

为了进一步说明本发明所提出的监测因子和系统的有效性。可使该结构发生6种不同程度的刚度损伤，刚度损伤在0.2秒发生，发生于结构的第一层。6种损伤程度分别为1％，2％，5％，10％，20％和40％的刚度损伤。表一给出了无刚度损伤结构和6种不同损伤工况下前5阶自振频率的比较。显然表中结果表明：在结构较小刚度损伤下，结构的自振频率改变非常小。结构频率的变化率甚至还低于常规的信号噪声的强度(约2％～5％)。因此，采用常规的基于振动的健康监测方法是无法识别结构的小程度刚度损伤事件。

图5a-图5f给出了不同损伤程度下的IMF1分量中的信号不连续状况。

(不同信号分解完毕有不同数量的IMF，可能有5个，可能有10几个，这取决于信号本身，但是对所有情况，只有第一个IMF包含有损伤信息，所以本文中只给出了第一个IMF分量)

由图中结果可知，在不同程度损伤下，在加速度信号的IMF分量的0.2s时刻均可观察到由于刚度损伤所引起的信号不连续。但是在损伤程度较小时(5％以下)，这种信号不连续已经非常微弱以至于难以用肉眼直观察觉。在损伤为1％时，仅仅从IMF分量本身根本无法判别发生了损伤事件。

图6a-图6e给出了采用本发明所提出的监测因子的监测结果。由图6a-图6e结果可知，第一层的监测因子在0.2秒时出现了一个非常明显的峰值，而其他时刻则没有峰值。这表明第一层在0.2s时发生了刚度损伤事故，因此可以判定损伤发生的时刻。进一步的对比不同楼层的监测因子的结果可知：只有第一层的监测因子存在峰值，而其他楼层没有峰值。因此通过对比监测因子在结构不同了楼层的分布，可以容易的确定刚度损伤发生的位置。

图7a-图7f给出了不同刚度损伤程度下的监测因子的识别结果。由图中结果可知在损伤程度5％以上时可以准确的识别出刚度损伤发生的时间，对比不同楼层的识别结果可以准确的监测到损伤发生的位置。但是在刚度损伤很小的时候(1％～2％)，并不能准确的识别损伤事件。这是因为在冲击荷载作用下，结构动力响应很快衰减并且具有很宽的频谱范围。同时由于损伤程度很小，因此损伤信号的能量很小，在原始加速度信号中很微弱，因此造成了损伤识别的困难。这是目前其他方法难以识别冲击荷载作用下的微小刚度损伤的原因。进一步的，联合利用本发明建立的监测因子和约束条件，可以成功地识别结构的微小损伤，如图8a-图8b所示。

图9示出了冲击荷载作用下损伤程度与监测因子之间的相互关系。图中结果表明，损伤程度与监测因子是明显的线性关系，这也符合前述的理论推导结果。图中还给出了通过最小二乘回归得到了损伤程度的数值模型。利用该模型可以识别结构的刚度损伤程度。

图10a–图10b显示了噪声对于损伤识别的影响(此时损伤程度为20％)。信号噪声通常认为是随机白噪声，噪声强度为2％和5％。由图9所示的结果可知，本发明所提出的监测因子在噪声干扰下依然能够有效地识别刚度损伤事件。

案例二

本案例二以前述建筑结构为例，考察在地震激励作用下的监测方法的有效性。本案例二中刚度损伤选用的六种损伤程度均与案例一相同，但损伤发生的时间为6.0s，地震地面运动选用El Centro地震波。

图11a-图11d给出了地震作用下的结构第一层的加速度时程曲线及其前三个IMF曲线。与冲击荷载作用下的情况类似，直接从加速度时程曲线中并不能直接发现损伤事件。但是具有最高频信号成分的第一个IMF分量中可以观察到由于刚度损伤引起的信号不连续。

图12a-图12e给出了不同楼层监测因子的结果。显然可以通过监测因子确定损伤发生的时间为6.0s，同时由监测因子的空间分布可以容易的确定刚度损伤发生的位置为第一层。

图13a-图13f给出了不同刚度损伤程度下的监测因子的识别结果。显然，在各种不同损伤程度下，监测因子均可有效的识别损伤事件。与冲击荷载作用下的情况不同，地震激励作用下结构动力响应较大，因此即使在微小损伤下，加速度响应中的损伤信号依然具有一定的能量，因此直接采用监测因子即可识别微小损伤时间，而不需要联合使用约束条件。这一点是与冲击荷载作用下的情况不同的。

图14显示了损伤程度与监测因子的相互关系。与冲击荷载作用下的情况相同，损伤程度与监测因子符合线性模型，利用该模型可以识别结构的刚度损伤程度。图15a-图15b结果则表明，本发明所提出的监测因子在噪声干扰下依然能够有效地识别地震作用下的结构刚度损伤事件。

案例三

本案例三对比了基于小波的监测识别方法与本发明的识别方法的效果比较。所采用的结构和损伤工况均与案例一相同。基于小波分析的监测方法主要是选择合适的小波基函数对原始加速度信号进行多层小波分解，通过观察高频小波系数的时程曲线来识别刚度损伤事件。在目前小波刚度损伤识别中db小波一类较为广泛采用的小波基函数。图16a-图16f给出了地震激励作用下的具有典型意义的db2和db4小波的分析结果。db4小波相比于db2小波具有更好的支撑长度，因此具有更加的识别分辨率。因此db4小波基函数在刚度损伤是别中被广泛采用。由图15a-图15b结果可知，在较小的刚度损伤下，db4小波均可有效的识别损伤事件，而db2小波的刚度损伤则为失败。

图17a-图17f给出了在冲击荷载作用下采用db4小波的识别结果。显然，在大损伤程度下，db4小波可以有效地识别刚度损伤事件。但是小损伤情况下(5％以下)，特别是微小损伤下(1％)，db4小波已经无法有效的识别损伤事件，这一点与本发明提出的监测因子具有明显的差异。

本发明还提供了一种基于EMD的结构刚度损伤监测系统，如图18所示，包括：

振动监测模块1，通过在结构上安装加速度传感器，实时监测结构的振动状况，获取结构各层的振动响应号。

基于EMD的监测信号处理模块2，用于对所述加速度响应进行EMD分解，并将各信号的第一个IMF分量提出出来。利用本发明提出的监测因子和约束条件，确定监测因子的具体数值。

刚度损伤时间判别模块3，根据监测因子时程曲线，并联合采用约束条件，判别刚度损伤事件发生的时间；

刚度损伤位置判别模块4，根据结构不同位置的监测因子的空间部分，并联合采用约束条件，判别刚度损伤事件发生的位置；

刚度损伤程度判别模块5，用于根据损伤发生时刻对应的监测因子的幅值确定结构的损伤程度。

在一种具体实施方式中，模块2通过以下公式对所述加速度响应信号进行EMD分解，并确定各个IMF分量：

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i\left(t\right)+r_n\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i\left(\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)\right)+r_n\left(t\right)$

其中：为加速度响应信号；加速度信号经过EMD分解后的第i个IMF分量；r_n(t)为加速度信号经过EMD分解后的残余趋势项。

在一种具体实施方式中，信号处理模块2根据如下公式计算i时刻的监测因子和约束条件：

MI_i＝|(D_i-D_i-1)+(D_i-D_i+1)|＝|2D_i-D_i-1-D_i+1| (i＝2,3,...,t_max-1)

${\mathrm{MI}}_{i-1}+{\mathrm{MI}}_{i+1}\stackrel{\&CenterDot;}{=}{\mathrm{MI}}_i$

其中，D_i表示i时刻的结构加速度响应第一个IMF分量的变化率；D_i-1和D_i+1表示i-1和i+1时刻的加速度响应第一个IMF分量的变化率。

在一种具体实施方式中，如图18所示。本发明的基于EMD的结构刚度损伤监测系统在实施过程中可以在预设荷载作用下使结构发生不同程度的损伤，通过模块2确定对应的监测因子，进一步的建立损伤程度与监测因子之间的数值模型关系，并通过统计回归的方法确定监测因子与损伤程度之间的模型参数α和β的数值，从而确定监测因子与损伤程度之间的定量模型；损伤程度判别模块5根据损伤发生时刻对应的监测因子大小，通过定量模型确定结构的刚度损伤程度。

本发明的基于EMD的结构刚度损伤监测系统中各个模块其具体功能的实现可采用上述的方法。

振动监测模块，通过在结构上安装加速度传感器，实时监测结构的振动状况，获取结构各层的振动响应号。

监测信号处理模块通过以下公式对所述加速度响应信号进行EMD分解并确定各个IMF分量：

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i\left(t\right)+r_n\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i\left(\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)\right)+r_n\left(t\right);$

其中：为加速度响应信号；加速度信号经过EMD分解后的第i个IMF分量；r_n(t)为加速度信号经过EMD分解后的残余趋势项；

计算模块根据如下公式计算i时刻的监测因子和约束条件：

MI_i＝|(D_i-D_i-1)+(D_i-D_i+1)|＝|2D_i-D_i-1-D_i+1| (i＝2,3,...,t_max-1)；

${\mathrm{MI}}_{i-1}+{\mathrm{MI}}_{i+1}\stackrel{\&CenterDot;}{=}{\mathrm{MI}}_i;$

其中，D_i表示i时刻的结构加速度响应第一个IMF分量的变化率；D_i-1和D_i+1表示i-1和i+1时刻的加速度响应第一个IMF分量的变化率。t_max为加速度响应信号的最大时间长度；

刚度损伤时间判别模块、刚度损伤位置判别模块和刚度损伤程度判别模块，在预设荷载作用下使结构发生不同程度的损伤，通过监测信号处理模块确定对应的监测因子，进一步的建立损伤程度与监测因子之间的数值模型关系，并通过统计回归的方法确定监测因子与损伤程度之间的模型参数α和β的数值，从而确定监测因子与损伤程度之间的定量模型；损伤程度判别模块根据损伤发生时刻对应的监测因子大小，通过定量模型确定结构的刚度损伤程度。

Claims (6)

1.一种基于EMD的结构刚度损伤监测方法，其特征是包括如下步骤：

步骤S1，在结构不同位置上安装多个加速度传感器实时监测所述结构的动力响应，获取结构不同位置的加速度响应信号；

步骤S2，对所述结构不同位置的加速度响应信号进行EMD分解，并获取具有最高频信号成分的第一个IMF分量；

步骤S3，计算各个时刻IMF分量的斜率，并计算结构不同位置的一种用以判断结构刚度损伤的监测因子；

步骤S4，利用约束条件剔除虚假的监测因子，所述监测因子与所述结构的刚度损伤程度呈线性关系；

步骤S5，根据所述监测因子随时间的变化特征确定结构刚度损伤发生的时刻，通过比较结构不同位置的监测因子的分布确定损伤发生的位置；

步骤S6，根据刚度损伤发生时刻对应的监测因子的幅值确定结构的刚度损伤程度。

2.根据权利要求1所述的基于EMD的结构刚度损伤监测方法，其特征是：所述的步骤S2具体为：通过以下公式对所述结构不同位置的加速度响应信号进行EMD分解，并确定各个不同的IMF分量：

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i\left(\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)\right)r_n\left(t\right);$

其中：为加速度响应信号；为加速度信号经过EMD分解后的第i个IMF分量；r_n(t)为加速度信号经过EMD分解后的残余趋势项。

3.根据权利要求1所述的基于EMD的结构刚度损伤监测方法，其特征是：所述的步骤S3计算结构不同位置的监测因子具体通过以下公式进行：

MI_i＝|(D_i-D_i-1)+(D_i-D_i+1)|＝|2D_i-D_i-1-D_i+1| (i＝2,3,...,t_max-1)；

其中，D_i表示i时刻的结构加速度响应第一个IMF分量的变化率；D_i-1和D_i+1表示i-1和i+1时刻的加速度响应第一个IMF分量的变化率；t_max为加速度响应信号的最大时间长度。

4.根据权利要求3所述的基于EMD的结构刚度损伤监测方法，其特征是：所述的步骤S3计算结构不同位置的监测因子的步骤如下：

S3-1，建立发生和不发生刚度损伤结构的频率之间的相互关系

若等效单自由度结构体系由于杆件屈曲或失稳，其结构的刚度在时刻t_i发生了损伤和减小，则结构刚度由正常值K₀减小了ΔK，变为K_s：

$K=\left\{\begin{array}{cc}{K}_0& (0\&le;t\&le;t_i)\\ K_s& (t_i<t)\end{array}\right.;$

ΔK＝K₀-K_s；

则无刚度损伤结构的频率f₀和有刚度损伤结构的频率f_s分别表示为：

$f_0=\frac{1}{2\mathrm{\&pi;}}\sqrt{\frac{K_0}{M}};$

$f_s=\frac{1}{2\mathrm{\&pi;}}\sqrt{\frac{K_0}{M}};$

式中，M为结构的质量；

结构的刚度变化用无刚度损伤结构和有刚度损伤结构的频率表示：

ΔK＝K₀-K_s＝4π²M(f₀ ²-f_s ²)；

S3-2，建立无失稳损伤的原始结构的动力响应计算方法

没有发生失稳损伤的结构不存在刚度损伤，因此其等效单自由度体系的运动方程表示为：

$\stackrel{\&CenterDot;\&CenterDot;}{x}+4\mathrm{\&pi;\&xi;}{f}_0\stackrel{\&CenterDot;}{x}+4{\mathrm{\&pi;}}^2{f}_0^{2}x=0;$

式中，ξ为结构体系的阻尼比；

结构在外荷载引起的脉冲作用下将发生具有初始速度为v₀的振动，则结构的位移、速度和加速度响应分别为：

$x\left(t\right)=\frac{v_0\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)}{2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d}\&CenterDot;e^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}};$

$\stackrel{\&CenterDot;}{x}\left(t\right)=v_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}(\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)-\frac{\mathrm{\&xi;}}{{\mathrm{\&xi;}}_d}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right));$

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=-2\mathrm{\&pi;}{f}_0{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}\&CenterDot;\frac{[\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)(2{\mathrm{\&xi;}}_d^2-1)+2{\mathrm{\&xi;\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)]}{{\mathrm{\&xi;}}_d};$

式中：

${\mathrm{\&xi;}}_d=\sqrt{1-{\mathrm{\&xi;}}^2};$

S3-3，建立发生失稳及刚度损伤的结构的动力响应计算方法

在结构在外荷载作用下发生振动中，假设在t_i时刻发生了构件失稳，由于构件失稳过程具有突然性，因此发生过程很短，这将导致构件的刚度在很短的时间内发生减小；为描述该刚度变化过程，采用一个新的时间坐标轴t₁＝t-t_i来描述发生失稳的结构的振动状况；因此，有刚度损伤结构的运动方程表示为：

${\stackrel{\&CenterDot;\&CenterDot;}{x}}_s+4\mathrm{\&pi;\&xi;}{f}_s{\stackrel{\&CenterDot;}{x}}_s+{4\mathrm{\&pi;}}^2{f_s}^2{x}_s=0,(t>t_i);$

则，有损伤结构振动的初始条件由无损结构在时刻t_i的位移和速度响应确定：

$x_s\left(0\right)=x\left(t_i\right)=\frac{\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)}{2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d}\&CenterDot;v_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i};$

${\stackrel{\&CenterDot;}{x}}_s\left(0\right)=\stackrel{\&CenterDot;}{x}\left(t_i\right)=v_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}(\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-\frac{\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)}{{\mathrm{\&xi;}}_d});$

由此计算得到在t₁时刻有损伤结构的加速度响应为：

${\stackrel{\&CenterDot;\&CenterDot;}{x}}_s\left(t_1\right)=-\frac{2\mathrm{\&pi;}{v}_0{f}_s(E_1+E_2+E_3)}{{\mathrm{\&xi;}}_d^2{f}_0}\&CenterDot;e^{-2\mathrm{\&pi;\&xi;}(f_s{t}_1+f_0{t}_i);}$

式中：

E₁＝f_ssin(2πf₀ξ_dt_i)[ξ_dcos(2πf_sξ_dt₁)-ξsin(2πf_sξ_dt₁)]；

$E_2=f_0\sin \left(2\mathrm{\&pi;}{f}_s{\mathrm{\&xi;}}_d{t}_1\right)(1-2{\mathrm{\&xi;}}_d^2)(\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right));$

$E_3=-2f_0\mathrm{\&xi;}\cos \left(2\mathrm{\&pi;}{f}_s{\mathrm{\&xi;}}_d{t}_1\right)[\mathrm{\&xi;}{\mathrm{\&xi;}}_d\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d^2\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)];$

由于工程结构的阻尼比往往很小，因此构件发生失稳破坏而导致的刚度损伤发生的时间很短，则有：

t₁＝t_i+1-t_i＝Δt≈0；

sin(2πf_sξ_dt₁)≈0；

cos(2πf_sξ_dt₁)≈1；

由此得结构在发生刚度损伤后后t_i+1的时刻的加速度响应为：

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t_{i+1}\right)={\stackrel{\&CenterDot;\&CenterDot;}{x}}_s\left(\mathrm{\&Delta;t}\right)=-\frac{2\mathrm{\&pi;}{f}_s{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}}{f_0{\mathrm{\&xi;}}_d}\{f_s\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-2f_0\mathrm{\&xi;}[\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)\left]\right\};$

S3-4，进行结构振动响应的EMD分解

为了建立损伤监测因子，则需对结构的加速度响应进行经验模态分解；

具体处理方法是：首先，确定加速度响应的多个局部极大值和局部极小值：采用多次样条函数将的局部极大值点与局部极小值点分别拟和得到其峰值的上包络曲线与下包络曲线然后计算两包络线的均值m₁(t)：

$m_1\left(t\right)=\frac{\stackrel{\&CenterDot;\&CenterDot;}{x}{\left(t\right)}_{\mathrm{envelope}}^u+\stackrel{\&CenterDot;\&CenterDot;}{x}{\left(t\right)}_{\mathrm{envelope}}^l}{2};$

将原加速度序列减去该平均包络m₁(t)后即得一个去掉低频的新加速度时程序列h₁(t)：

$h_1\left(t\right)=\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)-m_1\left(t\right);$

对得到的h₁(t)重复以上数据过程，重复k次直至所得到的平均包络趋于零为止：

h_1k(t)＝h_1(k-1)(t)-m_1k(t)；

其中：h_1k(t)为第k次处理所得加速度数据；h_1(k-1)(t)为第k-1次处理所得加速度数据；m_1k(t)为h_1(k-1)(t)上下包络线的均值；

由此得到该加速度响应的第一个内敛模函数分量c₁(t)：

c₁(t)＝h_1k(t)；

第一个IMF分量c₁(t)代表了原始加速度信号中的最高频成分；

将原始加速度响应减去第一个IMF分量c₁(t)，得去掉高频成分的加速度响应时程r₁(t)；将r₁(t)再作为要分解的信号重复上述过程，直至所剩余信号r₁(t)已是一单调函数时停止此分解过程；

此时的参与量r_n(t)代表原始加速度响应的低频趋势项；由此确定加速度响应的一组IMF分量c₁(t),c₂(t)…c_n(t)；原始的加速度响应由全部IMF分量和一个趋势项的叠加来表示：

$\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i\left(t\right)+r_n\left(t\right);$

将结构振动信号的EMD分解过程采用一个隐函数来表示，则无损伤结构加速度相应的EMD分解信号表示为：

$\begin{array}{c}\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)\text{=}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i\left(t\right)+r_n\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i\left(\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)\right)+r_n\left(t\right)\\ =\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i(-2\mathrm{\&pi;}{f}_0{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}\&CenterDot;\frac{[\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)(2{\mathrm{\&xi;}}_d^2-1)+2\mathrm{\&xi;}{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)]}{{\mathrm{\&xi;}}_d})+r_n\left(t\right)\end{array}$

其中：第i个IMF分量c_i(t)表示为：

$c_i\left(t\right)=f_{\mathrm{EMD}}^i(-2\mathrm{\&pi;}{f}_0{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;t}}\&CenterDot;\frac{[\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)(2{\mathrm{\&xi;}}_d^2-1)+2{\mathrm{\&xi;\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_{d}t\right)]}{{\mathrm{\&xi;}}_d});$

同理，发生刚度损伤后结构的加速度响应表示为：

$\begin{array}{c}{\stackrel{\&CenterDot;\&CenterDot;}{x}}_s\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i^s\left(t\right)+r_n^s\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i\left(\stackrel{\&CenterDot;\&CenterDot;}{x}\left(t\right)\right)+r_n^s\left(t\right)\\ =\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{f}_{\mathrm{EMD}}^i(-\frac{2\mathrm{\&pi;}{f}_s{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}}{f_0{\mathrm{\&xi;}}_d}\{f_s\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-2f_0\mathrm{\&xi;}[\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)]\left\}\right)+r_n^s\left(t\right)\end{array}$

；其中：第i个IMF分量表示为：

$c_i^s\left(t\right)=f_{\mathrm{EMD}}^i(-\frac{2\mathrm{\&pi;}{f}_s{v}_0{e}^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}}{f_0{\mathrm{\&xi;}}_d}\{f_s\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-2f_0\mathrm{\&xi;}[\mathrm{\&xi;}\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)-{\mathrm{\&xi;}}_d\cos \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)]\left\}\right);$

S3-5，确定损伤监测因子

EMD分解基于加速度响应局部特征时间尺度，从原加速度时程中提取固有模态函数，其本质是将加速度信号中不同频率和尺度的波动或趋势逐级分解开来；所分解出的各IMF分量分别包含了原加速度信号的不同时间尺度和频率特征的局部特征信息；

突变损伤时刻的信号不连续具有两个明显的特点：(1)信号的幅值在损伤时刻t_i到时刻t_i+1发生了很大的跳跃；(2)在时刻t_i-1和时刻t_i+1信号的斜率远小于损伤时刻t_i的信号斜率；事实上研究表明，发生刚度突然损伤时的结构加速度响应的第一个IMF分量也具有上述的两个相同的特点；

显然可知，发生瞬时失稳事故时，结构构件刚度突然减小，结构的加速度响应出现了一个突然的跳跃；这个突然的跳跃信号具有明显的高频特征和大振幅特点；由于突变刚度损伤具有高频特性，因此，其加速度响应的突变信号只保留在具有最高频成分的第一个IMF分量中；

损伤前后的结构加速度响应第一个IMF的变化率D_i表示为：

$\begin{array}{c}{D}_i=\frac{{\mathrm{\&Delta;c}}_1}{\mathrm{\&Delta;t}}=\frac{c_{1,i+1}-c_{1,i}}{\mathrm{\&Delta;t}}=f_{\mathrm{EMD}}^1({\stackrel{\&CenterDot;\&CenterDot;}{x}}_{i+1}-{\stackrel{\&CenterDot;\&CenterDot;}{x}}_i)\\ \stackrel{}{\stackrel{\&CenterDot;}{=}}{f}_{\mathrm{EMD}}^1(-\frac{2\mathrm{\&pi;}{\mathrm{\&xi;}}_d{v}_0({f_s}^2-{f_0}^2)}{f_0}\&CenterDot;e^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}\&CenterDot;\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right))(i=\mathrm{1,2},...,t_{\max }-1)\end{array};$

式中：Δt为加速度响应信号的时刻间距，t_max为加速度响应信号的最大时间长度；由于：

${f_s}^2-{f_0}^2=\frac{\mathrm{\&Delta;K}}{{4\mathrm{\&pi;}}^{2}M};$

则损伤前后的结构加速度响应第一个IMF的变化率D_i表示为：

$D_i\stackrel{\&CenterDot;}{=}{f}_{\mathrm{EMD}}^1(-\frac{\mathrm{\&Delta;}{\mathrm{K\&xi;}}_d{v}_0}{2\mathrm{\&pi;}{f}_{0}M}\&CenterDot;e^{-2\mathrm{\&pi;}{f}_0\mathrm{\&xi;}{t}_i}\&CenterDot;\sin \left(2\mathrm{\&pi;}{f}_0{\mathrm{\&xi;}}_d{t}_i\right)),(i=\mathrm{1,2},...,t_{\max }-1);$

由于EMD分解过程是一个线性过程，原始加速度信号表示为所有IMF分量和残余分量的线性叠加，因此存在如下关系：

|D_i|∝|ΔK|；

前述的加速度响应信号不连续的第二个特点在数学表示为：

$\left\{\begin{array}{c}\left|D_i\right|>>\left|D_{i-j}\right|\\ \left|D_i\right|>>\left|D_{i+j}\right|\end{array}\right.,(j=\mathrm{1,2},...,t_{\max });$

由此，得到一种结构由于瞬时失稳事故引起的刚度损伤的监测因子(Monitoring Index)MI_i：

MI_i＝|(D_i-D_i-1)+(D_i-D_i+1)|＝|2D_i-D_i-1-D_i+1| (i＝2,3,...,t_max-1)；

由于存在如下关系：

$\left\{\begin{array}{c}\left|D_i\right|>>\left|D_{i-1}\right|\\ \left|D_i\right|>>\left|D_{i+1}\right|\end{array}\right.,(j=\mathrm{2,3},...,t_{\max }-1);$

因此有：

${\mathrm{MI}}_i=|{2D}_i-D_{i-1}-D_{i+1}|\stackrel{\&CenterDot;}{=}2\left|D_i\right|;$

由上述推导过程可知，所采用的基于EMD的监测因子与结构刚度损伤的程度呈正比关系，即：

MI_i∝|ΔK|；

基于该损伤指标，对应于时刻t_i-1和t_i+1的损伤指标MI_i-1和MI_i+1表示为：

MI_i-1＝|2D_i-1-D_i-2-D_i|；

MI_i+1＝|2D_i+1-D_i-D_i+2|；

考虑i时刻前后的第一个IMF分量斜率均小于D_i,则时刻t_i-1和t_i+1的损伤指标MI_i-1和MI_i+1之和近似等于损伤时刻t_i,的损伤指标值MI_i，即得出监测因子MI_i的约束条件：

${\mathrm{MI}}_{i-1}+{\mathrm{MI}}_{i+1}\stackrel{\&CenterDot;}{=}{\mathrm{MI}}_i;$

S3-6，确定损伤监测因子与损伤程度的对应关系

建立线性模型：

MI＝α*S+β；

式中：MI为损伤指标的幅值；S为损伤程度大小；α、β均为线性模型中的常数参数；

实际应用过程中首先确定待监测结构的基本信息：质量、刚度和阻尼比；在预先设定的荷载作用下，采用数值分析方法或模型试验模拟结构发生多种不同程度的损伤，然后建立损伤程度与损伤指标的之间的数值关系，并采用数值回归的方法确定参数α和β的数值；由此得到监测因子与损伤程度之间的定量关系。

5.根据权利要求1所述的基于EMD的结构刚度损伤监测方法，其特征是：所述的步骤S4的监测因子约束条件为：

${\mathrm{MI}}_{i-1}+{\mathrm{MI}}_{i+1}\stackrel{\&CenterDot;}{=}{\mathrm{MI}}_i;$

将监测因子与约束条件联合使用剔除虚假的监测因子。

6.一种专用于如权利要求1至5任意一项所述方法的基于EMD的结构刚度损伤监测系统，其特征是：包括依次连接的下面五个模块：

振动监测模块，用于实时监测所述结构的振动状况，获取所述结构各层的加速度响应信号；

基于EMD的监测信号处理模块，用于对所述加速度响应信号进行EMD分解、并将所述加速度响应信号的第一个IMF分量提出出来；利用监测因子和约束条件，确定监测因子的具体数值并剔除虚假的监测因子；

刚度损伤时间判别模块，用于根据监测因子时程曲线，并联合采用约束条件，判别刚度损伤事件发生的时间；

刚度损伤位置判别模块，用于根据结构不同位置的监测因子的空间部分，并联合采用约束条件，判别刚度损伤事件发生的位置；

刚度损伤程度判别模块，用于根据损伤发生时刻对应的监测因子的幅值确定结构的损伤程度。