0.   引言

参数确定是岩土工程设计的基础。土体性质不确定性显著，工程实践中受经济和技术等因素制约，单个场地的土工测试数量有限，基于有限数据确定设计参数存在可靠性低的问题^[1-2]。以海上风电大直径单桩工程为例，砂土内摩擦角$\varphi$和黏土不排水抗剪强度${s_{\text{u}}}$是水平承载力设计p-y曲线法的关键参数^[3]。由于海洋工程勘察取样和测试成本高，数据量有限，一般根据有限的实测数据结合工程经验确定设计值，造成设计偏保守。利用有限数据对土体强度参数进行科学估计，降低大直径单桩设计参数的不确定性，实现设计优化，是当前海上风电工程提质增效的迫切需要。

目前，国内外在岩土体参数估计方面已开展了大量研究，提出了一些有效的方法，具体可分为3类。第一类方法将单个岩土体参数假定为随机变量，采用参数估计或非参数估计确定其概率分布^[4-7]。然而，这类方法在单个场地、有限小样本数据条件下可靠度往往较低，且仅能利用待估计岩土体参数的测试数据，忽略了不同类型岩土体参数之间的相关性，因此不能利用不同类型的数据。第二类方法通过经验公式^[8-12]将其它岩土体参数的实测值转换为待估参数。但经验公式一般通过全球或地区范围内实测数据拟合得到，具有显著的不确定性^[9]，不同经验公式的估计差异较大^[13]。第三类方法将多种类型的岩土体参数视为多元随机变量，以多元随机变量概率分布表示参数之间的相互关系，利用实测数据对概率分布的模型参数进行估计，代表性的工作有，Ching等^[14]针对9种岩体参数建立了岩体参数的多维联合分布；李典庆等^[1]和Tang等^[15]基于多维Gaussian Copula建立了多维土体参数的联合分布；汪海林等^[16]针对黏土的9个土性参数，基于贝叶斯理论对土性参数的联合后验分布进行了研究。上述研究主要利用单个场地的数据，当数据量有限时，结果存在较大的不确定性，为此一些学者通过融合不同场地的数据来降低不确定性。张广文等^[2]利用多场地数据建立黏聚力的先验分布，再利用某场地的实测数据得到后验分布；Bozorgzadeh等^[17]和Xiao等^[18]利用多场地数据分别对岩石力学特性参数和剪切波速参数进行了参数估计；Ching等^[19-20]提出了多元土体参数的混合贝叶斯方法，可综合利用场地数据；Ching等^[21]采用层次贝叶斯方法和混合贝叶斯方法对不排水抗剪强度比${s_{\text{u}}}/{\sigma '_{\text{v}}}$和归一化锥尖阻力${\text{(}}{q_{\text{t}}} - {\sigma _{\text{v}}}{\text{)}}/{\sigma '_{\text{v}}}$进行了二元正态分布的参数估计。整体而言，当前针对黏土不排水抗剪强度的研究主要集中在前两种方法，采用第三类方法的还较为有限^[21]，也较少涉及多场地数据的融合，而且现有多场地的研究主要考虑场地间的相似性，忽略了场地间的差异性。

针对上述不足，本文将采用基于单个场地数据进行多维参数估计的特定场地贝叶斯方法（Site-specific Bayesian Method，SBM）和可利用多场地数据进行多维参数估计的层次贝叶斯方法（Hierarchical Bayesian Method，HBM）对不排水抗剪强度的概率分布进行估计，并对这两种方法的准确性和适用性进行定量评价，为海上风电桩基工程设计优化和参数确定提供理论依据。

1.   多维土体参数概率分布的贝叶斯估计方法

1.1   贝叶斯参数估计的基本原理

岩土工程勘察往往涉及多个土体参数^[1]，如重度、含水率、稠度指标、不排水抗剪强度等。假定土体参数$\boldsymbol{X}=\left(X_1, X_2, \cdots, X_d\right)^{\mathrm{T}}$为d维随机变量，其概率密度函数的模型参数设为$\boldsymbol{\theta}$。例如，当X服从多元正态分布时，记$\boldsymbol{X} \sim {\text{N}}(\boldsymbol{\mu }, \boldsymbol{C})$，其中$\boldsymbol{\mu }$为均值列向量，$\boldsymbol{C}$为协方差矩阵，则模型参数$\boldsymbol{\theta} = (\boldsymbol{\mu }, \boldsymbol{C})$。

由贝叶斯理论，可由X的观测数据得到$\theta$的后验分布：

式中：Data为X的观测数据；$f(\left. {\bf{ Data}} \right|\mathit{\boldsymbol{\theta }} )$为似然函数，即给定$\theta$条件下得到Data发生的概率；$f(\mathit{\boldsymbol{\theta }} )$为$\theta$的先验分布；$f({\bf{ Data}})$是使后验概率总和为1的归一化常数。

需要指出的是，式（1）为通式。在不同数据条件下，需采用不同方法得到具体形式。特定场地贝叶斯方法利用本场地的数据对模型参数进行概率估计。层次贝叶斯方法通过构建多层级的贝叶斯模型，考虑不同场地的数据之间的相似性和差异性，利用多个场地的数据进行模型参数的概率估计。下面介绍特定场地贝叶斯方法和层次贝叶斯方法的基本原理。

1.2   特定场地贝叶斯方法

