CN1997914B - 三维反虚反射 - Google Patents

Abstract

本发明提出的方法在对海上地震数据进行反虚反射处理时考虑三维效应。该方法依赖于能得到横向测线方向的二阶空间导数。该二阶横向测线导数能被直接估计出来，或通过间接测量其他波场量并使用波动方程技术计算所希望的项。优选地，该方法利用在水面附近拖曳的多分量拖缆、在海面附近的双拖缆配置或者沿垂直或水平分开的三拖缆配置来估计压力的二阶垂直导数。

Description

三维反虚反射

技术领域

本发明涉及对获取的海上地震数据进行声波场分解，或者说“反虚反射(deghosting)”，更具体地说，涉及反虚反射过程中对三维效应的估算。

背景技术

图1是海上地震测量示意图，其中地震能量由水体表面6下深度h处的被拖曳震源(例如空气枪串阵列)发射并由被拖曳传感器(如以距离d分开的一对拖缆阵列2、2’，每个阵列有多个水中检波器S₁、S₂......S_N，并悬挂在浮筒/浮标8下面)。震源1向水中传送声波，造成一个波场，该波场相干地传播到水下的地层中。当波场撞击地层间的界面4时，波场被反射并沿路径5穿过地层和水到达传感器，在传感器中它被转换成电信号并被记录下来。

在其他海上测量方法中，传感器和/或震源放在海底3或靠近海底3处，或者放在穿过地层的井中(也称作井筒或井孔)。通过分析这些检测到的信号，有可能确定海下地层的形状、位置和岩性。

在海上测量以及反演垂直地震剖面(VSP)中遇到的一个问题是水柱的混响问题。这个问题是水面和海底(以及海下面各边界)固有的反射性质造成的，对该问题可做如下解释。由海底或海下面地层反射的地震波通常沿向上方向进入水中。这个波称作“初至波”，它穿过海水并到达传感器。传感器或在海底或在被拖曳的阵列中，它记录下初至波的存在(即该初至波的特征)。该波场继续向上前进，例如沿路径7到达水的表面，在那里被反射向下传播。这一反射波场，或称作“虚反射(ghost)”波场，也通过海水传播并到达传感器，在那里被再次记录。取决于在水底地球物质的性质，虚反射波场本身可能被反射向上穿过海水，产生一个或多个次生虚反射构成的序列，或称作“多重波”。

在海底地球物质特别硬的情况中，由震源产生的过剩声能或噪声还能陷入水柱，以与被反射地震波本身相同的方式形成交混回响，这一噪声往往有大振幅，其结果是趋向于覆盖供研究用的较弱地震反射波。在水中的这种地震波场混响模糊了地震数据，放大了某些频率成分而衰减了其他频率成分，从而使得难于分析地下地层。所以，为准确表征地层，反虚反射(或者说去掉虚反射波场)是重要的。本领域技术人员将会理解，单是反虚反射并不能完全解决多重波问题(尽管有其他已知方法解决多重波问题)，因为每个多重波会有上行部分以及下行部分(它的虚反射)。

迄今提议的大多数反虚反射解决方案(例如J.O.A.Robertsson，J.E.Kragh和J.Martin，1999，在地震数据中减小海面虚反射污染的方法和系统，GB专利2363459号；J.O.A.Robertsson和J.E.Kragh，2002，使用单个拖缆和压力梯度近似进行强浪反虚反射，Geophysics，67，2005-2011；以及J.O.A.Robertsson，L.Amundsen，T.Roesten以及J.E.Kragh，2003，使用垂直质点速度近似进行地震数据的强浪反虚反射，国际专利申请PCT/GB2003/002305，2003年5月27日提出)忽略了三维(“3D”)效应。数据被认为是利用与拖缆成一线发生的震源事件(或称“放炮”)获取的，否则则进行预处理以满足这一判据。

然而，实际上3D效应可能是显著的，这是由于下述若干不同的理由：

1、获取数据时的几何布局是3D的，在一些拖缆和震源之间存在显著的横向偏移。

2、2D处理过程假定波前在空间呈柱状扩散，而3D中假定波前的空间扩散是球状的。

3、海面有3D结构，造成平面外的散射；以及

4、表面下可能存在显著的横向变化，造成平面外的反射和散射。

在以往提议的那些“2D”解决方案中，我们能成功地处置上述第1和第2个理由。事实上，对于理由1，一种小型的反虚反射滤波器(例如参见Robertsson和Kragh，2002；Robertsson等，2003；以及T.Roesten，L.Amundsen，J.O.A.Robertsson和E.Kragh，使用垂直质点速度场近似进行强浪反虚反射，第64届EAGE大会，Florens，2002；或L.Amundsen，T.Roesten，J.O.A.Robertsson和E.Kragh，使用压力梯度近似对单拖缆地震数据进行强浪反虚反射，已提交Geophysics，2003)理想地适用于将实际传播平面投影到包含拖缆数据记录位置的垂直平面，无需进行不规则空间处理。Robertsson和Kragh(2002)说明了如何补偿理由2，并得出结论：由假定柱状扩散所造成的误差大部分可以忽略。列出的最后两项(理由3和4)不能使用“2D”方案来解决。其中列出的第3项，即强浪表面的3D结构多半不像第4项那样重要。在有复杂成像任务的区域(盐、断层块等)或多重波问题(例如绕射多重波)中，至关重要的是能够适当地估算脱离平面的波传播。

所以，需要对伴随2D方案的缺点提出解决方案。

在已知的2D反虚反射技术(例如Robertsson和Kragh，2002；Robertsson等，2003；Amundsen等，2003)中与3D效应有关的最大误差项对应于压力的二阶横向空间导数。如果以三条传统的拖缆并排分布在奈奎斯特波数的几分之一内(间距几米)，则能使用三点滤波器来实现这一点。但由于实际操作上的原因，这难于实现。

所以需要实际可行的方法和/或装置用于在进行海上地震数据反虚反射时估算3D效应。

定义

