基于完全贝叶斯估计的在役地下结构荷载随机反演方法

摘要

提出一种在役地下结构荷载随机反演方法。采用样条函数，将结构上的未知荷载参数化为一系列插值变量；基于贝叶斯框架，融入结构变形观测数据，并构建荷载参数的后验概率分布（PDF）；引入高效采样算法DREAM（differential evolution adaptive metropolis）实现对未知荷载的完全贝叶斯估计。实际案例结果表明：传统确定性反演方法表现出解不稳定的“病态反演”问题，而本方法的荷载期望值则与实测结果吻合较好；本方法得到了未知荷载的完整后验分布，体现出处理反演结果不唯一的优势。

关键词

地下结构; 荷载反演; 贝叶斯估计; 结构健康评估

识别荷载状态是对服役结构进行健康状态评估与行为预测的重要基础，如对在役结构进行数字建模、内力估算或残余承载力评估时均依赖于准确获知其当前荷载状态^［1］。特别地，由于地下服役环境的复杂性，许多在役地下结构荷载超出设计预期，引发超限变形与严重病害，服役安全性不容乐观^［2-4］。因此，感知在役地下结构的荷载状态显得尤为重要。然而，结构埋置于地下，由于监测技术、操作空间有限，直接测量结构外部水土荷载相对困难。特别针对在役结构，如未在结构浇筑前埋设高精度传感器，服役期将不具备相关元件安装与测量条件^［5］。在此情况下，基于其他易观测的结构响应信息，如结构变形数据，反演结构的当前外荷载状态的思路则更为可行。

目前，地下结构荷载反演研究大多针对设计荷载问题。朱合华等^［6］、周济民等^［7］、李策等^［8］、Yan等^［9］以预埋的传感器获取结构内力信息，基于内力与外荷载的关系反演了服役初期地下结构的设计荷载状态。这类研究预先假定荷载分布符合预设的荷载模式，将未知荷载参数化为对应模式下的特定参数，在此受限的参数空间内搜寻与观测数据最吻合的解答，一般可得到唯一且最优反演结果。然而，实测数据表明，大量在役结构的荷载状态早已超出设计预期，分布杂乱^［10-11］。针对未知的杂乱分布的荷载，Gioda等^［12］用多项式逼近未知荷载，将未知荷载分布参数化为多项式上的待定系数，得到观测数据与预测值拟合度最小的解答，将其视为反演结果。实际上，摒弃相关假设（如荷载模式）对参数空间放宽^［12］，反演结果可能会不唯一或不稳定^［13］。许多截然不同的荷载分布可能对应同一组位移观测值，使得反演结果没有唯一解；由于位移-荷载映射矩阵条件数过大，微小的位移观测误差将引起反演荷载的剧烈波动，又称之为“病态反演”。为解决上述问题，学者们通常引入正则化技术，施加一组人为强制性约束（如“平滑”）来惩罚参数空间内的部分成分，从而获得唯一最优解^［14］。张中生等^［15］基于截断奇异值分解正则化技术，反演了算例支架上的荷载，获得较好的反演结果。Liu等^［16］和Liu等^［17］分别采用截断奇异值分解和Tikhonov正则化技术反演了盾构机外壳上的荷载分布，并通过算例说明了该方法对观测误差的免疫效果。实际上，这类正则化约束被认为是一种主观先验信息，相关正则化参数选取通常依赖于研究人员的经验^［14，18］，强制要求反演结果在数学上的“平滑”可能并不适用于所有案例。此外，上述方法均基于确定性反演获得唯一最优解，未对解的不唯一性做出讨论。随机反演适合对解的不唯一性进行描述，朱永全等^［19］基于确定性反演的寻优框架，直接采用Monte Carlo模拟，近似获得了隧道结构上反演荷载的统计值，但该方法有效性有待进一步验证，而且Monte Carlo抽样计算成本过高。何涛^［20］在直接求逆法的基础上，采用概率密度演化理论求得结构上荷载的概率密度函数，但直接求逆法对观测数据要求较高，极易出现“病态反演”问题^［12，15］。随机反演中，贝叶斯方法凭借其完备的统计理论基础，在地下工程中的土参反演^［21-23］、结构损伤识别^［24-25］等方面得到了广泛的应用与发展。基于贝叶斯公式，待演参数的后验分布可以通过先验信息和似然函数求出。先验分布通常基于已知信息被预设为高斯分布等特殊分布形式，构造与似然函数同分布的函数，以获得后验概率的解析解，反演过程高效、简洁；通过高斯先验的强约束作用，解决“病态”反演问题。有学者已经证明，当先验信息取为拉普拉斯分布或高斯分布时，贝叶斯方法的最大后验解与确定性方法中的Lasso正则化或Tikhonov正则化反演结果完全等价^［26］。然而，对于在役地下结构的荷载反演问题，由于其荷载状态通常超出设计预期，因此要求荷载参数服从高斯分布等强先验分布过于主观。

