第 31 卷 第 7 期
岩石力学与工程学报 Vol.31
No.7
2012 年 7 月 Chinese Journal of Rock Mechanics and Engineering July，2012 地质雷达隧道超前地质预报检测信号的
HHT 分析法
凌同华，张 胜，李升冉
(长沙理工大学 土木与建筑学院，湖南 长沙 410114)
摘要：隧道工程建设因开挖而诱发的各种地质灾害时有发生，为确保隧道施工安全，目前普遍运用超前地质预报 探明掌子面前方的地质情况，如何提高超前地质预报的判释水平显得尤为重要。针对隧道超前地质预报地质雷达 信号的特点，采用 HHT 分析法提取原始图像各 IMF 分量的瞬时剖面图及其单道信号的瞬时幅值、瞬时相位、瞬 时频率。结果表明，HHT 分析法分解出的瞬时参数突显了原始信号的特征信息，能较好地定量识别不良地质体位 置与形状，提高地质雷达原始图像解释的精度和可信程度。 关键词：隧道工程；地质预报；地质雷达；EMD 方法；HHT 分析
中图分类号：U 45 文献标识码：A 文章编号：1000–6915(2012)07–1422–07
HILBERT-HUANG TRANSFORM METHOD FOR DETECTION SIGNAL OF TUNNEL GEOLOGICAL PREDICTION USING GROUND PENETRATING
RADAR
LING Tonghua，ZHANG Sheng，LI Shengran
(School of Civil Engineering and Architecture，Changsha University of Science and Technology，Changsha，Hunan 410114，China)
Abstract ： Geological hazard accidents induced by excavation occur frequently in the course of tunnel constructions. In order to ensure tunnel construction safety， geological prediction for exploring potential geological conditions in front of tunnel face is widely adopted. Therefore，it is of vital importance to explore how to improve the interpretation level of geological prediction. Based on the characteristics of ground penetrating radar(GPR) signals in tunnel geological prediction，the Hilbert-Huang transform(HHT) method is used to extract instantaneous profiles of each IMF(intrinsic mode function) component in the original image and instantaneous parameters of the single trace signal，such as instantaneous amplitude，instantaneous phase and instantaneous frequency. The results show that instantaneous parameters decomposed by HHT method highlight the characteristics of the original signal，which can qualitatively identify the position and shape of catastrophic geological body，and improve the accuracy and credibility of the original image interpretation of ground penetrating radar.
Key words：tunnelling engineering；geological prediction；ground penetrating radar；EMD method；HHT analysis
杂的地质区域，地质条件对山岭隧道和城市地下工
1 引 言
随着我国高速公路及铁路建设逐渐进入更加复
程施工的影响日益突出。在隧道施工过程中，不良 地质和施工地质灾害常常成为制约隧道建设的主要 因素[1]。隧道施工掌子面前方的地质情况预报是隧
收稿日期：2012–03–26；修回日期：2012–05–09 基金项目：国家自然科学基金资助项目(51078043)
作者简介：凌同华(1968–)，男，博士，2004 年于中南大学资源与安全工程学院采矿工程专业获博士学位，现任教授、博士生导师，主要从事岩土工 程及灾害控制技术方面的教学与研究工作。E-mail：lingtonghua@163.com
第 31 卷 第 7 期 凌同华等：地质雷达隧道超前地质预报检测信号的 HHT 分析法 · 1423 ·
道工程界普遍关注而又没有得到很好解决的难题[2]。 开展隧道超前地质预报，进行隧道信息化施工，对 减少施工盲目性、确保工程安全有着重要意义 。
[3]
decomposition)方法和 Hilbert 变换 2 部分组成。HHT
方法与建立在先验性谐波基函数上的傅里叶变换和 建立在小波基函数上的小波变换不同，它依据数据 本身的时间尺度特性进行分解，是当前国际上公认 的最新的时频局部化分析方法之一[13]。该方法消除 了人为因素的影响，克服了传统方法中用毫无意义 的谐波分量来表示非平稳、非线性信号的缺陷，并 可得到极高的时频分辨率，具有良好的时频聚集性。 自 1998 年 HHT 法被提出以来，已被广泛应用于潮 汐、地震信号、结构分析、语音信号、故障诊断、 图像处理等领域的信号分析与处理。 2.1 EMD 方法
EMD 方法具有多尺度分析的特点，信号分解 时，不需要构造基函数，其基函数由数据自身构造， 是一种特殊的自适应小波分解方法。EMD 分解的基 本原理是将信号分解为一系列 IMF 分量。IMF 分 量是满足单分量信号物理解释的一类信号，和简单 的单调函数相比，一个 IMF 分量代表了一个简单的 振动模态。EMD 分解方法利用局部极大值和极小 值的包络来进行分解，一旦获得所有极值点，所有 的局部极大值用三次样条插值函数连接起来形成上 包络线，同样，所有的局部极小值连接起来形成下 包络线。
最终原始信号 x(t)可由 n 个 IMF 分量 ci(t)和残 余函数 rn(t)构成[14]，即
x(t)??2.2 Hilbert 变换
原始信号 x(t)经 EMD 方法分解后得到多个 IMF 分量的组合，分别对每个 IMF 分量 ci(t)进行 Hilbert 变换[15]，ci(t)变换后可表示为
1?? ci (? )
c?i (t)?? d? (2)
π?? t?????
近年来，超前地质预报在我国隧道施工中已经开始 发挥作用，特别在一些长、大隧道施工中，超前地 质预报已被列为必不可少的技术环节[4]。
隧道施工超前地质预报方法有地质分析法、超 前钻探法、地震反射波法、地质雷达法、瞬变电磁 法、陆地声纳法、地震波层析成像法等，每种方法 都有各自的优、缺点及适用范围[5-6]。地质雷达法是 近年来应用于短距离超前地质预报较为普遍的一种 物探手段。地质雷达法具有对隧道掌子面实施快速、 连续扫描、结果显示直观的优点。目前，地质雷达 信号处理方法普遍采用的是基于传统的傅里叶分析 理论的仪器自带的分析软件，而傅里叶分析理论没 有时频局部化分析功能，无法对掌子面前方的地质 缺陷所处的空间位置及其几何尺寸进行精确计算， 严重影响了地质雷达法进行超前地质预报的实际应 用效果，探索与之相适应的科学信号解释技术成为 广大研究者共同关注的问题。
D. Gabor[7]在 1946 年首先提出复信号概念，之 后，许多学者试图将复信号分析技术引入到地质雷 达处理领域。该方法在将地质雷达记录的实信号分 解出多个相互独立的参数的过程中，放大了记录中 的噪声，使得多参数剖面中的大量有用信息被干扰， 从而无法从根本上克服地质雷达抗干扰能力弱的
[8]
?? c (t)?? r (t) ini??1
n
(1)
缺点 。小波(包)分析具备时域–频域双重良好的局 部性特点，是信号处理和分析的好工具，然而，不 同的小波基具有不同的时频特征，用不同的小波基 分析同一个信号会产生不同的结果[9]。相对小波(包) 等分析方法而言，Hilbert-Huang 变换是一种全新而 优越的具有自适应的时频分析方法[10-12]。本文针对隧 道施工超前地质预报地质雷达时间–深度剖面图及 单道信号的特点，采用 HHT(Hilbert-Huang transform) 方法提取各 IMF(intrinsic mode function)分量的瞬时 剖面图，并结合单道信号的瞬时幅值、瞬时相位、 瞬时频率等参数，对地质雷达信号进行处理与分 析，从而达到对地质雷达原始图像进行定量解释和 表征的目的。
??构造解析信号 zi(t)：
zi (t)?? ci (t)?? jci? ( t)?? ai (t)e j?i (t )
(3)
显然有：
ai (t)?? c 2 (it(4) ) ?? c?i 2 (t) ?i (t)?? arctan c?i (t) ci (t)
(5)
2 HHT 方法
HHT 方法是近年提出的适合处理稳态信号和 非平稳信号的有效方法，它由 EMD(empirical mode
式中：ai(t)为信号 ci(t)的瞬时幅值，φi(t)为信号 ci(t) 的瞬时相位。根据式(4)可知 ai(t)≥ci(t)，这就意味 着|ci(t)|的曲线包络于 ai(t)的曲线。
在此基础上，将信号 x(t)的瞬时频率定义为
· 1424 · 岩石力学与工程学报 2012 年 1 1 d?i (t)
wi (t)?? (6) fi (t)???2π 2π dt
采用一条横线的布线方式，如图 1(a)所示。根据现 场情况灵活布置测线，其原则就是尽可能靠近掌子 面轴心位置，使测线距离尽可能长，尽可能多地采 集数据，以备后期数据的分析处理。具体测试过程 见图 1(b)。
式中： wi (t) 为圆频率。
对分解后得到的各 IMF 分量进行 Hilbert 变换 后，则可把原始信号表示成：
x(t)?? Re? a (t) e?? Re? a (t) e???nn
i j?i (t )ij wi (t )dt
i??1 i??1
(7)
测线
式(7)中省略了残余函数 rn(t)，Re 表示实部。 式(7)称为 HilbertHilbert 谱，记为 H (w， t)时频谱，简称?? Re ai (t) e ???n
j wi (t )dt
i??1
?
(8)
测线
掘进方向
3 隧道掌子面超前地质预报
3.1 工程概况
将军山隧道为清(远)连(州)一级公路升级改造 (高速)工程，位于连州—凤埠段的连南改线段(左半 幅)，进口位于广东省连南县龙口村东、开丫山北侧 斜坡地带，出口位于老广村，为单洞隧道。隧道起 点里程桩号为 ZK2112+285，地面标高 161.18 m， 设计高程 157.27 m，底板埋深 3.91 m；终点里程桩 号为 ZK2113+615，地面标高 199.10 m，设计高程 192.65 m，底板埋深 6.45 m，全长 1 330 m。隧道区 大的地貌属碳酸盐岩溶蚀低山丘陵地貌，微地貌属 于较坚硬碳酸盐岩类丘陵地貌和松散岩类沟谷地 貌，总体呈南东向展布。隧道出口端桩号 ZK2113+ 515 附近围岩为强中风化灰岩、白云质灰岩，岩体 破碎，岩溶、裂隙发育，地表溶蚀现象强烈，岩体 无自稳能力，拱部和侧壁易坍塌，围岩等级为 IV 级。
3.2 仪器参数的设置
在本次隧道施工超前地质预报过程中，现场采 用 RAMAC/GPR 系列地质雷达，配以 100 MHz 天 线。按照奈圭斯特采样定理，采样频率至少应为反 射波最高频率的 2 倍，但在频率比仅仅为 2 时，雷 达信号失真很明显。地质雷达在实际应用过程中， 为了达到良好的探测效果，其采样频率应该是天线 主频率的 10 倍以上，即 100 MHz 天线应采取的采 样频率为 1 024 MHz [16]。记录时间和叠加次数根据 天线中心频率相应选取，还可参照地质雷达野外工 作参数选择的基本原则[17]。 3.3 测线布置
由于隧道掌子面采用上、下台阶法开挖，工作 面很狭小，根据掌子面的具体情况，在检测过程中
(a) 测线布置 (b) 现场测试
图 1 地质雷达测线布置及现场测试图 Fig.1 Layout of GPR lines and field test
4 超前地质预报信号分析的 HHT 法
以将军山隧道掌子面桩号 ZK2113+515 的地质 雷达探测结果分析为例，对地质雷达自带分析软件 与 HHT 分析方法得出的结果进行比较，分析了自 带分析软件存在的问题，说明了 HHT 分析的特性。 4.1 采用自带软件对雷达图像进行分析
采用地质雷达分析软件 REFLEXW5.0 对原始 雷达图像(见图 2(a))进行一系列的数字化信号处理， 包括：漂移去除、零线设定、背景去噪、自动增益、 谱值平衡、道间平衡、滑动平均、混波处理、单道 漂移去除等，其滤波后的地质雷达图像解释结果见 图 2(b)。图 2 中，左半部分为剖面图，右半部分为 单道信号。从图 2(b)中可以看出，掌子面左半部分 前方从左侧深 6 m 至掌子面中心处深 1 m 范围有一 个斜分界面，反射强烈，反射波强度大，其中，前 方深 4～6 m 处有较大溶沟或破碎带，溶沟或破碎 带内有泥或软弱岩体，溶沟或破碎带至 20 m 范围 内岩体很破碎；掌子面右半部分前方 20 m 范围内 岩体破碎，岩溶、裂隙发育。然而，地质雷达自带 分析软件普遍基于传统的傅里叶分析理论，而傅里 叶分析理论没有时频局部化分析功能，无法对掌子 面前方的地质缺陷所处空间位置及其几何尺寸进行 精确计算。
4.2 超前地质预报信号的 HHT 分析
HHT 方法能对非平稳、非线性信号进行平稳 化、线性化处理，并在分解的过程中保留原始信号 的固有特性。因此，在 Matlab 语言平台上采用 EMD
第 31 卷 第 7 期 凌同华等：地质雷达隧道超前地质预报检测信号的 HHT 分析法 · 1425 ·
距离/m
瞬时振幅无关这一特点，对追踪连续性差的弱反射 波及极性变化的反射波独具优点。图 3，4 中，椭圆 内为溶洞的具体位置。从图 3，4 可以看出，随着探
测剖面上时间或探测深度的增加，地质雷达信号第 一阶 IMF 分量剖面衰减缓慢，而其他低频 IMF 分 量剖面的有效成份所剩无几，因此，原始雷达图像 的特征信号主要集中在第一阶 IMF 分量内。
拾取图 2(a)中第 30，50 道信号，参照图 3 和 4， 对探测结果作进一步分析。图 2(a)中，第 30，50 道 信号的瞬时幅值、瞬时相位、瞬时频率分别见图 5， 6。图 5，6 中，箭头 B 所示为发射和接受天线直接 耦合波以及通过地表传播的电磁波叠加位置，箭头 A 为溶洞散射信号位置，X 和 Y 分别为拾取信号的 横、纵坐标值。从图 5，6 中可以看出：(1) 在竖向 采样点 105～158 ns 内，反射界面明显，反射面瞬
时间/ns 深度/m
(a) 地质雷达探测剖面图
距离/m
时幅值显著增强，能量衰减快；(2) 图 5(b)，6(b) 中，瞬时相位波形受掌子面前方溶洞的物性变化较 为敏感，表现得较为紊乱；(3) 第 30，50 道信号初 至值分别位于竖向采样点 105，84 ns (见图 5(c)， 6(c))，该值即为溶洞与隧道围岩的第一界面；(4) 溶 洞中心对应的瞬时幅值最大，向两边的瞬时幅值逐 渐减小。以不同的比例显示，瞬时幅值剖面上显示 出不同宽度(成像宽度)的溶洞，这实际上是由于各 道信号瞬时幅值变化所引起的，即成像宽度与地质
深度/m 时间/ns 雷达瞬时剖面上各道信号的瞬时幅值有关，这与 溶洞成像宽度物理模型试验结果相吻合[18]，这也 说明溶洞信号瞬时幅值成像宽度可以有效表征溶洞 的真实宽度。
根据上面的分析，初步判断溶洞所属道号的大 致位置(第 20～70 道)后，再从地质雷达信号第一阶 IMF 分量瞬时参数谱中拾取各道信号初至值、瞬时 幅值成像宽度，如表 1 和图 7 所示。图 7(a)为采样 时间、道号与信号瞬时幅值的三维图，即溶洞位置 信号能量衰减图；图 7(b)为溶洞信号瞬时幅值的成 像宽度，其横坐标拾取信号第 20～70 道，纵坐标拾 取信号竖向采样点 60～160 ns。从表 1 和图 7 可以 看出：(1) 随着地质雷达信号道号的增加，各道信 号初至值离隧道掌子面的距离越近，即掌子面开挖 越先暴露出溶洞；(2) 溶洞信号瞬时幅值成像宽度 均值约 52 ns，即溶洞瞬时幅值成像宽度为 2.5 m。 据此可以判断溶洞为掌子面从左至右 1.5～5.3 m(即 信号第 20～70 道所在位置)、水平向前宽度为 2.5 m。 因此，相比建立在傅里叶变换基础上的地质雷达自 带分析软件，HHT 方法能对地质雷达原始图像进行 定量解释和有效表征。
(b) 解释结果
图 2 地质雷达探测剖面图及解释结果
Fig.2 GPR detection profiles and explanation results
方法对隧道超前地质预报地质雷达信号进行分解得 到 n 个 IMF 分量，每个 IMF 分量都包含了不同的 特征尺度信息，通过 EMD 分解，原始信号的特征 完全可以由这 n 个 IMF 分量来表征。图 3(a)，(b) 和(c)分别为原始图像分解得到的前三阶 IMF 分量 的瞬时幅值剖面图，从图 3 中可以看出地质雷达信 号能量的瞬时变化情况，可以用来判断地质体有关 岩性状况。图 4(a)，4(b)和(c)分别为原始图像分解 得到的前三阶 IMF 分量的瞬时相位剖面图。利用与
