基于广义柔度矩阵的井架钢结构损伤识别

摘要

钻机井架结构在生产过程中会因各种原因而出现不同程度的损伤及破坏，为解决复杂工作环境下石油钻井井架结构健康监测中损伤识别问题，提出了一种基于广义柔度矩阵的井架结构损伤识别方法。相比于柔度矩阵方法，广义柔度矩阵仅需前几阶低阶固有频率及相应的模态振型便可确保计算的精度。在构建损伤控制方程的过程中，采用矩阵分块模式求解方法，与传统方法中需将单元刚度矩阵扩充整体刚度矩阵维数相比较，所提出的求解方法只需针对单元刚度矩阵的行数进行相应扩充，而列数不变，因此从计算量的统计角度分析可大大降低构建损伤识别控制方程的计算量。最后，以实际工程中的K型井架钢结构为例，结果显示对于不同的损伤位置及程度，该方法均具有良好的识别效果，可用于其他大型结构的损伤诊断或健康监测。

关键词

K型井架钢结构; 广义柔度矩阵; 损伤识别; 健康监测; 矩阵分块模式; 石油钻井

0 引言

石油钻井井架结构作为一种大型承重结构，在工程实际中起着重要的作用，尤其对于超深科学钻探的井架结构在长期服役过程中，由于长期作业活动和恶劣工作环境等因素的影响^［1-2］，常常会出现不同程度的损伤例如裂纹、变形和腐蚀等，进而降低井架结构承载能力，也同时影响着钻机的整体使用寿命^［3-4］，甚至出现累积损坏导致的灾难性事故以及巨大损失，因此早期对井架结构进行健康检测分析显得十分重要，并在近些年受到了学者的广泛关注。

近些年来，结构健康监测系统的研究大多是基于在线结构模态参数识别、有限元模型修正与仿真计算为核心的技术手段，并结合测试与传感技术、网络通信技术、信号处理技术等多领域，对工程结构的健康状态实现在线监测与评估^［5］。结构的健康状态通常可用特征参数进行描述，而特征参数通常包括结构的模态参数和物理参数。模态参数为结构的固有频率、阻尼比和振型，它们是结构最基本的动力学特性。当模态参数发生异常变动时，即判定结构处于损伤状态。因此发展一种快速、准确的结构模态特征参数识别技术至关重要。利用结构模态特征参数能够整体评估工程结构，检测费用低，且无需中断工程结构的正常使用，因此成为国内外结构健康监测领域的研究焦点。

基于结构模态参数的损伤识别方法充分利用了结构损伤前后动力特性的差异，即考虑结构的固有频率、模态振型及其曲率、模态柔度及其导数、模态应变能和频响函数在发生损伤前后的变化。这类方法可以对工程结构进行全面的损坏评估，并能够检测出结构内部是否发生损坏。值得注意的是，工程结构对应的高阶模态参数的测量难度远高于低阶模态参数，且实际测量过程中结构低阶模态参数的测量精度往往更容易控制。基于上述问题，Pandey等^［6］首先提出了一种利用结构损伤前后其柔度矩阵改变量的方法来识别结构的损伤程度。柔度矩阵指的是结构刚度矩阵的逆矩阵，ZHAO等^［7］对柔度矩阵进行了有关固有频率和模态振型的灵敏度分析，也证明柔度矩阵相比于刚度矩阵对于低阶模态参数的反应更加灵敏。

