多链MCMC有限元模型修正收敛判定方法研究

Convergence Determination Method for Multi-chain MCMC Finite Element Model Modification

摘要

针对多链MCMC算法主要依靠设置较多迭代步以生成足量后验样本的问题，提出了算法收敛及所需后验样本数量的自动判定方法 .在MCMC算法迭代初期引入最优竞争策略，通过定向差分代替原有随机差分，从而使多链样本快速向目标方向移动以提升计算效率；将当前迭代步多链带宽与预设精度指标相比较，自动判定由定向差分向随机差分的转换时机；基于抽样分布定理，采用一段区间内的样本构造t分布判定指标，对多链MCMC算法是否收敛进行自动判断，并在样本量满足统计要求后自动终止算法，以大幅减少算法平稳期的计算工作量.数值算例及实桥修正结果表明：在相同的计算精度下，所提判定方法可使多链MCMC算法燃烧期计算效率提高30%，考虑预设步长的情况下整个迭代过程可提速约50%，为基于贝叶斯的有限元模型修正在大型土木工程中的应用提供了方法支撑.

Abstract

An automatic determination method for the convergence and required posterior samples is proposed to address the problem that the multi-chain MCMC algorithm mainly relies on setting a large number of iteration steps to generate a sufficient number of posterior samples. The optimal competition strategy was introduced at the beginning of the MCMC algorithm iteration， replacing the random difference by the directional difference， so that the samples could be moved to the target direction quickly to accelerate the computational efficiency. Based on the sampling distribution theorem， a t-distribution determination index was constructed using samples within a period of time to automatically determine whether the multi-chain MCMC converged and automatically terminated the algorithm after the sample size meets the statistical requirements， so as to reduce the computational workload of the smooth period. The numerical examples and the results of the real bridge correction show that the proposed determination method can improve the computational efficiency of the multi-chain MCMC algorithm by 30% with the same calculation accuracy， and the entire iteration process can be accelerated by about 50% considering the preset step size， which provides a method to support the application of Bayesian-based finite element model correction in large civil engineering.

Graphical abstract

关键词

建筑结构 / 有限元模型修正 / 贝叶斯估计 / 马尔可夫链 / 判别指标

有限元分析已成为土木、机械及航空航天等领域进行数值模拟的主流手段^［1］. 大量实践表明，依靠设计或竣工图纸建立的初始有限元模型与实际工程结构并不完全相符，这主要是因为有限元模型建立过程中存在连接及约束简化、材料参数偏差、施工误差等因素，而模型修正是消除有限元模型与实际工程结构不一致的有效途径^［2］. 基于贝叶斯的有限元模型修正将修正参数视为随机变量，结合实际经验推导并计算其后验概率分布，成功解决了有限元模型中的不确定性问题，近年来逐渐成为有限元模型修正的重要分支^［3-4］.

在土木工程中，为解决贝叶斯方法无法直接建立似然函数显式表达式的问题，学界引入马尔可夫链蒙特卡罗（Markov chain Monte Carlo，MCMC）算法，通过成千上万次抽样逐步迭代以逼近修正参数的真实分布^［5-6］，并将从初始样本至样本总体平稳的波动阶段称为燃烧期，将样本在固定范围内游走的阶段称为平稳期. Gong等^［7］通过数学推导证明当有效样本量足够大时，马尔可夫链必然达到收敛状态. 因此，长期以来MCMC算法主要根据经验人为预先设定迭代次数进行大量寻优计算，随后采用人工方式对样本是否收敛以及燃烧期、平稳期的界限进行判定，存在计算时间长、客观性不足等缺点^［8］. 为解决这些问题，Geyer^［9］提出了平均收敛图法，通过计算每个样本与之前样本的平均距离，根据平均值图像是否稳定来判定收敛；Yu等^［10］构建了消费路径图法，通过计算单个样本与后验样本均值的累计误差，根据计算值与0刻度的接近程度以及计算曲线平滑度来判定收敛；Geweke^［11］将样本划分为数量相等的两部分，将其均值是否相等作为MCMC算法收敛的判定准则.

