基于EEMD-HHT變換的爆破損傷分析方法

宋肖龍，高文學，季金銘，葉明班，張登杰

(北京工業(yè)大學城市建筑學部，北京，100124)

大多數(shù)工程巖體包含多種尺度的節(jié)理、斷層、裂隙和破碎帶的多級結構面，會影響巖體的完整性和穩(wěn)定性，特別是在隧道爆破掘進過程中，爆破振動作用會對開挖輪廓以外的圍巖造成一定程度的擾動，加劇其損傷演化過程，進而導致圍巖力學性質的劣化[1]。因此，在隧道施工過程中定位與識別圍巖中的裂隙等薄弱結構面，對圍巖損傷進行實時監(jiān)測并研究爆破開挖過程中圍巖的損傷演化特征至關重要，有助于保障隧道的施工安全[2]。地質雷達(GPR)探測是一種高分辨率的無損探測技術，對巖體進行掃描與識別的精度達到厘米級別，可有效甄別隧道圍巖的異常區(qū)域應用[3]。

地質雷達在探測過程中向巖體深處發(fā)射高頻電磁波，通過接收天線接收反射回來的電磁波得到雷達的回波信號，根據(jù)回波信號的波形、振幅等信息來推斷巖體中結構面的位置與形態(tài)分布[4]。但電磁波在傳播過程中會發(fā)生繞射、多次反射等現(xiàn)象，同時，雷達的回波信號也會遇到各種雜波的干擾，需要對其進行相應的數(shù)據(jù)處理從而得到可反映巖體內部損傷特征的雷達圖像。一般的地質雷達自帶分析軟件，經過多種數(shù)學運算得到巖體損傷圖像，但過于依賴操作者的經驗，雷達自帶分析軟件一般通過數(shù)字濾波來去除干擾信號，會丟失一定的信息。很多研究者開始使用一些新的方法來對地質雷達信號進行編譯和解讀，如神經網格算法[5?6]，BP成像算法[7?8]和復信號分析方法[9?10]等，極大地豐富了地質雷達信號分析方法，凌同華等[11]通過HHT 變換提取地質雷達原始信號的瞬時剖面特征，獲得了瞬時振幅、瞬時頻率等信息，并認為這些瞬時信息能更真實地反映地質雷達信號的原始特征。

本文作者采用EEMD-HHT 變換方法，對雷達信號進行處理。取其一階瞬時幅值作為特征參量來分析圍巖損傷，利用數(shù)值算例對多種形態(tài)的損傷進行模擬，得到了很好的效果，并在實際隧道中進行了探測，可有效評估圍巖的損傷狀態(tài)，為爆破作用下圍巖的損傷演化分析提供參考，也有助于優(yōu)化隧道掘進爆破參數(shù)設計，提高圍巖的安全穩(wěn)定性。

1 基于EEMD-HHT 方法的地質雷達損傷探測技術

電磁波在傳播過程中存在繞射和多次反射等問題，導致地質雷達的回波信號很難保證與目標體的原有形態(tài)一致[12?13]，而且由于地質雷達采用寬頻帶記錄方式，在實際的監(jiān)測過程中，地質雷達會將處于記錄頻帶范圍內的信號全部采集，其中包含真實反射信號以及干擾信號[14]，需要使用一定的信號處理方法獲得真實的傳播信息。地質雷達的回波信號經過EEMD-HHT 變換可濾除大量雜波干擾，且得到的瞬時信號參量可以很好地反映原信號的核心特征，為巖層損傷探測提供了理論參考。

1.1 地質雷達信號處理方法

1.1.1 EEMD-HHT分析方法

HHT 變換是一種自適應性較強的時頻局部化分析方法[15]，在對非平穩(wěn)信號的處理和分析上有很好的效果。HHT 變換由EMD 分解和Hilbert 變換2部分組成，其核心部分就是EMD 分解。但EMD分解存在模態(tài)混疊問題，導致在處理部分信號時會出現(xiàn)端點處嚴重失真，且對包含奇異事件的非平穩(wěn)信號分解效果較差，在處理信號時受到了限制，為此，李丹丹等[16-17]對EMD 分解進行了改進，提出了EEMD分解(總體評價經驗模態(tài)分解)，通過加入小幅度的白噪聲來阻止低頻模態(tài)分量的擴散，在一定程度上緩解了模態(tài)混疊和端點效應，進而得到信號的真實分量。