基于柔度矩阵的识别方法，Ghosh等^［8］利用结构的高阶振型及其导数进行损伤定位，并以三维钢制抗弯矩框架模型为例进行探究，结果表明虽然利用结构的高阶振型及其导数能够有效地定位损伤，但它们的效率可能会显著降低，特别是对于较矮的建筑物中损伤的定位。YANG等^［9］通过将柔度矩阵进行分解，提出一种新的损伤检测方法，并通过算例进行验证，结果表明该法无需进行高阶灵敏度分析或迭代，即可准确地计算出刚度扰动参数。LI等^［10］通过广义柔度矩阵检测结构损伤，将广义柔度矩阵函数做关于损伤参数等于零处的一阶泰勒展开，再利用损伤前后广义柔度矩阵之差作为右端项，最终构建损伤控制方程并利用最小二乘法进行求解，该方法相比于柔度矩阵，所用的模态数据较少且识别精度更高。LIU等^［11］提出了改进的广义柔度矩阵方法用于确定损伤部位及程度，推导出了在无噪声和有噪声两种情况下，考虑了损伤程度的非负性的求解方法。Aghaeidoost等^［12］利用广义柔度矩阵和最优遗传算法构建了一种结构损伤识别方法，该方法基于结构在完整和损伤状态的模态信息建立目标函数，并结合特定的约束条件进行优化，最后通过以导管架式海洋平台为例，结果表明了在不同损伤级别下能够准确识别损伤位置和损伤程度。LIU等^［13］基于广义柔度矩阵提出了一种处理具有不完整模态振型数据的损伤检测方法。Katebi等^［14］针对平面桁架和框架结构提出了一种改进的模态灵敏度损伤识别方法，该方法利用灵敏度矩阵的变化和用实测的柔度数据，能够准确判断损伤程度及位置。谢少鹏等^［15］提出一种损伤定位方法，该方法将基于广义模态柔度损伤识别问题转化为最小二乘问题，并且利用最小平方正交三角分解（LSQR）法进行求解。然而上述研究的数值案例均是基于简单的力学模型结构，并未应用于实际大型复杂的工程结构。

本文在广义柔度矩阵方法的基础上，优化了控制方程的构建，在求解控制方程中采取分块计算的模式，并以工程实际大型K型井架结构作为数值模型，验证方法的有效性和工程适用性，结果表明该改进的广义柔度矩阵法在保证计算效率和计算精度的同时，可显著减少运算所需的计算量。

1 广义柔度矩阵灵敏度方法

由于工程结构的模态特征参数如频率、振型等是本身具备的物理特性，即均是关于质量、刚度等的函数，因此一旦工程结构发生损伤，其物理特性就会发生改变，因此利用结构模态特性的这一特点，即可实现判断结构是否发生损伤、损伤位置以及评估损伤程度。而基于结构模态特性进行结构损伤识别，在理论上就是将结构的模态特性关于结构设计变量或参数进行泰勒级数展开，从而获取结构模态特性与相关物理参数之间的变化关系^［16］，也就是对结构的模态特性进行灵敏度分析。因此通过结构的灵敏度分析即可实现对损伤位置及程度的评估，得到了广泛地应用^［17-19］。

由于结构的自身刚度作为反应结构强度的典型物理量，因此本文假定工程结构的损伤仅由结构本身刚度的折减所导致，因此结构的质量矩阵和自由度数目n可视为保持不变。在此前提下，工程结构的损伤问题可以下面数学模型表明^［13］：

$$K_{\mathrm{d}}=K_{\mathrm{u}}-{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{e}}}}{\alpha }_i{K}_{\mathrm{u}i}$$

(1)

式中：$K_{\mathrm{u}}$、$K_{\mathrm{d}}$——分别为结构未损伤及发生损伤后的整体刚度矩阵，均是$n\times n$维实对称矩阵；$K_{\mathrm{u}i}$——扩充为结构总体刚度矩阵形式后的单元刚度矩阵；${\alpha }_i$——第$i$个单元发生损伤的程度；$N_{\mathrm{e}}$——结构划分的单元个数。

将式（1）两端同时关于损伤程度${\alpha }_i$求导可得：

$$\frac{\partial K_{\mathrm{d}}}{\partial {\alpha }_i}=-K_{\mathrm{u}i}$$

(2)

将结构发生损伤后的$n\times n$维柔度矩阵表示为$F_{\mathrm{u}}$，结构未发生损伤的$n\times n$维柔度矩阵表示为$F_{\mathrm{d}}$，因此满足如下关系式^［17］：

