CPU/GPU协同并行计算的混合域全波形反演方法

摘要

本发明公开了一种混合域全波形反演方法，相对于传统的方法，本发明的方法通过CPU/GPU协同并行计算，显著提高了计算效率。本发明所公开的方法将全波形反演的正演部分放到时间域，即在时间域做正演，利用DFT转到频率域做反演，即采用离散傅里叶变换提取相应反演频率的波场成分，在频率域由低频到高频进行反演。本发明所公开的方法，有效地解决了常规时间域方法的收敛性问题，并且避免了常规三维频率域和拉普拉斯域波形反演内存占用量无法满足的问题，本发明的方法步骤少，减少了计算开销；并且通过GPU和CPU的协同并行计算，大幅度地提升了计算方法的加速比。

说明书

技术领域

本发明涉及一种全波形反演法，尤其涉及通过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 frequencydomain waveform inversion based on 3D time-domain forward modeling”，该发明所公开的技术方案中，提供了一种波形反演方法：首先将炮域数据利用DFT转换到频率域，然后利用时间域传播算子和DFT计算频率域震源正传波场，在频率域计算残差波场，并进一步利用逆离散傅里叶变换(Inverse DiscreteFourier Transform，简称IDFT)求取时间域的残差波场，之后，用时间域传播算子对残差波场进行反传，并同样利用DFT抽取频率域残差反传波场，最后，求取当前反演频率的梯度场。

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

发明内容

本发明公开了一种CPU/GPU协同并行计算的混合域全波形反演方法，可进行二维和三维的波形反演数据处理。本方法结合了时间域波形反演和频率域波形反演的优势，同时借助于GPU与CPU协同运算克服了时间域波形反演方法由于全频带数据的反演策略导致目标函数易陷入局部极小值的缺陷，并且避免了频率域和拉普拉斯域波形反演巨大的内存占用量，三维处理难以实施的问题。

通过本发明的技术方案，实现了下述技术效果：(1)针对时间域波形反演的数据处理加以改进，不仅保证反演的收敛性，而且相对频率域和拉普拉斯域的方法，极大地减小内存占用量；(2)采用GPU/CPU协同并行加速计算方法，最大程度使混合域波形反演的运算步骤得到并行处理，从而显著提高计算效率，大量节约计算时间。

在本发明中，使用了CPU、GPU等术语，并使用了函数、DFT等词，对于计算机领域技术人员而言，这些术语是显而易见的，CPU代表中央处理器，是一种具有通用计算能力的电子元件，常见的此类产品由Intel或AMD公司出品；GPU代表图形处理器，又称显卡，相对于CPU而言通用计算能力差，但是GPU专注于图像、图形等三维计算能力，在单位GPU内具有远高于CPU数目的并行的具有计算能力的计算单元，即使市场上的低端图形处理器所具有的计算单元(又称流处理器)也数以百计，由于GPU的功能是屏幕上合成显示数百万个像素的图像，也就是同时拥有几百万个任务需要并行处理，因此GPU在电路设计中即被设计成可并行处理很多任务。与CPU、GPU相关的是，CPU获得数据是通过内存实现数据传输的，GPU的数据传输则是通过显存；其中，DFT代表的是一种数学处理方法：Discrete Fourier Transform，也就是离散傅里叶变换，是连续傅里叶变换在时域和频域上都离散的形式，将时域信号的采样变换为在离散时间傅里叶变换(DTFT)频域的采样。

为达到本发明所述的技术效果，本发明是通过下述技术方案实现的：

CPU/GPU协同并行计算的混合域全波形反演方法，包括下述步骤：1，将波形反演的正演部分放到时间域；2，采用离散傅里叶变换提取频率域波场，并转到频率域进行反演。其中，上述步骤1、2是在GPU中完成计算的。

具体的，上述步骤由下述反演方法实现：

第一步，选取最优反演频率，并读取初始模型，反演参数及其他反演所需参数；对于本领域技术人员，这些参数是在反演开始前就确定的；

第二步，按照设定好的反演频率进行频率组循环或单频点循环，循环方式采用迭代循环；

第三步，进入炮点循环，读取单炮数据，确定单炮的成像孔径范围，并将相应的速度场由内存传输到GPU显存中；

