﻿ 叠前衰减补偿时间偏移及GPU实现
 石油地球物理勘探  2019, Vol. 54 Issue (1): 84-92  DOI: 10.13810/j.cnki.issn.1000-7210.2019.01.010 0
 文章快速检索 高级检索

### 引用本文

WU Jizhong, ZUO Hu. Attenuation compensation in prestack time migration and its GPU implementation. Oil Geophysical Prospecting, 2019, 54(1): 84-92. DOI: 10.13810/j.cnki.issn.1000-7210.2019.01.010.

### 文章历史

Attenuation compensation in prestack time migration and its GPU implementation
WU Jizhong , ZUO Hu
Jidong Oilfield Company, PetroChina, Tangshan, Hebei 063004, China
Abstract: The intrinsic absorption in subsurface media causes dissipation of seismic energy, thus decreases the amplitude and modifies the phase, and gives worse data resolution.We propose a prestack time migration in this paper, which compensates the attenuation along the seismic wave propagating path and enhance the resolution.For Q estimation, we calculate Q value with the joint use of VSP and seismic data.In general, we judge a right Q value based on visco-acoustic seismic synthetic seismogram and seismic traces at the well.For noise eli-mination, we use adaptive aperture according to dip angles, which can effectively remove migration noise.To improve computation efficiency, we adopt GPU calculating scheme.Model and real data results demonstrate that the proposed approach can not only compensate seismic energy but also improve data resolution.
Keywords: Q factor    prestack time migration    attenuation    adaptive aperture
0 引言

1 方法原理 1.1 井震联合求取Q

 $\frac{1}{{{Q_{{\rm{eff}}}}}} = \frac{1}{T}\sum\limits_{i = 1}^n {\frac{{\Delta {T_i}}}{{{Q_i}}}}$ (1)

 ${Q_{{\rm{well}}}} = \frac{{{\rm{ \mathsf{ π} }}\tau {\sigma ^2}}}{{{f_{{\rm{shot}}}} - {f_{{\rm{geo}}}}}}$ (2)

 $\delta = \sum\limits_\omega {\frac{{\rm{d}}}{{{\rm{d}}\omega }}\left[ {\ln F\left( \omega \right) - \ln {F_{\rm{Q}}}\left( \omega \right)} \right]}$ (3)

δ绝对值为零或最小时，表明FQ(ω)与F(ω)最接近，便认为对应的Q是最合理的。这种Q值求取方法与以往相比有一定改进与不同[11]。前人是对地下样本点进行“Q扫描+衰减补偿”处理，需要对叠前数据进行多轮次成像，计算量大，又由于增益控制因子的作用，低信噪比数据Q值求取稳定性差，另外在求取Q值过程中没有用到井信息，只利用了单一的叠前地震数据。而本文方法只对参考点进行“Q扫描+衰减”处理，对于叠后地震剖面而言，参考点选取若干即可，计算量大大降低。虽然该方法是在叠后开展的，但基于地震资料求取的Q值在空间上的趋势是有保证的，最后利用VSP数据计算的Q值对地面地震数据生成的Q场进行标定融合，从而又保证了Q值的数值精度。