为此，将贝叶斯反演方法引入在役地下结构的荷载反演问题，但对于在役结构荷载反演的具体背景，避免对先验分布进行预设，而采用更为合理的弱先验形式。针对由强先验预设解除带来的无解析解问题，引入自适应差分进化算法DREAM（differential evolution adaptive metropolis）^［27］高效采样，得到荷载后验分布的完全贝叶斯估计；针对由强先验预设解除带来的“病态反演”问题，在所得完整后验分布上积分得到荷载期望值，并依据所有可行解的“概率加权”平滑某一特定解不稳定的特征。

1 地下结构荷载随机反演方法

在役地下结构荷载分布通常杂乱无章，因此在荷载反演前需对其进行参数表征；基于贝叶斯公式融入结构变形观测数据与先验信息，构建荷载参数后验概率分布表达式；在后验概率分布表达式基础上对荷载参数进行完全贝叶斯估计，即可得到反演荷载的完整后验概率分布，进一步积分可得到荷载期望值。

1.1　荷载参数化

针对在役地下结构未知的荷载分布，摒弃相关荷载模式假设，采用简明而灵活的样条插值函数逼近任意荷载分布，即：将未知荷载分布参数化为结构上一系列等间距分布的插值参数。具体地，以广义坐标z表示结构上的区域，区域内以n个均布的插值参数x=（x₁，…，x_n）进行样条插值逼近未知荷载q，借助插值关系表达式即可将未知荷载反演问题转变为中间变量x的反演问题，如下所示：

$$\mathit{q}={\mathit{S}}_{\mathrm{p}}\mathit{x}$$

（1）

式中：S_p为插值算子，其推导见文献［28］。图1为条形结构的参数化示意图（圆弧形结构的插值可在极坐标下进行），未知的荷载分布由插值参数的不同取值而控制。值得一提的是，如图1所示，插值函数的参数化能力由参数个数n决定，即参数个数越多，参数化能力越强。只要n足够大，任何荷载形态均能囊括进此参数空间。实际上，为了数学上处理方便，正如奥卡姆剃刀原理，n的取值应当在囊括真实荷载形态的基础上越小越好。这里，进一步提出一种确定n值的思路，即不断增大n进行试算，当反演的结果几乎收敛不变时，可认为参数化能力无需再提升，真实的荷载已经反演得到，详细流程在2节中讨论。

图1  基于样条插值函数的荷载参数化方法示意图

Fig.1  Illustration example of load parameterization method based on spline function

1.2　贝叶斯反演框架

1.2.1　正演分析模型

预测某一组特定荷载参数x下的结构响应，即正分析过程（“荷载-结构”模型），可由最常见的有限单元法格式表示，如下所示：

