基于Wasserstein距离的近似贝叶斯模型修正

Approximate Bayesian model updating based on Wasserstein distance

摘要

针对传统贝叶斯方法在计算似然函数时归一化因子的高维积分难以求解的问题， 提出一种基于Wasserstein距离的近似贝叶斯有限元模型修正方法. 首先， 构建径向基函数模型替代有限元模型， 以减小计算量；其次， 选取Wasserstein距离作为相似性度量指标， 衡量试验观测数据与有限元模型模拟数据之间的整体分布差异， 以克服传统度量方式对局部误差的敏感性；最后， 采用近似贝叶斯-序贯蒙特卡罗采样方法进行似然求解， 经过少量迭代快速收敛， 从而高效识别参数的后验分布. 通过三自由度弹簧、简支梁和钢桁架数值算例以及 3 kW小型风力机叶片试验对所提方法的有效性进行验证.

Abstract

To address the issue of the high-dimensional integral of the normalization factor in the calculation of the likelihood function in traditional Bayesian methods， an approximate Bayesian finite element model updating method based on Wasserstein distance is proposed. Firstly， a radial basis function model is constructed to replace the finite element model to reduce the computational load. Secondly， the Wasserstein distance is selected as the similarity measurement index to measure the overall distribution difference between the experimental observation data and the finite element model simulation data， overcoming the sensitivity of the traditional measurement method to local errors. Finally， the approximate Bayesian-sequential Monte Carlo sampling method is adopted to solve the likelihood， which converges with only a few iterations， thereby efficiently identifying the posterior distribution of the parameters. The validity of proposed method is verified through numerical examples of a three-degree-of-freedom spring， a simply supported beam and a steel truss， as well as a 3 kW small wind blade test.

Graphical abstract

关键词

模型修正 / Wasserstein距离 / 贝叶斯方法 / 蒙特卡罗方法 / 径向基函数网络

Key words

model updating / Wasserstein distance / Bayesian methods / Monte Carlo methods / radial basis function networks

传统确定性模型修正方法已较为成熟， 在基于单次试验的确定性优化问题中得到广泛应用^［1］.然而， 在实际工程应用中，结构的材料参数、几何尺寸、简化建模以及传感器的测试数据带来的不确定性因素普遍存在^［2-5］， 导致传统的确定性模型修正方法难以满足工程精度要求. 如何准确识别参数的不确定性， 并对其进行量化表征， 已成为不确定性模型修正领域研究的关键焦点.

Beck等^［6］最早提出基于贝叶斯理论的不确定性模型修正方法， 突破了传统确定性方法难以表征不确定性的局限. 传统贝叶斯方法应用于复杂模型时难以获得后验分布的解析表达^［7］. Pritchard等^［8］提出近似贝叶斯计算（approximate Bayesian computation， ABC）， 通过计算观测与模拟数据之间的统计距离间接推断后验分布， 此方法的关键在于选择合适的统计距离^［9］. 欧氏距离因计算简单被广泛采用， Khodaparast等^［10］结合欧氏距离与Frobenius范数来同时修正均值和方差. 对于高斯分布的数据， 巴氏距离^［11］通过正态假设实现均值与方差的联合修正， 能更好地衡量数据分布的相似性. 张亚峰等^［12］以最小化巴氏距离为目标， 提出了参数均值和标准差的两步和同步求解. 然而， 上述传统距离在处理低重叠或复杂分布时均存在局限性. Wasserstein距离凭借其几何直观性和对分布支撑集偏移的鲁棒性， 在模型修正中展现出优势^［13-14］. 该距离基于最优传输理论， 通过衡量最小传输代价量化分布差异^［15］. Tao等^［16］将其引入变分贝叶斯框架中， 实现高效的参数后验识别， 优于基于马尔可夫链的随机采样方法.

此外， 在近似贝叶斯方法中， 频繁调用有限元模型增加了计算负担. 代理模型通过学习输入与响应之间的映射， 构建低成本替代模型， 有效降低对原始模型的依赖， 提升了计算效率. 常见的代理模型方法包括径向基函数^［17］（radial basis function， RBF）模型（简称径向基模型）、Kriging模型^［18］和神经网络模型^［19］等. 径向基模型因其结构简单、计算效率高且具备以任意精度逼近连续函数的优势， 在土木、航天和机械等领域得到了广泛应用.

基于上述理论背景， 提出一种基于Wasserstein距离的近似贝叶斯模型修正方法. 首先， 将待修正参数作为输入，通过粒子群优化算法（particle swarm optimization， PSO）优化径向基模型的相关系数， 构建高保真的径向基模型，以替代高计算成本的有限元模型进行迭代计算，减少计算量；其次， 利用Wasserstein距离构造近似似然函数；最后， 采用近似贝叶斯-序贯蒙特卡罗方法进行极大似然估计， 从而得到参数的后验分布. 通过三自由度弹簧-质量系统、简支梁和钢桁架数值算例以及3 kW小型风力机叶片试验对所提方法的有效性进行了验证.

1 径向基模型的构建

1.1 径向基模型原理