如图 1所示，假定土体参数X为d维随机变量，服从多元正态分布${\rm{N}}\left( {{\mathit{\boldsymbol{\mu }}_{\bf{s}}}, {\mathit{\boldsymbol{C}}_{\bf{s}}}} \right)$，模型参数$\mathit{\boldsymbol{\theta }} = {\mathit{\boldsymbol{\theta }}_{\bf{s}}} = \left( {{\mathit{\boldsymbol{\mu }}_{\bf{s}}}, {\mathit{\boldsymbol{C}}_{\bf{s}}}} \right)$，下标s表示SBM。当在某一测点得到数据$\boldsymbol{x}=\left(x_1, x_2, \cdots, x_d\right)^{\mathrm{T}}$时，式（1）中$f({\rm{ }}{\bf{Data}}{\rm{ }}\mid \mathit{\boldsymbol{\theta }})$可表示为

图  1  特定场地贝叶斯方法示意图

Figure  1.  Diagram of site-specific Bayesian method

下载: 全尺寸图片 幻灯片

式中：${\rm{N}}\left( {\mathit{\boldsymbol{x}}\mid {\mathit{\boldsymbol{\mu }}_{\bf{s}}}, {\mathit{\boldsymbol{C}}_{\bf{s}}}} \right)$为多元正态分布；${( \cdot )^T}$表示转置。

如在m个测点获得观测数据${\mathit{\boldsymbol{x}}_{\bf{s}}} = \left( {{\mathit{\boldsymbol{x}}_1}, {\mathit{\boldsymbol{x}}_2}, \cdots , {\mathit{\boldsymbol{x}}_m}} \right)$时，${\mathit{\boldsymbol{x}}_{\bf{s}}}$为$d \times m$矩阵。假设各测点土体参数独立同分布，由式（2）可得$f(\left. {{\bf{Data}}} \right|\mathit{\boldsymbol{\theta }})$为

假定${\mathit{\boldsymbol{\mu }}_{\bf{s}}}$和${\mathit{\boldsymbol{C}}_{\bf{s}}}$相互独立，则由式（1）可得${\mathit{\boldsymbol{\theta }}_{\bf{s}}}$的后验分布为

式中：$f(\left. {{\mathit{\boldsymbol{\mu }}_{\bf{s}}}, {\mathit{\boldsymbol{C}}_{\bf{s}}}} \right|{\mathit{\boldsymbol{x}}_{\bf{s}}})$为${\mathit{\boldsymbol{\mu }}_{\bf{s}}}$和${\mathit{\boldsymbol{C}}_{\bf{s}}}$的联合后验分布；$f({\mathit{\boldsymbol{\mu }}_{\bf{s}}})$和$f({\mathit{\boldsymbol{C}}_{\bf{s}}})$分别为${\mathit{\boldsymbol{\mu }}_{\bf{s}}}$和${\mathit{\boldsymbol{C}}_{\bf{s}}}$的先验分布。这里分别取二者的共轭先验，使后验分布与先验分布属于同一分布族。其中${\mathit{\boldsymbol{\mu }}_{\bf{s}}}$的共轭先验为多元正态分布${\rm{N}}({\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}, {\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}})$，${\mathit{\boldsymbol{C}}_{\bf{s}}}$的共轭先验为逆威希特分布^[22]，记为${\rm{IW}}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\bf{C}}}, {{\nu }_{\bf{C}}}} \right)$，则${\mathit{\boldsymbol{\mu }}_{\bf{s}}}$和${\mathit{\boldsymbol{C}}_{\bf{s}}}$的先验概率密度函数分别为

式中：${\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}$，${\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}}$分别为${\mathit{\boldsymbol{\mu }}_{\bf{s}}}$的均值列向量和协方差矩阵；${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}}$，${\nu _{\text{C}}}$分别为逆威希特分布的尺度矩阵和自由度；${\mathit{\Gamma }_{\rm{d}}}( \cdot )$表示自由度为d的多元伽马分布；${\text{tr(}} \cdot {\text{)}}$为矩阵的迹。一般称$({\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}, {\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}})$和$({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\rm{C}}}, {\nu _{\rm{C}}})$为超参数。

对式（4）中${\mathit{\boldsymbol{\theta }}_{\bf{s}}}$的联合后验分布积分可得${\mathit{\boldsymbol{\mu }}_{\bf{s}}}$和${\mathit{\boldsymbol{C}}_{\bf{s}}}$的边缘后验分布分别为

1.3   层次贝叶斯方法

