煤层双巷透射地震层析技术中，在经过不连续性地层时，透射波的速度、振幅、相位、频率等参数会产生变化。透射地震波层析技术以此来探查回采工作面内的地质异常体，如断层、陷落柱、煤厚变薄区、破碎带等。从反演参数区分，地震层析成像分为速度层析成像和衰减层析成像[1-2];从震源上来区分，可以分为主动源与被动源[3-4]。

主动源震波CT，使用可控的爆炸源、锤击等作为震源。主动源震波CT，把震源与检波器安装在工作面两侧巷道的煤壁上，地震波直接穿透工作面内的地质构造。该方法勘探的分辨率高，减少了隐伏小构造探测的遗漏，保障采煤工作面的安全生产[5-7]。主动层析成像技术中震源位置已知(如锤击、爆破等)，射线覆盖几何形状及其覆盖密度人为可控，从而保证了该技术的反演精度。

但是，较大倾角的煤层采用二维透射地震速度反演方法会导致结果出现较大偏差[8]。为克服煤层倾角带来的影响，三维透射波速度层析是较为可行的方法，该方法在煤矿透射波CT反演方面还存在一些尚未克服的困难。JACKSON M J等[9]提出煤层三维层析透射波速度方程欠定的问题。MJ Jackson等利用地面施工的钻孔，进行跨孔速度层析成像，反演煤层的三维构造特征。由于钻孔数量的限制，造成水平方向的数据密度不够，产生速度方程的欠定问题。同样的问题也会出现在双巷道间的三维地震波速度成像中。由于震源和检波器均布置在巷道内，在深度方向上的数据不够，速度方程的欠定问题更加严重。

1 透射地震拟2.5维速度层析反演原理

双巷间的二维透射地震速度层析，未考虑高程对反演结果的影响，而三维透射地震波速度层析又要求垂直方向上有足够多的数据。基于目前的技术条件，提出拟2.5维的透射折射波速度层析的技术方案(图1)。图1(b)为该方案的原理示意，其实质是在二维反演模型的基础上，在垂直方向上设置一个大网格，把平面模型变成一个板状模型。

该解决方案加入了高程参数，使得射线追踪的旅行路径更加准确;同时Z轴方向的网格数为1，使得速度模型本质上还是二维，进而解决了三维速度方程欠定的问题。

1.1 拟2.5维射线追踪

地震波层析成像中正演方法通常可分为两类：一类射线追踪法;另一类是波场数值模拟法。射线追踪方法是一种常见的正演方法，它通过对波场近似的方式达到确定震波射线路径和理论走时的目的，它的计算也非常简明直观，常用来指示复杂地质情况下的传播路径。经典的射线追踪方法有试射法和弯曲法[10]。试射法是最早提出的一种射线追踪方法，该方法基于Snell 定律，通过不断修改震源点处射线的方向，找到能传播到检波器的合适入射角。弯曲法利用费马原理，通过不断迭代修改从震源到接收点的初始路径，直到满足走时最小准则。上述算法实际过程中往往存在一些问题。

为此，有一些学者提出了近似算法以简化计算，提升计算效率，主要理念是采用近似等价的波前来描述地震波场的特征。由费马原理和惠更斯原理，相继提出了有限差分法、走时插值法、最短路径法及波前构建法等[11-14]。这几种方法各有优点，且广泛应用于地震波射线追踪领域[15]。有限差分法由VIDALE[14,16]首先提出，以处理射线追踪法在比较复杂的速度结构体中初至波旅行时求取困难的状况。有限差分法可一次性计算波前到所有节点的最小走时，效率较高。经过诸多学者的不断完善，提出了多种改进方法，有效改善了原始算法的稳定性[17-20]。

这些方法多基于三维相邻网格点间的最少旅行时，如三维有限差分射线追踪法，计算并选取相邻7个网格点的最小旅行时。拟2.5维模型在Z方向上只有一个网格，三维有限差分射线追踪明显不适用。为此，笔者先在XY平面进行二维射线追踪，并把炮点到检波点之间的高程差均匀代入，计算出每个三维网格内的射线长度。

(1)有限差分法首先要求将数值模型进行网格离散，通过已知网格点走时来构造程函方程，求解程函方程可获取网格中其他节点初至走时，具体的计算过程如图2所示。已知点M,N的地震波初至时间t1和t2，计算P点的旅行时。设两个网格内的慢度(OPNM和PQLN)分别为S,S′，网格宽度为L。但此处要考虑到折射波和绕射波出现时，即传播路径为M-N-P,以及由M直接达到P点时。

