梯度类迭代算法更新速度模型，对更新后的速度模型计算超级炮记录并与第 一步的超级炮记录做残差，当残差满足终止条件：
时该频段的全波形反演迭代结束，并把反演的模型做为下一个频 带全波形反演的初始模型，重复第二步到第五步，当满足给定的频率终止条件即：给定的反 演频带个数时完成最终的多尺度全波形反演流程，输出当前模型做为全波形反演的最终结 果。
[0098] 采用基于震源编码的超级震源全波形反演时，反演结果中会存在串扰噪音，并且 超级震源的引入使反演问题的非线性程度进一步加深，为压制随机串扰噪音，减小反演的 非线性性，本发明采用基于多重网格法的多尺度时间域全波形反演算法（即第五步中的反 演依次从低频段到高频反演的思想。），反演从大网格逐步精细化，从背景异常体逐渐精细 到目标体。
[0099] 图2为当显卡采用C2060时，采用单个GPU卡对二维不同网格尺度下的声波正演 模拟进行加速的一个加速比对比图，菱形折线表示采用CPU模拟声波传播的计算时间，方 块折线表示采用GPU模拟声波传播的计算时间，三角折线表示加速比。当采用多卡时，加速 比成线性关系增加，图3为当显卡采用C2060时，采用单个GPU卡对三维不同网格尺度下的 声波正演模拟进行加速的一个加速比对比图。图4-1至图4-4这为采用该实例对Marmousi2 模型采用声波全波形反演得到的速度剖面。图4-1为准确的Marmousi2速度模型，图4-2 为经过高斯大尺度平滑后的初始速度模型，图4-3是主频为5Hz的超炮集反演获得的反演 结果，从结果中可以看出模型中存在着较强的串扰噪音，把图4-3的结果作为初始模型，并 以此对高频带超炮集进行全波形反演，图4-4是主频为15Hz的超炮集反演所得速度模型， 从模型中可以看出，串扰噪音得到较好压制。图5-1至图5-8是采用该实例对Overthrust 模型采用弹性波全波形反演同时反演得到的拉梅常数剖面，其中图5-1，5-2为初始的拉梅 常数，图5-3,图5-4为采用主频为5. OHz的超炮集反演所得结果，把该结果作为初始模型对 10.0 Hz的超炮集反演，图5-5,图5-6是10.0 Hz的超炮集反演所得结果，把该结果作为初始 模型对15Hz的超炮集进行反演，图5-7,图5-8为其相应反演结果。
[0100] 巨大的计算量和存储量一直是制约3D弹性波动方程全波形反演发展的重要因 素。针对时间域波动方程全波形反演计算量大和存储量大的问题，本发明分别采用基于 CUDA的多GPU并行计算技术和有效边界存储策略予以解决，并结合三随机相位编码的多震 源技术对时间域弹性波动方程全波形反演加速处理，为全波形反演的准确实施提供了一套 高效可行的技术方案。
[0101] 上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本 发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本 发明上述【具体实施方式】所描述的方法，因此前面描述的方式只是优选的，而并不具有限制 性的意义。
【主权项】
1. 一种高效的时间域全波形反演方法，其特征在于：所述方法包括： 第一步、采用源极性随机、源位置随机、源个数随机Η随机方案对原始观测地震记录随 机形成不同频带超炮集； 第二步、采用基于CUDA的多GPU并行计算技术加速时间域一阶形式的声波方程系统和 弹性波动方程系统的高阶交错网格有限差分正演模拟，得到GPU模拟第i化eq频带超炮集； 在正演模拟时采用CPU和GPU协同技术，用多GPU核模拟不同的震源传播，把震源波场的边 界数据和最后一个时刻所有波场存储于计算机内存；在计算正向传播的震源波场时把各个 时刻内部区域与PML区域交界的质点振动速度分量的波场值及最大时刻所有波场值存储 在GPU显存中； 第Η步、计算波场残差及目标泛函值，然后判断是否满足终止条件，如果是，表示该频 段的全波形反演迭代结束，并把反演的模型作为下一个频带全波形反演的初始速度模型， 转入第六步，如果否，则转入第四步； 第四步、采用多GPU并行计算技术计算检波器反传波场，同时采用有效存储边界策略 和多GPU并行计算技术重构震源波场并计算目标泛函相对于模型参数的梯度； 第五步、采用梯度类迭代算法更新速度模型，然后返回第二步； 第六步、判断是否满足频率终止条件，如果是，则转入第走步，如果否，则返回第一步； 第走步、输出当前模型。2. 根据权利要求1所述的高效的时间域全波形反演方法，其特征在于；所述第一步是 送样实现的： 采用MPI并行输入原始单炮记录和原始雷克子波，利用维纳滤波器对原始单炮记录和 原始雷克子波滤波，得到不同频带的地震单炮记录和地震子波； 利用C语言库函数中产生随机种子的函数、随机函数W及符号函数对所有地震子波赋 予[-1,1]之间的一个数，即赋予每个地震子波不同的颜色，根据地震子波颜色的正负性， 赋予地震子波相应的极性，此过程为源极性随机； 采用堆排序算法或快速排序法对改变源极性的地震子波的颜色排序，颜色值逐渐增 大，并采用C语言中的随机函数根据总超级震源子波个数及总单个地震子波个数随机出每 个超级震源中单个地震子波的个数，此过程为源个数随机； 按照每个超级震源的单个地震子波的个数对排序后改变源极性的单个地震子波叠加， 此时叠加的每个地震子波的位置也相应随机，此过程为源位置随机，由此完成Η随机方案 形成不同频带的超级震源； 用同样的方式对地震单炮记录完成Η随机方案形成不同频带的超炮集。3. 根据权利要求2所述的高效的时间域全波形反演方法，其特征在于；所述第二步包 括： 由真实模型大尺度高斯平滑后得到初始速度模型； 采用高阶交错网格有限差分模拟波动方程在地下介质中的传播，具体如下： 采用GPU计算时把二维的整个求解区域视为一个核函数，把二维的整个求解区域分割 成多个子区域，每个子区域的大小为16X16网格节点，每个子区域映射到一个CUDA的线程 块，每个网格点采用线程块的单个进程独立更新，在每个时间步采用GPU中的并行线程同 时更新所有节点值，从而完成整个区域的波场计算； 对于Η维数据，把Η维数据依照最慢维的方向排列成二维片，然后按照与二维数据处 理相同的方式，把Η维节点映射到相应的进程，完成整个区域的正演模拟。4. 根据权利要求3所述的高效的时间域全波形反演方法，其特征在于；所述第Η步中 的终止条件如下：其中，dk表示第k次迭代的波场残差的目标泛函值，dk1表示第k-1次迭代的波场残差 的目标泛函值，ε为一个很小的量。5. 根据权利要求4所述的高效的时间域全波形反演方法，其特征在于；所述第四步包 括： 采用伴随方程计算检波器出的反传波场； 在计算正向传播的震源波场时把各个时刻内部区域与PML区域交界的质点振动速度 分量的波场值及最大时刻所有波场值存储在GPU显存中，在正向传播的最大时刻把第Η步 中的存储数据由GPU显存传递到CPU内存；在逆时反传检波残差波场的过程中，把CPU存储 波场拷贝到GPU显存，根据第Η步计算的超炮集的残差基于惠更斯原理对震源波场近似重 构，重构震源波场和计算检波波场时均采用多GPU并行计算技术加速。
【专利摘要】本发明提供了一种高效的时间域全波形反演方法，属于石油地震勘探数据高效率、低存储、高精度速度建模领域。第一步、采用源极性随机、源位置随机、源个数随机三随机方案对原始观测地震记录随机形成不同频带超炮集；第二步、采用基于CUDA的多GPU并行计算技术加速时间域一阶形式的声波方程系统和弹性波动方程系统的高阶交错网格有限差分正演模拟，得到GPU模拟第ifreq频带超炮集；在正演模拟时采用CPU和GPU协同技术，用多GPU核模拟不同的震源传播，把震源波场的边界数据和最后一个时刻所有波场存储于计算机内存。
【IPC分类】G01V1/28
【公开号】CN105319581
【申请号】CN201410373265
【发明人】王杰, 王立歆, 方伍宝, 胡光辉, 孙晶梅, 王振宇, 尹力
【申请人】中国石油化工股份有限公司, 中国石油化工股份有限公司石油物探技术研究院
【公开日】2016年2月10日
【申请日】2014年7月31日