镜像法表面多次波正演模拟研究

巩长春,刘　韬,李　琳,季玉新

(1.中国石油化工集团公司石油工程地球物理有限公司,北京100020;2.中国石油化工股份有限公司石油勘探开发研究院,北京100083;3.中国石油化工股份有限公司多波地震技术重点实验室,北京100083;4.中国地震局地球物理研究所,北京100081)



镜像法表面多次波正演模拟研究

巩长春1,刘韬2,3,李琳4,季玉新2,3

(1.中国石油化工集团公司石油工程地球物理有限公司,北京100020;2.中国石油化工股份有限公司石油勘探开发研究院,北京100083;3.中国石油化工股份有限公司多波地震技术重点实验室,北京100083;4.中国地震局地球物理研究所,北京100081)

针对海洋地震勘探的海底电缆(OBC)采集系统,提出了一种直接正演表面多次波的方法,其核心思想是采用镜像延拓的思路修改模型,将海底对称映射到海平面之上作为新边界。将原始模型的表面多次波传播路径等效成新模型下一次反射波的传播路径,因此,只需将新模型的表面边界设置成吸收边界就可以方便地计算出多次波的传播特征。将该方法应用于多个地质模型数据,正演出的多次波波场和实际多次波信息较为吻合,证明了方法的有效性和便捷性。

多次波;正演;海底电缆;镜像延拓;起伏边界

伴随着一次反射信号的多次波被地震记录仪接收。在海上地震勘探中,上行波场经海水面反射产生的表面多次波能量非常强,而且往往和有效波叠合在一起,严重影响地震成像精度和解释结果的可靠性。人们为此提出了多种去除多次波的方法,包括几何滤波类方法,如F-K方法、聚束滤波类方法[1-5];基于波动方程类的方法,如自由表面相关的多次波衰减方法(Surface related multiple elimination,SRME),逆散射级数方法和逆数据域方法等[6-9]。然而随着技术的发展,多次波自身携带的一些地质元素逐渐被挖掘出来,在特定情况下它甚至可以提供常规有效信号所缺乏的信息,作为一次反射波场的补充。例如近些年的多次波成像技术[10-12],就是利用了多次波的这种特性去改善成像质量和扩大照明范围。

因为多次波传播路径复杂,其波场特性相对于一次反射波更加复杂。为了更好地研究多次波的传播特性,开展多次波的正演模拟研究尤为重要。正演模拟中波动方程类方法在工业界应用最广,但是常规的波场正演涵括了所有的波场,包括一次波场、多次波场、转换波场等,很难自动地分离出多次波场。因此,有必要寻求一种直接正演多次波场的方法。在这方面前人已经开展过一些研究：HRON[13],顾汉明等[14]采用求取多次波射线路径来分析多次波的振幅特性;姚姚等[15]引入非线性褶积模型,实现了一维情况下层间多次波的正演;朱振宇等[16]基于波前重构的思路,利用修改的程函方程正演出了鬼波和层间多次波。以上方法基本上都是基于高频射线来进行多次波正演模拟,实现起来比较复杂。我们针对海上勘探的海底电缆(OBC)观测系统,利用多次波的几何学特征,采用

镜像延拓的方式将多次波场巧妙地转换成一次反射波场,从而直接用波动方程进行正演,可以方便地得到多次波场的传播特征。在介绍其基本原理之后,我们将该方法应用于对几个地质模型进行正演实验,并针对不同的地质模型介绍了镜像延拓法的实施步骤以及正演方法的选择,比较了计算的多次波结果和正演的全波场数据。结果显示该方法能够直接有效地模拟出多次波场信息,并且适用于任意复杂介质模型。

1　基本原理

图1和图2分别为海底水平情况下和起伏情况下镜像法计算多次波的示意图。图1中假设震源置于海面A处,检波器放置于海底D处,震源激发信号,穿过海水入射到地层反射点B处,经反射回来之后到达海面C处,反射形成多次波被D处检波器接收。这是OBC数据中常见的多次波传播路径,这种多次波可以通过镜像的方法被等效成一次反射波的传播路径：将放置于海底的检波器相对于海面进行镜像对称,得到镜像层如图1b中虚线所示,镜像层与海面的速度设置成海水的速度即可。假设海面是水平的,根据几何地震学原理,从海底之下反射回来的波场经海面C点反射之后到D点接收的这段路径可以等效于C点激发、E点接收的路径。其中E点和D点是关于海面镜像对称的。也就是说,在使用了镜像层之后,原先的

