0.   引言

根据勘探资料合理识别场地特征是岩土工程分析设计的基础，包括土体参数、土层剖面信息等^[1]。静力触探试验（cone penetration test，CPT）因其具有操作误差小、可重复性高、快速提供土层沿深度分布的力学特性等优点，在岩土工程场地勘察中广泛应用^[2-3]。比如，利用CPT数据与土体力学性质相关的特征，基于CPT数据的变化趋势解译土层信息，根据土体沿深度分布的力学特性识别土层剖面，与基于现场取样的室内土工试验土层划分结果对比分析，为岩土工程设计提供依据^[4]。

为此，文献中基于静力触探试验数据提出了土体分类指数I_c和基于I_c的土体经验分类图表^[5-10]。其中，刘松玉等^[10]通过分析江苏省不同地质成因的14个试验场地中95套CPTU及其邻近钻孔数据信息，发现Robertson等^[8]提出的土体分类图表对中国土体分类标准较为适用，并对其进行修正，使之与中国《岩土工程勘察规范（GB 50021-2001）》^[11]中土体分类标准一致。因此，在利用CPT数据计算得到I_c值后，可以根据修正后的土体经验分类图表划分土层剖面。

需要指出，CPT无法直接提供土样，需要根据I_c值推断土体类别。然而，基于I_c值的力学土体划分标准与根据钻孔取样的土体物性分类标准有所不同，导致基于CPT和钻孔的土层划分结果存在差异^[12]。然而，岩土体作为一种天然材料，在形成过程中受到搬运介质、沉积环境等因素的影响，具有碎散性、多相性和天然性，导致其物理力学性质复杂，相应的I_c值沿深度方向存在离散性和空间变异性。若仅利用沿深度变化的I_c值及土体经验分类图表划分土层剖面，可能存在过多薄层，与工程实际相悖。因此，工程应用中往往根据工程经验将薄层合并，以减少土层数目，获得合理的土层剖面^[14]。然而，不同工程师的判断标准及经验不同，导致土层剖面划分结果存在差异，存在主观不确定性。

针对上述问题，文献中提出了基于机器学习模型的土层划分方法，包括模糊集分析^[14]、聚类分析^[15-17]、小波分析^[18]、贝叶斯系统识别^{[13, 19]}和隐式马尔科夫随机场^[20-21]等。Cao等^[19]提出了土层剖面自动识别与不确定性表征贝叶斯方法，采用二阶平稳随机场模型表征I_c数据的空间变异性。相对于其它方法，Cao等^[19]所提方法不仅能够识别统计均匀的土层剖面，而且能够定量表征土层边界的不确定性。然而，由于需要通过高维积分计算不同土层数目的模型证据确定最可能土层剖面，当CPT探测深度较深、数据量较大时，该方法的计算效率显著降低；另一方面，由于模型证据的求解受到数据量的影响，导致该方法对于测试深度较短的CPT数据（如滨海滩地勘察）的适用性较差。此外，利用该方法识别土层剖面时，需要提前假定土体参数的空间相关结构，但是实际土层中土体力学参数的空间相关结构往往未知。因此，如何在考虑空间变异性条件下根据CPT数据高效地划分土层力学剖面，合理地反映不同土层土体力学响应的差异性仍是一个关键难点。

在贝叶斯学习框架下，提出了一种基于静力触探试验和变维联合后验分布的土层剖面高效优化识别方法。基于模拟退火优化算法，所提方法将基于CPT数据的土层力学分层问题定义为贝叶斯后验分布优化问题，直接计算变维联合后验分布，快速确定最可能土层剖面（包括土层数目和土层厚度）。此外，所提方法采用全高斯概率模型表征土体空间变异性，无需假定土体参数空间相关结构。为此，首先介绍了基于静力触探试验和变维联合后验分布的土层剖面高效优化识别方法及其计算流程，通过杭州某地铁区间现场CPT数据和两套模拟I_c数据说明了所提方法的有效性和合理性，并结合土层识别结果说明了所提方法的分层原理和特点。

