CN102759573B - 基于频率变化的结构损伤定位与损伤程度的评估方法 - Google Patents

Abstract

一种基于频率变化的结构损伤定位与损伤程度的评估方法，其特征在于：采用以下步骤：第一步：建立未损伤结构的有限元模型，作为基准有限元模型；第二步：对损伤结构进行振动测试，识别损伤结构的m阶模态频率第三步：对损伤结构进行损伤程度评估；第四步：对损伤结构进行损伤位置确定。本发明不需要结构的振型信息，仅利用结构损伤前后的有限的低阶频率信息即可准确地进行损伤定位和损伤程度评估；适用于单损伤、多种损伤工况，具有一定的实际应用价值。

Description

基于频率变化的结构损伤定位与损伤程度的评估方法

技术领域

本发明涉及结构损伤定位与损伤程度评估的方法，尤其涉及一种可以利用损伤结构实测的低阶频率来确定结构的损伤位置，并评估损伤的严重程度的基于频率变化的结构损伤定位与损伤程度的评估方法。属于海洋石油工程领域。

背景技术

海洋平台等大型工程结构在其服役期间损伤是不可避免的，而确保结构安全的唯一方法是尽早确定出结构的损伤位置，并对损伤的程度进行评估，以便进行及时的维修加固。

目前，对于结构损伤的检测方法主要包括：目测和超声、磁粉、声发射等局部无损检测方法等。然而，由于较弱的视觉观测条件以及损伤部位有可能被海生物覆盖等原因，阻碍了目测的效果。此外，上述方法中技术要求结构的损伤区域是以已知作为先决条件，要求配备特殊额外的测试设备和专业人员。因此，上述方法对海洋平台的检测不太方便，并且检测费用昂贵。

与上述方法相比较，基于振动测试的结构损伤诊断技术是相对简单、成本较低的，被公认为是较有发展前景的全局性损伤检测方法。这种方法的基本原理是：损伤将导致结构的系统刚度发生改变，因而导致结构的动力特性参数(如结构的频率、振型等)的变化，所以，结构的动力特性参数能够作为结构损伤诊断的指标，用于判定结构是否有损伤发生，并进而评估损伤的严重程度。这类方法最突出的优点是利用环境激励下的动力响应测试进行损伤诊断，整个损伤诊断操作过程不会影响结构的正常工作，能方便地完成结构损伤的在线检测和评估。

结构的损伤诊断包括四个层次：

(1)判断结构是否发生损伤(损伤识别)；

(2)确定结构的损伤位置(损伤定位)；

(3)评估结构的损伤程度(损伤程度评估)；

(4)结构剩余寿命的预测。

基于结构动力特性变化的损伤诊断方法经过多年发展，研究人员提出并发展了许多方法，但大都集中在一、二层次。可以同时进行损伤定位和损伤程度评估的研究成果较少，而且，这些方法在实际应用中，尤其是应用于大型结构存在一定的局限性如下：

(1)需要高阶的模态信息：由于激振荷载的频率范围和识别技术等方面的原因，一般很难得到高阶模态信息；

(2)需要完备的模态信息：由于传感器数量的限制和某些自由度信息(如转动自由度)难以量测，实际测试得到的结构模态参数是不完备的；

(3)需要质量归一化的阵型信息：对环境荷载激励下仅利用输出响应进行模态参数识别，其识别得到的振型无法质量归一化。

如何克服这些缺点，仅利用识别得到的有限的低阶模态信息进行损伤检测，是目前研究的重点和难点。

发明内容

本发明的主要目的在于克服现有技术存在的上述缺点，而提供一种基于频率变化的结构损伤定位与损伤程度的评估方法，其不需要结构的振型信息，仅利用结构损伤前后的有限的低阶频率信息即可准确地进行损伤定位和损伤程度评估；适用于单损伤、多种损伤工况，具有一定的实际应用价值。

