基于二阶系统的NCPML吸收边界三维声波逆时偏移方法

李青阳,吴国忱,杨凌云,王玉梅

(1.中国石油大学(华东)地球科学与技术学院,山东青岛266580;2.海洋国家实验室海洋矿产资源评价与探测技术功能实验室,山东青岛266071;3中国石油化工股份有限公司胜利油田分公司物探研究院,山东东营257022)

基于双程波动方程的逆时偏移方法可以处理不同种类的波动,如反射波、折射波、绕射波等,且不受成像倾角的限制,具备对高陡构造进行高精度成像的能力[1-5]。逆时偏移在实际应用中仍然受到一些因素的限制,由于波场延拓采用双程波动方程,故对于存储需求及I/O带来很大的负担。关于计算效率和内存优化的方法主要有两大类,一类是高性能计算,另一类是改进算法。FOLTINEK等[6]给出了高阶有限差分逆时偏移在GPU上的实现方法;SUH等[7]研究了基于集群的逆时偏移并行计算实施策略;MCGARRY等[8]通过计算域边界上检波点记录的波场来近似重建震源波场。

正演模拟是逆时偏移的基础,正演算法的优化能够有效提高计算效率,采用波动方程方法进行地震波正演时,一般都需要增加人工边界[9],将无限区域截断成有限区域。目前有多种吸收边界方法,例如基于傍轴近似的单程波动方程吸收边界条件[9-10]、海绵吸收边界条件[11-12]以及完全匹配层(perfectly matched layer,PML)吸收边界条件[13]等。在以上吸收边界条件中,完全匹配层吸收边界条件具有更好的吸收效果。经典PML是由复坐标伸展变换(complex coordinate stretching,CCS)推导给出的[14-15]。王守东[16]给出了声波模拟中PML吸收边界条件的基本原理。殷文等[17]将PML吸收边界应用于频率域弹性波正演。王永刚等[18]研究推导了二维声波方程PML吸收边界条件,并与旁轴近似吸收边界条件进行了比较,证明PML吸收边界条件更好。陈可洋[19]对传统的声波完全匹配层进行改进,提出了改进的声波分裂PML吸收边界条件。COLLINO等[20]对一阶方程实现PML的应用,KOMATITSCH等[21]将其扩展到二阶方程,并推导出分裂的PML公式(SPML)。然而CCS存在一定的缺陷,即无法吸收高角度的入射波和瞬逝波。KUZUOGLU等[22]提出复频移(complex frequency shift,CFS)PML,能够对高角度入射波进行衰减,BÉRENGER[23-24]同样验证了CFS-PML吸收边界条件对瞬逝波具有显著的吸收效果。但是CFS-PML边界条件需要对微分系统引入过多辅助变量,求解该系统需要计算大量卷积或者对波场进行分裂,计算代价巨大。秦臻等[25]在频率域弹性波模拟中给出辅助变量推导了微分形式下的CFS-PML吸收边界条件。熊章强等[26]提出了针对一阶速度-应力声波方程数值模拟CFS-PML吸收边界条件。RODEN等[27]提出卷积PML(CPML)吸收边界条件并利用递归方法求解CFS-PML中的卷积项[28]。KOMATITSCH等[29]用CPML解决一阶速度-应力系统弹性波方程。杨凌云等[30]提出了一种新的二阶系统下声波方程的CPML(NCPML)吸收边界条件,大大减少了内存使用,提高了计算效率。张鲁新等[31]将完全匹配层应用于孔隙介质中。PML吸收边界条件被引入到多孔弹性介质[32-33],各向异性介质[34-37]以及黏弹性介质[38-39]。PASALIC等[40]给出基于CPML吸收边界条件的二阶各向同性和各向异性声波方程数值模拟方法,DROSSAERT等[41]提出递归积分法(RI-PML),以积分计算替代卷积计算。除了常用的有限差分法数值模拟,李青阳等[42]研究了高阶伪谱法弹性波PML吸收边界条件,提高正演计算精度。覃发兵等[43]研究了完全匹配层在时域有限元弹性波数值模拟中的应用。此外,近似PML[44](nearly PML,NPML)与上述PML对坐标进行拉伸不同,NPML直接对波场进行拉伸变换无需卷积运算。刘守伟等[45]在三维逆时偏移GPU/CPU机群实现中引入NPML吸收边界条件。罗玉钦等[46]提出多轴复频移NPML吸收边界条件,提升了NPML对大角度入射波的吸收能力与NPML的稳定性。

