基于贝叶斯理论及MCMC-MH算法推演地基土材料阻尼比的概率分布模型

曹艳梅， 李东伟， 张玉玉， 杨 林

(北京交通大学 土木建筑工程学院，北京 100044)

近年来，交通系统对周围结构物及其内部仪器、设备等引发的微振动问题越来越受到重视。由于交通荷载引起的振动在土体中以波动形式进行传播和衰减，因此振动的预测和评估需要首先确定场地土的动力特性[1]，而其中土体动参数的准确与否对振动预测和评估结果影响非常大[2]。根据既有文献[3]，对地面振动预测影响最大的土体参数为土的剪切模量和材料阻尼比。Auersch[4]通过理论分析和试验研究证明了场地振动时波动幅值随距离的衰减指数主要受到材料阻尼和几何阻尼的影响。土体的材料阻尼是由土骨架中固体颗粒之间的摩擦以及固态骨架与孔隙流体之间的相对运动而产生的，体现为波动能量的耗散，而几何阻尼体现为波动的扩散，因此材料阻尼比的大小对振动能量的传播与衰减起到至关重要的作用[5-6]。

目前，材料阻尼比主要是通过测试场地中表面波的空间衰减特性而反演得到，反演过程中场地土通常被模拟为各向同性均质分层弹性半空间[7]。Lai等[8-9]开发了一种基于相位和幅值的回归技术，同时确定土壤的剪切模量和材料阻尼比。该方法通过几何衰减因子的迭代更新来计算几何阻尼，再由衰减系数反演材料阻尼比分布。Rix等[10]通过解耦剪切波波速和材料阻尼比的方法简化了反演地基土动参数的技术，但是仍然需要计算几何衰减因子，后来采用与Lai等开发的技术类似的方式获得衰减曲线，进而获得了地基土材料阻尼比分布。Badsar等基于半功率带宽法提出了一种确定材料阻尼比的新方法，不需要计算几何衰减因子，直接将半功率带宽法应用于土壤响应的波数域，有效地避开了不确定性剪切模量分布对几何衰减因子计算准确性的影响。

然而，由于衰减曲线对小空间尺度和较大深度土壤性质的变化不敏感，以及试验条件的限制和试验误差等不确定性因素，地基土的动参数在一定范围内不是确定的，而是随机分布的，如果采用确定的地基土参数对振动进行预测会出现一定的偏差。Schevenels等[11]研究了地基土动参数的不确定性对地面振动预测的影响，通过贝叶斯理论及与试验频散曲线构造似然函数，引入马尔科夫链得到了后验地基土剪切模量概率分布模型。肖张波[12]结合地质统计学和地球物理反演理论，在贝叶斯理论下基于马尔科夫链的扰动对先验空间随机搜索，利用地震数据约束得到地球储层物理参数的后验概率模型。

为了能够更加准确地预测交通系统引起的场地振动，本文建立地基土阻尼比的先验概率分布模型，根据试验衰减曲线构建似然函数，并将其与正演算法获得的理论衰减曲线进行比较，利用贝叶斯理论和MCMC-MH(Monte Carlo Markov chain-Metropolis)算法得到多组土体参数后验样本数据，进行独立性和收敛性检验，得到有效的土体阻尼比后验概率分布模型。进而应用独立的后验样本计算出自由场地的振动响应，利用核密度估计得到一定置信度的置信区间，通过和试验数据进行对比，验证了土体阻尼比非确定性概率模型的合理性和可靠性。

1 地基土材料阻尼比先验概率模型的建立

1.1 地基土材料阻尼比非高斯概率分布模型

基于Badsar等研究中的多通道表面波分析法(MASW)现场测试及其提出的半功率带宽法，可知场地土振动传递函数的频率-波数谱可表示为

(1)

(2)

图1 某场地土表面波衰减曲线

图2 反演的地基土材料阻尼比分布

