基于吉布斯采樣與壓縮感知的二維非平穩(wěn)CPT數(shù)據(jù)快速插值方法

朱文清 趙騰遠(yuǎn) 宋超 王宇 許領(lǐng)

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

關(guān)鍵詞：場地概率勘察;空間變異性;機(jī)器學(xué)習(xí);數(shù)據(jù)驅(qū)動;馬爾科夫鏈蒙特卡洛

中圖分類號：TU413.3 文獻(xiàn)標(biāo)志碼：A 文章編號：2096-6717（2022）05-0098-11

收稿日期：2021-06-02

基金項(xiàng)目：中央高?；究蒲袠I(yè)務(wù)費(fèi)（xjh012020046）

作者簡介：朱文清（1995- ），男，主要從事地質(zhì)災(zāi)害防治研究，E-mail：zhuwenqing0304@163.com。

趙騰遠(yuǎn)（通信作者），男，副教授，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

巖土工程場地勘察的目標(biāo)是通過室內(nèi)或現(xiàn)場原位試驗(yàn)合理地描述地下土層界面并表征土層內(nèi)土體物理、力學(xué)參數(shù)的空間變異性。與室內(nèi)試驗(yàn)相比，原位測試方法，如靜力觸探試驗(yàn)（Cone Penetration Test，CPT）成本更加低廉、測試更加快捷。此外，靜力觸探試驗(yàn)（CPT）能在深度方向獲得近乎連續(xù)性的錐尖阻力（q）和側(cè)摩阻力（f）數(shù)據(jù)，并以此作為土體的力學(xué)響應(yīng)。因此，CPT廣泛應(yīng)用于表征巖土場地參數(shù)、地下土體分層、砂土液化評估、確定隨機(jī)場的相關(guān)函數(shù)和參數(shù)等。

盡管CPT廣泛應(yīng)用于巖土場地勘察并沿深度方向獲得近乎連續(xù)的土體力學(xué)響應(yīng)，但由于工期、項(xiàng)目投入及技術(shù)等條件限制，在特定的工程場地或巖土工程項(xiàng)目中，沿水平方向的CPT鉆孔數(shù)量通常很少。為解決該問題，學(xué)者們提出了眾多方法來估計未采樣位置處的CPT數(shù)據(jù)。例如，F(xiàn)enton提出了采用隨機(jī)場的方法估計未采樣位置的CPT數(shù)據(jù)，但該方法需要大量的CPT數(shù)據(jù)標(biāo)定自相關(guān)函數(shù)。此外，該方法并未考慮CPT數(shù)據(jù)在水平方向的空間自相關(guān)性。Cai等雖然考慮了水平、深度方向的相關(guān)性，并利用條件隨機(jī)場來估計未采樣位置的CPT數(shù)據(jù)，卻依然難以回避平穩(wěn)隨機(jī)場假定以及參數(shù)估計問題。Juang等提出了利用神經(jīng)網(wǎng)絡(luò)估計三維場地未采樣位置的CPT數(shù)據(jù)。雖然該方法在數(shù)據(jù)擬合及外插上效果顯著，近年來在多個領(lǐng)域得以應(yīng)用，但由于其可解釋性差，較難合理量化插值過程中產(chǎn)生的不確定性。量化的不確定性可以有效反映插值結(jié)果的可靠性大小。地理統(tǒng)計方法，如克里金法可以有效解決這一問題，但該方法僅適用于滿足平穩(wěn)性假設(shè)的特定土層，很難應(yīng)用于地下存在多個土層且空間邊界不確定的巖土工程場地。Wang等提出了在CPT鉆孔數(shù)據(jù)較少時利用二維貝葉斯壓縮感知理論（Bayesian Compressive Sensing，BCS）進(jìn)行巖土分區(qū)。然而，由于CPT鉆孔深度方向數(shù)據(jù)較多，該方法計算量大，不能高效估計未采樣位置的二維CPT數(shù)據(jù)。對于具有空間變異性且邊界未知的地下土層，如何正確、高效地估計未采樣位置的二維非平穩(wěn)CPT數(shù)據(jù)，仍然是一個重大課題。

