CN103902834A - 一种基于岭估计和l曲线法的结构损伤识别方法 - Google Patents

Abstract

本发明公开了一种基于岭估计和L曲线法的结构损伤识别方法，首先根据结构模态应变能的灵敏度建立基于模态应变能灵敏度的结构损伤方程；然后根据Tikhonov正则化原理，确定岭估计的基本估计准则；再利用L曲线法确定岭估计的最优岭参数，基于最优岭参数计算出结构损伤方程的基本解，最后对结果进行修正，即可得到损伤单元的损伤系数，从而识别结构损伤。本发明采用岭估计可以大大降低病态问题对损伤识别的影响，采用L曲线法确定了近似最优的岭参数，同时提出了岭估计的修正策略，从而获得更优的损伤估计值。

Description

一种基于岭估计和L曲线法的结构损伤识别方法

技术领域

本发明涉及结构损伤识别，特别涉及一种基于岭估计和L曲线法的结构损伤识别方法，属于结构损伤检测技术。

背景技术

工程结构在投入使用后，由于环境的作用，常会产生各种损伤。结构的严重损伤会引起建筑物的倒塌或者失效，从而造成巨大的经济损失和人员的生命安全。因此，对结构损伤的识别研究一直是国际上的研究热点。利用损伤前后动力特性的改变，建立相应的损伤方程，从而求解出损伤系数是较为有效的方法，根据损伤系数即可判断结构的损伤程度。但是，在工程实际中，测试数据不可避免的会受到测量噪声、环境因素等的干扰，常会使损伤方程产生病态问题，从而给直接求解带来困难，造成损伤的解极不稳定，与真实值相差较远。这样，损伤方程组的求解和精度成为问题的症结。采用截断奇异值分解算法(TSVD)能够减少方程的参数空间，使损伤方程组变得数值稳定。但是，该方法需要人为的选择截断误差或者截断水平，不同的截断误差或截断水平会产生不同的识别结果，故选择恰当的截断误差或者截断水平决定了方程求解的精度。而结构的改变，甚至同一结构单元划分不同，都会影响截断误差或者截断水平的选取，从而带来不同的识别结果。

总之，现有的基于应变能灵敏度的损伤识别技术存在两种缺陷：

1）直接求逆的应变能灵敏度损伤方程求解方法难于处理病态损伤问题；

2）当采用单纯的截断奇异值分解算法处理应变能灵敏度损伤方程问题时，难于确定恰当的截断误差，从而容易造成识别精度欠佳的问题。

发明内容

针对现有技术存在的上述不足，本发明的目的在于提供一种可以大大降低病态问题对损伤识别的影响、从而提高损伤识别精度的基于岭估计和L曲线法的结构损伤识别方法。

本发明的技术方案是这样实现的：

一种基于岭估计和L曲线法的结构损伤识别方法，步骤如下，

1）首先根据结构模态应变能的灵敏度建立基于模态应变能灵敏度的结构损伤方程；

2）然后根据Tikhonov正则化原理，确定采用岭估计法求解结构损伤方程时岭估计的基本估计准则；

3）再利用L曲线法确定岭估计的最优岭参数，基于最优岭参数计算出结构损伤方程的基本解，该基本解即为结构损伤估计值，从而识别结构损伤。

在第3）步获得结构损伤方程的基本解后根据无损单元估计值的偏移量来确定修正的幅值，再结合岭估计的范数特性，对损伤单元的损伤系数进行修正。本发明中基本解、损伤估计值和损伤系数实质上是一个概念，基本解是从方程角度而言，损伤估计值是从岭估计角度而言，在本方法中实际上都是指损伤系数。

其中，第1）步建立基于模态应变能灵敏度的结构损伤方程具体过程如下，

第j个单元和第r阶模态的模态应变能

$E_{\mathrm{rj}}=\frac{1}{2}{\mathrm{\Φ}}_r^T{K}_j{\mathrm{\Φ}}_r---\left(1\right)$

式中，E_rj是模态应变能；Φ_r为第r阶质量归一化位移模态；K_j为第j个单元的刚度矩阵；则第j个单元和第r阶模态的模态应变能对任意设计参数p的敏感性为

