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

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

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

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

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

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

0 引 言

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

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

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

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

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

但在計算信號IPDI值時，HHT方法中各個關鍵參數因影響規律不清而隨意選取以及IPDI圖像兩端可能出現的端點效應將會導致溶洞頂底板位置識別不夠精準，亟待優化處理。

本文根據數值仿真模擬信號測試結果，探究了雙通道樁底溶洞探測方法中IPDI分析過程的關鍵參數影響規律和端點效應抑制方法，提出優化建議，并結合工程現場實測信號驗證了優化方法的有效性，提升了溶洞頂底板位置識別的精準度。

1 溶洞聲吶探測方法及數值仿真

1.1 雙通道樁底溶洞聲吶探測系統

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

圖2 雙通道樁底溶洞聲吶探測系統示意圖

1.2 數值仿真

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

(1)

(2)

(3)

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

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

圖3 樁底溶洞聲吶探測數值模擬模型

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

表1 數值模擬的材料參數

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

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

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

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

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

1.3 IPDI計算方法

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

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

(4)

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

(5)

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

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

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

(6)

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

(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組一定幅值的白噪聲來消除可能出現的模態混疊現象，而所加入白噪聲的具體組數以及其相應幅值是人為決定的。這將可能影響提取出來相位信息的真實程度。關于這一方面的討論將會在下一節詳細展開。

2 溶洞探測精度影響因素分析及優化

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

2.1 采樣率分析及參數優化

現實生活中我們很難記錄下連續信號，因此只能通過一定的采樣率將連續信號時間序列離散化(Flandrin et al.，2004)。采樣率越高，越能使記錄下的離散信號逼近原始的連續信號。而EEMD分解通過對信號的極值點進行插值計算來形成上下包絡線，并以上下包絡線的均值作為各階IMF，因此過小的采樣率可能導致信號極值點的位置和取值發生誤差，最終使獲得的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。結合基巖波速、黏土填充波速以及式(2)可以算得溶洞頂板距樁孔距離2.928m，溶洞直徑2.02m，誤差均小于5%。但當采樣率為2MHz時，其計算所得溶洞頂底板位置精度為未得到明顯提升但運算所需時間大幅增長。

圖5 信號采樣率對IPDI計算結果影響

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

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

2.2 總體平均次數以及標準差參數分析及優化

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

圖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時，其計算結果最準確(溶洞直徑2m，溶洞頂板距樁孔底距離2.912m)。

圖6 EEMD參數取值對IPDI計算結果影響

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

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

(11)

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

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

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

圖7 端點效應抑制

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

本節嘗試在構造完成的重構信號基礎上，在其數據兩端分別拾取3m、3n個數據點，建立灰色預測模型，再各自向兩端外預測m以及n個數據點，得到新的重構信號，經HT變換后再選擇與原始信號相同時間段的重構信號數據點，最終計算得到IPDI值。如圖7，可以看到經灰色預測模型延拓原信號數據點后，其IPDI曲線末端的“飛逸”得到明顯改善，證明了此方法對改善端點效應的有效性。但同時可以發現經端點效應優化處理后其峰值位置出現了一定偏移，可能影響最終溶洞頂底板位置的判定。

3 IPDI優化方法應用

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

3.1 數值模擬工況

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

圖8 數值模擬工況示意圖

圖9 數值模擬工況信號IPDI值優化結果

3.2 現場實際信號

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

在樁基施工之前，該場地利用3個鉆孔 Y63-1、Y64-2和 Y65-3進行了兩次跨孔地震波CT成像(圖10)。其反演的縱波速度剖面與鉆孔柱狀圖基本吻合，共同表明了該場地巖溶發育且分布不均，灰巖與上覆第四紀沉積物不整合接觸。僅就鉆孔Y64-2，其在鉆進過程中遇到串珠狀溶洞，溶洞發育于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中所揭示的結果相近。而48.7m所代表的曲線沒有明顯峰值，表明48.7m下為完整基巖。雖然34.9m處進行的溶洞探測能夠揭示溶洞頂板，但其溶洞底板無法直接從探測結果曲線中讀出。

圖11 實測信號IPDI優化

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

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

4 結 論

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

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

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

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