如图 2所示，假设共有r个场地。其中第i个场地的土体参数为多维随机变量${\mathit{\boldsymbol{X}}_i}$，${\mathit{\boldsymbol{X}}_i}$服从多元正态分布，记为${\mathit{\boldsymbol{X}}_i} \sim {\text{N}}({\mathit{\boldsymbol{\mu }}_i}, {\mathit{\boldsymbol{C}}_i})$，$i = 1, 2, \cdots , r$。${\mathit{\boldsymbol{\mu }}_i}$和${\mathit{\boldsymbol{C}}_i}$分别为${\mathit{\boldsymbol{X}}_i}$的均值列向量和协方差矩阵，记${\mathit{\boldsymbol{\theta }}_{\bf{h}}} = {\text{(}}{\mathit{\boldsymbol{\mu }}_{\bf{h}}}, {\mathit{\boldsymbol{C}}_{\bf{h}}}{\text{)}}$，其中${\mathit{\boldsymbol{\mu }}_{\bf{h}}} = ({\mathit{\boldsymbol{\mu }}_{\text{1}}}, {\mathit{\boldsymbol{\mu }}_{\text{2}}}, \cdots , {\mathit{\boldsymbol{\mu }}_r})$，${\mathit{\boldsymbol{C}}_{\bf{h}}} = {\text{(}}{\mathit{\boldsymbol{C}}_{\text{1}}}, {\mathit{\boldsymbol{C}}_{\text{2}}}, \cdots , {\mathit{\boldsymbol{C}}_r}{\text{)}}$，下标h表示HBM。假设各场地的${\mathit{\boldsymbol{\mu }}_i}$和${\mathit{\boldsymbol{C}}_i}$相互独立，但具有相同的统计特性，其中${\mathit{\boldsymbol{\mu }}_i}$服从多元正态分布${\text{N}}({\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}, {\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}})$，${\mathit{\boldsymbol{C}}_i}$服从逆威希特分布${\text{IW}}({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_C}, {\nu _C})$。假定超参数${\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}$，${\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}}$和${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{C}}$为随机变量，超参数${\nu _{C}}$为常数，记超参数随机变量${\mathit{\boldsymbol{\phi}} _{\;{\bf{h}}}} = ({\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}, {\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}}, {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\rm{C}}})$。

图  2  层次贝叶斯方法示意图

Figure  2.  Diagram of hierarchical Bayesian method

下载: 全尺寸图片 幻灯片

HBM方法中的模型参数$\mathit{\boldsymbol{\theta }}$既包含${\mathit{\boldsymbol{\theta }}_{\bf{h}}}$，也包含${\mathit{\boldsymbol{\phi}} _{\;{\bf{h}}}}$。由式（1），模型参数的联合后验分布可改写为

式中：$f(\left. {{\bf{Data}}} \right|{\mathit{\boldsymbol{\theta }}_{\bf{h}}})$为似然函数；$f({\mathit{\boldsymbol{\theta }}_{\bf{h}}}, {\mathit{\boldsymbol{\phi}} _{\;{\bf{h}}}})$为联合先验分布，

式中：$f(\left. {{\mathit{\boldsymbol{\theta }}_{\bf{h}}}} \right|{\mathit{\boldsymbol{\phi}} _{\;{\bf{h}}}})$为${\mathit{\boldsymbol{\theta }}_{\bf{h}}}$的先验分布；$f({\mathit{\boldsymbol{\phi}} _{\;{\bf{h}}}})$为超参数的先验分布。式（9）中第二行等式成立是因为与${\bf{Data}}$直接相关的参数只有${\mathit{\boldsymbol{\theta }}_{\bf{h}}} = ({\mathit{\boldsymbol{\mu }}_{\bf{h}}}, {\mathit{\boldsymbol{C}}_{\bf{h}}})$，而${\mathit{\boldsymbol{\phi}} _{\;{\bf{h}}}}$只能通过${\mathit{\boldsymbol{\theta }}_{\bf{h}}}$间接地影响${\bf{Data}}$^[22]。由此体现了HBM方法的层级特性。

当第i个场地有${m_i}$个测点，且各测点都有完整的d维数据时，则${\bf{Data}}$可表示为${\mathit{\boldsymbol{x}}_{\bf{h}}} = ({\mathit{\boldsymbol{x}}_1}, {\mathit{\boldsymbol{x}}_2}, \cdots , {\mathit{\boldsymbol{x}}_r})$，为$d \times \left( {\sum\nolimits_{i = 1}^r {{m_i}} } \right)$矩阵，其中${\mathit{\boldsymbol{x}}_i} = ({\mathit{\boldsymbol{x}}_{i1}}, {\mathit{\boldsymbol{x}}_{i2}}, \cdots , {\mathit{\boldsymbol{x}}_{i{m_i}}})$，为$d \times {m_i}$维矩阵，$i = 1, 2, \cdots , r$。此时，式（9）中$f(\left. {{\bf{Data}}} \right|{\mathit{\boldsymbol{\theta }}_{\bf{h}}})$对应为

与式（4）的推导类似，假定${\mathit{\boldsymbol{\mu }}_{\bf{h}}}$，${\mathit{\boldsymbol{C}}_{\bf{h}}}$相互独立，且超参数${\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}$，${\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}}$，${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}}$相互独立，则可得${\mathit{\boldsymbol{\theta }}_{\bf{h}}}$，${\mathit{\boldsymbol{\phi}} _{\;{\bf{h}}}}$的联合后验分布为

式中：$f(\left. {{\mathit{\boldsymbol{\mu }}_{\bf{h}}}} \right|{\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}, {\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}})$为${\mathit{\boldsymbol{\mu }}_{\bf{h}}}$的先验分布，${\mathit{\boldsymbol{\mu }}_i} \sim$ ${\text{N}}({\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}, {\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}})$；$f(\left. {{\mathit{\boldsymbol{C}}_{\bf{h}}}} \right|{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}}, {\nu _{\text{C}}})$为${\mathit{\boldsymbol{C}}_{\bf{h}}}$的先验分布，${\mathit{\boldsymbol{C}}_i} \sim$ ${\text{IW}}({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}}, {\nu _{\text{C}}})$；$f{\text{(}}{\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}{\text{)}}$，$f({\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}})$，$f({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}})$分别为超参数${\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}$，${\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}}$，${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}}$的先验分布。这里取${\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}$和${\mathit{\boldsymbol{C}}_{\mu }}$的共轭先验分别为${\text{N}}({\mathit{\boldsymbol{\mu }}_{\bf{0}}}, {\mathit{\boldsymbol{C}}_{\bf{0}}})$，${\text{IW}}({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\bf{0}}}, {\nu _{\bf{0}}})$，${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}}$的共轭先验为威希特分布^[22]，记为$W({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\Sigma }}, {\nu _{\Sigma }})$，则${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}}$的先验分布为

