CN102037379B - 匹配第一和第二地震反射数据集的反射时移的方法 - Google Patents

Abstract

本发明为一种用于使第一和第二地震反射数据集(10，30)匹配的方法，所述第一和第二地震反射数据集包括具有总体对应的地震反射序列(11，31)的第一和第二地震反射道(1，3)。间隔持续数月或数年的时间获取第一和第二地震数据集(10，30)。第二地震数据集(30)包括在第一地震数据集(10)中没出现过的新地震事件(4)的至少一个横向延伸序列(40)。按需要计算反射时移(22)，以用于将第二反射道(3)的地震反射(31)与第一反射道(1)的第二地震反射日期(11)相匹配。对所述第二反射道(3)执行计算出的时移。通过计算基函数的系数来进行时移(22)的计算。

Description

匹配第一和第二地震反射数据集的反射时移的方法

技术领域

本发明涉及一种匹配第一和第二地震反射数据集的反射时移的方法。

本发明总体涉及在石油油田的寿命期间，在不同时间获取的地震数据集之间的时移校正。在一个实施例中，本发明提供一种用于在新的地震事件出现在数据中的进一步约束条件下，校正两个地震数据集之间的时移差的方法，新的地震事件可能在先前测出的地震数据集中不具有配对物。

背景技术

为了发现地质层中的石油液体，地震法是主要的方法。地震信号在表面产生并向下传播并由每个地震阻抗差(impedance contrast)部分地反射。地震阻抗为地震声速与密度的乘积。在被反射后，地震信号由一系列地震传感器获取，且对于每个来自于地震源的地震传输，地震传感器处采集的时间序列被称为地震道。为了在石油油田生产期间监测和控制地质层的流体含量的发展情况，在石油油田的寿命期间，处理所谓的时间推移(timelapse)地震数据。地质层内的物质(material)变化可引起局部地震阻抗的变化，且可看作石油生产期间在不同时间获取的地震数据之间的时移。关于地质层的物质变化的已知参数可为如何控制石油液体生产提供关键信息，诸如调整气体或石油的生产率、调整产生石油液体的深度、或者确定气体或液体的喷射(injection)率，以支持石油液体生产。

US-B-6574563描述了一种处理从相同地下区域获取的第一和第二地震数据集的非刚性(non-rigid)方法。该方法被称为“NRM方法”。NRM方法包括：将第一和第二地震数据集设置成样本集，为单独来自于一个数据集的每个样本，生成表示方向和数量的位移矢量，该一个数据集可被移动以便改善与来自于另一个样本集的对应样本的匹配。该处理通过处理样本集之一的所建议的移动来完成。该方法的优点在于，可削弱可由噪声解释的第一和第二地震时间集之间的差异，诸如由于不同源特性所引起的噪声、在所使用的拖缆中的声传感器之间的差异、源和地震拖缆的定位和深度的差异、数据获取差异、以及不同的处理。在一个基本的实施例中，该方法可被限制成仅仅表明和处理样本的垂直移动，因为这样在试图沿着地震反射体关联样本时，几乎总是可以获得良好的匹配。

所谓的NRM方法的有效工作的条件是通常对于要被比较的两个地震数据时间集的所有样本实际上存在配对物。没有显著的(significant)配对物的地震事件可使得该方法试图移动由一组样本组成的地震事件，比如说，地震数据时间集之一中的新的地震层位，以更好拟合实际上在另一地震数据时间集中未发生的地震事件，并由此使得地震数据时间集中的其它地震事件以不适当的方式移位。

先前形成的用于匹配时移地震数据的方法可被显著地改善，以便于提供地震道的甚至更好的匹配。此外，背景技术的方法几乎不允许新的地震事件，这是因为，匹配过程可将不相关的匹配强加于地震数据上，同时局部地迫使时移地震数据的其它匹配部分的移位。

发明内容

根据本发明，提供了一种用于计算时移(22)以便在石油生产活动期间检测石油承载层中的物质变化的方法，所述方法包括：

-在第一时间(t₀)处，获取第一反射道(1)的第一地震反射数据集(10)，所述第一反射道包括第一反射序列(11)，

-执行石油生产活动，

-在较晚时间(t)处，获取第二反射道(3)的第二地震反射数据集(30)，所述第二反射道包括总体对应于所述第一反射序列(11)的第二反射序列(31)，

-计算所述反射时移(22)以使所述第二道(3)的所述第二反射序列(31)与对应的所述第一反射道(1)的所述第一反射(11)匹配，

-在允许所述第一和第二反射序列(11，31)中的时不变噪声的同时，计算所述时移(22)的基函数估计。

根据本发明的优选实施例的匹配过程的主要优点如下。首先，优点在于计算时移本身的过程的改进，以更好地将在第一时间t₀处获取的第一地震反射数据集与总体时移了的在第二时间t处获取的第二地震数据集、第二地震反射数据集相匹配，其中通过计算基函数估计的系数来计算所述时移。所计算出的时移随后可被应用于其中一个地震数据集或用作用于显示变化的参数本身。与根据背景技术的执行计算时移以便于沿着每个地震道单独移动时间样本所需要的较高量的运算相比，计算时移的基函数估计可显著减少计算的投入。

其次，优点在于：执行时移的过程的改进，以在时移了的第二地震数据集包括一个或多个新地震事件时，更好地将在第一时间t₀处获取的第一地震数据集与总体时移了的第二地震数据集匹配。根据本发明的方法可包括识别新的地震事件，并将所述新事件解释为随时(或空间)变的方差和所述新地震事件以外的时不变噪声。根据本发明的此优选实施例的方法，提供了在不同时间处获取的时移了的地震反射数据集之间的改善的匹配，并允许在至少其中一个地震数据集中发生新的地震反射事件。

在本发明的有利实施例中，时移(22)的计算包括计算样条函数的系数。优点在于，样条函数中的局部变化将不会不利地全面影响震反射道。

时移的计算可有利地包括计算样条函数、勒让德多项式、泰勒级数或傅立叶级数的系数。

