一种基于非迭代有限元模型修正的桥梁损伤识别方法

1.本发明涉及于桥梁安全监测和检测领域，具体地说是一种基于非迭代有限元模型修正的桥梁损伤识别方法，更新结果可用于桥梁结构的损伤识别和安全评估等。


背景技术：

2.随着有限元法的发展，使用有限元模型可以计算桥梁各种荷载、各种边界条件下的动力响应，分析速度快，结构设计周期短，与结构动力试验比较效率高且费用低，而且还可以广泛应用于桥梁结构健康监测。因此，一个可靠比较准确的有限元模型是进行桥梁结构损伤识别、健康诊断以及桥梁工作状态评估与预估的基础，即损伤识别的前提是修正初始有限元模型。在初始有限元模型基础上更新桥梁刚度和阻尼矩阵方法可以用来进行结构损伤识别和健康状态评估，因为结构的损伤必然会带来结构参数的变化，结构参数的变化必然导致试验数据的变化，根据桥梁更新后的刚度和阻尼矩阵矩阵元素的数据变化识别结构损伤。在此结构健康监测技术中，应该使建立的桥梁有限元模型中的刚度和阻尼矩阵能够全面、正确地反映结构的真实状况。
3.目前为止，大量研究工作的开展形成了基于不同种类的桥梁有限元模型修正的识别损伤方法。但大多数有限元模型修正方法使用的优化技术是迭代的，计算桥梁更新后的刚度和阻尼矩阵需要耗费大量的时间和成本，且基于频率和振型的迭代模型修正需要事先确定频率和振型权重。


技术实现要素：