图1　海底水平情况下镜像法计算多次波示意图解a 实际多次波路径示意; b 镜像多次波路径示意

图2　海底起伏情况下镜像法计算多次波示意图解a 实际多次波路径示意; b 镜像多次波路径示意

A-B-C-D的多次波反射路径可以等效成A-B-C-E的一次波反射路径。因此,我们只需设置一个镜像层,并在计算波场的时候将镜像层的边界设置为吸收边界,就可以简单地通过一次波的求解方式求解多次波的信息。对于海底起伏的情况,同样可以采用镜像延拓的方式,将海面多次波的传播路径映射成一次波的传播路径,然后再通过波动方程进行求解。唯一的区别在于：对于起伏的海底地形,镜像对称之后的边界是起伏的界面,相当于需要正演一个起伏地表的问题(图2)。

2　正演方法

我们选取基于交错网格的有限差分方法来进行正演计算。二维声波方程的速度-应力表达式满足：

(1)

式中：p代表声压;K代表弹性参数;ρ代表密度;vx和vz分别代表速度在x方向和z方向的分量。为了避免计算中带来的奇偶失联问题,速度-应力方程一般采用交错网格进行计算,其物理量的格点配置情况如图3所示。具体计算时采用的差分格式为：

(2)

(3)

式中：cm代表差分系数;l代表网格步长。

图3　交错网格的物理量配置方式

对于正演计算,除去核心区域的算子外,另一个非常重要的因素则是边界问题。设置不同的边界条件可以得到不同的波场信息。对于声波方程,常用的边界有两种情况：吸收边界和刚性边界。其中吸收边界应用最为广泛,可以吸收掉有效计算区域带来的人工反射问题,从而获得一次反射波的信息;而刚性边界也称为反射边界,例如在模拟海面多次波的情况下,直接将控制方程中的波场函数设置为0,即可在海面产生自由表面多次波。刚性边界条件是自由地表边界条件的一种特殊情况,其在计算的时候比较简单,直接置0即可。

本文应用的是PML吸收边界,也是近些年被工业界认为吸收效果较好的一种边界条件,具体的原理如下所示[17-18]。

考虑二维弹性波场的速度分量：

(4)

在x方向上采用吸收边界,可以等效于求解：

(5)

式中：d(x)代表衰减因子,当d(x)=0时即为正常的波动方程形式。同样地可以将PML边界条件下的方程扩展到其它变量的求解。

3　镜像法多次波正演

我们引入几个海底OBC采集的正演模型,采用镜像法对模型进行海面多次波的正演模拟,并将模拟结果与常规正演模拟的多次波进行对比。

3.1平层模型

引入如图4所示的平层模型。该模型深度为2km,分为4层介质：第1层为水层,速度设置为1500m/s;第2层速度设为2500m/s;第3层速度设为3500m/s;第4层速度设为4500m/s。震源放置于海面,检波器放置于海床。根据镜像法的原理,将检波器根据海面对称映射到海面之上,则可以进行海面多次波的正演模拟。模拟时,震源使用主频为25Hz的Ricker子波,60道单边接收,道间距40m,记录时长2500ms,4ms采样。

图4　平层模型a 原始OBC地震资料采集模型; b 应用镜像对称之后的模型

正演计算得到的波场信息如图5所示。其中图5a为将海面设置成刚性边界条件时计算出的波场,在单炮记录中出现了直达波、一次反射波以及较强的海面多次反射波场;图5b为将海面设置成吸收边界条件时计算出的波场,图中显示了较为明显的直达波和一次反射波;图5c是根据镜像法延拓之后计算得到的单炮记录,图中显示了直达波以及计算出的海面多次波信息。将图5c与图5a,图5b 对比,可以看到,图5c计算出的多次波场和实际产生的多次波场相吻合。

3.2海底为平层的复杂速度模型

我们引入海底为平层的复杂速度模型,如图6a所示。该模型横向10km,纵向2km,海床被设置

