地震层析成像概要

地震层析成像
摘要：层析成像方法是一种公认的基于地震数据的有效方法，近20年来，层析成像方法发展迅速。

从原理上讲，层析成像方法可分为两大类，一是基于射线理论走时层析成像，二是基于波动方程的散射层析成像。

本文介绍新的层析成像方法及其技术，包括各向异性介质的2D立体层析成像；时移层析成像的超声数据试验；绕射层析成像的迭代方法：真振幅偏移的本质；用于速度模型构建的下行波折封层析成像和反射层析成像；多尺度波动方程反射层析成像，并在后面展开层析成像方法应用于构造速度模型的分析和实例。

关键字：层析成像；偏移成像；速度模型；克希霍夫偏移。

一、引言
偏移成像在地震勘探和开发过程中，已经成为一种关键的地震数据处理技术。

成像的精度和可靠性依赖于速度模型的准确与否。

速度分析历来都是地震资料处理的基础工作，从均方根速度、层速度以及叠加速度等，贯穿于地震资料处理的方方面面，速度分析方法丰富多样。

迄今，层析成像方法是一种公认的基于地震数据的有效方法，近20年来，层析成像方法发展迅速。

从原理上讲，层析成像方法可分为两大类，一是基于射线理论走时层析成像，二是基于波动方程的散射层析成像。

后一种层析成像很复杂，正处于理论研究阶段。

尽管其实际应用不多，但却是层析成像的发展方向。

走时层析成像比较成熟，有很多的实际应用。

它又可细分为初至走时层析成像和反射走时层析成像。

初至走时层析成像方法简单直观，稳定性较好，主要应用于井间地震以及近地表的速度分析，但是，初至走时层析成像由于只利用初至走时，所以，得到的速度模型比较粗糙，分辨率也较低。

反射层析成像主要应用于地下速度和反射层深度的反演，以及叠前或叠后偏移的速度分析之中。

前者由于速度和深度之间的藕合关系，以及反射波到达时间及其层位难于拾取等，制约了它的广泛应用，但是，这是一种极具价值和潜力的反演方法。

后者则是利用经过叠前或叠后CRI道集中同相轴未被拉平的剩余时差，经过层析成像来修正用于偏移的速度模型。

这种构建速度模型的方法，目前正广泛应用于叠前深度或时间偏移中。

值得关注的还有，地震资料与其他地球物理资料间的联合反演，其反演结果互为验证、相得益彰，为我们提供了更为可靠的反演结果。

二、新的层析成像方法及其技术
1.各向异性介质的2D立体层析成像
立体层析成像是一种利用局部相关同相轴作为输人的斜率层析成像方
法。

本文首次将立体层析成像推广于各向异性介质中，因此不仅使得方程组更为复杂，而且保持反演的稳定性也困难得多。

利用射线扰动理论来计算雅可比矩阵，包括所有数据参数对所有模型参数的导数。

数值试验结果表明:在2DVTI 介质中，各向异性立体层析成像能很好地收敛于4个Thomsen 参数，尽管在反演时每次仅估算一个Thomsen 参数。

后续的工作将在算法中增加诸如co-depthing 的约束条件。

2.时移层析成像的超声数据试验
本文通过在实验室得到的井间超声宽频带的波形数据，来进行时移层析成像方法的有效性研究。

分别利用基于射线理论和散射理论的时移层析成像方法来估算速度差异，关于时移层析成像的方法描述如下：
利用一阶Rytov 近似来模拟波场，仅考虑P 波速度的变化，而忽略S 波以及密度的变化。

