地震成像现状存在问题及发展趋势分解

地震成像技术的发展现状存在问题及发展趋势
（杜炳毅地球探测与信息技术S1*******）
随着地震勘探难度的逐渐的增加和油气藏复杂性的增加，油气勘探开发对地震勘探精度的要求越来越高。

为了实现高精度的地震资料在油气勘探中的应用，近年来地震方法和技术重点发展了两个方向：一是地震成像技术，二是开发地震技术。

地震成像技术发展现状
地震成像是现代地震勘探数据处理中的重要组成部分，分为叠加成像和偏移成像。

随着油气勘探难度的增加，地震成像技术得到了迅速的发展，并且成为高精度地震勘探数据处理的关键技术。

地震偏移成像可以分为地震叠后偏移方法和地震叠前偏移方法。

叠后偏移是在共中心点叠加数据上进行零炮检距偏移，主要有叠后时间偏移和叠后深度偏移，叠后时间偏移主要包括射线偏移和波动方程偏移。

而叠后深度偏移可以有效的结果构造不太复杂，横向速度变化比较大的地质体的地震成像问题，并且能够提高地震成像的计算效率，常用的叠后深度偏移有Kirchhoff积分法，分步傅里叶法，有限差分法以及逆时偏（RTM）法。

叠前偏移是把共炮点道集记录或者共偏移距道集记录中的反射波归位到产生他的反射界面上，并使绕射波收敛到产生它的的绕射点上。

也分为叠前时间偏移和叠前深度偏移。

叠前时间偏移是基于绕射
叠加或者Claerbout发射波成像原则，是一种成像射线，能够解决叠后时间偏移存在的问题，叠前时间偏移的方法主要有Kirchhoff积分法叠前时间偏移，波动方程法叠前时间偏移（包括平面波分解法叠前时间偏移和F-K域法叠前时间偏移）；叠前深度偏移方法可以分为两类：第一类是基于射线理论的叠前深度偏移方法，另一类是基于波动方程理论的叠前深度偏移方法。

射线法叠前深度偏移方法主要有Kirchhoff积分法叠前深度偏移，高斯波束叠前深度偏移；波动方程叠前深度偏移主要有F-X域有限差分叠前深度偏移，SSF法波动方程叠前深度偏移，Fourier有限差分（FFD）法波动方程叠前深度偏移，广义屏近似波动方程叠前深度偏移，基于双平方根方程的波动方程的叠前深度偏移，基于波动方程的真振幅偏移，逆时叠前深度偏移。

地震偏移是一种将地震信息进行重排的反演运算,以便使地震波能量归位到其空间的真实位置,获取地下真实构造图像。

除了深度域构造成像外,地震偏移还为其它特殊处理提供振幅、相位等信息,用于速度估计和属性分析,建立在波动方程基础上的地震偏移成像技术代表了地震处理的极致。

地震偏移最初是在水平迭加基础上进行的,目的是使倾斜界面共深度映像聚焦,使绕射波归位,即将能量还原到它们正确位置上.早期人工偏移是按照偏移空间的时距关系作图；若将共深度点剖面看作一系列绕射点组成的源反射,可用计算机实现对这些绕射点的偏移,即建立在射线理论基础上的绕射扫描迭加方法以及后来的Kirchhoff 偏移.20 世纪70 年代初美国斯坦福大学以J. F. Claerbout 为首的
SEP 研究小组第一个对标量波动方程提出了有限差分近似解法, 实现了地震偏移.此后建立在波动方程基础上的地震偏移成像方法如有限差分法、Kirchhoff 积分法、F - K 方法及其各种变形等方法广泛应用爆炸反射面模型(Loewenthal et al. ,1976)为波动方程偏移成像条件的建立奠定了理论基础.由于波动方程描述地震波地下传播规律,因而波动方程偏移一方面可以解决复杂介质条件下成像问题,另一方面保持了波场的动力学特征。

地震偏移各种方法最初是作为时间偏移方法出现的, 目的是满足二维时间叠加剖面成图需要,后来为满足横向变速情况下成像精度需要, 发展了深度偏移方法(Hubral 1977,Larner 1981)。

