2. 中国地质科学院矿产资源研究所, 北京 100037
2. Institute of Mineral Resource, Chinese Academy of Geological Sciences, Beijing 100037, China
逆时偏移(Baysal et al., 1983)是全波场成像方法,近年来,取得了很大的发展,已成功应用到某些复杂构造地区,尤其是那些盐丘和陡倾角地区.地球介质各向异性普遍存在,Thomsen(1986)提出弱各向异性理论.基于各向异性的假设,各向异性逆时偏移能够有效提高复杂区的成像精度,被视为高精度成像的突破点.近年来,偏移方法由各向同性逆时偏移发展到各向异性逆时偏移.而TTI介质作为各向异性介质的一种典型近似,成为应用和研究热点.
鉴于三维TTI介质弹性波方程求解困难,计算量较大的问题,相对简单的拟声波方程被提出以实现各向异性逆时偏移.Alkhalifah(1998, 2000)基于精确频散关系提出VTI介质的“声波假设”近似,即假设对称轴方向的剪切横波速度为零.然而,精确频散关系是一个四阶方程,为了避免四阶方程的求解,Zhou等(2006a)、Hestholm(2007)、Duveneck等(2008)基于声波假设提出引入不同辅波场形式将四阶方程转换为二阶方程组,以便于求解. Zhang和Zhang(2008)指出在某些复杂地区,VTI介质近似会导致波场的几何与运动学特征错误. Zhou等(2006b)、Zhang和Zhang(2008)、Fletcher等(2008, 2009)通过VTI介质方程的旋转变换提出一系列TTI拟声波方程.Fowler等(2010)对这些二阶方程组进行了总结分析,认为通过矩阵相似变换,这些方程之间可以进行相互转化,彼此是等价的.本文以Fletcher等(2009)提出的二阶方程组研究,在实际应用中发现,TTI介质逆时偏移有几个关键性技术.首先,伪横波干扰的压制.二阶方程组中qP波与qSV波是耦合在一起的,尽管设定了对称轴方向横波速度为零,但是仍然能够产生伪横波,这是方程近似导致假象(Grechk et al., 2004).针对这一问题,Chu等(2013)、Zhan等(2013)、Kim等(2013)以及Xu和Zhou(2014)等推导得到纯P波方程,Zhang和Zhang(2009)基于特征值分析提出了新的方程组,减少了伪横波干扰,这些纯波场的推导和计算求解复杂,计算量大.Duveneck等(2008)针对VTI介质,提出通过设置锥形过渡区域,实现各向异性参数匹配,该方法计算简单方便.其次,当假设对称轴方向速度为零时,应用到对称轴剧烈变化的TTI介质中会导致数值计算的不稳定问题,Fletcher等(2009)引入较小的横波速度可以很好改善这种波场传播不稳定性问题,李博等(2012)分析了TTI介质波场模拟的稳定性条件,认为选择合适的SV波速度可以使SV波的全区各向异性和反射系数达到极小,并可有效的抑制SV波对纵波勘探的影响.最后,计算效率问题,为了提高计算效率,高性能计算领域中GPU并行计算技术加速效果较好,基于Micikevicius(2009)提出的基于CUDA平台GPU上三维有限差分算法的计算策略,目前GPU并行加速计算技术在地震偏移成像中广泛应用(李博等, 2010, 2012;刘红伟等, 2010, 2011;Liu et al., 2013).
针对上述问题,本文基于多GPUs给出了一个有效实用的TTI逆时偏移的计算策略,其算法是基于Fletcher等(2009)年提出的TTI拟声波方程.通过分析各向异性参数匹配方法,确定了一种伪横波压制策略,压制效果较好.对比分析了随机边界值填充策略,有效实现了TTI介质随机边界以计算换存储的计算策略.本文针对GPU上TTI逆时偏移提出了多异步流的思路来提高计算并发度,节省了数据传输的时间,实现了单卡和双卡的21计算效率.最后,用模型验证了方法的有效性.
1 基本理论本文采用Fletcher等(2009)提出的拟声波方程进行TTI介质的逆时偏移计算,其理论公式为:
(1) |
其中:
p是主波场,q是辅助波场,ε和δ是Thomsen各向异性参数,θ和φ分别是对称轴的倾角和方位角,H1和H2是导数算子,Vpz是对称轴方向纵波速度,Vsz是对称轴方向横波速度,Vsz由参数σ和Vpz来控制,σ表达式为
(2) |
Fletcher等(2009)研究表明,采用声学假设,即Vsz=0,会导致波场在遇到对称轴变化剧烈时,传播不稳定,给定有限大小的横波分量(
有效实施3D TTI介质逆时偏移需要梳理和解决一些关键技术问题.本文旨在讨论以下三个关键核心技术问题,分别是伪横波压制,随机边界处理以及计算效率优化核心.
2.1 应用“各向同性层”压制伪横波伪横波干扰是制约TTI介质逆时偏移应用的一个关键问题.在Duveneck等(2008)研究的基础上采用了各向异性参数匹配法来压制伪横波干扰.需要指出的是,这里提到的伪横波指的是源产生的伪横波,其只产生于震源位置在非椭圆各向异性区域(ε≠δ)的情况.由Takanashi和Tsvankin(2011)提出的公式(3)可以看出,ε与δ相等(椭圆各向异性介质)则不会产生这种形式的伪横波.公式(3)为
(3) |
基于式(3),采用各向异性参数匹配方法压制伪横波.在震源位置周围设置一个区域,文中称之为“各向同性层”区域,在该区域内让各向异性参数匹配相等以压制伪横波.具体有多种匹配方式,如ε:=δ与δ:=ε(:=为赋值符号,注意二者差异,前者是ε值被改变;后者是δ值被改变),采用均匀TTI模型测试,模型参数分别为Vpz=2500 m/s,ε=0.24,δ=0.1,θ=45°,φ=0°.测试中选择了四种典型的匹配方式进行研究:(1)完全各向同性ε=δ=θ=φ=0;(2) ε:=δ, θ=φ=0;(3)ε:=δ;(4)δ:=ε.图 1展示了四种方式处理下的波场快照, 图 2b中紫色、黄色、绿色及红色依次对应上述四种匹配方式.由图可以得到,四种各向异性参数匹配方式对伪横波均有一定的压制效果,但是波前面形态却差异较大,波前面产生了畸变,也表明了各向异性参数匹配方式处理可能会影响波场的传播,为了验证其影响大小,与未做各向异性参数匹配处理的原始结果(图 2a)进行波前面对比,如图 2b所示,对比发现,红色曲线与黑线几乎重合,即方式(4)的匹配方式对波场传播几乎没有影响.
进一步分析造成方式(3)与方式(4)差别的原因,对波动方程进行分析,在“各向同性层”区域中,无论δ:= ε还是ε:= δ,都有Vsz=0,则式(1)可简化为:
(4) |
对于主波场P而言,主要由方程(4a)控制,即ε对主波场P的影响大,故在加载各向同性层时,需要尽量保持ε恒定,这也是为什么“δ:=ε”的加载方式对波前影响最小.进一步对比图 1,在各向同性层分界面附近产生了较强的反射干扰,分析其产生的原因,提出对各向异性参数δ各向同性层区域增加一个薄层过渡带.如图 3a所示,让δ由匹配值平滑缓慢的过渡到其真实值,各向同性层边界的红色虚线之间为过渡带,由图 3b可见,取得了很好的压制效果,并且保留了波前的正确位置.
为了分析上述各向异性参数匹配方法压制伪横波对成像结果的影响,图 4分别对比了压制前后的成像结果.采用的模型参数如下:模型网格大384×384×400,间隔dx=dy=dz=10,主频15 Hz,采模拟时间间隔0.5 ms,Vpz=2500 m/s,ε=0.24,δ=0.1,θ=45,φ=30.对比结果表明,引入过渡带的各向异性参数匹配方法,几乎不影响成像结果.
在TTI介质逆时偏移过程中,因为互相关成像条件要求,检波点波场与炮点波场需要在同一时刻提取成像值.炮点波场是从t=0递推到t=tmax,而检波点波场逆时反传,从t=tmax递推到t=0,这就要求存储其中一个递推过程的波场,在应用成像条件时,读取对应时刻波场.该方法需要大量的磁盘存储量以及频繁的磁盘读写访问,效率较低.对此,许多研究者选择源波场重建的方法以避免存储和读写源波场.
Clapp(2009)基于声波逆时偏移提出随机速度边界的思想,仅需要保存最后两时刻波场,利用保存两时刻炮点波场反传恢复正传波场,实现与检波点波场的同步反传外推,避免存储大量波场和频繁磁盘读写,几乎是波场的零存储.引入随机边界的计算方式能充分利用GPU并行计算的优势,以计算换存储,提高计算效率.Liu等(2013)等文章说明了其在各向同性RTM中的有效性.TTI方程可以作类似处理,速度模型在边界区域进行随机化.然而与常规各向同性声波方程比,TTI方程还包含其他四个模型参数,如何对速度外的这四个模型参数的外扩边界区域处理是一个值得考虑的问题.大体有如下3种方式:(1)边界外扩区域做类似速度参数处理,设置随机值;(2)外扩区域设置为各向同性区域,即边界层区域的模型参数填充为0;(3)外扩区域设置为均匀各向异性区域,即边界层区域的模型参数填充为边界常数值.方式(1)的随机处理需要非常谨慎,这四个各向异性模型参数的微小波动,会引起各向异性视速度的很大波动,使方程在模拟过程中发散.这里主要分析方式2与3对波场模拟的影响.以均匀TTI模型Vpz=2500 m/s、ε=0.24、δ=0.1、θ=45°、φ=0°进行分析.图 5展示了这两种填充方式的波场对比结果.图 5a与5b是方式(2)与(3)处理的ε模型,红色虚线框内是实际模型大小,红色虚线框外是外扩填充区域.图 5c和d是方式(2)、(3)处理的波场快照.从图中明显可以看出,方式(3)的处理效果最好,方式(2)在填充界面会出现干扰反射波.图 6展示了应用方式(3)的源波场随机边界的应用效果,图 6f与6l说明了随机边界条件可以保证炮点波场反传基本可以恢复其正传波场.
GPU加速计算技术在保证波场模拟的精度的同时,提高了计算效率,而且成本低、简单方便,在地震波场模拟和成像中广泛应用.然而,对于大规模的三维数据,单GPU的显存无法满足精细的波场模拟,多GPU计算则很好的解决了这一问题,既拓宽了计算显存,又不影响波场模拟的精度.
多GPU编程的实现基础是沿着慢维方向基于模型空间域的分解.多GPU计算一个很核心的问题是GPU间的通信及其效率.如何有效隐藏GPU间的通信延迟,优化线程调度,最大化GPU的占用带宽,以提高计算效率.
GPU间的实时通信是通过GPU的P2P技术实现的.如图 7所示,以双卡为例进行说明.整个成像区域沿着深度域分为两子模型域,分别对应GPU0与GPU1上的计算.其中每个子模型域沿深度分为三个区域,对应图中黄、绿、蓝三色区域.GPU的实时通讯基于绿色区域与蓝色区域的数据拷贝流实现,蓝色区域0-1(1-1)用来接收对应绿色区域1-0(0-0)的波场值,即图中箭头所示.
波场递推计算过程中,无论对于炮点波场正传,反传还是检波点波场的反传,其四个区域(GPU0与GPU1上黄色区域与绿色区域)的波场递推是互相独立的,利用GPU的异步流的技术可以实现四者的并发执行,在四者并发执行过程中,一旦绿色区域计算完毕,立即采用流技术控制实现P2P的拷贝,将GPU间的拷贝时间隐藏在黄色区域波场的波场递推中.由于采用了随机边界处理,在反传过程中,检波点波场与炮点波场实现了同步递推,每一时刻的炮点波场反传与检波点波场反传计算彼此独立,二者的递推也是可以并发执行的,极大的提高了并行度.
成像计算过程中,常规的互相关成像算法是:等待当前时刻波场递推结束,利用最新得到的下一时刻(it-1)的波场进行互相关提取成像值.改进后的成像算法:直接利用当前时刻(it)未被更新的波场进行互相关成像,而不是等待递推结束才利用最新的波场提取成像值.这样成像过程与炮点波场以及检波点递推三者就可以并发执行,极大的提高了并行计算的并行度.值得注意的是,提高并行度的目的是最大化GPU占用带宽,提高计算效率.如果带宽已经达到极限,再提高计算并行度也没有意义.因此在实际计算过程中,需要根据GPU带宽,适当的提高计算并行度,带宽允许,并行度越高越有利于计算效率的提高.图 8所示展示了反传过程中,波场递推及成像的并行计算过程.GPU间的数据拷贝延迟被隐藏,高并行度,提高了计算效率.
按照上述计算方法进行模型实验,测试算法的效率.使用的GPU是Tesla K10,其显存为4 GB,CPU为Xeon X5570,8核,主频2.93 GHz,16 G内存.对比分析了未优化CPU代码,单GPU代码以及多异步流技术应用前后的双GPU代码的性能,图 9展示了不同偏移孔径下50个时间步长下各代码的耗时结果,其中单GPU相比单CPU性能明显提高,与无异步流相比,基于异步流技术的高并发策略双GPU加速效果更好.图 10展示了应用异步流前后的加速比对比图,从图中可以直观的看出,不应用异步流技术,加速比在1.5上下波动,采用本文的异步流技术后,高并发策略下的双卡加速比接近理想的2.总之,基于多卡GPU计算的TTI介质逆时偏移成像方法实现了很好的加速效果.
采用正演的三维TTI盐丘模型进行实验,图 11仅展示了三维各向异性参数中的速度Vpz和对称轴倾角θ参数.该模型成像网格901×901×501,对应的网格间距15 m×15 m×10 m.炮间距120 m,检波点距30 m,记录时间8 s,采样时间间隔4 ms.图 12展示了两个主测线方向的TTI逆时偏移成像结果.成像结果表明,对目标盐丘成像归位准确,成像效果较好.
本文具体分析了实现3D TTI RTM的一些关键技术,体现在3个方面,多种参数匹配方式中,通过引入过渡带的δ:=ε的各向异性参数匹配方式,伪横波压制效果最好,并且不影响成像效果;针对常规逆时偏移算法存储瓶颈,对速度外其他参数的模型外扩区域进行均匀边界值填充,采用基于GPU的TTI RTM随机边界条件,减少了大量存储,避免了频繁的磁盘读写,有效提高了计算效率.提出的基于多GPU计算的多异步流控制的高并行策略,大大提高了GPU的占用率,很大程度上提高了计算效率,使双卡计算效率相比单卡的加速比接近理想的21,降低了计算成本,最后通过正演的复杂TTI模型验证了该算法的有效性.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持![] | Alkhalifah T. 1998. Acoustic approximations for processing in transversely isotropic media[J]. Geophysics, 63(2): 623–631. DOI:10.1190/1.1444361 |
[] | Alkhalifah T. 2000. An acoustic wave equation for anisotropic media[J]. Geophysics, 65(4): 1239–1250. DOI:10.1190/1.1444815 |
[] | Baysal E, Kosloff D D, Sherwood J W C. 1983. Reverse time migration[J]. Geophysics, 48(11): 1514–1524. DOI:10.1190/1.1441434 |
[] | Chu C L, Macy B K, Anno P D. 2013. Pure acoustic wave propagation in transversely isotropic media by the pseudospectral method[J]. Geophysical Prospecting, 61(3): 556–567. DOI:10.1111/gpr.2013.61.issue-3 |
[] | Clapp R G. 2009. Reverse time migration with random boundaries[C].//79th Annual International Meeting, SEG, Expanded Abstracts, 2809-2813. http://scitation.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=SEGEAB000028000001002809000001&idtype=cvips&gifs=Yes |
[] | Duveneck E, Milcik P, Bakker P M, et al. 2008. Acoustic VTI wave equations and their application for anisotropic reverse-time migration[C].//78th Annual International Meeting, SEG, Expanded Abstracts, 2186-2190. https://www.researchgate.net/publication/249862537_Acoustic_VTI_wave_equations_and_their_application_for_anisotropic_reverse-time_migration |
[] | Fletcher R, Du X, Fowler P J. 2008. A new pseudo-acoustic wave equation for TI media[C].//78th Annual International Meeting, SEG, Expanded Abstracts, 2082-2086. http://www.researchgate.net/publication/249862440_A_new_pseudo-acoustic_wave_equation_for_TI_media |
[] | Fletcher R P, Du X, Fowler P J. 2009. Reverse time migration in tilted transversely isotropic (TTI) media[J]. Geophysics, 74(6): WCA179–WCA187. DOI:10.1190/1.3269902 |
[] | Fowler P J, Du X, Fletcher R P. 2010. Coupled equations for reverse time migration in transversely isotropic media[J]. Geophysics, 75(1): S11–S22. DOI:10.1190/1.3294572 |
[] | Grechka V, Zhang L B, Rector J W. 2004. Shear waves in acoustic anisotropic media[J]. Geophysics, 69(2): 576–582. DOI:10.1190/1.1707077 |
[] | Hestholm S. 2007. Acoustic VTI modeling using high-order finite-differences[C].//77th Annual International Meeting, SEG, Expanded Abstracts, 139-143. https://academic.oup.com/gji |
[] | Kim Y, Cho Y, Jang U, et al. 2013. Acceleration of stable TTI P-wave reverse-time migration with GPUs[J]. Computers & Geosciences, 52: 204–217. |
[] | Li B, Li M, Liu H W, et al. 2012. Stability of reverse time migration in TTI media[J]. Chinese Journal of Geophysics, 55(4): 1366–1375. DOI:10.6038/j.issn.0001-5733.2012.04.032 |
[] | Li B, Liu H W, Liu G F, et al. 2010. Computational strategy of seismic pre-stack reverse time migration on CPU/GPU[J]. Chinese Journal of Geophysics, 53(12): 2938–2943. DOI:10.3969/j.issn.0001-5733.2010.12.017 |
[] | Liu G F, Liu Y N, Ren L, et al. 2013. 3D seismic reverse time migration on GPGPU[J]. Computers & Geosciences, 59: 17–23. |
[] | Liu H W, Li B, Liu H, et al. 2010. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation[J]. Chinese Journal of Geophysics, 53(7): 1725–1733. DOI:10.3969/j.issn.0001-5733.2010.07.024 |
[] | Liu H W, Liu H, Li B, et al. 2011. Pre-stack reverse time migration for rugged topography and GPU acceleration technology[J]. Chinese Journal of Geophysics, 54(7): 1883–1892. DOI:10.3969/j.issn.0001-5733.2011.07.022 |
[] | Micikevicius P. 2009. 3D finite difference computation on GPUs using CUDA[C].//Proceedings of the 2nd Workshop on General Purpose Processing on Graphics Processing Units, Expanded Abstracts, 79-84. http://dl.acm.org/citation.cfm?id=1513905 |
[] | Takanashi M, Tsvankin I. 2011. Correction for the influence of velocity lenses on nonhyperbolic moveout inversion for VTI media[J]. Geophysics: WA13–WA21. |
[] | Thomsen L. 1986. Weak elastic anisotropy[J]. Geophysics, 51(10): 1954–1966. DOI:10.1190/1.1442051 |
[] | Xu S, Zhou H B. 2014. Accurate simulations of pure quasi-P-waves in complex anisotropic media[J]. Geophysics, 79(6): T341–T348. DOI:10.1190/geo2014-0242.1 |
[] | Zhan G, Pestana R C, Stoffa P L. 2013. An efficient hybrid pseudospectral/finite-difference scheme for solving the TTI pure P-wave equation[J]. Journal of Geophysics and Engineering, 10(2): 025004. DOI:10.1088/1742-2132/10/2/025004 |
[] | Zhang H Z, Zhang Y. 2008. Reverse time migration in 3D heterogeneous TTI media[C].//78th Annual International Meeting, SEG, Expanded Abstracts, 2196-2200. http://www.researchgate.net/publication/240736933_Reverse_time_migration_in_3D_heterogeneous_TTI_media |
[] | Zhang Y, Zhang H Z, Zhang G. 2009. A stable TTI reverse time migration and its implementation[C].//79th Annual International Meeting, SEG, Expanded Abstract, 2794-2798. http://adsabs.harvard.edu/abs/2011Geop...76A...3Z |
[] | Zhou H, Zhang G, Bloor R. 2006a. An anisotropic acoustic wave equation for VTI Media[C].//68th EAGE Conference and Exhibition, EAGE, Extended Abstracts, H033. https://www.researchgate.net/publication/288944646_An_Anisotropic_Acoustic_Wave_Equation_for_VTI_Media |
[] | Zhou H B, Zhang G Q, Bloor R. 2006b. An anisotropic acoustic wave equation for modeling and migration in 2D TTI media[C].//76th Annual International Meeting, SEG, Expanded Abstract, 194-198. https://www.onepetro.org/conference-paper/SEG-2006-0194 |
[] | 李博, 刘红伟, 刘国峰, 等. 2010. 地震叠前逆时偏移算法的CPU/GPU实施对策[J]. 地球物理学报, 53(12): 2938–2943. DOI:10.3969/j.issn.0001-5733.2010.12.017 |
[] | 李博, 李敏, 刘红伟, 等. 2012. TTI介质有限差分逆时偏移的稳定性探讨[J]. 地球物理学报, 55(4): 1366–1375. DOI:10.6038/j.issn.0001-5733.2012.04.032 |
[] | 刘红伟, 李博, 刘洪, 等. 2010. 地震叠前逆时偏移高阶有限差分算法及GPU实现[J]. 地球物理学报, 53(7): 1725–1733. DOI:10.3969/j.issn.0001-5733.2010.07.024 |
[] | 刘红伟, 刘洪, 李博, 等. 2011. 起伏地表叠前逆时偏移理论及GPU加速技术[J]. 地球物理学报, 54(7): 1883–1892. DOI:10.3969/j.issn.0001-5733.2011.07.022 |