地球物理学报  2013, Vol. 56 Issue (1): 262-268   PDF    
引用本文
王保利, 高静怀, 陈文超 等. 逆时偏移中用Poynting矢量高效地提取角道集. 地球物理学报, 2013, 56(1): 262-268, doi: 10.6038/cjg20130127.  
Wang B L, Gao J H, Chen W C, et al. Extracting efficiently angle gathers using Poynting vector during reverse time migration. Chinese J. Geophys. (in Chinese), 2013, 56(1): 262-268, doi: 10.6038/cjg20130127.  
逆时偏移中用Poynting矢量高效地提取角道集
王保利1,2 , 高静怀1 , 陈文超1 , 张唤兰3     
1. 西安交通大学电子与信息工程学院, 西安 710049;
2. 中煤科工集团西安研究院, 西安 710077;
3. 西安科技大学地质与环境学院, 西安 710054
收稿日期 2011-08-31, 2012-11-20 收修定稿
基金项目: 国家自然科学基金(40730424);国家科技重大专项(2011ZX05023-005, 2011ZX05040-002)和国家科技支撑计划(2012BAK04B04-01)资助
作者简介: 王保利, 男, 1981年生, 博士, 主要从事地震偏移成像和高性能计算方法研究.E-mail:pooly1981@163.com
通讯作者: 高静怀, 男, 1960年生, 教授, 博士生导师, 主要从事地震波传播、成像、地震资料处理以及信号处理方面的研究.E-mail:jhgao@mail.xjtu.edu.cn
摘要: 逆时偏移在提高复杂介质的成像质量方面表现出了优越的性能, 但逆时偏移对速度精度的要求比较高.共成像点道集是一种非常重要的叠前深度偏移输出, 它除了能为深度偏移处理提供速度信息外, 还能够提供振幅和相位等信息, 为后续的属性解释提供依据.本文提出一种在逆时偏移成像过程中提取角度域共成像点道集的方法, 该方法采用矢量波动方程进行波场传播, 并用能流密度矢量(Poynting vector)计算反射角, 最后应用互相关成像条件输出角度域共成像点道集.该方法简单易于实现, 且几乎不需要额外的计算量和存储量, 非常适合于进行逆时偏移速度分析, 同时提出的角道集也能用于进行AVA等分析.最后通过模型算例和实际数据检验了方法的有效性和优越性.
关键词: 逆时偏移      角道集      速度分析      Poynting矢量      AVA     
Extracting efficiently angle gathers using Poynting vector during reverse time migration
WANG Bao-Li1,2, GAO Jing-Huai1, CHEN Wen-Chao1, ZHANG Huan-Lan3     
1. The School of Electronic and Information Engineering, Xi'an Jiaotong University, Xi'an 710049, China;
2. Xi'an Research Institute of China Technology & Engineering Group Corp, Xi'an 710077, China;
3. College of Geology & Environment, Xi'an University of Science and Technology, Xi'an 710054, China
Abstract: Reverse time migration has shown its superior performance for imaging the complex subsurface structures, but it needs a more accurate velocity model. Common-image gathers are an important output of prestack depth migration. They not only can provide velocity model for depth migration, but also can provide amplitude and phase information for subsequent subsurface attribute interpretation. In this paper, we have developed a method for extracting angle domain common-image gathers during reverse time migration. In our method, vector wave equation is used to calculate the wavefield, energy flux density vector (Poynting vector) is used to calculate the reflection angle, and finally cross-correlation imaging condition is used to obtain the angle domain common-image gathers. Finally, a model and a field marine data are done to show the effectiveness and advantages of this new method..
Key words: RTM      ADCIG      Velocity analysis      Poynting vector      AVA     
1 引言