筆者提出一種計算效率較高的無參方法，可以根據(jù)數(shù)量有限的非平穩(wěn)CPT鉆孔數(shù)據(jù)估計二維剖面中未采樣位置的CPT數(shù)據(jù)。該方法結(jié)合壓縮感知（Compressive Sensing，CS）、貝葉斯框架、吉布斯采樣，并引入可快速更新計算結(jié)果的克羅內(nèi)克積，以提高計算效率。與已有文獻(xiàn)相比，該方法在估計二維CPT數(shù)據(jù)方面計算效率顯著提高。此外，相比于隨機(jī)場模型，該方法回避了自相關(guān)函數(shù)模型選擇與參數(shù)估計、數(shù)據(jù)平穩(wěn)性、高斯分布等假設(shè);相比于神經(jīng)網(wǎng)絡(luò)，該方法可以合理確定CPT鉆孔數(shù)目少帶來的不確定性。值得注意的是，該方法要求不同CPT鉆孔沿深度方向的力學(xué)響應(yīng)數(shù)據(jù)一致。

1 二維CPT數(shù)據(jù)的貝葉斯壓縮感知

1.1 二維貝葉斯壓縮感知框架

壓縮感知（CS）是信號處理領(lǐng)域中新提出的采樣方法?；谛盘柋旧淼目蓧嚎s性，CS可以將隨著時間或空間變化且自相關(guān)的信號從它的小部分?jǐn)?shù)據(jù)中重建出來?！翱蓧嚎s性”意味著信號可表示為少量基函數(shù)的加權(quán)求和。

以一組二維的歸一化CPTU數(shù)據(jù)Q（見圖1）說明CS的基本概念。

Q =（q-σ）/σ′且q=q+（1-a）u，其中，σ和σ′分別為土體的垂直總應(yīng)力和有效應(yīng)力，a和u分別為圓錐面積比和孔隙水壓力。用F表示上述二維Q數(shù)據(jù)，其為N×N的二維矩陣。數(shù)學(xué)上F可表示為

F=∑Nt=1Bω（1）

式中：N=N×N;B是與F同等大小的二維矩陣，表示一系列具有不同特征的二維基函數(shù)，其構(gòu)建與F本身無關(guān);ω是與B對應(yīng)的權(quán)重系數(shù)，表示B對F貢獻(xiàn)的大小。式（1）中，ω大多數(shù)元素的值都接近于零，這表明可以使用少量CPT鉆孔數(shù)據(jù)（見圖2）確定ω中極少數(shù)重要元素，進(jìn)而合理估計出二維CPT數(shù)據(jù)F。令表示估算的ω。重建的二維CPT數(shù)據(jù)則表示為

=∑Nt=1B（2）

由式（2）可以看出，重建的關(guān)鍵是根據(jù)有限數(shù)量的CPT鉆孔數(shù)據(jù)估算。

1.2 二維CPT數(shù)據(jù)與的關(guān)系

令矩陣Y表示從場地中n（n>1）個CPT鉆孔中獲得的數(shù)據(jù)集，并用它來估計。Y與ω的關(guān)系可以通過Y和F之間的關(guān)系推導(dǎo)出來。如圖1、圖2所示，CPT數(shù)據(jù)Y是F的子集，為大小N×n的矩陣。根據(jù)式（1），F(xiàn)是一個與ω有關(guān)的函數(shù)，故而Y亦是ω的函數(shù)。數(shù)學(xué)上，Y表示為Y=FΨ。其中，Ψ反映了CPT沿水平方向（即圖2（a）中x軸）的勘測位置。上標(biāo)“T”表示轉(zhuǎn)置運(yùn)算。

