引用本文
陈飞旭, 李振春, 张凯, 孙琦, 尚江伟. 基于波动方程偏移的宽方位三维角道集提取[J]. 物探与化探, 2015,39(4): 797-804.
CHEN Fei-Xu, LI Zhen-Chun, ZHANG Kai, SUN Qi, SHANG Jiang-Wei. Wide-azimuth 3D angle gathers based on wave equation migration[J]. Geophysical and Geochemical Exploration, 2015,39(4): 797-804.  
Doi:10.11720/wtyht.2015.4.23
Permissions
Copyright©2015, 物探与化探编辑部
? 《物探与化探》编辑部
基于波动方程偏移的宽方位三维角道集提取
陈飞旭1, 李振春1, 张凯1, 孙琦2, 尚江伟2
1.中国石油大学(华东) 地球科学与技术学院,山东 青岛 266580
2.中国石油塔里木油田分公司勘探开发研究院,新疆 库尔勒 841000

作者简介: 陈飞旭(1988-),男,现为中国石油大学(华东)地球科学与技术学院地球物理系在读硕士研究生,主要从事地震波偏移成像以及速度分析方面的研究工作。

收稿日期: 2014-08-13
基金: 国家自然科学基金项目(41204086、41374122)
摘要

基于三维波动方程偏移的宽方位三维角度道集不仅包含了局部反射角度的信息,还包含反射平面在地下的方位角信息,它充分反映了宽方位采集观测系统下勘探目标不同方位的照明情况。这一特点在复杂构造的三维偏移速度分析中是非常重要的。笔者从角度分解的思想出发,讨论了射线参数、反射角以及局部方位角的关系,实现了宽方位情况下局部偏移距域共成像点道集向角度域共成像点道集的映射方法,并给出具体流程。模型测试和实际资料应用结果表明,用该方法提取的角度道集可以真正反映局部反射处的不同方位角的反射信息,可以为全三维勘探地震成像中的深度偏移速度分析研究提供道集依据,也可以用于三维的振幅随角度关系分析。

关键词: 波动方程偏移; 宽方位; 三维角道集提取; 偏移速度分析
中图分类号:P631.4 文献标志码:A 文章编号:1000-8918(2015)04-0797-08
Wide-azimuth 3D angle gathers based on wave equation migration
CHEN Fei-Xu1, LI Zhen-Chun1, ZHANG Kai1, SUN Qi2, SHANG Jiang-Wei2
1.School of Geosciences,China University of Petroleum(East China),Qingdao 266580,China
2.Research Institution of Petroleum Exploration and Development,Tarim Oilfield,PetroChina,Korla 841000,China
Abstract

Wide-azimuth angle gathers based on 3D wave equation migration contain not only the information of local reflection angle but also that of subsurface azimuth angle.It fully reflects the illumination of the target from different azimuths under the wide-azimuth acquisition,which makes itself important in 3D migration velocity analysis in complex structure.Starting with the viewpoint of angle decomposition,this paper discussed the relationship of ray parameter,reflection angle and local azimuth angle,implemented the mapping of the local offset-domain common image gather to the angle-domain common image gather,and gave the work flow of the method.Synthetic model tests and field data application demonstrate that the angle gathers extracted with the proposed method can truly reflect the reflection information from different azimuths in the subsurface,provide gather basis for 3D seismic pre-stack depth migration, and thus can be used to analyze amplitude versus angle in 3D.

Keyword: wave equation migration; wide-azimuth; 3D angle gather extraction; migration velocity analysis

在速度强烈变化的复杂地质构造区域, 为了能得到精确的成像结果, 采用能适应强横向变速的偏移成像方法是必要的。目前主要的偏移方法有基于射线理论的Kirchhoff积分法[1, 2]和高斯波束法[3, 4]; 另一类是基于波动方程的偏移方法(如有限差分法逆时偏移、Fourier偏移方法等)。射线类方法计算效率高且灵活, 但是这类方法也有明显缺陷, 例如射线追踪前需要对速度场进行平滑, 在速度分布过于复杂的区域, 会出现焦散或阴影区, 这时计算出的旅行时场也就不准确。由于波动方程不存在这些困难, 因此, 近年来人们对波动方程法叠前深度偏移进行了大量的研究并做了部分应用。同时, 并行巨型机和高性能多节点集群域GPU集群的出现以及它们在地震数据处理中的应用, 为波动方程法叠前深度偏移提供了硬件条件和高的计算效率[5]。

