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

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

(1. 西安科技大學 地質與環(huán)境學院，西安 710054；2.西安交通大學 人居環(huán)境與建筑工程學院，西安 710049；3. 香港城市大學 建筑學及土木工程學系，香港 999077)

巖土工程場地勘察的目標是通過室內(nèi)或現(xiàn)場原位試驗合理地描述地下土層界面并表征土層內(nèi)土體物理、力學參數(shù)的空間變異性[1-2]。與室內(nèi)試驗[3]相比，原位測試方法，如靜力觸探試驗[4](Cone Penetration Test，CPT)成本更加低廉、測試更加快捷[5-6]。此外，靜力觸探試驗(CPT)能在深度方向獲得近乎連續(xù)性的錐尖阻力(qc)和側摩阻力(fs)數(shù)據(jù)，并以此作為土體的力學響應[7]。因此，CPT廣泛應用于表征巖土場地參數(shù)、地下土體分層[8-11]、砂土液化評估[12-13]、確定隨機場的相關函數(shù)和參數(shù)等[14-15]。

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

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

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

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

(1)

(2)

圖1 有兩個土層邊界的二維非平穩(wěn)Qt數(shù)據(jù)Fig.1 A set of 2D non-stationary Qt data with two

圖2 CPT測深位置及孔內(nèi)Qt隨深度的變化曲線Fig.2 CPT sounding locations and variation curves

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

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11a)

(11b)

2.1 吉布斯采樣(phthalic anhydride)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

圖3 生成二維CPT數(shù)據(jù)的流程圖Fig.3 A flow chart for simulating 2D CPT

圖4 引入克羅內(nèi)克積吉布斯采樣的偽代碼Fig.4 Pseudo code for fast simulation of using the Gibbs sampling method with Kronecker

3 數(shù)值模擬案例

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

(19)

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

圖5 二維CPT數(shù)據(jù)模擬結果示例Fig.5 Two examples of simulated 2D CPT

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

表1 不同CPT鉆孔數(shù)目下，使用序列更新技術(A)和不使用序列更新技術(B)生成100組二維Qt樣本的計算時間

*注：因運算所需內(nèi)存過大，方法B在筆者電腦中無法運行。

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

圖6 MCMC模擬Qt數(shù)據(jù)的均值Fig.6 Averaged Qt from the simulated MCMC

圖7 MCMC模擬Qt數(shù)據(jù)的標準差Fig.7 SD of Qt from the simulated MCMC

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

圖8 位置BH1到BH4處的Qt均值與原始數(shù)據(jù)的對比Fig.8 Comparison between averaged Qt profile and the original one at locations BH1 to

值得注意的是，圖8中虛線、灰線和粗實線之間存在一些差異。這是因為使用提出的方法時，輸入的CPT鉆孔數(shù)量較少。當nb增加時，二維Qt樣本會更加接近CPT真實數(shù)據(jù)。

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

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

圖9 nb對二維Qt樣本均值的影響Fig.9 Effect of nb on 2D Qt averaged

圖10 nb對MCMC樣本二維Qt數(shù)據(jù)標準差的影響Fig.10 Effect of nb on SD of 2D Qt estimated from

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

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

5 工程應用

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

圖11 CPT鉆孔201到207分布圖Fig.11 Spatial distribution of CPT soundings ranging

圖12 日本岡山河堤Qt剖面估計結果Fig.12 Estimated Qt profile of the Okayama Riverbank in

圖13 CPT-204估計值與實測值的對比Fig.13 Comparison between predicted and measured

6 結論

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