此外，式（1）中F可以分解為矩陣形式F=BΩ（B），B與B分別代表為N×N和N×N的一維正交矩陣。F=BΩ（B）是式（1）的矩陣形式表達(dá)。將F=BΩ（B）代入Y=FΨ，可得

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

式中：A=ΨB。為運(yùn)算方便，通過依次疊加（Y）的列，將Y向量化為長度為M=（n× N） 的列向量y

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

式中：vec（·） 表示向量化;ω=vec（Ω），是式（1）中ω的向量化表示;“”表示克羅內(nèi)克積;A定義為

A=（BA）=bA…bA

bAbA（5）

有關(guān)克羅內(nèi)克積的更多細(xì)節(jié)問題可參考文獻(xiàn)[31]。根據(jù)式（4）可在貝葉斯框架下推導(dǎo)并估計。

1.3 的貝葉斯估計

因?yàn)樨惾~斯方法可以有效處理各種不確定性，如模型不確定性與空間變異性，因此采用該方法估計D。在貝葉斯框架下，的估計信息通過其后驗(yàn)概率密度函數(shù)（PDF）來反映，即P（|y）。

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

式中：P（y|）為似然函數(shù)，反映了給定條件下得到y(tǒng)的合理性;P（）是不考慮y時的先驗(yàn)PDF;P（y）是歸一化常數(shù)，保證P（|y） 的積分為常數(shù)1。

顯然，構(gòu)建式（6）中P（y|）和P（）是進(jìn)行貝葉斯估計的核心，分別構(gòu)建為高斯似然函數(shù)和多層先驗(yàn)分布。為了推導(dǎo)P（y|），以均值為零且方差未知的高斯隨機(jī)變量ε表示y與Aω之間的殘差。另外，為方便推導(dǎo)，將ε的方差的倒數(shù)表示為τ。假定殘差之間相互獨(dú)立，P（y|）可表示為

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

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

式中：τ未知，可將其作為隨機(jī)變量并與同步更新。因此，需構(gòu)建P（，τ）=P（）P（τ）。根據(jù)Zhao等、Zhao等的研究，的先驗(yàn)信息模型構(gòu)建如下：首先，令的每個元素均服從均值為零、方差為α的高斯分布，則P（|α）可表示為P（|α）=∏Nt=1α（2π）-1/2exp[-α（）/2]，其中，α是表示α的矢量形式;其次，假定α服從參數(shù)為1和γ/2 （γ > 0） 的逆伽瑪分布，即P（α |γ） = γ/2αexp[-γ/2·α] （t = 1， 2…N）且P（α|γ）=∏Nt=1P（α|γ）;最后，假定γ服從伽馬分布，即P（γ）=（b）aγexp（-bγ）/Γ（a），其中a和b為極小非負(fù)數(shù)（如a = b = 10），這對構(gòu)建^2D的無信息先驗(yàn)分布至關(guān)重要。此外，假定τ亦服從伽瑪分布，且P（τ）=dcτexp（-dτ）/Γ（c），c和d分別取值為1和10，可以使得P（τ）成為τ的無信息先驗(yàn)分布。

將P（|α）、P（α|γ）和P（τ）代入P（，α，γ，τ），可得到聯(lián)合先驗(yàn)PDF

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

由于將α和γ均視為隨機(jī)變量，因此在式（8）中使用了P（，α， γ， τ），而非P（， τ）?；谑剑?）中的貝葉斯定理、式（7）中的似然函數(shù)以及式（8）中的先驗(yàn)PDF，后驗(yàn)PDF可表示為P（，α，γ，τ|y）=P（y|，τ）×P（，α，γ，τ）/P（y）。由于P（y）無解析表達(dá)式，故P（，α，γ，τ|y）亦無解析解。然而，筆者發(fā)現(xiàn)，在給定其他3個參數(shù)的條件下，（，α， γ， τ） 中剩余參數(shù)的后驗(yàn)概率密度函數(shù)均有解析解。例如，以α、γ、τ為條件的的后驗(yàn)PDF服從多元高斯分布，均值和協(xié)方差矩陣為

