近年来，随着我国地下基础设施建设加快，基坑开挖加深，作业环境与施工条件也随之复杂，工程项目不可避免的会遇到不良地质构造，诸如地下溶洞、空洞、底层软弱不均等，对施工安全和施工进度造成严重的影响。为了使地下工程顺利施工、安全通过不良地质分布区域，提前对不良地质构造进行探测并根据探测结果进行数值分析，充分评估其对施工进程造成的不利影响尤为重要^[1]。

弹性波CT和地质雷达作为常见的两种物探手段，在工程地质检测中应用广泛。弹性波CT法具有探测范围大、探测精度高、抗干扰能力强等优势，在溶洞探测中有良好的应用效果^[2-3]。Hu等^[4]使用电磁波CT法对有岩溶存在的高速公路路基的注浆加固效果进行了检测，根据所采集电磁波吸收系数的差异，得出了溶洞、含卵石土层、粘土层等不同路基介质的注浆加固效果，对高速公路路基的填土、修筑及加固修补都有重要指导意义；刘少虹等^[5]通过电磁波CT技术传播与探测理论，依靠电磁波吸收系数与掘进过程中冲击地压间的关系，建立了电磁波CT技术对冲击地压危险等级的评价模型，并据此给出危险等级评价指标；徐智勇等^[6]利用井间电磁波CT技术，对昆明地铁4号线苏家塘站专项线路工程地下岩溶发育地段进行测试，通过对岩石电磁波吸收系数差异的研究，依据吸收系数特征，探测到地铁线路工程地下岩溶及裂缝发育，划分出4个层次的岩溶区，与钻孔岩心资料对比验证了结果的正确性。弹性波CT利用速度反演可得出二维地质速度剖面图，通过速度剖面进行地质填图解释。

为了方便非物探人员与不具有相关经验的学者，满足对工程区域内进行快速全面的分析与评价，近年来，许多学者对弹性波CT的三维地质建模进行了研究。蔡连初等^[7]将多个弹性波CT二维剖面图作平面内最短射线路径搜索计算波速，建立三维空间波速阵列数据并导入专用立体绘图软件绘制了等值面图，得到了清晰的空间立体波速三维分布图；刘超杨等^[8]采用了克里金插值的方法，获得了岩体空间三维波速的表达模型，通过空间波速模型进行地质解释，提高了解释精度；张琦等^[9]利用现场钻孔数据与二维地质剖面数据表，将数据导入建模软件建立三维空间溶洞模型，从空间上给出桩基地下溶洞的三维空间形态；王启明等^[10-11]利用弹性波CT数值模拟将网格进行了分区克里金插值，获得了波速三维结构模型，根据波速进行峰值法区域分割，形成了研究区域内的三维网格模型。最后将三维重构技术应用到实际工程中取得了良好效果，精确地反应了地下岩溶分布空间状况，为地质稳定性支护与评价提供了重要意义。

虽然弹性波CT技术在常规不良地质的探测中获得了较大的发展，但相对于医学CT技术，无论是采集效率还是异常体的边界分析精度都有很大差距^[11]。分水岭算法是一种图像分割技术，主要用于边界提取和图像分割，在医学、数学等学科中应用广泛。分水岭算法最初是由Digabel等^[12]引入到图像处理邻域，后来由Beucher等^[13]对其做进一步的研究，进而得到在灰度图像上进行分水岭分割算法的结论。分水岭算法具有计算速度快和对图像的边缘能够精准定位的优点，同时对图像中微弱边缘具有良好的响应。

综上所述，弹性波CT的二维岩体病害识别和三维地质建模发展迅速，但是对于异常体边界的分析精度有待提高。本文引用分水岭算法对弹性波CT的速度剖面梯度图进行分割，识别溶洞；基于波速数据库，实现复杂岩体综合地质建模，并辅以地质雷达、窥孔成像等手段综合解释和验证。研究结果为基坑支护设计与开挖过程中的安全性和稳定性提供重要指导与借鉴，也为地下工程设计及防治水提供依据。

1.   数据采集与方法介绍

本文首先使用弹性波CT和地质雷达两种手段进行现场数据采集，然后利用分水岭算法进行岩体病害识别，最后对目标区域进行地质建模。其过程分为3个关键步骤：二维岩体病害识别、三维地质建模、岩体病害钻孔摄像验证。总体流程如图1所示。