在本申请中始终以被拖曳的海上地震测量为参考，其中一条或若干条拖缆被拖曳在具有一个或若干个震源的测量船后面。我们将指纵向测线方向为水平面内平行于拖缆的方向。还将使用直角坐标系，其中纵向测线方向是指x方向。相反地，我们将指横向测线方向为水平面内垂直于拖缆的方向。在直角坐标系中，横向测线方向是指y方向。

术语“多分量”是指使用一个或多个地震检波器(最好是三个正交的地震检波器以补偿方向灵敏度)检测质点速度信号，与检测压力梯度信号的水中检波器组合。为海底(也称海床)应用设计的一个商品化多分量系统实例是WesternGeco的Q-Seabed^TM系统。

发明内容

本发明涉及当对海上地震数据进行反虚反射处理时对3D效应的估算方法。这些方法依赖于能得到横向测线方向的二阶空间导数。这些导数能被直接测量(例如使三条拖缆很靠近并在水平面内彼此平行)，但更优选的作法是通过对其他波场量的间接测量和使用波动方程技术计算所希望的项来估计二阶纵向测线导数。

本发明方法优选利用在海面附近被拖曳的一条多分量拖缆，或者在海面附近的一对拖缆配置，或者是在海面下任意深度拖曳的三条拖缆配置，彼此沿垂直向上下排列。

对于三条拖缆的情况，或者直接测量或者用波动方程估计压力的二阶横向测线导数。

对于双拖缆配置的情况，压力的二阶横向测线导数是用波动方程估计的。

对于多分量拖缆或者双拖缆配置的情况，当然是对照不存在2D假设时直接测量的质点速度垂直分量或者使用双拖缆配置估计的质点速度垂直分量来校准只由单个拖缆压力数据得到的质点速度垂直分量，从而得到压力的二阶横向测线导数——这里假定波传播沿拖缆的纵向测线方向(Robertsson和Kragh，2002；Amundsen等，2003)。

根据本发明的另一方面，使用流体介质内的多分量传感器电缆获取地震波场数据。多分量传感器电缆可以由拖缆拖曳在流体介质表面附近，或者可以是与流体介质底床耦合的海底电缆。在拖缆的情况中(这里在其他地方提到)，多分量传感器电缆优选被拖曳在流体介质表面以下处于4-50米深度范围内。

附图说明

为了能详细理解本发明的上述特点和优点，通过参考附图中显示的本发明实施例，可以对上文概要说明的本发明有一个更具体的描述。然而，应该指出，附图显示的只是这一发明的典型实施例，所以不能被认为是对本发明范围的限制，因为本发明可以容许其他同样有效的实施例。

图1是根据本发明的一个方面所进行的海上地震测量的示意性展示。

图2显示一个解析3D反虚反射滤波器的图示，以及将三个不同的小型滤波器近似与该解析滤波器比较所得到的误差图。

图3是一个合成数据图的举例，代表使用2D有限差分模型产生的质点速度的垂直分量。

图4显示图3所示图形的细节部分，聚焦于一个特定的数据窗。

图5显示根据本发明使用双拖缆方案(配置1)得到的3D反虚反射基准结果与一个真实的上行垂直质点速度之间的差异图。

图6显示在图5的下部两图片中每道波形曲线的最大值被真实上行质点速度垂直分量最大值归一化后绘出的图形。

图7和图8显示与图5和图6相同的结果，但针对的是不同双拖缆(配置2)的情况。

图9、10、11和12显示与图5、6、7和8中所示图形类似的结果，但使用波动方程估计二阶水平导数。

图13以与图9和11类似的方式显示双拖缆配置1和2的3D反虚反射结果，但在对下拖缆进行滤波计算质点速度垂直分量时使用垂直入射近似，以避免单边差分近似。

图14和15显示使用双拖缆解决方案(分别为配置1和2)时每道波形曲线的反虚反射后数据与真实上行垂直质点速度之差的3D(归一化)最大值。

图16的图形代表使用双拖缆配置(配置1和2)时3D反虚反射的灵敏度，显示出当下拖缆相对于上拖缆沿横向测线方向横向偏移位移为3.125m时的效果。

图17是与图16相似的图形，但显示出下拖缆垂直偏移位移为0.3m(更深)时的效果。

图18显示根据本发明的多分量拖缆方案得到的3D反虚反射基准结果与一个真实的上行垂直质点速度之间的查意图。

图19显示在图18的下部两图片中每道波形曲线的最大值被真实上行质点速度垂直分量最大值归一化后绘出的图形。

图20、21、22和23显示与图18和19中所示结果类似的结果，但使用波动方程估计二阶水平导数，其中图20和21显示拖缆深度为4m时的结果，图22和23显示拖缆深度为6m时的结果

图24、25、26和27显示与图20、21、22和23中所示结果类似的结果，但显示出该算法对假定深度的扰动的灵敏性。

图28、29、30和31显示与图20、21、22和23中所示结果类似的结果，但使用显著浪高(SWH)为4m的强浪近似，以验证我们的确能处置强浪。

具体实施方式

3D效应和基本符号

本发明导出近似的三维反虚反射表达式，其中使用下述任何一种配置：

-三条拖缆并排分布在奈奎斯特波数几分之一内(间距几米)的一种配置，以直接测量压力的二阶横向测线导数；

-三条拖缆在海面下任意深度被拖曳，该配置中三条拖缆彼此垂直向上下排列；

-一条多分量拖缆；或

-双拖缆(即上/下一对拖缆)。

三条拖缆能以垂直分开配置在任意深度被拖曳，通过使用波动方程间接估计二阶横向测线导数。在使用多分量和双拖缆两种配置时，优选方式是在海面附近拖曳这些配置，因为它们利用海面近似来估计所希望的压力二阶横向测线空间导数。

声波场分解

可进行声波场分解(或反虚反射)，从而使结果量代表垂直质点速度的上行和下行分量或者压力。我们注意到，如果我们计算质点速度的上行垂直分量而不是压力，那么我们得到若干好处。最重要的一点是空间滤波器作用在压力记录上而不是在质点速度的垂直分量上。我们试图实现的正是对这个滤波器实现3D近似，在下面几节中我们将看到，由于该滤波器是作用在压力上，我们才能实现这一近似。应该指出，在下面对质点速度垂直分量进行的反虚反射中，反演到等效压力多少是不重要的，因为我们只处置上行波。其反演只是以一个常数除以入射角的余弦来进行标度，这能在fk域中容易地实现。