式中：${\text{W}}(\left. {{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}}} \right|{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\mathit{\Sigma }}, {\nu _\mathit{\Sigma }})$为威希特分布，其中${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\mathit{\Sigma }}$，${\nu _\mathit{\Sigma }}$分别为尺度矩阵和自由度。

需要指出的是，式（12）所得为模型参数${\mathit{\boldsymbol{\theta }}_{\bf{h}}}$，${\mathit{\boldsymbol{\phi}} _{\;{\bf{h}}}}$的联合后验分布，各边缘后验分布可通过类似式（7），（8）的积分得到。

1.4   土体参数的估计

上述两种方法除了可以得到各模型参数的后验分布，还可对土体参数本身进行估计。

（1）测点缺失数据的估计

由于条件限制，在某一测点处不一定能获得土体参数$\mathit{\boldsymbol{X}} = {({X_1}, {X_2}, \cdots , {X_d})^T}$的所有d维实测数据^{[9, 21]}。例如，对某测点的土样只开展了含水率等基本物性试验，但并未开展不排水抗剪强度试验，则该测点缺失${s_{\text{u}}}$的实测数据。此时可利用已知的数据，对缺失数据的土体参数进行估计。

假设在某测点d个土体参数中前k个有实测值，记为$\mathit{\boldsymbol{x}}_k^{\text{o}}$，对应式（1）中的${\bf{Data}}$；后$d - k$个参数无实测数据，记为$\mathit{\boldsymbol{x}}_{d - k}^{\text{u}}$，是待估计的参数。下文简化为实测值${\mathit{\boldsymbol{x}}^{\text{o}}}$和缺失项${\mathit{\boldsymbol{x}}^{\text{u}}}$。

基于观测数据得到模型参数$\mathit{\boldsymbol{\theta }}$的后验分布$f(\left. \mathit{\boldsymbol{\theta }} \right|{\bf{Data}})$后，通过积分可得${\mathit{\boldsymbol{x}}^{\text{u}}}$的后验分布为

式中：$f(\left. \mathit{\boldsymbol{\theta }} \right|{\bf{Data}})$由式（4）或（12）得到$f(\left. {{\mathit{\boldsymbol{x}}^{\text{u}}}} \right|\mathit{\boldsymbol{\theta }})$为${\mathit{\boldsymbol{x}}^{\text{u}}}$的似然函数：

式中：$\boldsymbol{\mu}$分块记为$\boldsymbol{\mu} = {({\boldsymbol{\mu}^{\text{o}}}, {\boldsymbol{\mu}^{\text{u}}})^T}$，$\boldsymbol{C}$分块记为

式中：${\boldsymbol{\mu }^{\text{o}}}$为k维列向量，${\boldsymbol{\mu }^{\text{u}}}$为$d - k$维列向量；${\boldsymbol{C}^{\text{o}}}$为$k \times k$矩阵，${\boldsymbol{C}^{{\text{ou}}}}$为$k \times (d - k)$矩阵，${\boldsymbol{C}^{{\text{uo}}}}$为$(d - k) \times k$维矩阵，${\boldsymbol{C}^{{\text{ou}}}} = {({\boldsymbol{C}^{{\text{uo}}}})^{\text{T}}}$，${\boldsymbol{C}^{\text{u}}}$为$(d - k) \times (d - k)$矩阵。

（2）未知点处土体参数的估计

工程中有时还需要对未知点位（无任何实测数据的位置）的土体参数${\boldsymbol{x}^{{\text{new}}}}$进行估计。类似式（14），对模型参数进行积分可得${\boldsymbol{x}^{{\text{new}}}}$的后验分布为

式中：$f(\left. {{\boldsymbol{x}^{{\text{new}}}}} \right|\boldsymbol{\theta })$为似然函数，即多元正态分布。

式（7），（8），（14），（17）的积分可通过马尔科夫链-蒙特卡洛模拟方法（Markov Chain Monte Carlo，MCMC）进行求解。MCMC的核心思想是构造一条平稳分布为待估参数后验分布的马尔科夫链，再利用马尔科夫链达到平稳分布时的样本进行蒙特卡洛积分^[22]。为提高计算效率，本文采用MCMC的吉布斯抽样算法，具体可见文献[22]。

2.   工程应用：不排水抗剪强度的参数估计

不排水抗剪强度${s_{\text{u}}}$是风电大直径单桩设计的重要参数。由于含水率等物理性质试验较强度试验更易开展，实践中常利用液性指数I_L来估计${s_{\text{u}}}$^[13]。本文以${X_1} = {\text{LI}}$和${X_2} = \ln ({s_{\text{u}}}/{\sigma '_{\text{v}}})$为正态随机变量，假定$\boldsymbol{X} = {({X_1}, {X_2})^T}$服从二元正态分布，利用LI对$\ln ({s_{\text{u}}}/{\sigma '_{\text{v}}})$进行估计。

2.1   数据集的建立