信號經過EEMD 分解之后可得到1 組IMF 分量，對每一階IMF 分量進行Hilbert 變換，用以構造解析信號z(t)：

式中：H[c(t)]為經過Hilbert變換后的IMF分量。

解析信號的極坐標形式可以反映Hilbert 變換的物理含義，它是通過1條正弦曲線的幅值調整獲得信號局部的最佳逼近[18]。可以進一步得到信號的瞬時幅值函數(shù)a(t)=和瞬時頻率函數(shù)ω(t)=

1.1.2 瞬時參量的選取及其物理意義

地質雷達的回波信號實質上是電磁波傳播沿途介質的反射系數(shù)R與雷達子波的卷積[19]，而R與介質相對介電常數(shù)有關。

式中：ε1為電磁波入射介質的相對介電常數(shù)；ε2為出射介質的相對介電常數(shù)。由此可以看出，當電磁波從高介電常數(shù)介質向低介電常數(shù)介質傳播時，反射系數(shù)為正，反之則為負。當高頻電磁波在致密巖體中傳播時，其相對介電常數(shù)差異較小，波形變化也就小，若反射回波發(fā)生較大波動，則說明電磁波在巖體內部遇到了相對介電常數(shù)差異較大的介質分界面。

圖1所示為一次地質雷達正演模擬的回波信號(已去除直達波影響)的幅值。由圖1可以看到：35 ns 前后信號幅值發(fā)生較大變化，說明電磁波在此時穿越了相對介電常數(shù)差異較大的不同介質。對該道信號進行EEMD-HHT變換后，可以得到多階瞬時數(shù)據(jù)。由于3階以上的數(shù)據(jù)失真嚴重，故提取1階和2階信號的瞬時數(shù)據(jù)，如圖2所示。

圖1 原雷達信號Fig.1 Original signal

由圖2可以看出：1 階IMF 分量較好地體現(xiàn)了原信號的變化形態(tài)；同樣，1 階瞬時振幅也比2 階瞬時振幅更能反映原信號的變化，并消除了原信號中小幅振蕩的變化波形，振幅最大值位置出現(xiàn)在變化區(qū)域的中心，1 階瞬時振幅的5 個極大值正好對應原信號的5個波峰或波谷，這說明1階瞬時振幅對應信號變化有很好的識別效果；對于瞬時頻率參量，原信號中變化區(qū)域頻率變化相對緩和，也有一定的識別效果，但在其他位置瞬時頻率振蕩嚴重，多道數(shù)據(jù)合并時彼此干擾嚴重，不利于后續(xù)分析。

圖2 EEMD-HHT變換瞬時參量圖Fig.2 Instantaneous parameter maps of EEMD-HHT transform

由此可以看出，瞬時振幅的實際意義更大，因為瞬時振幅越大，意味著雷達回波在該處經歷了較大幅度的振蕩，高頻電磁波可穿越相對介電常數(shù)差異明顯的介質，故可以使用1階瞬時振幅來識別巖體中的異常區(qū)域。

1.2 損傷判斷及閾值的選取

1.2.1 損傷評價指標

天然巖體中一般存在一定的先天缺陷，如裂隙、破碎帶等，相較完整巖體而言，這部分薄弱結構面會削弱巖體的承載能力，被視為損傷，在爆破振動作用的影響下，這類損傷會進一步演化和擴展。從微觀角度上看，巖石單元的損傷源于其內部的原生裂紋在拉應力作用下被激活并逐漸形成明顯的貫通型裂縫；從宏觀角度看，巖體內部裂隙以及破碎區(qū)域的擴展是無數(shù)微觀裂縫的發(fā)育和匯聚所導致的。裂縫影響區(qū)可視為巖體內部裂縫所占據(jù)的空間[20]。這些裂隙或破碎區(qū)域所形成的薄弱結構面決定了巖體的承載能力，因此，通過識別爆破掘進過程中圍巖的裂隙與破碎區(qū)域的演化規(guī)律，可獲得圍巖的爆破損傷演化特征。