有限差分逆时偏移方法[6]是典型的波动方程有限差分法, 它借助于差分计算, 把速度、密度等介质参数的影响体现在差分计算的矩阵中, 因而能自动适应速度场的任意变化。该方法基于双程波动方程, 将震源正传波场和地表记录反传的接收波场作零延时的互相关得到成像结果, 不存在地层倾角的限制, 但是计算量大, 并对计算机内存要求较高。

基于波场延拓的波动方程叠前深度偏移[7, 8, 9, 10]也是一种比较精确、稳定的成像方法。这类方法的思路是进行速度场分裂, 把复杂的介质速度场分裂为“ 常速背景+层内变速扰动” , 然后针对分裂后的速度场分别进行延拓处理。其中, 基于“ 观测沉降” 延拓法的双平方根(DSR)偏移只对上行波场进行延拓, 相当于向地下延拓(沉降)震源与接收点, 当二者重合时(零偏移距), 零时间的波场值就作为该空间点的成像值。相对于炮集偏移而言, DSR偏移的计算效率更高, 便于为速度估计、振幅与偏移距关系(AVO)振幅与角度关系(AVA)分析输出偏移距或角度域的共成像道集, 尤其在海上拖缆观测地震数据的成像中, DSR偏移的优势更明显[11]。将二维DSR扩展到三维, 容易得到全三维DSR偏移算子, 该算子在偏移过程中不对crossline作近似, 在成像道集中完整保留了方位角的信息, 因此全三维DSR偏移算子特别适合宽方位数据的偏移处理。

叠前深度偏移效果的优劣不仅与偏移方法有关系, 更重要的是还与偏移速度模型有关。对于波动方程叠前深度偏移, 由于偏移对速度模型很敏感, 所以速度模型对偏移成像结果有决定性的作用。反过来讲, 可以利用波动方程偏移对速度模型的敏感性来进行速度分析, 以成像质量最佳为速度场正确的判断准则。

偏移速度分析主要是利用偏移产生的成像道集对速度误差的敏感性不断迭代更新速度和成像结果, 因而道集的提取方法是速度分析的重要环节。目前, 角度域共成点像道集(ADCIG)是应用最广泛的共成像点道集, 它反映了偏移深度与反射角的关系, 速度不正确时会表现一定的剩余曲率, 且这种道集不受多路径对成像的影响, 最适合于偏移速度分析、层析速度反演和AVA分析, 因此许多著名学者对三维角度道集的提取方法作了深入的研究[12, 13, 14, 15]。二维情况下, ADCIG可以由局部偏移距共成像点道集(ODCIG)通过倾斜叠加映射得到。三维情况下, 尤其在地下构造复杂的地区, 射线很容易发生扭曲, 导致局部射线平面的法线发生改变, 即射线平面的方位角发生改变。目前三维角度道集的提取主要是采用方位角校正后的共方位角偏移方法[16, 17, 18], 然后提取三维角道集, 其方位角表征的是射线平面在地表的方位角, 而非局部反射处的方位角, 这种方法没有考虑由于地下构造倾角的信息引起的射线平面扭曲, 因此不能完全反映复杂构造处的偏移速度信息。

笔者从角度分解的思想出发, 讨论了射线参数、反射角以及局部方位角的关系, 实现了宽方位情况下局部偏移道集向局部角度道集映射, 并给出具体流程。通过模型测试和实际资料应用加以验证, 说明该方法提取的角度道集能真正反映局部反射处的不同方位角反射信息。

1 方法原理
1.1 扩展局部偏移距的成像条件

一般地, 偏移的实现过程就是对炮、检点波场进行延拓, 然后应用互相关成像条件[19]

I(x)=∫S* (x, ω)R(x, ω)dω, (1)