本发明的目的是由以下技术方案实现的：

一种基于频率变化的结构损伤定位与损伤程度的评估方法，其特征在于：采用以下步骤：

第一步：建立末损伤结构的有限元模型，作为基准有限元模型；

第二步：对损伤结构进行振动测试，识别损伤结构的m阶模态频率 ${\mathrm{\ω}}_i^{*}(i=1,...,m);$

第三步：对损伤结构进行损伤程度评估；

第四步：对损伤结构进行损伤位置确定。

所述第一步中，基准有限元模型以K和M表示，即：未损伤结构的刚度矩阵和质量矩阵，其特征值表示为

$K{\mathrm{\Φ}}_i={\mathrm{\ω}}_i^{2}M{\mathrm{\Φ}}_i---\left(1\right)$

其中，ω_i和Φ_i分别为基准有限元模型的第i阶频率和振型。

所述第二步中，识别损伤结构的m阶模态频率具体步骤如下：

(1)、利用传感器测量结构物损伤后的结构动力响应数据，并将数据存储入存储器中；

(2)、从存储器中读取存储的参数信息，利用模态参数识别方法识别其有限的低阶模态频率作为损伤结构模态频率

所述损伤结构是以未损伤结构的刚度矩阵和质量矩阵特征值公式：为基准，以K^*和M^*表示损伤结构的整体刚度矩阵和整体质量矩阵，结构损伤一般只引起结构刚度的变化，而对结构的质量影响很小，设定M^*＝M，则其特征值为

$K^{*}{\mathrm{\Φ}}_i^{*}={\mathrm{\ω}}_i^{*2}M{\mathrm{\Φ}}_i^{*}---\left(2\right)$

其中，和分别为损伤结构的第i阶频率和振型；利用公式(1)和(2)，得到

${\left({\mathrm{\Φ}}_i\right)}^T{K}^{*}{\mathrm{\Φ}}_i^{*}=\frac{{\mathrm{\ω}}_i^{*2}}{{\mathrm{\ω}}_i^2}{\left({\mathrm{\Φ}}_i\right)}^{T}K{\mathrm{\Φ}}_i^{*}---\left(3\right)$

其中，上标T代表矩阵的转置，即行列互换；

设定结构有N_d个单元发生了损伤，且损伤位置已知，则损伤结构的整体刚度矩阵表示为

$K^{*}=K+\underset{n=1}{\overset{N_d}{\mathrm{\Σ}}}{\mathrm{\α}}_n{K}_{l_n}---\left(4\right)$

其中，α_n和l_n分别为第n个损伤单元的损伤程度和第n个损伤单元的单元号；将公式(4)代入公式(3)，可得

$\underset{n=1}{\overset{N_d}{\mathrm{\Σ}}}{\mathrm{\α}}_n{\left({\mathrm{\Φ}}_i\right)}^T{K}_{l_n}{\mathrm{\Φ}}_i^{*}=\frac{{\mathrm{\ω}}_i^{*2}-{\mathrm{\ω}}_i^2}{{\mathrm{\ω}}_i^2}{\left({\mathrm{\Φ}}_i\right)}^{T}K{\mathrm{\Φ}}_i^{*}---\left(5\right)$

定义 $C_i={\left({\mathrm{\Φ}}_i\right)}^{T}K{\mathrm{\Φ}}_i^{*},$ $C_{n,i}={\left({\mathrm{\Φ}}_i\right)}^T{K}_{l_n}{\mathrm{\Φ}}_i^{*}$ 和 $b_i=\frac{{\mathrm{\ω}}_i^{*2}-{\mathrm{\ω}}_i^2}{{\mathrm{\ω}}_i^2},$ 则公式(5)可以写为

$\underset{n=1}{\overset{N_d}{\mathrm{\Σ}}}{\mathrm{\α}}_n{c}_{n,i}=b_i---\left(6\right)$

写成矩阵形式，即