近年来，为解决大型工程结构高维修正参数问题，多链MCMC算法已逐渐成为主流手段^［12］. 该算法初始抽取多个样本，各样本同步迭代形成多条马氏链，通过广泛搜索解空间解决高维修正参数的收敛问题. 然而多链运算导致计算量大幅攀升，特别是针对大型有限元模型，往往需数日才能完成修正计算；同时，早期提出的单链收敛判别方法已不再适用. 针对这一问题，Jones等^［13］根据样本平均值的置信区间宽度与指定阈值的比值大小进行收敛判定，通过自动终止迭代进程计算减少计算量；Gelman等^［14］提出了方差比法，通过计算各链样本的方差平均值与平均值方差，根据加权组合得到后验分布方差的无偏估计进行收敛判定.

上述方法仅能判定MCMC算法是否收敛，且计算复杂导致算法效率不高；同时，未对燃烧期、平稳期做出严谨划分，依靠人工手段无法消除参数统计的主观性. 鉴于此，本文提出一种多链MCMC算法收敛的判定方法. 通过构造t分布模型实现多链收敛自动判定，并在收敛后基于样本信息自动计算迭代次数以获取适量后验样本，避免因预设步数过大导致的计算资源浪费. 同时，根据定向差分具有快速收敛优势，在燃烧期采用定向差分替代随机差分以提升计算效率，而在平稳期保留原随机差分以保证MCMC算法的总体收敛性.

1 MCMC有限元模型修正基本原理

在土木工程中，常将有限元模型中单元的弹性模量、质量密度以及支座刚度等参数视为修正参数θ，并将所有修正参数组集为向量 θ =［θ₁，θ₂，…，θ_n ］^T， n为修正参数的个数. 当没有测试信息时，修正参数的概率分布记为π（ θ ），称为先验分布.

在测得结构响应信息或识别得到结构动力特性后，为便于计算，通常假定有限元模型的计算值 y （ θ ）与结构实测值 y_t之间的误差服从均值为0、协方差为Σ的高斯分布G（0， Σ）^［15］：

\begin{array}{l}p(\mathit{y}|\mathit{\theta })=\frac{\mathrm{1}}{{(\mathrm{2}\mathrm{\pi }\mathit{\Sigma })}^{\mathrm{1}/\mathrm{2}}}\\ \mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{[{\mathit{y}}_{\mathrm{t}}-\mathit{y}{(\mathit{\theta })]}^{\mathrm{T}}[{\mathit{y}}_{\mathrm{t}}-\mathit{y}(\mathit{\theta })]}{\mathrm{2}\mathit{\Sigma }}\right\}\end{array}（1）

式中：p(\mathit{y}|\mathit{\theta })表示 θ 给定下的条件分布，又称似然函数.

根据贝叶斯原理，通过似然函数对先验概率分布进行更新，可得修正参数的后验分布p（ θ | y ）^［16］：

p(\mathit{\theta }|\mathit{y})=\frac{p(\mathit{y}|\mathit{\theta })\pi (\mathit{\theta })}{{\int }_{\mathit{\theta }}p(\mathit{y}|\mathit{\theta })\pi (\mathit{\theta })\mathrm{d}\mathit{\theta }}\propto p(\mathit{y}|\mathit{\theta })\pi (\mathit{\theta })（2）

由式（2）可知，后验分布可直接通过积分进行求解，但由于土木工程结构尺寸大、构造复杂且修正参数维度高，往往无法直接建立式（1）中修正参数与有限元模型计算值 y （ θ ）之间的显式表达式，故常采用多链MCMC算法进行计算. 同时，因式（1）中协方差Σ往往未知，故在MCMC算法中定义目标函数F为：

F(\mathit{\theta })=[{\mathit{y}}_{\mathrm{t}}-\mathit{y}{(\mathit{\theta })]}^{\mathrm{T}}[{\mathit{y}}_{\mathrm{t}}-\mathit{y}(\mathit{\theta })]（3）

则MCMC算法的主要步骤可表述为：

1） 根据修正参数先验概率分布π（ θ ）生成m组初始样本{\mathit{x}}_{\mathrm{0}}^i（i=1， 2，…， m），其中m为马尔可夫链的链数，下标0表示初始样本.

