樁底溶洞探測中瞬時相位差強度計算方法的優(yōu)化及應用*

石振明 沈 健 劉 鎏 彭 銘 王登一 楊沛權 何 梟

(①同濟大學地下建筑與工程系， 上海 200092， 中國)

(②同濟大學巖土及地下工程教育部重點實驗室， 上海 200092， 中國)

(③中國科學院武漢巖土力學研究所， 巖土力學與工程國家重點實驗室， 武漢 430071， 中國)

(④廣東建科交通工程質(zhì)量檢測中心有限公司， 廣東 650506， 中國)

0 引 言

巖溶地貌在全球分布廣泛，而我國是世界上巖溶地質(zhì)發(fā)育最廣泛的國家之一，總面積達344.3萬平方千米(劉鎏， 2021)，約占國土面積的1/3(陳學軍等， 2019)。在巖溶地貌區(qū)開展地下工程建設通常充滿風險。因為時常存在無規(guī)則巖溶分布與填充情況(姚邦杰等， 2019)，不充分的樁底溶洞勘探可能導致護壁泥漿泄漏、鉆孔坍塌、周圍建筑開裂甚至周圍建筑結(jié)構(gòu)倒塌等嚴重工程事故(石振明等， 2016)。

傳統(tǒng)的探測方法提供了許多樁底溶洞勘探方法，但這些方法受灌注樁樁孔內(nèi)高壓泥漿環(huán)境影響，無法被應用于樁底溶洞探測。而彈性波CT法(朱文仲等， 2008)和管波測試法(李學文等， 2005)雖能被用于灌注樁樁底溶洞探測，但受限于其探測精度不足以及鉆孔布置繁瑣，并不常被用于實際工程的樁底溶洞探測?；谶@一現(xiàn)狀，相關研究團隊提出了一種鉆孔灌注樁樁底溶洞聲吶探測方法(劉鎏， 2021)，將樁孔中的泥漿作為聲波傳播和耦合的介質(zhì)，利用廣義S變換進行了聲波信號分析，揭示了樁底聲吶信號物理意義，實現(xiàn)了樁底溶洞的快速無損探測。但該分析方法是基于傳統(tǒng)Fourier變換提出的，在面對樁孔內(nèi)不規(guī)則多次反射面波的干擾時分析能力不足，對溶洞位置準確揭示能力受到影響。

相位是信號本質(zhì)特征之一，能夠表征信號任意時刻的狀態(tài)，是信號分析過程中重要的考察對象。傳統(tǒng)的信號分析方法包括相位反正切、DACM算法以及線性解調(diào)算法等。但這些方法只適用于單一分量信號，受噪聲影響較大，均無法應用于樁底聲吶探測的復雜信號分析。而希爾伯特-黃變換(Hilbert-Huang Transform， HHT)能夠?qū)碗s信號中的噪聲濾去，從而獲得信號準確的相位信息。如圖1a展示了經(jīng)HHT分析后得到兩個不同偏移距的樁底溶洞探測信號瞬時相位信息，從圖中不難發(fā)現(xiàn)信號相位均在-π到π間變化。圖1b則展示了兩信號間的瞬時相位差，可以看到瞬時相位差在2 ms及4 ms處接近于0，但這一特點仍無法直接清晰、定量地識別出來。

圖1 瞬時相位及瞬時相位差(Liu et al.， 2021)

基于瞬時相位信息特征，學者提出了一種雙通道樁底溶洞聲吶探測方法(Liu et al.， 2021)。該探測方法通過在樁底布置“一發(fā)雙收”(雙通道即指雙接收器)的聲吶接發(fā)換能器，獲取兩道聲吶應力波信號記錄，運用HHT分析方法從中獲得相應瞬時相位信息，并以定量指標-瞬時相位差強度(Instantaneous Phase Difference Intensity， IPDI)來衡量兩道信號間的相位差異，將溶洞頂?shù)装宓挠行w波反射信號從樁孔壁以及樁孔底多次反射面波信號中提取出來，從而實現(xiàn)溶洞頂?shù)装宓淖R別。

