隧道施工不良地质探地雷达超前探测全波形反演研究

张崇民， 张凤凯， 李 尧

(山东大学岩土与结构工程研究中心， 山东 济南 250061)

0 引言

随着我国经济建设的发展，大量隧道工程破土动工，隧道工程施工中需要面临的地质状况越来越复杂。溶洞、岩体破碎带、岩体裂隙等不良地质极易引发突水突泥、坍方等重大灾害，并造成严重的经济损失和人员伤亡[1]。因此，提前了解隧道掌子面前方赋存的不良地质体状况对于施工安全具有极为重要的意义。

探地雷达(ground penetrating radar，GPR)基于向地下地质体发射高频电磁信号并接收异常地质体反射信号来进行勘查，具有无伤、高效、高分辨率和对水体敏感的优点[2]，在隧道短距离不良地质超前探测中具有较好的应用效果。近年来，许多学者针对探地雷达隧道探测开展了广泛的研究，并通过实际工程资料分析[3]、正演模拟[4]等总结了典型隧道不良地质体的响应特征和辨识方法，对于提高隧道不良地质超前探测的准确性起到了关键性作用。随着探地雷达数据处理方法的发展，对于探测资料的解释也由传统的雷达剖面图分析向偏移[5-6]、反演[7]、复信号分析[8-9]等多种方法发展，进一步提高了探地雷达数据分析的精度，为隧道探地雷达资料处理提供了可靠的思路。总的来说，探地雷达在隧道不良地质超前探测中获得了广泛应用，基于实际工程和正演模拟总结得到的不良地质体的响应特征，实现了对不良地质体的简单辨识。由于电磁波传播过程中存在多次反射波和绕射波，导致探地雷达数据中的同相轴走向与实际不良地质体的形态差异较大[4,8,10]，特别是地质状况复杂区域的雷达探测结果中同相轴交错叠加，仅凭经验对雷达剖面图分析无法获得异常体的真实形态信息。

基于上述问题，本文提出将全波形反演运用到探地雷达数据处理中进行不良地质体判断的方法。全波形反演(full waveform inversion,FWI)是一种基于全波场信息实现地下介质参数反演的高分辨率成像方法，最早由A. Tarantola[11]提出并应用于地震勘探领域，后来逐渐在雷达信号处理领域中应用发展。S. Busch等[12]实现了地面雷达数据的相对介电常数和电导率的定量估计； 吴俊军[13]、孟旭[14]对时间域全波形反演方法做了系统推导及扩展； 李尧[15]将跨孔雷达运用到隧道施工不良地质超前探测中，并将全波形反演和逆时偏移成像方法相结合起来。总体来说，全波形反演方法能够较为准确地反演出地下介质的介电常数信息，在隧道掌子面前方不良地质体定位和识别方面显示出了良好的应用前景。

本文将全波形反演方法引入探地雷达探测数据处理流程，针对典型不良地质体进行探地雷达探测正演模拟和全波形反演研究，获得探地雷达探测典型不良地质体的响应特征和反演成像规律，以期为实际工程应用提供理论基础。

1 全波形反演基本原理及数据处理流程

全波形反演是一种基于全波场信息反演地下介质的反演方法，反演结果可以指示目标探测区域的较准确相对介电常数分布，成像结果更为直观。全波形反演的思想是寻找使目标函数最小时的地下介质物性参数分布(本文搜索的是相对介电常数和电导率)。笔者采用的目标函数为正演模拟数据与实际观测数据之间的二范数[13-16]。

(1)

式中：E(ε,σ),Eobs分别为正演模拟数据和实际观测数据；s为探地雷达探测位置；t为探测时间。

当目标函数的值最小时，正演模拟数据和实际观测数据最接近，此时得到的相对介电常数分布也就越接近真实模型。本文使用共轭梯度法进行目标函数寻优，目标函数梯度可由式(2)和式(3)计算[13-16]。

(2)

(3)

式(2)和式(3)中：I=GT·R；R=E-Eobs(R为残差，G为格林算子，I为逆时传播波场数据，均由正演模拟获得)。

介电常数和电导率更新的迭代公式如式(4)和式(5)[13-16]所示。

εk+1=εk-ζε,k·gε。

(4)

σk+1=σk-ζσ,k·gσ。

(5)

式(4)和式(5)中：k为迭代次数；ζε,k和ζσ,k为第k次迭代中相对介电常数和电导率的迭代步长。

ζε,k和ζσ,k可分别表示为[13-16]：