CFS-PML吸收边界条件在均匀各向同性介质下的三维二阶声波方程应用中,SPML吸收边界条件易于实现,但是波场分裂引入了9个波场变量[47](含3个中间变量),常规CPML吸收边界条件通过引入6个辅助变量[40]大大减少了内存消耗。赵茂强[48]通过忽略衰减因子空变特性消去了SPML的中间变量,将新增波场数量降为6个,但并未分析衰减因子忽略带来的影响。杨凌云等[30]进一步将赵茂强的处理思路[48]引入CPML研究中,提出了NCPML吸收边界条件,进一步减少了辅助变量的数量,同时文中分析指出牺牲衰减因子的空变特性后吸收效果完全可以接受。由于杨凌云等[30]提出的是二维声波方程的NCPML吸收边界条件,本文将二阶系统下的NCPML吸收边界推广至三维并成功应用于声波逆时偏移。首先回顾了二阶系统下的SPML吸收条件在声波方程中的应用,然后推导出含NCPML吸收边界条件的三维声波方程,最后进行三维均匀模型正演和SEG/EAGE推覆体模型和实际资料逆时偏移测试,通过对比NCPML吸收边界条件与常规二阶SPML吸收边界条件在数值模拟中的计算效率与内存占用情况,验证了NCPML在计算效率和稳定性上相比常规SPML所具有的优势。

1 方法原理

1.1 二阶系统下的SPML吸收边界条件

时间域二阶常密度声波方程表示为:

(1)

式中:x,y,z为笛卡尔坐标系中三个不同正交方向;P为标量波场;vP为速度参数;f为震源项。

忽略震源项,将方程(1)变换到频率域:

(2)

(3)

将方程(3)代入频率域声波方程(2)得到尺度变换后的声波方程:

(4)

式中:sx,sy,sz是x,y,z3个方向的衰减因子。

通过傅里叶反变换得到的时间域方程右端会产生卷积项,以x方向为例:

(5)

式中:F-1表示傅里叶反变换算子;*表示卷积符号。

传统SPML边界条件是将波场分裂为笛卡尔坐标系下的3个正交波场,其目的是避免对方程(4)进行傅里叶反变换时出现卷积项。忽略衰减因子1/sn空变特性(即令方程(5)右端第二项为0),应用于方程(4)和方程(5),经傅里叶反变换后得到常规SPML吸收边界条件:

(6)

如果不忽略衰减因子,将会出现波场对时间的三阶导数项,并因此引入3个中间变量。本文中SPML吸收边界条件采用赵茂强[48]的研究成果,不讨论含中间变量的情况,并参考了刘友山等[49]对弹性波PML吸收边界条件应用的研究。

1.2 二阶系统下的NCPML吸收边界条件

虽然SPML方法算法简单,计算效率高,但是该方法对于计算机内存要求较高,分裂的波场提高了正演的内存使用量。KOMATITSCH等[29]提出基于一阶系统方程的CFS-PML(简写为:CPML)吸收边界条件,通过直接计算卷积来节省内存。前文提到SPML提出的初衷即是避免计算卷积,显而易见,SPML和CPML本质上是在计算复杂度和计算内存取舍方面处于对立位置。PASALIC等[40]将CPML推广到二阶系统,其中新增6个辅助变量。杨凌云等[30]在二阶系统下忽略了衰减因子空变特性推导出NCPML吸收边界条件,辅助变量的数量进一步减少。以上研究均在二维条件下展开,本文将NCPML推广至三维。

杨凌云等[30]使用递归卷积方法有效地解决了方程(5)中的卷积如下:

(7)

(8)

式中:k表示当前时刻;其它衰减参数为:

(9)

将方程(7)至方程(9)代入方程(1)可以得到三维声波方程NCPML吸收边界条件:

(10)

