0.   引言

溶洞的存在不仅制约了隧道施工进度，还对人员财产安全构成严重威胁，勘察溶洞规模形态是制定相应处理方案的基础。因此，在隧道空间内准确探测溶洞的空间位置和形态是目前隧道建设工程中亟需解决的重要问题^[1-3]。

探地雷达因具有分辨率高、无损检测和快速成像等优点，目前已经成为隧道建设工程中浅埋岩溶不良地质主要探测手段^[4-6]。为提高探地雷达隧道岩溶不良地质的识别能力，含有不同物理意义的属性分析方法被逐渐引入到探地雷达数据解释中，如：刘东坤等^[7]研究了典型岩溶不良地质探地雷达反射信号的频谱特征；刘宗辉等^[8]研究了不同类型岩溶不良地质的属性特征，并提出了一种岩溶不良地质体反射信号探地雷达属性特征定量表示方法；Reis等^[9]应用不同探地雷达属性对地下岩溶塌陷进行了成像，研究发现部分探地雷达属性能够增强地层内部几何形状的成像效果。虽然目前在隧道探测方面应用探地雷达属性的研究成果较多，但大部分偏重于对二维图像和单道信号进行特征研究，而对隧道溶洞目标的三维几何形态研究相对较少。

事实上，已有许多学者对地下目标体的三维可视化展开了研究，并在不同探地雷达领域得到成功应用。李尧等^[10]通过改进BP算法提出了一种岩体裂隙的探地雷达三维成像方法。Zhu等^[11]利用三维探地雷达检测了树木的根系大小和分布，并通过自行设置振幅阈值实现了地下根系的三维可视化。邓海明等^[12]研究了基于半空间扫描测量模式的隧道坍腔探地雷达三维成像技术。胡群芳等^[13]应用三维探地雷达对城市市政管线渗漏情况进行了探测，并取得了较好的效果。上述研究结果表明探地雷达三维探测能够提高地下目标几何形态的识别能力，但大部分传统三维可视化方法主要是通过设置振幅阈值来显示目标反射信号，存在主观性强、理论性差等问题^[14]。

探地雷达三维属性是在二维属性的基础上发展起来的，主要通过二维数据插值或直接对三维数据进行计算获得。在本次研究中，提出了一种基于探地雷达三维属性分析的隧道溶洞三维可视化方法，并应用物理模型试验和现场案例数据对该方法的可靠性和适应性进行验证。

1.   方法原理

1.1   隧道内探地雷达三维探测方法

探地雷达采集系统通过发射天线向地下发射高频脉冲电磁波，当电磁波遇到不连续的地质体时会发生反射，并被接收天线记录。通过在测线方向上移动天线，探地雷达采集系统会根据采样步长自动记录每道反射信号，从而得到该测线上的探地雷达剖面。当需要获取地下目标体准确的空间几何形态特征时，进行三维探测则是非常有必要的。常用的三维探测方法有真三维和伪三维(2.5维)两种，真三维是指不同方向上步长增量相同，且小于Nyquist采样步长^[15]。伪三维是指采用多个平行测线或规则网格测线进行探测，并利用探地雷达文件中的坐标信息将不同测线上的雷达数据合成三维数据体。由于岩溶隧道探测环境复杂，探地雷达的三维探测通常采用伪三维探测方法，本次隧道内的探地雷达三维探测方法如图 1所示。

图  1  岩溶隧道探地雷达三维探测示意图

Figure  1.  Schematic diagram of GPR 3D detection in karst tunnel

下载: 全尺寸图片 幻灯片

传统探地雷达三维可视化方法是针对合成后的三维数据体设置振幅阈值，仅显示振幅较高的反射信号，从而提取出地下目标体的几何结构和空间分布特征。然而这种三维可视化方法的振幅阈值设置过度依赖于解译人员的经验判断，设置的振幅阈值不同，得到的结果也不相同^{[14, 16]}。

1.2   常用三维探地雷达属性