对质点速度垂直分量进行反虚反射的其他(次要)好处包括：

——能大大地减小地震干扰噪声，因为它大部分并不投射到质点速度的垂直分量上；以及

——声波分解中的空间滤波器包含一个零点替代一个极点，所以更容易实现。

下式用于分解数据(在频率-波数域中)：

$\left[1\right],V_Z^U=\frac{1}{2}(V_Z-\frac{k_Z}{\mathrm{\ρ\ω}}P),$

其中

是所希望的(反虚反射的)质点速度垂直分量的上行部分(大写字母代表波数-频率域中的表示)，V_Z是质点速度的垂直分量，P是压力，K_Z是垂直波数的绝对值，ω是角频率，ρ是水的密度。

已提出对式(1)的近似，其中假定kx＝0(垂直入射近似)，从而使式(1)简化为：

$\left[2\right],v_Z^U\≈\frac{1}{2}(v_Z-\frac{1}{\mathrm{\ρc}}P),$

其中c是水速，小些字母代表压力和质点速度的空间-频率域表示。该式将称作“垂直入射近似”，已经表明对于准确反虚反射不能令人满意。

式(1)能作为沿每条拖缆的空间滤波器很准确地实现。然而，实现式(1)要假定波只在2D平面内传播，波到达时不带有沿横向测线方向的分量。通过忽略这样的3D方面，我们引入了横向测线方向的“垂直入射近似”。

估计垂直波数的一种途径是利用水平波数近似

来计算一个空间小型滤波器(例如，Amundsen和Ikelle，1998；Roesten等，2002)。对水平波数的最简单的空间小型近似是在零波数附近的泰勒展开(对

有效)。在这里我们将重写对垂直波数的表述式并只在k_y方向展开它，因为已经能在纵向测线方向(即x方向)很准确地实现该滤波器：

$\left[3\right],k_z(\mathrm{\ω},k_x,k_y)\≈\sqrt{{(\mathrm{\ω}/c)}^2-{k_x}^2}(1-\frac{{k_y}^2}{2({(\mathrm{\ω}/c)}^2-{k_x}^2)}+O\left({k_y}^4\right))$

如上所述，式(3)使其自身作为空间小型滤波器来实现。这是因为水平波数k_y的每个方幂对应于与该波数方幂同阶的水平导数。然而，我们不必受限于垂直波数表达式的泰勒近似。Roesten等(2002)描述了也能造成小型滤波器并可能给出更好结果的其他更复杂的近似。我们关于小型滤波器的经验是：对于3点滤波器，最低阶泰勒近似难于有所改善。

式(3)包含两项。第一项与严格的2D近似完全相同，所以不包含k_y因子。第二项包含一个因子

它要求我们测量或估计横向测线方向的二阶空间导数或压力。

让我们用式(3)来检验对式(1)进行滤波器近似时好到什么程度[假定我们已知沿横向测线方向的二阶导数]。图2中的左上图显示解析性滤波器

的图形，而图2中其他图形显示的是误差图，表明三个不同的小型滤波器近似与该解析性滤波器比较所达到的近似程度。右上图描绘在最低阶近似(严格的2D近似)中的误差(与解析性滤波器比较)，该近似不包含k_y的任何因子。左下图描绘包括一个因子

的下一个较高阶近似。右下图描绘包括一个因子[在式(3)中未包括]的第二较高近似。黑圈表示与垂线成60度入射角。误差图是滤波器近似与真实解相比的差值相对于真实解归一化后得到的。当误差大于20％时，误差图已饱和。在黑圈内的区域对应于以与垂线相比入射角小于60度的入射能量。由图2清楚可见，通过包括一个有因子

的较高阶项，得到在横向测线方向对滤波器的好得多的近似。

多分量传感器电缆

Robertsson和Kragh(2002)提出一种单一拖缆反虚反射方法，该方法基于由在海面附近拖曳的一条传统拖缆记录压力来估计质点速度的垂直分量：

$\left[4\right],\frac{\∂p}{\∂z}=\frac{p}{h}-\frac{h}{2}{k}^{2}p-\frac{h}{2}(\frac{{\∂}^{2}p}{{\∂x}^2}+\frac{{\∂}^{2}p}{{\∂y}^2})+O\left(h^2\right)$

[Robertsson和Kragh，2002中的式(4)]，其中小写字母现在表示空间-频率域表达式，于是p是压力，k＝ω/c，h是在每个水中检波器位置的瞬时浪高作为时间的函数(需要在时间域中实现)，O(h2)表示展开式中的最大误差项。Amundsen等(2003)提出一种类似的方法，它对于稍大一些的拖缆深度是准确的：

$\left[5\right],\frac{\∂p}{\∂z}=k\cot \left(\mathrm{kh}\right)p-\frac{h}{3}[1+\frac{2}{15}{k}^2{h}^2](\frac{{\∂}^{2}p}{{\∂x}^2}+\frac{{\∂}^{2}p}{{\∂y}^2})+O\left(h^2\right)$

实现这两种方法时通常都假定相对于拖缆纵向测线2D波传播，从而忽略压力的二阶横向测线导数(见Robertsson和Kragh，2002，文章的附录)。通过将这一质点速度垂直分量估计值与在多分量拖缆中记录的实际值进行比较，我们能估计出在2D假定下的最大误差项。这里我们显示对于Amundsen等(2003)的方法如何进行这一误差估计，Amundsen等(2003)的方法给予我们在3D反虚反射近似中所需要的二阶横向测线导数。对Robertsson和Kragh(2002)的方法，也能导出类似的表达式。

