CN101556339A - 对不规则接收机位置海洋地震拖缆数据进行消重影的方法 - Google Patents

Abstract

本发明为对不规则接收机位置海洋地震拖缆数据进行消重影的方法，包括获取海洋拖缆中各地震源启动的地震数据以及在地震数据被变换至谱域之后获取各频率的地震数据。使用基于物理的预条件子的迭代共轭梯度方案被应用于所变换的地震数据以提供消重影等式系统的标准等式组的最小二乘解。该解被逆变换回空－时域以提供消重影地震数据。

Description

对不规则接收机位置海洋地震拖缆数据进行消重影的方法

技术领域

本发明大致涉及地球物理探测领域。更具体地，本发明涉及对海洋地震拖缆数据进行消重影的领域。

背景技术

在油气产业，地球物理探测通常被用于帮助地层的搜寻和评价。地球物理探测技术产生对地球表面下构造的认识，而这有助于发现和提取有价值的矿物资源，特别是诸如石油和天然气之类的碳水沉积物。地球物理探测的熟知技术是地震勘测。在以陆基的地震勘测中，在地球表面上或者在地球表面附近生成地震信号，并且随后向下传输入地球的表面下部分。在海洋地震勘测中，地震信号也可以向下传输通过覆在地球的表面下部分之上的水体。地震能量源被用于生成地震信号，而地震信号在传播入地球之后至少部分由表面下地震反射体反射。该地震反射体一般是具有不同弹性特性(具体是声波速度和岩石密度)的地层之间的分界面，而不同的弹性特性导致分界面处的声阻抗不同。在地球表面处或者在地球表面附近、在上覆的水体中或者在钻孔中的已知深度处，反射地震能量由地震传感器(也叫做地震接收机)来检测并且记录。

在执行地震勘测中获得的所得到的地震数据被处理以产生关于地质构造和所勘测的区域中的地层特性的信息。所处理的地震数据被处理用于显示和分析这些地层的潜在的碳氢化合物内容。地震数据处理的目的是为了从地震数据中提取关于地层的尽可能多的信息，以充分地反映地质表面下部分。为了识别地球表面下部分中可能找到石油聚集的方位，大量的金钱被花费在收集、处理和解释地震数据上。从所记录的地震数据中构建定义所关心的地下地层的反射体表面的处理提供了在深度或时间上的地球图像。

产生地球的表面下部分构造的图像以使解释器能够选择具有石油聚集的最大可能性的方位。为了验证石油的存在，必须钻孔打井。钻孔打井以确定石油矿床是否存在是非常昂贵且耗费时间的工作。由于那个缘故，不断需要改进地震数据的处理和显示，以产生地球的表面下部分构造的图像，从而无论是由计算机或者人来作出解释，都将提高解释器的能力以评估石油聚集存在于地球的表面下部分中的特定方位处的可能性。

用于在地面地震勘测中生成地震信号的适当的地震源可以包括爆炸物或者振动器。海洋地震勘测一般采用水下地震源，所述水下地震源由轮船牵引并且周期地启动以生成声波场(acoustic wavefield)。生成波场的地震源可以有数种类型，包括小炸药包、电火花或电弧、海洋振动器并且一般地包括枪(gun)。地震源枪可以是水枪、蒸汽枪以及最一般地，可以是空气枪。一般地，海洋地震源不是由单个源元件组成的，而是由源元件的空间分布阵列组成的。特别地，该设置对空气枪更是如此，是当前海洋地震源的最普通的形式。

地震传感器的适当的类型一般包括质点(particle)速度传感器(特别在地面勘测中)和水压传感器(特别在海洋勘测中)。有时质点位移传感器、质点加速传感器或者压力梯度传感器被用于代替质点速度传感器或者与质点速度传感器一起使用。质点速度传感器和水压传感器在现有技术中通常分别被认为是地震检波器和水中听音器。地震传感器可以单独配置，但更通常以传感器阵列配置。另外，在海洋勘测中压力传感器和质点速度传感器可以一起配置(成对地排列或者以阵列对部署)。

在一般的海洋地震勘测中，地震勘测船一般以大约5节行驶在水面上，并且包含地震采集设备，例如导航控制、地震源控制、地震传感器控制以及记录设备。地震源控制设备使由地震船在水体中牵引的地震源在所选择的时期启动。地震拖缆(seismic streamer)(也称为地震电缆)是伸长的类似电缆的构造，由牵引地震源的地震勘测船或者由另一艘地震勘测船在水体中牵引。一般地，多个地震拖缆被牵引在地震船之后。地震拖缆包含传感器以检测由地震源发起的并从反射分界面中反射的反射波场。按照惯例，地震拖缆包含诸如水中听音器之类的压力传感器，但建议地震拖缆除了包含水中听音器之外，还包含诸如地震检波器之类的水质点速度传感器或者诸如加速计之类的质点加速传感器。压力传感器和质点运动传感器可以配置在附近、沿地震电缆成对排列或者以阵列对排列。

在反射波到达拖缆电缆之后，波继续传播到在水表面处的水/空气分界面，波从分界面向下反射并且再一次由拖缆电缆中的水中听音器所检测。水表面是较好的反射体并且对于压力信号而言，水表面处的反射系数幅度是几乎统一的并且符号是负的。在表面处反射的波因此相对于向上传播的波相位移动180度。由接收机所记录的向下传播的波通常被称为表面反射或者“重影”信号(ghost signal)。由于表面反射，水表面作用就像滤波器，在所记录的信号中产生谱凹口(spectralnotch)，使在所选择的带宽之外难以记录数据。由于表面反射的影响，所记录的信号中的一些频率被放大而一些频率被衰减。

