基于吉布斯采样与压缩感知的二维非平稳CPT数据快速插值方法

朱文清 赵腾远 宋超 王宇 许领

摘 要：靜力触探试验（Cone Penetration Test，CPT）常被用于确定地下土体分层情况及层内土体的力学参数等。由于工期、工程投入、技术等条件限制，沿水平方向的CPT钻孔数目通常非常有限，有必要利用空间插值或随机模拟来估计未采样位置的CPT试验数据。提出一种有效的蒙特卡洛方法，可直接根据有限的CPT试验钻孔数据估计未采样位置的CPT数据，该方法将二维贝叶斯压缩感知框架与吉布斯采样相结合，并引入克罗内克积以提高其计算效率，然后用一系列数值及实际工程案例验证了所提方法的可靠性。结果表明：该插值方法合理，不仅能如实反映数据本身的非平稳特点，且采用序列更新技术后可显著降低时间成本，具有更强的适应能力。此外，插值结果的准确性、可靠性与已有CPT钻孔的距离成反比、与已有钻孔的数目成正比，反映出方法本身数据驱动的特点。

关键词：场地概率勘察;空间变异性;机器学习;数据驱动;马尔科夫链蒙特卡洛

中图分类号：TU413.3 文献标志码：A 文章编号：2096-6717（2022）05-0098-11

收稿日期：2021-06-02

基金项目：中央高校基本科研业务费（xjh012020046）

作者简介：朱文清（1995- ），男，主要从事地质灾害防治研究，E-mail：zhuwenqing0304@163.com。

赵腾远（通信作者），男，副教授，E-mail：tyzhao@xjtu.edu.cn。

Received：2021-06-02

Foundation item：Fundamental Research Funds for the Central Universities （No. xjh012020046）.

Author brief：ZHU Wenqing （1995- ）， main research interest： geological disaster prevention， E-mail： zhuwenqing0304@163.com.

ZHAO Tengyuan （corresponding author）， associate professor， E-mail： tyzhao@xjtu.edu.cn.

Efficient interpolation method for 2D non-stationary CPT data using Gibbs sampling and compressive sampling

ZHU Wenqing， ZHAO Tengyuan， SONG Chao， WANG Yu， XU Ling

（1. College of Geology and Environment， Xi'an University of Science and Technology， Xi'an 710054， P. R. China; 2. School of Human Settlements and Civil Engineering， Xi'an Jiaotong University， Xi'an 710049， P. R. China; 3. Department of Architecture and Civil Engineering， City University of Hong Kong， Hong Kong SAR 999077， P. R. China）

Abstract：Cone penetration test （CPT） is commonly used to determine the stratification of underground soil and  the mechanical parameters of soils in stratification. Due to time， resources and/or technical constraints， the number of CPT soundings along with a horizontal direction is generally limited. In such cases， spatial interpolation or stochastic simulation methods is a necessary choice to estimate CPT data at un-sampled locations. This paper proposes an efficient method for simulating CPT data at un-sampled locations directly from a limited number of CPT records. The approach couples the framework of 2D Bayesian compressive sensing with Gibbs sampling， where Kronecker product is introduced for facilitating its simulation efficiency. Both numerical simulations and case histories are used to illustrate the presented method.Results show that the proposed method is reasonable， which can not only reflect the non-stationary characteristics of the data， but also significantly reduce the time cost and have reasonable adaptability after using the sequential updating technique. In addition， the accuracy and reliability of interpolation are negatively and positively proportional to the distance from existing CPT soundings and the number of existing CPT soundings， which demonstrates the data-driven nature of the proposed method.

Keywords：probabilistic site investigation; spatial variability; machine learning; data-driven; Markov Chain Monte Carlo

岩土工程场地勘察的目标是通过室内或现场原位试验合理地描述地下土层界面并表征土层内土体物理、力学参数的空间变异性。与室内试验相比，原位测试方法，如静力触探试验（Cone Penetration Test，CPT）成本更加低廉、测试更加快捷。此外，静力触探试验（CPT）能在深度方向获得近乎连续性的锥尖阻力（q）和侧摩阻力（f）数据，并以此作为土体的力学响应。因此，CPT广泛应用于表征岩土场地参数、地下土体分层、砂土液化评估、确定随机场的相关函数和参数等。