径向基模型是一种基于局部加权回归的代理模型， 其核心原理是径向基函数插值法^［20］. 径向基模型按照信息流的方向分为输入层、隐藏层和输出层. 输入层接收原始数据， 隐藏层通过激活函数对输入数据进行处理， 以捕捉局部特征， 输出层则将隐藏层的输出进行线性组合. 隐藏层中的激活函数为高斯函数， 假设为{\alpha }_t(\mathit{X}),t=\mathrm{1,2},\dots ,k， t个隐藏神经元的输出值可以表示为：

{\alpha }_t\left(\mathit{X}\right)=\mathrm{E}\mathrm{x}\mathrm{p}\left(-‖\mathit{X}-{\mathit{\mu }}_t‖/{\sigma }_t^{\mathrm{2}}\right)（1）

式中：\mathit{X}=(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_m{)}^{\mathrm{T}}表示输入向量；{\mathit{\mu }}_t表示第t个隐藏神经元的中心向量；‖\mathit{X}-{\mathit{\mu }}_t‖表示\mathit{X}和{\mathit{\mu }}_t的统计学距离；{\sigma }_t表示第t个隐藏神经元的宽度.

1.2 径向基模型构建与验证

构建径向基模型步骤：首先，基于拉丁超立方抽样方法（Latin hypercube sampling， LHS）从待修正参数的\pm \mathrm{20}\mathrm{\%}区间内抽取样本点；其次，将所抽取的样本划分为互不重叠的训练集和测试集，基于训练集样本，构建初始径向基模型，并将测试集样本代入该模型，依据适应度函数（即模型的验证误差）评估其性能；然后，引入PSO算法优化高斯函数宽度{\sigma }_t， PSO模拟粒子的群体行为，粒子的状态通过迭代更新，依据适应度进行优化，通过自适应权重策略在参数空间内有效地探索，并最终收敛至最优{\sigma }_t值；最后，基于最优{\sigma }_t值和训练集样本，构造最终的径向基模型，通过将数据划分为多个测试子集进行交叉验证，进一步评估模型的性能， 增强模型的泛化能力.

在上述模型验证阶段， 采用均方根误差（root-mean-square error， RMSE）来评价模型的适应度， 其表达式为：

\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt[]{\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}({\mathit{y}}_{\mathrm{o}\mathrm{b}\mathrm{s},i}-{\mathit{y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d},i}{)}^{\mathrm{2}}}（2）

式中：{\mathit{y}}_{\mathrm{o}\mathrm{b}\mathrm{s},i}表示第i个真实观测值；{\mathit{y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d},i}表示第i个预测值；n表示样本总数.

2 近似贝叶斯-序贯蒙特卡罗原理

2.1 近似贝叶斯

基于贝叶斯的模型修正方法将先验知识与观测数据相融合， 通过贝叶斯理论推断后验概率密度函数. 待修正参数的后验分布可用贝叶斯公式表达， 具体为：

p(\mathit{\theta }|{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}})=\frac{L(\mathit{\theta })p(\mathit{\theta })}{\int L(\mathit{\theta })p(\mathit{\theta })\mathrm{d}\mathit{\theta }}（3）

式中：{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}}表示包含通过试验获得的n个观测值的数据集；\mathit{\theta }表示待修正模型的参数；L(\mathit{\theta })表示{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}}所对应的似然函数；p\left(\mathit{\theta }\right)表示模型参数的先验分布.

在处理复杂模型时， 由于参数分布存在不确定性， 很难确立式（3）中似然函数的解析表达式. 针对这一问题， 采用近似贝叶斯计算^［21］， 该方法通过构建近似似然函数来减少计算量. 基于统计学距离的近似似然函数为：

{\tilde{L}}_{Y_{\mathrm{o}\mathrm{b}\mathrm{s}}}\left(\mathit{\theta }\right)\propto \frac{\mathrm{1}}{\epsilon }I\left(\frac{d\left({\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}},{\mathit{Y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}\right)}{\epsilon }\right)（4）

式中：I\left(\cdot \right)表示返回0或1的指示函数；{\mathit{Y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}表示根据参数\mathit{\theta }生成的模拟数据；d\left({\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}},{\mathit{Y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}\right)表示试验观测数据与模拟数据之间的距离；\epsilon表示预设的容差阈值.

近似贝叶斯方法通过计算模拟数据间的距离， 将复杂的似然计算简化为接受-拒绝机制：若新样本满足预设容差， 则接受为近似后验， 否则舍弃. 通过大量抽样， 可获得参数的近似后验分布：

\left\{\begin{array}{l}p(\mathit{\theta }|{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}})\approx p_{\epsilon }(\mathit{\theta }|{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}})\\ (d({\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}},{\mathit{Y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})\le \epsilon ):{\mathit{\theta }}^{\left(t\right)}={\mathit{\theta }}^{\mathrm{*}}\end{array}\right.（5）

式中：{\mathit{\theta }}_{}^{(t)}表示满足容差并被接受的样本， t=\mathrm{1,2},\dots ,T， T表示被接受样本的总数；{\mathit{\theta }}^{\mathrm{*}}表示总体随机抽样中的一个样本.

2.2 序贯蒙特卡罗

