基于分层贝叶斯学习的滨海软土地层高效识别方法

Efficient Identification Method of Coastal Soft Soil Stratum Based on Hierarchical Bayesian Learning

摘要

基于静力触探试验数据划分土层依赖于经验图表和主观判断，划分的土层剖面不可避免地存在不确定性.提出了一种基于土体分类指数I_c 的土层界面快速识别和不确定性量化方法.在分层贝叶斯学习框架下，所提方法采用全高斯概率模型表征土体空间变异性，通过引入正态‒逆威沙特共轭分布实现似然函数的快速计算，高效求解模型证据，识别最可能土层数目和厚度.所提方法基于I_c 的统计特性自动划分土层，提高了识别结果的可靠性.

关键词

静力触探试验 / 分层贝叶斯学习 / 土层划分 / 不确定性 / 岩土工程

Key words

cone penetration test / hierarchical Bayesian learning / underground stratification / uncertainty / geotechnical engineering

0 引言

我国滨海地区软土（包括软黏性土、淤泥等）分布广泛，由于沉积环境特点和地质成因的特殊性，土层厚度和工程性质空间变异性显著.软土强度低、压缩性大，若处理不当，会给工程建设带来极大的风险（谬林昌，2012；陈宇航，2020）.岩土工程勘察的一项重要任务是利用勘探资料划分土层，为岩土工程分析和设计提供必要的土层剖面信息（Zhao and Wang，2020）.静力触探试验（CPT）因其具有快速、经济、可重复等优点，同时可以提供接近连续的原位土体圆锥贯入响应（包括锥尖阻力、侧壁摩阻力等）测量值，在岩土工程勘察中广泛应用（Robertson，2009；Wang et al.，2019；Guan and Wang，2022）.

国内外学者在基于CPT的土层划分方面进行了大量的研究，提出了基于CPT数据的土体经验分类图表和土体分类指数I _c（如Robertson，1990， 2009；Jefferies and Davies，1993；张诚厚等，1997；Robertson and Wride，1998）.其中，张诚厚等（1997）根据荷兰等地区试验场地50余套孔压静力触探CPTU数据及其邻近钻孔信息，提出了修正后的锥尖阻力q_t 和超孔隙压力系数B_p 的土体分类图，但并未考虑侧壁摩阻力f_s；Robertson（1990）基于归一化的静力触探试验参数（包括归一化锥尖阻力Q_t 、归一化摩阻比F_r 和孔压参数比B_q ）提出了考虑力学行为的土体分类系统（Soil Behavior Type，SBT）；Jefferies and Davies（1993）和Robertson and Wride（1998）在已有土体SBT分类图的基础上，更新了I _c的计算公式及其相应的土体分类标准.此后，Robertson（2009）基于以往学者研究和试验观测提出了I _c中应力指数n，进一步修正了I _c的计算公式如下：

I_c = \sqrt[]{{\left(\mathrm{3.47}-\mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{10}}{Q}_{tn}\right)}^{\mathrm{2}}+{\left(\mathrm{l}\mathrm{o}{\mathrm{g}}_{\mathrm{10}}{F}_r+\mathrm{1.22}\right)}^{\mathrm{2}}},(1)