图  1  总体流程图

Figure  1.  Overall flow chart

下载: 全尺寸图片 幻灯片

1.1   工程概况

徐州市地铁3号线三环南路站位于南三环路和北京路交叉口，沿北京路南北向布置，是轨道交通3号和4号线的换乘站，车站地理位置见图2。本标段土层以填土、第四系粘土及粉土为主，下伏基岩为寒武系灰岩及奥陶系灰岩。地貌类型为岩溶丘陵，标高46～55 m，由第四系全新统棕红、棕黄色黏土构成。场地范围内存在岩溶问题，溶洞多由黏土或黏土夹溶蚀碎石半充填，仅少量溶洞为空洞。岩溶发育程度为弱发育，一般发育在基岩顶面以下30 m范围内，局部地段较深。

图  2  项目地理位置

Figure  2.  Geographical location

下载: 全尺寸图片 幻灯片

根据规范^[14]中不良地质溶洞的勘察要求，对该项目场地进行岩溶的地质雷达和弹性波CT探测试验，从而获得溶洞的分布、形态和构造。由基坑边坡揭露的结构面发育特点及钻孔条件，在勘探区域中布置6个钻孔，钻孔布置见图3。

图  3  钻孔布置图

Figure  3.  Drilling layout

下载: 全尺寸图片 幻灯片

1.2   弹性波CT现场测试

弹性波CT是利用弹性波在探测介质中的传播规律，根据弹性波的走时和激发点与接收点之间的距离来计算弹性波波速与能量衰减等其他物性参数，然后借助计算机重建技术，重现物体内部结构^[15]。弹性波CT根据采集数据方式的不同又分为井间地震CT与井地地震CT。

采用井地和井间地震相结合的方式进行探测。探测仪器使用SUMMIT地震探测仪，井中水听器使用岳阳奥成公司的12道接收水听器和自制24道接收水听器。井地地震地面激发采用200 kg落重震源，井间地震井中激发采用TD-Sparker 20 KJ电火花震源。数据采集系统见图4。采集时，ZK1、ZK2、ZK3、ZK4、ZK5、ZK6共6个钻孔勘探区各采用东西向和南北向两条测线，共计12条测线，每条线20～40个激发点不等（激发点距1 m）的观测方式进行数据采集。每个激发点激发2～4次，选取最佳单炮记录，井中观测点设计为深度间隔1 m，共24道井中水听器接收。井地地震是根据GPS测量点在场地纵横方向布置炮点，炮距1 m，24道水听器放置井中，检波点距1 m。井间地震采集中6个钻孔分别激发，其余孔接收，激发点距1 m，接收点距1 m，激发、接收井的范围都为0～23 m，共形成13个井间地震CT断面。采集到的原始资料见图5。

图  4  弹性波CT数据采集系统

Figure  4.  Elastic wave CT data acquisition system

下载: 全尺寸图片 幻灯片

图  5  弹性波波形数据

Figure  5.  Elastic wave shape data

下载: 全尺寸图片 幻灯片

经过现场弹性波CT试验后使用专业商业软件Vista、Tomodel、Surfer进行反演成像处理，具体流程及模块设置参数如下：

（1）观测系统建立：首先根据地震数据野外采集班报建立处理时所需要的观测系统，包括观测点距、接受首段桩号、尾端桩号等。

（2）初至拾取：初至时间拾取是否精确，关系到能否建立可靠的时—深关系、平均速度和层速度的计算精度，为后续的排齐、提取子波等处理提供可靠的参数。初至拾取可以在Vista软件中拾取，也可以在绿山Picker模块中拾取。

（3）成像处理：成像处理是地震CT资料处理的关键步骤，对准确地进行构造分析、研究储层以及采空区的空间变化具有重要的意义。

（a）模型空间离散化。模型空间离散化包括，第一步将三维模型空间分成一系列的网格（Voxel）。网格大小可以人为选择，一般而言，网格应该足够小才能分辨感兴趣的最小地质体；同时，网格也应该足够大，才能获得充分的射线，增加成像的可信度，这里选择Inline方向为0.5 m，Crossline方向0.5 m，深度方向0.5 m。