1.2 衰减补偿叠前时间偏移

 $\begin{array}{l} P\left( {{k_x},\omega ,z + \Delta z} \right) = P\left( {{k_x},\omega ,z} \right)\exp \left( { - {\rm{i}}{k_z}\Delta z} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = f\left( \omega \right)\exp \left( { - {\rm{i}}{k_z}\Delta z} \right) \end{array}$ (4)

 ${k_z} = \sqrt {{{\left[ {\frac{\omega }{{c\left( \omega \right)}}} \right]}^2} - k_x^2}$ (5)

 $\begin{array}{l} P\left( {{k_x},\omega ,z + \Delta z} \right)\\ \;\;\;\;\;\;\; = f\left( \omega \right)\exp \left\{ { - {\rm{i}}\Delta z\sqrt {{{\left[ {\frac{\omega }{{c\left( \omega \right)}}} \right]}^2} - k_x^2} } \right\} \end{array}$ (6)

 $\left\{ \begin{array}{l} \frac{1}{{c\left( \omega \right)}} = \frac{1}{{v\left( \omega \right)}}\left( {1 - \frac{{\rm{i}}}{{2Q}}} \right)\\ v\left( \omega \right) = v\left( {{\omega _{\rm{c}}}} \right)\left( {1 + \frac{1}{{{\rm{ \mathsf{ π} }}Q}}\ln \frac{\omega }{{{\omega _{\rm{c}}}}}} \right) \end{array} \right.$ (7)

 $\frac{1}{{c\left( \omega \right)}} = \frac{1}{v}\left( {1 - \frac{{\rm{i}}}{{2Q}}} \right)\left( {1 - \frac{1}{{{\rm{ \mathsf{ π} }}Q}}\ln \frac{\omega }{{{\omega _{\rm{c}}}}}} \right)$ (8)

 $\begin{array}{l} P\left( {{k_x},\omega ,z + \Delta z} \right)\\ \approx f\left( \omega \right)\exp \left[ { - \frac{{{\rm{i}}\omega \Delta z}}{v}\sqrt {1 - \frac{2}{{{\rm{ \mathsf{ π} }}Q}}\ln \frac{\omega }{{{\omega _0}}} - {{\left( {\frac{v}{\omega }} \right)}^2}k_x^2 - \frac{{\rm{i}}}{Q}} } \right] \end{array}$ (9)

 $\begin{array}{l} P\left( {{k_x},\omega ,T} \right) \approx f\left( \omega \right) \times \\ \;\;\exp \left[ { - {\rm{i}}\omega \sum\limits_{i = 1}^n {\Delta {T_i}\sqrt {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_i}}}\ln \frac{\omega }{{{\omega _0}}} - {{\left( {\frac{{{v_i}}}{\omega }} \right)}^2}k_x^2 - \frac{{\rm{i}}}{{{Q_i}}}} } } \right] \end{array}$ (10)

 ${v_{{\rm{rms}}}} = \sqrt {\frac{1}{T}\sum\limits_{i = 1}^n {v_i^2\Delta {T_i}} }$ (11)

 $\begin{array}{l} \sum\limits_{i = 1}^n {\Delta {T_i}\sqrt {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_i}}}\ln \frac{\omega }{{{\omega _0}}} - {{\left( {\frac{{{v_i}}}{\omega }} \right)}^2}k_x^2 - \frac{{\rm{i}}}{{{Q_i}}}} } \\ \;\;\; \approx T\sqrt {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{{\rm{eff}}}}}}\ln \frac{\omega }{{{\omega _0}}} - {{\left( {\frac{{{v_{{\rm{rms}}}}}}{\omega }} \right)}^2}k_x^2 - \frac{{\rm{i}}}{{{Q_{{\rm{eff}}}}}}} - \\ \;\;\;\;\;\;\frac{{{\rm{i}}T}}{{2{Q_{{\rm{eff}}}}\sqrt {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{{\rm{eff}}}}}}\ln \frac{\omega }{{{\omega _0}}} - {{\left( {\frac{{{v_{{\rm{rms}}}}}}{\omega }} \right)}^2}k_x^2} }} \end{array}$ (12)

 $\begin{array}{l} P\left( {x,\omega ,T} \right) \approx \frac{\omega }{{2{\rm{ \mathsf{ π} }}}}\int {f\left( \omega \right)} \times \\ \;\;\;\exp \left[ {{\rm{i}}\omega \left( {T\sqrt {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{{\rm{eff}}}}}}\ln \frac{\omega }{{{\omega _0}}} - v_{{\rm{rms}}}^2p_x^2} + {p_x}x} \right)} \right] \times \\ \;\;\;\exp \left[ {\frac{{\omega t}}{{2{Q_{{\rm{eff}}}}}}{{\left( {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{{\rm{eff}}}}}}\ln \frac{\omega }{{{\omega _0}}} - v_{{\rm{rms}}}^2p_x^2} \right)}^{ - \frac{1}{2}}}} \right]{\rm{d}}{p_x} \end{array}$ (13)

 $\begin{array}{*{20}{c}} {P\left( {x,\omega ,T} \right) = f\left( \omega \right)\sqrt \omega \exp \left( { - \frac{{{\rm{i \mathsf{ π} }}}}{4}} \right)\frac{T}{{\tau _g^2v_{{\rm{rms}}}^2}} \times }\\ {\exp \left[ {{\rm{i}}\omega {\tau _g}\left( {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{{\rm{eff}}}}}}\ln \frac{\omega }{{{\omega _0}}}} \right)} \right]\exp \left( {\frac{{\omega {\tau _{\rm{g}}}}}{{2{Q_{{\rm{eff}}}}}}} \right)} \end{array}$ (14)

 $\begin{array}{l} {P^{\rm{D}}}\left( {x,\omega ,T} \right) = S\left( \omega \right)\sqrt \omega \exp \left( {\frac{{{\rm{i \mathsf{ π} }}}}{4}} \right)\frac{T}{{\tau _g^2v_{{\rm{rms}}}^2}} \times \\ \;\;\;\;\exp \left[ { - {\rm{i}}\omega {\tau _g}\left( {1 - \frac{2}{{{\rm{ \mathsf{ π} }}{Q_{{\rm{eff}}}}}}\ln \frac{\omega }{{{\omega _0}}}} \right)} \right]\exp \left( { - \frac{{\omega {\tau _{\rm{s}}}}}{{2{Q_{{\rm{eff}}}}}}} \right) \end{array}$ (15)

 $\begin{array}{l} I\left( {x,T} \right) = \frac{{{\tau _{\rm{s}}}}}{{{\tau _{\rm{g}}}}}\int {f\left( \omega \right)\sqrt \omega \exp \left( { - \frac{{{\rm{i \mathsf{ π} }}}}{4}} \right)} \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\exp \left[ {{\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{g}}}} \right)\left( {1 - \frac{1}{{{\rm{ \mathsf{ π} }}{Q_{{\rm{eff}}}}}}\ln \frac{\omega }{{{\omega _0}}}} \right)} \right] \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\exp \left[ {\frac{{\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{g}}}} \right)}}{{2{Q_{{\rm{eff}}}}}}} \right]{\rm{d}}\omega \end{array}$ (16)