Cα＝b    (7)

设定已经测得了损伤结构的前m阶频率信息，则C是维数为m×N_d的矩阵，α和b分别是维数为N_d×1和m×1的列向量；通过最小二乘法求解公式(7)，得到

α＝(C^TC)^-1C^Tb    (8)

α即为结构单元的损伤程度。

所述第三步中，对损伤结构进行损伤程度评估的具体步骤如下：

(1)设定损伤位置已知，利用基准有限元模型和损伤结构的模态频率构建m个方程，即：

$\underset{n=1}{\overset{N_d}{\mathrm{\Σ}}}{\mathrm{\α}}_n{\left({\mathrm{\Φ}}_i\right)}^T{K}_{l_n}{\mathrm{\Φ}}_i^{*}=\frac{{\mathrm{\ω}}_i^{*2}-{\mathrm{\ω}}_i^2}{{\mathrm{\ω}}_i^2}{\left({\mathrm{\Φ}}_i\right)}^{T}K{\mathrm{\Φ}}_i^{*}$

其中，K为基准有限元模型的总体刚度矩阵；和α_n分别为第n个损伤单元的单元刚度矩阵和损伤程度；N_d为损伤单元的个数；Φ_i和分别为基准有限元模型和损伤结构的第i阶振型；ω_i和分别为基准有限元模型和损伤结构的第i阶频率；上标T代表矩阵的转置；

定义 $C_i={\left({\mathrm{\Φ}}_i\right)}^{T}K{\mathrm{\Φ}}_i^{*},$ $C_{n,i}={\left({\mathrm{\Φ}}_i\right)}^T{K}_{l_n}{\mathrm{\Φ}}_i^{*}$ 和 $b_i=\frac{{\mathrm{\ω}}_i^{*2}-{\mathrm{\ω}}_i^2}{{\mathrm{\ω}}_i^2},$ 则上式可以写为矩阵形式Cα＝b；

(2)迭代求解损伤程度；

迭代求解损伤程度的具体计算步骤为：

①迭代初始赋值，k＝0，设定损伤程度为0，即α⁽⁰⁾＝0，上标“0”代表迭代初值；

②迭代开始，k＝k+1，由 $K^{*\left(k\right)}=K+\underset{n=1}{\overset{N_d}{\mathrm{\Σ}}}{\mathrm{\α}}_n^{(k-1)}{K}_{l_n}$ 和公式 $K^{*\left(k\right)}{\mathrm{\Φ}}_i^{*\left(k\right)}={\mathrm{\ω}}_i^{*2}M{\mathrm{\Φ}}_i^{*\left(k\right)}$ 计算

③利用②计算得到的计算矩阵C和向量b，最后，由α＝(C^TC)^-1C^Tb计算α^(k)；

④设定迭代终止条件，如果max{|α^(k)-α^(k-1)|}≤tol成立，则迭代停止，其中，tol为预先设定的容许误差，则α^(k)为损伤程度估计值；否则返回②继续。

所述第四步中，损伤位置确定的具体计算步骤为：

(1)设定共有N_k个损伤工况，利用N_m组频率变化对每一个损伤工况进行损伤程度评估，得到N_m个损伤程度估计值i＝1，…，N_m，下标k表示第k个设定的损伤工况；

(2)利用损伤程度估计值构建损伤定位因子DI_k，其定义为

${\mathrm{DI}}_k=\frac{\min \left(e_k\right)}{e_k}$

上式中 表示对应第k个设定损伤工况的N_m组估计值的均值，σ_k为其标准差，min(e_k)表示对e_k取最小值；则损伤定位因子等于1的单元为真正的损伤位置。

本发明的有益效果是：

(1)由于其不需要结构的振型信息，仅利用频率的变化来进行计算，而频率信息比振型信息更容易获取，因此，避开了实测振型完备性、质量归一化以及识别精度较低的缺点；