在近似贝叶斯计算中， 归一化因子的高维积分难以解析求解， 常借助马尔可夫链蒙特卡罗（Markov chain Monte Carlo，MCMC）求解其数值近似. 但传统MCMC采样效率低， 易陷入“采样停滞”^［22］. 序贯蒙特卡罗（sequential Monte Carlo， SMC）方法在此基础上进一步发展， 基于粒子滤波的思想， 引入顺序重要性抽样（sequential importance sampling， SIS）， 通过构建一系列中间分布并赋予样本重要性权值\omega \left(\mathit{\theta }\right)， 逐步逼近真实后验分布^［23］. 最终可获得难以直接计算的后验：

p_{{\epsilon }_p}(\mathit{\theta }|{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}})\approx \frac{\mathrm{1}}{T}{\displaystyle \sum _{i=\mathrm{1}}^T}{\delta }_{{\theta }_i}(\mathit{\theta })（6）

w(\mathit{\theta })=\frac{p(\mathit{\theta }|{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}})}{p_{{\epsilon }_p}(\mathit{\theta }|{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}})}（7）

p_{\epsilon }(\mathit{\theta }|{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}})\approx \frac{\mathrm{1}}{T}{\displaystyle \sum _{i=\mathrm{1}}^T}w({\mathit{\theta }}_i){\delta }_{{\theta }_i}(\mathit{\theta })（8）

式中：p_{{\epsilon }_p}\left(\mathit{\theta }|{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}}\right)表示中间分布；{\epsilon }_p表示对应总体的容差；{\delta }_{{\theta }_i}(\mathit{\theta })表示狄拉克函数， 用于分析ABC的收敛性.

2.3 近似贝叶斯-序贯蒙特卡罗

近似贝叶斯-序贯蒙特卡罗（approximate Bayesian computation-sequential Monte Carlo， ABC-SMC）方法基于顺序重要性抽样的思想， 通过近似贝叶斯计算生成一系列中间分布：