20世纪80年代就提出的逆时偏移[1-2]在对复杂地下结构进行成像时体现出了较强的优势, 因此越来越受到人们的重视[3-10].多路径、多次波、大倾角、回转波、棱镜波以及保幅等关键问题上, 逆时偏移几乎都可以较好地解决.但准确的逆时偏移成像比其它成像方法(如单程波偏移、Kirchhoff偏移和beam偏移等)对速度模型的准确度更加敏感.为了给逆时偏移建立一个较准确的速度模型, 需要一种迭代的方法来不断地更新速度.而共成像点道集(CIGs)是叠前深度偏移中的很重要的输出道集, 它可以提供速度建模所需要的信息, 以及反映地下岩性的振幅和相位等信息.正确的速度模型得到的共成像点道集中的同相轴是平的, 而不正确的速度得到的是弯曲的, 且弯曲的程度与速度校正量是有关系的, 通过这种关系就可以不断地更新速度模型[11].

通常, 有两种共成像点道集, 一种是偏移距域的(common-offset CIGs, 简称为COCIGs), 一种是角度域的(angle-domain CIGs, 简称为ADCIGs).在地下反射体是水平层状的情况, 这两种道集是可互换的.Kirchhoff偏移和beam偏移方法可以很方便快速地输出COCIGs[12].但多路径问题会影响到COCIGs的质量, 进而会影响速度分析和AVO分析的精度.而单程波偏移和逆时偏移中的COCIGs是很难提取的.Save和Formel(2003)[13]利用单程波偏移在二维情况下首先提取了地下局部偏移距域的CIGs, 然后再转换为ADCIGs.Fomel(2004)[14]将这种方法推广到三维情况, 并指出三维情况下该方法非常复杂.Save和Fomel(2006)[15]又提出一种新的方法, 该方法需要先生成延时CIGs, 然后再转换为ADCIGs.Xu等(2010, 2011)[16-17]提出一种直接从三维逆时偏移中提取ADCIGs的方法, 该方法在频率波数域求取传播方向, 进而计算出反射角和方位角信息, 最后使用给定的成像条件直接输出ADCIGs.这种方法可以得到质量较好的ADCIGs, 但需要保存所有时刻的炮点波场值和接收点波场值, 并进行四维的傅里叶变换, 计算量和存储量都非常得大, 不适合进行逆时偏移速度分析等处理.

Yoon和Marfurt(2003)[5]用坡印廷矢量来确定传播方向并用此修改成像条件来去除逆时偏移中的低频噪声, 取得了较好的效果.本文借助于这种思想, 用一阶应力速度方程(与常规标量应力波动方程相比, 增加了额外的存储量, 但减少了计算坡印廷矢量的耗时)进行波场延拓, 并用坡印廷矢量求取成像点处的传播方向, 最后使用给定的成像条件输出逆时偏移中的ADCIGs.

2 方法原理 2.1 波场传播

逆时偏移是通过波动方程在时间上对采集到的地震数据进行反向传播来实现的.其偏移过程使用的是双程波波动方程, 因此不受倾角限制, 且能正确地对回转波、多次波等进行成像.逆时偏移的原理本文不再赘述.

通常, 逆时偏移采用的声波波动方程为标量应力波动方程[18-22, 10], 其形式如下:

(1)

其中, υ为介质纵波速度, P为质点应力波场.该波动方程不利于使用坡印廷矢量(原因见下节).因此本文采用一阶应力-速度声波波动方程:

(2)

式中:ρ为介质密度, υ为介质纵波速度, v和w分别为x和z轴上的质点振动速度.本文对公式(2)采用交错网格技术[18]来构建波场传播算子.考虑到逆时偏移的计算量和为避免数值频散问题对逆时偏移结果的影响, 在时间方向采用二阶而在空间方向采用高阶差分格式(≥8阶).

2.2 Poynting矢量求取反射角

坡印廷矢量也称为能流密度, 1884年由坡印廷提出用来测量单位时间内穿过垂直于此矢量方向的单位表面的电磁场能量.Poynting矢量值等于电场强度与磁场强度的叉积, 矢量方向代表电磁波的传播方向.

Yoon和Marfurt(2003)[5]给出了地震波场中的坡印廷矢量计算式:

(3)

