CN106897717A - 基于环境激励数据的多次测试下贝叶斯模型修正方法 - Google Patents

Abstract

本发明目的在于克服传统方法的缺点，给出一种基于环境激励数据的多次测试下贝叶斯模型修正方法，可实现对多次测试数据进行直接的处理分析，可对多次测试得到的模态参数进行一次性输入，模型修正结果直接输出。本发明技术方案可用来解决基于实际测试数据的有限元模型修正问题。本发明分两个阶段，第一阶段是对多次测试下采集的环境激励下结构加速度数据进行分析，得到每次测试测得的结构的固有频率和振型，并计算这些模态参数的不确定性，用协方差矩阵来表示。第二个阶段是基于多次测试得到的结构模态参数及其协方差矩阵，基于贝叶斯理论构建目标函数，通过对目标函数的优化，得到需要修正的有限元模型的模型参数的最优值。

Description

基于环境激励数据的多次测试下贝叶斯模型修正方法

技术领域

本发明涉及结构有限元模型修正技术。

背景技术

在土木领域做结构设计的时候，通常都会建立有限元模型，而结构基于设计图纸建完之后，结构的一些固有模态参数包含固有频率，阻尼比，振型等与基于有限元建模计算得到的结果有明显的区别。模型修正技术就是基于实测数据识别得到的结构模态参数来对有限元模型进行修正，从而得到更为准确的有限元模型，为之后的结构健康监测、损伤识别等服务。

现有的技术有以下两个问题：

第一个问题是在实际振动测试过程中，传感器的数目往往少于需要测试的测点的数目，然而目前尚无直接基于多次测试数据直接进行模型修正的方法，传统方法需要进行多次独立的数据处理，因而易产生误差；或者只能进行单次测试的模型修正，但受到传感器数目的限制。

第二个问题是基于环境激励下的振动测试，输入激励是随机激励，因此输出的模态参数具有一定的误差并存在不确定性。传统技术只能直接利用识别出的模态参数的值本身，而对模态参数存在的误差及不确定性无法获得，从而也就不能进行利用。

另外，现有的基于振动数据的模型修正技术，通常是基于固有频率和振型两种主要参数建立目标函数，但是在目标函数中如何确定这两种参数的权重是目前尚无法合理解决的问题，传统方法往往通过经验来确定。

发明内容

本发明目的在于克服传统方法的缺点，给出一种基于环境激励数据的多次测试下贝叶斯模型修正方法，可实现对多次测试数据进行直接的处理分析，可对多次测试得到的模态参数进行一次性输入，模型修正结果直接输出。本发明可以基于计算得到的两种模态参数(固有频率和振型)的不确定性来得到它们在目标函数中的权重系数，从而从根本上得到了目标参数确定方法，无需人工经验。

本发明技术方案可用来解决基于实际测试数据的有限元模型修正问题。

为此，本发明需要保护的技术方案表征为：

一种基于环境激励数据的多次测试下贝叶斯模型修正方法，其特征在于，分两个阶段，

第一阶段是对多次测试下采集的环境激励下结构加速度数据进行分析，得到每次测试测得的结构的固有频率和振型，并计算这些模态参数的不确定性，用协方差矩阵来表示。

第二个阶段是基于多次测试得到的结构模态参数及其协方差矩阵，基于贝叶斯理论构建目标函数，通过对目标函数的优化，得到需要修正的有限元模型的模型参数的最优值。

所述的基于环境激励数据的多次测试下贝叶斯模型修正方法，其特征在于，总体构建方法如下：

从结构动力学的基本原理出发，考虑一个线弹性的结构满足以下的动力方程：

这里M,C,K分别表示结构的质量，阻尼和刚度矩阵，W是外力向量。假设该结构满足经典阻尼，结构的加速度可以从下式得到：

这里，u_i是第i阶全振型向量，是第i阶模态的模态加速度响应。刚度质量的关系可以通过以下特征方程得到：

这里ω_i表示结构的第i阶固有频率。让θ表示与结构的刚度矩阵K和质量矩阵M相关的结构参数。已知刚度和质量矩阵，结构的固有频率和全振型理论上通过特征值分解得到。因此，构建一个理论模型来进行模型修正从而确定θ。

让D＝{D_i:i＝1,...,n_s}表示多次测试得到的用来进行结构模态识别的数据，其中D_i表示第i次测试得到的数据。基于两阶段的模型修正公式和多次测试数据，得到结构参数θ的后验分布：

其中，p(θ)表示结构参数的先验分布；由固有频率和部分振型组成。由于可以通过有限元模型得到，其提供了在模型修正过程中第一阶段和第二阶段相互关联的以下信息。条件概率密度函数表示给定结构模型参数的条件下，结构模态参数的先验概率分布；表示综合了多次测试数据的的边缘后验分布，这里在第一阶段的先验分布被认为是均匀分布。假设有限元模型在预测结构模态参数的过程中不存在模型误差，那么条件概率密度函数可以通过一个Dirac-Delta方程得到：

这里，

其中，和分别表示固有频率和振型的理论解，它们可以通过解特征方程得到。

基于以上的推导的公式，当忽略模型误差时，p(θ|D)可以表示为只与有关。为了构建二阶段模型修正公式，接下来后验概率密度函数将通过利用环境激励下多次测试数据信息得到。

公式(4)中的后验概率密度函数公式即为该发明的总体框架公式，包含两个阶段，即第一阶段：贝叶斯模态识别；第二阶段：基于第一阶段得到的多次测试模态参数，进行贝叶斯模型修正。