式中:I(x)表示成像结果, S* 为震源波场S的复共轭, R是检波点的反向延拓波场, ω 为圆频率。

形如式(1)的成像条件可以称为传统的互相关成像条件, 偏移结果I(x)只是成像点x的函数。在偏移速度分析中, 为了保留与局部偏移距(或局部反射角度)有关的振幅信息, 成像时对偏移距轴进行扩展, 即应用扩展成像条件[20]

I(x, h)=∫S* (x-h, ω)R(x+h, ω)dω2

得到成像道集I(x, h), 称为局部偏移距共成像点道集。

从式(2)可知, I(x, h)是I(x)在局部偏移距h维度的扩展, 其包含的速度信息更丰富, 因此可以作为判断偏移正确与否的依据。当速度正确时, 道集能量向零偏移距(h=0)处聚焦, 速度不正确时, 道集向非零偏移距处发散。

在二维情况下或者三维均匀介质下, 射线平面始终在同一个垂直平面内, 如图1a所示。因此通过倾斜叠加可以将ODCIG映射到角度域, 得到ADCIG。当速度正确时, ADCIG是拉平的; 当速度不正确时, 道集具有曲率。ADCIG的拉平与否是判别偏移速度是否准确的标准, 偏移速度分析的过程就是不断更新速度, 使弯曲的ADCIG最终得到拉平[21]。

图1 射线路径在不同介质下的示意

在复杂介质情况下, 射线平面不再保持在同一个平面内, 如图1b所示。地下的盐丘、断层等复杂构造使得射线平面很容易发生扭曲, 地表的炮检点方位角已经不能用来表示地下反射点处的射线平面方位, 因此需要综合考虑反射点处的反射角和方位角的几何关系, 才能正确描述成像深度随反射角和方位角的变化规律。

1.2 ODCIG向ADCIG的映射

定义震源坐标s(sx, sy, sz)和检波点坐标r(rx, ry, rz), 对应的中心点坐标o(x, y, z)和半偏移距坐标h(hx, hy, hz), 有

s=o-h, (3)r=o+h。(4)

定义k和p分别为波数和射线参数, 则有

po=0.5(ps+pr), (5)ph=0.5(pr-ps)。(6)

图2为局部反射处的射线参数几何关系, 从中容易得出po与ph相互垂直, 则可以得到反射角γ 的表达式(Sava和Fomel[22])

|tanγ=|ph||po|=kh/ωko/ω=|kh||ko|=　　　khx2+khy2+khz2kx2+ky2+kz2, (7)

式中:ko和kh分别是中心点波数和局部偏移距波数, ω 为圆频率。

图2 局部反射处的射线几何关系

从图2中可以看出, 方位角θ 可以由测线方向单位矢量d0=(1, 0, 0)和与反射平面正交的矢量d=(dx, dy, dz)决定。d和d0之间的夹角表征了反射平面的法线相对于测线方向的旋转角度。这个关系可以写成

d=po×ph, (8)cosθ=d0·d|d0|d|。(9)

定义dh=(dx, dy, dz)为矢量d在x-y平面的投影, 则有

cosθ=khykx-khzky(khykz-khzky)2+(khzkx-khxkz)2。(10)

由po和ph相互垂直的关系可知, ko与kh也相互垂直, 则有

kh·ko=khxkx+khyky+khzkz, (11)khz=-khxkx-khykykz。(12)

由于全部成像点的偏移距道集I(x, y, hx, hy, z)是一个五维数组, 内存消耗非常大, 所以一般选取其中少量的中心点(x, y)应用扩展成像条件得到ODCIG, 因此傅里叶变换后的道集中未包含参数kx和ky, 可以通过引入局部倾角估计来化简方程。

定义 分别为沿着x和y方向的局部倾角, 可在成像剖面上利用平面波分解法估计得到。则式(12)可重新写成

khz=khxαx-khyαy。(13)

将式(13)代入式(7)和式(10)中, 就可以得到如下的最终映射公式

tanγ=(khx2+khy2)2+(-khxαx-khyαy)2(1+αx2+αy2)2kz2, 　(14)cosθ=khy(1+αy2)+αxαykhx[khy(1+αy2)+αxαykhx]2+[khx(1+αx2)+αxαykhy]2。(15)