2） 将初始样本{\mathit{x}}_{\mathrm{0}}^i代入有限元模型进行计算，得到m组有限元模型计算值{\mathit{y}}_{\mathrm{0}}^i\left({\mathit{x}}_{\mathrm{0}}^i\right).

3） 对第i条链样本进行交叉更新，得到迭代后的候选样本{\mathit{x}}^{i\mathrm{p}}：

{\mathit{x}}^{i\mathrm{p}}={\mathit{x}}_{\mathrm{0}}^i+\gamma \left({\mathit{x}}_{\mathrm{0}}^{R_{\mathrm{1}}}-{\mathit{x}}_{\mathrm{0}}^{R_{\mathrm{2}}}\right)+\mathit{e}（4）

式中：i表示第i条链的候选样本；R₁、R₂均为随机整数，分别表示第R₁、R₂条链，1≤R₁， R₂≤m且R₁≠R₂≠i；γ为比例因子，常取\gamma =\mathrm{2.38}/\sqrt[]{\mathrm{2}n}； e 为随机游走项， e\inG （0， c），c一般取10^-6.

4） 将候选样本{\mathit{x}}^{i\mathrm{p}}代入有限元模型进行计算，得到m组有限元模型计算值{\mathit{y}}^{i\mathrm{p}}\left({\mathit{x}}^{i\mathrm{p}}\right).

5） 判断是否接受候选样本，即计算接受概率{\mathit{r}}^i：

{\mathit{r}}^i=F\left({\mathit{y}}^{i\mathrm{p}}\right)/F\left({\mathit{y}}_{\mathrm{0}}^i\right)（5）

当{\mathit{r}}^i>1时，令{\mathit{x}}_{\mathrm{1}}^i={\mathit{x}}^{i\mathrm{p}}，否则{\mathit{x}}_{\mathrm{1}}^i={\mathit{x}}_{\mathrm{0}}^i，其中下标1表示第1步迭代.

6） 对所有的m条马尔可夫链重复以上步骤2）~5），完成本步内各链样本更新，得到1次更新后的后验样本{\mathit{x}}_{\mathrm{1}}^i（i=1，2，…，m）.

7） 将后验样本 x₁视为先验样本，重复步骤3）~6），直到达到设定的总迭代步数N. 将所得计算结果按迭代步数排列，则可得到马尔可夫链：{\mathit{x}}_{\mathrm{1}}\to {\mathit{x}}_{\mathrm{2}}\to \cdots \to {\mathit{x}}_t\to {\mathit{x}}_{t+\mathrm{1}}\to \cdots \to {\mathit{x}}_N.

8） 检验所得马尔可夫链的稳定性，若判定其从第t_b步开始平稳，则选取该步及以后的所有平稳样本进行后验分布的统计推断.

从MCMC计算的步骤可知，先验分布及初始样本往往是人为或随机选取的，与后验分布往往有较大差距，导致样本在迭代初期会出现大幅波动现象，即样本不平稳，称该阶段为燃烧期.若将燃烧期内产生的样本纳入计算势必造成较大误差，故业界通常将该阶段内样本抛弃，仅选取处于平稳期的样本进行统计. 但目前燃烧期与平稳期的判定通常凭借人工经验，导致所得后验分布客观性不强.更重要的是，样本何时达到平稳期无法事先设定，故为保证平稳期的样本数量满足统计推断要求，业界往往增加马尔可夫链计算步数N，导致计算量大幅攀升.

2 多链MCMC快速收敛判定算法

针对多链MCMC算法难以自动判定燃烧期和平稳期的问题，根据抽样分布定理，提出样本收敛判定方法及平稳期计算步数的自动确定方法；为提升燃烧期马尔可夫链计算效率，在燃烧期的初始计算期间采用定向差分算法代替随机差分算法，从而形成基于多链MCMC的快速有限元模型修正新方法.

2.1 马尔可夫链收敛自动判定算法