所述的基于环境激励数据的多次测试下贝叶斯模型修正方法，其特征在于，第一阶段-贝叶斯模态识别,具体实现方法如下：

2.1数据采集

采集数据时，将加速度或者速度传感器放在结构上，结构的激励可来自于周围的风荷载、交通荷载、环境噪音及结构中人员活动等。在传感器数目少于需要测试的测点数目时，通过多次测试完成。多次测试需要设置参考点，参考点位置以能采集到尽可能多的有效模态为原则。每次测试时间建议大于第一个周期长度的600倍。每次测试数据时间长度尽量保持一致。

2.2多次测试下贝叶斯模态识别目标函数构建

采集到的多次测试数据，分别进行单次测试数据的模态识别完成，模态识别分两部分进行，首先基于贝叶斯方法进行模态参数最优值的识别，然后进行模态参数后验不确定性的计算。将单次模态识别得到的模态参数进行收集，用于后期的模型修正。

贝叶斯模态识别方法基本原理是要识别的模态的快速傅里叶变换数据在对应的频域段内可以很好的近似为一个高斯概率密度函数。通过最大化这个高斯分布函数，从而可以将模态参数得到。该方法简单概述如下：

在第i个测试的加速度数据可以近似的模拟为：

其中是i次测试的理论加速度响应，该响应是通过将要识别的模态参数来构建的。这些模态参数包括固有频率，阻尼比，模态力的功率谱密度，预测误差的功率谱密度以及振型等。在公式(7)中，表示模型误差，N_i表示样本的数目，n_i表示单次测试自由度的数目。测试数据的快速傅里叶变化可以定义为：

这里，i²＝-1；Δt_i表示i次测试的样本时间间隔；k＝1,...,N_qi；N_qi＝int[N_i/2]+1是奈奎斯特频率的频率指标，int[.]表示整数部分。在i次测试中用来模态识别的数据D_i可以表示为

其中是在i次测试的快速傅里叶变换数据{F_ik}在第r个频域段数据的集合。n_B表示选择的频域段的数目。可以完全确定的概率分布的模态参数可以表示为：

其中

分别表示在r个频域段固有频率和阻尼比的集合；是模态力的功率谱密度，其可以在一个频域段内假设为一个常数；是预测误差的功率谱密度，其也可以在一个频域段内假设为一个常数。同时，

其中表示在第i次测试下第r个频域段的第j阶振型。

基于贝叶斯定理，给定第i次测试数据，的后验概率密度函数可以得到：

其中表示先验概率分布。假设先验信息满足均匀分布，先验概率密度函数可以认为是一个常数。因此后验概率密度函数可以认为直接跟似然函数成正比。当N_i足够大及Δt_i足够小时，不同频率的快速傅里叶变换可以证明其是近似独立的，同时他们的实部和虚部被证明满足高斯分布。因此似然函数可以写为：

其中表示负对数似然函数，其可以通过以下公式得到：

这里‘*’表示复数的共轭转置；

是在频率f_k理论时的理论功率谱密度矩阵；是一个单位矩阵；表示在r个频域段的模态正定转换矩阵，其(p,q)单元可以从下式得到：

其中

公式(16)是贝叶斯模态识别的目标函数，后续的最优值可以通过最小化负对数似然函数来实现。模态参数的协方差矩阵可以通过使其等于目标函数(16)的汉森矩阵的逆来实现。

2.3算法实现：

通过MATLAB编程来实现以上贝叶斯模态识别方法，程序收敛后，可以识别模态参数包括固有频率，振型，阻尼比，模态力的功率谱及预测误差的功率谱。若程序不收敛，需要重新选择频率段，进行重复循环。其中固有频率和振型以及其对应的参数的协方差矩阵将用到后期的第二阶段的贝叶斯模型修正。

总之，通过第一阶段的步骤可完成以下工作：

对于每一次测试，i＝1,...,n_s，对于每一个频域段，通过贝叶斯模态识别优化计算模态参数最优值：和其对应的后验协方差矩阵

所述的基于环境激励数据的多次测试下贝叶斯模型修正方法，其特征在于，第二阶段：贝叶斯模型修正，具体流程：

3.1构建目标结构的有限元模型

建立目标结构的有限元模型，在后续的模型修正过程中，直接进行调用。

3.2输入模态参数及输出模型参数的选择。

选定第一阶段获得的需要输入的多次测试得到的多组模态参数，尽可能将识别得到的所有模态用到，以提供最多的有效信息。同时根据有限元模型，选定需要修正的模型参数，模型参数数目需要根据测点信息及输入模态信息相对应，避免输出过多参数，导致后期优化过程中不收敛。

3.3构建目标函数并优化(贝叶斯模型修正后验概率密度函数)

在这个部分，定义选择矩阵，其可以将全局振型和在单次测试下得到的振型关联起来，从而基于多次测试数据构建模型参数的后验概率密度函数。这里首先介绍如何得到后验概率密度函数，然后介绍如何将负对数似然函数进行重构，从而方便进行优化。

3.3.1选择矩阵

全局振型Φ^(r)可以通过定义一个选择矩阵L_i来将其与i次测试时得到的振型关联。这个矩阵中，当自由度s在第r频道被测到，那么(r,s)对应的数值就等于1，其他值等于0.第i次测试的振型可以从以下公式得到：

为了方便，这个方法中，假设第i次测试的振型向量正则化为1。

3.3.2多次测试下模型参数的后验概率密度函数

让α＝{α_i,i＝1,...,n_s}表示所有测试下的模态参数。基于贝叶斯理论，给定所有测试的数据，α的后验概率密度函数可以通过下式得到：

