基于EEMD-HHT变换的爆破损伤分析方法

宋肖龙，高文学，季金铭，叶明班，张登杰

(北京工业大学城市建筑学部，北京，100124)

大多数工程岩体包含多种尺度的节理、断层、裂隙和破碎带的多级结构面，会影响岩体的完整性和稳定性，特别是在隧道爆破掘进过程中，爆破振动作用会对开挖轮廓以外的围岩造成一定程度的扰动，加剧其损伤演化过程，进而导致围岩力学性质的劣化[1]。因此，在隧道施工过程中定位与识别围岩中的裂隙等薄弱结构面，对围岩损伤进行实时监测并研究爆破开挖过程中围岩的损伤演化特征至关重要，有助于保障隧道的施工安全[2]。地质雷达(GPR)探测是一种高分辨率的无损探测技术，对岩体进行扫描与识别的精度达到厘米级别，可有效甄别隧道围岩的异常区域应用[3]。

地质雷达在探测过程中向岩体深处发射高频电磁波，通过接收天线接收反射回来的电磁波得到雷达的回波信号，根据回波信号的波形、振幅等信息来推断岩体中结构面的位置与形态分布[4]。但电磁波在传播过程中会发生绕射、多次反射等现象，同时，雷达的回波信号也会遇到各种杂波的干扰，需要对其进行相应的数据处理从而得到可反映岩体内部损伤特征的雷达图像。一般的地质雷达自带分析软件，经过多种数学运算得到岩体损伤图像，但过于依赖操作者的经验，雷达自带分析软件一般通过数字滤波来去除干扰信号，会丢失一定的信息。很多研究者开始使用一些新的方法来对地质雷达信号进行编译和解读，如神经网格算法[5−6]，BP成像算法[7−8]和复信号分析方法[9−10]等，极大地丰富了地质雷达信号分析方法，凌同华等[11]通过HHT 变换提取地质雷达原始信号的瞬时剖面特征，获得了瞬时振幅、瞬时频率等信息，并认为这些瞬时信息能更真实地反映地质雷达信号的原始特征。

本文作者采用EEMD-HHT 变换方法，对雷达信号进行处理。取其一阶瞬时幅值作为特征参量来分析围岩损伤，利用数值算例对多种形态的损伤进行模拟，得到了很好的效果，并在实际隧道中进行了探测，可有效评估围岩的损伤状态，为爆破作用下围岩的损伤演化分析提供参考，也有助于优化隧道掘进爆破参数设计，提高围岩的安全稳定性。

1 基于EEMD-HHT 方法的地质雷达损伤探测技术

电磁波在传播过程中存在绕射和多次反射等问题，导致地质雷达的回波信号很难保证与目标体的原有形态一致[12−13]，而且由于地质雷达采用宽频带记录方式，在实际的监测过程中，地质雷达会将处于记录频带范围内的信号全部采集，其中包含真实反射信号以及干扰信号[14]，需要使用一定的信号处理方法获得真实的传播信息。地质雷达的回波信号经过EEMD-HHT 变换可滤除大量杂波干扰，且得到的瞬时信号参量可以很好地反映原信号的核心特征，为岩层损伤探测提供了理论参考。

1.1 地质雷达信号处理方法

1.1.1 EEMD-HHT分析方法

HHT 变换是一种自适应性较强的时频局部化分析方法[15]，在对非平稳信号的处理和分析上有很好的效果。HHT 变换由EMD 分解和Hilbert 变换2部分组成，其核心部分就是EMD 分解。但EMD分解存在模态混叠问题，导致在处理部分信号时会出现端点处严重失真，且对包含奇异事件的非平稳信号分解效果较差，在处理信号时受到了限制，为此，李丹丹等[16-17]对EMD 分解进行了改进，提出了EEMD分解(总体评价经验模态分解)，通过加入小幅度的白噪声来阻止低频模态分量的扩散，在一定程度上缓解了模态混叠和端点效应，进而得到信号的真实分量。

信号经过EEMD 分解之后可得到1 组IMF 分量，对每一阶IMF 分量进行Hilbert 变换，用以构造解析信号z(t)：

式中：H[c(t)]为经过Hilbert变换后的IMF分量。

