0.   引言

在岩土地震工程学中,对于线性和非线性土性参数的估计一直是难题,尤其是原位土性参数的估计。工程上常采用取样进行室内试验来获得土性参数,但取样会造成土体原有结构的破坏,引起土样扰动、孔隙结构与含水量改变、应力卸载等,这些因素使得对土体动力参数的估计有较大偏差,特别是对中等到大变形水平下的动力参数影响更大。随着地震观测技术、电子信息技术的快速发展,越来越多的地震台站网在世界各地建立。这些地震台站网能够获得不同地点大量地震实测记录,这些记录相当于大自然做的大量天然实验,合理地运用这些数据,便可以估计特定场地的土性参数。一般工程钻孔能够得到一个较为粗略的土层剪切波速剖面,结合竖井台站井上井下地震记录,对土层剖面重新分层,能够反演一个更为精细的波速结构,这种方法可以被认为是一种获取原位土性的方法。

在地震场地动力响应分析中,频域中分析常采用等效线性法,Kausel等^[1]以及Huang等^[2]先后对等效线性法进行改进,提出了频率相关的等效线性法及改进方法。除了在频域中分析,也可直接使用数值积分方法在时域中进行分析,该方法可以结合不同的土体本构模型来表征土体的动力行为。Borja等^[3]针对黏土的循环加载行为,提出了边界面塑性模型。Wang等^[4]对该模型进行改进,简化了加卸载准则,提出了简化边界面塑性模型。简化边界面塑性模型能够较好地表征土体非线性动力行为,并且所需参数较少,应用简单。对于土体动力响应模型参数的反演,过去也有一些研究,如Asimaki等^[5]提出了一种基于小波域的全局优化的反演算法来反演土层的剪切波速结构；Tsai等^[6]提出了一种自学习模拟的反分析方法,将场地响应分析和实测数据相结合,建立了一种基于神经网络的土体本构模型。近年来,越来越多的研究开始关注贝叶斯方法在岩土参数估计的应用。Akeju ^[7]提出了一种贝叶斯方法来估计岩土参数之间的相关性,该方法能够应用较少的数据得到较为准确的参数相关性；Seylabi等^[8]提出了一种基于无迹卡尔曼滤波的贝叶斯反演方法,利用竖井台站地震记录反演土体动力模型参数。贝叶斯方法相比于传统方法,不仅能够获得参数的最优估计,还能对不确定性进行定量分析。

在Seylabi和Asimaki的贝叶斯反演方法基础上,发展了一套对深层土动力模型参数的反演方法,本文只介绍线性参数即剪切波速的反演。采用基于简化边界面塑性模型的一维场地动力响应分析方法,使用日本KiK-net竖井台站数据,对剪切波速进行反演。首先使用一个无模型误差的数值算例对算法在土体动力响应模型的适用性进行了验证,然后使用KiK-net实测地震动数据进行了剪切波速结构的反演计算。

1.   研究方法

1.1   一维场地动力响应分析

本研究采用OpenSees（Open System for Earthquake Engineering Simulation）^[9]建立土层动力响应有限元模型,假定满足一维场地响应,采用剪切梁建立土层模型。为表征土体小应变阻尼,对模型施加与刚度成正比的瑞利阻尼,

式中,C为阻尼矩阵,K为弹性刚度矩阵,$a_1$为阻尼系数。单元材料采用“Multiaxial Cyclic Plasticity”材料模型,该材料采用了简化边界面塑性模型。

简化边界面塑性模型的核心是在应力应变关系中没有纯弹性区域,土体的非线性通过模量的折减体现,模量从小应变下的初始模量渐变到完全塑性状态下的模量值。边界面${\rm B}$用于表征完全塑性状态：