地質雷達(GPR)使用高頻電磁波探測巖體內部，通過對雷達信號瞬時參量的提取可以有效識別密實巖體中的薄弱結構面的形態(tài)以及發(fā)育情況，通過對雷達圖像上瞬時振幅劇烈振蕩區(qū)域的數(shù)據(jù)進行提取，可以得到薄弱結構面的影響面積Sdam，完全致密巖體被視為無損傷，則有巖體損傷D=其中Sall為觀測區(qū)域面積。當觀測區(qū)域固定時，通過薄弱結構面影響面積的變化，可以得到巖體損傷的演化規(guī)律。

1.2.2 損傷閾值的確定

通過對雷達信號進行EEMD-HHT 變換所得到的瞬時振幅只能反映該處雷達回波的振蕩程度，不能準確體現(xiàn)巖體內部的損傷形態(tài)，需要通過設定損傷閾值得到較為精確的損傷圖像。

使用GprMax 軟件對單條裂縫形態(tài)損傷進行了正演模擬，模型長×寬為5 m×10 m，裂縫長度為2 m，寬度為0.2 m。模型以及模擬GPR 掃描該模型獲得的雷達信號經EEMD-HHT 變換后得到的瞬時幅值如圖3所示，設置損傷閾值k，當該處雷達信號的瞬時振幅大于k時即認為該處被損傷。

圖3 單裂隙正演模擬圖像Fig.3 Damage images of forward modeling for single fracture

由圖3可以看出：當損傷閾值k過小時，裂縫端點的反射回波也被計入損傷中，而當k過大時，大量實際損傷的反射信號被濾除，損傷形態(tài)失真嚴重；當k取4 000 左右時，能夠得到較為清晰、準確的損傷圖像，其形態(tài)、位置、大小都與實際相符。因此，選擇合適的損傷閾值對于巖體損傷圖像的分析至關重要。進行多次正演模擬發(fā)現(xiàn)，當損傷閾值設定為最大幅值的20%～40%時可以得到較為準確的損傷圖像，經多次試驗后，將損傷閾值設定為最大幅值的30%，既能有效濾除尖端的反射雜波，又能有效表征損傷形態(tài)。

一般來說，即便是完整的巖石區(qū)域其相對節(jié)點常數(shù)也不會完全一樣，在地質雷達的探測過程中也會產生小范圍振蕩，即巖石瞬時振幅小于閾值可以認為其較為完整，同時，瞬時振幅最大的位置應為損傷區(qū)域的中心位置，一般并不是損傷程度最嚴重的位置，為了獲得相對準確的損傷分布圖像，可將數(shù)據(jù)進行二值化處理，超過損傷閾值的即視為損傷，否則視為完好區(qū)域。

2 地質雷達正演模擬的數(shù)值算例

由于高頻電磁波的繞射現(xiàn)象，一維的正演模擬方法不能用來合成二維圖像，對地質雷達的二維正演模擬技術進行研究成為了地質雷達理論研究的主要內容。

時域有限差分法(FDTD)是近年來發(fā)展起來的以物理光學原理為基礎的高頻電磁場模擬方法，它具有直接時域計算、節(jié)約存儲空間和計算時間、計算程序通用性強等優(yōu)點，在地質雷達的正演模擬中得到了廣泛應用[21]。采用GprMax 軟件進行正演模擬。GprMax 是一款基于FDTD 的探地雷達正演模擬工具，該軟件可以設置Higdon 吸收邊界和理想匹配層(PML)作為邊界條件，在假設模擬介質為線性各向同質均勻媒介的前提下，通過引入電性參數(shù)輔助Maxwell 方程組求解。使用GprMax 進行探地雷達數(shù)據(jù)的正演模擬已在國內外得到廣泛應用[22]。

為了驗證EEMD-HHT 變換處理后的瞬時分量能否提取原雷達信號中的有效特征信息，基于GprMax軟件對多種形態(tài)的損傷類型分別進行正演模擬和分析。

