地震偏移成像基本原理.ppt

地震成像原理与方法,第一章 偏移成像第二章 三维叠前深度偏移成像理论与方法第三章 共聚集点偏移第四章 共反射面叠加第五章 偏移速度建模,第一章 偏移成像,1.1 偏移成像的基本原理1.2 波动方程偏移 1.3 叠前偏移 1.4 偏移速度分析 1.5 深度偏移 1.6 三维偏移 1.7 二维和三维叠前深度偏移,地震技术的发展趋势:,三维叠前深度偏移(3DPSDM)-地震成像(波动方程法3DPSDM,CRS叠加,CFP偏移)四维地震-开发地震(VSP技术,P-S技术,井间地震,3D_AVO技术,4D地震,弹性波阻抗反演,裂缝分析,岩石物理,地震相与地震属性分析,油藏描述等),Reflection point smearing,NMO_DMO_PostMig Correction,三大处理技术:反褶积、叠加、和偏移成像反褶积和叠加引自其它相关学科偏移成像基于古典技术偏移成像:1.具有地震勘探本身的特征。2.计算机使其研究由地震波运动学特征 过渡到地震波动力学特征 3.提高地震空间分辨率和保真度,一偏移成像的概念,偏移,反偏移,反射地震方法:1.激发弹性波，2.记录反射波,3.研究地质岩层结构和物性特征。是一种反散射问题。反射地震成像分做两步:1.记录反射波，2.处理反射波。地震偏移技术是使反射界面最佳成像的一种技术。,1.1 偏移成像的基本原理,地震偏移:叠前或/和叠后偏移叠前偏移:使CSP道集记录或COF道集记录中的反射波归位,绕射波收敛叠后偏移:基于水平叠加剖面，采用爆炸反射面概念实现倾斜反射层归位和绕射波收敛偏移原理和偏移效果见下图,1.偏移成像的基本概念,偏移原理图,偏移过程的定量分析图,2.发展史,1).古典的偏移技术(60年代前)-反射点的空间位置成像;2).早期的计算机偏移技术(6070年代)-定性和概念性地对反射波运动学特征成像;3).波动方程偏移技术(70年代后)-定性或/和定量地对反射波运动学或/和动力学特征成像.,1).有限差分法波动方程偏移：70年代初期，J.Claerbout教授首先提出了用有限差分法解单程波动方程的近似式，用地面观测的地震数据重建地震波在地下传播过程中的波场，从这些传播过程的波场中提取使地震界面成像的那些数据，组成地震偏移剖面。由于这种偏移方法在计算过程中要解波动方程或其近似式，所以被称为波动方程法偏移技术。2).Kirchhoff积分法波动方程法偏移：70年代中期，French和Schneider等在绕射偏移法的基础上使用了波动方程解的Kirchhoff积分公式，发展为地震偏移的波动方程积分法。使绕射偏移建立在可靠的波的基本原理上。因而改善了偏移剖面，取得了良好的效果。3).富里叶变换法波动方程法偏移：70年代后期，Stolt和Gazdag等又先后提出了在频率-波数域解波动方程，外推地震波场的方法。这种方法被称为F-K域偏移方法。由于该方法计算简单，效率高，因而很快得到了推广。,波动方程偏移技术的发展,3.偏移方法分类,二基于射线理论的叠后偏移与叠前偏移,经典的偏移方法和早期的计算机偏移方法 都是基于射线理论 经典的偏移方法只研究到达时间。叠后偏移有圆弧切线法和线段移动法；叠前偏移包括椭圆切线法和交会法等 早期的计算机偏移方法利用了波前、绕射等地震波传播的惠更斯原理，尽管只是定性的、概念性的，但与手工操作法相比偏移剖面除了归位精度提高外，还考虑了波形特征。叠后偏移有波前模糊法、绕射曲线叠加法;叠前偏移有Rockwell偏移叠加法和Paturet-Tariel偏移叠加法等。,1.叠后偏移,叠后偏移:即叠加偏移，是对叠加后的地震记录做偏移。下面介绍圆弧切线法、波前模糊法和绕射曲线(面)叠加法。1).圆弧切线法 一次反射波NMO后,得到时间叠加剖面 由（）得到视深度剖面 如果界面的倾角=0或者很小，例如只有 或更小，则视深度界面就是真深度界面。如果界面倾角不可忽略，则应当进行倾角校正，以求出反射界面的真实位置。,校正的做法是以地面各点为圆心，以各点下至视界面的垂直距离为半径做圆弧，其圆弧族的切线即为校正后的反射界面（v=cont）。,当速度是深度的函数时，例如 为常数时，则圆弧的圆心不位于地面上，而位于地面点的正下方某深度上。这时，圆心的深度和圆弧的半径由下式求出：,（）,2).波前模糊法,波前模糊法也可以称为波前切线法，它是对叠加后的地震剖面进行偏移的方法。这个方法是反推反射界面上的波场。,以地面接收点为中心，把相当于反射到达时间上的值送到以 的深度为半径的圆弧上去。如果我们把深度z仍以双程时间表示，就把反射数值送到以t为半径的圆弧上去（图1-4）。把各道上的所有反射波值都按这个原则去做，并把送到同一点的值叠加起来，就可以组成偏移剖面。把某道 上某时间t上的振幅值送到相邻各道上的时间 由下式算出：,（）,其中,用波前振幅叠加来求反射界面发出的波前实际上就是用这种方法做切线。,要求：较密的地震道和较高的信噪比，以得到满意的偏移剖面。,3).绕射曲线（面）叠加法,绕射曲线或绕射曲面叠加法是把地震剖面上的波场振幅值按绕射波时距曲线进行相加。因为绕射波时距曲线与所有反射波的时距曲线形状相比较，其凸率最大，故亦可称它为最大凸率法。,具体做法是，当要得到地震剖面上某个 点的偏移后的数据时，我们要计算一条以这点为顶点的绕射双曲线。它在各道上的时间t由下式算出：,（）,式中,在进行偏移时我们把各道上等于上式时间t的波场值取出来叠加在 点的波场值上，这就算完成了 点的偏移处理，如图1-5所示。,无论是波前模糊法还是绕射叠加法，其基本原理都是根据惠更斯原理提出来的。波前模糊法是把一个道上的波场值送到各个道上去叠加输出道法；而绕射叠加法是把各道上的相应值取来在一道上叠加输入道法。两者都符合反射波归位和绕射波收敛的要求，而且它们的叠加值也相等。,波前弧或绕射曲线在x方向上的范围L称为偏移孔径。L的范围是由最大实际倾角来决定的。倾角越大，L越大；有效波越深，L也越大。L的大小可用下式来估算：,（）,孔径的中心，原则上应当位于 处，但也可以是不对称的。,图1-7是用绕射叠加偏移法处理前后的地震剖面。从对比中可以看出，偏移后剖面上的地层层位关系得到了正确的反映。有利于地质解释。,2叠前偏移 叠前偏移:即偏移叠加，是对叠加前的多次覆盖的地震记录先偏移，再叠加。下面介绍椭圆切线法、Rockwell偏移叠加法和Paturet-Tariel偏移叠加法。,1）椭圆切线法 当给定CSP记录时，可用椭圆切线法(图1-8)。反射点(2D)位于以炮点和接收点为焦点的椭圆上，这个椭圆的方程可表示为：,（）,对每个炮检距的记录上的反射波画好椭圆弧。做椭圆弧族的切线即为偏移后的剖面。,2）Rockwell偏移叠加法 Rockwell偏移叠加法实际上是叠后偏移所使用的波前模糊法的一个扩展。,具体做法：把每个记录道上任一t时刻的采样值，在以炮检距中点的地面点为原点的直角坐标系中送到以 为长轴，为短轴的椭圆与各个地震记录道垂直线相交的各个点上去，并且与其它地震道送至该交点上的采样振幅值相加，即得偏移叠加剖面。偏移叠加实质上是用振幅叠加来做切线的。,3）Paturet-Tariel偏移叠加法 1971年Paturet-Tariel用相同炮检距的剖面进行叠前偏移，把所有相同炮检距的偏移后的剖面叠加得到偏移叠加剖面。叠前偏移的原理如图1-9所示。,绕射点M所产生的绕射波到达时曲线为：,（）,当炮检距 时，上式表现为：,（）,式中 为从M点到A点的双倍旅行时间。和 的曲线表示在图1-9的右图中。为了进行偏移，我们应当把 的曲线上的地震能量（即采样点振幅）送到零炮检距绕射双曲线的顶点M上去叠加。这样,把各个相同炮检距的剖面偏移后叠加在一起即得偏移叠加剖面。,图1-10.偏移叠加剖面与叠后偏移剖面对比图(a).水平叠加剖面,(b).叠后偏移剖面;(c).偏移叠加剖面,三基于波动方程的波场外推与地震成像原理 使用波动方程进行偏移，首先就是要重建反射波的原来波场。反射界面上刚刚产生的反射波，就认为是该反射面的像。为进行波场外推，把波动方程分解为上行波方程和下行波方程。1上行波和下行波,波动方程有两个解,一般表示为。在地震勘探中一般取深度方向向下为正z的方向。向正z方向传播的地震波称为下行波，即用 所表示的波。向负z方向传播的波为上行波，即用 代表的波。下行波即入射波，上行波为反射波。,只有在均匀各向同性完全弹性介质的情况下上行波和下行波才是分离的。分离过程如下:二维波动方程为：,（）,对(1.1.9)式相对x和t做二维付里叶正变换，并进行算子分解得到：,（）,其中利用了波散关系：,（）,由（）式得出：,（）,其中，正号代表上行波方程，负号代表下行波方程。,2波场外推 正向外推就是根据波在当前位置上的振动情况向波的自然传播方向用计算手段预测出波场。反向外推是向波的自然传播方向的反方向上重建原来的波场。对一个波场应是进行正向外推还是反向外推均有物理问题决定。,1）上行波的外推,（）,积分结果为：,（）,由此得出上行波的正、反向外推式。,（1）上行波正向外推公式 上行波的正向外推式就是向负z方向的外推公式。从（）式可求出为：,（）,根据这个公式可以计算模拟反射波的地震记录（地震图）。（2）上行波反向外推公式 上行波的反向外推式就是向正z方向的外推公式。从（）式可得出为：,（）,根据这个公式可以进行地震记录的向下半空间延拓，求出地下任何一点的波场，实现地震波偏移的目的。,2）下行波的外推,（）,积分结果为：,（）,据此可以得出下行波的正、反向外推公式。,（1）下行波正向外推公式 下行波的正向外推式是指沿正z方向的外推。其外推式为：,（）,这个方程可用来模拟下行波的地震记录。（2）下行波反向外推公式 下行波的反向外推是指沿负z方向的外推。其外推式为：,（）,上式可用来从下行波场进行反向求源的计算工作。,下面分析波场本身的条件对外推结果的影响,（）,当 时，为正或负的实数，这时所有外推公式中存在虚指数。说明在外推过程中波场发生相位变化。一般都能得出正确的结果。,当 时，值为虚数：,（）,（）,波场外推时只有振幅变化，而无相位变化。当指数项取负号时，外推的波场迅速衰减，称这种波为倏逝波。当指数项取正号时，外推波场迅速增大，这是一种实际不存在的波，只是进行波场计算时发生，我们称它为耗损波。在计算中要避免发生这种情况。,当 时，上行波的外推式可写为：,（）,此时反向外推遇到倏逝波，正向外推发生耗损波。分别表示为：,（）,（）,由此可见，用上行波方程进行向下波场外推永远是计算稳定的。而用上行波方程进行正向外推就可能遇到耗损波，因此有可能是不稳定的。除非在计算中不断地把 的波场滤除掉。同理可求出 时下行波的外推式为：,（）,此时也是反向外推遇到倏逝波，正向外推遇到耗损波。,3）波场外推的Kirchhoff积分法,Kirchhoff积分法并不直接解波动方程，而是用数学方法来描述关于波的传播的惠更斯原理，从而求出空间上任一点波场值的。,Kirchhoff早在1883年就证明了，从扰动区向外某点 传播的波的t时刻的波场 可以从扰动区封闭表面上的波场 以及该波场对时间和表面法线方向的导数通过积分式求出来。因此要假定 在封闭面上和封闭面内有直至二阶导数的连续性。Kirchhoff利用了格林定理：,（）,u取为波场函数，,当把观测点用包含有波前面在内的封闭曲面包围起来，如图1-11（a）,（b）那样的封闭时，这样的封闭面S和它所包围的体积V作为（）式的积分限，经过一定的推导后得出 点的正向外推波场为：,（）,这里 的方向取封闭表面的外法线方向。如果把观测点M移至封闭面外，则有：,（）,（）式中,（）式就是著名的Kirchhoff积分。它描述了物理波场传播的过程，也满足奇次波动方程，是它的积分形式解。对我们来说，也可以称它为正向外推公式。注意:Kirchhoff积分只满足均匀介质的情况。,是推迟场,下面讨论用Kirchhoff积分进行波场反向外推问题(地震偏移)。这时，所取的封闭体积V应在波前传播方向的反方向，计算点 就在这个封闭体内。根据格林定理同样可求出形式上相同的反向外推的Kirchhoff积分式：,（）,式中的u不再是推迟场，而是超前场。,（）式为用于波场反向外推的Kirchhoff积分式。它可用于上行波的反向外推，也可用于下行波的反向外推。当然，这种外推与正向外推不同，它不代表一个物理过程，而只是一种重建波场的计算过程。,3地震反射波场成像,从波动场的观点叙述反射波成像的一般原理。地震成像-地震偏移反射系数值-反映该反射点反射系数相对值的反射波振幅反射成像实际上就是把地面上观测到的反射波归位到产生它的反射点上去。地震偏移与地震成像在现阶段可以视为同一概念。,地震偏移成像:一是上行波场的反向外推;二是在外推波场中提取成像值。Claerbout提出下述反射波成像原则：反射面位于这些点上，其入射波的初至与反射波的产生时间相同。如图1-12所示,反射波成像的基本公式可写为：,（）,（）式没有考虑反射系数随着入射角变化的情况，它实质上是相位信息的公式。或者说，它对接近法线入射的情况时基本是正确的，能够反映反射系数在各点上的变化情况。应用（）式涉及到要选择下行波的初始时间。这是一个困难问题。我们通过假设下行波是最小相位而避开这个问题。我们把 作为初始时间，可推出如下的反射图象公式：,（）,当下行波是脉冲波时,(1.1.38)式很精确。但是，如果 是 一个短延续长度的子波时,它只是一个很好的近似成像公式。,爆炸反射面的概念:,水平叠加剖面后的自激自收剖面等价于在反射界面上同时爆炸产生地震波，并以半速度向外传播，在地面上观测到的上行波剖面。这就是Loewenthal等人首先提出来的爆炸反射面的概念（图1-14）。这个概念对于理解水平叠加剖面的偏移成像是很重要的。因为它比较直观地说明了这种剖面的成像原理。比前面所述的反射波成像的一般原理要容易理解得多。在这里我们引入了两种等价。一个是水平叠加剖面等价于自激自收剖面。另一个是自激自收剖面等价于界面上同时激发在地面记录的上行波剖面。这种等价只是概念上的，实际上只有一种水平叠加剖面，并没有三种剖面。,1.2 波动方程偏移,地震偏移成像技术发展到今天已经产生了各种形式的在各种域实现的方法。历史上曾经起过作用的根据几何光学原理的成像方法已经被淘汰。现在正在流行的是建立在波动方程基础上的三种方法，即Kirchhoff积分法，有限差分法和F-K法及其各种变形。这三种方法由于有相同的数理基础，因此它们的原理相同。同时，因计算方法不同，它们之间又有许多不同之处。下面讨论三种方法对水平叠加地震剖面的偏移。,一频率-波数域波动方程偏移,采用爆炸反射面的理论。为了成像，要求向地面以下反向外推地震波场。假定z轴垂直向下为正，测线沿x轴，则u(x,z,0)表示偏移后的真实剖面，而u(x,0,t)是未偏移的叠加剖面。,在均匀各向同性完全弹性介质中，用半速度代替地震波传播速度，则标量波动方程变为：,（）,（）,对（）式进行傅里叶变换并利用（）式有,（）,其中正号代表上行波，负号是下行波。,1Stolt偏移法,设 为 的二维傅里叶变换，对（）式进行上述变换得到：,将（）式代入上式有,按上行波求解，即取正值得,其中A与t无关。令t=0，上式变为：,从而，是待求的偏移剖面 的傅里叶变换。,下面讨论用水平叠加剖面 如何求出。对 做傅里叶逆变换得：,令z=0，上式变为：,（）,设水平叠加剖面 的二维傅里叶变换为，则,（）,其逆变换为：,（）,比较（）与（）有,这样,按上行波 取正号并对 微分得,（）,对 做二维傅里叶逆变换得到：,（）,就是要求取的偏移剖面。,输入零偏移距剖面,二维付氏变换,公式(1.2.5)得到,用公式（）把 映射成 并 标定振幅，得到,二维付氏逆变换,公式(1.2.8)得到,偏移剖面,图1-16 均速Stolt偏移流程图,上述偏移原理见图1-15。由图1-15和（）式可看出，在每个频率 移向新的频率 时，要乘上一个振幅比例。通过这个频率移动，把视倾角 转换为真倾角。其流程见图1-16。上述频率波数域的偏移方法称为Stolt偏移方法。Stolt法的偏移效果见图1-17，1-18和1-19。,2Gazdag相移法 对标量波动方程（）相对x和t做二维傅里叶变换得到：,（）,式中。求解（）式得出F-k域的向下外推公式,（）,偏移成像公式是把上式变换回到空间-时间域，并取t=0时刻的波场值为成像值。即,Gazdag相移法的流程见图1-20。偏移效果见图1-21。,输入零偏移距剖面,二维付氏变换,在新的z处计算（kx,z,）,用相移算子exp（-ikzz）将地面数据向地下延拓,相加所有频率（成像原理，t=0），（kx,z,t=0）,在x方向作付氏逆变换,偏移剖面,图1-20 Gazdag的相移偏移法流程图,二克希霍夫积分法波动方程偏移,前面导出了波动方程边值问题的Kirchhoff积分解。下面研究把它用于地震成像问题。现在转写反向外推的Kirchhoff积分如下：,（）,图1-22 求地震问题Kirchhoff积分解图示,因为在（）式中需知道，即波场在地面法向的导数值。但是，这个导数值目前是无法观测和计算的，因此，需想法去掉含 的项。为此不再用（）式，而从格林定理一般式：,（）,出发，设格林函数w为：,（）,来代替1/R，则可以达到目的。上式中：,（）,（）,把（）式代入（）式，得,（）,由此求出向下外推的Kirchhoff积分为：,（）,式中 A为地面的面积。,求下面的微分：,（）,（）,把以上二式代入（）式，得,（）,当我们把z取成地面上的点时，即z=0时，则有：,（）变为下列形式：,（）,式中,A为地面的面积。,（）式又可写为：,（）,（）式与下式等价：,（）,式中,现在我们来证明（）与（）式等价。我们先对 求导，再对 求积分。该过程如下：,（）,（）式与（）式完全相同，因此（）式也与（）完全相等。由（）式可知,（）,将（）式代入（）式，得到,（）,根据褶积的定义，我们把（）式写成三维褶积符号形式，则有,（）,其中,把（）式对x,y和t进行傅里叶变换，则可写为：,（）,式中褶积算子H为：,（）,对t积分得到：,（）,括号内积分是第一类Hankel函数，上式可写为,（）,利用圆柱函数间的关系：,式中 和 是Bessel函数和Neuman函数，把它们代入对y的积分式（）中，则最终得到H的表达式为：,（）,把（）式代入（）式中，得,（）,上式与频率-波数域的向下外推公式一致。因此，在常速介质中Kirchhoff积分法与频率-波数域的波场向下外推公式完全等价。,下面讨论利用Kirchhoff积分法对水平叠加剖面进行波动方程偏移的步骤。将水平叠加剖面看做是炮检距为零的自激自收地震剖面u(x,y,0,t)；利用爆炸反射面的思想将自激自收剖面等效为在反射界面上同时激发产生地震波,以半速度向外传播，在地面上观测到的上行波剖面u(x,y,0,t)；利用（）式将单程的上行波剖面u(x,y,0,t)向下延拓，得到深度为z的面上的波场值,（）,根据成像原理，对所有地下点（z0）取t=0时的波场值，即可实现三维偏移成像。此时，成像值为,（）,Kirchhoff积分法的偏移效果见图1-23。,三有限差分法波动方程偏移 下面讨论使用有限差分法对水平叠加地震剖面的偏移问题。为了把上行波方程表示为空间-时间域的表达式，需要把上行波方程表示为某种近似式。然后在空间-时间域研究其差分方程及求解问题。最后讨论一些计算方法和效果。,1上行波的空间-时间域方程 为了适应介质速度的空间变化，我们要在空间-时间域中进行偏移成像或地震图的模拟工作。首先就要把上行波方程表示在空间-时间域中，这需要用到某种根式展开。1）二项式展开,下面我们将用到 这样的二项式展开，在这里我们介绍几种展开式。,（1）Taylor展开 这是一个众所周知的显式展开式，它一般表达为：,（）,展开条件。,如果把这种展开式用于微分算子，在不进行辅助处理时将找不到稳定的有限差分方程来解相应的微分方程。因此我们在使用二级近似以上的展开式时不能用这种展开式。,（2）连续分式展开，或称为Pade展开 这个展开式表示为如下形式（）：,（）,这是一种隐式展开式。其各级展开式如下。,一级展开式：,（）,二级展开式：,（）,三级展开式：,（）,高级展开式可依此类推。,（3）迭代展开 这种隐式展开法，是把前一级的展开结果代入下一级的展开式中。设,则逐次迭代展开式可表示为（）：,（）,用这种展开方法得到的各级展开式如下。,一级展开式：,（）,二级展开式：,（）,三级展开式：,（）,高级近似式可依此类推出来。从（）和（）公式组可以看出，后两种展开是等价的。,2）上行波的空间-时间域方程 在第一节已经求出了频率-波数域的上行波方程(1.1.12)式：,用迭代展开法展开上行波方程：,（）,式中,由（）式求出各级近似式如下。,一级近似式：,（）,二级近似式：,（）,三级近似式：,（）,高级近似式可依次类推。现在，我们把（）式转换到空间-时间域，求出一级近似方程。上行波方程（）式可改写为：,（）,对（）式进行傅里叶反变换：,（）,根据傅里叶变换的微分性质：,（）,（）,（）,把(1.2.43a)、(1.2.43b)和(1.2.43c)式代入(1.2.42)式，得到：,（）,上式就是空间-时间域的一级近似的上行波方程，常常被称为 方程。,同理可求出空间-时间域的二级及二级以上近似的上行波方程。经推导，空间-时间域的二级近似的上行波方程为：,（）,上式常常被称为 波动方程。,3）浮动坐标系中的单程波方程上行波方程在一定的浮动坐标系中可以简化。我们对二维波动方程：,（）,进行如下的坐标转换：,坐标变换前后波场本身是不变的，因此存在：,（）,（）,从（）和（）导出下列导数等式：,（）,（）,（）,将（）各式代入（）式中得到新坐标系中的波动方程为：,（）,上式变换到频率-波数域为：,（）,（）式可改写为：,（）,从上式得到下列关系式：,它表示坐标变换前后的算子关系。因此，上行波方程可表示为：,（）,（）式的右端项可展开为各级近似式，便得到上行波各级近似方程。,一级近似式,（）,二级近似式,（）,三级近似式,（）,高级近似式可依次类推。由以上各式用前述方法可求出空间-时间域的各级近似方程。下面给出一级和二级近似方程。,一级近似方程 用推导（）式那样的方法可以求出一级近似方程为：,（）,与（）式相比，减少了一项。从而也减少了计算时间和差分时的时间层（少了一层）。另外，保持了计算的稳定性。,二级近似方程用推导（）式那样的方法可以求出二级近似方程为：,（）,这个方程与（）式相比也是少了一项。这也会减少计算工作量。,2有限差分法地震偏移技术 如前所述，水平叠加地震剖面可以看做是自激自收地震剖面；又可以看做是所有反射面同时爆炸产生波源向地面传播，被地面的接收器记录的上行波剖面。对于这种观测结果，为了成像，要求向地面以下反向外推地震波场。在外推过程中假设地震剖面上无任何多次波，也不存在任何规则干扰波，如折射波等。如果在剖面上存在这些波，在外推过程中也都按反射一次波处理，但它们是不能正确归位的，只能造成偏移成像剖面的干扰。因此，如果存在这些波，应当在偏移处理前把它们滤掉。,1）浮动坐标下的有限差分法地震偏移 采用浮动坐标系,只讨论一级近似的上行波二阶偏微分方程（）的有限差分偏移问题。考虑到爆炸反射面的概念，用v/2速度代替v。这样，（）式可重新写成：,（）,这里的速度，假设它是常数，在实用中它可以随 和 而变。,(1.2.59)式存在于空间 内，给定的混合条件为：,（）,（）,（）,（）,式中 即为在地面所观测的地震波场。,目的:通过解上述微分方程求出地面以下任何点（）上的曾经在该点出现过的上行波的波场值(位移或压力场振幅)。这样的问题是适定的，可解的。,下面用有限差分法来解这个方程。采用对称隐式(Crank-Nicolson)的差分格式（图1-24）。,实用中，常采用下列符号：,图1-24 差分格式,对图1-24格式的中心点进行差分，（）式可化为如下的差分方程：,（）,式中,由（）式整理得：,（）,为了求出以整采样点为依据的表达式，（）式可写成：,（）,从（）式可以解出求 的递推式。把求 做为递推结果值，不取 作为递推结果值是从物理条件（）式考虑的，这也是计算稳定性所要求的。（）式经过简单整理后得：,（）,式中,算子,（）式写成矩阵式为：,（）,式中、为三对角矩阵。为在 深度层和 时间层上沿 轴的波场值的列向量。，和 为相应层上沿x轴的波场值列向量。,这些矩阵和列向量表示如下：,（）,（）,（）,（）,（）,（）,如果把差分算子 取为，这相当于差分加权，可提高差分精度。这时差分方程（）可改写为：,（）,式中。,上式可以写成矩阵式(1.2.65),不过，矩阵A和B由下面的矩阵表示。,（）,（）,只要矩阵A可逆，（）可解。得到：,（）,矩阵A中的主对角元素是优势元素。因为 永远大于零。因此总有：,或者,因而矩阵A的行列式不为零，即,所以矩阵A是可逆的。,（）,矩阵方程（）式可用矩阵方法求解，也可以改写为代数方程组后用追赶法求解，要求：差分方程必须稳定。,在浮动坐标系中，成像时间是。,2）一般坐标下的有限差分法地震偏移 若按一般坐标系，叠加地震剖面的有限差分偏移成像过程可用图1-25说明。图1-25为（x,z,t）坐标系，地震剖面在（x,t）平面上，偏移剖面则在（x,z）平面上。有限差分偏移是按一定步长的z来外推地震剖面（x,t）,每外推一个步长，就将t=0的波场作为输出。这些输出结果就组成了偏移剖面（x,z,t=0）。3）有限差分法逆时偏移 Baysal等人在1983年提出了有限差分逆时偏移的方法，它是从一个波场为零的（x,z）起始平面，按时间反推，并以地震剖面资料u（x,z=0,t）作为每一步进时间的边界条件（z=0）,得出时间t=0的（x,z）平面就成为偏移结果 u（x,z,t=0）（见图1-26）。,图1-25.地表z=0的地震剖面（x,t），它向下延拓获得各个离散深度上的时间剖面。这里用粗黑箭头指出外推方向，偏移获得的剖面用t=0（根据成像原理）的（x,z）平面表示。,图1-26.逆时偏移从数据体底部的全零（x,z）平面开始，按时间向t=0反推，计算出不同时间的（x,z）平面切片；这些地下切片在图中用一系列水平面来表示，反推方向按粗黑箭头所示，每个时间平面（x,z）都包含有出自地震剖面的边界值（虚线表示的z=0平面上的x线），t=0的（x,z）平面即为偏移剖面（顶部的水平面）,3吸收边界条件仍然讨论二维情况。二维波动方程：,（）,一般求解域：,实际求解域：,在求解的过程中，一般给定如下的边界条件：,（）,（）,（）,（）,这样，就在两边和底界人为地造成了边界，称之为计算边界。这样的边界不可避免地会产生边界反射。这个边界反射是很强的。在现代偏移技术中为了避免在有用的地震剖面范围内出现边界反射波，常常要在地震剖面的两边进行扩边，扩边道有时要达到96道或更多。这就会使计算量和设备资源使用量增大。为了克服这个不足，提出了使用吸收边界条件的算法。首先，分析一下波入射到边界上的情况。考虑一个入射到右边边界上的简谐平面波：,（）,式中 为平面波波前与x轴间的夹角，k为波数，等于。,在 的区域上的反射波为：,（）,在 区域上总的波场为：,（）,把（）代入（）中，则得到反射系数：,说明反射系数是很强的。为了导出吸收边界，必须采用一种算子B，使得它作用在边界上的波场时，波场值等于零，即,（）,考虑到我们的边界条件是线性的，可以求出反射系数：,（）,（）,或,从上式可以看出，要选择这样一个算子，当它作用在波场上时，使界面上入射的波如同无边界那样，就不会产生边界反射了。为此要使，实际上就是要把波动方程分解为左行波方程和右行波方程。把右行波方程取为左边界的边界条件，把左行波方程取为右边界的边界条件。这与把波动方程分解为上行波和下行波方程的做法是一样的。,1）吸收边界条件的推导把波动方程（）式对z和t做二维傅氏变换，得,（）,进行算子分解，得出,（）,由此得出左行波方程为,（）,右行波方程为：,（）,把以上二式直接变换到空间-时间域得不到单程波方程，因此要把根式展开。,2）各级近似的边界条件（1）零级近似式零级近似式为：,（）,（）,和,经过反傅里叶变换到空间-时间域得到：,（）,（）,和,（）式即为左边界的吸收边界条件。（）式是右边界的吸收边界条件。这样简单的边界条件只能衰减向边界界面上垂直入射的波。这种波在地震剖面上是很少的。但是计算方法简单。为了衰减倾斜入射到边界界面上的波，可以使用高级近似式。（2）一级近似式一级近似式用(1.2.37)的方法可展成为：,（）,（）,和,以上二式经反傅里叶变换后的空间-时间域表达式为：,（）,（）,和,(1.2.94)式是左边界的边界条件。(1.2.95)为右边界的边界条件。同理可以求二级近似式的边界条件。（3）二级近似式二级近似展开式经过相应的计算表示为：,（）,（）,和,经过反傅里叶变换后，（），（）表示为：,(1.2.98)为左边界的边界条件。(1.2.99)式为右边界的边界条件。,（）,（）,和,由于在进行偏移时，使用的是浮动坐标系，所以应当把上面求出的吸收边界条件表示在浮动坐标系中。3）浮动坐标系中的吸收边界条件 根据浮动坐标系与原坐标系的关系式（）进行导数转换，可求出浮动坐标系中的吸收边界条件。（1）零级近似式,（）,（）,和,（）式为左边界条件，（）式为右边界条件。这个公式与原坐标系中的边界条件相同，未改变形式。（2）一级近似式,（）,（）,（）式为左边界条件，（）式为右边界条件。（3）二级近似式,（）,（）,（）为左边界条件，（）为右边界条件。,浮动坐标系中的吸收边界条件的公式可以进一步简化。简化的方法是，去掉含有对z轴的二阶导数项，然后把各项相同的导数从式中去掉。经过这样的简化后，浮动坐标系中的吸收边界条件可表示如下。（1）零级近似吸收边界条件,（）,（）,（2）一级近似吸收边界条件,（）,（）,（3）二级近似吸收边界条件,（）,（）,（）和（）只适合于衰减以直角入射到边界上的边界反射波，最多不应大于 的夹角（波前与边界面的夹角）。（）和（）适合于衰减夹角 及以下夹角的边界反射波，最多不超过 的夹角的边界反射波。对于大于这个夹角的边界入射波应当采用（）和（）式那样的边界条件。,4有限差分法偏移的效果 有限差分法波动方程偏移是地震成像技术上的一个飞跃。方法原理-基于波动方程处理效果-除反映相位关系外，还保持反射波的振幅特征与绕射叠加法相比，有限差分法具有以下几个方面的优点。基于波动方程，准确的偏移定量计算式，反射波正确归位，振幅相对保真。适于以研究波的特征为主的地质解释，特别有利于研究岩性、岩相变化，含气砂岩，油水接触面等。差分网格要小。小到,大约100 1000 取一个速度参数。而绕射叠加偏移要求在大范围内使用平均速度。且反射面越深，范围越大。要求大约 50 范围的速度值是不变的。对实际介质难于满足。因此，其归位效果不理想。,偏移综合效果，如S/N、波形特征、分辨力等都好于绕射法 因此，有限差分法从70年代中期就取代了绕射叠加法。今天，有限差分法偏移程序已成了一种常规的地震处理应用软件。有限差分法的偏移效果见图1-27和1-28。由图可看出偏移剖面的信噪比高，断层清楚，波形特征得到了保持。频率域有限差分法（频率空间域偏移）的偏移效果见图1-29。有限差分法偏移技术除了具有上述优点外，也有不足之处。具有一级近似方程的 有限差分法的使用有倾角限制，使用高阶方程偏移,克服对倾角的局限性（图1-29）。差分方法常常由于水平方向上采样不足会引起网格频散，即波的高频成分与低频成分偏移不到同一位置上去，高频成分变成一种干扰背景。克服办法:野外设计合适的水平采样密度;处理中进行空间道内插。空间道内插有几种实用程序。如上所述，有限差分法正在走向完善的过程中。,四.三种波动方程偏移方法的差异1偏移孔径不同 Kirchhoff积分法一般据偏移剖面上的倾角确定偏移孔径。理论上可取成满足 倾角的要求。但实际总是要小。浅层一般取 以内。深层要大些，但要以最大倾角为依据。否则，或者增加工作量，或者增强偏移噪声。F-K域偏移没有孔径限制，可自然满足 倾角的偏移。它可通过在F-K域中的二维滤波控制偏移孔径。有限差分法可通过数值粘滞性控制孔径，其实质也是一种二维滤波。另外，常用近似方程。实际偏移范围受方程限制。所用方程不同，偏移孔径的角度分别为 等。超过它们所允许角度的数据用数值粘滞性滤除，否则产生偏移噪声。,2对速度模型的适应性不同 给出的偏移速度模型要适合实际的地下地层的速度变化。总是希望偏移方法能适应速度的变化。有限差分法的速度适应性最强。它可在一个差分网格内取一个速度，另一个差分网格可用另外的速度来进行计算。例如，一般是在x方向两个道距，在z方向一个外推步长上用一个常速度值。如果 取20米，取60米，则差分网格面积为 平方米。在这个小面积内取常速值。向x方向移动一道，,形成新网格又可取另一个常速值。x方向每次要重叠一个道距，因此,在x方向上的速度变化不应是很剧烈的。z方向上差分网格向下外推时不重叠，速度变化可稍大些。总之适应速度的空间变化是有限差分法的显著优点之一。,图1-30 有限差分法速度变化图,F-K域法偏移有两种实现方法，即Stolt法和Gazdag相移法。Stolt法在波场外推时每次都是从地面的观测结果向下外推，所以它要求在每次外推时对全测线要使用一个平均常速度值。为适应速度变化，要在偏移前对未偏移的剖面进行时深转换。当然，这是一种近似办法，并不能正确地反映速度实际的空间变化。相移法可在深度方向上变化，但不能沿水平方向适应速度变化。因为它的每次向下外推是从上一次的外推结果计算的，所以每次可改变一个速度。关于适应横向变速的相移加内插法见参考文献1。Kirchhoff积分法需在一个大孔径内进行一次运算，而在这个孔径内只能用常速。因此，原则上不能变速。不过已有人提出速度可变的改进算法。以上讨论，或在x-t域，如有限差分法和Kirchhoff积分法，或在F-K域，如F-K域法。它们为三种基本偏移方法。另外，在()域的傅里叶法和横向滤波法，可在横向和垂向上实现变速。,3偏移成像的综合效果 最终效果受各种因素制约。影响因素有：所用方程的精确程度；方法对速度模型的适应性和计算方法与参数。方程越准确，原则上应越有好的偏移效果。若所用方程能满足实际介质倾角与绕射现象的成像，即使使用近似方程，也不会对偏移效果有显著影响。若所用方程不能满足上述要求，则偏移效果不佳。参考文献1对各阶方程的成像振幅的误差进行了分析，说明方程越精确，成像振幅的保真度越高。Kirchhoff积分法和F-K域方法使用准确方程。而有限差分法用近似方程。但若用高阶方程，如五阶方程，也能基本满足振幅保真的要求。方法对速度模型的适应性越好，偏移效果越佳。假定速度模型是真实的，方法对速度模型不能适应会引起全面的偏移误差，如位置不准，振幅不保真，波形特征不好，偏移噪声增强等。在速度适应性上，有限差分法最好，相移法等次之，,Kirchhoff积分法较差。当然，在（）域采用傅里叶变换求导数的方法和横向滤波方法都能比较好地适应速度的空间变化。应当指出，这种适应一般也是指速度变化不剧烈而言。当速度变化剧烈时还要采用深度偏移的方法。计算方法与所用参数不同将对偏移效果产生不同的影响。如有限差分法会产生频散，孔径选择的不适当会使Kirchhoff积分法偏移效果降低。水平方向采样不足是各种方法的偏移剖面上出现假频现象的主要原因。偏移剖面的最终效果是综合性原因造成的，不能以一种因素来决定使用何种方法进行偏移处理工作。偏移方法的效率是能否推广使用的一个重要因素。由于现代地震数据处理是一个工业性生产过程，它要日复一日、年复一年地进行大量的重复工作。因此，效率就显得大为重要。在同样偏移效果的前提下肯定是选用效率高的那种方法进行处理。,从全局看，Stolt法效率最高，其次是Kir