解析信号的极坐标形式可以反映Hilbert 变换的物理含义，它是通过1条正弦曲线的幅值调整获得信号局部的最佳逼近[18]。可以进一步得到信号的瞬时幅值函数a(t)=和瞬时频率函数ω(t)=

1.1.2 瞬时参量的选取及其物理意义

地质雷达的回波信号实质上是电磁波传播沿途介质的反射系数R与雷达子波的卷积[19]，而R与介质相对介电常数有关。

式中：ε1为电磁波入射介质的相对介电常数；ε2为出射介质的相对介电常数。由此可以看出，当电磁波从高介电常数介质向低介电常数介质传播时，反射系数为正，反之则为负。当高频电磁波在致密岩体中传播时，其相对介电常数差异较小，波形变化也就小，若反射回波发生较大波动，则说明电磁波在岩体内部遇到了相对介电常数差异较大的介质分界面。

图1所示为一次地质雷达正演模拟的回波信号(已去除直达波影响)的幅值。由图1可以看到：35 ns 前后信号幅值发生较大变化，说明电磁波在此时穿越了相对介电常数差异较大的不同介质。对该道信号进行EEMD-HHT变换后，可以得到多阶瞬时数据。由于3阶以上的数据失真严重，故提取1阶和2阶信号的瞬时数据，如图2所示。

图1 原雷达信号Fig.1 Original signal

由图2可以看出：1 阶IMF 分量较好地体现了原信号的变化形态；同样，1 阶瞬时振幅也比2 阶瞬时振幅更能反映原信号的变化，并消除了原信号中小幅振荡的变化波形，振幅最大值位置出现在变化区域的中心，1 阶瞬时振幅的5 个极大值正好对应原信号的5个波峰或波谷，这说明1阶瞬时振幅对应信号变化有很好的识别效果；对于瞬时频率参量，原信号中变化区域频率变化相对缓和，也有一定的识别效果，但在其他位置瞬时频率振荡严重，多道数据合并时彼此干扰严重，不利于后续分析。

图2 EEMD-HHT变换瞬时参量图Fig.2 Instantaneous parameter maps of EEMD-HHT transform

由此可以看出，瞬时振幅的实际意义更大，因为瞬时振幅越大，意味着雷达回波在该处经历了较大幅度的振荡，高频电磁波可穿越相对介电常数差异明显的介质，故可以使用1阶瞬时振幅来识别岩体中的异常区域。

1.2 损伤判断及阈值的选取

1.2.1 损伤评价指标

天然岩体中一般存在一定的先天缺陷，如裂隙、破碎带等，相较完整岩体而言，这部分薄弱结构面会削弱岩体的承载能力，被视为损伤，在爆破振动作用的影响下，这类损伤会进一步演化和扩展。从微观角度上看，岩石单元的损伤源于其内部的原生裂纹在拉应力作用下被激活并逐渐形成明显的贯通型裂缝；从宏观角度看，岩体内部裂隙以及破碎区域的扩展是无数微观裂缝的发育和汇聚所导致的。裂缝影响区可视为岩体内部裂缝所占据的空间[20]。这些裂隙或破碎区域所形成的薄弱结构面决定了岩体的承载能力，因此，通过识别爆破掘进过程中围岩的裂隙与破碎区域的演化规律，可获得围岩的爆破损伤演化特征。

地质雷达(GPR)使用高频电磁波探测岩体内部，通过对雷达信号瞬时参量的提取可以有效识别密实岩体中的薄弱结构面的形态以及发育情况，通过对雷达图像上瞬时振幅剧烈振荡区域的数据进行提取，可以得到薄弱结构面的影响面积Sdam，完全致密岩体被视为无损伤，则有岩体损伤D=其中Sall为观测区域面积。当观测区域固定时，通过薄弱结构面影响面积的变化，可以得到岩体损伤的演化规律。

1.2.2 损伤阈值的确定

通过对雷达信号进行EEMD-HHT 变换所得到的瞬时振幅只能反映该处雷达回波的振荡程度，不能准确体现岩体内部的损伤形态，需要通过设定损伤阈值得到较为精确的损伤图像。