2.1 正演模擬模型的建立

對于隧道圍巖來說，損傷主要表現(xiàn)為裂隙的貫通與擴展，進而形成松動面影響圍巖的穩(wěn)定性，其中常見的是宏觀尺寸較大的主裂縫區(qū)域，伴隨錯綜復雜的細小裂隙與破碎帶。為驗證雷達信號的瞬時振幅參量能否有效地描述損傷區(qū)域形態(tài)特征，使用GprMax軟件對不同情況下的損傷形態(tài)依次進行了正演模擬，并對得到的雷達回波信號進行EEMD-HHT 變換，提取其瞬時振幅數(shù)據(jù)。按設定的損傷閾值對瞬時振幅進行二值化處理，超過損傷閾值的位置認為其損傷，否則認為該處圍巖完整。

花崗巖巖石模型區(qū)域設置有效探測深度為5 m，雷達沿測量方向移動距離為10 m，模擬剖面法測量方式，其中。由于是單測線測量，模型厚度為最小離散步長，左下角為坐標起點，在其內部繪制出巖石介質與損傷形式，損傷區(qū)域填充介質為空氣。上表面為自由界面并放置發(fā)射與接收天線，其余3 面設置為PML 邊界，激勵源為400 MHz 的Rick 子波，發(fā)射天線的起點位置為(0.15 m，5 m)，接收天線位置為(0.25 m，5 m)，空間步長設置為0.01 m，天線移動步長0.01 m，采樣時間設為80 ns。模型材料的相關電性指標見表1。

表1 相關介質的電性參數(shù)表Table 1 Related dielectric electrical parameter table

2.2 交叉主裂縫形態(tài)損傷

交叉主裂縫模型如圖4(a)所示，裂縫寬度為0.1 m，為體現(xiàn)地質雷達識別的準確性，裂縫長度設置為一長一短。圖4(b)所示為經過二值化處理后的雷達信號瞬時振幅圖像，圖中藍色區(qū)域為完整巖體，紅色區(qū)域為強反射區(qū)域，即雷達信號的瞬時振幅出現(xiàn)劇烈振蕩的區(qū)域，可視為損傷處。

從圖4可以看出：經過處理的雷達信號瞬時幅值圖像基本定位了損傷區(qū)域，而且對裂縫位置也有較高的識別度。

圖4 交叉裂隙的正演模擬圖像Fig.4 Damage images of forward modeling for cross fractures

2.3 復雜損傷形態(tài)的正演模擬

2.3.1 網狀裂縫群形態(tài)損傷

圖5(a)所示為網狀裂縫群形態(tài)損傷的模型圖。為驗證地質雷達對復雜形態(tài)裂縫群的識別效果，設置了3 組平行裂縫彼此相交，構成了網狀裂縫群。裂縫寬度設置為0.1 m，損傷區(qū)域涉及模型的大部分區(qū)域。圖5(b)所示為經過處理后的瞬時振幅分量圖像。從圖5可以看出：反射雜波被基本濾除，基本可以反映出損傷區(qū)域的綜合分布和形態(tài)特征，但由于大量裂縫重疊交錯，雷達回波出現(xiàn)一定的信號淹沒現(xiàn)象，損傷圖像下半部分受到干擾較多，有一定失真。

圖5 網狀裂縫群的正演模擬圖像Fig.5 Damage images of forward modeling for map fractures

2.3.2 破碎帶區(qū)域形態(tài)損傷

為簡單起見，破碎帶形式設置為4排均勻分布且相同大小的36 個等邊三角形，三角形高為0.1 m，圖6所示為破碎帶形態(tài)損的正演模擬圖像。由圖6可以看出：經過EEMD-HHT 變換后得到瞬時幅值，再進行損傷閾值調節(jié)和二值化處理后，雜波基本上被去除，瞬時振幅基本上反映了破碎帶的分布與走向，與網狀裂縫群算例一樣，兩端的電磁波繞射現(xiàn)象明顯，淹沒了一部分反射信號。

圖6 破碎帶形態(tài)損的正演模擬圖像Fig.6 Damage images of forward modeling for fracture zone

