基于时域残余力向量特征分解的结构损伤识别方法研究

曹震雄，翁 顺，李佳靖，陈志丹，于 虹，闫俊锋，余兴胜; CAO Zhen-xiong; WENG Shun; LI Jia-jing; CHEN Zhi-dan; YU Hong; YAN Jun-feng; YU Xing-sheng

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

基于时域残余力向量特征分解的结构损伤识别方法研究 PDF

- ORCID：
曹震雄 ¹
✉
- ORCID：
翁顺 ¹
✉
- ORCID：
李佳靖 ¹
- ORCID：
陈志丹 ¹
- ORCID：
于虹 ²
- ORCID：
闫俊锋 ²
- ORCID：
余兴胜 ²

1. 华中科技大学土木与水利工程学院, 湖北武汉 430074； 2. 中铁第四勘察设计院集团有限公司, 湖北武汉 430063

中图分类号： TU312⁺.3

最近更新：2023-12-15

DOI：10.16385/j.cnki.issn.1004-4523.2023.06.006

摘要

提出基于时域残余力向量的结构损伤识别方法。从结构振动方程出发，推导强迫振动条件下由结构动力响应求解时域残余力向量的表达式。时域残余力向量包含与损伤过程相关的时间信息，可以用来识别损伤发生的时刻。通过结点时域残余力向量值的变异系数建立损伤位置指标，识别损伤位置。建立结构残余力矩阵特征分解方法，特征分解获得的特征向量在损伤前后不发生改变，而特征值与单元刚度相关。因此，可以从特征值中提取单元刚度参数构建损伤程度参数，直接识别损伤程度。将时域残余力向量方法应用于武汉军山大桥模型的损伤识别，结果表明：提出的方法在单点激励和多点激励下均可高精度地识别出损伤发生的时刻、位置和程度。提出的方法通过对时域残余力向量特征分解提取损伤时刻、位置和程度指标，计算过程直观且精度高，避免了繁琐的时频变换、优化迭代等计算过程。

关键词

结构健康监测; 结构损伤识别; 残余力向量; 矩阵特征分解

引言

残余力向量是将结构损伤状态下的模态数据代入到结构未损伤状态下的特征方程中所产生的误差。残余力向量法是一种基于结构动力学特性的损伤识别方法，结构的动力学特性会随其物理参数（质量、刚度和阻尼）变化，故而由此构建的损伤识别指标能够用于诊断结构的损伤。残余力向量作为一种损伤识别指标有着概念明确、计算简单的优点，能够反映结构的整体和局部信息，因此得到了国内外学者的关注，取得了许多研究成果。

Zimmerman等^［

1］提出了一种名为子空间转角法的损伤检测算法，用于分析结构损伤的位置和程度，后来该算法被称为残余力向量法。将该方法用于弹簧模型和二维桁架的试验中，成功地检测出了损伤，但该方法易受噪声的影响。周先雁等^{［参考文献 2

百度学术}2］将残余力向量法与灵敏度分析相结合，实现对结构损伤的识别。Fukunaga等^{［参考文献 3

百度学术}3］通过计算结构的残余力向量确定损伤的大致位置，依据残余力向量误差最小化确定实际损伤位置并得到结构的损伤程度，通过对堆成叠合板结构的数值模拟，验证了该方法的有效性，不过其识别结果受有限元模型精度和测量数据噪声的影响。Ge等^{［参考文献 4

百度学术}4］在使用残余力向量定位损伤的基础上，采用矩阵凝聚方法和比例损伤模型，通过计算结构损伤前后刚度和质量特性的变化来量化损伤。Amiri等^{［参考文献 5

百度学术}5］使用对角化方法优化模态残余力向量，以避免该方法对框架中间层损伤情况的误判。何伟等^{［参考文献 6

百度学术}6］针对测试模态信息不完备的情况以及自由度缩聚带来的误差影响，提出了一种改进的残余力向量法，将损伤前后结构同阶的残余力向量差值作为新的损伤指标进行损伤识别，并用简支梁模型数值模拟验证了该方法的有效性以及良好的抗噪性。Behtani等^{［参考文献 7

百度学术}7］将残余力向量法应用到复合材料层状梁结构的损伤识别，并用石墨/环氧复合梁的数值模拟验证了该方法的有效性。Wu等^{［参考文献 8

百度学术}8］提出了一种实用的残余力分解方法来识别网壳结构的损伤。张干等^{［参考文献 9

百度学术}9］将残余力向量法应用到高桩码头基桩的损伤识别中，并将残余力绝对值作为损伤指标，拓宽了残余力向量法的应用范围。

传统的残余力向量法都是依据模态参数，通过计算残余力向量判断损伤位置。在此基础上，于德介等^［