从图2可看出，土层材料阻尼比在不同深度具有随机分布的特性，通过地震锥穿透原位测试也显示出材料阻尼比随深度而变化的特点。本文忽略土体参数在水平方向上的空间变异性，仅对土体材料阻尼比随深度变化的不确定性进行研究。由于表面波的振动能量主要集中在地表，因此文中取地表以下6 m的深度范围来考虑材料阻尼比分布的随机性，而6 m以下半空间土体的阻尼比参数与6 m深度处的取值相同。计算时，假定除阻尼比以外的其它土参数均已确定：土体密度在整个地层中均为1 800 kg·m-3，泊松比υ为1/3；剪切波波速Cs从地表至1.5 m深度地层为128 m/s，1.5～4 m深度为175 m/s，4～6 m及半空间范围内为355 m/s。

由于阻尼比随深度具有一定的分散性，且所有土层材料的阻尼比参数不能均用简单的一个随机变量来表示，为了更加合理地模拟土参数分布，本文采用非高斯随机过程X(h,θ)来模拟[13]，其中θ是随机维度中的坐标，h是表征地基土深度的竖向坐标。文中采用24个均匀薄层+底部半空间对土体进行建模，每一薄层的厚度取为lc=0.25 m。

基于随机过程模拟理论[14]，非高斯随机过程X(h,θ)可以离散化为概率密度函数pX(θ)和协方差函数CX(h)。由图2所反演出的材料阻尼比可看出，阻尼比的取值范围在0.03～0.05。由于严格的先验随机土体模型会引起后验土体模型的偏差，且较广的先验概率模型能包含更多可能的土层剖面，所以本文先验概率分布模型的概率密度函数pX(θ)在深度0～6 m内采用0.01～0.08的均匀分布，以便能更好地模拟土层参数随深度的变化。协方差函数与两个取样点间的距离有关，如图3所示，它能较好地模拟土层参数随深度的变化，在地质空间统计中有广泛应用。

图3 非高斯过程的协方差函数

协方差函数CX(h)可以采用平稳的马特恩协方差函数来表示[15]

(3)

式中：KX为第二类修正贝塞尔函数；Γ为伽马函数；lc和v分别为确定随机过程X(h,θ)变化范围和平滑度的严格正参数,文中取lc=0.25 m和v=1；σ为概率密度函数pX(θ)的均方差。

1.2 地基土材料阻尼比先验概率分布模型的生成

对非高斯随机过程X(h,θ)直接进行Karhunen-Loeve分解，得到的随机变量ξk(θ)为互不相关的非高斯随机变量。由于概率分布较难确定不便于MATLAB软件程序生成实现，故文中首先将非高斯随机过程基于Nataf变换理论[16-17]和Gauss-Hermite多项式变换[3]转换为标准高斯随机过程η(h,θ)，其相关函数Rη如图4所示。

图4 标准高斯随机过程的相关函数

根据Karhunen-Loeve分解理论，标准高斯概率分布可以表示为具有随机系数的确定性函数的乘积和的形式，即

(4)

式中：ξk(θ)为标准高斯随机变量；λk和fk(h)分别为相关函数Rη矩阵分解的特征值和特征向量，如图5所示。

图5 相关函数Rη矩阵分解的特征值及特征向量

式(4)中，随着阶数的增高，相应的能量贡献不断降低，为了降低先验模型转为后验模型过程中的计算量，本文先验模型中取截断阶数M=18，这些多元随机变量的概率密度函数为

(5)

通过MATLAB随机生成M个服从标准高斯分布的随机数θi，经式(4)及等概率转换原则[18]获得地基土材料阻尼比的非高斯概率分布模型，即先验概率模型。从先验概率模型中提取4组样本数据，并根据TML-PLM理论及频率波数域-半功率带宽法正演出样本数据对应的4组表面波衰减曲线，如图6所示。因为还未考虑试验数据的约束，从图6中可以看出先验材料阻尼比样本以及对应的先验衰减曲线与图1、图2的试验数据都具有较大的差异性。

图6 土体材料阻尼比先验概率模型的4组样本数据及先验衰减曲线

2 材料阻尼比后验概率模型的建立

2.1 试验数据约束的似然函数

根据贝叶斯理论[19-20]，土体材料阻尼比的后验概率分布模型可以表示为

(6)

