基于Qt的Walkaway VSP逆时偏移软件开发研究

王维红,柯 璇,郭雪豹,张莹莹,于 鹏

(1.东北石油大学地球科学学院,黑龙江大庆163318;2.中国石油化工股份有限公司东北油气分公司,吉林长春130062)

基于Qt的Walkaway VSP逆时偏移软件开发研究

王维红1,柯 璇1,郭雪豹1,张莹莹1,于 鹏2

(1.东北石油大学地球科学学院,黑龙江大庆163318;2.中国石油化工股份有限公司东北油气分公司,吉林长春130062)

VSP观测系统的检波器位于井筒内,地震资料的信噪比高,波场信息丰富,利用基于双程波的叠前逆时偏移方法可实现VSP资料的高精度成像,也为井筒附近复杂介质和目的层流体性质的研究提供基础资料。在Walkaway VSP资料逆时偏移算法研究的基础上,在Linux操作系统下,以Qt环境为开发平台,设计了一套较为完整的VSP资料逆时偏移处理流程,初步开发了Walkaway VSP资料逆时偏移计算软件。该软件具有CPU和GPU两套核心算法模块,可适应不同的硬件系统,理论模型数据测试结果表明,该软件平台具有操作性强、移植性好和成像精度高的特点。

VSP;逆时偏移;Qt;软件界面

逆时偏移(RTM)理论自Whitmore等1983年提出之后[1-3],许多地球物理工作者对其进行了深入的研究,但是由于其计算量巨大,而且在计算过程中产生大量中间数据,对计算机的计算能力和存储性能都有很高的要求,所以在问世之初未能得到广泛应用。近年来,随着计算机硬件技术的发展,逆时偏移算法也得到很好的发展和应用。逆时偏移本身基于双程波方程,同时由于直接对波动方程进行离散,因此能够准确模拟地下复杂构造的波场特征,从而大幅度提高地震资料的成像精度。在波场数值模拟过程中,常用的方法包括有限差分、有限元和伪谱法[4-6],其中有限差分法由于其实现相对简单而得到普遍应用。在边界条件方面,主要有随机边界和吸收边界条件等。随机边界的引入主要是为了解决逆时偏移过程中的波场数据存储问题,但实际上也会将随机噪声引入到成像过程中[7];吸收边界中最常采用的是完全匹配层吸收边界条件,其吸收效果好,理论上能够吸收来自任何方向的入射波,但是计算效率较低,同时在偏移计算过程中产生的临时波场数据量较大,一般的硬件系统限制该方法的有效应用。

随着石油勘探和开发目标的复杂化,加之从勘探走向开发是当今地震技术的重要发展趋势之一,井中地震资料的高精度成像处理成为地震数据处理的必然要求,VSP地震资料以其精度高和直接测量目的层的特点,得到开发人员和油藏地球物理工作者的认可。针对VSP地震资料的成像方法,许多学者进行了深入研究。Wyatt[8]提出了用时间域正演计算合成VSP地震记录的方法,为后续的VSP地震资料成像研究提供数据基础;Lee等[9]提出了针对VSP地震数据的处理方法,并指出利用VSP地震数据可以有效进行子波整形、估算岩石波阻抗;Dillon[10]提出了对VSP数据进行Kirchhoff积分偏移计算,使得VSP地震资料的偏移归位更为准确,为资料的应用奠定基础。

前已述及,VSP资料对井周地层的精细刻画及深部构造高精度成像具有重要的意义，为进一步提高井周微幅构造的成像效果,需要研究基于VSP资料的逆时偏移方法。刘清林等[11]和朱金明等[12-13]提出了VSP资料的逆时偏移处理方法,但是包括上述作者在内的诸多学者对VSP资料逆时偏移的研究都仅限于理论与方法的研究。目前地震资料的逆时偏移软件多为国外产品。业界针对VSP地震资料处理的商业软件仅局限于速度分析等方面,尚无VSP资料的逆时偏移算法平台。本文在软件开发技术和VSP逆时偏移方法研究的基础上,在Linux操作系统环境下,应用Qt开发了图形用户界面,形成了一套基于Qt的Walkaway VSP逆时偏移处理软件,从而填补了VSP地震资料逆时偏移处理软件的空白,为VSP资料高精度成像方法研究和实际应用奠定了基础。

1 VSP逆时偏移原理