其中, v表示位移速度矢量, v=(v, w), P为应力波场.通过(3)式可知, 如果使用(1)式的标量声波波动方程, 需计算应力波场值的时间和空间导数, 计算比较复杂.而如果采用(2)式的一阶应力速度方程, 则应力P和位移速度矢量v都是逆时偏移过程中已有的, 不需任何额外计算.图 1是均匀介质中的某时刻计算的坡印廷矢量快照图, 可以很明显地判断出波的大致传播方向:在水平分量图(a)中蓝色和红色分别代表向右传和向左传, 而在垂直分量图(b)中蓝色和红色则分别代表向下传和向上传.

图 1 坡印廷矢量快照图 (a)水平分量; (b)垂直分量. Fig. 1 The snapshots of Poynting vector (a) Horizontal component; (b) Vertical component.

用SS和SR分别表示炮点波场的坡印廷矢量和接收点波场的坡印廷矢量, 则根据余弦定理, 可以求出两矢量的夹角θ,

(4)

由此可以得到波场的反射角为.为了增加抗噪性能, 夹角θ可在一小时窗内进行计算, 如

(5)

2.3 提取ADCIGs

获得反射角β后, 便可以进行ADCIGs的提取.常规的逆时偏移成像条件是对所有角度道集的叠加, 因此不能直接用来提取ADCIGs, 需要进行修改.考虑到计算得到的反射角通常与我们给定的采样角度(一般是等间隔的)不一致, 我们用高斯采样函数将常规的成像条件

(6)

修改为

(7)

其中S和R分别表示炮点波场和接收点波场(本文中为应力波场值), 而β表示当前时刻该成像点处计算得到的反射角, βk表示ADCIGs中的离散角, σ表示高斯函数的方差.根据(7)便能容易地提取出ADCGIs.

需特别说明的是, (5)式中, 分母在很小甚至为零情况下, 计算得到的角度β误差较大, 但此时, 其能量贡献也非常小, 因此不会影响最终角道集的质量.

具体的计算流程如图 2.

图 2 Poynting矢量提取角道集的流程图 Fig. 2 The flow chart of extracting ADCIGs using Poynting vector
3 模型试验 3.1 Marmousi模型的ADCIGs

用本文提出的方法提取出的Marmousi模型的ADCIGs如图 3所示.角度从0°~90°, 角度间隔为2°, 每100个CMP显示一个ADCIG, 一共23个ADCIGs.从图上可以看出每个ADCIG上的同相轴都是拉平的, 且低频噪声主要分布在大角度的道上. 图 4为第500个CMP位置处的ADCIG.

图 3 Marmousi模型的ADCIGs Fig. 3 The ADCIGs of Marmousi model
图 4 Marmousi模型第500个CMP位置处的ADCIG Fig. 4 A ADCIG located at 500th CMP of Marmousi model
3.2 低频噪声

逆时偏移的低频噪声主要是由透射波成像造成的, 而入射波和透射波的夹角在模型比较光滑情况下接近于180°(传播路径一微小段范围内), 因此可通过在角道集上对大角度的道进行衰减或直接置零来对逆时偏移成像剖面进行去噪.理论上这种去噪方法等同于拉普拉斯滤波[7, 11].图 5、图 6和图 7分别为用0°~40°、40°~70°和70°~90°等角道集叠加后的结果, 可明显看出在0°~40°内几乎没有任何低频噪声, 而在70°~90°的偏移剖面几乎全部是低频噪声.

图 5 Marmousi模型逆时偏移剖面(0°~40°) Fig. 5 The RTM image of Marmousi model(0°~40°)
图 6 Marmousi模型逆时偏移剖面(40°~70°) Fig. 6 The RTM image of Marmousi model(40°~70°)
图 7 Marmousi模型逆时偏移剖面(70°~90°) Fig. 7 The RTM image of Marmousi model(70°~90°)
3.3 AVA

角道集的另一个用途就是用于AVA分析[23-25], 本文使用单层模型分析了角道集上振幅的特征.图 8为用振幅随角度递增的模型参数(上下层速度分别为1500m/s和1800 m/s, 密度均为2.1g/cm3)得到的角道集AVA曲线与理论AVA曲线的对比.图 9为用振幅随角度递减的模型参数(上下层速度分别为1500m/s和1300 m/s, 密度均为2.1g/cm3)得到的角道集AVA曲线与理论AVA曲线的对比.对比可以看出, 两种模型得到的角道集的AVA曲线和理论AVA曲线基本一致.