给定α，假设在多次测试下数据在统计上是独立的，因此

这里应该注意到p(D_i|α)与其他测试时的参数无关，因此

p(D_i|α)＝p(D_i|α_i) (22)

从而

其中，

这里由i次测试得到的固有频率和部分振型组成

其中f_i和Φ_i分别表示在i次测试下所有选择的频率段内所有频率和阻尼比。参数υ_i由i次测试下剩下的其他模态参数组成，

υ_i＝{ζ_i,S_i,S_ei} (26)

其中ζ_i,S_i和S_ei分别表示在i次测试下所有选择的频率段内阻尼比，模态力的功率谱密度和预测误差的功率谱密度。

因此，基于贝叶斯定理，公式(23)可以由下式得到：

因为p(D)和p(D_i)可以认为是常数，所以公式(27)可以重新写为：

假设先验信息为均匀分布，可以得到：

因此，在第i次测试时，第一阶段的模态参数的后验概率密度函数p₀(α_i|D_i)可以从下式得到：

其中可以通过公式(16)得到。

假设每个是在全局范围内可识别的，在i次测试下，在公式(30)中的每个的后验概率密度函数可以很好的近似为一个高斯分布，其均值为最大可能值协方差矩阵为识别的模态参数协方差矩阵其分布可以写为：

在i次测试下，的边缘后验概率分布函数仍然是一个高斯分布，因此

其中和分别为的最优值和协方差矩阵，其可以从对应的和中的某一部分提取。

考虑多次测试下，基于公式(29)，我们可以得到：

其中

同时

假设固有频率和振型可以完全由结构模型参数决定，将(5)和(33)代入(4)，后验概率密度函数p(θ|D)可以表示为：

其中

这里表示在第r个频域段由有限元模型计算得到的固有频率，其中表示由有限元模型计算得到的对应测试自由度的振型。

3.3.3负对数似然函数的重构

由于振型存在着范数约束，在公式(37)中计算时会出现数值计算问题，因此非常有必要在计算过程中通过计算矩阵的特征基来克服这个问题。经过重构，公式(37)可以写为：

这里和分别是在i次测试下的r频域段中的汉森矩阵的特征值和特征向量。通过重构，不需要计算任何矩阵的逆。

基于目标函数(39)，通过输入模态参数及其协方差矩阵编写程序，优化使其达到最小值。若程序收敛，可以得到模型修正参数θ的最优值；若程序不收敛，那么需要回到开始的地方，调整有限元模型及选择模型修正参数进行循环计算，直至程序收敛。

3.4结构模型参数不确定性计算

在二次泰勒近似的情况下，当θ达到最优值时，后验协方差矩阵可以通过计算负对数似然函数的汉森矩阵的逆来得到，该汉森矩阵可以通过有限差分法来得到。从而我们可以实现评估得到的模型参数的不确定性。

与现有技术相比，该方法主要有以下几个优点：

1)本发明的技术比传统方法更加便捷，可以实现对多次测试数据进行整合，从而直接进行处理分析，实现多次测试的模态参数一次性输入，模型修正结果直接输出，克服了传统方法步骤繁琐的缺点。同时，由于减少了操作步骤，从而减少计算过程中的误差。

2)本发明目标函数中，关于频率和振型两种参数的权重可以通过识别的模态参数的协方差矩阵来确定，从而克服了传统方法确定目标函数权重指数需要人工经验确定的缺点，从理论根本上解决了这个问题。识别的模型参数比传统方法更加准确。

3)本发明与传统技术相比，可以用较少的传感器，通过设置参考点，来测量在实际模型修正过程中需要的大量测点，由于传感器往往由于价格昂贵，而本发明节约了传感器使用数目，因而其带来了一定的经济效益。

4)本发明目标函数构建中考虑了模型误差和模态参数识别误差的双重影响，而传统方法往往只能考虑其中之一的影响。

附图说明

图1：方法的总体框架图

图2：第一阶段贝叶斯模态识别流程图

图3：第二阶段贝叶斯模型修正流程图

具体实施方式

本发明是基于结构动力学的基本原理，将实测数据与有限元得到的结构模态参数经过推导，构建后验概率密度函数(目标函数)。该函数包含多次测试的模态分析结果及识别模态参数的不确定性等，可以从理论上推导得到在目标函数中固有频率和振型两种模态参数之间的关系。

该方法的总体框架图如图1所示，其分为两个阶段。本发明方法分两个阶段，

第一阶段是对多次测试下采集的环境激励下结构加速度数据进行分析，得到每次测试测得的结构的固有频率和振型，并计算这些模态参数的不确定性，用协方差矩阵来表示。

第二个阶段是基于多次测试得到的结构模态参数及其协方差矩阵，基于贝叶斯理论构建目标函数，通过对目标函数的优化，得到需要修正的有限元模型的模型参数的最优值，同时可基于有限差分法计算模型参数的不确定性

以下结合附图，对本发明方法过程详述如下。

1.两阶段模型修正方法总体框架构建如图1所示，总体构建方法如下：

该方法从结构动力学的基本原理出发，考虑一个线弹性的结构满足以下的动力方程：

这里M,C,K分别表示结构的质量，阻尼和刚度矩阵，W是外力向量。假设该结构满足经典阻尼，结构的加速度可以从下式得到：

这里，u_i是第i阶全振型向量，是第i阶模态的模态加速度响应。刚度质量的关系可以通过以下特征方程得到：