在本发明的优选实施例中，将应用样条函数，这是因为样条函数可允许局部新事件而不会沿着整个双程反射时间全面影响。最优选地，将应用所谓的b样条函数。

附图说明

将参照附图通过示例的方式描述本发明的实施例，其中：

图1示例了示出瞬时幅度的三个竖直设置的合成地震道。由“1”标示的左道为随机基道。由“2”标示的中间道为相同基道，但是假设有合成时移。命名为“3”的右道为具有如用于由“2”标示的道的合成时移且同时提供有新事件的基道；

图2示出了如根据背景技术所计算的用于将时移了的道(3)与图1中所示的合成基道(1)匹配的4毫秒间隔的样本的独立时移序列。横坐标为以毫秒表示的双程时间；

图3示出了正弦波形合成时移以及拟合至合成时移的对应5次勒让德多项式。横坐标为以秒表示的双程时间；

图4示例了根据本发明的方法的一个实施例，在假设有时移但没有加入新事件的合成基道上的结果，以及随后校正的结果，以及在校正前的原始差值及校正后的差值的结果；

图5示出了施加到合成地震道上的真实正弦形时移与基于勒让德多项式的合成道的估计时移的比较的笛卡尔图；

图6显示了根据本发明的方法的一个实施例的结果。从左到右示出的为：具有时移及给定的新地震事件的合成基道、合成基道、仅仅校正时移而不考虑新事件的道、原始基道和时移及新事件道之间的差值道、具有不成功地计算的基道和时移校正道之间的差值的道、以及具有理想时移和基道和时移及新事件道之间的新事件校正的差值道。后者为新事件的负像；

图7为类似于图5的笛卡尔图，其示出了，没有新地震事件的、很好匹配的、真实正弦形时移与仅仅基于勒让德多项式的估计时移的比较，并示出了对如图6中的具有新地震事件的地震道计算的附加估计时移；

图8示出了，显示时变噪声以及附加特定时间经验选择的噪声的笛卡尔图，该时变噪声包括地震道的整个双程反射时间的普通背景噪声，该附加特定时间经验选择的噪声基于双程反射时间100ms附近的新事件；

图9在左列显示了：根据本发明的方法的实施例在具有新事件的时移的合成基道上的结果；然后示出了合成基道的结果；然后是根据本发明的实施例的关于新事件校正了时移的道的结果；还示出了原始基道和时移及新事件道之间的差值道的结果；然后较成功地计算出的基道和校正了时移的道之间的差值的道的结果；而且，在右侧，示出了具有理想时移和基道与时移及新事件道之间的新事件校正差值的道的结果，即差值为就图6而言的新事件的负像；

图10示出了具有类似于图7的时移估计的笛卡尔图，同时在允许包括高幅度新事件作为噪声的时变变化形式的同时，还显示了用于长波正弦形时移的时移校正的曲线，并显示了根据本发明的方法所计算的附加的曲线；

图11示出了在North Sea的Grane石油油田，在第一时间处获取的第一地震数据集和在较晚的时间处获取的总体上时移的第二地震数据集之间的原始差值。第二地震道包括反射层，即所谓的“平点”，其为地质学解释的主要指标之一并且为本发明的主要问题；

图12和图11相同，具有由在180m到440m的偏移之间延伸的虚线标示的上述“平点”；

图13示出了被解释成如图8中所示例的时变方差的图11和12的地震剖面差值的平点序列；

图14示例了，根据本发明的利用多项式近似法-不包括使用将时移校正与归因于新事件的时变方差相匹配的方法，所计算的地震道集的估计时移；以及，

图15示例与图14相同的内容，除了根据本发明的实施例所计算出的时移与归因于总体横向延伸的新事件序列的时变方差相匹配。

具体实施方式

通过制作假设的或真实的地质柱状图的随机声阻抗日志(log)形成基反射道1。每层的声阻抗为声速乘以密度。此处所用的随机阻抗日志已经与雷克脉冲进行了卷积运算并在图1的左边第三个中由1示例。假设如此产生的随机反射道1已在第一时间t₀获取。由参考数字2标示的中间道为同样随机产生的基道，但是假设有合成时移，此处为正弦函数时间位移。假设时移的反射道2已在第二、较晚的时间t处获取。(第二时间t也可为在前的时间。)从道1、2将会看到，在由时间指示线5所标记的大约0.10秒处，对应纹波在时移的道2中发生的时间比基道1晚，并且在由另一个时间指示线6所标记的大约0.30s处，时移的道2的脉冲出现的比基道1地震图中的要早。根据US专利6574563的NRM方法能将差不多仅仅时移的道2的地震反射数据集20与基道1的第一地震反射数据集10单独地逐位相匹配。

右侧反射道3为诸如由2标示的具有合成时移的基道1，并附加地提供新事件4，新事件4在此为强的负反射，其不具有基道1中的对应特征也不具有纯粹合成时移反射曲线2中的对应特征。在“时移及新事件”道3中的这样的新地震事件4不易适用于背景技术的自动方法，这是因为，现有的自动道匹配方法会试图将这样的不匹配新地震事件4与经受匹配处理的第一道1的其它不相关的地震反射相匹配。该新事件4与基道1的不相关反射的强制匹配可将不适当的时移强加到新事件以下(或以上)的区域，可能错误地突出了道1、3的两个集之间的非物质的不真实的差异。

根据本发明的优选实施例的匹配处理的主要目的如下。首先，该方法用于实施时移以便更好的将在第一时间t₀处获取的第一地震数据集10与总体时移了的第二地震数据集20、30相匹配，第二地震反射数据集20、30在第二时间t处获取。在将计算出的时移实际应用于地震数据集之一之前，即使如此计算和显示时移，也会提供关于地震阻抗中的变化的有价值的信息。其次，该方法用于实施时移以便在时移了的第二地震数据集30包括一个或多个新地震事件4时，更好地将在第一时间t₀处获取的第一地震反射数据集10与总体时移了的第二地震数据反射集30相匹配。该第二地震数据集30在第二时间t处获取。

