跨孔地震CT技术在地铁工程施工溶洞探测方面的应用

2013-04-15段成龙阎长虹许宝田吴焕然邹明阳

地质论评 2013年6期

段成龙，阎长虹，许宝田，吴焕然，邹明阳

南京大学地球科学与工程学院，南京，210046

内容提要：本文介绍了在南京地铁4号线某区间应用跨孔地震CT技术进行现场数据采集的方法，并运用Matlab语言编制程序，通过最小二乘法对大型线性稀疏矩阵方程进行求解，由生成的三维网格数据绘制波速等值线图，从而实现了孔间介质波速层析成像的功能。根据成像图中低速带范围及波形分析确定了溶洞位置和充填情况，经checkerboard测试，并与钻探资料和地质雷达影像对比分析，结果表明跨孔地震CT技术用于地下岩溶探测时可靠性高，能够实现对孔间岩溶发育状况的超前预报。

岩溶是威胁地下工程施工安全的一大隐患。在南京地铁4号线工程勘探中发现，部分灰岩分布地段有大规模岩溶发育，洞径小则几十厘米，大至4～5m，岩溶类型为覆盖型(曾玉等，2011)，覆盖层为第四系沉积物。为避免施工过程中发生坍塌、涌水等工程事故，必须采取有效的手段查明岩溶的具体位置和规模。传统钻探方法成本高，钻孔间距和数量受到较大限制，因此单凭钻探资料并不能对测区地质异常体的分布情况做出整体判断(李世民，2004)。若采用地面物探方法，如垂直地震剖面(VSP)法和高频电磁波法，上覆土层可大量吸收弹性波和电磁波能量，探测深度达不到要求；直流电法因受到表面高阻层(水泥或沥青路面)屏蔽传导电流的影响，探测效果也不理想(陈易玖，1999)。相比较而言，跨孔地震CT技术在基岩中激发弹性纵波，探测深度不受影响，勘探准确性也得到了保证；此外，该技术可对孔间岩土体进行逐层扫描，并对目标探测区内介质的波速成像，使得孔间异常体能在图中直观地表现出来(李世民，2004)，因此地震CT技术用于探测孔间岩溶分布不失为一种较理想的方法。

对跨孔地震CT技术来说，成像分辨率是制约解释精度的重要因素。它既受激发频率、观测系统等硬件性能的影响，也受反演软件的影响(裴正林等，2002)。目前针对石油物探的反演软件已得到广泛应用，以走时反演为例，其核心算法从分辨率极低的反投影技术(BPT)，过渡到以联合迭代重建技术(SIRT)和正交分解(LSQR)为代表的迭代法(杨文采等，1993；杨文采，1997；罗炬等，2011)，到90年代初又提出了高分辨率级联法(杜剑渊等，1991；杨文采等，1993)。不足的是，这些方法或针对大型反演而设计(如迭代法)，或需要较高的计算成本(如级联法)，若盲目应用于规模不大的工程浅层勘探很可能事倍功半。因此，本文基于广义最小二乘法反演的思想(杨文采，1997)，在Matlab中编程反演计算介质波速，并将反演结果在二维剖面内成像。根据波速图中相对低速区的范围，对测区岩溶发育位置和大小作出判断。

图 1 南京地铁4号线一期工程桦墅站—仙林东站区间探测孔相对位置及跨孔测线分布情况Fig. 1 Plan of the relative locations of boreholes and survey lines for the cross-hole seismic CT test on metro line 4, Nanjing City

1 现场跨孔地震CT试验

1.1 场地工程地质与水文地质条件

南京地铁4号线一期工程桦墅站—仙林东站区间(矿山及明挖段)在钻探过程中发现规模较大的溶洞，溶洞为充填或半充填，充填物主要为粉质粘土和碎石土。钻孔深度范围内揭示的地层为上土下岩结构：上覆土层以杂填土、粉质粘土为主，厚度从0.5～20m不等；下伏基岩类型为中三叠统周冲组灰色中厚层灰岩和泥质灰岩。研究区地下水类型随地貌条件而异：在地势相对较高的构造剥蚀低山丘陵地貌单元内，地下水类型为基岩裂隙水和岩溶水，二者存在水力联系；在山体坡脚以下地段内，地下水类型为孔隙潜水、基岩裂隙水和岩溶水。