目前探地雷达属性分析方法较多，不同属性对不同地下目标反射信号的表征效果也存在差异^[17-18]。由于岩溶不良地质体与围岩之间的介电常数差异较大，电磁波反射能量较强，使得有效反射数据的振幅强度较高。因此与振幅相关的属性在实际工程中较为常用，如最大谱振幅、瞬时振幅和能量等常见探地雷达属性。

(1) 最大谱振幅

最大谱振幅是一种可以描述目标信号幅值和频域特征的常见频谱类属性。由于电磁波在不同介质中传播时频谱和反射信号幅值均会发生改变，因此利用最大谱振幅分析地下岩性变化处的探地雷达数据时具有较好的应用效果^[9]。

最大谱振幅属性的计算方法一般是通过短时傅里叶变换或小波变换把指定时窗内实信号转换到频域进行分析，然后获取频域内的最大幅值^[19]。计算时窗的大小一般接近天线主频的周期。图 2显示了不同频谱类属性的关系。

图  2  频谱类属性示意图

Figure  2.  Schematic diagram of frequency attributes

下载: 全尺寸图片 幻灯片

(2) 瞬时振幅

瞬时振幅是一种描述信号振幅瞬时特征的分析方法。在探地雷达应用中，通常采用希尔伯特变换或小波变换将分析信号构造复信号后求取^[20]。瞬时振幅的幅值始终为正，该属性对信号振幅较为敏感，但对相位敏感程度较低。

对于实信号x(t)的希尔伯特变换可表示为

瞬时振幅为

(3) 能量属性

能量属性是探地雷达信号分析中较为常用的振幅类属性分析方法，能够显著增强有效反射信号的强度^[19]。该属性计算方法是指定时窗内样本的振幅平方和除以样本数，可表示为

式中，${x_i}$为样点的幅值，n为计算的窗口大小。该属性对于低振幅反射具有较好增强效果，在应用该属性绘制三维图像时，具有主观性低和稳定性好等优点。

1.3   K-Means聚类原理

K-Means聚类算法是一种无监督的聚类算法，在信号处理中，该算法主要是对数据点进行有规律的分类。为提高分类精度，需要选取合适的相似度评估函数来进行分类，以期将各数据点准确分配到相应的数据集中。误差平方和准则函数因其理论简单、计算速度快和分类效果好，已经成为目前K-Means聚类算法中最常用的数据点相似度评估函数之一^[21]。误差平方和准则函数可表示如下：

式中：r_i为每个类别的中心点；b为第i类数据集X_i中的样点数据。

对于输入的三维探地雷达数据，K-Means聚类算法主要步骤如下：①格式转换。将三维探地雷达数据转换为一维格式，降低各数据点之间的横向联系；②设置初始聚类中心。设置背景数据和有效反射数据两种类别，并随机生成两个聚类中心点；③样点分类。计算需要分类的样本与中心点的最小欧几里得距离，并分配至最近的数据集中；④修正聚类中心点。计算各类数据集的均值作为新的聚类中心；⑤迭代聚类。重复步骤③和④，直到达到初始设定的最小误差平方和限定条件，或迭代次数；⑥聚类结束。

1.4   隧道内溶洞三维可视化方法

针对目前大部分探地雷达三维可视化方法振幅阈值设置存在主观性强的问题，本文提出了一种基于属性分析和K-Means聚类分析的隧道溶洞探地雷达三维可视化方法。该方法首先通过计算合适的属性来增强有效反射数据的对比度，然后通过聚类方法提取有效反射数据，进而利用数理统计方法获取有效反射数据的振幅阈值参数，最后基于该参数进行三维可视化和溶洞上反射面提取。

