CN106909738B - 一种模型参数辨识方法 - Google Patents

一种模型参数辨识方法 Download PDF

Info

Publication number
CN106909738B
CN106909738B CN201710104641.5A CN201710104641A CN106909738B CN 106909738 B CN106909738 B CN 106909738B CN 201710104641 A CN201710104641 A CN 201710104641A CN 106909738 B CN106909738 B CN 106909738B
Authority
CN
China
Prior art keywords
state
value
calculating
time
restoring force
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
CN201710104641.5A
Other languages
English (en)
Other versions
CN106909738A (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.)
Beijing University of Technology
Beijing Technology and Business University
Original Assignee
Beijing University of Technology
Beijing Technology and Business 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 Beijing University of Technology, Beijing Technology and Business University filed Critical Beijing University of Technology
Priority to CN201710104641.5A priority Critical patent/CN106909738B/zh
Publication of CN106909738A publication Critical patent/CN106909738A/zh
Application granted granted Critical
Publication of CN106909738B publication Critical patent/CN106909738B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • 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/06—Power analysis or power optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及一种模型参数辨识方法，包括：选取研究对象变化过程的中间变量和待辨识参数作为系统状态变量并对相关参数进行初始化；建立系统的状态变化模型和恢复力测量模型；根据状态变化模型和各状态变量的初始值计算下一时间点的系统状态预测值；计算每个sigma点的恢复力测量预测值和系统恢复力测量预测值；计算不敏卡尔曼滤波器增益；计算系统状态优化估计值；更新所述状态变化模型。本发明将Bouc‑Wen模型从具体的土木结构模型中抽离，以数据为驱动解决模型参数的准确辨识问题，具有通用性，可拓展应用到各类相关变量间满足或近似满足Bouc‑Wen模型的辨识问题中去。

Description