表1 试验仪器Table 1 Testing instruments

图 2 南京地铁4号线Q16G95—Q16Z73跨孔地震CT现场观测系统(两孔间距26m)Fig. 2 The field observation system of cross-hole ( in this case, boreholes are Q16G95 and Q16Z73 on metro line 4, Nanjing City) seismic CT (the distance between the two boreholes is 26m)

1.2 现场数据采集

测试区间内共12个勘探孔，见图1，分成11组进行跨孔地震CT试验。根据勘察要求，目标是对10～35m深度处的溶洞进行探测，由于孔间距不超过30m，因此为取得尽可能完全的投影数据，勘探孔深度应大于40m。

为具体说明跨孔地震CT的试验步骤，这里以Q16G95-Q16Z73两孔间的试验为例展开讨论，其余测线的处理办法与之类似。

试验开始前，首先应收集钻探资料，以便确定激发深度并检验成像结果的可靠性。钻孔资料显示，Q16G95孔基岩埋深7.4m，地下水位距地表8.8m，在9.1～12.2m段发现溶洞，26.7～29.8m段为破碎灰岩；Q16Z73孔基岩埋深0.5m，地下水位距地表11.4m，灰岩完整性较好。根据现场测线布置示意图(图2)，激发区间在Q16G95孔下8.5～20.5m。完成仪器连接后，检查RAS-24数字地震仪(表1)与计算机间的通信是否顺畅，系统触发状态是否正常，若不存在问题，即可将电火花震源(表1)置于Q16G95孔内20.5m深度处激发纵波；为了增加数据区域覆盖率，将12道三分量检波器以1m的道间距排列成串(表1)，底部固定在Q16Z73孔24.5m深度处，以便同时接收同一震源点激发的纵波信号；最后通过RAS-24地震仪和计算机将这12道信号记录下来。为保证探测分辨率，设定炮间距为1m，移动Q16G95孔中的震源至19.5m深度后再次激发地震纵波，并由地震仪和计算机记录下纵波信号。依上述流程逐次完成激—接工作，直至地震波覆盖目标探测区域为止。

图 3 直射线追踪与区域网格剖分示意图Fig. 3 Schematic illustration of an example straight ray path between the source and receiver and the mesh generation of the study area

图 4 反演计算流程图Fig. 4 Flowchart showing the steps of the inversion calculation计算公式参考杨文采(1997)Thecalculation formulas refer to Yang Wencai (1997)

值得注意的是，试验过程中为了使获取的信号有较高的信噪比，激发探头和检波器必须处于地下水或泥浆环境中，以保持传感器与井壁良好的耦合(杨文采，1993)。对耦合条件不佳的，可适当提高充电电压来加大激发能量。应当指出，加大震源能量不如改善耦合条件，这是CT采集的主要经验(杨文采等，1993)。此外，为保证地震仪能最大程度地记录高频信号，在记录系统中设置前放增益值为12dB来使波形有较大的振幅动态范围(章建宏等，2005)。

2 基于走时反演的波速层析成像

2.1 反演模型

如图3所示，在S处发射地震波，在R处接收。目标是通过第i个地震波的初至走时ti来反演断面内速度分布C(r)或慢度S(r) = 1/C(r)(i= 1, 2 ,…,N，r为空间位置矢量)。设地震波的传播路径为Li，则

(1)

式中，dl为弧长微元。

测区地质资料表明该区成像介质均一性好，弹性波主要穿透灰岩地层，由于不同深度的灰岩弹性波速变化很小(陆基孟等，1991)，因此波的衍射和散射效应可忽略不计。根据Fermat原理，路径Li可由图3中的直线唯一确定，这一方法称为直射线追踪法(Gu Hanming et al.，2006)。以矩形网格离散成像区域(图3)，并假设网格波速cj为常量(下标j为网格号，j= 1,2,…,M)，则该网格的慢度sj也是常量，代入(1)式得