不同于传统的地面地震勘探,VSP观测系统是在地面激发,在井筒内的一定深度范围布置检波器,接收来自地下界面多次覆盖的地震数据,这在一定程度上降低了地面噪声的干扰,井中检波器记录的地震反射波无需传播至地面,旅行路径更短。因此,VSP地震资料具有更高的信噪比和更小的能量衰减。VSP数据偏移是资料有效应用的基础,目前常规VSP偏移处理技术,如VSP-CDP转换、射线法等,难以获得理想的成像效果,逆时偏移方法基于双程波方程,结合VSP数据资料特点,可对井旁复杂构造精确成像。

1.1 VSP逆时偏移的基本原理

根据VSP资料和逆时偏移算法的特点,在实际计算时首先进行震源波场的正演模拟,再在接收点位置载入检波点波场,并进行反传播模拟,然后将互相关成像条件应用于上述两种波场,进行同时刻成像,最后叠加生成偏移数据体。对Walkaway VSP地震资料而言,正、反传播波场的计算公式表述如下：

(1)

(2)

(3)

应用有限差分法对波场进行数值模拟,其时间二阶、空间2M阶精度的有限差分格式如下：

(4)

式中：i,j和m分别为波场在x,z和t方向上的离散序号;Δx,Δz和Δt分别为x,z和t方向的采样间隔;Cl为有限差分系数。本文应用12阶的高阶差分格式进行正演模拟和逆时偏移计算,同时为进一步提高波动方程数值解的精度,应用了优化差分系数的方法[14]。

1.2 随机边界条件

随机边界通过在速度模型边界设置随机速度层,从而降低边界反射的相干性,减小其对成像结果的影响。由于该方法仅在速度模型的外围区域设置随机速度层,并未对波动方程做出改变,因此不会造成波场信息的任何损失,对于震源波场正传来说,可以通过保存正传波场最后几个时刻(一般为两个)的波场信息,再进行反传达到节省存储空间的目的。其在边界部位的随机函数构造如下：

(5)

式中：v(x,z)为构造边界的随机速度;v1(x,z)为构造边界的原始速度;r为随机数;d为内边界与速度点之间的距离,其保证了速度的随机性随距离的增加而增加[15]。

逆时计算时,首先将震源波场正传至最大时刻,并存储最后两个时刻的波场信息,再利用这两个时刻的波场信息逆推之前各个时刻的波场,这一过程与检波点波场反传同时进行,即可同时得到每一时刻的震源波场和检波点波场,应用一定的成像条件(互相关成像条件或者震源归一化成像条件)进行成像,直至零时刻,即可得到逆时偏移结果。随机边界方法可以极大地节省存储空间,但计算量需增加一倍,即随机边界是一种以计算换存储的计算方案,其计算效率的问题可通过后续给出的GPU加速技术得到有效解决。

1.3 成像条件

逆时偏移中常用的成像条件有激发时刻成像条件、反褶积成像条件和互相关类成像条件等,由于互相关类成像条件实施简易且计算稳定而得到普遍采用,互相关成像条件可表示为：

(6)

式中：s(x,z,t)代表t时刻的震源波场;r(x,z,t)代表t时刻的检波点波场;s(x,z,t)与r(x,z,t)的乘积表示在t时刻对波场进行一次成像运算;求和计算代表最终所成的像是各时间步长成像的叠加。由此可知,互相关成像条件充分利用了双程波场进行成像,是地面和井中地震资料逆时偏移中广泛采用的成像条件。

针对VSP资料的特点,为保证目的层的成像效果,提高目的层的成像精度,本文采用的成像条件为震源归一化互相关成像条件,其计算公式如下：

(7)

式中：ε为稳定系数,保证成像过程计算稳定。

与(6)式所示的传统互相关成像条件相比,震源归一化互相关成像在成像数据的振幅上保持相对均衡,并且可在一定程度上补偿深层的能量,提高深部的成像效果。

1.4 Laplacian去噪

由逆时偏移相关成像条件的原理可知,在震源与检波点之间的射线路径上均存在满足相关条件的成像点,在对有效构造进行成像的同时,也会产生低频噪声,可采用Laplacian噪声压制方法压制该噪声,从而提高叠加数据体的信噪比[16-17]。Laplacian算子可表示为：

(8)

式中：kx,ky,kz为成像域波数在空间方向的分量;kI为成像域的波数矢量;ω为角频率;v为介质速度;θ为入射角。此处对Laplacian低频噪声压制方法的实现步骤不做重点阐述。

2 GPU/CPU并行加速

与地面地震资料的逆时偏移方法相同,为降低Walkaway VSP逆时偏移算法的计算成本,采用基于CUDA平台的GPU并行加速算法提高计算效率。因GPU拥有更多的计算核心,因此,多线程调用可实现大规模并行计算,提高Walkaway VSP逆时偏移方法的计算效率。