$$\left\{\begin{array}{c}{F}_{\mathrm{d}}=\underset{}{{\displaystyle \sum _{j=\mathrm{1}}^k}\frac{\mathrm{1}}{{\lambda }_{\mathrm{d}j}}{\phi }_{\mathrm{d}j}{\phi }_{\mathrm{d}j}^T}\\ F_{\mathrm{u}}=\underset{}{{\displaystyle \sum _{j=\mathrm{1}}^k}\frac{\mathrm{1}}{{\lambda }_{\mathrm{u}j}}{\phi }_{\mathrm{u}j}{\phi }_{\mathrm{u}j}^T}\end{array}\right.$$

(3)

式中：$\left({\phi }_{\mathrm{d}j},{\lambda }_{dj}\right)$、$\left({\phi }_{\mathrm{u}j},{\lambda }_{\mathrm{u}j}\right)$——分别为结构损伤以及未损伤情况时的第$j$阶特征向量与特征值；$k$——截断的模态个数。

且$F_{\mathrm{u}}$和$F_{\mathrm{d}}$同时满足：

$$F_{\mathrm{d}}{K}_{\mathrm{d}}=I$$

(4)

$$F_{\mathrm{u}}{K}_{\mathrm{u}}=I$$

(5)

对式（4）的两端同时关于损伤程度${\alpha }_i$求导，并结合式（2）可得：

$$\frac{\partial F_{\mathrm{d}}}{\partial {\alpha }_i}=F_{\mathrm{d}}{K}_{\mathrm{u}i}{F}_{\mathrm{d}}$$

(6)

利用结构损伤后的模态满足关于质量矩阵归一化的条件，一般广义柔度矩阵^［20］可表示为：

$$f_{\mathrm{d}}^g\left(\alpha \right)=F_{\mathrm{d}}(M{F}_{\mathrm{d}}{)}^L={\Phi }_{\mathrm{d}}{\Lambda }_{\mathrm{d}}^{-L-\mathrm{1}}{\Phi }_{\mathrm{d}}^T$$

(7)

其中，$M$为$n\times n$维质量矩阵，本文中假设在损伤前后未发生变化；$\alpha ={\left({\alpha }_{\mathrm{1}},{\alpha }_{\mathrm{2}},\cdots ,{\alpha }_{Ne}\right)}^T$是由所有单元损伤程度组成的一个向量。${\Phi }_{\mathrm{d}}=\left({\phi }_{\mathrm{d}\mathrm{1}},{\phi }_{\mathrm{d}\mathrm{2}},\cdots ,{\phi }_{\mathrm{d}k}\right)$、${\Lambda }_{\mathrm{d}}=diag\left({\lambda }_{\mathrm{d}\mathrm{1}},{\lambda }_{\mathrm{d}\mathrm{2}},\cdots ,{\lambda }_{\mathrm{d}k}\right)$分别代表结构损伤后的前$k$个低阶特征向量矩阵和特征值对角矩阵。由式（7）可以得出，当$L$值越大，高阶模态在$f_{\mathrm{d}}^g\left(\alpha \right)$中所占比例越小。值得注意的是，当$L=\mathrm{0}$时广义柔度矩阵即可表示成一般柔度矩阵，在本文中，仅考虑$L=\mathrm{1}$的情况，如下式所示：

$$f_{\mathrm{d}}^g\left(\alpha \right)=F_{\mathrm{d}}M{F}_{\mathrm{d}}={\Phi }_{\mathrm{d}}{\Lambda }_{\mathrm{d}}^{-\mathrm{2}}{\Phi }_{\mathrm{d}}^T$$

(8)

将式（8）两端同时对损伤参数${\alpha }_i$求导可得：

$$\frac{\partial f_{\mathrm{d}}^g}{\partial {\alpha }_i}=\frac{\partial F_{\mathrm{d}}}{\partial {\alpha }_i}M{F}_{\mathrm{d}}+F_{\mathrm{d}}M\frac{\partial F_{\mathrm{d}}}{\partial {\alpha }_i}$$

(9)

将式（6）代入式（9）中得到损伤程度${\alpha }_i=\mathrm{0}$时的广义柔度矩阵灵敏度：

$${\left.\frac{\partial f_{\mathrm{d}}^g}{\partial {\alpha }_i}\right|}_{{\alpha }_i=\mathrm{0}}=F_{\mathrm{u}}{K}_{\mathrm{u}i}{F}_{\mathrm{u}}M{F}_{\mathrm{u}}+F_{\mathrm{u}}M{F}_{\mathrm{u}}{K}_{\mathrm{u}i}{F}_u$$

(10)

同时，将式（8）中广义柔度矩阵在损伤程度${\alpha }_i=\mathrm{0}$处进行泰勒级数展开，并仅保留一阶项，可得：

$$f_{\mathrm{d}}^g\approx f_{\mathrm{u}}^g+{\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{e}}}}{\alpha }_i{\left.\frac{\partial f_{\mathrm{d}}^g}{\partial {\alpha }_i}\right|}_{{\alpha }_i=\mathrm{0}}$$