(2)

式中，aij是第i条射线沿路径Lij经过第j个网格时所截线段长度(图3)。

考虑所有N条射线的情况，就得到了一个关于网格慢度的M元线性方程组，写成矩阵方程的形式，即

AS=T

(3)

2.2 数据处理

依据上述思想，以南京地铁4号线Q16G95—Q16Z73两孔间的现场测试结果为例说明数据处理方法，分为以下3个步骤：

图 5 实测信号经40～200Hz带通滤波后的波形记录Fig. 5 Waveform record processed by 40～200Hz band-pass filtering to the original signal

(1) 滤波。以Q16G95 孔内1处激发(图2)所获实测信号为例，进行40～200Hz的带通滤波处理。滤波后的地震记录如图5所示，其中横轴为接收通道编号，纵轴为时间(单位 s)，信号采样间隔0.125ms，数字信号记录格式为SEG-D 8038，模拟信号前放增益为12dB。

表 2 初至波走时拾取结果(单位ms)Table 2 First arrival picking results (unit ms)

图 6 南京地铁4号线Q16G95—Q16Z73断面反演计算点Fig. 6 The inversion calculation points on the cross section of Q16G95-Q16Z73 boreholes on metro line 4, Nanjing City

(2) 拾取初至波走时。由于初至波能量很强，有明显的波峰或波谷，因此正常情况下，初至时间取在初至波信号起跳位置。但对噪声较大以致初至波难以辨识的情形，应拾取第一个波峰(或波谷)对应的时间(李世民，2004)。因初至波走时拾取往往会引入观测误差，从而使方程(3)的解偏离真实情况，所以要对走时数据进行调试和校正，以确保反演结果的合理性。具体方法：①在单炮记录中，对初至走时无法确定的记录道，利用已确定的初至波走时值，通过样条插值重新确定，得到一连续的走时曲线；②将插值后的走时值映射到记录道中，在此附近寻找振幅最大的点作为最终初至波走时(潘树林等，2006)。按上述方法共拾取了156个初至波走时值，将它们合成初至波走时矢量T= (ti)156×1，如表2所示。

图 7 南京地铁4号线Q16G95—Q16Z73跨孔二维剖面纵波波速等值线图(图中所标数字为等值线对应的波速值)Fig. 7 Contoured interval P-wave velocities calculated by seismic CT in the cross-hole plane between boreholes Q16G95 and Q16Z73 on metro line 4, Nanjing City (the number marked in the figure represents velocity of the contour line)

(3) 成像平面网格离散。为使溶洞在成像时有较高的分辨率，必须指定合适的网格尺寸。一般来说，由于炮间距和道间距均为1m，当成像区被网格离散后，必须保证一个网格至少有一条射线经过(杨文采等，1993；张平松等，2012)。其次，成像物理分辨率的估计对网格尺寸的划分也具有参考意义。以图5为例，信号主频取为174Hz，由于纵波在溶洞内的速度不超过1km/s，为得到稳定的分辨率估计，波速取上限1km/s，求得波长为5.7m。因此，沿地震波传播方向溶洞的分辨率必大于1/4个波长，即0.7m(Fowler C M R，1990；Gu Hanming et al.，2006)。最终将网格尺寸取为1m×1m。

杨文采等(1993)提出，走时反演得到的图像分辨率约为3～5个像元宽度，且像元宽度最好小于探测目标特征尺度的1/3(像元在这里等同于网格)。而实际需查明对工程有影响的溶洞平均特征尺度大于3m，因此上述网格划分办法是合理的。令各网格中心点为反演计算点，则成像平面内312个计算点如图6所示，图中横轴代表以Q16G95为原点，跨孔断面内任一点的水平距离，纵轴代表埋深(注：断面内任一点埋深已经过地形校正，零点均设在地表处)。

表 3 由波速成像结果得到的溶洞方位及其充填情况Table 3 The position of karst caves and filled conditions by the results of velocity tomography