GPU和CPU协同并行计算中,将需要GPU加速计算的部分打包成函数形式,称之为kernel函数,可在GPU上直接运行加速算法。在调用kernel函数前,需划分计算量,并将其分配给不同的线程(Thread)计算,主要以Thread为基本单位。但是,CPU和GPU各具优势,因此,在编程过程中应根据其特性予以考虑,以充分利用两者优点达到节省计算资源和提高计算速度的目的。

当函数在GPU端计算时,线程(Thread)为基本单位,若干线程组成块(Block),若干块组成栅格(Grid),如图1所示。Block的大小具有一定限制,主要由显卡的型号决定,CUDA有对应变量对Block和Grid进行设置,并且每个线程都有编号与之对应。GPU上具有多种存储器,其结构如图2所示。图中仅显示常用的存储器,在调用kernel函数时,需先将数据传输至全局存储器内,每个Block对应一个共享存储器,共享存储器的大小固定,每个Block内的线程均可从共享存储器中读取数据,Block之间的共享存储器互不通信,其数据交互通过全局存储器进行,每个线程对应一个寄存器,并且只能从其对应寄存器内存取数据,但寄存器本身很小,在线程进行运算时,从共享存储器和寄存器进行数据读取较全局存储器更快,但由于这两种存储器本身存在限制,因此需充分考虑实际计算的需求。

图1 GPU编程模式

图2 GPU存储器结构

波场数值模拟是逆时偏移的核心,也是最耗时的部分,其本身包含大量的并行运算。对震源波场正传和检波点波场反传进行加速计算,可以明显提高逆时偏移的计算效率。实际编程时,首先利用CPU在串行计算方面的优势,从文件中读取相应的数据,包括震源子波、速度模型等参数,将其传输至全局存储器中,然后对计算任务进行优化分配,以充分利用计算资源,调用kernel函数进行波场正传的计算,并采用相应的存储策略。之后进行检波点处波场反传,在反传过程中,利用成像条件进行成像。最后,将所需要的数据从全局存储器中传回主机端进行保存。利用GPU并行计算,在很大程度上降低了逆时偏移算法的计算成本。

3 软件设计

3.1 Qt环境

Qt是基于C++的跨平台图形界面程序开发软件,具有跨系统平台移植性,可提供“一次编写,随处编译”的开发框架[18-19]。Qt通过信号和槽机制完成任意Qt对象之间的通信操作,整个过程由开发者自行定义,信号和槽的匹配关系不受数量限制,可提高对象处理的灵活性与多样性。创建软件界面主要有两种方式：一是通过Qt平台设计工具绘制图形;二是编写程序代码设计图形界面。本文采用后者,这样可为后续界面组件的编辑和调试提供方便,同时形成的软件系统也具有良好的可移植性。

3.2 架构设计

软件的架构设计与软件的执行效率密切相关,合理的架构设计可提高软件的执行效率,大幅缩短地震资料的处理周期。软件架构设计思路如图3所示。对逆时偏移算法模块进行预编译,生成可执行程序,软件界面向可执行程序发送参数和命令,反过来可执行程序给软件界面提供功能支持。在界面代码执行过程中,仅进行参数传递等操作,也就是说直接调用前期编译完成的算法可执行程序。这样做的优点主要体现在：形成的软件平台避免了界面代码与算法代码的混合编译,使编译更容易,计算效率更高,代码算法也更易于实现。另外需要指出的是,按照上述方法设计软件具有很好的可移植性和扩展性,只要提供可执行软件和必要的参数接口,就可以完成相关计算程序的界面化移植。

图3 软件执行路线

3.3 流程设计

软件平台的操作流程和技术流程对于软件的执行和应用效果至关重要,作业流程的设计将直接影响软件的执行效率。VSP逆时偏移软件平台的基本流程为工区建立、数据载入、参数设置、作业发送等4个关键步骤,根据软件系统不同子模块的结构和VSP逆时偏移的基本原理,设计了基于GPU并行计算的技术流程图(图4)。

图4 VSP地震数据逆时偏移处理技术流程

4 软件测试