式中：π(ξ|E)为受试验数据约束的后验概率分布模型；p(ξ)为不受试验数据约束的先验概率分布模型；L(E|ξ)为似然函数，它是建立在土体模型参数与试验数据的匹配度之上的。基于先验模型中的理论衰减曲线(见图6(b))与试验测得的多组衰减曲线下界的相关性，构建似然函数L(E|ξ)

图7 基于试验数据构建的频散曲线范围

2.2 MCMC-MH算法生成后验概率分布模型

由于似然函数L(E|ξ)是一个与试验数据匹配的函数，因此式(6)不能直接计算后验概率分布模型，本文联合应用马尔科夫链蒙特卡洛方法和Metropolis-Hastings算法[21]来获得符合后验概率π(ξ|E)的土体剖面。根据MCMC方法的原理，首先将图2所示的土体材料阻尼比分布通过NATAF变换生成马尔科夫链的初始状态ξ1，其次通过概率转移函数q(ξj+1|ξj)产生新的状态ξj+1,随后生成均匀分布随机变量U～Uniform[0,1]，若满足式(8)的条件则接受转移状态ξj+1，否则拒绝该转移状态。定义转移函数q(ξj+1|ξj)为

(8)

(9)

式中，均方差σ用来确定马尔可夫链的步长，σ越大马尔可夫链的扰动越大，收敛越快，但状态转移中候选状态的接受率会越低。本文取σ=0.12。

对于式(9)，由于π(ξ)无明确的概率计算公式，故在计算时分两个阶段来判断是否接受候选状态ξj+1。首先使用先验概率p(ξ)来考虑扰动后的新状态ξj+1是否通过第一阶段，定义

(10)

若第一阶段U<α1再进行第二阶段，即使用似然函数L(E|ξ)来考虑是否接受扰动后的新状态ξj+1，此时定义

(11)

两阶段过程中，第一阶段控制马尔科夫链的扰动达到平稳收敛，第二阶α2段等于0或1进行试验数据的约束得到后验概率分布模型。将以上的MCMC-MH算法采用MATLAB编程实现，由此计算出的马尔科夫链前十组的后验材料阻尼比分布，如图8所示。

图8 材料阻尼比的10组后验样本数据

2.3 收敛性和独立性检验

在MCMC-MH算法下，分别从变量ξi对应不同样本组数时的均值和均方差两个方面来判断推演所得的后验概率分布模型的马尔可夫链是否收敛。

(12)

(13)

式中，n为状态ξi的组数，每组ξi包含M个相关变量(M=18)。根据式(12)和式(13)计算出马尔科夫链的均值mξ和均方差σξ，分别如图9(a)和图9(b)所示。从图中可以看出，随着样本组数n的增多，均值mξ和均方差σξ随着马尔科夫链的扰动不断变小，最终逐渐趋于平稳。

平稳随机过程的马尔科夫链可通过相关系数ρξ来描述样本之间的相关性，相关系数ρξ=0的样本即认为是相互独立的后验概率模型。相关系数和相关函数的定义为

(14)

(15)

式中：τ为马尔科夫链上两组样本的间隔步数；Rξ为两组样本的相关函数；N为样本总组数。

图9 变量ξi不同样本组数的均值、均方差以及不同间隔步长的相关系数

根据式(14)和式(15)计算得到相关系数，如图9(c)所示。可以看出，对于本文的算例，样本间隔步数大约在500步时相关系数接近0，此时即认为样本之间相互独立。取马尔科夫链间隔步数大于500步的状态为相互独立的后验概率模型，其中10组相互独立的土体阻尼比后验概率分布如图10所示。

3 地基土材料阻尼比后验概率模型的应用

对于N组独立的土体动参数的后验样本，通过TLM-PML理论可以得到N组地基土自由场竖向位移的传递函数H(x,ω)，即土体表面的竖直位移u(x,ω)与锤击力p(ω)的比值。本文将距振源不同距离D和不同频率f的振动位移传递函数通过核密度估计[22]表示为具有概率密度函数的经验分布，进而用具有一定置信度水平的置信区间来表示。

图10 10组相互独立的后验材料阻尼比样本