10］结合神经网络方法，将结构残余力向量作为网络的输入，使用正交空间格点法对网络的训练样本进行预处理，解决了训练样本在数据空间分布不均匀的问题，从而提高识别精度和收敛速度。袁颖等^{［参考文献 11

百度学术}11］将残余力向量法与改进的遗传算法相结合，基于常规模态分析，将结点的残余力向量构造成遗传算法目标函数，在噪声条件下利用改进的遗传算法进行结构损伤定位和定量研究。Yun等^{［参考文献 12

百度学术}12］提出一种将残余力向量的灵敏度用于子集选择过程的损伤检测方法，首先进行一个子集选择过程来识别最可能损坏的单元，然后采用稳态遗传算法进行优化以获取与损伤单元相关联的参数的值，在识别多损伤工况时表现良好。Saberi^{［参考文献 13

百度学术}13］将残余力向量法与电荷系统搜索算法结合，实现对结构损伤的识别。赵一霖等^{［参考文献 14

百度学术}14］将残余力向量法与树种算法结合，通过残余力向量法找出结构的疑似损伤单元，再利用改进的树种算法反演和识别损伤程度，简支梁数值模拟结果表明该方法有较高的求解精度和较好的鲁棒性。上述文献中使用的方法皆基于动态测试的模态参数，可统称为模态残余力向量，该方法概念明确，相关理论较为成熟，但是测试数据需要进行模态分析等时频变换，且该方法对局部损伤不敏感，不含与损伤相关的时间信息。

杨秋伟等^［

15］提出了一种静力残余力向量法，利用静力测试位移数据来进行损伤评估，精度更高且更易实现，并进一步提出了静力缩聚残余力向量法，使该方法的应用范围更广，通过对桁架和梁结构模型的数值模拟验证了该方法的有效性。但该方法不含与损伤相关的时间信息，且实际工程中难以施加较大静力荷载。与模态参数相比，动力响应数据能够直接测得，不需进行模态分析等复杂操作，相较于静力测试数据更加容易采集。

因此，本文基于动力响应数据，推导强迫振动下时域残余力向量的计算公式，根据其值在时间上的变化确定损伤时刻；根据时域残余力向量的定义，使用其变异系数构建损伤位置指标，判定损伤发生的位置；在此基础上，对结构刚度变化矩阵进行特征分解，定义了损伤程度参数，以此识别结构的损伤程度。利用时域残余力向量法对武汉军山大桥模型进行损伤识别，数值算例表明，该方法能准确识别损伤的时间、位置和程度，且计算过程简单，无需繁琐的时频变换和迭代运算，具备较好的噪声鲁棒性。

1 时域残余力向量法

1.1 时域残余力向量

对于一个多自由度的线性体系，其结构动力平衡方程为：

M \ddot{x} (t) + C \dot{x} (t) + K x (t) = F (t)

（1）

式中 M，C和K分别为结构的质量、阻尼和刚度矩阵；F（t）为外荷载； $x (t)$ ， $\dot{x} (t)$ 和 $\ddot{x} (t)$ 分别为结构的位移、速度和加速度响应。

一般认为结构发生损伤时刚度发生改变，质量不改变，则有如下关系：

M_{d} = M_{u}, K_{d} = K_{u} - Δ K

（2）

式中 $M_{u}$ ， $K_{u}$ 分别为未损结构的质量和刚度矩阵； $M_{d}$ ， $K_{d}$ 分别为损伤结构的质量和刚度矩阵； $Δ K$ 为结构发生损伤前后刚度矩阵的变化。

由式（1）和（2）可得损伤结构的动力平衡方程为：

M_{u} \ddot{x} (t) + C \dot{x} (t) + (K_{u} - Δ K) x (t) = F (t)

（3）

将式（3）进行变量分离可得：

Δ K x (t) + F (t) = M_{u} \ddot{x} (t) + C \dot{x} (t) + K_{u} x (t)

（4）

式（4）左边含有参量 $Δ K$ 和 $x (t)$ 所组成的项以及外加荷载项，其中 $Δ K$ 是与损伤相关的未知量。式（4）右边由参量 $M_{u}$ ， $C$ ， $K_{u}$ 和动力响应 $x (t)$ ， $\dot{x} (t)$ ， $\ddot{x} (t)$ 所组成，其中 $M_{u}$ ， $C$ 和 $K_{u}$ 是已知参量，动力响应可以通过传感器测得，故通过求解式（4）可识别结构损伤状况。

定义结构的时域残余力向量（Time Domain Residual Force Vector，TRF）为：

T R F = Δ K x (t) = M_{u} \ddot{x} (t) + C \dot{x} (t) + K_{u} x (t) - F (t)

（5）

通过前面的分析可知，时域残余力向量是将结构损伤状态下的动力响应代入到结构未损伤状态下的动力平衡方程所产生的误差，故而可以通过结构损伤前的物理参数和损伤后的动力响应求得结构的时域残余力向量。

由式（5）可知，当结构没有损伤时， $Δ K = 0$ ，此时时域残余力向量也为零。结构发生损伤后， $Δ K \neq 0$ ，此时时域残余力向量不为零，且其值与结构的损伤程度相关。因此，可以通过时域残余力向量的突变点判断损伤发生的时刻。

1.2 损伤定位

矩阵TRF中的非零行与损伤位置相关，考虑到模型误差与噪声的影响，TRF中元素离散程度较大的行表明了潜在的损伤位置，为排除结构不同位置结点数据尺度差异的影响，取各结点自由度TRF的变异系数作为损伤定位的指标。实际工程中遇到的噪声多为高斯白噪声，其均值和方差为常数，而噪声影响下未损伤单元相关结点的TRF主要由噪声引起，故而未损伤单元相关结点自由度的TRF变异系数相近，可以通过减去损伤无关结点自由度的TRF变异系数来避免噪声可能造成的损伤位置的误判。为降低噪声影响，定义结构损伤位置指标为：