使用GprMax 软件对单条裂缝形态损伤进行了正演模拟，模型长×宽为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过大时，大量实际损伤的反射信号被滤除，损伤形态失真严重；当k取4 000 左右时，能够得到较为清晰、准确的损伤图像，其形态、位置、大小都与实际相符。因此，选择合适的损伤阈值对于岩体损伤图像的分析至关重要。进行多次正演模拟发现，当损伤阈值设定为最大幅值的20%～40%时可以得到较为准确的损伤图像，经多次试验后，将损伤阈值设定为最大幅值的30%，既能有效滤除尖端的反射杂波，又能有效表征损伤形态。

一般来说，即便是完整的岩石区域其相对节点常数也不会完全一样，在地质雷达的探测过程中也会产生小范围振荡，即岩石瞬时振幅小于阈值可以认为其较为完整，同时，瞬时振幅最大的位置应为损伤区域的中心位置，一般并不是损伤程度最严重的位置，为了获得相对准确的损伤分布图像，可将数据进行二值化处理，超过损伤阈值的即视为损伤，否则视为完好区域。

2 地质雷达正演模拟的数值算例

由于高频电磁波的绕射现象，一维的正演模拟方法不能用来合成二维图像，对地质雷达的二维正演模拟技术进行研究成为了地质雷达理论研究的主要内容。

时域有限差分法(FDTD)是近年来发展起来的以物理光学原理为基础的高频电磁场模拟方法，它具有直接时域计算、节约存储空间和计算时间、计算程序通用性强等优点，在地质雷达的正演模拟中得到了广泛应用[21]。采用GprMax 软件进行正演模拟。GprMax 是一款基于FDTD 的探地雷达正演模拟工具，该软件可以设置Higdon 吸收边界和理想匹配层(PML)作为边界条件，在假设模拟介质为线性各向同质均匀媒介的前提下，通过引入电性参数辅助Maxwell 方程组求解。使用GprMax 进行探地雷达数据的正演模拟已在国内外得到广泛应用[22]。

为了验证EEMD-HHT 变换处理后的瞬时分量能否提取原雷达信号中的有效特征信息，基于GprMax软件对多种形态的损伤类型分别进行正演模拟和分析。

2.1 正演模拟模型的建立

对于隧道围岩来说，损伤主要表现为裂隙的贯通与扩展，进而形成松动面影响围岩的稳定性，其中常见的是宏观尺寸较大的主裂缝区域，伴随错综复杂的细小裂隙与破碎带。为验证雷达信号的瞬时振幅参量能否有效地描述损伤区域形态特征，使用GprMax软件对不同情况下的损伤形态依次进行了正演模拟，并对得到的雷达回波信号进行EEMD-HHT 变换，提取其瞬时振幅数据。按设定的损伤阈值对瞬时振幅进行二值化处理，超过损伤阈值的位置认为其损伤，否则认为该处围岩完整。

花岗岩岩石模型区域设置有效探测深度为5 m，雷达沿测量方向移动距离为10 m，模拟剖面法测量方式，其中。由于是单测线测量，模型厚度为最小离散步长，左下角为坐标起点，在其内部绘制出岩石介质与损伤形式，损伤区域填充介质为空气。上表面为自由界面并放置发射与接收天线，其余3 面设置为PML 边界，激励源为400 MHz 的Rick 子波，发射天线的起点位置为(0.15 m，5 m)，接收天线位置为(0.25 m，5 m)，空间步长设置为0.01 m，天线移动步长0.01 m，采样时间设为80 ns。模型材料的相关电性指标见表1。

表1 相关介质的电性参数表Table 1 Related dielectric electrical parameter table

2.2 交叉主裂缝形态损伤

交叉主裂缝模型如图4(a)所示，裂缝宽度为0.1 m，为体现地质雷达识别的准确性，裂缝长度设置为一长一短。图4(b)所示为经过二值化处理后的雷达信号瞬时振幅图像，图中蓝色区域为完整岩体，红色区域为强反射区域，即雷达信号的瞬时振幅出现剧烈振荡的区域，可视为损伤处。

从图4可以看出：经过处理的雷达信号瞬时幅值图像基本定位了损伤区域，而且对裂缝位置也有较高的识别度。

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