大量研究结果表明，土木工程中单元弹性模量、质量密度、支座刚度等修正参数的后验分布满足正态分布，即所有马尔可夫链上第k维（即第k个修正参数，k=1， 2， …， n）收敛后的样本分布服从均值为{\mathit{\mu }}_k、方差为{\mathit{\sigma }}_k^{}的高斯分布\mathit{G}({\mathit{\mu }}_k,{\mathit{\sigma }}_k^{\mathrm{2}}). 根据统计假设检验原理^［17］，此时可采用区间分析对样本是否满足高斯分布进行判定. 对于第i条链上的第k维参数，当马尔可夫链计算到第t步时，截取链上第t-t_L 步到第t步共t_L （偶数）个样本为分析区间，并将t-t_L 步到t-t_L /2步内样本记为子区间{\mathit{I}}_k^{\mathrm{1}i}，将t-t_L /2+1步到t步的样本记为子区间{\mathit{I}}_k^{\mathrm{2}i}，计算两子区间之差：

{\mathit{I}}_k^i={\mathit{I}}_k^{\mathrm{1}i}-{\mathit{I}}_k^{\mathrm{2}i}（6）

计算所有马尔可夫链第k维参数的子区间之差，并将结果随机组成集合{\mathit{I}}_k^{}. 此时{\mathit{I}}_k^{}中共有L=t_Lm/2个样本，其方差估计值为：

S_k^{\mathrm{2}}=\frac{\mathrm{1}}{L-\mathrm{1}}{\displaystyle \sum _{j=\mathrm{1}}^L}(I_{k,j}^{}{)}^{\mathrm{2}}（7）

由抽样分布定理可知{\mathit{I}}_k^{}/{\sigma }_k^{}~G(\mathrm{0,1})，当用样本方差S_k^{\mathrm{2}}代替总体方差{\sigma }_k^{\mathrm{2}}时，其中L个随机变量的平方和构成的随机变量服从自由度为L-1的卡方分布：

Q_k=\frac{(L-\mathrm{1})S_k^{\mathrm{2}}}{{\sigma }_k^{\mathrm{2}}}\sim {\chi }^{\mathrm{2}}(L-\mathrm{1})（8）

通过标准高斯分布{\mathit{I}}_k^{}/{\sigma }_k^{}和卡方分布Q_k可以构造自由度为L-1的t分布，即

{\mathit{Z}}_k=\frac{{\mathit{I}}_k^{}/{\sigma }_k^{}}{\sqrt[]{Q_k/\left(L-\mathrm{1}\right)}}=\frac{{\mathit{I}}_k^{}}{S_k^{}}~T(L-\mathrm{1})（9）

当马尔可夫链收敛时，Z_k 序列将会围绕0值波动，即Z_k 值满足^［18］：

{\mathit{{\rm Z}}}_k\in (-u_{\alpha /\mathrm{2}},+u_{\alpha /\mathrm{2}})（10）

式中：α为显著性水平，假设检验中通常取5%；u表示检验临界值，当t分布自由度大于120时，u_{\mathrm{0.05}}=\mathrm{1.96}. 此时Z_k 中的样本数量L也应大于120，即分析区间宽度t_L 的取值为：

t_L\ge ⌈\mathrm{240}/m⌉（11）

式中：⌈·⌉表示向上取整的偶数.如判定未收敛，则令 t=t+t_L，进入下一次判定流程.

2.2 平稳期计算步数自动确定算法

马尔可夫链稳定后仍需迭代计算，以满足从样本中统计推断得出修正参数后验分布的要求. 根据统计学原理，对于方差为{\sigma }^{\mathrm{2}}的高斯分布样本，当样本数量达到N_c时，即可保证均值误差在95%置信水平下小于1%：

N_{\mathrm{c}}=(\mathrm{100}{u}_{\mathrm{0.05}}{\times \sigma )}^{\mathrm{2}}（12）

考虑到多链MCMC算法有m条链，同一修正参数在收敛后分布特性相同，故各链样本均可纳入统计计算，则第k个修正参数判定收敛后，平稳期需迭代计算的步数t_ck 为：

t_{\mathrm{c}k}^{}=\mathrm{i}\mathrm{n}\mathrm{t}\left(N_{\mathrm{c}}/m\right)（13）

