城市地下孤石的地震波场响应特征及其在地震层析成像中的应用研究

朱智清，王道坤，汪伟民，吴殿伟，施卫力，张 赛

(1.机械工业勘察设计研究院有限公司，陕西西安 710043；2.安徽省皖北煤电集团有限责任公司，安徽宿州 234000；3.中铁上海工程局集团有限公司，上海 201900；4.安徽省建设工程测试研究院有限公司，安徽合肥 230051)

0 引言

随着我国经济的高速发展，各种类型的城市建筑物数量逐年增长，城市土地资源紧缺，城市地下空间利用与开发是缓解这一状况的重要措施(赵镨等，2017；杨洋等，2019；阎浩等，2020)。然而城市浅层地质条件复杂，地下孤石等地质灾害体较为常见且危害巨大，为降低城市地下空间开发过程中的安全隐患，必须详细查明地下孤石的分布情况(施有志等，2019)，目前孤石探测仍处于发展阶段(吴成平和胡祥云，2007；党如姣，2012；徐佩芬等，2012)。通常情况下孤石分布在特定的地层中(冯涛，2007；谢壮，2010；杨雄飞和袁鸿，2011；杨书江，2014)，在构造条件简单的情况下可以借助钻探查明一部分建构筑物基础区域的地质情况，但大多数情况下大量采用钻孔勘察会大大提高成本和工期，而充分利用工程物探手段，往往能取得较好的效果(李万伦等，2018)。

在工程物探领域中，地震层析成像技术有着不可替代的重要作用(王妙月，1997；赵永贵等，2000；孙宇等，2003；肖宽怀等，2003；郭慧丽和徐佩芬，2011；张凯等，2013)。这种方法借鉴了医学上X射线断面扫描诊断的原理，利用地震波在不同方向投射的波场信息，得出被测区域内岩体波速分布规律，从而对地下介质内部结构进行精细成像(刘畅等，2020)。与常规工程物探方法相比，地震层析具有很高的分辨率，有助于全面细致地对岩体进行质量评价，查明地质异常体(谭显江等，2013；舒森和胡金星，2016；秦晶晶等，2018；张雨飞，2019；何耀京，2019)。

目前，地震层析成像技术已在孤石探测问题上引起重视并得到了研究(李红立等，2016)，但针对地震波在含孤石介质中的传播特征及探测效果的研究较少。本文以正演模拟为基础，构建含孤石的地质模型，分析地震波在含孤石介质中的传播特征，根据地震波的传播特征利用层析成像技术对孤石进行速度成像，并通过厦门某轨道工程实例进行了验证。

1 基本原理

地震层析成像通常使用两个或多个钻孔，经过多点激发、多点接收，在被测区域内形成致密的射线交叉网络(成谷等，2002；胡明顺等，2009；吕子强，2009；孙嘉妤，2019)。这些射线的传播时间包含了被测区域内介质的波速信息。为了求解岩体的波速分布规律，首先需要将被测区域划分成许多规则的矩形网格，网格化方法及观测方式如图1所示。

图1 地震层析成像测区网格化方法

不同接收关系下地震波的到达时间用数学方程可表示为：

(1)

式(1)中，i为射线条数，i=1，2，…，n，n为射线总条数；j为划分成规则矩形单元个数，j=1，2，…，m，m为矩形网格总个数；Ti为第i条射线接收的时间；vj(x,y)为坐标x、y的第j个矩形单元的速度值(待求)，Sj(x,y)为坐标x、y的第j个矩形单元的慢度值；lij为第i条射线通过第j个单元矩形的长度。

当矩形网格划分得足够小时，可以认为每个网格内的介质是均匀的，即波速单一。因此公式(1)可写成离散化形式：

(2)

公式(2)实际上为线性方程组，可表示为：

(3)

其矩阵表达式为：

Tn×1=Ln×mSm×1

(4)