为一个平层,海底之下存在一个起伏地层。我们采用有限差分方法进行正演模拟,其中震源放置于海面,检波器放置于海床,60个检波器记录,道间距40m,记录时间2500ms,采样间隔4ms,震源采用主频为25Hz的Ricker子波。

同样的,采用镜像延拓方法,将检波器从海床对称延拓到海面以上,如图6b所示,进行正演模拟。图7和图8分别给出了震源在不同位置时采用不同边界条件以及镜像法计算得到的单炮记录。其中,图7的震源放置于横向3km处。图7a为将海面设置为刚性边界条件下的正演结果,其中,直达波、一次反射波、表面多次波以及起伏边界产生的绕射波都清晰可见;图7b为将海面设置成吸收边界条件的正演波场,其中表面多次波的信息已经消除,只剩下直达波场、一次波场;图7c 是采用镜像法正演得到的单炮记录。将图7c与图7a,图7b对比可以看出,镜像法多次波正演在复杂介质中也取得了良好的效果。图8的震源放置于横向1km处,可以看到,因为海底构造的变化导致了单炮记录上反射同相轴和多次波同相轴的变化。同样的,对比图8a和图8c可以发现,镜像法在复杂介质的多次波正演中依然方便、有效。

图5　平层模型单炮模拟记录a 将海面设为刚性边界条件时的单炮记录; b 将海面设为吸收边界条件时的单炮记录; c 镜像法计算出的单炮记录

图6　海底为平层的复杂速度模型a 原始OBC地震资料采集模型; b 应用镜像对称之后的模型

图7　海底为平层的复杂速度模型震源在3km处的单炮记录a 将海面设为刚性边界条件时的单炮记录; b 将海面设为吸收边界条件时的单炮记录; c 采用镜像法计算出的单炮记录

3.3海底起伏的速度模型

引入海底起伏的速度模型,如图9a所示。该模型横向10km,纵向2km,检波器放置于起伏的海底之上。海底起伏的速度模型比平层模型复杂,但是基本原理一致：将起伏的海底相对于海面进行镜像对称,反转到海面之上,如图9b所示。在镜像对称之后,将检波器的位置同样镜像对称到海面之上,再进行正演计算即可得到海底起伏情况下的表面多次波。

需要注意的是,在镜像延拓后转换成正演计算一个起伏地表模型。有限差分方法在使用矩形网格刻画边界时容易引起强烈绕射波场,从而影响整体计算效果;同时,起伏边界的吸收条件非常复杂,实现起来也很不方便。我们将起伏地表之外的部分边界填充成为一个规则边界,如图9b所示。将填充区域的速度设置成和海水速度一致,从而可以减少填充区域和海水的速度差异导致的强散射问题。扩充后的表面是平的,因此也可以使用PML吸收边界条件计算。在最终正演计算的时候,仍然选取60个检波器进行单边接收,其中道间距为40m,记录时间为2500ms,采样间隔为4ms,震源是主频为25Hz的Ricker子波。

图8　海底为平层的复杂速度模型震源在1km处的单炮记录a 将海面设为刚性边界条件时的单炮记录; b 将海面设为吸收边界条件时的单炮记录; c 采用镜像法计算出的单炮记录

图9　海底起伏的速度模型a 原始OBC地震资料采集模型; b 应用镜像对称之后的模型

正演计算的单炮记录如图10和图11所示,其中,图10中的单炮记录震源位置在1km处,图11中单炮记录的震源位置在3km处。图10a为将海面设置成刚性边界条件的单炮记录,可以看到比较明显的直达波、一次波和多次波信息。因为海底起伏不平,因此计算出的反射同相轴不是双曲线形状,发生了一定的扭曲,同样的情况也体现在海面多次波的同相轴中。图10b为将海面设为吸收边界条件计算的单炮记录,该单炮记录上无表面多次波,可以看到直达波和3条一次反射波的同相轴。图10c为采用镜像方法计算出的单炮记录。对比图10a和图10c可以发现,直接正演计算出的多次波走时信息和全波场计算出的多次波是吻合的。图11是震源位于横向3km处的计算结果。类似地,反射同相轴与图10产生了差异,而直接正演出的多次波场和全波场正演中的多次波走时信息也是一致的。计算结果证明了起伏地表正演的有效性以及镜像法多次波正演的准确性。