上述ODCIG向ADCIG映射的原理可以用图3所示的流程来实现。

图3 ODCIG向ADCIG映射流程

2 模型试算

首先用一个简单模型进一步说明宽方位角道集提取方法的实现原理和过程。将一个有限带宽的子波置于道集的零偏移距(hx=0, hy=0, hz=1.5 km)处, 模拟合成理想情况下聚焦的ODCIG道集I(hx, hy, z); 对道集作三维傅里叶正变换, 得到I(khx, khy, kz); 利用式(14)、式(15)所示的映射关系可以得到I(θ , γ , kz); 最后, 沿kz轴作一维傅里叶逆变换, 得到ADCIG点道集I(θ , γ , z)。图4所示的三维图像依次代表了四个域中的成像道集。

图4 映射过程中道集在不同域的表示

由于角度道集是偏移距道集经过映射得到的, 偏移距域的采样和范围决定着角度域的采样和范围, 因此, 理解两个域中采样和范围的关系是角度道集提取中重要的一部分。利用扩展成像条件得到的偏移距道集, 其采样和偏移距范围已经给定, 则需要分析哪个范围内的角度可以得到准确地恢复以及如何选取参数以消除空间假频的干扰。

从图4c中可以看到, 在方位角θ 等于45° 、135° 、225° 、315° 的情况下(红色箭头所示), 道集所能映射到的信息比其他方位角更多, 这是因为对于某一个kz, 当方位角为这四个角度时, 可以同时取到x、y方向的最大偏移距波数(khx, max, khy, max), 此时能映射到的反射角度为最大; 当方位角θ 从这四个角度分别向0° 、90° 、180° 、270° (绿色箭头所示)变化时, 映射的反射角角度范围逐渐减小。当反射角度|γ |小于25° 时, 不论哪个方位角, 角度道集的能量都能很完整地得到映射。

接下来对SEG/EAGE三维盐丘模型进行试算。采用C3-NA的观测系统, 所用的数据参数为:50条炮线(x坐标方向), 每条测线上激发97炮。每一炮有8条测线接收, 每条测线上最多有68个检波器。crossline方向(x坐标方向)的偏移距范围是(-700 m, 700 m), 和inline方向(y坐标方向)的偏移距范围是(-2 683 m, 140 m)。炮点的crossline方向和inline方向的间隔是160 m和80 m。接收点crossline和inline方向的间隔均为40 m。记录长度为4 912 ms, 8 ms采样。盐丘速度模型网格数大小为676× 676, x、y和z方向网格间隔均为20 m。

图5为盐丘模型的速度场。可以看出东北区盐丘发育, 并伴有多处断裂。图6为图5模型下z=0.48 km的切片, 其中红点区域为观测系统对应的炮点范围, 白色框为偏移成像数据体的范围(中心点范围)。

图5 三维盐丘模型的速度场

图6 炮点区域和成像数据体范围

采用裂步傅里叶法波动方程全三维DSR偏移算子对数据进行偏移。图7是从盐丘模型中截取的与成像结果对应的速度模型。图8是偏移成像结果, x、y方向分别为crossline和inline方向。从inline剖面和顶部水平切片可以看出, 盐丘边界和盐丘上部的断层都得到清晰的成像。从crossline剖面可以看出, 盐丘的冲顶周围成像效果欠佳, 其他部位成像较好。观察顶部水平切片剖面, 对比图7中的顶部水平切片剖面可以看出, 断层边界刻画的也非常精确。值得注意的是, 由于盐丘覆盖范围大, 厚度比较大, 会对一次反射波场产生屏蔽作用, 这造成盐丘下部构造成像非常困难, 但除此之外的其他部分均能清晰地反映构造成像。整个成像结果反映了偏移方法的正确性和有效性。

图7 盐丘成像体对应的速度模型

图8 盐丘偏移成像结果