2.3 复杂损伤形态的正演模拟

2.3.1 网状裂缝群形态损伤

图5(a)所示为网状裂缝群形态损伤的模型图。为验证地质雷达对复杂形态裂缝群的识别效果，设置了3 组平行裂缝彼此相交，构成了网状裂缝群。裂缝宽度设置为0.1 m，损伤区域涉及模型的大部分区域。图5(b)所示为经过处理后的瞬时振幅分量图像。从图5可以看出：反射杂波被基本滤除，基本可以反映出损伤区域的综合分布和形态特征，但由于大量裂缝重叠交错，雷达回波出现一定的信号淹没现象，损伤图像下半部分受到干扰较多，有一定失真。

图5 网状裂缝群的正演模拟图像Fig.5 Damage images of forward modeling for map fractures

2.3.2 破碎带区域形态损伤

为简单起见，破碎带形式设置为4排均匀分布且相同大小的36 个等边三角形，三角形高为0.1 m，图6所示为破碎带形态损的正演模拟图像。由图6可以看出：经过EEMD-HHT 变换后得到瞬时幅值，再进行损伤阈值调节和二值化处理后，杂波基本上被去除，瞬时振幅基本上反映了破碎带的分布与走向，与网状裂缝群算例一样，两端的电磁波绕射现象明显，淹没了一部分反射信号。

图6 破碎带形态损的正演模拟图像Fig.6 Damage images of forward modeling for fracture zone

从以上正演模拟算例可以看出，地质雷达信号在通过EEMD-HHT 变换后得到的瞬时振幅图像，经阈值调节后，可以准确反映简单形式的损伤分布形态；对于复杂形态的损伤，存在一定的信号淹没现象，但对于损伤的识别和定位也有比较高的吻合度。

3 工程应用

为了进一步研究和验证地质雷达损伤识别技术，在正演模拟的基础上还进行了场地试验，对正在开挖过程的隧道进行现场探测。

3.1 试验场地基本情况

爆破试验场地位于北京市延庆区，属于北京冬奥会造雪引水系统工程，进行爆破开挖的隧道为一泄洪洞，隧道洞口海拔1 013.310 m，全长113 m，高6 m，宽5 m，采用全断面开挖方式，药量为96 kg，炸药单耗为0.71 kg/m3。

使用地质雷达在爆破掘进的间隙对侧壁围岩进行探测，现场使用的是LTD 系列探地雷达，配以400 MHz天线，采样频率为1 024 MHz，测量方式为测距轮测距，可以有效探测围岩40 m 范围内的岩体损伤特征，采集时窗设为200 ns，探测深度约为10 m。每一次爆破前后使用地质雷达对右边墙进行扫描，探测过程中采用1 条纵线的布线方式，测线位置距离地面1.5 m。探测区域自距掌子面5 m 处起始，向洞口方向测量，测线长度40 m，为便于后续的分析而将其分为宽度为10 m 的4 个区域，如图7所示，炮孔布置图如图8所示。

图7 地质雷达测线布置平面图Fig.7 Line plan of GPR measurement

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

3.2 围岩损伤分析

对地质雷达探测的数据进行处理后，得到围岩的爆破损伤图像。图9所示为围岩中3个区域的损伤图像，图像左端靠近掌子面，右端靠近洞口，探测宽度×深度为10 m×10 m。使用图形处理软件可以识别到图像中高亮区域(即损伤区域)并得到相关的面积与宽度，所有相互独立且未联通的高亮区域均可视为一处损伤，选取总损伤面积、平均损伤面积、最大损伤宽度和平均损伤宽度共4个特征参数作为研究指标。其中损伤总面积为所有高亮区域的面积的总和，平均损伤面积和平均损伤宽度分布表示所有未联通区域的平均面积和平均直径，最大损伤宽度为未联通区域的最大直径，损伤总面积反映了观测范围内裂缝影响区域的变化，平均损伤面积和平均损伤宽度可以反映裂缝及破碎区域扩展情况，最大损伤宽度可以反映贯通性裂缝的演化程度，具体数据如表2～3所示。