4.本发明是为避免上述现有基于迭代有限元模型修正的桥梁损伤识别方法所存在的不足之处，提供一种基于非迭代有限元模型修正的桥梁损伤识别方法，只需更新一次状态矩阵，就能得到与实际桥梁模态参数匹配的刚度和阻尼矩阵，从而提高识别损伤的效率，并减少计算成本，进而能有效解决迭代的有限元模型修正计算量大和收敛速率较低等问题。
5.本发明为达到上述发明目的，采用如下技术方案：
6.本发明一种非迭代有限元模型修正的桥梁损伤识别方法的特点包括以下步骤：
7.步骤1：获取实际桥梁的基本参数，用于确定初始桥梁有限元模型；利用式(1)计算初始桥梁有限元模型的状态空间矩阵ac，从而得到所述状态空间矩阵ac的第i个特征值和第i个特征向量所构成的第i个特征对
[0008][0009]
式(1)中，c表示连续时间域上初始桥梁的状态空间模型，n1表示初始桥梁有限元模型的自由度，z(n1)，i(n1)分别是n1×
n1阶零矩阵和单位矩阵，m，c，k分别是初始桥梁有限元模型的n1×
n1阶时不变质量矩阵、阻尼和刚度矩阵；
[0010]
步骤2：识别连续时间域上实际桥梁有限元模型的状态空间矩阵：
[0011]
步骤2.1：将初始桥梁有限元模型相应的自由度节点作为加速度测点，利用加速度传感器获取实际桥梁加速度响应s(t)，通过随机子空间算法对加速度响应s(t)识别状态空间矩阵并对其进行标准化之后，再进行平均化，从而得到实际桥梁有限元模型离散时间域上不同阶数的状态空间矩阵a
′r；
[0012]
步骤2.2：通过零阶保持器方法将实际桥梁有限元模型离散时间域上不同阶数的状态空间矩阵a
′r变换到连续时间域上不同阶数的状态空间矩阵ar；
[0013]
步骤3：确定状态空间矩阵的可赋值的目标特征值和特征向量对
[0014]
步骤3.1：对不同阶数的ar计算其第2i个特征值及其相应的特征向量从而构成第2i个特征对
[0015]
步骤3.2：利用式(2)计算不同阶数的ar的第2i阶特征值从而根据不同阶数的ar的第2i阶特征值的虚部，计算出频率和阻尼比所构成的参数对再通过稳定图和相关的实际约束去除所述参数对的非物理根；
[0016][0017]
式(2)和式(3)中，分别是实际桥梁在物理空间上第i阶模态的固有频率和阻尼比，k为复数单位；
[0018]
步骤3.3：利用式(3)和式(4)计算第i阶目标特征值
[0019][0020][0021]
式(3)和式(4)中，分别是实际桥梁第i阶模态稳定的频率和阻尼比；
[0022]
步骤3.4：利用式(5)计算2n1阶数的状态空间矩阵ar的第2i阶特征值所对应的特征向量从而根据特征向量计算实际桥梁第i阶振型
[0023][0024]
步骤3.5：当初始桥梁有限元模型自由度小于等于阈值α时，利用式(6)和式(7)计算与目标特征值对应的目标特征向量当初始桥梁有限元模型自由度大于阈值α时，令目标特征向量为
[0025][0026][0027]
步骤4：计算增益矩阵gc：
[0028]
步骤4.1：当i＝1，2，...，2m1时，初始化第i个特征值初始化第i个特征向量当i＝2m1 1，...，2n1时，初始化第i个特征值初始化第i个特征向量
[0029]
步骤4.2：定义第i个组合矩阵并计算其零空间矩阵其中，b表示满足rank([λii(ni)-a
c b])≥n1且矩阵的上半部分元素为“0”的输入矩阵，rank表示矩阵的秩；和分别表示第i个零空间矩阵的上半部分和下半部分并具有与ac相同的阶数；
[0030]
步骤4.3：利用式(8)计算第i个vi的关联向量zi：
[0031][0032]
式(7)中，表示moore-penrose伪逆；
[0033]
步骤4.4：利用式(9)计算第i个zi的关联向量wi：
[0034][0035]
步骤4.5：利用式(10)计算增益矩阵gc：
[0036][0037]
式(10)中，r和i分别表示相应复数元素的实部和虚部，gc表示具有赋值状态空间矩阵的前2m1阶目标特征值和特征向量功能的增益矩阵；
[0038]
步骤5：利用增益矩阵gc计算更新后的状态矩阵
[0039]
步骤5.1：利用式(11)选取输入矩阵b：
[0040][0041]
式(11)中，b1和b2分别表示矩阵b的左下分块和右下分块；
[0042]
步骤5.2：利用式(12)计算更新后的桥梁的状态矩阵
[0043][0044]
式(12)中，g1，g2，g3，g4表示顺时针划分增益矩阵gc后的四个分块；
[0045]
步骤6：利用(13)式计算更新后的桥梁的刚度矩阵和阻尼矩阵
[0046][0047]
步骤7：识别桥梁损伤单元j及其损伤程度γj，j＝1，2，...，n1：
[0048]
步骤7.1：计算刚度矩阵在对角线上的第i个元素相对于刚度矩阵k在对角线上的第i个元素的数值变化率δi，i＝1，2，...，n1；当δi大于阈值η时，则表示节点i为损伤节点，然后利用损伤节点的连通性信息确定损伤单元j，j＝1，2，...，n1；
[0049]
步骤7.2：确定与损伤单元j相关联的所有损伤节点的数值变化率并取平均值，即
为损伤单元j的损伤程度γj。
[0050]
与现有技术相比，本发明的有益效果体现在：
[0051]
1、本发明大大降低了桥梁更新刚度和阻尼矩阵过程中的计算量和计算成本。它采用随机子空间识别算法(ssi)建立状态空间模型，并利用其连续时间状态矩阵的特征值和特征向量对有限元生成的系统矩阵建立的初始状态空间模型进行更新，避免了将状态空间映射回物理空间。该优化过程只需更新一次的非迭代特性在很大程度上降低了计算量和计算成本，提高了计算效率，更新结果可以更快地实时用于桥梁损伤识别和健康状态评估。
[0052]
2、本发明避免了以往基于频率和振型的迭代的有限元模型修正识别损伤需要事先确定频率和振型权重的困境，并且不需要模态序列信息，对桥梁低强度损伤也能很好地识别。
附图说明
[0053]
图1为本发明的说明书摘要流程图；
[0054]
图2为本发明的数值模拟等截面简支混凝土梁图；
[0055]
图3为本发明的脉冲激励下等截面简支混凝土梁的加速度响应图；
[0056]
图4为本发明的脉冲激励下等截面简支混凝土梁的加速度响应的频谱图；
[0057]
图5为本发明的等截面简支混凝土梁的随机子空间识别(ssi)的稳定图；
[0058]
图6为本发明的等截面简支混凝土梁的识别损伤与实际损伤的对比图；
[0059]
图7为本发明的数值模拟等截面简支梁桥图；
[0060]
图8为本发明的脉冲激励下等截面简支梁桥的加速度响应图；
[0061]
图9为本发明的脉冲激励下等截面简支梁桥的加速度响应的频谱图；
[0062]
图10为本发明的等截面简支梁桥的随机子空间识别(ssi)的稳定图；
[0063]
图11为本发明的等截面简支梁桥的识别损伤与实际损伤的对比图；
具体实施方式
[0064]
算例1：图2所示的等矩形截面简支混凝土梁，梁跨长为6m，弹性模量为3.2gpa，密度为2500kg/m3，截面宽0.2m，高0.25m，惯性矩为2.604
×
10-4
m4。将混凝土梁等间距划分为6个平面欧拉梁单元，用matlab软件建立未损伤下的初始梁有限元模型。假设单元4处发生损伤，相应的刚度降低40％，阻尼比为1％rayleigh阻尼矩阵，用matlab软件建立损伤下的实际梁有限元模型。在节点5处沿-y方向作用500n的脉冲激励，采用newmark-β法计算梁的动力响应。采样频率为500hz且在加速度数据中加入5％的均方根噪声。利用加速度响应数据和未损伤下的初始梁有限元模型更新损伤下的混凝土梁的实际有限元模型，如图1所示，一种非迭代的桥梁刚度和阻尼矩阵更新方法步骤如下：
[0065]
步骤1：获取实际桥梁的基本参数，梁跨长为6m，弹性模量为3.2gpa，密度为2500kg/m3，截面宽0.2m，高0.25m，惯性矩为2.604
×
10-4
m4，将混凝土梁等间距划分为6个平面欧拉梁单元，用于确定未损伤下的初始桥梁有限元模型；利用式(1)计算初始桥梁有限元模型的状态空间矩阵ac，从而得到所述状态空间矩阵ac的第i个特征值和第i个特征向量所构成的第i个特征对
[0066][0067]
式(1)中，c表示连续时间域上初始桥梁的状态空间模型，n1表示初始桥梁有限元模型的自由度，z(n1)，i(n1)分别是n1×
n1阶零矩阵和单位矩阵，m，c，k分别是初始桥梁有限元模型的n1×
n1阶时不变质量矩阵、阻尼和刚度矩阵；
[0068]
步骤2：识别连续时间域上实际桥梁有限元模型的状态空间矩阵：
[0069]
步骤2.1：在节点5处沿-y方向作用500n的脉冲激励，将初始桥梁有限元模型相应的自由度节点作为加速度测点，利用加速度传感器获取实际桥梁加速度响应s(t)，采样频率为500hz且在加速度数据中加入5％的均方根噪声；加速度响应图如图3所示，加速度响应的频谱图如图4所示；通过随机子空间算法(ssi)对加速度响应s(t)识别状态空间矩阵并对其进行标准化之后，再进行平均化，减小噪声对信号的影响，从而得到实际桥梁有限元模型离散时间域上5至30不同阶数的状态空间矩阵a
′r；
[0070]
步骤2.2：通过零阶保持器方法(zoh)将实际桥梁有限元模型离散时间域上5至30不同阶数的状态空间矩阵a
′r变换到连续时间域上5至30不同阶数的状态空间矩阵ar；
[0071]
步骤3：确定状态空间矩阵的可赋值的目标特征值和特征向量对m1＝6：
[0072]
步骤3.1：对5至30不同阶数的ar计算其第2i个特征值及其相应的特征向量从而构成第2i个特征对
[0073]
步骤3.2：利用式(2)计算5至30不同阶数的ar的第2i阶特征值从而通过5至30不同阶数的ar的第2i阶特征值的虚部计算出频率和阻尼比所构成的参数对再通过稳定图和相关的实际约束去除所述参数对的非物理根；稳定图如图5所示；
[0074][0075]
式(2)和式(3)中，分别是实际桥梁在物理空间上第i阶模态的固有频率和阻尼比，k为复数单位。
[0076]
步骤3.3：利用式(3)和式(4)计算目标特征值
[0077][0078][0079]
式(3)和式(4)中，分别是实际桥梁第i阶模态稳定的频率和阻尼比；
[0080]
步骤3.4：利用式(5)计算2n1阶数的ar的第2i阶特征值所对应的特征向量从而根据特征向量计算实际桥梁第i阶振型
[0081][0082]
步骤3.5：当初始桥梁有限元模型自由度小于等于阈值α时，α＝20，利用式(6)和式
(7)计算与目标特征值对应的目标特征向量当初始桥梁有限元模型自由度大于阈值α时，令目标特征向量为
[0083][0084][0085]
步骤4：计算增益矩阵gc：
[0086]
步骤4.1：当i＝1，2，...，2m1时，初始化第i个特征值初始化第i个特征向量当i＝2m1 1，...，2n1时，初始化第i个特征值初始化第i个特征向量
[0087]
步骤4.2：定义第i个组合矩阵并计算其零空间矩阵其中，b表示满足rank([λii(ni)-a
c b])≥n1且矩阵的上半部分元素为“0”的输入矩阵，rank表示矩阵的秩；和分别表示第i个零空间矩阵的上半部分和下半部分并具有与ac相同的阶数；
[0088]
步骤4.3：利用式(8)计算第i个vi的关联向量zi：
[0089][0090]
式(7)中，表示moore-penrose伪逆。
[0091]
步骤4.4：利用式(9)计算第i个zi的关联向量wi：
[0092][0093]
步骤4.5：利用式(10)计算增益矩阵gc：
[0094][0095]
式(10)中，r和i分别表示相应复数元素的实部和虚部；gc表示具有赋值状态空间矩阵的前2m1阶目标特征值和特征向量功能的增益矩阵；
[0096]
步骤5：利用增益矩阵gc计算更新后的状态矩阵
[0097]
步骤5.1：利用式(11)选取输入矩阵b：
[0098][0099]
式(11)中，b1和b2分别表示矩阵b的左下分块和右下分块；b1＝i(n1)，b2＝i(n1)；
[0100]
步骤5.2：利用式(12)计算更新后的桥梁的状态矩阵
[0101][0102]
式(12)中，g1，g2，g3，g4表示顺时针划分增益矩阵gc后的四个分块；
[0103]
步骤6：利用(13)式计算更新后的桥梁的质量矩阵和阻尼矩阵
[0104][0105]
步骤7：识别桥梁损伤单元j及其损伤程度γj，j＝1，2，...，n1：
[0106]
步骤7.1：计算刚度矩阵在对角线上的第i个元素相对于刚度矩阵k在对角线上的第i个元素的数值变化率δi，i＝1，2，...，n1；当δi大于阈值η时，则表示节点i为损伤节点，然后利用损伤节点的连通性信息确定损伤单元j，j＝1，2，...，n1；
[0107]
步骤7.2：确定与损伤单元j相关联的所有损伤节点的数值变化率并取平均值，即为损伤单元j的损伤程度γj。
[0108]
利用修正模型更新后的桥梁的质量矩阵和阻尼矩阵计算前三阶频率、阻尼比修正结果与实际模型和初始模型的对比，修正模型和初始模型与实际模型的前三阶振型之间的mac如表1所示。
[0109]
表1
[0110][0111]
表1数据说明修正后的桥梁有限元模型能够较为准确地代表实际桥梁有限元模型，由此，可用更新后的桥梁质量矩阵识别桥梁损伤位置及其程度，如图6所示。
[0112]
算例2：图7所示的两跨等截面简支梁桥，桥梁跨长为30m，弹性模量27.5gpa，密度为2400kg/m3，截面宽2.4m，高1m，惯性矩为0.2m4。将桥梁等间距划分为30个平面欧拉梁单元，用matlab软件建立未损伤下的初始梁有限元模型。假设单元4和单元26处发生损伤，相应的刚度降低40％，阻尼比为1％rayleigh阻尼矩阵，节点7处沿-y方向作用5000kn的脉冲激励，采用newmar
·
k-β法计算梁动力响应。采样频率为500hz且在加速度数据中加入10％的均方根噪声。计算步骤与算例1一致，加速度响应如图8所示，加速度响应的频谱如图9所示，稳定图如图10所示。利用修正模型更新后的桥梁的质量矩阵和阻尼矩阵计算的前五阶频率、阻尼比修正结果与实际模型和初始模型的对比，修正模型和初始模型与实际模型的前五阶振型之间的mac如表2所示。
[0113]
表2
[0114][0115]
表2数据说明修正后的桥梁有限元模型能够较为准确地代表实际桥梁有限元模型，由此，可用更新后的桥梁质量矩阵识别桥梁损伤位置及其程度，如图11所示。
[0116]
算例1和算例2充分说明了本发明的方法能够准确的更新实际桥梁的质量矩阵和阻尼矩阵用该方法能够大大降低更新过程中的计算量和计算成本，解决了迭代的更新过程计算量大和收敛速率较低等问题，更新结果的矩阵元素的数值变化可以更快地用于结构损伤识别和健康状态评估。