娜，朱培民

娜，朱培民_大人地球物理科学

徐智勇，赵　娜，朱培民

（中国地质大学地球物理与空间信息学院，武汉，430074）

摘　要：煤层内的断层与陷落柱等地质异常体的勘探非常重要。常规井下煤层地震勘探和煤层的成像主要利用的是纵波。如果震源在煤层中，检波器在巷道中时，因为煤层中的速度较围岩的速度低很多，接收到的地震初至波主要是在围岩中传播的折射波。由于这些折射波没有穿过煤层，因此无法利用初至波进行煤层层析成像，而检波器接收到的槽波主要是在煤层中传播，但对煤层有效的槽波信息却在层析成像中没有利用。为了解决煤层及其顶底板的成像问题，可以考虑利用在煤层中传播的槽波和在围岩中传播的初至波进行联合层析成像。基于波前快速推进法计算旅行时，采用联合迭代重建技术进行层析成像。通过理论模型的层析成像试验，证明该方法是可行的。

关键词：槽波，初至波，快速推进法，联合层析成像

1　引　言

在煤矿井下勘探中，煤层中小断层、陷落柱等不连续构造会严重影响采煤生产的正常进行，同时这些不连续构造还常常是诱发矿井安全事故的主要根源，在采煤前或采煤过程中及时准确地查明这些地质构造非常重要^［1］。在采煤前或采煤过程中进行地震勘查时，震源和检波器只能布置在巷道中。由于受井下巷道中地震勘查观测系统的布局限制，将地震层析成像技术应用于煤层成像，探测异常体或断层的位置和大小，是一种自然而然的选择。然而由于煤层相对于上下围岩为低速层，常规层析成像所利用的初至波信息主要是在围岩中传播的折射波，穿过煤层的射线很少，难以得到煤层的速度。而煤层与上下围岩之间的分界面相当于一个良好的反射界面，煤层是一个典型的低速夹层，在物理上构成一个“波导”。当煤层中激发了体波，激发的部分能量由于顶底界面的多次全反射被禁锢在煤层及其邻近的岩石中，不向围岩辐射，在煤槽中相互叠加、相互干涉，形成沿煤层传播的槽波^［2］。有效利用槽波进行煤层的勘查或层析成像是一种必然的选择。

槽波和初至波是两种不同的波类，直接进行联合层析成像是非常困难的。因为槽波主要是纵波和横波干涉产生的一种波动，具有频散特性，没有稳定的传播能量，其传播速度、频散特征与能量变化都受煤层的纵、横波速度和煤层厚度的影响。初至波主要是在高速围岩中传播的纵波，速度较高，而槽波主要是在煤层中传播，速度低，不能用常规的射线追踪方法来同时完成层析成像算法要求的纵波和槽波的射线路径与旅行时的计算。解决的方案是设计同时进行纵波和槽波射线追踪、同时进行层析成像的算法，但是这在目前存在着很大的困难。

笔者提出了一种利用常规射线追踪和层析成像算法进行初至波和槽波联合层析成像的方法，其基本思想是：考虑到槽波只在中间煤层中传播且传播速度低于在围岩中传播的初至波（折射波）的特点，先进行煤层槽波层析成像，然后在此成像基础上，进行围岩顶底板的初至波层析成像，获得整个煤层及围岩的成像结果。

为此，可先建立一个上下围岩速度低于煤层速度的初始模型进行煤层槽波层析成像，结果得到一个中间为高速、上下为低速的模型；然后，将上下的低速层置换为比煤层速度还要高的速度层，这样得到的就是一个上下为高速、中间为低速的模型；第三步，以第二步置换后的模型作为初至波层析成像的初始模型，进行初至波层析成像，得到围岩速度信息。在第三步中，由于射线基本上是沿着高速围岩传播，低速煤层的速度没有修改。在上述三步全部完成后，就得到了煤层、围岩的速度，并成像出其中的异常体。