这一解决方案不限于多分量拖缆电缆，也可通过其他多分量传感器电缆来实现，例如：与流体介质底床耦合的海底电缆(OBC)以及在流体介质表面附近被拖曳的多分量拖缆。在转让给WesternGeco的美国专利申请2004/0068373号中公开了一种适用的OBC电缆配置，特别是其中包括一个多分量传感器容器。在拖缆的情况中(这里在其他地方提到)，多分量传感器电缆最好是被拖曳在流体介质表面之下，处在4-50m的深度范围。根据相关领域普通技术人员的水平，这种拖缆最好是可由已知的多分量拖缆方案修改而成。在这方面参考下列多分量拖缆方案：转让给WesternGeco的美国专利6061302号；转让给I/O Exploration Products的美国专利5475652号；转让给Lister的美国专利4618949号；以及转让给Petroleum Geo-Services的专利申请公开20040042341号。

在多分量拖缆数据中，所记录的质点速度垂直分量v_z能通过下列运动方程转换成垂直压力梯度：

$\left[6\right],\frac{\∂p}{\∂z}={\mathrm{i\ω\ρv}}_z$

式(5)和式(6)给予我们的压力二阶横向测线导数是：

$\left[7\right],\frac{{\∂}^{2}p}{{\∂y}^2}=\frac{3}{1+\frac{2}{15}{k}^2{h}^2}[\frac{k\cot \left(\mathrm{kh}\right)}{h}p-\frac{\mathrm{i\ω\ρ}}{h}{v}_z]-\frac{{\∂}^{2}p}{{\∂x}^2}+O\left(h\right)$

式(7)是在空间-频率域中表示的。因为海面随时间变化，式(7)应利用空间中和时间中的小型滤波器在空间-时间域中实现(Amundsen等，2003)。这能通过对不包含纯因子k的项进行级数展开来实现(转换成时间导数)。这也将有助于解决由于在第一虚反射凹槽(反虚反射算子有效区之外)余切项趋于无限大所造成的稳定性问题。在下面的测试中，我们对照其泰勒展开中的前三项来替换式(7)中的余切项：

$\frac{k}{h}\cot \left(\mathrm{kh}\right)=\frac{1}{h^2}-\frac{k^2}{3}-\frac{k^4{h}^2}{45}$

最后，请注意，式(7)中最大量级的误差项与h的一次幂成正比。下面我们将测试这是否足够准确。

双拖缆配置

如我们在上文中看到的那样，为了对地震数据进行反虚反射处理，即把数据分解成上行和下行组分，必须知道质点速度的垂直分量。在多分量拖缆中这是直接测量的，然而在双拖缆配置中(例如见图1中的2、2’)，它是通过对这对拖缆的水中检波器数据进行滤波估计出来的(例如Robertsson等，2003)。这一滤波过程引入第二个2D假设，其中，在全3D滤波器的2D近似中的最大误差也是压力的一个二阶横向测线空间导数。

Robertsson等(2003)导出使用小型反虚反射滤波器，应用于在多个深度记录的拖缆数据，从而估计质点速度垂直分量的表达式：

[8]

$\mathrm{Vz}(\mathrm{\ω},k_x,k_y,z_1)=$

$\frac{i}{\mathrm{\ρ\ω\Δz}}\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}{(-1)}^m{F}_m^{\left(1\right)}{\mathrm{\κ}}^{2m}P(\mathrm{\ω},k_x,k_y,z_1)-\frac{i}{\mathrm{\ρ\ω\Δz}}\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}{(-1)}^m{F}_m^{\left(2\right)}{\mathrm{\κ}}^{2m}P(\mathrm{\ω},k_x,k_y,z_2)$

在式(8)中，P(ω，k_x，k_y，z₂)和P(ω，k_x，k_y，z₁)分别代表最深的和最浅的拖缆记录的压力，ρ是水的密度，Δz是拖缆垂直间距，ω是角频率，是垂直波数，k＝ω/c，c是水中的速度，κ²＝k_x ²+k_y ²，k_x和k_y是水平波数。在式(8)展开式中的最低阶系数是：

$\left[9\right],F_0^{\left(1\right)}=\mathrm{k\Δ}z\cot \left(\mathrm{k\Δz}\right)$

$\left[10\right],F_1^{\left(1\right)}=\frac{\mathrm{\Δz}}{2k}\frac{d}{d\left(\mathrm{k\Δz}\right)}\mathrm{k\Δ}z\cot \left(\mathrm{k\Δz}\right)$

${\left[11\right],F}_0^{\left(2\right)}=\frac{\mathrm{k\Δz}}{\sin \left(\mathrm{k\Δz}\right)}$

$\left[12\right],F_1^{\left(2\right)}=\frac{\mathrm{\Δz}}{k}\frac{\sin \left(\mathrm{k\Δz}\right)-\mathrm{k\Δ}z\cos \left(\mathrm{k\Δz}\right)}{{\sin }^2\left(\mathrm{k\Δz}\right)}$

在式(8)的级数展开中只保留最低阶项[式(9)和式(11)]所得到的表达式对所有拖缆间距都是有效的，而且只对垂直入射波是严格准确的。所以它等效于早期的双拖缆反虚反射中使用的所谓“移位-求和”技术。

在级数展开中保留较高阶项能显著改善其结果。下一个较高阶近似将包括如式(9)-(12)给出的四项。这能用3点滑动空间滤波器来实现。

请注意，式(8)对于3D情况是严格准确的。在本申请的其余部分中，我们将不再关注沿拖缆纵向测线所发生的事情。我们能在式(8)中使用必须的任意多项，因为在纵向测线(x)方向的高阶空间滤波器不是一个问题。反之，我们将焦点放在沿横向测线能得到什么结果。我们观察到，如果我们能测量或估计沿横向测线方向的压力二阶导数，我们就能使用上文提到的3点滤波器[系数(9)-(12)]。我们受到这一观察结果的鼓舞，因为我们知道这一近似对于处理纵向测线数据也是很准确和足够的。