（b）初始速度模型定义。定义初始模型速度，在理想情况下，可以根据该区域已知的地质和地球物理信息来定义初始模型速度。如测井资料来建立，或利用梯度变化的速度模型。通常而言，利用合适的初始速度模型开始迭代能够提高最终成像的可信度。另外，初始速度模型应当保证地震能量能够从炮点传播到检波点。

（c）射线追踪和分割。射线追踪和分割包括计算每个炮检对的旅行时间和射线路径，以及根据射线路径获得通过每一个网格的射线段。射线段的长度在反演问题中用于计算加权和更新速度。在绿山的层析方法中，采用Um和Thurber于1987年提出的最大速度梯度射线追踪三维算法，这是一种两点射线追踪方法，基于计算效率以及射线追踪方法的优点，可以避免内插。

（d）速度反演。在层析中设置速度反演的批量处理参数，如射线追踪参数，包括最小的射线路径反射点、最大的扰动误差等；速度反演参数，包括速度平滑、最小速度和最大速度的边界值；收敛参数，包括定义不满足射线数量的百分比、速度的变化百分比、超出模型边界的射线条数百分比来设置停止运算标准。层析方法也设置了使用初至门槛值，拒绝质量差的初至参与计算，提高速度反演的精度。在批量参数设置中，用户可以选择下一个迭代开始的迭代次数，也可是输出中间速度结果，进行对比。

处理过程包括观测系统建立、模型空间离散化、初始速度模型定义、射线追踪和分割、速度反演最终得到弹性波CT速度成像切片如图6所示。

图  6  南北方向CT速度剖面图

Figure  6.  CT velocity profile in north-south direction

下载: 全尺寸图片 幻灯片

1.3   地质雷达现场测试

地质雷达是一种利用人工激发的电磁波在具有电性差异的介质中传播，通过研究和观测电磁场的分布和变化特征，实现介质内部异常体和不同介质层面的识别和定位的物探方法^[16]。结合现场实际情况，采集仪器选择Pulse EKKO雷达，测试系统见图7。

图  7  EKKO探地雷达测试系统

Figure  7.  EKKO GPR test system

下载: 全尺寸图片 幻灯片

采用高低频天线联合探测，综合解释，提高解释准确度。具体采集参数为：低频天线频率100 MHz，时窗400 ns，采样间隔0.8 ns，天线间距1 m，采样间距0.25 m；高频天线频率250 MHz，时窗150 ns，天线步距0.25 m，采样间距0.05 m。实际探测过程中，通过皮尺现场标定测线，确定钻孔与测线间的距离。安装好设备之后，通过等距离移动天线进行现场采集。经过初步滤波、降噪以及增益等数据处理后，得到反射剖面图如图8所示。

图  8  雷达探测剖面图

Figure  8.  Radar detection profile

下载: 全尺寸图片 幻灯片

由于电磁波在传播过程中会出现绕射现象，所以对于溶洞、结构面以及地层界面等不同地质结构的反射特征将会和二维反射剖面的同相轴形状相近。因此，从雷达探测剖面图上，可以直接判断和读取相应的地质构造情况。从图8雷达反射剖面图中，可以清楚地看到同相轴呈现出双曲线形状，据此可以初步判读此处为溶洞构造。

2.   岩体病害的判断与识别

按照前期地质资料显示，场地范围内存在岩溶问题，溶洞多由黏土或黏土夹溶蚀碎石半充填，仅少量溶洞为空洞，岩溶发育程度为弱发育。下面从弹性波CT、地质雷达手段进行病害识别并用分水岭算法进行精准识别，同时获取岩体病害的三维空间坐标。

2.1   岩体病害的弹性波CT和地质雷达初步判断

深基坑岩体内部存在溶洞，结合地质资料可知溶洞为岩体内部受到地下水、化学等侵蚀形成的充填或半充填洞穴现象，其溶洞周边岩体完整性较好。从弹性波CT速度剖面图上的判断依据主要是剖面上形成的低速圈闭现象，而低速圈周边速度较高，形态上表现为不规则圆圈状，与溶洞形态相对应。地质雷达对于溶洞的判断，主要的依据为反射波同相轴的形态特征，即同相轴形态上表现为双曲线特征。对于深基坑岩质边坡中存在的断层裂隙带而言，其形成原因主要为岩体或岩层顺着破裂面发生位移的构造，其在弹性波CT速度剖面图上表现为条带状形态，而在地质雷达剖面图上表现为同相轴错断特征。据此，对断层裂隙带进行识别判断。