图 8 单层模型得到的角道集AVA曲线与理论AVA曲线的对比(上下层速度分别为1500 m/s和180 0m/s, 密度均为2.1 g/cm3) Fig. 8 The comparison of the oretical AVA curve and calculated AVA curve from single-layer model (The velocity of the upper and lower layers are 1500m/s and 1800m/s respectively, the density was 2.1 g/cm3)
图 9 单层模型得到的角道集AVA曲线与理论AVA曲线的对比(上下层速度分别为1500 m/s和1300 m/s, 密度均为2.1 g/cm3) Fig. 9 The comparison of the oretical AVA curve and calculated AVA curve from single-layer model (The velocity of the upper and lower layers are 1500m/s and 1300m/s respectively, the density was 2.1g/cm3)
4 实际数据试验

图 10为从某地区的一条测线中提取出的部分ADCIGs.可以看出1.3~1.8km段中同相轴基本是平的, 而在1.8~2.3km断层同相轴是向下弯曲的, 说明偏移过程中所用的速度偏高.通过不断迭代地进行速度调整即可实现逆时偏移速度分析.

图 10 实际数据中提取的部分角道集 Fig. 10 Part of the ADCIGs from field data
5 结论

本文在逆时偏移过程中采用一阶应力速度波动方程进行波场传播，并用Poynting矢量来确定传播方向并计算反射角，进而通过给定的成像条件进行逆时偏移角道集的提取.整个计算过程中几乎不额外增加计算量和存储量, 提取的共角度道集可用于进行速度分析和AVA分析等.该方法可用于三维逆时偏移中角道集等的提取(包括方位角道集), 实现简单, 易于进行并行加速(如GPU等).

后续工作将在此工作基础上, 使用提取的角道集进行逆时偏移速度分析, 以得到较好的偏移速度剖面.