当满足0≤t2-t1≤时，可由程函方程求取P点到时

(x,y)

(1)

由有限差分可得

，

若S≥S′时，tP=t2+LS′;t2-t1≥时，故可以得到

(2)

根据费马原理，P点的到时应该取这3种情况的最小值，即

,,

(3)

(2)根据上述方法求取相邻点旅行时，并逐步扩展直至完成所有网格节点上的走时计算。

然后依据平面波理论，沿着走时数据的最大梯度方向反推，就可以得到炮点到检波器的最小路径。根据已知的炮检对关系，即可得到二维情况下，第i条射线上第j个网格中射线的长度dij。已知第i条射线的高程差hi，把总高程差平均的加入到经过的M个网格中，即可得到三维情况下的射线长度Dij。

(4)

1.2 拟2.5维透射成像原理

地震层析成像中，透射波的旅行时从数学角度可以认为是慢度沿地震波传播方向上的线积分。实际反演计算时需要用到离散图像的重建技术，其表达式为(LUXBACHER,2008):

→VT=L

(5)

,

(i=1,…,N)

(6)

式中，Ti为地震波的总旅行时;Li为第i条射线的路径长度;V(x,y)为2.5维空间下的地震波传播速度;S(x,y)为对应的慢度;Dij为第 i条射线上第j个网格中射线的长度;N为射线总数;M为网格数。

此时可在探测区域建立一个网格矩阵A(N×M)，可形成

T=AS→S=A-1T

(7)

这时，可以得到一系列像元的地震波走时Ti，速度的重建问题就可转化为求取离散像元内的慢度S。地震波的传播路径只经过很少的像元，因此方程(7)中的系数矩阵A为一个大型的稀疏矩阵，大部分元素为0。对方程(7)的系数矩阵进行求解，即可得到探测区域内的速度分布信息。对于系数矩阵的求取方法研究较多，常见的方法有代数重建法(ART)、联合迭代重建技术(SIRT)、最小二乘法(LSQR)、奇异值分解法(SVD)等方法。联合迭代重建技术由于具有较好地收敛速度和计算效率，求解的稳定性等优点，被笔者所采用。

2 数值模拟实验

地震波微分方程数值模拟方法实质是对计算模型网格化后,通过数值方法求解地震波微分方程在正演模型中的数值解,模拟地震波在正演模型中的传播。笔者应用高阶交错网格有限差分数值模拟技术，用数值模拟得到的t透射波信号进行技术验证。

在进行弹性波传播数值模拟时，利用位移-应力方程表达三维各向同性介质中弹性波方程组，即

(8)

(9)

其中,u,v,w分别为质点位移在直角坐标系中的x,y,z分量;σxx,σyy,σzz为正应力分量;σxy,σxz,σyz为切应力分量;ρ为质点的介质密度;λ,μ为拉梅常数;Vx,Vy,Vz分别为质点在X,Y,Z方向上的振动速度分量，,,

构造不同倾角的单斜煤层三维速度模型。测线和炮点分别布置在两个平行巷道内的煤壁上，巷道的走向为平行于X轴方向。3个模型的巷道A和巷道B之间的高程差分别为0,10,25 m。模型的尺寸分别为100 m×100 m×50 m，煤层厚度为5 m。顶底板岩石纵波速度设为3 000 m/s,煤层速度设为2 000 m/s,巷道声波速度为340 m/s。观测系统采用一发11收的采集方式,检波器个数为11个，震源个数为11个。炮点间距和道间距均为10 m。数值模拟的模型及观测系统如图3所示。

基于上述模型及观测系统，可以得到相应的数值模拟信号。图4为选取的3个模型中，XY坐标一致的3个共炮点道集(x分量)。由图4可知，其首至波起跳点清晰，波形完整。根据正演地震波数据可拾取首至波走时并以此为依据，进行二维透射波速度反演。按照二维震波层析的反演流程，可以得到上述3个模型的二维速度结果(图5)。

图5(a)为水平煤层透射纵波二维速度反演结果，反演模型的速度主要集中在2 900 m/s;说明煤层直达纵波主要沿顶或底板传播，主要反映顶底板的速度特征。