(11)

上式中$f_{\mathrm{u}}^g$对应结构未发生损伤时的广义柔度矩阵。将式（10）代入（11）并结合式（8）可得最终的损伤识别控制方程，即：

$${\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{e}}}}{\alpha }_i\left(F_{\mathrm{u}}{K}_{\mathrm{u}i}{F}_{\mathrm{u}}M{F}_{\mathrm{u}}+F_{\mathrm{u}}M{F}_{\mathrm{u}}{K}_{\mathrm{u}i}{F}_{\mathrm{u}}\right)=\mathrm{\Delta }f$$

(12)

其中：

$$\mathrm{\Delta }f=f_{\mathrm{d}}^g-f_{\mathrm{u}}^g={\displaystyle \sum _{j=\mathrm{1}}^k}\frac{\mathrm{1}}{{\lambda }_{\mathrm{d}j}^{\mathrm{2}}}{\phi }_{\mathrm{d}j}{\phi }_{\mathrm{d}j}^T-{\displaystyle \sum _{j=\mathrm{1}}^k}\frac{\mathrm{1}}{{\lambda }_{\mathrm{u}j}^{\mathrm{2}}}{\phi }_{\mathrm{u}j}{\phi }_{\mathrm{u}j}^T$$

(13)

上式中$\left({\lambda }_{\mathrm{d}j},{\phi }_{\mathrm{d}j}\right)$和$\left({\lambda }_{\mathrm{u}j},{\phi }_{\mathrm{u}j}\right)$分别代表损伤结构与未损伤结构的第j个低阶特征值和特征向量。

在构建方程（12）的过程中需要大量的数据计算，很大程度上降低了方法的求解效率。因此本文在式（12）的基础上，利用矩阵分块计算的方法，实现了计算效率上的提升。

2 改进的损伤识别控制方程构建及求解

为提高式（12）的计算效率，首先利用式（5）这一关系，可将式（12）左右两端均乘以矩阵，进行如下变形：

$${\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{e}}}}{\alpha }_i\left(K_{\mathrm{u}i}{F}_{\mathrm{u}}M+M{F}_{\mathrm{u}}{K}_{\mathrm{u}i}\right)=K_{\mathrm{u}}\mathrm{\Delta }f{K}_{\mathrm{u}}$$

(14)

令$\mathrm{\Delta }F=K_{\mathrm{u}}\mathrm{\Delta }f{K}_{\mathrm{u}}$，$C_i=M{F}_{\mathrm{u}}{K}_{\mathrm{u}i}$，即式（14）可写成：

$${\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{e}}}}{\alpha }_i\left(C_i^T+C_i\right)=\mathrm{\Delta }F$$

(15)

将$C_i$中的K_ui做如下分块处理：

$$K_{\mathrm{u}i}=\left[\mathrm{0},{\tilde{K}}_{\mathrm{u}i},\mathrm{0}\right]$$

(16)

其中${\tilde{K}}_{\mathrm{u}i}$表示由第$i$个单元刚度矩阵进行扩充后的$n\times m$维矩阵，$m$代表单元刚度矩阵的维数。对具有约束的工程结构而言，其有限元模型结构对应的刚度矩阵通常可视为对称非奇异矩阵，因此结合式（16），也可将$C_i$进行如下分块：

$$C_i=\left[\mathrm{0},M{F}_{\mathrm{u}}{\tilde{K}}_{\mathrm{u}i},\mathrm{0}\right]$$

(17)

在采取分块模式求解$C_i$的过程中，由于${\tilde{K}}_{\mathrm{u}i}$的列数$m$远小于整体刚度矩阵的维数$n$（即$n-m\approx n$），因此相比于直接求解$C_i$，经分块处理后所需要的计算量减少了大约$\mathrm{2}{n}^{\mathrm{3}}$次运算量，而整个控制方程需要求解$N_e$个系数项$C_i$，且网格划分越密集，相应的计算量也将大大减少，具体计算量如表1所示。