1.3 倾角域自适应偏移孔径

 $\tan \left( {\beta + \alpha } \right) = \frac{{L + h}}{{{v_{{\rm{rms}}}}t}}$ (17)
 $\tan \left( {\beta - \alpha } \right) = \frac{{L - h}}{{{v_{{\rm{rms}}}}t}}$ (18)
 图 1 偏移孔径与地层倾角的关系 m为炮点s与接收点r的中点，设定输入地震数据的半炮检距为h，β为深度t处的地层倾角，P为深度t处最远成像道对应的成像点，p为成像点P在地面的投影，M为过点m的垂线与过P点的水平线的交点，α是地震波在成像点处的最大入射角

 $L = \frac{{{v_{{\rm{rms}}}}t\left( {{{\tan }^2}\beta - 1} \right) + \sqrt {v_{{\rm{rms}}}^2{t^2}\left( {1 + {{\tan }^2}\beta } \right) + 4{h^2}{{\tan }^2}\beta } }}{{2\tan \beta }}$ (19)

 图 2 地质模型(a)及其对应的倾角成像道集(b)
1.4 处理流程及GPU加速实现 1.4.1 计算流程

(1) 常规预处理；

(2) 在炮检距初叠剖面上选取若干CDP样本；

(3) 生成样本CDP对应的倾角成像道集，并确定样本CDP对应的倾角成像区；

(4) 利用样本CDP对应的倾角成像区插值出所有CDP对应的倾角成像区；

(5) 在CDP对应的倾角成像区内，对叠前数据进行衰减补偿成像处理；

(6) 重复步骤(5)，得到所有CDP最终的成像结果。

1.4.2 GPU加速实现

GPU设备适合并行非逻辑性运算，比如常规的加减乘除，而选择、判断等逻辑运算更适合用CPU设备去处理，因此利用GPU设备实现程序加速，不是简单的将程序加载到GPU设备上运算，而是将程序中计算量最大、耗时最长的部分，在这里称之为计算热点，放到GPU设备上进行加速，因此计算热点的选取直接决定了整体的加速效果。本文将频率点的衰减补偿运算作为计算热点，这是因为每个成像点都要在有效频带内逐频点做衰减补偿运算，计算量与耗时在整个程序中是最大的，并且各频率点的运算互不影响，没有先后顺序，更适合并行运算。计算粒度是用来描述计算单元大小的指标，一个线程可以计算一个频率点的衰减补偿，也可以计算整个成像空间所有频率点的衰减补偿，本文将前者称之为细粒度计算，后者称之为粗粒度计算。在细粒度计算中，线程并发量很大，每个线程的计算量却很小，线程的数据读取与交换用时比有效计算时间要大得多，总体加速效果不理想。在粗粒度计算中，线程的并发量太小，GPU中处于计算激活状态的单元少，设备中很多计算单元没有充分利用起来，负载率低，最终的加速效果也不好。借鉴这两种计算粒度的不足，本文对成像空间进行区块划分，每个线程只负责单个成像点的衰减补偿运算，既兼顾了线程计算量，又考虑了GPU的线程并发量，可以获得最佳的加速效果。

 图 3 GPU计算架构