参考文献
[1] Whitmore N D. Iterative depth migration by backward time propagation. 53rd Annual International Meeting, SEG Expanded Abstracts, 1983, 2: 827-830.
[2] Baysal E, Kosloff D D, Sherwood J W C. Reverse time migration. Geophysics , 1983, 48(11): 1514-1524. DOI:10.1190/1.1441434
[3] 张美根, 王妙月. 各向异性弹性波有限元叠前逆时偏移. 地球物理学报 , 2011, 44(5): 711–719. Zhang M G, Wang M Y. Prestack finite element reverse-time migration for anisotropic elastic waves. Chinese J. Geophys (in Chinese) , 2011, 44(5): 711-719.
[4] Yoon K, Shin C, Hong S. 3D reverse-time migration using acoustic wave equation. An experience with the SEG/EAGE data set. The Leading Edge , 2003, 22(1): 38-41. DOI:10.1190/1.1542754
[5] Yoon K, Marfurt Kurt J, Starr W. Challenges in reverse-time migration. 74th Annual International Meeting, SEG Expanded Abstracts, 2003: 1057-1060.
[6] Bednar J B, Bednar C J, Shin C. Two-way versus one-way: A case study style comparison. 76th Annual International Meeting, SEG Expanded Abstracts, 2006, 25: 2343-2347.
[7] Zhang Yu, Xu S, Bleistein N, et al. True-amplitude, angle-domain, common-image gathers from one-way wave-equation migrations. Geophysics , 2007, 72(1): S49-S58. DOI:10.1190/1.2399371
[8] Zhang Y, James S. Practical issues of reverse time migration: True amplitude gathers, noise removal and harmonic-source encoding. First Break , 2009, 27(1): 53-60.
[9] Liu F Q, Zhang G Q, Morton S A, et al. Reverse time migration using one way wavefield imaging condition. 77th Annual International Meeting, SEG Expanded Abstracts, 2007, 26(1): 2170-2174
[10] Liu F Q, Zhang G Q, Morton S A, et al. An effective imaging condition for reverse-time migration using wavefield decomposition. Geophysics , 2011, 76(1): S29-S39. DOI:10.1190/1.3533914
[11] Zhou H B, Gray S H, Young J, et al. Tomographic residual curvature analysis: The process and its components. 73rd Annual International Meeting, SEG Expanded Abstracts, 2003, 22: 666-669.
[12] Xu S, Chauris H, Lambaré G, et al. Common-angle migration: A strategy for imaging complex media. Geophysics , 2001, 66(6): 1877-1894. DOI:10.1190/1.1487131
[13] Sava P C, Formel S. Angle-domain common-image gathers by wavefield continuation methods. Geophysics , 2003, 68(3): 1065-1074. DOI:10.1190/1.1581078
[14] Fomel S. Theory of 3-D angle gathers in wave-equation imaging. 74th Annual International Meeting, SEG Expanded Abstracts, 2004, 1(1): 1053-1056.
[15] Sava P, Formel S. Time-shift imaging condition in seismic migration. Geophysics , 2006, 71(6): S209-S217. DOI:10.1190/1.2338824
[16] Xu S, Zhang Y, Tang B. 3D common image gathers from Reverse time migration. 80th Annual International Meeting, SEG Expanded Abstracts, 2010: 3257-3260.
[17] Xu S, Zhang Y, Tang B. 3D angle gathers from Reverse time migration. Geophysics , 2011, 76(2): S77-S92. DOI:10.1190/1.3536527
[18] 李博, 刘国峰, 刘洪. 地震叠前时间偏移的一种图形处理器提速实现方法. 地球物理学报 , 2009, 52(1): 245–252. Li B, Liu G F, Liu H. A method of using GPU to accelerate seismic pre-stack time migration. Chinese J. Geophys (in Chinese) , 2009, 52(1): 245-252.
[19] 李博, 刘红伟, 刘国峰, 等. 地震叠前逆时偏移算法的CPU/GPU实施对策. 地球物理学报 , 2010, 53(12): 2938–2943. Li B, Liu H W, Liu G F, et al. Computational strategy of seismic pre-stack reverse time migration on CPU/GPU. Chinese J. Geophys (in Chinese) , 2010, 53(12): 2938-2943.
[20] 刘红伟, 李博, 刘洪, 等. 地震叠前逆时偏移高阶有限差分算法及GPU实现. 地球物理学报 , 2010, 53(7): 1725–1733. Liu H W, Li B, Liu H, et al. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation. Chinese J. Geophys (in Chinese) , 2010, 53(7): 1725-1733.
[21] 刘红伟, 刘洪, 邹振, 等. 地震叠前逆时偏移中的去噪与存储. 地球物理学报 , 2010, 53(7): 2171–2180. Liu H W, Liu H, Zou Z, et al. The problems of denoise and storage in seismic reverse time migration. Chinese J. Geophys (in Chinese) , 2010, 53(7): 2171-2180.
[22] 刘红伟, 刘洪, 李博, 等. 起伏地表叠前逆时偏移理论及GPU加速技术. 地球物理学报 , 2011, 54(7): 1883–1892. Liu H W, Liu H, Li B, et al. Pre-stack reverse time migration for rugged topography and GPU acceleration technology. Chinese J. Geophys (in Chinese) , 2011, 54(7): 1883-1892.
[23] 杨午阳. 模拟退火叠前AVA同步反演方法. 地球物理学进展 , 2010, 25(1): 219–224. Yang W Y. Pre-stack simultaneous AVA inversion algorithm with simulated annealing. Progress in Geophys (in Chinese) , 2010, 25(1): 219-224.
[24] 王保利. 弹性阻抗反演及应用研究. 地球物理学进展 , 2005, 20(1): 89–92. Wang B L. Elastic impedance inversion and its application. Progress in Geophys (in Chinese) , 2005, 20(1): 89-92.
[25] 喻岳钰, 杨长春, 王彦飞, 等. 叠前弹性阻抗反演及其在含气储层预测中的应用. 地球物理学进展 , 2009, 24(2): 574–580. Yu Y Y, Yang C C, Wang Y F, et al. Application of pre-stack seismic elastic impedance inversion to gas reservoir. Progress in Geophys (in Chinese) , 2009, 24(2): 574-580.