最大衰减将发生在检测水中听音器和水表面之间的传播距离等于二分之一波长的频率处。最大放大将发生在检测水中听音器和水表面之间的传播距离是四分之一波长的频率处。声波的波长等于速度除以频率，而水中声波的速度大约为1500米/秒。因此，容易确定所产生的谱凹口的频谱中的方位。例如，对于在7米深度处的地震拖缆并且波垂直入射而言，最大衰减将发生在大约107Hz的频率处，而最大放大将发生在大约54Hz的频率处。

诸如地震检波器之类的质点运动传感器具有方向敏感性，但诸如水中听音器之类的压力传感器不具有方向敏感性。因此，由相互靠近的地震检波器和水中听音器所检测的上行波场信号将是同相的，而下行波场信号将被记录异相180度。已经提出各种技术用于使用该相差来减少由表面反射所引起的谱凹口，并且如果记录是在海底进行，还可以衰减水传播倍数(water borne multiple)。应该指出，使地震检波器和水中听音器共处在一地的可供选择的办法是具有足够的传感器空间密度，从而由水中听音器和地震检波器所记录的相应波场可以被内插或外推以产生处在同一方位的两个波场。

当前，在地球物理海洋地震探测的现有技术中，船牵引着非常长的拖缆，所述拖缆具有许多附接的地震接收机，一般为水中听音器。这些接收机注册源自地震源测深的散射声波场的一部分。由地震源所生成的声波场通过地球中的反射和折射而散射。因为这些拖缆非常长、具有许多接收机并且被牵引在移动的船之后，所以在航行方向或者在同线方向上的覆盖区域非常大。

在常规的海洋地震采集中，拖缆的接收机位于海表面下特定深度位置处的所谓阵列配置中。在固定水平面的该实施通过标准方法来完成。由于该设置，所谓的一次反射(从源到表面下部分并随后到接收机的直接响应)由所谓的虚反射(ghost reflection)而掩盖(cloud)，所述虚反射来自从源传输到表面下部分并随后经由海面而传输到接收机的波。

从海洋地震数据中去除虚反射是预处理数据以增加分辨力的第一阶段。该步骤被称为“消重影”。在当代的海洋地震采集中，记录了接收机的位置。对于平行于水面的拖缆中的接收机而言，可以由在所记录数据的谱域表示中所实现的相对简单的代数运算而实现消重影。此时，地震拖缆中的接收机的深度是唯一的物理输入参数。

标准消重影步骤是以接收机被放置在平行于海平面的水平面上为条件。在实践中，难以将具有接收机的长拖缆保持在恒定的垂直深度面上。进一步，有目的的非水平拖缆的例子就是所谓的倾斜拖缆。然而，当接收机被不规则地放置在水平方向和垂直方向上时，标准处理就会失效或者最多是给出近似的解。

因此，存在能够对在任意垂直和水平接收机位置处记录的海洋地震数据进行消重影的方法的需要。优选地，该方法应当能够在没有表面下部分(subsurface)的任何先验知识的情况下工作。

发明内容

本发明是对从具有不规则接收机位置的海洋拖缆所记录的地震数据进行消重影的方法。获取海洋拖缆中各地震源启动的地震数据和在地震数据被变换到谱域之后各频率的地震数据。使用基于物理的预条件子(physically-based preconditione)的迭代共轭梯度(iterativeconjugate gradient)方案被应用于所变换的地震数据，以提供消重影等式系统(system of equations)的标准等式组(normal set of equations)的最小二乘解。该解被逆变换回空-时域(space-time domain)以提供消重影地震数据。

附图说明

参考下面详细说明和附图，本发明及其优点可以更容易理解，在附图中：

图1是示例用于对在任意接收机位置处记录的海洋地震拖缆数据进行消重影的本发明第一实施例的流程图；

图2是示例用于对在任意接收机位置处记录的海洋地震拖缆数据进行消重影的本发明第二实施例的初始部分的流程图；

图3是示例起始于图2的用于对在任意接收机位置处记录的海洋地震拖缆数据进行消重影的本发明第二实施例的最后部分的流程图；以及

图4是示例用于迭代前承条件(preconditioned)共轭梯度方案的本发明实施例的流程图。

虽然本发明将结合其优选实施例来描述，但当然本发明不受限于此。相反，本发明旨在覆盖可以被包括在由所附上的权利要求所限定的本发明范围之内的所有备选方案、变型和等价物。

具体实施方式

本发明是对在任意接收机位置处记录的海洋地震拖缆数据进行消重影(deghosting marine seismic streamer data)的方法，其前提是深度坐标唯一相关于水平坐标。该前提排除了垂直接收机电缆。在现代海洋采集系统中，记录了该接收机位置信息。把消重影过程看作傅立叶类型等式系统的解说明该解来源于具有定义明确的傅立叶(kernel)的等式系统的逆(inversion)。事实上对于水平电缆的情形而言，等式系统减少为代数算子的简单情形。另一方面，对于具有不同深度面(在各水平位置为单值)的拖缆的例子而言，用等式系统的鲁棒的最小二乘解可以实现消重影，其中现在所有的接收机深度起物理输入参数的作用以建立该等式系统的核。

