[发明专利]CPU/GPU协同并行计算的混合域全波形反演方法

说明书

技术领域

本发明涉及一种全波形反演法，尤其涉及通过CPU/GPU协同并行计算在混合域上执行的全波形反演方法，属于计算机领域。

背景技术

目前我国西部及海外老油区开发的难度越来越大，深度不断加深，对高精度的地球物理成像及反演方法需求很大；而高精度的成像方法最终还要依赖高精度速度模型的构建。速度模型的反演方法根据反演信息使用的不同分为走时反演方法，振幅反演方法和全波形反演三类方法。由于全波形反演方法综合利用叠前地震波场的动力学和运动学信息，能够高精度地重建地下速度场，成为目前国内外勘探地球物理研究的重要方向。

已有的全波形反演(Full Waveform Inversion，简称FWI)方法包括Tarantola在1984年提出的基于广义最小二乘反演理论的时间域全波形反演，随后Pratt，Shin等人在其基础上分别发展了频率域和拉普拉斯域的波形反演。全波形反演直接利用地震记录反演地下介质参数，充分地利用反射、折射等波形信息，在一定条件下，可以同时反演纵波速度、密度和横波速度的相对变化，从而得到精度较高的速度模型，是目前国内外勘探地球物理研究的重要方向。虽然全波形反演具有诸多的技术优势，但是由于全波形反演计算量巨大(约为逆时偏移的30-100倍)，难以实现三维处理，导致这种方法无法适应工业生产的需要。

在已有的研究中，对于全波形反演方法从计算方式上分为三类：第一类是时间域波形反演(Tarantola 1984，1986)，该方法利用全频带数据计算梯度场，由于目标函数的极度非线性导致其易陷入局部极值，从而使反演不收敛；第二类是频率域波形反演(Pratt 1998，1999)，此方法在频率域做正演，之后在频率域从低频到高频进行逐频(组)反演，相对高频数据，低频数据与速度模型具有更好的线性关系，该方法可以有效避免时间域反演过程中的局部极小问题；在频率域正演的方法分为直接法(如：LU分解等)和迭代法两种，由于迭代法效率低，频率域正演多采用LU分解等直接法来求解，这种方法的优势是对阻抗矩阵进行一次分解后，分解后的结果可以用于所有炮点的计算，计算效率高，该种方法也存在着难以克服的弊端，在三维情形下，LU分解内存占用量和计算量均十分巨大，导致三维波形反演实现起来不现实。第三类是拉普拉斯域波形反演(Shin2008，2009)，该方法在频率域波形反演的基础上，对反演参数增加了衰减因子，所以反演不仅是从低频到高频，还从浅层到深层，从而减小了反演结果对初始模型的依赖，但是此方法也面临与频率域方法一样的内存瓶颈，也就是在其正演部分需要LU分解，致使三维计算实现困难。

目前在国内，随着勘探复杂度的增加，许多专家学者也开展了波形反演的研究及应用工作，但是一般限于二维，难以用于三维。总体而言，波形反演技术从其诞生之日起，尤其是近几年来一直是学术界研究的热点问题；国际上波形反演的研究开展的较多，尤其是国外工业界各个环节的研究都在紧张进行。相反，国内波形反演的研究应用基本上处于起步阶段，相对比较落后，因此，研究实用性强的波形反演方法对填补我国的技术空白具有重要的意义。

2007年的美国专利US11/756，384，“System and method for 3D frequency domain waveform inversion based on 3D time-domain forward modeling”，该发明所公开的技术方案中，提供了一种波形反演方法：首先将炮域数据利用DFT转换到频率域，然后利用时间域传播算子和DFT计算频率域震源正传波场，在频率域计算残差波场，并进一步利用逆离散傅里叶变换(Inverse DiscreteFourier Transform，简称IDFT)求取时间域的残差波场，之后，用时间域传播算子对残差波场进行反传，并同样利用DFT抽取频率域残差反传波场，最后，求取当前反演频率的梯度场。

与本领域早期采用的计算方案相比较，该技术方案虽然取得了一定的进步，但其流程步骤相对复杂，执行效率受到一定影响；而且，混合域FWI中主要的耗时模块(如DFT，时间域传播算子等)并行粒度很细，仅采用CPU运算，未能最大程度地使程序并行化，计算效率会受到很大影响。本发明的反演策略与其不同，在此反演流程的基础上少了一次DFT和一次IDFT；并进一步，采用GPU/CPU协同并行运算，最大程度地使程序并行，明显地增加了相对于CPU程序的加速比。

发明内容