该方法主要意在匹配时移了的地震数据以便检测引起新地震事件4的地质层中的物质变化。假设物质变化发生在第一时间t₀和第二时间t之间的间隔中。这样的新地震事件4，即，声阻抗的显著局部变化可归因于，油贮层中的气体喷射或气/液界面的形成(或消失)、含气层压力的显著变化、在地质形成中新的油/水界面的引入、地质层中化学沉淀物的沉积、或者已稳定在地质层中的含液沉积粒子的人工喷射、或任何其它引起经受地震研究的地质柱的部分的声阻抗显著变化的物理变化。

将两个地震数据集10、30相匹配的其它目的可为比较在两个不同地震数据集10、30的现场获取期间所使用的地震数据采集设备和方法。

出现在时移数据集20和时移新事件数据集30中的差不多线性的时移可归因于多个地震数据采集参数，诸如：

-不同源特性，

-所使用的拖缆中的声传感器之间的差异，

-采样率中的差异，

-放大、声学测量的前置滤波、所采集的道的后置滤波中的差异，

-堆栈过程中的差异，即，满堆栈与近堆栈或远堆栈，

-迁移算法或参数中的差异，

-由航行错误以及拖缆漂移引起的横向或直线源或者拖缆定位差异，

-以及源和地震拖缆的深度差异。

通常，可以持续数月或数年的时间间隔获取第一和第二地震数据集。由于上述原因中的一个或多个，在具有诸如按月或年衡量的较长时间延迟的不同时间处获取的第一数据集10和第二或连续的地震数据集20、30的时移很少是恒定的。如上所述的，除了仅仅时移，也可能存在幅度差异以及波谱中的差异。期望在由时移引起的数据集间的差异和由幅度变化引起的差异之间进行区分。可能存在由石油液体生产、液体迁移或液体喷射引起的显著的新事件。为了能与已知的基数据集比较，基于合成数据提出了本发明的方法的一个实施例。提供了对来自于North Sea的Grane油田的地震数据的非商业测试，以便证明该方法对于真实地震测量的可行性。

如图1所示，基于来自于一系列性质相似的声阻抗的随机序列的并与雷克脉冲进行卷积的反射来计算合成(synthetic)基地震道1，并且在第一时间t₀处获取该合成基地震道。由参考数字2标示的中间道为同样的道，但假设有合成时移，此处为低频正弦函数时间位移，请参考图3的数字22t。假设在第二、较后的时间t处获取第二地震道2。

右侧道3根据基道1形成，具有诸如由2标示的合成时移，并附加地提供新事件4，这里新事件4为强的负反射，其不与基道1或合成地时移的曲线2相一致。

4D地震数据可由下式表示：

d_m(t)＝d_b(t+Δt)+d_4D(t)+ε           (等式I)其中，d_m(t)为作为在时间t处获取的随后的第二地震数据集的时间(或深度)的函数的道的地震数据。在等式右侧的第一项d_b(t+Δt)为作为在初始时间t₀处获取的初始第一地震数据集10的时间(或深度)加上时移Δt的时间的函数的道的地震数据，其中，时移Δt为对于第一地震道1和第二地震道2或3的每个对应测量(形成峰或谷)的小时间差。第二地震道可在第一地震道之后的数月或数年获取。在右侧的d_4D(t)项为新地震事件4，以及项ε为噪声。

上述等式I可被泰勒展开到一阶时间导数，结果如下：

$d_m\left(t\right)=d_b\left(t\right)+\left(\frac{\∂}{\∂t}{d}_b\right)\mathrm{\Δt}+(d_{4D}\left(t\right)+\mathrm{\ϵ})$ (等式II)

可重新排列该方程式以便得到所测量的地震数据和基数据之间的差：

$d_m\left(t\right)-d_b\left(t\right)=\left(\frac{\∂}{\∂t}{d}_b\right)\mathrm{\Δt}+(d_{4D}\left(t\right)+\mathrm{\ϵ})$ (等式III)

在该方程中除时移Δt外所有事情都是已知的，且最后项(d_4D(t)+ε)是已知的。如果以进行未处理的(raw)除法，则得到

$\frac{d_m\left(t\right)-d_b\left(t\right)}{\left(\frac{\∂}{\∂}{d}_b\right)}=\mathrm{\Δt}+\frac{(d_{4D}\left(t\right)+\mathrm{\ϵ})}{\left(\frac{\∂}{\∂t}{d}_b\right)}$ (等式IV)

图2中示例了该结果：除了使时移及新事件道3更好拟合至基曲线1所需要的局部时移之外，曲线中还出现了时移及新事件道3的低频正弦曲线时移(其在合成道中是已知的)。这基本上与非刚性方法相同，但是难点在于使该方案稳定。

估计时移的基函数的融合(inclusion)

根据背景技术的时移计算需要滤波并单独计算数据像素或立体像素的位移，以便更好地使第二道集20中的第二道2拟合至第一基道集10中的第一基道1相。

不是为每个位单独计算时移，而是通过计算用于将第二道2拟合至第一道1的时移函数，可节省计算投入，第二道2拟合至第一道的通过利用多项式拟合来近似时移实现。图3为多项式拟合给定曲线的能力的示例，此处为正弦形时移和表示拟合至该正弦曲线的五次勒让德多项式的对应曲线的两个几乎重叠的曲线。

Δt(t)＝c₁L₁(t)+c₂L₂(t)+...Λ

此处，c₁L₁(t)、c₂L₂(t)为第一和第二勒让德系数和函数，以及Λ为剩余误差。图3的横坐标为以秒表示的双程(two-way)反射时间。纵坐标为以毫秒表示的时移。很大程度上拟合的5次勒让德多项式与其试图表示的原始正弦形位移曲线重叠。估计的勒让德多项式仅仅由计算出的5个系数表示，而原始的正弦形位移由具有4毫秒的间隔的100个采样表示。因此，多项式函数可用于加强时移函数的平滑性，同时减小为成对的道计算满意的时移所需要的计算的数字复杂性。由此，典型地，128时移计算可由第一五次勒让德系数c₁、c₂、..、c₅的计算所代替。允许非常好的多项式拟合的简化特征实际上在于合成正弦形曲线是没有噪声的。这不是真实数据的情况。对于好的曲线近似，可以使用数种类型的基函数；勒让德多项式近似仅仅为一个好的范例。可以用于构建该曲线的目的的其它相关的基函数为泰勒级数、傅立叶级数以及样条函数。