一种模型参数辨识方法
技术领域
本发明涉及信息科学与土木工程的交叉技术领域，具体涉及一种模型参数辨识方法。
背景技术
由于地震对人类生命和财产安全造成严重威胁，土木结构的抗震性能得到了越来越高的重视。结构混合试验是评估大型复杂土木结构抗震性能的有效试验技术，该试验综合运用了物理加载试验和数值模拟技术，其中数值模拟所依赖的结构数学模型的准确性直接影响着试验结果。与此同时，结构数学模型的一大重要作用是反映大量土木结构存在的迟滞非线性恢复力。
恢复力-形变Bouc-Wen模型通过具有不确定参数的非线性微分方程将结构恢复力与形变相联系，通过参数控制滞回环的形状，以模拟实际工程问题中遇到的非线性现象。Bouc-Wen模型因其光滑的滞回特性和较强的通用性在结构数学模型中得到了广泛的应用。但是，Bouc-Wen模型具有较为复杂的数学关系表达和物理意义不明确的参数，由于难以对其进行高效准确辨识，因此不利于土木结构数学模型的准确表达，进而不利于有效分析结构的地震反应及抗震性能。
另外，现有研究成果和技术方法中，多数将Bouc-Wen模型用于刻画特定结构运动方程中的恢复力，参数与具体结构特征存在极大关联性，导致所运用的参数辨识方法通用性较弱。
发明内容
本发明针对恢复力-形变(Bouc-Wen)模型，在获取不同土木结构形变位移及对应恢复力试验测量数据的基础上，运用不敏卡尔曼滤波器(Unscented Kalman Filter,UKF)对模型中的多个时变参数进行辨识。本发明考虑描述恢复力与形变非线性关系的Bouc-Wen模型：
r＝αku+(1-α)kz (1)
Figure BDA0001232740040000021
其中r为结构对象为抵抗外力作用而产生的恢复力，u、v分别为结构对象在外力作用下产生的位移及相应的速度，z为滞变位移，k为初始结构刚度，α为第二刚度系数(即屈服后刚度与初始刚度之比)，β控制滞回环形状和大小的参数，n控制滞回环光滑度。
本发明解决上述技术问题的技术方案如下：一种模型参数辨识方法，包括：
步骤1、选取研究对象变化过程的中间变量和待辨识参数作为系统状态变量并对相关参数进行初始化；
步骤2、建立系统的状态变化模型和恢复力测量模型；
步骤3、根据状态变化模型和各状态变量的初始值计算下一时间点的系统状态预测值；
步骤4、根据恢复力测量模型和每个sigma点的状态预测值计算每个sigma点的恢复力测量预测值和系统恢复力测量预测值；
步骤5、根据系统状态预测值和系统恢复力测量预测值计算不敏卡尔曼滤波器增益；
步骤6、读取试验测量数据，根据研究对象的恢复力实验测量值、系统恢复力测量预测值和不敏卡尔曼滤波器增益计算系统状态优化估计值；
步骤7、利用结构对象的状态优化估计值中的参数值更新所述状态变化模型。
本发明的有益效果是：将恢复力-形变(Bouc-Wen)模型从具体的土木结构模型中进行抽离，以数据为驱动解决模型参数的准确辨识问题，该方法具有通用性，可拓展应用到各类相关变量间满足或近似满足Bouc-Wen模型的辨识问题中去。
在上述技术方案的基础上，本发明还可以做如下改进。
进一步，所述步骤1具体包括：
选定系统状态变量x＝[z,K,β,n,α]T，其中，z为中间变量，K、β、n和α为待辨识参数，将状态变量进行初始化，初始化方法为：
x(0|0)＝x0+q×rand(0,1)
其中，x0为自定义的状态变量初始值，q为过程噪声的均值，rand(0,1)用于产生取值范围(0,1)内的随机数；
状态估计值的协方差初始值为P(0|0)＝P0，用于产生初始sigma点。
进一步，所述步骤2中，系统状态变化模型为：
x(k)＝f(x(k-1),vel(k-1))+w(k-1)
其中，k为步数，x(k)为tk时刻的系统状态向量，tk为采样时刻，x(k)＝[z(k),K(k),β(k),n(k),α(k)]T为5维状态列向量；x(k-1)为tk-1时刻的系统状态向量；vel(k-1)为tk-1时刻到tk时刻之间结构对象的速度均值；f(·,·)为系统的状态变化方程；w(k-1)为过程噪声，其均值为0，方差为Q；
所述恢复力测量模型为：
r(k)＝h[x(k),u(k)]+v(k)
其中，r(k)为通过监测得到的结构对象在实验过程中下一时间点的恢复力实验测量值，h[x(k),u(k)]为测量方程，x(k)为tk时刻的系统状态向量；u(k)为经测量得到的研究对象在tk时刻由于加载的外力而产生的位移；v(k)为实验中测量恢复力时的高斯测量白噪声，其均值为0，方差为R，且与过程噪声w(k)相互独立。
进一步，所述步骤3具体包括：
步骤3.1、根据各状态变量的初始值计算每个sigma点的状态估计值
Figure BDA0001232740040000041
Figure BDA0001232740040000042
Figure BDA0001232740040000043
其中，
Figure BDA0001232740040000044
为第i个sigma点的状态估计值,i＝0,1,…,2nx，nx代表状态变量的维数；
Figure BDA0001232740040000045
表示在tk-1时刻的系统状态估计值；P(k-1|k-1)表示在tk-1时刻的状态估计值的协方差；
Figure BDA0001232740040000046
为第i个sigma点的测量输出权值(i＝0,1,…,2nx)；
Figure BDA0001232740040000047
为第i个sigma点的测量输出方差权值(i＝0,1,…,2nx)；
Figure BDA0001232740040000048
表示
Figure BDA0001232740040000049
的均方差矩阵的第i列，κ、χ、o为UKF的相关尺度参数，在具体迭代过程中保持不变；
步骤3.2、根据所述系统状态变化模型和每个sigma点的状态估计值
Figure BDA00012327400400000410
计算每一个sigma点的状态预测值：
Figure BDA00012327400400000411
其中，
Figure BDA00012327400400000412
为第i个sigma点基于其在tk-1时刻的状态估计计算得出的tk时刻的状态预测值，
Figure BDA00012327400400000413
为第i个sigma点在tk-1时刻的状态估计值；w(i)为第i个sigma点的过程噪声模拟量，来重构过程噪声对状态预测的影响；vel(k-1)为tk-1时刻到tk时刻之间结构对象的速度均值；
步骤3.3、计算所有sigma点状态预测值的加权值，即系统状态预测值
Figure BDA0001232740040000051
Figure BDA0001232740040000052
的计算公式为：
Figure BDA0001232740040000053
其中，
Figure BDA0001232740040000054
表示基于tk-1时刻结构对象的状态估计得出的tk时刻的状态预测值；
Figure BDA0001232740040000055
为第i个sigma点在tk时刻的状态预测值；
Figure BDA0001232740040000056
为第i个sigma点的状态估计权值(i＝0,1,…,2nx)；
步骤3.4、根据每个sigma点的状态估计值
Figure BDA0001232740040000057
和系统状态预测值
Figure BDA0001232740040000058
计算向前一步预测的估计协方差P(k|k-1)；
P(k|k-1)的计算公式为：
Figure BDA0001232740040000059
其中，P(k|k-1)表示tk时刻基于结构对象的状态预测值得到的状态协方差预测值；
Figure BDA00012327400400000510
为第i个sigma点在tk时刻的状态预测值；
Figure BDA00012327400400000511
表示tk时刻的结构对象状态预测值；
Figure BDA00012327400400000512
为第i个sigma点的状态估计方差权值(i＝0,1,…,2nx)；Q(k-1)为过程噪声协方差。
进一步，所述步骤4具体包括：
步骤4.1、计算每个sigma点的恢复力测量预测值
Figure BDA00012327400400000513
Figure BDA00012327400400000514
其中，
Figure BDA00012327400400000515
为第i个sigma点的恢复力测量预测值；
Figure BDA00012327400400000516
为第i个sigma点的状态预测值；
Figure BDA00012327400400000517
为每个sigma点的状态预测值
Figure BDA00012327400400000518
及结构对象在tk时刻的位移值u(k)代入测量方程得到的结果；v(i)为第i个sigma点的测量噪声模拟量，用来重构测量噪声对测量预测的影响；
步骤4.2、计算所有sigma点的测量预测值的加权值，即系统恢复力测量预测值
Figure BDA00012327400400000622
Figure BDA0001232740040000062
其中，
Figure BDA0001232740040000063
为结构对象在tk时刻的测量预测值；
Figure BDA0001232740040000064
为第i个sigma点的测量预测值；
Figure BDA0001232740040000065
为第i个sigma点的测量输出权值(i＝0,1,…,2nx)。
进一步，所述步骤5具体包括：
步骤5.1、计算目标状态测量预测值的协方差：
Figure BDA0001232740040000066
其中，
Figure BDA0001232740040000067
表示在tk时刻结构对象的测量预测值的协方差；
Figure BDA0001232740040000068
为第i个sigma点的测量预测值；
Figure BDA0001232740040000069
为目标状态测量预测值；
Figure BDA00012327400400000610
为第i个sigma点的测量输出方差权值(i＝0,1,…,2nx)；R(k)为系统的测量方差；
步骤5.2、计算目标状态预测值与目标状态测量预测值的交叉协方差：
Figure BDA00012327400400000611
其中，Pxr(k)为交叉协方差矩阵；
Figure BDA00012327400400000612
为第i个sigma点的状态估计方差权值(i＝0,1,…,2nx)；
Figure BDA00012327400400000613
为第i个sigma点的状态预测值；
Figure BDA00012327400400000614
表示tk时刻的结构对象状态预测值；
Figure BDA00012327400400000615
为第i个sigma点的测量预测值；
Figure BDA00012327400400000616
为目标测量预测值；
步骤5.3、计算不敏卡尔曼滤波器增益：
Figure BDA00012327400400000617
其中，Kf(k)为tk时刻的滤波器增益；Pxr(k)为交叉协方差矩阵；
Figure BDA00012327400400000618
为结构对象的测量预测值的协方差，
Figure BDA00012327400400000619
表示协方差的逆。
进一步，所述步骤6具体包括：
步骤6.1、计算目标状态的修正值X(k)：
Figure BDA00012327400400000620
其中，X(k)为目标状态修正值；Kf(k)为tk时刻滤波器增益，r(k)为tk时刻的测量数据，
Figure BDA00012327400400000621
为tk时刻的测量的估计值；
步骤6.2、计算系统状态的优化估计值
Figure BDA0001232740040000071
Figure BDA0001232740040000072
其中，
Figure BDA0001232740040000073
表示tk时刻的结构对象状态优化估计值；
Figure BDA0001232740040000074
表示基于tk-1时刻结构对象的状态估计得出的tk时刻的对象状态预测值，；X(k)为目标状态修正值；
步骤6.3、计算此时的状态估计值的协方差，计算公式为
Figure BDA0001232740040000075
其中，P(k|k)表示tk时刻结构对象的状态估计值的协方差，P(k|k-1)表示tk时刻的结构对象状态预测值的协方差，Kf(k)为tk时刻的滤波器增益；
Figure BDA0001232740040000076
为交叉协方差矩阵。
进一步，所述步骤7具体包括：利用结构对象的状态优化估计值
Figure BDA0001232740040000077
中的状态分量K、β、n和α代入所述状态变化模型。
进一步，所述方法还包括：
步骤8、判断恢复力测量系统中是否还有未处理的原始测量数据，若是则返回步骤3，否则结束。
附图说明
图1是本发明实施例提供的一种模型参数辨识方法的流程图；
图2是本发明实施例提供的另一种模型参数辨识方法的流程图；
图3是实施例中的参数k辨识结果；
图4是实施例中的参数α辨识结果；
图5是实施例中的参数β辨识结果；
图6是实施例中的参数n辨识结果；
图7是实施例中的中间变量z的辨识结果；
图8是实施例中的恢复力与形变的非线性关系。
具体实施方式
以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。
如图1所示为本发明提供一种模型参数辨识方法100的流程图，该方法100包括：
101、选取研究对象变化过程的中间变量和待辨识参数作为系统状态变量并对相关参数进行初始化；
具体的，选定系统状态变量x＝[z,K,β,n,α]T，其中，为区分后续步骤中的步数k，将式(1)中的k记为K,z为中间变量，K、β、n和α为待辨识参数，将状态变量进行初始化，初始化方法为：
x(0|0)＝x0+q×rand(0,1) (3)
其中，x0为自定义的状态变量初始值，q为过程噪声的均值，rand(0,1)用于产生取值范围(0,1)内的随机数。
状态估计值的协方差初始值为P(0|0)＝P0，用于产生初始sigma点。
102、建立系统的状态变化模型和恢复力测量模型；
具体的，系统状态变化模型为：
x(k)＝f(x(k-1),vel(k-1))+w(k-1) (4)
其中，k为步数，x(k)为tk时刻的系统状态向量，tk为采样时刻，x(k)＝[z(k),K(k),β(k),n(k),α(k)]T为5维状态列向量；x(k-1)为tk-1时刻的系统状态向量；vel(k-1)为tk-1时刻到tk时刻之间结构对象的速度均值；f(·,·)为系统的状态变化方程；w(k-1)为过程噪声，其均值为0，方差为Q；
所述恢复力测量模型为：
r(k)＝h[x(k),u(k)]+v(k) (5)
其中，r(k)为通过监测得到的结构对象在实验过程中下一时间点的恢复力实验测量值，h[x(k),u(k)]为测量方程，x(k)为tk时刻的系统状态向量；u(k)为经测量得到的研究对象在tk时刻由于加载的外力而产生的位移；v(k)为实验中测量恢复力时的高斯测量白噪声，其均值为0，方差为R，且与过程噪声w(k)相互独立。
103、根据状态变化模型和各状态变量的初始值计算下一时间点的系统状态预测值
Figure BDA0001232740040000091
具体的，步骤3具体包括：
步骤3.1、根据各状态变量的初始值计算每个sigma点的状态估计值
Figure BDA0001232740040000092
Figure BDA0001232740040000093
Figure BDA0001232740040000094
其中，
Figure BDA0001232740040000095
为第i个sigma点的状态估计值,i＝0,1,…,2nx，nx代表状态变量的维数；
Figure BDA0001232740040000096
表示在tk-1时刻的系统状态估计值；P(k-1|k-1)表示在tk-1时刻的状态估计值的协方差；
Figure BDA0001232740040000097
为第i个sigma点的测量输出权值(i＝0,1,…,2nx)；
Figure BDA0001232740040000098
为第i个sigma点的测量输出方差权值(i＝0,1,…,2nx)；
Figure BDA0001232740040000099
表示
Figure BDA00012327400400000910
的均方差矩阵的第i列，κ、χ、o为UKF的相关尺度参数，在具体迭代过程中保持不变。
步骤3.2、根据式(4)的系统状态变化模型和每个sigma点的状态估计值
Figure BDA0001232740040000101
计算每一个sigma点的状态预测值：
Figure BDA0001232740040000102
其中，
Figure BDA0001232740040000103
为第i个sigma点基于其在tk-1时刻的状态估计计算得出的tk时刻的状态预测值，
Figure BDA0001232740040000104
为第i个sigma点在tk-1时刻的状态估计值；w(i)为第i个sigma点的过程噪声模拟量，来重构过程噪声对状态预测的影响；vel(k-1)为tk-1时刻到tk时刻之间结构对象的速度均值。
步骤3.3、计算所有sigma点状态预测值的加权值，即系统状态预测值
Figure BDA0001232740040000105
具体的，
Figure BDA0001232740040000106
的计算公式为：
Figure BDA0001232740040000107
其中，
Figure BDA0001232740040000108
表示基于tk-1时刻结构对象的状态估计得出的tk时刻的状态预测值；
Figure BDA0001232740040000109
为第i个sigma点在tk时刻的状态预测值；
Figure BDA00012327400400001010
为第i个sigma点的状态估计权值(i＝0,1,…,2nx)。
步骤3.4、为了衡量向前一步预测的过程中存在的偏差，根据每个sigma点的状态估计值
Figure BDA00012327400400001011
和系统状态预测值
Figure BDA00012327400400001016
计算向前一步预测的估计协方差P(k|k-1)；
具体的，P(k|k-1)的计算公式为：
Figure BDA00012327400400001012
其中，P(k|k-1)表示tk时刻基于结构对象的状态预测值得到的状态协方差预测值；
Figure BDA00012327400400001013
为第i个sigma点在tk时刻的状态预测值；
Figure BDA00012327400400001014
表示tk时刻的结构对象状态预测值；
Figure BDA00012327400400001015
为第i个sigma点的状态估计方差权值(i＝0,1,…,2nx)；Q(k-1)为过程噪声协方差。
104、根据恢复力测量模型和每个sigma点的状态预测值
Figure BDA00012327400400001017
计算每个sigma点的恢复力测量预测值
Figure BDA0001232740040000111
和系统恢复力测量预测值
Figure BDA0001232740040000112
具体的，
Figure BDA0001232740040000113
的计算公式为：
Figure BDA0001232740040000114
其中，
Figure BDA0001232740040000115
为第i个sigma点的恢复力测量预测值；
Figure BDA0001232740040000116
为第i个sigma点的状态预测值；
Figure BDA0001232740040000117
为每个sigma点的状态预测值
Figure BDA0001232740040000118
及结构对象在tk时刻的位移值u(k)代入测量方程得到的结果；v(i)为第i个sigma点的测量噪声模拟量，用来重构测量噪声对测量预测的影响。
计算所有sigma点的测量预测值的加权值，即为系统恢复力测量预测值，计算公式如下：
Figure BDA0001232740040000119
其中，
Figure BDA00012327400400001110
为结构对象在tk时刻的测量预测值；
Figure BDA00012327400400001111
为第i个sigma点的测量预测值；
Figure BDA00012327400400001112
为第i个sigma点的测量输出权值(i＝0,1,…,2nx)。
105、根据系统状态预测值
Figure BDA00012327400400001113
和系统恢复力测量预测值
Figure BDA00012327400400001114
计算不敏卡尔曼滤波器增益Kf(k)。
具体的，为了衡量在预测目标状态的测量的过程中出现的偏差，可计算目标状态测量预测值的协方差为：
Figure BDA00012327400400001115
其中，
Figure BDA00012327400400001116
表示在tk时刻结构对象的测量预测值的协方差；
Figure BDA00012327400400001117
为第i个sigma点的测量预测值；
Figure BDA00012327400400001118
为目标状态测量预测值；
Figure BDA00012327400400001119
为第i个sigma点的测量输出方差权值(i＝0,1,…,2nx)；R(k)为系统的测量方差。
目标状态预测值与目标状态测量预测值的交叉协方差为：
Figure BDA00012327400400001120
其中，Pxr(k)为交叉协方差矩阵；
Figure BDA00012327400400001121
为第i个sigma点的状态估计方差权值(i＝0,1,…,2nx)；
Figure BDA00012327400400001122
为第i个sigma点的状态预测值；
Figure BDA00012327400400001123
表示tk时刻的结构对象状态预测值；
Figure BDA00012327400400001124
为第i个sigma点的测量预测值；
Figure BDA00012327400400001125
为目标测量预测值。
则不敏卡尔曼滤波器增益为
Figure BDA0001232740040000121
其中，Kf(k)为tk时刻的滤波器增益；Pxr(k)为交叉协方差矩阵；
Figure BDA0001232740040000122
为结构对象的测量预测值的协方差，
Figure BDA0001232740040000123
表示协方差的逆。
106、读取试验测量数据，根据研究对象的恢复力实验测量值r(k)、系统恢复力测量预测值
Figure BDA0001232740040000124
和滤波器增益Kf(k)计算系统状态优化估计值
Figure BDA0001232740040000125
具体的，计算目标状态的修正值X(k)：
Figure BDA0001232740040000126
其中，X(k)为目标状态修正值；Kf(k)为tk时刻滤波器增益，r(k)为tk时刻的测量数据，
Figure BDA0001232740040000127
为tk时刻的测量的估计值。
则系统状态的优化估计值为：
Figure BDA0001232740040000128
其中，
Figure BDA0001232740040000129
表示tk时刻的结构对象状态优化估计值；
Figure BDA00012327400400001210
表示基于tk-1时刻结构对象的状态估计得出的tk时刻的对象状态预测值，；X(k)为目标状态修正值。
为了求取在下一时刻用于向前一步预测的sigma点，还需计算此时的状态估计值的协方差。计算公式为
Figure BDA00012327400400001211
其中，P(k|k)表示tk时刻结构对象的状态估计值的协方差，P(k|k-1)表示tk时刻的结构对象状态预测值的协方差，Kf(k)为tk时刻的滤波器增益；
Figure BDA00012327400400001212
为交叉协方差矩阵。
107、利用结构对象的状态优化估计值中的参数值更新步骤2中的状态变化模型；
具体的，利用结构对象的状态优化估计值
Figure BDA00012327400400001213
中的状态分量K、β、n和α代入步骤2中的状态变化模型，即得到在下一时刻用于预测计算的新的状态变化模型；
可选地，在该实施例中，如图2所示，该方法还包括：
108、判断恢复力测量系统中是否还有未处理的原始测量数据，如果有，则返回步骤103；否则结束。
具体的，若此时k的值小于原始测量数据的长度，说明恢复力测量系统中还有未处理的原始测量数据，则返回继续计算系统状态修正值和状态优化估计值，否则结束并输出参数辨识结果。
实施例1：
本发明针对恢复力-形变Bouc-Wen模型进行参数辨识，选用25456组试验测量数据对本发明方法进行仿真计算，本次试验中，时间步长为0.1s，数据格式见表1.
表1实验1测量数据
序号 t(s) u(mm) F(kN) 序号 t(s) u(mm) F(kN)
1 0 0 0 12226 1222.5 -11.19 442.15
2 0.1 2.5 247.2 12227 1222.6 -11.09 444.23
3 0.2 2.5 248.66 12228 1222.7 -10.96 446.3
4 0.3 2.44 249.03 12229 1222.8 -10.8 453.14
5 0.4 2.41 249.15 12230 1222.9 -10.64 460.83
6 0.5 2.5 249.76 12231 1223 -10.47 462.54
7 0.6 2.49 250.86 … … … …
8 0.7 2.47 252.69 25452 2545.1 3.08 -344.49
… … … … 25453 2545.2 3.11 -347.3
12223 1222.2 -11.55 431.78 25454 2545.3 3.1 -352.31
12224 1222.3 -11.37 436.54 25455 2545.4 3.14 -353.77
12225 1222.4 -11.25 440.08 25456 2545.5 3.11 -352.79
考虑试验结构自身特点，将状态变量初始值设置为x0＝[0,276.97,0.55,1,0.025]T，x＝x0+q×rand(0,1)，过程噪声的均值q＝0.001，协方差为
Figure BDA0001232740040000131
测量噪声的均值为1.2，协方差为1.5。估计协方差的初始值设为
Figure BDA0001232740040000141
按具体实施方式中的步骤依次计算系统状态预测值
Figure BDA0001232740040000142
系统状态更新值
Figure BDA0001232740040000143
恢复力测量预测值
Figure BDA0001232740040000144
滤波器增益Kf(k)等，将测量数据进行完全遍历，并在迭代过程中对恢复力-形变Bouc-Wen模型的参数辨识结果进行输出，参数k、α、β、n的辨识结果分别见附图3—附图6，中间变量z的辨识结果见附图7，本实验中恢复力与形变的非线性关系见附图8。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种模型参数辨识方法，用于结构混合试验领域，其特征在于，包括以下步骤：
步骤1、选取研究对象变化过程的中间变量和待辨识参数作为系统状态变量并对相关参数进行初始化；
步骤2、建立系统的状态变化模型和恢复力测量模型；
步骤3、根据状态变化模型和各状态变量的初始值计算下一时间点的系统状态预测值；
步骤4、根据恢复力测量模型和每个sigma点的状态预测值计算每个sigma点的恢复力测量预测值和系统恢复力测量预测值；
步骤5、根据系统状态预测值和系统恢复力测量预测值计算不敏卡尔曼滤波器增益；
步骤6、读取试验测量数据，根据研究对象的恢复力实验测量值、系统恢复力测量预测值和不敏卡尔曼滤波器增益计算系统状态优化估计值；
步骤7、利用结构对象的状态优化估计值中的参数值更新所述状态变化模型；
所述步骤1具体包括：
选定系统状态变量x＝[z,K,β,n,α]T，其中，z为中间变量，K、β、n和α为待辨识参数，K代表结构对象的刚度值，α为第二刚度系数，即屈服后刚度与初始刚度之比，在结构对象的位移与恢复力形成的滞回曲线中，β控制滞回曲线的形状和大小，n决定滞回曲线的光滑度，将状态变量进行初始化，初始化方法为：
x(0|0)＝x0+q×rand(0,1)
其中，x0为自定义的状态变量初始值，q为过程噪声的均值，rand(0,1)用于产生取值范围(0,1)内的随机数；
状态估计值的协方差初始值为P(0|0)＝P0，用于产生初始sigma点；
所述步骤2中，系统状态变化模型为：
x(k)＝f(x(k-1),vel(k-1))+w(k-1)
其中，k为步数，x(k)为tk时刻的系统状态向量，tk为采样时刻，x(k)＝[z(k),K(k),β(k),n(k),α(k)]T为5维状态列向量；x(k-1)为tk-1时刻的系统状态向量；vel(k-1)为tk-1时刻到tk时刻之间结构对象的速度均值；f(·,·)为系统的状态变化方程；w(k-1)为过程噪声，其均值为0，方差为Q；
所述恢复力测量模型为：
r(k)＝h[x(k),u(k)]+v(k)
其中，r(k)为通过监测得到的结构对象在实验过程中下一时间点的恢复力实验测量值，h[x(k),u(k)]为测量方程，x(k)为tk时刻的系统状态向量；u(k)为经测量得到的研究对象在tk时刻由于加载的外力而产生的位移；v(k)为实验中测量恢复力时的高斯测量白噪声，其均值为0，方差为R，且与过程噪声w(k)相互独立；
所述步骤3具体包括：
步骤3.1、根据各状态变量的初始值计算第i个sigma点的状态估计值
Figure FDA0002507723950000021
nx代表状态变量的维数：
Figure FDA0002507723950000022
Figure FDA0002507723950000023
Figure FDA0002507723950000031
其中，
Figure FDA0002507723950000032
表示在tk-1时刻的系统状态估计值；P(k-1|k-1)表示在tk-1时刻的状态估计值的协方差；
Figure FDA0002507723950000033
为第i个sigma点的状态估计权值，i＝0，1，…，2nx；
Figure FDA0002507723950000034
为第i个sigma点的状态估计方差权值，i＝0，1，…，2nx；
Figure FDA0002507723950000035
表示
Figure FDA0002507723950000036
的均方差矩阵的第i列，κ、χ、ο为不敏卡尔曼滤波器的相关尺度参数，在具体迭代过程中保持不变；
步骤3.2、根据状态估计值
Figure FDA0002507723950000037
计算第i个sigma点基于其在tk-1时刻的状态估计计算得出的tk时刻的状态预测值：
Figure FDA0002507723950000038
其中，w(i)为第i个sigma点的过程噪声模拟量，来体现过程噪声对状态预测的影响；vel(k-1)为tk-1时刻到tk时刻之间结构对象的速度均值；
步骤3.3、计算所有sigma点状态预测值的加权值，即基于tk-1时刻结构对象的状态估计得出的tk时刻的系统状态预测值
Figure FDA0002507723950000039
Figure FDA00025077239500000310
的计算公式为：
Figure FDA00025077239500000311
步骤3.4、根据tk时刻所有sigma点的状态估计值和系统状态预测值计算tk时刻基于结构对象的状态预测值得到的状态协方差预测值P(k|k-1)；
P(k|k-1)的计算公式为：
Figure FDA00025077239500000312
其中，Q(k-1)为过程噪声协方差。
2.根据权利要求1所述的方法，其特征在于，所述步骤4具体包括：
步骤4.1、计算第i个sigma点的恢复力测量预测值
Figure FDA00025077239500000313
Figure FDA00025077239500000314
其中，
Figure FDA0002507723950000041
为每个sigma点的状态预测值
Figure FDA0002507723950000042
及结构对象在tk时刻的位移值u(k)代入测量方程得到的结果；v(i)为第i个sigma点的测量噪声模拟量，用来体现测量噪声对测量预测的影响；
步骤4.2、计算所有sigma点的测量预测值的加权值，即结构对象在tk时刻的测量预测值
Figure FDA0002507723950000043
Figure FDA0002507723950000044
3.根据权利要求2所述的方法，其特征在于，所述步骤5具体包括：
步骤5.1、计算在tk时刻结构对象的测量预测值的协方差：
Figure FDA0002507723950000045
其中，R(k)为系统的测量方差；
步骤5.2、计算目标状态预测值与目标状态测量预测值的交叉协方差：
Figure FDA0002507723950000046
步骤5.3、计算tk时刻的不敏卡尔曼滤波器增益：
Figure FDA0002507723950000047
其中，
Figure FDA0002507723950000048
表示协方差的逆。
4.根据权利要求3所述的方法，其特征在于，所述步骤6具体包括：
步骤6.1、计算目标状态的修正值X(k)：
Figure FDA0002507723950000049
其中，X(k)为目标状态修正值；r(k)为tk时刻的测量数据；
步骤6.2、计算tk时刻的结构对象状态优化估计值
Figure FDA00025077239500000413
Figure FDA00025077239500000410
其中，X(k)为目标状态修正值；
步骤6.3、计算tk时刻结构对象的状态估计值的协方差，计算公式为：
Figure FDA00025077239500000411
其中，
Figure FDA00025077239500000412
为交叉协方差矩阵Pxr(k)的转置。
5.根据权利要求4所述的方法，其特征在于，所述步骤7具体包括：利用tk时刻的结构对象的状态优化估计值
Figure FDA0002507723950000051
中的状态分量K、β、n和α代入所述状态变化模型。
6.根据权利要求1-5任一项所述的方法，其特征在于，所述方法还包括：
步骤8、判断恢复力测量系统中是否还有未处理的原始测量数据，若是则返回步骤3，否则结束。
7.根据权利要求6所述的方法，其特征在于，所述步骤8具体包括：判断此时k的值是否小于原始测量数据的长度，若是则返回步骤3，否则结束并输出参数辨识结果。
CN201710104641.5A 2017-02-24 2017-02-24 一种模型参数辨识方法 Active CN106909738B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710104641.5A CN106909738B (zh) 2017-02-24 2017-02-24 一种模型参数辨识方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710104641.5A CN106909738B (zh) 2017-02-24 2017-02-24 一种模型参数辨识方法

