什么是波场延拓 波动方程偏移
反射地震资料的偏移处理方法，波动方程偏移，基于二维地质建模的两种地震数值模拟方法的应用及分析。
本文导航
反射地震资料的偏移处理方法
反射地震资料的偏移校正、射线偏移和波动方程偏移等方法统称偏移处理。偏移处理可使倾斜界面的反射波，断层面上的断面波，弯曲界面上的回转波以及断点、尖灭点上的绕射波收敛和归位，得到地下反射界面的真实位置和构造形态，得到清晰可辨的断点和尖灭点。因此，偏移处理对提高地震勘探的横向分辨率具有重要的作用。
偏移处理通常又可称为偏移归位、偏移成像、波场延拓成像等。从偏移原理考虑，可将偏移处理分为射线偏移和波动方程偏移两大类。若以偏移处理的流程，也可分为叠前偏移和叠后偏移。一般对三维资料，应用三维偏移处理，对二维资料则用二维偏移方法处理。以下讨论主要以二维偏移方法为主，对三维偏移处理方法，可在二维方法的基础上扩展即可。
10.4.1 偏移的概念
反射波水平叠加剖面相当于自激自收记录剖面，在叠加剖面上的反射波同相轴与地下的反射界面有关。当反射界面水平时，反射波同相轴与地下界面形态一致，如图10-11（a）所示。当反射界面倾斜时，反射波同相轴则与反射界面形态不一致，若直接将反射时间作时深转换，所得视界面为，如图10-11（b）。与地下真实反射界面 R1 R2 比较，不论是界面长度、界面位置及界面倾角两者均不一致。视界面相对界面 R1 R2，向界面下倾方向偏移，而且倾角变小。我们称这种现象为偏移现象，R1 至的水平距离称为偏移距。偏移现象随反射界面的埋深和陡度增加而越严重。由于偏移现象的存在，当地下构造复杂时，自激自收剖面上反映的视界面因位置不正确可能产生界面交叉重叠或出现空白带。如图10-12，背斜界面出现空白带，而向斜界面出现界面交叉重叠。
图10-11 偏移的反射分析
另外，根据绕射理论，在断点，尖灭点等岩性突变点还会产生绕射波，这些绕射波再与偏移后的反射波叠加，就使得水平叠加地震剖面上的反射波变得很复杂。若直接用水平叠加剖面解释地下界面，很难得出正确的结论。可见偏移现象使地震剖面的横向分辨率降低。若能使偏移后的波场归位，绕射波收敛到绕射点，就可恢复反射界面的真实形态，因此偏移处理就是针对偏移现象的反偏移方法。偏移处理经常也简称为偏移，但含义和前面提到偏移是不一样的。
图10-12 反射界面位置不正确造成空白或干涉
对非零炮检距地震记录的偏移处理（叠前偏移），可理解为实现地震波传播的逆过程，让波场反向传播（称为延拓）。相当于将激发点和接收点平面逐渐向地下移动，随之炮检距也变小。当激发点和接收点移至反射点时，炮检距为零，这时的激发点和接收点位置为偏移处理后的反射点真实位置。
10.4.2 射线偏移方法
射线偏移是建立在几何地震学基础上的一类偏移方法。可实现叠后偏移，也可实现叠前偏移。偏移的基本原理可用地震脉冲的偏移响应来说明，偏移方法可分为圆法偏移、绕射扫描叠加偏移以及椭圆法偏移。
10.4.2.1 偏移脉冲响应
偏移脉冲响应可分为输入剖面脉冲响应和输出剖面的脉冲响应。输入脉冲响应指在输入的时间剖面（水平叠加剖面）或共炮集记录上的一个脉冲，在深度剖面上可能存在的界面位置轨迹。例如在均匀介质中，水平叠加时间剖面上的一个脉冲在深度剖面上的偏移脉冲响应是一个圆；而炮集记录中某一道中的一个脉冲对应的偏移响应是一个椭圆。输出脉冲响应指在输出的深度剖面上的一个脉冲，对应时间剖面上的时间轨迹。例如在均匀介质，目标空间有一个脉冲（或绕射点），在自激自收时间剖面上脉冲响应为一绕射双曲线；在非零炮检距的炮集记录上时距曲线也为双曲线，但两个双曲线的计算公式不同。根据偏移脉冲响应可形成一系列射线偏移方法。
10.4.2.2 绕射扫描叠加偏移
在多种射线偏移方法中，以绕射扫描叠加偏移为例，介绍如下。
利用这一方法作偏移处理时，只考虑几何关系，将绕射双曲线上的能量会聚于其顶点。首先，将地下空间划分为网格，认为每个网格点都是绕射点。根据网格点坐标计算出它的绕射波时距曲线
勘查技术工程学
式中：（x，z）是绕射点R（地下网格点）的坐标；而（xi，0）是接收点坐标；t0为绕射点至地表的双程垂直传播时间。然后按此绕射双曲线的时距关系ti-xi在实际记录道上取对应的振幅值，将它们相加后放置在绕射点R处，作为偏移后该点的输出振幅（图10-13）。依次对每个网格点都作如上处理就完成了绕射扫描叠加偏移的工作。
如果R点是真正的绕射点（界面点），则按绕射双曲线取出的各道记录振幅应当是同相的，它们相加是同相叠加，能量增强，偏移后R点处振幅突出。若R点不是真正的绕射点（非界面点），则参与叠加的幅值是随机的，叠加结果必然会相互完全抵消或部分抵消，从而使R点处振幅相对较小。因此，偏移后的剖面上，绕射波自动收敛到其绕射点处，在有反射界面处振幅变大，无界面处振幅自然相对减小，显示出了真实反射界面的位置（绕射双曲线顶点连线）。
图10-13 绕射扫描叠加的双曲线示意图
10.4.3 波动方程偏移方法
射线偏移是一种近似的几何偏移，虽然地震波的运动学特点得以恢复，但波的动力学特点（如振幅、波形、相位等）却受到畸变。因此，射线偏移已逐渐被高精度的波动方程偏移所代替。波动方程偏移是以波动理论为基础的偏移处理方法，其基本思路是，当地表产生弹性波向下传播（称为下行波），遇到反射界面时将产生反射，这时可将反射界面看作新的波源，又有新的波以波动理论向上传播（称为上行波），在地表接收到的地震记录就可看作反射界面产生的波场效应。偏移就是将地表接收到的波场按波动方程的传播规律反向向下传播，通常称为波场反向延拓，当波场反向延拓到反射界面时成像（成像剖面为偏移剖面），从而找到了真实反射界面，达到了偏移处理的目的。可见波动方程偏移主要由波场延拓和成像两部分组成。波场延拓可用多种不同的方法实现，随之形成了多种不同的波动方程偏移方法。成像也有成像的原理，叠前和叠后偏移各有不同的成像条件。
10.4.3.1 波动方程偏移的成像原理
（1）爆炸反射界面成像原理
该原理属叠后偏移成像原理。叠加剖面相当自激自收剖面，若将剖面中时间除2，或将传播速度减一半，就可将自激自收剖面看作在反射界上同时激发的地震波沿界面法线传播到地表所接收的记录，即可将界面看作爆炸源，称为爆炸反射界面。若用波动方程将地表接收的波场（叠加剖面）作反时间方向传播（向下延拓），当波场延拓到时间t为零（t=0）时，该波场的所在位置就是反射界面位置。因此，t=0成为叠后波动方程偏移的成像条件。从延拓的结果（地下各点的波场）中取出地下各点处零时刻的波场值组成的剖面就为成像剖面，该剖面为叠后波动方程偏移结果。
（2）波场延拓的时间一致性成像原理
时间一致性成像原理适用于叠前偏移。此成像原理可描述为：在地下某一深度存在一反射界面R（如图10-14（a）），在地面S点激发的下行波D到达界面R时产生反射上行波U，到达G点被接收。下行波D到达界R面的时间（或空间位置）与上行波U产生的时间（或空间位置）是一致的，即称为时间（或空间位置）一致性。设波从S点到R的传播时间为ts，从R至G的传播时间为tg，从S到G的总时间为tsg=ts+tg。在叠前偏移中，若模拟一震源函数D自S点正向（向下）延拓，而将G点接收到的上行波U反向延拓，当D和U延拓深度为Z1时，D的正向传播时间和U的反向传播时间分别为ts1和tg1。因Z1ts1说明上行波和下行波所在的时间（或空间位置）不一致（如图10-14（b）），当D和U延拓深度为zz=ZR时，下行波正向传播时间为ts1=ts，上行波反向传播时间为tg2=tg，即有tsg-tg2=ts2，或tsg-tg=ts。这时上、下行波所在的时间（或空间位置）是一致的。再将D、U延拓到Z3，Z3>ZR，即当延拓深度Z>ZR以后，不会再出现时间（或深度位置）一致的现象。在上、下行波延拓过程中，若求下行波场D和上行波场U的零移位互相关，在满足时间（或空间位置）一致性条件时，相关值最大，而在其他情况相关值很小或为零，延拓过程中的相关结果就为叠前偏移成像剖面。
图10-14 时间一致性成像原理示意图
10.4.3.2 有限差分法波动方程偏移
有限差分法波动方程偏移是以地面上获得的水平叠加时间剖面作为边界条件，用差分代替微分，对只包含上行波的近似波动方程求解以得到地下界面的真实图像。这也是一个延拓和成像的过程。
（1）延拓方程
由二维波动方程出发，
勘查技术工程学
经数学推导，即可得到下面方程：
勘查技术工程学
式中：uxx，uττ，uτt分别表示波场u（x，z，t）的二次导数。注意，此方程仍然包含了上行波和下行波，仍不能用来进行延拓。
当上行波的传播方向与垂直方向之间的夹角较小时（小于15°），uττ可以忽略，而对下行波来说，uττ不能忽略。忽略掉uττ项，就得到只包含上行波的近似方程
勘查技术工程学
此即15°近似方程（因为它只适用于夹角小于15°的上行波，或者只有倾角小于15°的界面形成的上行波才能满足它），为常用的延拓方程。
为了求解此方程还必须给出定解条件。由于震源强度有限，可给出如下定解条件：
1）测线两端外侧的波场为零，即
勘查技术工程学
2）记录最大时间以外的波场为零，即
勘查技术工程学
3）自激自收记录（水平叠加剖面）为给定的边界条件，即时间深度τ=0处的波场值u（x，0，t）已知。
有了这些定解条件就可对方程（10.4-3）求解得到地下任意深度处的波场值u（x，τ，t），这是延拓过程。再根据前述成像原理，取时间t=0时刻时的波场值，即为时间t=τ时刻的波场值u（x，τ，t）就组成了偏移后的输出剖面。
（2）差分方程
为了求解微分方程（10.4-3），用差分近似微分，采用如图10-15所示的12点差分格式，将uxx、uτt表示为差分表达式，可得差分方程
勘查技术工程学
图10-15 12点差分格式
式中：I和T为向量。
勘查技术工程学
α、β为标量
勘查技术工程学
（3）计算关系和偏移结果
图10-16画出了偏移时的计算关系及结果取值位置。A 表示地面观测到的叠加剖面。由 A 计算下一个深度Δτ处的波场值B，计算 B 时先算第1′排的数值（只用到 A 中第1排数值），再算第2′排数值（要用 A 中第1、2排和 B 中第1′排数值），依此类推，直到 t=τ为止。再由 B 算下一个深度2Δτ处波场值C……在二维空间（x，t=τ）上呈现出需要的结果剖面信息。
图10-16 偏移结果取值位置图
当延拓计算步长Δτ与地震记录的采样间隔Δt一样时，由图10-16的几何关系可以看到，偏移剖面是该图中45°对角线上的值。实际工作中Δτ不一定要与Δt相等，可根据界面倾角大小确定Δτ，倾角较大时应取较小的Δτ，倾角较小时Δτ可取的较大些，以减少计算工作量。中间值可用插值求得。
与其他波动方程偏移方法相比，有限差分法有能适应横向速度变化，偏移噪声小，在剖面信噪比低的情况下也能很好地工作等优点。但15°有限差分法对倾角太大的情况不能得到好的偏移效果。因此，相继又研究发展了45°、60°有限差分偏移方法和适应更大倾角的高阶有限差分分裂算法。
10.4.3.3 克希霍夫积分偏移
克希霍夫积分偏移是一种基于波动方程克希霍夫积分解的偏移方法。
三维纵波波动方程的克希霍夫积分解（可见原理部分）为
勘查技术工程学
式中：Q为包围点（x，y，z）的闭曲面；n为Q的外法线；r为由（x，y，z）点至Q面上各点的距离；〔 〕表示延迟位；〔u〕=u（t-r/v）。
此解的实质是由已知的闭曲面Q上各点波场值计算面内任一点处的波场值。它正是惠更斯原理的严格数学形式。
选择闭曲面Q由一个无限大的平面Q0和一个无限大的半球面Q1所组成。Q1面上各点波场值的面积分对面内一点波场函数的贡献为零。因此，仅由平地面Q0上各点的波场值计算地下各点的波场值
勘查技术工程学
此时，原公式中的项消失，积分号前的负号也因 z 轴正向与n 相反而变为正。
以上是正问题的克希霍夫积分计算公式。偏移处理的是反问题，是将反射界面的各点看作为同时激发上行波的源点，将地面接收点看作为二次震源，将时间“倒退”到t=0时刻，寻找反射界面的源波场函数，从而确定反射界面。反问题也能用上式求解，差别仅在于〔 〕
不再是延迟位而是超前位，〔u〕=u（t+）。根据这种理解，克希霍夫积分延拓公式应为
勘查技术工程学
按照成像原理，此时t=0时刻的波场值即为偏移结果。只考虑二维偏移，忽略掉y坐标，将空间深度z转换为时间深度t0=2z/v，得到克希霍夫积分偏移公式
勘查技术工程学
式中：τ=[]1/2，xl 为地面记录道横坐标，x 为偏移后剖面道横坐标，r=〔z2+（x-xl）2〕1/2（见图10-17）。
由=-cosθ，得
勘查技术工程学
由此可见，克希霍夫积分偏移与绕射扫描叠加十分相似，都是按双曲线取值叠加后放在双曲线顶点处。不同之处在于：①不仅要取各道的幅值，还要取各道的幅值对时间的导数值参加叠加。②各道相应幅值叠加时不是简单相加，而是按（10.4-10）式的加权叠加。
图10-17 克希霍夫偏移公式中各量示意图
正因如此，虽然形式上克希霍夫积分法与绕射扫描叠加类似，但二者有着本质区别。前者的基础是波动方程，可保留波的动力学特性，后者属几何地震学范畴，只保留波的运动学特征。
与其他波动方程偏移法相比，克希霍夫积分法具有容易理解，能适应大倾角地层等优点。它在速度横向变化较大的地区难以使用，且偏移噪声较大。
波动方程偏移
波动方程偏移与绕射扫描叠加偏移相比有本质上的重大改进，是目前实际生产中使用的主要偏移方法，其中又以15°有限差分偏移最为典型。
1.15°度有限差分法波动方程偏移
15°度有限差分法波动方程偏移是以地面上获得的水平叠加时间剖面作为边界条件，用差分代替微分，对只包含上行波的近似波动方程求解得到地下各点波场值，并进而获得地下界面真实图像的一种偏移方法。其偏移的过程也是一个延拓和成像的过程。
1）延拓方程的推导
由下述二维波动方程出发
地震波场与地震勘探
根据爆炸反射面模型，将速度缩小一半，即用v/2代替v，可得：
地震波场与地震勘探
此方程有二个解，分别对应于上行波和下行波。地震记录是单纯的上行波记录，故不能用此方程进行延拓，必须将它化为单纯的上行波方程才能利用。通常采用的方法是进行坐标变换后取近似。第一步是坐标变换，令
地震波场与地震勘探
上式中第一个变换无任何改变；第二个变换只是将空间深度z换成时间深度τ，也无实质性变化。关键是第三个变换，它表示不再用传统的旧时钟计时，而是用一个运行速度与旧钟一样，但起始时刻各深度不同的新时钟计时。采用新时钟计时时，上、下行波就表现出差异。
因为坐标变换并不会改变实际波场，故原坐标系中的波场u（x，z，t）与新坐标系中的波场
是完全一样的，即
地震波场与地震勘探
由复合函数微分法，得：
地震波场与地震勘探
将上述二阶偏微分结果代入方程（4-4-3），整理后得：
地震波场与地震勘探
为书写方便，以u、x、t分别代替
，则（4-4-5）式可写为
地震波场与地震勘探
式中uxx、uττ、uτt分别表示u的二次导数。注意。此方程仍然包含了上行波和下行波，仍不能用来进行延拓，故还有第二步。
经过了坐标变换，虽然波场不变，但在新的坐标系下上、下行波表现出差异，此差异主要表现为uττ的大小不同：当上行波的传播方向与垂直方向之间的夹角较小时（小于15°），uττ可以忽略；而对下行波来说，uττ不能忽略。忽略掉uττ项，就得到只包含上行波的近似方程
地震波场与地震勘探
此即15°上行波近似方程（因为它只适用于运行方向与垂直方向间的夹角小于15°的上行波，或曰只有倾角小于15°的界面形成的上行反射波才能满足它），为常用的延拓方程。
为了求解此方程还必须给出定解条件。由于震源强度有限，可以给出如下定解条件：
a.测线两端外侧的波场为零，即
u（x，τ，t）≡0 当 x＞xmax或 x ＜ xmin时
b.记录最大时间以外的波场为零，即
u（x，τ，t）≡0 当 t＞tmax时
c.自激自收记录（水平叠加剖面）为给定的边界条件，即时间深度τ＝0处的波场值u（x，0，t）已知。
图4-4-5 12点差分格式
有了这些定解条件就可以对方程（4-4-7）式求解得到地下任意深度处的波场值u（x，τ，t），这是延拓过程。再根据前述成像原则，取传统旧时钟零时刻时的波场值，即新时钟时间t＝τ时刻的波场值u（x，τ，τ）就组成了偏移后的输出剖面。
2）差分方程的建立
为了求解微分方程（4-4-7）式，用差分近似微分，采用如图4-4-5所示的12点差分格式，可得：
地震波场与地震勘探
将（4-4-8）式和（4-4-9）式代入（4-4-7）式中得：
地震波场与地震勘探
地震波场与地震勘探
定义向量I、T：
I＝［0，1，0］ T＝［-1，2，-1］
令向量u （x，j，l）为
u（x，j，l）＝［u（i-1，j，l），u（i，j，l），u（i＋1，j，l）］
则（4-4-10）式可简写为
地震波场与地震勘探
又令
地震波场与地震勘探
则（4-4-11）式可写成如下形式：
［I-（α＋β）T］u（x，j＋1，l＋1）-［I＋（α-β）T］u（x，j，l＋1）＋［I-（α＋β）T］u（x，j，l）＝［I＋（α-β）T］u（x，j＋1，l）
因此有：
地震波场与地震勘探
此即适合计算机计算的差分方程。
3）计算步骤和偏移结果
差分方程（4-4-12）形式上是一个隐式方程，即时间深度τ＝（j＋1）Δτ处的波场值不能单独地用时间深度τ＝jΔτ处的波场值组合得到，方程右边仍然有τ＝（j＋1）Δτ的项。如图4-4-6所示，为了求得一排数据u（x，j＋1，l），必须用到三排数据u（x，j＋1，l＋1），u（x，j，l＋1）和u（x，j，l）。一般来说，隐式方程的求解必须用求解联立方程的方法进行，比较麻烦，但这里可以利用有利的定解条件，无须复杂的联立运算。
利用定解条件b，在计算新的深度τ＝（j＋1）Δτ处的波场值时，由最大时间开始，首先计算t＝tmax的那一排值。因u（x，j＋1，tmax＋Δt）≡0和u（x，j，tmax＋Δt）≡0，有：
地震波场与地震勘探
计算u（x，j＋1，tmax）只用到已知的u（x，j，tmax）值，十分容易。然后再利用（4-4-12）式递推地求τ＝（j＋1）Δτ深度处任何时刻的波场值就没有任何困难了。
具体计算时由地面向下延拓，计算深度Δτ处的波场值：首先计算此深度处在t＝tmax时的波场，然后向t减小的方向进行计算直至本深度处的全部波场值计算完。一个深度的波场值计算结束后，再向下延拓一个步长Δτ继续计算。依此类推，可以得到地下所有点在不同时刻的波场值。
如前所述，在新时钟t＝τ时刻的波场值正是所欲求的“像”。因此，每次递推计算某一深度τ处的波场值时，由t＝tmax向t减小的方向计算至t＝τ时就可以结束了，u （x，τ，τ）为该深度处的“像”。不同深度处的“像”组成偏移后的输出剖面。
图4-4-6 有限差分法偏移求解中的一步
①u（x，j，l＋1），②u（x，j，l），③u（x，j＋1，l＋1），④u（x，j＋1，l）
图4-4-7 偏移结果取值位置图
图4-4-7画出了偏移时的计算关系及结果取值位置。A表示地面观测到的叠加剖面，由A计算下一个深度Δτ处的波场值B，计算B时先算第1′排的数值（只用到A中第1排数值），再算第2′排数值（要用到A中第1、2排和B中第1′排的数值），依此类推进行计算，直到算出t＝Δτ的值为止，再由B计算下一个深度2Δτ处的波场值C，到算出t＝2Δτ的值为止……在二维空间（x，t＝τ）上呈现出需要的结果剖面信息。
当延拓计算步长Δτ与地震记录的采样间隔Δt一样时，由图4-4-7的几何关系可以看到，偏移剖面是该图中45°对角线上的值。实际工作中Δτ不一定要与Δt相等，应当根据界面倾角大小确定Δτ，倾角较大时应取较小的Δτ，倾角较小时Δτ可取得大一些，以减少计算工作量，中间值用插值方法求得。
与其他波动方程偏移方法相比，有限差分法有能适应横向速度变化、偏移噪声小、在剖面信噪比低的情况下也能很好地工作等优点，但15°有限差分法在界面倾角太大时不能得到好的偏移效果。因此，又发展了45°、60°甚至90°的有限差分偏移方法，有兴趣的读者可参阅有关文献。
2.频率波数域波动方程偏移
有限差分偏移方法是在时间空间域中进行计算的。利用傅里叶变换也可以使偏移在频率波数域中实现。
与有限差分法偏移的思想完全一样，认为水平叠加剖面是由界面上无数震源同时向上发出的上行波在地面处的波场值u（x，0，t），用它反求地下任一点的波场值u（x，z，t）是延拓过程。再根据成像原理，取其在t＝0时刻的值u（x，z，0），就组成了偏移后的输出剖面。
仍由速度减半后的波动方程（4-4-3）出发，对方程两边作关于x和t的二维傅里叶变换，得到一个常微分方程：
地震波场与地震勘探
式中：U＝U（kx，z，ω）是波场函数u（x，z，t）的二维傅里叶变换，ω＝2πf为角频率，kx为x方向上的空间波数。
（4-4-13）式是常微分方程，很容易求解。其解有二个，分别对应于上行波和下行波。偏移研究的是上行波的向下延拓问题，故只考虑上行波解
地震波场与地震勘探
其中U（kx，0，ω）是解的初值，即上行波在地面（z＝0）处记录的傅里叶变换。因此，式（4-4-14）表示由z＝0 处波场的傅里叶变换求出地下任何深度处波场傅里叶变换的过程，是频率波数域中的波场延拓。
通过傅里叶反变换可以由U（kx，z，ω）求出地下任何深度处的波场值：
地震波场与地震勘探
根据成像原理，偏移结果应是该深度处t＝0时刻的波场值：
地震波场与地震勘探
这就是频率波数域偏移的数学模型。其具体实现步骤就不赘述了。
如果求解常微分方程（4-4-13）时初值不取z＝0处波场值的傅里叶变换，而取任一较浅处的波场傅里叶变换值，则可得到：
地震波场与地震勘探
从而得到相移法偏移的数学模型：
地震波场与地震勘探
利用此式可逐步向下延拓成像，每延拓一次所用的速度均可改变，所以相移法能够适应速度的纵向变化。
由于快速傅里叶变换的应用，频率波数域偏移法效率十分高，运行时间少，是波动方程偏移算法中最经济的方法，且适用于大倾角地区。因为计算在频率波数域中进行，需要注意假频问题，且此法对横向速度变化的地区不太适应。
3.克希霍夫积分偏移
克希霍夫积分偏移是一种基于波动方程克希霍夫积分解的偏移方法。
三维纵波波动方程的克希霍夫积分解（见第一章）为
地震波场与地震勘探
式中Q为包围点（x，y，z）的闭曲面，n为Q的外法线，r为由（x，y，z）点至Q面上各点的距离，［ ］表示延迟位，
。
此解的实质是由已知的闭曲面Q上各点波场值计算面内任一点处的波场值，它正是惠更斯原理的严格数学形式。
选择闭曲面Q由一个无限大的平地面Q0 和一个无限大的半球面Q1 所组成。Q1 面上各点波场值的面积分对面内一点波场函数的贡献为零，因此仅由平地面Q0 上各点的波场值计算地下各点的波场值。在此条件下，地下任一点的波场值为
地震波场与地震勘探
此时，原公式中的
项消失，积分号前的负号也因z轴正向与n相反而变为正。
已知源函数，求取波传到某点的波场值是正问题。以上是正问题的克希霍夫积分计算公式。偏移处理的是反问题，是将地面接收到的波场值看作为二次震源，将时间“倒退”寻找地下波场值，取t＝0 时刻的波场值确定反射界面的问题。反问题也能用上式求解，差别仅在于［ ］不再是延迟位而是超前位，
。根据这种理解，克希霍夫积分延拓公式应为
地震波场与地震勘探
按照成像原理，t＝0时刻的波场值即为偏移结果。只考虑二维偏移，忽略掉y坐标，将空间深度z转换为时间深度t0＝2z /v，得到克希霍夫积分偏移公式
地震波场与地震勘探
式中：
；xl为地面记录道横坐标；x为偏移后剖面道横坐标（见图4-4-8）。
图4-4-8 克希霍夫偏移公式中各量示意图
由
，得：
地震波场与地震勘探
由此可见，克希霍夫积分偏移与绕射扫描叠加十分相似，都是按照绕射双曲线取值叠加后放在双曲线顶点处。不同之处在于：
a.不仅要取各道的幅值，还要取各道幅值对时间的导数值
参加叠加；
b.各道相应幅值叠加时不是简单的相加，而是按（4-4-22）式的加权叠加。
正因如此，所以虽然形式上克希霍夫积分偏移法与绕射扫描叠加偏移类似，但二者有着本质的区别。前者的基础是波动方程，可保留波的动力学特性；后者属几何地震学范畴，只保留波的运动学特征。
与其他波动方程偏移法相比，克希霍夫积分法具有容易理解，能适应大倾角地层等优点。它在速度横向变化较大的地区难以使用，且偏移噪声较大。
地震波实际上是在三维空间中传播的，故要实现完全的偏移必须是三维偏移。目前三维偏移方法已经得到了极大的发展，从二步法到分裂法，到目前已经实现了没有近似的一步法完全三维偏移。
上面介绍的叠后偏移有一个基本假设，即水平叠加剖面是自激自收剖面，实际在地下界面较复杂时这一假设是不成立的。为了实现真正的叠加和偏移，发展了叠前偏移方法，它将偏移与叠加同时进行，保证了能达到真正的共反射点叠加。但是，一次就完成偏移和叠加的任务对于求取速度参数是不利的。为此又发展了叠前部分偏移方法，它只进行部分偏移，使共中心点道集成为真正的共反射点道集，以保证实现共反射点叠加。叠前部分偏移方法加叠后偏移就等于叠前偏移。有了共反射点道集，就好进行速度分析。
目前大部分偏移还是时间偏移，即偏移后得到的是时间剖面。时间偏移没有考虑地震波穿过界面时的弯折现象。如果考虑这一现象，偏移后得到的就是深度剖面，这种偏移称为深度偏移。目前，三维叠前深度偏移是最先进的偏移方法。
地震波实际上是弹性波。目前的偏移方法均使用声波方程，它只是一种声波偏移，要实现真正完全的偏移应当是弹性波偏移。因为需要多波地震资料，目前弹性波偏移还很少使用。
基于二维地质建模的两种地震数值模拟方法的应用及分析
赵忠泉
（广州海洋地质调查局 广州 510760）
作者简介：赵忠泉，男，（1983—），硕士，主要从事海洋油气资源调查研究工作，E-mail：zzqhello＠163.com。
摘要 利用地震数值模拟技术结合实际资料，可以建立各种地质体的地震识别模型，有效地避免地震现象的多解性，从而可以提高解释的精度。本文介绍了二维地质建模的方法流程及两种模拟方法-褶积法和PSPI波动方程法，前者无边界条件约束和频率域中的信号损失，简洁易行，计算稳定，应用广泛，是最早的地震波场模拟方法；后者通过求解波动方程，包含丰富的波场信息，能够充分反映地震波的动力学和运动学特征。实际应用中利用褶积法对三维潮道模型及简化的碳酸盐岩多旋回倾斜薄互层沉积模型进行了模拟；利用零炮检距的频率波数域的波动方程法模拟了生物礁的地震响应，结果对于碳酸盐岩生物礁识别有一定指导意义。
关键词 地质建模 数值模拟 褶积法PSPI法
不同地质体由于其岩性、物性、含油气性、内部结构和岩石组合等的差异，在地震上具有不同的反射特征，包括内部结构、外部形态、振幅、频率等参数。由于地震波在地下地质体中传播的复杂性，加上各种干扰，造成了地震剖面中的各种反射现象存在多解性，大大增加了地震解释的难度。利用地震数值模拟技术结合实际资料，在建立不同地质体的地震识别模型的同时也有效地避免地震现象的多解性，从而可以提高解释的精度。
1 地质建模
地震数值模拟技术的基础是地质地球物理模型的建立，可归结为对地质及地球物理模型结构的数学描述。
二维封闭结构模型用于建立复杂地质模型。二维封闭结构模型就是定义相同地质属性为一独立封闭的地质单元，按照地质属性将地质模型划分成多个独立封闭的地质单元，把所有独立封闭地质单元按照空间分别有序地排列起来，这样组成的集合体就构建了一个二维地质模型。封闭结构模型是以积木方式定义地下地质结构，可以描述非常复杂的地质体。二维封闭结构模型被描述为具有相同地质属性（速度、密度等），并被地层界面、断层界面或模型边界所围成的地质单元的有机组合。对封闭结构模型的描述，实际上就是描述封闭地质单元和封闭地质单元之间的关系，前者包括对封闭地质单元属性和封闭地质单元边界的描述；后者是对地质单元空间关系的描述，也就是描述封闭地质单元边界相接关系及地层属性［1］。
在进行数值模拟过程中，为了验证某些复杂地质体的波场特征，需要绘制多种不同的地质模型，通常可借助常规绘图软件（绘图板、Photoshop、CorelDraw，AutoCAD等）绘制好二维封闭结构面，再根据图像处理中的区域填充算法（种子填充和扫描转换填充），对不同二维封闭结构面进行不同颜色的填充。其中不同颜色代表不同的二维封闭结构面属性（速度、密度等）；合并相同属性的封闭面，形成最终的二维封闭结构模型［1］。为了得到二维封闭结构模型的属性（速度、密度等）模型，需要对二维封闭结构模型的彩色图进行速度像素空间和属性空间转换，根据颜色空间和属性空间的相互映射，就可以得到复杂地质体的属性（速度、密度等）模型，如图1为模型创建流程图。
图1 二维封闭结构模型建立流程图
2 两种数值模拟方法
2.1 褶积模型
在褶积模型中，我们把地震反射信号s（t）看作是地震子波w（t）与地下反射率r（t）的褶积。地震子波w（t），使用实际地震系统记录到的地下一个单独的反射界面反射的波形（如图2，理想的无噪声褶积过程）。反射率r（t）则代表理想的无噪声地震记录。记录到的地震道s（t）可看作是地震信号w（t）* r（t）与可加噪声n（t）之和，因此可以把地震道看作是一种有噪声干扰的，经过了滤波的地下反射率的变形。
在无噪声褶积模型中，我们把地震信号S（t）看作是地震子波w（t）和地下反射系数r（t）的褶积：
南海地质研究.2012
式中：s（t）——合成地震记录；
r（t）——反射系数；
w（t）——地震子波。
图2 褶积过程
2.2 PSPI波动方程法
通过求解波动方程的数值模拟方法，能够充分反映地震波的动力学和运动学特征，波场信息丰富，模拟结果较为准确。这里仅介绍适合横向速度剧烈变化的频率-波数域相移加插值的波场延拓方法［2］。
相位移加插值的波场延拓方法，简称PSPI法，基本思想是在波场向下延拓的每个深度步长Δz之内，将波场的延拓分成两部进行，首先用L个参考速度V1，V2，…VL，将位于深度zi处的波场p（x，zi，ω）延拓到zi＋1＝zi＋Δz处，得到L个参考波场p1（x，zi＋1，ω），p2（x，zi＋1，ω），…，PL（x，zi＋1，ω）。第二步，按实际的偏移速度V（x，z）同参考速度V1，V2…，VL的关系，用波场插值的方法求出zi＋1处的波场p（x，zi＋1，ω），按同样的步骤，可将zi＋1处的波场值p（x，zi＋1，ω）延拓到深度zi＋2，得p（x，zi＋2，ω），直到延拓到最大的深度zmax为止。
对于各向同性介质，取二维标量声波方程作为延拓的基本方程：
南海地质研究.2012
式中，p＝p（x，z，t）为二维地震波场值；x，z分别为水平方向和垂直方向坐标轴；t为时间轴；v（x，z）为纵、横向都可变的地震波传播速度。将式（2）分别对x、t作傅氏变换，考虑到并考虑到ə2/əx2和与（-ikx）2和（iw）2的对应关系，可得：
南海地质研究.2012
式中， 是p（x，z，t）的二维傅氏变换；v为地震波速度；w为圆频率；kx为水平波数；kz为垂直波数。零炮检距情况下的地震记录模拟只考虑单程波，因此可得到相位移波场延拓公式如下：
南海地质研究.2012
式中， （kx，zi，w）为频率波数域波场值；Δz为深度延拓步长；kx为测线方向波数；kz为深度方向波数。式（4）为二维波场正演公式，其延拓方向为由地下向地面延拓；式（5）为二维波场偏移公式，其延拓方向为由地面向地下延拓。
为了适应地下地震波场速度在纵横向均可变的要求，在同一延拓深度内用几个不同地震波速度分别作相移，再用拉格朗日插值公式进行插值，就可求出所有的以不同速度传播的延拓波场值P（x，zi＋1，t），从而近似地解决了横向变速时的波场延拓问题［3］。
3 模拟实例
3.1 三维潮道数值模拟
运用褶积原理建立了一个简单三维潮道模型，此三维潮道事实上为多个（128）二维剖面排列而成，三维模型的采样点为128×128×128，利用MATLAB实现。选用子波为雷克（Ricker）子波，其公式为：
南海地质研究.2012
其中fp为主频。在处理过程中选用主频为fp＝40 Hz、采样间隔2 ms，对称采样点数为24，子波波形如图3。
图3 雷克子波
图4 潮道平面图
图4为潮道平面图，该图仅反映了潮道的平面形态，作为计算机实现三维建模的边界控制，横坐标代表inline线，纵坐标代表xline（crossline）线，图5为三维地质模型示意图，模型较简单，整体由三个水平层叠置而成，在第二层和第三层之间镶嵌了形如图4的潮道，此潮道没有考虑进水方向，根据此地质模型进行计算机地震正演模拟，可得到相应三维地震数据体，从图中可以看到，黄色虚线（上）和蓝色虚线（下）位置上，分别横跨了三个潮道分支和两个潮道分支，就是说在相应两条虚线位置上的两条测线应该分别有三个和两个潮道显示，提取相应的两条剖面如下图6和图7：
图5 三维地质模型
图6 xline＝100（黄色虚线）剖面
图7 xline＝100（蓝色虚线）剖面
再在三维数据体中沿水平方向做切片，即提取时间切片。图8为时间切片在地震剖面上的位置示意图，图中五条标示线从上到下依次为白色实线、黄色虚线、白色实线、红色虚线和白色实线，与之对应的时间分别为70 ms、85 ms、95 ms、99 ms和110 ms（时间范围是0～128 ms），图9～图13为相应切片，从图中可以看出，随着所做切片时间的增大（深度的增加），潮道的展布范围逐渐减小，由于地层是水平层状的，使得时间切片等同于地层切片和沿层切片，其切片效果非常明显，切片中潮道形态得到了很好的展示，但是在多个切片中发现，从可以见到潮道形态一直到潮道消失的时间范围是在70～110 ms之间，而潮道的真实范围是在80～100 ms之间，显然依据切片所圈定的潮道的范围相比真实的范围扩大了，究其原因是由于不管选取哪一种子波，子波都有一定的延续长度和有限频宽，这就限制了合成地震记录本身的分辨率并不能达到等时厚度反射系数序列的分辨率。因此在对实际地震资料进行解释的时候，对地质异常体边界的识别应该考虑地震子波并非脉冲波所带来的影响。
图8 剖面示意图
图9 切片t＝70 ms
图10 切片t＝85 ms
图11 切片t＝95 ms
图12 切片t＝99 ms
图13 切片t＝110 ms
3.2 薄互层沉积模型
图14为简化的碳酸盐岩多旋回倾斜薄互层沉积模型（Zeng，2003），模型简化是为了更好地突出由岩相控制的波阻抗结构和地震信号之间的相互关系。该模型所有倾斜的倾角都相同，每层都有相同的垂直时间厚度（5 ms或15 m，速度为6000 m/s），泥岩与低孔隙度颗粒灰岩的波阻抗差，以及低孔隙度颗粒灰岩与高孔隙度颗粒灰岩的波阻抗差都相同，所有高孔隙度颗粒灰岩具有相同的深度范围，综合起来形成了一个水平的岩性地层单元。
其时间域地震响应（图15）中，高频情况下（60 Hz雷克子波），地震反射被建设性地调谐到时间地层单元，因此地震同相轴沿着时间地层单元分布（图15a）。当子波频率减到40 Hz时，地震反射对时间地层单元和岩性地层单元都有响应（图15b）。当用30 Hz雷克子波时（图15c），地震同相轴破坏性地调谐到时间地层单元和建设性地调谐到岩性地层单元，因而时间地层单元的反射进一步变弱，地震同相轴被岩相反射所控制［4］。
这个模拟过程强调了了解地质格架和时间地层单元以及岩性地层相带厚度尺度的重要性。时间地层（图15a）和岩性地层（图15c）成像都是有用的，前者用于对比，后者用于粗略的储层评价。然而，这两种响应不能混淆在一起。图15b中的两组相互矛盾的地震同相轴会造成地震假象［4］。
图14 简化的碳酸盐岩多旋回倾斜薄互层沉积模型
3.3 生物礁数值模拟［5～7］
频率—波数域的相移加插值偏移（PSPI）在每一个深度间隔内使用多个参考速度进行偏移，由多个偏移结果插值生成最终的偏移剖面，所用插值的速度越多，越能反映实际介质的速度变化情况，此方法在成像精度及横向变速适应性上具有很大的优越性，但处理所需的时间稍长，鉴于本文的二维叠后建模对处理时间没有过高要求，因此应用PSPI方法做正演、偏移。
图16为某区块过生物礁的原始地震剖面，图17为根据此剖面建立的生物礁速度模型：模型速度变化范围是5600 m/s到5980 m/s，从图16中可以看出生物礁的底界面清晰可辨，围岩有披覆现象，内部呈杂乱反射。为了检验该地质建模的正确性，先采用PSPI方法对该模型进行了波场正演模拟计算，其模拟剖面如图18所示。由于生物礁埋藏深，生物礁顶底反射的弧度较大，不规则点的绕射波杂乱，因此用图15的速度模型对其进行叠后时间偏移，得到了偏移剖面（图19），横向表示256个地震道，纵向表示零偏移距反射时间，礁体最大时间厚度约40 ms。从图19可以看出，模拟记录中的礁体顶界与原始剖面有一定差距，但是生物礁底界反射和内幕反射以及侧翼反射与原始剖面基本一致，其他的地层界面形态与原始剖面也吻合较好，在一定程度上验证了地质模型的正确性，说明当生物礁与围岩之间存在一定波阻抗差异时，在地震剖面上必然出现异常反映，经过有效的构造和参数反演，能够将其分辨出来。相信通过模型改进以及算法中参数的调整，能够与原始剖面更好地吻合，从而为生物礁的地震解释提供一种有力的验证工具。
图15 图14模型时间域地震响应
4 结论
地震数值模拟（正演）技术基于地球物理模型的建立，运用概念二维封闭结构地质模型的建立方法，得到复杂地质体的数学模型，结合各种算法对其进行模拟从而可以验证相应地质体的地震波场特征；结合实际资料建立不同地质体的地震识别模型，可以有效地减少地震现象的多解性，从而提高解释的精度；褶积法无边界条件约束和频率域中的信号损失，简洁易行，计算稳定，应用广泛，本文用此方法模拟的伪三维潮道模型及倾斜薄互层模型取得了较好的效果；通过求解波动方程的数值模拟方法，包含丰富的波场信息，能够充分反映地震波的动力学和运动学特征，PSPI波场沿拓方法为其中之一，利用正演与偏移相结合的流程模拟了生物礁的地震响应特征，检验解释成果的正确性，为生物礁的地震解释提供了一种有力的检验工具。
图16 原始剖面
图17 生物礁地质速度模型（256×256）
图18 正演记录（子波主频30Hz）
图19 偏移剖面（子波主频30Hz）
参考文献
［1］刘远志.碳酸盐岩地震相分析与数值模拟［D］.成都：成都理工大学，2009.
［2］韩建彦.复杂地质体地震正演与偏移［D］.成都：成都理工大学，2008.
［3］贺振华，王才经等.反射地震资料偏移处理与反演方法［M］.重庆：重庆大学出版社，1989.
［4］Zeng Hongliu ＆Kerans，C.Seismic frequency control on carbonate seismic stratigraphy；a case study of the Kingdom Abosequence，West Texas，American Association of Petroleum Geologists Bulletin，2003.87，273～293.
［5］贺振华，黄德济，文晓涛，等.碳酸盐岩礁滩储层多尺度高精度地震识别技术［R］.成都：成都理工大学地球探测与信息技术教育部重点实验室，2009.
［6］熊晓军，贺振华，黄德济.生物礁地震响应特征的数值模拟［J］.石油学报，2009，30（1）：7～65.
［7］熊忠，贺振华，黄德济.生物礁储层的地震数值模拟与响应特征分析［J］.石油天然气学报，2008，30（1）：75～78
The application and analysis of two kinds of seismic numerical simulation method based on the2D-geological modeling
Zhao Zhongquan
（Guangzhou Marine Geological Survey，Guangzhou，510760）
Abstract：Pick to using seismic numerical simulation technology combined with the actual seismicdata，we can build all kinds of seismic recognition model of geologic body and effectively avoidthe multiple solutions of seismic phenomenon，which can improve the precision of the explana-tion.This paper describes the method of the process of 2D geological modeling and two simulationmethods，seismic convolution method and PSPI wave equation method，the former has no bounda-ry condition and the signal loss in frequency domain，is concise and easy，it can be calculatedsteadily and be applied widely，is the earliest simulation method in seismic wave field，the latterbased on the wave equation，it contains the rich information in wave field，can fully reflect thedynamics and kinematics characteristics of seismic wave.In the practical application，we use theconvolution model in 3D-tidal channel model and the multi-cyclic simplified deposition model oftilt thin interbed layer of carbonate；We simulate the seismic response of reefs using the method ofzero-offset wave equation in frequency and wave number domain，it is confirmed that the resulthas definite significance for the identify of the reef.
Key words：Geological modeling Numerical simulation Convolution PSPI method
