因为贝叶斯方法可以有效处理各种不确定性，如模型不确定性与空间变异性，因此采用该方法估计$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$。在贝叶斯框架下，$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$的估计信息通过其后验概率密度函数(PDF)来反映，即$ P\left( {{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}\mid {\mathit{\boldsymbol{y}}^{2{\rm{D}}}}} \right)$^[32]。

式中：$ P\left( {{\mathit{\boldsymbol{y}}^{2{\rm{D}}}}\mid {{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}} \right)$为似然函数，反映了给定$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$条件下得到y^2D的合理性；$ P\left( {{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}} \right)$是不考虑y^2D时$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$的先验PDF；P(y^2D)是归一化常数，保证$ P\left( {{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}\mid {\mathit{\boldsymbol{y}}^{2{\rm{D}}}}} \right)$的积分为常数1。

显然，构建式(6)中$ P\left( {{\mathit{\boldsymbol{y}}^{2{\rm{D}}}}\mid {{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}} \right)$和$ P\left( {{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}} \right)$是进行贝叶斯估计的核心，分别构建为高斯似然函数和多层先验分布。为了推导$ P\left( {{\mathit{\boldsymbol{y}}^{2{\rm{D}}}}\mid {{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}} \right)$，以均值为零且方差未知的高斯随机变量ε表示y^2D与Aω^2D之间的残差。另外，为方便推导，将ε的方差的倒数表示为τ。假定残差之间相互独立，$ P\left( {{\mathit{\boldsymbol{y}}^{2{\rm{D}}}}\mid {{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}} \right)$可表示为

式中：τ未知，可将其作为随机变量并与$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$同步更新。因此，需构建P($ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$, τ)=P($ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$)P(τ)。根据Zhao等^[33]、Zhao等^[34]的研究，$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$的先验信息模型构建如下：首先，令$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$的每个元素均服从均值为零、方差为α_t^-1的高斯分布，则$ P\left( {{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}\mid \mathit{\boldsymbol{\alpha }}} \right)$可表示为$ P\left( {{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}\mid \mathit{\boldsymbol{\alpha }}} \right) = \prod\limits_{t = 1}^N {\alpha _t^{1/2}} (2{\rm{ \mathsf{ π} }} ) - 1/2\exp \left[ { - {\alpha _t}{{\left( {\mathit{\boldsymbol{\widehat \omega }}_t^{2{\rm{D}}}} \right)}^2}/2} \right]$，其中，α是表示α_t的矢量形式；其次，假定α_t服从参数为1和γ/2 (γ>0) 的逆伽玛分布，即P(α_t|γ) = γ/2α_t^-2exp[-γ/2·α_t^-1] (t = 1, 2…N)且$ P(\mathit{\boldsymbol{\alpha }}\mid \gamma ) = \prod\limits_{t = 1}^N P \left( {{\alpha _t}\mid \gamma } \right)$；最后，假定γ服从伽马分布，即P(γ)=(b₀)a₀γ^a₀-1exp(-b₀γ)/Γ(a₀)，其中a₀和b₀为极小非负数(如a₀ = b₀ = 10^-4)，这对构建$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$的无信息先验分布至关重要。此外，假定τ亦服从伽玛分布，且P(τ)=d₀c₀τ^c₀-1exp(-d₀τ)/Γ(c₀)，c₀和d₀分别取值为1和10^-4，可以使得P(τ)成为τ的无信息先验分布。

将$ P\left( {{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}\mid \mathit{\boldsymbol{\alpha }}} \right)$、P(α|γ)和P(τ)代入P($ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$, α, γ, τ)，可得到联合先验PDF

由于将α和γ均视为随机变量，因此在式(8)中使用了$ P\left( {{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}{\rm{, }}\mathit{\boldsymbol{\alpha }}{\rm{, }}\gamma {\rm{, }}\tau } \right)$, 而非$ P\left( {{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}{\rm{, }}\tau } \right)$。基于式(6)中的贝叶斯定理、式(7)中的似然函数以及式(8)中的先验PDF，后验PDF可表示为$ P\left( {{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}, \mathit{\boldsymbol{\alpha }}, \gamma , \tau |} \right.\mathit{\boldsymbol{y}}) = P\left( {{\mathit{\boldsymbol{y}}^{2{\rm{D}}}}\mid {{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}, \tau } \right) \times P\left( {{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}, \mathit{\boldsymbol{\alpha }}, \gamma , \tau } \right)/P\left( {{\mathit{\boldsymbol{y}}^{2{\rm{D}}}}} \right)$。由于P(y^2D)无解析表达式，故P($ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$, α, γ, τ|y)亦无解析解。然而，笔者发现，在给定其他3个参数的条件下，($ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$, α, γ, τ) 中剩余参数的后验概率密度函数均有解析解。例如，以α、γ、τ为条件的$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$的后验PDF服从多元高斯分布，均值和协方差矩阵为

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

式中：a_t=$ \mathit{\boldsymbol{\hat \omega }}_{t}^{2{\rm{D}}}$；b_t=γ；p=-1/2；K_p(·)表示参数为p的第二类修正Bessel函数。P(τ|$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$, α, γ, y²)和P(γ|$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$, α, τ, y^2D)可分别为

式中：c_n=M/2+c₀；d_n=d₀+1/2(y^2D)^T(y^2D)-($ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$)^TA^T y^2D+1/2($ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$)^TA^TA$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$；γ_a=N+a₀和γ_b=b₀+$ \sum\limits_{t = 1}^N {\alpha _t^{ - 1}} $。由于($ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$, α, γ, τ)这些变量相互依存，难以得出$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$的解析解。鉴于上述条件，概率密度函数的解析性即式(9)~式(11)，笔者采用MCMC模拟中的吉布斯采样算法最终估计获得$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$的统计特征。

2 $ {{\hat \omega }}_{}^{2{\rm{D}}}$的高效随机模拟
跳转到：本文顶部跳转到：本文顶部　1 二维CPT数据的贝叶斯压缩感知　    1.1 二维贝叶斯压缩感知框架　    1.2 二维CPT数据与$ \hat \omega _t^{2{\rm{D}}}$的关系　    1.3 $ \hat \omega _{}^{2{\rm{D}}}$的贝叶斯估计　2 $ {{\hat \omega }}_{}^{2{\rm{D}}}$的高效随机模拟　    2.1 吉布斯采样(phthalic anhydride)　    2.2 $ {{{\rm{CO}}}}{{{{\rm{V}}}}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$和$ {{{\mu }}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$的序列更新　3 数值模拟案例　4 CPT测深数量的影响　5 工程应用　6 结论　参考文献

2.1 吉布斯采样(phthalic anhydride)
跳转到：本文顶部跳转到：本文顶部　1 二维CPT数据的贝叶斯压缩感知　    1.1 二维贝叶斯压缩感知框架　    1.2 二维CPT数据与$ \hat \omega _t^{2{\rm{D}}}$的关系　    1.3 $ \hat \omega _{}^{2{\rm{D}}}$的贝叶斯估计　2 $ {{\hat \omega }}_{}^{2{\rm{D}}}$的高效随机模拟　    2.1 吉布斯采样(phthalic anhydride)　    2.2 $ {{{\rm{CO}}}}{{{{\rm{V}}}}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$和$ {{{\mu }}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$的序列更新　3 数值模拟案例　4 CPT测深数量的影响　5 工程应用　6 结论　参考文献

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

由于$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$服从多元高斯分布，可在步骤1)中用Cholesky分解方法随机产生$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$样本，表示为

$ \hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}=\boldsymbol{\mu}_{\hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}}+\boldsymbol{L}_{\hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}} \boldsymbol{z} $

(12)

式中：$ {\bf{CO}}{{\bf{V}}_{{{\widehat {\bf{ \pmb{\mathsf{ ω}} }}}^{2{\rm{D}}}}}} = \left( {{\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}}} \right){\left( {{\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}}} \right)^{\rm{T}}}{\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}}$是下三角矩阵；z是长度为N=N₁×N₂的列向量，其中每个元素都是均值为零、方差为1的标准高斯随机变量。式(12)表明，可通过独立分布的z来随机生成$ {\mathit{\boldsymbol{\widehat \omega }}^{2{\rm{D}}}}$样本。然而，由于$ {\mathit{\boldsymbol{\widehat \omega }}^{2{\rm{D}}}}$的协方差矩阵$ {\bf{CO}}{{\bf{V}}_{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}}$过于庞大，基于式(12)的分解无法在普通个人电脑上执行。此外，吉布斯采样流程涉及$ {\bf{CO}}{{\bf{V}}_{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}}$的反复更新与分解，这就使得直接通过式(12)生成$ {\mathit{\boldsymbol{\widehat \omega }}^{2{\rm{D}}}}$的样本极其困难，进而导致难以高效随机生成非平稳的二维CPT数据。

2.2 $ {{{\rm{CO}}}}{{{{\rm{V}}}}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$和$ {{{\mu }}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$的序列更新
跳转到：本文顶部跳转到：本文顶部　1 二维CPT数据的贝叶斯压缩感知　    1.1 二维贝叶斯压缩感知框架　    1.2 二维CPT数据与$ \hat \omega _t^{2{\rm{D}}}$的关系　    1.3 $ \hat \omega _{}^{2{\rm{D}}}$的贝叶斯估计　2 $ {{\hat \omega }}_{}^{2{\rm{D}}}$的高效随机模拟　    2.1 吉布斯采样(phthalic anhydride)　    2.2 $ {{{\rm{CO}}}}{{{{\rm{V}}}}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$和$ {{{\mu }}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$的序列更新　3 数值模拟案例　4 CPT测深数量的影响　5 工程应用　6 结论　参考文献

由式(9)可知，$ {\bf{CO}}{{\bf{V}}_{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}}$是关于式(5)中已定义的A的函数。根据式(5)，计算$ {\bf{CO}}{{\bf{V}}_{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}}$的因子A^TA可推导为

$ \boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}=\left(\boldsymbol{B}_{z}^{\mathrm{T}} \boldsymbol{B}_{z}\right) \otimes\left(\boldsymbol{A}_{x}^{\mathrm{T}} \boldsymbol{A}_{x}\right)=\boldsymbol{I}_{z} \otimes\left(\boldsymbol{A}_{x}^{\mathrm{T}} \boldsymbol{A}_{x}\right) $

(13)

式中：B_z为单位正交矩阵。因此，I_z= B_z^TB_z为单位矩阵。将式(13)代入$ {\bf{CO}}{{\bf{V}}_{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}}$中，得到$ {\bf{CO}}{{\bf{V}}_{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}}$=[(I_z⊗(τA_x^TA_x))+ D^α]^-1。根据克罗内克积的定义，I_z⊗(τA_x^TA_x)可表示为

$ \boldsymbol{I}_{z} \otimes\left(\tau \boldsymbol{A}_{x}^{\mathrm{T}} \boldsymbol{A}_{x}\right)=\left[\begin{array}{ccc} \left(\tau \boldsymbol{A}_{x}^{\mathrm{T}} \boldsymbol{A}_{x}\right) & \cdots & {\bf{0}} \\ \vdots & \ddots & \vdots \\ {\bf{0}} & \cdots & \left(\tau \boldsymbol{A}_{x}^{\mathrm{T}} \boldsymbol{A}_{x}\right) \end{array}\right] $

(14)

式中“0”是一个N_x×N_x的零元素矩阵。将式(14)代入$ {\bf{CO}}{{\bf{V}}_{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}}$，得

$ {\bf{C O V}}_{\hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}}=\left[\begin{array}{ccc} \left(\tau \boldsymbol{A}_{x}^{\mathrm{T}} \boldsymbol{A}_{x}\right)+\boldsymbol{D}_{1}^{\alpha} & \cdots & {\bf{0}} \\ \vdots & \ddots & \vdots \\ {\bf{0}} & \cdots & \left(\tau \boldsymbol{A}_{x}^{\mathrm{T}} \boldsymbol{A}_{x}\right)+\boldsymbol{D}_{N_{1}}^{\alpha} \end{array}\right]^{-1} $

(15)

式中：D₁^α、D₂^α… D_N1^α均为大小N₂×N₂的对角矩阵，分别表示D^α中N₂个连续且互斥的对角线元素。由式(15)可知，$ {\bf{CO}}{{\bf{V}}_{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}}$为分块对角矩阵，可按式(16)计算。

$ {\bf{COV}}_{\hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}}=\left[\begin{array}{ccc} {\bf{COV}}_{\hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}}^{1} & & {\bf{0}} \\ \vdots & \ddots & \vdots \\ {\bf{0}} & \cdots & {\bf{C O V}}_{\hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}}^{N_{1}} \end{array}\right] $

(16)

式中：$ {\bf{COV}}_{{{\mathit{\boldsymbol{\widehat \omega }}}^{{\rm{2D}}}}}^i = {\left[ {\left( {\tau \mathit{\boldsymbol{A}}_x^{\rm{T}}{\mathit{\boldsymbol{A}}_x}} \right) + \mathit{\boldsymbol{D}}_i^\alpha } \right]^{ - 1}}$。式(16)将(N₁×N₂) × (N₁×N₂) 的巨型矩阵的逆矩阵运算简化为一系列尺寸为N₂×N₂小矩阵的逆矩阵的运算。此外，式(16)中无B_z项，这进一步提高了$ {\bf{CO}}{{\bf{V}}_{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}}$的分解效率。

利用A^T=(B_z^1D⊗A_x)^T=(B_z^1D)^T⊗(A_x)^T与式(16)可得到$ {\mathit{\boldsymbol{\mu }}_{\mathit{\boldsymbol{\hat \omega }}\mathit{D}_{}^{2{\rm{D}}}}}$

$ \boldsymbol{\mu}_{\hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}}=\left[\begin{array}{c} \boldsymbol{\mu}_{\hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}}^{1} \\ \vdots \\ \boldsymbol{\mu}_{\hat{\boldsymbol{\omega}}{}^{\mathrm{2D}}}^{N_{1}} \end{array}\right] $