式中,${\widehat{\sigma }}^{\prime }$为边界面上的应力像点,应力应变符号中带有记号（$\text{'}$）则表示偏应力或偏应变,$\beta$为边界面应力反向点,$R$为边界面半径,可通过不排水抗剪强度来确定,$R=\sqrt{8/3}{S}_u$。对于不排水的情况下,边界面和屈服面在偏应力平面（$\pi$-平面）的投影为圆,边界面上的应力像点${\widehat{\sigma }}^{\prime }$由当前应力点${\sigma }^{\prime }$和应力反向点${{\sigma }^{\prime }}_{\text{0}}$的映射关系确定,如图1所示。本构方程为^[4]

图  1  边界面应力映射关系

Figure  1.  Mapping rules of stress on bounding surface

下载: 全尺寸图片 幻灯片

式中,$G_{\max }$为最大剪切模量,可由土体密度和剪切波速确定,$G_{\max }=\rho V_s{}^2$,${\dot{\sigma }}^{\prime }$,${\dot{\epsilon }}^{\prime }$为偏应力应变变化率,$H^{\prime }$为硬化模量。

硬化模量可通过以下方程确定：

式中,$H_0$为动力硬化模量,h,m为硬化参数,控制剪切模量的折减速率。

1.2   土体动力模型参数的贝叶斯反演方法

将土体在地震作用下的动力响应模型视为概率状态空间模型,第k个时窗$[t_1^k,t_2^k]$的地表响应预测$y_k$与土层性质参数${\theta }_k$的关系（即量测模型）为

式中,$T$为非线性函数,将土层性质参数${\theta }_k\in {ℝ}^n$和输入加速度序列${\ddot{u}}_{g,t_1^0:t_2^k}$映射到地表响应向量${\text{y}}_k\in$${ℝ}^{n_y}$。假设${\theta }_k$不随时间变化,可表示为高斯随机游走模型,地表响应的实测值和预测值的误差假定满足高斯分布,

式中,$w\in R^n$、$v\in R^{n_y}$均为均值为0的高斯噪声,$w$的协方差为$Q\in R^{n\times n}$,$v$的协方差为$R\in R^{n_y\times n_y}$。为了找到一组最优土层参数${\widehat{\theta }}^{est}$,得到地表响应较为准确的预测,本文采用基于无迹卡尔曼滤波（UKF）技术的贝叶斯估计方法,算法整体框架如图2所示。

图  2  无迹卡尔曼滤波算法

Figure  2.  Algorithm of unscented Kalman filtering

下载: 全尺寸图片 幻灯片

为了避免不符合实际物理意义的估计出现,需要给每个参数加上一定的约束,

对于有约束的随机变量,在使用无迹变换生成的样本点（也称sigma点）也需要满足该约束,本文采用Mandela等 ^[10]提出的带约束的无迹变换,其中sigma点${\widehat{\theta }}_k^{(i)}$和权值$W^{(i)}$的定义如下：

式中,$s_i$为矩阵$S=[\sqrt{P_k^{-}}-\sqrt{P_k^{-}}]\in R^{n\times 2n}$的第i列,${\alpha }_k=\min (X)$,矩阵$X$的每个元素由下式确定：

若${\alpha }_k=\sqrt{n+\kappa }$则为无约束无迹变换,$\kappa$为比例系数,决定了sigma点在均值附近的分布,这里取$\kappa$ = 0。在预测步骤中,参数的估计${\widehat{\theta }}_k^{+}$有可能超出约束,因此,在进入到下一步迭代前,需要进行调整,这里采用Simon^[11]提出的一种概率密度函数截断法。

2.   数据来源及处理

本文中,地表响应计算模型采用一维模型,因此,需要选择较符合一维波传播假定的测站进行反演计算。Thompson等^[12]对KiK-net 100个测站的场地响应建模的复杂性进行了研究,提出了一种测站分类方法。他们根据测站不同地震记录的传递函数的变异性和实测记录传递函数与一维假定的理论传递函数的匹配度,选出16个变异性小、匹配度高的测站,经验证这16个测站能够较好地满足一维波传播假定,使用一维计算模型能够较为准确地估计出地表响应。在本文中,选择其中IWTH08测站进行土体动力模型参数的反演计算,该测站的位置及土层剖面如图3所示。该测站位于日本近海岸,土层30 m内的等效剪切波速V_s30为305 m/s,土层岩性较单一,岩基上有约30 m左右的软土覆盖层。

图  3  IWTH08测站位置及其土层剖面

Figure  3.  Location and soil profile of IWTH08

下载: 全尺寸图片 幻灯片

KiK-net中包含了大量地震记录,每个测站均有井上和井下三向加速度时间序列^[13]。从该测站中获取井上地震动PGA小于0.1g的地震动数据用于反演剪切波速结构,以确保土体动力响应处于线性阶段。对于每个地震记录,根据震源和测站的方位角,以及地震记录的东西分量和南北分量获取井上井下地震动SH 分量,然后将井下地震动加速度时间序列作为地底地震动输入,井上地震动加速度时间序列作为地表加速度响应进行反演计算。

3.   剪切波速反演算法验证研究

假定KiK-net提供的IWTH08测站的剪切波速结构为真实波速结构,并假定刚度瑞利阻尼系数a₁= 0.001,则土层参数向量的真实值为

式中,V_s为剪切波速,反演初始值设定为${\widehat{\theta }}_0^{\text{+}}\text{=}1.3{\theta }^{\text{true}}$,参数的初始协方差$P_0^{+}=\text{diag{(}0.15{\widehat{\theta }}_0^{+}{)}^2\}$,过程噪声协方差$Q=\text{diag}\left\{{({10}^{-4}{\widehat{\theta }}_0^{+})}^2\right\}$,量测噪声协方差$R=$${(0.001g)}^{2}I$。其中,$\text{diag}\left\{\text{x}\right\}$表示以向量$\text{x}$为对角线元素的对角矩阵,$g$为重力加速度（g$=9.81$ m/s²）,$\text{I}$为单位矩阵。

该算例的地底地震动输入如图4所示,该地震波的PGA为0.0028g,Arias强度5%～95%的时段为26～36 s。使用${\theta }^{\text{true}}$作为土层参数计算得到的地表加速度响应作为反演的匹配目标,以消除模型误差的影响。时窗起点$[t_1^1,t_2^1]=[24\text{s},26\text{s}]$,每次迭代时窗移动0.1 s,地震波取点时间间隔dt= 0.01 s,因此地表加速度响应向量${\text{y}}_k\in {ℝ}^{200}$。

图  4  地底地震动加速度时间序列

Figure  4.  Time histories of acceleration of down hole ground motion

下载: 全尺寸图片 幻灯片

图5为土层各个参数后验估计的均值比真实值随着迭代过程的变化曲线,灰色区域为一个标准差内的区域。从该图可以看到,随着迭代的进行,每个参数不断地向真值（红色水平线）逼近,并且标准差也不断减小（灰色条带逐渐变窄）,说明不确定性不断减小。最终各个参数均收敛到真值。

图  5  反演参数随迭代的变化曲线

Figure  5.  Evolution of parameters with iteration

下载: 全尺寸图片 幻灯片

图6为剪切波速剖面及传递函数的估计值和初始值与真实值的对比,估计结果和真实结果几乎完全一样。图7为分别用剪切波速的估计值、初始值和真实值计算得到的地表加速度响应,从图中可以看出,使用剪切波速的估计值能够得到非常准确的地表响应估计结果。根据该数值算例的反演结果,可以得到在没有模型误差的情况下,UKF算法能够准确地估计出土体动力模型参数,参数的估计值能够完全收敛于真实值,并且能够得到非常准确的地表响应估计结果。

图  6  剪切波速及传递函数估计结果

Figure  6.  Estimation of shear wave velocity and transfer function

下载: 全尺寸图片 幻灯片

图  7  地表加速度响应对比

Figure  7.  Comparison of time histories of ground surface acceleration

下载: 全尺寸图片 幻灯片

4.   剪切波速反演算法应用实例

4.1   剪切波速反演算例分析

图8为地震记录IWTH080508161146加速度时间序列,该地震发生于2005年8月16日,震源深度为42 km,震中距239.06 km,震级为7.2级,井上地震动PGA为0.0491g。井下地震动Arias强度5%～95%的时段为32～68 s,以时窗加速度序列作为反演目标${\text{y}}_k$,初始时窗$[t_1^1,t_2^1]=[28\text{ s},30\text{ s}]$,每次迭代时窗移动0.2 s。

图  8  地震记录IWTH080508161146的加速度时间序列

Figure  8.  Time histories of acceleration of event IWTH080508161146

下载: 全尺寸图片 幻灯片

尽管KiK-net提供了测站的剪切波速剖面,但这些波速剖面过于粗略且不够精准,无法准确表征场地的放大效应。以KiK-net提供的剪切波速数据作为反演初始值,瑞利阻尼系数$a_1$初始值假定为0.001,则${\widehat{\theta }}_0^{+}=[150,280,280,280,280,680,680,900,2120,0.001]$,参数的初始协方差$P_0^{+}=\text{diag{(}0.15{\widehat{\theta }}_0^{+}{)}^2\}$,$P_0^{+}$越大,意味着对初始参数的置信度越小,太大的$P_0^{+}$有可能造成算法前期迭代不稳定。过程噪声协方差$Q=$$\text{diag}\left\{{({10}^{-4}{\widehat{\theta }}_0^{+})}^2\right\}$。与第2节不同,这里使用实测地表响应作为反演的匹配目标。由于地表响应计算模型为理想的一维模型,与实际情况存在较大的误差,同时,地震动的记录也存在噪声。用量测误差的协方差矩阵$R$来表征量测误差大小,$R$越大,意味着量测误差越大,参数更新时对匹配目标的估计误差更不敏感,这有可能使得估计结果不准确；$R$越小,意味着量测误差越小,参数更新时对匹配目标的估计误差更敏感,这有可能使得算法前期迭代不稳定。经过大量的试验计算,量测误差取地表加速度时间序列PGA的3%～10%能得到相对较好的反演结果,在该算例中取$R={(0.004g)}^{2}I$。

图9为反演参数随迭代的变化曲线,在前50次迭代,参数变化比较剧烈,随后参数慢慢收敛,标准差也逐渐减小。图10为该算例反演得到的剪切波速结构以及根据该波速结构计算得到的传递函数与实测结果的对比,图11为地表加速度响应的估计结果和实测结果的对比。

图  9  反演参数随迭代的变化曲线

Figure  9.  Evolution of parameters with iteration

下载: 全尺寸图片 幻灯片

图  10  剪切波速及传递函数估计结果

Figure  10.  Estimation of shear wave velocity and transfer function

下载: 全尺寸图片 幻灯片

图  11  地表加速度响应估计与实测的对比

Figure  11.  Comparison of time histories of ground surface acceleration

下载: 全尺寸图片 幻灯片

从加速度时程上看,反演得到的估计结果在每个相位上都和实测结果非常接近,峰值也基本吻合。反演得到传递函数估计结果也和实测结果较接近,在前3个主频上和实测结果基本吻合。这些结果说明,反演得到的剪切波速结构能够较好地反映实际的土层性质,运用该波速结构能较好地预测出在该地震激励下的地表响应。

4.2   多条地震记录反演结果

在IWTH08测站中选取10条以上井上地震动PGA在0.1g以下的地震记录进行剪切波速的反演计算,得到的剪切波速剖面反演结果及传递函数估计结果如图12所示。剪切波速剖面中,黑粗线为反演初始值,每一条灰线代表每条地震记录反演得到的波速结构。对这些剪切波速结构求平均,可以得到每个测站的剪切波速反演结果（图中红线所示）。传递函数中,红线为根据剪切波速反演结果得到的理论传递函数,每条灰线分别为根据每个地震记录井上井下地震动数据得到的传递函数,黑线为所有灰线的平均值。

图  12  多条地震记录剪切波速剖面和传递函数估计结果

Figure  12.  Estimation of shear wave velocity profile and transfer function using different ground motions

下载: 全尺寸图片 幻灯片

从图中可以看到,传递函数的估计结果和实测结果比较接近,每个主频频率基本上能够匹配准确,说明反演得到的剪切波速结构能够基本表征土层的线性性质。由于场地存在一定的空间变异性以及地震波传播的复杂性,实测传递函数呈现带状分布,不同地震记录的传递函数存在一定的差异性,这也使得不同地震记录反演得到的剪切波速结构存在较大的差异性,如图12所示。此外,反演也可能存在一定程度上的多解性。增大用于反演的地震记录数据量,或者采用一种多条地震记录联合反演的方法,能够在一定程度上减小多解性的影响,得到更接近实际情况的剪切波速结构。

5.   结论

本文探索了一种基于无迹卡尔曼滤波的贝叶斯估计方法,使用日本KiK-net台站数据对测站的土体动力响应模型参数进行反演。首先使用一个数值算例验证了该算法的适用性,然后使用IWTH08测站实测地震动数据进行了剪切波速反演计算。

（1）基于无迹卡尔曼滤波的贝叶斯反演算法在土体动力响应模型的适用性较好,在没有模型误差的情况下,参数的估计值能完全收敛于真实值,得到非常准确的地表响应估计结果。

（2）使用实测地震数据能够反演得到更精细的剪切波速结构,运用反演得到的剪切波速结构能够较好地表征土层对地震动的放大效应,得到较为准确的地表响应估计结果。

本文基于KiK-net地震记录数据,运用贝叶斯反演方法对剪切波速结构进行了反演计算研究。理论上来说,该算法还适用于其他参数及其他土体动力响应模型的反演计算。KiK-net还记录了很多大震（井上地震动PGA>0.1g）数据,运用这些大震数据,结合本文的反演算法,还可对土体的非线性特性进行反演分析,获得土体在大变形下的性质特征,这些还有待未来进一步的研究。