尽管CPT广泛应用于岩土场地勘察并沿深度方向获得近乎连续的土体力学响应，但由于工期、项目投入及技术等条件限制，在特定的工程场地或岩土工程项目中，沿水平方向的CPT钻孔数量通常很少。为解决该问题，学者们提出了众多方法来估计未采样位置处的CPT数据。例如，Fenton提出了采用随机场的方法估计未采样位置的CPT数据，但该方法需要大量的CPT数据标定自相关函数。此外，该方法并未考虑CPT数据在水平方向的空间自相关性。Cai等虽然考虑了水平、深度方向的相关性，并利用条件随机场来估计未采样位置的CPT數据，却依然难以回避平稳随机场假定以及参数估计问题。Juang等提出了利用神经网络估计三维场地未采样位置的CPT数据。虽然该方法在数据拟合及外插上效果显著，近年来在多个领域得以应用，但由于其可解释性差，较难合理量化插值过程中产生的不确定性。量化的不确定性可以有效反映插值结果的可靠性大小。地理统计方法，如克里金法可以有效解决这一问题，但该方法仅适用于满足平稳性假设的特定土层，很难应用于地下存在多个土层且空间边界不确定的岩土工程场地。Wang等提出了在CPT钻孔数据较少时利用二维贝叶斯压缩感知理论（Bayesian Compressive Sensing，BCS）进行岩土分区。然而，由于CPT钻孔深度方向数据较多，该方法计算量大，不能高效估计未采样位置的二维CPT数据。对于具有空间变异性且边界未知的地下土层，如何正确、高效地估计未采样位置的二维非平稳CPT数据，仍然是一个重大课题。

笔者提出一种计算效率较高的无参方法，可以根据数量有限的非平稳CPT钻孔数据估计二维剖面中未采样位置的CPT数据。该方法结合压缩感知（Compressive Sensing，CS）、贝叶斯框架、吉布斯采样，并引入可快速更新计算结果的克罗内克积，以提高计算效率。与已有文献相比，该方法在估计二维CPT数据方面计算效率显著提高。此外，相比于随机场模型，该方法回避了自相关函数模型选择与参数估计、数据平稳性、高斯分布等假设;相比于神经网络，该方法可以合理确定CPT钻孔数目少带来的不确定性。值得注意的是，该方法要求不同CPT钻孔沿深度方向的力学响应数据一致。

1 二维CPT数据的贝叶斯压缩感知

1.1 二维贝叶斯压缩感知框架

压缩感知（CS）是信号处理领域中新提出的采样方法。基于信号本身的可压缩性，CS可以将随着时间或空间变化且自相关的信号从它的小部分数据中重建出来。“可压缩性”意味着信号可表示为少量基函数的加权求和。

以一组二维的归一化CPTU数据Q（见图1）说明CS的基本概念。

Q =（q-σ）/σ′且q=q+（1-a）u，其中，σ和σ′分别为土体的垂直总应力和有效应力，a和u分别为圆锥面积比和孔隙水压力。用F表示上述二维Q数据，其为N×N的二维矩阵。数学上F可表示为

F=∑Nt=1Bω（1）

式中：N=N×N;B是与F同等大小的二维矩阵，表示一系列具有不同特征的二维基函数，其构建与F本身无关;ω是与B对应的权重系数，表示B对F贡献的大小。式（1）中，ω大多数元素的值都接近于零，这表明可以使用少量CPT钻孔数据（见图2）确定ω中极少数重要元素，进而合理估计出二维CPT数据F。令表示估算的ω。重建的二维CPT数据则表示为

=∑Nt=1B（2）

由式（2）可以看出，重建的关键是根据有限数量的CPT钻孔数据估算。

1.2 二维CPT数据与的关系

令矩阵Y表示从场地中n（n>1）个CPT钻孔中获得的数据集，并用它来估计。Y与ω的关系可以通过Y和F之间的关系推导出来。如图1、图2所示，CPT数据Y是F的子集，为大小N×n的矩阵。根据式（1），F是一个与ω有关的函数，故而Y亦是ω的函数。数学上，Y表示为Y=FΨ。其中，Ψ反映了CPT沿水平方向（即图2（a）中x轴）的勘测位置。上标“T”表示转置运算。

此外，式（1）中F可以分解为矩阵形式F=BΩ（B），B与B分别代表为N×N和N×N的一维正交矩阵。F=BΩ（B）是式（1）的矩阵形式表达。将F=BΩ（B）代入Y=FΨ，可得

Y=BΩ（ΨB）=BΩA（3）

式中：A=ΨB。为运算方便，通过依次叠加（Y）的列，将Y向量化为长度为M=（n× N） 的列向量y

y=vec（Y）=[BA]vec（Ω）=Aω （4）

式中：vec（·） 表示向量化;ω=vec（Ω），是式（1）中ω的向量化表示;“”表示克罗内克积;A定义为

A=（BA）=bA…bA

bAbA（5）

有關克罗内克积的更多细节问题可参考文献[31]。根据式（4）可在贝叶斯框架下推导并估计。

1.3 的贝叶斯估计