如图9中所示，钻孔1-2速度剖面图、钻孔2-3速度剖面图中红圈部位为低速的圈闭现象，形态上为不规则圆形且周围部位为高速区域，据此推断为R1、R2、R3、R4、R5溶洞所在部位。同样的，从地质雷达在此处的测线中可以看到，波形图上相应部位同相轴呈现出双曲线特征，推断为溶洞反射。图9（a）地质雷达图中R1与R2两处溶洞的判断与速度剖面图判断相吻合；而图10（b）地质雷达中R3与R4两处溶洞的判断与图9（b）速度剖面图中R3溶洞位置相对应。因此，通过地质雷达与弹性波CT手段可识别出相同位置溶洞4处，分别为R1、R2、R3、R4。弹性波CT速度剖面图上R5处溶洞无法在地质雷达波形图上识别对应，主要原因地质雷达在现场探测过程中信号受到干扰，对于深部岩体溶洞识别电磁波能量衰减严重以及雷达分辨率不足等原因导致无法识别。

图  9  弹性波CT病害识别图

Figure  9.  Elastic wave CT disease identification diagram

下载: 全尺寸图片 幻灯片

图  10  地质雷达病害识别图

Figure  10.  Geological radar disease identification

下载: 全尺寸图片 幻灯片

下面利用地质雷达与弹性波CT手段，对断层裂隙带进行判断。如图9（a）所示，在钻孔1-2速度剖面图中存在中间速度低两边速度高的带状低速区域，初步推断此处可能存在裂隙破裂带，其走向为东西方向。在图10（a）地质雷达剖面图上，可以看到F1处同相轴错断现象，其深度位置与图9（a）中F1处断层位置一致。因此，此处可判断为断层现象。

2.2   基于分水岭算法的岩体病害的弹性波CT识别

地质雷达和弹性波CT可用于岩体隐患识别，为准确圈定异常体范围，本文尝试采用分水岭算法实现之。分水岭算法是基于形态学中对图像进行梯度分割的算法，它是将空间位置与灰度值相近的像素点连接起来构成封闭轮廓线。图像的灰度空间就如地表结构，不同像素的灰度值代表不同的高度，其中灰度值较大的像素连成的线可以视为山脊，也就是分水岭。如图11所示为深基坑岩体从边坡到基坑处南北方向上的速度切片图，切片位置依次为X=2.4 m、5.9 m、9.4 m、12.9 m。

图  11  南北方向速度剖面切片图

Figure  11.  North-South velocity section

下载: 全尺寸图片 幻灯片

依据弹性波CT速度剖面图，利用MATLAB软件中的Sobel算子对速度图像进行梯度计算，图像的梯度相当于图像上相邻两点像素的差值。当图像存在边缘时，就会存在较大梯度值；当图像较为平滑时，其灰度值变化较小，则相应的梯度变化较小。通过梯度计算，观察速度剖面梯度幅值可以知道，速度变化较大的图像部位，其梯度值差异较大；速度变化较小的图像部位，其梯度值差异较小。

对于溶洞与断层裂隙带而言，由于其与周围岩体的速度差异较大，因此其在速度剖面梯度幅值上表现为梯度存在差异，于是据此可以对岩体不良地质体进行识别。但是直接利用基于梯度进行分水岭计算，会对图像造成过度分割现象，难以给出有用信息。本文采用梯度加掩模的三次分水岭算法对图像进行分割，其效果大大改善，处理结果见图12所示。

图  12  速度剖面梯度图与分水岭图

Figure  12.  Velocity profile gradient and watershed

下载: 全尺寸图片 幻灯片