这里ω_i表示结构的第i阶固有频率。让θ表示与结构的刚度矩阵K和质量矩阵M相关的结构参数。如果刚度和质量矩阵都知道，那么结构的固有频率和全振型就都可以理论上通过特征值分解得到。因此，通过这个基本原理，我们可以构建一个理论模型来进行模型修正从而确定θ。

让D＝{D_i:i＝1,...,n_s}表示多次测试得到的用来进行结构模态识别的数据，其中D_i表示第i次测试得到的数据。基于两阶段的模型修正公式和多次测试数据，我们可以得到结构参数θ的后验分布：

其中，p(θ)表示结构参数的先验分布；由固有频率和部分振型组成。由于可以通过有限元模型得到，其提供了在模型修正过程中第一阶段和第二阶段相互关联的信息。条件概率密度函数表示给定结构模型参数的条件下，结构模态参数的先验概率分布；表示综合了多次测试数据的的边缘后验分布，这里在第一阶段的先验分布被认为是均匀分布。假设有限元模型在预测结构模态参数的过程中不存在模型误差，那么条件概率密度函数可以通过一个Dirac-Delta方程得到：

这里，

其中，和分别表示固有频率和振型的理论解，它们可以通过解特征方程得到。

基于以上的推导的公式，当忽略模型误差时，p(θ|D)可以表示为只与有关。为了构建二阶段模型修正公式，接下来后验概率密度函数将通过利用环境激励下多次测试数据信息得到。

公式(4)中的后验概率密度函数公式即为该发明的总体框架公式，包含两个阶段，即第一阶段：贝叶斯模态识别；第二阶段：基于第一阶段得到的多次测试模态参数，进行贝叶斯模型修正。具体实现方法如下：

2.第一阶段-贝叶斯模态识别,如图2所示

2.1数据采集

采集数据时，将加速度或者速度传感器放在结构上，结构的激励可来自于周围的风荷载、交通荷载、环境噪音及结构中人员活动等。在传感器数目少于需要测试的测点数目时，通过多次测试完成。多次测试需要设置参考点，参考点位置以能采集到尽可能多的有效模态为原则。每次测试时间建议大于第一个周期长度的600倍。每次测试数据时间长度尽量保持一致。

2.2多次测试下贝叶斯模态识别目标函数构建

采集到的多次测试数据，分别进行单次测试数据的模态识别完成，模态识别分两部分进行，首先基于贝叶斯方法进行模态参数最优值的识别，然后进行模态参数后验不确定性的计算。将单次模态识别得到的模态参数进行收集，用于后期的模型修正。

贝叶斯模态识别方法基本原理是要识别的模态的快速傅里叶变换数据在对应的频域段内可以很好的近似为一个高斯概率密度函数。通过最大化这个高斯分布函数，从而可以将模态参数得到。该方法简单概述如下：

在第i个测试的加速度数据可以近似的模拟为：

其中是i次测试的理论加速度响应，该响应是通过将要识别的模态参数来构建的。这些模态参数包括固有频率，阻尼比，模态力的功率谱密度，预测误差的功率谱密度以及振型等。在公式(7)中，表示模型误差，N_i表示样本的数目，n_i表示单次测试自由度的数目。测试数据的快速傅里叶变化可以定义为：

这里，i²＝-1；Δt_i表示i次测试的样本时间间隔；k＝1,...,N_qi；N_qi＝int[N_i/2]+1是奈奎斯特频率的频率指标，int[.]表示整数部分。在i次测试中用来模态识别的数据D_i可以表示为

其中是在i次测试的快速傅里叶变换数据{F_ik}在第r个频域段数据的集合。n_B表示选择的频域段的数目。可以完全确定的概率分布的模态参数可以表示为：

其中

分别表示在r个频域段固有频率和阻尼比的集合；是模态力的功率谱密度，其可以在一个频域段内假设为一个常数；是预测误差的功率谱密度，其也可以在一个频域段内假设为一个常数。同时，

其中表示在第i次测试下第r个频域段的第j阶振型。

基于贝叶斯定理，给定第i次测试数据，的后验概率密度函数可以得到：

其中表示先验概率分布。假设先验信息满足均匀分布，先验概率密度函数可以认为是一个常数。因此后验概率密度函数可以认为直接跟似然函数成正比。当N_i足够大及Δt_i足够小时，不同频率的快速傅里叶变换可以证明其是近似独立的，同时他们的实部和虚部被证明满足高斯分布。因此似然函数可以写为：

其中表示负对数似然函数，其可以通过以下公式得到：

这里‘*’表示复数的共轭转置；

是在频率f_k理论时的理论功率谱密度矩阵；是一个单位矩阵；表示在r个频域段的模态正定转换矩阵，其(p,q)单元可以从下式得到：

其中

公式(16)是贝叶斯模态识别的目标函数，后续的最优值可以通过最小化负对数似然函数来实现。模态参数的协方差矩阵可以通过使其等于目标函数(16)的汉森矩阵的逆来实现。

2.3算法实现：

通过MATLAB编程来实现以上贝叶斯模态识别方法，程序收敛后，可以识别模态参数包括固有频率，振型，阻尼比，模态力的功率谱及预测误差的功率谱。若程序不收敛，需要重新选择频率段，进行重复循环。其中固有频率和振型以及其对应的参数的协方差矩阵将用到后期的第二阶段的贝叶斯模型修正。

总之，通过第一阶段的步骤可完成以下工作：

对于每一次测试，i＝1,...,n_s，对于每一个频域段，通过贝叶斯模态识别优化计算模态参数最优值：和其对应的后验协方差矩阵

3第二阶段：贝叶斯模型修正，具体流程如图3所示