(6)

(7)

式(6)和式(7)中：κε和κσ为小稳定因子，反演过程中必须选择适当的值，并且随着迭代更新。

根据上述梯度和迭代步长不断更新，直至达到终止条件，完成全波形反演。

探地雷达探测数据全波形反演处理过程，具体分为以下8个步骤[11,13]。

1) 建立初始模型；

2) 将源子波加载到发射天线，进行正演模拟，并记录全时刻正传波场和接收天线处的正演波形；

3) 将实际接收数据与正演数据相减得到残差E-Eobs和目标函数S；

4) 将残差E-Eobs在时间上倒置，并以此为源加载到接收天线上，计算逆时残场；

5) 通过正传波场和逆时残场计算相对介电常数迭代梯度gε和电导率迭代梯度gσ；

6) 计算对应于梯度方向的相对介电常数迭代步长ζε和电导率迭代步长ζσ；

7) 更新模型参数；

8) 计算目标函数并判断是否达到终止条件，如未达到终止条件则跳转至步骤2)，如已达到终止条件则完成全波形反演过程。

本文在基于全波形反演方法处理探地雷达正演模拟数据时，首先假设初始模型为均匀模型，使用全波形反演方法得到不良地质区域的相对介电常数。全波形反演流程如图1所示。

图1 全波形反演流程

2 隧道施工典型不良地质体探地雷达探测数值算例

隧道围岩地质情况复杂多变，各种不良地质造成的安全事故也屡见不鲜，隧道施工不良地质常见于溶洞、岩体裂隙和岩体破碎带等[17]。基于上述不良地质类型，本文以空气、灰岩、潮湿黏土、水4种介质设计了4个数值算例,进行正演模拟和全波形反演结果分析。介质电性参数设置如表1所示。

表1 介质电性参数设置[18]

本文探地雷达采用单发射单接收型平板天线，发射天线和接收天线收发距固定。发射电磁波中心频率为100 MHz，采样时间为365 ns。图2为探地雷达隧道检测示意图，左侧为探地雷达测线，平板天线在掌子面上横向移动，间隔0.1 m进行1次探测，并得到1道数据。当电磁波在掌子面前方岩体内传播至不良地质体界面时发生发射，反射信号被接收天线接收。本文数值算例中采用的源子波类型为雷克子波，运用Matlab编制的时域有限差分程序对宽10 m、长20 m的模型进行模拟，正演网格空间步长0.05 m，整个程序运算均在1台CPU配置为Intel Xeon E5-2603的小型服务器上进行。通过Matlab内置的parfor函数，在1次正演过程中同时进行12线程的运算，耗时约为81 s。单次迭代过程共包括正传波场计算、逆时残场计算和2次迭代步长计算等4次正演过程。全波形反演过程中对于内存需求主要集中在几千个时间步长的正传波场场量存储上。在程序计算中，模型网格数为200×400，单个网格存储类型为single单精度类型，大小为4字节，并且需要同时存储12个线程的场量，这就需要产生大约10 GB的内存。

图2 探地雷达隧道检测示意图

2.1 溶洞模型数值算例

溶洞成因于岩石的岩溶作用，其物性参数与周围岩体存在明显差异。溶洞内部常常会充填碎黏土、水、空气等物质，使溶洞与周围岩体之间的物性差异更加明显[19]。溶洞模型如图3所示。设置溶洞模型Ⅰ、模型Ⅱ、模型Ⅲ如图3(a)、3(d)、3(g)所示，背景介质为灰岩，在距离掌子面5，7，13 m处分别设置3个大小约为1 m×1 m，2 m×2 m，4 m×4 m的不规则溶洞，溶洞内部分别充填空气、潮湿黏土和水。在图3(b)、3(e)、3(h)所示的正演结果中，黑色虚线为溶洞位置。溶洞响应均呈现出连续的双曲线同相轴，弧顶处为溶洞的近掌子面边界，远端界面及溶洞形状均无法获得。模型Ⅱ和模型Ⅲ正演结果中的反射信号出现了明显的同相轴反转。由于电磁波信号在岩层中衰减，溶洞③响应并不强烈，反射信号勉强能够辨识出其近掌子面位置。在如图3(c)所示的模型Ⅰ全波形反演结果中，溶洞①和溶洞②的反演结果较为理想，溶洞形状基本被反演出来，结果中指示的溶洞位置分别位于掌子面前方5 m和7 m处，与实际位置一致，相对介电常数值为0～2，也与预设模型相符合。模型Ⅰ正演结果中双曲线弧面较长，仅能判断近掌子面界面，而无法分辨其横向宽度。溶洞③的反演结果并不理想，但是仍能辨识出溶洞③的前界面在掌子面前方13 m处及大致横向宽度为4 m。在如图3(f)和图3(i)所示的模型Ⅱ和模型Ⅲ的全波形反演结果中，溶洞①和溶洞②的形状与模型吻合较好，模型Ⅱ反演的溶洞①和溶洞②的介电常数值与预设模型接近，但是模型Ⅲ中溶洞①和溶洞②的反演结果在12～15，与预设模型中水的相对介电常数81差距较大。模型Ⅱ和模型Ⅲ的全波形反演结果中，均没有反映出溶洞③的介电常数范围,只能对其近掌子面位置进行判断,考虑原因为接收到的溶洞③反射信号微弱，同时介质电性参数差异过大,影响反演结果。

