桥梁损伤位置的无损识别方法、存储介质及设备

：
1.本发明属于桥梁工程健康监测技术领域，具体涉及桥梁损伤位置的无损识别方法、存储介质及设备。


背景技术：

2.桥梁作为交通运输生命线工程中的节点和咽喉，在促进区域链接、保障物资运输、推动经济发展等方面有着不可替代的作用。随着桥梁服役年限的增长，重载交通作用、环境侵蚀、材料老化、疲劳效应等对桥梁结构的性能影响越来越大，致使桥梁局部抗力下降、损伤累积，诱发破坏甚至坍塌。如何尽早获取桥梁损伤位置，对于桥梁及时养护维修、预防桥梁大范围破坏、保障人员安全及节能环保方面意义重大。
3.目前，按照是否破坏原结构，桥梁损伤位置检测方法主要分为有损检测法和无损检测法两类，其中有损检测法主要是借助工具获取桥梁某位置的实体样本，如通过打孔取样等，并进一步观察或试验分析，这无疑会对原结构产生影响，且检测速度慢，效率低；无损检测法包括超声波探测技术、声发射探测技术、探地雷达检测技术等，但这都需要专门的设备，检测成本比较高，且需要专业的技术人员操作。因此，亟需发展适用于桥梁结构且易于操作的无损检测方法。


技术实现要素：