(10)式中辅助变量仅为3个,进一步降低了二阶声波波动方程数值模拟中内存的使用。表1展示了SPML和NCPML吸收边界条件下一次波场延拓所需要的的内存。

1.3 逆时偏移

在各向同性介质中,逆时偏移的理论基础是Claerbout提出的时间一致性原理[50],即反射面存在于地层内下行波波至时间与某一上行波波至时间相一致之处。成像条件可以表示为震源波场(下行波场)和检波点波场(上行波场)之间的零延迟互相关。

(11)

式中:R,U和D分别表示成像剖面、检波点波场和震源波场;smax和tmax分别表示最大炮数和最大记录时间;ti和sj分别表示初至时间和震源。

表1 内存分配情况统计

逆时偏移包含3个步骤[51]:①震源波场正向延拓;②检波点波场逆向延拓;③应用成像条件。互相关成像的量纲为波场振幅的平方,并且该方法对所有波进行成像,偏移剖面低频噪声较严重,KAELIN等[52]提出了震源归一化互相关成像条件,利用震源波场的照明度归一化互相关成像剖面,不仅可以消除部分低频噪声,而且可以消除震源强度对剖面的影响,获得的偏移剖面数值接近于真实的反射系数。震源归一化互相关成像条件表示为:

(12)

三维逆时偏移涉及大量的波场延拓,其对计算机的计算效率、硬件性能有很高的要求,而单一的PC机不能满足三维RTM计算需求,需要借助集群并行计算和图形处理器(graphic processing unit,GPU)来实现三维RTM。MPI是最常用的编程工具,为多节点并行计算提供了一组良好定义的库函数,CPU集群解决了粗粒度并行问题,可以用于不同炮的正、反传运算。GPU则解决了细粒度并行问题,波场进行有限差分运算中,空间域中每个网格点都是独立的,所有网格点可以细粒度并行计算。本文中逆时偏移利用MPI控制的CPU和GPU协同并行策略实现,如图1所示。第一级为MPI控制命令,控制第二级多个CPU节点并行运算,CPU节点控制第三级多个GPU运算。

图1 MPI控制下CPU/GPU协同并行策略

综上所述,本文采用的三维逆时偏移CPU/GPU集群实现方案如图2所示。NCPML边界条件应用于波场的正向和反向延拓,由于二者延拓方向相反,互相关成像条件在同一时刻互相关,因此本文采用FENG等[53]提出的逐阶递增模拟精度的源波场重建方法。其中,三维RTM框架最外层为炮循环,采用MPI控制CPU集群进行多炮偏移并行运算,每一次循环所计算的炮数为GPU卡数总和。假设每个节点包含b块GPU卡,共有a台CPU机器,则每次循环共算出a×b个炮集数据的成像结果。因为GPU显存较小,因此每一块GPU卡只负责某一炮的波场正反延拓中的差分运算。如果存在GPU型号不同的情况,也需要保证GPU最小显存能容纳一次波场重构和反向延拓所需要的的数组内存之和。所有炮集数据偏移过程中的I/O均在一块共享硬盘上完成,并在主节点进行最后多炮成像结果的叠加。

图2 三维逆时偏移CPU/GPU机群实现方案

2 数值测试

2.1 3D均匀介质模型

为了测试NCPML完全匹配层的吸收衰减效果,建立了一个均匀介质模型。模型大小为4000m×4000m×4000m,地震波速度为3000m/s。空间步长为10m,震源采用主频为15Hz的雷克子波,震源位于坐标(2000m,2000m,2000m)处,时间采样间隔为0.8ms;PML吸收层数为60层;采用时间二阶差分精度和空间十阶差分精度,在扩展坐标变量中,方向的参数如下[54]:

(13)

其中,

(14)

式中:x0和d分别表示NCPML边界的起始位置坐标和厚度;R0是理论上的反射系数[55],这里取1×10-6;n1和n2是控制NCPML层衰减变化的指数因子,通常2≤n1≤5且0≤n2≤1,n3是NCPML层的频移变换尺度的指数因子,通常取n1=2,n2=0,n3=1。在y和z方向上参数选择与x方向类似。