图10　海底起伏速度模型震源在1km处正演的单炮记录a 将海面设为刚性边界条件时计算的单炮记录; b 将海面设为吸收边界条件时计算的单炮记录; c 镜像法计算出的单炮记录

图11　海底起伏速度模型震源在3km处正演的单炮记录a 将海面设为刚性边界条件时计算的单炮记录; b 将海面设为吸收边界条件时计算的单炮记录; c 镜像法计算出的单炮记录

4　结束语

本文采用镜像延拓方法对海底OBC地震采集模型中的表面多次波进行了正演研究。相对于其它多次波正演方法,本文方法充分利用了地震波传播的几何特性,将多次波的传播路径转换成一次波的传播路径,从而可以方便地使用常规的正演算子对海底电缆(OBC)采集模型进行多次波正演模拟。在镜像延拓后模拟起伏地表模型中,本文采取的填充方法,将起伏边界规则化,并填充相应的海水速度,从而能够有效规避矩形网格在刻画边界过程中产生的强绕射波场,使得正演结果更加精确可靠。文中引入的几个速度模型的表面多次波正演结果和常规全波场正演结果的对比显示,镜像延拓方法能够方便、有效地模拟表面多次波场的传播特征。

[1]CLAERBOUT J F.Fundamentals of geophysical data processing [M].London:Blackwell Scientific Publications,1976:1-100

[2]朱宪怀.应用预测反滤波迭代法消除多次波[J].石油物探,1983,22(1)：54-61

ZHU X H.Elimination of multiples by predicative inverse filter iteration[J].Geophysical Prospecting for Petroleum,1983,22(1):54-61

[3]李凌云.应用Radon变换消除多次波[J].石油天然气学报,2006,28(4)：264-266

LI L Y.Elimination of multiples with Radon transform method[J].Journal of Oil and Gas Technology,2006,28(4):264-266

[4]王维红,崔宝文.双曲Radon变换法多次波衰减[J].新疆石油地质,2007,28(3)：363-365

WANG W H,CUI B W.Multiple attenuation by hyperbolic Radon transform[J].Xinjiang Petroleum Geology,2007,28(3):363-365

[5]刘喜武,刘洪,李幼铭.高分辨率Radon变换方法及其在地震信号处理中的应用[J].地球物理学进展,2004,19(1)：8-15

LIU X W,LIU H,LI Y M.High resolution radon transform and its application in seismic signal processing[J].Progress in Geophysics,2004,19(1):8-15

[6]SEN M K,LIU F,STOFFA P L,et al.A unified treatment to free surface multiple elimination methods[J].Journal of Seismological Exploration,1998,7:129-143

[7]田继强,胡天跃.反馈迭代法在自由表面多次波压制中的应用[J].石油物探,2008,47(5)：449-454

TIAN J Q,HU T Y.Application of feedback iteration method in free surface multiple attenuation [J].Geophysical Prospecting for Petroleum,2008,47(5):449-454

[8]金德刚,常旭,刘伊克,等.逆散射级数法预测层间多次波的算法改进及其策略[J].地球物理学报,2008,51(4):1209-1217

JIN D G,CHANG X,LIU Y K,et al.Algorithm improvement and strategy of internal multiples prediction based on inverse scattering series method[J].Chinese Journal of Geophysics,2008,51(4):1209-1217

[9]马继涛,陈小宏,黄小宁.反数据域压制多次波方法研究[J].石油地球物理勘探,2009,44(5):537-542

MA J T,CHEN X H,HUANG X N.Studying on multiple elimination in inverse data domain[J].Oil

Geophysical Prospecting,2009,44(5):537-542

[10]SCHUSTER G T,YU J,SHENG J,et al.Interferometric daylight seismic imaging[J].Geophysical Journal International,2004,157(2):838-852

[11]SNIEDER R,WAPENAAR K,LARNER K.Spurious multiples in seismic interferometry of primaries [J].Geophysics,2006,71(4):SI111-SI124

[12]吴世萍,黄录忠,胡天跃.Walkaway VSP多次波成像技术研究[J].石油物探,2011,50(2):115-123