本发明工作在不具有表面下部分的任何先验知识的情况。在每次发射所记录的数据上分别执行消重影操作。对于每一个频率而言，其包括求解等式系统的方法。使用前承条件共轭梯度迭代方法而获得等式系统的解。预条件子(proconditioner)是基于当前问题的物理性质，而不是如常规的数学。作为有效的预条件子，水平拖缆(horizontalsteamer)的逆被采用，因为不规则非水平拖缆可以以规则水平拖缆的扰动来对待。不规则拖缆的特征是消重影步骤对频域中凹口处的基本零点信息的较低的不敏感性。

图1是示例用于对在任意接收机位置处所记录的海洋地震拖缆数据进行消重影的本发明第一实施例的流程图。

在方框11中，获取各地震源启动及被变换至谱域之后各频率的地震数据(seismic data are obtained for each seismic source activiationand for each frequency after being transformed to a spectral domain)。

在方框12中，使用基于物理的预条件子的迭代共轭梯度方案被应用于所变换的地震数据，以提供最小二乘解到消重影等式系统的标准等式组。

在方框13中，地震数据被逆变换回空-时域以提供消重影地震数据。

图2是示例用于对在任意接收机位置处记录的海洋地震拖缆数据进行消重影的本发明第二实施例的起始部分的流程图。参考图3中的流程图，最后部分将在下文中讨论。

在方框21中，获取用海洋拖缆记录的地震数据。各发射(shot)(地震源启动)位置的接收机数据和所有记录时间t的接收机数据被作为输入地震数据。笛卡尔坐标结构中接收机的空间位置给定为 $x^R=(x_1^R,x_2^R,x_3^R)$ 。在本发明中，接收机深度 $x_3^R=x_3^R(x_1^R,x_2^R)$ 是水平坐标x₁ ^R和x₂ ^R的单值函数。因此，拖缆不是垂直的，但不必定是水平的。进一步，在本发明中，水平坐标x₁ ^R和x₂ ^R可以不规则放置。因此，拖缆中接收机的水平位置不必定以规则的网格图形排列。

在方框22中，从方框21中获取的地震数据中获取发射的所记录的地震数据。

在方框23中，来自方框22的发射的所记录的地震数据被暂时从空-时域(space-time domain)变换为空-频域(space-frequency domain)。在本发明的实施例中，接收机方位处的散射波场p^sct被从空时(x₁ ^R，x₂ ^R，x₃ ^R，t)域变换至空频(x₁ ^R，x₂ ^R，x₃ ^R，s)域，因此：

$p^{\mathrm{sct}}=p^{\mathrm{sct}}(x_{1,q}^R,x_{2,r}^R,x_3^R;s)---\left(1\right)$

其中，

s＝jω＝j2πf

并且p^sct是空-频域中的散射声波场，x_1，q ^R是同线(in-line)接收机坐标，q是同线接收机数量，x_2，r ^R是交叉线(cross-line)接收机坐标，r是交叉线接收机数量，x₃ ^R是作为x_1，q ^R和x_2，r ^R函数的接收机深度，s是复拉普拉斯变换参数，j是虚部

，ω是角频率而f是频率。

对于各发射位置而言，完成了接收机方位处的散射波场p^sct的此变换，并且接收机区域处的散射波场p^sct的该变换可以通过应用熟知的变换(例如标准傅立叶变换)而实现。

在方框24中，对于频率f获取来自方框23的所变换的地震数据。

在方框25中，来自方框24的频率的所变换的地震数据被从空-频域变换至谱域。在本发明的实施例中，散射波场p^sct被从空-频域(x₁ ^R，x₂ ^R，x₃ ^R，s)变换至谱域(sα_n，sβ_m，x₃ ^R，s)，通过：

$P_{n,m}^{\mathrm{sct}}={\mathrm{\Δx}}_1^R{\mathrm{\Δx}}_2^R\underset{q=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\Σ}}}\underset{r=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\Σ}}}\exp ({\mathrm{js\α}}_n{x}_{1,q}^R+{\mathrm{js\β}}_m{x}_{2,r}^R)p^{\mathrm{sct}},---\left(2\right)$

其中，

sα_n＝nΔ(sα)，sβ_m＝mΔ(sβ)，    (3)

并且：

$\mathrm{\Δ}\left(\mathrm{s\α}\right)=\frac{2\mathrm{\π}}{{\mathrm{N\Δx}}_1^R},\mathrm{\Δ}\left(\mathrm{s\β}\right)=\frac{2\mathrm{\π}}{{\mathrm{M\Δx}}_2^R}.---\left(4\right)$

此时，P_n，m ^sct是谱域中的散射声波场，n是同线谱数量，m是交叉线谱数量，Δx₁ ^R是同线接收机采样距离，Δx₂ ^R是交叉线接收机采样距离，jsα_n是同线谱傅立叶参数，Δ(sα)是同线谱采样距离。jsβ_m是交叉线谱傅立叶参数，Δ(sβ)是交叉线谱采样距离，N是同线接收机的总数，而M是交叉线接收机的总数。