因为贝叶斯方法可以有效处理各种不确定性，如模型不确定性与空间变异性，因此采用该方法估计D。在贝叶斯框架下，的估计信息通过其后验概率密度函数（PDF）来反映，即P（|y）。

P（|y）=P（y|）P（）P（y）（6）

式中：P（y|）为似然函数，反映了给定条件下得到y的合理性;P（）是不考虑y时的先验PDF;P（y）是归一化常数，保证P（|y） 的积分为常数1。

显然，构建式（6）中P（y|）和P（）是进行贝叶斯估计的核心，分别构建为高斯似然函数和多层先验分布。为了推导P（y|），以均值为零且方差未知的高斯随机变量ε表示y与Aω之间的残差。另外，为方便推导，将ε的方差的倒数表示为τ。假定残差之间相互独立，P（y|）可表示为

P（y|，τ）=τ/（2π）×

exp-τ（y-A）（y-A）/2（7）

式中：τ未知，可将其作为随机变量并与同步更新。因此，需构建P（，τ）=P（）P（τ）。根据Zhao等、Zhao等的研究，的先验信息模型构建如下：首先，令的每个元素均服从均值为零、方差为α的高斯分布，则P（|α）可表示为P（|α）=∏Nt=1α（2π）-1/2exp[-α（）/2]，其中，α是表示α的矢量形式;其次，假定α服从参数为1和γ/2 （γ > 0） 的逆伽玛分布，即P（α |γ） = γ/2αexp[-γ/2·α] （t = 1， 2…N）且P（α|γ）=∏Nt=1P（α|γ）;最后，假定γ服从伽马分布，即P（γ）=（b）aγexp（-bγ）/Γ（a），其中a和b为极小非负数（如a = b = 10），这对构建^2D的无信息先验分布至关重要。此外，假定τ亦服从伽玛分布，且P（τ）=dcτexp（-dτ）/Γ（c），c和d分别取值为1和10，可以使得P（τ）成为τ的无信息先验分布。

将P（|α）、P（α|γ）和P（τ）代入P（，α，γ，τ），可得到联合先验PDF

P（，α，γ，τ）=P（|α）p（α|γ）p（γ）p（τ）（8）

由于將α和γ均视为随机变量，因此在式（8）中使用了P（，α， γ， τ），而非P（， τ）。基于式（6）中的贝叶斯定理、式（7）中的似然函数以及式（8）中的先验PDF，后验PDF可表示为P（，α，γ，τ|y）=P（y|，τ）×P（，α，γ，τ）/P（y）。由于P（y）无解析表达式，故P（，α，γ，τ|y）亦无解析解。然而，笔者发现，在给定其他3个参数的条件下，（，α， γ， τ） 中剩余参数的后验概率密度函数均有解析解。例如，以α、γ、τ为条件的的后验PDF服从多元高斯分布，均值和协方差矩阵为

μ=COVAyτ

COV=（AAτ+D）（9）

式中：D是N×N的对角矩阵，且对角线上元素为D（t，t）=α。其他几个分布P（α|，τ，γ，y）、P（τ|，α，γ，y）和P（γ|，α，τ，y）分别服从广义逆高斯分布和两个伽玛分布。更多细节可以参考文献[32]，此处只将结果展示如下，其中P（α|，τ，γ，y）推导为

P（α|，τ，γ，y）=∏Nt=1exp-aα+bα2×

（α）p-1（a/b）2K（ab）（10）

式中：a=;b=γ;p=-1/2;K（·）表示参数为p的第二类修正Bessel函数。P（τ|，α，γ，y）和P（γ|，α，τ，y）可分别为

P（τ|，α，γ，y）=（d）cτexp（-dτ）/Γ（c）（11a）

P（γ|，α，τ，y）=（γ）γγexp（-γγ）/Γ（γ）（11b）

式中：c=M/2+c;d=d₀+1/2（y）（y）-（）Ay+1/2（）AA;γ=N+a和γ=b+∑Nt=1α。由于（， α， γ， τ）这些变量相互依存，难以得出的解析解。鉴于上述条件，概率密度函数的解析性即式（9）～式（11），笔者采用MCMC模拟中的吉布斯采样算法最终估计获得的统计特征。

2 的高效随机模拟

2.1 吉布斯采样（phthalic anhydride）