p_{{\epsilon }_p}(\mathit{\theta }|{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}})=\frac{p(\mathit{\theta })}{T}{\displaystyle \sum _{r=\mathrm{1}}^{T^{\mathrm{\text{'}}}}}I\left(d({\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}},{\mathit{Y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}})\right)（9）

式中：p_{{\epsilon }_p}(\mathit{\theta }|{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}})表示参数在第p代总体中的分布；p\left(\mathit{\theta }\right)表示参数的先验分布；T^{\mathrm{\text{'}}}表示样本数量.

新的中间分布是由粒子扰动分布k\left(·\right)产生的， 它的作用是扰动上一代种群中的一个粒子， 从而产生一个新粒子. 在这一阶段， 第p代种群中粒子的采样重要性权重由式（10）给出：

w_p^{(\mathrm{0})}=\frac{p({\mathit{\theta }}_p^{(\mathrm{0})})}{{\displaystyle \sum _{j=\mathrm{1}}^r}{w}_{p-\mathrm{1}}^{j}k({\mathit{\theta }}_p^{(\mathrm{0})}|{\mathit{\theta }}_{p-\mathrm{1}}^{(\mathrm{0})})}（10）

式中：w_{p-\mathrm{1}}^j表示p-\mathrm{1}代粒子的权重；r表示第p-\mathrm{1}代粒子的总数；j是p-\mathrm{1}代粒子的索引， 表示第p-\mathrm{1}代种群中的第j个粒子；k({\mathit{\theta }}_p^{(\mathrm{0})}|{\mathit{\theta }}_{p-\mathrm{1}}^{(\mathrm{0})})是核函数， 表示从p-\mathrm{1}代粒子{\mathit{\theta }}_{p-\mathrm{1}}^j生成新粒子{\mathit{\theta }}_p^{(\mathrm{0})}的概率密度.

在ABC-SMC中， 微扰核的选择对参数空间的探索效果和计算效率有重要影响. 为此， 本文采用均匀分布摄动核， 即以0为中心， 确保新的采样点不会偏离上一代粒子太远， 同时又能覆盖一定的探索空间. 摄动范围为上一代粒子的参数变化幅度：

\left\{\begin{array}{l}{k}_p({\mathit{\theta }}_p^{(i)}|{\mathit{\theta }}_{p-\mathrm{1}}^{(i)})~\mathrm{u}\mathrm{n}\mathrm{i}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}(-\frac{\mathrm{1}}{\mathrm{2}}{L}_{p-\mathrm{1}},\frac{\mathrm{1}}{\mathrm{2}}{L}_{p-\mathrm{1}})\\ L_{p-\mathrm{1}}=\mathrm{m}\mathrm{a}{\mathrm{x}}_{\mathrm{1}\le t\le T}\left\{{\mathit{\theta }}_{p-\mathrm{1}}^{(t)}\right\}-\mathrm{m}\mathrm{i}{\mathrm{n}}_{\mathrm{1}\le t\le T}\left\{{\mathit{\theta }}_{p-\mathrm{1}}^{(t)}\right\}\end{array}\right.（11）

式中：L_{p-\mathrm{1}}表示上一代粒子的参数变化幅度.

ABC-SMC方法的实施依赖于一系列递减的容差阈值， 以逐步逼近参数的真实后验分布. 然而， 手动设置容差阈值可能会导致过早收敛或低粒子接受率. 为解决这一问题， 采用自适应容差阈值策略， 在保证算法有效运行的同时避免低粒子接受率. 容差阈值序列在预设区间内按迭代线性递减. 此外， 为避免手动设置容差阈值可能导致的算法收敛问题， 引入了自适应终止机制， 当某次迭代的接受率低于1%时， 采样过程将自动停止. 这种自适应容差递减方案在保持算法稳定性的同时， 有效避免了过度收敛导致的低采样率问题， 从而保证了算法的可靠性和计算效率.

3 基于Wasserstein距离的ABC-SMC模型修正

3.1 Wasserstein距离

在ABC-SMC推断过程中， 距离测度是构造近似似然的关键. Wasserstein距离基于最优传输理论， 提供更精确直观的相似性评估^［24］. 将其嵌入ABC-SMC框架， 可以更准确地评估模拟数据与观测数据之间的相似性， 从而在容差阈值的动态调整过程中提供更可靠的依据.

假设存在两个概率分布\mu和\nu， 以及一个成本函数c\left(X,Y\right)， 传统Wasserstein距离定义为：

W_{\lambda }\left(\mu ,\nu \right)={\left[\underset{\gamma \in \Gamma (\mu ,\nu )}{\mathrm{i}\mathrm{n}\mathrm{f}}{{\int }_{X\times Y}c\left(X,Y\right)}^{\lambda }\mathrm{d}\gamma (X,Y)\right]}^{\mathrm{1}/\lambda }（12）

式中：X和Y是分别对应于概率分布\mu和\nu的样本空间；\mathrm{i}\mathrm{n}\mathrm{f}表示在所有可能的联合分布\gamma \left(X,Y\right)中寻找使整体运输代价最小的方案；c\left(X,Y\right)表示从X运输到Y的代价；\lambda是一个实数， 通常取1或2；\Gamma (\mu ,\nu )表示所有联合分布\gamma (X,Y)的集合， 满足其边缘分布分别为概率分布\mu和\nu， 约束条件为：

{\int }_Y\gamma (X,Y)\mathrm{d}Y=P(X)（13）

{\int }_X\gamma (X,Y)\mathrm{d}X=Q(Y)（14）

当概率分布是一维的情况时， Wasserstein距离可以表示为：

W_{\lambda }(X,Y)=({\int }_{\mathrm{0}}^{\mathrm{1}}{\left|F_X^{-\mathrm{1}}(\alpha )-F_Y^{-\mathrm{1}}(\alpha )\right|}^{\lambda }{\mathrm{d}\alpha )}^{\mathrm{1}/\lambda }（15）

式中：F_X和F_Y表示变量X和Y累积分布函数；F_Y^{-\mathrm{1}}表示从一个均匀随机变量到Y的最优映射. 当实数\lambda =\mathrm{1}时， 式（15）可简化为：

W_{\mathrm{1}}(X,Y)={\int }_{\mathrm{0}}^{\mathrm{1}}\left|F_X^{-\mathrm{1}}\left(t\right)-F_Y^{-\mathrm{1}}\left(t\right)\right|\mathrm{d}t（16）

3.2 模型修正

本文提出的基于Wasserstein距离的ABC-SMC模型修正流程如图1所示.

1）先验抽样：从参数的先验分布p\left(\mathit{\theta }\right)中随机抽取N_{\mathrm{a}}个样本， 记为{\mathit{\theta }}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(N_{\mathrm{a}})}， 并为每个粒子赋予相同的初始权重值{\omega }_i^{\mathrm{0}}.

2）径向基模型的构建：以步骤1）得到的样本作为输入， 构建径向基模型， 并进行精度评估.

3）生成模拟数据：运行径向基模型生成对应的模拟数据， 记为{\mathit{Y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(N_{\mathrm{a}})}.

4）计算Wasserstein距离：计算每个模拟数据{\mathit{\theta }}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(N_a)}生成的模拟样本集{\mathit{Y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(N_{\mathrm{a}})}与试验观测数据{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}}之间的Wasserstein距离d_W=W_{\mathrm{1}}({\mathit{Y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(N_{\mathrm{a}})},{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}}).

5）计算近似似然函数：将步骤4）中计算的Wasserstein距离代入式（4）， 得到{\tilde{L}}_{{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}}}\left({\mathit{\theta }}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(N_{\mathrm{a}})}\right).

6）更新权重与重采样：使用序贯蒙特卡罗抽样算法，结合步骤4）中计算的Wasserstein距离更新粒子权重，当距离值小于当前容差阈值时， 保留粒子权重；否则，对粒子进行重采样， 直到满足终止条件，停止迭代.

7）后验分布评估：基于最终粒子合集的后验分布p\left({\mathit{\theta }}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}^{(N_{\mathrm{a}})}|{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}}\right)进行评估， 得到最终参数修正值.

4 数值分析

4.1 三自由度弹簧-质量系统数值分析

三自由度（three-degree-of-freedom， 3-DOF）弹簧-质量系统是一种典型的动力学模型，其广泛应用于悬挂系统、减震器设计等方面，采用该系统仿真对所提方法进行验证， 如图2所示.

对待修正参数进行不确定性建模， 其中均值\mu和标准差\sigma是给定区间内的未知数， 如表1所示.

为验证Wasserstein距离在模型修正中的效果， 本文对比分析了欧氏距离与巴氏距离两种方法， 在相同条件下进行参数修正， 结果见表2.

结果表明， 3种方法在参数均值估计上性能相近，但Wasserstein距离在标准差估计方面表现更优，显著降低了相对误差，体现出其在量化数据离散性方面的优势.同时，Wasserstein方法得到的95%置信区间更窄，表明其不确定性较小， 预测更集中， 展现出优于欧氏和巴氏距离的稳健性与准确性.

4.2 简支梁数值分析

在工程实际应用中，许多工程结构都可以简化为简支梁形式.为了进一步验证所提方法的有效性， 以钢结构简支梁为对象，进行数值分析. 简支梁结构如图3所示.尺寸为\mathrm{2}\mathrm{000}\mathrm{m}\mathrm{m}\times \mathrm{100}\mathrm{m}\mathrm{m}\times \mathrm{10}\mathrm{m}\mathrm{m}，将该外伸梁划分为20个单元，21个节点，每单元尺寸为\mathrm{100}\mathrm{m}\mathrm{m}\times \mathrm{100}\mathrm{m}\mathrm{m}\times \mathrm{10}\mathrm{m}\mathrm{m}，只考虑每个节点的垂向自由度.在梁两端设置约束方式，结构1节点为固定铰支座，21节点为简支支座.简支梁材料为Q235 钢，材料密度为7 800 kg/m³，弹性模量为210 GPa. 结构频率分别为5.88 Hz、23.53 Hz、52.97 Hz. 初始有限元频率为6.09 Hz、22.79 Hz、51.29 Hz.

4.2.1 简支梁径向基模型的构建

采用拉丁超立方抽样生成2 000个样本， 80%用于训练， 20%用于验证， 按1.2节方法构建径向基模型. 一阶频率预测与真实值对比如图4所示，数据分布表明模型预测精度较高.

径向基模型在不同阶数下的误差如表3所示， 除了RMSE外， 还列出了平均绝对误差（mean absolute error， MAE）、平均偏差（mean bias error， MBE）指标作为评判标准， 从多个角度综合评估模型的预测性能， 以确保模型的可靠性和准确性. 表中结果表明， 径向基模型误差普遍处于10^-4数量级， 说明其在精度上能够较好地逼近有限元模型的计算结果， 可代替有限元模型进行迭代计算.

4.2.2 简支梁模型修正

简支梁的动态特性受多种不确定性因素的影响， 梁的几何尺寸相对稳定. 由于材料性质和制造工艺的差异， 弹性模量E和材料密度d可能会在梁的横截面内呈现出不均匀分布， 将其作为待修正参数. 对待修正参数进行不确定性建模， 如表4所示.

通过从不确定性参数中抽取先验样本， 将这些样本输入径向基模型， 得到1 000个模拟样本{\mathit{Y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}， 计算观测数据{\mathit{Y}}_{\mathrm{o}\mathrm{b}\mathrm{s}}与模拟数据{\mathit{Y}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}之间的Wasserstein距离构建近似似然函数， 采用SMC方法迭代获取参数后验样本. 均值{\mu }_d和{\mu }_E的先验样本与后验样本的结果分别如图5和图6所示. 可以看出， 修正前参数存在较大的不确定性. 修正后的后验分布趋于收敛， 且接近真实样本的分布. 表明Wasserstein距离修正效果良好.

在相同条件下， 比较欧氏距离与巴氏距离的修正结果， 如表5所示.

Wasserstein距离在参数的标准差估计方面表现最优， 修正误差分别为0.14%、2.26%， 显著优于欧氏距离的1.46%和6.1%， 以及巴氏距离的39.99%和17.89%.

修正前后模态频率均值与真实值的对比如图7所示.Wasserstein距离法修正结果最优，频率均值为5.87、23.49、52.81 Hz， 优于欧氏距离法（5.79、23.17、52.16 Hz）和巴氏距离法（5.86、23.43、52.75 Hz）. 整体来看，Wasserstein距离法有效降低了修正误差.

4.3 钢桁架试验数值分析

本节基于工程实践需求，针对实验室钢桁架模型进行数值仿真研究.该下承式钢桁架模型采用Q235钢材，总质量为54 kg.该结构体系由两个跨度2 800 mm的主桁架与连接系统共同构成，沿x方向划分为7个400 mm等分节段.主桁架截面尺寸为270 mm（高）×390 mm（宽），跨高比为10∶1，其中弦杆选用等边热轧角钢， 竖向腹杆采用热轧T型钢， 斜向腹杆则选用矩形钢， 其结构如图8所示.

在ANSYS软件中建立该钢桁架的有限元模型， 结构杆件根据传力特性分为三类：上弦杆、中部杆和下弦杆. 所有构件采用BEAM188单元建模， 梁单元通过直接连接节点连接， 将每个构件离散为单个梁单元，完整模型共包含28个节点与87个单元. 边界条件方面，固定结构底部四个角节点的全部自由度，用以模拟支座约束条件.在完成建模后，设置分析类型为模态分析，选用Block Lanczos算法求解前三阶频率，为后续模型修正提供模态观测基准.结构前三阶频率分别为18.129 Hz、45.493 Hz、57.968 Hz.将初始弹性模量均值赋予各构件， 设定初始弹性模量E_{\mathrm{0}}=\mathrm{190}\mathrm{G}\mathrm{P}\mathrm{a}，初始有限元模型的前三阶固有频率分别为17.924 Hz、42.232 Hz、48.340 Hz.

4.3.1 钢桁架径向基模型的构建

本文采用LHS生成2 000组样本， 通过有限元软件响应面模块完成样本生成， 并基于1.2节的方法构建径向基模型. 模型精度从局部误差和全局统计指标两方面评估.

各测试组一阶预测频率与真实频率对比及对应的RMSE如图9所示. 误差总体处于10^-3数量级.

从整体层面对代理模型精度进行评估， 结果如表6所示. 结果表明该径向基模型可替代有限元模型进行迭代计算.

4.3.2 钢桁架模型修正

鉴于结构损伤通常表现为刚度退化的特征， 本文选取三类杆组的弹性模量作为待修正参数. 对待修正参数进行不确定性建模， 其中均值\mu和标准差\sigma是给定区间内的未知数， 如表7所示.

3种方法下，参数的修正后验直方图如图10、图11和图12所示.不同方法下的参数修正值、与真实值的误差以及95%置信区间的具体对比结果如表8所示.

采用Wasserstein距离时所得的后验分布形状相对更为平滑， 对角线上的直方图呈现出近似高斯分布的特征， 样本集中度较高， 所得到的修正值最接近真实参数，误差最小.而欧氏距离所对应的后验分布则相对分散，部分参数的直方图偏离高斯形态，显示出较强的不确定性，置信区间显著扩大， 修正值误差也随之增加.巴氏距离生成的后验分布较集中，但均值位置相较真实值存在一定偏移.

综上，根据图示的后验分布形态及表中定量指标综合分析，Wasserstein距离在参数修正任务中展现出分布集中性良好、修正值精确度高及不确定性适中等优势.

图13对比了修正前后模态频率均值. Wasserstein距离修正后均值为18.372 Hz、45.470 Hz、53.617 Hz，二阶频率预测误差最小，绝对误差为0.023 Hz，相对误差为0.05%.欧氏距离修正后的频率均值为18.191 Hz、46.178 Hz、54.633 Hz，一阶误差较小，但高阶精度下降；巴氏距离修正后频率均值为17.755 Hz、42.231 Hz、50.528 2 Hz，各阶误差较大，高阶相对误差达13.26%.

5 试验验证

为了进一步评估所提方法的有效性，本文选取3 kW小型风力机叶片作为研究对象.该叶片的总长度为2.4 m，其材料为环氧树脂基体及玻璃纤维增强体.如图14所示，试验中采用悬臂梁式安装方法，将叶片根部利用螺栓孔固定于试验台， 最大限度地复现风机叶片的实际工作条件.

在试验过程中， 采用INV3062-C2（L）型网络分布式DASP采集仪进行数据采集， 并结合DASP-V11信号分析软件进行数据处理和模态参数识别. 设置采样频率为2 000 Hz， 同时在叶片表面布置加速度传感器. 采用MSC-3型弹性力锤对叶片施加垂向激励， 在锤头激励下获取结构的低阶模态.采用频响函数采峰法获取模态频率.

由于缺乏叶片的设计图纸，参考文献［25］中提出的建模策略，采用逆向工程方法构建叶片的有限元模型.采用3DSHandy-7XB Plus手持式三维扫描仪对叶片进行外形扫描，以获取其精确的几何数据. 该设备的扫描精度高达0.02 mm，确保了所建立三维模型的高保真度.由于三维扫描仪获取的数据为激光点云数据，无法直接用于有限元建模，因此，后续采用Design X软件对点云数据进行处理，并转换为实体模型.在此基础上，利用ANSYS软件建立叶片的有限元模型，使用Solid 186实体单元进行网格划分.网格尺寸设定为50 mm.最终，所建立的有限元模型包含 13 561个单元和6 089个节点.但是该模型在参数和边界条件等方面难以完全反映真实结构特性，因此需进行模型修正以提高其预测精度.前三阶模态频率以及对应的初始有限元分析结果如表9所示.

5.1 径向基模型的构建及评估

在初始有限元计算值的±20%范围内，采用LHS生成3 000个样本构造径向基模型. 对训练后的径向基模型进行精度评估， 如图15所示.RMSE值整体保持在较低的水平，即使在预测值波动幅度较大的点也能够很好地捕捉真实值的变化趋势.

径向基模型在不同阶数下的误差如表10所示， 包括MAE、MBE和RMSE. 从表中数据可以看出， MAE 数值较小， 表明径向基模型预测的误差整体较低. MBE绝对值较小， 表明没有显著的系统性偏差. 所构建的径向基模型可以替代有限元模型.

5.2 风力机叶片模型修正

选取对动力学响应的敏感度最高的质量密度{\rho }_{\mathrm{b}}和弹性模量E_{\mathrm{b}}作为待修正参数.为了有效描述待修正参数的不确定性，对待修正参数进行不确定性建模. 待修正参数的不确定性的建模如表11所示， \mu表示均值，\sigma表示标准差，均为给定区间内的未知数.

三组后验分布的直方图如图16~图18所示，每组直方图对应不同的距离度量方法，分别对应Wasserstein距离、欧氏距离和巴氏距离.其中，各子图横坐标分别为材料密度与弹性模量，纵坐标为概率密度函数（PDF， probability density function），用于量化参数后验分布的概率密度水平.综合来看，基于Wasserstein距离修正的后验分布与真实值匹配较好.

为进一步分析不同距离的效果，表12总结了基于不同距离度量， 包括Wasserstein距离、欧氏距离和巴氏距离，计算的修正参数均值及其对应的修正误差.根据表中数据可以看出，Wasserstein距离在所有参数上的误差最小，显示出较高的修正精度和稳定性.欧氏距离的误差略高于Wasserstein距离，但仍然保持在较低水平.巴氏距离的误差显著高于前两者，特别是在{\mu }_{E_{\mathrm{b}}}和{\sigma }_{E_{\mathrm{b}}}上，显示出较大的不确定性.

修正前后的模态频率均值以及其与真实模态频率的对比如图19所示.基于Wasserstein距离的方法修正后的频率均值为8.260、12.527、30.719 Hz，基于欧氏距离修正后的频率均值为8.149、15.712、31.696 Hz，基于巴氏距离方法修正后的频率均值为7.802、14.024、29.246 Hz.采用Wasserstein距离的方法在修正后， 其误差明显小于其他方法，且修正后的参数值更接近真实值，显示出更高的修正精度.相比之下，欧氏距离和巴氏距离方法的修正误差较大，尤其是巴氏距离方法在三阶模态频率上的误差最为显著.结果表明，Wasserstein距离方法在模态频率修正方面表现最优，能够有效提高参数的准确性.

6 结 论

本文提出了一种结合Wasserstein距离的近似贝叶斯有限元模型修正方法. 结论如下：

1）构建了径向基模型替代有限元模型， 通过粒子群优化算法确定径向基模型的超参数配置， 构建的代理模型在保持有限元模型特征的同时， 兼具良好的拟合精度与泛化能力. 相较在每次采样时使用有限元模型进行迭代计算的传统方法， 显著提升了计算效率.

2）引入Wasserstein距离构建近似似然函数，实现无须显式计算似然函数的贝叶斯模型修正.由于Wasserstein距离的对称性以及无须支撑集重叠的特性，避免了欧氏距离因非对称性导致的先验分布偏移敏感性问题，从而提升了模型修正的稳定性与可靠性.

3）使用近似贝叶斯-序贯蒙特卡罗进行似然求解， 通过动态调整提议分布和采用自适应重采样机制， 实现对参数后验分布的高效探索与快速收敛.

参考文献

[1] EREIZ S， DUVNJAK I， FERNANDO JIMÉNEZ-ALONSO J ．Review of finite element model updating methods for structural applications［J］．Structures， 2022， 41：684-723．

[2] 刘纲， 谭帅帅， 邹春蓉， 等 ．多链MCMC有限元模型修正收敛判定方法研究［J］. 湖南大学学报（自然科学版）， 2024， 51（11）： 85-93．

[3] LIU G， TAN S S， ZOU C R， et al ．Convergence determination method for multi-chain MCMC finite element model modification［J］．Journal of Hunan University （Natural Sciences）， 2024， 51（11）：85-93．（in Chinese）

[4] 陶言和， 郭勤涛， 周瑾， 等 ．观测不确定性下变分贝叶斯高效模型修正［J］．航空学报， 2024， 45（19）：187-201．

[5] TAO Y H， GUO Q T， ZHOU J， et al ．Efficient variational Bayesian model updating under observation uncertainty［J］．Acta Aeronautica et Astronautica Sinica， 2024， 45（19）：187-201．（in Chinese）

[6] 王未寅， 王佐才， 辛宇， 等 ．基于模块化贝叶斯推理的随机非线性模型修正［J］．振动与冲击， 2023， 42（2）：79-88．

[7] WANG W Y， WANG Z C， XIN Y， et al ．Stochastic nonlinear model updating based on modular Bayesian inference［J］．Journal of Vibration and Shock， 2023， 42（2）：79-88．（in Chinese）

[8] 姚凯. 基于高斯过程的结构动力不确定性量化及有限元模型修正研究［D］.南昌：南昌大学， 2022.

[9] YAO K. Study on quantification of structural dynamic uncertainty and finite element model correction based on Gaussian process［D］．Nanchang：Nanchang University，2022. （in Chinese）

[10] BECK J L， KATAFYGIOTIS L S ．Updating models and their uncertainties．I：Bayesian statistical framework［J］．Journal of Engineering Mechanics， 1998， 124（4）：455-461．

[11] BI S F， BEER M， COGAN S， et al ．Stochastic model updating with uncertainty quantification：an overview and tutorial［J］．Mechanical Systems and Signal Processing， 2023， 204：110784．

[12] PRITCHARD J K， SEIELSTAD M T， PEREZ-LEZAUN A， et al ．Population growth of human Y chromosomes：a study of Y chromosome microsatellites［J］. Molecular Biology and Evolution， 1999， 16（12）： 1791-1798．

[13] NUNES M A， BALDING D J ．On optimal selection of summary statistics for approximate Bayesian computation［J］．Statistical Applications in Genetics and Molecular Biology，2010，9：34.

[14] KHODAPARAST H H， MOTTERSHEAD J E. Efficient methods in stochastic model updating ［C］//Proceedings of ISMA2008. Leuven： ISMA， 2008： 1855-1869.

[15] BI S F， BROGGI M， BEER M ．The role of the Bhattacharyya distance in stochastic model updating［J］．Mechanical Systems and Signal Processing， 2019， 117：437-452．

[16] 张亚峰， 彭珍瑞， 张雪萍， 等 ．基于径向基模型和巴氏距离的随机有限元模型修正［J］．振动与冲击， 2021， 40（19）：221-229．

[17] ZHANG Y F， PENG Z R， ZHANG X P， et al ．Stochastic finite element model updating based on radial basis model and Bhattacharyya distance［J］．Journal of Vibration and Shock， 2021， 40（19）：221-229．（in Chinese）

[18] PANARETOS V M， ZEMEL Y ．Statistical aspects of Wasserstein distances［J］．Annual Review of Statistics and Its Application， 2019， 6：405-431．

[19] 杨乐昌， 韩东旭， 王丕东 ．基于Wasserstein距离测度的非精确概率模型修正方法［J］．机械工程学报， 2022， 58（24）：300-311．

[20] YANG L C， HAN D X， WANG P D ．Imprecise probabilistic model updating using a Wasserstein distance-based uncertainty quantification metric［J］. Journal of Mechanical Engineering， 2022， 58（24）：300-311．（in Chinese）

[21] BERNTON E， JACOB P E， GERBER M， et al ．On parameter estimation with the Wasserstein distance［J］. Information and Inference：A Journal of the IMA， 2019， 8（4）： 657-676．

[22] TAO Y H， GUO Q T， ZHOU J， et al ．An efficient model updating method based on variational Bayesian inference with Wasserstein distance metric［J］．Journal of Mechanical Science and Technology， 2024， 38（11）：5949-5959．

[23] 郁胜， 周林仁， 欧进萍 ．基于径向基函数响应面方法的超大跨悬索桥有限元模型修正［J］．铁道科学与工程学报， 2014， 11（1）：1-9．

[24] YU S， ZHOU L R， OU J P ．Finite element model updating of large suspension bridge based on radial basis function response surface［J］．Journal of Railway Science and Engineering， 2014， 11（1）：1-9．（in Chinese）

[25] 王增辉， 殷红， 彭珍瑞， 等 ．基于加速度频响函数小波变换的贝叶斯模型修正［J］．振动与冲击， 2022， 41（10）：30-39．

[26] WANG Z H， YIN H， PENG Z R， et al ．Bayesian model updating method based on the wavelet transform of acceleration frequency response functions［J］．Journal of Vibration and Shock， 2022， 41（10）： 30-39．（in Chinese）

[27] 翁顺， 雷奥琦， 陈志丹， 等 ．基于模型嵌入循环神经网络的损伤识别方法［J］．湖南大学学报（自然科学版）， 2024， 51（7）：21-29．

[28] WENG S， LEI A Q， CHEN Z D， et al ．Model-embedding based damage detection method for recurrent neural network［J］. Journal of Hunan University （Natural Sciences）， 2024， 51（7）：21-29．（in Chinese）

[29] JIANG Q H， ZHU L L， SHU C， et al ．An efficient multilayer RBF neural network and its application to regression problems［J］．Neural Computing and Applications， 2022， 34（6）：4133-4150．

[30] BEAUMONT M A ．Approximate Bayesian computation in evolution and ecology［J］．Annual Review of Ecology， Evolution， and Systematics， 2010， 41：379-406．

[31] 彭珍瑞， 郑捷， 白钰， 等 ．一种基于改进MCMC算法的模型修正方法［J］．振动与冲击， 2020， 39（4）： 236-245．

[32] PENG Z R， ZHENG J， BAI Y， et al ．A model updating method based on an improved MCMC algorithm［J］．Journal of Vibration and Shock， 2020， 39（4）：236-245．（in Chinese）

[33] SISSON S A， FAN Y， TANAKA M M ．Sequential Monte Carlo without likelihoods［J］．Proceedings of the National Academy of Sciences of the United States of America，2007，104（6）：1760-1765．

[34] RUBNER Y， TOMASI C， GUIBAS L J ．The Earth mover’s distance as a metric for image retrieval［J］．International Journal of Computer Vision， 2000， 40（2）：99-121．

[35] WANG B K， SUN W L， WANG H W， et al ．Research on rapid calculation method of wind turbine blade strain for digital twin［J］．Renewable Energy， 2024， 221：119783．

基金资助

国家自然科学基金资助项目(51768035)

National Natural ScienceFoundation of China(51768035)

甘肃省自然科学基金资助项目(24JRRA172)

Natural Science Foundation of Gansu Province(24JRRA172)