三维逆时偏移中的随机边界条件

三维逆时偏移中的随机边界条件

0多波地震资料偏移处理的意涵基于标量波理论的反射纵波勘探技术在过去的能源开发过程中发挥了重要作用。随着油气开发对地震勘探精度的要求的不断提高，勘探目标结构和岩性复杂性的不断增加，单因素扩散的纵波采集精度的提高越来越受到其理论上的局限性和波场信息的个体化。碳酸盐岩裂缝储层、煤层和岩浆岩储层的勘探难以以德容忍。近年来,以矢量波动理论为基础的多波地震勘探技术引起了业界的广泛重视,同常规纵波勘探技术相比,多波地震技术的理论假设与实际情况更为接近,且横波在复杂介质中具有独特的传播机理,故理论上多波地震技术解决实际问题的能力更强,更有利于提高勘探精度或实现某些特殊目标的精确勘探。叠前深度偏移是多波地震技术的重要研究内容,对多波地震资料进行叠前深度偏移处理的意义在于:一方面,地下构造的精确成像需要借助高精度的叠前偏移技术来实现;另一方面,多波偏移剖面中所包含的地震波动力学特征还可用于岩性和流体识别;同时偏移得到的共成像点道集还可为振幅随偏移距变化(AVO)、振幅随入射角变化(AVA)或其它叠前反演工作提供输入数据。因此,多波地震资料偏移应当实现两个目标:(1)地下构造的纵、横波准确成像;(2)获取保真的纵横波叠前反演道集和叠加剖面。基于弹性波理论的多分量地震资料联合逆时偏移技术是实现多波地震纵横波偏移成像的有效工具逆时偏移思想提出之初,在叠后地震资料偏移成像中获得了良好的成像剖面目前的野外采集普遍向高密度、三维方向发展,为确保复杂构造观测波场的连续性,多采用小炮距、小接收线距、长排列接收,部分探区单炮记录过万道本文针对有效边界存储在三维情况下存储过大的情况,提出应用随机边界条件,在不保存边界波场的情况下,实现正向延拓波场的重构;考虑三维逆时偏移情况下,开辟数组空间内存较大,提出应用区域分解技术将计算数据分配到不同节点,实现基于MPI+CUDA的协同并行。1纵波场叠加响应的基本原理弹性波逆时偏移将多分量地震记录当作弹性波场的边值问题弹性波逆时偏移的基本原理为:假定地下介质由一系列绕射点组成,地面、海底或井中接收到的三分量记录为各绕射点产生的纵横绕射波场在接收点处的叠加响应,其中各绕射点纵横绕射波的产生时间为炮点波场的主能量到达该点的时间1.1时间上1阶差分根据各向同性介质中,位移-应变关系、应力-应变关系、应力-位移关系其中,υ在交错网格空间中对式(1)进行差分离散可得到时间2阶,空间任意偶数阶的差分格式,以式(1)中第一个方程为例:其中,Δx、Δy、Δz分别表示坐标系中x、y、z方向上的网格间距;Δt表示时间采样间隔;i,j,k分别代表x、y、z方向上离散点的序号;n表示时间上的离散点的序号。且ρ、μ半网格点上的值由同一方向上相邻整网格点求取算术平均得到。2N为差分阶数,D1.2无分裂pml边界弹性波方程交错网格有限差分格式属于显示差分型递推格式,计算当前时刻波场值时,需利用上一时刻波场值。由于数值计算截断误差的影响,随着时间递推误差会累积,破坏数值稳定性甚至引起计算溢出通过对三维各向同性介质一阶双曲型弹性波方程进行平面谐波分析式(4)中,Δx、Δy、Δz分别表示坐标系中x、y、z方向上的网格间距;Δt表示时间采样间隔;υ采用无分裂式PML边界条件解决阶段边界的伪反射问题上式中Ω其中,d(x)、d(y)、d(z)为三个方向的衰减因子;根据公式(5)、(6)可依此类推基于无分裂PML边界条件的弹性波控制方程,其计算步骤为:(1)对包含PML边界的整个区域进行常规的弹性波波场求解;(2)在PML边界区域将得到的常规波场减去衰减量,然后按照时间轴方向更新辅助变量。由控制方程和计算步骤可知,无分裂的PML边界条件避免了对弹性波波场进行分裂,降低了存储。由于采用同一套方程也不会出现在边界处两套方程耦合的编程难题。1.3偏移剖面的基本原理成像条件决定着逆时偏移成像剖面的质量。目前业界广泛应用互相关成像条件,其基本原理为:利用炮点波场和接收点逆时延拓波场做零延时互相关,在反射点位置的波场将得到加强,非反射点处的波场将削弱,最后所有时刻波场相加得到偏移剖面。互相关成像条件的实现算法为:其中,I2横波场的重构上述弹性波逆时偏移基本流程中,在利用互相关成像条件时需同时应用相同时刻的检波点波场和炮点波场,上述两个波场的时间序列相反,因此需要重构炮点波场以实现纵横波成像;本文采用随机边界条件散射人工边界的伪反射,实现波场的重构。为了解决三维情况下单卡GPU显存不足的难题,利用区域分解技术将数据划分,实现大规模数据的三维弹性波逆时偏移。2.1机速度无线性的影响在解决炮点、检波点波场时间序列不一致的问题中,国内外学者先后提出Checkpointing方法、有效边界存储策略随机边界条件实质为在计算区域外部随机赋速度值,利用随机速度的无规律性,将波场值散射,这样在进行互相关成像和叠加时,无规律的波场值能量得到减弱,反射点成像值相对得到增强。所以,随机边界条件的关键之处在于随机边界层内随机速度的选择,常用随机速度赋值公式为:式中,V为了验证随机边界条件的适用性,本文采用二维均匀速度模型,验证随机边界的重构和散射效果。二维均匀速度模型大小1500m×1500m,纵波速度为2500m/s,横波速度1800m/s,炮点位置位于(750m,750m),纵波源激发,记录不同时刻v2.2区域分解和多卡并行数据传递GPU中有大量晶体管,是专为执行复杂的数学和几何计算而设计的。当今的GPU已不再将这些计算局限于图形渲染,其通用计算技术的发展和天生的并行处理机制在计算方面提供了高于CPU数十倍乃至上百倍于CPU的性能。业界已经利用GPU的并行特性在二维逆时偏移取得了一系列突破在对弹性波方程的差分格式,进行有限差分计算编程时,将一个维度上(以y方向为例)的区域进行划分,按照实际资料的大小和显卡显存的限制来决定划分的数量。利用区域分解进行差分计算,需要解决的核心问题在于如何隐藏卡与卡间的数据交换。如图5所示,为了保证差分计算的正确性,在每一个时间步需要在虚线区域进行数据交换。以其中的两块卡为例,说明分解区域后每块卡的计算与数据传递。如图6所示,每块卡可以分为三部分(图中分别用A、B和C表示),C区域为两块卡之间的重叠部分,厚度由差分格式来决定,B区域为传递部分与C区域保持一样的厚度,A区域为显卡本身的计算区域。在正演和延拓计算的每一个时间步内,C区域的网格处在边界而无法进行差分运算,所以需要将B区域的数据传递到C。这样计算过程就分为了如下步骤:首先是每块卡计算B区域;然后,每块卡计算A区域的同时,传递B区域数据到相邻显卡。计算能力高于2.0的NVIDIA显卡支持在同节点上进行显卡与显卡之间的数据交换。而本文提出的多卡并行是基于多节点的并行,在数据传递时需要将传递数据通过内存中转传递,由MPI的每一个进程控制每块显卡,在MPI进程组内进行非阻塞的数据交换。整个的实现流程图如图7所示:(1)由主进程读取到炮记录和速度模型并分发于进程;(2)各进程控制GPU显卡实现每一时间步的图4B区域波场计算;(3)GPU回传数据并在内存中交换,同时实现图4A区域的波场计算;(4)循环迭代波场,最后由主进程规约得到最终的偏移结果。3盐丘模型三维弹性波逆时偏移试验SEG/EAEG盐丘模型是一个国际上通用测试三维复杂构造成像效果的地质模型。本文采用SEG/EAEG盐丘速度模型作为纵波模型,通过固定的纵横波速度比构造横波模型,炮集记录由一阶弹性波方程有限差分正演获得。模型参数如下:模型大小3380m×3380m×1050m,网格大小Δx=Δy=Δz=5m,时间采样间隔0.35ms,记录长度2.8s,震源为35Hz雷克子波。观测系统参数为:总炮数4690炮,67条测线,炮线间距50m,炮间距60m,接收线间距10m,接收道间距5m,每条测线70炮,每炮30400道接收;偏移算法采用本文所提出的基于MPI+CUDA的三维弹性波逆时偏移并行算法。图8为三维盐丘模型的纵横波速度模型以及多分量逆时偏移剖面。另外,PP、PS的成像效果较SP、SS的成像效果更好。原因为:PP,PS波是在纵波倾斜入射到弹性界面时产生能量较强,而SP、SS是由纵波倾斜入射到弹性界面产生的转换横波再次入射得到的地震记录,该情况涉及波型的多次转换,情况复杂。因此PP成像结果主要代表了P波的反射系数,PS主要成像结果代表了PS波的反射系数;而SP、SS成像结果并不对应一种波的反射系数,多次转换后并不清楚它的能量分布。PP波和PS波均可以对盐丘模型构造得到准确的偏移成像结果。PP波主频比PS波主频低,其原因在于:模型资料的纵横波记录是由同一纵波源激发得到,且转换横波在空间域中的波长小于纵波波