$$\mathit{g}(\mathit{x})={\mathit{K}}^{-\mathrm{1}}\mathit{f}={\mathit{K}}^{-\mathrm{1}}{\mathit{L}}_{\mathrm{E}}\mathit{q}={\mathit{K}}^{-\mathrm{1}}{\mathit{L}}_{\mathrm{E}}{\mathit{S}}_{\mathrm{p}}\mathit{x}$$

（2）

式中：$\mathit{g}(\mathit{x})$为结构响应预测值，如结构变形；K为结构离散后的刚度矩阵；f为离散单元节点上的等效荷载，可由荷载分布q及其算子L_E根据虚功原理等效而得。线弹性问题可直接由式（2）求得结构响应，非线性问题则可基于牛顿迭代等方法求解。常见“荷载-结构”模型的有限元法格式推导已十分成熟^［29］，不再赘述。

若实际观测一组结构变形数据d，常规确定性方法通过寻求预测响应与实测变形差值的二阶范数最小解来反演荷载参数，从而获得荷载，如下所示：

$${\mathit{x}}_{\mathrm{M}\mathrm{A}\mathrm{P}}=\underset{\mathit{x}}{\mathrm{m}\mathrm{i}\mathrm{n}}\{(\mathit{g}{(\mathit{x})-\mathit{d})}^{\mathrm{T}}(\mathit{g}(\mathit{x})-\mathit{d})\}$$

（3）

结构工程中，由于K的条件数十分大，因此d上微小的观测误差就会引起反演结果的剧烈波动，此外若K欠定，则反演结果没有唯一解。这时，学者们通常引入正则化技术，人为修正刚度矩阵^［15］或在式（3）右侧加入惩罚因子^［17］来限定x的可行解。与之不同，本研究基于随机分析的思想，在贝叶斯反演框架下反演荷载。

1.2.2　贝叶斯公式

对参数x的随机反演结果用给定观测数据d下的条件概率表示，贝叶斯公式表达式如下所示：

$$P(\mathit{x}\left|\mathit{d}\right.)=\frac{P(\mathit{x})P(\mathit{d}\left|\mathit{x}\right.)}{\int P(\mathit{x})P(\mathit{d}\left|\mathit{x}\right.)\mathrm{d}\mathit{x}}$$

（4）

式中：$P(\mathit{x}\left|\mathit{d}\right.)$是参数x融入后观测数据d（本研究主要讨论结构变形观测数据）的后验概率分布（PDF）；$P(\mathit{x})$是先验概率；$P(\mathit{d}\left|\mathit{x}\right.)$是似然函数；分母又称“边际似然”，是使PDF积分始终为1的归一化常数，在PDF的数值估计流程中通常可以忽略。

1.2.3　先验概率

先验概率以概率密度函数的形式表征参数x的先验信息。任何合理的工程判断均可作为先验信息纳入反演过程。通常而言，对于在役结构的荷载反演问题，由于其荷载状态通常超出设计预期，因此要求荷载参数服从高斯分布等强先验的分布过于主观。

采用一种弱先验形式，即均匀分布。均匀分布在边界内不提供任何信息，认为在融入观测数据前，参数的所有取值均具有相等概率，对其取值没有任何偏好，即x_i~U（q_min，q_max）（i=1，…，n）。据此尽可能依据观测数据推算真实的荷载情况。其中，均匀分布的边界q_min与q_max可在具体案例依据工程经验给出，以过滤掉明显不合理的取值（如某个数量级明显不正常的结果）。

1.2.4　似然函数

似然函数用于表征结构响应预测值与实测位移数据的近似程度，两者差值构造如下所示：

$$\mathit{e}=\mathit{g}(\mathit{x})-\mathit{d}$$

（5）

式中：e为误差向量。对于一些简单的结构而言，由于正演分析模型误差较小，e主要由数据观测误差构成。由中心极限定理，可以假设其为服从均值为零的高斯分布，于是似然函数可构造为