如我们所看到的那样，在对双拖缆数据进行反虚反射过程中的两个步骤(估计质点速度的垂直分量和将这一估计与压力记录进行上/下分解组合)造成最大量级误差项，这些误差项包含横向测线方向的二阶空间导数。在双拖缆配置中，我们优选两条拖缆垂直分开几米之内。此外，我们还优选双拖缆配置位于海面附近(例如拖缆在海面下5m和8m深)。这将给我们三个已知压力的垂直深度(在海面的压力为零)，使我们在原则上能使用有限差分近似来估计压力的二阶垂直导数。然后，波动方程告诉我们如何计算横向测线的压力二阶导数(在上部拖缆的位置)：

$\left[13\right],\frac{{\∂}^{2}p}{{\∂y}^2}=\frac{{\∂}^{2}p}{{\text{\∂t}}^2}-\frac{{\∂}^{2}p}{{\∂x}^2}-\frac{{\∂}^{2}p}{{\∂z}^2}$

在另一个双拖缆配置实施例中，通过参照在上部拖缆处(在深度z₁)质点速度垂直分量v_z的式(4)或(5)来校准式(8)，可得到压力的二阶横向测线导数。这造成一个方程，但有两个未知数，即在上拖缆和下拖缆两个位置(在深度z₁、z₂)的压力二阶横向测线导数。然而，通过参照式(4)或式(5)重复校准式(8)能得到关于这同样两个未知数的第二个方程式，但这一次是对下拖缆(在z₂)处的v_z。所以我们是在使用方案2的获取配置(双拖缆)，但使用方案1(多分量拖缆)的处理方法，即不使用式(13)。

合成数据

为了评估本申请中描述的3D反虚反射技术，实际上使用2D合成就足够了。我们假定有在垂直于拖缆的2D平面内的模拟数据，而且假定在纵向测线方向已圆满完成反虚反射。所以，在2D合成中接收器相对于震源的不同位置将告诉我们，对于不同的入射角3D反虚反射在横向测线方向的作用如何。该反虚反射技术中唯一未被充分评估的部分是在式(7)和式(13)中减去含有压力二阶纵向测线导数

的项，我们假定该项是被准确估计的。这不是对测试的重要限制，因为这一项能以几乎与同一表达式中的时间导数同样高的精确度计算出来。如我们将会看到的那样，在式(13)中的最不准确项是沿垂直方向的二阶导数。至于式(7)，引入误差的原因是该表达式只准确到海面下拖缆深度h的一阶。

合成数据是使用Robertsson等(1994)的2D有限差分(FD)程序代码产生的。这一程序代码允许模拟强浪海面的效应(Robertsson，1996)。在本申请的主要部分，我们使用来自两套模拟的合成数据，这两组模拟过程是使用含有一层均匀水体的模型进行的。在第一模拟中，在模型的全部边缘使用吸收边界，从而使这些数据将只包含沿离开震源位置的一个方向传播的波(我们的参考数据代表“完美的”反虚反射)。在第二模拟中，模型顶面以平的海面为边界。到本申请的末尾，我们还使用当存在有效浪高(SWH)为4m的强浪海面时产生的合成数据。

在全部模拟中使用位于海面下150m的50Hz Ricker子波点源。

双拖缆合成

在本申请中我们将考虑两个不同的拖缆配置：

——配置1：浅拖缆在海面下3m，较深拖缆在海面下6m，垂直位于浅拖缆下方。

——配置2：浅拖缆在海面下6m，较深拖缆在海面下9m，垂直位于浅拖缆下方。

沿每条拖缆以3.125m采样间隔记录压力数据。该记录的水平延伸距离是相对于震源的水平位置(在本申请的图中位于x＝1860m处)+/-1500m。

多分量拖缆合成

在海面下4m和6m的拖缆以3.125m采样间隔记录压力和质点速度垂直分量数据。该记录的水平延伸距离是相对于震源的水平位置(在本申请的图中位于x＝1860m处)+/-1500m。

合成数据(质点速度的垂直分量)的实例示于图3和图4。图3描绘在FD模拟中记录的质点速度垂直分量。图4是图3中感兴趣数据窗的“放大”图，显示出在本申请中我们聚焦的特定数据窗(对双拖缆和多分量拖缆两种情况)。请注意跟随在直达波到达之后的来自下吸收边界的轻微反射。

结果：双拖缆方案

我们提出的用于3D反虚反射的两种方法都依赖于估计式(3)中需要的压力二阶横向测线空间导数。如我们在上文中看到的那样，双拖缆方案还需要压力的二阶横向测线空间导数来估计质点速度的垂直分量[式(8)-(12)]。我们首先通过直接由合成数据计算二阶水平导数而不是估计该导数来检验对于双拖缆方案我们能做到的最好程度。这将用作我们的基准方案以在评估反虚反射结果时判定我们能达到何种程度。

回顾一下我们有两种双拖缆配置要评估。在配置1中，浅拖缆是在海面下3m，而较深拖缆是在海面下6m，垂直位于浅拖缆下方。在配置2中，浅拖缆是在海面下6m，而较深拖缆是在海面下9m，垂直位于浅拖缆下方。

图5显示使用双拖缆配置1时利用反虚反射处理后的数据与真实质点速度上行垂直分量之间的差值给出的基准反虚反射结果。在该基准方案中，我们在FD模拟中使用直接计算出的压力二阶水平导数(这是我们可能的最好做法)。图5中的左手一列显示使用本申请中展现的3D反虚反射方案(即式(1)和(3))得到的结果。左下图是在左上图的放大。图5中的右手一列显示使用式(2)的“垂直入射近似”得到的结果，其右下图是右上图的放大。这是当应用2D纵向测线反虚反射时沿横向测线方向间接达到的结果。所以，该结果代表当使用2D反虚反射方案时沿横向测线方向的误差。

在图6中我们画出了图5下部两图中每道波形曲线的最大值被真实上行质点速度垂直分量最大值归一化后的结果。所以，这一图形代表反虚反射误差作为入射角的函数。实线曲线代表由式(2)中的“垂直入射近似”导出的2D反虚反射结果，而虚线曲线代表由使用式(1)和式(3)的反虚反射导出的3D反虚反射结果。应该指出，通过引入3D方案，我们已显著增大了成功进行反虚反射的区域。对于2D反虚反射，在10度入射角时已达到1.5％相对误差，而对于3D反虚反射方案，在40度入射角时才达到这一相对误差。