在VSP逆时偏移算法的基础上,基于Qt环境开发了垂直地震剖面数据逆时偏移的软件平台,为测试算法的准确性和软件平台的实用性,设计如图5 所示的复杂断陷盆地模型。该模型包括了高陡构造、小断层和尖灭等复杂地质构造,模型在x方向长度为2700m,z方向深度为2800m。在利用该模型进行VSP逆时偏移软件测试前,首先进行VSP数据的正演模拟,在模型沿x方向0和2700m处共布置两口VSP井,正演观测系统参数如下：共270炮激发,炮间距10m;井中共280级检波器接收,检波器间距10m;空间网格采用边长10m的正方形;时间采样间隔1ms;震源子波为主频25Hz的Ricker子波;记录长度2s。图6为正演模拟的单炮记录。其中图6a为模型x=0处VSP井接收到的炮点坐标为0的单炮记录;图6b 为模型x=2700m处VSP井接收到的炮点坐标为1000m的单炮记录。分析可知,模拟资料的上、下行波波场特征清晰。由于正演模拟所用的模型较为复杂,导致炮记录上的时距曲线特征较为复杂,该模拟数据将用于测试VSP逆时偏移算法软件的成像效果。

将正演模拟得到的270炮记录作为数据输入,基于图5所示的速度模型,完成了逆时偏移试算,经过叠加和低频噪声压制后,得到的成像剖面如图7a 所示。为进行成像效果的对比分析,图7b给出了基于同样速度模型的地面地震数据的逆时偏移成像剖面,需要说明的是该数据同样是采取270炮激发,270个检波器布置于地面0～2700m处,检波点间距10m,除检波点位置不同外,其它的所有参数都和VSP资料的参数一致。对比VSP资料与地面地震资料的成像剖面可知,VSP逆时偏移算法对模型底部成像更清晰,主要原因是检波器在井下,可实现对目的层的直接照明和成像,从而使得深层介质及陡倾角地层的成像效果比地面地震更清晰。

图5 复杂断陷盆地模型

在成像精度测试的基础上,对软件平台所包括的CPU和GPU并行方式的Walkaway VSP逆时偏移模块进行了计算效率对比。利用前述模拟的270炮VSP地震记录进行了逆时偏移计算时间测试，测试结果如表1所示。基于CUDA平台GPU的Walkaway VSP逆时偏移方法可获得较高的计算效率,相比传统的CPU计算方式,可加速10倍左右。

图6 正演模拟的单炮记录

图7 VSP资料(a)与地面地震资料(b)逆时偏移成像效果对比

表1 CPU和GPU不同计算方式Walkaway VSP逆时偏移计算时间对比

计算方式单炮计算时间/s270炮计算时间/sCPU12634020GPU123240

5 结束语

逆时偏移是当前地震资料高精度成像的重要算法,Walkaway VSP地震资料的逆时偏移可实现复杂勘探和开发研究区构造的高精度成像。基于Qt开发工具集成Walkaway VSP逆时偏移算法,进而形成了软件系统平台,为实际Walkaway VSP资料的高精度成像计算提供了方便。结合CPU和GPU两种计算架构,实现了Walkaway VSP的逆时偏移计算,扩展了软件的应用范围,GPU并行加速算法,可在传统CPU计算的基础上,大幅提高计算效率,缩短成像计算的周期。复杂断陷盆地模型的成像试算结果表明,所形成的软件平台具有计算精度高、计算效率高和实用性强的优点。该软件开发的研究实践对同类软件的设计和开发具有一定的借鉴意义。

[1] Whitmore N D.Iterative depth migration by backward time propagation[J].Expanded Abstracts of 53rdAnnual Internat SEG Mtg,1983,382-385

[2] Baysal E,Kosloff D D,Sherwood J W C.Reverse time migration[J].Geophysics,1983,48(11):1514-1524

[3] McMechan G A.Migration by extrapolation of time-dependent boundary values[J].Geophysical Prospecting,1983,31(3):412-420

[4] 薛东川,王尚旭.波动方程有限元叠前逆时偏移[J].石油地球物理勘探,2008,43(1):17-21 Xue D C,Wang S X.Finite element wave equation prestack reverse time migration[J].Oil Geophysical Prospecting,2008,43(1):17-21

[5] 刘红伟,李博,刘洪.地震叠前逆时偏移高阶有限差分算法及GPU实现[J].地球物理学报,2010,53(7):1725-1733 Liu H W,Li B,Liu H.The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation[J].Chinese Journal of Geophysics(in Chinese),2010,53(7):1725-1733

[6] 郭念民,吴国忱.基于PML边界的变网格高阶有限差分声波方程逆时偏移[J].石油地球物理勘探,2012,47(2):256-265 Guo N M,Wu G C.High-order finite difference method in reverse-time migration with variable grids based on PML boundary condition[J].Oil Geophysical Prospecting,2012,47(2):256-265