图3a至图3c分别表示未加吸收边界、常规二阶SPML吸收边界和二阶NCPML吸收边界进行波场模拟计算得到的1100ms时刻的波场快照,从图3a可以看出,不加吸收边界时,外行波在人工边界处完全反射,而加载SPML和NCPML吸收边界后的波场(如图3b和3c所示)边界处外行波明显被衰减和吸收。由此可见,常规SPML吸收边界与NCPML吸收边界对外行波都具有良好的吸收效果,在边界处不会产生明显的虚假反射。

图3 1100ms时刻三维波场快照a 未加吸收边界; b 常规SPML吸收边界; c 二阶NCPML吸收边界

图4a显示了分别利用常规二阶SPML吸收边界和本文NCPML吸收边界算法数值模拟检波点(1000m,1000m,0)处地震记录的数值解及其与理论值的对比,该理论值是通过大幅度扩大计算范围以消除边界影响得到数值解近似而来的,并不是通过给定一个具有一定频带范围的子波,然后在三维均匀空间中求取的理论解。两种吸收边界模拟的结果与理论值很接近,误差主要出现在980～1200ms的时间段,我们截取该时间段地震记录单独显示(图4b),NCPML和SPML吸收边界模拟结果均存在误差,二者波形振幅在理论值附近波动,而SPML的误差比NCPML的误差更大,为更加直观地定量分析二者的误差,将常规二阶SPML和二阶NCPML吸收边界模拟的地震记录与理论值求差,结果如图4c所示。由图4c可见,SPML吸收边界模拟记录的误差最大值达到5.2×10-3,而NCPML吸收边界模拟记录的误差最大值仅为2.7×10-3,大约是常规二阶SPML的1/2。因此,NCPML吸收边界对于边界处的吸收效果要优于常规二阶SPML。表2是两种不同边界得到的地震记录的累计绝对误差和累计误差的统计表,即对误差值取平方求和,从表中可以看到,本文提出的NCPML吸收边界比传统的SPML吸收边界的计算精度更高。

图5a和图5b显示了理论值以及利用两种吸收边界条件模拟(对应于图3中红色竖线)得到的1100ms和1200ms时刻的波场值。如图所示,边界附近仍然存在一定程度的反射波,但NCPML吸收边界下的反射波形比SPML吸收边界下的反射波形在振幅上与理论值更接近。因此本文推导的NCPML吸收边界的计算精度比常规SPML更高。

图4 SPML和NCPML吸收边界条件在接收点(2000m,2000m,0)位置处的地震记录(a)、地震记录局部显示(b)以及SPML和NCPML吸收边界条件的地震记录与理论值之差(c)

表2 误差统计表

2.2 三维SEG/EAGE推覆体模型

利用三维SEG/EAGE推覆体模型进行逆时偏移测试,速度模型如图6所示,图6a为真实速度模型,用以合成地震记录,图6b为偏移速度模型。模型大小为16000m×16000m×4800m,x,y和z方向网格间距均为40m。模型测试整体共设置676炮(26×26),炮点均匀分布于地表,起始炮点坐标(0,0),相邻两炮640m间隔向地表两个正交方向延拓。检波点布置于地表,间距为40m,即道距和线距均为40m。观测系统最大偏移距为±4000m,当炮点位置靠近模型边缘时,例如炮点坐标为(0,0),该炮在偏移过程中仅100条接收线,每条线有100道,而炮点位于模型中央时则接收线数达到最大的200条,同时每条接收线包含了200道。模型数据震源子波为主频8Hz的雷克子波。地震记录的总时长为1.6s,采样点数为1600,时间采样间隔1ms。

图5 SPML和NCPML边界条件下1100ms(a)和1200ms(b)时刻波场快照抽道对比

图6 三维SEG/EAGE推覆体模型a 真实速度模型; b 偏移速度模型

三维推覆体模型逆时偏移使用CPU集群共18个节点,其中共66块GPU显卡,GPU型号是Tesla K10C,GPU并行环境为CUDA,MPICH2版本。表3是三维推覆体模型进行676炮的逆时偏移内存和计算耗时结果,包括所有的计算与I/O读写时间,内存对比是指RTM过程中单个CPU节点和单个GPU卡的对比。可以看出使用NCPML吸收边界进行三维逆时偏移时在内存和计算效率上均得到很大优化。