尽管P（|y）无解析解，但可采用吉布斯采样方法生成大量后验样本对其进行表征。吉布斯采样过程为：1）提供α、τ、和γ的初始值，将初始值代入式（9），计算的均值和协方差矩阵并由此生成一组样本。2）将随机产生的样本与已知的τ、γ值同时代入式（10），产生一组α样本。3）将随机产生的、α样本以及已知的γ，代入式（11a），τ为P（τ|，α，γ，y）中的唯一变量。4）在式（11b）中代入α、τ和，即γ为P（γ|，α，τ，y）中的唯一变量。5）收集α、τ和γ的样本，并在步骤1）中替换α、τ和γ的值，然后重复该过程，直到获得指定数量的MCMC样本。值得注意的是，吉布斯采样是马尔科夫链蒙特卡洛（Markov Chain Monte Carlo， MCMC）的一種特殊形式，与传统的MCMC方法不同，吉布斯采样可以有效处理高维参数问题。这是由于吉布斯采样需要已知待处理参数的条件概率分布，如式（9）～式（11）。相比之下，传统MCMC方法中的算法，如Metroplis-Hasting更加泛化，在未知参数维度较高的时候，由于难以找到合适的高维提议分布（proposal probability density function），进而导致生成样本的接受率过低，计算效率因此受到极大影响。

由于服从多元高斯分布，可在步骤1）中用Cholesky分解方法随机产生样本，表示为

=μ+Lz（12）

式中：COV=（L）（L）L是下三角矩阵;z是长度为N=N×N的列向量，其中每个元素都是均值为零、方差为1的标准高斯随机变量。式（12）表明，可通过独立分布的z来随机生成样本。然而，由于的协方差矩阵COV过于庞大，基于式（12）的分解无法在普通个人电脑上执行。此外，吉布斯采样流程涉及COV的反复更新与分解，这就使得直接通过式（12）生成的样本极其困难，进而导致难以高效随机生成非平稳的二维CPT数据。

2.2 COV和μ的序列更新

由式（9）可知，COV是关于式（5）中已定义的A的函数。根据式（5），计算COV的因子AA可推导为

AA=（BB）（AA）=I（AA）（13）

式中：B为单位正交矩阵。因此，I=BB为单位矩阵。将式（13）代入COV中，得到COV=[（I（τAA））+D]。根据克罗内克积的定义，I（τAA）可表示为

I（τAA）=（τAA）…00…（τAA）（14）

式中“0”是一个N×N的零元素矩阵。将式（14）代入COV，得

COV=（τAA）+D…00…（τAA）+D（15）

式中：D、D…D均为大小N×N的对角矩阵，分别表示D中N个连续且互斥的对角线元素。由式（15）可知，COV为分块对角矩阵，可按式（16）计算。

式中：COVi=[（τAA）+D]。式（16）将（N×N） × （N×N） 的巨型矩阵的逆矩阵运算简化为一系列尺寸为N×N小矩阵的逆矩阵的运算。此外，式（16）中无B项，这进一步提高了COV的分解效率。

利用A=（BA）=（B）（A）与式（16）可得到μ

μ=μμN（17）

式中：μ=τ∑Nj=1bCOVAy，为μ的第i项第N个连续元素;y为y的第j项第n个连续元素;b为位于矩阵B第j行与第i列的元素。将式（16）、式（17）代入式（12），得

式中：=μ+Lz，表示的第N个连续元素;L是对应于COV的下三角矩阵，且COV=（L）（L）;z是长度为N的列向量，表示N个独立的、标准的高斯随机变量。式（18）表明随机样本的产生转变为N组标准独立高斯随机变量z的随机实现。相比式（12），式（18）中进行的随机生成时仅涉及一系列大小为（N×N）的矩阵运算，避免了（N×N）×（N×N）巨型矩阵的求逆与分解。Xiao等、Ching等、Yang等也用到了類似方法。因此，该方法极大提高了样本生成的计算效率。

将生成的随机样本代入式（2），生成大量非稳态二维CPT数据。值得注意的是，该方法不需要对F的平稳性进行假设，因此，可直接用于估计具有空间变异性和多个未知边界的土层的二维非稳态CPT数据。虽然上述公式推导看似复杂，但通过商业软件，如MATLAB 可以轻易地将数学表达式编译为用户函数。为了进一步说明该方法的逻辑过程，这里给出了文中方法的流程图，主要包含5个步骤。步骤1～3属于方法的准备阶段，步骤4属于方法的核心，步骤5属于方法的结尾，详见图3。此外，为了便于读者重复该工作，给出了步骤4，即方法核心的伪代码，详见图4。

3 数值模拟案例

为了验证该方法，选取一组分布在3层土层内的二维CPT锥尖阻力（Q）数据，如图1所示。此例中，采用高斯平稳随机场生成器（如循环嵌入算法）在垂直方向和水平方向分别以η=0.02 m和η=0.5 m的分辨率生成每层内的Q数据。随机场模拟中涉及的参数包括Q的均值、标准差SD、垂直方向以及水平方向相关长度λ和λ。第1～3层的均值分别取为12、40和12;SD分别为2、5和2;相关长度分别为λ=2 m和λ=20 m，与文献[40]一致。在生成每层Q数据时采用单指数自相关函数