3.1构建目标结构的有限元模型

利用计算软件MATLAB或者有限元分析软件ANASYS等，建立目标结构的有限元模型，在后续的模型修正过程中，方便直接进行调用。

3.2输入模态参数及输出模型参数的选择。

选定第一阶段获得的需要输入的多次测试得到的多组模态参数，尽可能将识别得到的所有模态用到，以提供最多的有效信息。同时根据有限元模型，选定需要修正的模型参数，模型参数数目需要根据测点信息及输入模态信息相对应，避免输出过多参数，导致后期优化过程中不收敛。

3.3构建目标函数并优化(贝叶斯模型修正后验概率密度函数)

在这个部分，定义选择矩阵，其可以将全局振型和在单次测试下得到的振型关联起来，从而基于多次测试数据构建模型参数的后验概率密度函数。这里首先介绍如何得到后验概率密度函数，然后介绍如何将负对数似然函数进行重构，从而方便进行优化。

3.3.1选择矩阵

全局振型Φ^(r)可以通过定义一个选择矩阵L_i来将其与i次测试时得到的振型关联。这个矩阵中，当自由度s在第r频道被测到，那么(r,s)对应的数值就等于1，其他值等于0.第i次测试的振型可以从以下公式得到：

为了方便，这个方法中，假设第i次测试的振型向量正则化为1。

3.3.2多次测试下模型参数的后验概率密度函数

让α＝{α_i,i＝1,...,n_s}表示所有测试下的模态参数。基于贝叶斯理论，给定所有测试的数据，α的后验概率密度函数可以通过下式得到：

给定α，假设在多次测试下数据在统计上是独立的，因此

这里应该注意到p(D_i|α)与其他测试时的参数无关，因此

p(D_i|α)＝p(D_i|α_i) (22)

从而

其中，

这里由i次测试得到的固有频率和部分振型组成

其中f_i和Φ_i分别表示在i次测试下所有选择的频率段内所有频率和阻尼比。参数υ_i由i次测试下剩下的其他模态参数组成，

υ_i＝{ζ_i,S_i,S_ei} (26)

其中ζ_i,S_i和S_ei分别表示在i次测试下所有选择的频率段内阻尼比，模态力的功率谱密度和预测误差的功率谱密度。

因此，基于贝叶斯定理，公式(23)可以由下式得到：

因为p(D)和p(D_i)可以认为是常数，所以公式(27)可以重新写为：

假设先验信息为均匀分布，可以得到：

因此，在第i次测试时，第一阶段的模态参数的后验概率密度函数p₀(α_i|D_i)可以从下式得到：

其中可以通过公式(16)得到。

假设每个是在全局范围内可识别的，在i次测试下，在公式(30)中的每个的后验概率密度函数可以很好的近似为一个高斯分布，其均值为最大可能值协方差矩阵为识别的模态参数协方差矩阵其分布可以写为：

在i次测试下，的边缘后验概率分布函数仍然是一个高斯分布，因此

其中和分别为的最优值和协方差矩阵，其可以从对应的和中的某一部分提取。

考虑多次测试下，基于公式(29)，我们可以得到：

其中

同时

假设固有频率和振型可以完全由结构模型参数决定，将(5)和(33)代入(4)，后验概率密度函数p(θ|D)可以表示为：

其中

这里表示在第r个频域段由有限元模型计算得到的固有频率，其中表示由有限元模型计算得到的对应测试自由度的振型。

3.3.3负对数似然函数的重构

由于振型存在着范数约束，在公式(37)中计算时会出现数值计算问题，因此非常有必要在计算过程中通过计算矩阵的特征基来克服这个问题。经过重构，公式(37)可以写为：

这里和分别是在i次测试下的r频域段中的汉森矩阵的特征值和特征向量。通过重构，不需要计算任何矩阵的逆。

基于目标函数(39)，通过输入模态参数及其协方差矩阵编写程序，优化使其达到最小值。若程序收敛，可以得到模型修正参数θ的最优值；若程序不收敛，那么需要回到开始的地方，调整有限元模型及选择模型修正参数进行循环计算，直至程序收敛。

3.4结构模型参数不确定性计算

在二次泰勒近似的情况下，当θ达到最优值时，后验协方差矩阵可以通过计算负对数似然函数的汉森矩阵的逆来得到，该汉森矩阵可以通过有限差分法来得到。从而我们可以实现评估得到的模型参数的不确定性。

本发明为本领域贡献的关键技术：

1)贝叶斯模型修正目标函数的推导过程。

2)基于多次测试数据用有限传感器覆盖尽可能多测点并进行模态识别和模型修正的理念；

3)基于模态参数最优值和其不确定性的多次测试贝叶斯模型修正目标函数。

Claims (6)

1.一种基于环境激励数据的多次测试下贝叶斯模型修正方法，其特征在于，分两个阶段，

第一阶段是对多次测试下采集的环境激励下结构加速度数据进行分析，得到每次测试测得的结构的固有频率和振型，并计算这些模态参数的不确定性，用协方差矩阵来表示；

第二个阶段是基于多次测试得到的结构模态参数及其协方差矩阵，基于贝叶斯理论构建目标函数，通过对目标函数的优化，得到需要修正的有限元模型的模型参数的最优值。

2.如权利要求1所述的基于环境激励数据的多次测试下贝叶斯模型修正方法，其特征在于，总体构建方法如下：

从结构动力学的基本原理出发，考虑一个线弹性的结构满足以下的动力方程：

$M\stackrel{\·\·}{x}\left(t\right)+C\stackrel{\·}{x}\left(t\right)+Kx\left(t\right)=W\left(t\right)---\left(1\right)$