(2)仅利用损伤前后有限的低阶频率信息即可进行损伤位置确定和损伤程度评估，因此，频率的识别精度高于振型的识别精度，可以把模态参数识别误差降至最低。

(3)适用于单损伤、多种损伤工况，能够准确地确定结构的损伤位置并评估其损伤程度，具有一定的实际应用价值。

附图说明：

图1为本发明平面钢架结构有限元模型示意图。

图2为本发明平面钢架结构的损伤工况示意图。

图3为本发明平面钢架结构的损伤定位效果图(工况1)。

图4为本发明平面钢架结构的损伤定位效果图(工况2)。

图5为本发明平面钢架结构的损伤定位效果图(工况3)。

图6为本发明平面钢架结构的损伤定位效果图(工况4)。

图7为本发明平面钢架结构的损伤定位效果图(工况5)。

图中主要标号说明：

1.竖直杆单元、2.水平杆件单元、3.斜杆单元、4.水平杆件单元、

5.竖直杆单元、6.水平杆件单元、7.斜杆单元、8.水平杆件单元、

9.竖直杆件单元、10.水平杆件单元11.斜杆单元、12.水平杆件单元、

13.竖直杆单元、14.水平杆件单元、15.斜杆单元、16.水平杆件单元、

17.竖直杆单元、18.水平杆件单元、19.斜杆单元、20.水平杆件单元、

21.竖直杆单元。

具体实施方式

如图1，图2所示，本发明采用以下步骤：

第一步：建立未损伤结构的有限元模型，作为基准有限元模型；

第二步：对损伤结构进行振动测试，识别损伤结构的m阶模态频率 ${\mathrm{\ω}}_i^{*}(i=1,...,m);$

第三步：对损伤结构进行损伤程度评估；

第四步：对损伤结构进行损伤位置确定。

以上步骤具体算法推导如下：

一、以K和M表示基准有限元模型，即：未损伤结构的刚度矩阵和质量矩阵，其特征值表示为

$K{\mathrm{\Φ}}_i={\mathrm{\ω}}_i^{2}M{\mathrm{\Φ}}_i---\left(1\right)$

其中，ω_i和Φ_i分别为基准有限元模型的第i阶频率和振型；

损伤结构是以未损伤结构的刚度矩阵和质量矩阵特征值公式：为基准，以K^*和M^*表示损伤结构的整体刚度矩阵和整体质量矩阵，结构的损伤一般只引起结构刚度的变化，而对结构的质量影响很小，设定M^*＝M，则其特征值问题为

$K^{*}{\mathrm{\Φ}}_i^{*}={\mathrm{\ω}}_i^{*2}M{\mathrm{\Φ}}_i^{*}---\left(2\right)$

其中，和分别为损伤结构的第i阶频率和振型。利用公式(1)和(2)，可以得到

${\left({\mathrm{\Φ}}_i\right)}^T{K}^{*}{\mathrm{\Φ}}_i^{*}=\frac{{\mathrm{\ω}}_i^{*2}}{{\mathrm{\ω}}_i^2}{\left({\mathrm{\Φ}}_i\right)}^{T}K{\mathrm{\Φ}}_i^{*}---\left(3\right)$

其中，上标T代表矩阵的转置，即行列互换；

设定结构有N_d个单元发生了损伤，且损伤位置已知，则损伤结构的整体刚度矩阵可以表示为

$K^{*}=K+\underset{n=1}{\overset{N_d}{\mathrm{\Σ}}}{\mathrm{\α}}_n{K}_{l_n}---\left(4\right)$

其中，α_n和l_n分别为第n个损伤单元的损伤程度和第n个损伤单元的单元号。将公式(4)代入公式(3)，可得