ρ=exp-2（Δz）λ+（Δx）λ（19）

式中：Δz=z-z和Δx=x-x分别表示位置（z，x）和（z，x）沿z和x方向的相对距离。利用上述随机场参数和随空间变化的岩土层边界，可获得3层内的非稳态二维Q数据，如图1所示。尽管每层中的Q数据是稳态的，但由于不同土层中Q的统计量不同，图1所示的二维数据集是非稳态的。

如图2所示，当n=5时的CPT钻孔锥尖阻力数据可作为所提出方法的输入Y （步骤1）。首先，令y=vec（Y） （步骤2）。之后，记录Q沿x方向的位置，以此构造矩阵Ψ;确定N和N，N代表Q数据点沿深度方向的数量，本例中N=500。N表示分辨率为η的插值后的CPT数据沿x方向的数量，计算公式为N=h/η+1，式中h代表场地的长度。例如，令h=127.5 m、η=0.5 m，则N=256。在MATLAB中输入使用“dctmtx”以及N、N构造一维正交矩阵B与B，且A=ΨB（步骤3）。

接下来按照图4引入克罗内克积的吉布斯采样方法，生成大量样本（步骤4）。此例中，首先随机产生了2 100个样本。然后每隔20步收集一次样本，从而保证在丢弃前100个不可信的样本后，样本之间的统计相关性较弱。如此，共生成n = 100个相互独立的随机样本。将每次产生的样本代入式（2），产生一个二维的Q样本。100组样本产生了n = 100组二维且独立的Q样本（步骤5）。图5给出了上述两个CPT样本。值得注意的是，虽然这里只生成了100个相互独立的Q样本，但其统计特征基本收敛，如每一点的均值及标准差并不随着MCMC样本的增多而发生变化。

为了说明引入克罗内克积对计算效率的提升，记录生成上述100个独立CPT样本的时间，见表1。由表1可见，使用64位Windows 10操作系统的IntelCorei7-9700 3.0 GHz CPU和16.0 GB RAM的计算机，生成100组独立二维CPT样本大约需40.6 s。然而，若不采用序列更新方法，即直接采用式（12）生成样本时，同一台计算机上大约需1 421.2 s。相比于原始方法，该算法在计算效率上提高了约35倍。随着CPT勘测钻孔数量n的增加，使用序列更新技术的计算效率提高会更为明显。

利用上述100个二维Q样本可得其均值，如图6所示。结果表明，即使在未知地下边界处，Q均值的分布也与图1所示较为一致。表明利用该方法生成的二维非稳态Q样本在统计上保留了图1中Q原始数据的非稳态模式。此外，根据生成的二维Q样本还可计算每一点的标准差，见图7。结果表明，CPT钻孔位置处的标准差非常小，说明在这些位置估计的Q可靠性和可信度較高。因为这些位置的Q数据已知，并将其作为建议方法的输入。相反，因为远离测量点的位置有效信息极少，导致远离CPT钻孔位置的样本标准差相对较大。

为了进一步验证该方法插值的合理性，比较图2（a）中4个未测量钻孔（即BH、BH、BH、BH）的Q数据。图8（a）～（c）分别以虚线绘制了该方法在这4个钻孔位置插值的Q数据。为进行比较，图8以粗线给出了BH至BH处的原始Q变化曲线。图8显示，每个子图中虚线的变化趋势与实线一致，表明利用提出的方法生成的Q样本是合理、可靠的，并较为合理地刻画了图1中CPT数据的非平稳特点。

图8 位置BH到BH处的Q均值与原始数据的对比

Fig.8 Comparison between averaged Q profile and the original one at locations BH to BH[]

值得注意的是，图8中虚线、灰线和粗实线之间存在一些差异。这是因为使用提出的方法时，输入的CPT钻孔数量较少。当n增加时，二维Q样本会更加接近CPT真实数据。

4 CPT测深数量的影响

采用n=5、15、25和50的案例探讨该方法的性能。对于每个n案例，前述方法随机生成100个独立的二维Q样本，然后利用这100个Q样本计算每个位置的Q均值及标准差，详见图9、图10。由图9可以看出，随着n的增加，每个位置的Q样本均值越来越接近于图1。当n=25时，Q样本均值与实测Q几乎相同。此外，图10显示，随着n的增加，每个位置估算的Q标准差不断减小，表明每个位置Q估算值的可靠性和置信度都有所提高。该计算结果反映了本文所提方法的数据驱动本质。