1～4号区域中心位置距离爆源越来越大，距离爆源最近的1号区域中心位置距掌子面10 m，每个区域宽10 m。从表2可以看出：经历1次爆破循环后，4 个区域的损伤均有一定程度增加，2 号和3号围岩损伤影响区域较小，说明该区域岩体较为致密，其完整性较好，受到爆破振动效应影响后损伤演化现象并不严重，相较而言，1 号与4 号都有明显的损伤扩展现象。

表2 图像损伤面积特征数据Table 2 Feature data of damage area

由图9可见：靠近围岩外轮廓线的区域损伤区域分布明显，越靠近围岩内部其完整性越好，损伤越小；但1 号和4 号区域在距围岩表面4～6 m 范围内都出现了大面积的损伤区域，并在爆破振动效应的扰动下损伤进一步扩展，但两者又表现出不同的损伤演化特征。

表3 图像损伤宽度特征数据Table 3 Feature data of damage width

1号区域距离爆源最近，该区域中心在爆破前距离掌子面10 m，该区域岩体受到爆破荷载影响最深，从损伤图像和损伤特征来看，在经历了1次循环爆破作用后，产生了较大的新生损伤区域，损伤平均宽度和最大宽度都有一定增加，说明损伤区域演化明显，但损伤总面积和平均面积都减小。这是由于在爆破振动效应影响下，1号区域产生了大量新生的微小损伤区域，这些损伤区域并未连接到一起，因此平均损伤减小，而受到岩体自重效应和岩体内部应力重分布影响，图9(a)中的裂缝影响区域被压紧、密实，在1 次爆破扰动后，损伤面积出现暂时的减小现象。

图9 某次爆破循环前、后围岩损伤图像Fig.9 Images of surrounding rock damage area before and after one cyclic blasting

4号区域中心位置距离爆源30 m，比1号区域经历过更多次的循环爆破作用，因此，该区域靠近隧道外轮廓面的位置，即图像上部位置，其爆前的损伤分布明显。在1 次爆破振动效应影响下，损伤区域有着明显扩展，结合表2的数据来看，损伤面积和宽度都有着不同程度的增加，但爆前图像上的细小损伤在爆破后有部分贯通、连接在一起，有部分细小损伤同样受到岩体自重效应影响被压紧，未被地质雷达探测到，在爆破作用影响下，本身存在的主裂缝群和破碎区域在这个过程中被拉断、破坏，有了明显的扩展与贯通迹象，需要在后续爆破开挖过程中对此区域进行技术监测与支护。

值得注意的是，1号区域在爆破振动效应影响下，产生了很多微小的损伤区域，而4号区出现更多的是贯通型损伤，即在距离爆源较远区域，在爆破振动效应影响下，已有的损伤进一步演化，但不会产生新生损伤。

可以看出，地质雷达的探测结果在经过EEMD-HHT 变换后提取的瞬时振幅数据能很好地反映出隧道围岩受爆破振动效应影响下的损伤演化趋势，基于EEMD-HHT 变换的地质雷达损伤识别技术可以在隧道的开挖过程中对围岩进行实时探测和损伤评价，对振动损伤演化剧烈的围岩进行及时支护，对保障隧道施工安全有着重要意义。

4 结论

1)雷达信号经EEMD 分解后获得的瞬时振幅分量可以有效地反映出电磁波在岩体传播过程中穿越裂隙、破碎带等形式的结构面，进而可对岩体内部的损伤区域进行识别与定位。

2)基于时域有限差分法(FDTD)对多种形式的围岩损伤进行了正演模拟，雷达信号的一阶瞬时振幅分量对损伤区域有很好的定位与识别效果。

3)在爆破振动效应影响下，观测区围岩的损伤演化受其完整性影响较大，内部含有一定破碎和裂缝区域损伤会发生进一步扩展，而完整性较好区域损伤演化不明显，且距爆源较近的区域会产生新损伤，较远区域多出现已有损伤连接与贯通。

4)利用EEMD-HHT 变换方法对围岩的地质雷达损伤信号进行信号处理，操作简单，可执行度高，能有效去除干扰信息，形成较为直观、准确的损伤图像，丰富了岩体爆破损伤探测技术的理论体系，对于隧道施工安全进行有一定的参考意义。