在波场沿拓过程中采用扩展偏移距成像条件得到ODCIG道集, 如图8蓝色标线所示, 选择坐标为(6.16 km, 5.32 km)的CDP点M, 对其进行扩展成像, 提取ODCIG。图9和图10分别是M所在的x和y方向的垂直成像剖面局部倾角, 利用文中提出的道集映射方法就可以提取得到该点的ADCIG道集。图11是M点对应的ODCIG道集, 局部偏移距hx和hy的范围分别是(0.72 km, 0.72 km)和(0.72 km, 0.05 km), 可以看出inline方向的道集聚焦于零偏移距; 由于盐丘数据采的集观测系统不是真正的宽方位角采集, 而导致crossline方向的覆盖次数较少, 所以crossline方向的道集呈现一定程度的不聚焦。图12为M点对应的ADCIG, 前侧剖面表示方位角为0° 的角度道集, 右侧剖面表示反射角γ =0° 时不同方位角的成像深度。

图9 x=6.16 km的垂直成像剖面局部倾角

图10 y=5.32 km的垂直成像剖面局部倾角

图11 CDP点M对应的ODCIG

图12 CDP点M对应的ADCIG

图13 从图12中抽取θ =0° , 45° , …, 315° 的反射角度道集, 其中反射角γ ∈ (-45° , 45° )

为了分析不同方位角对应的角度道集, 从图12中沿着方位角θ 轴等间隔抽取θ =0° , 45° , …, 315° 的反射角度道集, 其中反射角γ ∈ (-45° , 45° ), 如图13所示。从图13可以看出, 0° 、180° 方位的角道集拉平程度较好, 且效角度范围较大(红色方框所示), 而90° 、270° 方位的角度道集会呈现一定的弯曲, 有效角度范围较小(红色方框所示)。这也说明crossline和inline方向角道集的拉平程度与偏移距道集的聚焦性相一致。

3 实际资料应用

图14是某地区的全三维DSR偏移结果。可以看出, 探区整体构造成像清晰, 右侧砂砾岩体的边界和左侧细小断层都得到精细地刻画。选择CDP点N(5.75 km, 3.9 km), 抽取其对应的ODCIG道集(图15所示), 可以看出道集在inline方向聚焦较好, 而在crossline方向由于偏移速度不准确而呈现不聚焦。图16是N点对应的ADCIG道集, 前侧剖面表示方位角为0° 的角度道集, 右侧剖面表示反射角γ =0° 时不同方位角的成像深度。从图16中沿着方位θ 角方向等间隔抽取θ =0° , 30° , …, 330° 的反射角度道集, 其中反射角γ ∈ (-45° , 45° ), 如图17所

示。从图17可以看出, 方位的角道集拉平程度较好, 有效角度范围较大, 而90° 、270° 方位的角度道集会呈现一定的弯曲, 有效角度范围较小(红色框所示)。这与盐丘模型分析结果一致, 说明该方法同样适用于实际资料。

图14 实际资料三维偏移成像体

图15 CDP点N对应的ODCIG

图16 CDP点N对应的ADCIG

图17 从图16中抽取θ =0° , 30° , …, 330° 的反射角度道集, 其中反射角γ ∈ (-45° , 45° )

4 结论

三维复杂构造情况下射线平面容易发生扭曲, 使得提取三维角度道集变得复杂。针对这一问题, 通过讨论反射点处射线参数矢量之间的几何关系, 实现了从偏移距共成像点道集(ODCIG)到角度域共成像点道集(ADCIG)的映射, 给出了提取角度道集的具体流程。理论推导和模型试算、实际资料应用可以得出以下结论。

(1)宽方位角度域共成像点道集表示的是偏移深度与局部方位角和反射角的关系, 速度不正确时会表现一定的剩余曲率, 且这种道集不受多路径对成像的影响, 可以为复杂条件下的速度反演提供道集依据。

(2)宽方位偏移速度分析的重要前提是成像算子也要体现宽方位, 即在偏移产生的偏移距道集中要充分体现不同方位的成像信息。笔者采用全三维DSR偏移, 偏移产生的偏移距道集包含了局部方位角的信息, 很适合向宽方位角道集映射。