从CLAY/10/7490数据库^{[20, 23]}中提取141个场地共1679组LI和$\ln ({s_{\text{u}}}/{\sigma '_{\text{v}}})$数据建立数据集，各场地的数据量至少5组。液性指数I_L和竖向有效应力${\sigma '_{\text{v}}}$分别由室内试验和现场孔压静力触探测定，${s_{\text{u}}}$为不排水抗剪强度修正值，由各类不排水剪切试验的实测值转换得到，转换公式详见文献[23]。原数据库^{[20, 23]}中部分场地仅有编号，无具体位置信息，故本文也仅用编号代表各场地。本文数据集中有13个场地含20组及以上数据，数据量较大，故将这13个场地作为目标场地进行参数估计，分别标记为场地1~场地13。如图 3所示为各场地土体的塑性图，其中高液限黏土CH和低液限黏土CL分别占49.1%，34.7%，高液限粉土MH和低液限粉土ML分别占13.1%，2.6%。

图  3  各场地土体的类型

注：数据集中少量LL或PI数据缺失。

Figure  3.  Plasticity chart

下载: 全尺寸图片 幻灯片

经Shapiro-Wilk检验^[24]，I_L和$\ln ({s_{\text{u}}}/{\sigma '_{\text{v}}})$的p值分别为1.33和4.85，均大于0.05，即在置信度为95%的条件下，I_L和$\ln ({s_{\text{u}}}/{\sigma '_{\text{v}}})$的边缘分布均满足正态性假设。采用表 1所示4种Copula函数对数据进行拟合，得到Copula参数${\theta _{\text{C}}}$，再利用AIC指标和BIC指标比较各Copula函数的拟合优度^{[15, 25]}，最终确定数据的最优拟合Copula函数为Gaussian Copula函数。因此，本文数据满足二元正态分布假设。

表  1  Copula函数拟合

Table  1.  Copula function fitting

Copula	$C({u_1}, {u_2})$	${\theta _{\text{C}}}$	AIC	BIC
Gaussian	${\mathit{\Phi} _2}\left[ {{\mathit{\Phi} ^{ - 1}}({u_1}), {\mathit{\Phi} ^{ - 1}}({u_2});{\theta _{\text{C}}}} \right]$	-0.05	-3.50	1.93
Clayton	${(u_1^{ - {\theta _{\text{C}}}} + u_2^{ - {\theta _{\text{C}}}} - 1)^{ - 1/{\theta _{\text{C}}}}}$	1.45×10-6	2.00	7.43
Frank	$- \frac{1}{{{\theta _{\text{C}}}}}\ln \left[ {1 + \frac{{({{\text{e}}^{ - {\theta _{\text{C}}}{u_1}}} - 1)({{\text{e}}^{ - {\theta _{\text{C}}}{u_2}}} - 1)}}{{{{\text{e}}^{ - {\theta _{\text{C}}}}} - 1}}} \right]$	-0.07	1.80	7.22
Gumbel	$\exp \left\{ { - {{\left[ {{{( - \ln {u_1})}^{{\theta _{\text{C}}}}} + {{( - \ln {u_2})}^{{\theta _{\text{C}}}}}} \right]}^{^{1/{\theta _{\text{C}}}}}}} \right\}$	1.00	2.00	7.43
注：${u_1}$和${u_2}$分别为变量${X_1}$和${X_2}$的边缘累积分布；${\mathit{\Phi} _2}$为标准二元正态分布的累积分布函数；${\mathit{\Phi} ^{ - 1}}( \cdot )$为标准正态分布累积分布的逆函数；加粗值为AIC和BIC值的最小值，对应最优拟合Copula函数。

下载: 导出CSV 

| 显示表格

各目标场地的实测数据如图 4所示，${X_1} = {I_{\text{L}}}$的分布区间为$\left[ { - {\text{0}}{\text{.44}}, 3.25} \right]$，${X_2} = \ln ({s_{\text{u}}}/{\sigma '_{\text{v}}})$的分布区间为$\left[ { - 3.68, 0.91} \right]$。场地的均值$\overline {{x_1}}$、$\overline {{x_2}}$和协方差${S_{11}}$、${S_{12}}( = {S_{21}})$、${S_{22}}$汇总于表 2。

图  4  目标场地实测数据的散点图

Figure  4.  Scatter plots of measured data from target sites

下载: 全尺寸图片 幻灯片

表  2  目标场地的样本统计量

Table  2.  Sample statistics of measured data

编号	${m_i}$	${\bar x_1}$	${\bar x_2}$	${S_{11}}$	${S_{12}}$	${S_{22}}$
1	22	0.131	-0.632	0.015	0.027	0.163
2	39	0.839	-0.666	0.603	-0.021	1.222
3	30	0.861	-0.761	0.637	-0.005	0.806
4	32	0.932	-1.752	0.648	-0.008	0.031
5	20	0.526	-1.176	0.036	0.011	0.028
6	21	0.728	-1.306	0.139	0.015	0.065
7	20	0.598	-1.067	0.032	-0.072	0.464
8	25	0.847	-1.256	0.007	0.006	0.265
9	21	1.236	-1.213	0.093	-0.017	0.066
10	34	1.019	-1.096	0.055	0.091	0.268
11	22	0.742	-1.969	0.037	0.002	0.105
12	20	0.827	-2.253	0.092	0.009	0.120
13	22	2.368	-0.526	0.246	-0.026	0.018
注：${m_i}$为第i个场地的数据量（测点数量）。

下载: 导出CSV 

| 显示表格

2.2   先验分布的参数确定

先验分布是贝叶斯估计的重要组成部分，前文已设定各先验均为共轭先验，本节确定先验分布中的参数。Cao等^[26]基于文献总结了LI均值和标准差的先验范围分别为0.5~2.5，0.025~2.2，${s_{\text{u}}}/{\sigma '_{\text{v}}}$均值和标准差的先验范围分别为0.23~1.4，0.01~1.26。为降低先验分布对后验分布影响，本文参考Cao等^[26]并采用更为平坦的先验分布：假定${X_1} = {I_{\text{L}}}$和${X_2} = \ln ({s_{\text{u}}}/{\sigma '_{\text{v}}})$的均值都为0，标准差都为5。

SBM方法中需要确定${\mathit{\boldsymbol{\mu }}_{\bf{s}}}$和${\mathit{\boldsymbol{C}}_{\bf{s}}}$的先验分布${\text{N}}({\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}, {\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}})$，${\text{IW}}({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}}, {\nu _{\text{C}}})$中的参数取值。根据上述假定和相关研究^{[22, 27]}，确定如表 3所示的参数值。HBM方法中需要确定超参数${\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}$，${\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}}$，${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}}$的先验分布${\text{N}}({\mathit{\boldsymbol{\mu }}_{\bf{0}}}, {\mathit{\boldsymbol{C}}_{\bf{0}}})$，${\text{IW}}({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\bf{0}}}, {\nu _{\bf{0}}})$，${\text{W}}({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\mathit{\Sigma }}, {\nu _\mathit{\Sigma }})$中的参数，以及超参数${\nu _{\text{C}}}$。这里假定${\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}$的标准差为${\mathit{\boldsymbol{\mu }}_i}$标准差的一半，结合威希特分布期望的性质，确定先验分布参数的取值，汇总在表 3。

