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

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

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

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

(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ù)與的關(guān)系

(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ù)值模擬案例

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

(19)

式中：Δz=zi-zm和Δx=xj-xn分別表示位置(zi,xj)和(zm,xn)沿z和x方向的相對距離。利用上述隨機(jī)場參數(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ù)模擬結(jié)果示例Fig.5 Two examples of simulated 2D CPT

為了說明引入克羅內(nèi)克積對計算效率的提升，記錄生成上述100個獨(dú)立CPT樣本的時間，見表1。由表1可見，使用64位Windows 10操作系統(tǒng)的Intel?Core?i7-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ù)量nb的增加，使用序列更新技術(shù)的計算效率提高會更為明顯。

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

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

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

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

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

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

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

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

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

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

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

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

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

值得強(qiáng)調(diào)的是，該方法可以較為合理地考慮CPT鉆孔數(shù)量導(dǎo)致的不確定性，為了定量說明不確定性的大小，根據(jù)圖9、圖10分別計算出了nb=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)化Qt數(shù)據(jù)當(dāng)作輸入。所有的數(shù)據(jù)均來自TC304免費(fèi)數(shù)據(jù)庫(http://140.112.12.21/issmge/tc304.htm)。按照該方法，可以得到Qt的二維剖面分布，其均值與標(biāo)準(zhǔn)差見圖12。由圖12可知，估計的Qt數(shù)據(jù)離CPT鉆孔位置愈接近，對應(yīng)的標(biāo)準(zhǔn)差愈小，反之，則愈大。由于這是真實(shí)原位試驗(yàn)數(shù)據(jù)，無法得知CPT-204以外的真實(shí)數(shù)據(jù)。這里以CPT-204對應(yīng)的Qt數(shù)據(jù)與此位置的估計數(shù)據(jù)進(jìn)行對比，見圖13。圖13顯示，雖然估計值與試驗(yàn)值之間存在一定差距，但兩者趨勢基本一致；此外，CPT-204的大部分實(shí)測試驗(yàn)值落在估計值的95%置信區(qū)間里，間接說明該方法的準(zhǔn)確性以及魯棒性。

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

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

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

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)的魯棒性。