$$P(\mathit{d}\left|\mathit{x}\right.)=\frac{\mathrm{1}}{(\mathrm{2}\mathrm{\pi }{\overline{\sigma }}_{\mathit{e}}^{\mathrm{2}}{)}^{H/\mathrm{2}}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\mathit{e}}^{\mathrm{T}}\mathit{e}}{\mathrm{2}{\overline{\sigma }}_{\mathit{e}}^{\mathrm{2}}}\right)$$

（6）

式中：H为观测数据向量长度；${\overline{\sigma }}_{\mathit{e}}^{}$为观测误差的方差，由数据误差水平决定，可基于数据测量精度预估。

1.3　后验概率分布的完全贝叶斯估计

通常而言，在表征完先验概率与似然函数后，将两者代入式（4）对参数求导，则可以直接求出最大后验解。然而，最大后验解本质上是基于确定性反演的思路，如果缺少强先验信息的支持，就会面临解不稳定的“病态反演”问题^［14］。因此，这里对参数的后验概率分布进行完全贝叶斯估计，采用统计推断的思想，依据完整后验概率分布进一步做出工程决策。特别地，在完整后验分布上积分可进一步得到反演荷载的期望值$\overline{\mathit{q}}$，计算式如下所示：

$$\overline{\mathit{q}}=\int P(\mathit{x}\left|\mathit{d}\right.){\mathit{S}}_{\mathrm{p}}\mathit{x}\mathrm{d}\mathit{x}$$

（7）

根据式（7），可认为荷载期望值是综合所有可行解并按其概率大小“加权”的结果。根据所有可行解的综合推断可以平滑某一特定解大幅波动的特征，解决“病态反演”问题，在地球物理、生命科学等领域已有应用与证明，被描述为“天然正则化特性”^［30-31］。这一点也将直接体现在2节的反演结果中。

对于本问题，由于似然函数与先验概率非同分布，并且考虑到高维参数、反演分析模型可能非线性等问题，后验概率（见式（4））的解析求解几乎不可能。在此情况下，基于马尔可夫链蒙特卡洛方法（MCMC）对后验分布进行完全贝叶斯估计则是一种较为现实的数值方法。该方法的基本思路为：构造一条各态遍历的马尔可夫链^［32］，在马尔可夫链上迭代采样，样本最终将收敛到参数的后验分布。通常而言，“荷载-结构”模型的正分析过程较为耗时，为适应MCMC高频采样的要求，引入更为高效的自适应差分进化算法DREAM^［27］构造马尔可夫链。DREAM并行多条马尔可夫链，通过链间样本相互“学习”，自适应调整采样步长与方向，可快速逼近后验分布。具体步骤如下所示（t代表迭代步数）：

步骤1   t=1时，根据先验概率在N条马尔可夫链上随机产生初始样本${\mathit{x}}^{\mathrm{1},j}=(x_{\mathrm{1}}^{\mathrm{1},j},\dots ,x_n^{\mathrm{1},j}) (j=\mathrm{1},\dots ,N)$。

步骤2   t≥2时，开始进化采样过程，过程如下：

（1）变异产生候选样本，计算式如下所示：

$${\mathit{z}}^{t,j}={\mathit{x}}^{t-\mathrm{1},j}+\epsilon +(\mathrm{1}+\chi )\gamma (\delta ,n_{\mathrm{d}\mathrm{e}\mathrm{f}})•$$

$$\left({\displaystyle \sum _{k=\mathrm{1}}^{\delta }}{\mathit{x}}^{t-\mathrm{1},r_{\mathrm{1}}(k)}-{\displaystyle \sum _{l=\mathrm{1}}^{\delta }}{\mathit{x}}^{t-\mathrm{1},r_{\mathrm{2}}(l)}\right)$$

（8）