第四步，用空间高阶时间域传播算子对震源进行正传，同时消除边界反射波；在正传的每个传播时刻内，记录地表接受波场，同时利用DFT抽取相应反演频率的频率域正传波场，即对每个时刻的波场快照计算该时刻波场对当前频率域波场的贡献，对这些贡献进行累加；

第五步，在时间域，对地表接受波场和实际炮记录做差，同时计算该反演频率对应的目标函数值；

第六步，以残差波场作为震源，利用时间域传播算子对其进行反传，利用DFT抽取相应反演频率的频率域波场，即累加所有时刻时间域波场对当前频率域波场的贡献；

第七步，计算当前炮记录在当前频率或频率组下对应的梯度场；

第八步，对下一炮集，重复第三步到第七步，同时累加该共炮点道集在当前频率下的梯度场，直到所有炮点循环完毕，然后求取此梯度场中的最大值用于优化步长的选取；

第九步，利用抛物线拟合法求取优化步长；

第十步，用优化步长和梯度场更新速度场，将此速度场继续存放在GPU显存中，用于下次迭代，重复上述步骤二到九，直到满足迭代终止条件或者达到最大迭代次数，然后进行下个频率或频率组的反演；

第十一步，当所有的频率或频率组均反演完毕，将最终更新的速度场由GPU显存传输到内存，并写入存储器中，完成反演过程；

其中，所述第四步至第十步的过程由GPU进行并行加速运算。

通过上述过程，本发明的方法显著提高了全波形反演的效率，具体表现在：

时间域传播算子采用空间高阶、时间二阶差分格式，波场传播和梯度场求取等过程中每个网格点都是相互独立的，所有网格点可以并行计算，即每个线程只计算一个或多个网格点的值，GPU并行计算核心多，能够高效率地进行细粒度并行运算，使得本发明的方法得以实现对混合域波形反演的加速处理。

另一方面，针对混合域算法的特点，本发明采用矩阵分片的思想。因为高阶有限差分算法，每个网格点波场值的计算都需要用到周围多个网格点的数据，内存读取冗余度很高，借助GPU共享内存(share memory)重用内存数据，可明显降低内存读取冗余度(read redundancy)，所以本发明的方法使用GPU共享内存有效提升了加速比；同时，与内存相比，share memory在GPU芯片内部，具有非常高的带宽，其速度接近于寄存器的速度，相对于内存具有很高的数据传输效率，通过共享内存对数据实现复用，显著改善了反演执行效率。

其中，所述第四步中消除边界反射波可以通过多种方式实现，这些方法为本领域广为采用，包括但不限于采用最佳匹配层条件、海绵吸收边界条件、随机边界条件、透射边界条件等。优选的，本发明采用最佳匹配层条件，即PML边界条件消除边界反射波。

在本发明中，所述第七步中求取梯度场的公式为

$>\mathrm{Grad}=\frac{1}{v^3\left(x\right)}\underset{\mathrm{shot}}{\Sigma }\underset{\mathrm{recev}}{\Sigma }\underset{o}{\overset{w}{\int }}(-{\omega }^2)P_f(x,\omega ;x_s)P_b{(x,\omega ;x_s,x_r)}^{*}\mathrm{dw}$>

其中，P_f为震源正传波场，P_b为残差反传波场，v是地下介质的速度，x_s和x_r分别是震源和检波器的水平位置，ω为当前反演的频率，*表示对变量取共轭。

在本发明中，第九步构建抛物线可以采用数学上已有的抛物线拟合算法实现抛物线拟合，优选采用下述方法构建抛物线，并求取极小值点：将抛物线坐标系的纵轴定义目标函数的值，横轴为步长，选取三个步长，和相应的目标函数值；已知α₀＝0时的目标函数值，之后通过迭代搜索α₁和α₂的值，其目标函数值F(α)满足F(α₁)＜F(α₀)且F(α₂)＞F(α₁)的条件；用(α₀，F(α₀))，(α₁，F(α₁))，(α₂，F(α₂))三组点确定一条抛物线，将其极小值点作为优化步长。

在本发明的混合域全波形反演法中，第四步还可以包括如下步骤：在波场正传的时间循环内，利用DFT算子求取每个时刻对频率域波场的贡献，直到所有的时刻都传播完毕，得到当前反演频率或频率组对应的频率域正传波场。