式中： Q_{tn} = \left(\left(q_t-{\sigma }_{v\mathrm{0}}\right)/P_a\right) \times {\left(P_a/\left({\sigma }_{v\mathrm{0}}^{\text{'}}\right)\right)}^n为考虑应力水平的归一化锥尖阻力； F_r = \left[f_s/\left(q_t-{\sigma }_{v\mathrm{0}}\right)\right] \mathrm{100} \mathrm{\%}为归一化摩阻比；q_t 为修正后的锥尖阻力； {\sigma }_{v\mathrm{0}}为竖向总应力；P_a 为标准大气压； {\sigma }_{v\mathrm{0}}^{\text{'}}为竖向有效应力；f_s 为侧壁摩阻力.对于净砂土和黏土而言，应力指数n分别等于0.5和1.0，其他土体类型的应力指数则可以通过公式 n = \mathrm{0.381} \left(I_c\right) + \mathrm{0.05} \left({\sigma }_{v\mathrm{0}}^{\text{'}}/P_a\right) - \mathrm{0.15}迭代计算.刘松玉等（2013）分析了江苏省不同典型地质成因的14个试验场地中95套CPTU及其邻近钻孔数据，表明Robertson and Wride（1998）提出的土体分类图表对我国土体分类较为适用，并基于土类指数提出了修正的土体分类表，使之与我国《岩土工程勘察规范（GB 50021-2001）》中土体分类一致，如表1所示.

基于I _c的土体分类图表（如表1）是半经验分类标准.由于数据有限以及其离散性，该分类标准在识别特定场地土体类型时不可避免地存在不确定性（Ku et al.，2010）.另一方面，由于CPT无法提供土样，只能根据每个测点的土体分类推断土层界面位置，每个测点土体分类的不确定性将对土层识别产生直接影响（Wang et al.，2013；Cao et al.，2019）.为合理考虑基于CPT经验图表土体分类不确定性对土层划分的影响，国内外学者进行了大量的研究，提出了模糊集分析（Zhang and Tumay，1999）、聚类分析（Hegazy and Mayne，2002；蔡国军等，2009）、贝叶斯系统识别（Wang et al.，2013）以及小波分析（Ching et al.，2015）等土层划分方法.尽管上述方法均能识别土层边界位置，但无法合理表征边界识别结果的不确定性.在工程实际中，相邻土层界面处的土体物理组成复杂，CPT锥尖阻力受锥尖位置处上覆土层和下卧土层的影响（Robertson and Campanella，1983），在密砂层该“影响区”范围甚至可达10~20倍锥径（Ahmadi and Robertson，2005）.此外，土体是一种天然材料，同一土层内不同位置的I _c必然具有离散性，这些因素都导致了土层界面位置存在不确定性.

近年来，不同机器学习模型如随机森林、贝叶斯学习、深度卷积网络等在岩土工程数据分析、反演、解译等方面得到广泛应用（Chen et al.， 2021；黄康等，2022；王德涛和陈国雄，2022；Wu et al.，2022）.为了量化土层边界的不确定性，Cao et al.（2019）基于贝叶斯学习，采用平稳随机场模型表征I _c的空间变异性，可识别最可能土层数目，并定量表征所识别土层边界的不确定性，但该方法需要假设不同土层具有相同的自相关函数，以量化不同测点之间的空间相关性，且计算效率较低.Zheng et al.（2019）通过优化算法，提出了基于CPT的快速贝叶斯土层划分方法，但该方法需要假定I_c 的自相关函数为单指数型.然而，对于实际土层，土层剖面划定前土体参数的空间相关结构未知，亟需提出一种更具普适性的高效土层划分方法，避免在土层划分前假设自相关函数类型.

本文提出了一种基于分层贝叶斯学习的高效土层划分方法，采用全高斯概率模型表征I _c的空间变异性，在无需假设空间相关结构的条件下高效识别土层界面并量化其不确定性.为此，本文将首先介绍所提分层贝叶斯学习框架及其计算流程，然后通过温州市仙湖调蓄工程试验区现场CPT数据和模拟I _c数据说明所提方法的有效性和合理性.

1 基于分层贝叶斯学习的土层划分方法

1.1　分层贝叶斯学习

令 \underline{\xi } = [ {\underline{\xi }}_{\mathrm{1}} , {\underline{\xi }}_{\mathrm{2}} , \dots , {\underline{\xi }}_N ]代表一套CPT数据， 其中 {\underline{\xi }}_n = [ {\xi }_{n,\mathrm{1}} = \mathrm{l} \mathrm{n} I_{c,\mathrm{1}} , {\xi }_{n,\mathrm{2}} = \mathrm{l} \mathrm{n} I_{c,\mathrm{2}} , \dots , {\xi }_{n,m} = \mathrm{l} \mathrm{n} I_{c,m} ]为第n层的m个测点I _c的对数值.根据给定的一套CPT数据划分土层，需要确定最可能土层数目N和每个土层的厚度 {\underline{H}}_N = \left\{H_n,n=\mathrm{1,2},\dots ,N\right\}（或相邻土层界面深度 {\underline{D}}_N = \left\{D_n,n=\mathrm{1,2},\dots ,N-\mathrm{1}\right\}）.根据贝叶斯定理，在给定一套CPT数据条件下，土层厚度和土层数目的联合后验分布 P \left({\underline{H}}_N,N|\underline{\xi }\right)为：

P \left({\underline{H}}_N,N|\underline{\xi }\right) = \frac{P\left(\underline{\xi }|{\underline{H}}_N,N\right)P\left({\underline{H}}_N,N\right)}{P\left(\underline{\xi }\right)}，(2)

式中： P \left({\underline{H}}_N,N\right)为土层厚度 {\underline{H}}_N和土层数目N的先验分布； P \left(\underline{\xi }|{\underline{H}}_N,N\right)为似然函数，反映了CPT数据提供的关于 {\underline{H}}_N和N的信息； P \left(\underline{\xi }\right)为归一化常数，与 {\underline{H}}_N和N无关.

式（2）中的后验分布 P \left({\underline{H}}_N,N|\underline{\xi }\right)是一个随着N变化的变维度概率密度函数，难以同时识别N和 {\underline{H}}_N.注意到，在求解似然函数 P \left(\underline{\xi }|{\underline{H}}_N,N\right)时，涉及高维积分问题.如图1所示，第n个土层似然函数为： P \left({\underline{\xi }}_n|{\underline{H}}_N,N\right) = \int \int P \left({\underline{\xi }}_n|{\underline{\mu }}_n,{\mathbf{\Sigma }}_n,{\underline{H}}_N,N\right) P \left({\underline{\mu }}_n,{\mathbf{\Sigma }}_n|{\underline{H}}_N,N\right) \mathrm{d} {\underline{\mu }}_n \mathrm{d} {\mathbf{\Sigma }}_n.为了解决变维度贝叶斯更新的求解问题，同时提高计算效率，本文基于分层贝叶斯学习框架提出了土层高效自动划分方法，具体表述如下.

为了表征土体的空间变异性，需要引入I_c 的概率模型，例如Cao et al.（2019）中采用的稳态随机场模型.对于稳态随机场而言，通常假定同一土层内的I_c 的均值、标准差和波动范围均为定值，同一土层内不同测点的空间相关性由自相关函数定量表征.因此，利用稳态随机场模型进行分析时，需要假设 I_c的自相关函数类型.然而，对于实际土层，土层剖面划定前土体参数的空间相关结构未知.本文采用的全高斯概率模型，采用均值向量 \underline{\mu }和协方差矩阵 \mathbf{\Sigma }来表征土体的空间变异性，无需事先进行任何相关结构假设.由于 I_c必为正值，本文采用对数正态随机场作为I_c 的概率模型 M_p \left(\underline{\theta }\right)，其中 \underline{\theta } = \left\{{\underline{\mu }}_n,{\mathbf{\Sigma }}_n,n=\mathrm{1,2},\dots ,N\right\}为模型参数.对于第n层土体而言， {\underline{\xi }}_n = \mathrm{l} \mathrm{n} {\underline{I}}_{cn} ~ N \left({\underline{\mu }}_n,{\mathbf{\Sigma }}_n\right).

根据条件概率，在给定一套CPT数据条件下，模型参数 \underline{\theta }、土层数目 N和土层厚度 {\underline{H}}_N的联合后验分布 P \left(\underline{\theta },{\underline{H}}_N,N|\underline{\xi }\right)可表达如下：

P \left(\underline{\theta },N,{\underline{H}}_N|\underline{\xi }\right) = P \left(N|\underline{\xi }\right) P \left({\underline{H}}_N|N,\underline{\xi }\right) P \left(\underline{\theta }|N,{\underline{H}}_N,\underline{\xi }\right)，(3)

式（3）中 P \left(N|\underline{\xi }\right)为给定 \underline{\xi }时，N的条件分布，表征了不同N值的发生概率； P \left({\underline{H}}_N|N,\underline{\xi }\right)为 {\underline{H}}_N的条件分布，量化了给定 \underline{\xi }、N时不同 {\underline{H}}_N值的可能性； P \left(\underline{\theta }|N,{\underline{H}}_N,\underline{\xi }\right)为 \underline{\theta }的条件分布，代表了给定 \underline{\xi }、N、 {\underline{H}}_N（即给定土层剖面）时不同 \underline{\theta }值的可能性.

式（3）可以将目标分布转换为3个条件概率分布，利用分层贝叶斯方法逐层学习目标变量（Lunn et al.，2012；Bozorgzadeh et al.，2019）.图2所示为所提土层划分方法的图模型，展示了变量及参数之间的关系，图中方形阴影填充标记代表模型变量的超参数，圆形阴影填充标记代表观测的CPT数据 \underline{\xi }.本节将对其计算求解过程进行详细介绍.

1.2　土层数目的后验概率

由于CPT测试深度有限，且其力学响应受到锥尖前后阻力的影响，测试深度范围内土层数目有限.令N _max为CPT测试深度范围内最大可能土层数目，每个可能的土层数目 N发生的概率为 P \left(N|\underline{\xi }\right)，其最大值 P \left(N^{\mathrm{*}}|\underline{\xi }\right)代表了最可能土层数目为 N^{\mathrm{*}}层.根据贝叶斯定理，土层数目 N的后验分布 P \left(N|\underline{\xi }\right)计算如下（Cao and Wang，2013；Wang et al.，2013， 2014）：

P \left(N|\underline{\xi }\right) = \frac{P\left(\underline{\xi }|N\right)P\left(N\right)}{P\left(\underline{\xi }\right)}，(4)

式中： P \left(\underline{\xi }|N\right)为给定土层数目N条件下 \underline{\xi }的概率密度函数，反映了CPT数据提供的关于含有N个土层的分层模型的信息，通常被称为“模型证据”； P \left(\underline{\xi }\right) = {\mathrm{\Sigma }}_{N=\mathrm{1}}^{N_{\mathrm{m}\mathrm{a}\mathrm{x}}} P \left(\underline{\xi }|N\right) P \left(N\right)为与土层数目 N无关的归一化常数.当缺乏N的先验信息时，可假设每个可能的 N值发生的概率相等，即 P \left(N\right)为均匀分布.此时， P \left(N\right)和 P \left(\underline{\xi }\right)均为常数， P \left(N|\underline{\xi }\right)与 P \left(\underline{\xi }|N\right)成正比.因此，最可能土层数目 N^{\mathrm{*}}对应 P \left(N^{\mathrm{*}}|\underline{\xi }\right)和 P \left(\underline{\xi }|N^{\mathrm{*}}\right)的最大值，可通过对比不同土层数目的模型证据 P \left(\underline{\xi }|N\right)确定最可能土层数目 N^{\mathrm{*}}.

1.3　土层厚度的后验概率分布

在识别土层厚度 {\underline{H}}_N时，可首先给定土层数目N.根据贝叶斯定理，土层厚度的后验分布 P \left({\underline{H}}_N|N,\underline{\xi }\right)为：

P \left({\underline{H}}_N|\underline{\xi },N\right) = \frac{P\left(\underline{\xi }|{\underline{H}}_N,N\right)P\left({\underline{H}}_N|N\right)}{P\left(\underline{\xi }|N\right)}，(5)

式中 P \left(\underline{\xi }|{\underline{H}}_N,N\right)为似然函数，反映了CPT数据提供的关于土层厚度的信息； P \left({\underline{H}}_N|N\right)为土层厚度的先验分布，反映了获取CPT数据之前关于土层厚度的先验信息； P \left(\underline{\xi }|N\right)为归一化常数，与土层厚度 N无关，即式（4）中的模型证据.

当土层数目为 N时，土层厚度 {\underline{H}}_N有无数种组合，其中任意一种组合之和为土层厚度H.记 \underline{h} = \left[h_{\mathrm{1}},h_{\mathrm{2}},\dots ,h_N\right]，其中 h_n = H_n / H为第n层土层厚度H_n 在土层总厚度H中所占比例，故 h_{\mathrm{1}} + h_{\mathrm{2}} + \dots + h_N = \mathrm{1}.假设 \underline{h}服从Dirichlet分布，其概率密度函数为（Bishop，2006）：

f \left(h_{\mathrm{1}},\dots ,h_N;{\alpha }_{\mathrm{1}},\dots ,{\alpha }_N\right) = \frac{\mathrm{1}}{B\left(\underline{\alpha }\right)} {\displaystyle \prod _{n=\mathrm{1}}^N} h_n^{{\alpha }_n-\mathrm{1}} = \frac{\mathrm{\Gamma }\left(\stackrel{N}{\underset{n=\mathrm{1}}{\mathrm{\Sigma }}}{\alpha }_n\right)}{{\displaystyle \prod _{n=\mathrm{1}}^N}\mathrm{\Gamma }\left({\alpha }_n\right)} {\displaystyle \prod _{n=\mathrm{1}}^N} h_n^{{\alpha }_n-\mathrm{1}}，(6)

其中， \underline{\alpha } = \left({\alpha }_{\mathrm{1}},\dots ,{\alpha }_N\right)为Dirichlet分布的参数. \underline{h} = \left[h_{\mathrm{1}},h_{\mathrm{2}},\dots ,h_N\right]分布在 N - \mathrm{1}标准单纯形 {\mathrm{\Omega }}_x中，且对于所有的 n \in \left\{\mathrm{1},\dots ,N\right\}均有 h_n \ge \mathrm{0}.当缺乏土层厚度的先验信息时，常采用对称Dirichlet分布作为土层层厚先验分布，此时 \underline{\alpha }中所有元素取值相等，式（6）表示为：

f \left(h_{\mathrm{1}},\dots ,h_N;\alpha \right) = \frac{\mathrm{\Gamma }\left(\alpha N\right)}{\mathrm{\Gamma }{\left(\alpha \right)}^N} {\displaystyle \prod _{n=\mathrm{1}}^N} h_n^{\alpha -\mathrm{1}}，(7)

在实际层厚样本空间 \mathrm{\Omega }中：

P \left({\underline{H}}_N|N\right) = \frac{\mathrm{\Gamma }\left(\alpha N\right)}{\mathrm{\Gamma }{\left(\alpha \right)}^N} {\displaystyle \prod _{n=\mathrm{1}}^N} H_n^{\alpha -\mathrm{1}} \cdot \frac{\mathrm{1}}{H^{N\alpha -\mathrm{1}}}，(8)

式中： H^{N\alpha -\mathrm{1}}为归一化常数.

1.4　基于正态‒逆威沙特分布的土层厚度似然函数

当土层数目为N时，式（5）中的似然函数 P \left(\underline{\xi }|{\underline{H}}_N,N\right)定量反映了CPT数据提供的关于土层厚度 {\underline{H}}_N的信息，假设不同土层的CPT数据 {\underline{\xi }}_n相互独立，则 P \left(\underline{\xi }|{\underline{H}}_N,N\right)可表达为：

P \left(\underline{\xi }|{\underline{H}}_N,N\right) = {\displaystyle \prod _{n=\mathrm{1}}^N} P \left({\underline{\xi }}_n|{\underline{H}}_N,N\right)，(9)

式中： P \left({\underline{\xi }}_n|{\underline{H}}_N,N\right)为第n个土层的似然函数，依赖于第n层土体空间变异性的概率模型及其模型参数 {\underline{\theta }}_n = \left\{{\underline{\mu }}_n,{\mathbf{\Sigma }}_n\right\}.根据贝叶斯理论，给定N、 {\underline{H}}_N、 {\underline{\xi }}_n条件下，第n个土层模型参数 {\underline{\theta }}_n = \left\{{\underline{\mu }}_n,{\mathbf{\Sigma }}_n\right\}的后验分布可表示为：

P \left({\underline{\mu }}_n,{\mathbf{\Sigma }}_n|{\underline{\xi }}_n,{\underline{H}}_N,N\right) = \frac{P\left({\underline{\xi }}_n|{\underline{\mu }}_n,{\mathbf{\Sigma }}_n,{\underline{H}}_N,N\right)P\left({\underline{\mu }}_n,{\mathbf{\Sigma }}_n|{\underline{H}}_N,N\right)}{P\left({\underline{\xi }}_n|{\underline{H}}_N,N\right)}，(10)

式中： P \left({\underline{\xi }}_n|{\underline{\mu }}_n,{\mathbf{\Sigma }}_n,{\underline{H}}_N,N\right)为似然函数， P \left({\underline{\mu }}_n,{\mathbf{\Sigma }}_n|{\underline{H}}_N,N\right)为模型参数的先验分布， P \left({\underline{\xi }}_n|{\underline{H}}_N,N\right) = \int \int P \left({\underline{\xi }}_n|{\underline{\mu }}_n,{\mathbf{\Sigma }}_n,{\underline{H}}_N,N\right) P \left({\underline{\mu }}_n,{\mathbf{\Sigma }}_n|{\underline{H}}_N,N\right) \mathrm{d} {\underline{\mu }}_n \mathrm{d} {\mathbf{\Sigma }}_n为与 {\underline{\mu }}_n、 {\mathbf{\Sigma }}_n无关的归一化常数，同时也是式（9）中第n个土层的似然函数.

在贝叶斯框架中，恰当的先验分布的选取至关重要，当先验分布为任意分布时，后验分布类型未知.当模型参数维度较高时，将会面临高维积分问题（Tsionas，2000）.本文采用全高斯概率模型表征土体的空间变异性， P \left({\underline{\xi }}_n|{\underline{\mu }}_n,{\mathbf{\Sigma }}_n,{\underline{H}}_N,N\right)服从多维正态分布.设模型参数的先验分布 P \left({\underline{\mu }}_n,{\mathbf{\Sigma }}_n|{\underline{H}}_N,N\right)采用正态‒逆威沙特（Normal-Inverse Wishart，NIW）分布（Murphy，2012）：

\begin{array}{l}\mathrm{N}\mathrm{I}\mathrm{W}\left({\underline{\mu }}_n,{\mathbf{\Sigma }}_n|{\underline{m}}_n,k_n,v_n,{\mathbf{S}}_n,{\underline{H}}_N,N\right)\triangleq \\ N\left({\underline{\mu }}_n|{\underline{m}}_n,\frac{\mathrm{1}}{k_n}{\mathbf{\Sigma }}_n,{\underline{H}}_N,N\right)\times IW\left({\mathbf{\Sigma }}_n|{\mathbf{S}}_n,v_n,{\underline{H}}_N,N\right)\end{array}，(11)

式中： N（·）为多维高斯分布， I W（·）为Inverse Wishart分布.本文中NIW先验分布参数取值如下： {\mathit{S}}_n = \mathrm{d} \mathrm{i} \mathrm{a} \mathrm{g} {\left\{{\sigma }_n^{\mathrm{2}},{\sigma }_n^{\mathrm{2}},\dots ,{\sigma }_n^{\mathrm{2}}\right\}}_{m\times m}为对角阵， {\sigma }_n为第n层m个CPT测点的标准差， {\underline{m}}_n = {\underline{\xi }}_n， v_n = m + \mathrm{2}，k_n 可取一个较小值（比如0.01），其中参数v_n 和k_n 的取值保证了所取先验为弱信息先验（Murphy， 2012）.由于NIW分布为似然函数 P \left({\underline{\xi }}_n|{\underline{\mu }}_n,{\mathbf{\Sigma }}_n,{\underline{H}}_N,N\right)的共轭先验， {\underline{\theta }}_n的后验分布亦为NIW分布：

P \left({\underline{\mu }}_n,{\mathbf{\Sigma }}_n|{\underline{\xi }}_n,{\underline{H}}_N,N\right) = \mathrm{N} \mathrm{I} \mathrm{W} \left({\underline{\mu }}_n,{\mathbf{\Sigma }}_n|{\underline{m}}_n^{\mathrm{\text{'}}},k_n^{\text{'}},v_n^{\text{'}},{\mathit{S}}_n^{\text{'}}\right)，(12)

式中： {\underline{m}}_n^{\mathrm{\text{'}}} = {\underline{m}}_n， k_n^{\text{'}} = k_n + \mathrm{1}， v_n^{\text{'}} = v_n + \mathrm{1}， {\mathit{S}}_n^{\text{'}} = {\mathit{S}}_n.

如前文所述，式（10）中 P \left({\underline{\xi }}_n|{\underline{H}}_N,N\right)为归一化常数，与模型参数 {\underline{\theta }}_n = \left\{{\underline{\mu }}_n,{\mathbf{\Sigma }}_n\right\}无关，根据全概率定理可计算如下：

P \left({\underline{\xi }}_n|H_N,N\right) = \int \int P \left({\underline{\xi }}_n|{\underline{\mu }}_n,{\mathbf{\Sigma }}_n,{\underline{H}}_N,N\right) P \left({\underline{\mu }}_n,{\mathbf{\Sigma }}_n|{\underline{H}}_N,N\right) \mathrm{d} {\underline{\mu }}_n \mathrm{d} {\mathbf{\Sigma }}_n.(13)

由于 P \left({\underline{\mu }}_n,{\mathbf{\Sigma }}_n|{\underline{H}}_N,N\right)为 P \left({\underline{\xi }}_n|{\underline{\mu }}_n,{\mathbf{\Sigma }}_n,{\underline{H}}_N,N\right)的共轭先验，式（13）中积分项存在解析解， P \left({\underline{\xi }}_n|{\underline{H}}_N,N\right)计算如下（Murphy，2012）：

P \left({\underline{\xi }}_n|{\underline{H}}_N,N\right) = \frac{\mathrm{1}}{{\mathrm{\pi }}^{m/\mathrm{2}}} {\left(\frac{k_n}{k_n^{\text{'}}}\right)}^{m/\mathrm{2}} \frac{{\left|{\mathit{S}}_n\right|}^{v_n/\mathrm{2}}}{{\left|{\mathit{S}}_n^{\text{'}}\right|}^{v_n^{\text{'}}/\mathrm{2}}} \frac{{\mathrm{\Gamma }}_m\left(v_n^{\text{'}}/\mathrm{2}\right)}{{\mathrm{\Gamma }}_m\left(v_n/\mathrm{2}\right)}，(14)

基于式（14）的计算结果，代入式（9）可得到整个土层剖面的似然函数值，避免了 P \left(\underline{\xi }|{\underline{H}}_N,N\right)求解时的积分计算，有助于提升计算效率.根据本节的推导过程，虽然在建模过程中需要随机场模型参数 {\underline{\theta }}_n = \left\{{\underline{\mu }}_n,{\mathbf{\Sigma }}_n\right\}，在土层识别过程中并不关心它们的具体值，通过式（13）中的积分计算避免了识别随机场模型参数，直接识别土层数目和土层厚度.虽然所提分层贝叶斯分析方法图模型中包含了土层数目、土层厚度和随机场模型参数（见图2），但是本研究的识别目标是土层数目和厚度，无需显式识别随机场模型参数.

由上述介绍可知，求解后验分布 P \left({\underline{H}}_N|\underline{\xi },N\right)和模型证据 P \left(\underline{\xi }|N\right)是定量表征土层剖面不确定性的关键.本文采用基于子集模拟的贝叶斯更新方法求解 P \left({\underline{H}}_N|\underline{\xi },N\right)和 P \left(\underline{\xi }|N\right)，介绍如下.

2 基于子集模拟的贝叶斯更新方法（aBUS_SuS）

基于结构可靠度的贝叶斯更新方法（BUS）（Straub and Papaioannou，2015；DiazDelaO et al.，2017；Betz et al.，2018）利用似然函数 P \left(\underline{\xi }|{\underline{H}}_N,N\right)和均匀随机变量U构建失效域F，将贝叶斯更新问题转换为结构可靠度问题.因此，可靠度分析方法可用来求解后验分布 P \left({\underline{H}}_N|\underline{\xi },N\right).利用BUS算法，构建的失效域F定义如下：

，(15)

其中：U为服从（0，1）均匀分布的随机变量； P \left(\underline{\xi }|{\underline{H}}_N,N\right)为式（9）中定义的似然函数.在给定土层数目N下，c为保证所有 {\underline{H}}_N满足 c^{-\mathrm{1}} \ge P \left(\underline{\xi }|{\underline{H}}_N,N\right)的常数.根据抽样拒绝准则（Straub and Papaioannou，2015；DiazDelaO et al.，2017；Betz et al.，2018），从 P \left({\underline{H}}_N|N\right)中抽取的随机样本，当其落入失效域 F时，即符合后验分布 P \left({\underline{H}}_N|\underline{\xi },N\right).由式（15）可知，执行BUS算法前，需要事先确定乘数c，通常十分困难.

为了解决该问题，Betz et al.（2018）提出了基于子集模拟（SuS）的自适应BUS算法（aBUS_SuS），在子集模拟执行过程中，自适应地确定c.对于高维贝叶斯更新问题，aBUS_SuS提供了一种高效求解方法.在aBUS_SuS中，中间失效域 F_i定义如下：

，(16)

式中：h_i 为中间失效事件的阈值.为自适应地确定未知的乘数 c，式（16）中的中间失效域F_i 可等效表达为（Betz et al.，2018）：

，(17)

式中： h_i^{\mathrm{*}} = h_i + \mathrm{l} \mathrm{n} \left(c/c^{\mathrm{*}}\right).根据式（16）和式（17）可知，F_i 可由常数 c^{\mathrm{*}}、中间阈值水平 h_i^{\mathrm{*}}定义.

在执行过程中，首先将c ^-1初始值设置为采用直接蒙特卡洛模拟方法从先验分布 P \left({\underline{H}}_N|N\right)中抽取的样本的似然函数最大值；然后，根据给定的c ^-1和条件概率p ₀，确定第一层的中间阈值水平h ₁；接着，在失效域F ₁中通过马尔科夫链蒙特卡洛模拟产生N _cs个条件样本.令 {\left(c^{\mathrm{*}}\right)}^{-\mathrm{1}}为N _cs个条件样本似然函数最大值，若 ，将 c^{-\mathrm{1}}更新为 {\left(c^{\mathrm{*}}\right)}^{-\mathrm{1}}，此时 c^{-\mathrm{1}} \ge P \left(\underline{\xi }|{\underline{H}}_N,N\right)仍然成立；相应地，用 h_i^{\mathrm{*}}代替h_i 确保中间失效域F_i 保持不变（见式（16）和式（17））.在每一个中间失效层，c ^-1和h_i 都采用这样的方式自适应地进行调整.aBUS_SuS一直执行到 {h_i}^{\mathrm{*}} = \mathrm{0}，此时产生的样本已全部落入目标失效域，即为服从后验分布的样本，此时子集模拟自动终止运算，同时可计算失效概率 P_{\mathrm{\Omega }}，对应的模型证据计算如下：

P \left(\underline{\xi }|N\right) = c^{-\mathrm{1}} \cdot P_{\mathrm{\Omega }}，(18)

根据模型证据 P \left(\underline{\xi }|N\right)和后验样本分别可以确定最可能土层数目和厚度（或界面深度）.在给定土层数目N时，根据式（5）可确定 {\underline{H}}_N后验样本概率密度函数 P \left({\underline{H}}_N|\underline{\xi },N\right)最大值，其对应的土层剖面，即为最可能土层剖面.

3 计算流程

本文所提方法的计算流程，主要包括以下6个步骤：

（1） 通过静力触探试验获取一套CPT数据 \underline{\xi }.

（2） 根据现场勘探信息及工程经验，确定最大可能土层数目N _max.

（3） 当土层数目N=1时，土层厚度为触探深度H，模型证据对数值 \mathrm{l} \mathrm{n} P \left(\underline{\xi }|N=\mathrm{1}\right)通过式（14）计算.

（4） 当土层数目N>1时，根据式（9）构建似然函数及式（15）中的失效域F，采用aBUS_SuS产生土层厚度 {\underline{H}}_N的后验样本，并计算模型证据对数值 \mathrm{l} \mathrm{n} P \left(\underline{\xi }|N\right).

（5） 重复步骤（4）N _max－1次，比较不同土层数目N的模型证据对数值，确定最可能土层数目N ^*.

（6） 根据N ^*对应的土层厚度的后验样本计算N ^*-1个界面深度，确定最可能界面深度 {\underline{D}}_N^{\mathrm{*}}，并统计出界面深度的均值和标准差，界面深度的标准差反映了其不确定性.

4 算例分析

4.1　工程概况

本节通过温州市仙湖调蓄工程试验区现场CPT数据（见图3a和3b）说明所提方法的有效性.图3c为根据式（1）计算得到的一套I_c 数据，图3d为根据钻孔取样绘制的试验区土层剖面图.如图3所示，CPT测试深度范围内土体共分为6层，主要为滨海相软土.根据钻孔，第3~5层均为淤泥层（Q₄ m），但其平均含水量逐渐降低，分别为71.5%、65.0%和56.4%，3层淤泥层均夹杂少量粉细砂，同时第4层局部夹少量贝壳碎屑，第5层局部夹少量腐殖质.由于CPT本身无法提供土样，因此无法直接绘制土层剖面图.本文所提方法主要应用于没有钻探土样而只有CPT勘探数据的现场，应用CPT数据 \underline{\xi }进行土体的力学分层.

本算例选取N _max=10（即土层数目N的取值范围为1~10），对于每个可能的N值，采用aBUS_SuS产生土层厚度 {\underline{H}}_N的后验样本，并计算模型证据对数值 \mathrm{l} \mathrm{n} P \left(\underline{\xi }|N\right).在aBUS_SuS执行过程中，子集模拟每层样本数量为10 000，p ₀设置为0.1，子集模拟执行终止时，共产生10 000个后验样本.同时，Dirichlet分布中α取4，关于α取值及其对模型选择的影响将在第5节讨论.

4.2　土层剖面识别结果

图4所示为不同土层数目N对应的模型证据对数值 \mathrm{l} \mathrm{n} P \left(\underline{\xi }|N\right).当N＜4时， \mathrm{l} \mathrm{n} P \left(\underline{\xi }|N\right)逐渐增大；当N>4时， \mathrm{l} \mathrm{n} P \left(\underline{\xi }|N\right)逐渐减小.随着N的增大，模型证据对数值 \mathrm{l} \mathrm{n} P \left(\underline{\xi }|N\right)先增大后减小，与文献中（如Cao and Wang，2013；Wang et al.，2013）观察的现象一致，故本算例中最可能土层数目N^* =4.当土层数目N=4时，根据aBUS_SuS产生的土层厚度 {\underline{H}}_N的后验样本，可估计每个土层厚度的概率分布和统计特征（平均值、标准差、最可能值）以及相应的土层界面深度 {\underline{D}}_N的后验分布.本文采用界面深度最可能值划分土层，使用土层界面深度后验样本标准差定量反映土层剖面的不确定性.

图5a所示为本文所提方法划分土层剖面，共4个土层，图中红色实线代表划分土层界面位置，灰色竖虚线代表不同土类的I_c 边界（见表1）.如图5a所示，所划分土层剖面同一土层内I_c 值表现出相似的统计特性，而相邻土层内I_c 的统计特性存在显著差异，如第2层和第3层内I_c 值均未表现出明显的趋势性，但第2层内I_c 的空间变异性强于第3层，由此可见本文所提方法划分的土层剖面合理地反映了土体力学响应的统计相似性.

4.3　对比分析

根据Cao et al.（2019）所提方法，使用单指数型自相关函数表征不同测点之间的空间相关性，识别土层剖面如图5b所示，与图5a划分的结果基本一致，但如图5c所示，本文所提方法划分土层边界的不确定性相对较小，主要原因在于Cao et al.（2019）使用均匀分布作为模型参数的先验分布，而本文采用相对更有弱信息的NIW分布.在计算效率方面，使用16.0 GB内存、3.20 GHz主频的Intel Core i7-8700 CPU处理器配置的台式计算机，采用Cao et al.（2019）中提出的算法耗时约为5.73 h，而所提方法计算时间约为7.56 min，显著提高了计算效率.

现场钻孔取土分层结果如图3d所示，与所提方法划分土层剖面存在差异，可能原因如下：钻孔取土分层结果受现场取样、样品运输、室内试验结果以及工程师的工程经验等因素影响；同时，CPT在进行测试过程中，其力学响应会受到锥尖前后阻力的影响，如所提方法划分的一个土层界面位置为9.79 m，而现场钻孔划分界面位置约为12 m，可能就是由于“过渡区”的存在造成的差异.另外，现场钻孔取土分层只能考虑土体的物理特性（比如粒径分布和含水率），而工程师往往更关注于土体的力学响应，所提方法可在考虑土体力学响应统计特征条件下划分土层，为岩土工程分析和设计提供了合理的土层信息.

5 方法验证

5.1　模拟I_c 数据

在基于静力触探的岩土工程勘察中，工程师只能根据CPT数据推断土层数目和土层厚度，而实际土层界面位置未知.为验证本文所提方法，本节采用对数正态随机场模型，使用单指数型自相关函数，生成一套已知土层界面位置的模拟数据.令 {\underline{I}}_{cn}表示第n个土层I_c 的对数正态随机场，随机场模型参数 {\underline{\theta }}^{\mathrm{\text{'}}} = \left\{{\mu }_n^{\text{'}},{\sigma }_n^{\text{'}},{\lambda }_n^{\text{'}},n=\mathrm{1,2},\dots ,N\right\}，其中 {\mu }_n^{\text{'}}、 {\sigma }_n^{\text{'}}、 {\lambda }_n^{\text{'}}分别为第 n个土层的均值、标准差和波动范围.

如图6所示，该模拟场地土层数目N=5，土层厚度分别为2 m、3 m、5 m、10 m和15 m，对应的土层界面深度为2 m、5 m、10 m和20 m.每个土层内模拟数据点间隔为0.05 m，随机场参数取值见图6a，采用相关矩阵分解法模拟I_c 数据.如图6b所示，红色实线代表实际土层界面位置，在识别过程中假设未知.根据图6b所示模拟I_c 数据采用所提方法对该模拟场地进行土层划分，将识别土层结果与实际土层剖面进行对比，以验证所提方法的合理性.需要注意的是，用于生成模拟I_c 数据的概率模型及随机场参数仅用于生成各层具有统计特性差异的模拟I_c 数据，而在土层识别过程中未使用任何相关信息.

5.2　模拟数据识别结果

本节中仍取最大土层数目N _max=10（即N的取值范围为1~10），采用aBUS_SuS产生土层厚度 {\underline{H}}_N的后验样本，并计算模型证据对数值 \mathrm{l} \mathrm{n} P \left(\underline{\xi }|N\right).如图7所示，不同土层数目N下的模型证据对数值 \mathrm{l} \mathrm{n} P \left(\underline{\xi }|N\right)先增大后减小，与第4节中算例结果类似.当N=5时， \mathrm{l} \mathrm{n} P \left(\underline{\xi }|N\right)最大，因此N^* =5为最可能土层数目，与实际土层数目一致，验证了本文所提方法.

图8a所示为不同土层数目N下最可能土层界面深度的识别结果，所识别的土层剖面随 N的变化趋势.当N=2时，根据所提方法识别出的土层界面深度为20.00 m；当N=3时识别出2个土层界面，分别在10.02 m和20.04 m处，包含了 N=2时识别出的土层界面位置.所提方法首先识别出钻探深度内相邻两层I_c 统计特征存在显著差异的土层界面，如本例中模拟场地的实际土层边界，随着N值的增大，更多的土层界面位置可逐渐被识别出来.如图8b所示，当N=5时，根据所提方法识别出的土层界面深度分别为： 2.04 m、5.04 m、10.04 m和20.03 m，识别结果与实际土层剖面一致，验证了所划分土层的合理性.

图8c为根据所提方法识别的土层界面深度的标准差，第1个边界深度 D_{\mathrm{1}}^{\mathrm{*}}（2.04 m）的标准差最大，第4个边界深度 D_{\mathrm{4}}^{\mathrm{*}}（20.03 m）的标准差次之，而第2个边界深度 D_{\mathrm{2}}^{\mathrm{*}}（5.04 m）和第3个边界深度 D_{\mathrm{3}}^{\mathrm{*}}（10.04 m）的标准差相对较小.该结果与图8b所示的I_c 数据一致，如第2层与第3层I_c 统计特征差异性较大，类似第3层与第4层之间I_c 统计特征差异性也相对较大，因此 D_{\mathrm{2}}^{\mathrm{*}}、 D_{\mathrm{3}}^{\mathrm{*}}的标准差相对较小.另外，所提方法识别土层剖面 D_{\mathrm{1}}^{\mathrm{*}}的标准差最大为 0.044 m，小于模拟I_c 数据点间隔0.05 m，说明所提方法识别结果可靠性较高.所提方法根据I_c 数据统计特征识别土层剖面，合理地表征了界面深度的不确定性，且能定量反映识别结果的可靠性.

5.3　关于Dirichlet分布 \alpha取值的讨论

为探讨Dirichlet分布 \alpha取值对土层识别结果的影响，采用5.1节模拟I_c 数据进行了20次重复计算，采用不同 \alpha正确识别土层数目的次数和正确率如图9所示.

当 \alpha = \mathrm{1}时，土层数目正确识别率为10%；随着 \alpha增大，所提方法正确识别率逐渐提高，当 \alpha = \mathrm{4}时，正确识别率达到100%.值得注意的是，当 \alpha = \mathrm{1}时，对称Dirichlet分布等效于N-1标准单纯形上的均匀分布，即在样本空间内数据点均匀分布，随着 \alpha的增大，样本空间边缘数据点发生的概率逐渐降低，即样本逐渐向样本空间中心聚拢.Nobile（2004）指出Dirichlet先验分布的取值对于模型后验分布有较大的影响.Frühwirth-Schnatter（2006）通过模拟数据试验发现，当 \alpha =1时，真实模型的后验概率远小于 \alpha =4时的后验概率，同时，在避免过拟合能力方面， \alpha =4优于其他取值.

本文中 \alpha的取值决定了含有I_c 数据点较少甚至不含数据点的薄层发生的可能性，同时会对模型选择（即最可能土层数目识别）产生影响.由于样本数量有限，当 \alpha超过一定值时，可能造成真实模型样本的丢失，产生“漏层”现象.因此，在使用所提方法进行土层划分时，为了有效提升正确识别率，同时避免潜在的“漏层”现象，建议 \alpha的取值为4.

6 结论

本文提出了一种基于分层贝叶斯学习的滨海软土复杂地层快速识别方法，包括分层贝叶斯学习框架及分层学习方法.本文使用全高斯概率模型打破了计算协方差矩阵时需要给定自相关函数的限制.同时，引入NIW分布显著提升了计算效率.通过温州市仙湖调蓄工程试验区现场CPT数据和模拟I_c 数据说明了所提方法的有效性和合理性.主要结论如下：

（1） 本文从分层贝叶斯学习的角度提出了基于土类指数I_c 的土层划分方法，并给出了相应的图模型.采用均值向量 \underline{\mu }和协方差矩阵 \mathbf{\Sigma }来表征土体的空间变异性，无需事先假设相关结构.

（2） 本文引入NIW分布作为全高斯随机场模型参数的先验分布，避免了对随机场参数进行高维积分，有效地提升了土层厚度的似然函数和后验分布的计算效率.

（3） 所提方法不仅能在考虑CPT数据空间变异性的前提下合理划分土层，同时可定量表征土层识别不确定性.通过引入全高斯随机场模型以及随机场模型参数的弱信息先验分布，所提方法识别的土层剖面不确定性较小、可靠性较高.

（4） 钻孔取样分层主要考虑土体物性特征（比如粒径分布和含水率），无法考虑其力学参数，同时划分的土层剖面与工程师的经验相关，其可靠性未知.相反，所提方法根据I_c 统计特性划分土层，当相邻两层I_c 数据统计特性差异较大时，不确定性较小.

参考文献

[1] Ahmadi, M. M., Robertson, P. K., 2005. Thin-Layer Effects on the CPT qc Measurement. Canadian Geotechnical Journal, 42(5): 1302-1317.https://doi.org/10.1139/t05-036

[2] Betz, W., Papaioannou, I., Beck, J. L., et al., 2018. Bayesian Inference with Subset Simulation: Strategies and Improvements. Computer Methods in Applied Mechanics and Engineering, 331: 72-93. https://doi.org/10.1016/j.cma.2017.11.021

[3] Bishop, C. M., 2006. Pattern Recognition and Machine Learning. Springer, New York, 76-78.

[4] Bozorgzadeh, N., Harrison, J. P., Escobar, M. D., 2019. Hierarchical Bayesian Modelling of Geotechnical Data: Application to Rock Strength. Géotechnique, 69(12): 1056-1070. https://doi.org/10.1680/jgeot.17.P.282

[5] Cai, G.J., Liu, S.Y., Tong, L.Y., et al., 2009. Soil Classification Using CPTU Data Based upon Cluster Analysis Theory. Chinese Journal of Geotechnical Engineering, 31(3): 416-424 (in Chinese with English abstract).

[6] Cao, Z. J., Wang, Y., 2013. Bayesian Approach for Probabilistic Site Characterization Using Cone Penetration Tests. Journal of Geotechnical and Geoenvironmental Engineering, 139(2): 267-276. https://doi.org/10.1061/(ASCE)GT.1943-5606.0000765

[7] Cao, Z. J., Zheng, S., Li, D. Q., 2019. Bayesian Identification of Soil Stratigraphy Based on Soil Behaviour Type Index. Canadian Geotechnical Journal, 56(4): 570-586. https://doi.org/10.1139/cgj-2017-0714

[8] Chen, L., Lin, W. B., Chen, P., et al., 2021. Porosity Prediction from Well Logs Using Back Propagation Neural Network Optimized by Genetic Algorithm in One Heterogeneous Oil Reservoirs of Ordos Basin, China. Journal of Earth Science, 32(4): 828-838. https://doi.org/10.1007/s12583-020-1396-5

[9] Chen, Y. H., 2020. Variability Analysis of Undrained Shear Strength of Coastal Soft Soil (Dissertation). Southeast University, Nanjing (in Chinese with English abstract).

[10] Ching, J., Wang, J. S., Juang, C. H., et al., 2015. Cone Penetration Test (CPT)-Based Stratigraphic Profiling Using the Wavelet Transform Modulus Maxima Method. Canadian Geotechnical Journal, 52(12): 1993-2007. https://doi.org/10.1139/cgj-2015-0027

[11] DiazDelaO, F. A., Garbuno-Inigo, A., Au, S. K., 2017. Bayesian Updating and Model Class Selection with Subset Simulation. Computer Methods in Applied Mechanics and Engineering, 317: 1102-1121. https://doi.org/10.1016/j.cma.2017.01.006

[12] Frühwirth-Schnatter, S., 2006. Finite Mixture and Markov Switching Models. Springer, New York, 141-143.

[13] Guan, Z., Wang, Y., 2022. CPT-Based Probabilistic Liquefaction Assessment Considering Soil Spatial Variability,Interpolation Uncertainty and Model Uncertainty. Computers and Geotechnics, 141: 104504. https://doi.org/10.1016/j.compgeo.2021.104504

[14] Hegazy, Y. A., Mayne, P. W., 2002. Objective Site Characterization Using Clustering of Piezocone Data. Journal of Geotechnical and Geoenvironmental Engineering, 128(12): 986-996. https://doi.org/10.1061/(ASCE)1090-0241(2002)128:12(986)

[15] Huang, K., Sun, R.L., Yuan, S.Q., et al., 2022. Effect of Number of Pumping Tests and Prior Information on Hydraulic Conductivity Estimation of Three-Dimensional Heterogeneous Aquifer. Earth Science, 47(2): 689-699 (in Chinese with English abstract).

[16] Jefferies, M. G., Davies, M. P., 1993. Use of CPTu to Estimate Equivalent SPT N₆₀. Geotechnical Testing Journal, 16(4): 458-468. https://doi.org/10.1520/gtj10286j

[17] Ku, C. S., Juang, C. H., Ou, C. Y., 2010. Reliability of CPT Ic as an Index for Mechanical Behaviour Classification of Soils. Géotechnique, 60(11): 861-875. https://doi.org/10.1680/geot.09.P.097

[18] Liu, S.Y., Cai, G.J., Zou, H.F., 2013. Practical Soil Classification Methods in China Based on Piezocone Penetration Tests. Chinese Journal of Geotechnical Engineering, 35(10): 1765-1776 (in Chinese with English abstract).

[19] Lunn, D., Jackson, C., Best, N., et al., 2012. The BUGS Book: A Practical Introduction to Bayesian Analysis. CRC Press, New York, 219-252.

[20] Miu, L. C., 2012. Mechanical Properties and Engineering Practice of Soft Soil. Science Press, Beijing, 3-7 (in Chinese).

[21] Murphy, K. P., 2012. Machine Learning: A Probabilistic Perspective. MIT Press, London, 125-161.

[22] Nobile, A., 2004. On the Posterior Distribution of the Number of Components in a Finite Mixture. The Annals of Statistics, 32(5): 2044-2073. https://doi.org/10.1214/009053604000000788

[23] Robertson, P. K., 1990. Soil Classification Using the Cone Penetration Test. Canadian Geotechnical Journal, 27(1): 151-158. https://doi.org/10.1139/t90-014

[24] Robertson, P. K., 2009. Interpretation of Cone Penetration Tests—A Unified Approach. Canadian Geotechnical Journal, 46(11): 1337-1355. https://doi.org/10.1139/T09-065

[25] Robertson, P. K., Campanella, R. G., 1983. Interpretation of Cone Penetration Tests. Part I: Sand. Canadian Geotechnical Journal, 20(4): 718-733. https://doi.org/10.1139/t83-078

[26] Robertson, P. K., Wride, C. E., 1998. Evaluating Cyclic Liquefaction Potential Using the Cone Penetration Test. Canadian Geotechnical Journal, 35(3): 442-459. https://doi.org/10.1139/t98-017

[27] Straub, D., Papaioannou, I., 2015. Bayesian Updating with Structural Reliability Methods. Journal of Engineering Mechanics, 141(3): 04014134. https://doi.org/10.1061/(ASCE)EM.1943-7889.0000839

[28] Tsionas, E. G., 2000. Numerical Bayesian Inference with Arbitrary Prior. Statistical Papers, 41(4): 437-451. https://doi.org/10.1007/BF02925762

[29] Wang, D. T., Chen, G. X., 2022. Seismic Wave Impedance Inversion Based on Temporal Convolutional Network. Earth Science, 47(4): 1492-1506 (in Chinese with English abstract).

[30] Wang, X., Wang, H., Liang, R. Y., et al., 2019. A Semi-supervised Clustering-Based Approach for Stratification Identification Using Borehole and Cone Penetration Test Data. Engineering Geology, 248: 102-116. https://doi.org/10.1016/j.enggeo.2018.11.014

[31] Wang, Y., Huang, K., Cao, Z. J., 2013. Probabilistic Identification of Underground Soil Stratification Using Cone Penetration Tests. Canadian Geotechnical Journal, 50(7): 766-776. https://doi.org/10.1139/cgj-2013-0004

[32] Wang, Y., Huang, K., Cao, Z. J., 2014. Bayesian Identification of Soil Strata in London Clay. Géotechnique, 64(3): 239-246. https://doi.org/10.1680/geot.13.T.018

[33] Wu, S., Ching, J., Phoon, K. K., 2022. Quasi-Site-Specific Soil Property Prediction Using a Cluster-Based Hierarchical Bayesian Model. Structural Safety, 99: 102253. https://doi.org/10.1016/j.strusafe.2022.102253

[34] Zhang, C. H., Shi, J., Dai, J. Q., 1997. The Application of Piezocone Tests in China. Chinese Journal of Geotechnical Engineering, 19(1): 50-57 (in Chinese with English abstract).

[35] Zhang, Z., Tumay, M. T., 1999. Statistical to Fuzzy Approach toward CPT Soil Classification. Journal of Geotechnical and Geoenvironmental Engineering, 125(3): 179-186. https://doi.org/10.1061/(ASCE)1090-0241(1999)125:3(179)

[36] Zhao, T. Y., Wang, Y., 2020. Interpolation and Stratification of Multilayer Soil Property Profile from Sparse Measurements Using Machine Learning Methods. Engineering Geology, 265: 105430. https://doi.org/10.1016/j.enggeo.2019.105430

[37] Zheng, S., Cao, Z. J., Li, D. Q., 2019. A Fast Bayesian Approach for CPT-Based Soil Stratification with BUS. In: Ching, J., Li, D. Q., Zhang, J., eds, Proceedings, 7th International Symposium on Geotechnical Safety and Risk, Research Publishing, Singapore, 364-369.

[38] 蔡国军, 刘松玉, 童立元, 等, 2009. 基于聚类分析理论的CPTU土分类方法研究. 岩土工程学报, 31(3): 416-424.

[39] 陈宇航, 2020. 滨海软弱土不排水抗剪强度变异性分析（硕士学位论文）. 南京: 东南大学.

[40] 黄康, 孙蓉琳, 袁淑卿, 等, 2022. 抽水组数和先验信息对估算三维非均质含水层渗透系数的影响. 地球科学, 47(2): 689-699.

[41] 刘松玉, 蔡国军, 邹海峰, 2013. 基于CPTU的中国实用土分类方法研究. 岩土工程学报, 35(10): 1765-1776.

[42] 谬林昌，2012.软土力学特性与工程实践.北京：科学出版社，3-7.

[43] 王德涛, 陈国雄, 2022. 基于时间卷积网络的地震波阻抗反演. 地球科学, 47(4): 1492-1506.

[44] 张诚厚, 施健, 戴济群, 1997. 孔压静力触探试验的应用. 岩土工程学报, 19(1): 50-57.

基金资助

国家自然科学基金项目(52278368;52025094)