表1  系数矩阵$C_i$计算量分析统计

Table 1  Analysis and statistics of the computation amount of coefficient matrix $C_i$

计算系数矩阵$C_i$	计算步骤	计算量	总的计 算量
直接法	对$K_{\mathrm{u}}$进行Choloschy分解	$\frac{n^3}{6}$	$\frac{13n^3}{6}$
	由$R=K_{\mathrm{u}}^{-1}{K}_{\mathrm{u}i}$回带求解矩阵$R$	$n^3$
	计算$C_i=MR$	$n^3$
分块求解方法	对$K_{\mathrm{u}}$进行Choloschy分解	$\frac{n^3}{6}$	$\frac{n^3}{6}+2m{n}^2$
	由${\displaystyle \tilde{R}}=K_{\mathrm{u}}^{-1}{{\displaystyle \tilde{K}}}_{\mathrm{u}i}$回带求解矩阵${\displaystyle \tilde{R}}$	$m{n}^2$
	计算$C_i=\left[0,M{\displaystyle \tilde{R}},0\right]$	$m{n}^2$

最后，基于分块模式求解矩阵$C_i$后，式（13）可写成如下形式进行求解：

$$A\alpha =b$$

(18)

上式中$\alpha =\left({\alpha }_{\mathrm{1}},{\alpha }_{\mathrm{2}},\cdots ,{\alpha }_{N_{\mathrm{e}}}\right)$表示各单元损伤程度的$N_{\mathrm{e}}$维列向量，$b=vec\left(∆F\right)$表示将$n\times n$维的$∆F$矩阵拉直成一个$n^{\mathrm{2}}$维的列向量，同理$A=\left[vec\left(C_{\mathrm{1}}\right),vec\left(C_{\mathrm{2}}\right),\cdots ,vec\left(C_{N_{\mathrm{e}}}\right)\right]$是$C_i\left(i=\mathrm{1,2},\cdots ,N_{\mathrm{e}}\right)$依次进行拉直处理后由左向右排列的$n^{\mathrm{2}}\times N_{\mathrm{e}}$维矩阵，最后利用最小二乘法求解此方程组即可确定结构的损伤位置及损伤程度。

3 数值实验

本节中，仅考虑由于弹性模量的减少或出现裂纹导致的横截面积减少引起的损伤情况，因此可统一归结为结构有限元模型的刚度折减损伤情况，并且当单元损伤程度<5%时，则认为单元未发生损伤破坏^［20］，而质量矩阵在结构损伤前后不会发生改变。所有运算程序均使用Intel（R） Fortran Composer XE 2019，并基于Math Kernel Library函数数据库编译并计算。为验证本文所提出的损伤识别方法在实际工程结构中的有效性与适用性，故选用某科研钻探工程K型井架实体结构建立相应的有限元模型，井架整体结构高约57 m，宽约10 m，实体结构见图1（a）。由于石油井架在钻采过程中会伴随自身结构的前后摆动，且变形较大处大多发生在结构的中下部位置，故假设在该位置处随机产生多种损伤情况（见表2），并以模拟计算结果与之比对，以验证算法的有效性。

图1  井架实体结构及有限元模型

Fig.1  Entity structure and finite element model of

drilling rig derrick structure

表2  井架结构损伤情况

Table 2  Structural damage to the derrick

损伤情形	损伤单元编号	损伤程度/%
Ⅰ多重损伤	104 106 145	25 25 20
Ⅱ双重损伤	159 165	25 25
Ⅲ单重损伤	160	25

考虑如图1（b）所示井架有限元模型，该模型共有219个单元，每个单元具有2个节点，共计100个节点，每个节点具有6个自由度，井架模型结构最底层的4个节点自由度完全约束，模型总的自由度数为$n=\mathrm{576}$。在该数值算例中，单元矩阵维数$m=\mathrm{12}$，故根据表1的浮点计算量统计可知，通过直接法计算单个系数矩阵$C_i$的计算量约为$\mathrm{4.14}\times {\mathrm{10}}^{\mathrm{8}}$，而利用本文中所提出的矩阵分块技术计算单个系数矩阵$C_i$的计算量约为3.98×10⁷，是直接法计算量的10%，因此大大提高计算效率。材料几何参数如下：弹性模量$E=\mathrm{210}\mathrm{G}\mathrm{P}\mathrm{a}$，密度$\rho =\mathrm{7850}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}$，泊松比$\upsilon =\mathrm{0.288}$。