2 数据应用 2.1 模拟数据

 图 4 二维模型

 图 5 常规叠前时间偏移剖面

 图 6 衰减补偿叠前时间偏移剖面

 图 7 x=690m处常规(a)和衰减补偿(b)叠前时间偏移单道时频谱对比
2.2 实际数据

 图 8 Q值求取过程 (a)常规合成记录；(b)衰减合成记录；(c)衰减合成记录与地震剖面的标定结果；(d)优化调整后的Q值曲线

 图 9 地震资料求取的Q场VSP标定结果

 图 10 常规叠前时间偏移剖面(a)及其反Q滤波结果(b)与衰减补偿叠前时间偏移剖面(c)对比

 图 11 常规叠前时间偏移剖面(a)及其反Q滤波结果(b)与衰减补偿叠前时间偏移剖面(c)局部细节对比

 图 12 常规叠前时间偏移剖面(蓝色)及其反Q滤波结果(绿色)与衰减补偿叠前时间偏移剖面(红色)的频谱对比
3 结论

 [1] Futterman W I. Dispersive body waves[J]. Journal of Geophysical Research, 1962, 67(13): 5279-5291. DOI:10.1029/JZ067i013p05279 [2] Kjartansson E. Constant Q wave propagation and attenuation[J]. Journal of Geophysics Research, 1979, 84(1): 4737-4748. [3] Hargreaves N D, Calvert A J. Inverse Q filtering by fourier transform[J]. Geophysics, 1991, 56(1): 519-527. [4] Wang Y. A stable and efficient approach of inverse Q filtering[J]. Geophysics, 2002, 67(1): 657-663. [5] 陈增保, 陈小宏, 李景叶, 等. 一种带限稳定的反Q滤波算法[J]. 石油地球物理勘探, 2014, 49(1): 68-74. CHEN Zengbao, CHEN Xiaohong, LI Jingye, et al. A band-limited and robust inverse Q filtering algorithm[J]. Oil Geophysical Prospecting, 2014, 49(1): 68-74. [6] Mittet R, Sollie R, Hokstak K. Prestack depth migration with compensation for absorption and dispersion[J]. Geophysics, 1995, 60(1): 1485-1494. [7] Zhang J, Wapenaar K. Wavefield extrapolation and prestack depth migration in anelastic inhomogeneous media[J]. Geophysics Prospecting, 2002, 50(1): 629-643. [8] Tonn R. The determination of the seismic quality factor Q from VSP data:A comparison of different computational methods[J]. Geophysical Prospecting, 1991, 39(1): 1-27. DOI:10.1111/gpr.1991.39.issue-1 [9] Quan Y, Harris J M. Seismic attenuation tomography using the frequency shift method[J]. Geophysics, 1997, 62(1): 895-905. [10] Xin K, Hung B.3-D tomographic Q inversion for compensation frequency dependent attenuation and dispersion[C]. SEG Technical Program Expanded Abstracts, 2009, 28: 4014-4017. http://library.seg.org/doi/abs/10.1190/1.3255707 [11] Zhang C, Ulrych T J. Estimation of quality factors from CMP records[J]. Geophysics, 2002, 67(1): 1542-1547. [12] Zhang J F, Wu J Z, and Li X Y. Compensation for absorption and dispersion in prestack migration:An effective Q approach[J]. Geophysics, 2013, 78(1): 1-14. [13] 吴吉忠, 刘成权, 张建坤, 等. 倾角域自适应孔径叠前时间偏移[J]. 石油地球物理勘探, 2017, 52(3): 502-508. WU Jizhong, LIU Chengquan, ZHANG Jiankun, et al. Adaptive aperture prestack time migration in the dip angle domain[J]. Oil Geophysical Prospecting, 2017, 52(3): 502-508. [14] Klokov A, Fomel S.Optimal migration aperture for conflicting dips[C]. SEG Technical Program Expan-ded Abstracts, 2012, 31: 504-507. [15] 陈志德, 初海红. 黏滞介质吸收补偿叠前时间偏移的倾角道集孔径优选[J]. 石油地球物理勘探, 2016, 51(3): 444-451. CHEN Zhide, CHU Haihong. Optimization of vari-able migration aperture in dip angle gathers for pre-stack time migration with absorption compensation in viscous elastic media[J]. Oil Geophysical Prospecting, 2016, 51(3): 444-451. [16] Audebert F, Froidevaux P, Racotoarisoa H, et al.Insights into migration in the angle domain[C]. SEG Technical Program Expanded Abstracts, 2002, 21: 1188-1191. https://www.onepetro.org/conference-paper/SEG-2002-1188 [17] 刘和年, 王建立, 吴蕾, 等. 地层倾角约束自适应孔径叠前时间偏移[J]. 石油地球物理勘探, 2014, 49(5): 899-903. LIU Henian, WANG Jianli, WU Lei, et al. Adaptive aperture prestack time migration with dip constrain[J]. Oil Geophysical Prospecting, 2014, 49(5): 899-903. [18] 李博, 刘国峰, 刘洪. 地震叠前时间偏移的一种图形处理器提速实现方法[J]. 地球物理学报, 2009, 52(1): 245-252. LI Bo, LIU Guofeng, LIU Hong. A method of using GPU to accelerate seismic pre-stack time migration[J]. Chinese Journal of Geophysics, 2009, 52(1): 245-252. [19] 李肯立, 彭俊杰, 周仕勇. 基于CUDA的Kirchhoff叠前时间偏移算法设计与实现[J]. 计算机应用研究, 2009, 26(12): 4474-4477. LI Kenli, PENG Junjie, ZHOU Shiyong. Implement Kirchhoff prestack time migration algorithm on CUDA architecture[J]. Application Research of Compu-ters, 2009, 26(12): 4474-4477. DOI:10.3969/j.issn.1001-3695.2009.12.020 [20] 马召贵, 赵改善, 武港山, 等. Kirchhoff叠前时间偏移的GPU移植与性能优化技术[J]. 石油学报, 2014, 35(4): 700-705. MA Zhaogui, ZHAO Gaishan, WU Gangshan, et al. GPU-based porting and optimization of Kirchhoff pre-stack time migration[J]. Acta Petrolei Sinica, 2014, 35(4): 700-705. [21] 陈志德, 赵忠华, 王成. 黏滞声学介质地震波吸收补偿叠前时间偏移方法[J]. 石油地球物理勘探, 2016, 51(2): 325-333. CHEN Zhide, ZHAO Zhonghua, WANG Cheng. Prestack time migration with compensation for absorption of viscous-elastic media[J]. Oil Geophysical Prospecting, 2016, 51(2): 325-333. [22] 吴吉忠, 杨晓利, 龙洋. 一种稳定高效的等效Q值反Q滤波算法及应用[J]. 石油地球物理勘探, 2016, 51(1): 63-70. WU Jizhong, YANG Xiaoli, LONG Yang. A roubust approach of inverse Q filtering with equivalent Q[J]. Oil Geophysical Prospecting, 2016, 51(1): 63-70. [23] Bleistein N, Cohen J K, Stockwell J W. Mathematics of Multidimensional Seismic Inversion[M]. Springer, 2001. [24] 卢宝坤, 张剑锋, 蒲泊伶. 角度域叠前时间偏移:非均匀介质中的单程波方法[J]. 地球物理学报, 2012, 55(11): 3795-3804. LU Baokun, ZHANG Jianfeng, PU Boling. Pre-stack time migration in angle domain:a one-way propagator approach in inhomogeneous media[J]. Chinese Journal of Geophysics, 2012, 55(11): 3795-3804. DOI:10.6038/j.issn.0001-5733.2012.11.026 [25] 吴吉忠, 李君, 石文武. 南堡1号构造火成岩发育区稳健保幅处理技术[J]. 石油地球物理勘探, 2017, 52(增刊1): 1-9. WU Jizhong, LI Jun, SHI Wenwu. Robust amplitude-preserving processing strategy for well-developed igneous rocks in the structure Nanpu 1[J]. Oil Geophysical Prospecting, 2017, 52(S1): 1-9.