WU S P,HUANG L Z,HU T Y.Multiple reflection imaging by using walkaway VSP data[J].Geophysical Prospecting for Petroleum,2011,50(2):115-123

[13]HRON F.Efficient computation of the amplitude distance curves for multiples and converted phase in horizontally layered media[J].Expanded Abstracts of 68thAnnual Internat SEG Mtg,1998:1727-1729

[14]顾汉明,王纬,陈国俊.复杂介质中地震多次反射波快速正演模拟[J].地球科学,2001,26(5):541-544

GU H M,WANG W,CHEN G J.Fast synthetic seismograms for multiples in any complex stratified media[J].Earth Science,2001,26(5):541-544

[15]姚姚,蔡其新,张果,等.微屈多次波正演模拟及其在中原油田深层地震资料解释中的应用[J].物探与化探,2004,28(5):443-446

YAO Y,CAI Q X,ZHANG G,et al.The forward simulation of peg-leg multiple wave and its application to the deep seismic data interpretation in the zhongyuan oilfield[J].Geophysical and Geochemical Exploration,2004,28(5):443-446

[16]朱振宇,赵伟.基于地质建模的波前重构正演方法模拟多次波[J].天然气工业,2007,27(增刊A): 155-156

ZHU Z Y,ZHAO W.Modeling of multiples with wave reconstruction based on geological model building[J].Natural Gas Industry,2007,27(S1):155-156

[17]COLLINO F,TSOQKA C.Application of the perfectly matched absorbing layer model to the linear elastic dynamic problem in anisotropic heterogeneous media[J].Geophysics,2001,66(1): 294-307

[18]单启铜,乐友善.PML边界条件下二维粘弹性介质波场模拟[J].石油物探,2007,46(2):126-130

SHAN Q T,LE Y S.Wavefield simulation of 2-D viscoelastic medium in perfectly matched layer boundary[J].Geophysical Prospecting for Petroleum,2007,46(2):126-130

(编辑：陈杰)

Free surface multiples modeling with mirror image method

GONG Changchun1,LIU Tao2,3,LI Lin4,JI Yuxin2,3

(1.SinopecGeophysicalCorporation,Beijing100020,China;2.SinopecPetroleumExplorationandProductionResearchInstitute,Beijing100083,China;3.SinopecKeyLaboratoryofMultiComponentsSeismicTechnology,Beijing100083,China;4.InstituteofGeophysics,ChinaEarthquakeAdministration,Beijing100081,China)

Multiples always exist with the primary seismic reflection for offshore seismic acquisition.Most of the time,the multiples are regarded as noises,which may destroy the seismic imaging results and seismic interpretation.Thence,the removal of seismic multiples is very important to seismic exploration,and seismic modeling helps to know more about the performance of multiples.In this paper,a new method to directly simulate free surface multiples for ocean bottom cable (OBC) data acquisition is proposed.For the OBC data acquisition,the original geophones are placed at the position symmetrically to the sea surface.In that way,the traveling path of the seismic data recorded by those virtual geophones equal to free surface multiples in original acquisition geometry.Thus the free surface multiples can be computed after modifying the velocity model by setting a new boundary,which is the mirror image of the seabed to the sea surface.The free surface multiples modeling with mirror image method was applied on several OBC data acquisition system,and the modeling results indicate that it can simulate the free surface multiples correctly and efficiently.

multiples,seismic modeling,ocean bottom cable (OBC),mirror extension,irregular boundary

2015-09-20;改回日期：2016-03-05。

巩长春(1983—),女,工程师,研究方向为可控震源高效采集方法及弹性波地震勘探等。

国家高技术研究发展计划(863计划)(2013AA064201)、国家重点基础研究发展计划(973计划)(2012CB214802)和国家科技重大专项(2011ZX05005-005)共同资助。

P631

A

1000-1441(2016)04-0475-08DOI:10.3969/j.issn.1000-1441.2016.04.002

This research is financially supported by the National High-tech R&D Program (863 Program) (Grant No.2013AA064201), the National Key Basic Research and Development Program of China (973 Program) (Grant No.2012CB214802) and the National Science and Technology Major Project of China (Grant No.2011ZX05005-005).