2010, Vol. 53
2. 北京吉星吉达科技有限公司, 北京 100029;
3. 中石油大庆油田有限责任公司勘探开发研究院, 大庆 163712
2. Beijing Geostar Science and Technology Ltd. Beijing 100029, China;
3. Exploration and Development Institute, Daqing Oil Field Ltd, PetroChina, Daqing 163712, China
逆时偏移的思想早已由Hemon在1978年提出[1], 继而, 由Baysal、Whitemore、McMechan等[2~4]将之应用于纵波资料的叠后偏移, 他们在获得好效果的同时, 明确告诫尚存在艰难的数值计算效率和存储等关键问题待解决.当时的逆时偏移方法虽然只能限于在2D或者叠后偏移中应用, 却引发了业界的巨大兴趣.逆时偏移方法是通过双程波波动方程在时间域上对人工给予的震源子波正向传播和接收到的地震资料进行反向传播, 并结合成像条件实现偏移[5].一方面逆时偏移需同时计算震源在正向传播和地震资料的反向传播过程中的各个时刻的波场, 因其各自时间的延拓方向不同, 致使实际应用中必须存储其中一个方向的传播过程, 这就是何以需要巨大额外存储空间.如今的地震资料采集已经发展到了三维甚至四维, 其数据量更是十分庞大, 逆时偏移在推向实际应用时尚需面对此类困难.另一方面, 相较于单程波方程的波场延拓而言, 逆时偏移运用双程波波动方程进行波场延拓, 避免了上下行波的分离处理, 因而成为最准确的成像算法, 且不受倾角的限制, 并能实现回转波和多次波的成像.正基于此, 研究的焦点也集中于提高逆时偏移的计算效率和减少偏移中的存储I/O问题[6].近年来, 随着计算机技术的发展和计算能力的大幅度提高, 业界在逆时偏移研究工作中不断取得突破性的成果.例如各向异性的逆时偏移、弹性波逆时偏移以及逆时偏移的噪音衰减技术等[7, 8].2008年以来, 快速发展的GPU (Graphic Processing Unit)以及其编程架构CUDA (Computing Unified Device Architecture)已为深度偏移的快速计算提供更加有力的工具.本文着重由声波方程的高阶有限差分近似为切入点, 分析现行多种逆时偏移的实现策略, 提出一种利用随机边界的CPU/GPU为编程环境的地震逆时偏移算法实施的对策, 并经数值模型和实际资料验证其有效性.
2 波动方程的高阶有限差分近似三维声波方程的表达如下:
|
(1) |
其中t代表时间坐标, x, y, z代表空间坐标, P表示位移函数, V表示速度函数.利用时间二阶中心有限差分近似
|
(2) |
其中,
|
利用N阶中心差分近似得:
|
(3) |
其中,
|
分别用Δx, Δy, Δz表示差分网格的间距, h表示最大网格间距, n表示一个波形中的采样点数, Δt表示延拓步长, fmax表示子波频率的最大值, Vmin表示速度模型的最小值, 那么差分格式的频散条件可表示为
|
(4) |
稳定性条件:
|
(5) |
成像条件的制定是地震偏移成像算法的关键之一, 它直接影响成像的效果和计算成本, 本文应用Claerbout的互相关成像原理.
互相关成像条件的成像公式为
|
(6) |
其中, U(x, z, t)为上行波波场, D(x, z, t)为下行波波场, dt为延拓时间步长.双程波逆时偏移中以外推观测波场PR(x, z, t)取代上行波, 以震源外推波场PS(x, z, t)取代下行波, (6)式变为
|
(7) |
(7)式中的被积函数PR(x, z, t)PS(x, z, t)表示t时刻对波场做一次成像运算, 积分说明如空间I(x, z)中的像是多个时间步成像的叠加.所以, 互相关成像条件充分利用了成像信息, 在增强成像信号的同时也有效压制了成像噪音.逆时偏移互相关成像条件的成像步骤是:
(1) 首先要完成一次波动方程正演计算, 并保存每一时间步的波场信息;
(2) 逆时外推记录波场, 存储每一时间步的波场信息;
(3) 分别读取保存的同一时刻的震源波场和记录波场做成像运算, 即PR(x, z, t)PS(x, z, t), 再累加入成像空间.
4 算法实现策略分析逆时偏移的成像条件需要使用在同一时刻的震源波场(经过激发和正演模拟波的传播)和记录波场(逆时反传播回地下的波场).因为一个是正传波场, 另一个是反传波场, 如果想同时得到相同时刻的两个波场, 则必须存储其中一个波场的传播过程, 即每一时刻的波场分布, 需要消耗甚大的存储资源, 这个要求在实际操作中是难以满足的.为此, 许多学者找到一些解决办法.其中包括:
策略一 震源和记录同时从Tn(表示n*dt时刻)延拓到T0, 这样就不需要额外的存储空间, 并且可以达到最小计算量O(2N)的量级, 然而Tn时刻的震源波场并不是已知数据, 所以不能直接实现.
策略二 存储T0到Tn的波场, 然后反传记录波场, 当传播的某个时刻时读取那个时刻的激发点波场.此方法也可以达到O(2N)计算量的要求, 但是需要巨大的磁盘空间来存储波场, 即使磁盘空间可以满足需求, 但是因此产生的I/O时间消耗十分庞大.
策略三 只存储少数几个时刻的震源波场, 在记录波场反向延拓的过程中, 利用这些存储的震源波场插值来近似当前时刻的波场然后成像.这种方法对波场的插值首先是不准确的, 并且计算量也随之增加.
策略四 首先把震源波场正向延拓到Tn和Tn-1两个时刻, 把这两个波场作为初始条件, 与记录波场同时反向延拓, 并且随延拓过程中应用成像条件, 不需要额外的存储空间.这种策略只可以应用在密度为常数的介质中, 并且激发点边界条件是Dirichlet边界条件.但是这种波场的延拓假设是在无限大空间内传播的, 而通常的计算方法是加入了人工边界的, 包括自由边界条件、吸收边界条件、指数衰减边界条件和完美匹配层边界条件(PML)等[9~11].其中除自由边界条件外, 其他三种边界条件都会破坏波场的完整性, 使得波场的延拓不可逆.为此只能在适当的时刻在适当的空间位置上插入波场能量以补充波场的完整性[12, 13], 这样就需要存储每个时刻的位于人工边界处的波场值, 当计算规模增大后, 存储量也会快速上升, 尽管小于策略二中需要的存储量, 但是也是通常不能满足实际的要求.
策略五 首先正向延拓波场, 并且记录多组检查点(Pn, Pn-1)作为中间时刻的波场初始条件, 然后反向延拓检波点波场到某时刻, 用最近的检查点来正向传播得到该时刻的震源波场, 然后应用成像条件, 当应用PML边界条件的时候波场还需要额外的存储空间, 并且需要大量的重复计算过程, 最大的重复计算可能会达到O(mN)的量级, 其中m为检查点的个数.Griewank利用这个思路提出优化的检查点的选取方法, 可以使得重复计算的次数减少到O(Nlogm)的量级[14].尽管如此, 逆时偏移的消耗时间和存储的矛盾也并未得到完全解决, 但是该方法是目前最经济的计算策略.
在上述5种逆时偏移的计算策略中可以看出计算量级与存储规模是一对矛盾, 若要得到最小的计算量就会达到最大的存储量, 反之亦然.然而在策略四中不需要额外的存储空间, 而其他策略都要涉及到存储规模与计算量级的平衡问题.在策略四中, 人工边界是主要的影响因素, 然而自由边界条件又会引入很多的成像假象, 其他的边界条件又不符合该策略的要求只能舍弃.针对策略四, 本文依据2009年Robert提出的随机速度边界模型提出一种利用GPU优化的逆时偏移计算方法.
5 随机边界的建立2009年Robert[15]提出随机边界模型, 其思想是消除人工边界自由边界条件反射波的相干性, 使边界反射不能成像.
具体实现过程为将有限空间外扩一定的距离, 然后在外扩的空间内填充随机速度, 从而形成随机边界速度模型(如图 1b), 当波传到随机速度区域时波前面将被随机化, 使波场变成随机噪音反传回真实速度区域, 破坏了边界反射的相干性使得边界反射不能成像.
|
图 1 原始速度模型(a)及随机速度边界模型(b) Fig. 1 Original velocity model (a) and random velocity boundary model (b) |
本文构造随机边界函数如下:
|
(8) |
其中Ṽ(x, z)为边界点的随机速度函数, V(x, z)为边界点的原始速度函数, r是随机数, d为速度点与内层边界的空间距离.
运用CPU/GPU作为偏移计算核心, 采用随机速度边界模型, 将激发点波场传播至Tn时刻, 然后利用Tn时刻的波场作为初始条件, 同时反传激发点波场和检波点波场, 并且同时利用成像条件进行成像.进而避免了额外的存储空间, 虽然激发点的波场重复计算一次, 但利用GPU作为波场延拓的数值计算核心其耗时与大规模的磁盘I/O相比还是非常经济的.另外, 从文后的实际测试中也可表明, 利用了GPU作为波场延拓函数的计算核心, 延拓过程的计算已经不是主要瓶颈.
存储检查点的计算策略五是目前认为最经济的计算策略, 所以利用该策略与本文的计算策略做对比.
6 数据试算为了验证上述方法的有效性, 我们设计了一个空间采样为256×c的网格(如图 1所示), 取边界宽度NL=30, 网格间距为10 m, 在网格的三个边界生成随机速度模型, 在模型的中间放置震源, 然后记录不同时刻的波前快照.
计算结果如图 2所示, 波前面正传到120×4 ms时的波场尚未到达模型的边界(图 2a), 正传到180×4 ms时的波场已经到达随机边界并且有部分反射能量(图 2b), 但是与自由边界条件(图 2e)相比在同一时刻的反射能量的相干性已经削弱很多.当正传到400×4 ms时的波场已经完全成为随机噪音.我们将400×4 ms时的波场作为初始值再反传到120×4 ms时刻, 其波场示于图 2d, 利用正传波场与反传波场在120×4 ms时刻的差示于图 2f, 其误差值量级在10-7以内, 与单精度浮点运算的精度范围相同, 所以我们认为这样的误差是可以接受的.
|
图 2 随机边界模型的波场正演快照 (a) 正传到120×4ms时的波场; (b) 正传到180×4ms时的波场; (c) 正传到400×4ms时的波场; (d) 反传到120×4ms时的波场; (e) 应用自由边界条件时正传到180×4ms时的波场; (f) 正传波场(a) 与反传波场(d) 在同一时刻的差. Fig. 2 random boundary model forward propagation wave field snapshot (a) Forward wave field at 120×4 ms; (b) Forward wave field at 180×4 ms; (c) Forward wave field at 400×4 ms; (d) Backward wave field at 120×4 ms; (e) Forward wave field at 120×4 ms with free reflection boundary condition; (f) Difference between (a) and (d). |
利用常用的BP2004模型作为实例, 运用本文提供的方法做数值试验, 分别对该模型利用自由边界条件(FBC)、吸收边界条件(ABC)和本文的随机边界条件(RBC)做偏移测试, 其中自由边界采用成像策略四, 吸收边界采用成像策略五, 随机边界采用本文的成像策略.使用吸收边界条件时的时间延拓步长与其他两种方法相同, 但成像时间步长数据的采样率为4 ms, 这样使得内存的使用量可以达到最小.为了不失公平性, 偏移其他参数完全相同.偏移结果示于图 3, 计算所需时间和存储需求示于表 1.
|
图 3 三种边界条件的偏移结果 (a) BP2004速度模型; (b) 自由边界偏移结果; (c) 随机边界偏移结果; (d) 吸收边界偏移结果. Fig. 3 Three boundary condition migration result (a) BP2004 velocity model; (b) FBC migration result; (c) RBC migration result; (d) ABC migration result. |
|
表 1 计算用时和存储量 Table 1 Computing time and memory requirements |
从图 3中的计算结果来看, 随机边界的成像质量大大优于自由边界条件的成像效果, 其原因是显而易见的.随机边界的成像质量略优于吸收边界, 笔者认为, 本文采用的吸收边界为1977年Clayton & Engquist提出的单程波方程的二阶旁轴近似, 其吸收效果并不是十分理想, 仍然有边界反射能量反传回来而造成成像结果较差, 若采用完美匹配边界条件(PML)可能效果会有所改善.随机边界条件使得反射波成为随机噪音反传回来, 经过多炮的随机叠加有消除噪音的作用.
从表 1中可见, 若利用CPU作为计算核心, 完成逆时偏移的耗时非常巨大, 而利用CPU/GPU联合作为计算核心, 计算耗时已经大幅度降低, 这样更加方便应用于实际当中.若利用CPU/GPU联合作为计算核心, 不同边界条件的耗时基本相当, 但是吸收边界仍然需要大量的存储空间, 而随机边界条件和自由边界条件都只需要很少的内存空间即可完成, 并且而随机边界用时要少于吸收边界的计算时间, 略大于自由边界条件的计算用时.笔者认为, 随机边界尽管要重复计算一次震源波场的逆向传播, 但是相比于吸收边界的存储IO时间来说还是值得的.
8 结语本文利用CPU/GPU作为计算核心, 并将计算量最大的波场延拓放入GPU中计算.在基本保证成像质量的前提下, 凭借GPU的计算优势, 大大减少了计算周期, 并且利用随机速度边界的思想解决了大规模的存储I/O问题.不过在模型很大的情况下, 还有可能遇到GPU显存不足的问题, 还有待继续研究.
致谢感谢北京吉星吉达科技有限公司允许发表本文; 感谢BP公司提供模型数据.
| [1] | Hemon C. Equations d'onde et modeles. Geophysical Prospecting , 1978, 26: 790-821. DOI:10.1111/gpr.1978.26.issue-4 |
| [2] | Whitmore D W. Iterative depth migration by backward time propagation. 53rd Annual International Meeting, SEG, Expanded Abstracts, 1983. 382~385 http://www.oalib.com/references/18988173 |
| [3] | Baysal E, Kosloff D D, Sherwood J W C. Reverse time migration. Geophysics , 1983, 48: 1514-1524. DOI:10.1190/1.1441434 |
| [4] | McMechan G A. Migration by extrapolation of time-dependent boundary values. Geophysical Prospecting , 1983, 31: 413-420. DOI:10.1111/gpr.1983.31.issue-3 |
| [5] | Claerbout J. Toward a unified theory of reflector mapping. Geophysics , 1971, 36: 467-481. DOI:10.1190/1.1440185 |
| [6] | Yoon K, Marfurt K J, Starr W. Challenges in reverse-time migration:74th Annual International Meeting, SEG, Expanded Abstracts, 2004.1057~1060 http://www.oalib.com/references/18986628 |
| [7] | Yu Zhang, James Sun, Samuel Gray. Reverse-time migration:Amplitude and implementation issues. 78th Annual International Meeting, SEG, Expanded Abstracts, 2009.2145~2149 |
| [8] | Houzhu (James) Zhang, Yu Zhang. Reverse time migration in 3D heterogeneous TTI media. 78th Annual International Meeting, SEG, Expanded Abstracts, 2009. 2196~2200 |
| [9] | Robert W. Clayton and Engquist, Absorbing boundary conditions for wave equation migration. Geophysics , 1980, 45: 95-101. |
| [10] | Virieux J. P-sv wave propagation in heterogeneous media:Velocity stress finite difference method. Geophysics , 1986, 51: 889-901. DOI:10.1190/1.1442147 |
| [11] | Marfurt K J. Accuracy of finite difference and finite element modeling of the scalar and elastic wave equations. Geophysics , 1984, 49: 533-549. DOI:10.1190/1.1441689 |
| [12] | Symes W W. Mathematical Foundations of Reflection Seismology. Technical Report, Rice University. 1995 http://www.oalib.com/references/19002042 |
| [13] | Symes W W. Reverse time migration with optimal check-pointing. Geophysics , 2007, 72: 213-221. DOI:10.1190/1.2742686 |
| [14] | Griewank A. Achieving logarithmic growth of temporal and spatial complexity in reverse automatic differentiation. Optimization Methods and Software , 1992, 1: 35-54. DOI:10.1080/10556789208805505 |
| [15] | Robert G. Clapp, Reverse time migration with random boundaries. 79th Annual International Meeting, SEG, Expanded Abstracts, 2009.2809~2813 http://www.oalib.com/references/18988276 |