μ=COVAyτ

COV=（AAτ+D）（9）

式中：D是N×N的對角矩陣，且對角線上元素為D（t，t）=α。其他幾個分布P（α|，τ，γ，y）、P（τ|，α，γ，y）和P（γ|，α，τ，y）分別服從廣義逆高斯分布和兩個伽瑪分布。更多細(xì)節(jié)可以參考文獻(xiàn)[32]，此處只將結(jié)果展示如下，其中P（α|，τ，γ，y）推導(dǎo)為

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

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

式中：a=;b=γ;p=-1/2;K（·）表示參數(shù)為p的第二類修正Bessel函數(shù)。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α。由于（， α， γ， τ）這些變量相互依存，難以得出的解析解。鑒于上述條件，概率密度函數(shù)的解析性即式（9）～式（11），筆者采用MCMC模擬中的吉布斯采樣算法最終估計獲得的統(tǒng)計特征。

2 的高效隨機(jī)模擬

2.1 吉布斯采樣（phthalic anhydride）

盡管P（|y）無解析解，但可采用吉布斯采樣方法生成大量后驗(yàn)樣本對其進(jìn)行表征。吉布斯采樣過程為：1）提供α、τ、和γ的初始值，將初始值代入式（9），計算的均值和協(xié)方差矩陣并由此生成一組樣本。2）將隨機(jī)產(chǎn)生的樣本與已知的τ、γ值同時代入式（10），產(chǎn)生一組α樣本。3）將隨機(jī)產(chǎn)生的、α樣本以及已知的γ，代入式（11a），τ為P（τ|，α，γ，y）中的唯一變量。4）在式（11b）中代入α、τ和，即γ為P（γ|，α，τ，y）中的唯一變量。5）收集α、τ和γ的樣本，并在步驟1）中替換α、τ和γ的值，然后重復(fù)該過程，直到獲得指定數(shù)量的MCMC樣本。值得注意的是，吉布斯采樣是馬爾科夫鏈蒙特卡洛（Markov Chain Monte Carlo， MCMC）的一種特殊形式，與傳統(tǒng)的MCMC方法不同，吉布斯采樣可以有效處理高維參數(shù)問題。這是由于吉布斯采樣需要已知待處理參數(shù)的條件概率分布，如式（9）～式（11）。相比之下，傳統(tǒng)MCMC方法中的算法，如Metroplis-Hasting更加泛化，在未知參數(shù)維度較高的時候，由于難以找到合適的高維提議分布（proposal probability density function），進(jìn)而導(dǎo)致生成樣本的接受率過低，計算效率因此受到極大影響。

由于服從多元高斯分布，可在步驟1）中用Cholesky分解方法隨機(jī)產(chǎn)生樣本，表示為

=μ+Lz（12）

式中：COV=（L）（L）L是下三角矩陣;z是長度為N=N×N的列向量，其中每個元素都是均值為零、方差為1的標(biāo)準(zhǔn)高斯隨機(jī)變量。式（12）表明，可通過獨(dú)立分布的z來隨機(jī)生成樣本。然而，由于的協(xié)方差矩陣COV過于龐大，基于式（12）的分解無法在普通個人電腦上執(zhí)行。此外，吉布斯采樣流程涉及COV的反復(fù)更新與分解，這就使得直接通過式（12）生成的樣本極其困難，進(jìn)而導(dǎo)致難以高效隨機(jī)生成非平穩(wěn)的二維CPT數(shù)據(jù)。

2.2 COV和μ的序列更新

由式（9）可知，COV是關(guān)于式（5）中已定義的A的函數(shù)。根據(jù)式（5），計算COV的因子AA可推導(dǎo)為

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

式中：B為單位正交矩陣。因此，I=BB為單位矩陣。將式（13）代入COV中，得到COV=[（I（τAA））+D]。根據(jù)克羅內(nèi)克積的定義，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個連續(xù)且互斥的對角線元素。由式（15）可知，COV為分塊對角矩陣，可按式（16）計算。