式中：x^t-1，j为第j条马尔可夫链上的第（t-1）代样本；ε和χ为保证链条随机遍历的极小随机数； γ（δ，n_def）为跳动步长缩放系数，根据δ和n_def调整；δ为生成候选样本的并行链条对数；n_def为候选样本中被更新参数的数目，若全部更新，则n_def=n；$r_{\mathrm{1}}(k)$和 $r_{\mathrm{2}}(l)$为随机抽取的并行链条编号，且满足$r_{\mathrm{1}}(k)$≠ $r_{\mathrm{2}}(l)$≠j（k，l=1，…，δ）。根据Vrugt等^［27］，相关参数建议取值为ε~N（0，10^-24）；χ~U（-0.05，0.05）；每次迭代时在｛1，2，3｝中随机抽取δ；$\gamma (\delta ,n_{\mathrm{d}\mathrm{e}\mathrm{f}})=\mathrm{2.38}/\sqrt[]{\mathrm{2}\delta n_{\mathrm{d}\mathrm{e}\mathrm{f}}}$。

（2）交叉确定是否接受候选样本分量$z_i^{t,j} (i=\mathrm{1},\dots ,n)$，若u<1-C_R，则不接受（$z_i^{t,j}=x_i^{t-\mathrm{1},j}$，n_def更新为n_def=-1）；反之，接受新值。u是［0，1］均匀分布产生的随机数，定义交叉概率C_R∈［0，1］，$z_i^{t,j}$的表达式如下所示：