(a) 溶洞充填空气模型(模型Ⅰ)

(b) 模型Ⅰ正演结果

(c) 模型Ⅰ全波形反演所得相对介电常数分布

(d) 溶洞充填潮湿黏土模型(模型Ⅱ)

(e) 模型Ⅱ正演结果

(g) 溶洞充填水模型(模型Ⅲ)

(i) 模型Ⅲ全波形反演所得相对介电常数分布

2.2 岩体裂隙模型数值算例

岩体裂隙主要成因于岩体内部应力作用或地壳运动，其往往会造成岩体应力失稳[20]。岩体裂隙模型如图4所示。设置岩体裂隙模型Ⅳ、模型Ⅴ、模型Ⅵ如图4(a)、4(d)、4(g)所示，背景介质为灰岩，在距离掌子面7～12 m处设置2个倾角(以掌子面为基准)分别为30°、40°的裂隙，其长度分别为6 m和5 m，内部分别填充空气、潮湿黏土、水。在图4(b)、4(e)、4(h)所示的模型正演模拟结果中，黑色虚线指示裂隙长度走向。从电磁波反射信号中可以清楚地观察到2组反射非常强烈的同相轴，方向与设置模型中的裂隙走向呈交叉状，裂隙尖端呈现抛物线状绕射信号，图4(e)和图4(h)所示模型Ⅴ和模型Ⅵ的正演结果中反射信号同相轴反转。由于同相轴走向与真实裂隙形态差异较大，单纯基于正演剖面图分析无法获知前方岩体裂隙形态。图4(c)和图4(f)分别为模型Ⅳ和模型Ⅴ的全波形反演结果。预制模型中灰岩的相对介电常数为7，空气的相对介电常数为1，潮湿黏土的相对介电常数为15，反演结果中指示的背景围岩相对介电常数值为6～8，图4(c)和图4(f)中异常体相对介电常数分别为0～2和14～15，均与模型较为吻合，反演得到的裂隙形态走向、长度与模型相比也较为相符，基本实现了对模型的复现，能够通过全波形反演结果简单地辨识出掌子面前方岩体裂隙的实际形态。图4(i)所示模型Ⅵ的反演结果中，裂隙形态反映较好，介电常数值在18～20，与设置模型参数相差较大，但是与图3(i)所示的模型Ⅲ反演结果相比，模型Ⅵ反演得到的高介电常数分布更加理想，其可能原因为裂隙狭长形态使得在探测中能够反射更多的电磁波信息。

2.3 岩体破碎带模型数值算例