[7] 刘红伟,刘洪,邹振,等.地震叠前逆时偏移中的去噪与存储[J].地球物理学报,2010,53(9):2171-2180 Liu H W,Liu H,Zou Z,et al.The problems of denoise storage in seismic reverse time migration[J].Chinese Journal of Geophysics(in Chinese),2010,53(9):2171-2180

[8] Wyatt K D.Synthetic vertical seismic profile[J].Geophysics,1981,46(6):880-891

[9] Lee M W,Balch A H.Computer processing of vertical seismic profile data[J].Geophysics,1983,48(3):272-287

[10] Dillon P B.Vertical seismic profile migration using the Kirchhoff integral[J].Geophysics,1988,53(6):786-799

[11] 刘清林,骆毅.VSP资料的叠前波动方程逆时偏移[J].石油地球物理勘探,1989,24(3):281-289 Liu Q L,Luo Y.Prestack reverse time migration of VSP data[J].Oil Geophysical Prospecting,1989,24(3):281-289

[12] 朱金明,严俊华.VSP逆时偏移处理[J].石油地球物理勘探,1991,26(5):564-570 Zhu J M,Yan J H.VSP reverse-time migration[J].Oil Geophysical Prospecting,1991,26(5):564-570

[13] 朱金明,董敏煜,李承楚.VSP的双程无反射波动方程逆时偏移[J].石油地球物理勘探,1989,24(3):256-270 Zhu J M,Dong M Y,Li C C.VSP reverse-time migration using two-way nonreflection wave equation[J].Oil Geophysical Prospecting,1989,24(3):256-270

[14] Zhang J H,Yao Z X.Optimized finite-difference operator for broadband seismic wave modeling[J].Geophysics,2012,77(1):A13-A18

[15] Clapp R G.Reverse time migration with random boundaries[J].Expanded Abstracts of 79thAnnual Internat SEG Mtg,2009,2809-2813

[16] 郭念民,冯雪梅,李海山.高阶拉普拉斯算子逆时偏移低频噪声去除方法[J].石油物探,2013,52(6):642-649 Guo N M,Feng X M,Li H S.Research on higher order Laplacian operator denoising method in reverse-time migration[J].Geophysical Prospecting for Petroleum,2013,52(6):642-649

[17] 刘诗竹,石颖,郭雪豹,等.VSP数据逆时偏移方法研究[J].地球物理学进展,2014,29(5):2211-2218 Liu S Z,Shi Y,Guo X B,et al.Investigation on reverse-time migration of VSP data[J].Progress in Geophysics(in Chinese),2014,29(5):2211-2218

[18] 郑阿奇.QT4开发实践[M].北京：电子工业出版社,2011:30-48 Zheng A Q.QT4 development practices[M].Beijing:Publishing House of Electronics Industry,2011:30-48

[19] 丁仁伟,李振春,仝兆岐.基于Qt的共聚焦点技术地震资料处理系统研发[J].石油物探,2009,48(2):199-205 Ding R W,Li Z C,Tong Z Q.Research and development of seismic data processing system based on CFP technology by utilizing Qt graphical user interface library[J].Geophysical Prospecting for Petroleum,2009,48(2):199-205

(编辑：陈 杰)

Qt based software development of reverse time migration for walkaway VSP seismic data

Wang Weihong1,Ke Xuan1,Guo Xuebao1,Zhang Yingying1,Yu Peng2

(1.CollegeofEarthSciences,NortheastPetroleumUniversity,Daqing163318,China;2.Oil&GasExplorationDepartment,NortheastOil&GasBranch,SINOPEC,Changchun130062,China)

The geophones in VSP geometry are located down in the borehole,therefore VSP seismic data has high S/N and with rich wavefield information.Based on two-way acoustic wave raypath,the prestack reverse time migration (RTM) has high precision for imaging,and provides important data for the study of complex subsurface structure near the borehole and fluid properties in targeted formation.According to the algorithm of reverse time migration for walkaway VSP seismic data,the complete RTM processing flow of VSP data was designed and the VSP RTM software was preliminarily developed,which is developed by Qt language as well as in Linux operating system.Two types of core algorithm systems,CPU and GPU can be installed in different hardware.The test of model data demonstrated that the developed software is characterized by easy operating,good transferring and high imaging precision.

VSP,reverse time migration,Qt,software interface

2015-01-23;改回日期：2015-05-21。

王维红(1975—),男,副教授,主要研究方向为地震资料数字处理。

国家自然科学基金面上项目(41474118)、国家高技术研究发展计划(863计划)项目(2012AA061202)联合资助。

P631

A

1000-1441(2015)04-0452-07

10.3969/j.issn.1000-1441.2015.04.012