式中：int（•）表示向上取整.

因有n个修正参数，故MCMC算法稳定后需要迭代的次数t_c应取马尔可夫链各维所需迭代的最大值，即t_{\mathrm{c}}=\mathrm{m}\mathrm{a}\mathrm{x}(t_{\mathrm{c}k}).

2.3 燃烧期高效寻优迭代算法

燃烧期迭代计算的样本不会用于修正参数的统计推断，故减少燃烧期内迭代步数可有效提高计算效率且不影响MCMC算法的计算精度，特别当修正参数维数较高时，可有效解决过长燃烧期导致的计算效率低下问题. 鉴于此，在迭代初期引入最优竞争策略，采用定向差分代替式（4）中的随机差分，以确保每次更新均向较优位置迭代：

{\mathit{x}}_{t+\mathrm{1}}^{i\mathrm{p}}={\mathit{x}}_t^i+{\gamma }^{\mathrm{*}}({\mathit{x}}_t^{\mathrm{B}\mathrm{E}\mathrm{S}\mathrm{T}}-{\mathit{x}}_t^i)+{\mathit{e}}^{\mathrm{*}}（14）

式中：γ^* 为随机比例因子，当迭代步数为偶数时从均匀分布U（0，1）中随机抽取，为奇数时与式（4）中γ取值相同；{\mathit{x}}_t^{\mathrm{B}\mathrm{E}\mathrm{S}\mathrm{T}}表示第t步迭代产生样本中的最优样本，为式（3）计算最小值所对应的样本； e^*为自适应游走项， e^*\inG（0， c^*），c^{\mathrm{*}}={\mathrm{10}}^{-\mathrm{6}t/t_N}.

采用定向差分进行迭代，各链将更快向目标空间搜索，从而大幅提高计算效率. 但该策略将打破细致平衡条件，所形成的马尔可夫链不具备收敛性^［19］. 因此，当样本波动幅度较小时，将定向差分重新转变为原有随机差分. 针对第k个修正参数，设当前步采用定向差分迭代共产生m个样本，选取其中最大值、最小值并计算平均值，分别记为{\mathit{x}}_k^{\mathrm{m}\mathrm{a}\mathrm{x}}、{\mathit{x}}_k^{\mathrm{m}\mathrm{i}\mathrm{n}}、{\mathit{x}}_k^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}，则第k个修正参数的带宽R_k 归一化后可表示为：

R_k=\frac{{\mathit{x}}_k^{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathit{x}}_k^{\mathrm{m}\mathrm{i}\mathrm{n}}}{{\mathit{x}}_k^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}}\le \epsilon（15）

式中：ε表示预设精度指标. 该指标的选取将在第3节进行讨论.

因有n个修正参数，故马尔可夫链当前总带宽应取各维参数带宽的最大值，即R=\mathrm{m}\mathrm{a}\mathrm{x}(R_k^{}).

本文所提算法与传统算法的流程对比如图1所示.

3 数值模拟算例

3.1 搜索性能检验

为了检验算法改进后的全局搜索性能，采用 式（16）、式（17）表示的常用基准测试函数进行检验. 其中Rosenbrock函数（f₁）最优解隐藏在狭长的通道中不易获得，Ackley函数（f₂）通过大量局部极值检验算法跳出局部最优解的能力，函数图像见图2.

f_{\mathrm{1}}{={\displaystyle \sum _{\mathrm{1}}^{d-\mathrm{1}}}[\mathrm{100}(x_{i+\mathrm{1}}-x_i^{\mathrm{2}})}^{\mathrm{2}}+(x_i{-\mathrm{1})}^{\mathrm{2}}（16）

f_{\mathrm{2}}=-\mathrm{20}{\mathrm{e}}^{\left(-\mathrm{0.2}\sqrt[]{\frac{\mathrm{1}}{d}{\displaystyle \sum _{i=\mathrm{1}}^d}{x}_i^{\mathrm{2}}}\right)}+\mathrm{20}-{\mathrm{e}}^{\left(\frac{\mathrm{1}}{d}{\displaystyle \sum _{i=\mathrm{1}}^d}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\mathrm{\pi }{x}_i)\right)}+\mathrm{e}（17）