在方框26中，等式系统被确定用于对来自方框25的频率的所变换的地震数据进行消重影。因为接收机方位(location)是任意安置的，所以谱域中的波场不能明确地被确定。因此，本发明的消重影步骤包括求解等式系统。对于谱域中的任意接收机位置 $x^R=(x_{1,q}^R,x_{2,r}^R,x_{3,q,r}^R)$ 而言，在非水平拖缆上的接收机的深度由水平接收机坐标的单值函数来给定。具体地，接收机深度被写为 $x_{3,q,r}^R=x_3^R(x_{1,q}^R,x_{2,r}^R)$ 。在现代海洋采集系统中，记录了该接收机位置信息。消重影操作随后可以简洁地写为等式系统的解：

$\underset{n=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\Σ}}}\underset{m=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\Σ}}}{K}_{q,r;n,m}{P}_{n,m}^{\mathrm{dgh}}=p_{q,r}^{\mathrm{sct}},q=-\frac{1}{2}N+1,...,\frac{1}{2}N,r=-\frac{1}{2}M+1,...,\frac{1}{2}M,---\left(5\right)$

其中，系统矩阵(system matrix)K_{q，r；n，m}被给定为

$K_{q,r;n,m}=\frac{\mathrm{\Δ}\left(\mathrm{s\α}\right)\mathrm{\Δ}\left(\mathrm{s\β}\right)}{{4\mathrm{\π}}^2}\underset{n=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\Σ}}}\underset{m=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\Σ}}}\exp ({-\mathrm{js\α}}_n{x}_{1,q}^R-\mathrm{js}{\mathrm{\β}}_m{x}_{2,r}^R)2\sinh \left({\mathrm{s\Γ}}_{n,m}{x}_3^R\right).---\left(6\right)$

此处，未知矢量P_n，m ^dgh是零深度处的消重影的声波场，并且在谱域中给定为：

$p_{n,m}^{\mathrm{dgh}}=p_{n,m}^{\mathrm{dgh}}({\mathrm{s\α}}_n,s{\mathrm{\β}}_m,0;s),---\left(7\right)$

已知的矢量p_q，r ^sct是在空-频域中所记录的地震数据，给定为：

$p_{q,r}^{\mathrm{sct}}=p^{\mathrm{sct}}(x_{1,q}^R,x_{2,r}^R,x_{3,q,r}^R;s)---\left(8\right)$

并且垂直传播系数sΓ_n，m是：

${\mathrm{s\Γ}}_{n,m}=j\sqrt{{\left(\frac{2\mathrm{\πf}}{c}\right)}^2-{\left({\mathrm{s\α}}_n\right)}^2-{\left({\mathrm{s\β}}_m\right)}^2},---\left(9\right)$

其中，c是在水中的声波速度。

在方框27中，该过程继续到图3的方框31以求解方框26中所确定的等式系统。

图3是示例用于对在任意接收机位置处所记录的海洋地震拖缆数据进行消重影的本发明第二实施例最后部分的流程图。参考图2中的流程图，最后部分继续上面所讨论的起始过程。

在方框31中，标准等式组被从图2的方框26中的消重影联立方程中确定。因为已知的矢量p_q，r ^sct在空-频域中定义，并且未知的矢量P_n，m ^dgh在谱域中定义，所以在本发明的实施例中采用空-频域中的最小二乘解。等式(5)的两边都乘以系统矩阵K_{q，r；n′，m′}的复共轭，而该结果在同线和交叉线接收机数量(分别是q和r)的所有值上求和。在q和r上的较后的求和与n和m上的求和相互交换之后，获得了标准等式系统，如下：

$\underset{n=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\Σ}}}\underset{m=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\Σ}}}{A}_{n^{\′},m^{\′};n,m}{P}_{n,m}^{\mathrm{dgh}}=B_{n^{\′},m^{\′}},n^{\′}=-\frac{1}{2}N+1,...,\frac{1}{2}N,m^{\′}=-\frac{1}{2}M+1,...,\frac{1}{2}M.---\left(10\right)$

此处，矩阵A_{n′，m；n，m}给定为：

$A_{n^{\′},m^{\′};n,m}=\underset{q=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\Σ}}}\underset{r=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\Σ}}}\stackrel{\‾}{K_{q,r;n^{\′},m^{\′}}}{K}_{q,r;n,m},---\left(11\right)$

并且已知矢量B_n′，m′给定为

$B_{n^{\′},m^{\′}}=\underset{q=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\Σ}}}\underset{r=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\Σ}}}\stackrel{\‾}{K_{q,r;n^{\′},m^{\′}}}{p}_{q,r}^{\mathrm{sct}},---\left(12\right)$

其中n和n’是同线谱数量，m和m’是交叉线谱数量，并且系统矩阵K_{q，r；n′，m′}上的跨线(overline)指明复共轭。

本发明中的主要处理包括等式(10)中该线性标准等式组的求解。该求解具有现有拖缆配置的一般本质，但可以适应新的拖缆配置，例如倾斜拖缆。

需要注意，常规的海洋拖缆(严格水平的拖缆配置)只是由本发明所处理的一般情形的特殊情形。在该常规的(水平的)情形中，接收机深度位置x₃ ^R是常数并且等式系统的求逆减少为标准消重影步骤。等式系统可以通过逆离散傅立叶变换来求解。等式(5)和(6)成为等式系统：