$\underset{n=1}{\overset{N_d}{\mathrm{\Σ}}}{\mathrm{\α}}_n{\left({\mathrm{\Φ}}_i\right)}^T{K}_{l_n}{\mathrm{\Φ}}_i^{*}=\frac{{\mathrm{\ω}}_i^{*2}-{\mathrm{\ω}}_i^2}{{\mathrm{\ω}}_i^2}{\left({\mathrm{\Φ}}_i\right)}^{T}K{\mathrm{\Φ}}_i^{*}---\left(5\right)$

定义 $C_i={\left({\mathrm{\Φ}}_i\right)}^{T}K{\mathrm{\Φ}}_i^{*},$ $C_{n,i}={\left({\mathrm{\Φ}}_i\right)}^T{K}_{l_n}{\mathrm{\Φ}}_i^{*}$ 和 $b_i=\frac{{\mathrm{\ω}}_i^{*2}-{\mathrm{\ω}}_i^2}{{\mathrm{\ω}}_i^2},$ 则公式(5)可以写为

$\underset{n=1}{\overset{N_d}{\mathrm{\Σ}}}{\mathrm{\α}}_n{c}_{n,i}=b_i---\left(6\right)$

写成矩阵形式，即

Cα＝b    (7)

设定已经测得了损伤结构的前m阶频率信息，则C是维数为m×N_d的矩阵，α和b分别是维数为N_d×1和m×1的列向量。通过最小二乘法求解公式(7)，得到

α＝(C^TC)^-1C^Tb    (8)

α即为结构单元的损伤程度。

需要说明的是，利用公式(8)求解结构的损伤程度需要用到损伤结构完备的振型信息，这在实际应用中是不可能的，为了克服该限制条件，采用迭代技术来计算结构的损伤程度，具体步骤如下：

第1步：迭代初始赋值；k＝0，设定损伤程度为0，即α⁽⁰⁾＝0，上标“0”代表迭代初值。

第2步：迭代开始，k＝k+1；由 $K^{*\left(k\right)}=K+\underset{n=1}{\overset{N_d}{\mathrm{\Σ}}}{\mathrm{\α}}_n^{(k-1)}{K}_{l_n}$ 和公式 $K^{*\left(k\right)}{\mathrm{\Φ}}_i^{*\left(k\right)}={\mathrm{\ω}}_i^{*2}M{\mathrm{\Φ}}_i^{*\left(k\right)}$ 计算

第3步：利用第2步计算得到的计算矩阵C和向量b，最后由α＝(C^TC)^-1C^Tb计算α^(k)；

第4步：设定迭代终止条件，如果max{|α^(k)-α^(k-1)|}≤tol成立，则迭代停止，其中tol为预先设定的容许误差；否则返回第2步继续。

在上面的求解中，设定损伤构件的位置是已知的，当损伤构件的位置未知时，必须先寻求一种方法来确定结构的损伤位置。本发明提出采用多组频率的变化来进行损伤定位。即：只有真正的损伤位置才能造成同实际测试完全相符的多组频率的变化。为此，对每个“设定”的损伤位置，采用上述的方法进行损伤程度评估，结果记为i＝1，…，Nm，其中下标k表示第k个设定的损伤位置，上标i表示采用第i组频率变化来进行损伤程度评估，Nm表示可用的频率变化组数。

对真正的损伤位置，利用Nm组频率变化进行损伤程度评估时，理论上得到的估计值应该是完全相等的。但由于计算误差或实测频率误差，估计值不会完全相等，但应该近似相等。于是引入变异系数来表示其差值的大小，定义为

$e_k=\frac{{\mathrm{\σ}}_k}{\stackrel{\‾}{{\mathrm{\α}}_k^i}},i=1,\·\·\·,\mathrm{Nm}$

上式中，表示对应第k个“设定”损伤工况的Nm组估计值的均值；σ_k为其标准差。于是e_k中最小值对应的单元为损伤单元，定义损伤定位因子

${\mathrm{DI}}_k=\frac{\min \left(e_k\right)}{e_k}$