1.   基于CPT数据的土层力学剖面参数变维联合后验分布

1.1   给定CPT数据条件下土层识别贝叶斯方法

在给定一套CPT数据条件下，划分土层剖面需要确定最可能土层数目N和土层厚度H_N = [H_n, n =1, 2, 3, …, N]（或相邻土层界面深度D_N = [D_n, n = 1, 2, 3, …, N-1]）。根据CPT数据划分土层可以定义为求解土层厚度H_N和土层数目N的联合条件概率分布P(H_N, N|$\underline \xi$)。其中，$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi } = [{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi } _n}$, n = 1, 2, 3, …, N]为根据CPT数据计算的I_c的对数值（即$\underline \xi$ = lnI_c）；${\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi } _n}$ = [ξ_{n, i} = lnI_{c, i}, i = 1, 2, 3, …, m_n]表示第n个土层中m_n个I_c的对数值。在贝叶斯框架下，该分布可以表示为^[19]

式中：P(H_N, N)为土层厚度H_N和土层数目N的联合先验概率分布；P($\underline \xi$ | H_N, N)为似然函数；P($\underline \xi$) = ∫∫P($\underline \xi$ | H_N, N)P(H_N, N)dH_N dN为归一化常数，与土层厚度H_N和土层数目N无关。

需要指出，P(H_N, N)为变维度概率密度函数，即H_N和N同时发生改变，通常需要先确定最可能土层数目，再划分最可能土层界面^[19]。为避免同时识别参数H_N和N带来的求解难度，根据条件概率公式对P(H_N, N)展开：

将式（2）代入式（1）可得

在给定$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$条件下，通过设置一系列可能的N值（即N = 1, 2, 3, …, N_max），其中N_max为CPT触探深度范围内可能的最大土层数目，计算所有候选土层剖面P(H_N, N| $\underline \xi$)的最大值，即可确定最可能的土层数目和对应的土层剖面，分别用N^*和H_N^*表示。由式（3）可知，求解H_N和N的联合后验分布P(H_N, N| $\underline \xi$)需要确定先验概率分布P(H_N | N)和P(N)以及似然函数P($\underline \xi$| H_N, N)，分述如下。

1.2   先验分布

先验分布P(H_N, N)反映了获取CPT数据之前的土层剖面信息，包括土层数目和土层厚度。通过条件概率公式对P(H_N, N)展开后（式（2）），P(N)为N的先验概率分布，P(H_N |N)为给定N值条件下H_N的先验概率密度函数，反映了获取现场勘察信息之前关于H_N的信息。当CPT触探深度内最大可能土层数目为N_max且无其它关于N的先验信息时，可以假定每个可能的N值发生的概率相等^[19]，即取P(N) = 1/N_max。

当土层数目为N时，土层厚度H_N存在无数种组合，但其深度之和始终等于CPT触探深度H。令h = [h_n, n = 1, 2, 3, …, N]，其中h_n = H_n / H，则h₁ + h₂ + … + h_N = 1。假设h服从Dirichlet分布，其概率密度函数为^[22]

式中：α = (α₁, α₂, α₃, …, α_N)为Dirichlet分布的参数；B(α)为归一化常数，由多元Beta函数表示；Γ(·)为伽马函数。h分布在N-1标准单纯形中，且对于所有的h_n均满足h_n ≥0。如果对于土层厚度H_N无更多先验信息时，P(H_N | N)通常采用对称Dirichlet分布^[19]，此时α中所有元素取值相等，则式（4）改写为

由于h_n = H_n / H，H_N的先验分布可以表示为

注意到，在静力触探试验过程中，上覆土层和下卧土层会对CPT探头锥尖位置处的测量值产生影响^[23]，该“影响区”范围可达（10~20）倍锥径^[24]。为合理考虑CPT触探过程中该影响区的存在，同时避免划分不符工程实际的过小“薄层”，工程师可以结合探头尺寸和工程地质报告设置最小土层厚度。根据最小土层厚度Δh和式（6），H_N的先验分布可改写为