$\frac{\mathrm{\Δ}\left(\mathrm{s\α}\right)\mathrm{\Δ}\left(\mathrm{s\β}\right)}{{4\mathrm{\π}}^2}\underset{n=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\Σ}}}\underset{m=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\Σ}}}\exp ({-\mathrm{js\α}}_n{x}_{1,q}^R-\mathrm{js}{\mathrm{\β}}_m{x}_{2,r}^R)2\sinh \left({\mathrm{s\Γ}}_{n,m}{x}_3^R\right)P_{n,m}^{\mathrm{dgh}}$

${=p}^{\mathrm{sct}}(x_{1,q}^R,x_{2,r}^R,x_3^R;s),q=-\frac{1}{2}N+1,...,\frac{1}{2}N,r=-\frac{1}{2}M+1,...,\frac{1}{2}M.---\left(13\right)$

常数x₃ ^R的消重影可以在谱域中显式地实现。零深度面处的谱域中的消重影域(如等式(7)中定义的p_n，m ^dgh)可以直接获得为：

$P_{n,m}^{\mathrm{dgh}}=\frac{P_{n,m}^{\mathrm{sct}}}{2\sinh \left({\mathrm{s\Γ}}_{n,m}{x}_3^R\right)}.---\left(14\right)$

对于等式(10)中一般线性等式组的求解而言，许多有效的求解器(solver)可用，例如直接求解器和迭代求解器。非常有效的迭代求解器是基于最小化所测量数据和建模响应之间的差异的前承条件共轭梯度迭代(preconditioned conjugate gradient iterative)方法。在本发明的进一步实施例中，考虑等式(5)并且谱消重影波场的解P_n，m ^dgh被采用共轭梯度迭代方案而展开，从而误差ERR给定为：

$\mathrm{ERR}=\underset{q=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\Σ}}}\underset{r=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\Σ}}}|p_{q,r}^{\mathrm{sct}}-\underset{n=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\Σ}}}\underset{m=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\Σ}}}{K}_{q,r;n,m}{P}_{n,m}^{\mathrm{dgh}}{|}^2---\left(15\right)$

被最小化。

返回图3，在方框32中，来自方框31的标准等式被用基于物理的预条件子而进行预处理(precondition)。预条件子被经常使用在共轭梯度迭代方法中，以使等式较好地经调整处理(well-conditioned)。例如，诸如多级ILU类型(multi-level ILU type)之类的标准预条件子通常基于系统矩阵的数学结构。不幸的是，对于当前问题而言，这些基于数学的预条件子是低效的。在本发明的实施例中，使用了预条件子，而该预条件子是基于所解决问题的物理性质(非水平拖缆)。

因为不规则非水平拖缆可以被认为是规则水平拖缆的扰动(perturbation)，所以不规则拖缆的系统矩阵的逆由水平拖缆的系统矩阵的逆来近似。因此，等效的水平拖缆的系统矩阵的逆被当作预条件子。需要注意，如果拖缆是水平的，那么共轭梯度迭代方案在具有该预条件子的单次迭代之内终止。对于水平拖缆而言，等式(10)的标准等式系统是对角的并且其解可以直接获得为：

$P_{n,m}^{\mathrm{dgh}}={\left[A_{n,m;n,m}\right]}^{-1}{B}_{n,m}---\left(16\right)$

在本发明的实施例中，所采用的预条件子为：

${\tilde{P}}_{n,m}={\left[A_{n,m;n,m}\right]}^{-\frac{1}{2}}={\left(\underset{q=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\Σ}}}\underset{r=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\Σ}}}\right|K_{q,r;n,m}{|}^2)}^{-\frac{1}{2}}---\left(17\right)$

为了方便起见，等式(5)可以以简洁的方式写出：

$\underset{n}{\mathrm{\Σ}}{K}_{q,n}{P}_n^{\mathrm{dgh}}=p_q^{\mathrm{sct}},---\left(18\right)$

其中，现在n代表(n，m)，而q代表(q，r)。

在方框33中，前承条件共轭梯度迭代方案被应用于求解等式(18)中的等式系统(来自方框32的前承条件标准等式组)。在各迭代中，用数字k来指明，获得具有残差(residual error)r_q ^(k)的谱消重影波场P_n ^dgh(k)的近似值。因此，在第k次迭代中，第k次残差被给定为：

$r_q^{\left(k\right)}=p_q^{\mathrm{sct}}-\underset{n}{\mathrm{\Σ}}{K}_{q,n}{P}_n^{\mathrm{dgh}\left(k\right)}.---\left(19\right)$

参考图4中的流程图的讨论，下文描述了用于前承条件共轭梯度迭代方案的本发明的特定实施例。

在方框34中，确定在图2的方框24中所选择的任何频率是否继续存在(remain)。如果频率继续存在，那么过程返回至图2的方框24。如果没有频率继续存在，那么过程继续至步骤35。因此，对于各频率f而言，通过共轭梯度迭代的求解被重复。

在方框35中，确定在图2的方框22中所选择的任何发射是否继续存在。如果发射继续存在，该过程返回到图2的方框22。如果没有发射继续存在，那么该过程继续至步骤36。因此，对各发射而言，通过共轭梯度迭代的求解被重复。

在方框36中，来自方框33的消重影解被从谱域变换回至空-时域。在获取谱域中的消重影波场之后，空-时域中零面(zero level)处的消重影波场 $p^{\mathrm{dgh}}=p^{\mathrm{dgh}}(x_{1,q}^R,x_{2,q}^R,0;s)$ 被获得为：