近20 年来, 偏移方法又发展到了三维和叠前偏移,三维叠前深度偏移代表地震偏移的发展水平. 当今各种各样的偏移技术方法极为丰富, 如时间偏移、深度偏移、二维、三维、迭前、迭后; 如使用共炮集、共方位角道集、面炮方法实现等等.偏移算法也多种多样.实际应用中根据具体情况和要求选取相适应的方法. 考虑如此复杂而庞大的内容, 本文只对波动方程地震偏移中核心算法。

波场延拓和成像的现状与进展进行阐述, 并对其发展趋势进行展望.
波动方程地震偏移成像方法的研究现状
波动方程地震偏移成像的各种方法都是建立在波场反向外推基础上,按照算法实现的原理可以分为两大类: 基于射线理论的偏移方法和基于波场延拓的偏移方法[10]。

基于射线的Kirchhoff 积分类方法,
依靠射线追踪获得成像所需的旅行时,不受反射界面顷角限制, 计算效率高,灵活,但在复杂地质条件下,多值走时使射线追踪难于获得正确旅行时,导致成像效果较差.此外,基于射线的方法缺少动力学信息如振幅等.基于波场延拓的方法,如有限差分法、F-K 方法等,物理概念清晰,自然解决了多值走时问题,能够更为精确成像.这类方法包括由双程波动方程导出的逆时偏移,由单程波动方程导出的各种方法.实质上,建立在波动方程基础上的Kirchhoff 积分方法与波场外推的F-K方法、有限差分方法数理基础相同。

已证明在常速介质中Kirchhoff 积分方法与F-K方法的波场外推公式完全等价,而有限差分方法使用波动方程的各种近似,其波场外推公式除相位精度外,形式与前两者基本相同,是它们的近似式.然而,由于波场反向外推的实现算法不同,导致它们各自不同的特点.
1. 1 Kirchhoff 积分类偏移方法
Kirchhoff 积分法波动方程偏移建立在波动方程Kirchhoff 积分解[ 11] 的基础上, 把Kirchhoff 积分中的格林函数用它的高频渐进解( 即射线理论解) 来代替.其基本过程包括从震源和接收点同时向成像点进行射线追踪,然后按照相应的走时从地震记录中拾取子波并进行叠加,如果对所有的路径计算出的走时都正确,对所有记录数据的叠加结果会在某些部位产生相对较大的值, 这些值给出地下界面( 即反射体位置)。

Schneider( 1978)建立了Kirchhoff积分偏移的波动方程理论基础,Bleistein( 1987) 将Kirchhoff积分方法拓展
到求解反射系数,进一步推进到偏移后的参数估计.
虽然Kirchhoff 积分公式是严格的波动方程分解, 但它的实现是利用波动方程的零阶高频渐近近似(射线方程),这种近似只有在t 时刻的圆频率 大时才合理,因而源点或接收点的几个波长以内的绕射点不能正确成像,存在焦散区.其次,在复杂介质中,由于速度不均匀和高频近似,绕射点与源点和接收点之间的传播距离要远大于几个波长,这种大距离反向外推波场,就存在多重路径问题,造成旅行时求取困难,近年来人们不断改进旅行时求取方法Audebert et al ( 1997) 对这些方法进行了总结和对比。

由Kirchhoff 积分法成像时对绕射面扫描到的数据,没有考虑频率成份, 绕射面较陡部分在它扫过未偏移数据平坦部分时,将对地震子波进行重采样,造成算子假频.Gray (1992)和Lumleyet al. (1994) 提出减少绕射面段陡部分扫描到子波频率成份,克服了这个问题。

归纳起来, Kirchhoff积分法偏移成像局限为:(1) 分辨率随着
深度的增加而变差,从而导致对深部结构分辨率降低,这一现象源于
利用射线解近似格林函数时对菲涅尔带的影响;(2)成像信息中缺乏
正确的振幅信息,这一现象源于射线近似在复杂介质中存在焦散、多重路径和干涉等问题.
为了保持灵活性，同时提高Kirchhoff积分偏移的精度,Hill( 1990, 2001)提出了高斯束偏移方法。