这里M,C,K分别表示结构的质量，阻尼和刚度矩阵，W是外力向量。假设该结构满足经典阻尼，结构的加速度可以从下式得到：

$\stackrel{\·\·}{x}\left(t\right)=\underset{i}{\Σ}{u}_i{\stackrel{\·\·}{\mathrm{\η}}}_i\left(t\right)---\left(2\right)$

这里，u_i是第i阶全振型向量，是第i阶模态的模态加速度响应。刚度质量的关系可以通过以下特征方程得到：

${\mathrm{Ku}}_i={\mathrm{\ω}}_i^2{\mathrm{Mu}}_i---\left(3\right)$

这里ω_i表示结构的第i阶固有频率。让θ表示与结构的刚度矩阵K和质量矩阵M相关的结构参数。已知刚度和质量矩阵，结构的固有频率和全振型理论上通过特征值分解得到。因此，构建一个理论模型来进行模型修正从而确定θ。

让D＝{D_i:i＝1,...,n_s}表示多次测试得到的用来进行结构模态识别的数据，其中D_i表示第i次测试得到的数据。基于两阶段的模型修正公式和多次测试数据，得到结构参数θ的后验分布：

其中，p(θ)表示结构参数的先验分布；由固有频率和部分振型组成。由于可以通过有限元模型得到，其提供了在模型修正过程中第一阶段和第二阶段相互关联的以下信息。条件概率密度函数表示给定结构模型参数的条件下，结构模态参数的先验概率分布；表示综合了多次测试数据的的边缘后验分布，这里在第一阶段的先验分布被认为是均匀分布。假设有限元模型在预测结构模态参数的过程中不存在模型误差，那么条件概率密度函数可以通过一个Dirac-Delta方程得到：

这里，

其中，和分别表示固有频率和振型的理论解，它们可以通过解特征方程得到。

3.如权利要求2所述的基于环境激励数据的多次测试下贝叶斯模型修正方法，其特征在于，基于以上的推导，当忽略模型误差时，p(θ|D)可以表示为只与有关。所述后验概率密度函数将通过利用环境激励下多次测试数据信息得到，为总体框架公式。

4.如权利要求1或者2所述的基于环境激励数据的多次测试下贝叶斯模型修正方法，其特征在于，第一阶段-贝叶斯模态识别,具体实现方法如下：

2.1数据采集

采集数据时，将加速度或者速度传感器放在结构上，结构的激励可来自于周围的风荷载、交通荷载、环境噪音及结构中人员活动。在传感器数目少于需要测试的测点数目时，通过多次测试完成。

2.2多次测试下贝叶斯模态识别目标函数构建

采集到的多次测试数据，分别进行单次测试数据的模态识别完成，模态识别分两部分进行，首先基于贝叶斯方法进行模态参数最优值的识别，然后进行模态参数后验不确定性的计算。将单次模态识别得到的模态参数进行收集，用于后期的模型修正。

贝叶斯模态识别方法基本原理是要识别的模态的快速傅里叶变换数据在对应的频域段内可以很好的近似为一个高斯概率密度函数。通过最大化这个高斯分布函数，从而可以将模态参数得到。该方法如下：

在第i个测试的加速度数据可以近似的模拟为：

${\stackrel{\·\·}{y}}_{ij}={\stackrel{\·\·}{x}}_{ij}+e_{ij}---\left(7\right)$

其中是i次测试的理论加速度响应，该响应是通过将要识别的模态参数来构建的。这些模态参数包括固有频率，阻尼比，模态力的功率谱密度，预测误差的功率谱密度以及振型等。在公式(7)中，表示模型误差，N_i表示样本的数目，n_i表示单次测试自由度的数目。测试数据的快速傅里叶变化可以定义为：

$F_{ik}=\sqrt{\frac{2{\mathrm{\Δt}}_i}{N_i}}\underset{j=1}{\overset{N_1}{\Σ}}{\stackrel{\·\·}{y}}_k\exp \[-2\mathrm{\π}i\frac{(k-1)(j-1)}{N_i}\]---\left(8\right)$

这里，i²＝-1；Δt_i表示i次测试的样本时间间隔；k＝1,...,N_qi；N_qi＝int[N_i/2]+1是奈奎斯特频率的频率指标，int[.]表示整数部分。在i次测试中用来模态识别的数据D_i可以表示为

$D_i=\{D_i^{\left(r\right)}:r=1,\mathrm{...},n_B\}---\left(9\right)$

其中是在i次测试的快速傅里叶变换数据{F_ik}在第r个频域段数据的集合。n_B表示选择的频域段的数目。可以完全确定的概率分布的模态参数可以表示为：

${\mathrm{\α}}_i^{\left(r\right)}=\[f_i^{\left(r\right)},{\mathrm{\Φ}}_i^{\left(r\right)},{\mathrm{\ζ}}_i^{\left(r\right)},S_i^{\left(r\right)},S_{ei}^{\left(r\right)}\]---\left(10\right)$

其中

$f_i^{\left(r\right)}=\[f_{i1}^{\left(r\right)},\mathrm{...},f_{{\mathrm{im}}_r}^{\left(r\right)}\]\∈R^{m_r}---\left(11\right)$

${\mathrm{\ζ}}_i^{\left(r\right)}=\[{\mathrm{\ζ}}_{i1}^{\left(r\right)},\mathrm{...},{\mathrm{\ζ}}_{{\mathrm{im}}_r}^{\left(r\right)}\]\∈R^{m_r}---\left(12\right)$