另外，随着n的增加，采用序列更新技术的计算时间会略有延长，如表1所示。但与不采用序列更新技术的方法相比，该方法所增加的计算时间几乎可以忽略。当n=15时，该方法仅需约79.6 s，而未采用序列更新技术的方法需9 h以上（见表1）。两种方法计算时间的差距表明，所提出的方法可显著提高计算效率。

值得强调的是，该方法可以较为合理地考虑CPT钻孔数量导致的不确定性，为了定量说明不确定性的大小，根据图9、图10分别计算出了n=5、15、25、50不同情况下变异系数 （δ=σ/μ×100%）的中位数，分别是11.80%、10.86%、9.36%与8.10%。此外，根据计算经验，该方法能适用的最少CPT钻孔数在4～5个左右，如果CPT钻孔更少，应谨慎使用该方法。

5 工程应用

为进一步验证方法的有效性，选取位于日本冈山河堤的某工程场地进行验证研究。该场地自上而下主要由砂土、黏土或粉土、砂土组成。其中上层砂土厚约2 m，黏土最厚约3 m，粉土厚约1～2 m，粉土下面仍主要以砂土为主。为了探明该河堤的软土空间分布，日本工程师在此进行了一系列的CPT。选取其中的CPT-201～CPT-207来验证该方法。图11给出了上述7个CPT钻孔的位置分布图，其中，CPT-204用来验证，其余CPT标准化Q数据当作输入。所有的数据均来自TC304免费数据库（http：//140.112.12.21/issmge/tc304.htm）。按照该方法，可以得到Q的二维剖面分布，其均值与标准差见图12。由图12可知，估计的Q数据离CPT钻孔位置愈接近，对应的标准差愈小，反之，则愈大。由于这是真实原位试验数据，无法得知CPT-204以外的真实数据。这里以CPT-204对应的Q数据与此位置的估计数据进行对比，见图13。图13显示，虽然估计值与试验值之间存在一定差距，但两者趋势基本一致;此外，CPT-204的大部分实测试验值落在估计值的95%置信区间里，间接说明该方法的准确性以及鲁棒性。

6 结论

基于吉布斯采样与贝叶斯压缩感知方法，提出了一种快速估计未采样位置的非平稳静力触探测试数据（CPT）的方法。該方法属于无参估计范畴，能够自动考虑静力触探数据沿水平方向的相关性，同时回避了水平、深度方向自相关函数的采用以及数据平稳性等假设。生成的CPT数据可用于解决各种岩土工程和地质工程问题，如土壤分层和分区、土壤液化潜力评估及其空间变异性表征、岩土设计参数的间接估计等。该方法为解决实际工程问题提供了一种概率工具，尤其是针对在实践中经常遇到的沿水平方向的CPT勘测钻孔数量较少的情况。最后，通过数值与实际工程案例对提出的方法进行验证，结果表明，该方法较为准确，并且采用序列更新技术后可显著提高计算效率，具有较强的鲁棒性。

参考文献：

[1] MATTHEWS M C， SIMONS N E. Site investigation： A handbook for engineers [M]. Hoboken， New Jersey， U.S： Wiley-Blackwell， 1995.

[2] 张继周， 缪林昌， 王华敬. 土性参数不确定性描述方法的探讨[J]. 岩土工程学报， 2009， 31（12）： 1936-1940.

ZHANG J Z， MIAO L C， WANG H J. Methods for characterizing variability of soil parameters [J]. Chinese Journal of Geotechnical Engineering， 2009， 31（12）： 1936-1940. （in Chinese）

[3] 刘松玉， 蔡正银. 土工测试技术发展综述[J]. 土木工程学报， 2012， 45（3）： 151-165.

LIU S Y， CAI Z Y. Review of the geotechnical testing [J]. China Civil Engineering Journal， 2012， 45（3）： 151-165. （in Chinese）

[4] 刘松玉， 吴燕开. 论我国静力触探技术（CPT）现状与发展[J]. 岩土工程学报， 2004， 26（4）： 553-556.

LIU S Y， WU Y K. On the state-of-art and development of CPT in China [J]. Chinese Journal of Geotechnical Engineering， 2004， 26（4）： 553-556. （in Chinese）

[5] 沈小克， 蔡正银， 蔡国军. 原位测试技术与工程勘察应用[J]. 土木工程学报， 2016， 49（2）： 98-120.

SHEN X K， CAI Z Y， CAI G J. Applications of in situ tests in site characterization and evaluation [J]. China Civil Engineering Journal， 2016， 49（2）： 98-120. （in Chinese）

[6] 孟高头， 张德波， 刘事莲， 等. 推广孔压静力触探技术的意义[J]. 岩土工程学报， 2000， 22（3）： 314-318.