如图12（a）所示，从梯度幅值图像上可以看出，速度剖面图上梯度存在较大差异，边缘线清晰且梯度图出现了较为清晰的“梯度圈闭”现象。依据梯度物理意义，梯度差异大，则速度变化率大，且形态上与溶洞相似，据此可推断R1、R2、R3、R4、R5五处位置为溶洞现象。通过分水岭计算可以看到，图像上出现了封闭的轮廓线，即分水岭线。将分水岭线与梯度图像进行对比，发现R1、R2、R3、R4、R5五处位置溶洞位置与梯度图像相对应，据此初步推断共有5处溶洞现象。浅层位置封闭轮廓线，推断其为表层杂填土；深部边缘处轮廓线，由于其所处位置为图像边缘，形成的封闭线无法判断其为溶洞边界线，且其位置处于高速岩体位置，与低速异常现象不符合，因此此处不是溶洞边界线。其余南北方向速度剖面梯度图像与分水岭图像如图12（d）和图12（f）所示。其中部分水岭图像中轮廓线较杂乱，难以做出有效判断，推断是由于该部位岩体较破碎所导致的。

下面对断层裂隙带进行判断识别。如上图12（a）所示，从梯度幅值图像上可以看到一条梯度差异较大的带状形态F1，带状形态两边梯度较大，则其速度较大，带状形态处速度较小，符合断层裂隙带的速度分布与形态特征。另外，通过分水岭计算，图像上同一位置同样出现了F1带状形态。据此推断此处为断层裂隙带现象。同理，在图12（c）～图12（f）中梯度图像与分水岭图像上同样出现了带状形态F2和F3。根据前文所述，判断F2与F3为断层。另外，从空间坐标关系中可以看出，F1、F2、F3在深度与Y轴方向上的位置一致，X轴方向上由X=2.4 m到X=9.4 m，跨度7 m。

综上所述，通过梯度图像与分水岭图像判断岩体病害，识别出溶洞5处，分别为R1、R2、R3、R4、R5；断层裂隙带1处，为F1位置。

2.3   测区岩体病害成果统计

通过地质雷达、弹性波CT以及图像处理技术，共识别出5处溶洞与1处断层。依据上述对溶洞的识别，下面利用弹性波CT速度剖面图、地质雷达反射波分析以及分水岭图像对溶洞的位置坐标、深度、走向进行分析判断。空间坐标X与Y表示溶洞在相应两个方向上的跨度，深度表示溶洞和断层在竖直方向上下界面的位置。空间坐标X与Y、深度由分水岭图像上直接读出，岩体病害空间描述如下表1所示。

表  1  岩体病害三维空间位置

Table  1.  3D spatial location of rock mass diseases

序号	X/m	Y/m	深度/m
R1	7.6～10	13～15	6.8～8.5
R2	7.8～10	6～8	10.8～11.2
R3	8～10	21～23	7～9
R4	9～10	18～22	14.3～15.5
R5	8～10	23～26	22～24
F1(A) F1(B)	2.4	11	10
	2.4	9	15

下载: 导出CSV 

| 显示表格

3.   岩体病害的三维地质建模

通过分水岭算法手段对异常体进行识别后，给出了溶洞与断层裂隙面的空间信息。建模的数据来源主要是以弹性波CT的三维空间波速信息为基础，结合分水岭算法得到的空间坐标信息获得建模数据库，将上述数据导入GOCAD软件，经过插值、平滑以及中值滤波等手段形成三维模型，如图13所示为溶洞的三维地质解释图。通过不同方位的展示，给出溶洞与断层的空间分布。其中R表示溶洞，F表示断层。

图  13  溶洞三维地质解释图

Figure  13.  3D geological interpretation of karst cave

下载: 全尺寸图片 幻灯片

从图13溶洞三维地质解释图可知该深基坑岩体中溶洞共有5处，溶洞分布区域较广，从浅部7 m到深部24 m均存在，纵向长度多集中于2 m左右。结合现场地质勘探资料可知，目标探测区域内揭露的溶洞大多半充填，多以硬～可塑黏性土为主，夹灰岩碎块，结构致密性差，透水性较强，有少量溶洞无充填，需加强支护和防水措施，或者采取必要的加固防护措施；对于跨度过大的溶洞，可采用钢筋混凝土梁跨越等方式；为防止洞穴坍塌和加强洞穴顶板岩层的稳定，可采用注浆加固等措施并加强衬砌。

同时，现场有断层裂隙面共有1处，分布在深度Y=9～16 m之间，研究区的岩体内部结构面较为发育，岩质边坡存在较大滑移风险。另外，通过图13可以看到，岩溶和裂隙面（潜在滑移面）存在贯通现象，进一步会增加深基坑边坡滑移的风险。

