一种高效的时间域全波形反演方法
【技术领域】
[0001] 本发明属于石油地震勘探数据高效率、低存储、高精度速度建模领域，具体涉及一 种高效的时间域全波形反演方法，针对时间域全波场速度建模所需巨大计算量和海量存储 问题，采用震源编码、多GPU结合的并行计算技术策略和无数据I/O的有效边界存储策略完 成高效、高分辨率速度建模。
【背景技术】
[0002] 利用叠前地震波数据运动学和动力学信息，全波形反演具有揭示复杂地质条件下 的地质构造、岩性参数的能力，目前其被认为是速度建模精度最高的技术。然而三维全波形 反演的实现还受到种种限制，其中两个主要难点是巨大的计算量和存储量一直未得到较好 地解决。
[0003] 由于全波形反演的正确实施需成百上千次正演模拟，因而正演模拟的效率是制约 全波形反演发展的本质所在。对二维反演问题，采用直接解法的频率域波动方程正演模拟 技术在处理多震源模拟时效率较高，然而对三维问题，即使采用效率较高的嵌套剖分直接 解法，频率域正演的计算复杂度也将达到n snw0 (n4l〇gn)，所需内存0 (η5)，其计算量和存储 量巨大，目前的计算机条件很难实现。对三维正演，时间域有限差分正演的时间计算复杂度 为〇(η 4)，所需内存复杂度为0(η3)，相对于三维频率域正演，时间域正演具有计算量小，内 存占用低等优点。目前对正演模拟的加速研究主要有两种思路：一、通过计算机硬件的加速 以提高反演的计算效率（如基于CUDA平台的GPU并行计算等），二、提高算法本身的计算效 率，基于随机相位编码的多震源反演技术思路较好，其基本思想是把所有的单个震源数据 叠加成单个或多个超级震源数据，并对超级震源数据进行处理。由于全波形反演的计算效 率与震源的个数近似成正比，因而采用超级震源可以大大提高全波形反演的计算效率。
[0004] 常规时间域全波形反演面临的另外一个难题是需要海量数据输入和存储，海量数 据的输入输出（I/O)不仅对计算机硬件提出了较高的挑战，而且严重影响了反演的计算效 率。即使采用GPU并行计算技术对需输出数据的常规时间域全波形反演加速时，效率提升 有限。Moussa (2009)采用GPUPU对逆时偏移研究时表明91 %的运行时间用于CPU内存和 GHJ显存的数据传递，这是由于各时刻的震源波场需存储在硬盘上，限制了 GPU的并行效率 所致。
【发明内容】
[0005] 本发明针对时间域全波形反演所需计算量和存储量巨大的问题，提供一种高效的 时间域全波形反演方法，为大规模数据的高精度速度建模提供一种高效的处理手段。
[0006] 本发明是通过以下技术方案实现的：
[0007] 第一步、采用源极性随机、源位置随机、源个数随机三随机方案对原始观测地震记 录随机形成不同频带超炮集；
[0008] 第二步、采用基于CUDA的多GPU并行计算技术加速时间域一阶形式的声波方程系 统和弹性波动方程系统的高阶交错网格有限差分正演模拟，得到GPU模拟的第ifreq频带 超炮集；在正演模拟时采用CPU和GPU协同技术，用多GPU核模拟不同的震源传播，把震源 波场的边界数据和最后一个时刻所有波场存储于计算机内存；在计算正向传播的震源波场 时把各个时刻内部区域与PML区域交界的质点振动速度分量的波场值及最大时刻所有波 场值存储在GPU显存中；
[0009] 第三步、计算波场残差及目标泛函值（用模拟的超炮集和原始超炮集相减得波场 残差叠加求和得目标泛函值），然后判断是否满足终止条件，如果是，表示该频段的全波形 反演迭代结束，并把反演的模型作为下一个频带全波形反演的初始速度模型，转入第六步， 如果否，则转入第四步；
[0010] 第四步、采用多GPU并行计算技术计算检波器反传波场，同时采用有效存储边界 策略和多GPU并行计算技术重构震源波场并计算目标泛函相对于模型参数的梯度；
[0011] 第五步、采用梯度类迭代算法更新速度模型，对更新后的速度模型计算超级炮记 录并与第一步的超级炮记录做残差；
[0012] 第六步、判断是否满足频率终止条件，如果是，则转入第七步，如果否，则返回第一 I K 少；
[0013] 第七步、输出当前模型。
[0014] 所述第一步是这样实现的：
[0015] 采用MPI并行输入原始单炮记录和原始雷克子波，利用维纳滤波器对原始单炮记 录和原始雷克子波滤波，得到不同频带的地震单炮记录和地震子波；
[0016] 利用C语言库函数中产生随机种子的函数、随机函数以及符号函数对所有地震子 波赋予[-1，1]之间的一个数，即赋予每个地震子波不同的颜色，根据地震子波颜色的正负 性，赋予地震子波相应的极性，此过程为源极性随机；
[0017] 采用堆排序算法或快速排序法对改变源极性的地震子波的颜色排序，颜色值逐渐 增大，并采用C语言中的随机函数根据总超级震源子波个数及总单个地震子波个数随机出 每个超级震源中单个地震子波的个数，此过程为源个数随机；
[0018] 按照每个超级震源的单个地震子波的个数对排序后改变源极性的单个地震子波 叠加，此时叠加的每个地震子波的位置也相应随机，此过程为源位置随机，由此完成三随机 方案形成不同频带的超级震源；
[0019] 用同样的方式对地震单炮记录完成三随机方案形成不同频带的超炮集。
[0020] 所采用的维纳滤波器为：
[0022] Fkctict是维纳滤波器，Wtoiginal是原始彳目号，WTa_为期望得到的彳目号，ω是角频率， ε选取10 6。
[0023] 所述第二步包括：
[0024] 由真实模型大尺度高斯平滑后得到初始速度模型；
[0025] 采用高阶交错网格有限差分模拟波动方程在地下介质中的传播，具体如下：
[0026] 采用GPU计算时把二维的整个求解区域视为一个核函数，把二维的整个求解区域 分割成多个子区域，每个子区域的大小为16 X 16网格节点，每个子区域映射到一个CUDA的 线程块，每个网格点采用线程块的单个进程独立更新，在每个时间步采用GPU中的并行线 程同时更新所有节点值，从而完成整个区域的波场计算；
[0027] 对于三维数据，把三维数据依照最慢维的方向排列成二维片，然后按照与二维数 据处理相同的方式，把三维节点映射到相应的进程，完成整个区域的正演模拟。
[0028] 所述第三步中的终止条件如下：
[0030] 其中，dk表示第k次迭代的波场残差的目标泛函值，dk i表示第k-Ι次迭代的波场 残差的目标泛函值，ε为一个很小的量。
[0031] 所述第四步包括：
[0032] 采用伴随方程计算检波器出的反传波场；
[0033] 在计算正向传播的震源波场时把各个时刻内部区域与PML区域交界的质点振动 速度分量的波场值及最大时刻所有波场值存储在GPU显存中，在正向传播的最大时刻把第 三步中的存储数据由GPU显存传递到CPU内存；在逆时反传检波残差波场的过程中，把CPU 存储波场拷贝到GPU显存，根据第三步计算的超炮集的残差基于惠更斯原理对震源波场近 似重构，重构震源波场和计算检波波场时均采用多GPU并行计算技术加速。
[0034] 与现有技术相比，本发明的有益效果是：利用多GPU并行计算技术加速时间域波 动方程系统的数值模拟，采用三随机方案的震源编码技术减少震源的模拟个数，使时间域 全波形反演的效率大大提升，针对传统利用GPU进行全波形反演的数据输入输出引起的计 算效率低下问题，本发明采用有效存储边界策略，最大程度降低了数据的输入输出，从而使 全波形反演的计算效率得到进一步提升，最终为时间域全波形反演的正确实施提供了一种 高效可行的技术方案。
【附图说明】
[0035] 图1 :本发明方法的步骤框图
[0036] 图2 :单核GPU对二维弹性波动方程进行有限差分正演模拟的加速效果对比图
[0037] 图3 :单核GPU对三维弹性波动方程进行时间域有限差分正演模拟的加速效