(3)实现了真正的宽方位角度道集提取, 可以为全三维勘探地震成像中的深度偏移速度分析研究提供道集依据, 也可以用于三维的振幅与角度关系(AVA)分析。

The authors have declared that no competing interests exist.

参考文献
[1] Keho T H, Beydoun W B. Paraxial ray Kirchhoff migration[J]. Geophysics, 1988, 53(12): 1540-1546. [本文引用:1]
[2] Stolk C C, Symes W. Artifacts in Kirchhoff common image gathers[C]//Expand ed Abstracts of the 72nd SEG Annual International Meeting, 2002, 21: 1129-1132. [本文引用:1]
[3] Hill N R. Gaussian beam migration[J]. Geophysics, 1990, 55(11): 1416-1428. [本文引用:1]
[4] Gray S H. Gaussian beam migration of common-shot records[J]. Geophysics, 2005, 70(4): S71-S77. [本文引用:1]
[5] 李振春. 地震叠前成像理论与方法[M]. 东营: 中国石油大学出版社, 2011: 114. [本文引用:1]
[6] Baysal E, Kosloff D, Sherwood J W C. Reverse time migration[J]. Geophysics, 1983, 48: 1514-1524. [本文引用:1]
[7] Gazdag J. Wave equation migration with the phase-shift method[J]. Geophysics, 1978, 43: 1342-1351. [本文引用:1]
[8] Gazdag J, Sguazzero P. Migration of seismic data by phase shift plus interpolation[J]. Geophysics, 1984, 49: 124-131. [本文引用:1]
[9] Stoffa P L, Fokkema J T, Freir R M, et al. Split-step Fourier migration[J]. Geophysics, 1990, 55: 410-421. [本文引用:1]
[10] Popovici. Prestack migration by split-step DSR[J]. Geophysics, 1996, 61(5): 1412-1416. [本文引用:1]
[11] 程玖兵, 王华忠, 马在田, 等. 双平方根方程三维叠前深度偏移[J]. 地球物理学报, 2003, 46(5): 676-683. [本文引用:1]
[12] Biondi B, Symes W. Angle-domain common image gathers for migration velocity analysis by wavefield continuation imaging[J]. Geophysics, 2004, 69(5): 1283-1298. [本文引用:1]
[13] Biondi B, Tisserant T, Symes W. Wavefield-continuation angle-domain common-image gathers for migration velocity analysis[C]//Dallas: Expand ed Abstracts of the 73rd SEG Annual International Meeting, 2003: 2104-2107. [本文引用:1]
[14] Xu S, Zhang Y, Tang B. 3D common image gathers from reverse-time migration[C]//Expand ed Abstracts of the 80th SEG Annual International Meeting, 2010: 3257-3262. [本文引用:1]
[15] Sava P, Fomel S. Angle-domain common image gathers by wavefield continuation methods[J]. Geophysics, 2003, 68(3): 1065-1074. [本文引用:1]
[16] 王华忠, 程玖兵, 马在田, 等. 炮点—全偏移距域共方位角道集3D叠前深度偏移[C]//CPS/SEG 2004国际地球物理会议论文集, 2004: 325-328. [本文引用:1]
[17] 程玖兵, 王华忠, 耿建华, 等. 波动方程共方位角三维叠前偏移方法[J]. 石油地球物理勘探, 2006, 41(6): 629-634. [本文引用:1]
[18] Biondi B, Palacharla G. 3-D prestack migration of common-azimuth data[J]. Geophysics, 1996, 61: 1822-1832. [本文引用:1]
[19] Claerbout J F. Toward a unified theory of reflector mapping[J]. Geophysics, 1971, 36: 467-481. [本文引用:1]
[20] Rickett J, Sava P. Offset and angle-domain common image-point gathers for shot-profile migration[J]. Geophysics, 2002, 67: 883-889 [本文引用:1]
[21] 张凯. 叠前偏移速度分析方法研究[D]. 上海: 同济大学, 2008. [本文引用:1]
[22] Sava P, Fomel S. Coordinate-independent angle-gathers for wave equation migration[C]//SEG Technical Program Expand ed Abstracts, 2005: 2052-2055. [本文引用:1]