该方法将源点和接收点的波场局部分解为“束”导回地下.几个束可能来自不同地面位置, 且不同束指向不同的初始传播方向.每一束由各自射线管引导,独立
于其他束传播.射线管可以重叠,这样能量可以通过多个路径在像点位置与源点和接收点之间旅行,解决了多路径问题Bevc(1997)也提出一种解决多路径方法.该方法先向下应用标准的非递归Kirchhoff偏移方法将记录偏移到地表以下多个波长深度处.在这个深度范围内认为多值走时不很严重,Kirchhoff积分方法可以较精确进行.接着在该深度处采用Kirchhoff基准面方法(Berryhill,1984)计算一个向下延拓的波场.然后将这个波场用于下一个有限范围深度内Kirchhoff 偏移.经过这些Kirchhoff偏移和向下延拓的结合,完成偏移.这种方法通过分级办法处理多路径问题, 但只对二维实用。

1.2 基于波场延拓的波动方程偏移方法
Kirchhoff 积分偏移采用地表所有记录数据单一整体空间褶积计
算每一个点的像,并且地表位置与像点位置之间通常只采用一个路径, 基于波场外推的波动方程偏移方法递推地从前一个深度Z的波场计
算深度Z Z
+∆的波场,自然考虑到每个深度可能的绕射点与每个源点
或接收点的多重路径,因而可以在较为复杂介质条件精确成像。

随着
偏移技术的发展,波场外推偏移成像从最早单一的T - X域实现,发展
为诸如T-K 、F-K 、F-X 、P
τ-域、小波变换域等多种域实现偏移方法.偏移方法包括由双程波方程导出偏移方法,如T-X 域有限差分法, 逆时偏移方法；以及基于单程波方程近似解的各种方法.
1.2.1 T-X域或波数域有限差分方法( F-D)
有限差分法波动方程偏移是最早提出的一种波场延拓波动方程偏移方法。

这种方法直接对T-X 波动方程进行坐标变换并略去二阶导数项, 得到变换后的简化的波动方程,然后再利用有限差分方法求解波动方程, 进行成像.当地面已知波场()
P x t向地面下延拓至反射
,0,
点Q( x’, z’)时,波场函数P( x’, z’, t’)的旅行时t#为波从震源到反射点Q ( x’,z’)的下行传播时间[23].F - X域有限差分偏移方法存在一些固有的困难[ 18] :①由于空间离散化造成的数值频散, 导致不同频率的波以不同的速度传播,从而造成成像的误差和人为假象；②由于矩形网格划分导致的三维数值各向异性,造成沿不同方向波传播速度不同.
用差分方程近似替代偏微分方程, 可分为隐格式和显格式两大类. 隐格式求解困难, 但精度高,对倾角较大情况偏移效果好. 隐格式有限差分偏移自Claerbout(1971) 提出15°有限差分偏移以来得到充分发展,马在田( 1982)和张关泉( 1986)分别提出高阶方程降阶方法, 有效地解决了隐式有限差分方法难以高角度成像问题, 但是没有考虑到降低计算量的问题[ 35, 37] .李志明的双线性变换和三维隐式有限差分多方向分裂算法中各项异性的补偿是隐式有限差分偏移方法的出色工作[ 27] .
显格求解较容易,快速,但精度较低,对大倾角适应性差.显式有限差分法由Berkout提出[ 26],Holberg和Hale( 1999)在二维情况下有效解决了该算法的稳定性[ 2],Blacquiere进一步将其推广到三维[2] .
对近似方程中x 做Fourier 变换到波数域,可进一步得到近似方程,再利用有限差分求解.
基于波场外推的波动方程偏移方法(尤其是F - D方法)与Kirchhoff 积分偏移方法在偏移孔径处理上存在显著差异[ 2]。

Kirchhoff 偏移方法通过绕射曲线对输入道叠加或将输入采样散布到输出道.波场外推偏移从记录地表将整个波场向下延拓,如果需要较大孔径,这个孔径必须作为一部分包含在整个计算内,换句话说即使记录地表这些道位置处没有能量,在所有深度偏移运算也必须包含所有的输出道.这就导致额外巨大计算量,在三维叠前深度偏移中, 难于经济适用。