图7a和图7b分别为SEG/EAGE推覆体模型应用SPML和NCPML边界条件的逆时偏移三维剖面。三维推覆体模型的目标层位于中深层的河道处,从水平切片中可以看出两种边界条件下的三维逆时偏移方法均能够对河道处的构造细节准确归位,具有较高的精度。从速度纵向剖面与对应的应用SPML和NCPML边界条件的逆时偏移纵向切片对比,可以看出中深层成像能量均衡,浅层由于震源因素影响成像结果,有稍许噪声,中深层成像精度较高。

表3 三维推覆体模型逆时偏移所需内存和耗时

综上所述,在本节的三维SEG/EAGE推覆体模型逆时偏移模型测试中,应用NCPML边界条件与常规SPML边界条件计算得到的偏移结果精度几乎一致。但是NCPML在内存占用上比SPML边界条件减少了约35%,计算效率也有小幅提高。需指出该数据仅在本文所用测试模型中有效,对于大型三维逆时偏移,随着正演模型增大,NCPML的优势更加明显。

图7 不同吸收边界条件下三维推覆体模型逆时偏移剖面a SPML; b NCPML

2.3 实际应用

选取中国东部某三维工区的叠前地震资料进行试算。数据覆盖区域约120km2,满覆盖区域约为40km2。该区地表起伏平缓,地震资料信噪比较高。所选区块共1944炮,最大偏移距5464m,面元网格25m×25m。炮集采样间隔2ms,通过对炮集滤波处理,地震数据主频集中在20Hz。采用偏移速度分析所得的速度模型如图8所示,最大深度为4000m,网格间距为10m×10 m×10m。测试中所用的偏移距为0～2000m,即单一炮的网格大小为400m×400m×400m。时间采样间隔0.8ms,针对炮集网格稀疏与偏移中模拟网格规则化的不符,我们采用多项式插值法进行观测系统自适应。偏移中使用了所有炮以增加覆盖次数,压制噪声。

数据预处理阶段,由于三维逆时偏移只用到资料中反射波的信息,因此在实施之前去除资料中的面波和异常噪声,在尽量保留有效信号的前提下压制噪声,提供高保真的道集数据。同时,我们对所选工区实际资料做能量补偿以消除目的层纵、横向能量差异,有利于目的层弱信号的成像。

图8三维初始速度模型

三维叠前资料逆时偏移利用GPU集群8个节点,所用GPU型号是M2090,共12个GPU,GPU并行环境为CUDA,节点间的并行环境为MPICH2。表4为利用上述运行环境进行三维逆时偏移所需的内存消耗及时间消耗,这里包括所有的I/O读写时间。由于计算设备数量和运算性能的限制,三维逆时偏移测试的运算效率较慢,因此只进行了NCPML边界条件的应用。

表4 三维实际资料逆时偏移所需内存与耗时

图9是应用NCPML边界条件的三维逆时偏移测试结果,分别为抽取In-line测线和Cross-line测线的逆时偏移结果,可以看出逆时偏移结果深浅层能量均衡,成像精度高,这是使用归一化成像条件的原因,使得中、深层能量得到照明补偿。边界处覆盖次数较低,因而成像结果中噪声较多。浅层由于震源等因素的影响,有一定的噪声。综上所述,NCPML在实际资料应用中计算效率和内存占用都比较令人满意。

图9 不同测线的三维逆时偏移效果a In-Line测线;b Cross-Line测线

3 结论

本文将二阶系统下的NCPML吸收边界条件推广至三维并成功应用于声波逆时偏移。均匀模型正演以及三维SEG/EAGE推覆体模型逆时偏移测试结果表明:NCPML应用效果优于常规SPML,同时计算效率和内存使用方面NCPML比SPML更具优势。实际资料的逆时偏移测试证明了NCPML吸收边界条件在实际应用中的稳定性。此外,针对二阶声波方程提出的NCPML可以推广到更复杂的介质波动方程模拟,例如声学各向异性建模和弹性介质。