但在計算信號IPDI值時，HHT方法中各個關鍵參數(shù)因影響規(guī)律不清而隨意選取以及IPDI圖像兩端可能出現(xiàn)的端點效應將會導致溶洞頂?shù)装逦恢米R別不夠精準，亟待優(yōu)化處理。

本文根據(jù)數(shù)值仿真模擬信號測試結(jié)果，探究了雙通道樁底溶洞探測方法中IPDI分析過程的關鍵參數(shù)影響規(guī)律和端點效應抑制方法，提出優(yōu)化建議，并結(jié)合工程現(xiàn)場實測信號驗證了優(yōu)化方法的有效性，提升了溶洞頂?shù)装逦恢米R別的精準度。

1 溶洞聲吶探測方法及數(shù)值仿真

1.1 雙通道樁底溶洞聲吶探測系統(tǒng)

雙通道樁底溶洞聲吶探測系統(tǒng)示意圖如圖2所示。圖中O點位于樁孔底中心處，代表了聲吶應力波發(fā)射換能器。而A、B兩點分別代表了兩個相隔一定距離的接收換能器，接收沿樁孔壁及樁孔底傳播的多次反射面波和溶洞頂?shù)装宓奈⑷躞w波反射信號。C、D點分別代表溶洞頂板和底板。圖中OA、AB間距離(即發(fā)射器至接收器距離以及兩接收器間距離)都是人為可調(diào)的。當根據(jù)面波波速以及接收信號主頻設置AB間的距離為半面波波長時，信號記錄中面波的相位差近似為±π。同時當溶洞離開樁孔底距離較遠時，兩接收換能器中P波反射信號的相位差就應當接近于0。反射體波相位差與多次反射面波相位差的明顯不同為從信號記錄中提取出所需的微弱反射體波信號并由此計算溶洞頂?shù)装逦恢锰峁┝丝赡堋?/p>

圖2 雙通道樁底溶洞聲吶探測系統(tǒng)示意圖

1.2 數(shù)值仿真

為進行數(shù)值模擬測試，首先建立了二維樁底溶洞聲吶探測的數(shù)值模型，利用二維的黏彈性波方程(式(1)～式(3)，式中變量均以二階張量形式表示)和交叉網(wǎng)格時域有限差分(Ke et al.，2016)獲取數(shù)值模擬結(jié)果。

(1)

(2)

(3)

式中：vi為速度分量；τij為應力分量；ρ為介質(zhì)密度；λ和G是拉梅常量；εL和εT分別為縱波阻尼系數(shù)和橫波阻尼系數(shù)。

分析探測波場和波剖面包括了樁孔、基巖和樁底溶洞(圖3)。模擬計算區(qū)域的大小為8m×8m，樁孔模型位于計算區(qū)域的上部中心，其直徑Dp為1.6m，深度為2m。溶洞頂板距離樁孔底距離為Hc，其為直徑Dc。且通過設置不同波速、密度以及阻尼系數(shù)可以模擬不同溶洞填充物情況。豎直方向的入射波將在樁底孔的中心向下發(fā)射，該點即為發(fā)射換能器，其坐標(0， 0)。同時在計算邊界外加入了完美匹配邊界(Perfectly matched layer， PML)(Berenger， 1994)來吸收并壓制在計算邊界上產(chǎn)生的反射以模擬波場實際傳播過程。圖3中PML邊界厚度為2m，分布于模型四周邊界。其中紅色區(qū)域(見電子版文章)表示邊界重合處。

圖3 樁底溶洞聲吶探測數(shù)值模擬模型

數(shù)值模擬中不同介質(zhì)的材料模擬參數(shù)設定見表1，其具體過程如下：

表1 數(shù)值模擬的材料參數(shù)

(1)根據(jù)樁孔和樁底溶洞的地質(zhì)結(jié)構(gòu)建立交叉網(wǎng)格時域有限差分模型，在考慮計算精度和計算成本后，選取合適空間和時間精度。此處選取空間10階，時間2階的差分形式。