式中：COVi=[（τAA）+D]。式（16）將（N×N） × （N×N） 的巨型矩陣的逆矩陣運(yùn)算簡化為一系列尺寸為N×N小矩陣的逆矩陣的運(yùn)算。此外，式（16）中無B項(xiàng)，這進(jìn)一步提高了COV的分解效率。

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

μ=μμN(yùn)（17）

式中：μ=τ∑Nj=1bCOVAy，為μ的第i項(xiàng)第N個連續(xù)元素;y為y的第j項(xiàng)第n個連續(xù)元素;b為位于矩陣B第j行與第i列的元素。將式（16）、式（17）代入式（12），得

式中：=μ+Lz，表示的第N個連續(xù)元素;L是對應(yīng)于COV的下三角矩陣，且COV=（L）（L）;z是長度為N的列向量，表示N個獨(dú)立的、標(biāo)準(zhǔn)的高斯隨機(jī)變量。式（18）表明隨機(jī)樣本的產(chǎn)生轉(zhuǎn)變?yōu)镹組標(biāo)準(zhǔn)獨(dú)立高斯隨機(jī)變量z的隨機(jī)實(shí)現(xiàn)。相比式（12），式（18）中進(jìn)行的隨機(jī)生成時僅涉及一系列大小為（N×N）的矩陣運(yùn)算，避免了（N×N）×（N×N）巨型矩陣的求逆與分解。Xiao等、Ching等、Yang等也用到了類似方法。因此，該方法極大提高了樣本生成的計算效率。

將生成的隨機(jī)樣本代入式（2），生成大量非穩(wěn)態(tài)二維CPT數(shù)據(jù)。值得注意的是，該方法不需要對F的平穩(wěn)性進(jìn)行假設(shè)，因此，可直接用于估計具有空間變異性和多個未知邊界的土層的二維非穩(wěn)態(tài)CPT數(shù)據(jù)。雖然上述公式推導(dǎo)看似復(fù)雜，但通過商業(yè)軟件，如MATLAB 可以輕易地將數(shù)學(xué)表達(dá)式編譯為用戶函數(shù)。為了進(jìn)一步說明該方法的邏輯過程，這里給出了文中方法的流程圖，主要包含5個步驟。步驟1～3屬于方法的準(zhǔn)備階段，步驟4屬于方法的核心，步驟5屬于方法的結(jié)尾，詳見圖3。此外，為了便于讀者重復(fù)該工作，給出了步驟4，即方法核心的偽代碼，詳見圖4。

3 數(shù)值模擬案例

為了驗(yàn)證該方法，選取一組分布在3層土層內(nèi)的二維CPT錐尖阻力（Q）數(shù)據(jù)，如圖1所示。此例中，采用高斯平穩(wěn)隨機(jī)場生成器（如循環(huán)嵌入算法）在垂直方向和水平方向分別以η=0.02 m和η=0.5 m的分辨率生成每層內(nèi)的Q數(shù)據(jù)。隨機(jī)場模擬中涉及的參數(shù)包括Q的均值、標(biāo)準(zhǔn)差SD、垂直方向以及水平方向相關(guān)長度λ和λ。第1～3層的均值分別取為12、40和12;SD分別為2、5和2;相關(guān)長度分別為λ=2 m和λ=20 m，與文獻(xiàn)[40]一致。在生成每層Q數(shù)據(jù)時采用單指數(shù)自相關(guān)函數(shù)

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

式中：Δz=z-z和Δx=x-x分別表示位置（z，x）和（z，x）沿z和x方向的相對距離。利用上述隨機(jī)場參數(shù)和隨空間變化的巖土層邊界，可獲得3層內(nèi)的非穩(wěn)態(tài)二維Q數(shù)據(jù)，如圖1所示。盡管每層中的Q數(shù)據(jù)是穩(wěn)態(tài)的，但由于不同土層中Q的統(tǒng)計量不同，圖1所示的二維數(shù)據(jù)集是非穩(wěn)態(tài)的。