L o c (d) = \frac{σ_{d}}{μ_{d}} - m e d i a n (\frac{σ}{μ})

（6）

式中 $d$ 为自由度编号； $σ_{d}$ 为 $d$ 自由度对应的时域残余力向量的标准差； $μ_{d}$ 为 $d$ 自由度对应的时域残余力向量的平均值； $m e d i a n (σ / μ)$ 为各结点自由度对应的时域残余力向量的变异系数的中位数。由时域残余力向量的性质可知，损伤位置对应自由度的 $L o c$ 值会远大于其他位置的 $L o c$ 值。

1.3 损伤定量

为了确定结构损伤的程度，通过矩阵分解技术^［

16］，对时域残余力向量中的

Δ K

进行分解，从中提取出刚度损伤参数，即为损伤程度参数，以此来衡量结构的损伤程度。

整体刚度矩阵为单元刚度矩阵经过转换后的叠加：

K = \sum_{i = 1}^{N} T_{i}^{T} K_{i} T_{i}

（7）

式中 $T_{i}$ 为将局部坐标下 $n \times n$ 单元刚度矩阵转换为整体坐标下与自由度对应的 $m \times m$ 单元刚度矩阵的变换矩阵；n为单元自由度个数；m为结构整体自由度个数； $K_{i}$ 表示第i个单元的刚度矩阵；N为结构的单元个数。

对单元刚度矩阵 $K_{i}$ 进行特征值分解，得到：

K_{i} = u_{i} p_{i} u_{i}^{T}

（8）

式中 $u_{i}$ 为第 $i$ 个单元的 $n \times r$ 特征向量矩阵； $p_{i}$ 为第 $i$ 个单元的非零特征值｛ $p$ ｝ $_{i}$ 的对角矩阵；r为矩阵的秩。

以一个12自由度的梁单元为例，可以得到其非零特征值和相应特征向量分别为：

\{p\} = [\begin{matrix} \frac{2 E A}{L} & \frac{2 G J}{L} & \frac{2 E I_{z z}}{L} & \frac{6 E I_{z z} (L^{2} + 4)}{L^{3}} & \frac{2 E I_{y y}}{L} & \frac{6 E I_{y y} (L^{2} + 4)}{L^{3}} \end{matrix}]

（9）