图7和图8显示与图5和图6相同的结果，但图7和图8是关于双拖缆配置2的(拖缆在水面下6m和9m)。这样，图7描绘使用双拖缆方案(配置2)时3D反虚反射的基准方案。在该基准方案中我们使用在FD模拟中直接计算出的压力二阶水平导数(这是我们可能做的最好方案)。全部四张图反映反虚反射处理后的数据与真实的上行垂直质点速度之间的差值。左上图描绘使用式(1)和式(3)的“3D”反虚反射，而在上图描绘使用式(2)中的“垂直入射近似”的“2D”反虚反射。左下图和右下图分别为左上图和右上图的放大。在图8中，两条曲线显示每道反虚反射处理后的数据与真实上行垂直质点速度之差的(归一化)最大值。实线描绘使用式(2)中的“垂直入射近似”进行的反虚反射。虚线描绘使用式(1)和式(3)进行的反虚反射。对于这一配置，我们要指出，通过引入3D方案，我们已显著增大了成功进行反虚反射的区域。配置1和配置2的结果总体上没有显著差别。

图9、10、11、12显示与图5、6、7、8中所示类似的结果，但这一次我们使用式(13)估计二阶水平导数。这些结果并不像基准方案的结果(图5-8)那样令人鼓舞(在基准方案中我们使用在FD模拟中直接计算出的二阶水平导数值)，而且不会为3D反虚反射提供足够准确的结果。

在反虚反射过程中，有三个地方使用压力的二阶水平导数。首先是在浅拖缆和较深拖缆处当估计质点速度垂直分量时需要它[式(8)-(12)]。当使用式(3)进行上/下分离时在浅拖缆位置也需要它。式(13)能用于准确地计算在较浅拖缆位置的二阶水平导数。然而，在较深拖缆位置其估计值将不那样准确，因为我们不能使用二阶垂直导数的中心有限差分估计(这将需要在下面有另一条拖缆)。替代地，我们使用不那么准确的单边FD近似计算图9-12所示反虚反射结果中的二阶垂直导数。通过多次测试，我们发现正是这一近似降低了3D反虚反射的准确度。

在图13中我们显示与图9和11完全相同的双拖缆配置1和2得到的3D反虚反射结果(即我们使用了式(13)估计二阶水平导数)，但这一次在对较低拖缆滤波以计算质点速度垂直分量时使用了垂直入射近似，以避免单边差分近似。在所有其他地方都使用了这一完整表达式。虽然这是偏离最准确反虚反射方案的一个折衷，但我们能看到现在的结果已得到显著改善。图14和15显示误差曲线作为入射角的函数。我们现在得到的方案是对标准的2D反虚反射方案的改进-至少是离开了垂直入射，而在垂直入射处2D方案似乎更好些。在垂直入射时的欠佳表现多少有些令人迷惑，但可能由合成或滤波过程中的数值计算假象造成。压力的二阶水平导数是用式(13)计算出来的，但在对较低拖缆滤波以计算质点速度垂直分量时使用垂直入射近似。上述两条曲线显示每道反虚反射处理后的数据与真实的上行垂直质点速度之间差值的(归一化)最大值。实线曲线描绘使用式(2)中的“垂直入射近似”进行的反虚反射。虚线曲线描绘使用式(1)和式(3)进行的反虚反射。在图15中，上述两条曲线显示每道反虚反射处理后的数据与真实上行垂直质点速度之间差值的(归一化)最大值。

接下来我们研究基于双拖缆配置的方法对扰动的敏感程度。首先我们研究两条拖缆沿横向侧线方向的横向偏移的效应。图16显示在配置1和2的情况下当下拖缆相对于上拖缆沿横向侧线方向有一个未被计入的3.125m横向偏移时得到的结果。请注意，这是一个相当显著的偏移，因为这两条拖缆的间距只有3m。图中左列描绘的是配置1，而右列描绘的是配置2。上面一行图形是下面一行图形的放大图。虚线和实线曲线显示每道反虚反射处理后的数据与真实的上行垂直质点速度之间差值的最大值。如这里的其他图形那样，实线曲线描绘的是使用式(2)中的“垂直入射近似”进行的反虚反射。虚线曲线描绘的是使用式(1)和式(3)进行的反虚反射。

在离开垂直入射时这一方案迅速失灵。然而，值得指出的是，在横向侧线方向，3D反虚反射方案决不比2D反虚反射方案差。尽管如此，我们的结论是，为了使用双拖缆配置在3D中对数据进行反虚反射处理，我们必须小心地关注两条拖缆沿横向侧线方向彼此相对位置。尽管两条拖缆可能不需要一条拖缆总是正好在另一条的上方，但希望准确地知道两条拖缆彼此相对位置，从而能对其予以考虑。

最后，图17显示在配置1和2的情况下当下部拖缆有比预期位置深0.3m的未被计入的垂直偏移(相对于拖缆间距有10％的扰动)时得到的结果。尽管该扰动显然对结果有负面影响，但我们注意到这一效应不像横向测线扰动的效应那样强烈。这可能也与引入的扰动幅度较小有关。

结果：多分量拖缆方案

3D反虚反射的多分量拖缆方案只依赖于估计压力的二阶横向测线空间导数供式(3)中使用。首先，我们再次通过直接由合成数据计算二阶水平导数而不是估计该导数来验证对于该反虚反射方案我们能做到的最好程度。这将用作我们的基准方案以在评估反虚反射结果时判定我们能达到何种程度。

图18显示使用多分量拖缆方案时利用作为基准的反虚反射处理后的数据与真实的上行质点速度垂直分量之间的差值给出的反虚反射结果。如前文所述，在该基准方案中我们在FD模拟中使用直接计算出的压力二阶水平导数(这是我们可能的最好做法)，图18中的左手一列显示使用本申请中展示的3D反虚反射方案得到的结果。左下图是左上图的放大。图18中的右手一列显示使用“垂直入射近似”得到的结果，或者说是当应用2D纵向侧线反虚反射时沿横向侧线方向间接达到的结果。其右下图是右上图的放大。

