[发明专利]基于EMD的结构刚度损伤监测方法及系统

权利要求书

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

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

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

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

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

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

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

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

$\stackrel{\·\·}{x}\left(t\right)=\underset{i=1}{\overset{n}{\Σ}}{f}_{EMD}^i\left(\stackrel{\·\·}{x}\right(t\left)\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为加速度响应信号的最大时间长度；

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

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

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

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

ΔK＝K₀-K_s；

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

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

$f_s=\frac{1}{2\mathrm{\π}}\sqrt{\frac{Ks}{M}};$

式中，M为结构的质量；

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

$\mathrm{\Δ}K=K_0-K_s=4{\mathrm{\π}}^{2}M(f_0^2-f_s^2);$

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

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

$\stackrel{\·\·}{x}+4{\mathrm{\π\ξf}}_0\stackrel{\·}{x}+4{\mathrm{\π}}^2{f}_0^{2}x=0;$

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

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

$x\left(t\right)=\frac{v_{0}sin\left(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_{d}t\right)}{2{\mathrm{\πf}}_0{\mathrm{\ξ}}_d}\·e^{-2{\mathrm{\πf}}_0\mathrm{\ξ}t};$

$\stackrel{\·}{x}\left(t\right)=v_0{e}^{-2{\mathrm{\πf}}_0\mathrm{\ξ}t}\left(cos\right(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_{d}t)-\frac{\mathrm{\ξ}}{{\mathrm{\ξ}}_d}sin(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_{d}t\left)\right);$

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

式中：

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

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

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

${\stackrel{\·\·}{x}}_s+4{\mathrm{\π\ξf}}_s{\stackrel{\·}{x}}_s+4{\mathrm{\π}}^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{\πf}}_0{\mathrm{\ξ}}_d{t}_i\right)}{2{\mathrm{\πf}}_0{\mathrm{\ξ}}_d}\·v_0{e}^{-2{\mathrm{\πf}}_0{\mathrm{\ξt}}_i};$

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

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

${\stackrel{\·\·}{x}}_s\left(t_1\right)=-\frac{2{\mathrm{\πv}}_0{f}_s(E_1+E_2+E_3)}{{\mathrm{\ξ}}_d^2{f}_0}\·e^{-2\mathrm{\π}\mathrm{\ξ}(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{\πf}}_s{\mathrm{\ξ}}_d{t}_1\right)(1-2{\mathrm{\ξ}}_d^2)\left(\mathrm{\ξ}sin\right(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_d{t}_i)-{\mathrm{\ξ}}_{d}cos(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_d{t}_i\left)\right);$

$E_3=-2f_0\mathrm{\ξ}cos\left(2{\mathrm{\πf}}_s{\mathrm{\ξ}}_d{t}_1\right)\[{\mathrm{\ξ\ξ}}_{d}sin\left(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_d{t}_i\right)-{\mathrm{\ξ}}_d^{2}cos\left(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_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{\·\·}{x}\left(t_{i+1}\right)={\stackrel{\·\·}{x}}_s\left(\mathrm{\Δ}t\right)=-\frac{2{\mathrm{\πf}}_s{v}_0{e}^{-2{\mathrm{\πf}}_0{\mathrm{\ξt}}_i}}{f_0{\mathrm{\ξ}}_d}\{f_s\sin \left(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_d{t}_i\right)-2f_0\mathrm{\ξ}\[\mathrm{\ξ}\sin \left(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_d{t}_i\right)-{\mathrm{\ξ}}_d\cos \left(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_d{t}_i\right)\]\};$

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

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

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

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

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

$h_1\left(t\right)=\stackrel{\·\·}{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{\·\·}{x}\left(t\right)=\underset{i=1}{\overset{n}{\Σ}}{c}_i\left(t\right)+r_n\left(t\right);$

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

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

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

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

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

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

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

$c_i^s\left(t\right)=f_{EMD}^i(-\frac{2{\mathrm{\πf}}_s{v}_0{e}^{-2{\mathrm{\πf}}_0{\mathrm{\ξt}}_i}}{f_0{\mathrm{\ξ}}_d}\{f_s\sin \left(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_d{t}_i\right)-2f_0\mathrm{\ξ}\[\mathrm{\ξ}\sin \left(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_d{t}_i\right)-{\mathrm{\ξ}}_d\cos \left(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_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{\Δc}}_1}{\mathrm{\Δ}t}=\frac{c_{1,i+1}-c_{1,i}}{\mathrm{\Δ}t}=f_{EMD}^1({\stackrel{\·\·}{x}}_{i+1}-{\stackrel{\·\·}{x}}_i)\\ \stackrel{\·}{=}{f}_{EMD}^1(-\frac{2{\mathrm{\π\ξ}}_d{v}_0(f_s^2-f_0^2)}{f_0}\·e^{-2{\mathrm{\πf}}_0{\mathrm{\ξt}}_i}\·\sin \left(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_d{t}_i\right)),(i=1,2,\mathrm{...},t_{\max }-1)\end{array};$

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

$f_s^2-f_0^2=\frac{\mathrm{\Δ}K}{4{\mathrm{\π}}^{2}M};$

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

$D_i\stackrel{\·}{=}{f}_{EMD}^1(-\frac{{\mathrm{\ΔK\ξ}}_d{v}_0}{2{\mathrm{\πf}}_{0}M}\·e^{-2{\mathrm{\πf}}_0{\mathrm{\ξt}}_i}\·sin(2{\mathrm{\πf}}_0{\mathrm{\ξ}}_d{t}_i\left)\right),(i=1,2,\mathrm{...},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=1,2,\mathrm{...},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=2,3,\mathrm{...},t_{max}-1);$

因此有：

${\mathrm{MI}}_i=|2D_i-D_{i-1}-D_{i+1}|\stackrel{\·}{=}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{\·}{=}{\mathrm{MI}}_i;$

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

建立线性模型：

MI＝α*S+β；

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

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