(2)設置所有參數(shù)，包括不同材料介質(zhì)物性參數(shù)以及發(fā)射、接收裝置位置坐標。這里設置兩接收器間距離為半面波波長。

(3)設置PML邊界條件，輸入PML邊界參數(shù)。

(4)進行時域有限差分，求解黏彈性波動方程。

(5)記錄模擬信號并使用IPDI方法對其進行分析。

1.3 IPDI計算方法

瞬時相位差強度(IPDI)其計算流程如圖4所示，可以簡單概括為：(1)布置探測系統(tǒng)，并使兩個接收器之間的距離為半面波波長；(2)運用EEMD將原始信號分解為一系列固有經(jīng)驗模態(tài)函數(shù)(IMF)；(3)選擇合適的IMF組合成為重構(gòu)信號；(4)對重構(gòu)信號作用希爾伯特變換得到瞬時相位，并計算瞬時相位差；(5)計算瞬時相位差強度(IPDI)，其具體定義如式(4)：

圖4 瞬時相位差強度計算流程圖

(4)

式中：t為任意時刻；T為信號周期，數(shù)值上等于信號主頻的倒數(shù)。當t為反射波到達時刻后一個周期的中點時刻t0，I(t)達到最大值，在IPDI圖像中表現(xiàn)為峰值。因此，反射波到達時刻ta可以定義為：

(5)

此過程中的總體平均經(jīng)驗模態(tài)分解(Ensemble Empirical Mode Decomposition， EEMD)過程意義在于將接受到的信號中多余噪聲濾去而盡可能保留真實相位信息。其具體計算步驟如下：

(1)給待處理的信號x(t)加上一組高斯白噪聲w(t)，得到新的總體信號X(t)

X1(t)=x(t)+w1(t)

(6)

(2)對X1(t)進行EMD處理，得到其分解結(jié)果

(7)

式中：c1k為分解所得的各IMF分量；r1n為分解所得殘差。

(3)再給原信號x(t)添加N-1組不同的白噪聲，重復EMD過程

Xi(t)=x(t)+wi(t)，i=2，3，…，N

(8)

(9)

(4)取相應的IMF均值為最終的IMF

(10)

從以上步驟可以看到，EEMD關鍵在于加入了N組一定幅值的白噪聲來消除可能出現(xiàn)的模態(tài)混疊現(xiàn)象，而所加入白噪聲的具體組數(shù)以及其相應幅值是人為決定的。這將可能影響提取出來相位信息的真實程度。關于這一方面的討論將會在下一節(jié)詳細展開。

2 溶洞探測精度影響因素分析及優(yōu)化

這一部分將結(jié)合樁底溶洞聲吶探測仿真信號測試采樣率、EEMD中總體平均次數(shù)以及標準差對IPDI分析結(jié)果的影響規(guī)律，并嘗試就這些參數(shù)取值以及端點效應的抑制給出優(yōu)化意見。

2.1 采樣率分析及參數(shù)優(yōu)化

現(xiàn)實生活中我們很難記錄下連續(xù)信號，因此只能通過一定的采樣率將連續(xù)信號時間序列離散化(Flandrin et al.，2004)。采樣率越高，越能使記錄下的離散信號逼近原始的連續(xù)信號。而EEMD分解通過對信號的極值點進行插值計算來形成上下包絡線，并以上下包絡線的均值作為各階IMF，因此過小的采樣率可能導致信號極值點的位置和取值發(fā)生誤差，最終使獲得的IMF不夠準確，降低了IPDI值對信號異常相位信息的提取能力(Liu et al.，2021)。

圖5a展示了在采樣頻率分別為0.5MHz、1MHz和2MHz時計算所得的IPDI值。該算例模擬了溶洞直徑為2m、溶洞頂板距樁孔底距離3m、填充物為黏土的情況。圖中2 ms處和4 ms處附近的峰值分別代表了溶洞頂板以及底板的反射信號。而3 ms處的峰值代表了溶洞頂板反射PS波的垂直分量。從圖中可以看到，采樣率為0.5MHz和1MHz時，IPDI峰值重合，其所對應的時刻分別為：C-1.904 ms， D-3.924 ms。結(jié)合基巖波速、黏土填充波速以及式(2)可以算得溶洞頂板距樁孔距離2.928m，溶洞直徑2.02m，誤差均小于5%。但當采樣率為2MHz時，其計算所得溶洞頂?shù)装逦恢镁葹槲吹玫矫黠@提升但運算所需時間大幅增長。