u = [\begin{matrix} \frac{1}{\sqrt[]{2}} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{\sqrt[]{2}}{\sqrt[]{L^{2} + 4}} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{\sqrt[]{2}}{\sqrt[]{L^{2} + 4}} \\ 0 & \frac{1}{\sqrt[]{2}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{\sqrt[]{2}} & - \frac{L}{\sqrt[]{2} \sqrt[]{L^{2} + 4}} \\ 0 & 0 & - \frac{1}{\sqrt[]{2}} & \frac{L}{\sqrt[]{2} \sqrt[]{L^{2} + 4}} & 0 & 0 \\ - \frac{1}{\sqrt[]{2}} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & - \frac{\sqrt[]{2}}{\sqrt[]{L^{2} + 4}} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & - \frac{\sqrt[]{2}}{\sqrt[]{L^{2} + 4}} \\ 0 & - \frac{1}{\sqrt[]{2}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & - \frac{1}{\sqrt[]{2}} & - \frac{L}{\sqrt[]{2} \sqrt[]{L^{2} + 4}} \\ 0 & 0 & \frac{1}{\sqrt[]{2}} & \frac{L}{\sqrt[]{2} \sqrt[]{L^{2} + 4}} & 0 & 0 \end{matrix}]

（10）

式中 E为单元的弹性模量；A为单元的截面面积；G为剪切模量； $I_{z z}$ 和 $I_{y y}$ 为截面惯性矩；L为单元长度。

由式（9）和（10）可知，采用式（8）特征分解后，梁单元刚度矩阵特征向量只与单元长度相关，而一般认为损伤前后单元长度不发生改变，因此单元特征向量矩阵在损伤前后不发生改变，而单元刚度非零特征值则与弹性模量和剪切模量相关，即为刚度参数。

整体刚度矩阵可写为：

K = \sum_{i = 1}^{N} T_{i}^{T} u_{i} p_{i} u_{i}^{T} T_{i}

（11）

定义连接矩阵：

A = [T_{1}^{T} u_{1} T_{2}^{T} u_{2} \dots T_{N}^{T} u_{N}]

（12）

进一步整体刚度矩阵可表示为：

K = A P A^{T}

（13）

式中 P为各单元刚度矩阵非零特征值 ${P}$ 的对角矩阵， ${P}$ 为：

{P} = [{p}_{1} {p}_{2} \dots {p}_{N}]^{T} = [P_{1} P_{2} \dots P_{n r}]^{T}

（14）

其中：

n r = \sum_{i = 1}^{N} r_{i}

（15）

若各单元刚度参数的折减程度为 $α_{i}$ ，则损伤前后的单元刚度参数有如下关系：

P_{i, d} = (1 - α_{i}) P_{i, u}

（16）

将式（16）代入式（13）中，可得如下关系式：

K_{u} = A P_{u} A^{T}, K_{d} = A (1 - α) P_{u} A^{T}

（17）

其中：

α = [\begin{matrix} α_{1} \\ α_{2} \\ ⋱ \\ α_{n r} \end{matrix}], α_{i} \in [0,1]

（18）

将式（17）代入式（2）中，可得：

Δ K = K_{u} - K_{d} = A α P A^{T}

（19）

则时域残余力向量TRF的表达式可写为：

T R F = A α P A^{T} x (t)

（20）

式（20）中只有一个未知量 $α$ ，对上式方程两边进行变换可得：

α P A^{T} x (t) = A^{+} T R F

（21）

式中 $A^{+}$ 为 $A$ 的广义逆矩阵，可通过下式求解：

A^{+} = A^{T} (A A^{T})^{- 1}

（22）

为了方便求解，定义参数：

γ = P A^{T} x (t), χ = A^{+} T R F

（23）

则式（21）可写作：

α γ = χ

（24）

将矩阵展开，可得：

[\begin{matrix} α_{1} \\ α_{2} \\ ⋱ \\ α_{n r} \end{matrix}] \{\begin{matrix} γ_{1} \\ γ_{2} \\ ⋮ \\ γ_{n r} \end{matrix}\} = \{\begin{matrix} χ_{1} \\ χ_{2} \\ ⋮ \\ χ_{n r} \end{matrix}\}

（25）

由上述分析可知，结构动力响应是测试数据，则 $γ_{i}$ 和 $χ_{i}$ 都可以通过式（23）计算得到。用单元对应刚度参数的减少评估单元损伤程度，其表达式为：

α_{i} = χ_{i} / γ_{i}

（26）

式中 $α_{i}$ 即为单元损伤程度参数，其值为单元损伤程度，下标 $i$ 表示单元刚度参数编号。当某单元损伤程度参数为0，说明损伤前后该单元刚度参数相同，即单元未发生损伤；当损伤程度参数为1，说明该单元损伤后刚度参数为0，即完全损坏。

2 数值试验

2.1 武汉军山大桥算例

将上述基于时域残余力向量的损伤识别方法应用于武汉军山大桥模型的损伤识别。军山大桥全长4881.178 m，桥梁长2847 m，引道长2034 m。军山大桥包括主桥、过渡孔桥、引桥、引道等部分，是一座五跨连续半漂浮体系双塔双索面钢箱梁斜拉桥。主桥长为964 m（48 m+204 m+460 m+204 m+48 m），主梁连续长度为964 m，宽为38.8 m，过渡孔桥主孔跨径为56 m，引桥跨径均为30 m。桥梁全景图及有限元模型如图1所示。

图1 武汉军山大桥及其有限元模型

Fig.1 Wuhan Junshan Bridge and its finite element model

根据施工图纸建立该工程结构有限元模型，包含758个单元、611个结点、3634个自由度，模型总质量为 $2.16 \times 10^{9}$ kg。整个模型的单元类型共有“BEAM4”和“LINK8”两种，其中与拉索对应的451~594号单元采用“LINK8”单元，单元结点有x，y，z三个方向的平动自由度。与索塔对应的1~284号单元和747~758号单元采用“BEAM4”单元、与主梁对应的285~450号单元、595~746号单元均采用“BEAM4”单元，单元结点有x，y，z三个方向的平动自由度和绕x轴、y轴、z轴的转动自由度。以单元刚度折减模拟结构的损伤，为模拟该方法在结构受到单点激励和多点激励下的表现，外荷载分别采用简谐荷载和地震荷载两种。设置不同损伤位置、损伤程度的9种损伤工况，工况设置如表1所示。其中，简谐荷载只作用在结构的176号结点处，为单点激励；地震荷载作用在结构整体上，为多点激励。在简谐荷载和地震荷载两种外荷载作用下，工况1和4为单损伤工况，工况2和5为多损伤工况，工况3和6为添加10%噪声条件下的多损伤工况，工况7~9为结构发生较小损伤时的工况。各工况损伤位置如图2所示，16号单元为桥墩处的单元，单元两边的结点分别为175和176号结点，64号单元为桥塔处的单元，单元两边的结点分别为184和185号结点。各损伤工况中，损伤均设置为第10 s发生。采用本文提出的时域残余力向量方法，分别对不同工况下的结构损伤进行识别。

表1 不同激励和损伤程度的工况

Tab.1 Cases with different excitations and damage levels

工况	外荷载	损伤模式
工况1	简谐荷载	16号单元发生30%程度损伤
工况2	简谐荷载	16号单元发生30%程度损伤，64号单元发生20%程度损伤
工况3	简谐荷载	16号单元发生30%程度损伤，64号单元发生20%程度损伤，加速度响应中含10%噪声
工况4	地震荷载	16号单元发生30%程度损伤
工况5	地震荷载	16号单元发生30%程度损伤，64号单元发生20%程度损伤
工况6	地震荷载	16号单元发生30%程度损伤，64号单元发生20%程度损伤，加速度响应中含10%噪声
工况7	地震荷载	16号单元发生3%程度损伤
工况8	地震荷载	16号单元发生3%程度损伤，加速度响应中含10%噪声
工况9	地震荷载	16号单元发生3%程度损伤，加速度响应中含2%噪声

图2 各工况损伤示意图

Fig.2 Damage location of each case

在数值试验中，使用高斯白噪声模拟测量噪声，则测量加速度为计算加速度与测量噪声的和^［

17］：

{\ddot{x}}_{m} = {\ddot{x}}_{c} + β σ ({\ddot{x}}_{c}) N

（27）

式中 ${\ddot{x}}_{m}$ 为测量加速度； ${\ddot{x}}_{c}$ 为计算的加速度； $β$ 为噪声水平； $σ ({\ddot{x}}_{c})$ 为加速度的标准差； $N$ 为均值为0、标准差为1的正态分布随机向量。

考虑阻尼影响，以黏滞阻尼理论为基础，根据下式^［

18］：

[\begin{matrix} a_{0} \\ a_{1} \end{matrix}] = 2 \frac{ω_{m} ω_{n}}{ω_{n}^{2} - ω_{m}^{2}} [\begin{matrix} ω_{n} & - ω_{m} \\ - \frac{1}{ω_{n}} & \frac{1}{ω_{m}} \end{matrix}] [\begin{matrix} ξ_{m} \\ ξ_{n} \end{matrix}]

（28）

由结构固有频率计算出瑞利阻尼系数 $a_{0}$ 和 $a_{1}$ ，其中 $ω_{m}$ 和 $ω_{n}$ 一般为结构的基频和对结构动力响应贡献显著的高阶频率， $ξ_{m}$ 和 $ξ_{n}$ 为对应的阻尼比。

本文取结构的第1阶和第5阶固有频率进行计算，对应阻尼比均取0.02。

2.2 简谐荷载作用下的损伤识别

模拟简谐荷载作用在军山大桥的176号结点处（对应桥墩处的16号单元）。采样频率为200 Hz。施加荷载形式如图3所示，为两个简谐荷载叠加，如下式所式：

F = 20 s i n (2 t) + 25 s i n (5 t)

（29）

图3 简谐荷载

Fig.3 The harmonic load

将Newmark‑ $β$ 法计算损伤结构的动力响应作为现场试验数据。根据式（5）将损伤结构的动力响应数据代入到未损伤结构的动力平衡方程中，求出结构各结点的时域残余力向量。利用求出的时域残余力向量矩阵计算损伤位置指标，从而确定损伤单元。最后根据式（23）~（26）计算损伤单元对应的损伤程度。

2.2.1 损伤时刻识别结果

由前述对时域残余力向量的分析可知，当结构未发生损伤时，时域残余力向量理论值为0；结构发生损伤后，时域残余力向量非零，考虑到计算误差和测量噪声的影响，可通过分析相关单元结点处时域残余力向量值的变化来判断损伤发生的时刻。以16号单元对应的176号结点为例，损伤单元对应结点的时域残余力向量值与时间的关系如图4所示。

图4 简谐荷载作用下工况1~3的损伤时刻识别结果

Fig.4 Identification of damage existence for cases 1~3 under harmonic loading

图4（a），（b）表明，在不含噪声的工况1和工况2中，176号结点的时域残余力向量值在第10 s之前都约为0，在第10 s突然增大，并开始变化。由此可以判断两工况下16号单元损伤的发生时刻皆为第10 s，与预设的损伤时刻相符。图4（c）给出加速度响应含10%测量噪声下的损伤识别结果（工况3），图中176号结点的残余力向量值自第10 s开始突变，且此后的值远大于第10 s之前的残余力向量值，由此可以判断该工况下，损伤时刻为第10 s。综合比较图4（a）~（c）可以得出，简谐荷载作用时，单处损伤、多处损伤以及噪声工况下，模型的损伤时刻都可通过时域残余力向量的突变来识别。

2.2.2 损伤位置与程度识别结果

由时域残余力向量的定义可知，未损伤处结点的时域残余力值会远小于损伤处的时域残余力值。在噪声影响下，未损伤处结点的时域残余力在整个时间段中的变化是相对平稳的，而损伤处结点的时域残余力在损伤时刻则会发生突变。因此用各结点的损伤位置识别指标先确定损伤位置，再计算损伤单元的损伤程度，可以减小矩阵的计算量。各结点的残余力向量在其每个自由度上都存在分量，在识别损伤位置的过程中，只需要选取结点的一个自由度方向上的时域残余力向量计算Loc值来代表此结点即可，这里选取的是y自由度方向上的时域残余力向量。

简谐荷载作用下的3个工况损伤位置识别结果如图5所示。

图5 简谐荷载作用下工况1~3的损伤位置识别结果

Fig.5 Identification of damage location in cases 1~3 under harmonic loading

图5（a）表明，175号与176号结点的损伤位置指标要远大于其他结点，可知与这两结点相关的16号单元发生损伤。图5（b），（c）分别是无噪声和有噪声下结构发生多处损伤的工况，两图中175号、176号和184号、185号结点的损伤位置指标都明显大于其他结点，由此可以确定与之相关的16号单元和64号单元发生损伤。对比图5（b）和（c）发现，在简谐荷载作用下，噪声对损伤位置的识别结果影响不大，该方法在识别结构损伤位置方面具有较强的噪声鲁棒性。

获取了损伤位置之后，根据式（23）~（26）计算损伤单元的损伤程度，得到损伤程度识别结果如图6所示，表2给出了各工况损伤程度识别结果的相对误差。

图6 简谐荷载作用下工况1~3的损伤程度识别结果

Fig.6 Damage identification results for cases 1~3 under harmonic loading

表2 简谐荷载作用下工况1~3损伤程度识别结果

Tab.2 Damage identification results for cases 1~3 under harmonic loading

工况	识别损伤	预设损伤	识别误差/%
1	0.29992	0.3	0.03
2	0.30033	0.3	0.11
2	0.21391	0.2	6.96
3	0.30037	0.3	0.12
3	0.21409	0.2	7.04

由图6和表2所示，工况1简谐荷载作用下，识别出16号单元损伤，其损伤程度为0.29992，相比于实际设置的损伤程度0.3，误差为0.03%，此单处损伤工况下，损伤程度识别精度较高；工况2是结构处于简谐荷载作用下，桥墩处的16号单元和桥塔处的64号单元发生损伤，两处的损伤程度识别结果分别为0.30033和0.21391，相对于实际的损伤程度0.3和0.2，误差分别为0.11%和6.96%；在噪声条件下，两单元的损伤程度识别结果分别为0.30037和0.21409，误差分别为0.12%和7.04%，相较于不含噪声的工况2误差略大。因为结构的时域残余力向量矩阵实质上是病态矩阵，对扰动比较敏感，而且当简谐荷载作用在16号单元时，64号单元的动力响应较小，容易受到噪声影响，但对比工况2和3，可以发现噪声对损伤单元的损伤程度识别结果影响较小。

2.3 地震荷载作用下的损伤识别

为了研究时域残余力向量方法在地震荷载下的损伤识别效果，设置了工况4~6。地震荷载采样频率为200 Hz，总的时间步为7851，地震波加速度时程曲线如图7所示。

图7 地震波加速度时程曲线

Fig.7 Seismic acceleration time history curve

2.3.1 损伤时刻识别结果

计算过程与简谐荷载作用工况类似，同样以176号结点为例，得到地震荷载作用下工况4~6的损伤时刻识别结果如图8所示。

图8 地震荷载作用下工况4~6损伤时刻识别结果

Fig.8 Identification of damage existence for cases 4~6 under seismic loading

图8表明，在地震作用下的3个工况中，176号结点的时域残余力向量皆在第10 s时突变，与预设的损伤发生时刻相符。比较图8（a），（b）的突变时刻可知，在结构发生单处损伤和两处损伤时，通过时域残余力向量都能较为准确地识别出损伤时刻，结合图8（c）可知，噪声的影响较小。

2.3.2 损伤位置和程度识别结果

类似地，求出结构各结点损伤位置指标如图9所示。图9（a）表明，第175和176号结点的损伤位置指标要远大于其他结点，说明与这两结点相关的16号单元发生了损伤，与预设的单损伤工况的损伤位置相符。图9（b）则证明了在地震荷载作用下，该方法在结构发生多处损伤时的损伤位置识别的有效性。对比分析图9（b）和（c），可以发现噪声对损伤位置的识别影响较小，在加速度含10%水平的测量噪声的情况下，该方法仍能识别出损伤的位置。

图9 地震荷载作用下工况4~6损伤位置识别结果

Fig.9 Identification of damage location in cases 4~6 under seismic loading

根据上述分析确定的损伤位置，计算损伤单元的损伤程度参数，结果如图10所示，表3给出各工况识别结果的相对误差。

图10 地震荷载作用下工况4~6损伤程度识别结果

Fig.10 Damage identification results for cases 4~6 under seismic loading

表3 地震荷载作用下工况4~6损伤程度识别结果

Tab.3 Damage identification results for cases 4~6 under seismic loading

工况	识别损伤	预设损伤	识别误差/%
4	0.29995	0.3000	0.02
5	0.29995	0.3000	0.02
5	0.20000	0.2000	0.00
6	0.29869	0.3000	0.44
6	0.20478	0.2000	2.39

如图10和表3所示，工况4地震荷载作用下，识别出16号单元损伤，其损伤程度为0.29995，相比于实际设置的损伤程度0.3，误差为0.02%，此单处损伤工况下，损伤程度识别误差精度较高；工况5是结构处于地震荷载作用下，桥墩处的16号单元和桥塔处的64号单元发生损伤，两处的损伤程度识别结果分别为0.29995和0.2000，相对于实际的损伤程度0.3和0.2，误差分别为0.11%和0.00%，识别结果精度较高；在结构加速度含有噪声的工况下，16号单元和64号单元的损伤程度识别结果分别为0.29869和0.20478，误差分别为0.44%和2.39%，相较于无噪声工况误差较大。由表2可知，简谐荷载作用下的多损伤工况2和3的损伤单元64的损伤程度识别误差分别为6.96%和7.04%，要大于地震荷载作用下的对应工况5和6该单元的识别误差，这是因为简谐荷载为单点激励，这种工况下64号单元的响应较小，损伤识别结果受干扰程度较大，而地震荷载为多点激励，结构整体都受到了作用，误差和噪声的影响较小。结果很好地说明了本文方法在多点激励下的识别效果要好于单点激励。

2.4 小损伤工况下的损伤识别

为了研究时域残余力向量方法在结构所受损伤较小情况下的损伤识别效果，设置了工况7~9。外荷载为地震荷载。

2.4.1 损伤时刻识别结果

以176号结点为例，由时域残余力向量方法得到地震荷载作用下的小损伤工况7~9的损伤时刻识别结果，如图11所示。

图11 地震荷载作用下工况7~9损伤时刻识别结果

Fig.11 Identification of damage existence for cases 7~9 under seismic loading

图11表明，在结构发生小损伤的工况下，176号结点的时域残余力向量皆在第10 s时发生突变，与预设的损伤时刻相符。比较图11（a）~（c），在结构发生小损伤时，对于不含噪声的工况7和包含噪声影响的工况8，9，通过时域残余力向量都能较为准确地识别出结构的损伤时刻，噪声对识别损伤时刻的影响较小。

2.4.2 损伤位置和程度识别结果

分别求出小损伤工况7~9结构各结点损伤位置指标如图12（a）~（c）所示。

图12 地震荷载作用下工况7~9损伤位置识别结果

Fig.12 Identification of damage location in cases 7~9 under seismic loading

图12（a）表明，175和176号结点的损伤位置指标远大于其他结点，说明工况7下这两个结点对应的16号单元发生了损伤，与预设的损伤位置相符。工况8和9分别在工况7的基础上添加了10%和2%水平的测量噪声，对比分析图12（b），（c）和图9（a），可以发现在噪声影响下仍能较为准确地识别出结构发生小损伤时的损伤位置，说明噪声的影响较小。

根据损伤位置的识别结果，可以确定小损伤工况7~9皆为与175和176号结点对应的第16号单元发生损伤，计算其损伤程度和识别误差如图13和表4所示。

图13 地震荷载作用下工况7~9损伤程度识别结果

Fig.13 Damage identification results for cases 7~9 under seismic loading

表4 地震荷载作用下工况7~9损伤程度识别结果

Tab.4 Damage identification results for cases 7~9 under seismic loading

工况	识别损伤	预设损伤	识别误差/%
7	0.029995	0.03	0.02
8	0.033338	0.03	11.13
9	0.030858	0.03	2.86

由图13和表4所示，工况7识别出结构16号单元损伤，其损伤程度为0.029995，相比于预设的损伤程度0.03，误差为0.02%。在不含噪声的情况下，可以较为准确地识别出结构发生小损伤时的损伤程度。工况8在工况7的基础上在结构加速度响应中添加了10%水平的测量噪声，此工况下识别出16号单元的损伤程度为0.033338，相较于实际的损伤程度0.03，识别误差为11.13%，此时误差超过10%；工况9在工况7的基础上在结构加速度响应中添加了2%水平的测量噪声，此工况下识别出16号单元的损伤程度为0.030858，识别误差为2.86%，此时误差较小。对比分析工况7~9的损伤程度识别结果可知，对于结构发生小损伤的工况，其损伤程度的识别结果误差受噪声影响较大，在噪声水平较高时，其识别误差较大，在没有噪声或者噪声水平较低时，仍能较为准确地识别出损伤单元的损伤程度。

3 结论

本文在传统模态残余力向量的基础上，利用时域动力响应数据，推导了强迫振动下时域残余力向量的表达式。根据损伤位置和未损伤位置在时域残余力向量上表现出来的差异，利用结点时域残余力向量值确定损伤时刻，利用其变异系数定义损伤位置指标用来确定损伤位置。在此基础上，进一步推导了刚度矩阵的特征分解公式，由分解公式可知，损伤前后单元刚度矩阵特征向量不发生变化而其非零特征值随损伤程度同步变化，由此可以提取反映损伤程度的单元刚度特征值的变化系数，作为结构的损伤程度参数，用以确定单元损伤的程度。与模态残余力向量法不同，时域残余力向量包含与损伤相关的时间信息，可以通过其判断损伤发生的时刻。将此方法用于武汉军山大桥的损伤识别，识别结果表明：所推导的时域残余力向量可用于单点激励和多点激励下的结构损伤识别，能够准确地识别损伤的发生时刻，定位损伤单元，且能够准确识别单元的损伤程度，具备较强的噪声鲁棒性，多点激励下的损伤识别结果比单点激励下的结果更为准确。目前为保障损伤识别的精度，残余力向量方法需要大量测点的动力响应来计算残余力向量，在之后的研究中，结合动态缩聚^［

19］、子结构方法^{［参考文献 20‑21}20‑21］等有望有效解决测点需求多的难题。

参考文献

Zimmerman D C， Kaouk M. Structural damage detection using a subspace rotation algorithm［C］//33rd Structures， Structural Dynamics and Materials Conference. Dallas， TX， USA， 1992： 2341-2350. [百度学术]

周先雁，沈蒲生，程翔云. 用振动参数识别技术对混凝土框架进行破损评估［J］. 土木工程学报， 1998， 31（2）： 39-45. [百度学术]

Zhou Xianyan， Shen Pusheng， Cheng Xiangyun. Damage assessment of reinforced concreter structures by vibration parameter identification technique［J］. China Civil Engineering Journal， 1998， 31（2）： 39-45. [百度学术]

Fukunaga H， Kameyama M， Ogi Y. Damage identification of laminated composite structures based on dynamic residual forces［J］. Advanced Composite Materials， 2001， 10（2-3）： 209-218. [百度学术]

Ge M， Lui E M. Structural damage identification using system dynamic properties［J］. Computers and Structures， 2005， 83（27）： 2185-2196. [百度学术]

Amiri Gholamreza Ghodrati， Hosseinzadeh Ali Zare， Abdollah Bagheri， et al. Damage prognosis by means of modal residual force and static deflections obtained by modal flexibility based on the diagonalization method［J］. Smart Materials and Structures， 2013， 22（7）： 075032. [百度学术]

何伟，陈淮，王博，等. 运用改进残余力向量法的结构损伤识别研究［J］. 振动、测试与诊断， 2009， 29（4）： 379-382. [百度学术]

He Wei， Chen Huai， Wang Bo， et al. Structural damage identification by using improved residual force vector［J］. Journal of Vibration， Measurement and Diagnosis， 2009， 29（4）： 379-382. [百度学术]

Behtani A， Bouazzouni A， Khatir S， et al. Damage localization and quantification of composite stratified beam structures using residual force method［J］. Journal of Physics： Conference Series， 2017， 842： 012028. [百度学术]

Wu J， Luo Y F， Wang L. Practical residual force decomposition method for damage identification of existing reticulated shells［J］. International Journal of Steel Structures， 2019， 19（6）： 2024-2038. [百度学术]

张干，孙熙平，王倩，等. 基于残余力向量的高桩码头基桩损伤诊断研究［J］. 振动与冲击， 2016， 35（4）： 188-192. [百度学术]

Zhang Gan， Sun Xiping， Wang Qian， et al. Damage diagnosis for high-piled wharf piles based on residual force vector［J］. Journal of Vibration and Shock， 2016， 35（4）： 188-192. [百度学术]

于德介，雷慧. 一种基于神经网络的结构参数识别方法［J］. 湖南大学学报（自然科学版）， 1999， 26（4）： 39-43. [百度学术]

Yu Dejie， Lei Hui. A method for structural parameters identification using neural networks［J］. Journal of Hunan University（Natural Sciences Edition）， 1999， 26（4）： 39-43. [百度学术]

袁颖，林皋，闫东明，等. 基于残余力向量法和改进遗传算法的结构损伤识别研究［J］. 计算力学学报， 2007， 24（2）： 224-230. [百度学术]

Yuan Ying， Lin Gao， Yan Dongming， et al. Study on structural damage identification based on residual force method and improved genetic algorithm［J］. Chinese Journal of Computational Mechanics， 2007， 24（2）： 224-230. [百度学术]

Yun G J， Ogorzalek K， Dyke S J， et al. A structural damage detection method based on subset selection and evolutionary computation［C］//Proceedings of the 2008 Structures Congress. Vancouver， British Columbia， Canada， 2008： 24-26. [百度学术]

Saberi M. Damage detection of space structures using charged system search algorithm and residual force method［J］. Iranian Journal of Science and Technology- Transactions of Civil Engineering， 2015， 39（C2）： 215-229. [百度学术]

赵一霖，刘济科，吕中荣，等. 基于残余力向量法和树种算法的结构损伤识别二步法［J］. 中山大学学报（自然科学版）， 2017， 56（4）： 46-50. [百度学术]

Zhao Yilin， Liu Jike， Lyu Zhongrong， et al. Structural damage identification based on residual vectors and Tree-seed algorithm［J］. Acta Scientiarum Naturalium Universitatis Sunyatseni， 2017， 56（4）： 46-50. [百度学术]

杨秋伟，周聪，李翠红，等. 基于静力残余力向量的结构损伤评估方法［J］. 计算力学学报， 2021， 38（5）： 625-630. [百度学术]

Yang Qiuwei， Zhou Cong， Li Cuihong， et al. Structural damage assessment based on static residual force vector［J］. Chinese Journal of Computational Mechanics， 2021， 38（5）： 625-630. [百度学术]

Doebling S W， Peterson L D， Alvin K F. Experimental determination of local structural stiffness by disassembly of measured flexibility matrices［J］. Journal of Vibration and Acoustics， 1998， 120（4）： 949-957 [百度学术]

和禹含，刘坤. 基于奇异值分解和时域灵敏度分析的框架结构损伤识别方法［J］. 声学与振动， 2018， 6（1）： 8-15. [百度学术]

He Yuhan， Liu Kun. Damage detection method of frame structure based on singular value decomposition and time domain sensitivity analysis［J］. Open Journal of Acoustics and Vibration， 2018， 6（1）： 8-15. [百度学术]

Clough R W， Penzien J. Dynamics of Structures［M］. New York： McGraw-Hill， 1975. [百度学术]

Weng Shun， Tian Wei， Zhu Hongping， et al. Dynamic condensation approach to calculation of structural responses and response sensitivities［J］. Mechanical Systems and Signal Processing， 2017， 88： 302-317. [百度学术]

Weng Shun， Xia Yong， Xu Youlin， et al. Substructure based approach to finite element model updating［J］. Computers and Structures， 2011， 89（9-10）： 772-782. [百度学术]

翁顺，朱宏平. 基于有限元模型修正的土木结构损伤识别方法［J］. 工程力学， 2021， 38（3）： 1-16. [百度学术]

Weng Shun， Zhu Hongping. Damage identification of civil structures based on finite element model updating［J］. Engineering Mechanics， 2021， 38（3）： 1-16. [百度学术]