$p^{\mathrm{dgh}}=\frac{\mathrm{\Δ}\left(\mathrm{s\α}\right)\mathrm{\Δ}\left(\mathrm{s\β}\right)}{{4\mathrm{\π}}^2}\underset{n=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\Σ}}}\underset{m=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\Σ}}}\exp (-\mathrm{js}{\mathrm{\α}}_n{x}_{1,q}^R-\mathrm{js}{\mathrm{\β}}_m{x}_{2,r}^R)P_{n,m}^{\mathrm{dgh}}.---\left(20\right)$

图4是示例用于迭代前承条件共轭梯度方案的本发明实施例的流程图。这参照于上面图3的方框33中的方案。

在方框41中，选择了迭代计数器k＝0的零起始估计。在一个实施例中，该估计是：

$P_n^{\mathrm{dgh}\left(0\right)}=0,r_q^{\left(0\right)}=p_q^{\mathrm{sct}}.---\left(21\right)$

在其它实施例中，可以使用共轭梯度方案的起始估计的不同选择。例如，不同选择可以包括但不限定于前一频率的所获得的解或者不同发射位置的解。

在方框42中，计算了迭代计数器k＝1的估计。在一个实施例中，这些估计是：

$g_n^{\left(1\right)}={\tilde{P}}_n\underset{q}{\mathrm{\Σ}}\stackrel{\‾}{K_{q,n}}{r}_q^{\left(0\right)},$ $v_n^{\left(1\right)}={\tilde{P}}_n{g}_n^{\left(1\right)},$

${\mathrm{Kv}}_n^{\left(1\right)}=\underset{n}{\mathrm{\Σ}}{K}_{q,n}{v}_n^{\left(1\right)},$ ${\mathrm{\α}}^{\left(1\right)}=\frac{{\mathrm{\Σ}}_n|g_n^{\left(1\right)}{|}^2}{{\mathrm{\Σ}}_q|{\mathrm{Kv}}_n^{\left(1\right)}{|}^2},$

$P_n^{\mathrm{dgh}\left(1\right)}=P_n^{\mathrm{dgh}\left(0\right)}+{\mathrm{\α}}^{\left(1\right)}{v}_n^{\left(1\right)},$ $r_q^{\left(1\right)}=r_q^{\left(0\right)}-{\mathrm{\α}}^{\left(1\right)}{\mathrm{Kv}}_n^{\left(1\right)}.---\left(22\right)$

在方框43中，迭代计数器k增加了1。因此，k＝k+1，使得按顺序为k＝2，3，...。

在方框44中，使用了迭代计数器k-1的估计而计算迭代计数器k的估计。在一个实施例中，该估计是：

$g_n^{\left(k\right)}={\tilde{P}}_n\underset{q}{\mathrm{\Σ}}\stackrel{\‾}{K_{q,n}}{r}_q^{(k-1)},$ $v_n^{\left(k\right)}={\tilde{P}}_n{g}_n^{\left(k\right)}+\frac{{\mathrm{\Σ}}_n|g_n^{\left(k\right)}{|}^2}{{\mathrm{\Σ}}_n|g_n^{(k-1)}{|}^2}{v}_n^{\left(k\right)},$

$K{v}_n^{\left(k\right)}=\underset{n}{\mathrm{\Σ}}{K}_{q,n}{v}_n^{\left(k\right)},$ ${\mathrm{\α}}^{\left(k\right)}=\frac{{\mathrm{\Σ}}_n|g_n^{\left(k\right)}{|}^2}{{\mathrm{\Σ}}_q|{\mathrm{Kv}}_n^{\left(k\right)}{|}^2},$

$P_n^{\mathrm{dgh}\left(k\right)}=P_n^{\mathrm{dgh}(k-1)}+{\mathrm{\α}}^{\left(k\right)}{v}_n^{\left(k\right)},$ $r_q^{\left(k\right)}=r_q^{(k-1)}-{\mathrm{\α}}^{\left(k\right)}K{v}_n^{\left(k\right)}.---\left(23\right)$

在方框45中，计算了来自方框44的残差r_q ^(k)的当前第k次迭代估计的误差准则。在一个实施例中，该误差准则是：

$\frac{{\mathrm{\&Sigma;}}_q|r_q^{\left(k\right)}{|}^2}{{\mathrm{\&Sigma;}}_q|p_q^{\mathrm{sct}}{|}^2}<{10}^{-3}.---\left(24\right)$

在方框46中，确定方框45中的误差准则是否被满足。如果不满足误差准则，那么迭代过程返回至方框43以继续另一估计组。如果满足误差准则，那么迭代过程结束。

应该明白，前述仅是本发明特定实施例的详细说明并且根据这里所公开的内容在不背离本发明范围的基础上，可以作出对所公开的实施例的许多变化、变型和替代。因此，前述说明不是为了限定本发明的范围。相反，本发明的范围将只由附上的权利要求及其等价物所确定。

Claims (12)

1.一种用于对从具有不规则接收机位置的海洋拖缆中记录的地震数据进行消重影的方法，包括：

获取各地震源启动的以及被变换至谱域之后的各频率的地震数据；

将使用基于物理的预条件子的迭代共轭梯度方案应用于所变换的地震数据，以提供消重影等式系统的标准等式组的最小二乘解；以及

将所述解逆变换回空-时域以提供消重影地震数据。

2.根据权利要求1所述的方法，其中所述获取地震数据包括：

获取用海洋拖缆记录的地震数据；