1. 2. 2 逆时偏移( Reverse time migration)
逆时偏移也是一种应用有限差分求解波动方程实现波场延拓的方法,它不是深度域外推,而是进行时间外推,求解双程波声波或弹性波方程, 并且允许波向各个方向传播.这种方法没有倾角限制, 精度较高. 它的计算方法正好与地震正演模拟的计算顺序相反,以最大时间开始向最小时间计算。

Baysal et al( 1983)和Mcmechan( 1983) 给出这种方法的详细内容并描述了高精度全倾角成像能力[ 24, 25]。

之后许多学者提出了不同的差分格式,不同计算方法实现逆时偏移方法, 如双线性变换逆时偏移[ 27] ( 李志明, 1991) 。

逆时偏移法完全遵循全波波动方程,是最为精确的方法,但同有限差分模拟一样,存在稳定性和数值频散问题,同时计算量非常大。

1.2.3 基于单程波方程的偏移方法
从理论上讲,波场延拓应该用波动方程的边值问题来解决,但由于地面上的观测范围有限,这种边值问题是不稳定的,为此,地球物理学家提出用单程波方程来做波场延拓.通常的叠后偏移(也就是零炮检距剖面偏移)就要用到单程波方程波场延拓.叠前深度偏移的共炮集偏移,是通过炮点波场和接收点波场向下延拓实现的。

在叠后偏移中的成像运算是取零时刻的波场值。

炮集叠前深度偏移的成像运算是炮点波场与接收点波场做互相关。