此外，Frühwirth-Schnatter^[25]通过模拟试验发现，当Dirichlet参数α的取值为4时，正确模型后验概率远大于α = 1，同时在避免过拟合能力方面，α = 4优于其他值，故式（7）中取α = 4。除了上述先验分布，似然函数的推导也是求解贝叶斯后验分布的关键步骤，如下所述。

1.3   似然函数

当土层数目为N时，式（1）中的似然函数P($\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$| H_N, N)定量反映了CPT数据$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$提供的关于N层土体剖面的信息。假设土层剖面中不同土层之间的ξ_n相互独立，则P($\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$ | H_N, N)可以表示为

式中：P($\underline{\xi}_n$ | H_N, N)为第n个土层的似然函数。前面提到，由于岩土体具有复杂的物理力学性质，导致I_c数据剖面异常复杂，存在离散性和空间变异性。为了合理地考虑岩土体的空间变异性，需要在似然函数P($\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$ |H_N, N)建立过程中引入描述I_c空间变异性的概率模型，如Cao等^[19]采用的稳态随机场模型。对于稳态随机场，通常需要假定同一土层内I_c的均值、标准差和波动范围，且同一土层内不同测点的空间相关性由自相关函数表征。常见的自相关函数类型包括单指数型、三角函数型、余弦指数型、二阶马尔科夫型、平方指数型等^[26-27]。然而，实际土层的自相关函数往往未知，如何在划分土层剖面过程中合理选择自相关函数类型表征空间变异性仍不为所知。因此，本文采用全高斯概率模型，通过均值向量μ和协方差矩阵Σ直接表征CPT数据的空间变异性，无需假定其空间相关结构^[28]。根据I_c计算公式^[8-9]，本文采用对数正态随机场作为I_c的概率模型M_P(θ)，其中θ = [μ_n, Σ_n, n = 1, 2, 3, …, N]为模型参数，则$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$ = lnI_c为正态随机场。对于第n层土体，$\underline{\xi}_n=\ln \underline{I}_{c n} \sim N\left(\underline{\mu}_n, \mathbf\Sigma_n\right)$。图 1所示为土层划分示意图及本文采用的全高斯随机场模型。

图  1  图土层划分示意图及全高斯随机场模型^[28]

Figure  1.  Illustration of full Gaussian random field model for soil stratification^[28]

下载: 全尺寸图片 幻灯片

为了得到似然函数的解析表达式，本文采用正态-逆威沙特分布（Normal-Inverse Wishart, NIW）作为模型参数θ_n= {μ_n, Σ_n}的先验分布P(μ_n, Σ_n| H_N, N)，第n个土层的似然函数可表示为^[28-29]

式中：k_n为缩放常数，可取较小值（如0.01），以保证本文所取NIW分布为弱信息先验分布；${k'_n}$ = k_n + 1，${v_n}$ = m_n+2，${v'_n}$= ${v_n}$+1。关于式（9）中参数取值和详细说明参考相关文献[28，29]。

由上可知，基于CPT数据识别最可能土层剖面等价于计算后验分布P(H_N，N| $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$)的最大值。因此，本文采用后验分布P(H_N，N| $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$)定义目标函数，通过模拟退火优化算法对其进行优化，搜索P(H_N，N| $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$)的最大值，确定最可能土层剖面的N^*和H_N^*值，具体介绍如下。

2.   土层识别优化目标函数

根据前文所述，识别最可能土层剖面的关键在于计算后验分布P(H_N，N| $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$)的最大值。本文将计算该问题转化为优化问题，采用模拟退火算法快速确定目标函数的最大值，高效识别最可能土层剖面。模拟退火算法在给定初始温度、降温函数等参数条件下引入了与温度相关的概率性接受准则，从可行解随机出发，持续进行“备选解-判断-接受/舍弃”的迭代过程，具有跳出局部最优解，逐步逼近全局最优解的优势^[30]。此外，Goffe等^[30]对比了模拟退火算法与3种传统优化方法，结果表明模拟退火算法不仅能找到全局最优解，而且鲁棒性较强。因此，本文采用基于Metropolis准则的模拟退火算法优化土层剖面参数，计算目标函数，识别最可能土层剖面，关于模拟退化算法详见文献[30，31]。