在图19中我们画出了图18下部两图中每道波形曲线的最大值被真实上行质点速度垂直分量最大值归一化后的结果。所以，这一图形代表反虚反射误差作为入射角的函数。实线曲线代表2D反虚反射结果，而虚线曲线代表3D反虚反射结果。应该指出，通过引入3D方案，我们已显著增大了成功进行反虚反射的区域。对于2D反虚反射，在10度入射角时已达到1.5％相对误差，而对于3D反虚反射方案，在40度入射角时才达到这一相对误差。

图20、21、22、23显示与图18和19中所示类似的结果，但这一次我们使用式(7)估计压力的二阶水平导数。在图20和21中我们显示拖缆深度为4m时的结果。这些结果极好，几乎如图18和19的基准方案一样好。图22和23显示拖缆深度为6m时的结果。这些结果也显示出比标准的2D反虚反射有显著改善，尽管与4m深拖缆的结果相比质量有所下降。在这些测试中使用50Hz Ricker子波，在高达110Hz时仍有显著的能量。如果使用频率成分稍低的子波，我们当然预期对于6m深的拖缆其结果会有所改善。

多分量拖缆对于在DC或较高凹槽处的虚反射不存在问题，因为压力和质点速度分量在各自的频率中有互补的凹槽。所以在什么深度拖曳拖缆是无关紧要的。然而，为了减小噪声，希望在尽可能大的深度拖曳拖缆。在平静天气，在例如4m深处拖曳拖缆是可能的，特别是如果这样做能使得以高频实施本申请中所示3D反虚反射方案的话。

图24、25、26、27显示与图20、21、22、23中所示类似的结果，但这一次我们在3D反虚反射处理中使用错误的拖缆深度，以研究该算法对扰动的敏感性。在图24和25中我们显示拖缆深度为4m的结果，但这里我们代之以假定深度为3.9m。在图26和27中我们显示拖缆深度为6m的结果，但这里我们代之以假定深度为5.9m。对于这一参数的0.1m小扰动，其结果质量却迅速下降，由此我们得出结论：拖缆在海面下的深度必须在这一距离内已知，以提供准确的3D反虚反射处理后的结果。再次使用式(7)计算压力的二阶水平导数。在图25和27中，曲线显示每道反虚反射处理后的数据与真实的上行垂直质点速度之间差值的(归一化)最大值，实线曲线描绘使用式(2)中的“垂直入射近似”进行的反虚反射，而虚线曲线描绘使用式(1)和式(3)进行的反虚反射。

已设计出以小型反虚反射滤波器实现的3D反虚反射算法，以去掉强浪效应。图28、29、30、31显示与图20、21、22、23所示类似的结果，但这一次我们使用有效浪高(SWH)为4m的强浪，以验证我们的确能处置强浪。我们的结论是：在存在强浪的情况下该算法也能工作得很好。图28和29描绘使用在平均海面下4m拖曳的多分量拖缆方案进行的3D反虚反射强浪测试(4m SWH)。图30和31描绘使用在平均海面下6m拖曳的多分量拖缆方案进行的3D反虚反射强浪测试(4m SWH)。

结论

在本申请中我们讨论了反虚反射的3D方面，并提出两种解决方案。通常所谓“3D”效应可细分为四组：

-在一些拖缆和震源之间存在显著横向测线偏移的获取3D数据几何布局。

-2D方法假定波前在空间中呈柱状扩散，而在3D中为球状扩散。

-海面有造成非平面散射的3D结构。

-在表面下可能存在显著的横向测线变化，造成非平面反射和散射。

在2D方法中，如在过去提出的那些方法(Robertsson和Kragh，2002；Robertsson等，2003；Amundsen等，2003)中，前两组能被成功地处置。所列举的后两组不能用“2D”方法解决。其中，所列举的第三组，即强浪海面的3D结构多半不像第四组那样重要。在有复杂成像任务(盐、断层块体等)或多重问题(如多重绕射)的领域中，能够适当地考虑非平面的波传播是至关重要的。本申请中提出的方法还可以提供用于这些领域的解决方案。

我们已经展现了用于3D反虚反射的三种解决方案：

-第一种技术要求的数据是由在海面附近(如在6m和9m深度)被拖曳的双拖缆(上/下)配置记录的。

-第二种技术要求使用也是在海面附近(如在6m深度)被拖曳的一条多分量拖缆记录的压力和质点速度垂直分量。

-第三种技术要求在海面下任何深度水平或垂直分开拖曳的三条拖缆，以便直接地或间接地(使用波动方程)估计压力的二阶横向测线导数。

这些技术包括完全强浪反虚反射。

这些方法依赖于可得到横向测线方向的二阶导数。虽然在原理上这些量能被直接测量(例如让三条拖缆在水平面内彼此平行地十分靠近)，但这种途径可能在操作上、后勤上或经济上并没有吸引力。替代地，我们提出通过间接测量其他波场量并使用波动方程技术计算所希望的项来估计二阶横向测线导数。

双拖缆技术给出比较好的结果，但发现该方法对扰动多少有些敏感。然而，虽然引入扰动会迅速降低3D方案的质量，但发现在相似条件下该3D方案总是如2D方案一样好(或优于2D方案)。另一方面，对于海面下4m和6m的拖缆深度，以及存在强浪的情况下，多分量拖缆技术都能给出极好的结果。粗略地说，我们得到在3D算法中对40度入射角的相对误差类似于2D算法在离开垂线10度时的相对误差。这是一个很显著的改进。我们还发现，为了得到高质量结果，必须在10cm误差范围内已知每个水中检波器在海面下的深度。

3D反虚反射方案对有效地衰减绕射多重波会特别有意义。Robertsson和Kostov(2003)概述了一种3D多重波压缩技术，它不要求获取3D数据的几何布局，而是基于逐道波形。然而，该方法假定数据已首先经过适当的反虚反射处理(3D)。

