基于gpu加速和随机边界的三维弹性波逆时偏移算法

基于gpu加速和随机边界的三维弹性波逆时偏移算法

1.gpu并行计算技术逆时偏移需要精确地求解双向声波或弹性波的波动方程，反向地向外推地震数据，并选择合适的成像条件来精确成像。目前逆时偏移是公认的最能精确成像的偏移技术,原因是其适用于任意复杂介质、不存在任何倾角限制、能够对任意类型的波进行数值模拟成像相比CPU,GPU旨在增强计算机图形处理显示能力,天生具备高并行性和高密集度数值计算的巨大潜能。目前,GPU硬件平台和CUDA编程日渐成熟,越来越多的业内人士开始寻求利用GPU超强的并行计算能力解决三维逆时偏移、三维模型数值模拟、速度建模以及全波形反演等地球物理勘探的前沿问题。目前将GPU并行计算应用于海量地球物理数据运算已取得大规模提速,为大规模实际地球物理资料的逆时偏移实现搭建了平台基础。针对逆时偏移需要巨大临时文件存储量的问题,业内提出了Checkpointing方法、边界存储策略等方法来降低硬盘存储,但仍需存储相当大的数据。贯彻“用计算换存储”策略的随机边界存储方法目前三维逆时偏移仍受限于计算量和存储量巨大等问题而难以推广,基于GPU并行加速和随机边界存储策略的三维弹性波逆时偏移算法研究,提高逆时偏移效率,是一个有很强现实意义的课题。2.三维弹性波方程的逆偏移原理2.1d交错网格数值离散3D各向同性介质中的一阶速度—应力弹性波方程如下:式(1)中,x、y、z表示空间直角空间坐标方向,t代表时间,ρ代表密度。v本文采用3D交错网格对方程(1)进行数值离散,交错网格中速度分量v表1、图1展示了3D交错网格中各弹性波波场分量、弹性参数的空间分布。在3D交错网格空间中对式(1)进行有限差分,得到时间2阶,空间上任意偶数阶的有限差分格式,以式(1)中第1个方程举例:式(2)中,Δx、Δy、Δz代表空间x、y、z方向离散间隔;Δt代表时间离散间隔;i、j、k、n表示空间x、y、z方向和时间上的离散点序号；i2.2受内部弹性介质的传播方向本文通过傅立叶变换和矩阵特征值分析可得稳定性条件在三维各向同性介质中,弹性波的传播速度是与传播方向无关的,因此x,y,z三方向应满足稳定条件是相同的。2.3随机速度随机扩散波场模拟本文利用随机边界存储策略进行逆时偏移,基本思路是:选择合适的随机方程,生成的随机边界,消除弹性波方程延拓中人工边界和自由边界条件下反射波的相干性,使得边界反射不能进行有效成像。具体实现实现过程如下:将有限计算空间向外扩充一定的区域,在外扩的区域内填充随机速度,当弹性波场越过模型边界向外继续传播到随机速度区域时,内充随机速度对波前面进行随机化,使波场信息变成随机噪声,反传回真实的速度区域,破坏模型边界反射的相干性,进而使模型边界反射不能进行有效成像。注意:随机边界机速度也要满足稳定性条件公式(3)。本文使用3D各向同性的地质模型来验证随机边界的效果。三维均匀介质模型大小为500m×500m×500m,采样网格Δx、Δy、Δz为5m,纵波速度v2.4矢量横波标量化据马德堂等,纵波为无旋场,横波为无散场该方法所得纵波为标量,横波则为矢量,需要进行矢量横波标量化,这样获得的标量横波可以直接用于逆时偏移的成像。2.5成像剖面的选取弹性波逆时偏移的另一个关键是偏移成像条件的选取,这将直接影响成像剖面的最终质量。本文采用归一化互相关成像条件,通过计算两个波场的零延迟互相关,并做归一化处理,进行偏移成像。3.gpu并行性分析大型地震资料的逆时偏移因其巨大的计算量和数据存储需求而备受限制,上文分析可得逆时偏移的延拓和成像算法是可分解和可线性叠加的,因而具备大规模并行计算的可能。待移植程序的计算部分是否具有并行性以及计算量大小是程序移植之前必须考虑的2个问题。若不能同时满足这两点要求,即是程序成功移植到GPU平台,也会出现并行性不高,效率没有明显提升的现象,有时甚至会降低。串行计算中三维弹性波逆时延拓由CPU单一线程顺序实现,GPU平台可通过多线程并行执行。GPU内没有(或有少量)高速缓存,事实证明访问共享存储器的速度和访问寄存器的速度相差不多,因此可以通过访问共享存储器来实现不同线程间的最小通信延迟。龙桂华等分析发现GTX295显卡上利用片内共享存储计算1024×1024两矩阵乘法时速度大约为利用片内共享内存时计算速度的4.5倍4.空间大时延拓的随机边界存储策略逆时偏移所采用的成像条件需要对同一时刻的炮点和检波点波场进行相关计算,而炮点和检波点波场是相互反向的过程,不能同时由计算,因此常规逆时偏移通常选择在炮点波场的正向延拓时记录整个波场的值,在检波点波场的逆时延拓过程中直接读入所需同一时刻的震源波场值。这种需要存储整个模型波场值的方法硬盘存储量是巨大的,大约x×y×z×t。而巨大硬盘存储要求是制约3D资料弹性波逆时偏移推广的又一瓶颈。另外,震源波场正向延拓记录数值、检波点波场逆时延拓读取数值等操作都将导致I/O访问过大,耗费大量时间,严重影响计算效率。对于CUDA这种异构编程模式而言,CPU/GPU通信是必须考虑的问题,需要最大限度地减少I/O访问、读取操作。因此,本文采用随机边界存储策略,在炮点波场正向延拓过程中,使用随机速度构建随机边界,只需保存最后t时刻的炮点波场值(包括随机边界范围在内的全部炮点波场值),在检波点波场进行逆时延拓的同时,通过读取保存的最后t时刻随机波场快照和随机边界同步逆时重构炮点波场。这样,相同时刻的炮点波场值和检波点波场值能随时计算出来,进而进行偏移成像。虽然这种方法对炮点波场延拓进行了2次计算,但最大限度地降低了硬盘空间存储量和文件读写、拷贝的工作量,相比于GPU极高的并行计算能力,增加的计算量可忽略不计。5.地表剖面位置及三维倾斜模型三维倾斜模型,底部有一高速圆形隆起,模型1500m×1500m×1500m,采样网格Δx、Δy、Δz为5m,时间间隔为0.5ms,记录1.5s。一共布设10条测线,测间距为150m,测线内炮间距为100m,每条测线布设15炮,每炮30×30道接收,处理面元50m×50m,在地表处(z=0)接收。三维倾斜模型的纵横波速度模型以及逆时偏移剖面如图3所示。结果显示3D弹性波逆时偏移算法能够比较准确对模型进行成像,识别出其中的水平地层、倾斜地层和底部穹形隆起。本次3D倾斜模型一共正演150炮,逆时偏移150炮。比较一炮数据计算时间,CPU算法:6301.62s,GPU并行+随机边界算法(GTX560Ti,核心448):137.41s,加速比达45.86,计算效率