如圖2所示，當(dāng)n=5時的CPT鉆孔錐尖阻力數(shù)據(jù)可作為所提出方法的輸入Y （步驟1）。首先，令y=vec（Y） （步驟2）。之后，記錄Q沿x方向的位置，以此構(gòu)造矩陣Ψ;確定N和N，N代表Q數(shù)據(jù)點(diǎn)沿深度方向的數(shù)量，本例中N=500。N表示分辨率為η的插值后的CPT數(shù)據(jù)沿x方向的數(shù)量，計算公式為N=h/η+1，式中h代表場地的長度。例如，令h=127.5 m、η=0.5 m，則N=256。在MATLAB中輸入使用“dctmtx”以及N、N構(gòu)造一維正交矩陣B與B，且A=ΨB（步驟3）。

接下來按照圖4引入克羅內(nèi)克積的吉布斯采樣方法，生成大量樣本（步驟4）。此例中，首先隨機(jī)產(chǎn)生了2 100個樣本。然后每隔20步收集一次樣本，從而保證在丟棄前100個不可信的樣本后，樣本之間的統(tǒng)計相關(guān)性較弱。如此，共生成n = 100個相互獨(dú)立的隨機(jī)樣本。將每次產(chǎn)生的樣本代入式（2），產(chǎn)生一個二維的Q樣本。100組樣本產(chǎn)生了n = 100組二維且獨(dú)立的Q樣本（步驟5）。圖5給出了上述兩個CPT樣本。值得注意的是，雖然這里只生成了100個相互獨(dú)立的Q樣本，但其統(tǒng)計特征基本收斂，如每一點(diǎn)的均值及標(biāo)準(zhǔn)差并不隨著MCMC樣本的增多而發(fā)生變化。

為了說明引入克羅內(nèi)克積對計算效率的提升，記錄生成上述100個獨(dú)立CPT樣本的時間，見表1。由表1可見，使用64位Windows 10操作系統(tǒng)的IntelCorei7-9700 3.0 GHz CPU和16.0 GB RAM的計算機(jī)，生成100組獨(dú)立二維CPT樣本大約需40.6 s。然而，若不采用序列更新方法，即直接采用式（12）生成樣本時，同一臺計算機(jī)上大約需1 421.2 s。相比于原始方法，該算法在計算效率上提高了約35倍。隨著CPT勘測鉆孔數(shù)量n的增加，使用序列更新技術(shù)的計算效率提高會更為明顯。

利用上述100個二維Q樣本可得其均值，如圖6所示。結(jié)果表明，即使在未知地下邊界處，Q均值的分布也與圖1所示較為一致。表明利用該方法生成的二維非穩(wěn)態(tài)Q樣本在統(tǒng)計上保留了圖1中Q原始數(shù)據(jù)的非穩(wěn)態(tài)模式。此外，根據(jù)生成的二維Q樣本還可計算每一點(diǎn)的標(biāo)準(zhǔn)差，見圖7。結(jié)果表明，CPT鉆孔位置處的標(biāo)準(zhǔn)差非常小，說明在這些位置估計的Q可靠性和可信度較高。因?yàn)檫@些位置的Q數(shù)據(jù)已知，并將其作為建議方法的輸入。相反，因?yàn)檫h(yuǎn)離測量點(diǎn)的位置有效信息極少，導(dǎo)致遠(yuǎn)離CPT鉆孔位置的樣本標(biāo)準(zhǔn)差相對較大。

為了進(jìn)一步驗(yàn)證該方法插值的合理性，比較圖2（a）中4個未測量鉆孔（即BH、BH、BH、BH）的Q數(shù)據(jù)。圖8（a）～（c）分別以虛線繪制了該方法在這4個鉆孔位置插值的Q數(shù)據(jù)。為進(jìn)行比較，圖8以粗線給出了BH至BH處的原始Q變化曲線。圖8顯示，每個子圖中虛線的變化趨勢與實(shí)線一致，表明利用提出的方法生成的Q樣本是合理、可靠的，并較為合理地刻畫了圖1中CPT數(shù)據(jù)的非平穩(wěn)特點(diǎn)。