4.本发明为了解决现有的检测技术不能简便识别桥梁结构损伤位置的问题，进而提出一种桥梁损伤位置的无损识别方法。
5.一种桥梁损伤位置的无损识别方法，包括以下步骤：
6.针对桥梁结构，确定桥梁结构对应状态的初始值，并组成初始状态量χ0，并依据卡尔曼滤波原理确定初始状态量的协方差矩阵，简称初始状态量协方差p0；其中χ0和p0分别称作第0时间步的状态量和状态量协方差；
7.基于无迹卡尔曼滤波器算法进行初步识别，基于无迹卡尔曼滤波器算法进行初步识别的过程中，需要基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的观测误差εk和第k时间步的量测预测协方差p
yy,k
，并基于εk和p
yy,k
计算输出每步对应的灵敏参数
8.然后绘制输出灵敏参数ηk的时程曲线，如果ηk时程曲线有峰值脉冲出现，则基于自适应无迹卡尔曼滤波器算法识别损伤位置，在基于自适应无迹卡尔曼滤波器算法进行识别的过程中，需要判断计算的ηk与灵敏参数阈值η0的大小，如果ηk＜η0，则基于无迹卡尔曼滤波器算法继续识别；如果ηk≥η0，则继续执行以下步骤：
9.设置一个初始值为0且维度等于初始状态量χ0的向量，用字母l表示；取状态量协方差的全部主对角元素组成一个新的对角方阵且保持原对角元素的位置不变；依次扩大中每一个弹性模量参数对应的协方差值，且每次只对一个协方差值进行扩大，并
执行一步完整的无迹卡尔曼滤波运算，针对协方差值对应的行号或列号为z的扩大，执行一步完整的无迹卡尔曼滤波运算得到灵敏参数ηk，并令l(z)＝ηk；找到最小的l对应的位置号z
min
，该位置就是损伤位置；然后令扩大20倍，并继续执行第(k 1)时间步，并再次判断η
k 1
与灵敏参数阈值η0的大小，直到全部循环结束得出所有的损伤位置。
10.进一步地，利用无迹卡尔曼滤波器算法进行识别的过程中和利用自适应无迹卡尔曼滤波器算法进行识别的过程中，选择各梁单元节点的竖向位移响应作为观测值y。
11.进一步地，所述的桥梁结构对应状态包括桥梁各梁单元的弹性模量、桥梁结构的位移和速度。
12.进一步地，在基于自适应无迹卡尔曼滤波器算法进行识别的过程中，计算灵敏参数ηk的过程包括以下步骤：
13.步骤7.1、基于无迹卡尔曼滤波器算法的ut变换原理，使用第(k-1)时间步的状态量χ
k-1
和状态量协方差p
k-1
生成(2n 1)个sigma点，并通过状态方程求解出每个sigma点对应的状态量其中k从1开始，且k∈[1,n]，n为总的时间步数，n为状态量的维数，i为第i个sigma点，且i∈[1,2n 1]；
[0014]
步骤7.2、基于无迹卡尔曼滤波器算法的时间更新步完成从第(k-1)时间步到第k时间步的状态量和状态量协方差的更新，分别记作和公式如下所述：
[0015][0016][0017]
式中，和分别为第k时间步第i个sigma点的权重值，为第k时间步第i个sigma点对应的状态量估计值，qk为第k时间步的噪声；
[0018]
步骤7.3、基于无迹卡尔曼滤波器算法的ut变换原理，使用步骤7.2中更新的和生成(2n 1)个sigma点，并通过观测方程求解出每个sigma点对应的观测估计值
[0019]
步骤7.4、基于无迹卡尔曼滤波器算法的量测更新步计算输出第k时间步的量测预测值且
[0020]
式中，为第k时间步第i个sigma点的权重值，为第k时间步第i个sigma点对应的观测估计值；
[0021]
步骤7.5、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的观测误差εk，且
[0022]
式中，yk为第k时间步的观测值，为第k时间步的量测预测值；
[0023]
步骤7.6、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的量测预测协方差p
yy,k
，且
[0024]
式中，为第k时间步第i个sigma点的权重值，为第k时间步第i个sigma点对应的观测估计值，为第k时间步的量测预测值，rk为第k时间步的噪声；
[0025]
步骤7.7、基于步骤7.5及步骤7.6计算的εk和p
yy,k
构造灵敏参数ηk，且并计算输出每步的ηk值。
[0026]
进一步地，在基于自适应无迹卡尔曼滤波器算法进行识别的过程中，在得到计算灵敏参数ηk后还要计算第k时间步与的互协方差p
xy,k
并进行数据更新，具体过程包括以下步骤：
[0027]
步骤7.8、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步与的互协方差p
xy,k
，
[0028]
步骤7.9、更新第k时间步的卡尔曼增益矩阵：
[0029]
步骤7.10、更新并输出第k时间步的状态量：
[0030]
步骤7.11、更新并输出第k时间步的状态量协方差：
[0031]
进一步地，找到最小的l对应的位置号z
min
的具体过程包括以下步骤：
[0032]
步骤7.14、将第k时间步状态量协方差的主对角线元素中第一个弹性模量参数对应的协方差值的行号或列号记为m，将状态量协方差的主对角线元素中最后一个弹性模量参数对应的协方差值的行号或列号记为l，则总的弹性模量参数个数为(l-m 1)；
[0033]
取当前时间步中的所有主对角元素组成一个新的对角方阵且保持原对角元素的位置不变；依次扩大中每一个弹性模量参数对应的协方差值，令且每次只对一个协方差值进行扩大，的其余元素值保持不变，其中和分别代表行数列数均为z位置处的协方差值，为一个标量；λ为扩大倍数；
[0034]
z起始值为m，且z∈[m,l-m 1]；依据z＝m时计算的和第k时间步的状态量执行一步完整的无迹卡尔曼滤波运算，即执行步骤7.1～步骤7.11，并输出ηk，并令l(z)＝ηk；
[0035]
步骤7.15、令z＝m 1，继续执行步骤7.14，直到z＝l-m 1时结束；
[0036]
步骤7.16、不计零值，找到l中最小值对应的位置号，记为z
min
，即z
min
为损伤位置。
[0037]
进一步地，所述扩大倍数λ＝1
×
10
ω
，10
ω
等于初始状态量协方差p0中最小数量级的倒数。
[0038]
进一步地，所述桥梁结构对应的运动控制微分方程或有限元模型是基于欧拉-伯努利梁单元建立的。
[0039]
一种存储介质，所述存储介质中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现所述的一种桥梁损伤位置的无损识别方法。
[0040]
一种桥梁损伤位置的无损识别设备，所述设备包括处理器和存储器，所述存储器中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现所述的一种桥梁损
伤位置的无损识别方法。
[0041]
有益效果：
[0042]
本发明仅基于桥梁局部竖向位移响应作为观测值即可识别出桥梁的损伤位置，操作简单，且不破坏原有结构。同时，利用本发明对桥梁损伤位置识别，有益于桥梁结构的及时养护维修，可以有效防止大面积破坏，符合绿色、环保、节能、低碳的发展理念。
附图说明：
[0043]
为了易于说明，本发明由下述的具体实施及附图作以详细描述。
[0044]
图1为灵敏参数时程曲线图。
[0045]
图2为灵敏参数时程曲线有峰值脉冲的效果图。
[0046]
图3为图2的局部视图。
[0047]
图4为实施例的车-桥耦合系统图，其中车以一定速度从桥的一端行驶到另一端，图中，1-四分之一车辆模型的车体质量，用字母m1表示；2-四分之一车辆模型的车体与轮胎之间的悬挂刚度，用字母k1表示；3-四分之一车辆模型的车体与轮胎之间的悬挂阻尼，用字母c1表示；4-四分之一车辆模型的轮胎质量，用字母m2表示；5-四分之一车辆模型的轮胎与桥梁之间的接触刚度，用字母k2表示；6-轮胎与桥梁接触点；7-简支边界约束的固定端；8-梁单元；9-梁单元节点；10-简支边界约束的滑动端。
具体实施方式：
[0048]
为使本发明的目的、技术方案和优点更加清楚明了，下面通过附图中示出的具体实施例来描述本发明。但是应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。
[0049]
在此，还需要说明的是，为了避免因不必要的细节而模糊了本发明，在附图中仅仅示出了与根据本发明的方案密切相关的结构和/或处理步骤，而省略了与本发明关系不大的其它细节。
[0050]
损伤位置识别时需要先制定损伤指标，并基于此判断桥梁结构是否发生损伤。一般认为桥梁结构的质量恒定，阻尼为质量和刚度的函数，且为了方便分析，一般都假定为瑞利阻尼，由此，桥梁损伤指标一般通过刚度参数反映。进一步，刚度由弹性模量和截面惯性矩组成，且一般忽略桥梁横截面的变化。因此，更深层次的损伤指标由弹性模量这一微观参数表示。此外，考虑到损伤位置识别的时效性及桥梁的尺寸特征，一般基于欧拉-伯努利梁单元构建桥梁的有限元模型，或者基于欧拉-伯努利梁单元建立桥梁结构的运动控制微分方程；并将此作为真实桥梁结构的合理简化。桥梁结构的位移、速度等响应可通过梁单元节点自由度的对应值表征，且每个梁单元的弹性模量属性都可独立赋予。由此，可根据不同梁单元的属性判断桥梁具体的损伤位置。
[0051]
具体实施方式一：
[0052]
本实施方式为一种桥梁损伤位置的无损识别方法，由于本发明是基于桥梁局部竖向位移响应进行识别的，所以本发明是一种基于桥梁局部竖向位移响应识别桥梁损伤位置的无损识别方法，本质上是一种基于自适应无迹卡尔曼滤波器算法的识别方法。
[0053]
本实施方式所述的一种桥梁损伤位置的无损识别方法，包括以下步骤：
[0054]
步骤1、根据达朗贝尔原理、结构动力学及有限元理论，推导结构的运动控制微分方程或建立结构有限元模型，所述结构为桥梁结构，为方便表述，简称结构；
[0055]
在此过程中，基于欧拉-伯努利梁单元建立桥梁结构的运动控制微分方程或有限元模型，并将各个梁单元的弹性模量作为损伤指标；
[0056]
步骤2、将结构位移x、结构速度以及弹性模量e组成的向量称作状态量，用符号表示；
[0057]
在一些实施例中，根据状态量χ，依据步骤1中结构的运动控制微分方程，并基于线性代数矩阵运算及数值分析的数值微分和积分运算，或，基于线性代数矩阵运算及结构动力学的newmark-β法等数值积分方法推导出求解状态量的方程关系；
[0058]
或者，
[0059]
在另一些实施例中，根据状态量χ，依据步骤1中结构有限元模型进行状态量的输出设置，并将其作为状态方程；该过程中考虑噪声影响；
[0060]
步骤3、将传感器采集的各梁单元节点的竖向位移响应作为观测值y。
[0061]
观测值一般是指包含位移、速度、加速度、应变、应力、力等一系列能通过传感器测量的物理量，对于桥梁结构，本实施方式中使用位移测量值作为观测值。
[0062]
本发明中，状态方程是状态量的微分运算，状态量写成位移和速度的形式，是方便状态方程的推导，因为对位移求导是速度，对速度求导是加速度，而根据结构的运动控制微分方程很容易就能推出状态方程。另外，当观测值为位移时，因为状态量中的位移和观测值中的位移含义相同，那么观测方程中的位移关系也很容易得到。概括说，状态量中的位移和速度是为状态方程推导以及观测方程中的位移关系服务的，观测值中的位移是为修正观测方程计算的位移服务的。由于此过程的具体内容为公知常识，所以本发明中不再赘述。需要说明的是，本发明中的速度相当于中间量，不需要输出，它的存在只是为状态方程的推导服务的，它的计算由算法递推迭代完成。
[0063]
根据观测值类型，依据步骤1中结构的运动控制微分方程并基于线性代数矩阵运算、数学移项、数学合并同类项等知识推导出求解观测值的方程关系；
[0064]
或者，
[0065]
根据观测值类型，依据步骤1中结构有限元模型进行相应观测值的输出设置，并将其作为观测方程；该过程中考虑噪声影响；
[0066]
步骤4、将结构的初始位移、结构的初始速度以及各弹性模量的初始值组成的向量称作初始状态量，用符号χ0表示，并根据卡尔曼滤波器算法原理得出初始状态量的协方差矩阵，简称初始状态量协方差，用符号p0表示，其中χ0和p0分别称作第0时间步(启始步)的状态量和状态量协方差；
[0067]
步骤5、基于无迹卡尔曼滤波器算法进行初步识别，过程如下：
[0068]
步骤5.1、基于无迹卡尔曼滤波器算法的ut变换原理，使用第(k-1)时间步的状态量χ
k-1
和状态量协方差p
k-1
生成(2n 1)个sigma点，并通过状态方程求解出每个sigma点对应的状态量其中k从1开始，且k∈[1,n]，n为总的时间步数，n为状态量的维数，i为第i个sigma点，且i∈[1,2n 1]；
[0069]
步骤5.2、基于无迹卡尔曼滤波器算法的时间更新步完成从第(k-1)时间步到第k时间步的状态量和状态量协方差的更新，分别记作和公式如下所述：
[0070][0071][0072]
式中，和分别为第k时间步第i个sigma点的权重值，为第k时间步第i个sigma点对应的状态量估计值，qk为第k时间步的噪声；
[0073]
步骤5.3、基于无迹卡尔曼滤波器算法的ut变换原理，使用步骤5.2中更新的和生成(2n 1)个sigma点，并通过观测方程求解出每个sigma点对应的观测估计值步骤5.4、基于无迹卡尔曼滤波器算法的量测更新步计算输出第k时间步的量测预测值且
[0074]
式中，为第k时间步第i个sigma点的权重值，为第k时间步第i个sigma点对应的观测估计值；
[0075]
步骤5.5、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的观测误差εk，且
[0076]
式中，yk为第k时间步的观测值，为第k时间步的量测预测值；
[0077]
步骤5.6、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的量测预测协方差p
yy,k
，且
[0078]
式中，为第k时间步第i个sigma点的权重值，为第k时间步第i个sigma点对应的观测估计值，为第k时间步的量测预测值，rk为第k时间步的噪声；
[0079]
步骤5.7、基于步骤5.5及步骤5.6计算的εk和p
yy,k
构造灵敏参数ηk，且并计算输出每步的ηk值；
[0080]
步骤5.8、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步与的互协方差p
xy,k
，
[0081]
步骤5.9、更新第k时间步的卡尔曼增益矩阵：
[0082]
步骤5.10、更新并输出第k时间步的状态量：
[0083]
步骤5.11、更新并输出第k时间步的状态量协方差：
[0084]
步骤5.12、时间步变为(k 1)，重复步骤5.1～步骤5.11直到最大时间步n完成，即直到循环结束。
[0085]
步骤6、绘制步骤5.7中输出的灵敏参数时程曲线，如果曲线整体平稳，无脉冲响应
出现(参见图1)，则无需调用自适应无迹卡尔曼滤波器算法，按常规无迹卡尔曼滤波器算法，即步骤5(但忽略步骤5.5和步骤5.7)识别即可；如果ηk时程曲线有峰值脉冲出现(参见图2)，则需调用自适应无迹卡尔曼滤波器算法(步骤7)识别损伤位置，并令峰值脉冲前出现的最大灵敏参数值等于灵敏参数阈值η0(参见图3，其中图3为图2的局部视图)。
[0086]
步骤7、基于自适应无迹卡尔曼滤波器算法进行识别，过程如下：
[0087]
步骤7.1、同步骤5.1；
[0088]
步骤7.2、同步骤5.2；
[0089]
步骤7.3、同步骤5.3；
[0090]
步骤7.4、同步骤5.4；
[0091]
步骤7.5、同步骤5.5；
[0092]
步骤7.6、同步骤5.6；
[0093]
步骤7.7、同步骤5.7；
[0094]
步骤7.8、同步骤5.8；
[0095]
步骤7.9、同步骤5.9；
[0096]
步骤7.10、同步骤5.10；
[0097]
步骤7.11、同步骤5.11；
[0098]
步骤7.12、判断步骤7.7中计算的ηk与灵敏参数阈值η0的大小，如果ηk＜η0，则时间步变为(k 1)，重复步骤7.1～步骤7.11继续计算；如果ηk≥η0，则继续执行步骤7.13～步骤7.18；
[0099]
步骤7.13、设置一个初始值为0且维度为n的向量，用字母l表示，需要说明的是代表该向量的字符可自选，此处使用字母l表示。
[0100]
步骤7.14、将状态量协方差的主对角线元素中第一个弹性模量参数对应的协方差值的行号或列号记为m，将状态量协方差的主对角线元素中最后一个弹性模量参数对应的协方差值的行号或列号记为l，则总的弹性模量参数个数为(l-m 1)。
[0101]
取步骤7.11计算的的所有主对角元素组成一个新的对角方阵且保持原对角元素的位置不变。依次扩大中每一个弹性模量参数对应的协方差值，令且每次只对一个协方差值进行扩大，的其余元素值保持不变，其中和分别代表行数列数均为z位置处的协方差值，为一个标量，并且式中λ＝1
×
10
ω
，且10
ω
等于初始状态量协方差p0中最小数量级的倒数；z起始值为m，且z∈[m,l-m 1]。依据z＝m时计算的和步骤7.10计算的执行一步完整的无迹卡尔曼滤波运算，即执行步骤7.1～步骤7.11，并输出ηk，并令l(z)＝ηk。
[0102]
步骤7.15、令z＝m 1，继续执行步骤7.14，直到z＝l-m 1时结束。
[0103]
步骤7.16、不计零值，找到l中最小值对应的位置号，记为z
min
，即z
min
为损伤位置。
[0104]
步骤7.17、仅令中的元素扩大20倍，其中代表行数列数均为z位置处的协方差值，且倍数20来源于数值模拟统计结果。
[0105]
步骤7.18、时间步变为(k 1)，重复步骤7.1～步骤7.17直到最大时间步n完成，即
直到循环结束，找到所有损伤位置并输出。
[0106]
实施例
[0107]
为了充分说明本发明，本发明以车-桥耦合系统进行实施例的说明。
[0108]
为了充分说明本发明，本实施例先对图4所示的车-桥耦合系统进行说明：
[0109]
所述的车-桥耦合系统，包括四分之一车辆模型的车体质量1、四分之一车辆模型的车体与轮胎之间的悬挂刚度2、四分之一车辆模型的车体与轮胎之间的悬挂阻尼3、四分之一车辆模型的轮胎质量4、四分之一车辆模型的轮胎与桥梁之间的接触刚度5、轮胎与桥梁接触点6、简支边界约束的固定端7、梁单元8、梁单元节点9和简支边界约束的滑动端10；
[0110]
所述的轮胎与桥梁接触点6是指在车辆行驶过程中轮胎与桥梁始终密接，不发生分离；
[0111]
所述的梁单元8不局限于图中示意的位置，图中一共有6个梁单元；
[0112]
所述的梁单元节点9不局限于图中示意的位置，其余具有相同形状符号的均是梁单元节点，包括桥梁两端位置也都含有梁单元节点。
[0113]
方法实施过程如下：
[0114]
1、基于车辆与结构动力相互作用理论，并结合有限元理论，推导车-桥耦合系统的运动控制微分方程，并基于车辆部分的运动控制微分方程推导接触力关系，即轮胎与桥面板间的相互作用力方程，由此当车辆通过桥梁时可计算出桥梁所受的外荷载激励。
[0115]
2、基于有限元理论，桥梁结构的位移和速度状态可通过各梁单元节点自由度的位移和速度表征，动态荷载作用前，认为各节点自由度的位移和速度都为0，而各梁单元的弹性模量初始值可基于桥梁结构的材料组成求得，由此初始状态量χ0已知，再基于卡尔曼滤波器原理，获得初始状态量协方差p0。
[0116]
3、荷载作用过程中，桥梁各梁单元节点的竖向位移可通过传感器或数值模拟手段采集或计算得到，由此观测值已知。
[0117]
4、基于步骤1～7及上述初始信息可进行桥梁结构损伤位置的识别。
[0118]
为便于本发明的应用说明，此处对算法重点部分做进一步阐释，分别选择初始状态量χ0和中间向量l作进一步说明，重点选择初始状态量是因为它与初始状态量协方差p0的设置、中间向量l的设置以及最后识别的状态量都密切相关。根据前文阐述，本实施例基于欧拉-伯努利梁单元构建桥梁有限元模型，由于每个梁单元有2个节点，每个节点有2个自由度，而根据图4显示，本实施例共有6个梁单元，故共有7个节点，包含14个自由度。因此，本桥梁的位移和速度状态分别通过14个参数表征，再考虑6个弹性模量参数，故本实施例的状态量维度为(14 14 6＝34)。根据方法介绍，l的维度同样等于34，且其初始值都为0。同样，为方便表达识别效果，假设第30和第31位置处的弹性模量参数发生折减损伤，且状态量中第30和第31位置处的弹性模量分别对应桥梁的两个梁单元，具体编号可事先定义好。通过方法计算，可优先识别出位置30和31两处的状态量参数异常，进而通过梁单元编号可迅速确定损伤位置。
[0119]
具体实施方式二：
[0120]
本实施方式为一种存储介质，所述存储介质中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现一种桥梁损伤位置的无损识别方法。
[0121]
应当理解为本实施方式所述的存储介质包括但不限于磁存储介质和光存储介质；
所述磁存储介质包括但不限于ram、rom，以及其他硬盘、u盘等存储介质。
[0122]
具体实施方式三：
[0123]
本实施方式为一种桥梁损伤位置的无损识别设备，所述设备包括处理器和存储器，所述存储器中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现一种桥梁损伤位置的无损识别方法。
[0124]
应当理解为本实施方式所述的设备包括但不限于包括处理器和存储器的设备，还可以包括其他具有信息采集、信息交互、控制能功能的单元或模块所对应的设备，例如所述设备还可以包括信号采集装置等。所述设备包括但不限于pc机、工作站、移动设备等。
[0125]
以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。