MENG G T， ZHANG D B， LIU S L， et al. The significance of piezocone penetration test [J]. Chinese Journal of Geotechnical Engineering， 2000， 22（3）： 314-318. （in Chinese）

[7] LUNNE T， POWELL J J M， ROBERTSON P K. Cone penetration testing in geotechnical practice [M]. London， UK： Taylor & Francis： CRC Press， 2002.

[8] 曹子君， 郑硕， 李典庆， 等. 基于静力触探的土层自动划分方法与不确定性表征[J]. 岩土工程学报， 2018， 40（2）： 336-345.

CAO Z J， ZHENG S， LI D Q， et al. Probabilistic characterization of underground stratigraphy and its uncertainty based on cone penetration test [J]. Chinese Journal of Geotechnical Engineering， 2018， 40（2）： 336-345. （in Chinese）

[9] CAO Z J， ZHENG S， LI D Q， et al. Bayesian identification of soil stratigraphy based on soil behaviour type index [J]. Canadian Geotechnical Journal， 2019， 56（4）： 570-586.

[10] 刘松玉， 邹海峰， 蔡国军， 等. 基于CPTU的土分类方法在港珠澳大桥中的应用[J]. 岩土工程学报， 2017， 39（Sup2）： 1-4.

LIU S Y， ZOU H F， CAI G J， et al. Application of CPTU-based soil classification methods in Hong Kong-Zhuhai-Macao Bridge [J]. Chinese Journal of Geotechnical Engineering， 2017， 39（Sup2）： 1-4. （in Chinese）

[11] 林軍， 蔡国军， 刘松玉， 等. 基于孔压静力触探力学分层的土体边界识别方法研究[J]. 岩土力学， 2017， 38（5）： 1413-1423.

LIN J， CAI G J， LIU S Y， et al. Identification of soil layer boundaries using mechanical layered method base on piezocone penetration test data [J]. Rock and Soil Mechanics， 2017， 38（5）： 1413-1423. （in Chinese）

[12] WANG Y， FU C， HUANG K. Probabilistic assessment of liquefiable soil thickness considering spatial variability and model and parameter uncertainties [J]. Géotechnique， 2017， 67（3）： 228-241.

[13] 邹海峰， 刘松玉， 蔡国军， 等. 基于电阻率CPTU的饱和砂土液化势评价研究[J]. 岩土工程学报， 2013， 35（7）： 1280-1288.

ZOU H F， LIU S Y， CAI G J， et al. Evaluation of liquefaction potential of saturated sands based on piezocome penetration tests on resistivity [J]. Chinese Journal of Geotechnical Engineering， 2013， 35（7）： 1280-1288. （in Chinese）

[14] STUEDLEIN A W， KRAMER S L， ARDUINO P， et al. Geotechnical characterization and random field modeling of desiccated clay [J]. Journal of Geotechnical and Geoenvironmental Engineering， 2012， 138（11）： 1301-1313.

[15] CAO Z J， WANG Y. Bayesian model comparison and selection of spatial correlation functions for soil parameters [J]. Structural Safety， 2014， 49： 10-17.

[16] 鄭栋， 李典庆， 黄劲松. 基于CPTU和MASW勘察信息融合的二维土性参数剖面贝叶斯表征方法[J]. 应用基础与工程科学学报， 2021， 29（2）： 337-354.

ZHENG D， LI D Q， HUANG J S. A Bayesian characterization approach for 2D profiles of soil properties via integrating information from CPTU and MASW in site investigation [J]. Journal of Basic Science and Engineering， 2021， 29（2）： 337-354. （in Chinese）

[17] FENTON G A. Random field modeling of CPT data [J]. Journal of Geotechnical and Geoenvironmental Engineering， 1999， 125（6）： 486-498.

[18] CAI Y M， LI J H， LI X Y， et al. Estimating soil resistance at unsampled locations based on limited CPT data [J]. Bulletin of Engineering Geology and the Environment， 2019， 78（5）： 3637-3648.

[19] JUANG C H， JIANG T， CHRISTOPHER R A. Three-dimensional site characterization： Neural network approach [J]. Géotechnique， 2001， 51（9）： 799-809.

[20] 王长虹， 朱合华， 钱七虎. 克里金算法与多重分形理论在岩土参数随机场分析中的应用[J]. 岩土力学， 2014， 35（Sup2）： 386-392.

WANG C H， ZHU H H， QIAN Q H. Application of Kriging methods and multi-fractal theory to estimate of geotechnical parameters spatial distribution [J]. Rock and Soil Mechanics， 2014， 35（Sup2）： 386-392. （in Chinese）

[21] 刘志平， 何秀凤， 张淑辉. 多测度加权克里金法在高边坡变形稳定性分析中的应用[J]. 水利学报， 2009， 40（6）： 709-715.