4.   岩体病害钻孔摄像验证

弹性波CT现场探测过程中人为操作与存在噪声等因素，会对信号造成干扰，进而影响反演获得的数据真实性。另一方面，地质雷达现场探测过程电磁波受环境因素干扰较强，造成数据的失真。因此，对于三维地质模型建立的准确性，需要依靠更加直观准确的方法进行验证。钻孔摄像手段是对现场存在的不良地质体的直接体现，具有清晰直观特点。下面依据不同钻孔不同深度处的摄像图，对三维地质模型的建立进行验证。

如图14所示，通过钻孔摄像给出了溶洞的直观图像与位置。根据弹性波CT现场探测可知，1号钻孔坐标为（6,10）。另外1号钻孔在深度10.8～11.2 m处存在溶洞，此处钻孔摄像下的溶洞位置与图13溶洞三维地质解释图中R1处溶洞相对应。

图  14  溶洞钻孔摄像图

Figure  14.  Video of karst cave drilling

下载: 全尺寸图片 幻灯片

在2号钻孔摄像图中，在深度15.3～16 m处存在溶洞，2号钻孔的坐标为（6,15）。此处钻孔摄像下的溶洞位置与图13溶洞三维地质解释图中R2处溶洞相对应。

在3号钻孔摄像图中，在深度8.7～9.2 m处、13.4～13.7 m处存在两处溶洞，3号钻孔的坐标为（6,24）。此处钻孔摄像下的溶洞位置分别与图13溶洞三维地质解释图中R4和R3处溶洞相对应。

由于溶洞三维地质解释图中R5溶洞所处位置未有钻孔穿过，因此无法利用钻孔摄像对R5溶洞进行验证。

总之，在有钻孔穿过的4处溶洞R1、R2、R3、R4，均在钻孔摄像中有明确的发现；R5溶洞因无钻孔穿过未能获得直接证据，但雨季施工时R5溶洞处有基坑坑壁出水现象，可为该处溶洞的存在提供佐证。

5.   结论

本文以地铁深基坑岩体边坡为研究对象，通过弹性波CT的数据处理，结合分水岭算法对探测区域的不良地质进行识别。同时依据弹性波CT和地质雷达手段对不良地质体的空间位置进行定量计算，减少了单一建模带来的误差性，综合识别不良地质体并进行相互解释与验证，给出整体岩体边坡三维地质模型。通过本文工作，得到了以下认识和结论。

（1）多种探测手段相结合，对于探测岩体病害有重要价值。通过弹性波CT与地质雷达的现场实测，获得目标区域的波速数据信息。分别对弹性波CT速度剖面图、地质雷达图进行初步病害识别。由于弹性波CT在准确圈定异常体边界上存在困难且地质雷达对深部岩体病害识别较差，因此引用分水岭算法对弹性波速度图进行精准识别。实践表明，采用分水岭算法可将弹性波速度梯度相近的像素点所构成的封闭轮廓线作为溶洞的边界，对于断层裂隙带等面状异常体则表现为梯度差异较大的带状形态，整体提高了识别精度。在研究区域识别溶洞5处和断层1处，同时获取了岩体病害的三维空间坐标。

（2）建立的基坑边坡岩体病害模型具备高精度、高准确度的特点，能够满足工程岩体稳定性评价的需要。通过弹性波CT手段对异常体进行识别后，给出了溶洞与断层裂隙面的空间信息，结合分水岭算法得到的空间坐标信息获得建模数据库，将数据库导入GOCAD软件，经过数据处理后形成以波速为属性的深基坑岩体模型。

（3）使用钻孔摄像对三维地质解释模型进行验证，判断二维岩体病害识别和三维地质模型的准确性。结果表明，所建立的岩体病害模型，与地质钻探、窥孔成像、现场观察等高度吻合，验证情况与勘察报告相符合。

本文从三维可视化角度给出了岩体病害的定量判断，为深基坑的稳定支护提供帮助，提出的岩体病害识别方法和建模技术可为岩体支护设计提供重要参考，同时也为城市地下工程水害防治提供地质依据。

致谢：中国矿业大学潘冬明教授、李元海教授等对本文的弹性波CT的数据采集处理和窥孔成像等工作提供了帮助。