圖8 位置BH到BH處的Q均值與原始數(shù)據(jù)的對比

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

值得注意的是，圖8中虛線、灰線和粗實(shí)線之間存在一些差異。這是因?yàn)槭褂锰岢龅姆椒〞r，輸入的CPT鉆孔數(shù)量較少。當(dāng)n增加時，二維Q樣本會更加接近CPT真實(shí)數(shù)據(jù)。

4 CPT測深數(shù)量的影響

采用n=5、15、25和50的案例探討該方法的性能。對于每個n案例，前述方法隨機(jī)生成100個獨(dú)立的二維Q樣本，然后利用這100個Q樣本計算每個位置的Q均值及標(biāo)準(zhǔn)差，詳見圖9、圖10。由圖9可以看出，隨著n的增加，每個位置的Q樣本均值越來越接近于圖1。當(dāng)n=25時，Q樣本均值與實(shí)測Q幾乎相同。此外，圖10顯示，隨著n的增加，每個位置估算的Q標(biāo)準(zhǔn)差不斷減小，表明每個位置Q估算值的可靠性和置信度都有所提高。該計算結(jié)果反映了本文所提方法的數(shù)據(jù)驅(qū)動本質(zhì)。

另外，隨著n的增加，采用序列更新技術(shù)的計算時間會略有延長，如表1所示。但與不采用序列更新技術(shù)的方法相比，該方法所增加的計算時間幾乎可以忽略。當(dāng)n=15時，該方法僅需約79.6 s，而未采用序列更新技術(shù)的方法需9 h以上（見表1）。兩種方法計算時間的差距表明，所提出的方法可顯著提高計算效率。

值得強(qiáng)調(diào)的是，該方法可以較為合理地考慮CPT鉆孔數(shù)量導(dǎo)致的不確定性，為了定量說明不確定性的大小，根據(jù)圖9、圖10分別計算出了n=5、15、25、50不同情況下變異系數(shù) （δ=σ/μ×100%）的中位數(shù)，分別是11.80%、10.86%、9.36%與8.10%。此外，根據(jù)計算經(jīng)驗(yàn)，該方法能適用的最少CPT鉆孔數(shù)在4～5個左右，如果CPT鉆孔更少，應(yīng)謹(jǐn)慎使用該方法。

5 工程應(yīng)用

為進(jìn)一步驗(yàn)證方法的有效性，選取位于日本岡山河堤的某工程場地進(jìn)行驗(yàn)證研究。該場地自上而下主要由砂土、黏土或粉土、砂土組成。其中上層砂土厚約2 m，黏土最厚約3 m，粉土厚約1～2 m，粉土下面仍主要以砂土為主。為了探明該河堤的軟土空間分布，日本工程師在此進(jìn)行了一系列的CPT。選取其中的CPT-201～CPT-207來驗(yàn)證該方法。圖11給出了上述7個CPT鉆孔的位置分布圖，其中，CPT-204用來驗(yàn)證，其余CPT標(biāo)準(zhǔn)化Q數(shù)據(jù)當(dāng)作輸入。所有的數(shù)據(jù)均來自TC304免費(fèi)數(shù)據(jù)庫（http：//140.112.12.21/issmge/tc304.htm）。按照該方法，可以得到Q的二維剖面分布，其均值與標(biāo)準(zhǔn)差見圖12。由圖12可知，估計的Q數(shù)據(jù)離CPT鉆孔位置愈接近，對應(yīng)的標(biāo)準(zhǔn)差愈小，反之，則愈大。由于這是真實(shí)原位試驗(yàn)數(shù)據(jù)，無法得知CPT-204以外的真實(shí)數(shù)據(jù)。這里以CPT-204對應(yīng)的Q數(shù)據(jù)與此位置的估計數(shù)據(jù)進(jìn)行對比，見圖13。圖13顯示，雖然估計值與試驗(yàn)值之間存在一定差距，但兩者趨勢基本一致;此外，CPT-204的大部分實(shí)測試驗(yàn)值落在估計值的95%置信區(qū)間里，間接說明該方法的準(zhǔn)確性以及魯棒性。

