基于吉布斯采樣與壓縮感知的二維非平穩CPT數據快速插值方法

朱文清 趙騰遠 宋超 王宇 許領

摘 要：靜力觸探試驗（Cone Penetration Test，CPT）常被用于確定地下土體分層情況及層內土體的力學參數等。由于工期、工程投入、技術等條件限制，沿水平方向的CPT鉆孔數目通常非常有限，有必要利用空間插值或隨機模擬來估計未采樣位置的CPT試驗數據。提出一種有效的蒙特卡洛方法，可直接根據有限的CPT試驗鉆孔數據估計未采樣位置的CPT數據，該方法將二維貝葉斯壓縮感知框架與吉布斯采樣相結合，并引入克羅內克積以提高其計算效率，然后用一系列數值及實際工程案例驗證了所提方法的可靠性。結果表明：該插值方法合理，不僅能如實反映數據本身的非平穩特點，且采用序列更新技術后可顯著降低時間成本，具有更強的適應能力。此外，插值結果的準確性、可靠性與已有CPT鉆孔的距離成反比、與已有鉆孔的數目成正比，反映出方法本身數據驅動的特點。

關鍵詞：場地概率勘察;空間變異性;機器學習;數據驅動;馬爾科夫鏈蒙特卡洛

中圖分類號：TU413.3 文獻標志碼：A 文章編號：2096-6717（2022）05-0098-11

收稿日期：2021-06-02

基金項目：中央高?；究蒲袠I務費（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）是信號處理領域中新提出的采樣方法?；谛盘柋旧淼目蓧嚎s性，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.

（編輯 黃廷）