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

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

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

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

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

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

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

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

(1)

(2)

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

圖2 CPT測(cè)深位置及孔內(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ī)場(chǎng)生成器(如循環(huán)嵌入算法)在垂直方向和水平方向分別以ηz=0.02 m和ηx=0.5 m的分辨率生成每層內(nèi)的Qt數(shù)據(jù)[39]。隨機(jī)場(chǎng)模擬中涉及的參數(shù)包括Qt的均值、標(biāo)準(zhǔn)差SD、垂直方向以及水平方向相關(guān)長(zhǎng)度λz和λx。第1～3層的均值分別取為12、40和12；SD分別為2、5和2；相關(guān)長(zhǎng)度分別為λz=2 m和λx=20 m，與文獻(xiàn)[40]一致。在生成每層Qt數(shù)據(jù)時(shí)采用單指數(shù)自相關(guān)函數(shù)

(19)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 工程應(yīng)用

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

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

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

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

6 結(jié)論

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