为了计算方便，根据式（3）将后验分布P(H_N, N| $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$)展开，并将其对数化后作为模拟退火算法的目标函数：

由式（10）可知，似然函数P($\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$ | H_N, N)对计算Y(H_N, N)至关重要。根据式（8），（9），P($\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$ | H_N, N)可表示为

式中：M为CPT触探深度范围内数据点数目；σ_n则为第n层CPT数据$\underline{\xi}_n$的标准差。在给定一套数据条件下，式（11）中第一项1/${(101{\text{π }})^{M/2}}$为常数；第三项$\prod\nolimits_{n = 1}^N {\Gamma [({m_n} + 3)/2]}$/${[\Gamma (3/2)]^N}$只与第n个土层内的数据量有关，而与I_c值无关；第二项$\prod\nolimits_{n = 1}^N {\sigma _n^{ - {m_n}}}$表明每一层土中数据的离散性对土层识别结果的影响显著，即所提方法通过最大化目标函数倾向于首先识别数据离散性较小（即${\sigma _n}$较小）且较厚（即m_n较大）的土层。

在给定土层数目N值条件下，本文通过模拟退火算法搜寻目标优化函数Y(H_N, N)的最大值，然后比较不同N对应的Y(H_N, N)的最大值，确定最可能土层剖面，即N^*和H_N^*。值得注意的是，尽管本文采用模拟退化算法搜寻给定N条件下的最可能土层剖面，其它全局优化算法也可以用于本文所提方法。

3.   计算流程

所提方法的计算流程主要包括以下4个步骤：①通过岩土工程现场勘探获取一套CPT数据，计算得到I_c，并将其转换为$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$；②确定CPT勘探深度范围内最大可能土层数目N_max和最小土层厚度Δh；③在给定土层数N条件下，采用模拟退火算法确定Y(H_N, N)的最大值及其对应的土层剖面；④重复步骤③N_max次，比较不同N对应的目标函数Y(H_N, N)的最大值，确定所有候选土层模型中Y(H_N, N)的最大值，其相应的土层数目N^*和土层厚度H_N^*，即为最可能土层剖面。

4.   算例分析

4.1   场地概况

本节采用杭州某地铁线路区间勘探CPT数据说明所提方法，所采用的CPT数据如图 2所示。本文采用Robertson等^[8]方法计算I_c值，其主要计算公式如下：

式中：Q为归一化锥尖阻力，由锥尖阻力q_c（见图 2（a））计算；F为归一化摩阻比，利用q_c和侧壁摩阻力f_s（见图 2（b））计算。需要指出，通过式（12）计算I_c时涉及迭代计算，具体计算过程参考文献[8]。根据式（12）以及图 2（a），（b）所示的CPT数据（即q_c和f_s），计算得到I_c值，如图 2（c）所示。图 2（d）为根据临近钻孔绘制的土层剖面图，钻探深度范围内土层根据物性共分为5层。根据钻探取土分层信息，0~12.76 m为砂质粉土，12.76~15.80 m为粉砂夹砂质粉土，15.80~22.18 m为砂质粉土，22.18~26.38 m为淤泥质粉质黏土夹黏质粉土，26.38~34.00 m为淤泥质黏土。值得注意的是，CPT无法直接提供土样，本文所提方法仅需利用图 2（c）所示的I_c数据划分土层，图 2（d）所示钻孔剖面仅用来验证所提方法的合理性，具体如下所述。

图  2  杭州某地铁区间内一套CPT数据

Figure  2.  CPT data obtained from some section of Hangzhou metro line

下载: 全尺寸图片 幻灯片