在本文中，槽波和初至波的射线追踪采用快速推进算法（fast marching method，FMM）。笔者比较了FMM分别采用一阶差分与二阶差分的旅行时计算精度，发现二阶差分算法精度较高，而且采用二阶差分计算时，误差主要集中在炮点周围，不会随离炮点距离的增加而增加。因此，在层析成像时我们采用二阶差分快速推进法来计算旅行时和射线路径。

我们首先介绍了射线追踪所使用的快速推进算法，然后重点讨论了槽波和初至波联合层析成像算法的思想和方法，最后用一个理论模型说明该算法的效果。

2　快速推进算法

2.1　旅行时场的差分方程

地震波从炮点激发逐步向外传播，其传播旅行时逐渐增大，且不可能出现负值。快速推进法是由Sethian等^［3，4］基于求解哈密顿方程提出的算法，被用于程函方程来计算初至旅行时、反射旅行时、多次反射和透射旅行时等^［5，6］。FMM基于迎风有限差分格式求解程函方程，采用窄带延拓重建旅行时波前，向外扩展求解程函方程的方法。由于建立了窄带环绕波前，可以大大减少无用的计算工作量，提高计算效率。

计算旅行时场的程函方程为：

式中：（x，y，z）表示空间位置坐标；T和S分别是旅行时函数和介质的慢度分布函数；▽T（x，y，z）表示旅行时的梯度。计算地震波旅行时的时侯，用上风格式表示式（1）的梯度函数，可得到如下的表达式：

式中：D^-x、D^＋x分别是x方向的向前、向后差分算子；D^-y、D^＋y分别是y方向的向前、向后差分算子；D^-z、D^＋z分别是z方向的向前、向后差分算子；max和min表示取两个数中最大值和最小值的运算符号。

x方向的向前、向后差分算子一阶差分格式定义为：

x方向的向前、向后差分算子二阶差分格式为：

y和z方向的差分格式与x方向的差分格式相同，仅需将上式中的x改成y或z。

2.2　计算最小旅行时的步骤

图1是快速推进法的射线追踪计算过程的示意图。算法的基本思想是将全区的点分为下风点、窄带点和上风点。上风区中点的旅行时已经确定，在以后的计算中不改变；窄带区的点的旅行时虽然已经计算出来，但不一定是最小旅行时，是可以改变的；下风区的点的旅行时还没有计算，是未知的。

射线追踪的具体步骤是：①计算开始时，先确定一个震源点O，并将旅行时设置为0；②在窄带区（灰色点组成的区域）搜索旅行时最小的点（例如，点B）；③将旅行时最小的点从窄带区移出并移入到上风区（黑色点组成的区域）；④将旅行时最小点B周围下风区（空心点组成的区域）的点移入窄带区；⑤根据方程式（2）重新计算该点（点B）邻近非上风区点的旅行时，如果比原值大，则保留原值，如果比原值小，用新值取代之；⑥重复②到⑤，直到全部点进入上风区。

图1　快速推进法示意图

（a）首先确定震源O点的值；（b）将震源O点周围几个点移入窄带区（灰色的点A、B、C、D），用迎风法计算窄带里面点的可能值；（c）从窄带区中找出值最小的点（假设是B），将B点移入上风区；（d）将B点周围的点（E，F，G）移入窄带区，并计算其可能值；（e）从窄带区中找出值最小的点（假设是A），移入上风区；（f）将A点周围的点（E，H，I）移入窄带区，并用A点的值计算其可能值

由于该算法总是在上风区的邻域点计算旅行时，不可能产生比任何一个上风区中点旅行时还小的值，并且总可以选择窄带区中旅行时最小的点，调整周围点的值，不断向前推进。

该算法的计算量为O（NlogN_B），其中N为网格点数，N_B为窄带区的点数。由于窄带区的点数在计算的过程中是不断变化的，因此该算法的计算量上限为O（NlogN）^［7］。

2.3　计算精度分析

在用快速推进法计算时，影响旅行时计算精度的因素主要有差分阶数、网格间距等。下面以速度均匀分布的三维模型为例讨论计算误差。图2（a）是一个速度均匀介质模型，x、y、z三个方向的网格数均为80，3个方向上的网格间距均为1m，速度大小为2 000m/s，震源点位置坐标为（40m，40m，0m）；图2（b）是理论旅行时；图2（c）是采用一阶差分计算得到的旅行时；图2（d）是采用二阶差分计算得到的旅行时。