基于单程波方程的方法有w-K域stolt法、w-K域相移法(phase shift )和相移加插值法( phase shift plus interpolation),w-X域有限差分算法(finite difference)、w-K域与w-X域交替的裂步Fourier算法(split step Fourier),w- K域与w-X域Fourier有限差分算法(Fourier finite difference),K域与w-X域广义相位(generalized screen propagator)方法.
(1) Stolt 偏移算法[28]
Stolt( 1978) 提出的w- K 域偏算法.该算法不采用波场延拓途径而是通过已知叠加数据p( x , t , z = 0)的谱p(kx ,w,z = 0) 一步求出偏移剖面p ( x , t= 0, z ) 的谱p(( kx , t= 0, kz ),再经过二维Fourier反变换形成偏移剖面.该算法运算极快,效率高,无倾角限制和频散.该算法要求每次外推时全部使用一个平均速度,不能适应横向变速, 是一个时间偏移方法。

(2) 相移法( PS) 和相加插值算法( PSPI)[29 -31]
Gazdag( 1978)提出相位移方法.( 3)式是其基本公式.该方法无倾角限制,无频散,但由于在w-kx 域中延拓,kx与V( x )关系无法确定, 因而不能适应横向变速情况.Gazdag(1984)提出相移加插值算法克服这个问题,在每一个深度用多个不变的参考速度计算几个外推波场,再变换到w-x域,真正波场根据V( x )与参数速度的关系通过线性插值得到.该方法纯粹是一种数学手段,无物理本质意义,计算量较大,且不能完全改善对横向变速的适应.
(3) 裂步Fourier 方法( SSF)[32,33]
针对相移法不能适应横向变速情况,Stoffa et al(1990)提出一种分裂步Fourier方法.该方法基于小扰的理论,将速度场分为背景速度和扰动项之和,推导出波场延拓公式,交替在w-kx域和w-x域进行, 实现对横向速度变化的处理.该方法无倾角限制和频散现象.
对强烈的横向变速情况,扰动理论不成立.为解决这一问题,类似于相移加插值思想,Kessinger( 1992)将多参考速度对数(MRVL)引入SSF方法,对两参考速度之间重叠区进行插值,但增加了较大计算量. (4) w-X域有限差分方法(F-D)[34—37]
也是针对(3)式不能处理横向变速,想办法将其变换到w-x域通过求解差分方程,实现波场延拓.这种方法的最大好处是可以处理任意方向任意变化的速度.缺点是特别陡倾角的轴无法成像.但通过对平方根算子使用高阶近似( Glaerbout,1985)可以改善倾角限制。

马在田( 1982) 提出的平方根算子连分式展开,将高阶近似方程分裂为多个低阶方程,在w-x域用有限差分方法求解这些低阶方程实现波场延
拓,从而进一步解决了大倾角有限差分偏移问题.
( 5)Fourier有限差分方法(FFD)[38,39]
Ristow,et al(1994)提出结合相移法和有限差分方法各自优点在双域操作的一种混合算法.通过加入有限差分校正项对Stoffa et al ( 1990) 的SSP 方法进行一般化并进行改善,以适应剧烈横向变速. 该方法将速度分裂为常速的背景场和横向变化速度场,将向下延拓算子分裂为适于常速度的相位移算子和适于变速度的有限差分算子.
单程波方程双域(波数-间)解法,首先将问题变换到频率域,然后利用快速Fourier变换( FFT)在空间域和波数域之间往返转换,实现双域操作。

波的自由传播是在波数域通过一个具有某种参考速度的均匀介质进行的,参考速度可以纵向变化；而波和介质非均匀部分的相互作用是在空间域中依靠修改波前来进行的.双域操作的方法没有离散化造成的格点频散,计算量比时间-间域的有限差分小得多.
( 6) 广义屏方法(GSP)[40-42]
广义屏方法(GSP)同FFD方法一样,是一个结合双域操作得混合算法,两者的思想和做法一致.目前GSP传播算子主要基于两条途径导出:一是基于散射理论和Born近似[43],另一途径是基于Hamilton 路径积分[44].对双域传播算子的改进提高了大角度波响应得精度,并相应给出这些算子的命名,如拟屏传播算子(Pesudo screen propagator),复屏传播算子(Complex creen propagator),加窗传播算子(Windowed creen propagator),Pade屏传播算子(Pade screen propagator)，高阶广义屏传播算子(High order screen propagator).这些算子在实
际应用中得到很大得发展,能够在非常复杂的地区成像.
三维波动方程叠前深度偏移的基础是三维波场延拓
在将二维单程波方程波场延拓推广到三维时,分裂步方法和相移加插值方法都没有什么限制,FD和FFD则遇到矩阵求逆的困难.提高大角度成像效果,人们探索了两个方向,一个是绕开FD和FFD,找出快速实现的显格式短算子;另一个是探索矩阵求逆的快速算法设法实现三维FD 和FFD.在短算子构造方面张剑峰提出一种优化方法[45],可以较快构造出二维短算子,然后引入Mcclellan变换进行一个二维到三维的变换,求出三维短算子.Clear bout1997年提出了螺旋坐标系谱分解[46]的方法,把多维数据体看成一维数据体,引入螺旋坐标系使得拉普拉斯算子的矩阵表述具有良好性质,采用谱因式分解方法可以比较快地求出矩阵的LU(下三角和上三角)分解,从而找到一种至少在均匀介质中最准确的求逆方法.Ristow和Ruhl则提出用多方向分裂近似求算子的逆[ 39] .
1.3基于波场延拓的波动方程偏移方法的结合
有学者结合射线追踪和有限差分各自的优点提出射线+ 有限差分偏移方法[47]。

该方法利用射线追踪速度快的优点,替代炮记录上波动方程偏移中源子波的正向传播过程,确定地震子波在各个时刻的位置,又利用有限差分法保持振幅的优点对原始记录作反向传播成像.
波动方程地震偏移成像方法研究的进展
现代地震偏移技术的发展已从叠后偏移到叠前偏移,从二维偏移到三维偏移,从时间偏移到深度偏移,但就偏移成像的核心技术仍是波场延拓和速度估计.在围绕继续解决陡倾角和变速问题,构造稳定、快速、精确算法的同时,保结构和保幅波动方程偏移算法成为重点研究方向.
在人们致力于对传统波场表示(傅立叶变换)的波场延拓算子进行分析同时,数学和物理领域新兴的理论和工具为波动方程偏移提供了新的思路和方法,如前文提到的高斯束方法,以及小波分析和现代群论逐渐被应用于地震波场研究中,取得重要进展.
2. 1 相空间小波分析地震波场偏移成像的步进算法[48-53]
研究波传播与成像的常用方法是利用特定的基函数( 如Fourier 谐波,Green 函数)等对波场作分解.这条途径是否有效其关键在于: 变换后的基函数所满足的方程是否易解,边界条件和初始条件是否容易处理.若对波场的展开能满足这个条件,称这种展开为有效表示,例如求均匀介质中的Helmholz方程解时,波场对u(x ,z)对x 做Fourier 展开, 这时在变换域u)(kx ,z)满足的方程易于求解. 但是在非均匀介质中,u)( kx,z)所满足的方程与u( x , z )所满足的方程一样难于求解.在解决非均匀介质中波传播问题时, Fourier变换之所以失效,其根源是Fourier变换中基函数没有局部化性质.应用具有良好时频局部化性质的基函数对波场进行展开,将对波场求解有益.Weyle
eisenberg框架(W - H框架)、短时Fourier变换(WFT) 、小波变换(WT)等近年来发展起来的数学分析方法,由于其基函数具有时空局部化特性,成为函数和算子表示的有力工具.
Steinberg(1993)在频率-间域从标量波动方程出发研究了非均匀缓变介质中波场的步进算法,使用短时Fourier变换和扰动方法, 通过对非均匀介质作部分均匀近似的途径,推导出相空间地震波传播算子,给出步进算法,但计算量较大.在理论上Steinberg(1993)提出的方法开辟了研究非均匀介质中波传播的新的重要途径.高静怀(2000)对Steinberg步进算法进行改进,提出基于W-H框架的波场步进算法,并将这种算法用于地震成像,减少了计算量.吴如山等(2000)将非均匀介质分解为局部背景介质和局部扰动介质,推导出了离散表示的局部化传播算子,将波场分解为一个紧标架(Garbor aubechies 标架)或局部余弦基的表示,并且均有快速算法.用于叠后SEG/ EAEG 盐丘模型成像,取得好的效果.
相空间小波分析类地震偏移试图用具有局部化特征的小波类的基函数取代传统的Fourier变换基函数,从根源上解决非均匀介质波传播和成像问题，值得深入研究和探讨.
2. 2 哈密顿体系下地震波场延拓的辛群算法和李群算法[54-62]
在地震波传播理论和偏移成像研究中,应用辛群和李群算法是中国科学院地质与地球物理研究所以刘洪为首的研究小组探索的一个方向,获得了比较准确、快速的二维和三维波场延拓算法.
牛顿体系、拉格朗体系与哈密顿体系是力学的三个等价的体系. 牛顿体系将运动表示为二阶微分方程;拉格朗体系将力学体系表示为一个变分极值问题; 哈密顿体系则给出在相空间(位置和动量) 的一阶微分方程并指出相空间物体运动满足辛几何.当将它们表示成数值算法时,计算过程和结果是不等价的,有限差分方法是牛顿力学体系表述下的算法,有限元是拉格朗日体系表述下的算法, 哈密顿表述下的算法是辛几何算法.辛几何就是体积不变的几何,而欧氏几何是长度不变. 冯康院士在研究拉格朗日体系后,独立地发现了有限元方法冯康又进一步发现,当把连续的Hamilton体系离散化时,存在一种离散化方法,可以得到离散的Hamilton 系统,这种系统的运动方程满足辛几何性质,这种离散化的方法,就是辛几何方法.地震波传播过程从本质上说是一个具有无穷维自由度的哈密顿体系随时间演化过程,系统在任意时刻状态由初始状态通过单参数辛变换群变换得到.为正确计算地震波传播和偏移成像过程必须采用保持系统基本性质的算法.
对(3)式中平方根算子的不同处理方法是波动方程偏移不同算法的基本原理.对于横向均匀介质在w-x域可对单程波方程算子e 指数进行准确计算,为处理横向变速介质,相继开展了对单程波方程平方根算子的不同近似的方法研究.单程波深度延拓算子是由辛算子不对称分解得到的李群算子,将其表示成一个辛算子与一个窗算子的乘积, 用辛几何方法改进辛算子数值计算精度,用短算子非平稳滤波来加快窗算子计算.已有方法(如FD、FFD)中辛算子的原有计算方法与冯康辛几何方法中的一阶Pade近似相同.将冯康辛几何方法中二阶Pade近。