获取地震源启动用的所记录的地震数据；

将所述地震源启动的所记录的地震数据从所述空-时域变换至空-频域；

获取频率的所变换的地震数据；以及

将频率的所变换的地震数据从所述空-频域变换至所述谱域。

3.根据权利要求2所述的方法，其中所述变换所记录的地震数据包括：

使用傅立叶变换将散射声波场从所述空-时(x₁ ^R，x₂ ^R，x₃ ^R，t)域变换至所述空-频(x₁ ^R，x₂ ^R，x₃ ^R，s)域，因此：

$p^{\mathrm{sct}}=p^{\mathrm{sct}}(x_{1,q}^R,x_{2,r}^R,x_3^R;s)$

其中，

s＝jω＝j2πf，

并且p^sct是所述空-频域中的所述散射声波场，x_1，q ^R是同线接收机

并且p^sct是所述空-频域中的所述散射声波场，x_1，q ^R是同线接收机坐标，q是同线接收机数量，x_2，r ^R是交叉线接收机坐标，r是交叉线接收机数量，x₃ ^R是作为x_1，q ^R和x_2，r ^R函数的接收机深度，s是复拉普拉斯变换参数，j是虚部

ω是角频率，而f是频率。

4.根据权利要求3所述的方法，其中所述变换所变换的地震数据包括：

将所述散射波场p^sct从所述空-频域(x₁ ^R，x₂ ^R，x₃ ^R，s)变换至所述谱域(sα_n，sβ_m，x₃ ^R，s)，这通过如下实现：

$P_{n,m}^{\mathrm{sct}}={\mathrm{\&Delta;x}}_1^R{\mathrm{\&Delta;x}}_2^R\underset{q=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\&Sigma;}}}\underset{r=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\&Sigma;}}}\exp (\mathrm{js}{\mathrm{\&alpha;}}_n{x}_{1,q}^R+\mathrm{js}{\mathrm{\&beta;}}_m{x}_{2,r}^R)p^{\mathrm{sct}},$

其中：

sα_n＝nΔ(sα)，sβ_m＝mΔ(sβ)，

并且：

$\mathrm{\&Delta;}\left(\mathrm{s\&alpha;}\right)=\frac{2\mathrm{\&pi;}}{\mathrm{N\&Delta;}{x}_1^R},$ $\mathrm{\&Delta;}\left(\mathrm{s\&beta;}\right)=\frac{2\mathrm{\&pi;}}{\mathrm{M\&Delta;}{x}_2^R},$

并且P_n，m ^sct是所述谱域中的所述散射声波场，n是同线谱数量，m是交叉线谱数量，Δx₁ ^R是同线接收机采样距离，Δx₂ ^R是交叉线接收机采样距离，jsα_n是同线谱傅立叶参数，Δ(sα)是同线谱采样距离，jsβ_m是交叉线谱傅立叶参数，Δ(sβ)是交叉线谱采样距离，N是同线接收机的总数，而M是交叉线接收机的总数。

5.根据权利要求4所述的方法，还包括：

确定用于对频率的所变换的地震数据进行消重影的等式系统；

确定所述消重影等式系统的标准等式组；以及

将基于物理的预条件子应用于所述标准等式组。

6.根据权利要求5所述的方法，其中所述确定用于消重影的等式系统包括应用如下等式：

$\underset{n=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\&Sigma;}}}\underset{m=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\&Sigma;}}}{k}_{q,r;n,m}{P}_{n,m}^{\mathrm{dgh}}=p_{q,r}^{\mathrm{sct}},$ $q=-\frac{1}{2}N+1,...,\frac{1}{2}N,$ $r=-\frac{1}{2}M+1,...,\frac{1}{2}M,$

其中，所述系统矩阵K_{q，r′，n，m}给定为：

$K_{q,r;n,m}=\frac{\mathrm{\&Delta;}\left(\mathrm{s\&alpha;}\right)\mathrm{\&Delta;}\left(\mathrm{s\&beta;}\right)}{{4\mathrm{\&pi;}}^2}\underset{n=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\&Sigma;}}}\underset{m=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\&Sigma;}}}\exp (-\mathrm{js}{\mathrm{\&alpha;}}_n{x}_{1,q}^R-\mathrm{js}{\mathrm{\&beta;}}_m{x}_{2,r}^R)2\sinh \left(s{\mathrm{\&Gamma;}}_{n,m}{x}_3^R\right).$

7.根据权利要求6所述的方法，其中所述确定标准等式组包括应用下面等式：

$\underset{n=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\&Sigma;}}}\underset{m=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\&Sigma;}}}{A}_{n^{\&prime;},m^{\&prime;};n,m}{P}_{n,m}^{\mathrm{dgh}}=B_{n^{\&prime;},m^{\&prime;}},$ $n^{\&prime;}=-\frac{1}{2}N+1,...,\frac{1}{2}N,$ $m^{\&prime;}=-\frac{1}{2}M+1,...,\frac{1}{2}M,$

其中，所述矩阵A_{n′，m′，n，m}给定为：