Publications (2)

Publication Number Publication Date
CN106909738A CN106909738A (zh) 2017-06-30
CN106909738B true CN106909738B (zh) 2020-07-24

Family

ID=59208126

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710104641.5A Active CN106909738B (zh) 2017-02-24 2017-02-24 一种模型参数辨识方法

Country Status (1)

Country Link
CN (1) CN106909738B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109068138B (zh) * 2018-08-07 2021-12-24 北京市商汤科技开发有限公司 视频图像的处理方法及装置、电子设备和存储介质
CN109885916B (zh) * 2019-02-02 2020-06-16 东南大学 一种基于lssvm的混合试验在线模型更新方法
CN110531708A (zh) * 2019-07-02 2019-12-03 上海大学 一种三自由度压电作动平台的辨识与补偿控制方法
CN112907860B (zh) * 2021-01-18 2022-03-18 南京大学 一种光纤周界安防系统入侵点检测的方法、系统及其装置
CN115685746B (zh) * 2022-09-20 2023-09-01 四川大学 一种机床工作台的离线和在线结合的系统辨识方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5459659A (en) * 1992-05-29 1995-10-17 Honda Giken Kogyo Kabushiki Kaisha Attitude stabilization control system for a legged mobile robot
CN103308896A (zh) * 2013-05-07 2013-09-18 北京工商大学 一种适于非引擎机动目标的高精度跟踪方法
CN103529424A (zh) * 2013-10-23 2014-01-22 北京工商大学 一种基于rfid及ukf实现室内目标快速跟踪的方法
CN103678937A (zh) * 2013-12-29 2014-03-26 中国地震局工程力学研究所 基于等效单自由度体系的钢筋混凝土框架结构整体地震损伤水平评估方法
CN104808662A (zh) * 2015-03-13 2015-07-29 哈尔滨工程大学 一种基于数据驱动的抑制船舶航向扰动的控制方法
US9619491B2 (en) * 2015-04-02 2017-04-11 Sas Institute Inc. Streamlined system to restore an analytic model state for training and scoring

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5459659A (en) * 1992-05-29 1995-10-17 Honda Giken Kogyo Kabushiki Kaisha Attitude stabilization control system for a legged mobile robot
CN103308896A (zh) * 2013-05-07 2013-09-18 北京工商大学 一种适于非引擎机动目标的高精度跟踪方法
CN103529424A (zh) * 2013-10-23 2014-01-22 北京工商大学 一种基于rfid及ukf实现室内目标快速跟踪的方法
CN103678937A (zh) * 2013-12-29 2014-03-26 中国地震局工程力学研究所 基于等效单自由度体系的钢筋混凝土框架结构整体地震损伤水平评估方法
CN104808662A (zh) * 2015-03-13 2015-07-29 哈尔滨工程大学 一种基于数据驱动的抑制船舶航向扰动的控制方法
US9619491B2 (en) * 2015-04-02 2017-04-11 Sas Institute Inc. Streamlined system to restore an analytic model state for training and scoring

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
约束UKF初始参数对Bouc-Wen模型参数识别的影响;王涛等;《黑龙江科技大学学报》;20141130;第24卷(第6期);第651-666页 *

