果对 比图
[0038] 图 4-1 :真实 Marmousi2 模型
[0039] 图4-2 :初始纵波速度模型
[0040] 图4-3 :声波方程5. OHz多尺度反演结果，串扰噪音严重
[0041] 图4-4 :声波方程15. OHz多尺度反演结果串扰噪音得以压制
[0042] 图5-1 :初始的Overthurst模型（拉梅常数lambda)
[0043] 图5_2 :初始的Overthurst模型（拉梅常数mu)
[0044] 图5-3 :弹性波动方程5. OHz多尺度反演结果（拉梅常数lambda)
[0045] 图5-4 :弹性波动方程5. OHz多尺度反演结果（拉梅常数mu)
[0046] 图5-5 :弹性波动方程10.0 Hz多尺度反演结果（拉梅常数lambda)
[0047] 图5-6 :弹性波动方程10.0 Hz多尺度反演结果（拉梅常数mu)
[0048] 图5-7 :弹性波动方程15. OHz多尺度反演结果（拉梅常数lambda)
[0049] 图5-8 :弹性波动方程15. OHz多尺度反演结果（拉梅常数mu)
【具体实施方式】
[0050] 下面结合附图对本发明作进一步详细描述：
[0051] 本发明一种基于震源相位编码和多GPU并行计算技术且无数据I/O的高效时间域 波动方程全波形反演方法，为大规模数据的高精度速度建模提供一种高效的处理手段。
[0052] 如图1所示，本发明方法包括：
[0053] 第一步、采用源极性随机、源位置随机、源个数随机三随机方案对原始观测地震记 录随机形成不同频带超炮集；
[0054] 具体如下：
[0055] 采用MPI并行输入原始单炮记录和原始雷克子波，利用维纳滤波器对原始单炮记 录和原始雷克子波滤波，得到不同频带的地震单炮记录和地震子波，所采用的维纳滤波器 为：
[0057] 其中Fkotct是维纳滤波器，Wtoiginal是原始彳目号，W lallif3t为期望得到的彳目号，ω是角 频率，ε是一个很小的值，为了避免不稳定现象的发生，本发明选取1〇 6。利用C语言库函 数中产生随机种子的srandO函数、随机函数randO以及符号函数对所有地震子波赋予 [-1，1]之间的一个数，即赋予每个地震子波不同的颜色，根据地震子波颜色的正负性，赋予 地震子波相应的极性，此时完成源极性随机；
[0058] 采用堆排序算法或快速排序法对改变源极性的地震子波的颜色排序，颜色值逐渐 增大，并采用C语言中的随机函数rand ()根据总超级震源子波个数及总单个地震子波个数 随机出每个超级震源中单个地震子波的个数，此过程为源个数随机；
[0059] 此后按照每个超级震源的单个地震子波的个数对排序后改变源极性的单个地震 子波叠加，此时叠加的每个地震子波的位置也相应随机，此过程为源位置随机，由此完成 三随机方案形成不同频带超级震源，用同样的实现方式对地震单炮记录形成不同频带超炮 集。
[0060] 第二步、采用基于CUDA的多GPU并行计算技术加速时间域一阶形式的声波方程和 弹性波动方程系统的高阶交错网格有限差分正演模拟。
[0061] 第二步为正演模拟人工合成超炮集并采用第一步得到的超炮集与模拟的超炮集 做残差，其中采用的子波为超级震源子波，初始模型由真实模型大尺度高斯平滑后得到。
[0062] 正演问题是反演问题的基础，本发明采用的一阶形式的声波方程系统为：
[0064] 其写成算子形式为：
[0068] 本发明采用的一阶形式的弹性波动方程系统为：
[0072] 其写成算子形式为：
[0078] 本发明采用高阶交错网格有限差分模拟波动方程在地下介质中的传播，在波场外 推时，每个网格节点的波场在更新时用到该节点前一个时刻的波场值和周围节点前半个时 刻的波场值，该问题具有很好的并行性，在采用GPU计算时把二维的整个求解区域（第一 步得到超炮集，第二步是用GPU模拟超炮集，目的是求两者的残差）视为一个核函数，把二 维的整个求解区域分割成多个子区域，每个子区域的大小为16X16网格节点，每个子区域 映射到一个CUDA的线程块，每个网格点采用线程块的单个进程独立更新，在每个时间步采 用GPU中的并行线程同时更新所有节点值，从而完成整个区域的波场计算。对于三维问题， 把三维数据依照最慢维的方向（即笛卡尔坐标系的y轴方向）排列成二维片，基于二维问 题的处理方式划分进程块和进程，把三维节点映射到相应的进程，完成整个区域的正演模 拟。由于各个震源（一个震源对应一个节点）的数值模拟相互独立，因而在正演模拟时可 以采用CPU/GPU协同技术，用多GPU核模拟不同的震源传播。
[0079] 第三步、利用GPU加速计算震源波场时（即第二步）把震源波场的边界数据和最 后一个时刻所有波场存储于计算机内存（即图1中的"正演（CPU)存边界波场"）；
[0080] 第四步、采用多GPU并行计算技术计算检波器反传波场，同时采用有效存储边界 策略（有效边界存储策略是一种重构震源波场的办法，该办法需要用到存储的数据，常规 的时间域全波形反演中重构震源波场采用CPU，本发明采用GPU。）和多GPU并行计算技术 重构震源波场（即把正常传播过程中存储的每个时刻的边界波场作为边界条件，存储的最 后一个时刻的所有波场作为初始条件，然后解反向传播的波动方程，这是一个数学物理方 程的求解问题，解完即可重构波场。）并计算目标泛函相对于模型参数的梯度（梯度可以看 做震源波场和伴随波场的互相关，当两个波场解出后，梯度也相应可解出。）；
[0081] 本发明采用有效存储边界策略重构震源波场，在计算正向传播的震源波场时把各 个时刻内部区域与PML(最佳匹配层）区域交界的质点振动速度分量的波场值（存储的边 界层数为空间精度差分的一半即可，在重构波场时采用低阶差分向高阶差分过度的思想） 及最大时刻所有波场值存储在GPU显存中（这个存储是采用GPU计算的时候需要申请的 GPU内存，当得到最大时刻的波场后，需要把GPU内存中需要的数据传递到CPU内存中，可以 看做是第三步的中间过程），在正向传播的最大时刻把上述存储数据由GPU显存传递到CPU 内存。在逆时反传检波残差波场的过程中，把CPU存储波场拷贝到GPU显存，根据第二步计 算的超炮集的残差基于惠更斯原理对震源波场近似重构，重构震源波场和计算检波波场时 均可采用多GPU并行计算技术加速。该过程，仅需要对GPU的显存和CPU的内存做两次数 据传递，且不需要对数据存储和读取，相较于传统需要数据存储和读取的逆时偏移和全波 形反演，仅多了对内部区域震源波场重构这一过程，该过程可采用不加 CPML边界条件（前 边的方程都是采用基于CPML (卷积最佳匹配层-)的声波方程和弹性波方程，以压制边界的 反射波，而重构震源波场时，采用的是不加 CPML的声波方程和弹性波方程）的声波方程和 弹性波动方程模拟，当采用多GPU时其计算效率大大提升。重构波场后，与伴随波场互相关 可以求解梯度，梯度求解出后可以采用优化算法完成速度模型迭代更新，即全波形反演。
[0082] 本发明采用伴随方程计算检波器出的反传波场（即图1中的伴随波场），对一阶速 度-应力声波方程,其表达式为：
[0083] Lt λ = g
[0084] 其中：
[0085] 对一阶速度-应力弹性波动方程，采用的伴随方程为：
[0086] Lt λ = g
[0087] 其中：
[0089] 本发明采用的声波方程目标泛函关于模型参数的梯度的表达式为：
[0091] 其中NS表示采用三随机方案后得到的超级震源的个数。
[0092] 本发明采用的弹性波动方程目标泛函关于模型参数的梯度表达式为：
[0095] 第五步、基于三随机方案、多重网格法和维纳滤波器技术实现高效时间域波动方 程多尺度全波形反演。
[0096] 采用