(17)

式中：$ \mathit{\boldsymbol{\mu }}_{\mathit{\boldsymbol{\hat \omega }}^{2{\rm{D}}}}^i = \tau \sum\limits_{j = 1}^{{N_1}} {{b_{j, i}}} {\bf{COV}}_{\mathit{\boldsymbol{\hat \omega }}^{2{\rm{D}}}}^{i }\mathit{\boldsymbol{A}}_x^{\rm{T}}\mathit{\boldsymbol{y}}_j^{2{\rm{D}}}$，为$ {\mathit{\boldsymbol{\mu }}_{\mathit{\hat \omega }_{}^{2{\rm{D}}}}}$的第i项第N_x个连续元素；y_j^2D为y的第j项第n_b个连续元素；b_{j, i}为位于矩阵B₁第j行与第i列的元素。将式(16)、式(17)代入式(12), 得

$ \hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}=\left[\begin{array}{c} \hat{\boldsymbol{\omega}}{}_{1}^{2 \mathrm{D}} \\ \vdots \\ \hat{\boldsymbol{\omega}}{}_{N_{1}}^{2 \mathrm{D}} \end{array}\right]=\left[\begin{array}{c} \boldsymbol{\mu}_{\hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}}^{1} \\ \vdots \\ \boldsymbol{\mu}_{\hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}}^{N_{1}} \end{array}\right]+\left[\begin{array}{ccc} \boldsymbol{L}_{\hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}}^{1} & & {\bf{0}} \\ \vdots & \ddots & \vdots \\ {\bf{0}} & \cdots & \boldsymbol{L}_{\hat{\boldsymbol{\omega}}{}^{2 \mathrm{D}}}^{N_{1}} \end{array}\right]\left[\begin{array}{c} \boldsymbol{z}_{1} \\ \vdots \\ \boldsymbol{z}_{{N}_{1}} \end{array}\right] $