图5(b)和(c)为不同倾角煤层透射纵波二维速度反演结果，图5(b)和(c)中速度反演结果整体较图5(a)低。图5(b)的平均速度约为2 860 m/s,图5(c)的平均速度约为2 750 m/s。该现象表明:煤层倾斜会造成透射波二维速度反演结果偏小;随倾角增大，对于反演结果的影响越明显。同时，反演结果随煤层倾角的变大，会出现“X”型高速异常，即反演模型的对角线出现高异常值。

通过对上述现象分析可得:二维透射波速度反演结果与炮检对之间高程差有紧密关系;煤层直达纵波主要沿顶或底板传播，主要反映顶底板的速度特征;二维速度反演中，炮检间距离为其水平距离，忽略了高程对传播距离的影响。高程差的存在，使得实际传播距离大于其炮检间水平距离，造成反演速度偏小于实际速度。若高程差不变，炮检水平距离越小，对于反演结果的影响就越大;故对角线位置相对于其他位置，高程差影响较小，出现高速“异常”:炮检水平距离越近，计算的速度结果就越低，越偏离真实值。

根据数值模拟结果，对于图5的透射波二维反演结果进行了高程校正(图6)。拟2.5维透射波速度层析反演代入了高程，计算的传播距离更加符合真实值，使得“X”现象得到有效校正且整体速度接近于2 900 m/s，跟水平煤层反演基本一致。该结果说明:拟2.5维速度层析反演方法基本消除了煤层倾角对速度层析反演带来的影响。数值模拟实验验证了拟2.5维速度层析反演理论上的可行性。

3 应用实例

3.1 地质概况

2612工作面位于河南省永城市车集矿区，该工作面目标煤层底板标高在-789～-831 m，施工巷道的高程差较大。工作面呈一宽缓背斜构造，未发现陷落柱、煤层变薄区、煤层急剧变化带、冲刷带等影响工作面回采的地质因素异常情况(图7)。目标煤层平均厚度2.76 m，属中厚煤层。地质条件非常适合本次测试。

3.2 探测技术方案

本次探测采用分布式地震仪和主频60 Hz检波器进行数据采集。结合2612工作面的实际情况，本次工作面探测观测系统如下:

(1)激发点参数。激发点测线布置于2612工作面回风巷中，测线总长约440 m，炮间距均为10 m，总炮数45个。巷道平均高程约为-800 m，呈东高西低的走势。

(2)接收点参数。接收测线布置于2612工作面进风巷中，共布置35个接收点，道间距10 m，接收测线总长约340 m。巷道平均高程约为-820 m，巷道高程变化微小。每一接收点要求布置在煤层中心部位，平行于顶底板，并将传感器与煤层良好耦合。使用250 g乳胶炸药作为震源，炮点孔径φ42 mm，并用黄泥进行封孔。观测系统如图8所示。

对比分析可得，二维透射波速度层析结果呈现较明显的交叉装型高速异常且整体速度偏低。拟2.5维透射波速度层析结果整体速度值相较于二维结果增大且高速异常区集中，交叉型高速异常分布现象减弱。总体而言，拟2.5维透射波速度层析减弱了高程影响，“X”型异常影响减弱，速度异常区与地质构造空间形状和位置更加吻合，取得了较好的应用效果。现场实验验证了拟2.5维透射波速度层析方法的可行性和有效性。

4 结 论

(1)倾斜煤层会造成透射波二维速度反演结果偏小且随倾角增大，对于反演结果的影响越明显。同时，反演结果随煤层倾角的变大，会出现交叉型高速异常，即反演模型的对角线出现高异常值。

(2)数值模拟验证了拟2.5维速度层析反演方法在理论上的可行性;拟2.5维速度层析反演方法基本消除了双巷高程差对透射波速度层析反演带来的影响。

(3)在现场实验中，相较于二维透射波速度层析反演结果，拟2.5维透射波速度层析方法基本消除了双巷高程的影响，交叉行高速速度异常影响减弱，速度异常区与地质构造空间形状和位置更加吻合，验证了拟2.5维透射波速度层析方法在实践中的可行性和有效性。

参考文献(References) :

[1] WESTMAN E C,HARAMY K Y,ROCK A D.Seismic tomography for longwall stress analysis[J].Rock Mechanics Tools and Techniques,1996:397-403.

[2] WESTMAN E C.Use of tomography for inference of stress redistribution in rock[J].IEEE Transactions on Industry Applications,2004,40(5):1413-1417.