与上述类似，第六步还可以包括如下步骤：以残差波场作为震源反传，在该时间循环内，利用DFT算子求取每个时刻残差反传波场对频率域波场的贡献，直到所有的时刻都传播完毕，得到当前反演频率或频率组对应的频率域反传波场。

本发明的方法，CPU与GPU是协同运算的，并非所有运算均在GPU完成，具体的，二者的协同计算方式为CPU中调用GPU执行空间高阶有限差分GPU模块，PML吸收边界GPU模块，地表接受波场提取的GPU模块，波场正传时的离散傅立叶变化GPU模块，时间域残差求取的GPU模块，目标函数值求取的GPU模块，残差波场加载的GPU模块，残差波场反传时的离散傅里叶变换GPU模块，梯度场求取的GPU模块，求取梯度场最大值的GPU模块，优化步长求取的GPU模块，更新速度模型的GPU模块。

上述的模块，容易为本领域技术人员所理解，是执行所述功能的程序(可以采用任何编程语言实现，通常是C或C++)或多个程序封装好的程序包，通过利用GPU研发厂商，例如NVIDIA或AMD公司提供的已经封装好的函数库调用及修改来实现，CPU调用这些模块使得GPU执行模块中的指令以实现相应的功能。

通过上述技术改进，本发明的全波形反演方法可进行二维和三维的混合域波形反演，即在时间域做正演，之后利用DFT转到频率域做反演。本发明的方法技术方案基于GPU和CPU的协同并行计算有效地解决了常规时间域方法的收敛性问题，并且回避了常规三维频率域和拉普拉斯域波形反演内存占用量无法满足的问题；同时，相对前人的方法，本发明的方法计算步骤减少了一次DFT和一次IDFT。

附图说明

图1是本发明技术方案的流程图。

图2是技术方案第四步和第六步中两次DFT的流程图：(a)DFT抽取频域震源正传波场；(b)DFT抽取频域残差反传波场。

图3是技术方案第九步中抛物线拟合法求取优化步长的函数曲线示意图。

图4、5、6分别是实施例中所公开的初始速度模型、真实速度模型和反演速度模型。

其中底色为灰色的图框是采用GPU实现了加速的计算环节。

具体实施方式

下面结合附图对本发明的实施方式进行详细阐述，这些附图作为本发明的一部分，表示出了本发明的一种具体实现。在理解本发明的附图和优选实施例及本发明实质精神的基础上可知本发明允许多种不同形式的实施方式，因此在不脱离本发明的范围内，也可以有其它实施例被使用和构建，并且可能做出其它可选的改变。

如图1所示，本发明的混合域全波形反演方法，具有包括：

第一步，选取最优反演频率，并读取初始模型，反演参数及其他参数。

第二步，按照设定好的反演频率进行频率(组)循环；频率(组)循环内是迭代循环。

第三步，进入炮点循环，并读取单炮数据，确定成像孔径范围，并将速度场又CPU内存传输到GPU显存中。

第四步，如图2a所示，用时间域传播算子对震源进行正传，记录地表接受记录，同时利用DFT抽取相应反演频率的频率域正传波场。对每个时刻的波场快照计算该时刻波场对频率波场的贡献，并累加之。由于各个空间点是相互独立的，所以DFT的实现非常适合GPU细粒度并行计算。

第五步，在时间域，对地表接受记录和实际炮记录做差，同时计算该反演频率对应的目标函数值。

第六步，如图2b所示，以残差波场为震源，利用时间域传播算子对其进行反传，同样利用DFT抽取相应反演频率的频率域波场，累加所有时刻时间域波场对频率波场的贡献。

第七步，利用公式：

$>\mathrm{Grad}=\frac{1}{v^3\left(x\right)}\underset{\mathrm{shot}}{\Sigma }\underset{\mathrm{recev}}{\Sigma }\underset{o}{\overset{w}{\int }}(-{\omega }^2)P_f(x,\omega ;x_s)P_b{(x,\omega ;x_s,x_r)}^{*}\mathrm{dw}$>

计算当前炮记录在当前频率(组)下相应的梯度场，并求取梯度场的最大值。

第八步，对下一炮集重复第三步到第七步，同时累加该共炮点道集在当前频率下的梯度场，直到所有炮点循环完毕。