從以上正演模擬算例可以看出，地質雷達信號在通過EEMD-HHT 變換后得到的瞬時振幅圖像，經閾值調節(jié)后，可以準確反映簡單形式的損傷分布形態(tài)；對于復雜形態(tài)的損傷，存在一定的信號淹沒現(xiàn)象，但對于損傷的識別和定位也有比較高的吻合度。

3 工程應用

為了進一步研究和驗證地質雷達損傷識別技術，在正演模擬的基礎上還進行了場地試驗，對正在開挖過程的隧道進行現(xiàn)場探測。

3.1 試驗場地基本情況

爆破試驗場地位于北京市延慶區(qū)，屬于北京冬奧會造雪引水系統(tǒng)工程，進行爆破開挖的隧道為一泄洪洞，隧道洞口海拔1 013.310 m，全長113 m，高6 m，寬5 m，采用全斷面開挖方式，藥量為96 kg，炸藥單耗為0.71 kg/m3。

使用地質雷達在爆破掘進的間隙對側壁圍巖進行探測，現(xiàn)場使用的是LTD 系列探地雷達，配以400 MHz天線，采樣頻率為1 024 MHz，測量方式為測距輪測距，可以有效探測圍巖40 m 范圍內的巖體損傷特征，采集時窗設為200 ns，探測深度約為10 m。每一次爆破前后使用地質雷達對右邊墻進行掃描，探測過程中采用1 條縱線的布線方式，測線位置距離地面1.5 m。探測區(qū)域自距掌子面5 m 處起始，向洞口方向測量，測線長度40 m，為便于后續(xù)的分析而將其分為寬度為10 m 的4 個區(qū)域，如圖7所示，炮孔布置圖如圖8所示。

圖7 地質雷達測線布置平面圖Fig.7 Line plan of GPR measurement

圖8 炮孔布置圖Fig.8 Blast hole distribution map

3.2 圍巖損傷分析

對地質雷達探測的數(shù)據(jù)進行處理后，得到圍巖的爆破損傷圖像。圖9所示為圍巖中3個區(qū)域的損傷圖像，圖像左端靠近掌子面，右端靠近洞口，探測寬度×深度為10 m×10 m。使用圖形處理軟件可以識別到圖像中高亮區(qū)域(即損傷區(qū)域)并得到相關的面積與寬度，所有相互獨立且未聯(lián)通的高亮區(qū)域均可視為一處損傷，選取總損傷面積、平均損傷面積、最大損傷寬度和平均損傷寬度共4個特征參數(shù)作為研究指標。其中損傷總面積為所有高亮區(qū)域的面積的總和，平均損傷面積和平均損傷寬度分布表示所有未聯(lián)通區(qū)域的平均面積和平均直徑，最大損傷寬度為未聯(lián)通區(qū)域的最大直徑，損傷總面積反映了觀測范圍內裂縫影響區(qū)域的變化，平均損傷面積和平均損傷寬度可以反映裂縫及破碎區(qū)域擴展情況，最大損傷寬度可以反映貫通性裂縫的演化程度，具體數(shù)據(jù)如表2～3所示。

1～4號區(qū)域中心位置距離爆源越來越大，距離爆源最近的1號區(qū)域中心位置距掌子面10 m，每個區(qū)域寬10 m。從表2可以看出：經歷1次爆破循環(huán)后，4 個區(qū)域的損傷均有一定程度增加，2 號和3號圍巖損傷影響區(qū)域較小，說明該區(qū)域巖體較為致密，其完整性較好，受到爆破振動效應影響后損傷演化現(xiàn)象并不嚴重，相較而言，1 號與4 號都有明顯的損傷擴展現(xiàn)象。

表2 圖像損傷面積特征數(shù)據(jù)Table 2 Feature data of damage area

由圖9可見：靠近圍巖外輪廓線的區(qū)域損傷區(qū)域分布明顯，越靠近圍巖內部其完整性越好，損傷越?。坏? 號和4 號區(qū)域在距圍巖表面4～6 m 范圍內都出現(xiàn)了大面積的損傷區(qū)域，并在爆破振動效應的擾動下?lián)p傷進一步擴展，但兩者又表現(xiàn)出不同的損傷演化特征。