$\frac{{\∂E}_{\mathrm{rj}}}{\∂p}={\mathrm{\Φ}}_r^T\left[\begin{array}{cc}{K}_j& 0\end{array}\right]{K_{U/D}}^{-1}{F}_p{\mathrm{\Φ}}_r+\frac{1}{2}{\mathrm{\Φ}}_r^T\frac{{\∂K}_j}{\∂p}{\mathrm{\Φ}}_r---\left(2\right)$

其中， $K_{U/D}=\left[\begin{array}{cc}K-{\mathrm{\λ}}_{r}M& -M{\mathrm{\Φ}}_r\\ -{\mathrm{\Φ}}_r^{T}M& 0\end{array}\right]$

$F_p=\left[\begin{array}{c}-(\frac{\∂K}{\∂p}-{\mathrm{\λ}}_r\frac{\∂M}{\∂p})\\ \frac{1}{2}{\mathrm{\Φ}}_r^T\frac{\∂M}{\∂p}\end{array}\right]$

式中，λ_r是第r阶特征值；K和M是整体刚度矩阵和质量矩阵；p为设计参数；

当结构损伤时，其第j个单元损伤后的刚度矩阵

K_j ^d=K_j-c_jK_j；              (3)

式中，c_j为第j个单元的损伤系数，其取值范围为(0≤c_j≤1)；由单元模态应变能的灵敏度表达式可以推导出基于模态应变能变化的损伤方程组

S△C=△R                (4)

式中，S是系数矩阵；△R是残余向量；△C是损伤系数向量；S矩阵的组成具体如下：