本算例选取N_max = 10（即可能的土层数目取值范围为1到10）。由于现场原位测试所用CPT探头截面积为10 cm²，锥径为35.7 mm。根据前文所述，最小土层厚度Δh设为0.7 m。采用模拟退火算法计算每个可能N值对应的目标函数最大值，识别最可能的土层数目和土层边界。

4.2   土层剖面识别结果

图 3所示为每个N值对应的目标函数最大值。由图可知，随着N值增加，目标函数呈现先增加后降低的趋势，当N = 4时，达到最大值516.12，即最可能土层数目为4层，存在3个土层边界。图 4（a）采用红色横实线展示了3个土层界面的位置。为了便于比较，图中用黑色圆圈表示I_c数据，灰色竖虚线表示不同土类的I_c边界^[10]。可以发现，所提方法可以根据I_c的统计特征合理地划分土层剖面，同一土层内的I_c数据具有相似的力学统计特性，而不同土层间的统计特性则存在显著差异。比如，图 4（a）中第1层和第2层I_c数据的趋势性并不明显，均未表现出明显的趋势性，但第2层的空间变异性更大。如前所述，基于I_c和钻孔信息的土体分类标准有所差异。对比所识别剖面（见图 4（a））与钻孔剖面（见图 4（c））可知，虽然根据钻孔信息划分的砂质粉土层与粉砂夹砂质粉土层在物性方面不同，但是它们的I_c统计特征相似，导致所提方法将其识别为一层，如图 4（a）所示；此外，钻孔剖面第1层和第3层均为砂质粉土层，然而其I_c值显著不同，所提方法根据I_c数据识别出了深度位于15.8 m处的土层界面。

图  3  目标函数最大值

Figure  3.  Maximum values of objective optimization function

下载: 全尺寸图片 幻灯片

图  4  最可能土层剖面识别结果

Figure  4.  Most probable soil stratification

下载: 全尺寸图片 幻灯片

此外，为了进一步验证所提方法的有效性，基于所提方法识别了杭州该地铁区间内其它CPT数据的土层剖面，识别结果类似。篇幅所限，根据所提方法对这些CPT数据的土层划分结果在此不做介绍。结果表明，基于I_c的力学分层和基于钻孔的物性分层存在一定差异，尽管物性相似，土体力学响应也可能存在不同。例如，所提方法将钻孔剖面上部的砂质粉土层和粉砂夹砂质粉土层识别为一层。所提方法能够基于I_c数据的统计特性较合理地识别力学分层剖面。在本算例中，基于I_c数据识别的淤泥质粉质黏土夹黏质粉土层和淤泥质黏土层位置与钻孔剖面比较接近，如图 4（a），（c）所示。表明所提方法在该场地能够有效地识别上述两个土层，对滨海地区浅层勘探具有一定的应用价值。

4.3   对比分析

根据文献[19]所提方法划分土层剖面的结果如图 4（b）所示，其中表征I_c数据的自相关函数假定为单指数型。对比图 4（a），（b）可知，两者划分的土层剖面结果较为类似，均为4层。在计算效率方面，使用16.0 GB内存、3.20 GHz主频的Intel Core i7-8700 CPU处理器的台式计算机，采用Cao等^[19]所提方法耗时约为4.86 h，而本文所提方法仅需要1.45 min，显著提高基于CPT数据的土层识别，说明了所提方法用于识别土层剖面的高效性。

5.   方法验证及讨论

5.1   模拟数据Ⅰ

为了进一步说明所提方法，本节基于对数正态随机场模型生成已知土层界面位置的I_c数据（即模拟数据Ⅰ），模拟数据时采用单指数型自相关函数表征空间相关性。令I_cn表示第n层土中I_c的对数正态随机场，随机场模型参数θ^'={μ_n^'，σ_n^'，λ_n^'，n = 1，2，3，…，N}，其中μ_n^'，σ_n^'，λ_n^'分别为第n个土层的均值、标准差和波动范围。根据I_c对数随机场模型参数的取值范围^[19]以及Robertson等^[8]的I_c土分类表，模拟数据Ⅰ的随机场参数取值如图 5（a）所示。