图2　快速推进法计算旅行时试验

（a）三维速度模型，速度大小为2 000m/s；（b）速度模型的理论旅行时；（c）采用一阶差分计算得到的旅行时；（d）采用二阶差分计算得到的旅行时

比较一阶差分和二阶差分计算得到的旅行时，可以发现二阶差分的精度较高。图3（a）、（b）分别是一阶差分和二阶差分计算得到的旅行时误差（旅行时误差＝计算旅行时—理论旅行时）分布图，图3（c）、（d）分别是z＝0平面的一阶差分和二阶差分旅行时误差分布图。其中，一阶差分最大误差为0.000 985s，二阶差分的最大误差为0.000 311s。从z＝0平面的旅行时误差分布图可以看出：①一阶差分格式，其较大的误差主要分布在远离炮点的位置；②二阶差分格式，其较大的误差主要分布在炮点周围，炮点周围的误差大小决定了整个区域的误差大小，而且误差并没有随着离炮点距离的增大而变大。

图3　一阶差分和二阶差分旅行时误差

（a）一阶差分的旅行时误差；（b）二阶差分的旅行时误差；（c）z＝0平面，一阶差分的旅行时误差；（d）z＝0平面，二阶差分的旅行时误差

3　三维层析成像

3.1　层析成像方程离散化

由于层析成像是一个非线性问题，反演时需线性化。设慢度场为S（r），地震波沿射线路径R［s（r）］传播的时间t可表示为^［8］：

上式离散后变为线性方程组：

式中：A_i是第i（i＝1，2，…，I）条射线穿过各个网格的长度所形成的向量；S是慢度向量，S＝（S₁，S₂，…，S_n），实际计算时对所有I条射线的数据方程进行组合，计算得到慢度向量。

上述方程组的核矩阵A＝［A₁，A₂，…，A_I］，具有超定、稀疏、病态的特点。我们采用联合迭代重建技术（simultaneous iterative reconstruction techniques，SIRT）算法求解方程组（6）。

3.2　联合迭代重建技术（SIRT）

SIRT通过被重建图像中某一个像素的所有射线修正值来确定这个像素的平均修正值，这样的平均值可以消除某些干扰因素。从数学上考虑，我们不一定要求某个方程的残差消失，而是要使残差呈递减趋势。从物理上说，是对穿过图像重建中某一个像素的射线取平均修正量，即

式中：i表示第i条射线，j表示第j个网格（像素）；A为核矩阵，A_ij表示第i条射线所通过第j个网格的射线长度；r_i为第i条射线的旅行走时残差；q为迭代次数；N_j表示核矩阵A的第j列中非零元素的个数，即通过第j个网格的射线条数。在效果上，利用上式可以消除某些干扰和随机测量误差；在计算上，该方法比较稳定_［9］。

3.3　SIRT实现步骤

（1）选定一组初值s^（0）_j，j＝1，2，…，J；

（2）计算估计值，i＝1，2，…，I；

（3）计算、实际值与估算值的差r_i＝t_i--t_i；

（4）求出第j个像素内平均修正值Δs^（q）_j，假设在像素内共有Nj条射线通过，则有：

（5）用平均修正值Δs^（q）对第j个像素的慢度值s进行修正：

_jj

s^（q＋1）＝s^（q）＋Δs^（q）

_jj

j

同时要加约束条件：

s≤s^（q＋1）≤s

_minjmax

当s_j＜s_min，令s_j＝s_min；s_j＞s_max，令s_j＝s_max；

（6）设ε≥0，当

｜s^（q）-s^（q＋1）｜≤ε，j＝1，2，…，J

_jj

时停止迭代，否则进入下一轮迭代，重复（2）～（6）的步骤。

3.4　层析成像模型试验