分别表示在r个频域段固有频率和阻尼比的集合；是模态力的功率谱密度，其可以在一个频域段内假设为一个常数；是预测误差的功率谱密度，其也可以在一个频域段内假设为一个常数。同时，

${\mathrm{\Φ}}_i^{\left(r\right)}=\[{\mathrm{\Φ}}_{i1}^{\left(r\right)},\mathrm{...},{\mathrm{\Φ}}_{{\mathrm{im}}_r}^{\left(r\right)}\]\∈R^{n_i\×m_r}---\left(13\right)$

其中表示在第i次测试下第r个频域段的第j阶振型。

基于贝叶斯定理，给定第i次测试数据，的后验概率密度函数可以得到：

$p\left({\mathrm{\α}}_i^{\left(r\right)}\right|D_i^{\left(r\right)})\∝p\left({\mathrm{\α}}_i^{\left(r\right)}\right)p\left(D_i^{\left(r\right)}\right|{\mathrm{\α}}_i^{\left(r\right)})---\left(14\right)$

其中表示先验概率分布。假设先验信息满足均匀分布，先验概率密度函数可以认为是一个常数。因此后验概率密度函数可以认为直接跟似然函数成正比。当N_i足够大及Δt_i足够小时，不同频率的快速傅里叶变换可以证明其是近似独立的，同时他们的实部和虚部被证明满足高斯分布。因此似然函数可以写为：

$p\left(D_i^{\left(r\right)}\right|{\mathrm{\α}}_i^{\left(r\right)})\∝\exp \[-L_i^{\left(r\right)}\left({\mathrm{\α}}_i^{\left(r\right)}\right)\]---\left(15\right)$

其中表示负对数似然函数，其可以通过以下公式得到：

$L_i^{\left(r\right)}\left({\mathrm{\α}}_i^{\left(r\right)}\right)=\underset{k\∈I_r}{\Σ}\[ln\left|E_{ik}^{\left(r\right)}\right|+F_{ik}^{*}{E}_{ik}^{\left(r\right)-1}{F}_{ik}\]---\left(16\right)$

这里‘*’表示复数的共轭转置；

$E_{ik}^{\left(r\right)}={\mathrm{\Φ}}_i^{\left(r\right)}{H}_{ik}^{\left(r\right)}{\mathrm{\Φ}}_i^{\left(r\right)T}+S_{ei}^{\left(r\right)}{I}_{n_i}\∈C^{n_i\×n_i}---\left(17\right)$

是在频率f_k理论时的理论功率谱密度矩阵；是一个单位矩阵；表示在r个频域段的模态正定转换矩阵，其(p,q)单元可以从下式得到：

$H_{ik}^{\left(r\right)}(p,q)=S_i^{\left(r\right)}(p,q){\[({\mathrm{\β}}_{ipk}^{\left(r\right)2}-1)+i\left(2{\mathrm{\ζ}}_{ip}^{\left(r\right)}{\mathrm{\β}}_{ipk}^{\left(r\right)}\right)\]}^{-1}{\[({\mathrm{\β}}_{iqk}^{\left(r\right)2}-1)-i\left(2{\mathrm{\ζ}}_{iq}^{\left(r\right)}{\mathrm{\β}}_{iqk}^{\left(r\right)}\right)\]}^{-1}---\left(18\right)$

其中

2.3算法实现：

通过MATLAB编程来实现以上贝叶斯模态识别方法，程序收敛后，可以识别模态参数包括固有频率，振型，阻尼比，模态力的功率谱及预测误差的功率谱。若程序不收敛，需要重新选择频率段，进行重复循环。其中固有频率和振型以及其对应的参数的协方差矩阵将用到后期的第二阶段的贝叶斯模型修正。

5.如权利要求4所述的基于环境激励数据的多次测试下贝叶斯模型修正方法，其特征在于，通过第一阶段的步骤完成以下工作：

对于每一次测试，i＝1,...,n_s，对于每一个频域段，通过贝叶斯模态识别优化计算模态参数最优值：和其对应的后验协方差矩阵

6.如权利要求1或者2所述的基于环境激励数据的多次测试下贝叶斯模型修正方法，其特征在于，第二阶段：贝叶斯模型修正，具体流程：

3.1构建目标结构的有限元模型

建立目标结构的有限元模型，在后续的模型修正过程中，直接进行调用。

3.2输入模态参数及输出模型参数的选择。

选定第一阶段获得的需要输入的多次测试得到的多组模态参数，同时根据有限元模型，选定需要修正的模型参数，模型参数数目需要根据测点信息及输入模态信息相对应。

3.3构建目标函数并优化(贝叶斯模型修正后验概率密度函数)

定义选择矩阵，将全局振型和在单次测试下得到的振型关联起来，从而基于多次测试数据构建模型参数的后验概率密度函数。

3.3.1选择矩阵

全局振型Φ^(r)可以通过定义一个选择矩阵L_i来将其与i次测试时得到的振型关联。这个矩阵中，当自由度s在第r频道被测到，那么(r,s)对应的数值就等于1，其他值等于0.第i次测试的振型可以从以下公式得到：

${\mathrm{\Φ}}_i^{\left(r\right)}=L_i{\mathrm{\Φ}}^{\left(r\right)}---\left(19\right)$

假设第i次测试的振型向量正则化为1。

3.3.2多次测试下模型参数的后验概率密度函数

让α＝{α_i,i＝1,...,n_s}表示所有测试下的模态参数。基于贝叶斯理论，给定所有测试的数据，α的后验概率密度函数可以通过下式得到：