$A_{n^{\&prime;},m^{\&prime;};n,m}=\underset{q=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\&Sigma;}}}\underset{r=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{K_{q,r;n^{\&prime;},m^{\&prime;}}}{K}_{q,r;n,m},$

已知的矢量B_n′，m′给定为：

$B_{n^{\&prime;},m^{\&prime;}}=\underset{q=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\&Sigma;}}}\underset{r=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{K_{q,r;n^{\&prime;},m^{\&prime;}}}{p}_{q,r}^{\mathrm{sct}},$

并且n和n’是同线谱数量，m和m’是交叉线谱数量，并且跨线指明复共轭。

8.根据权利要求6所述的方法，其中所述预条件子包括水平拖缆的系统矩阵的逆。

9.根据权利要求8所述的方法，其中所述应用预条件子包括应用如下等式：

${\tilde{P}}_{n,m}={\left[A_{n,m;n,m}\right]}^{-\frac{1}{2}}={\left(\underset{q=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\&Sigma;}}}\underset{r=-\frac{1}{2}M+1}{\overset{\frac{1}{2}N}{\mathrm{\&Sigma;}}}{\left|K_{q,r;n,m}\right|}^2\right)}^{-\frac{1}{2}}$

10.根据权利要求1所述的方法，其中所述应用迭代共轭梯度方案包括：

选择迭代计数器k＝0的起始估计，如下：

$P_n^{\mathrm{dgh}\left(0\right)}=0,$ $r_q^{\left(0\right)}=p_q^{\mathrm{sct}}.$

计算迭代计数器k＝1的估计，如下：

$g_n^{\left(1\right)}={\tilde{P}}_n\underset{q}{\mathrm{\&Sigma;}}\stackrel{\&OverBar;}{K_{q,n}}{r}_q^{\left(0\right)},$ $v_n^{\left(1\right)}={\tilde{P}}_n{g}_n^{\left(1\right)},$

${\mathrm{Kv}}_n^{\left(1\right)}=\underset{n}{\mathrm{\&Sigma;}}{K}_{q,n}{v}_n^{\left(1\right)},$ ${\mathrm{\&alpha;}}^{\left(1\right)}=\frac{{\mathrm{\&Sigma;}}_n{\left|g_n^{\left(1\right)}\right|}^2}{{\mathrm{\&Sigma;}}_q{\left|{\mathrm{Kv}}_n^{\left(1\right)}\right|}^2},$

$P_n^{\mathrm{dgh}\left(1\right)}=P_n^{\mathrm{dgh}\left(0\right)}+{\mathrm{\&alpha;}}^{\left(1\right)}{v}_n^{\left(1\right)},$ $r_q^{\left(1\right)}=r_q^{\left(0\right)}-{\mathrm{\&alpha;}}^{\left(1\right)}{\mathrm{Kv}}_n^{\left(1\right)}.$

执行下面步骤直到满足误差准则为止：

通过k＝k+1而增加所述迭代计数器；以及

使用为k-1计算的先前估计而计算迭代计数器k的估计，如下：

$g_n^{\left(k\right)}={\tilde{P}}_n\underset{q}{\mathrm{\&Sigma;}}\stackrel{\&OverBar;}{K_{q,n}}{r}_q^{(k-1)},$ $v_n^{\left(k\right)}={\tilde{P}}_n{g}_n^{\left(k\right)}+\frac{{\mathrm{\&Sigma;}}_n{\left|g_n^{\left(k\right)}\right|}^2}{{\mathrm{\&Sigma;}}_n{\left|g_n^{(k-1)}\right|}^2}{v}_n^{\left(k\right)},$

${\mathrm{Kv}}_n^{\left(k\right)}=\underset{n}{\mathrm{\&Sigma;}}{K}_{q,n}{v}_n^{\left(k\right)},$ ${\mathrm{\&alpha;}}^{\left(k\right)}=\frac{{\mathrm{\&Sigma;}}_n{\left|g_n^{\left(k\right)}\right|}^2}{{\mathrm{\&Sigma;}}_q{\left|{\mathrm{Kv}}_n^{\left(k\right)}\right|}^2},$

$P_n^{\mathrm{dgh}\left(k\right)}=P_n^{\mathrm{dgh}(k-1)}+{\mathrm{\&alpha;}}^{\left(k\right)}{v}_n^{\left(k\right)},$ $r_q^{\left(k\right)}=r_q^{(k-1)}-{\mathrm{\&alpha;}}^{\left(k\right)}{\mathrm{Kv}}_n^{\left(k\right)}.$

11.根据权利要求5所述的方法，其中所述误差准则给定为：

$\frac{{\mathrm{\&Sigma;}}_q{\left|r_q^{\left(k\right)}\right|}^2}{{\mathrm{\&Sigma;}}_q{\left|p_q^{\mathrm{sct}}\right|}^2}<{10}^{-3}.$

12.根据权利要求4所述的方法，其中对所述解进行逆变换包括应用下面的等式：

$p^{\mathrm{dgh}}=\frac{\mathrm{\&Delta;}\left(\mathrm{s\&alpha;}\right)\mathrm{\&Delta;}\left(\mathrm{s\&beta;}\right)}{{4\mathrm{\&pi;}}^2}\underset{n=-\frac{1}{2}N+1}{\overset{\frac{1}{2}N}{\mathrm{\&Sigma;}}}\underset{m=-\frac{1}{2}M+1}{\overset{\frac{1}{2}M}{\mathrm{\&Sigma;}}}\exp (-\mathrm{js}{\mathrm{\&alpha;}}_n{x}_{1,q}^R-\mathrm{js}{\mathrm{\&beta;}}_m{x}_{2,r}^R)P_{n,m}^{\mathrm{dgh}}.$

Images