设有一个三层水平层状介质模型［图4（a）］，模型网格数为20×20×20，x、y、z方向网格间距都为10m，三层速度分别是4 000m/s、3 000m/s、2 000m/s，为了达到最佳射线覆盖观测效果，炮点分别布置在模型的x＝0m、y＝0m、z＝0m三个面上，各面有66炮激发，检波器布置在模型的x＝200m、y＝200m、z＝200m3个面上，各面有66道接收，共设计有198炮激发、198道接收。初始速度模型采用均匀介质模型，速度为1 000m/s，经过400次迭代得到如图4（b）所示的层析成像结果，从图中可以看出层析成像结果的分辨率较高，这说明笔者设计的算法是可行的。

图4　层析成像模型试验

（a）理论模型；（b）层析成像结果

4　槽波和初至波联合层析成像

煤层相对于上下围岩为低速层。根据费马原理，当炮点放置在巷道或煤层中，接收点放置于巷道中，进行井下地震勘探时，射线基本上会绕过低速煤层，以折射波在高速的围岩中传播。例如，一个如图5所示的地震地质模型，上下围岩速度为4 000m/s，中间煤层速度为2 000m/s，图中左边红色点代表炮点，右边绿色点代表检波点，炮点和检波点都设置在煤层中，射线穿过低速煤层进入高速岩层中传播，中间煤层大部分区域中没有射线通过。

根据槽波只在低速层中传播和初至波只经过高速层的特点，设计以下槽波和初至波联合层析成像流程（图6）。

图5　煤层地质模型射线轨迹图

第一步，建立一个比理论煤层速度低的均匀速度模型，FMM射线追踪得到槽波旅行时，利用观测记录中的槽波信息，层析成像得到煤层的槽波速度。

第二步，利用槽波速度与纵波速度之间的关系，将槽波速度换算为纵波速度，将模型中纵波速度低于煤层速度的模型值置换为高于煤层速度的值。

第三步，以第二步置换后的速度模型作为初至波层析成像的初始模型，利用观测记录中的初至波信息，层析成像得到围岩的速度。

图6　槽波和初至波联合层析成像流程图

由于在第一步层析成像中得到的是煤层的速度，在第三步层析成像中，只修改了围岩的速度，因此，最后得到的速度就是整个模型的速度。

5　槽波和初至波联合层析成像模型试验

5.1　理论模型

为了验证槽波和初至波联合层析成像的可行性，笔者设计了一个三层模型（围岩-煤层-围岩），上下围岩的速度均为4 000m/s，中间煤层的速度为2 000m/s，划分成20×20×20个网格，网格间距0.4m，模型长、宽、高各为8m，煤厚2m，煤层中有一高速异常体，异常体大小为4×4×4个网格，即1.6m×1.6m×1.6m，异常体速度为4 000m/s，如图7所示。炮点和检波器均设置在煤层的4个表面上（图8），采用四面放炮、四面接收，每面各有36个震源、36个检波器。共144炮激发，144道接收。

图7　联合层析成像理论模型

（a）理论模型；（b）模型切片

图8　理论观测系统

为了得到槽波理论观测数据，通过将围岩速度值替换为比煤层速度低的值，以此射线追踪得到的旅行时作为槽波观测数据。设计的模型如图9所示。

5.2　槽波层析成像

槽波层析成像采用的初始模型为均匀介质模型，速度为1 000m/s［图10（a）］，按照槽波与初至波联合层析成像流程的第一步，经过400次迭代，得到的层析结果如图10（b）所示，中间的煤层速度与所设计模型的煤层速度非常接近，但是反演的煤层厚度比真实模型要厚。

图9　槽波层析成像理论模型

（a）理论速度模型；（b）模型切片

图10　槽波层析成像

（a）槽波层析成像初始模型；（b）槽波层析成像结果

5.3　速度置换和初至波层析成像

速度置换：设定一个煤层速度的上下限范围（1 800～2 200m/s），以槽波层析成像结果为基础，将低于煤层速度下限的速度值替换为一个高于煤层速度上限的速度值（3 000m/s），中间煤层的槽波速度乘以一个比例因子（在文中所用比例因子是1.2），将之换算为煤层纵波速度，以此作为速度置换后的初至波层析成像的初始模型［图11（a）］。

