中国的岩溶地貌发育广泛，约占国土面积的三分之一^[1]，主要分布在广西、广东、云南和贵州等地区。岩溶发育对地铁工程的不良影响主要体现在地面塌陷^[2]、隧道涌水突泥^[3]、盾构机沉陷^[4]等，同时严重危害既有建筑物和人民生命财产安全。因此，在岩溶地区地铁建设过程中，岩溶专项勘察工作必不可少。对岩溶勘察而言，地面物探成果因城市各种噪音干扰，精度不够。最直观的方法是钻探，但是钻探受城市场地条件和成本限制，孔距往往较大，具有一定的偶然性，不能完整刻画孔间的岩溶分布及连通等情况^[5]，造成勘察空白区。因此需要结合孔间地球物理勘探方法，进一步查明岩溶发育空间分布情况。

随着计算机技术及医学CT技术的发展，20世纪80年代CT技术（层析成像）也逐渐应用到地球物理领域。1984年SEG年会上首次公布地震层析成像的研究成果后，国内外开展了大量的试验和生产，该方法在油气勘探、裂缝监测、生产动态监测、矿产等领域得到广泛应用^[6-7]。近年来，国内在地震CT技术的方法、算法、软件和仪器等方面不断取得新进展。在实际工作中，地震CT技术以其施工方便、分辨率高、效果直观等优点，近几年被广泛用于城市地铁岩溶勘察中，取得了较好的探测效果^[8-11]。

根据算法的不同，地震CT主要分为反投影技术（back projection technique，BPT），代数重建技术（algebra reconstruction technique，ART），联合迭代重建技术（simultaneous iterative reconstruction technique，SIRT），共轭梯度最小二乘法（conjugate gradient least squares method，CGLS），正交分解最小二乘法（least squares with QR factorization，LSQR）以及最大熵法（maximum entropy method，ME）等^[7]，无论哪种算法，其本质都是求解一个系数矩阵为大型稀疏的方程组。

BPT算法虽然计算简单、快捷，但是分辨率低。ART和SIRT是与直射线相适应的，考虑到介质波速扰动对射线的作用，射线呈弯曲状态并且分布不均匀，系数矩阵常常病态，有时不能收敛。对于病态的方程组，基于外部迭代的CGLS和LSQR可以不断修正系数矩阵，得到更好的成像，尤其是LSQR具有更好的收敛性和稳定性^[12-13]。

根据地震波的传播理论，地震CT分为射线层析和波动方程层析两大类，射线层析理论基础稳定牢固，应用效果良好^[14]。当介质比较均匀时，地震波可以近似地看作沿直线传播，其反演结果与实际情况更接近。当介质非均匀分布时，波的传播路径就会发生弯曲，如果继续采用直线传播，反演误差必然很大。岩溶地区地质情况复杂，岩溶和围岩性质差异很大，属于波速差异较大的非均匀介质，采用弯曲射线追踪方法更符合实际情况。

1.   地质背景和地球物理条件

1.1   地质背景

研究区位于贵阳市云岩区，地形相对较为平坦，主要为溶蚀-侵蚀类的丘林沟谷地貌，地下水位埋深约6～15 m。根据勘察钻孔资料揭露，研究区的地层主要分为：①覆盖层主要为第四系杂填土（Q₄^ml）和红黏土（Q₄^el+dl），其中杂填土广泛分布于场地表层范围内，红黏土局部分布；②下伏基岩以可溶性灰岩为主，泥岩和砂岩局部分布，分别为二叠系茅山口组（P₁m）石灰岩，二叠系栖霞组（P₁q）石灰岩，石炭系摆佐、黄龙组（C₂b+h）石灰岩、泥岩，石炭系大塘组（C₁d）石灰岩，泥盆系蟒山群（D₂m）砂岩。

根据钻探成果显示，研究区内岩溶率高，岩溶发育情况复杂，形态主要以溶洞、溶沟（槽）、裂隙为主，溶洞一般呈全充填状态，充填物多为软塑粉质黏土或粉质黏土夹砂、夹灰岩碎屑等，具体示例详见图1。钻探在钻进过程中存在漏水现象，推测场地内岩溶可能具有连通性，但是需要采取其他勘察手段进一步判断岩溶连通情况。

图  1  全充填型溶洞示例

Figure  1.  Fully filled Karst caves

下载: 全尺寸图片 幻灯片

1.2   地球物理条件

工区位于城市繁华地段，地铁线位沿城市主干道敷设，人车流量大，地下电力、通讯、热力、排水等各种管道错综复杂。传统的地面物探方法受浅地表电磁、震动干扰，探测精度或深度往往不够，效果不佳。因此，选用以钻孔约束、地震CT技术为主的岩溶勘察方法，以详细查明盾构区间隧道影响范围内岩溶发育的埋深、规模和延伸等空间展布情况。