(18)

式中：$ \mathit{\boldsymbol{\widehat \omega }}_i^{2{\rm{D}}} = \mathit{\boldsymbol{\mu }}_{{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}}^i + \mathit{\boldsymbol{L}}_{{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}}^i{\mathit{\boldsymbol{z}}_i}$，表示$ {\mathit{\boldsymbol{\widehat \omega }}^{2{\rm{D}}}}$的第N₂个连续元素；$ \mathit{\boldsymbol{L}}_{{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}}^i$是对应于$ {\bf{COV}}_{{{\mathit{\boldsymbol{\widehat \omega }}}^{{\rm{2D}}}}}^i$的下三角矩阵，且$ {\bf{COV}}_{{{\mathit{\boldsymbol{\widehat \omega }}}^{{\rm{2D}}}}}^i = \left( {\mathit{\boldsymbol{L}}_{{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}}^i} \right){\left( {\mathit{\boldsymbol{L}}_{{{\mathit{\boldsymbol{\widehat \omega }}}^{2{\rm{D}}}}}^i} \right)^{\rm{T}}}$；z_i是长度为N₂的列向量，表示N_x个独立的、标准的高斯随机变量。式(18)表明$ \mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}$随机样本的产生转变为N₁组标准独立高斯随机变量z_i的随机实现。相比式(12)，式(18)中进行$ {\mathit{\boldsymbol{\widehat \omega }}^{2{\rm{D}}}}$的随机生成时仅涉及一系列大小为(N_x×N_x)的矩阵运算，避免了(N_z×N_x)×(N_z×N_x) 巨型矩阵的求逆与分解。Xiao等^[35]、Ching等^[36]、Yang等^[37]也用到了类似方法。因此，该方法极大提高了$ _{\mathit{\boldsymbol{\hat \omega }}_{}^{2{\rm{D}}}}$样本生成的计算效率。

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

