Inversion Imaging Research of Cross Hole Seismic CT in Complex Karst Areas: Guiyang Metro Line 3 under Construction Presented as an Example
-
摘要:
本文以贵阳在建地铁3号线为例,详细研究岩溶复杂地区勘察过程中跨孔地震CT技术的应用,揭露研究区岩溶发育特征。通过对不同算法的对比分析表明,相比于CGLS算法,LSQR算法对于低速异常的识别更灵敏和精确,具有更好的收敛性和可靠性。钻探验证表明,采用LSQR算法和弯曲射线追踪法进行地震层析成像反演,跨孔地震CT探测可达米级精度,获得满意的速度图像,能够查明钻孔之间岩溶埋深、规模、延伸等情况,弥补钻探的局限性,说明基于LSQR算法的地震CT技术探测岩溶是切实可行的。
Abstract:This article used Guiyang Metro Line 3 under construction as an example to study the application of cross hole seismic computed tomography (CT) technology in the exploration process of Karst complex areas, thereby revealing the characteristics of Karst development in the study area. Through comparative analysis of various algorithms, it was determined that compared to the conjugate gradient least squares (CGLS) algorithm, the sparse linear equations and least square problems (LSQR) algorithm is more sensitive and accurate in identifying low-speed anomalies, and has better convergence and reliability. Drilling verification shows that using the LSQR algorithm and bending ray tracing method for seismic tomography inversion can achieve meter level accuracy in cross hole seismic CT detection, obtain high-quality velocity images, and identify the depth, scale, and extension of Karst between boreholes, neutralizing the limitations of drilling. This indicates that seismic CT technology based on the LSQR algorithm is practical and feasible for detecting Karst.
-
Keywords:
- seismic CT technique /
- LSQR algorithm /
- wave velocity /
- drilling /
- Karst survey
-
中国的岩溶地貌发育广泛,约占国土面积的三分之一[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。根据勘察钻孔资料揭露,研究区的地层主要分为:①覆盖层主要为第四系杂填土(Q4ml)和红黏土(Q4el+dl),其中杂填土广泛分布于场地表层范围内,红黏土局部分布;②下伏基岩以可溶性灰岩为主,泥岩和砂岩局部分布,分别为二叠系茅山口组(P1m)石灰岩,二叠系栖霞组(P1q)石灰岩,石炭系摆佐、黄龙组(C2b+h)石灰岩、泥岩,石炭系大塘组(C1d)石灰岩,泥盆系蟒山群(D2m)砂岩。
根据钻探成果显示,研究区内岩溶率高,岩溶发育情况复杂,形态主要以溶洞、溶沟(槽)、裂隙为主,溶洞一般呈全充填状态,充填物多为软塑粉质黏土或粉质黏土夹砂、夹灰岩碎屑等,具体示例详见图1。钻探在钻进过程中存在漏水现象,推测场地内岩溶可能具有连通性,但是需要采取其他勘察手段进一步判断岩溶连通情况。
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以上地层无地下水,数据质量差,未采集该深度范围的数据。
钻孔分别沿隧道左右中线布置,受场地条件限制,部分钻孔稍作偏移,孔深根据隧道埋深确定,钻孔间距约20 m,CT测线布设方式如图3所示。
2.2 数据预处理
采集到原始数据后,先对每个CT剖面数据进行道编辑、建立观测系统、滤波去噪等预处理,然后根据地震信号的振幅、频率和相位的变化拾取初至获得走时文件,准确拾取地震波初至时间是最关键的一步,在很大程度上制约整个剖面的最终处理结果。
图4为预处理后的典型单炮波形记录图,从图中可以看出,记录深度范围内浅部3 m能量较小,其余深度剖面波形完整、光滑,初至明显。
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为待求未知量。设v1,v2,···,vm为n维空间的m个无关向量,令Vm=(v1,v2,···,vm),然后寻找近似解xm,使得(VmTAVm)xm=VmTb,同时利用Lanczos方法构造向量(v1,v2,···,vm),使得矩阵VmTAVm具有三角对称形式。当A为不对称方阵时,需将A对角化,转化为对称方阵,然后再应用Lanczos方法。该算法在迭代过程中只涉及到非零元素,所需存储空间小,运算速度快,最适用于求解奇异或病态方程。程序设计主要步骤如下:
两种算法类似,但对不同的数据拟合会表现出不同的精确性[15],下面通过岩溶区实例对比选取合适的算法。以钻孔ZK-85和ZK-86地震剖面数据为例,采用弯曲射线追踪方法,以钻孔边界条件为约束,选取相同的初始模型,网格大小为0.5 m×0.5 m。通过反演试算,二者迭代1次耗时都在22 s左右,迭代速度接近。二者速度拟合误差都在8% 左右,CGLS算法需迭代8次才能趋于稳定,LSQR算法迭代5次即可稳定,说明LSQR算法具有更快的收敛性(图5)。
两种算法反演结果如图6所示,从图中反演结果可以看出,二者波速剖面形态基本一致。
对于岩溶低速区,LSQR算法反映更明显,速度值更小,其揭示的小规模低速异常R8也得到了钻探验证(图7),说明LSQR算法对于异常的识别更灵敏,反演结果更加可靠和精细。因此本文最终选用LSQR算法,以求得到更好的反演结果。
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 反演结果偏小 从验证情况可以看出,反演相对误差绝对值随岩溶规模增大而减小,说明岩溶区越大,反演越精确;对于小规模岩溶异常R8仍然有较准确反映,可达米级探测精度,说明LSQR算法具有较高的分辨率和可靠性,可以获得满意的速度图像。
5. 结论
本文通过贵阳市地铁3号线勘察中跨孔地震CT技术的应用研究,得到如下结论:
(1)研究区岩溶发育在纵向埋深上分布具有随机性,在横向上具备连通性和延展发育的特点,易形成地下岩溶水渗流通道。
(2)通过对不同算法分析,尤其是CGLS和LSQR算法的实例对比研究,表明在相同初始模型和反演条件下,LSQR算法对于低速异常的识别更灵敏和精确,具有更好的收敛性和可靠性。
(3)钻探验证表明,采用LSQR算法和弯曲射线追踪法,跨孔地震CT探测可达米级精度,能够查明钻孔之间岩溶埋深、规模、延伸等情况,说明基于LSQR算法的地震CT技术探测岩溶是切实可行的。但是反演结果的相对误差仍然较大,有待于探索更有效的反演方法或约束条件来改善反演精度。
-
表 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 反演结果偏小 -
[1] 袁道先. 现代岩溶学在中国的发展[J]. 地质论评, 2006, (6): 733−736. DOI: 10.3321/j.issn:0371-5736.2006.06.002. YUAN D X. The development of modern Karstology in China[J]. Geological Review, 2006, (6): 733−736. DOI: 10.3321/j.issn:0371-5736.2006.06.002. (in Chinese).
[2] 蒙彦, 雷明堂. 岩溶塌陷研究现状及趋势分析[J]. 中国岩溶, 2019, 38(3): 411−417. MENG Y, LEI M T. Analysis of situation and trend of sinkhole collapse[J]. Carsologica Sinica, 2019, 38(3): 411−417. (in Chinese).
[3] 张可能, 张岳, 廖阳, 等. 贵阳某地铁车站岩溶发育特征及突水模式分析[J]. 中国岩溶, 2018, 37(2): 300−306. ZHANG K N, ZHANG Y, LIAO Y, et al. Analysis on Karst development and water burst in a subway station[J]. Carsologica Sinica, 2018, 37(2): 300−306. (in Chinese).
[4] 郭红梅, 马培贤. 济南地铁建设中的岩土工程问题分析[J]. 西部探矿工程, 2012, 24(6): 9−12. [5] 张雨飞. 井间地震层析成像技术在地铁岩溶勘察中的应用[J]. 工程地球物理学报, 2019, 16(5): 591−595. ZHANG Y F. The application of cross-well seismic tomography technique to Karst exploration of metro engineering[J]. Chinese Journal of Engineering Geophysics, 2019, 16(5): 591−595. (in Chinese).
[6] 陈仲侯, 傅唯一. 浅层地震勘探[M]. 成都: 成都地质学院, 1986. [7] 杨文采, 李幼铭. 应用地震层析成像[M]. 北京: 地质出版社, 1993. [8] 王玉海. 井地地震CT技术在岩溶勘察中的应用研究[D]. 南京: 南京大学, 2011. WANG Y H. Application and research of well-ground seismic CT technology to Karst caves survey[D]. Nanjing: Nanjing University, 2011. (in Chinese).
[9] 史晓忠. 跨孔地震CT在地铁岩溶勘察中的应用[J]. 城市道桥与防洪, 2018, (7): 306−309. [10] 张腾飞. 跨孔地震CT技术在大连地铁五号线岩溶段探测的应用研究[D]. 成都: 西南交通大学, 2019. ZHANG T F. Application research of cross-hole seismic CT technology in the detection of Karst section of Dalian Metro Line 5[D]. Chengdu: Southwest Jiaotong University, 2019. (in Chinese).
[11] 林松, 王薇, 金聪, 等. 地震CT在岩溶精细探测中的应用与探讨—以深圳地铁14号线为例[J]. 科学技术与工程, 2019, 19(25): 18-23. LIN S, WANG W, JIN C, et al. Application and discussion of seismic CT in detailed Karst detection: A case of shenzhen metro line 14[J]. Science Technology and Engineering, 2019, 19(25): 18-23. (in Chinese).
[12] 张东, 乔友峰, 姜麟舜, 等. 地震层析成像中LSQR算法的快速求解[J]. 物探化探计算技术, 2011, 33(6): 632−635. [13] 杨薇, 刘四新, 冯彦谦. 跨孔层析成像LSQR算法研究[J]. 物探与化探, 2008, 32(2): 199−202. YANG W, LIU S X, FENG Y Q. A study of the LSQR algorithm for cross-hole tomography[J]. Geophysical & Geochemical Exploration, 2008, 32(2): 199−202. (in Chinese).
[14] 胡刚, 段宝平, 何正勤, 等. 利用井间地震层析成像方法探测隐伏断层[J]. CT理论与应用研究, 2015, 24(3): 345−355. DOI: 10.15953/j.1004-4140.2015.24.03.03. HU G, DUAN B P, HE Z Q, et al. Application of cross-well seismic tomography in detecting buried fault[J]. CT Theory and Applications, 2015, 24(3): 345−355. DOI: 10.15953/j.1004-4140.2015.24.03.03. (in Chinese).
[15] 陈善雄, 熊海灵, 廖建伟, 等. 一种基于CGLS和LSQR的联合优化的匹配追踪算法[J]. 自动化学报, 2018, 44(7): 1293−1303. CHEN S X, XIONG H L, LIAO J W, et al. A matching pursuit algorithm of jointing optimization based on CGLS and LSQR[J]. Acta Automatica Sinica, 2018, 44(7): 1293−1303. (in Chinese).
[16] HESTENES M R, STIEFEL E. Methods of conjugate gradients for solving linear systems[J]. Journal of Research of the National Bureau of Standards. 1952, 49: 409-436.
[17] PAIGE C C, SAUNDERS M A. LSQR: An algorithm for sparse linear equations and sparse least squares[J]. ACM Transactions on Mathematical Software, 1982, 8(1): 43−71. DOI: 10.1145/355984.355989.