其中，min(e_k)表示对e_k取最小值；则损伤定位因子等于1的单元为真正的损伤位置。损伤位置确定后，损伤的严重程度可以通过前述的方法来确定。

上述第二步中，识别损伤结构的m阶模态频率具体步骤如下：

(1)、利用传感器测量结构物损伤后的结构动力响应数据，并将数据存储入存储器中；

(2)、从存储器中读取存储的参数信息，利用模态参数识别方法识别其有限的低阶模态频率作为损伤结构模态频率

二、建立平面钢架有限元模型

如图1所示，本算例模拟研究的平面钢架结构由竖直杆、水平杆、斜杆等构件组成，共21个单元。其中，竖直杆单元为(1、5、9、13、17、21)，水平杆件单元为(2、6、10、14、18)，斜杆单元为(3、7、11、15、19)。利用MATLAB软件编写有限元程序，通过计算机建立有限元模型，作为未损伤结构的基准有限元模型。然后，再模拟不同的损伤工况，得出模拟实测的低阶模态频率。本算例模拟了五种损伤工况，包括不同位置的单个构件损伤、两个构件损伤以及不同程度的损伤。损伤位置(如图2所示)。

三、损伤定位与损伤程度评估分析

基于前三阶实测频率，利用本发明的方法进行损伤定位与损伤程度评估，各工况说明如下：

损伤工况一，水平杆件单元6，单元刚度损失20％，用本发明方法可以准确地进行损伤定位，损伤定位因子(如图3所示)，准确地指出损伤构件为水平杆件单元6，损伤程度估计值为20.007％。

损伤工况二，竖直杆件单元9，单元刚度损失10％，用本发明方法可以准确地进行损伤定位，损伤定位因子(如图4所示)，准确地指出损伤构件为竖直杆件单元9，损伤程度估计值为10.005％。

损伤工况三，斜杆单元11，单元刚度损失15％，用本发明方法可以准确地进行损伤定位，损伤定位因子(如图5所示)，准确地指出损伤构件为斜杆单元11，损伤程度估计值为15.038％。

损伤工况四，水平杆件单元6、斜杆单元11同时发生损伤，其刚度分别损失20％和10％，用本发明方法可以准确地进行损伤定位，损伤定位因子(如图5所示)，其x轴准确地指出损伤构件为水平杆件单元6，y轴准确地指出损伤构件为斜杆单元11；损伤水平杆件单元6和损伤斜杆单元11的损伤程度估计值分别为19.992％和10.004％。验证了该方法对发生两处损伤，即：水平杆件和斜杆，可以很好的进行损伤定位和损伤程度评估。

损伤工况五，斜杆单元7与竖直杆单元13同时发生损伤，其刚度损失分别为20％和10％，用本发明方法可以准确地进行损伤定位，损伤定位因子(如图6所示)，其x轴准确地指出损伤构件为斜杆单元7，y轴准确地指出损伤构件为竖直杆单元13；损伤斜杆单元7和损伤竖直杆单元13的损伤程度估计值分别为19.992％和10.004％。验证了该方法对发生两处损伤(斜杆单元7和竖直杆单元13)可以很好的进行损伤定位和损伤程度评估。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (1)

1.一种基于频率变化的结构损伤定位与损伤程度的评估方法，其特征在于：采用以下步骤：

第一步：建立未损伤结构的有限元模型，作为基准有限元模型；基准有限元模型以K和M表示，即：未损伤结构的刚度矩阵和质量矩阵，其特征值表示为

$K{\mathrm{\Φ}}_i={\mathrm{\ω}}_i^2{\mathrm{M\Φ}}_i---\left(1\right)$

其中，ω_i和Φ_i分别为基准有限元模型的第i阶频率和振型；

第二步：对损伤结构进行振动测试，识别损伤结构的m阶模态频率识别损伤结构的m阶模态频率具体步骤如下：