根据钻孔波测试结果：中风化石灰岩的纵波波速一般在4000 m/s以上，岩石破碎会导致波速下降，但一般不低于3300 m/s；溶蚀裂隙发育区的纵波波速一般在2300～3100 m/s；溶洞发育区（充泥）的纵波波速一般在1500～2500 m/s。岩溶发育区与基岩之间弹性波波速差异明显，这是开展地震CT法探测工作的地球物理前提。

2.   数据采集和预处理

2.1   现场采集

跨孔地震CT是分别在两个钻孔中完成地震的激发和接收，利用观测到的地震波走时、振幅和波形等信息，通过求解非线性方程组，重新构建孔间介质速度、刻画地质体构造的跨孔物探方法。

现场数据采集使用Geometrics公司生产的NZXP高精度地震仪、湘潭无线电有限责任公司生产的XW5512B电火花震源、广东物探院与中科院声学所共同研制的CH3型高灵敏度多道声波探头。发射电压为5 kV，采样间隔为20.833 μs，采样时长为100 ms，孔中炮间距和道间距均为1 m，图2为采集装置示意图。数据采集之前进行水听器一致性检测，采集一般选择噪音干扰较小的夜间进行，采集过程中采用单炮信号重复叠加提高信噪比。由于场地内标高1 111 m以上地层无地下水，数据质量差，未采集该深度范围的数据。

图  2  采集装置示意图

Figure  2.  Schematic map of the acquisition device

下载: 全尺寸图片 幻灯片

钻孔分别沿隧道左右中线布置，受场地条件限制，部分钻孔稍作偏移，孔深根据隧道埋深确定，钻孔间距约20 m，CT测线布设方式如图3所示。

图  3  测线布置示意图

Figure  3.  Schematic map of the survey line layout

下载: 全尺寸图片 幻灯片

2.2   数据预处理

采集到原始数据后，先对每个CT剖面数据进行道编辑、建立观测系统、滤波去噪等预处理，然后根据地震信号的振幅、频率和相位的变化拾取初至获得走时文件，准确拾取地震波初至时间是最关键的一步，在很大程度上制约整个剖面的最终处理结果。

图4为预处理后的典型单炮波形记录图，从图中可以看出，记录深度范围内浅部3 m能量较小，其余深度剖面波形完整、光滑，初至明显。

图  4  典型单炮波形记录图

Figure  4.  Recording map of a typical single shot waveform

下载: 全尺寸图片 幻灯片

3.   反演方法对比

根据上文所述，BPT、ART和SIRT三种算法是与直射线追踪相适应的，我们经过反演试算，3种方法对于岩溶地层确实得不到理想结果，本文不再研究。本文主要对CGLS和LSQR算法进行研究对比，两种算法本质都是求方程${\boldsymbol{Ax}}={\boldsymbol{b}} $的最小二乘问题$\min\big\|{\boldsymbol{Ax}}-{\boldsymbol{b}}\big\|_2 $，CGLS算法是一种解决不对称线性方程和最小二乘问题的共轭梯度法，LSQR是基于双对角化的Golub和Kahan法^[15]。

共轭梯度法最早由Hestenes等^[16]在1952年提出，共轭梯度法作为迭代法，首先需要计算误差向量和梯度向量，并初始化一个搜索方向向量。然后根据搜索方向向量，利用线性方程组的解法求出一个步长，将参数向量更新为新的值。接着计算新的误差向量和梯度向量，并重新搜索方向向量，重复这个过程，直到达到收敛条件为止。该方法介于最速下降法与牛顿法之间，利用1阶导数信息，克服了最速下降法收敛慢的缺点，同时避免了牛顿法需要存储和计算海森矩阵并求逆的缺点。

LSQR算法最早由Paige等^[17]在1982年提出，利用Lanczos方法求解最小二乘问题，在求解过程中会用到QR因子分解，因此称为LSQR方法（least squares QR-factorization）。在方程Ax=b中，系数矩阵A为n×n阶实数对称方阵，x和b为n维向量，x为待求未知量。设v₁，v₂，···，v_m为n维空间的m个无关向量，令V_m=(v₁，v₂，···，v_m)，然后寻找近似解x_m，使得(V_m^TAV_m)x_m=V_m^Tb，同时利用Lanczos方法构造向量(v₁，v₂，···，v_m)，使得矩阵V_m^TAV_m具有三角对称形式。当A为不对称方阵时，需将A对角化，转化为对称方阵，然后再应用Lanczos方法。该算法在迭代过程中只涉及到非零元素，所需存储空间小，运算速度快，最适用于求解奇异或病态方程。程序设计主要步骤如下：

两种算法类似，但对不同的数据拟合会表现出不同的精确性^[15]，下面通过岩溶区实例对比选取合适的算法。以钻孔ZK-85和ZK-86地震剖面数据为例，采用弯曲射线追踪方法，以钻孔边界条件为约束，选取相同的初始模型，网格大小为0.5 m×0.5 m。通过反演试算，二者迭代1次耗时都在22 s左右，迭代速度接近。二者速度拟合误差都在8％ 左右，CGLS算法需迭代8次才能趋于稳定，LSQR算法迭代5次即可稳定，说明LSQR算法具有更快的收敛性（图5）。

图  5  LSQR算法程序设计主要步骤

Figure  5.  The main steps of LSQR algorithm program design

下载: 全尺寸图片 幻灯片

两种算法反演结果如图6所示，从图中反演结果可以看出，二者波速剖面形态基本一致。

图  6  ZK-85～ZK-86地震剖面反演结果对比

Figure  6.  Comparison of the ZK-85～ZK-86 seismic profile inversion results

下载: 全尺寸图片 幻灯片

对于岩溶低速区，LSQR算法反映更明显，速度值更小，其揭示的小规模低速异常R8也得到了钻探验证（图7），说明LSQR算法对于异常的识别更灵敏，反演结果更加可靠和精细。因此本文最终选用LSQR算法，以求得到更好的反演结果。

图  7  CT-82～CT-86地震剖面解释图

Figure  7.  Interpretation map of the CT-82～CT-86 seismic profile

下载: 全尺寸图片 幻灯片

4.   地质解释及钻探验证

4.1   地质解释

解释原则：①根据钻孔资料约束，对速度剖面图进行对比分析，确定不同地层的波速范围及特征；②在地震剖面图中，出现较明显低速闭合或半闭合区域，一般推断解释为岩溶发育区。

以CT-82～CT-86地震CT剖面为例，该剖面位于钻孔ZK-82、ZK-83、ZK-84、ZK-85、ZK-86之间，剖面长度约80 m，反演结果如图7所示。根据CT-82～CT-86地震剖面波速分布特征分析：岩溶发育区与完整基岩波速差异明显，共揭示9处岩溶发育区。其中R1、R2、R4、R8、R9发育规模较小，呈闭合不规则形状，波速范围为2300～3 100 m/s（黄、绿色），解释为溶蚀裂隙发育区；R3、R5、R6、R7发育规模较大，以横向条带状为主，波速范围为1500～3 100 m/s，解释为溶洞发育区，溶洞发育区与基岩存在波速过渡区，推测为溶蚀裂隙，溶洞发育区内部充泥，波速小于2 300 m/s（蓝色）。

地震CT结果揭露工区岩溶发育区在纵向埋深上均有分布，水平方向比垂直方向发育规模大，表明岩溶异常在横向上较垂向更具备连通性和延展发育的特点，易形成地下岩溶水渗流通道，为下一步设计和施工提供可靠的地质依据。

4.2   钻探验证

为了验证地震CT剖面对岩溶发育区刻画的准确性，在ZK-85和ZK-86之间布置验证孔YZK-1。钻探揭示的3处岩溶发育区（即图6中的黄色粗线部分）和地震CT剖面揭示的R6、R7、R8岩溶区基本一致，验证结果对比分析详见表1。

表  1  YZK-1钻孔验证结果对比分析

Table  1.  Comparative analysis of YZK-1 drilling verification results

序号	岩溶编号	反演高度/m	钻探验证高度/m	反演相对误差/％	结果分析
1	R6	3.3	4.6	-28.2	反演结果偏小
2	R7	3.2	1.9	68.4	反演结果偏大
3	R8	1.2	2.3	-47.8	反演结果偏小

下载: 导出CSV 

| 显示表格

从验证情况可以看出，反演相对误差绝对值随岩溶规模增大而减小，说明岩溶区越大，反演越精确；对于小规模岩溶异常R8仍然有较准确反映，可达米级探测精度，说明LSQR算法具有较高的分辨率和可靠性，可以获得满意的速度图像。

5.   结论

本文通过贵阳市地铁3号线勘察中跨孔地震CT技术的应用研究，得到如下结论：

（1）研究区岩溶发育在纵向埋深上分布具有随机性，在横向上具备连通性和延展发育的特点，易形成地下岩溶水渗流通道。

（2）通过对不同算法分析，尤其是CGLS和LSQR算法的实例对比研究，表明在相同初始模型和反演条件下，LSQR算法对于低速异常的识别更灵敏和精确，具有更好的收敛性和可靠性。

（3）钻探验证表明，采用LSQR算法和弯曲射线追踪法，跨孔地震CT探测可达米级精度，能够查明钻孔之间岩溶埋深、规模、延伸等情况，说明基于LSQR算法的地震CT技术探测岩溶是切实可行的。但是反演结果的相对误差仍然较大，有待于探索更有效的反演方法或约束条件来改善反演精度。