采用DEMC算法^［20］与改进算法分别对测试函 数进行20次优化计算. 计算时设定目标函数为 min（f），其最优值为0，链数取为60. 每次优化均选择距离函数最优解空间距离最远的坐标点作为初始样本，f₁对应取（-10，-10），f₂对应取（-30，-30），两种算法的平均迭代过程如图3所示.

由图3可知，DEMC算法寻优速率在计算初期便急剧降低，很快便陷入局部最优解；而改进算法的计算精度相较于前者有明显提升，证明该算法能够有效跳出局部最优解并收敛到真实值.

3.2 简支梁模型参数

建立长、宽、高分别为10、0.4、0.6 m的C40混凝土简支梁有限元模型，并沿梁长方向等长度划分为10个单元. 混凝土材料密度为2 500 kg/m³. 假设各单元的弹性模量真实值E_0，i 相互独立，且服从正态分布G_i （E_0，i， 0.1²），其中E_0，i 见表1，E₀=3.25×10⁴ MPa.

根据各单元弹性模量的分布进行抽样，将计算得到的前6阶频率及前3阶振型作为“实测值”，各阶频率见表2，振型见图4. 取各单元弹性模量E_i 与E_0，i 之比作为修正参数，即θ_i =E_i /E_0，i .

3.3 算法效果对比

马尔可夫链数取为8，根据式（11）得区间宽度 t_L =30，取精度指标ε=0.05. 因多链MCMC算法具有一定随机性，故采用DEMC算法、DREAM^［21］算法及本文算法分别进行10次修正运算. 假设先验分布服从均值分布U［0， 3］，即每次修正的初始样本均取0~3之间的随机数. 本文算法某次修正过程中的带宽及检验值变化情况如图5所示.

从图5可以看出，经过374次迭代，总带宽R满足精度要求，算法由定向差分转换为随机差分；随后经过22次收敛判定，即再迭代22×30=660次，Z_k 收缩到检验临界范围内，马尔可夫链判定收敛. 三种算法的平均迭代步数及修正效果对比结果见表3.

由于各参数的迭代收敛规律一致，为节约篇幅，以第5个参数为例，其收敛路径及参数分布对比如图6所示.

结合表3、图6可以看出： DEMC算法在多维修正参数情况下易采样失败导致算法不收敛；DREAM算法虽然能够保证收敛，但迭代收敛成本较高；本文所提算法相较于DEMC算法、DREAM算法能够快速将多链迭代至真实解，并自动判定收敛，算法可减少燃烧期迭代步数30%以上，考虑预设步长的情况下整个迭代过程可提速约50%.

3.4 精度指标ε取值分析

由于有限元模型计算值与结构实测值之间的误差服从均值为0的高斯分布，故各修正参数所满足的分布与结构形式本身并无直接关联，因此本文选择对典型结构开展大量模拟计算，得到ε的建议取值.

分别取不同的ε值各进行10次独立修正，所需燃烧期迭代次数平均值如表4所示. 定向差分、随机差分及燃烧期总步数的平均值如图7所示.

结合表4、图7可知，预设精度指标ε取值过大时，由于指标宽松易满足判定条件，定向差分所需迭代步数较少，但此时各链之间距离较远，无法有效临近后验分布，故随机差分所需迭代步数大幅增加；ε取值过小时，定向差分所需迭代步数增加，且因该阶段马尔可夫链不满足细致平衡条件，各链可能被限制在距真实解一定范围的狭小空间内，导致链间差分值过小，需要一定步数扩展搜索范围，故随机差分所需步数也会增加. 从有效提升计算效率角度出发，要求定向差分、随机差分总迭代步数较少，经过大量模拟计算，本文建议预设精度指标ε为0.05.

4 实桥算例

4.1 桥梁测试信息及初始有限元模型

以重庆市高家花园大桥为例，该桥为预应力混凝土连续刚构桥，主跨径组合为140 m+240 m+ 140 m，主梁形式为三向预应力钢筋混凝土单室箱型梁. 测试主要依靠随机车辆荷载施加激励，利用高灵敏度加速度传感器（941-B型低频传感器）及采集仪记录加速度响应数据. 由于设备限制，测试只对桥梁竖向平动加速度信号进行采集，传感器布置如图8所示.