实际计算中，网格很小，穿过每个网格的射线条数很多，但每条射线穿过的总网格数都小于2 m，所以公式(4)中的L是一个大型的稀疏矩阵。为了求解这个稀疏矩阵，首先将被测区域中的v或s赋予一个初值，然后将所得到的投影值残差沿射线方向均匀地反投影过去，并不断对网格中的数值进行校正，一直到满足精度条件为止。上述过程用数学公式可表示为：

(5)

公式(5)可计算出一系列向量S(1)，S(2)，…，当k足够大时，这个向量收敛于所求的值。这个值即为被测区域内波速分布，常用等值线图进行表示。

2 数值模拟

地震波场数值模拟是通过计算机模拟地震波在给定模型中的传播规律，具有方便、灵活成本低等特点。随着计算机技术的发展，地震波场的数值模拟技术已经成为地球物理勘探的重要手段(陈建，2015)。

2.1 单个孤石体的数值模拟

由于地震层析成像的研究对象为直达波走时，为了简化模型提高计算效率，采用声波方程模拟含孤石均匀地层中的地震波场特征。三维正演模型如图2a所示，模型在X、Y、Z方向上的大小分别是50 m、50 m、50 m。其中XY方向为水平地面方向，Z方向为垂直方向。孤石为直径4 m的球形体，球心位于X、Y方向上25 m，Z方向上15 m处。模型在X、Y和Z方向上空间网格大小均为0.5 m。模型的弹性参数如表1。震源采用频率为120 Hz的零相位雷克子波爆炸震源，模型周围添加PML吸收边界，减小模型边界对有效波场的影响。观测系统布置如图2b，采用单孔激发，单孔接收观测系统。震源位置从(5，5，0)到(5，5，30)，炮间距为1 m。检波器位置从(45，45，0)到(45，45，30)，道间距为1 m。

表1 模型介质参数表Table 1 Parameter of model media

图2 三维模型与观测系统

对地震波场三维空间传播特征进行分析。图3为第1炮不同时刻波场快照。其中，左边为三维显示，右边为对应的剖面图。通过分析波场快照可以得出，在20 ms时刻，震源激发产生的直达波在均匀介质中呈球形扩散；在40 ms时刻，直达纵波在球形扩散过程中遇到孤石异常，在接触点发生绕射现象，绕射波相位与原波场相反，且能量减弱。

图3 不同时刻三维地震波场快照与切片

图4为前6炮地震记录。其中，同相轴1为直达波；同相轴2为孤石产生的绕射波。二者速度均为800 m/s。通过对比信号的能量可知，直达波的能量强，孤石绕射波能量弱于直达波，但频率略高于直达波，这与波场快照的分析吻合。前12道检波器可接收到分离的直达波与绕射波，后18道接收到的是直达波与绕射波叠加后的信号。

图4 地震记录

建立好观测系统后对地震信号进行初至时间拾取。初至时间是地震层析成像反演工作中最重要的基础数据，其精度直接影响层析的精度，因此初至拾取工作需要反复多次的进行。图5为第1、5、10、15、25炮每一道信号绘制对应的初至时间曲线。

图5 不同炮地震记录初至时间曲线

将探测区域划分成规则正方形网格，网格边长为0.5 m，对到时数据按公式(5)进行迭代成像，结果如图6。成像结果中，黑色虚线表示数值模拟设置的孤石位置。通过层析成像结果可知，探测区域内存在明显的高速异常，异常范围较大。异常位置在X方向23～27 m，Z方向-13～-16 m处，异常速度为2300 m/s以上，与建立的模型基本吻合，但在高速异常周围存在较大的速度渐变。说明层析成像能很好地反映孤石的尺寸与位置，但成像精度与反演网格大小有关。

图6 单个异常成像结果

2.2 两个孤石体的数值模拟

为研究分析反演网格大小对地震层析成像的精度影响和孤石模型可以区分的纵向距离，设置垂直方向组合的孤石模型。