$S_{\mathrm{jir}}=\left\{\begin{array}{cc}{\mathrm{\Φ}}_r^T\mathrm{\Λ}{\mathrm{\Φ}}_r,& i\≠j\\ {\mathrm{\Φ}}_r^T\mathrm{\Λ}{\mathrm{\Φ}}_r+\frac{1}{2}{\mathrm{\Φ}}_r^T{K}_j{\mathrm{\Φ}}_r-\frac{1}{2}{\mathrm{\Φ}}_r^{\mathrm{dT}}{K}_j{\mathrm{\Φ}}_r^d,& i=j\end{array}\right.---\left(5\right)$

$\mathrm{\Δ}{R}_{\mathrm{jr}}=\frac{1}{2}{\mathrm{\Φ}}_r^T{K}_j{\mathrm{\Φ}}_r-\frac{1}{2}{\mathrm{\Φ}}_r^{\mathrm{dT}}{K}_j{\mathrm{\Φ}}_r^d---\left(6\right)$

其中， $\mathrm{\Λ}=\left[\begin{array}{cc}{K}_j& 0\end{array}\right]{K_{U/D}}^{-1}\left\{\begin{array}{c}-K_i\\ 0\end{array}\right\};$

为损伤后第r阶位移模态；

式（4）即为基于模态应变能灵敏度的结构损伤方程。

其中，第2）步对应于结构损伤方程岭估计的基本估计准则为：

$\begin{array}{c}{\left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|}^2+\mathrm{\α\Ω}\left(\mathrm{\Δ}\hat{C}\right)={\left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|}^2+\mathrm{\α}{\left(\mathrm{\Δ}\hat{C}\right)}^T\left(\mathrm{\ΔC}\right)\\ ={\left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|}^2+\mathrm{\α}{\left|\right|\mathrm{\Δ}\widehat{C\left|\right|}}^2=\min \end{array}---\left(7\right)$

式中，

是损伤系数向量△C的估计值；α为岭参数；||·||为欧氏2-范数；Ω(△C)为稳定泛函；通过上述岭估计的基本估计准则可得到岭估计的解为

$\mathrm{\Δ}\hat{C}={(S^{T}S+\mathrm{\αI})}^{-1}{S}^T\mathrm{\ΔR}---\left(8\right)$

式中，I为单位阵。

其中，第3）步利用L曲线法确定岭估计的最优岭参数方法为，

令 $\mathrm{\η}={\left|\right|\mathrm{\Δ}\hat{C}\left|\right|}^2,\mathrm{\ρ}={\left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|}^2---\left(9\right)$

两边取对数，可得

$\widehat{\mathrm{\η}}=\mathrm{lg\η}=2\lg \left|\right|\mathrm{\Δ}\hat{C}\left|\right|,\widehat{\mathrm{\ρ}}=\mathrm{lg\ρ}=2\lg \left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|---\left(10\right)$

则L曲线由许多点

拟合而成；L曲线上点的曲率κ的计算公式为

$\mathrm{\κ}=2\frac{{\widehat{\mathrm{\ρ}}}^{\′}{\widehat{\mathrm{\η}}}^{\′\′}-{{\widehat{\mathrm{\ρ}}}^{\′\′}\widehat{\mathrm{\η}}}^{\′}}{{({\left({\widehat{\mathrm{\ρ}}}^{\′}\right)}^2+{\left({\widehat{\mathrm{\η}}}^{\′}\right)}^2)}^{3/2}}---\left(11\right)$

在式(11)中，

分别为

的一阶和二阶导数，并且都是岭参数α的函数，对式(11)求最大值，就可以得到L曲线的最大曲率κ_max，所对应点的α值就是所求的最优岭参数。

其中，对损伤单元的损伤系数进行修正的修正公式如下

${\hat{c}}_j={\hat{c}}_j+\frac{({\hat{c}}_j-d)}{d}\mathrm{mean}\left(\mathrm{\Δ}\stackrel{\‾}{C}\right)---\left(12\right)$

式中，是第j个单元损伤系数的估计值；

是损伤向量

的偏移估计量，是由

中所有小于0的无损单元估计值缩减而成的向量；

表示无损单元偏移估计量的均值；d是基准参数。

该基于岭估计和L曲线的损伤识别方法可以较好的进行损伤识别。其具体优点如下：

1、对于应变能敏感度损伤方程受测量噪声等产生的病态问题，采用本方法的岭估计方法可以大大降低该类病态问题对求解的影响。

2、相对于简单的岭估计，本方法采用了L曲线法确定了近似最优的岭参数，从而获得更优的损伤估计值。

3、相对于基本的岭估计和L曲线相结合方法，本方法提出了岭估计的修正策略，该修正策略可以减少损伤估计值的线性偏差，从而进一步提高了损伤的识别精度。

附图说明

图1-本发明损伤识别流程图。

图2-实施例Euler-Bernoulli简支梁结构示意图。

图3-第一种损伤情况下的四种损伤识别方案的识别结果。

图4-第一种损伤情况下的四种损伤识别方案的识别结果（考虑3%的测量噪声）。

图5-第二种损伤情况下的四种损伤识别方案的识别结果。

图6-第二种损伤情况下的四种损伤识别方案的识别结果（考虑3%的测量噪声）。

图7-第三种损伤情况下的四种损伤识别方案的识别结果。

图8-第三种损伤情况下的四种损伤识别方案的识别结果（考虑3%的测量噪声）。

具体实施方式

本发明的目的是提出一种新的应变能灵敏度损伤方程求解方法，并且提供一种改进策略来提高损伤识别的精度。本方法不仅可以处理损伤方程的病态问题，而且还可以通过对估计值偏差的修正策略来提高损伤识别的精度。由于岭估计（RE）法可以较好的处理病态问题，而L曲线方法可以较好地处理优化参数的选取问题，故本发明采用岭估计和L曲线法相结合，来解决损伤方程的病态问题以及优化参数的选取问题。具体实现思路可以参见图1，从图1可以看出，本方法主要是采用基于模态应变能灵敏度损伤方程和岭估计方法进行损伤识别研究，首先提取结构的加速度响应数据，生成基本的刚度矩阵，基于加速度响应数据和刚度矩阵得到测量结构的频率和振型模态，从而得到模态应变能的表达方式，并根据模态应变能的灵敏度建立相应的损伤方程，然后根据Tikhonov正则化原理，确定岭估计的基本估计准则，再利用L曲线法确定岭估计的最优岭参数，从而可以计算出结构损伤的基本解，最后采用岭估计的修正策略来获取更精确的损伤识别结果。

具体如下：

1基于模态应变能灵敏度的结构损伤方程

基于模态应变能的损伤识别方法是较为有效的一种方法，可以根据模态应变能的灵敏度建立相应的损伤方程。

第j个单元和第r阶模态的模态应变能

$E_{\mathrm{rj}}=\frac{1}{2}{\mathrm{\Φ}}_r^T{K}_j{\mathrm{\Φ}}_r---\left(1\right)$

式中，E_rj是模态应变能；Φ_r为第r阶质量归一化位移模态；K_j为第j个单元的刚度矩阵。则第j个单元和第r阶模态的模态应变能对任意设计参数p的敏感性为

$\frac{{\∂E}_{\mathrm{rj}}}{\∂p}={\mathrm{\Φ}}_r^T\left[\begin{array}{cc}{K}_j& 0\end{array}\right]{K_{U/D}}^{-1}{F}_p{\mathrm{\Φ}}_r+\frac{1}{2}{\mathrm{\Φ}}_r^T\frac{{\∂K}_j}{\∂p}{\mathrm{\Φ}}_r---\left(2\right)$

其中， $K_{U/D}=\left[\begin{array}{cc}K-{\mathrm{\λ}}_{r}M& -M{\mathrm{\Φ}}_r\\ -{\mathrm{\Φ}}_r^{T}M& 0\end{array}\right]$

$F_p=\left[\begin{array}{c}-(\frac{\∂K}{\∂p}-{\mathrm{\λ}}_r\frac{\∂M}{\∂p})\\ \frac{1}{2}{\mathrm{\Φ}}_r^T\frac{\∂M}{\∂p}\end{array}\right]$

式中，λ_r是第r阶特征值；K和M是整体刚度矩阵和质量矩阵；p为设计参数。

当结构损伤时，其第j个单元损伤后的刚度矩阵

K_j ^d=K_j-c_jK_j;           (3)

式中，c_j为第j个单元的损伤系数，其取值范围为(0≤c_j≤1)。由单元模态应变能的灵敏度表达式可以推导出基于模态应变能变化的损伤方程组

S△C=△R        (4)

式中，S是系数矩阵；△R是残余向量；△C是损伤系数向量。S矩阵的组成具体如下：

$S_{\mathrm{jir}}=\left\{\begin{array}{cc}{\mathrm{\Φ}}_r^T\mathrm{\Λ}{\mathrm{\Φ}}_r,& i\≠j\\ {\mathrm{\Φ}}_r^T\mathrm{\Λ}{\mathrm{\Φ}}_r+\frac{1}{2}{\mathrm{\Φ}}_r^T{K}_j{\mathrm{\Φ}}_r-\frac{1}{2}{\mathrm{\Φ}}_r^{\mathrm{dT}}{K}_j{\mathrm{\Φ}}_r^d,& i=j\end{array}\right.---\left(5\right)$

$\mathrm{\Δ}{R}_{\mathrm{jr}}=\frac{1}{2}{\mathrm{\Φ}}_r^T{K}_j{\mathrm{\Φ}}_r-\frac{1}{2}{\mathrm{\Φ}}_r^{\mathrm{dT}}{K}_j{\mathrm{\Φ}}_r^d---\left(6\right)$

其中， $\mathrm{\Λ}=\left[\begin{array}{cc}{K}_j& 0\end{array}\right]{K_{U/D}}^{-1}\left\{\begin{array}{c}-K_i\\ 0\end{array}\right\};$

为损伤后第r阶位移模态。

通过求解该方程可以同时识别结构的损伤位置和程度。但是，由于结构损伤常常出现在结构的局部某些位置或单元，而大部分结构单元未发生损伤，因此，在系数矩阵中常会出现几乎线性相关的列，特别是模态参数在测试时会出现噪声的干扰，故方程组的求解和精度往往成为问题的关键症结。采用截断奇异值分解法可以较为有效的处理该问题，但是，由于在现实中，损伤的位置和程度是未知的，截断奇异值分解法选择不同的截断误差或者截断水平，所得计算结果往往差别很大，从而对损伤识别结果造成较大的干扰。故本发明采取岭估计和L曲线法相结合，来求解损伤方程，以获取更具有可靠性的计算结果。

2岭估计和L曲线法及其改进

2.1损伤的岭估计

可以采用岭估计法求解损伤方程，以避免因测量噪声等引起的损伤方程病态问题。岭估计方法是一种在均方误差意义下优于最小二乘估计的有偏估计方法，是从Tikhonov正则化原理推导出来，该方法可以克服矩阵方程的病态性。根据Tikhonov正则化原理，对应于损伤方程岭估计的估计准则为：

$\begin{array}{c}{\left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|}^2+\mathrm{\α\Ω}\left(\mathrm{\Δ}\hat{C}\right)={\left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|}^2+\mathrm{\α}{\left(\mathrm{\Δ}\hat{C}\right)}^T\left(\mathrm{\ΔC}\right)\\ ={\left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|}^2+\mathrm{\α}{\left|\right|\mathrm{\Δ}\widehat{C\left|\right|}}^2=\min \end{array}---\left(7\right)$

式中，是损伤系数向量△C的估计值；α为岭参数；||·||为欧氏2-范数；Ω(△C)为稳定泛函。则可得到岭估计的解为

$\mathrm{\Δ}\hat{C}={(S^{T}S+\mathrm{\αI})}^{-1}{S}^T\mathrm{\ΔR}---\left(8\right)$

式中，I为单位阵。αI项的参与，可以有效的抑制方程的病态性，从而可以得到可靠的损伤估计值。岭估计的核心是如何确定岭参数，本发明采用L曲线法来选取恰当的岭参数。

2.2L曲线法

在式(7)中，

和都是岭参数α的函数，通过选择不同的α值，可绘出残余范数

和正则化解范数

之间的二维曲线，由于该曲线通常呈L形，故称为L曲线，其拐角为岭参数的最优值，在此处可求出相应的最优岭参数α。该方法的关键是定位L曲线上的曲率最大点，可通过对数形式来求解。具体如下：

令 $\mathrm{\η}={\left|\right|\mathrm{\Δ}\hat{C}\left|\right|}^2,\mathrm{\ρ}={\left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|}^2---\left(9\right)$ 两边取对数，可得

$\widehat{\mathrm{\η}}=\mathrm{lg\η}=2\lg \left|\right|\mathrm{\Δ}\hat{C}\left|\right|,\widehat{\mathrm{\ρ}}=\mathrm{lg\ρ}=2\lg \left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|---\left(10\right)$

则L曲线由许多点拟合而成。则L曲线上点的曲率κ的计算公式为[9]

$\mathrm{\κ}=2\frac{{\widehat{\mathrm{\ρ}}}^{\′}{\widehat{\mathrm{\η}}}^{\′\′}-{{\widehat{\mathrm{\ρ}}}^{\′\′}\widehat{\mathrm{\η}}}^{\′}}{{({\left({\widehat{\mathrm{\ρ}}}^{\′}\right)}^2+{\left({\widehat{\mathrm{\η}}}^{\′}\right)}^2)}^{3/2}}---\left(11\right)$

在式(11)中，

分别为

的一阶和二阶导数，并且都是岭参数α的函数，对式(11)求最大值，就可以得到L曲线的最大曲率κ_max,所对应点的α值就是所求的岭参数。

2.3估计值的修正

由于岭估计是有偏估计，而L曲线法求解的也是近似最优解，故损伤的估计值需要修正。由岭估计的估计准则方程式(7)可知，该准则是用来确定一次损伤系数的2-范数项最小，这样，该岭估计的解会造成近似线性的偏差，故本文采用了线性修正的方式来提高识别精度。损伤系数c_j应大于0，当估计值小于0时，则认为是无损伤，而该项本身是由于有偏估计造成的，故可根据无损单元估计值的偏移量来确定修正的幅值，再结合岭估计的范数特性，可以对可能损伤单元的损伤系数进行修正，修正公式如下

${\hat{c}}_j={\hat{c}}_j+\frac{({\hat{c}}_j-d)}{d}\mathrm{mean}\left(\mathrm{\Δ}\stackrel{\‾}{C}\right)---\left(12\right)$

式中，

是第j个单元损伤系数的估计值，左边是修正后的损伤估计值，右边是修正前的损伤估计值；

是损伤向量的偏移估计量，是由

中所有小于0的无损单元估计值缩减而成的向量；

表示无损单元偏移估计量的均值；d是基准参数。采用该修正公式可以减少估计值的线性偏差，从而提高损伤的识别精度。

下面结合具体实施例对本发明做进一步详细说明。

假设有如图2所示的一个Euler-Bernoulli简支梁。梁总长为6m，划分为20个单元，假设未损伤结构的几何参数和物理参数如下：梁的横截面积A为0.005m²，惯性矩I=1.67m⁴，弹性模量E=32GPa，密度ρ=2500kg/m³，基准参数为0.25。假设有三种损伤情况，第一种情况，在单元3和12发生损伤，刚度分别降低15%和20%；第二种情况，在单元2、10、14发生损伤，刚度分别降低25%，15%和20%；第三种情况，在单元5、11、19发生损伤，刚度分别降低20%，20%和30%。分别采用截断奇异值分解(TSVD)法和岭估计(RE)法求解损伤方程。由于低阶模态比高阶模态更能反映结构的损伤，其抗噪能力更强，而且更易被测量和获取。因此，本发明仅采用第1阶位移模态进行结构的损伤识别研究。岭参数的确定主要是通过对公式(11)求最大值κ_max，κ_max所对应的点就是L曲线的所求的点。这样就定位了L曲线上曲率最大的点，这个点所对应的α值就是所求的岭参数。

第一种情况

当单元3和12发生损伤时，其损伤系数的计算结果如图3所示。采用了四种方案（其它两种情况相同），即截断误差0.01的TSVD法、截断误差0.05的TSVD法、基本RE和L曲线法，修正RE和L曲线法。从图3中可以发现，截断误差0.01的TSVD法的损伤定位结果不理想，易对单元1和20产生误判，其余三种方法定位识别较好，而从定量角度分析，可以较为明显地发现修正RE和L曲线法的定量识别结果最好。将单元3和12的损伤系数单独提取分析，结果如表1所示。从表1可以观察到，修正RE和L曲线法的平均误差最小，识别结果最好。

表1.第一种情况的损伤定量识别结果

当考虑位移模态有3%的测量噪声干扰时，则四种损伤识别方案的一次计算结果如图4所示。首先从损伤定位角度分析，截断误差0.01的TSVD法已经丧失了损伤定位功能，截断误差0.05的TSVD法则具有一定的损伤定位效果，但是易对单元2和19产生误判，而两种RE和L曲线法的损伤定位效果较好。然后从损伤定量角度分析，从图4中可以发现，本发明的修正RE和L曲线法相对较好。故修正RE和L曲线法不仅具有较好的定位能力，具有一定的定量分析能力。由于测量噪声具有一定的随机性，则每次考虑随机噪声的计算结果都会产生一些小的变化，故未列出定量结果的平均误差表。

第二种情况

当单元2、10和14发生损伤时，其损伤系数的计算结果如图5所示。从图5中可以发现，截断误差0.01的TSVD法的损伤定位结果不理想，易对单元1、19、20产生误判，其余三种方法的定位识别结果较好。再从定量角度分析，可以较为明显地发现修正RE和L曲线法的定量识别结果最好。将单元2、10和14的损伤系数单独提取分析，结果如表2所示。从表2可以观察到，修正RE和L曲线法的平均误差最小，识别结果最好。

表2.第二种情况的损伤定量识别结果

当考虑位移模态有3%的测量噪声干扰时，计算结果如图6所示。首先从损伤定位角度分析，截断误差0.01的TSVD法已经丧失了损伤定位功能，截断误差0.05的TSVD法则具有一定的损伤定位效果，但是易对单元3、18和19产生误判，而两种RE和L曲线法的损伤定位效果较好。然后从损伤定量角度分析，从图6中可以发现，本发明的修正RE和L曲线法相对较好。由于测量噪声具有一定的随机性，故未列出定量结果的平均误差表。

第三种情况

当单元5、11和19发生损伤时，其损伤系数的计算结果如图7所示。从图7中依然可以发现，截断误差0.01的TSVD法的损伤定位结果仍不理想，易对单元1、2和20产生误判，其余三种方法的定位识别结果较好。再从定量角度分析，可以较为明显地发现修正RE和L曲线法的定量识别结果最好。将单元5、11和19的损伤系数单独提取分析，结果如表3所示。从表3可以观察到，本文建议的修正RE和L曲线法的平均误差最小，识别结果最好。

表3.第三种情况的损伤定量识别结果

当考虑位移模态有3%的测量噪声干扰时，计算结果如图8所示。首先从损伤定位角度分析，截断误差0.01的TSVD法已经基本丧失了损伤定位功能，截断误差0.05的TSVD法则具有一定的损伤定位效果，但是易对单元2、3和18产生误判，而基本和修正的RE和L曲线法的损伤定位效果较好。然后从损伤定量角度分析，从图8中可以发现，本发明的修正RE和L曲线法相对较好。由于测量噪声具有一定的随机性，故未列出定量结果的平均误差表。

从以上Euler-Bernoulli简支梁三个不同损伤情况的例子可以观察到，当仅采用第一阶位移模态进行损伤识别时，截断误差为0.01时，TSVD法的识别效果较差，截断误差为0.05时，TSVD法的识别效果相对较好，故TSVD法的关键是选取恰当的截断误差值。由于结构损伤是未知的，故截断误差的选择也应存在一个优化问题。当采用基本岭估计和L曲线法，其定位和定量结果相对更好，这是因为L曲线可以选择近似最优的岭参数，从而解决了参数优化问题。修正岭估计和L曲线法的识别效果最好，这是由于该方法采用了未损单元的有偏估计均值对损伤系数估计值进行了修正和改进。综上所述，基本岭估计和L曲线法可以较可靠地求解出损伤方程，并较好地识别出结构的损伤，而修正的岭估计和L曲线法可以更可靠地进行损伤的定性和定量识别，该方法优于基本岭估计和L曲线法以及未经过优化的截断奇异值分解法。

本发明的上述实施例仅仅是为说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其他不同形式的变化和变动。这里无法对所有的实施方式予以穷举。凡是属于本发明的技术方案所引申出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (6)

1.一种基于岭估计和L曲线法的结构损伤识别方法，其特征在于：步骤如下，

1）首先根据结构模态应变能的灵敏度建立基于模态应变能灵敏度的结构损伤方程；

2）然后根据Tikhonov正则化原理，确定采用岭估计法求解结构损伤方程时岭估计的基本估计准则；

3）再利用L曲线法确定岭估计的最优岭参数，基于最优岭参数计算出结构损伤方程的基本解，该基本解即为结构损伤估计值，从而识别结构损伤。

2.根据权利要求1所述的基于岭估计和L曲线法的结构损伤识别方法，其特征在于：在第3）步获得结构损伤方程的基本解后根据无损单元估计值的偏移量来确定修正的幅值，再结合岭估计的范数特性，对损伤单元的损伤系数进行修正。

3.根据权利要求1所述的基于岭估计和L曲线法的结构损伤识别方法，其特征在于：第1）步建立基于模态应变能灵敏度的结构损伤方程具体过程如下，

第j个单元和第r阶模态的模态应变能

$E_{\mathrm{rj}}=\frac{1}{2}{\mathrm{\Φ}}_r^T{K}_j{\mathrm{\Φ}}_r---\left(1\right)$

式中，E_rj是模态应变能；Φ_r为第r阶质量归一化位移模态；K_j为第j个单元的刚度矩阵；则第j个单元和第r阶模态的模态应变能对任意设计参数p的敏感性为

$\frac{{\∂E}_{\mathrm{rj}}}{\∂p}={\mathrm{\Φ}}_r^T\left[\begin{array}{cc}{K}_j& 0\end{array}\right]{K_{U/D}}^{-1}{F}_p{\mathrm{\Φ}}_r+\frac{1}{2}{\mathrm{\Φ}}_r^T\frac{{\∂K}_j}{\∂p}{\mathrm{\Φ}}_r---\left(2\right)$

其中， $K_{U/D}=\left[\begin{array}{cc}K-{\mathrm{\λ}}_{r}M& -M{\mathrm{\Φ}}_r\\ -{\mathrm{\Φ}}_r^{T}M& 0\end{array}\right]$

$F_p=\left(\begin{array}{c}-(\frac{\∂K}{\∂p}-{\mathrm{\λ}}_r\frac{\∂M}{\∂p})\\ \frac{1}{2}{\mathrm{\Φ}}_r^T\frac{\∂M}{\∂p}\end{array}\right)$

式中，λ_r是第r阶特征值；K和M是整体刚度矩阵和质量矩阵；p为设计参数；

当结构损伤时，其第j个单元损伤后的刚度矩阵

K_j ^d=K_j-c_jK_j；              (3)

式中，c_j为第j个单元的损伤系数，其取值范围为(0≤c_j≤1)；由单元模态应变能的灵敏度表达式可以推导出基于模态应变能变化的损伤方程组

S△C=△R               (4)

式中，S是系数矩阵；△R是残余向量；△C是损伤系数向量；S矩阵的组成具体如下：

$S_{\mathrm{jir}}=\left\{\begin{array}{cc}{\mathrm{\Φ}}_r^T\mathrm{\Λ}{\mathrm{\Φ}}_r,& i\≠j\\ {\mathrm{\Φ}}_r^T\mathrm{\Λ}{\mathrm{\Φ}}_r+\frac{1}{2}{\mathrm{\Φ}}_r^T{K}_j{\mathrm{\Φ}}_r-\frac{1}{2}{\mathrm{\Φ}}_r^{\mathrm{dT}}{K}_j{\mathrm{\Φ}}_r^d,& i=j\end{array}\right.---\left(5\right)$

$\mathrm{\Δ}{R}_{\mathrm{jr}}=\frac{1}{2}{\mathrm{\Φ}}_r^T{K}_j{\mathrm{\Φ}}_r-\frac{1}{2}{\mathrm{\Φ}}_r^{\mathrm{dT}}{K}_j{\mathrm{\Φ}}_r^d---\left(6\right)$

其中， $\mathrm{\Λ}=\left[\begin{array}{cc}{K}_j& 0\end{array}\right]{K_{U/D}}^{-1}\left\{\begin{array}{c}-K_i\\ 0\end{array}\right\};$

为损伤后第r阶位移模态；

式（4）即为基于模态应变能灵敏度的结构损伤方程。

4.根据权利要求3所述的基于岭估计和L曲线法的结构损伤识别方法，其特征在于：第2）步对应于结构损伤方程岭估计的基本估计准则为：

$\begin{array}{c}{\left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|}^2+\mathrm{\α\Ω}\left(\mathrm{\Δ}\hat{C}\right)={\left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|}^2+\mathrm{\α}{\left(\mathrm{\Δ}\hat{C}\right)}^T\left(\mathrm{\ΔC}\right)\\ ={\left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|}^2+\mathrm{\α}{\left|\right|\mathrm{\Δ}\widehat{C\left|\right|}}^2=\min \end{array}---\left(7\right)$

式中，

是损伤系数向量△C的估计值；α为岭参数；||·||为欧氏2-范数；Ω(△C)为稳定泛函；通过上述岭估计的基本估计准则可得到岭估计的解为

$\mathrm{\Δ}\hat{C}={(S^{T}S+\mathrm{\αI})}^{-1}{S}^T\mathrm{\ΔR}---\left(8\right)$

式中，I为单位阵。

5.根据权利要求4所述的基于岭估计和L曲线法的结构损伤识别方法，其特征在于：第3）步利用L曲线法确定岭估计的最优岭参数方法为，

令 $\mathrm{\η}={\left|\right|\mathrm{\Δ}\hat{C}\left|\right|}^2,\mathrm{\ρ}={\left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|}^2---\left(9\right)$

两边取对数，可得

$\widehat{\mathrm{\η}}=\mathrm{lg\η}=2\lg \left|\right|\mathrm{\Δ}\hat{C}\left|\right|,\widehat{\mathrm{\ρ}}=\mathrm{lg\ρ}=2\lg \left|\right|\mathrm{S\Δ}\hat{C}-\mathrm{\ΔR}\left|\right|---\left(10\right)$

则L曲线由许多点

拟合而成；L曲线上点的曲率κ的计算公式为

$\mathrm{\κ}=2\frac{{\widehat{\mathrm{\ρ}}}^{\′}{\widehat{\mathrm{\η}}}^{\′\′}-{{\widehat{\mathrm{\ρ}}}^{\′\′}\widehat{\mathrm{\η}}}^{\′}}{{({\left({\widehat{\mathrm{\ρ}}}^{\′}\right)}^2+{\left({\widehat{\mathrm{\η}}}^{\′}\right)}^2)}^{3/2}}---\left(11\right)$

在式(11)中，

分别为

的一阶和二阶导数，并且都是岭参数α的函数，对式(11)求最大值，就可以得到L曲线的最大曲率κ_max，所对应点的α值就是所求的最优岭参数。

6.根据权利要求2所述的基于岭估计和L曲线法的结构损伤识别方法，其特征在于：对损伤单元的损伤系数进行修正的修正公式如下

${\hat{c}}_j={\hat{c}}_j+\frac{({\hat{c}}_j-d)}{d}\mathrm{mean}\left(\mathrm{\Δ}\stackrel{\‾}{C}\right)---\left(12\right)$

式中，

是第j个单元损伤系数的估计值；

是损伤向量

的偏移估计量，是由

中所有小于0的无损单元估计值缩减而成的向量；

表示无损单元偏移估计量的均值；d是基准参数。

Images