没有新事件的合成数据的结果

无论对勒让德多项式作出了多好的近似，近似被证明了，对于如上所述的一部分纯粹的长波正弦波，找到表示道的时移的好的基函数集仍旧不理想。首先，如图4中所示例的，演示了合成基数据道的勒让德多项式拟合。第二个左侧栏是作为可以秒记录的双程时间的函数的合成基道1。左侧栏示出了图5中示例的假设有4ms正弦波形时移22t的时移道2，其中对基曲线1引起的以形成时移曲线2的真实4ms正弦形时移由灰色虚正弦曲线22t示出。用于校正时移道2以便更好地拟合基道1的5次勒让德多项式估计时移由低幅度线22ts给出。二次迭代可能会改善拟合，但发明人认为，对具有几个月或几年时间间隔的两个地震行程期间获取的两个实时数据集，获得非常好的拟合可能是不可能的。使用估计的时移22ts以此方式计算的勒让德多项式估计校正的道2c示于图4的中部。时移曲线2和基道1之间的差值由差值道2d给出。计算勒让德多项式估计时移曲线2c和原始基道1之间的差值并在图4的右侧以命名为“校正后差值”道2cd的道中绘制出。从校正道2c和基道1可以看出，匹配良好，其进一步由校正后差值道2cd的小幅度所证实。因此，为了使没有新事件出现的地震道1、2的两个时移地震集10、20相匹配，根据优选方法的通过计算多项式函数的系数来估计时移的方法被证明对这些人造道起作用。

关于新地震事件的问题

因此与背景技术的方法相比，上述改进方法将减少使两个时移地震道集20、10匹配所需要的数字投入。改进的方法也对时移道2与基道1的匹配有效，这是因为校正后差值道2cd具有小的幅度，如图4中所示。然而，如果在时移道2中引入新事件4，因此形成图6的左栏中所示的“时移及新事件”道3，则上述的方法可能试图将新事件4当作应当在基道1中的不存在的对应事件中具有其配对物的地震事件。在图6中，上述方法用于计算时移及新事件道3的校正时移以产生“校正时移但不校正新事件”道3c0。不非常成功的“校正时移但不校正新事件”道3c0与基道1之间的差值示为时移校正后初始差值道3d0，其具有显著的高于和低于新事件4的移位的负像4’的错误差异。看样子似乎，该算法试图将道3中的新事件与基道1的较低的不相关的波动(wiggle)相匹配，导致了上和下(above and below)的波动的强制移位。时移的新事件道3和基道1之间的校正差值可为图6的右手侧曲线中示出的所谓校正新事件差值3cd。如果时移道3还包括基道1中未出现的新事件4，则根据背景技术且根据上述多项式途径的方法同样地表明实施了不符合要求的匹配。进一步，在大约0到100ms的双程反射时间处，实施了如图7中所示的曲线22n所示的太大并且随后又太小的时移校正，高于图6中所示的新事件4的实际反射时间100ms，以及低于新事件4的实际反射时间大约150ms的太低的时移校正，对非匹配新事件4的适应进一步引起错误的但是较不显著的时移，如在350ms和400ms的双程时间之间的部分曲线22n的小凸起。注意与图5相比的图7中的估计时移的比例的变化。

新事件的融合

在作为测量出的地震数据d_m(t)和基数据d_b(t)之间的4D信号差方程式表述中，描写成：

$d_m\left(t\right)-d_b\left(t\right)=\left(\frac{\∂}{\∂t}{d}_b\right)\mathrm{\Δt}+(d_{4D}\left(t\right)+\mathrm{\ϵ})$ (等式III)

在该等式中，除时移Δt外的每件事都是已知的，并且最后项(d_4D(t)+ε)是已知的。以上，在假定新事件为高斯噪声的前提下，已经证明了时移及新事件道3与基道1之间的逼近(approaching)匹配方法是不成功的。知道由新事件4表示的新反射体的出现，d_4D(t)项可作为时变噪声ε_4D(t)并入如下算法中：

$d_m\left(t\right)-d_b\left(t\right)=\left(\frac{\∂}{\∂t}{d}_b\right)\mathrm{\Δt}+({\mathrm{\ϵ}}_{4D}\left(t\right)+\mathrm{\ϵ})$ (等式V)

在该等式中，噪声为时变噪声，该时变噪声可用于使得，对在已知其出现的道中的时间间隔，噪声局部地并入新地震事件4，并对于沿着道时间轴的其它时间，抑制所允许的噪声至较低的水平。图8中示例了如此精选的时间方差的范例。

噪声通常以方差的形式表示。新事件4的合理的时变方差描绘可以希尔伯特变换来表示：

σ_4d(t)＝|H(d_4D(t))|                 (等式VI)

总方差σ(t)于是以下述等式表示：

σ(t)²＝|H(d_4D(t))|²+σ_ε            (等式VII)

其中，附加项σ_ε为与时间无关的噪声的方差。如上所述，图3以及图5、7和10中的正弦形时移曲线22t不含有噪声。在根据本发明的算法中这是不可能的，所以噪声水平应被给定一个值，且被设定成如图8中所示例的新事件4的峰值的1/10。

在图9中，示例了用于计算“时移及新事件”道3的校正时移以产生“校正的时移及新事件”道3cn的上述改进方法。更成功的“校正的时移及新事件”道3cn和基道1之间的差值以时移及新事件校正后的初始差值道3dn的形式示出，其已经高度削减了超出以及低于现在更恰当地设置的新事件4的负像4cn的错误差值。看起来，根据优选方法的新算法更成功地将新事件4并入到道3中，而没有试图将新事件4与基道1的任何上或下的不相关的波动相匹配。这具有的优点是，显著减少了时移校正的波动以上和以下的错误时移位移的量。恰好在新事件4的实际时间(或实际深度，如果道被转换成深度的话)以上和以下的基道1的波动和“校正的时移及新事件”道3cn的波动现在高度对应。对于图4，时移新事件道3与基道1之间的校正差值为图9中右部所示的称为“校正新事件差值3cd的道。从而证明了根据本发明的包括时变噪声的实施例并根据多项式逼近的方法提供了也包括在基道中没出现的新事件4的时移道3的大大改进的匹配。