岩体破碎带是完整岩体在构造运动作用下产生的破碎区域，主要表现为区域岩体分块破碎，块体空隙内部可能充填空气或水等[21]。岩体破碎带模型如图5所示。设置岩体破碎带模型Ⅶ、模型Ⅷ、模型Ⅸ如图5(a)、5(d)、5(g)所示，背景介质为灰岩，在距离掌子面8～12 m处为破碎区域，破碎岩体空隙分别填充空气、潮湿黏土和水。图5(b)为模型Ⅶ正演模拟得到的雷达剖面图，可以看出异常体区域的电磁波反射信号同相轴交叉错乱，因此考虑异常体内部介质性质变化较大； 同时异常信号反射区域的前后边界面较为清晰，指示的异常区域也与设置模型基本吻合。然而，仅依靠正演结果分析无法对异常信号区域介质性质进行判断。图5(c)为模型Ⅶ的全波形反演结果，相较于正演结果，破碎带前界面走向在图5(c)得到了良好的重现； 异常体区域的相对介电常数主要分布在1～2和6～8，考虑到空气和灰岩的相对介电常数，可以推断不良地质区域存在大量破碎，破碎区域内部充填空气。反演得到的低介电常数位置，基本实现了对岩体内部小空隙的归位，基于全波形反演结果，可以更加直观地了解破碎区域内部的岩体碎块和空隙分布状态。图5(e)和图5(h)分别为模型Ⅷ和模型Ⅸ正演结果，可以看出电磁波反射信号主要集中在破碎带浅层区域，同相轴呈现明显反转。与图5(b)所示模型Ⅶ正演结果相比，电磁波在含水地层中衰减比在空气中更加严重，破碎带深部区域反射信号微弱，特别是模型Ⅷ设置的潮湿黏土中，深层电磁波信号几乎完全衰减。图5(f)和图5(i)分别为模型Ⅷ和模型Ⅸ的全波形反演结果，可以看出掌子面前方8 m处出现有明显的相对介电常数高值区域，其范围分别为12～13和9～10，且高值区域走向和预设模型吻合。模型Ⅸ的反演结果较差，反演结果与实际模型参数差值在70左右，但是考虑其较高的介电常数，仍可判断异常区域含水。由于电磁波信号在含水地层衰减，模型Ⅷ和模型Ⅸ中破碎带深层区域反演结果均不理想。

(a) 裂隙充填空气模型(模型Ⅳ)

(b) 模型Ⅳ正演结果

(c) 模型Ⅳ全波形反演所得相对介电常数分布

(d) 裂隙充填潮湿黏土模型(模型Ⅴ)

(e) 模型Ⅴ正演结果

(f) 模型Ⅴ全波形反演所得相对介电常数分布

(g) 裂隙充填水模型(模型Ⅵ)

(h) 模型Ⅵ正演结果

(i) 模型Ⅵ全波形反演所得相对介电常数分布

(a) 岩体破碎带充填空气模型(模型Ⅶ)

(b) 模型Ⅶ正演结果

(c) 模型Ⅶ全波形反演所得相对介电常数分布

(d) 岩体破碎带充填潮湿黏土模型(模型Ⅷ)

(e) 模型Ⅷ正演结果

(f) 模型Ⅷ全波形反演所得相对介电常数分布

(g) 岩体破碎带充填水模型(模型Ⅸ)

(h) 模型Ⅸ正演结果

(i) 模型Ⅸ全波形反演所得相对介电常数分布

2.4 不良地质组合模型数值算例

在实际工程中，地质状况更加复杂，多类不良地质体常常相伴出现，因此，本文在上述不良地质类型的基础上设计了组合地质模型。不良地质组合模型如图6所示。设置模型Ⅹ如图6(a)所示，背景介质为灰岩。在掌子面前方9 m处设置1个大小为2 m×4 m的椭圆形溶洞，溶洞两端连接裂隙，内部充填水。在模型Ⅹ的正演结果(图6(b))中，可以看到3组不同走向的同相轴，中间1组双曲线同相轴为溶洞响应，弧顶处为溶洞的近掌子面位置； 两侧同相轴为裂隙响应，其位置与裂隙的实际位置存在偏差。在模型Ⅹ的全波形反演结果中，裂隙起于掌子面前方8 m处，终止于掌子面前方12 m处，呈120°夹角，与实际模型对应良好； 反演结果中指示的内部介质相对介电常数值在12～15，与设置模型相差较大，但仍可依据反演结果中裂隙区域较高介电常数值推测地层含水； 溶洞前界面的反演结果较好，其位置位于掌子面前方9.5 m处，未能反演出溶洞的后界面。设置模型Ⅺ如图6(d)所示，背景介质为灰岩。在掌子面前方8～9 m处设置1条破碎带，破碎带后方距掌子面13 m处有1个大小约为3 m×3 m的圆形溶洞，内部充填空气。在模型Ⅺ的正演结果(图6(e))中，可以看到明显的杂乱反射区域，反射信号同相轴交错叠加； 在模型中后部能观察到一近似平行于掌子面的同相轴，但是反射信号响应十分微弱。在模型Ⅺ的全波形反演结果中，破碎位置与设置模型吻合，反演得到的介电常数值在2～3，与设置模型相近。溶洞反演结果较差，在图6(f)中能辨识到其近掌子面界面，且介电常数较低，在3～5。设置模型Ⅻ如图6(g)所示，背景介质为灰岩。在掌子面右侧前方8 m处存在4 m×4 m的岩体破碎区域，破碎区域一端连接裂隙，内部充填空气。模型Ⅻ的正演结果如图6(h)所示，能明显辨识出一倾斜同相轴，为裂隙位置。破碎区域电磁波反射信号较为杂乱，且多条同相轴与裂隙反射信号交错。通过正演结果仅能识别裂隙位置，无法分辨破碎区域。在模型Ⅻ的全波形反演结果中将裂隙的位置很好地呈现出来，指示的内部介质相对介电常数值在1～2，与设置模型一致。反演结果中岩体破碎位置和范围与设置模型吻合较好，对于破碎区域的相对介电常数反演结果，近掌子面结果与模型较为一致，深部区域的介电常数反演较差。