图 8 南京地铁4号线Q16G95—Q16Z73跨孔剖面checkerboard分辨率检测结果Fig. 8 Result of the checkerboard resolution test for the cross-hole plane between boreholes Q16G95 and Q16Z73 on metro line 4, Nanjing City

2.3 反演结果和讨论

确定了成像平面的尺寸、网格大小及初至波走时矢量后，调用事先编制好的解译程序，得到图7所示的二维剖面纵波波速等值线图，其坐标轴含义同图6，纵波波速在图中通过不同颜色予以区分。由图可知，波速分布于0.8～1.3km/s的区域为整个成像区的低速带。应当指出，成像平面经网格离散后，每一网格所反映的是该区域的平均波速，当地质异常体(如岩溶发育区)的尺寸小于网格尺寸时，其波速会受到周围介质的影响，因此图7所示低速带的边界处波速值偏大。

进一步分析，可从图7中获得下述探测结果：①根据低速带的边界确定出溶洞位置和平均洞径(见表3)；②低速带波速范围大于纵波在空气中的速度0.43km/s(陆基孟，1991)，因此解释为充填型溶洞。表3给出了包括Q16G95—Q16Z73跨孔在内的5组测线上探测到的溶洞的各项具体信息。

表 4 波速成像结果与钻探资料对比Table 4 The comparisonbetween the results of velocity tomography and drilling data

一般地球物理反演问题存在多解性，跨孔地震CT问题也不例外。为了评价反演结果的可靠性，拟综合checkerboard测试、钻探资料和地质雷达探测结果来说明该问题。

图 9 南京地铁4号线Q16G95—Q16Z73跨孔剖面地质雷达探测结果Fig. 9 Interpreted radar-gram showing distribution of karst caves for the cross-hole plane between boreholes Q16G95 and Q16Z73 on Metro 4, Nanjing City图中横轴为水平距离(单位m)；左侧纵轴代表反射波的回声时间(单位ns)；右侧纵轴代表探测深度(单位m)。计算探测深度时假定电磁波传播速度为0.1m/ns，并采用电磁波单程走时，即回声时间的一半The horizontal axis in the figure represents horizontal distance (unit m). The left and right vertical axes represent echo time of reflected wave (unit ns) and detecting depth (unit m), respectively. Assume the velocity of electromagnetic wave is0.1m/ns and use one way travel time (i.e. half of the echo time) when calculating detecting depth

根据checkerboard(或称检查板、棋盘格)法检测成像分辨率的原理(Humphreys E R et al.，1988；Zhao Dapeng et al.，1992)，在初始速度模型中加入绝对值为3%的正负相间的扰动量，对Q16G95—Q16Z73跨孔层析成像结果的分辨率进行分析。检测结果如图8所示，从图中直观地看出，成像结果中的主要特征都是可信的。

表 5 地质雷达探测到的溶洞方位及充填情况Table 5 The position of karst caves and filled conditions byground penetrating radar (GPR)

将钻孔周围解译出溶洞的4组数据与相应的钻探资料对比，对比结果如表4所示。从中不难看出，在容许误差内，解译结果与钻探资料有较高的一致性。

同时，还采用地质雷达方法对CT探测结果进行了验证。本次使用的地质雷达是瑞典Mala Geoscience公司生产的第三代ProEx雷达主机及50MHz超强地面耦合天线，测点间距0.2m，雷达有效探测深度为22m(图9)。根据反射波能量大小及同相轴特征，得到Q16G95—Q16Z73跨孔断面内溶洞分布如图9所示。表5详细列出了除该断面外其他4组跨孔测线的雷达解译结果，对比于由波速成像得到的解译结果(表3)，二者基本一致。

综上所述，利用介质波速层析成像方法可以有效给出溶洞的大小规模和发育位置，能够达到预期的分辨率要求，因而解译结果可信度高。若对测区内所有符合跨孔地震CT技术测试要求的场地进行探测，即可对测区岩溶发育情况作出准确、系统的评价。