第九步，利用图3中的抛物线拟合求取优化步长。其中纵轴为目标函数的值，横轴为步长，选取三个步长和相应的目标函数值；已知α0＝0时的目标函数值，之后通过迭代搜索α1和α2的值，其目标函数值F(α)满足F(a1)＜F(α0)且F(α2)＞F(α1)的条件；进一步地，用这三组点来构造一个抛物线，并求取极小值点，以此作为优化步长。

第十步，用优化步长和梯度场来更新速度场，并将此速度场继续存放在GPU显存中，用于下次迭代，直到满足迭代终止条件或者达到最大迭代次数，然后进行下个频率(组)的反演，即重复第二步到第十步。

第十一步，当所有的频率均反演完毕，将最终更新的速度场由GPU显存传输到CPU，并写入存储器中。此时，反演完毕。

如图4-图6中所示，本实施例采用CAS(Chinese Academy of Sciences简写)字母模型来验证混合域方法的有效性。

所选用模型大小为2000m×500m，纵横采样间距均为5m，所选炮的数据为炮点个数为65炮，每炮400道，记录长度为1.6s，时间采样间隔为1ms，炮间距为30m。

反演参数设置为反演频率从2Hz到35Hz以1.5Hz为间隔选取23个频率，反演按照频率组进行循环，每个频率组选取两个频率。

迭代终止条件为：

(E(m_k+1)-E(m_k))/E(m_k+1)＜0.001

其中，E(mk)和E(mk+1)分别是第k次和第k+1次的目标函数值。

本实施例中用于反演的电脑采用了型号为i7-2600的CPU，型号为GTX680的GPU。通过附图可以看出，采用本发明的方法，反演结果很好地刻画真实速度模型。

具体的执行过程如下：首先将一炮数据从硬盘读入到计算机内存中，并确定该炮所对应的孔径大小，然后将该孔径内的速度和当前炮数据从CPU内存拷贝到GPU显存中去；之后在CPU中依次调用tstep_kernel(高阶有限差分GPU执行程序，用于在时间域进行波传播)，pml_kernel(PML吸收边界GPU执行程序，用于吸收来自边界的反射)，record_kernel(地表接受波场提取的GPU程序，用于从波场快照中提取相应的地表接受波场)，DFT_kernel(离散傅立叶变化GPU执行程序，用于求取当前时刻时间域波场对当前频率域波场的贡献)，residual_kernel(时间域残差求取的GPU程序，用于求取预测时域波场和真实波场的差值)，obj_kernel(目标函数值求取的GPU程序，用于计算当前速度模型下的目标函数值)，addresi_kernel(残差波场加载的GPU执行程序，用于在残差波场反传时，将残差波场以震源形式加载)，tstep_kernel和pml_kernel(用于传播残差波场)，DFT2_kernel(用于计算当前时域波场对当前频率域波场的贡献)，grad_kernel(梯度场求取的GPU程序，用于计算当前输入炮集对应的梯度场)，step_cuda(步长求取的GPU程序，此模块中包含多个kernel程序模块，用于以抛物线拟合的方法求取优化步长)，updatevel_kernel(更新速度模型的GPU程序，利用梯度场和优化步长来更新速度场)。所有GPU程序调用完毕以后，一个完整的迭代过程就结束了，之后进入到下一个迭代，直到满足迭代终止条件或者达到最大迭代次数为止。迭代结束以后，即进入下一个反演频率组，直至所有的反演频率执行完毕。

在该过程中，GPU承担了大部分的计算量，从而使计算效率大幅提升。

在上述具体实施中，使用了多个函数名称，如pml_kernel等，均可通过调用及修改NVIDIA公司(即GTX680研发商)公开提供的CUDA函数库中的函数实现。

而将上述CAS模型在单纯的i7-2600的CPU上进行的实验显示，CPU在与上述协同计算执行全波形反演相同的时间内无法完成此反演流程，得不到结果；实际上即使采用上述同样算法单纯使用CPU进行全波形反演，反演是无法在可接受时间内完成的，基于上述硬件平台，采用相同的算法，单纯的基于CPU的运算过程单次迭代需要耗时5034s，而本发明的方法单次迭代仅仅耗时241s，由此可见本发明的方法显著提升了计算效率，具有很好的加速比。

上述包括附图的实施例仅是本发明的示例，并且仅是提供了对于本发明的图解。其它变换形式的实施方式对于本技术领域中的技术人员而言是显而易见的。因此，本发明的保护范围由权利要求及说明书中公开的内容为准，而不局限于所给出的实施例。