初至波层析成像：按照槽波与初至波联合层析成像流程的第三步，经过400次迭代，得到的层析成像结果如图11（b）所示。可以看出中间煤层的速度和厚度与所设计的煤层的速度和厚度非常接近，煤层围岩分界面比较明显。

图11　初至波层析成像

（a）初至波层析成像初始模型；（b）初至波层析成像结果

5.4　成像结果分析

图12（a）为图10（b）槽波层析成像结果的切片图，图12（b）为图11（b）初至波层析成像结果的切片图，这也是联合层析成像的最终结果。从图12（b）中可以清晰地看出煤层、煤岩分界面以及高速异常体，说明利用该方法对陷落柱异常体等的成像是基本可行的。

从图10至图12中，可以明显看出，煤层、煤层中间的高速异常体和紧挨着煤层的围岩得到了有效的成像，而离煤层较远的模型速度与真实模型速度不一致。原因是观测系统设计在煤层中射线从震源到检波器只通过少量的离煤层较近的岩层，所以有射线通过的岩层其速度得到了较好的成像，而较远的没有射线通过的岩层则保持给定的初始速度不变。

从图10（b）、图12（a）中可以看出，第一步利用槽波层析成像，已基本上把煤层、异常体、煤岩分界面成像出来了，但是煤岩分界面与真实模型并不一致。通过第三步层析成像，煤岩分界面回复到了正确位置。

图12　槽波和初至波层析成像结果

（a）槽波层析成像结果的切片；（b）纵波初至波层析成像结果的切片

6　结论和问题

针对煤层含陷落柱特殊地质情况下的成像问题，我们提出的分步槽波和初至波联合层析成像算法，经过数值模拟试验，验证了该方法的有效性，但还需要实际数据的验证。

在此文中射线追踪采用快速推进算法，比较了在相同网格间距情况下，一阶差分与二阶差分旅行时计算结果的精度，通过对比，得出二阶差分旅行时计算误差明显低于一阶差分误差，且二阶差分计算时，误差主要集中在炮点周围，不随着离炮点距离的增加而增加，因此，在层析成像中，我们采用二阶差分，有效减少了旅行时计算误差，提高了计算精度。

目前这种算法还不够成熟，槽波和初至波的射线追踪是分开进行的，而且槽波的射线追踪还没有考虑频散和模式的变化。下一步的研究工作，重点将放在槽波纵波的射线追踪上，通过一次计算，同时得到槽波和初至波信息，从而为煤层及异常体提供更精确的成像。特别是，由于槽波射线追踪需要横波速度，这种考虑了槽波层析成像的方法还可以同时得到横波的信息，有利于含水层、含气层和陷落柱的探测。

参考文献

［1］张守恩，姜克富.地震槽波方法试验与研究［J］.地球物理学报，1983，26（2）：198-203.

［2］刘天放，潘冬明，李德春，等.槽波地震勘探［M］.北京：中国矿业大学出版社，1994.

［3］Sethian J A.Fast Marching Methods［J］.SIAM Review，1999，41（2）：199-235.

［4］Kimmel R，Sethian J A.Computing Geodesic Paths on Manifolds［J］.Proceedings of National Academy of Sciences，1998，95（15）：8 431-8 435.

［5］Popovici A M，Sethian J A.3-D imaging using higher order fast marching traveltimes［J］.Geophysics，2002，67（2）：604-609.

［6］Rawlinson N，Sambridge M.The fast marching method：an effective tool for tomographic imaging and tracking multiple phase in complex layered media［J］.Exploration Geophysics，2005，36（4）：341-350.

［7］单联瑜，刘连升.三维高阶快速步进波前重构旅行时计算精度分析［J］.石油地球物理勘探，2007，42（4）：413-417.

［8］王守东.三维表层有限差分层析成像［J］.勘探地球物理进展，2003，26（4）：252-254.

［9］吴律.层析基础及其在井间地震中的应用［M］.北京：石油工业出版社，1997.