[3] LUXBACHER K,WESTMAN E,SWANSON P,et al.Three-dimensional time-lapse velocity tomography of an underground longwall panel[J].International Journal of Rock Mechanics & Mining Sciences,2008,45(4):478-485.

[4] DOU Linming,MU Zonglong,LI Zhenlei,et al.Research progress of monitoring,forecasting,and prevention of rockburst in underground coal mining in China[J].International Journal of Coal Science & Technology,2014,1(3):278-288.

[5] 彭苏萍,凌标灿,刘盛东.综采放顶煤工作面地震CT探测技术应用[J].岩石力学与工程学报,2002,21(12):1786-1790.

PENG Suping,LING Biaocan,LIU Shengdong.Application of seismic tomography in longwall top-coal caving face[J].Chinese Journal of Rock Mechanics and Engineering,2002,21(12):1786-1790.

[6] 潘乐荀,任川,兰鹏波.地震波CT技术在煤矿工作面煤层变薄带探测中的应用[J].安徽理工大学学报(自科版),2016,36(5):72-77.

PAN Lexun,REN Chuan,LAN Pengbo.Application of seismic wave CT technology to detect the thinning zone of coal seam[J].Journal of Anhui University of Science and Technology(Natural Science),2016,36(5):72-77.

[7] SALNIKOV A S,KANAREIKIN B A,DOLGOVA S V,et al.Transmitted-wave seismic tomography for Kuznetsk coal basin:Technology and results[J].Seismic Technology,2012(2):74-88.

[8] 李娜,邢立杰,张建智,等.地震CT反演方法在井下起伏厚煤层中的应用[J].中国煤炭地质,2015,27(12):79-82.

LI Na,XING Lijie,ZHANG Jianzhi,et al.Application of seismic CT inversion in minefield with undulating coal seams[J].Coal Geology of China,2015,27(12):79-82.

[9] JACKSON M J,FRIEDEL M J,TWEETON D R,et al.Three-dimensional imaging of underground mine structures using seismic tomography[J].Digital Library Home,1995:221-230.

[10] PEREYRA V,KELLER H B,LEE W H K.Computational methods for inverse problems in geophysics:inversion of travel time observations[J].Physics of the Earth & Planetary Interiors,1980,21(2-3):120-125.

[11] ASAKAWA E,KAWANAKA T.Seismic ray tracing using linear travel time interpolation[C].Eaeg Meeting,1991.

[12] MOSER T J.Shortest path calculation of seismic rays[J].Geophysics,1991,56(1):59-67.

[13] VINJE V,IVERSEN E,GJOYSTDAL H.Traveltime and amplitude estimation using wavefront construction[J].Geophysics,1993,58(8):1157-1166.

[14] VIDALE J E.Finite and difference calculations of travel times[J].Bulletin of Seismological Society of America,1988,78(6):2062-2076.

[15] 王东鹤,陈祖斌,刘昕,等.地震波射线追踪方法研究综述[J].地球物理学进展,2016(1):344-353.

WANG Donghe,CHEN Zubin,LIU Xin,et al.Review of the seismic ray tracing method[J].Progress in Geophysics,2016,31(1):344-353.

[16] VIDALE J E.Finite-difference calculation of travel time in three dimensions[J].Geophysics,1990,55:521-526.

[17] PODVIN P,LECOMTE I.Finite difference computation of traveltimes in very contrasted velocity models:A massively parallel approach and its associated tools[J].Geophysical Journal International,1991,105(1):271-284.

[18] SUN J G,YANG H,HAN F.A finite-difference scheme for solving the eikonal equation with varying grid spacing[A].77th SEG Technical Program Expanded Abstracts[C].San Antonio:SEG,2007,2120-2124.

[19] 孙章庆,孙建国,韩复兴.三维起伏地表条件下地震波走时计算的不等距迎风差分法[J].地球物理学报,2012,55(7):2441-2449.

SUN Zhangqing,SUN Jianguo,HAN Fuxing.Traveltime computation using the upwind finite difference method with nonuniform grid spacing in a 3D undulating surface condition[J].Chinese Journal of Geophysics,2012,55(7):2441-2449.

[20] 孙章庆,孙建国,韩复兴.针对复杂地形的三种地震波走时算法及对比[J].地球物理学报,2012,55(2):560-568.

SUN Zhangqing,SUN Jianguo,HAN Fuxing.The comparison of three schemes for computing seismic wave traveltimes in complex topographical conditions[J].Chinese Journal of Geophysics,2012,55(2):560-568.