Also Published As

Publication number Publication date
CN106909738A (zh) 2017-06-30

Similar Documents

Publication Publication Date Title
CN106909738B (zh) 一种模型参数辨识方法
KR101958674B1 (ko) 시퀀스 재귀 필터링 3차원 변분(3d-var) 기반의 실측 해양 환경 데이터 동화방법
Sudret Meta-models for structural reliability and uncertainty quantification
US11709979B1 (en) Bridge damage identification method considering uncertainty
CN108875178B (zh) 用于减小结构模态识别不确定性的传感器布置方法
KR100816269B1 (ko) 언센티드 필터를 적용한 강인한 동시 위치 추정 및 지도작성 방법
Seyedpoor et al. A two-step method for damage identification in moment frame connections using support vector machine and differential evolution algorithm
Abudu et al. Modeling of daily pan evaporation using partial least squares regression
CN109583100B (zh) 一种基于ago-rvm的陀螺仪故障预测方法
US9659122B2 (en) Aerodynamic design optimization using information extracted from analysis of unstructured surface meshes
Uchida et al. Parameter identification for multibody systems expressed in differential-algebraic form
CN116956428A (zh) 基于物理信息神经网络的移动荷载下桥梁响应计算方法
Solari The role of analytical methods for evaluating the wind-induced response of structures
CN110489800B (zh) 一种基于矩阵正则化的结构动荷载稀疏识别方法
Friedman et al. A worked-out example of surrogate-based bayesian parameter and field identification methods
CN114692529B (zh) 一种cfd高维响应的不确定度量化方法、装置、计算机设备
CN115809575A (zh) 一种工况传递路径分析方法
CN113158297B (zh) 一种水工弧形钢闸门参数荷载识别方法
Zhu et al. A two-step approach for structural damage localization and quantification using static and dynamic response data
Park et al. Updating structural parameters with spatially incomplete measurements using subspace system identification
CN114580246A (zh) 一种基于非迭代有限元模型修正的桥梁损伤识别方法
Gubarev et al. Identification in cognitive maps in the impulse processes mode with full information
Li et al. Prediction error method‐based second‐order structural identification algorithm in stochastic state space formulation
CN107145665B (zh) 一种巷道围岩应力建模与预测方法
CN112016956A (zh) 基于bp神经网络的矿石品位估值方法及装置

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