图2所示是假定多重损伤Ⅰ对应的识别情况。从图中可看出，假设损伤部位出现在中下部左侧部位，单元编号分别为104、106以及145。在该工况当中，广义柔度矩阵损伤识别方法仅考虑前两阶频率和相应的振型，计算结果如图3所示，从图中可见，损伤程度较为突出的位置发生在104号单元、106号单元以及145号单元，其对应的损伤程度分别为0.23、0.23、0.16，可判定以上3个单元为损伤部位。损伤程度与假定损伤情况${\alpha }_{\mathrm{104}}=\mathrm{0.25}$，${\alpha }_{\mathrm{106}}=\mathrm{0.25}$以及${\alpha }_{\mathrm{145}}=\mathrm{0.20}$相比，较为接近，可见，所提出的方法能够准确地判断出结构发生破坏的位置，在实际复杂工程结构中应用具有适用性及有效性。

图2  井架结构损伤工况Ⅰ

Fig.2  Damage scenario Ⅰ for derrick structure

图3  井架结构多重损伤情形Ⅰ

Fig.3  Derrick structure multiple damage scenario Ⅰ

双重损伤工况的识别结果如图4所示。在该工况当中，假设损伤部位出现在中下部两侧部位，损伤单元编号为159以及165，假定的损伤程度均为0.25。在该工况中，基于广义柔度矩阵损伤识别方法仅考虑前两阶频率和相应的振型，计算结果如图5所示，可以确定出159、165两个单元对应的位置发生了较大的损伤情况，且损伤程度分别为0.2475和0.2454，与预先假定的损伤工况相比，损伤位置判断准确，损伤程度计算接近于预先设定的损伤程度。在一些其他单元出现了轻微损伤的误判，如单元编号为131和158位置，但损伤程度仅为0.0786和0.0807。虽然有误判，但是已损伤的单元均已突出，因此可以认为对结构的损伤部位实现了较为准确的损伤判断。

图4  井架结构损伤预设Ⅱ

Fig.4  Damage preset Ⅱ for derrick structure

图5  井架结构损伤工况Ⅱ

Fig.5  Damage scenario Ⅱ of derrick structure

图6所示为单个损伤的识别情况，由图中可以看出，预先假定的损伤单元出现在下部后方的斜直梁，单元编号为160，损伤程度设为${\alpha }_{\mathrm{160}}=\mathrm{0.25}$。由于结构损伤后相比结构损伤前仅发生单个单元的损伤情况，结构改动较小，因此为了精确捕捉损伤单元及其损伤程度，在识别方法当中，考虑前四阶频率和相应的振型，计算结果如图7所示，显示损伤的部位同样为第160号单元，损伤程度为0.2864，与假定工况相比，对损伤位置的定位精确，损伤程度识别接近预设损伤数值。

图6  井架结构损伤预设Ⅲ

Fig.6  Damage scenario Ⅲ for derrick structure

图7  井架结构单重损伤情形Ⅲ

Fig.7  Single damage situation Ⅲ of derrick structure

基于上述的数值算例结果可知，本文中所提出的基于广义柔度矩阵的改进方法，在处理实际的大型复杂K型井架结构时，能够准确判断出发生损伤破坏的位置，并具备一定的损伤程度评估功能，且由于对应有限元模型单元划分较多，整体控制方程组维数较大，因此能显著降低计算机的浮点运算次数，提高求解效率。因此该方法同样具备实际大型复杂工程结构健康预测及损伤判断的潜力。

4 结论

（1）基于广义柔度矩阵法，进行结构损伤识别研究，在该方法当中，仅需要利用结构前几个低阶频率以及对应的振型，而且在计算过程中采用分块求解模式，在确保损伤识别精度的同时减少了运算的计算量。

（2）以某科研钻探工程K型井架实体结构为依据，构建与之对应的有限元模型结构，并随机针对多种不同的结构损伤工况进行数值实验，结果表明，该方法能够较为准确地识别损伤位置及相应的损伤程度，具有一定的可靠性。