在单个孤石模型的基础上添加另外一个孤石模型，其它参数保持不变。三维正演模型如图7所示的。孤石1为直径4 m的球形体，球心位于X、Y方向上25 m，Z方向上15 m处。孤石2为直径4 m的球形体，球心位于X、Y方向上25 m，Z方向上21 m处。两个孤石中心在垂直方向上间距为6 m。

图7 三维模型与观测系统

将探测区域划分成规则正方形网格，网格边长为0.2 m，对采集的数据按上述方法进行迭代成像，结果如图8。成像结果中，黑色虚线表示数值模拟设置的孤石位置。通过层析成像结果可知，探测区域内存在两个明显的高速异常，异常区域相对于单个孤石模型的反演结果明显收敛。异常1在X方向23～27 m，Z方向-13～-17 m范围内，异常速度为2300 m/s以上；异常2在X方向23～27 m，Z方向-20～-23 m范围内，异常速度为1700 m/s左右。由于异常1和异常2距离近，波场相互干涉，导致反演结果与实际模型在速度特征与实际大小上存在一定误差。

图8 两个异常成像结果

数值模拟结果表明，地震层析成像能很好地反映孤石异常的位置与速度特征，通过减小反演网格的大小可以明显提升反演结果的精度；当存在相邻孤石异常时，孤石异常的纵向距离大于孤石半径时可以有效区分。

3 工程实例

本次选取实际工程实例来进一步说明地震层析成像的效果。厦门地质条件特殊，某轨道交通施工区间存在较多的孤石，对工程建设造成了较大的困难。根据前期的岩土勘测资料，浅层的岩土主要为碎石填土和花岗岩。碎石素填土主要由黏性土、碎石角砾组成，也有部分块石和残积黏性土，成分不均，纵波速度均值为850 m/s。花岗岩以中等风化花岗岩为主，土质较坚硬，纵波速度均值为2400 m/s；速度的差异作为此次地球物理探查的物性基础。震源采用电火花震源，16道孔中检波器作为信号接收仪器。图9为实测地震记录与拾取的初至。

图9 地震记录初至拾取

将震源与检波器所围成的区域采用0.2 m×0.2 m的网格进行划分。对提取到的时间初至，采用公式(5)中的算法进行迭代，当迭代次数足够多，反演结果基本不再变化时，即为地震层析成像得到的最终结果。图10为联合所有钻孔不同深度的平面反演结果图。对各异常区域解释如下：

图10 异常解释

异常区域1在水平方向上4～16 m处，竖直方向上-8～-14 m处，异常区范围较大。包含2～3个高速异常中心，高速异常中心呈椭球形，其地震波波速2100 m/s以上，解释为多个孤石组成的孤石群；异常区域2在水平方向上24～32 m处，竖直方向上-8～-12 m处，异常区范围较大。主要为1个高速异常中心，高速异常中心呈椭球型，其地震波波速2100m/s以上，解释为单个孤石异常；异常区域3在水平方向上38～43 m处，竖直方向上-10～-14 m处，异常区范围较大。主要为2个高速异常中心，高速异常中心呈椭球形，其地震波波速2100 m/s以上，解释为两个孤石异常。

经现场钻孔施工验证，在划分的异常区域内均有孤石存在，且孤石的大小与规模与反演结果基本吻合，说明地震层析成像法对被测区域内速度差异反映敏感，能有效用于城市地下工程建设的孤石探测，防范孤石地质灾害带来的生命财产安全损失。

4 结论

(1)文章通过数值模拟方法得到了不同孤石异常模型的地震波响应与层析成像结果，说明地震层析成像方法能很好地反映孤石异常的位置与速度特征，通过减小反演网格的大小可以明显提升反演结果的精度。

(2)由于实际孤石异常复杂，当存在相邻孤石异常时，孤石异常的纵向距离大于孤石半径时可以有效区分相邻孤石。

(3)现场工程实践有效查明了厦门某轨道交通施工区间内孤石的分布与规模，说明地震层析成像是一种经济有效的地质异常勘察方法。