本文采用矩阵分解法模拟I_c数据，数据点间隔为0.05 m，总深度为35 m，共包括5层土，其厚度分别为2，3，5，10，15 m，如图 5（a）所示。图 5（b）分别用黑色圆圈和红色实线展示了模拟数据Ⅰ的剖面和该模拟场地实际土层界面位置。和现场CPT数据（见第4节）相比，模拟场地的“实际”土层界面位置已知，用于比较和说明所提方法的有效性。需要指出，图 5（a）用于生成模拟场地的对数正态随机场模型及其参数仅用于模拟各层具有统计特性差异的I_c数据，应用所提方法识别最可能土层剖面过程时不需要这些信息。

图  5  模拟场地Ⅰ土层剖面与模拟数据

Figure  5.  Virtual site Ⅰ and simulated data

下载: 全尺寸图片 幻灯片

本节仍取最大土层数目N_max = 10，并设置最小土层厚度Δh = 0.7 m，采用模拟退火算法计算每个可能N值对应目标函数最大值。图 6为不同N值对应的目标函数最大值，当N = 5时，目标函数达到最大值，因此所识别的最可能土层剖面数目为5层，与真实模型相同。

图  6  模拟场地Ⅰ的目标函数最大值

Figure  6.  Maximum values of objective function for virtual site Ⅰ

下载: 全尺寸图片 幻灯片

使用上述相同配置的台式计算机，本文所提方法正确识别模拟数据Ⅰ的最可能土层剖面仅需2.03分钟，说明了所提方法在计算效率方面的优势。图 7为模拟数据Ⅰ的土层剖面划分结果。其中，图 7（a）为N取值不同时根据模拟退火算法搜索的土层界面。然后，比较不同土层数目对应的联合后验分布最大值（见图 6），识别最可能土层剖面，如图 7（b）所示。由图可知，当N = 2时，识别出土层界面深度为20.00 m，随着N值的增大，逐渐识别出所有真实土层界面位置；当N = 5时，识别的土层剖面与模拟场地真实模型土层剖面一致，共包括4个界面深度，如图 7（b）所示。需要指出，本文所提方法识别的土层界面是服从联合后验分布的随机样本，允许落在数据点之间。因此，当N = 5时，土层界面的识别顺序依次为20.04，10.04，5.04，2.04 m。如前文所述，所提方法首先识别I_c离散性较小（方差较小）且较厚的土层（如第5层），以上结果验证了所提方法的合理性。

图  7  模拟场地Ⅰ的土层划分结果

Figure  7.  Results of soil stratification for virtual site Ⅰ

下载: 全尺寸图片 幻灯片

5.2   模拟数据Ⅱ

本文通过模拟数据Ⅰ说明了CPT测试数据较深时，所提方法用于识别土层剖面的合理性。为了进一步说明所提方法对于浅层CPT测试数据的可靠性，本节基于对数正态随机场模型（见图 8（a））生成一套较短的I_c模拟数据（见图 8（b）），模拟过程和方法与前文相同，模拟深度为5 m，包括厚度分别为2，3 m的两个土层。由图 8（b）中灰色竖直虚线所示的I_c边界可知，模拟数据Ⅱ中两层土的土体类型不同。为了说明所提方法对于浅层CPT数据力学分层的有效性，分别采用本文所提方法和文献[19]所提方法根据模拟数据Ⅱ划分土层剖面。

图  8  模拟场地Ⅱ土层剖面与模拟数据Ⅱ

Figure  8.  Virtual site Ⅱ and simulated data Ⅱ

下载: 全尺寸图片 幻灯片

由于模拟场地Ⅱ深度较浅，仅为5 m，故本节取最大土层数目N_max = 5。本节首先采用文献[19]所提方法识别模拟场地Ⅱ土层剖面。需要指出，该方法需要近似计算不同N值对应的模型证据对数值lnP($\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\xi }$ | N)，具体算法流程可参见文献[19]，计算结果如图 9（a）所示。采用本文所提方法划分模拟场地Ⅱ的土层剖面时，首先设置最小土层厚度Δh = 0.7 m，然后通过模拟退火算法计算每个可能N值下目标函数最大值，如图 9（b）所示。