圖5 信號采樣率對IPDI計算結(jié)果影響

圖5b展示了在采樣率分別為0.25MHz和0.5MHz時計算所得的IPDI值。從圖中可以發(fā)現(xiàn)，當采樣率為0.25MHz時， 2 ms處IPDI峰值發(fā)生了偏離，其對應的時刻為1.896 ms，由此計算所得的距離偏小。并且4 ms處的峰值所對應的具體時刻辨認困難，為溶洞直徑的計算制造了困難。

綜合以上兩例，當進行IPDI計算時，為保證計算結(jié)果的精確，采樣率不宜過小，數(shù)值模擬結(jié)果表明采樣率通常取0.5～1MHz能夠使誤差足夠小。

2.2 總體平均次數(shù)以及標準差參數(shù)分析及優(yōu)化

雖然EEMD是一種自適應的信號時頻分析方法，但在其分解過程中仍需人為選定兩個參數(shù)：總體平均次數(shù)N以及加入白噪聲的標準差ε。其中：總體平均次數(shù)決定了加入白噪聲的組數(shù)，而標準差則決定了白噪聲幅值與原信號幅值之比。通?？傮w平均參數(shù)取100～200，而標準差大致取0.2能夠取得比較理想的分解結(jié)果(Wu et al.，2009)?，F(xiàn)就采樣率為0.25MHz時探索這兩個參數(shù)的何種取值方式能夠使IPDI計算結(jié)果得到提升。

圖6a展示了當標準差ε分別選取0.05、0.1、0.4時計算所得IPDI曲線。標準差為0.4及0.1時，溶洞頂板反射信號在IPDI曲線中對應的時刻均為1.896 ms； 標準差為0.05時，溶洞頂板反射信號在IPDI曲線中對應的時刻為1.88 ms，其計算所得溶洞頂板距樁孔底距離(2.88 ms)誤差最大。標準差為0.05及0.1時，溶洞頂板反射信號在IPDI曲線中對應的時刻均為3.896 ms； 標準差為0.4時，溶洞頂板反射信號在IPDI曲線中對應的時刻為3.948ms。綜合來看，當標準差取0.1時，其計算結(jié)果最準確(溶洞直徑2m，溶洞頂板距樁孔底距離2.912m)。

圖6 EEMD參數(shù)取值對IPDI計算結(jié)果影響

圖6b展示了當總體平均次數(shù)N分別選取100、300和600時計算所得IPDI曲線。溶洞頂板反射信號在IPDI曲線中對于3個總體平均次數(shù)取值均為1.896 ms； 而當總體平均次數(shù)取300和600時， 4 ms處溶洞底板反射信號峰值更清晰，更容易判斷其對應時刻。

綜合以上分析，在進行EEMD過程中，標準差取0.1且總體平均次數(shù)取較大值如300～600能夠提升最終溶洞識別準確度。這一建議與式(11)(Wu et al.，2009)相吻合。

(11)

式中：ε代表白噪聲標準差；N代表總體平均次數(shù)；εn代表最終誤差的標準差，即輸入信號與相應IMF間的誤差。

2.3 基于灰色預測模型的端點效應抑制

端點效應指溶洞聲吶探測方法所計算得到IPDI圖像中，其兩端出現(xiàn)異常“飛逸”現(xiàn)象，如圖7左右兩端。端點效應作為HHT公開的七大理論問題(Rilling et al.，2009)，已經(jīng)引起了學者們的廣泛關注。相關理論研究表明，端點效應可能同時存在于EMD以及HT兩個過程中。前者主要由左右端點插值精度差導致，后者發(fā)生的原因由傅里葉變換中“Gibbs現(xiàn)象”引起的頻率泄露導致。

圖7 端點效應抑制