如在图7中所观察到的，且其在图10中被重复，在大约100ms双程时间，通过忽略新事件4的多项式逼近在实际反射时间到紧邻区域施加了太大且随后太小的时移校正，如由以划线-点-点线标示的估计时移曲线22n所示。然而，利用时移算法中的时变方差的融合，对非匹配新事件4的新的调整，提供了由估计的“时移及新事件校正”曲线22tse示出的显著改进的时移校正，该曲线以点划线指示，该曲线与实时的正弦形时移曲线22t更相似。如图9中所示的，根据本发明的该实施例的改进算法更成功地在时移校正曲线22tse中并入新事件4，而没有强制将新事件(图9中的4)与基道(图9中的1)的任何上或下的不相关波动相匹配。这具有的优点是，显著减少了上和下时移校正的波动的错误时移位移的量，显然这是估计时移校正曲线22n的问题，该估计时移校正曲线为试图将新事件4与不相关反射相匹配的结果。从而证明了根据本发明的包括时变噪声的实施例并同样根据多项式逼近的方法，提供了如曲线22tse所示例的大大改进的时移校正并允许在未出现在基道1中的新事件4的时移过程中进行融合。

真实地震数据的时变方差的估计

对于真实数据，难处在于估计所需的时变方差。在图8中以简化的方式示出一个单个道的一个单个时间事件4的该时变方差。校正地震道的时移的动机是更好地看到新事件。即，我们想从用于将包括新地震事件的时移道集与先前获取的不包括后来出现的地震事件的基道相匹配的时移校正过程得到的信息，十分类似于我们输入的信息。幸运的是，至根据本发明的实施例的算法的输入允许所怀疑的新事件时变噪声描绘的十分粗糙的估计，并仍然以显著的程度改善了结果。下面给出来自于North Sea的Grane油田的范例。

图11为包括在Grane油田在第一时间t₀处获取的第一地震数据集10的第一基道集1和地震数据30的总体时移的第二地震道集3的反射道的地震剖面(section)之间的初始差值，第二地震数据集30在第二时间t处获取。第二地震道集30包括新事件4。该数据包括由喷射的气体引起的显著时移。新事件4出现在形成道3的所谓的“平点”序列40的数个道中，道3的所谓“平点”序列出现在被认为包含所喷射的气体的地层顶上。

图12和图11相同，具有沿着上部标示的上述的“平点”、中央定位的强反射40，强反射40由在大约65ms双程反射时间处并在沿着剖面的大约180m和大约440m之间的范围内的点线标示。

图13为解释成如图8中所示的时变方差的图11和12的地震剖面差值的平点序列。所解释的平点利用35Hz的雷克小波和从初始平点差值获得的幅度进行建模，所述差值为第一地震数据集10和包含新事件4的地震数据集30之间的差值。加权函数如下定义：