$$z_i^{t,j}=\left\{\begin{array}{l}{x}_i^{t-\mathrm{1},j},u<\mathrm{1}-C_R\\ z_i^{t,j},u\ge \mathrm{1}-C_R\end{array}\right.$$

（9）

（3）计算新旧样本后验概率，求得接受概率，以接受概率$\alpha ({\mathit{x}}^{t-\mathrm{1},j},{\mathit{z}}^{t,j})$接受候选样本，x^t，j=z^t，j；否则拒绝，x^t，j=x^t-1，j。$\alpha ({\mathit{x}}^{t-\mathrm{1},j},{\mathit{z}}^{t,j})$的表达式如下所示：

$$\alpha ({\mathit{x}}^{t-\mathrm{1},j},{\mathit{z}}^{t,j})=\left\{\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}\left(\frac{P({\mathit{z}}^{t,j}\left|\mathit{d}\right.)}{P({\mathit{x}}^{t-\mathrm{1},j}\left|\mathit{d}\right.)},\mathrm{1}\right),P({\mathit{z}}^{t,j})>\mathrm{0}\\ \mathrm{0},P({\mathit{z}}^{t,j})=\mathrm{0}\end{array}\right.$$

（10）

退出预烧期（burn-in period）后，链条收敛到后验分布。收敛性判断一般采用Gelman等^［33］给出的指标（收敛因子），如下所示：

$${\widehat{R}}_i=\sqrt[]{\frac{t_j-\mathrm{1}}{t_j}+\frac{N+\mathrm{1}}{N{t}_j}\frac{B}{W}}$$

（11）

式中：t_j为第j（j=1，…，N）条链上的迭代次数；N为总并行链条数；B为链内样本方差在链间的方差；W为链内样本方差在链间的平均值。一般认为${\widehat{R}}_i<\mathrm{1.2}$时，该链条收敛^［33］。

2 方法应用与讨论

2.1　案例背景

本节基于一个实际工程案例说明方法的应用流程并证明方法的有效性。Smethurst等^［34］记录了一个抗滑桩和滑体相互作用的案例。如图2a所示，路基边坡发生滑动，使得抗滑桩向坡脚方向挠曲。图2a中，γ₁、γ₂、γ₃、γ₄分别表示碎石、黏土、软弱夹层、硬黏土的土体重度，（EI）_u和（EI）_d分别表示抗滑桩上端3 m和下端7 m的抗弯刚度。抗滑桩的变形数据由测斜仪记录下来，如图2b所示。本节的目标是通过桩的挠曲变形数据反演桩上的未知净土压力（见图2c），并与文献中记录的实测净土压力进行对比，验证反演方法的有效性。

图2  Smethurst等^［34］记录的实测案例

Fig.2  A field case recorded by Smethurst et al.^[34]

2.2　反演过程演示

（1）参数化。桩上分布的未知净土压力参数化为样条插值函数上的n个均布插值参数x=（x₁，…，x_n）。正如1节讨论的，针对n的取值，本案例进行了一系列试算，如表1所示。

表1  荷载参数化试算

Tab.1  Trial calculations for parameterization

试算编号	n取值
S1	4
S2	8
S3	10
S4	12
S5	14
S6	16

（2）正演分析模型。根据Smethurst等^［34］，桩底土体牢牢将桩底固定住，桩底几乎没有任何位移产生。因此，本案例的正演分析模型可以抽象为一个悬臂梁（见图2c），其有限元法格式推导（见式（2））可参考文献［29］。

（3）先验概率与似然函数。先验概率取均匀分布，其边界值可根据工程经验给出一个弱约束：桩上的水平土压力值（大致估计为14.5 m（埋深）×20 kN•m^-3（土体重度）≈300 kPa）均小于桩最底部的土体自重。于是，先验概率确定为x_i~U（-300， 300）（i=1，…，n）；将实测变形数据（见图2a）与正演分析模型（见式（2））代入式（6），得到似然函数，其中${\overline{\sigma }}_{\mathit{e}}^{}$按数据测量精度取值为1^［34］。

（4）完全贝叶斯估计。在DREAM算法框架下开展荷载参数的完全贝叶斯估计。针对每一个试算，根据Vrugt等^［27］建议，并行链条数N取值与参数个数n保持一致（见表1）。每条链迭代次数取值为10 000，并在迭代中监测收敛性（见式（10））。

（5）反演效果评判指标。引入拟合优度R²，评价反演结果优劣^［12］。R²越接近1，表示反演结果和真实值越接近；当然，当反演结果极差时，R²可能为负数，这时修正为零，表示反演完全失败。R²计算式如下所示：

$$R^{\mathrm{2}}=\mathrm{1}-\frac{{‖{\mathit{q}}_{\mathrm{A}}-{\mathit{q}}_{\mathrm{I}}‖}_{\mathrm{2}}}{{‖{\mathit{q}}_{\mathrm{A}}-\frac{\mathrm{1}}{M_{\mathrm{p}}}‖{\mathit{q}}_{\mathrm{A}}‖‖}_{\mathrm{2}}}$$

（12）

式中：q_A为真实荷载分布向量；q_I为反演荷载分布向量；M_p为向量长度。

2.3　方法验证与讨论

以试算S1展示反演中间结果。如图3所示，在迭代约300步左右，4个参数的收敛因子已经快速收敛并小于阈值1.2，可认为马尔可夫链已经退出预烧期并收敛到后验分布。因此，取后50%样本（后 5 000步）用于参数后验分布的估计。

图3  试算S1中收敛因子演化过程

Fig.3  Evolution process of convergence factor in trail calculation S1

试算S1中4个荷载参数的后验边际分布如图4中对角线所示。图4中，μ和σ表示均值与方差。4个参数的后验边际分布具有峰值，意味着融入变形数据后对参数的定量判断确定性增大。值得注意的是，参数x₁分布相对较为平坦，而且标准差较大，意味着x₁的后验判断不确定性相对较大，这是因为x₁在“悬臂梁”的固定端，完全的位移约束使得结构变形对外部荷载变化不敏感，桩底的荷载反演不确定性很大。图4也给出了参数间的两两相关关系。从图4可以发现，x₂、x₃、x₄两两间体现了强相关关系，强相关关系意味着参数间可能存在多种取值组合来对应同一组观测数据，暗示了解的非唯一性^［35］。

图4  试算S1中荷载参数后验边际分布及相关关系

Fig.4  Posterior marginal distribution of parameters (diagonal) and their correlation in trail calculation S1

试算S1的最终反演结果如图5a所示（图中M表示反演荷载后验期望值）。由概率密度云图定量地给出了反演荷载的完整后验概率分布，即高亮区域为荷载取值高概率区域，反之为荷载取值低概率区域。进一步在后验概率分布上积分（见式（7）），得到反演荷载的后验期望值。试算S1-S6结果分别如图5a-f所示。可以发现，当n较小时（S1-S3），参数化能力较弱，反演荷载的后验期望值形态较为简单，而且其形态随着参数个数的增加而逐渐复杂。当n大于12后（S4之后），荷载后验期望值变化较小。具体地，将S1-S6反演得到的荷载期望值绘制于图6。可以发现，随着参数的增多，后验期望值趋于收敛，最终S5与S6的结果几乎完全重合。据此可以认为，试算S6的参数化能力已经足够囊括真实荷载，可视作最终的反演结果，如图5f所示。

图5  反演荷载后验概率分布云图

Fig.5  Color map of PDFs of inversion pressures

图6  试算S1-S6中的反演荷载后验期望值

Fig.6  Expectation of inversion pressures in trial calculations S1-S6

Smethurst等^［34］记录的真实荷载分布也绘制在图5f中（图中A表示）。总体而言，随机反演得到的荷载后验期望值与真实荷载吻合较好，但值得注意的是，两者在桩底处存在一定差距，这是由桩底位移对荷载变化不敏感所致，位移约束端的荷载反演具有一定困难。总体上，两者拟合优度R${}_M^{\mathrm{2}}$=0.75，证明了随机反演方法于本案例应用的有效性。

针对本案例开展了确定性反演，反演的条件与随机反演保持一致，反演结果如图5f中亮色实线（图中D表示）所示。可以发现，确定性反演解波动十分大，与真实荷载相差甚远，展示出典型的“病态反演”问题（或称“过拟合”现象），而R${}_D^{\mathrm{2}}$=0也表明了确定性反演于本案例中应用的失败。相较而言，随机反演不需要任何主观正则化信息，而凭借其后验期望值体现出“天然正则化”特性，具有较强的工程实用性。需要说明的是，随机反演的结果不限于后验期望值，而在于其完整的后验概率分布。如：基于后验分布（见图5f），进一步可推导桩上的内力分布（见图7a）。这意味着桩上任意位置的内力概率分布均可定量得到（如图7b中A、B、C 3个示意点）。总而言之，随机反演承认解的不唯一性，基于完整概率密度函数进行统计推断，而相关工程决策亦可根据相关概率分布或统计特征值等（如养护维修成本的期望值）进一步开展。

图7  桩上内力（弯矩）的概率分布

Fig.7  PDFs of internal force (bending moment) on pile

3 结论

（1）无需任何主观正则化技术，本方法能够解决确定性反演中常见的“病态反演”问题，反演得到的荷载期望值与实测荷载吻合较好，说明本方法具有较强的工程实用性。

（2）反演结果不限于荷载期望值，还得到了反演荷载的完整后验分布，并建议基于不确定性分析的思想进行进一步统计推断，如推算结构内力的完整概率分布，以相关统计特征值进行工程决策。

（3）本方法于结构固定端的反演结果存在较大不确定性，具体体现为此处的反演期望值与真实荷载存在一定误差，这是由完全的位移约束使得该处的变形响应对荷载变化不敏感所致。

限于篇幅，仅给出了抗滑桩的应用实例，但方法基于通用框架建立，意味着本方法可应用于更多工程案例，为在役地下结构的数字化维养技术提供基础。

作者贡献声明

田志尧：研究构思，算法实现。

宫全美：稿件撰写与审核。

赵 昱：数据采集与处理。

周顺华：方案设计，稿件修订。