对于经过常规处理后的多个平行测线探地雷达数据，本文三维可视化方法处理步骤如下：①偏移成像处理。利用偏移成像方法处理每个探地雷达数据，提高溶洞目标成像精度。因为探地雷达收发天线具有一定的波瓣和频带宽度，使得地下目标在雷达图像上的反射信号与真实位置存在偏差^[22]。②合成三维数据体。利用探地雷达文件中的坐标信息将偏移成像处理后的数据插值为三维数据体。③计算三维属性。通过选取合适的探地雷达属性分析方法，增强探地雷达数据体中有效反射数据的对比度。通过室内参数调试和物理模型试验测试，本文主要使用最大谱振幅属性作为输入属性。④归一化处理。为方便分析，对获取的探地雷达属性体进行振幅归一化处理。⑤聚类分析。利用K-Means聚类算法将该数据分为背景数据和有效反射数据两种类别。⑥计算振幅阈值参数。计算标记为有效反射数据的振幅平均值，该值即为归一化处理后探地雷达属性体的可视化振幅阈值。⑦三维可视化。结合振幅阈值和归一化处理后的属性体，利用相关软件(如Matlab软件中的体可视化功能)实现溶洞目标的三维可视化。⑧溶洞上反射数据提取。将探地雷达属性体中小于振幅阈值的数据置为0，通过编程获取溶洞目标区域深度方向上的第一个不为0的点，去除奇异点(去除一定体积内分布较少的点)之后即可得到溶洞上反射面的点云数据。⑨绘制溶洞顶板面。通过编程或相关软件对提取的点云数据进行三维面绘制。

该方法的具体实现流程如图 3所示。

图  3  探地雷达三维可视化流程图

Figure  3.  Work flow of method of GPR 3D visualization

下载: 全尺寸图片 幻灯片

值得注意的是由于电磁波的衍射会产生类似双曲线反射的干扰信号，可能会与溶洞顶板反射信号形成叠加干扰，不仅给地下目标成像造成干扰，也不利于溶洞顶板反射面提取。因此，偏移成像处理在本文中是非常重要的步骤。偏移成像方法不仅可以用来聚焦双曲线反射，降低衍射干扰对溶洞成像和顶板面提取的影响，还可以帮助评估波速，有利于后续溶洞顶板岩层厚度估算。

2.   物理模型试验

为确定最佳输入属性并验证本文方法的可靠性，进行了室内物理模型试验。试验材料选用完整的灰岩石块，构建的物理模型如图 4所示。模型大小为0.9 m×0.5 m×0.7 m，共设计两个溶洞模型，直径均为0.04 m，深度分别为0.25，0.50 m。测试仪器采用意大利IDS公司生产的K2型探地雷达，选用天线的中心频率为1600 MHz。采样时窗为20 ns，采样点数为1024，道间距0.4 cm。在石块侧面共布置了14条水平测线，测线长度0.7 m，间距0.03 m。本文探地雷达数据预处理流程包括：数据编辑、废道剔除、零值偏移、直达波去除、增益、带通滤波、背景去除、噪声压制(如中值滤波)和合成三维数据体等。

图  4  物理模型试验几何模型和现场照片

Figure  4.  Geometric model and field image of physical model

下载: 全尺寸图片 幻灯片

图 5(a)为物理模型试验第7条测线得到的探地雷达图像。从探地雷达图像中可以发现两个溶洞模型均存在明显的双曲线反射信号。利用Reflexw软件中的速度分析方法进行速度标定后，可以看出电磁波在该灰岩试块中的传播速度约为0.11 m/ns。为进一步提高地下目标的成像精度，利用F-K偏移成像方法^[23]处理了图 5(a)中的探地雷达图像，结果如图 5(b)所示。图中可以看出经过偏移成像处理后，溶洞模型反射信号主要集中在时间深度3 ns和8 ns附近。将每个偏移成像处理后的探地雷达数据根据坐标信息插值为三维数据体(图 5(c))，然后采用等间距时间切片进行显示，图中同样可见在3，8 ns处的反射信号最明显(图 5(d))。

图  5  探地雷达二维剖面和三维数据体

Figure  5.  GPR images and 3D data

下载: 全尺寸图片 幻灯片