通过频域分解法获取大桥的模态信息^［22］，经过稳定测试后得到结构的竖向前6阶频率及前3阶振型，频率实测值如表5所示，振型如图9所示.

采用商业有限元软件建立桥梁模型，其中，主梁、桥墩均采用梁单元进行模拟，桥墩底部根据现场实际情况设置为固结，桥梁主梁两端支座均设定为滑动铰支座. 根据设计图纸，初始模型材料参数特性取值见表6，由于整体配筋率较小，初始有限元模型建立时不考虑钢筋作用.

4.2 模型修正

有限元模型沿桥长方向每2 m划分为一个单元，共260个单元. 由于实桥体积庞大，若针对每个单元均设置修正参数，则修正参数过多，易出现不适定问题. 针对这一问题，将若干单元组集形成超单元，采用同一个修正参数进行表达，全桥各单元的组集方式如图10所示，共组集为10个超单元，选取各超单元的弹性模量及质量密度为修正参数，共20个修正参数.

采用本文所提方法对模型参数进行修正，取链数m=10，根据自由度取值要求，区间宽度t_L =24，初始样本取0~3之间的随机数值. 某次迭代过程中某参数收敛路径如图11所示.

由图11可知，该参数在923步时已满足5%预设精度指标，达到转换为随机差分的要求；通过t分布收敛诊断方法对随后区间进行收敛诊断，在1 643步时，马尔可夫链中所有参数均已收敛. 将估计参数代入模型计算获得修正后的频率响应及模态置信因子（MAC），结果如表7所示.

由表7可知，修正后频率的最大误差为3.01%，振型MAC最小值为0.968，各阶频率和振型与实测吻合度均有大幅提升，验证了所提方法在大型工程有限元模型修正中的有效性.

5 结 论

针对传统多链MCMC算法依靠人工判定或大量计算来生成足量后验参数样本、主观误差大、无法自动终止的问题，通过构造t分布指标判定多链收敛性，并在迭代初期通过定向差分代替原随机差分使样本快速向目标方向迭代，结论如下：

1）通过构造t分布指标，采用最新样本区间进行多链收敛判定，剔除迭代过程中陈旧和不稳定的样本，成功突破了传统算法一步一判定的局限性，提高了诊断效率；

2）迭代初期采用定向差分代替随机差分能够使样本快速向目标方向迭代，在5%精度指标下燃烧期能够缩短30%；

3）多链收敛判定后根据置信水平及精度要求进行固定次数迭代，能够有效避免预设迭代次数不合理导致的计算资源浪费，从而提高计算效率.

参考文献

[1] 封周权， 王文赞， 华旭刚， 等 ．基于模态参数与改进萤火虫算法的结构模型修正［J］. 湖南大学学报（自然科学版）， 2022， 49（11）：252-259．

[2] FENG Z Q， WANG W Z， HUA X G， et al. Structural model updating based on modal parameters and modified firefly algorithm［J］．Journal of Hunan University （Natural Sciences）， 2022， 49（11）： 252-259．（in Chinese）

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

[4] 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）

[5] 万华平， 任伟新， 黄天立. 基于贝叶斯推理的随机模型修正方法［J］. 中国公路学报，2016， 29（4）： 67-76．

[6] WAN H P， REN W X， HUANG T L. Stochastic model updating approach by using Bayesian inference［J］. China Journal of Highway and Transport， 2016， 29（4）： 67-76．（in Chinese）

[7] 杨朋超， 薛松涛， 谢丽宇 ．消能减震建筑结构的贝叶斯有限元模型修正［J］. 土木工程学报， 2021， 54（）：13-19．

[8] YANG P C， XUE S T， XIE L Y. Bayesian finite-element model updating of passively controlled building structures［J］. China Civil Engineering Journal， 2021， 54（Sup.1）： 13-19．（in Chinese）