(a)溶洞-裂隙充填水模型(模型Ⅹ)

(b) 模型Ⅹ正演结果

(c) 模型Ⅹ全波形反演所得相对介电常数分布

(d) 岩体破碎-溶洞充填空气模型(模型Ⅺ)

(e) 模型Ⅺ正演结果

(f) 模型Ⅺ全波形反演所得相对介电常数分布

(g) 岩体破碎-裂隙充填空气模型(模型Ⅻ)

(h) 模型Ⅻ正演结果

(i) 模型Ⅻ全波形反演所得相对介电常数分布

2.5 全波形反演适应性分析

影响全波形反演的主要因素有初始模型、源子波形态、反演目标体形态等。本文所采用的初始模型为均匀模型，得到了良好的全波形反演结果。若基于先验地质信息建立初始模型，可以进一步对反演结果进行优化。在数值模拟中，源子波形态是已知的，但是在处理实际数据时，并不知道实际的发射波形，这就需要对源子波进行估计，并在每次迭代过程中，更新源子波信号。全波形反演对于狭长的、平行于掌子面的不良地质体有更好的反演效果，其原因为狭长的界面能够反射更多的电磁波信号。基于本文展示的裂隙模型反演结果，均清晰地显示出了多条裂隙的走向和形态。在溶洞模型中，目标体相对于雷达测线的反射面较小，能够获得的反射信息较少。对于靠近掌子面的小型溶洞，反演结果较好，位置及相对介电常数均较为清晰，而深层溶洞仅能反演得到前界面位置； 特别是在设置的组合模型Ⅺ中，破碎带后方的溶洞响应极其微弱，反演结果得到的前界面也非常微弱。在破碎带模型的反演结果中，破碎区域能够反射更多的电磁波信号，浅层破碎区域走向较为清晰。在实际工程中，由于并不知道不良地质体相对于掌子面的走向和形态，因此可以采用2条测线(1条水平、1条垂直)，对多个方向的地质雷达数据进行反演分析。当反演区域内部的相对介电常数变化过大时，全波形反演过程中对于区域介电常数的迭代更新往往不足，反演结果与实际模型存在较大差值，像充填水模型反演得到的介电常数值与实际差值在50以上，但是能观察到介电常数高值区域。此时，应结合实际地质资料进行综合评判。

3 结论与讨论

针对传统隧道探地雷达数据处理中仅能依靠经验对不良地质体进行判断的局限性，将全波形反演方法应用到探地雷达数据处理中，针对典型不良地质体建立数值模型，通过系统的正演模拟和全波形反演研究，取得了较好的效果。

1)电磁波在围岩中的传播特征较为复杂，雷达剖面图中的同相轴与实际异常体形态往往存在较大差异。在岩体裂隙数值算例中，电磁波反射信号同相轴与实际裂隙走向完全不同，直接依据雷达剖面图无法获得异常体的形态特征。

2)采用全波形反演方法可以反演出掌子面前方异常体的电性参数和位置形态，4个数值算例的全波形反演结果均有助于对预设模型的判断。例如： 含水地层在反演结果中表现出较高的相对介电常数值，通常在15以上； 空洞在反演结果中表现出较低的相对介电常数，在0～2。

3)当不良地质体介电常数差距悬殊时，全波形反演结果不理想； 同时，电磁波在传播过程中的衰减会影响全波形反演的结果。例如，溶洞算例中后界面反射信息较少、岩体破碎带算例中深部破碎反射信号微弱，都影响了相关区域的全波形反演结果，但是综合先验地质资料分析，仍然可以对前方地质状况进行大致判断。

4)由于实际探测地质雷达数据信噪比较低，全波形反演在隧道超前探测中的应用效果并不理想。下一步将重点开展实际数据滤波降噪算法，以提高实际数据全波形反演效果。