表  3  先验分布的参数取值

Table  3.  Parameter values of prior distributions

SBM					HBM
${\mathit{\boldsymbol{\mu }}_{\mu }}$	${\mathit{\boldsymbol{C}}_{\mu }}$	$\mathit{\boldsymbol{{ \boldsymbol{\varSigma} }_{\text{C}}}}$	${\nu _{\text{C}}}$		${\mathit{\boldsymbol{\mu }}_{\bf{0}}}$	${\mathit{\boldsymbol{C}}_{\bf{0}}}$	${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\bf{0}}}$	${\nu _{\bf{0}}}$	${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_\mathit{\Sigma }}$	${\nu _\mathit{\Sigma }}$	${\nu _{\text{C}}}$
$\left( {\begin{array}{*{20}{c}} 0 \\ 0 \end{array}} \right)$	$\left( {\begin{array}{*{20}{c}} {25}&0 \\ 0&{25} \end{array}} \right)$	$\left( {\begin{array}{*{20}{c}} {25}&0 \\ 0&{25} \end{array}} \right)$	4		$\left( {\begin{array}{*{20}{c}} 0 \\ 0 \end{array}} \right)$	$\left( {\begin{array}{*{20}{c}} {6.25}&0 \\ 0&{6.25} \end{array}} \right)$	$\left( {\begin{array}{*{20}{c}} {25}&0 \\ 0&{25} \end{array}} \right)$	4	$\left( {\begin{array}{*{20}{c}} {6.25}&0 \\ 0&{6.25} \end{array}} \right)$	4	4

下载: 导出CSV 

| 显示表格

2.3   预测准确度的评价

采用留一交叉验证法（Leave-One-Out Cross- Validation，LOO-CV）对两种方法的预测准确度进行比较。该方法通过每次留出目标场地一个测点的${X_2}$实测值作为验证集，利用其他测点的实测数据和该测点的${X_1}$值对${X_2}$进行后验估计，循环直到场地内各测点都完成验证为止。

LOO-CV可结合点预测密度的对数值（log pointwise predictive density，lppd）^{[22, 28]}作为评价指标。假设目标场地是第i个场地，${m_i}$为该场地的数据量，则目标场地lppd_loo-cv指标的计算如下：

式中：${f_{{\text{post}}\left( {\backslash j} \right)}}( \cdot )$为基于训练集${{x}_{\backslash j}}$得到的后验概率密度函数，相应的验证集为$2 \times 1$维列向量${x_{ij}} = {({x_{ij, 1}}, {x_{ij, 2}})^{T} }$。SBM和HBM两种方法的训练集不同，SBM方法中$\boldsymbol{x}_{\backslash j}=\boldsymbol{x}_{\backslash j}^{\bf{s}}=\left(x_{i 1}, \ldots, x_{i, j-1}, x_{i, j+1}, \ldots, x_{i m_i}\right)$为$2 \times ({m_i} - 1)$维矩；HBM方法中$\boldsymbol{x}_{\backslash j}=\boldsymbol{x}_{\backslash j}^{\bf{h}}=\left(\boldsymbol{x}_1, \cdots, \boldsymbol{x}_{i \backslash j}, \cdots, \boldsymbol{x}_r\right)$为$2 \times \left( {\sum\nolimits_{i = 1}^r {{m_i}} - 1} \right)$维矩阵。其中r为场地总数；${\boldsymbol{x}_k}(k = 1, \cdots , i - 1, i + 1, \cdots , r)$为第k个场地所有测点数据的集合，为$2 \times {m_k}$维矩阵；${\boldsymbol{x}_{i\backslash j}} = ({x_{i1}}, \cdots , {x_{i, j - 1}}, {x_{i, j + 1}}, \cdots , {x_{i{m_i}}})$为第i个场地除验证集${x_{ij}}$以外测点数据的集合，为$2 \times ({m_i} - 1)$维矩阵。