那么，在震源位置s r 激发，在检波器位置r r 处接收的Rytov 波场),,(ωs r R r r R 表示为
)),,()(,,(),,(0
0ωωωs r B s r s r R r r P P r r P r r R = （1） 式中，),,(0ωs r r r P 表示参考速度模型)(0r v 时的参考波场),,(ωs r B r r P 是一阶Born 波场，表示为
dV r r G r r P r v r v r r P s r s r v s r ),,(),,()
()(2),,(0302
0ωωωω⎰∆= （2） 式中，),,(ωs r r r G 表示随机参考介质中远场的格林函数;△v(r)表示速度扰动量。

相对于参考波场，从方程(1)推导得到的散射场走时时移),(s r r r t ∆和振幅变化).(/0s r r r A A ∆取一阶近似，表示为
⎰∆∆=∆v t s r dV r K r v r r t )()(),( （3）
⎰∆∆=∆v A s r dV r K r v r r A A )()(),(0 （4）
式中，)(r K t ∆和)(r K A ∆是著名的Frechet 核函数。

对于速度为0v 的均匀参考介质中2D 的波传播，走时时移的Frechet 核函数为
dV x L x v Lz v v v A x L x v L z x K v v v v D
t ⎰∆+∆-∆+-⨯--=00)4
)(sin()()(),(502502ππ （5） 3D 的情况则为
dV x L x v z y L v v v A x L x v L z x K v v v v D
t ⎰∆+∆-∆-+⨯--=00))(sin()()(),(022303π （6） 超声试验的物理模型如图la 所示，图中阴影线部分表示油藏层的注水带，时移层析成像的目的在于确定注水带的范围以及检侧速度变化。

利用同样的模型，来模拟注水前后的情况，因此，采集了两套数据。

在模拟井间超声数据采集时，采用500HZ 的压电传感器作为展源和接收器(均为51个)，共有2601道记录，其频率介于200HZ:和500HZ 之间。

应用散射层析成像，其结果如图lb ，图lb 表示注水前后层析成像结果的差异，可以看到:噪声较大，且低估了注水带的速度，这是因为分别对两个数据进行层析成像，所用的参数不尽相同之故。

为了避免这些问题，直接从两个波场的差异来进行层析成像，其结果如图1c ，有明显的改善。

3.绕射层析成像的迭代方法:真振幅偏移的本质
本文提出了一种非均匀背景的绕射层析成像算法，绕射层析成像可以理解为滤波反向传播，如偏移+建波器。

在非均匀背景介质中显式计算滤波器的计算量太大，可以通过迭代方法来间接计算，而且，还能很容易地将先验信息作为约束条件考虑进算法中，并为单一的频率分量或多频率分量反演的选择提供了极大灵活性，开辟了规则化绕射层析成像的一种新途径。

为方便起见，这里讨论2D 目标函数的情况。

在弱散射假设条件下，正演散射模型的表达式为
r d r r G r P r O r k r P in SC '''''-=⎰),,(),()()(),(20ωωω
图la 为真实的物理模型;b 为注水前后的速度差，c 为采用本文的层析成像获得的注水前后的速度差
可重写为
dr r r G r r G r O r k r O F r r D g s s s g ),,(),,()()()}({),,(20ωωω⎰-==
式中，),,(ωs g r r D 表示震源位于s r 接收点位于g r 处的波场。

这一关系式表明：散射波场是目标函数的线性映射，可认为是算子F 作用于目标函数上。

绕射层析成像是利用记录的波场重构目标函数。

当),,(ωr r G '是震源位于r 在g r 处记录的格林函数),,(ωg r r G ，人射波场),(ωr P in 可表示为震源位于r ，在r 处记录的格林函数),,(ωr r G s s 。

类似于偏移，上式中的格林函数有几种计算方法，我们利用Lambare 等(1996)提出的动力学射线追踪方法。

Born 近似所隐含的物理假设是：目标和背景间的速度差小，与波长相比目标轮廓小，所以有；(l)绕射点的每一个无穷小部分可视为独立的次级震源；
(2)可忽略多绕射。

因此，如果速度模型用网格上的速度分布数值表示，有 )()}({)}({)}({)}({21ω++=+++=g s m r r D r O F r O F r O F r O F 式中，m 表示模型的网格数；m r r r ,,,21 表示所有网格的位置矢量。

以矩阵形式重写
),,()()()}({11111ωg s n m nm ni m r r D d d r O r O f f f f r O F =⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=
这里，n 表示观测n d d d ,,,21 的数量。

一旦以矩阵形式建立起资料(散射波场)和模型(速度扰动)间的关系式，就能迭代地解决反演问题。

数值试验如图2所示，可见其绕射层析成像结果是令人满意的。

4.用于速度模型构建的下行波折封层析成像和反射层析成像
在Trinidad 和Aserbaijan 近海岸地区，由于存在浅层气导致地震信号的信噪比差，因此，构造速度模型是非常困难的任务。

利用反射层析成像通常产生不合适的且分辨率差的速度模型。

在这种情况下，下行波折射层析成像更适合用于构造初始速度模型。

Bell 等(1994，2000)利用折射层析成像处理Mississippi 三角洲、西Texas 和其他地区的地震资料，并取得了浅层速度模型用于确定静校正量。

图2由字母“UT”组成的真实的目标函数(a)，绕射层析成像分别对8个(b)和12个(c)频率重建的图像;右边的3张图表示有先验信息约束时的情况，带有方框的目标函数(b)以及同样的8个(d)和12个(f)频率重建的图像
下行波层析成像:初至时间的折射层析成像是一个迭代过程，包括初始速度模型的建立、走时计算以及走时残差的最小化过程。

下行波层析成像的流程如下:
(1)在炮点/共偏移距道集拾取初至走时;
(2)建立初始的3D近地表速度模型;
(3)计算从炮点至检波点的理论走时;
(4)根据理论走时和观测走时间的差，求解线性方程组，得到速度扰动量;
(5)更新速度模型;
(6)循环步骤(2)至步骤(5)，直至走时残差最小。

值得注意的是，步骤2对于下行波层析成像是至关重要的，因为差的速度模型将导致算法不能收敛或不具有地质意义。

对Trinidad地区的实际资料的下行波层析成像结果如图3所示，同时用于静校正如图4所示。

有意思的是静校正量与另一个独立的浅层气解释的轮廓匹配得很好。

从图5的下行波层析成像结果，以及反射层析成像结果看，两者极为相似。

5.多尺度波动方程反射层析成像
本文利用频率域的公式来讨论了一种波动方程反射层析成像或偏移速度分析方法。

波动方程反射层析成像考虑整个波场，无须识别各个波至，尽管有时需要利用时窗或其他处理来分离单一的散射相位。

该方法包括了零化子(annihilator)对资料的应用，波场零化子是一
种算子，如果波速模型预测到了所有的观测波场，零化子对波场的作用将消失。

从这个意义上说，成功的零化能得到一个可接受的对介质模型估计。

零化准则与传统的不匹配准则相比，具有显著的有点，即能消除地下散射效应，并能有效利用散射波场中固有的信息冗余。

三、层析成像方法应用于构造速度模型
1.折封层析成像的实际应用
折射走时层析成像适用于获取在反射数据处理中静校正所需要的浅地表结构信息。

近年来，折射层析成像也用于获取波形反演和Kirchhoff偏移(对初始速度模型敏感)的初始速度模型.早期的折射层析成像利用射线追踪算法，射线追踪容易且便宜，但只适用于简单的结构而不是复杂的高速度差的模型。

为了得到更加精确的折射层析成像，提出了基于波场的算法和基于Fresnel带的方法。

尽管基于波场的算法对于复杂的高速度差的模型能得到精确的解，但其计算t却比传统的射线追踪方法大得多。

另一方面，基于Fresnel带的方法，与传统的射线追踪方法相比，不需要额外的计算，因此，比基于波场的算法更加有效。

然而，基于Fresnel带的方法在展源和接收器附近以低频率反演有时会失败。

所有的折射走时层析成像算法均有其各自的优点和缺点，在某种程度上大多数方法能给出数值模型的可靠结果。

一个好的折射走时层析成像算法，其必要条件之一是能得到实际资料的可靠速度结构。

对于基于波场的折射走时层析成像，我们选择由Pyun等(2005)和Min 等提出的方法，这两种方法相互类似，都是利用阻尼波场的对数来计算初至走时，区别在于Pyun方法直接计算Frechet导数，而Min方法则是利用反传播算子。

采用在非洲Congo地区采集的资料，共有399个共炮点道集，检波器408个，炮间距和道间距分别为50m和25m，最小偏移距150m。

由上述两种折射走时层析成像获得的速度模型，也可应用于叠前Kirchhoff偏移成像。

反演过程中，利用对0.6283H:的单色阻尼波场取对数来提取模拟数据的初至或相位。

选择低频的波场是为了避免周波跳跃效应。

初始速度模型的速度从 1.5km/s线性增加至4km深度处
的反演结果对比的4.5km/s。

代的反演结果(图6a)和Min方法第
201次及401次迭代的反演结果(图
6b)和(图6c)对比。

图7给出了这两
种方法的均方误差。

图8a和图8b
分别表示采用线性增加速度模型
和Min方法的反演结果作为初始
模型的叠前Kirchhoff偏移结果，
前者未能正确描述盐丘构造左侧
的层位，而后者却能较为精确地揭
示出来，而且还能看到清晰的盐丘
构造。

可见，利用折射走时层析成像获得的速度模型作为Kirchhoff偏移的
初始速度模型，能得到很好的叠前偏移成像结果。

2.利用多尺度层析成像建立深度速度模型
多尺度层析成像的基本思想是以多种尺度的网格来剖分模型(即多个子空间)，以期提高射线的空间彼盖程度，并各自成像，然后对各个子空间的成像结果加权求和，得到下一次迭代的速度模型。

影响层析成像质量的一个关键因素是走时拾取的工作量及其拾取精度，因此，如何进行高精度的走时自动拾取及其对拾取的质量控制，是近年来许多研究者的工作重点。

本文提出了一个新的指标，用于指导走时的自动拾取，即相似性参数，定义如下:
Aa Sc Sm *)0.1(*ωω-+=
式中，Sm （Similarity ）表示相似性；Sc （Semblance coefficient ）表示相似系数；Aa （Average amplitude ）表示平均振幅；ω（0.10≤〈ω）表示权重。

相似系数和平均振幅定义如下：
}]),([/{}]),([{22∑∑∑∑=j i if j M j i if j Sc
M j i if j Aa /]),([∑∑=
式中，),(j i f 表示第i 道第j 个样点，如图9所示。

多尺度层析成像的射线分布见图10所示。

检验的模型采用SEG/EAGE 的盐丘模型，分别采用数值模型的90%真实速度和多尺度层析成像建立的速度模型作为速度输人，进行叠前深度偏移成像，其结果如图11所示。

3.三维不均匀VTI 介质中P 波层析成像速度建立。

近年来，为了能更好地成像，各向异性模型尤其是VTI 模型已经应用于各向异性深度成像中。

Zhou 等(2003)提出了一种层析成像方法用于弱VTI
介质中的速度更新。

本文利用精确的群速度近似式(Fomel ，2004)，给出了一种强各向异性介质中更新速度的层析成像方法。

Fomel 群速度近似式可表示为
])21(cos sin )1(16)(1)()43([)1(41)(122224220
P P p v v E E v n ηφφμηφφηηφ++++++= 对于给定的群角，)(φE 表示椭圆速度，满足
)
21(sin cos )(1222022ηφφφ++=pn p v v E 式中，φ表示群角；p v 表示群速度；pn v 表示P 波时差速度；0p v 表示P 波垂直速度；η表示非椭圆参数。

那么，对于VTI 介质，其层析成像可表述为如下的线性方程组，即 []r v t v t pn pn =⎥⎦
⎤⎢⎣⎡∆∆∂∂∂∂ηη// 式中，r 表示走时残差矢量。

注意，由于更新时差速度和η之后，关于走时的导数矩阵将重新计算，这使得上述方程组的反演是高度非线性的，而且是病态的，原因有： 时差速度和η的数值标量不同，它们对走时的影响不同，以及射线分布的不均匀导致层析反演的收敛慢。

另外，残差的准确拾取对于任何层析成像反演都是至关重要的，因此，本文也提出了一种基于多项式曲线的自动拾取技术，具有高阶校正的功能。

数值模型和实际资料的应用结果表明，利用上述方法建立的各向异性速度模型是精确有效的(图12)。

4.利用层析成像反演进行高频静校正
本文通过在拾取初至走时时增加偏移距的范围，来达到层析成像反演的网格细分，从而提高反演的分辨率，获得合理的近地表速度模型，用于解决高频静校正问题。

其中，利用射线的分布程度，以及将层析成像得到的反演速度与井口速度进行对比，来进行层析成像的质t控制，以保障层析成像结果的可靠性.实际资料处理结果(图13)表明：(1)网格尺度影响刻画模型细节的精度，也决定了解决高频静校正问题的能力，因此，网格尺度是一个关键的参数，应该基于诸如筱盖次数、道间距、接受线距离等因素来进行选取；
(2)对于高密度采集的地震数据，层析成像反演能获得较高精度的速度模型，从而改善高频静校正的精度。

利用层析成像构造的速度模型也能有效应用于
传统探勘项目中的短波长静校正问题。

5.基于有限差分、图论和射线育曲的泥合射线追蛛和振幅计算本文利用
有限差分、图论和射线弯曲，提出了一种新的混合射线追踪方法，来计算走
时和射线路径。

利用有限差分方法来计算随机介质中的走时以避免盲区。

利用基于走时的图论来有效计算锯齿型的射线路径，并利用准射线弯曲方法来修正锯齿型射线路径，最后得到光滑和精确的射线路径。

同时，还提出了一种新的射线几何扩散方程，利用这一方程，能计算基于走时Laplacian 的射线振幅。

这一方法的主要特点是利用了射线弯曲方法，来修正锯齿型的射线路径，其方法如图14所示。

如果点1-k x 和1+k x 固定，调节k x ，那么，寻找点k
x '，使得走时最小。

通过迭代，调节射线路径来获得光滑的射线路径和精确的走时。

为了确定k x '，需要知道mid x 的方向矢量才和扰动量R 。

，根据费
马原理和泰勒展开式，分别表示为
2111111)
)](([→-→+→-→+→-→+→---∙∇-∇=k k k k k k x x x x x x v v n
]2)4()1([)4()
1(222
mid
mid mid mid mid c cv L v n c cv v n c cv R +∇∙++∇∙+-=→→ )11(21,1
111+-→-→++=
-=k k k k v v c x x L 振幅满足如下方程
)(12
22
v A s A t ∂∂-=∇ 利用对数函数偏导数的特征，得到新的射线几何扩散方程为
t v v
A s 22
)][log(∇=∂∂ 数值例子如图15所示，由于模型是水平均匀的，因此，走时和射线路径对称，这样能有些检验算法的可靠性。

混合算法的结果表明：不管是走时、射线路径还是振幅，均具有很好的对称性。

四、结束语
总结各地学者发表的论文等情况可以看出，层析成像的发展趋势可以归纳为以下几点：
(1)走时层析成像技术已经处于较为成熟的应用阶段，包括地面层析成像，其目标主要是为了得到一个较为精确的宏观速度模型，为其他地震资料处理服务，如叠前深度偏移。

(2)考虑各向异性的层析成像技术正不断涌现。

(3)大力发展基于声波方程的波形层析成像方法，其分辨率要其初始模型一般是利用走时层析成像结果。

但是，波形层析成像方法幅敏感，如何提取震源子波等等。

(4)发展基于弹性波方程的全波层析成像技术。

参考文献：
[1]单娜琳,程志平.地震映像方法及应用[J].桂林工学院学报2003,23(1):36~40.
[2] 张三宗,徐义贤.地震记录小波域高阶相关叠加技术[J].地球物理学报, 2006, 49(2): 554~560.
[3] 陈雨红,杨长春,曹齐放,等.几种时频分析方法比较[J].地球物理学进展,2006,21(4):1180~1185.
[4] 黄德济,赵宪生.时频域中主参数剖面的计算及应用[J].物探与化探计算技术,1994,16(3):197~204.
[5] 刘喜武,刘洪,李幼铭,年静波.局域波分解及其在地震信号时频分析中的应用[J].地球物理学进展,2007,22(2):364~375.
[6] 刘喜武,刘洪,李幼铭,等.基于广义S变换研究地震地层特征[J].地球物理学进展,2006,21(2):440~451.李琳，孙宝喜，当议地震映像勘探，黑龙江水利科技，2007年第5期
[7 ]Sengbush R L.Seismic exploration methods[M].Boston:International Human Resources DevelopmentCorp.1983.
[8] 王振东.浅层地震勘探应用技术[M].北京:地质出版社,1988.
[9] 何樵登.地震勘探原理和方法[M].北京:地质出版社,1986.
[10] 陆基孟.地震勘探原理[M].北京:石油工业出版社,1982.
[11] 郑海山,张中杰,田小波. VTI介质中非线性对地震波频率和频宽变化的约束[J].地球物理学报, 2006, 49(3): 885~894.
[12]陈文超,高静怀.基于改进的最佳匹配地震子波的地震资料衰减特性分析[J].地球物理学报, 2007, 50(3): 837~843.
[13] Geoffrey A D. Modern 3-D seismic interpretation. The Leading
Edge,1998,(9).
[14] ZOU De-yun, JIANG Fu-xing. Research of energy storing and gestation mechanism and forecasting of rock burst in the coal and rock mass[J]. Journal of China Coal Society, 2004,29(4):159-163
[15] QI Qing-xin, CHEN Shang-ben, W ANG Huai-xin,et al. Study on the relations among coal bump, Rock-burst and mining tremor with numrical simulation [J]. Chinese Journal of Rock Mechanics and Engineering, 2003,22(11):1852-1857
[16]王维红,崔宝文,刘洪.表面多次波衰减的研究现状与进展[J].地球物理学进展,2007,22(1):156~164.
[17]王大兴,辛可锋,李幼铭,高静怀,伍向阳.地层条件下砂岩含水饱和度对波速及衰减影响的实验研究[J].地球物理学报, 2006, 49(3): 908~914.
[18]殷文,印兴耀.基于MPI的时频分布的改进及应用[J].地球
物理学进展,2005,20(1):165~169.
[19]高静怀,杨森林.利用零偏移VSP资料估计介质品质因子方
法研究[J].地球物理学报, 2007, 50(4): 1198~1209.
[20]邱宁,何展翔,昌彦君.分析研究基于小波分析与谱分析提高重
力异常的分辨能力[J].地球物理学进展,2007,22(1):112~120.。