为对比常见探地雷达属性的成像效果，同时确定本文可视化方法的最佳输入属性，在计算探地雷达三维数据体的不同探地雷达属性(振幅归一化处理)后，提取不同属性体在8 ns处的切片进行对比分析，结果如图 6所示。从图 6中可以看出原始数据虽然可以大致显示溶洞模型的形状，但视觉效果较差(图 6(a))。而瞬时振幅、能量及最大谱振幅属性的视觉效果均较好，溶洞模型的几何形态特征识别度较高，其中最大谱振幅图像最清晰(图 6(b)~(d))。

图  6  不同探地雷达数据体在8 ns处切片对比

Figure  6.  Comparison of depth slices of different 3D volumes at 8 ns

下载: 全尺寸图片 幻灯片

为进一步验证本文基于K-Means聚类提取可视化振幅阈值方法的可靠性，通过聚类分析分别获取了原始数据体及其能量、瞬时振幅和最大谱振幅属性体的振幅阈值。通过聚类分析得到的各数据体振幅阈值参数分别为：原始数据24.6%，能量属性28.2%，瞬时振幅属性23.5%，最大谱振幅属性37.5%。基于振幅阈值参数分别提取了相应的等值面，结果如图 7所示。图中可以直观地看出最大谱振幅属性体的可视化效果最好，能够很好反映出溶洞模型的几何形态。三维可视化结果同时验证了最大振幅谱属性能够提高背景数据和有效反射数据的区分度。

图  7  原始数据和不同属性体三维可视化结果对比

Figure  7.  Comparison of 3D visualization results between original data and different attribute volumes

下载: 全尺寸图片 幻灯片

图 8为利用图 7(d)中最大谱振幅属性体提取得到的溶洞模型上反射面。从图 8中可以看出提取的反射面位置与实际溶洞模型基本一致，证明了本文溶洞上反射界面提取方法的有效性。

图  8  溶洞模型上反射面提取结果

Figure  8.  Extracted results of upper reflection interface of karst cave model

下载: 全尺寸图片 幻灯片

3.   现场案例

为进一步验证本文方法的适应性，利用在建高速公路隧道探地雷达实测数据来进行测试。所依托的高速公路隧道位于中国西南部广西壮族自治区红水河中段。同泽隧道所穿越山体主要为二叠纪下统茅口组灰岩，隧道进出口主要由微风化灰岩组成，表层覆盖层较薄，基岩大多裸露，岩质较硬，结构面发育。洞身段主要为微风化灰岩，岩体岩质较硬、结构完整，局部区域节理发育、破碎。隧道走向约为258°，在采用爆破开挖方式开挖至YK412+465附近，隧底左侧揭露出隐伏溶洞，揭露洞口长约7 m，宽约4 m，洞底最大深度约15 m(图 9)。经现场踏勘发现该溶洞存在两条支线，其中一条支线(走向约为262°)沿隧道走向延伸，另一条支线沿隧道边墙两侧延伸，与隧道走向近似垂直，受地形阻碍未能进行更详细的勘测。为保障隧道施工安全和后续治理方案提供依据，探明该区域隧道下伏溶洞的顶板厚度是非常有必要的。

图  9  隧道定位和现场照片

Figure  9.  Site case map and tunnel images

下载: 全尺寸图片 幻灯片

本文采用意大利IDS公司生产的RIS-K2型探地雷达，选用100 MHz收发分置的屏蔽天线。沿隧道轴线方向共布置了5条间隔2 m的纵向测线，测线长度为18 m，现场测线布置如图 10所示。探地雷达采样时窗为600 ns，采样间隔为0.05 m，采样点数为1024。探地雷达数据预处理流程与小模型试验数据一致，其中增益方法采用球面波指数补偿增益，带通滤波范围为40~240 MHz，使用中值滤波方法来减少奇异点噪声的影响。

图  10  现场网格测线布置示意图

Figure  10.  Schematic diagram of layout of site grid survey lines

下载: 全尺寸图片 幻灯片

