2. 中国科学院地质与地球物理研究所, 北京 100029;
3. 中国石油化工股份有限公司上海海洋油气分公司研究院, 上海 200120
2. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
3. Institute of SINOPEC Shanghai Offshore Oil & amp; Gas Company, Shanghai 200120, China
地震波正演技术在波形反演和逆时偏移中扮演重要的角色, 数据量大以及计算速度慢已成为制约这两项技术长期难以快速发展的瓶颈.尽管电子技术的进步加速了大规模并行计算机计算速度和计算尺度, 如现有并行集群系统已经能够解决二维的中、大规模或三维小规模逆时偏移和反问题[1], 但对于三维大尺度精细结构成像或岩性成像, 即便是在目前最先进的并行计算机上也依旧无能为力.
通用GPU计算有望改变这一局面, CUDA的出现更是减轻人们学习图形AP I的困难[2]. CUDA具有可利用类C等高级语言直接操作GPU硬件资源的特点, 使得GPU上数值算法实现和CPU上一样方便.作为CPU的一个协处理器, GPU具有更多的并行计算单元和更高的存储器带宽, 运行在主机上的并行数据和高密度计算通常被加载到这个设备上, 并分配到成百上千个不同线程上并行执行, 因此相对于串行CPU或者多核CPU来说具有更快的计算速度.
相对于传统的并行计算机, GPU计算机在性价比和功耗方面都具有明显优势, 使得GPU计算在HPC领域越来越流行.GPU计算涉及的研究领域包括生物物理、图像和信号处理、地球物理成像、分子动力学以及计算流体力学等方面[3~5].在地球物理方面, CGGVeritas的Deschizeaux等[4]在GPU上实现了地震资料处理的SRMIP算法, 该算法以频率域有限差分为核心, 在G80显卡上获得了比单线程CPU快8~15x的加速比.Micikevicius[6]通过研究三维空间优化型版的差分格式, 设计和实现了时间域的多阶有限差分算法, 以及多GPU之间多核异步执行和内存访问优化在不同区域及GPU/CPU之间通过消息传递的区域分解算法, 更高效地利用了GPU的计算能力和带宽, 并计算了规模较大的三维地震波正演问题.Jahanson[7]同样在GPU上实现了波动方程的有限差分正演方法, 与上面几位不同的是, 他分析和比较了该算法在不同GPU上的计算效率及带宽占用率, 为波形反演而非逆时偏移打下了基础.对于GPU在加速逆时偏移方面的应用, Abdelkhalek等[8]展示了GPU在地震波传播和逆时偏移方面的极大优势, 他通过充分利用线程块内的共享存储器带宽进一步加快了计算速率.Acceleware的Foltinek团队[9]以及GeoStar的刘钦、佟小龙团队[10]的研究成果更是代表了GPU逆时偏移的方向.李博等[11]在GPU上实现了非对称走时叠前时间偏移算法, 并将其用于实际地震资料的处理, 取得了较单个CPU快15x以上的加速比, 相当于10节点CPU集群计算速度的1.5~2.0倍.刘国峰等[12]通过GPU对叠前时间偏移算法进行了优化, 应用单颗Tesla C1060 GPU取得了较主频为2.5 GHz单核CPU快70x的加速比.Zhang等[13]利用cuFFT库在GPU上实现了三阶最优分裂步(OSP)深度偏移, 并比较了CPU上FFTW和GPU上cuFFT的执行效率, 通过对三维SEG/EDGE盐丘模型进行试算, 取得了较CPU深度偏移快25~40x的加速比.
在GPU上能够比较容易实现矩阵乘法运算, 而且还能直接套用现存的优化数学库cuBLAS[14].将一阶速度应力方程写成矩阵形式, 利用现存的GPU矩阵乘积算法和交错网格的Fourier伪谱微分矩阵算法[15], 应该可以加快波场模拟的计算速度和提高计算精度.对于Fourier伪谱方法, Fornberg[16]和Dault等[17]发现其精度可以看作高阶有限差分方法的阶数达到无穷时的极限情况.其交错网格形式不仅可以消除波场的谱在Nyquist波数处的不连续性, 还能消除"#"字现象和环形伪影[18], 这使得交错网格Fourier伪谱具有较常规方法更高的精度.
借助于伪谱微分矩阵算子精度高及GPU计算速度快等特点, 本文首先推导了交错网格的Fourier伪谱微分矩阵算子和一阶速度应力方程的矩阵形式, 然后利用GPU上优化的矩阵乘法, 对二维介质中的地震波场进行了快速计算, 为我们快速分析目标反射层在地震剖面中所对应的同相轴位置, 制定优化采集方案具有重要指导意义.
2 交错网格Fourier伪谱微分矩阵算子
交错网格的Fourier伪谱微分算子可以通过在常规网格算子中添加半个网格步长的相移算子
(1) |
其中, lj(x)为插值基函数, INu(x)在离散点xj处满足INu(xj)=u(xj).通过Fourier分析, lj(x)可以表示成
(2) |
又因为INu(x)的微分在Hilbert空间SN=span{eikx-N/2≤k≤N/2-1}上满足(INu)′=IN(u′).当N为偶数时, 存在如下的关系式:
(3) |
于是, u′(x)的插值基函数可以表示为
(4) |
如果在上式中令x=xi±1/2, 我们就可以得到插值基函数l′j(xi±1/2)的微分矩阵形式:
(5) |
该式即交错网格的Fourier微分矩阵算子.相对于普通格式的伪谱微分算子, 该算子的优势在于能使繁琐的微分计算转化成简单的矩阵乘积, 从而利于直接套用现存的矩阵乘积并行算法.交错网格的Fourier伪谱微分矩阵算子还能消除Nyquist波数处的不连续性, 从而使微分计算获得更高的精度[22, 23].
3 地震波场模拟的GPU实现 3.1 一阶速度应力方程及PML边界条件从弹性动力学本构方程出发, 可以推导出频率域各向同性介质中的一阶速度应力方程:
(6) |
其中, v为速度矢量, σ为二阶应力张量, c为弹性张量.将梯度算子沿∇平行和垂直于边界方向进行分解:
(7) |
其中, ∂n=
(8) |
通过引入扩展坐标系ñ(n)=n-
(9) |
利用关系∂n/∂ñ=iω/(iω+d(n)), 通过Fourier反变换可以将(9)式变回到时间域:
(10) |
考虑到纯声波介质内任一体积元表面只存在与法线方向相反的正应力而无剪切力, 弹性介质中的二维波动方程(10)便可以简化成与剪切力无关的声波方程:
(11) |
其中, 弹性模量K=ρV2, V为介质的纵波速度; σ在声波介质中为无向声压, σx和σz为声压沿两个坐标轴方向的分裂, vx和vz为两个坐标方向上的速度分量.在PML边界外物理模型区域内令d(x)=0和d(z)=0, 能确保(11)式中改进的带PML边界条件的波动方程退化成原无限介质中波动方程形式[26].
为了便于利用微分矩阵算子方法求解(11)式中关于应力和速度对空间的导数, 需要将波动方程(11)式离散并写成矩阵形式, 具体剖分方法可参考文献[21]关于黏弹性声波方程的离散方法.约定符号"⊗"表示矩阵对应元素之间的乘积, 于是声波方程(11)的矩阵形式可以表示成:
(12) |
其中, Vx和Vz分别表示沿两个坐标方向的速度场, Σx和Σz分别表示声压沿边界方向的分裂, 上标k为当前波场的模拟时刻; Bx和Bz分别表示浮度、Kx和Kz表示弹性参数; αx和αz分别为Δt/Δx和Δt/Δz; DxBH、DzBV、DxFH和DzFV分别为上节中推导的微分矩阵, 上标F表示前向微分, B表示后向微分, H表示沿着x方向, V表示沿z方向.Tx、Tz、T′x和T′z分别表示与PML边界阻尼函数d(n), n=x, z有关的矩阵.矩阵中元素分别表示为
(13) |
对于矩阵乘法的GPU实现, 很多文献做过详细的介绍和分析比较[14, 27].类似于基于消息传递通信MPI的矩阵乘法, GPU首先通过对矩阵进行分块, 将各块数据分配到各个线程块, 块内数据通过共享内存进行通信, 并利用块内各线程进行计算.GPU核心众多、存储器带宽高, 以及同时活动的线程数目多, 其计算能力远远高于CPU.另外, 高度并行的编程模式也是算法在GPU上获得高性能的基础.
对于矩阵相乘, 我们参考CUDA编程手册[27]中介绍的关于矩阵乘积的分块方法(如图 1所示), 将矩阵A和B根据设定线程块内线程数量BLOCK_SIZE×BLOCK_SIZE分解成Asub (A.width, BLOCK_SIZE)和矩阵Bsub (BLOCK_SIZE, B.height)大小子矩阵, 并保存在同一线程块内共享存储器中.受制于块内共享内存大小的限制, 通常我们还需要将子矩阵Asub和Bsub再次分解为BLOCK_SIZE×BLOCK_SIZE大小, 子矩阵的乘法在线程块内执行, 然后通过累加得到最终结果.由于从块内各线程在共享存储器中读取数据的带宽较从片外全局存储器中读取数据的带宽要宽很多, 而且从全局存储器中读取矩阵A和B的次数分别减小到B.width/BLOCK_SIZE和A.width/BLOCK_SIZE次, 因此通过将Asub和Bsub再次分解的子矩阵同时调入同一线程块共享内存上计算将有利于大大减小读取开销.从大小为1024×1024的两矩阵相乘看, 是否利用块内共享内存的两种方式在GTX295显卡上的加速比大约为4.5x.
在单机上我们分别用CPU (E5520)和GPU (GTX295)测试了不同大小矩阵的乘法, 计算时间如图 2所示.为显示清楚, GPU计算时间尺度扩大了20倍.GPU相对于CPU在计算矩阵乘积时的加速, 具体体现在矩阵规模为256×256其加速比为125, 而当矩阵规模增加到2048×2048其加速比增加到1691.这也就是说, 矩阵规模越大, 加速比越高, 十分适合计算区域范围较大的地震波场模拟.
波动方程的矩阵形式在GPU计算机上实现分Host (CPU)和Device (GPU)两部分:Host负责模型参数的读入、微分矩阵的生成、PML边界设置、各波场变量初始化及震源序列的生成, 然后通过cudaMemcpy或cudaMemcpy2D将数据从Host上传到Device上; Device负责GPU上的矩阵乘法和加法运算, 并将数据从Device下载回Host.整个波场模拟的流程如图 3所示.相比于CPU串行计算, GPU计算不同之处在于多了两个CPU到GPU以及GPU到CPU数据的拷贝过程.但最重要的是, GPU中矩阵乘法计算通过多线程来实现, 该并行化执行可以隐藏内存存取的等待时间.图 3中GPU计算部分对应于下面罗列计算流程的第2部分①~⑩.
GPU计算的具体流程描述如下:
(1) 在Host上读入模型参数ρ和V, 生成微分矩阵DxBH、DzBV、DxFH和DzFV, 设置PML边界Tx、Tz、T′x和T′z, 初始化各波场变量Vx、Vz、Σx和Σz, 生成震源时间序列s, 并将上述数据上传到Device上;
(2) 在Device上对时间t作循环, 并执行下面10个核运算:
①∂σ/∂x微分运算:
②vx速度场计算:
③∂σ/∂z微分运算:
④vz速度场计算:
⑤∂vx/∂x微分运算:DVxk+1/2<<<
⑥σx应力场计算:Σxk+1<<<Tx⊗Σxk+αxT′x⊗Kx⊗DVxk+1/2;
⑦∂vz/∂z微分运算:DVzk+1/2<<<DzFV
⑧σz应力场计算:Σzk+1<<<Tz⊗Σzk+αzT′z⊗Kz⊗DVzk+1/2;
⑨σ=σx+σz应力场求和:Σk+1<<<Σxk+1+Σzk+1;
⑩加入震源s项:Σk+1<<<Σk+1+s.
(3) 令t=t+Δt, 返回步骤(2), 直到t=tend, 然后从Device上下载数据回Host, 保存波场数据, 结束计算.
上述过程中, "<<<"表示加法或乘法运算可以写成核函数在D e v i c e上执行. D e v i c e上矩阵乘法采用共享内存策略减小数据计算时从块外全局内存的存取次数, 矩阵加法和对应元素相乘则采用全局存储器存取, 此类计算过程未作共享内存优化.
4 数值试验 4.1 Corner Edge模型对于交错网格Fourier伪谱微分矩阵算子的精度, 龙桂华等[21]在模拟弹性介质中波场传播的理论解和数值解时做过比较, 这里不再赘述.本节我们着重比较在GPU和CPU架构上模拟Corner Edge模型上地震波传播用时.测试硬件平台和3.2节一样, GPU采用2×240核GTX295显卡, 其着色器频率1.24 GHz, 核心频率为576 MHz; CPU采用4核Intel® Xeon® E5520, 主频为2.26 GHz.
设定Corner Edge模型密度恒定为1.5 g/cm3, 两层速度分别为1500 m/s和3500 m/s, 模型空间采样间隔为10 m.为较好地利用GPU线程块中的线程数以及内存访问优化, 模型大小被设置成线程块大小16×16的整数倍.模拟网格大小Lx×Lz分别为256×256、512×512、768×768、1024×1024和1280×1280, 四周包含20×20完美匹配层边界.震源位于地下(3/8Lx, 1/4Lz)处, 为主频25Hz的Ricker子波, 子波采样间隔为1 ms.
模型在0.5 s时的波场快照如图 4所示.地震波在界面处的反射、透射及角点绕射波形清楚, 侧面反射也能从波场快照中识别出来.顶边界PML吸收良好, 未有顶界面的虚假波场反射.从子波频率、介质速度以及网格空间采样长度之间的关系, 不难看出算子即便在每波长采样点数为3的情况下, 频散仍然抑制得非常好, 精度效果优于16阶交错网格的有限差分算子, 具体比较参考文献[21].
同一设置模型参数, 地震波场模拟在单GPU和单核单线程CPU上计算时间如表 1所示.从两者比较来看, 相对于CPU, GPU计算时间几乎可以忽略不计.即便在计算1280×1280大小模型时, 其消耗的时间也大约只有CPU上计算时间的0.43%, 约300 s.该时间包含Host生成微分矩阵以及预先设置参数以及边界条件的时间, 因此在核上的计算时间还可以进一步减小, 这就意味着波场计算本身在不考虑上述矩阵生成以及参数及边界条件预设时间时, 还能有比256×256时182x和1280×1280时231x更高的加速比.再则, 波场计算过程中矩阵之间的求和运算以及对应元素的乘积运算还可以通过类似矩阵乘积运算一样, 通过将全局存储器上的数据预取入共享存储器上减小读取开销等优化来进一步加快计算速度.
利用波场快照分析目标地质体在地震记录上对应的同相轴, 可以设置观测系统的最佳采集位置.为了利用波场快照指导观测系统设计, 我们选取非均匀介质模型如下图 5所示.模型大小为472×216, 空间采样间隔沿两个方向分别为25 m.介质为左边4层右边5层结构, 层界面起伏有致, 层中间有一右倾的高速盐丘.模型速度分布范围为2500 m/s到6500 m/s.由于模型比较简单, 各界面或目标层反射容易在地震记录中标示.
非均匀介质中地震波传播的波场快照如图 6所示.层界面和角点处反射和绕射清楚, 中间拱形高速区反射强、透射弱, 但遇到高阻层L4m后其反射信号明显再次增强, 这一现象在地表地震记录图 7中也可看到.由于高速层的屏蔽作用, 我们感兴趣的区域, 如L4r和L5r与L4m和底层高速岩层交汇处构造, 在炮点为(0 m, 0 m)和(2000 m, 0 m)时的地震记录Fig. 7a和Fig. 7b中并没有很好地显示出来.通过炮点平移, 根据上覆层岩性特点, 结合不同时刻地震波场快照图, 能够比较准确地锁定目标反射L5r在原始道集上的同相轴位置.根据上面计算的共炮记录剖面, 针对我们感兴趣区域构造反射, 应该能够指导我们设计最佳的采集窗口.
本文推导了交错网格的Fourier伪谱微分矩阵算子, 将矩阵乘积微分算法在GTX295GPU上进行了实现, 并比较了CPU和GPU构架的计算速度, 分析了制约GPU计算速度因素.通过利用线程块内共享存储技术减小线程从全局存储器上存取数据时间, 从而将GPU矩阵乘法计算速度提升了4.5倍.将速度应力方程中的空间微分转化成矩阵乘积形式, 并利用前面实施的GPU上微分矩阵乘积算法, 对不同规模的CornerEdge模型进行了实验, 分析比较了单CPU和单GPU的计算速率以及不同尺度下的加速比, 进一步描述了加速地震波模拟速度可采取的方案.在GPU上快速精确地进行地震波模拟, 可作为对正演要求高且计算密集的大尺度反演或逆时偏移的基础.通过模拟复杂非均匀介质中的地震波场, 结合不同时刻波场快照和地震剖面, 准确定位了目标反射在原始道集上的同相轴位置.根据剖面上特定地层反射的同相轴分析, 我们可以制定相应的采集方案, 这对复杂区域的地震数据采集具有指导意义.
致谢感谢两位匿名审稿人对本文所作的建议和有益讨论; 感谢中国科学院深圳先进技术研究院技术平台超算中心提供联想深腾7000G GPU集群服务器.
[1] | Vireux J, Operto S. An overview of full-waveform inversion in exploration geophysics. Geophysics , 2009, 74(6): 127-152. DOI:10.1190/1.3237087 |
[2] | http://developer.download.nvidia.com/compute/cuda/23/docs/CUDAGettingStarted2.3.Linux.pdf |
[3] | Owens J D, Houston M, Luebke D, et al. GPU Computing. Proceedings of the IEEE , 2008, 96(5): 879-899. DOI:10.1109/JPROC.2008.917757 |
[4] | Nguyen H, Zeller C, Hart E, et al. GPU Gems 3. Addison-Wesley Professional, 2007 |
[5] | Cohen J M, Molemaker M J. A fast double precision CFD code using CUDA. Proceedings of Parallel CFD, 2009 |
[6] | Micikevicius P. 3D finite difference computation on GPUs using CUDA. GPGPU-2:Proceeding of 2nd Workshop on General Purpose Processing on Graphics Processing Units (New York, USA), ACM, 2009. 79~84 |
[7] | Johansen O. Seismic Shot Processing on GPU. Norway:Norwegian University of Science and Technology, 2009 http://cn.bing.com/academic/profile?id=647157418&encoded=0&v=paper_preview&mkt=zh-cn |
[8] | Abdelkhalek R, Calandra H, Coulaud O, et al. Fast seismic modeling and reverse time migration on a GPU cluster. The 2009 High Performance Computation & Simulation-HPCS'09, 2009 |
[9] | http://www.acceleware.com/tasks/sites/default/assets/pdf/Acceleware_SEG_2009_RTMonGPU.pdf |
[10] | http://www.geostar.net.cn/show.php?id=85 |
[11] | 李博, 刘国峰, 刘洪. 地震叠前时间偏移的一种图形处理器提速实现方法. 地球物理学报 , 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. |
[12] | 刘国峰, 刘洪, 李博, 等. 山地地震资料叠前时间偏移方法及其GPU实现. 地球物理学报 , 2009, 52(12): 3101–3108. Liu G F, Liu H, Li B, et al. Method of prestack time migration of seismic data of mountainous regions and its GPU implementation. Chinese J. Geophys. (in Chinese) , 2009, 52(12): 3101-3108. |
[13] | Zhang J H, Wang S Q, Yao Z X. Accelerating 3D Fourier migration with graphics processing units. Geophysics , 2009, 74(6): 129-139. |
[14] | http://developer.download.nvidia.com/compute/cuda/2_3/toolkit/docs/CUBLAS_Library_2.3.pdf |
[15] | 龙桂华, 李小凡, 张美根. 错格傅里叶伪谱微分算子在波场模拟中的应用. 地球物理学报 , 2009, 52(1): 193–199. Long G H, Li X F, Zhang M G. The application of staggered-grid Fourier pseudospectral differentiation operator in wavefield modeling. Chinese J. Geophys. (in Chinese) , 2009, 52(1): 193-199. |
[16] | Fornberg B. The pseudospectral method:accurate representation of interfaces in elastic wave calculation. Geophysics , 1988, 53(5): 625-637. DOI:10.1190/1.1442497 |
[17] | Dault C R, Brail L W, Nowack R L, et al. A comparison of finite-difference and Fourier method calculations of synthetic seismograms. Bull. Seism. Am. Soc. , 1989, 79(4): 1210-1230. |
[18] | Özdenvar T, McMechan G A. Cause and reduction of numerical artifacts in pseudo-spectral wavefield extrapolation. Geophys. J. Int. , 1996, 126(3): 819-828. DOI:10.1111/gji.1996.126.issue-3 |
[19] | Bracewe R N. The Fourier transform and its application (2nd Ed.), McGraw-Hill, Inc., 1974 |
[20] | Fornberg B. High-order finite differences and the pseudospectral method on staggered grids. SIAM Journal on Numerical Analysis , 1990, 27(4): 904-918. DOI:10.1137/0727052 |
[21] | 龙桂华.粘弹性介质中的地震波传播与波形反演研究.北京:中国科学院地质与地球物理研究所, 2009 Long G H. The study of seismic wave propagation and waveform inversion in viscoelastic media (in Chinese). Beijing:Institute of Geology and Geophysics, Chinese Academy of Sciences, 2009 http://cn.bing.com/academic/profile?id=2670324508&encoded=0&v=paper_preview&mkt=zh-cn |
[22] | Corrêa G J, Spiegelman M, Carbotle S, et al. Centered and staggered Fourier derivative and Hilbert transforms. Geophysics , 2002, 67(5): 1558-1563. DOI:10.1190/1.1512801 |
[23] | Bale R A. Staggered grids for 3D pseudospectral modeling in anistropic elastic media. Calgary:CREWES Research Report, 2002 |
[24] | Chew W C, Liu Q. Perfectly matched layers for elastodynamics:a new modified Maxwell's equations with stretched coordinates. Microwave Opt. Technol. Lett. , 1996, 7: 599-604. |
[25] | Cillino F, Tsogka C. Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media. Geophysics , 2001, 66(1): 294-301. DOI:10.1190/1.1444908 |
[26] | Komatitsch D, Tromp J. A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation. Geophys. J. Int. , 2003, 154: 146-153. DOI:10.1046/j.1365-246X.2003.01950.x |
[27] | http://developer.download.nvidia.com/compute/cuda/2_3/toolkit/docs/NVIDIA_CUDA_Programming_Guide_2.3.pdf |