另一方面，3D反虚反射可能是使Amundsen去多重波方法(Amundsen，2001)对拖缆配置有效和处置绕射多重波的重要一步。我们还预期经2D反虚反射的数据和经3D反虚反射的数据之间的差别能用于标记具有横向测线分量的波至，从而能在成像过程中压缩它们。

由前文的描述将会理解，在本发明的优选的和替换的实施例中可以进行各种修改和变化而不脱离它的真正精神。

本描述只是为了以例说明，不应在限定的意义上加以解释。本发明的范围只应由跟随本描述的权利要求的语言来确定。在权利要求中的术语“包括”是指“至少包括”，从而在一个权利要求中给出的要素列表是一个开放组。“一(个)”或其他单数术语是要包括其复数形式，除非特别予以排除。

Claims (23)

1.一种对地震数据进行反虚反射处理的方法，包括如下步骤：

通过使用一个或多个传感器间接测量其他波场量而对海上地震波场数据确定压力二阶横向测线导数；以及

在分解地震波场数据的过程中使用所确定的压力二阶横向测线导数，其中所述确定步骤包括使用所述一个或多个传感器与海平面之间的距离。

2.权利要求1的方法，其中海上地震波场数据是使用放在流体介质中的至少一个地震源和至少一个地震接收器获取的。

3.权利要求1的方法，其中和确定的压力二阶横向测线导数一起使用波动方程技术，从而分解所获取的地震波场数据。

4.权利要求1的方法，其中地震波场数据是使用流体介质中的多分量传感器电缆获取的。

5.权利要求4的方法，其中多分量传感器电缆是在流体介质表面附近被拖曳的一条拖缆。

6.权利要求5的方法，其中多分量传感器电缆是在该流体介质表面下至少2米深度被拖曳。

7.权利要求6的方法，其中多分量传感器电缆是在该流体介质表面下不超过50米的深度被拖曳。

8.权利要求4的方法，其中多分量传感器电缆是与流体介质底床耦合的海底电缆。

9.权利要求4的方法，其中所获取的地震波场数据包括质点速度垂直分量。

10.权利要求9的方法，其中的确定步骤包括将所获取的质点速度垂直分量转换成垂直压力梯度。

11.权利要求10的方法，其中的转换是使用下列运动方程实现的：

$\frac{\∂p}{\∂z}={\mathrm{i\ω\ρv}}_z$

其中

ω是角频率，

ρ是水密度，以及

v_z是质点速度垂直分量。

12.权利要求11的方法，其中的确定步骤是根据下列方程实现的：

$\frac{{\∂}^{2}p}{{\∂y}^2}=\frac{3}{1+\frac{2}{15}{k}^2{h}^2}[\frac{k\cot \left(\mathrm{kh}\right)}{h}p-\frac{\mathrm{i\ω\ρ}}{h}{v}_z]-\frac{{\∂}^{2}p}{{\∂x}^2}+O\left(h\right)$

其中，

p是压力，

k＝ω/c，

ω是角频率，

c是水中速度，

ρ是水密度，

h是在每个水中检波器位置的作为时间的函数的瞬时浪高(需要在时间域实现)，以及

O(h²)代表展开式中的最大误差项。

13.权利要求1的方法，其中地震波场数据是使用至少3条彼此垂直分开的拖缆配置获取的。

14.权利要求13的方法，其中该至少3条拖缆彼此垂直分开至少2米。

15.权利要求13的方法，其中的确定步骤是使用下列波动方程实现的：

$\frac{{\∂}^{2}p}{{\∂y}^2}=\frac{{\∂}^{2}p}{{\∂t}^2}-\frac{{\∂}^{2}p}{{\∂x}^2}-\frac{{\∂}^{2}p}{{\∂z}^2}$

16.权利要求1的方法，其中地震波场数据是使用海面附近的双拖缆配置获取的。

17.权利要求16的方法，其中的确定步骤包括分别对双拖缆数据应用小型反虚反射滤波器，以确定在每条拖缆处的垂直压力梯度近似值。

18.权利要求17的方法，其中的确定步骤包括根据下列方程确定在每条拖缆位置的质点速度垂直分量：

$\mathrm{Vz}(\mathrm{\ω},k_x,k_y,z_1)=$

$\frac{i}{\mathrm{\ρ\ω\Δz}}\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}{(-1)}^m{F}_m^{\left(1\right)}{\mathrm{\κ}}^{2m}P(\mathrm{\ω},k_x,k_y,z_1)-\frac{i}{\mathrm{\ρ\ω\Δz}}\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}{(-1)}^m{F}_m^{\left(2\right)}{\mathrm{\κ}}^{2m}P(\mathrm{\ω},k_x,k_y,z_2)$

其中

P(ω，k_x，k_y，z₂)和P(ω，k_x，k_y，z₁)分别代表在最深拖缆和最浅拖缆处记录的压力，

ρ是水密度，

Δz是垂直拖缆间距，

ω是角频率，

是垂直波数，

k＝ω/c，

c是水中速度，

${\mathrm{\κ}}^2=k_x^2+k_y^2,$

以及k_x和k_y是水平波数。

19.权利要求16的方法，其中的确定步骤使用下列波动方程实现：

$\frac{{\∂}^{2}p}{{\∂y}^2}=\frac{{\∂}^{2}p}{{\∂t}^2}-\frac{{\∂}^{2}p}{{\∂x}^2}-\frac{{\∂}^{2}p}{{\∂z}^2}$

20.权利要求16的方法，其中双拖缆的深度相差不超过6米。

21.权利要求16的方法，其中，双拖缆中的最上部深度不大于6米。

22.权利要求1的方法，其中的确定步骤包括利用拖缆与海面的接近来确定所希望的压力二阶横向测线空间导数。

23.一种对地震数据进行反虚反射处理的方法，包括如下步骤：

使用3条拖缆确定海上地震波场数据的压力二阶横向测线导数，该3条拖缆并列排放并且空间间隔在奈奎斯特波数的几分之一内；以及

使用所确定的压力二阶横向测线导数和适当的技术来分解地震波场数据。

Images