$p\left(\mathrm{\α}\right|D)=\frac{p\left(\mathrm{\α}\right)}{p\left(D\right)}p\left(D\right|\mathrm{\α})=\frac{p\left(\mathrm{\α}\right)}{p\left(D\right)}p\left(\right\{D_1,D_2,\mathrm{...},D_{n_s}\left\}\right|\mathrm{\α})---\left(20\right)$

给定α，假设在多次测试下数据在统计上是独立的，因此

$p\left(\right\{D_1,D_2,\mathrm{...}{D}_{n_s}\left\}\right|\mathrm{\α})=p\left(D_1\right|\mathrm{\α})p\left(D_2\right|\mathrm{\α}),\mathrm{...},p\left(D_{n_s}\right|\mathrm{\α})---\left(21\right)$

这里应该注意到p(D_i|α)与其他测试时的参数无关，因此

p(D_i|α)＝p(D_i|α_i) (22)

从而

$p\left(\mathrm{\α}\right|D)=\frac{p\left(\mathrm{\α}\right)}{p\left(D\right)}\underset{i=1}{\overset{n_s}{\Π}}p\left(D_i\right|{\mathrm{\α}}_i)---\left(23\right)$

其中，

这里由i次测试得到的固有频率和部分振型组成

其中f_i和Φ_i分别表示在i次测试下所有选择的频率段内所有频率和阻尼比。参数υ_i由i次测试下剩下的其他模态参数组成，

υ_i＝{ζ_i,S_i,S_ei} (26)

其中ζ_i,S_i和S_ei分别表示在i次测试下所有选择的频率段内阻尼比，模态力的功率谱密度和预测误差的功率谱密度。

因此，基于贝叶斯定理，公式(23)可以由下式得到：

$p\left(\mathrm{\α}\right|D)=\frac{p\left(\mathrm{\α}\right)}{p\left(D\right)}\[\underset{i=1}{\overset{n_s}{\Π}}\frac{p\left(D_i\right)}{p\left({\mathrm{\α}}_i\right)}\]\[\underset{i=1}{\overset{n_s}{\Π}}p\left({\mathrm{\α}}_i\right|D_i)\]---\left(27\right)$

因为p(D)和p(D_i)可以认为是常数，所以公式(27)可以重新写为：

$p\left(\mathrm{\α}\right|D)\∝\frac{p\left(\mathrm{\α}\right)}{\underset{i=1}{\overset{n_s}{\Π}}p\left({\mathrm{\α}}_i\right)}\underset{i=1}{\overset{n_s}{\Π}}p\left({\mathrm{\α}}_i\right|D_i)---\left(28\right)$

假设先验信息为均匀分布，可以得到：

$p\left(\mathrm{\α}\right|D)\∝\underset{i=1}{\overset{n_s}{\Π}}p\left({\mathrm{\α}}_i\right|D_i)---\left(29\right)$

因此，在第i次测试时，第一阶段的模态参数的后验概率密度函数p₀(α_i|D_i)可以从下式得到：

$p_0\left({\mathrm{\α}}_r\right|D_i)\∝\underset{r=1}{\overset{n_B}{\Π}}{p}_0\left({\mathrm{\α}}_i^{\left(r\right)}\right|D_i^{\left(r\right)})\∝\exp \[-\underset{r=1}{\overset{n_B}{\Σ}}{L}_i^{\left(r\right)}\left({\mathrm{\α}}_i^{\left(r\right)}\right)\]---\left(30\right)$

其中可以通过公式(16)得到。

假设每个是在全局范围内可识别的，在i次测试下，在公式(30)中的每个的后验概率密度函数可以很好的近似为一个高斯分布，其均值为最大可能值协方差矩阵为识别的模态参数协方差矩阵其分布可以写为：

$p_0\left({\mathrm{\α}}_i^{\left(r\right)}\right|D_i^{\left(r\right)})=\mathrm{\φ}({\mathrm{\α}}_i^{\left(r\right)};{\widehat{\mathrm{\α}}}_i^{\left(r\right)},{\hat{C}}_i^{\left(r\right)})---\left(31\right)$

在i次测试下，的边缘后验概率分布函数仍然是一个高斯分布，因此

其中和分别为的最优值和协方差矩阵，其可以从对应的和中的某一部分提取。

考虑多次测试下，基于公式(29)，可以得到：

其中

同时

假设固有频率和振型可以完全由结构模型参数决定，将(5)和(33)代入(4)，后验概率密度函数p(θ|D)可以表示为：

其中

这里表示在第r个频域段由有限元模型计算得到的固有频率，其中表示由有限元模型计算得到的对应测试自由度的振型。

3.3.3负对数似然函数的重构

由于振型存在着范数约束，在公式(37)中计算时会出现数值计算问题，为此，在计算过程中通过计算矩阵的特征基。经过重构，公式(37)可以写为：

这里和分别是在i次测试下的r频域段中的汉森矩阵的特征值和特征向量。通过重构，不需要计算任何矩阵的逆。

基于目标函数(39)，通过输入模态参数及其协方差矩阵编写程序，优化使其达到最小值。若程序收敛，可以得到模型修正参数θ的最优值；若程序不收敛，那么需要回到开始的地方，调整有限元模型及选择模型修正参数进行循环计算，直至程序收敛。

3.4结构模型参数不确定性计算

在二次泰勒近似的情况下，当θ达到最优值时，后验协方差矩阵可以通过计算负对数似然函数的汉森矩阵的逆来得到，该汉森矩阵可以通过有限差分法来得到，实现评估得到的模型参数的不确定性。