3.   结果与讨论

3.1   数据缺失情况下的参数估计

为模拟实际工程中单个场地LI数据量大而$\ln ({s_{\text{u}}}/{\sigma '_{\text{v}}})$数据量少的情形，只取本文数据集中的部分数据作为贝叶斯参数估计的${\bf{Data}}$，利用SBM方法和HBM方法分别对$\ln ({s_{\text{u}}}/{\sigma '_{\text{v}}})$缺失数据$x_2^{\text{u}}$、未知点处的$x_2^{{\text{new}}}$和场地$\ln ({s_{\text{u}}}/{\sigma '_{\text{v}}})$均值${\mu _2}$进行后验估计。

以场地1为例，场地1共22个测点，即数据集中有22组$\boldsymbol{X} = {({X_1}, {X_2})^T}$实测数据。采用SBM方法时，取该场地22个${X_1}$数据和n个${X_2}$数据构成${\bf{Data}}$，其中${X_2}$数据对应n个${X_1}$分别（近似）为该场地${X_1}$数据的$n + 1$分位数。利用式（14），（17），（7）可分别得到$x_2^{\text{u}}$，$x_2^{{\text{new}}}$和${\mu _2}$的后验分布。采用HBM方法时，${\bf{Data}}$包括SBM方法中构造的本场地数据和数据集中其余140个场地的全部数据。

首先分析讨论给定数据量，即给定n值时，$x_2^{\text{u}}$的后验估计结果。这里以$n = 5$为例。如图 5（a）所示，场地9共有21个测点，图中黑色圆点表示Data中5个信息完备的测点，白色圆点表示16个假设${X_2}$值缺失的测点（待估点），线段表示$x_2^{\text{u}}$的95%后验置信区间（posterior confidence interval，简称PCI）。由图 5（a）可见，测点真实值均位于PCI中。SBM方法所得PCI的范围（红线）大，甚至超过了图 4中${X_2}$的实际分布区间$\left[ { - 3.68, 0.91} \right]$，表明SBM方法无法提供有效估计。HBM方法的PCI范围（蓝线）相比SBM的显著缩小，恰好覆盖真实值的范围，表明HBM方法可有效降低土体参数估计的不确定性。

图  5  $x_2^{\text{u}}$的95%后验置信区间

Figure  5.  95% posterior confidence interval of $x_2^{\text{u}}$

下载: 全尺寸图片 幻灯片

场地11和场地12的实测数据和PCI分别如图 5（b），5（c）所示，与上述结论基本一致，其余场地也得到了与图 5类似的结果，限于篇幅不再展示。各目标场地的估计中，MCMC抽样次数均为10000，由于抽样2000次后马尔科夫链即达到平稳，因此取后8000组样本作为后验分布的稳定样本，下文抽样过程与此一致。

下面讨论数据量n对场地未知点处$x_2^{{\text{new}}}$和${\mu _2}$后验估计的影响，并对两种方法进行比较。如图 6所示，x轴表示${\bf{Data}}$中有n组${X_2}$数据。以场地9（图 6（a））为例，当$n = 0$时，SBM方法的PCI范围很大，$x_2^{{\text{new}}}$的PCI范围为$\left[ { - 11.2, 13.4} \right]$，${\mu _2}$的PCI范围为$\left[ { - 9.1, 11.2} \right]$。这是因为SBM方法仅利用本场地的信息，当没有实测数据时，后验分布完全由先验分布决定。随着n的增大，PCI的范围大幅减小，表明SBM方法受数据量的影响显著。但即使${\bf{Data}}$使用了该场地的全部数据（$n = 21$），PCI的范围仍较大，$x_2^{{\text{new}}}$和${\mu _2}$的不确定性也较大。

图  6  数据量对参数估计的影响

Figure  6.  Influences of data volume on prediction of parameters

下载: 全尺寸图片 幻灯片

对于HBM方法，即使目标场地没有任何数据（$n = 0$），$x_2^{{\text{new}}}$和${\mu _2}$的PCI（$\left[ { - 2.8, 1.1} \right]$和$\left[ { - 2.5, 0.9} \right]$）已经明显小于SBM方法的，说明HBM参数估计的不确定性更低，这是因为HBM方法可利用其余140个场地的数据降低估计的不确定性。当数据量n增大到5之后，HBM所得PCI的范围已基本稳定，表明HBM方法的结果受目标场地数据量的影响较小。场地11和场地12的结果分别如图 6（b），6（c）所示，与上述结论一致。

3.2   不同场地后验估计和预测准确度的比较

下面分别利用13个目标场地的全部数据进行后验估计，得到各场地$x_2^{{\text{new}}}$和${\mu _2}$的PCI如图 7所示。结果表明，SBM方法仅利用了单个场地的数据，因此PCI的范围较大；而HBM方法可结合多个场地的数据进行估计，因此PCI的范围显著缩小，$x_2^{{\text{new}}}$和${\mu _2}$的不确定性大幅降低。

图  7  $x_2^{{\text{new}}}$和${\mu _2}$的95%后验置信区间

Figure  7.  95% posterior confidence intervals of $x_2^{{\text{new}}}$ and ${\mu _2}$

下载: 全尺寸图片 幻灯片

13个目标场地的lppd_loo-cv指标如表 4所示。lppd_loo-cv越大，表示预测准确度越高。由表 4可知，各场地中HBM方法对应的lppd_loo-cv基本都大于SBM方法的结果。其中场地2，3的数据离散、变异性较大（表 2），造成了两种方法的结果差异较小。将全部13个场地对应的评价指标求和，得到HBM方法的整体lppd_loo-cv指标为-76.7，SBM方法的整体lppd_loo-cv指标为-164.6，表明HBM方法的预测准确度更高，更适于有限数据下土体参数的概率估计。

表  4  两种方法的lppd_loo-cv评价指标比较

Table  4.  Comparison of lppd_loo-cv between SBM and HBM

场地编号	lppdloo-cv（SBM）	lppdloo-cv（HBM）
1	-10.6	-4.3
2	-26.9	-26.9
3	-18.7	-18.2
4	-11.6	3.4
5	-9.1	2.7
6	-9.6	-1.5
7	-11.5	-8.8
8	-12.4	-10.1
9	-9.6	-1.3
10	-15.4	-6.8
11	-10.2	-4.4
12	-9.7	-4.5
13	-9.5	4.2
整体求和	-164.6	-76.7

下载: 导出CSV 

| 显示表格

3.3   HBM方法统计量后验均值的工程应用

HBM方法所得13个目标场地统计量的后验均值汇总在表 5，与表 2对照可知，HBM方法所得统计量的后验均值与实测数据的统计量基本吻合。实际应用中，若需对这些目标场地中任一未知点的强度参数进行取值，则可直接根据${\mu _2}$和测点处的${\sigma '_{\text{v}}}$求得不排水抗剪强度。本文重要研究成果之一是利用HBM方法得到了超参数${\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}}$，${\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}}$，${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}}$的后验均值，分别为${\mathit{\boldsymbol{\mu }}_\mathit{\boldsymbol{\mu }}} = {\left( {1.040, - 0.931} \right)^T}$，${\mathit{\boldsymbol{C}}_\mathit{\boldsymbol{\mu }}} = \left( {\begin{array}{*{20}{c}} {0.539}&{ - 0.018} \\ { - 0.018}&{0.585} \end{array}} \right)$，${\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\text{C}}} = \left( {\begin{array}{*{20}{c}} {0.088}&{ - 0.008} \\ { - 0.008}&{0.313} \end{array}} \right)$。当需对一个新场地进行设计时，可直接利用上述结果确定新场地的不排水抗剪强度的设计值或者概率分布。

表  5  HBM所得场地统计量的后验均值

Table  5.  Posterior means of site statistics obtained by HBM

编号	${\mu _1}$	${\mu _2}$	${C_{11}}$	${C_{12}}$	${C_{22}}$
1	0.132	-0.635	0.018	0.026	0.170
2	0.843	-0.678	0.589	-0.021	1.198
3	0.871	-0.771	0.619	-0.007	0.790
4	0.936	-1.751	0.631	-0.008	0.040
5	0.529	-1.174	0.038	0.010	0.042
6	0.733	-1.304	0.137	0.014	0.077
7	0.598	-1.064	0.034	-0.069	0.457
8	0.847	-1.251	0.010	0.006	0.266
9	1.234	-1.211	0.093	-0.017	0.077
10	1.020	-1.093	0.056	0.087	0.269
11	0.743	-1.959	0.039	0.001	0.115
12	0.829	-2.238	0.092	0.008	0.130
13	2.342	-0.524	0.240	-0.025	0.032

下载: 导出CSV 

| 显示表格

4.   结论

针对场地有限数据条件下黏土不排水抗剪强度不确定性大的问题，建立了一个LI和$\ln ({s_{\text{u}}}/{\sigma '_{\text{v}}})$双变量参数数据集，应用特定场地贝叶斯SBM方法和层次贝叶斯HBM方法分别对13个目标场地进行参数估计，并对两种方法的估计结果进行了对比，得到2点结论。

（1）在测点缺失数据的情况下，SBM方法的PCI范围大，而HBM方法的PCI范围较小，刚好覆盖真实值的分布范围；随着数据量的增加，两种方法得到的PCI范围均有减小。SBM方法受测点数据量的影响更显著，且参数估计的不确定性较大。HBM方法受本场地数据量的影响较小，当数据量大于5时，PCI的范围就基本稳定，且参数估计的不确定性小。

（2）分别利用数据集中13个目标场地的全部实测数据对$\ln ({s_{\text{u}}}/{\sigma '_{\text{v}}})$进行估计，结果表明HBM方法得到$x_2^{{\text{new}}}$和${\mu _2}$的PCI范围均小于SBM方法对应的结果。采用留一交叉验证法对两种方法的预测准确度进行比较，得到HBM方法的整体lppd_loo-cv指标为-76.7，SBM方法的整体lppd_loo-cv指标为-164.6，表明HBM方法的预测准确度更高。

综上，HBM方法更适用于场地有限数据下土体参数的概率估计。此外，HBM方法得到的超参数${{\mu }_{\mu }}$，${{C}_{\mu }}$和${{\Sigma }_{\text{C}}}$的后验均值还可用于新场地的设计和可靠度计算。需要指出的是，上述结论依赖于本文数据集，在其他类型土体参数上的适用性有待后续研究。