为了验证较少数量独立样本的核密度估计是否准确，下面采用非独立的马尔科夫链上的所有后验样本与独立的100组后验样本的核密度估计进行对比，如图11所示，其中实线为独立的100组后验样本数据计算所得，虚线为非独立的25 000组后验样本数据计算所得。从图11中可以发现，分别用收敛的马尔科夫链中提取的100组独立样本与已经收敛的非独立的马尔科夫链后验样本所得的概率密度，数值的大小和形状趋势基本是一致。由于独立样本数量少，因而在两种图形的峰值附近存在较小的差异性，但是这种差异性对较大置信度的双侧置信区间基本无影响。

图11 自由场表面竖向位移传递函数H(x,ω)的核密度估计

根据概率密度函数计算出距振源不同距离D的传递函数H(x,ω)在95%置信度的双侧置信区间，并与Badsar等研究中的试验数据进行对比,如图12所示。

从图12中可以看出,在后验概率模型约束的频率范围15～70 Hz内，95%置信度的置信区间基本上涵盖了所有的试验数据。个别频率的不一致可能是由于试验中存在的不利影响因素和试验数据的处理误差等引起的。因此，通过随机的后验概率土体模型可以在一定频率范围内以一定置信度的置信区间预测地基土的振动传递，可以有效地解决地基土反演问题中存在的非唯一性问题，避免交通荷载环境振动预测中的各种不确定性因素的影响。

从图12中还可以看出，在较低频率下(约低于15 Hz)，预测的置信区间和试验数据差异较大，这是由于低频时，土体中振动波的波长较大，表面波测试无法测得较深的土层，而理论模型也未考虑土层深度大于6 m时材料阻尼比的空间变化，因此可以通过增加表面波测试的空间分辨率获得较低频率情况下的衰减曲线，进而得到低频范围内可靠的后验传递函数置信区间。

当距振源的距离D较远时，如图12(c)(D=100 m)，试验数据的高频阶段高于预测的置信区间，而距离D=10 m和D=50 m情况下均未发现高频阶段的差异性，因此可以判定这种远距离处的高频阶段的差异性主要是由周围环境的振动或噪声等干扰引起的，进一步说明了环境振动远场响应的预测可能会存在一定的误差。可以通过优化表面波测试或使用其他的测试方法来得到更高精度的衰减曲线，可以获得更大的预测频率范围，缩小一定置信度下的置信区间，进一步来减少轨道交通引起的振动预测中的不确定性。

图12 后验概率模型传递函数的95%置信区间(灰色阴影)及试验数据(黑色实线)

4 结 论

本文基于非高斯随机过程建立了地基土材料阻尼比的先验概率模型，通过贝叶斯理论和似然函数以MCMC-MH算法得到土体动参数的后验概率模型，进而以一定置信度的置信区间预测了地基土的振动传递，并与试验数据进行比较，验证了本文提出的土体阻尼比非确定性概率模型的合理性和可靠性。通过本文的研究，得到如下主要结论：

(1)地基土材料阻尼比的随机分布特征对场地振动的传播与衰减具有较大的影响，因此建议在场地振动的研究中考虑进地基土参数取值的非确定性。

(2)本文提出的地基土材料阻尼比后验概率模型的推演方法是一种非常有效的随机振动分析方法，可以有效地解决地基土反演问题中存在的非唯一性问题，能够应用于交通环境场地振动的预测、评估、振动传递和分布的随机性分析和研究中。

(3)基于多通道表面波分析的半功率带宽法可以避免不确定性的剪切模量概率分布对求解衰减系数的影响，能够准确有效地应用于地基土材料阻尼比的反演和随机分析中。

(4)可以通过扰动系数来控制马尔可夫链的步长及收敛速度，但由于马尔科夫链的收敛性需要大量的后验样本数据进行检验，所以建议采用平行马尔科夫链并行运算方法来减少计算机的运行时间，提高计算效率。

(5)后验材料阻尼比土体参数具有很大的随机性，基于多组相互独立的后验样本数据，从概率的角度以一定置信度的置信区间预测地基土的振动传递，避免了不确定因素的影响，对环境振动预测及控制具有重要意义。