图 3

图 3 生成二维CPT数据的流程图 Fig. 3 A flow chart for simulating 2D CPT data

图 4

图 4 引入克罗内克积吉布斯采样的伪代码 Fig. 4 Pseudo code for fast simulation of $ {\hat \omega }_{}^{2{\rm{D}}}$ using the Gibbs sampling method with Kronecker product

3 数值模拟案例
跳转到：本文顶部跳转到：本文顶部　1 二维CPT数据的贝叶斯压缩感知　    1.1 二维贝叶斯压缩感知框架　    1.2 二维CPT数据与$ \hat \omega _t^{2{\rm{D}}}$的关系　    1.3 $ \hat \omega _{}^{2{\rm{D}}}$的贝叶斯估计　2 $ {{\hat \omega }}_{}^{2{\rm{D}}}$的高效随机模拟　    2.1 吉布斯采样(phthalic anhydride)　    2.2 $ {{{\rm{CO}}}}{{{{\rm{V}}}}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$和$ {{{\mu }}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$的序列更新　3 数值模拟案例　4 CPT测深数量的影响　5 工程应用　6 结论　参考文献

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

$ \rho=\exp \left(-2 \sqrt{\frac{(\Delta z)^{2}}{\lambda_{z}^{2}}+\frac{(\Delta x)^{2}}{\lambda_{x}^{2}}}\right) $

(19)

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

如图 2所示，当n_b=5时的CPT钻孔锥尖阻力数据可作为所提出方法的输入Y (步骤1)。首先，令y =vec(Y^T) (步骤2)。之后，记录Q_t沿x方向的位置，以此构造矩阵Ψ；确定N_z和N_x，N_z代表Q_t数据点沿深度方向的数量，本例中N_z=500。N_x表示分辨率为η_x的插值后的CPT数据沿x₂方向的数量，计算公式为N_x=h_x/η_x+1，式中h_x代表场地的长度。例如，令h_x=127.5 m、η_x=0.5 m，则N_x=256。在MATLAB^[37]中输入使用“dctmtx”以及N_z、N_x构造一维正交矩阵B_z^1D与B_x^1D，且A_x= ΨB_x^1D(步骤3)。

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

图 5

图 5 二维CPT数据模拟结果示例 Fig. 5 Two examples of simulated 2D CPT data

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

表 1

表 1 不同CPT钻孔数目下，使用序列更新技术(A)和不使用序列更新技术(B)生成100组二维Qt样本的计算时间 Table 1 Computational time for simulating 100 sets of independent 2D Qt samples using the proposed method (method A) and that without sequential updating technique (method B) under different number nb of CPT soundings 方法 计算时间/s nb = 5 nb = 15 nb = 25 nb = 50 A 40.6 79.6 146.2 301.3 B 1 421.2 34 073.9 137 056.5 *无结果 *注：因运算所需内存过大，方法B在笔者电脑中无法运行。

表 1 不同CPT钻孔数目下，使用序列更新技术(A)和不使用序列更新技术(B)生成100组二维Qt样本的计算时间 Table 1 Computational time for simulating 100 sets of independent 2D Qt samples using the proposed method (method A) and that without sequential updating technique (method B) under different number nb of CPT soundings

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

图 6

图 6 MCMC模拟Qt数据的均值 Fig. 6 Averaged Qt from the simulated MCMC samples

图 7

图 7 MCMC模拟Qt数据的标准差 Fig. 7 SD of Qt from the simulated MCMC samples

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

图 8

图 8 位置BH1到BH4处的Qt均值与原始数据的对比 Fig. 8 Comparison between averaged Qt profile and the original one at locations BH1 to BH4

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

4 CPT测深数量的影响
跳转到：本文顶部跳转到：本文顶部　1 二维CPT数据的贝叶斯压缩感知　    1.1 二维贝叶斯压缩感知框架　    1.2 二维CPT数据与$ \hat \omega _t^{2{\rm{D}}}$的关系　    1.3 $ \hat \omega _{}^{2{\rm{D}}}$的贝叶斯估计　2 $ {{\hat \omega }}_{}^{2{\rm{D}}}$的高效随机模拟　    2.1 吉布斯采样(phthalic anhydride)　    2.2 $ {{{\rm{CO}}}}{{{{\rm{V}}}}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$和$ {{{\mu }}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$的序列更新　3 数值模拟案例　4 CPT测深数量的影响　5 工程应用　6 结论　参考文献

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

图 9

图 9 nb对二维Qt样本均值的影响 Fig. 9 Effect of nb on 2D Qt averaged from the MCMC samples

图 10

图 10 nb对MCMC样本二维Qt数据标准差的影响 Fig. 10 Effect of nb on SD of 2D Qt estimated from the MCMC samples

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

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

5 工程应用
跳转到：本文顶部跳转到：本文顶部　1 二维CPT数据的贝叶斯压缩感知　    1.1 二维贝叶斯压缩感知框架　    1.2 二维CPT数据与$ \hat \omega _t^{2{\rm{D}}}$的关系　    1.3 $ \hat \omega _{}^{2{\rm{D}}}$的贝叶斯估计　2 $ {{\hat \omega }}_{}^{2{\rm{D}}}$的高效随机模拟　    2.1 吉布斯采样(phthalic anhydride)　    2.2 $ {{{\rm{CO}}}}{{{{\rm{V}}}}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$和$ {{{\mu }}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$的序列更新　3 数值模拟案例　4 CPT测深数量的影响　5 工程应用　6 结论　参考文献

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

图 11

图 11 CPT钻孔201到207分布图 Fig. 11 Spatial distribution of CPT soundings ranging from CPT-201 to CPT-207

图 12

图 12 日本冈山河堤Qt剖面估计结果 Fig. 12 Estimated Qt profile of the Okayama Riverbank in Japan

图 13

图 13 CPT-204估计值与实测值的对比 Fig. 13 Comparison between predicted and measured Qt at CPT-204

6 结论
跳转到：本文顶部跳转到：本文顶部　1 二维CPT数据的贝叶斯压缩感知　    1.1 二维贝叶斯压缩感知框架　    1.2 二维CPT数据与$ \hat \omega _t^{2{\rm{D}}}$的关系　    1.3 $ \hat \omega _{}^{2{\rm{D}}}$的贝叶斯估计　2 $ {{\hat \omega }}_{}^{2{\rm{D}}}$的高效随机模拟　    2.1 吉布斯采样(phthalic anhydride)　    2.2 $ {{{\rm{CO}}}}{{{{\rm{V}}}}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$和$ {{{\mu }}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$的序列更新　3 数值模拟案例　4 CPT测深数量的影响　5 工程应用　6 结论　参考文献

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

参考文献跳转到：本文顶部跳转到：本文顶部　1 二维CPT数据的贝叶斯压缩感知　    1.1 二维贝叶斯压缩感知框架　    1.2 二维CPT数据与$ \hat \omega _t^{2{\rm{D}}}$的关系　    1.3 $ \hat \omega _{}^{2{\rm{D}}}$的贝叶斯估计　2 $ {{\hat \omega }}_{}^{2{\rm{D}}}$的高效随机模拟　    2.1 吉布斯采样(phthalic anhydride)　    2.2 $ {{{\rm{CO}}}}{{{{\rm{V}}}}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$和$ {{{\mu }}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$的序列更新　3 数值模拟案例　4 CPT测深数量的影响　5 工程应用　6 结论　参考文献跳转到：本文顶部跳转到：本文顶部　1 二维CPT数据的贝叶斯压缩感知　    1.1 二维贝叶斯压缩感知框架　    1.2 二维CPT数据与$ \hat \omega _t^{2{\rm{D}}}$的关系　    1.3 $ \hat \omega _{}^{2{\rm{D}}}$的贝叶斯估计　2 $ {{\hat \omega }}_{}^{2{\rm{D}}}$的高效随机模拟　    2.1 吉布斯采样(phthalic anhydride)　    2.2 $ {{{\rm{CO}}}}{{{{\rm{V}}}}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$和$ {{{\mu }}_{{{\hat \omega }}_{}^{2{\rm{D}}}}}$的序列更新　3 数值模拟案例　4 CPT测深数量的影响　5 工程应用　6 结论　参考文献

[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) DOI:10.3321/j.issn:1000-4548.2009.12.021
[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) DOI:10.3321/j.issn:1000-4548.2004.04.025
[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) DOI:10.3321/j.issn:1000-4548.2000.03.010
[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. DOI:10.1139/cgj-2017-0714
[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. DOI:10.1680/jgeot.15.P.219
[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. DOI:10.1061/(ASCE)GT.1943-5606.0000723
[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. DOI:10.1016/j.strusafe.2013.06.003
[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. DOI:10.1061/(ASCE)1090-0241(1999)125:6(486)
[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. DOI:10.1007/s10064-018-1318-2
[19]	JUANG C H, JIANG T, CHRISTOPHER R A. Three-dimensional site characterization: Neural network approach[J]. Géotechnique, 2001, 51(9): 799-809. DOI:10.1680/geot.2001.51.9.799
[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) DOI:10.3321/j.issn:0559-9350.2009.06.010
[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. DOI:10.1139/cgj-2019-0131
[23]	CANDES E J, WAKIN M B. An introduction to compressive sampling[J]. IEEE Signal Processing Magazine, 2008, 25(2): 21-30. DOI:10.1109/MSP.2007.914731
[24]	JI S H, XUE Y, CARIN L. Bayesian compressive sensing[J]. IEEE Transactions on Signal Processing, 2008, 56(6): 2346-2356. DOI:10.1109/TSP.2007.914345
[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. DOI:10.1109/TIT.2005.862083
[28]	DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306. DOI:10.1109/TIT.2006.871582
[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. DOI:10.3208/sandf.42.2_131
[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. DOI:10.1016/j.enggeo.2018.09.022
[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. DOI:10.1016/j.enggeo.2020.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. DOI:10.1061/(ASCE)GT.1943-5606.0001875
[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. DOI:10.1061/(ASCE)EM.1943-7889.0001859
[37]	YANG Z Y, CHING J. Simulation of three-dimensional random field conditioning on incomplete site data[J]. Engineering Geology, 2021, 281: 105987. DOI:10.1016/j.enggeo.2020.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. DOI:10.1029/93WR01070
[40]	PHOON K K, KULHAWY F H. Characterization of geotechnical variability[J]. Canadian Geotechnical Journal, 1999, 36(4): 612-624. DOI:10.1139/t99-038