CN114692279B - 抗边界约束干扰的时变刚度参数识别方法、存储介质及设备 - Google Patents

抗边界约束干扰的时变刚度参数识别方法、存储介质及设备 Download PDF

Info

Publication number
CN114692279B
CN114692279B CN202210398865.2A CN202210398865A CN114692279B CN 114692279 B CN114692279 B CN 114692279B CN 202210398865 A CN202210398865 A CN 202210398865A CN 114692279 B CN114692279 B CN 114692279B
Authority
CN
China
Prior art keywords
covariance
state quantity
time step
kalman filter
filter algorithm
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210398865.2A
Other languages
English (en)
Other versions
CN114692279A (zh
Inventor
张延哲
丁勇
张立平
卜建清
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Shijiazhuang Tiedao University
Original Assignee
Harbin Institute of Technology
Shijiazhuang Tiedao University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology, Shijiazhuang Tiedao University filed Critical Harbin Institute of Technology
Priority to CN202210398865.2A priority Critical patent/CN114692279B/zh
Publication of CN114692279A publication Critical patent/CN114692279A/zh
Application granted granted Critical
Publication of CN114692279B publication Critical patent/CN114692279B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/20—Administration of product repair or maintenance
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/08—Construction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Human Resources & Organizations (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Tourism & Hospitality (AREA)
  • Strategic Management (AREA)
  • Marketing (AREA)
  • Economics (AREA)
  • General Business, Economics & Management (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Computational Mathematics (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Structural Engineering (AREA)
  • Civil Engineering (AREA)
  • Architecture (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Primary Health Care (AREA)
  • Feedback Control In General (AREA)
  • Filters That Use Time-Delay Elements (AREA)

Abstract

抗边界约束干扰的时变刚度参数识别方法、存储介质及设备，属于桥梁工程健康监测技术领域。为了解决现有的识别方法不能有效识别桥梁结构的时变刚度参数问题。本发明首先基于桥梁损伤位置检测技术识别桥梁结构的损伤位置并编号。然后根据桥梁结构对应的初始状态量及初始状态量协方差，并基于无迹卡尔曼滤波器算法进行初步识别，计算出每步对应的灵敏参数并绘制时程曲线。根据时程曲线是否有峰值脉冲出现，进行进一步识别，当有峰值脉冲出现，调用抗边界约束干扰的自适应无迹卡尔曼滤波器算法进行识别，在此识别过程中针对结构状态量协方差进行二次更新，并更新卡尔曼增益矩阵；进一步更新第k时间步的状态量；直到循环结束获得待识别参数。

Description

抗边界约束干扰的时变刚度参数识别方法、存储介质及设备
技术领域：
本发明属于桥梁工程健康监测技术领域，涉及一种时变刚度参数识别方法、存储介质及设备。
背景技术：
桥梁作为重要的基础设施，其在长期服役过程中，由于重载交通作用、环境侵蚀、材料老化、疲劳效应等使结构损伤不断累积，设计、施工过程中的初始缺陷或日常养护不足更可能导致桥梁结构抗力衰减，致使结构参数具备时变特性，准确地掌握结构参数的时变规律对于结构优化设计、维修加固及灾后救援方案的快速制定等意义重大。
结构参数识别属于反问题，是结构健康监测的核心组成部分。常规参数识别算法一般都假定待识别参数恒定，故很难识别出时变参数的变化规律，且很难收敛到真实参数值。桥梁作为梁式结构，考虑到参数识别这一反问题的可执行性，一般认为其质量恒定，阻尼为质量和刚度的函数，且一般假定为瑞利阻尼。因此，桥梁损伤指标一般通过刚度参数反映。而刚度由弹性模量和截面惯性矩表征，且一般忽略桥梁横截面的变化。由此，更深层次的损伤指标由弹性模量这一微观参数表示。考虑到参数识别问题的时效性及桥梁的尺寸特征，一般基于欧拉-伯努利梁单元构建桥梁的有限元模型，并将此作为真实桥梁结构的合理简化。桥梁结构的位移、速度等响应可通过梁单元节点自由度的对应值表征，且每个梁单元的弹性模量属性都可独立赋予并识别。常规的参数识别算法一般都假定各个梁单元的弹性模量相同且恒定，识别时使用一个参数代表，故很难识别出每个梁单元弹性模量的真实情况。更为真实的识别方式应将所有梁单元的弹性模量同时作为待识别参数，并同步识别，这无疑会增加问题求解的复杂度和困难度，尤其还要面临时变参数这一更具挑战性的问题。
此外，结构的内力和响应只有在足够的边界条件下才能求解，且边界约束位置处的刚度明显大于结构的其它部分。例如，简支梁两端的平动刚度非常大，而其余不受约束部分的平动刚度则取决于结构自身，由此导致结构刚度矩阵的巨大数值差异，再考虑上时变参数的影响，使得常规识别算法在反问题识别时，很难识别出边界位置的桥梁刚度，边界约束条件对参数的正确识别干扰很大。
综上，常规识别算法在识别桥梁结构的损伤时会遇到三个问题，一是无法考虑参数的时变特性，二是无法克服边界条件带来的干扰，三是无法同步识别各个梁段的参数情况。
发明内容：
本发明为了解决现有的识别方法不能有效识别桥梁结构的时变刚度参数问题，进而提出一种抗边界约束干扰的时变刚度参数识别方法。
抗边界约束干扰的时变刚度参数识别方法，包括以下步骤：
针对桥梁结构，确定桥梁结构对应状态的初始值，并组成初始状态量χ0，并依据卡尔曼滤波原理确定初始状态量的协方差矩阵，简称初始状态量协方差P0；其中χ0和P0分别称作第0时间步的状态量和状态量协方差；桥梁结构对应状态包括待识别的参数；
基于无迹卡尔曼滤波器算法进行初步识别，基于无迹卡尔曼滤波器算法进行初步识别的过程中，需要基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的观测误差εk和第k时间步的量测预测协方差Pyy,k，并基于εk和Pyy,k计算输出每步对应的灵敏参数
然后绘制输出灵敏参数ηk的时程曲线，如果ηk时程曲线有峰值脉冲出现，则基于抗边界约束干扰的自适应无迹卡尔曼滤波器算法识别待识别参数；在基于抗边界约束干扰的自适应无迹卡尔曼滤波器算法进行识别的过程中，需要判断计算的ηk与灵敏参数阈值η0的大小，如果ηk＜η0，则基于无迹卡尔曼滤波器算法继续识别；如果ηk≥η0，则继续执行以下步骤：
如果tr为矩阵的迹，则遗忘因子/>否则αk＝1；
基于遗忘因子αk分别修正量测预测协方差Pyy,k、互协方差Pxy,k及状态量协方差再基于遗忘因子αk修正状态量协方差/>并更新卡尔曼增益矩阵；基于更新的卡尔曼增益矩阵更新第k时间步的状态量；同时更新时间步继续滤波，直到循环结束，最终获得识别过程状态量中待识别参数；所述待识别的参数为桥梁各梁单元的弹性模量。
进一步地，基于遗忘因子αk修正状态量协方差的过程包括以下步骤：
将桥梁结构的损伤位置编号用A1、A2、…、As表示，其中1≤s≤总的待识别弹性模量个数；
当s＝1时，修正公式为：
式中，代表行数列数均为A1位置处的协方差值，修正公式含义为：扩大行数列数均为A1位置处的协方差值；/>为根据αk确定的扩大倍数；
当s≥2时，有两种修正方式：
方式①：首先绘制损伤位置排布表，将损伤位置编号分别作为行号和列号，使损伤位置编号两两对应构成一个元素，实际对应状态量协方差中的元素位置；
然后，修正损伤位置排布表中除第一行和第一列外每个元素位置对应的协方差值，公式如下：
式中，x∈[2,s]，y∈[2,s]，代表行数为Ax、列数为Ay位置处的协方差值，修正公式含义为：扩大行数为Ax、列数为Ay位置处的协方差值；
方式②：首先绘制损伤位置排布表，将损伤位置编号分别作为行号和列号，使损伤位置编号两两对应构成一个元素，实际对应状态量协方差中的元素位置；
然后，基于损伤位置排布表选取相应主对角元素，修正公式如下：
式中，x∈[2,s]，代表行数列数均为Ax位置处的协方差值，修正公式含义为：扩大行数列数均为Ax位置处的协方差值。
进一步地，所述扩大倍数取值范围为
进一步地，当2≤s≤3时，优先选择方式①；当s≥4时，优先选择方式②。
进一步地，基于抗边界约束干扰的自适应无迹卡尔曼滤波器算法识别待识别参数的具体过程如下：
步骤7.1、基于有限元理论离散桥梁结构，将桥梁结构划分为单元形式，并编号，确定损伤位置并输出损伤单元的编号；
步骤7.2、基于无迹卡尔曼滤波器算法的UT变换原理，使用第(k-1)时间步的状态量χk-1和状态量协方差Pk-1生成(2n+1)个sigma点，并通过状态方程求解出每个sigma点对应的状态量其中k从1开始，且k∈[1,N]，N为总的时间步数，n为状态量的维数，i为第i个sigma点，且i∈[1,2n+1]；
步骤7.3、基于无迹卡尔曼滤波器算法的时间更新步完成从第(k-1)时间步到第k时间步的状态量和状态量协方差的更新，分别记作和/>公式如下所述：
式中，为第m时间步第i个sigma点的权重值，/>为第c时间步第i个sigma点的权重值，/>为第k时间步第i个sigma点对应的状态量估计值，Qk为第k时间步的噪声；
步骤7.4、基于无迹卡尔曼滤波器算法的UT变换原理，使用步骤7.3中更新的和生成(2n+1)个sigma点，并通过观测方程求解出每个sigma点对应的观测估计值/>
步骤7.5、基于无迹卡尔曼滤波器算法的量测更新步计算输出第k时间步的量测预测值且/>
式中，为第m时间步第i个sigma点的权重值，/>为第k时间步第i个sigma点对应的观测估计值；
步骤7.6、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的观测误差εk，且
式中，yk为第k时间步的观测值，为第k时间步的量测预测值；
步骤7.7、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的量测预测协方差Pyy,k，且
式中，为第c时间步第i个sigma点的权重值，/>为第k时间步第i个sigma点对应的观测估计值，/>为第k时间步的量测预测值；
步骤7.8、基于步骤7.6及步骤7.7计算的εk和Pyy,k构造灵敏参数ηk，且并计算输出每步的ηk值；
步骤7.9、判断步骤7.8中计算的ηk与灵敏参数阈值η0的大小，如果ηk＜η0，则基于步骤5.8～步骤5.12继续识别；如果ηk≥η0，则继续执行步骤7.10～步骤7.15；
步骤7.10、构造遗忘因子αk，如果则/>否则αk＝1；
式中，tr为矩阵的迹；
步骤7.11、基于遗忘因子αk分别修正量测预测协方差Pyy,k、互协方差Pxy,k及状态量协方差并输出状态量协方差，具体公式如下：
步骤7.12、基于遗忘因子αk修正状态量协方差
步骤7.13、基于步骤7.11计算的量测预测协方差Pyy,k和互协方差Pxy,k更新卡尔曼增益矩阵，即
步骤7.14、基于步骤7.13计算的卡尔曼增益矩阵更新并输出第k时间步的状态量，即：
步骤7.15、时间步长变为(k+1)，继续执行步骤7.1～7.14，直到循环结束。
进一步地，利用无迹卡尔曼滤波器算法进行识别的过程中和利用抗边界约束干扰的自适应无迹卡尔曼滤波器算法进行识别的过程中，选择各梁单元节点的竖向位移响应作为观测值y。
进一步地，所述的桥梁结构对应状态包括桥梁各梁单元的弹性模量、桥梁结构的位移和速度。
进一步地，所述桥梁结构对应的运动控制微分方程或有限元模型是基于欧拉-伯努利梁单元建立的。
一种存储介质，所述存储介质中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现所述的抗边界约束干扰的时变刚度参数识别方法。
一种抗边界约束干扰的时变刚度参数识别设备，所述设备包括处理器和存储器，所述存储器中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现所述的抗边界约束干扰的时变刚度参数识别方法。
有益效果：
1、通过引入遗忘因子参数单次修正量测预测协方差、互协方差及两次修正状态量协方差的方式克服了桥梁结构边界约束条件带来的干扰影响，能够准确地同步地识别出各个梁段的弹性模量值，进而准确估计出桥梁结构各梁段的刚度值。可以有效的解决常规参数识别方法无法考虑参数的时变特性，无法克服边界条件带来的干扰以及无法同步识别桥梁各个梁段的参数值的问题。
2、通过时变刚度参数识别方法的应用，能够准确得出桥梁结构时变参数发生的位置及时变参数演变规律，并可同步识别出桥梁其余位置的刚度值值，对于结构优化设计及维修加固等方面意义重大。
3、本发明对识别过程中遗忘因子参数的选取不敏感，可有效避免不同人操作带来的偶然误差。
附图说明：
为了易于说明，本发明由下述的具体实施及附图作以详细描述。
图1为灵敏参数时程曲线图。
图2为灵敏参数时程曲线有峰值脉冲的效果图。
图3为图2的局部视图。
图4为实施例的车-桥耦合系统图，其中车以一定速度从桥的一端行驶到另一端，图中，1-四分之一车辆模型的车体质量，用字母m1表示；2-四分之一车辆模型的车体与轮胎之间的悬挂刚度，用字母k1表示；3-四分之一车辆模型的车体与轮胎之间的悬挂阻尼，用字母c1表示；4-四分之一车辆模型的轮胎质量，用字母m2表示；5-四分之一车辆模型的轮胎与桥梁之间的接触刚度，用字母k2表示；6-轮胎与桥梁接触点；7-简支边界约束的固定端；8-梁单元；9-梁单元节点；10-简支边界约束的滑动端。
图5为图4中单一梁单元损伤时各梁单元的识别效果示例，其中图(a)-(f)分别对应六个梁单元的效果。
图6为图4中两个梁单元损伤时各梁单元的识别效果示例，其中图(a)-(f)分别对应六个梁单元的效果。
图7为图4中四个梁单元损伤时各梁单元的识别效果示例，其中图(a)-(f)分别对应六个梁单元的效果。
具体实施方式：
为使本发明的目的、技术方案和优点更加清楚明了，下面通过附图中示出的具体实施例来描述本发明。但是应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。
在此，还需要说明的是，为了避免因不必要的细节而模糊了本发明，在附图中仅仅示出了与根据本发明的方案密切相关的结构和/或处理步骤，而省略了与本发明关系不大的其它细节。
具体实施方式一：
抗边界约束干扰的时变刚度参数识别方法是一种针对桥梁结构的刚度参数识别方法，其本质上是一种基于抗边界约束干扰的自适应无迹卡尔曼滤波器算法的识别方法；
本实施方式所述的抗边界约束干扰的时变刚度参数识别方法，包括以下步骤：
步骤1、根据达朗贝尔原理、结构动力学及有限元理论，推导结构的运动控制微分方程或建立结构有限元模型，所述结构为桥梁结构，为方便表述，简称结构；
在此过程中，基于欧拉-伯努利梁单元建立桥梁结构的运动控制微分方程或有限元模型，并将各个梁单元的弹性模量作为待识别参数；
步骤2、将结构位移X、结构速度以及弹性模量E组成的向量称作状态量，用符号表示；
根据状态量χ，依据步骤1中结构的运动控制微分方程，并基于线性代数矩阵运算及数值分析的数值微分和积分运算，或者，基于线性代数矩阵运算及结构动力学的Newmark-β法等数值积分方法推导出求解状态量的方程关系；
或者，
根据状态量χ，依据步骤1中结构有限元模型进行状态量的输出设置，并将其作为状态方程；该过程中考虑噪声影响；
步骤3、将传感器采集的各梁单元节点的竖向位移响应作为观测值y。
观测值包含位移、速度、加速度、应变、应力、力、温度等一系列能通过传感器测量的物理量，对于本桥梁结构，使用位移测量值作为观测值。
本发明中，状态方程是状态量的微分运算，状态量写成位移和速度的形式，是方便状态方程的推导，因为对位移求导是速度，对速度求导是加速度，而根据结构的运动控制微分方程很容易就能推出状态方程。另外，当观测值为位移时，因为状态量中的位移和观测值中的位移含义相同，那么观测方程中的位移关系也很容易得到。概括说，状态量中的位移和速度是为状态方程推导以及观测方程中的位移关系服务的，观测值中的位移是为修正观测方程计算的位移服务的。由于此过程的具体内容为公知常识，所以本发明中不再赘述。需要说明的是，本发明中的速度相当于中间量，不需要输出，它的存在只是为状态方程的推导服务的，它的计算由算法递推迭代完成。
根据观测值类型，依据步骤1中结构的运动控制微分方程并基于线性代数矩阵运算、数学移项、数学合并同类项等知识推导出求解观测值的方程关系；
或者，
根据观测值类型，依据步骤1中结构有限元模型进行相应观测值的输出设置，并将其作为观测方程；该过程中考虑噪声影响；
步骤4、将结构的初始位移、结构的初始速度以及各弹性模量的初始值组成的向量称作初始状态量，用符号χ0表示，并根据卡尔曼滤波器算法原理得出初始状态量的协方差矩阵，简称初始状态量协方差，用符号P0表示，其中χ0和P0分别称作第0时间步(启始步)的状态量和状态量协方差；
步骤5、基于无迹卡尔曼滤波器算法计算每步的灵敏参数值，过程如下：
步骤5.1、基于无迹卡尔曼滤波器算法的UT变换原理，使用第(k-1)时间步的状态量χk-1和状态量协方差Pk-1生成(2n+1)个sigma点，并通过状态方程求解出每个sigma点对应的状态量其中k从1开始，且k∈[1,N]，N为总的时间步数，n为状态量的维数，i为第i个sigma点，且i∈[1,2n+1]；
步骤5.2、基于无迹卡尔曼滤波器算法的时间更新步完成从第(k-1)时间步到第k时间步的状态量和状态量协方差的更新，分别记作和/>公式如下所述：
式中，为第m时间步第i个sigma点的权重值，/>为第c时间步第i个sigma点的权重值，/>为第k时间步第i个sigma点对应的状态量估计值，Qk为第k时间步的噪声；
步骤5.3、基于无迹卡尔曼滤波器算法的UT变换原理，使用步骤5.2中更新的和生成(2n+1)个sigma点，并通过观测方程求解出每个sigma点对应的观测估计值/>
步骤5.4、基于无迹卡尔曼滤波器算法的量测更新步计算输出第k时间步的量测预测值且/>
式中，为第m时间步第i个sigma点的权重值，/>为第k时间步第i个sigma点对应的观测估计值；
步骤5.5、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的观测误差εk，且
式中，yk为第k时间步的观测值，为第k时间步的量测预测值；
步骤5.6、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的量测预测协方差Pyy,k，且
式中，为第c时间步第i个sigma点的权重值，/>为第k时间步第i个sigma点对应的观测估计值，/>为第k时间步的量测预测值；
步骤5.7、基于步骤5.5及步骤5.6计算的εk和Pyy,k构造灵敏参数ηk，且并计算输出每步的ηk值；
步骤5.8、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步与/>的互协方差Pxy,k，/>
步骤5.9、更新第k时间步的卡尔曼增益矩阵：Kk＝Pxy,k(Pyy,k+Rk)-1；
式中，Rk为第k时间步的噪声；
步骤5.10、更新并输出第k时间步的状态量：
步骤5.11、更新并输出第k时间步的状态量协方差：
步骤5.12、时间步变为(k+1)，重复步骤5.1～步骤5.11直到最大时间步N完成，即直到循环结束。
步骤6、绘制步骤5.7中输出的灵敏参数时程曲线，如果曲线整体平稳，无脉冲响应出现(参见图1)，则无需调用抗边界约束干扰的自适应无迹卡尔曼滤波器算法，按常规无迹卡尔曼滤波器算法，即步骤5(但忽略步骤5.5和步骤5.7)识别即可；如果ηk时程曲线有峰值脉冲出现(参见图2)，则需调用抗边界约束干扰的自适应无迹卡尔曼滤波器算法(步骤7)，并令峰值脉冲前出现的最大灵敏参数值等于灵敏参数阈值η0(参见图3，其中图3为图2的局部视图)。
步骤7、基于抗边界约束干扰的自适应无迹卡尔曼滤波器算法进行待识别参数的识别，过程如下：
步骤7.1、基于有限元理论离散桥梁结构，将桥梁结构划分为单元形式，并编号，再基于超声波探测技术、声发射探测技术、探地雷达检测技术、灵敏度算法或其它卡尔曼滤波算法等方法识别桥梁结构的损伤位置，并输出损伤单元的编号；
步骤7.2、同步骤5.1；
步骤7.3、同步骤5.2；
步骤7.4、同步骤5.3；
步骤7.5、同步骤5.4；
步骤7.6、同步骤5.5；
步骤7.7、同步骤5.6；
步骤7.8、同步骤5.7；
步骤7.9、判断步骤7.8中计算的ηk与灵敏参数阈值η0的大小，如果ηk＜η0，则基于步骤5.8～步骤5.12继续识别；如果ηk≥η0，则继续执行步骤7.10～步骤7.15；
步骤7.10、构造遗忘因子αk，如果则/>否则αk＝1；
式中，tr为矩阵的迹；
步骤7.11、基于遗忘因子αk分别修正量测预测协方差Pyy,k、互协方差Pxy,k及状态量协方差并输出状态量协方差，具体公式如下：
步骤7.12、提取步骤7.1输出的损伤位置编号，为便于公式书写和理解，用符号A1、A2、…、As表示，其中1≤s≤总的待识别弹性模量个数，基于遗忘因子αk再次修正状态量协方差并输出状态量协方差，具体操作如下，
当s＝1时，修正公式为：
式中，代表行数列数均为A1位置处的协方差值，修正公式含义为：扩大行数列数均为A1位置处的协方差值；/>为扩大倍数，建议取值范围为/>本实施方式中扩大倍数选择为/>
当s≥2时，有两种修正方式，且当2≤s≤3时，优先选择方式①；当s≥4时，优先选择方式②。因为当2≤s≤3时，方式①的精度更高且鲁棒性更好；当s≥4时，方式②的精度更高且便于书写，具体内容如下：
方式①：首先绘制损伤位置排布表，如表1所示，并基于表1确定欲修正元素，其中表1中除第一行和第一列外每个单元格均代表一个元素位置，实际对应状态量协方差中的元素位置；
表1损伤位置排布表
然后，修正表1中除第一行和第一列外每个元素位置对应的协方差值，公式如下：
式中，x∈[2,s]，y∈[2,s]，代表行数为Ax、列数为Ay位置处的协方差值，修正公式含义为：扩大行数为Ax、列数为Ay位置处的协方差值；/>为扩大倍数，建议取值范围为/>本实施方式中扩大倍数选择为/>
方式②：同样基于损伤位置排布表1选取相应主对角元素，修正公式如下：
式中，x∈[2,s]，代表行数列数均为Ax位置处的协方差值，修正公式含义为：扩大行数列数均为Ax位置处的协方差值；/>为扩大倍数，建议取值范围为 本实施方式中扩大倍数选择为/>
步骤7.13、基于步骤7.11计算的量测预测协方差Pyy,k和互协方差Pxy,k更新卡尔曼增益矩阵，即
步骤7.14、基于步骤7.13计算的卡尔曼增益矩阵更新并输出第k时间步的状态量，即：
步骤7.15、时间步长变为(k+1)，继续执行步骤7.1～7.14，直到循环结束。
步骤8、识别结束后，根据识别过程状态量中待识别参数的时程曲线可以判断各弹性模量的变化情况，并能从曲线中得出待识别参数的识别值，通过与桥梁结构相应初始参数进行比较可以及时发现刚度异常并预警。
实施例
为了充分说明本发明，本发明以车-桥耦合系统进行实施例的说明。
为了充分说明本发明，本实施例先对图4所示的车-桥耦合系统进行说明：
所述的车-桥耦合系统，包括四分之一车辆模型的车体质量1、四分之一车辆模型的车体与轮胎之间的悬挂刚度2、四分之一车辆模型的车体与轮胎之间的悬挂阻尼3、四分之一车辆模型的轮胎质量4、四分之一车辆模型的轮胎与桥梁之间的接触刚度5、轮胎与桥梁接触点6、简支边界约束的固定端7、梁单元8、梁单元节点9和简支边界约束的滑动端10；
所述的轮胎与桥梁接触点6是指在车辆行驶过程中轮胎与桥梁始终密接，不发生分离；
所述的梁单元8不局限于图中示意的位置，图中一共有6个梁单元；
所述的梁单元节点9不局限于图中示意的位置，其余具有相同形状符号的均是梁单元节点，包括桥梁两端位置也都含有梁单元节点。
方法实施过程如下：
1、基于车辆与结构动力相互作用理论，并结合有限元理论，可推导车-桥耦合系统的运动控制微分方程，并进一步基于车辆部分的运动控制微分方程可推导接触力的关系方程，即轮胎与桥面板之间的相互作用力，由此桥梁所受的外荷载激励已知。
2、基于有限元理论，桥梁结构的位移和速度状态可通过各梁单元节点自由度的位移和速度表征，动态荷载作用前，认为各节点自由度的位移和速度都为0，而各梁单元的弹性模量初始值可基于桥梁结构的材料组成求得，由此初始状态量χ0已知，再基于卡尔曼滤波原理，获得初始状态量协方差P0；
3、荷载作用过程中，桥梁各梁单元节点的竖向位移可通过传感器或数值模拟手段采集或计算得到，由此观测值已知；
4、基于步骤1～9及以上初始信息可进行桥梁刚度参数的识别，为了更好展示识别效果，此处给出三个应用结果示例，分别对应图5、图6和图7。
具体实施方式二：
本实施方式为一种存储介质，所述存储介质中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现所述的抗边界约束干扰的时变刚度参数识别方法。
应当理解为本实施方式所述的存储介质包括但不限于磁存储介质和光存储介质；所述磁存储介质包括但不限于RAM、ROM，以及其他硬盘、U盘等存储介质。
具体实施方式三：
本实施方式为一种抗边界约束干扰的时变刚度参数识别设备，所述设备包括处理器和存储器，所述存储器中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现所述的抗边界约束干扰的时变刚度参数识别方法。
应当理解为本实施方式所述的设备包括但不限于包括处理器和存储器的设备，还可以包括其他具有信息采集、信息交互、控制能功能的单元或模块所对应的设备，例如所述设备还可以包括信号采集装置等。所述设备包括但不限于PC机、工作站、移动设备等。
以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (8)

1.抗边界约束干扰的时变刚度参数识别方法，包括以下步骤：
针对桥梁结构，确定桥梁结构对应状态的初始值，并组成初始状态量χ0，并依据卡尔曼滤波原理确定初始状态量的协方差矩阵，简称初始状态量协方差P0；其中χ0和P0分别称作第0时间步的状态量和状态量协方差；桥梁结构对应状态包括待识别的参数；
基于无迹卡尔曼滤波器算法进行初步识别，基于无迹卡尔曼滤波器算法进行初步识别的过程中，需要基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的观测误差εk和第k时间步的量测预测协方差Pyy,k，并基于εk和Pyy,k计算输出每步对应的灵敏参数
然后绘制输出灵敏参数ηk的时程曲线，如果ηk时程曲线有峰值脉冲出现，则基于抗边界约束干扰的自适应无迹卡尔曼滤波器算法识别待识别参数；在基于抗边界约束干扰的自适应无迹卡尔曼滤波器算法进行识别的过程中，需要判断计算的ηk与灵敏参数阈值η0的大小，如果ηk＜η0，则基于无迹卡尔曼滤波器算法继续识别；如果ηk≥η0，则继续执行以下步骤：
如果tr为矩阵的迹，则遗忘因子/>否则αk＝1；
基于遗忘因子αk分别修正量测预测协方差Pyy,k、互协方差Pxy,k及状态量协方差
其特征在于：
基于遗忘因子αk分别修正量测预测协方差Pyy,k、互协方差Pxy,k及状态量协方差之后，再基于遗忘因子αk修正状态量协方差/>并更新卡尔曼增益矩阵；基于更新的卡尔曼增益矩阵更新第k时间步的状态量；同时更新时间步继续滤波，直到循环结束，最终获得识别过程状态量中待识别参数；所述待识别的参数为桥梁各梁单元的弹性模量；
其中，基于无迹卡尔曼滤波器算法计算每步的灵敏参数值，过程如下：
步骤5.1、基于无迹卡尔曼滤波器算法的UT变换原理，使用第(k-1)时间步的状态量χk-1和状态量协方差Pk-1生成(2n+1)个sigma点，并通过状态方程求解出每个sigma点对应的状态量其中k从1开始，且k∈[1,N]，N为总的时间步数，n为状态量的维数，i为第i个sigma点，且i∈[1,2n+1]；
步骤5.2、基于无迹卡尔曼滤波器算法的时间更新步完成从第(k-1)时间步到第k时间步的状态量和状态量协方差的更新，分别记作和/>公式如下所述：
式中，为第m时间步第i个sigma点的权重值，/>为第c时间步第i个sigma点的权重值，/>为第k时间步第i个sigma点对应的状态量估计值，Qk为第k时间步的噪声；
步骤5.3、基于无迹卡尔曼滤波器算法的UT变换原理，使用步骤5.2中更新的和/>生成(2n+1)个sigma点，并通过观测方程求解出每个sigma点对应的观测估计值/>
步骤5.4、基于无迹卡尔曼滤波器算法的量测更新步计算输出第k时间步的量测预测值且/>
式中，为第m时间步第i个sigma点的权重值，/>为第k时间步第i个sigma点对应的观测估计值；
步骤5.5、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的观测误差εk，且
式中，yk为第k时间步的观测值，为第k时间步的量测预测值；
步骤5.6、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的量测预测协方差Pyy,k，且
式中，为第c时间步第i个sigma点的权重值，/>为第k时间步第i个sigma点对应的观测估计值，/>为第k时间步的量测预测值；
步骤5.7、基于步骤5.5及步骤5.6计算的εk和Pyy,k构造灵敏参数ηk，且并计算输出每步的ηk值；
步骤5.8、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步与/>的互协方差Pxy,k，/>
步骤5.9、更新第k时间步的卡尔曼增益矩阵：Kk＝Pxy,k(Pyy,k+Rk)-1；
式中，Rk为第k时间步的噪声；
步骤5.10、更新并输出第k时间步的状态量：
步骤5.11、更新并输出第k时间步的状态量协方差：
步骤5.12、时间步变为(k+1)，重复步骤5.1～步骤5.11直到最大时间步N完成，即直到循环结束；
基于抗边界约束干扰的自适应无迹卡尔曼滤波器算法识别待识别参数的具体过程如下：
步骤7.1、基于有限元理论离散桥梁结构，将桥梁结构划分为单元形式，并编号，确定损伤位置并输出损伤单元的编号；
步骤7.2、基于无迹卡尔曼滤波器算法的UT变换原理，使用第(k-1)时间步的状态量χk-1和状态量协方差Pk-1生成(2n+1)个sigma点，并通过状态方程求解出每个sigma点对应的状态量其中k从1开始，且k∈[1,N]，N为总的时间步数，n为状态量的维数，i为第i个sigma点，且i∈[1,2n+1]；
步骤7.3、基于无迹卡尔曼滤波器算法的时间更新步完成从第(k-1)时间步到第k时间步的状态量和状态量协方差的更新，分别记作和/>公式如下所述：
式中，为第m时间步第i个sigma点的权重值，/>为第c时间步第i个sigma点的权重值，/>为第k时间步第i个sigma点对应的状态量估计值，Qk为第k时间步的噪声；
步骤7.4、基于无迹卡尔曼滤波器算法的UT变换原理，使用步骤7.3中更新的和/>生成(2n+1)个sigma点，并通过观测方程求解出每个sigma点对应的观测估计值/>
步骤7.5、基于无迹卡尔曼滤波器算法的量测更新步计算输出第k时间步的量测预测值且/>
式中，为第m时间步第i个sigma点的权重值，/>为第k时间步第i个sigma点对应的观测估计值；
步骤7.6、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的观测误差εk，且
式中，yk为第k时间步的观测值，为第k时间步的量测预测值；
步骤7.7、基于无迹卡尔曼滤波器算法的量测更新步计算第k时间步的量测预测协方差Pyy,k，且
式中，为第c时间步第i个sigma点的权重值，/>为第k时间步第i个sigma点对应的观测估计值，/>为第k时间步的量测预测值；
步骤7.8、基于步骤7.6及步骤7.7计算的εk和Pyy,k构造灵敏参数ηk，且并计算输出每步的ηk值；
步骤7.9、判断步骤7.8中计算的ηk与灵敏参数阈值η0的大小，如果ηk＜η0，则基于步骤5.8～步骤5.12继续识别；如果ηk≥η0，则继续执行步骤7.10～步骤7.15；
步骤7.10、构造遗忘因子αk，如果则/>否则αk＝1；
式中，tr为矩阵的迹；
步骤7.11、基于遗忘因子αk分别修正量测预测协方差Pyy,k、互协方差Pxy,k及状态量协方差并输出状态量协方差，具体公式如下：
步骤7.12、再基于遗忘因子αk修正状态量协方差
步骤7.13、基于步骤7.11计算的量测预测协方差Pyy,k和互协方差Pxy,k更新卡尔曼增益矩阵，即
步骤7.14、基于步骤7.13计算的卡尔曼增益矩阵更新并输出第k时间步的状态量，即：
步骤7.15、时间步长变为(k+1)，继续执行步骤7.1～7.14，直到循环结束；
再基于遗忘因子αk修正状态量协方差的过程包括以下步骤：
将桥梁结构的损伤位置编号用A1、A2、…、As表示，其中1≤s≤总的待识别弹性模量个数；
当s＝1时，修正公式为：
式中，代表行数列数均为A1位置处的协方差值，修正公式含义为：扩大行数列数均为A1位置处的协方差值；/>为根据αk确定的扩大倍数；
当s≥2时，有两种修正方式：
方式①：首先绘制损伤位置排布表，将损伤位置编号分别作为行号和列号，使损伤位置编号两两对应构成一个元素，实际对应状态量协方差中的元素位置；
然后，修正损伤位置排布表中除第一行和第一列外每个元素位置对应的协方差值，公式如下：
式中，x∈[2,s]，y∈[2,s]，代表行数为Ax列数为Ay位置处的协方差值，修正公式含义为：扩大行数为Ax、列数为Ay位置处的协方差值；
方式②：首先绘制损伤位置排布表，将损伤位置编号分别作为行号和列号，使损伤位置编号两两对应构成一个元素，实际对应状态量协方差中的元素位置；
然后，基于损伤位置排布表选取相应主对角元素，修正公式如下：
式中，x∈[2,s]，代表行数列数均为Ax位置处的协方差值，修正公式含义为：扩大行数列数均为Ax位置处的协方差值。
2.根据权利要求1所述的抗边界约束干扰的时变刚度参数识别方法，其特征在于：所述扩大倍数取值范围为
3.根据权利要求2所述的抗边界约束干扰的时变刚度参数识别方法，其特征在于：当2≤s≤3时，选择方式①；当s≥4时，选择方式②。
4.根据权利要求1所述的抗边界约束干扰的时变刚度参数识别方法，其特征在于：利用无迹卡尔曼滤波器算法进行识别的过程中和利用抗边界约束干扰的自适应无迹卡尔曼滤波器算法进行识别的过程中，选择各梁单元节点的竖向位移响应作为观测值y。
5.根据权利要求4所述的抗边界约束干扰的时变刚度参数识别方法，其特征在于：所述的桥梁结构对应状态包括桥梁各梁单元的弹性模量、桥梁结构的位移和速度。
6.根据权利要求5所述的抗边界约束干扰的时变刚度参数识别方法，其特征在于：所述桥梁结构对应的运动控制微分方程或有限元模型是基于欧拉-伯努利梁单元建立的。
7.一种存储介质，其特征在于：所述存储介质中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现权利要求1至6之一所述的抗边界约束干扰的时变刚度参数识别方法。
8.一种抗边界约束干扰的时变刚度参数识别设备，其特征在于：所述设备包括处理器和存储器，所述存储器中存储有至少一条指令，所述至少一条指令由处理器加载并执行以实现权利要求1至6之一所述的抗边界约束干扰的时变刚度参数识别方法。
CN202210398865.2A 2022-04-15 2022-04-15 抗边界约束干扰的时变刚度参数识别方法、存储介质及设备 Active CN114692279B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210398865.2A CN114692279B (zh) 2022-04-15 2022-04-15 抗边界约束干扰的时变刚度参数识别方法、存储介质及设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210398865.2A CN114692279B (zh) 2022-04-15 2022-04-15 抗边界约束干扰的时变刚度参数识别方法、存储介质及设备

Publications (2)

Publication Number Publication Date
CN114692279A CN114692279A (zh) 2022-07-01
CN114692279B true CN114692279B (zh) 2023-09-15

Family

ID=82142083

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210398865.2A Active CN114692279B (zh) 2022-04-15 2022-04-15 抗边界约束干扰的时变刚度参数识别方法、存储介质及设备

Country Status (1)

Country Link
CN (1) CN114692279B (zh)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20160110773A (ko) * 2015-03-12 2016-09-22 숙명여자대학교산학협력단 칼만 필터 모델을 이용한 이동 객체 추적 방법 및 장치
KR101907589B1 (ko) * 2018-01-22 2018-10-12 연세대학교 산학협력단 확장 칼만 필터와 유전자 알고리즘을 이용한 건축 구조물 시스템 식별 장치 및 그 방법
CN110874450A (zh) * 2019-11-20 2020-03-10 武汉理工大学 一种基于车载监测的铁路桥梁轨道不平顺计算方法
CN114036605A (zh) * 2021-10-29 2022-02-11 河海大学 基于自适应控制的卡尔曼滤波钢桁架桥梁结构参数监测方法
CN114186595A (zh) * 2021-12-14 2022-03-15 哈尔滨工业大学 时变结构参数识别方法、存储介质及设备
CN114322911A (zh) * 2021-12-31 2022-04-12 重庆大学 一种联合卡尔曼滤波的桥梁路面平整度间接精准识别方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20160110773A (ko) * 2015-03-12 2016-09-22 숙명여자대학교산학협력단 칼만 필터 모델을 이용한 이동 객체 추적 방법 및 장치
KR101907589B1 (ko) * 2018-01-22 2018-10-12 연세대학교 산학협력단 확장 칼만 필터와 유전자 알고리즘을 이용한 건축 구조물 시스템 식별 장치 및 그 방법
CN110874450A (zh) * 2019-11-20 2020-03-10 武汉理工大学 一种基于车载监测的铁路桥梁轨道不平顺计算方法
CN114036605A (zh) * 2021-10-29 2022-02-11 河海大学 基于自适应控制的卡尔曼滤波钢桁架桥梁结构参数监测方法
CN114186595A (zh) * 2021-12-14 2022-03-15 哈尔滨工业大学 时变结构参数识别方法、存储介质及设备
CN114322911A (zh) * 2021-12-31 2022-04-12 重庆大学 一种联合卡尔曼滤波的桥梁路面平整度间接精准识别方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Structural Parameters Identification Based on Extended Kalman Filter;Bin Zeng等;《2019 International Conference on Artificial Intelligence and Advanced Manufacturing (AIAM)》;第566-570页 *
基于车桥耦合振动信号和扩展卡尔曼滤波的桥梁结构损伤识别;黄进鹏等;《土木工程与管理学报》;第35卷(第05期);第140-144页 *
车辆—轨道系统参数识别方法;张延哲;《中国优秀硕士学位论文全文数据库信息科技辑》(第2期);第1-83页 *

Also Published As

Publication number Publication date
CN114692279A (zh) 2022-07-01

Similar Documents

Publication Publication Date Title
CN107729706B (zh) 一种非线性机械系统的动力学模型构建方法
CN110704801B (zh) 桥梁集群结构运营安全智能监测与快速检测成套方法
CN103162678B (zh) 批量mems陀螺信息融合方法
CN111324949A (zh) 一种考虑噪声影响的工程结构柔度识别方法
CN106503730A (zh) 一种基于级联字典与稀疏正则化的桥梁移动荷载识别方法
CN109902408B (zh) 一种基于数值运算和改进的正则化算法的载荷识别方法
CN103884776B (zh) 一种提高随机损伤定位向量法监测结果准确性的方法
CN107862170B (zh) 一种基于动态缩聚的有限元模型修正方法
CN112949131B (zh) 连续桥梁集群损伤诊断的概率损伤定位向量法
CN107679328A (zh) 一种系统参数识别的最优传感器布置方法
Hung et al. Structural damage detection using the optimal weights of the approximating artificial neural networks
CN114925526B (zh) 一种结合多工况响应的结构模态参数识别方法
CN111595541A (zh) 一种海量振动传递率数据卷积神经网络处理的多维结构损伤识别方法
CN108038315A (zh) 一种基于谱随机有限元模型的随机动载荷识别方法
CN110362902B (zh) 一种基于区间逐维分析的单源动载荷识别方法
CN112182697A (zh) 一种有阻尼吊杆系统张力的高精度动测法
CN114692279B (zh) 抗边界约束干扰的时变刚度参数识别方法、存储介质及设备
CN115128300A (zh) 传感器优化布置下的结构动态载荷/参数联合识别方法
CN114692465B (zh) 桥梁损伤位置的无损识别方法、存储介质及设备
CN111274630A (zh) 一种用于工程结构柔度识别的物理模态提取方法
CN106383003A (zh) 基于柔度识别的索结构索力的测量方法及测量系统
CN111027133B (zh) 一种基于径向基神经网络的结构动态分布载荷识别方法
CN105651537B (zh) 一种高损伤敏感性的桁架结构损伤实时监测系统
CN105808845A (zh) 基于萤火虫群智能算法的结构损伤识别方法
CN113688465B (zh) 一种基于载荷与状态结合的飞行器结构强度数字孪生方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant