跨孔地震CT探测基岩面附近岩溶研究

2022-02-21李卫卫蒙爱军

工程地球物理学报 2022年1期

李卫卫, 熊鑫, 蒙爱军

(1.浙江华东建设工程有限公司，浙江杭州 311202；2.浙江省工程物探勘察设计院有限公司,浙江杭州 315012)

1 引言

很多建筑物的持力层都落在基岩上，若基岩面附近存在岩溶则是危害建筑物基础的重大隐患。高密度电法[1-5]、跨孔地震 CT[6]、地质雷达[7]以及综合物探方法[8,9]广泛应用于岩溶探测。当基岩上部有覆盖层时，借助于勘探孔进行地震跨孔CT探测岩溶是一种很好的选择。谈顺佳等[10]和杨永龙等[11]在对一般建筑物基础的岩溶分布探测中采用了跨孔地震CT技术，取得了较好的效果。段成龙等[12]和郭书兰等[13]在地铁基础岩溶探测中采用了跨孔地震CT技术，根据成像图中低速带范围及波形分析确定了溶洞位置和充填情况。朱文仲等[14]对弹性波CT技术几个重要问题进行了正反演研究，就弹性波 CT技术的分辨率问题、“低速色盲”问题和像元划分问题，提出了各种观测系统的最高分辨率估算公式，指出了二侧透视观测系统的重大缺陷、弹性波CT技术的“低速色盲”本质以及划分像元的一般原则。石振明等[15]基于波动方程的时域有限差分对跨孔地震岩溶探测进行正演模拟，从模拟的地震记录中提取初至旅行时，最后以基于程函方程的初至层析方法反演得到钻孔间的地层速度模型。潘纪顺等[16]利用 Matlab 编程对堤防工程进行地震断面CT数值模拟，通过ART和SIRT重建方法进行反演，分析得出二者在堤防工程CT应用中的优缺点及适用条件。

为了更好、更精确地探测岩溶，本文从基岩面附近的多岩溶探测入手，依据射线理论进行正反演模拟研究，讨论岩溶与基岩面的距离、围岩与岩溶速度差异大小以及观测系统对反演结果的影响。本文还对将岩土介质一起加入反演的情况进行研究。

2 基本的算法原理

本文采用韩佩恩等[17]提出的改进的Moser方法[18]进行初至时间场计算，利用Vidale[19]和Asakawa和Kawanaka[20]提出的走时线性插值算法确定射线路径。

韩佩恩等在Moser算法的基础上，提出了加入动态节点在同样的网格间距下提高计算初至时间精度的方法，如图1所示，在加入动态节点情况下采用Moser算法计算所有节点的初至旅行时，满足根据已知网格化的速度模型正演计算激发点至接收点旅行时的要求。Moser算法难以提供反演要求的准确射线路径，走时线性插值算法则可以有效解决这个问题，该算法分为两步：第一步是向前算法，即计算从震源出发到达网格各个节点的地震波初至旅行时；第二步是向后算法，即根据各个节点的旅行时和接收点位置反向计算接收点至震源的确切射线路径。具体计算过程参考相关文献[19,20]。

图1 改进的Moser方法示意图Fig.1 Schematic diagram of the improved Moser method

反演采用联合迭代重建法(SIRT)，其基本原理为：将所有射线的投影差值计算出来，然后将各个射线的投影差值反投影到单元网格上；将网格单元排序，把所有穿过网格单元的射线的走时投影汇总起来，利用这些汇总数据对网格单元地震波慢度(速度的倒数)进行修正，直到获得的结果满足预设的精度或者达到设定的最大迭代次数。在做跨孔地震波CT过程中，网格化所需要勘察的区域，根据先验信息给网格化后的每个网格单元赋予一个慢度初值。再计算每条经过网格的初至波射线的走时差，然后将走时差分配到每个网格单元中，并对所有分配到的走时差进行统计，依据平均值来修正网格单元的地震波慢度。

设lkj为第k条射线穿过第j个网格(按照一维编序号)的射线长度，设Lk为第k条射线长度，即

(1)

式(1)中Mk为第k条射线所经过的模型网格个数。

设Δtk为第k条射线的观测值与理论旅行时之差，设βk为第k条射线的模型修正系数，即

(2)

则第k条射线对第j个网格的慢度修正量为(μ为微小的量)：

(3)

所以，所有射线对第j个网格的慢度修正量为：

(4)

式(4)中Nj表示第j个网格所穿过的射线条数。

模型慢度修正(反演)为：

(5)

其中q为迭代次数。

3 理论模型的正反演计算

本文理论模拟的研究思路是建立多岩溶不同排列组合的正方形网格化孔间地震纵波速度剖面模型，根据模型设计跨孔地震波初至CT成像技术的观测系统，如图2所示，采用改进的Moser方法[17]计算所有炮点(激发点)至检波点(接收点)的初至旅行时，以此观测系统和初至旅行时数据反演两孔之间的速度分布。反演的步骤是：①根据观测系统和初至旅行时建立初始反演模型；②采用改进的Moser方法[18]计算所有炮检点的初至旅行时；③采用Vidale[19]提出的走时线性插值向后算法计算所有炮检点的射线路径；④修改每个网格的速度(慢度)值；⑤判断是否满足反演终止条件；⑥如果不满足终止条件，则继续②～⑤的步骤。

图2 跨孔地震CT观测系统示意图Fig.2 Schematic diagram of cross-hole seismic CT observation system

本文建立初始模型的思路如下：设Ns个激发点，其纵坐标为ZSi(i=1, 2, 3,...,Ns)，与其对应纵坐标的接收点为ZPj(j=1, 2, 3,...,Np)，地震波初至的最小旅行时为tij，假设能够找到ZPk=ZSi，将其简称为Zi(i=1, 2, 3,...,Ns)，假设钻孔间距为dx，按照直射线传播假设，则传播速度为V(Zi)=dx/tik(i=1, 2, 3,...,Ns)，用速度序列V(Zi)进行垂向线性插值，建立横向均匀垂向连续介质模型。

本次理论模拟的孔间距固定为15 m，孔深为20 m，左边为激发孔，从孔口到孔底1 m一个激发点，右边为接收孔，从孔口到孔底1 m一个接收点，共计21个接收点。观测系统采用60度张望角，排列长度(接收段长度)L=2×15 m×tan60°≈52 m>2×20 m，所以21个接收点均可接收到每一个激发点出发的地震波。模型的网格化间距为0.25 m×0.25 m，反演模型的网格化间距均为0.5 m×0.5 m，速度剖面的色标值是动态的，未统一在同一个值域，通常情况下，围岩的速度设为3 000 m/s，岩溶填充物的速度为1 500 m/s。

3.1 岩溶与基岩面的距离对反演结果的影响

设计一组两个水平相邻的岩溶模型，半径均为R=1.5 m，岩溶中心间距均为L=4 m，岩溶中心埋深(与基岩面的距离)分别为H=4 m、5.5 m、7 m。正反演结果如图3～图5所示。可以看出：当岩溶中心与基岩面的垂直距离不大于5.5 m时,反演的结果只能看到一个岩溶体，埋深较准，水平位置在两个岩溶中间(图3d，图4d)；当岩溶中心与基岩面的垂直距离达到7 m时(图5d)，则看到两个岩溶体的准确中心位置，紧靠接收孔的假异常则可以通过接收钻孔排除其为岩溶体可能。当岩溶中心与基岩面的垂直距离为7 m半径缩为1.25 m时(图6d)，也能探测到两个岩溶体的存在。以上研究表明，当两个岩溶体距离基岩面太近时，无法解析多岩溶体的空间分布；当两个岩溶体与基岩面相隔一定距离时，则可以分辨多岩溶体的空间位置。

图3 模型1正反演结果(H=4 m, L=4 m, R=1.5 m)Fig.3 The forward and inverse result for model 1(H=4 m, L=4 m, R=1.5 m)

图4 模型2正反演结果(H=5.5 m, L=4 m, R=1.5 m)Fig.4 The forward and inverse result for model 2 (H=5.5 m, L=4 m, R=1.5 m)

图5 模型3正反演结果(H=7 m, L=4 m, R=1.5 m)Fig.5 The forward and inverse result for model 3 (H=7 m, L=4 m, R=1.5 m)

图6 模型4正反演结果(H=7 m, L=4 m, R=1.25 m)Fig.6 The forward and inverse result for model 4 (H=7 m, L=4 m, R=1.25 m)

3.2 岩溶横向分辨能力探讨(贯穿和非贯穿模型，岩溶间距)

为了探讨跨孔地震CT的横向分辨能力，本文设计了两类模型，一是两个水平相邻的岩溶，半径均为R=1.5 m，岩溶中心间距L和中心埋深H可变；二是一个水平矩形条状模型，高度h=1.5 m，其中心埋深H和宽度w可变。图7是将图3(a)中的岩溶中心间距从4 m增加到5 m的正反演结果，能明显看到两个岩溶体的存在，右边的岩溶体位置误差较大；如图8所示，将图7(a)的岩溶中心埋深从4 m减到2.5 m时，已经探测不到两个岩溶体的准确位置；图9是将图8(a)中的岩溶中心间距从5 m增加到9 m的正反演结果，能明显地看到两个岩溶体的存在，位置基本正确。以上研究表明：如果岩溶体距离基岩面较近，当岩溶体的中心距增加到一定程度时，跨孔地震CT也可以分辨多岩溶体的存在。

图7 模型5正反演结果(H=4 m, L=5 m, R=1.5 m)Fig.7 The forward and inverse result for model 5 (H=4 m, L=5 m, R=1.5 m)

图8 模型6正反演结果(H=2.5 m, L=5 m, R=1.5 m)Fig.8 The forward and inverse result for model 6 (H=2.5 m, L=5 m, R=1.5 m)

图9 模型7正反演结果(H=2.5 m, L=9 m, R=1.5 m)Fig.9 The forward and inverse result for model 7 (H=2.5 m, L=9 m, R=1.5 m)

如图10(a)和图11(a)所示，水平矩形条岩溶模型的宽度均为w=8 m，中心埋深分别为H=5 m、H=10 m，反演结果如图10(d)和图11(d)所示，中心埋深反演地都很准确，但左右边缘位置都相差太大。图12是将图10(a)的水平矩形条岩溶模型的宽度增至15 m(横向贯穿)时的正反演结果，可以看出跨孔地震CT能够对横向贯穿性岩溶进行准确探测。以上研究表明：采用跨孔地震CT技术对条带岩溶进行探测，当条带岩溶水平方向贯穿到两个钻孔时该方法可以准确探测，而对于水平方向非贯穿岩溶则只能探测其中心埋深而不能探测其左右边缘。

图10 模型8正反演结果(H=5 m, w=8 m, h=1.5 m)Fig.10 The forward and inverse result for model 8 (H=5 m, w=8 m, h=1.5 m)

图11 模型9正反演结果(H=10 m, w=8 m, h=1.5 m)Fig.11 The forward and inverse result for model 9 (H=10 m, w=8 m, h=1.5 m)

图12 模型10正反演结果(H=5 m, w=15 m, h=1.5 m)Fig.12 The forward and inverse result for model 10 (H=5 m, w=15 m, h=1.5 m)

3.3 围岩与岩溶充填物速度差异大小对反演的影响

本文设计一组两个水平相邻的岩溶体，半径均为R=1.5m，岩溶中心间距均为L=5 m，岩溶中心埋深为H=2.5 m，一共三个模型，通过调整围岩与岩溶填充物速度大小关系，探讨围岩与岩溶充填物速度差异大小对跨孔地震CT反演的影响。图13和图14反演结果一定程度上恢复了两个岩溶的空间位置，但图15中两个岩溶体则无法分辨出来。岩溶充填物与围岩的速度差相对于围岩速度分别为(2 500 m/s-1 500 m/s)/2 500 m/s=0.4、(2 000 m/s-1 500 m/s)/2 000=0.25 m/s、(2 500 m/s-2000 m/s)/2 500 m/s=0.2，依据反演结果可以看出，岩溶充填物与围岩的速度差与围岩速度的比值越大，跨孔地震CT技术越易于分辨两个水平相邻的岩溶体。

图13 模型11正反演结果(H=2.5 m, L=5 m, R=1.5 m)Fig.13 The forward and inverse result for model 11 (H=2.5 m, L=5 m, R=1.5 m)

图14 模型12正反演结果(H=2.5 m, L=5 m, R=1.5 m)Fig.14 The forward and inverse result for model 12 (H=2.5 m, L=5 m, R=1.5 m)

图15 模型13正反演结果(H=2.5 m, L=5 m, R=1.5 m)Fig.15 The forward and inverse result for model 13 (H=2.5 m, L=5 m, R=1.5 m)

图16 模型14正反演结果(H=2.5 m, L=5 m, R=1.5 m)Fig.16 The forward and inverse result for model 14 (H=2.5 m, L=5 m, R=1.5 m)

图17 模型15正反演结果(H=4 m, L=5 m, R=1.5 m)Fig.17 The forward and inverse result for model 15 (H=4 m, L=5 m, R=1.5 m)

3.4 覆盖层对反演的影响

前文研究表明，当两个岩溶体中心间距较小且又接近于基岩面时，在基岩中布置观测系统无法分辨这两个岩溶体。基于此，本文设计含覆盖层的模型，探讨在覆盖层中也布置激发点和接收点能否探测清楚基岩面附近的多个(大于等于2)岩溶体。参考前人研究结果[11, 13]，本文将覆盖层的速度参数亦作为反演参数。

图18 模型16正反演结果(H=5.5 m, L=5 m, R=1.5 m)Fig.18 The forward and inverse result for model 16 (H=5.5 m, L=5 m, R=1.5 m)

图19 模型17正反演结果Fig.19 The forward and inverse result for model 17

本文设计了4个模型，覆盖层厚度都是5 m，前三个模型都是基岩里含两个水平相邻的岩溶体，半径均为R=1.5 m，岩溶中心间距均为L=5 m，岩溶中心埋深(与基岩面的距离)分别为H=2.5 m、4 m、5.5 m，第4个模型是基岩中不含岩溶。基于图16(a)、图17(a)、图18(a)模型做正演，它们的反演结果如图16(d)、图17(d)、图18(d)所示，这3张速度剖面中都看不到两个岩溶体的存在，它们都与图19(d)相似，即都在基岩面下方约5 m处模型水平中心位置有一个低速假异常体。以上研究表明：如果将激发点和接收点也布置在覆盖层中，那么该观测系统严重影响了反演的结果；若将激发点和接收点布置在基岩中(图7d)则能准确探测到两个岩溶体。