[9] KATAFYGIOTIS L S， BECK J L. Updating models and their uncertainties.Ⅱ： Model identifiability［J］. Journal of Engineering Mechanics， 1998， 124（4）： 463-467．

[10] JONES G L， QIN Q. Markov chain Monte Carlo in practice［J］．Annual Review of Statistics and Its Application，2022，9：557-578．

[11] GONG L， FLEGAL J M ．A practical sequential stopping rule for high-dimensional Markov chain Monte Carlo［J］．Journal of Computational and Graphical Statistics，2016，25（3）：684-700．

[12] 张建新 ．基于贝叶斯方法的有限元模型修正研究［D］．重庆：重庆大学， 2014．

[13] ZHANG J X ．Research on the modification of finite element model based on Bayesian method［D］．Chongqing：Chongqing University，2014．（in Chinese）

[14] GEYER C J ．Practical Markov chain Monte Carlo［J］．Statistical Science，1992，7（4）： 473-511．

[15] YU B， MYKLAND P ．Looking at Markov samplers through cusum path plots：a simple diagnostic idea［J］．Statistics and Computing，1998，8（3）：275-286．

[16] GEWEKE J. Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments［J］. Bayesian Statistics， 1992， 4（1）： 641-649．

[17] LAM H F， HU J， ZHANG F L， et al ．Markov chain Monte Carlo-based Bayesian model updating of a sailboat-shaped building using a parallel technique［J］. Engineering Structures， 2019， 193：12-27．

[18] JONES G L， HARAN M， CAFFO B S， et al ．Fixed-width output analysis for Markov chain Monte Carlo［J］. Journal of the American Statistical Association， 2006， 101（476）： 1537-1547．

[19] GELMAN A， RUBIN D B ．Inference from iterative simulation using multiple sequences［J］．Statistical Science，1992，7：457-472．

[20] 周云， 贾凡丁， 奚树杭 ．基于贝叶斯理论的多模型结构识别的试验研究［J］．湖南大学学报（自然科学版），2018，45（5）：36-45．

[21] ZHOU Y， JIA F D， XI S H ．Experiment research on multi-model structural identification based on Bayesian theory［J］．Journal of Hunan University （Natural Sciences）， 2018， 45（5）： 36-45．（in Chinese）

[22] 方圣恩， 陈杉， 董照亮 ．结构概率损伤识别的改进近似贝叶斯计算［J］．振动工程学报， 2019， 32（2）： 224-233．

[23] FANG S E， CHEN S， DONG Z L ．Improved approximate Bayesian computation for probabilistic damage identification of structures［J］．Journal of Vibration Engineering，2019，32（2）：224-233．（in Chinese）

[24] 胡政发， 肖海霞 ．应用数理统计与随机过程［M］. 北京： 电子工业出版社， 2021： 102-133．

[25] HU Z F， XIAO H X ．Applied mathematical statistics and stochastic process［M］．Beijing：Publishing House of Electronics Industry， 2021： 102-133．（in Chinese）

[26] COSMA I A， ASGHARIAN M ．Principle of detailed balance and convergence assessment of Markov Chain Monte Carlo methods and simulated annealing［J］．Statistics，2008： 1-29．

[27] MEYN S， TWEEDIE R L， GLYNN P W ．Markov chains and stochastic stability［M］．Cambridge：Cambridge University Press，2009：11-15．

[28] BRAAK C J F TER ．A Markov chain Monte Carlo version of the genetic algorithm differential evolution：easy Bayesian computing for real parameter spaces［J］．Statistics and Computing， 2006， 16（3）： 239-249．

[29] 蒋伟， 刘纲. 基于多链差分进化的贝叶斯有限元模型修正方 法［J］. 工程力学， 2019， 36（6）： 101-108．

[30] JIANG W， LIU G. Bayesian finite element model updating method based on multi-chain differential evolution［J］. Engineering Mechanics， 2019， 36（6）： 101-108．（in Chinese）

[31] HIZAL Ç， AKTAŞ E. Probabilistic investigation of error propagation in frequency domain decomposition-based operational modal analysis［J］．Structural Control and Health Monitoring，2021，28（8）： 1-24．

基金资助

国家自然科学基金资助项目（52078084）