(1)、利用传感器测量结构物损伤后的结构动力响应数据，并将数据存储入存储器中；

(2)、从存储器中读取存储的参数信息，利用模态参数识别方法识别其有限的低阶模态频率作为损伤结构模态频率

第三步：对损伤结构进行损伤程度评估；对损伤结构进行损伤程度评估的具体步骤如下：

⑴设定损伤位置已知，利用基准有限元模型和损伤结构的模态频率构建m个方程，即：

$\underset{n=1}{\overset{N_d}{\mathrm{\Σ}}}{\mathrm{\α}}_n{\left({\mathrm{\Φ}}_i\right)}^T{K}_{l_n}{\mathrm{\Φ}}_i^{*}=\frac{{\mathrm{\ω}}_i^{{*}^2}-{\mathrm{\ω}}_i^2}{{\mathrm{\ω}}_i^2}{\left({\mathrm{\Φ}}_i\right)}^T{\mathrm{K\Φ}}_i^{*}$

其中，K为未损伤结构的刚度矩阵；和α_n分别为第n个损伤单元的单元刚度矩阵和损伤程度；N_d为损伤单元的个数；Φ_i和分别为基准有限元模型和损伤结构的第i阶振型；ω_i和分别为基准有限元模型和损伤结构的第i阶频率；上标T代表矩阵的转置；

定义 $C_i={\left({\mathrm{\Φ}}_i\right)}^T{\mathrm{K\Φ}}_i^{*},C_{n,i}={\left({\mathrm{\Φ}}_i\right)}^T{K}_{l_n}{\mathrm{\Φ}}_i^{*}$ 和 $b_i=\frac{{\mathrm{\ω}}_i^{{*}^2}-{\mathrm{\ω}}_i^2}{{\mathrm{\ω}}_i^2},$ 则上式可以写为矩阵形式：Cα＝b           (7)；

设定已经测得了损伤结构的前m阶频率信息，则C是维数为m×N_d的矩阵，α和b分别是维数为N_d×1和m×1的列向量；通过最小二乘法求解公式(7)，得到

α＝(C^TC)^-1C^Tb           (8)

α即为结构单元的损伤程度；

⑵迭代求解损伤程度；

迭代求解损伤程度的具体计算步骤为：

①迭代初始赋值，k＝0，设定损伤程度为0，即α⁽⁰⁾＝0，上标“0”代表迭代初值；

②迭代开始，k＝k+1，由 $K^{*\left(k\right)}=K+\underset{n=1}{\overset{N_d}{\mathrm{\Σ}}}{\mathrm{\α}}_n^{(k-1)}{K}_{l_n}$ 和公式 $K^{*\left(k\right)}{\mathrm{\Φ}}_i^{*\left(k\right)}={\mathrm{\ω}}_i^{*2}{\mathrm{M\Φ}}_i^{*\left(k\right)}$ 计算

③利用②计算得到的计算矩阵C和向量b，最后，由α＝(C^TC)^-1C^Tb计算α^(k)；

④设定迭代终止条件，如果max{|α^(k)-α^(k-1)|}≤tol成立，则迭代停止，其中，tol为预先设定的容许误差，则α^(k)为损伤程度估计值；否则返回②继续；

第四步：对损伤结构进行损伤位置确定，损伤位置确定的具体计算步骤为：

⑴设定共有N_k个损伤工况，利用N_m组频率变化对每一个损伤工况进行损伤程度评估，得到N_m个损伤程度估计值i＝1,…,N_m，下标k表示第k个设定的损伤工况；

⑵利用损伤程度估计值构建损伤定位因子DI_k，其定义为

${\mathrm{DI}}_k=\frac{\min \left(e_k\right)}{e_k}$

上式中表示对应第k个设定损伤工况的N_m组估计值的均值，σ_k为其标准差，min(e_k)表示对e_k取最小值；则损伤定位因子等于1的单元为真正的损伤位置。