$w=\frac{{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}{{\mathrm{\σ}}_{4D}{\left(t\right)}^2+{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}$ (等式VIII)

图14示例了根据本发明计算出的估计时移，估计真实地震数据基集10和随后利用多项式近似所获取的时移了的、包含新事件的、真实地震数据集30之间的时移，但是不包括将时移校正与由新事件4引起的时变方差相匹配的本发明的实施例。在大约250ms处标示出了解释成强的“新事件”4的反射层40。在200m到260m的偏移范围中，当没有考虑强的新事件4的序列40时，存在强加到随后的地震数据集30上的强力负时移，其尤其强的新事件4的序列40上但是也可在其下。在所述序列下，可假设该数据包括由所喷射的气体引起的显著时移。在大约100ms双程时间的范围处且在220m和270m之间的偏移范围中，也可在新事件的序列40以上的范围中看到强的正时移组。

图15示出了与图14相同的剖面，除了根据本发明的实施例来计算时移，本发明的该实施例包括与由新事件4引起的时变方差相匹配的方法的使用。此处，在200m到260m的相同偏移范围中，当考虑强的新事件4的序列40时，在强的新事件4的序列40以上但是也可以在其下，强加到随后的地震数据集30的时移被显著减小了。此外，在大约100ms双程时间处且在220m到270m之间的偏移范围内，出现在高于新事件的序列40中的强的正时移组明显较少地被在图15中表示出。这表明，由没有考虑新事件4的序列40出现的方法所强加的时移，对部分地震剖面实施了非实时时移校正，这是因为当考虑基于新事件序列的时变方差时，在新事件4的序列40的邻近区域中计算出的部分时移被显著减少。然而，出现在从大约310m到大约340m的偏移范围中的在大约300ms与大约360ms之间的时间段中的强的负时移的黑的、阴影状的组在图14至图15中占主要地位并被看作物质源。此外，在从大约400ms到大约460或者甚至接近500ms的较深时间剖面中，对于在210m与380m之间偏移的估计时移在两种方法之间具有较小的时移差。这进一步表明那些计算出的时移是真实的。穿过大约440ms双程反射时间的3-D“时间-水平”剖面可指示高至大约190m的偏移范围中的较小的估计时移值，随后为高至大约380或甚至420m的较长范围中的较高时移值，接着是超出600m偏移的总体较低的时移。图像总体与在图14和图15中的此范围的相同。在本例中，这可指示在标示深度和偏移范围中的气体的引入。

包括第一和第二地震数据集10、30的4-D地震集包含时移和幅度变化。本方法的一个目的是将时移对地震数据集10、30的影响与幅度变化对相同的地震数据集的影响分开，以便于研究这样的时移，并便于研究这样的地震幅度变化。如果可基于良好确立的幅度变换地震数据，可提供所喷射的气体量的改进的评估。

因此，利用需要比背景技术的方法少的计算投入的基函数计算，优选实施例可提供时移地震数据集之间的好的匹配，其中背景技术的方法对每个时间样本计算单独的时移。在新地震事件4的横向序列40的约束下，本发明的实施例还能够在在第二地震数据集中出现新事件4的条件下计算时移，而不会将错误的时移强加到地震数据中。

通过知道显著引发的时移及它们的深度，可计算出地震阻抗的变化。地震阻抗中的变化可以归因于所讨论的地质层的地震速度或密度中的变化，或者归因于该两者的变化。知道这样的关于地质层中的物质变化的参数可提供关键信息输入以便控制石油液体的生产，诸如调整气或石油的生产率、调整产生石油液体的深度或确定气体或液体的喷射率，以便维持石油液体的生产。

特别参照示例的范例描述了本发明的实施例。然而，应当理解，可在本发明的范围内对所描述的范例进行改变和更改。

附录：

                根据时间推移地震的时移估计

                     Espen Oen Lie

                     2008年2月25日

1 数学模型

利用十分简单的假设导出用于估计时移的模型。

d_m(t)＝d_b(t+Δt)+d_4D(t)+ε         (1)

对该等式以Δt进行麦克劳林级数展开

$d_m\left(t\right)=d_b\left(t\right)+\left({\∂}_t{d}_b\right)\mathrm{\Δt}+O\left({\mathrm{\Δt}}^2\right)+d_{4D}\left(t\right)+\mathrm{\ϵ}---\left(2\right)$

对其进行重新排列并忽略二次(以及高次)项

$d_m\left(t\right)-d_b\left(t\right)=\left({\∂}_t{d}_b\right)\mathrm{\Δt}+d_{4D}\left(t\right)+\mathrm{\ϵ}---\left(3\right)$

一些注释：

·如果去掉包含Δt的项，我们得到时间推移差立方的标准解释。即，那个差值仅仅与幅度效应和噪声有关

·如果去掉幅度效应(d_4D(t))，我们得到与能够用于导出在NRM中使用的Horn-Schunk法相同的等式

我们将暂时省略d_4D(t)项以导出该方法的基本原理。稍后将包括幅度项。我们能够写下简洁版的等式

Δd＝DΔt+ε                 (4)

其中，Δd为d_m(t)-d_b(t)，以及D为对角为的对角矩阵。该等式可直接逆转换成Δt＝D^-1Δd，但是那将在数学上和物理上都是不合理的。该结果将会非常嘈杂和非物理的。从物理学观点，时移应该在纵向和横向都十分平滑。也许除了误差周围。我们认为速度上的变化是不连续的，但是得到的时移仍旧会纵向连续。

该方法中适合的平滑约束是基于基函数的。代替直接以Δt计算，我们声称该函数可由一组基函数充分描绘。能够使用任意组的函数，但是目前b样条函数是优选的。基函数的引入能够表示为

Δt＝Fc                 (5)

F为包含n_f个基函数的的n_t×n_f矩阵。n_f典型地大大小于n_t。将其代入，得到

Δd＝DFc+ε             (6)

注意，由于D为对角矩阵，DF不需要作矩阵乘法。我们宁愿构建包含所有基函数次数的矩阵。该等式可容易地以最小二乘的方式求逆，且对于只包含时移的合成结果非常稳定。解为

c＝(F^TDF)^-1F^TΔd          (7)

Δt＝F(F^TDF)^-1F^TΔd       (8)

这可被称为问题的1D解。

1.1 横向约束

对于合成数据，1D解是足够的，但是真实地震包含许多不能用这些简单等式建模的效应。对于这些效应的通常途径是将所有其他事情看作噪声。此处将采用同样的途径，且假设时移比该“误差”平滑。可能的解决方案是引入2D基函数。然而，这将施加非常强的约束，其将忽视误差和地震中的其它较不平滑的特征。

我们将假设，从道到道的时移变化是高斯的。

该假设总体上是实际的，赞成理由如下：

·约束和求解是全局的(没有滑动窗)

·求解非常快

·横向时移变化非常灵活

最重要的反对理由在于假设从某种角度是靠不住的。我们没对问题进行统计分析，但是仍旧利用统计的方法对它求解。

通过施加统计假设，通常以贝叶斯定理求解。为了完整性，包括方程

p(c|d)∝p(c)p(d|c)           (9)

其中，d为测量的数据，以及c为期望的参数。在前的高斯函数可被表示成

$p\left(c_{j,k}\right)\∝\exp \left[{-\left(F(c_{j,k}-c_{j,k}^0)\right)}^T{C}_c^{-1}\left(F(c_{j,k}-c_{j,k}^0)\right)\right]---\left(10\right)$

$c_{j,k}^0=\frac{1}{4}(c_{j-1,k}+c_{j+1,k}+c_{j,k-1}+c_{j,k+1})---\left(11\right)$

其中，j和k与空间采样相关。注意到，包括以F表示的基函数。这是因为它是时移而不是被比较的基函数系数。C_c为先前的协方差矩阵。通常，似然性(likelihood)估计为

$p\left(d_{j,k}\right|c_{j,k})\∝\exp \left[{-(D_{j,k}{\mathrm{Fc}}_{j,k}-{\mathrm{\Δd}}_{j,k})}^T{C}_{\mathrm{\ϵ}}^{-1}(D_{j,k}{\mathrm{Fc}}_{j,k}-{\mathrm{\Δd}}_{j,k})\right]---\left(12\right)$

其中，C_ε为协方差矩阵。注意到，似然性是局部的。该问题的最大或然性(probability)解为

$c_{j,k}={[{\left(D_{j,k}F\right)}^T{C}_{\mathrm{\ϵ}}^{-1}{D}_{j,k}F+F^T{C}_c^{-1}F]}^{-1}---\left(13\right)$

$[{\left(D_{j,k}F\right)}^T{C}_{\mathrm{\ϵ}}^{-1}{\mathrm{\Δd}}_{j,k}+F^T{C}_c^{-1}{\mathrm{Fc}}_{j,k}^0]---\left(14\right)$

$={[F^T{D}_{j,k}^T{C}_{\mathrm{\ϵ}}^{-1}{D}_{j,k}F+F^T{C}_c^{-1}F]}^{-1}---\left(15\right)$

$[F^T{D}_{j,k}^T{C}_{\mathrm{\ϵ}}^{-1}{\mathrm{\Δd}}_{j,k}+F^T{C}_c^{-1}{\mathrm{Fc}}_{j,k}^0]---\left(16\right)$

$\≡{(L_{j,k}+\mathrm{\Γ})}^{-1}(b_{j,k}-{\mathrm{\Γc}}_{j,k}^0)---\left(17\right)$

然而，该解是非局部的，这是因为是非局部的。如此处所述的，不可能求解该问题，但是它是寻求求解该问题的启发方法。(实际上，它是雅可比迭代法方案，为了我们的目的要慢慢来)

1.2 作为线性方程系求解

存在多种求解该方程系的方法，但是我们应该指出问题的大小。我们具有n_f个基函数，n_i为受测线(inline)以及n_x为横测线(crossline)。这形成矩阵n_fn_in_x×n_fn_in_x矩阵，其典型地为20000000×20000000。所以我们需要利用非常稀疏(以及带状)的矩阵结构。我们最终利用共轭梯度法求解，但是首先我们将为此目的重新排列方程。

$(L_{j,k}+\mathrm{\Γ})c_{j,k}-\frac{1}{4}\mathrm{\Γ}(c_{j-1,k}+c_{j+1,k}+c_{j,k-1}+c_{j,k+1})=b_{j,k}---\left(18\right)$

这导致了带状系，其可被写成

$c_j=\left(\begin{array}{c}{c}_{j,1}\\ c_{j,2}\\ c_{j,3}\\ \·\\ \·\\ \·\end{array}\right),$ $b_j=\left(\begin{array}{c}{b}_{j,1}\\ b_{j,2}\\ b_{j,3}\\ \·\\ \·\\ \·\end{array}\right)---\left(20\right)$

其中，D_Γ是对角为Γ的块对角矩阵。

为了数值的目的，Γ可从对角消除：

${\hat{b}}_j=\left(\begin{array}{c}{\mathrm{\Γ}}^{-1}{b}_{j,1}\\ {\mathrm{\Γ}}^{-1}{b}_{j,2}\\ {\mathrm{\Γ}}^{-1}{b}_{j,3}\\ \·\\ \·\\ \·\end{array}\right)---\left(23\right)$

1.3 指定协方差-先验知识的融合

除了基函数(其必须被指定)之外，两个协方差矩阵加入方程。这些必须被指定，且具有各自的用途。最简单的途径为，除需要平滑的结果外，我们不知道任何其它事情。在这种情形下，协方差是对角的且为常量。于是仅仅加入了σ_ε/σ_c，其可被解释成噪声信号比。或者为加权函数。

到目前为止，我们已经忽视了d_4d，其为所期望的幅度效应。所有其它的方法将此当作噪声，我们也将其看作噪声，但是具有空间变化的方差。C_ε为协方差，其到目前为止被当作常量。我们将假设零时间相关，以使得C_ε为对角的，但是具有时间和空间变化的方差所建议的这个的实际实施在于，在时移估计之前，将4D幅度效应建模。这个通过利用所解释的层和具有等于地震的频率的雷克小波来进行。幅度可从初始的差立方来获取。注意到，基函数内插于建模区域上，以使得主要效应为这很重(它不会直接影响时移)。由于方差等效于能量，我们需要将建模的4D效应转变为能量。通过利用希耳伯特变换，可同时移除相位，使方差十分平滑并确保能量一致(在希耳伯特变换之后，信号的总能量是恒定的)。与时间有关的方差于是为

σ_4D(t)＝|H(d_4D(t))|              (24)

其中，H(·)为希耳伯特变换。其给出了总方差

σ_ε(t)²＝σ_4D(t)²+σ²             (25)

其中，σ²为经典意义上的噪声方差。这将为对角的的对角。

第二协方差C_c也能够用于包括现有知识。这与时移横向上多么相似有关。当位于均一区划的内部时，认为时移在横向上是相似的。即，位于地质单元内部。时移不必连续地超过断层(fault)和其它地质学边界。通过C_c可包括这种情况。注意到，这样的融合很可能导致时移与输入约束相一致，因此或许不采用地质模型。这起因于它们的边界受到一些臆测。断层属性或其它数据驱动的测量更为合适。再一次，方差应该从能量方面来说。

注意到C_ε和C_c应被归一化(normalized)，以使得C_c＝σ_cI_c，其中I_c在对角上大于1。这导致在对角上小于1，且减小或保持F的数量。这应该全局进行，以使得输入为I_ε，I_c和这也意味着将结构输入信息与平滑度分离。除了可行性方面，这也确保了数值的稳定性。

2.数值求解

最好的检测技术是共轭梯度法(CG)。与最初检测的高斯—赛德尔法相比，这提供了极大的加速。CGs的优势之一在于，它仅仅需要矩阵向量积而不需要矩阵求逆。归因于我们系统的结构，矩阵向量积非常容易得到：

q＝Ap           (26)

$q_{j,k}={\mathrm{\Γ}}^{-1}{L}_{j,k}{p}_{j,k}-\frac{1}{4}(p_{j+1,k}+p_{j-1,k}+p_{j,k+1}{p}_{j,k-1})---\left(27\right)$

无疑我们不在每次迭代进行Γ^-1L_j，k，这仅仅是使得符号简单。

实际的算法如下：

对一些初值估计x⁰计算r⁰＝b-Ax⁰

for i＝1:maxit

${\mathrm{\ρ}}^{i-1}=r^{{i-1}^T}{r}^{i-1}$

     if i＝1

        p¹＝r⁰

  else

     β^i-1＝ρ^i-1/ρ^i-2

     pⁱ＝r^i-1+β^i-1p^i-1

  endif

    qⁱ＝Apⁱ

${\mathrm{\α}}^i={\mathrm{\ρ}}^{i-1}/{p^i}^T{q}^i$

    xⁱ＝x^i-1+βⁱp(i)

    rⁱ＝r^i-1-βⁱq(i)

  check convergence；continue if necessary

end。

Claims (26)

1.一种用于计算反射时移（22）以便在石油生产活动期间检测石油承载层中的物质变化的方法，所述方法包括：

-提供在第一时间（t₀）处获取的第一反射道（1）的第一地震反射数据集（10），所述第一反射道包括第一反射序列（11），

-提供在较晚时间（t）处获取的第二反射道（3）的第二地震反射数据集（30），所述第二反射道包括总体对应于所述第一反射序列（11）的第二反射序列（31），

-计算所述反射时移（22）以使所述第二反射道（3）的所述第二反射序列（31）与对应的所述第一反射道（1）的所述第一反射序列（11）匹配，

其特征在于，

-在允许所述第一反射序列和所述第二反射序列（11，31）中的时不变噪声的同时，计算所述反射时移（22）的基函数估计。

2.根据权利要求1所述的方法，其中，所述反射时移（22）的计算包括计算样条函数的系数、勒让德多项式的系数、泰勒级数的系数或傅立叶级数的系数。

3.根据权利要求1或2所述的方法，还包括：对所述第二反射道（3）执行计算的所述反射时移（22），以与所述第一反射道（1）相匹配。

4.根据权利要求1或2所述的方法，其中，所述第二地震反射数据集（30）的所述第二反射序列（31）包括在所述第一地震反射数据集（10）中不存在的新地震事件（4）的至少一个横向延伸的序列（40），并且还包括：

-将所述新地震事件（4）解释为时变方差和所述新地震事件（4）以外的所述时不变噪声，并计算所述基函数估计，同时容许所述新地震事件（4）附近的所述时变方差和远离所述新地震事件（4）的所述时不变噪声。

5.根据权利要求1或2所述的方法，包括识别不存在于所述第一反射道中的至少一个横向延伸的系列的新地震事件并将所述新地震事件解释为时变方差和所述新地震事件以外的时不变噪声。

6.根据权利要求5所述的方法，其中，所述新地震事件（4）的识别包括：计算从所述第二地震反射数据集（30）减去所述第一地震反射数据集（10）的差的幅度，用取决于所述时不变噪声及所述时变方差的权重函数乘以所述幅度，并与地震脉冲函数进行卷积运算。

7.根据权利要求6所述的方法，其中，所述地震脉冲函数为雷克脉冲。

8.根据权利要求1或2所述的方法，其中，所述第一时间（t₀）和所述较晚时间（t）之间的间隔是两个月以上。

9.根据权利要求1或2所述的方法，包括沿地震剖面线显示作为反射时间函数的所述第二反射道（3）的所述第二反射系列（31）的所计算的所述反射时移（22）。

10.根据权利要求1或2所述的方法，包括从地球的地下区获取所述第一地震反射数据集（10）和所述第二地震反射数据集（30）。

11.根据权利要求1或2所述的方法，还包括执行所述石油生产活动。

12.根据权利要求1所述的方法，其中，计算所述反射时移的步骤包括比较所述第一反射道和所述第二反射道，以确定所述第一反射道和所述第二反射道中的事件之间的反射时移。

13.根据权利要求12所述的方法，其中，比较步骤包括形成每个第一反射道和对应的第二反射道之间的差并以所述第二反射道对时间的一阶导数除所述差。

14.根据权利要求12或13所述的方法，其中，计算基函数估计的步骤包括：通过计算基函数的系数来将所述基函数与所述反射时移匹配以实现匹配，以及将所述基函数应用于所述第一反射道和所述第二反射道中的至少一个以提供时间匹配的第一反射道和第二反射道。

15.根据权利要求1或2所述的方法，还包括使用所计算的反射时移来检测所述石油承载层中的物质变化。

16.根据权利要求1或2所述的方法，还包括在石油生产活动期间使用所计算的所述反射时移。

17.根据权利要求16所述的方法，包括使用所计算的所述反射时移来监测或控制石油液体生产。

18.根据权利要求17所述的方法，其中所述监测或控制石油液体生产是从包括下列操作的组中选择的以支持所述石油液体生产：调整气体或石油的生产率、调整生产石油液体的深度或确定气体或液体的喷射率。

19.一种确定或检测地球的地下区中的地球物理性质的变化的方法，所述方法包括执行根据权利要求14所述的方法，其中，所述第一地震反射数据集和所述第二地震反射数据集（10，30）是从所述地下区获取的，并且所述方法还包括：根据所述时间匹配的第一反射道和第二反射道确定所述第一时间和所述较晚时间之间所述地下区的所述地球物理性质的变化。

20.根据权利要求19所述的方法，其中所述地球的地下区中的地球物理性质的所述变化是地质层中的物质变化。

21.根据权利要求19或20所述的方法，还包括在石油生产活动期间，使用所确定的或所检测的所述地球物理性质的变化。

22.根据权利要求21所述的方法，包括使用所确定的或所检测的所述地球物理性质的变化来监测或控制石油液体生产。

23.根据权利要求22所述的方法，其中所述监测或控制石油液体生产是从包括下列操作的组中选择的以支持所述石油液体生产：调整气体或石油的生产率、调整生产石油液体的深度或确定气体或液体的喷射率。

24.一种支持石油生产活动的方法，其中，在石油生产活动期间，使用根据权利要求1至15中的任一项所述的方法计算的所述反射时移或根据权利要求19或20所述的方法确定的所述地球物理性质的所述变化。

25.根据权利要求24所述的方法，包括使用所述反射时移或所述地球物理性质的所述变化来监测或控制石油液体生产。

26.根据权利要求25所述的方法，其中所述监测或控制石油液体生产是从包括下列操作的组中选择的以支持所述石油液体生产：调整气体或石油的生产率、调整生产石油液体的深度或确定气体或液体的喷射率。