LIU Z P， HE X F， ZHANG S H. Multi-distance measures weighted Kriging method for deformation stability analysis of steep slopes [J]. Journal of Hydraulic Engineering， 2009， 40（6）： 709-715. （in Chinese）

[22] WANG Y， HU Y， ZHAO T Y. Cone penetration test （CPT）-based subsurface soil classification and zonation in two-dimensional vertical cross section using Bayesian compressive sampling [J]. Canadian Geotechnical Journal， 2020， 57（7）： 947-958.

[23] CANDES E J， WAKIN M B. An introduction to compressive sampling [J]. IEEE Signal Processing Magazine， 2008， 25（2）： 21-30.

[24] JI S H， XUE Y， CARIN L. Bayesian compressive sensing [J]. IEEE Transactions on Signal Processing， 2008， 56（6）： 2346-2356.

[25] 赵腾远， ALADEJARE Adeyemi Emman， 王宇. 基于贝叶斯方法的模型选择以及岩石性质概率表征[J]. 武汉大学学报（工学版）， 2016， 49（5）： 740-744.

ZHAO T Y， ALADEJARE A E， WANG Y. Bayesian model selection and characterization for rock properties [J]. Engineering Journal of Wuhan University， 2016， 49（5）： 740-744. （in Chinese）

[26] 曹子君， 赵腾远， 王宇， 等. 基于贝叶斯等效样本的土体杨氏模量的统计特征确定方法[J]. 防灾减灾工程学报， 2015， 35（5）： 581-585.

CAO Z J， ZHAO T Y， WANG Y， et al. Characterization of Young's modulus of soil using Bayesian equivalent samples [J]. Journal of Disaster Prevention and Mitigation Engineering， 2015， 35（5）： 581-585. （in Chinese）

[27] CANDES E J， ROMBERG J， TAO T. Robust uncertainty principles： Exact signal reconstruction from highly incomplete frequency information [J]. IEEE Transactions on Information Theory， 2006， 52（2）： 489-509.

[28] DONOHO D L. Compressed sensing [J]. IEEE Transactions on Information Theory， 2006， 52（4）： 1289-1306.

[29] BURNS S E， MAYNE P W. Analytical cavity expansion-critical state model for piezocone dissipation in fine-grained soils [J]. Soils and Foundations， 2002， 42（2）： 131-137.

[30] ZHAO T Y， HU Y， WANG Y. Statistical interpretation of spatially varying 2D geo-data from sparse measurements using Bayesian compressive sampling [J]. Engineering Geology， 2018， 246： 162-175.

[31] PETERSEN K， PEDERSEN M. The matrix cookbook [R]. Technical University Denmark， Kongens Lyngby， Denmark， 2012.

[32] ZHAO T Y， WANG Y. Non-parametric simulation of non-stationary non-Gaussian 3D random field samples directly from sparse measurements using signal decomposition and Markov Chain Monte Carlo （MCMC） simulation [J]. Reliability Engineering & System Safety， 2020， 203： 107087.

[33] ZHAO Q B， ZHANG L Q， CICHOCKI A. Bayesian sparse tucker models for dimension reduction and tensor completion [J/OL]. Computer Science， https：//arxiv.org/abs/1505.02343.

[34] ZHAO T Y， XU L， WANG Y. Fast non-parametric simulation of 2D multi-layer cone penetration test （CPT） data without pre-stratification using Markov Chain Monte Carlo simulation [J]. Engineering Geology， 2020， 273： 105670.

[35] XIAO T， LI D Q， CAO Z J， et al. CPT-based probabilistic characterization of three-dimensional spatial variability using MLE [J]. Journal of Geotechnical and Geoenvironmental Engineering， 2018， 144（5）： 04018023.

[36] CHING J， HUANG W H， PHOON K K. 3D probabilistic site characterization by sparse Bayesian learning [J]. Journal of Engineering Mechanics， 2020， 146（12）： 04020134.

[37] YANG Z Y， CHING J. Simulation of three-dimensional random field conditioning on incomplete site data [J]. Engineering Geology， 2021， 281： 105987.

[38] MATHWORKS I. MATLAB： The language of technical computing [EB/OL]. [2021-05-21]. http：//www.mathworks.com/products/matlab/.

[39] DIETRICH C R， NEWSAM G N. A fast and exact method for multidimensional Gaussian stochastic simulations [J]. Water Resources Research， 1993， 29（8）： 2861-2869.

[40] PHOON K K， KULHAWY F H. Characterization of geotechnical variability [J]. Canadian Geotechnical Journal， 1999， 36（4）： 612-624.

（編辑 黄廷）