图  9  模拟场地Ⅱ的计算结果

Figure  9.  Calculated results for virtual site Ⅱ

下载: 全尺寸图片 幻灯片

由图 9（a）可知，当N的取值范围为1~5时，采用文献[19]所提方法计算的模型证据对数值逐渐减小，当N = 1时，模型证据最大，即文献[19]所提方法识别的最可能土层剖面为1层，未能根据浅层CPT数据识别出正确土层界面。由于文献[19]所提方法在CPT数据深度较浅时（如5 m）无法准确计算模型证据，存在近似误差，难以有效识别土层数目。因此，当CPT数据较短时，应谨慎使用基于模型证据的贝叶斯土层识别方法。与之相反，本文所提方法利用模拟退火优化算法避免了耗时的模型证据近似求解，快速确定了最可能土层数目为2层（即N = 2时目标函数最大，如图 9（b）所示）。本文所提方法识别的最可能土层剖面数目与模拟场地Ⅱ真实模型土层数目一致，说明了所提方法在浅层勘探的有效性。

模拟场地Ⅱ土层划分结果如图 10所示。由图 10（a）可知，根据文献[19]所提方法，虽然当N = 2时识别的土层剖面与模拟场地Ⅱ较为稳合，但该方法未能识别出正确的土层数（见图 9（a）），因此无法根据浅层CPT数据正确识别土层界面。相对而言，本文所提方法不仅能够正确识别最可能的土层数目，所识别边界位置D₁^* = 1.84 m与真实边界深度（2.00 m）基本一致（见图 10（b）），说明了所提方法在划分浅层CPT数据土层剖面时的有效性。此外，在计算效率方面，使用相同配置的台式计算机，采用文献[19]所提方法耗时约为31 min，而本文所提方法仅需要11 s。综上可知，相对于基于模型证据的贝叶斯土层识别方法，本文所提方法在浅层CPT数据分析方面具有一定优势，适用于浅层CPT数据分析（如基于滨海滩地浅层CPT数据的力学分层）。

图  10  模拟场地Ⅱ的土层划分结果

Figure  10.  Results of soil stratification for virtual site Ⅱ

下载: 全尺寸图片 幻灯片

6.   结论

本文提出了一种基于静力触探试验和变维联合后验分布的土层剖面高效优化识别方法，通过杭州某地铁区间CPT实测数据和模拟数据说明了所提方法的合理性和有效性。

（1）所提方法通过模拟退火优化算法搜索土层剖面参数（土层数目和土层厚度）联合后验分布的最大值，确定最可能土层剖面。在优化过程中，采用基于全高斯概率模型推导的似然函数计算目标函数，显著缩短了基于贝叶斯学习和I_c数据的土体力学分层方法的计算时间，提高了土层识别计算效率，便于工程应用。此外，由于所提方法在土层识别过程无需计算模型证据，不仅适合勘探深度较深的CPT数据分析，也适用于较浅层CPT数据分析（如滨海滩地勘探等）。

（2）所提方法根据I_c数据的统计相似性进行力学分层，倾向于首先识别数据离散性较小且较厚的土层。基于I_c统计特征的力学分层结果与基于钻孔取样的物性分层结果存在一定差异，物性分层结果未必能体现土层力学响应的统计相似性，因此应该根据勘探目标选择合适的土体分类标准划分土层剖面。

（3）所提方法仅根据CPT数据划分土层力学剖面，未能考虑岩土工程勘探中可能获取的其它数据。近年来，不同机器学习方法在岩土工程数据分析、反演、解译等方面得到了广泛应用。如何基于机器学习方法拓展所提方法，合理地考虑多种勘探数据划分土层，是本研究需要进一步探讨的问题。