6 結(jié)論

基于吉布斯采樣與貝葉斯壓縮感知方法，提出了一種快速估計未采樣位置的非平穩(wěn)靜力觸探測試數(shù)據(jù)（CPT）的方法。該方法屬于無參估計范疇，能夠自動考慮靜力觸探數(shù)據(jù)沿水平方向的相關(guān)性，同時回避了水平、深度方向自相關(guān)函數(shù)的采用以及數(shù)據(jù)平穩(wěn)性等假設(shè)。生成的CPT數(shù)據(jù)可用于解決各種巖土工程和地質(zhì)工程問題，如土壤分層和分區(qū)、土壤液化潛力評估及其空間變異性表征、巖土設(shè)計參數(shù)的間接估計等。該方法為解決實(shí)際工程問題提供了一種概率工具，尤其是針對在實(shí)踐中經(jīng)常遇到的沿水平方向的CPT勘測鉆孔數(shù)量較少的情況。最后，通過數(shù)值與實(shí)際工程案例對提出的方法進(jìn)行驗(yàn)證，結(jié)果表明，該方法較為準(zhǔn)確，并且采用序列更新技術(shù)后可顯著提高計算效率，具有較強(qiáng)的魯棒性。

參考文獻(xiàn)：

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

[2] 張繼周， 繆林昌， 王華敬. 土性參數(shù)不確定性描述方法的探討[J]. 巖土工程學(xué)報， 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] 劉松玉， 蔡正銀. 土工測試技術(shù)發(fā)展綜述[J]. 土木工程學(xué)報， 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] 劉松玉， 吳燕開. 論我國靜力觸探技術(shù)（CPT）現(xiàn)狀與發(fā)展[J]. 巖土工程學(xué)報， 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] 沈小克， 蔡正銀， 蔡國軍. 原位測試技術(shù)與工程勘察應(yīng)用[J]. 土木工程學(xué)報， 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] 孟高頭， 張德波， 劉事蓮， 等. 推廣孔壓靜力觸探技術(shù)的意義[J]. 巖土工程學(xué)報， 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]. 巖土工程學(xué)報， 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的土分類方法在港珠澳大橋中的應(yīng)用[J]. 巖土工程學(xué)報， 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] 林軍， 蔡國軍， 劉松玉， 等. 基于孔壓靜力觸探力學(xué)分層的土體邊界識別方法研究[J]. 巖土力學(xué)， 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]. 巖土工程學(xué)報， 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勘察信息融合的二維土性參數(shù)剖面貝葉斯表征方法[J]. 應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報， 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] 王長虹， 朱合華， 錢七虎. 克里金算法與多重分形理論在巖土參數(shù)隨機(jī)場分析中的應(yīng)用[J]. 巖土力學(xué)， 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] 劉志平， 何秀鳳， 張淑輝. 多測度加權(quán)克里金法在高邊坡變形穩(wěn)定性分析中的應(yīng)用[J]. 水利學(xué)報， 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] 趙騰遠(yuǎn)， ALADEJARE Adeyemi Emman， 王宇. 基于貝葉斯方法的模型選擇以及巖石性質(zhì)概率表征[J]. 武漢大學(xué)學(xué)報（工學(xué)版）， 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] 曹子君， 趙騰遠(yuǎn)， 王宇， 等. 基于貝葉斯等效樣本的土體楊氏模量的統(tǒng)計特征確定方法[J]. 防災(zāi)減災(zāi)工程學(xué)報， 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.

（編輯 黃廷）