图 11为测区内第1条测线的探地雷达图像及其偏移成像处理结果。图中存在多个反射信号，溶洞顶板反射信号区分度较低，且局部存在同相轴叠加干扰(图 11(a))。此外，通过Reflexw软件中的速度分析功能发现该区域电磁波传播速度约为0.11 m/ns。经过偏移成像处理后，浅层的异常反射信号得到很好还原(红色箭头)，而且还可以清楚地看出该测线处的溶洞顶板反射信号主要集中在12 m附近(红色虚线)。偏移成像后溶洞顶板反射同相轴叠加干扰区域得到有效处理，偏移成像结果有利于后续溶洞顶板面提取(图 11(b))。

图  11  第1条测线探地雷达剖面及其偏移成像

Figure  11.  GPR image of line 1 and its migration imaging

下载: 全尺寸图片 幻灯片

图 12为探测区域探地雷达三维数据体(经过偏移成像处理)的深度切片结果。从图 12中可以大致看出下伏溶洞顶板反射界面相对平稳，溶洞顶板主要分布在地下10 m左右。相对于三维数据体，二维剖面分析仅能判断单个测线方向上异常体存在与否，而三维数据体能够更详细的展示测区内地下不同深度的反射信号情况。然而，二维剖面和三维切片分析均无法直观地看出溶洞的几何结构。

图  12  探地雷达深度切片

Figure  12.  Deep slices of GPR data

下载: 全尺寸图片 幻灯片

通过计算隧道溶洞探地雷达三维数据体的最大谱振幅属性，并利用本文方法实现了下伏溶洞探地雷达三维可视化，结果如图 13所示。从图 13中可以看出溶洞顶板反射信号(青色)主要分布在8~11 m，相对较为平缓。此外，从图 13中还可以看出该区域存在一些独立的异常体(红色)，分布在深度4.2，9.2 m附近，三维可视化振幅阈值为最大振幅的62.1%。

图  13  基于最大谱振幅属性的隧道下伏溶洞三维可视化

Figure  13.  3D visualization of the hidden karst cave based on maximum spectral amplitude attribute

下载: 全尺寸图片 幻灯片

为进一步确定隧道下伏溶洞的顶板分布，提取了图 13属性体中溶洞的顶板反射面，结果如图 14所示。

图  14  隧道下伏溶洞顶板面提取结果

Figure  14.  Top panel of hidden karst cave

下载: 全尺寸图片 幻灯片

从图 14中可以看出，所提取的反射面能够较为直观地反映该区域溶洞顶板的分布情况。现场案例分析结果表明本文方法不仅能够增强隧道溶洞三维可视化效果，还可以帮助区分独立的异常体和主要的地下层反射目标。

4.   结论

(1) 本文针对探地雷达隧道内溶洞目标检测中存在空间形状确定困难的问题，提出了一种基于探地雷达属性分析的隧道溶洞三维可视化方法，该方法利用探地雷达属性来增强背景数据和目标反射数据的对比度，然后基于K-Means聚类分析方法提取三维可视化的振幅阈值参数，并采用等值面提取技术实现隧道溶洞三维可视化，最后基于可视化结果提取溶洞顶板反射面。

(2) 物理模型试验和现场案例应用结果表明：三维可视化方法能够更细致地反映地下目标的大小、形态和位置信息，提高了基于探地雷达视图地质解释的直观性；能量、瞬时振幅和最大谱振幅三类属性均具有较好的视觉效果，但其中最大谱振幅属性对探地雷达数据中的背景数据和有效反射数据区分度更高，最大谱振幅是本文可视化方法的较优分析属性；所提出的基于K-Means聚类分析振幅阈值提取方法有效解决了传统可视化方法在选择振幅阈值时存在过度依赖解译人员经验判断的问题，为隧道溶洞三维探测提供了解释依据。

(3) 由于溶洞内部多次反射信号的干扰，导致本文三维可视化方法对隧道底板下方隐伏溶洞上界面三维成像精度优于下界面。此外，受隧道现场探测环境制约，布置的测线相对稀疏，且未能利用其他物探技术对隧道溶洞进行详尽检测。在类似工程研究中应适量密集探测，并综合多种物探方法进行相互验证，提高隧道围岩内部溶洞检测精度。