為解決這一問題，本節(jié)選擇了灰色模型預測模型(賀智， 2011)來進行端點效應的處理?；疑A測模型是基于灰色系統(tǒng)理論發(fā)展的(Kayacan et al.，2010； Li et al.，2010)，其本質(zhì)在于通過運用信號端點少量數(shù)據(jù)來對信號進行延拓，包括EMD過程中根據(jù)端點極值的延拓以及HT變換前對IMF的延拓，其具體原理如下：

本節(jié)嘗試在構(gòu)造完成的重構(gòu)信號基礎上，在其數(shù)據(jù)兩端分別拾取3m、3n個數(shù)據(jù)點，建立灰色預測模型，再各自向兩端外預測m以及n個數(shù)據(jù)點，得到新的重構(gòu)信號，經(jīng)HT變換后再選擇與原始信號相同時間段的重構(gòu)信號數(shù)據(jù)點，最終計算得到IPDI值。如圖7，可以看到經(jīng)灰色預測模型延拓原信號數(shù)據(jù)點后，其IPDI曲線末端的“飛逸”得到明顯改善，證明了此方法對改善端點效應的有效性。但同時可以發(fā)現(xiàn)經(jīng)端點效應優(yōu)化處理后其峰值位置出現(xiàn)了一定偏移，可能影響最終溶洞頂?shù)装逦恢玫呐卸ā?/p>

3 IPDI優(yōu)化方法應用

本節(jié)嘗試就1個數(shù)值模擬工況以及1個現(xiàn)場實際工程所得信號以測試上文所述的IPDI值優(yōu)化方法，并證明其有效性。下面介紹這些工況和其相應優(yōu)化結(jié)果。

3.1 數(shù)值模擬工況

該工況下溶洞直徑Dc設置為2m，溶洞頂板距樁孔底距離Hc設置為3m均未改變，只通過設置不同的波速、密度以及阻尼系數(shù)(表1)來模擬樁孔下部溶洞被3種不同材質(zhì)的填充材料填充——水、空氣以及黏土(圖8)。圖9a為未作優(yōu)化處理的該工況下IPDI值圖像， 2 ms左右處的峰值反映了溶洞頂板的反射波信號， 4 ms左右處的峰值則反映了溶洞底板的反射波信號。可以看到，只有黏土填充和水填充的IPDI圖像上在4 ms處能看到明顯峰值，而由于波在空氣中傳播速度慢導致其反射信號未在采樣時窗內(nèi)被捕獲。為優(yōu)化該工況下的IPDI值圖像，采用了提高采樣率至1MHz以及端點效應優(yōu)化的方法(圖9b)。從優(yōu)化結(jié)果中可以看到，黏土填充溶洞底板反射信號IPDI峰值更明顯，其對應時刻更容易判斷，為3.954 ms(對應可算得溶洞直徑1.938m)。而且，水填充溶洞頂板反射信號IPDI峰值對應時刻由1.898 ms變?yōu)?.961 ms，計算所得溶洞頂板距樁孔底距離的相對誤差由2.8%改進為1.4%。此外，本文所采取的端點效應優(yōu)化方法一定程度上壓制了信號兩端的“飛逸”現(xiàn)象。

圖8 數(shù)值模擬工況示意圖

圖9 數(shù)值模擬工況信號IPDI值優(yōu)化結(jié)果

3.2 現(xiàn)場實際信號

該信號取自廣清高速改擴建工程的一樁基建設場地。該工程沿線有29km的巖溶區(qū)，主要為下石炭統(tǒng)的淺海相夾含煤碎屑巖，由深灰、灰黑色灰?guī)r、白云質(zhì)灰?guī)r、泥灰?guī)r和硅質(zhì)灰?guī)r組成。舊路勘察資料、施工記錄及新路初勘資料表明，巖溶區(qū)巖溶分布密度大，見洞率高，巖溶發(fā)育深度正處在路基、橋基的影響深度。為保證工程安全進行，有必要在樁基施工階段對樁底溶洞進行詳細勘察。