表3 圖像損傷寬度特征數(shù)據(jù)Table 3 Feature data of damage width

1號區(qū)域距離爆源最近，該區(qū)域中心在爆破前距離掌子面10 m，該區(qū)域巖體受到爆破荷載影響最深，從損傷圖像和損傷特征來看，在經歷了1次循環(huán)爆破作用后，產生了較大的新生損傷區(qū)域，損傷平均寬度和最大寬度都有一定增加，說明損傷區(qū)域演化明顯，但損傷總面積和平均面積都減小。這是由于在爆破振動效應影響下，1號區(qū)域產生了大量新生的微小損傷區(qū)域，這些損傷區(qū)域并未連接到一起，因此平均損傷減小，而受到巖體自重效應和巖體內部應力重分布影響，圖9(a)中的裂縫影響區(qū)域被壓緊、密實，在1 次爆破擾動后，損傷面積出現(xiàn)暫時的減小現(xiàn)象。

圖9 某次爆破循環(huán)前、后圍巖損傷圖像Fig.9 Images of surrounding rock damage area before and after one cyclic blasting

4號區(qū)域中心位置距離爆源30 m，比1號區(qū)域經歷過更多次的循環(huán)爆破作用，因此，該區(qū)域靠近隧道外輪廓面的位置，即圖像上部位置，其爆前的損傷分布明顯。在1 次爆破振動效應影響下，損傷區(qū)域有著明顯擴展，結合表2的數(shù)據(jù)來看，損傷面積和寬度都有著不同程度的增加，但爆前圖像上的細小損傷在爆破后有部分貫通、連接在一起，有部分細小損傷同樣受到巖體自重效應影響被壓緊，未被地質雷達探測到，在爆破作用影響下，本身存在的主裂縫群和破碎區(qū)域在這個過程中被拉斷、破壞，有了明顯的擴展與貫通跡象，需要在后續(xù)爆破開挖過程中對此區(qū)域進行技術監(jiān)測與支護。

值得注意的是，1號區(qū)域在爆破振動效應影響下，產生了很多微小的損傷區(qū)域，而4號區(qū)出現(xiàn)更多的是貫通型損傷，即在距離爆源較遠區(qū)域，在爆破振動效應影響下，已有的損傷進一步演化，但不會產生新生損傷。

可以看出，地質雷達的探測結果在經過EEMD-HHT 變換后提取的瞬時振幅數(shù)據(jù)能很好地反映出隧道圍巖受爆破振動效應影響下的損傷演化趨勢，基于EEMD-HHT 變換的地質雷達損傷識別技術可以在隧道的開挖過程中對圍巖進行實時探測和損傷評價，對振動損傷演化劇烈的圍巖進行及時支護，對保障隧道施工安全有著重要意義。

4 結論

1)雷達信號經EEMD 分解后獲得的瞬時振幅分量可以有效地反映出電磁波在巖體傳播過程中穿越裂隙、破碎帶等形式的結構面，進而可對巖體內部的損傷區(qū)域進行識別與定位。

2)基于時域有限差分法(FDTD)對多種形式的圍巖損傷進行了正演模擬，雷達信號的一階瞬時振幅分量對損傷區(qū)域有很好的定位與識別效果。

3)在爆破振動效應影響下，觀測區(qū)圍巖的損傷演化受其完整性影響較大，內部含有一定破碎和裂縫區(qū)域損傷會發(fā)生進一步擴展，而完整性較好區(qū)域損傷演化不明顯，且距爆源較近的區(qū)域會產生新?lián)p傷，較遠區(qū)域多出現(xiàn)已有損傷連接與貫通。

4)利用EEMD-HHT 變換方法對圍巖的地質雷達損傷信號進行信號處理，操作簡單，可執(zhí)行度高，能有效去除干擾信息，形成較為直觀、準確的損傷圖像，豐富了巖體爆破損傷探測技術的理論體系，對于隧道施工安全進行有一定的參考意義。