在樁基施工之前，該場地利用3個鉆孔 Y63-1、Y64-2和 Y65-3進行了兩次跨孔地震波CT成像(圖10)。其反演的縱波速度剖面與鉆孔柱狀圖基本吻合，共同表明了該場地巖溶發(fā)育且分布不均，灰?guī)r與上覆第四紀沉積物不整合接觸。僅就鉆孔Y64-2，其在鉆進過程中遇到串珠狀溶洞，溶洞發(fā)育于32.9～34.9m、38.4～41.6m以及41.6～44.4m等多個深度處。

圖10 跨孔地震波CT成像(Liu， 2021)

圖11a為該工況下的IPDI值圖像，從圖中34.9m處進行的溶洞探測表示的IPDI曲線在2.42 ms處的峰值可以計算得到溶洞頂板距樁孔底3.7m，即位于38.6m深度處，這與圖10中所揭示的結(jié)果相近。而48.7m所代表的曲線沒有明顯峰值，表明48.7m下為完整基巖。雖然34.9m處進行的溶洞探測能夠揭示溶洞頂板，但其溶洞底板無法直接從探測結(jié)果曲線中讀出。

圖11 實測信號IPDI優(yōu)化

當我們仔細觀察信號記錄，不難發(fā)現(xiàn)在信號記錄尾端發(fā)生了零漂現(xiàn)象，即一段時間的信號記錄不再關于時間軸對稱(圖11b)。鑒于這一現(xiàn)象的存在，嘗試對經(jīng)過EEMD得到的重構(gòu)信號進行處理，即使重構(gòu)信號減去EEMD所得的最后一個殘余趨勢項，并再運用上文所述的灰色預測模型信號延拓方法，最終得到的結(jié)果如圖11c?？梢钥吹讲粌H左端端點的“飛逸”得到了抑制，而且右端原先表現(xiàn)為端點效應的4.09 ms處呈現(xiàn)了一個新的IPDI峰值。結(jié)合黏土填充的縱波波速進行計算，該峰值極可能對應圖10中41.6m處的溶洞底板。

通過這一實際工程的溶洞頂?shù)装逦恢锰綔y結(jié)果優(yōu)化分析，我們可以看到，當溶洞底板所對應的反射信號處于采樣時窗末端時，其可能被IPDI計算過程以及采集儀器導致的零漂現(xiàn)象端點效應所掩蓋。此時，將EEMD后得到的重構(gòu)信號減去分解過程中自然得到的殘余趨勢項(趨勢項表達了信號零漂的具體趨勢)，再運用灰色預測模型進行端點效應處理，能夠有效將被掩蓋的溶洞底板反射信號提取出來，從而揭示出溶洞底板位置。

4 結(jié) 論

本文基于一種新型的樁孔底雙通道聲吶探測技術，開展了對于樁底溶洞探測中瞬時相位差分析過程優(yōu)化方法的研究。通過對數(shù)值仿真模擬信號的測試以及實測信號驗證，本文得到結(jié)論如下：

(1)當利用IPDI分析方法進行樁底溶洞探測時，設置采樣率為0.5～1MHz為宜，過小的采樣會影響最終對溶洞位置的判定； 當采樣率較小時，總體平均次數(shù)取300～600且標準差取0.1，這能夠增強IPDI方法對溶洞頂?shù)装逦恢糜嬎銣蚀_度。

(2)本文運用灰色預測模型，對重構(gòu)信號序列兩端的數(shù)據(jù)進行了預測延拓，彌補了求取信號瞬時相位運算方法中存在的不足，一定程度上抑制了端點效應。但這一方法可能影響溶洞頂?shù)装逦恢锰綔y精度，建議可以人為加長采樣時窗來進一步處理。

(3)信號傳感器余震累計或工程現(xiàn)場電磁干擾都可能導致信號記錄出現(xiàn)零漂現(xiàn)象。溶洞探測IPDI計算結(jié)果尾端可能會因此出現(xiàn)端點效應，從而埋藏需識別的溶洞底板反射信號。建議將重構(gòu)信號減去EEMD過程得到的殘余趨勢項，這將幫助重新準確識別溶洞底板位置并計算溶洞大小。