CN1007555B - 由熵引导的地震信号的消褶积 - Google Patents

由熵引导的地震信号的消褶积

Info

Publication number
CN1007555B
CN1007555B CN85108418A CN85108418A CN1007555B CN 1007555 B CN1007555 B CN 1007555B CN 85108418 A CN85108418 A CN 85108418A CN 85108418 A CN85108418 A CN 85108418A CN 1007555 B CN1007555 B CN 1007555B
Authority
CN
China
Prior art keywords
signal
operator
observed
observed signal
seismic
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired
Application number
CN85108418A
Other languages
English (en)
Other versions
CN85108418A (zh
Inventor
拉尔夫·威金斯
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Schlumberger Overseas SA
Original Assignee
Schlumberger Overseas SA
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Schlumberger Overseas SA filed Critical Schlumberger Overseas SA
Publication of CN85108418A publication Critical patent/CN85108418A/zh
Publication of CN1007555B publication Critical patent/CN1007555B/zh
Expired legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364Seismic filtering
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/50Corrections or adjustments related to wave propagation
    • G01V2210/56De-ghosting; Reverberation compensation

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Geology (AREA)
  • Environmental & Geological Engineering (AREA)
  • Acoustics & Sound (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Electrophonic Musical Instruments (AREA)
  • Control Of Vending Devices And Auxiliary Devices For Vending Devices (AREA)

Abstract

本发明是一个处理过程和系统,用它来改善地震回波和其它信号,这些信号代表由于未知子波的相互作用而被降低质量的那些非时间变化系列。使用二项式算子序列,其中每个算子都从最新的信号中导出,并且每个算子都被应用于该信号以获得下一个更新的信号。这些算子是这样选择的,使得它们的应用趋于减小(原文有误—译者)信号的熵,以及去掉未知子波的作用,并使被观察信号变得更能揭示感兴趣的信号。

Description

本发明涉及对石油,天然气和其它地下资源的地震勘探,并涉及使用消褶积技术改善地下岩层的地震测绘,它甚至能在被观察的地震信号不是最小相位或离高斯分布不很远的情况下给出良好的结果。从广义上讲,对由一个未知子波与一个未知静态随机序列的物理相互作用(褶积)产生的信号进行改善,从而使子波的平滑效应可以被除掉或者至少基本被除掉,使随机序列可以被揭示出来,这只是本发明所涉及的一个方面。在这个意义上说,可用于地震勘探的爆炸所产生的子波就是未知子波的一个例子,地下地震反射体的位置则是静态随机序列成员的例子,而要改善的信号就是地震记录迹,这个记录迹是在对地震信号的地面测量的基础上得到的,其中,地震信号是由爆炸引起的子波与各个地下反射体(位于具有不同声学性质的地下岩层的界面上)的相互作用(褶积)产生的。本发明的其它方面将在下面详细叙述。
地下反射地震学是根据从地下岩层反射回来的信号的到达时间来测绘地下岩层的。在地表产生的地震能量穿过地球的分层介质,其中一部分能量在各岩层之间的界面上被反射回来。反射回地表的能量由地面上的接收机(探测器)测量。在陆地上,产生地震能量的典型方法是使用化学炸药、振动机械或撞击装置。预先放置的接收机也称做探测器或地音探测器,它们列成一排配置,用来探测反射回来的地震能量。在海上,当舰船沿预定路线走过时,每隔几秒种就启动诸如用一排喷枪做的地震源。返回的地震能量用船后边拖着的装在电缆(一个浮筒)里的多个探测器采集。
地下构造的复杂性,被观察信号的噪声及其它不足,使得人们必须通过一系列的操作来处理这些信号,才能得出一个有用的地下岩层的分布图,而上述每一步操作又包括严格的近似和简化。在这种信号处理中,一个主要操作是速度分析,它所以必要是因为:在地面观察到的地震信号表示成作为时间函数的反射信号,它必须被转换成深度的函数,以便导出更有用的地下构造的分布图。然而,地震波的速度对传播介质的依赖性很大,因而地震波在向地下传播过程中其速度变化很大。总的看来,速度随深度增加而提高,当然也有可能偶而遇到某些地层使其速度减小。就一个给定的地表点来说,作为深度的函数画出的速度曲线通常称做速度函数。因此,有两个重要的变量:反射信号的时间和速度。知道了这些变量,就可以确定反射层位的深度。因为速度有很大的横向变化,而且速度函数从一个地点到另一个地点也有变化,所以对整个勘探只用一个给定的速度函数那是不合适的,而是必须在勘探区域中从一处到另一处不断地对其进行修正。当在正确的地方能获得足够深的钻孔并对其记录试验结果时,那么可以从置于钻孔不同深度的地震探测器的输出信号中估计速度函数,或者通过垂直地震断面图(“VSP”)来获得速度函数。然而在大多数情况下,必须通过局限于地面的测量和对这些测量结果进行复杂的信号处理来估计速度函数。速度分析结果还用在下面马上要讲到的动态校正中。
对被观察信号处理中的另一个主要操作是对接收机的输出进行各种校正。所谓静态校正包括震源校正和探测器校正,它们的组合效应是把地震信号源和地震信号探测器置于假设的水平面上。实际上,这一点是通过对多个接收机输出进行适当的时移来实现的。动态校正把每个接收机的输出转换成这样的输出,它是假定震源和接收机处于同一横向点(即实际震源和接收机位置间的中点)时被接收的。在这种转换中,各轨迹均被认为是由主反射组成的,所谓主反射是由从震源向下到反射地 层,然后直接反射回接收机的路径构成的。如果各地层是水平的,那么所有的反射点(或称深度点)正如位于震源和接收机之间中点的下边。如果各层是倾斜的,则各深度点偏离这个中心点。因此动态校正既取决于速度函数,又取决于反射岩床的倾斜情况。由震源和接收机的分离所致的校正分量称为正常的外移校正,而由于倾斜所致的分量称为倾斜校正。因此,对记录的轨迹(探测器输出)的四个重要校正是:震源校正,接收机校正,正常外移校正和倾斜校正。
每个被观察信号(记录轨迹)可以看成由反射波与各种干扰波和噪声构成的时间系列。人们需要的反射波是主反射,它是那些直接向下传到某反射体,然后直接返回地面并在地面加以记录的波。人们不希望的一种主要的干扰波是多次反射,它所经历的路径是:向下奔向某一反射体,然后向上奔向另一个反射体,它再向下到一个反射体,再向上,或许再上下数次,直到多次反射的信号到达地表面被记录为止。在任何分层系统中都可以有无数个多次反射。这种多次反射的存在使得对主反射的识别变得困难。因此,必须尽可能把多次反射衰减掉。
为了便于取出有用的信号,通过布置一个震源和一排探测器,开动震源及记录探测器的输出信号,然后沿探测器行列横向移动整个布置,并且重复这个过程,来记录典型多倍的有效区域。当这种布置在爆破(shots)之间以足够小的增量移动时,每个深度点可被复盖若干次,从而导致相当大的信息重复度,它能使主反射波随之增强,而使多次反射受到抑制。如果将勘探中的所有轨迹进行分离,归类成组,并把每组叫作汇集,使得某一给定汇集中的所有轨迹具有共同的震源与接收机间的中点,则适当的外移校正把给定汇集中的每个轨迹转换成大致相同的等效轨迹,即假如震源和接收机正好位于上述共同中点时所接收的主反射轨迹。换句话说,在适当的外移校正之后,所有轨迹的主反射将趋于同相,而多次反射将趋于因轨迹而异,因此一般将异相。如果这种方法先 用于正常的外移校正,那么它也可以借助对不同的轨迹汇集的利用来进行适当的震源,接收机和倾斜校正。另一方面,同时校正方法,既可以是直接的,也可以通过采用各种迭代技术加以利用。如果把具有共同中点的汇集内的校正后的各种轨迹相加,则结果是多次反射波被大大衰减,因为它们之间相位不一致。而主反射波将得到放大,因为这种信号在适当选择多个汇集和校正后互相之间将同相。因此,对每个中点来说,均可以产生一个输出轨迹,这个轨迹叫作这个中点的层迭轨迹。
多次反射中有一类特殊情况叫做混响。一种看法认为混响是指这样的能量,它聚集在靠近地表的多个岩层中,并不断被上下反射。当这种能量通过靠近地表岩层时,可能附着在主反射上。结果,被观察到的反射不再那么明显清楚,而且时间分辩下降,它后边跟着长长的混响波列。这些波列再被后续的反射的波列重迭,因此整个地震轨迹(被观察信号)可能具有阻尼振荡或正弦曲线的特征,从而难以或不能分辩出主反射的起始点。消除每个混响波列能量,并保持主反射原封不动,可以提高所有反射信号的分辩率。这个消混响波列的过程一般被称为消褶积,可以把它看成地震源子波与地下岩层的褶积(即卷在一起或相互作用)的地球物理学过程的反演。然而,正如下面将要更清楚叙述的,这不是消褶积的数学过程。
消褶积的一个已知类型是基于这样的考虑:被观察信号的子波被看成是由主反射和附在它上面的混响组成的。根据这个假设,从它的能量大部分集中在其起始部分这个意义上讲,该子波的延迟应该最小。此外,从整个波形来看所有这些子波应该是相似的。由于主反射来自据信是厚度无规律分布的地质层,所以主反射波的到达时间实际上是随机分布的。因此轨迹的自相关函数将与子波的自相关函数大致相同,而且根据自相关函数,人们就可以估计要求的反演消褶积算子。把该算子运用于轨迹将得出消褶积轨迹,即得出其中子波的混响分量被除掉的轨迹, 因此提高了主反射的分辩率。消褶积处理也可以用来消除某些长周期的多次反射。实际上,消褶积用于地震轨迹既可在层迭之前也可在它之后,这取决经费和其它一些考虑。
消褶积之后,被处理的信号一般还要经过所谓转移或反传播这样的技术处理,这个技术可以被概括为这样一个过程:把在地表面观察到的波及时反向传播到地下的地层结构,以揭示地下各反射点在深度方面的实际空间位置。这个过程也可描述成把在地表观察到的信号转换成在地下深度上应能观察到的信号。
正如本发明人在一篇文章中所讨论的那样(R.A.威金斯(Wigyns)),“最小熵消褶积”,地质勘探,16(1978),21-25),如果震源子波至少是近似知道的,那么借助估计震源子波的最小二乘方或借助维纳(Wiener)变换对子波整形可直接得到一个良好的希望信号估值。如果震源子波是最小相位或最小延迟,而且希望信号几乎是全频谱的,那么可以使用预测消褶积。如果震源子波的振幅谱相对于希望信号的振幅谱来说是平滑的,则可以使用同态消褶积。当希望信号由几个大的尖脉冲组成,而且有这样的采样信号整体存在,对这个整体来说震源子波形状保持不变,而各尖脉冲间的时间间隔不同时,则这时可能按照该文章中详细叙述的方式使用最小熵消褶积。在威金斯的文章的第34和35页所引用的参考文献中还能找到进一步的背景信息。从1969年艾尔斯威尔(Elsevier)科学出版公司出版的E.A.鲁宾逊(Robinson)等人所写的“在石油和天然气勘探中的地球物理时间序列的消褶积”一书以及1981年古斯庞德(Gooes    Pond)出版社出版的E.A.鲁宾逊所写的“时序分析和应用”一书中还能找到更多的背景信息。威金斯文章与上述鲁宾逊等人所写的书中的内容在这个说明书中也被引用了。
人们确信,已知的用于地震轨迹的多种“消褶积”处理有许多缺点。例如,这些处理方法据信对不具有最小相位(延迟)或离高斯分布不甚远 的地震轨迹不能很好地进行处理,不同的轨迹需要选用不同的处理技术。此外,据说已有的各种处理术太复杂而且精度不高。
因此,本发明的一个重要目的是对子波和非高斯随机过程的褶积结果(诸如地震源子波与地下反射体的相互作用结果)提供一个改进的消褶积技术。其处理方式能对所有的或大多数典型地震轨迹很好地工作,不管轨迹的类型如何均使用同一方法,它处理轨迹的时间仅限于充分地对它消褶积所需要的时间,而且操作快速有效。从广义上说,一个重要目的是要找到一种实用的办法解决所谓盲目消褶积问题,这个问题是在为了某一实际目的,把观察和测量的信号看成是一个未知子波与一个静态随机时序的褶积(相互作用)时遇到。目的是将子波的影响消除到这样的程度,即允许在揭示和研究基本随机时序时不受或基本不受子波的平滑效应的影响。本发明特别适用于以实用方式解决地震勘探中长期存在的问题,但它也适用于其中观察到的信号不是有价值的而是子波与感兴趣的时序的褶积这样一些其它领域,诸如毫无限制地应用在某些声纳和雷达活动及在使用声音信号和其它信号和其它一些领域中,这其中的信号受了与子波相互作用而变得不纯净,而这种相互作用就实际目的来说可以处理成褶积。
用概念术语来说,正象它被应用到地震勘探时那样,该改进的消褶积技术包括:找到可以被表征成反演滤波器的东西,当该滤波器应用于被观察信号时,它能把该信号变成深度点的反射系数和位置(按时间)的良好估值,而这些深度点则产生引起被观察信号的主反射。在本发明的一个无限制的例子的实施过程中,这一点是采用二项算子对被观察的信号进行连续修正来完成的。第一个算子根据被观察轨迹及对预期反射系数系列的熵的估计加以选择,并把它运用到被观察的信号以产生一个修正的观察轨迹。第二个算子是基于该被修正的被观察轨迹及类似的熵方面的考虑来选择的,将它应用于已修正的轨迹产生一个新的修正轨迹。 第三个算子根据上一个修正的被观察信号和类似的熵考虑加以选择,然后应用于上一个修正后的被观察信号以产生另一个修正的被观察信号,以此类推直到没有进一步的明显改进为止。从而最后修正的被观察轨迹就是所关心的地下反射体的反射系数和位置(按时间)的良好估值。每个二项算子可以是〔1,0,0,……a〕形式或者〔a,0,0,……1〕的形式,其中1项和a项间的间隔是根据最新修正(既使有也极少)的被观察信号的自动相关的最大偏差的滞后来决定的,而a项的振幅是根据最新修正的被观察信号的其它特征(如它的自动相关)求出的。在该处理过程中的给定运算中使用各种可能算子中的哪一个,要根据哪个算子对修正的被观察信号给出的熵较小来决定,这可以用诸如Varimax值来估计。
当把本发明应用到地震勘探时,实施本发明的具体无限制例子的主要步骤如下:
1.产生一个被观察信号,如果把一个层迭的轨迹用于被观察信号,在第一次通过步骤1时,该被观察信号是层迭信号本身,而在以后每次通过步骤1时,则被观察信号是下面通过步骤5的最新结果。被观察信号可以是数字信号y(i),其中i=1,2,……I,而标号i表示在一个周期内进行的数字采样,这个周期相当长,足以包括所关心的最深反射体(一般是几秒钟)而采样频率足够的高,足以满足所关心的最高被观察信号频率所要求的尼奎斯特(Nyguest)判据。虽然为了简化只对使用层迭轨迹作为被观察信号的例子进行了详细讨论,但不是说非这样不可。作为无限制的例子,被观察信号可以是一个接收机的输出,层迭轨迹以外的各接收机输出间的某些组合,也可以是地面反射地震学中的被观察信号或者在其它类型的地震勘探中的被观察信号,诸如垂直地震剖面图(VSP),或者还可以是某些其它被测量信号,诸如声纳或雷达的信号;
2.求出i=0,1,2,……I时的被观察信号y(i)的自相关函数R(yy,
3.从步骤2求出的自    相关函数R中求出自    相关的被选定旁瓣滞后d(例如以步骤i表示的R(y,y,0)处的自    相关的中心瓣和在R(y,y,i)处的被选定旁瓣间的距离)。这个滞后d决定了该算子上述两项之间的零的个数(如果d=i,则有(d-1)个零)。
4.从步骤2中求出的同一自 相关R(y,y)求出a项的振幅。对求出a项振幅有用的关系的例子之一是a=〔-R(yy;d)〕/R(yy;0)〕。R(yy;d)是在被选择旁瓣上现有被观察轨迹的自 相关函数的振幅,而R(yy,0)是中心瓣的自 相关振幅。另外,也可以转而使用自 相关特征之间的其它关系,其条件是当a项为一个二项算子的第一项或最后一项用来对现有被观察信号运算时,最终的a项能把结果的熵减到最小(或者至少是大大减小),这可通过诸如核查Varimax量V是否是最大值(即核查新的V值是否比前一个V值加上阈值之和大,这个过程仅当回答是肯定时继续进行)来确定。确定了间隔d(因此1项和a项之间的零的个数也就确定了)和a项的振幅,则把两个可能的二项算子〔1,0,0,……a〕和〔a,0,0,……1〕的其中一个应用于当前被观察信号,可根据它们对被观察信号的熵的减小的效果大小来对这两者加以选择。例如,这一点可根据V=y(i)4/〔y(i)22这样的Varimax量V来实现,上式中y(i)是把被检查的二项算子应用到现有被观察信号而产生的修正的被观察信号。给出较大V值的二项算子被选用,因为它能更有效地减小被观察信号的熵;
5.把在步骤3和4中选定的二项算子应用到被观察的信号以产生一个修正的被观察信号。在第一次通过该步骤5时,这个二项算子加到实际观察到的信号例如层迭轨迹上;在随后每次通过步骤5时,新选择的二项算子被应用到以前的该步骤5的结果上去。
图1是以简图形式表示陆地反射地震探测系统的某些有关方面;
图2表示:(a)是一个混相子波的被观察信号曲线;(b)是被观察 信号的自    相关曲线;
图3说明:(a1)是通过把图2所示被观察信号与图3所示的具有〔1,0,0,……a〕形式的(c1)算子进行褶积运算所形成的修正的被观察信号的曲线,(a2)是图2的同一被观察信号被图3所示具有〔a,0,0,……1〕形式的算子(c2)进行褶积运算后的曲线;(b)是图3中修正的被观察信号((a1)或(a2),它们具有相同的自    相关)的自    相关的右半部的曲线;(c1)是从自    相关函数导出的算子的曲线,其中自    相关函数表示在图2的(b)中而算子具有〔1,0,0,……0,a〕的形式;(c2)也是从图2(b)所示的函数中导出的算子的曲线,但这个算子具有〔a,0,0,……0,1〕的形式。
图4表示:(a1)是修正的被观察信号的曲线,这个信号是通过把具有较大Varimax量(表示的图3的(a1)中)的修正被观察信号与一个算子进行褶积运算形成的,这个算子是从图3(b)中的自    相关函数导出的,取〔1,0,……0,a〕的形式;(a2)是把(a1)中的同样的信号与具有〔a,0,……0,1〕形式的算子进行褶积运算所形成信号的曲线;(b)是表示在图4的(a1)和(a2)的信号的自    相关曲线;(c1)是累积算子的曲线,这个累积算子是通过把图3所示的具有较大Varimax量(表示在图3的(cl)中)的算子与从表示在图3(b)中的自    相关函数导出并具有〔1,0,……0,a〕形式的算子进行褶积运算得到的。(c2)是一个累积算子的曲线,它与(c1)中的算子相对应,但具有〔a,0,……0,1〕的形式。
图5表示与图4相似类型的曲线,但它是在七次迭代以后在该过程的较晚阶段导出的。
图6是在处理被观察信号对它们进行消褶积并产生地震图的过程中所包括的主要步骤的流程图;
图7表示被观察信号可能类型的性质;
图8表示对一个几乎全通的子波的识别和处理;
图9表示对一个具有复数因式分解的子波的处理;
图10表示被观察信号位于某一高斯参数内的本发明的用处。
图1表示使用本发明的典型实施例在陆地上进行反射地震勘探的典型布置。在地面上放一个地震能源10,这个地震能源可以象在所谓“连续震动法(Vibroseis)”技术中那样是一辆卡车,它使用一个震动器对地球施加机械震动,它也可以是其它类型的适用的地震能源。地震能波向下传播,通过地下岩层12,14和16。本图中表示了这种波的几条路径:路径18a向下传播到岩层12和14间的界面上,然后沿反射路径18b(作为主反射)返回地面,被地音探测器20接收。路径22a向下传播到同一反射体,然后沿反射路径22b(另一主反射)返回到探测器24。同样作为主反射,图中还示出了向下传播到岩层14与16之间的反射体和到达岩层16底部然后反射回到接收机20和24的各条类似路径。此外,还表示出探测器26接收的一个多次反射:沿路径28a向下传播到达位于岩层16底部的反射体,然后向上沿28b反射,向下沿28c反射,再向上沿28d反射,再向下沿28e反射,又一次向上沿28f反射。为了简化起见,没有表示为路径从一个岩层进入另一个岩层时的折射。地音探测器的输出可以是被示意表示为30的一个或多个放大器放大,也可被示意表示的一个或几个滤波器32滤波,并在送往存储和处理系统36之前在34将其数字化。在地面特定位置向地下发送适当数量的地震能量脉冲并用适当的接收机记录,然后将能源10移到新的地点。各个接收机也可以移到各个新地点,这个过程可以重复多次,以便对所关心的地下岩层进行信号多重覆盖。
为了简化,假设各接收机的输出信号已运用已知技术通过层迭步骤被处理,假设层迭轨迹可以画成一个时间序列曲线,这个序列包括具有图2的40-0中所示形状的子波,图2中纵轴是振幅,横轴是向右递增的时间。由于噪声及各种各样的干扰因素,其中包括最上边几个岩层的地震能的混响,所以层迭轨迹本身并不是能够满足要求的按照地下反射体 的时间排列的反射波(或是反射系数和信置)的时间图。假设在这种情况下层迭轨迹40-0(被观察的信号或者更准确地说在这种情况下是其被表示出来的那一部分)是一个混合相位子波而不是最小相位(延迟)子波,而且它离高斯分布不太远,那么最小熵消褶积可以很好地奏效。按照本发明的一个典型实施例,诸如层迭轨迹40-0这样的被观察信号由于向它连续施加算子而被修正,每个后续的算子都加到随着每个前面的算子的应用而进行累积修正的被观察信号上,其中每个算子的选定要根据该算子的应用所能减小熵的程度来作出,只要它能有效地减小熵,这个修正就一直延续下去,因此最终得到的累积修正的被观察信号就可以作为对按时间排列的有关地下反射系数的振幅和位置的良好估计而加以应用。为此,为了导出零阶算子,图2线(a)上的被观察信号40-0被自    相关,以导出位于线(b)上的曲线42-0,曲线42-0仅表示了该自    相关结果的右半边。所关心的二项算子的a项是用上面简述过并且在下边还要详细叙述的方式从曲线42-0中导出来的。
在图3中,位于(a1)线上的曲线40-1是被观察信号40-0,它被施加了图3线(c1)上44-1所示的二项算子〔1,0,……0.a〕而被修正,而位于(a2)线上的曲线40-1′是修正的被观察信号,获得它的方式与前面提到的方式相同,只是其二项算子取图3线(c2)上的44-1′所示的〔a,0,……0,1〕的形式。图3线(b)上的曲线42-1是曲线40-1和40-1′的自    相关(只是它的右半边),40-1和40-1′有相同的自    相关。如图3所示,在40-1处的修正的被观察信号的Varimax量V比在40-1′处的修正的被观察信号的相应量要大,因此用于修正被观察信号40-0以产生被观察信号40-1的具有〔1,0,0,……a〕形式的算子被选择供进一步使用。
在图4中,(a1)线上的曲线40-2表示了这样的褶积运算结果,这个褶积是图3中所示具有较大Varimax量的曲线〔图3中(a1)线上的曲线40-1〕与从图3线(b)上的自    相关函数42-1导出的二项算子的褶积, 该算子以〔1,0,……0,a〕的形式被使用。图4中(a2)线上的曲线40-2′也是修正的被观察信号,它是通过把具有较大Varimax量的图3所示的同一曲线(曲线40-1)与从图3中的42-1曲线导出的同一二项算子进行褶积运算得到的,但算子具有〔a,0,……0,1〕的形式。
在图4中,(c1)线上的曲线44-2表示累积算子,这个算子是由图3中所示的二项算子44-1〔它导致了图3中(a1)线上的较大Varimax曲线40-1〕与从图3所示曲线42-1导出的并且具有〔1,0,……0,a〕形式的二项算子进行褶积运算得到的。图4中(c2)线上的曲线44-2′表示了一个类似的累积算子,它与44-2形成的方式相同,但它所用的二项式算子是从曲线42-1导出的,并且具有〔a,0,……0,1〕的形式。
在图5所表示的各条曲线,尖峰脉冲和Varimax量说明了该处理过程的较后的阶段(迭代),这一点下面还要详细叙述,它与前面阶段(迭代)相应量的关系与图4和图3间的关系相同。
从图2-5可以看出,在各个图中表示的具有较大Varimax量(低熵)的修正过的被观察信号,在加上了附加的二项算子后趋于给出更详细地信息,而且有关的累积算子44变得更加复杂正如从图5中看到的,七次迭代(通过该整个过程)以后的结果是一个累积的修正的被观察信号40-7′,它与图2中所示的最初的被观察信号40-0相比要详细得多,并且它具有更大的Varimax量,因此熵较低。该处理过程还可以继续再进行若干次迭代,只要新修正的被观察信号的Varimax量显出足够的改进即可,即如果新的Varimax量至少要比前一个大出一个选定的阈值。
就实施本发明的一个处理过程的无限制的例子而言,当把它运用到地震勘探时,其主要步骤表示在图6中。它从步骤50开始,在这里收集地震接收器1,2,……k的输出,并把这些输出送到步骤52上去,在52中如果需要的话,可对这些输出进行普通的预处理,比如经过层迭这一步。在步骤54采用一个层迭轨迹作为被观察的地震信号y(i),式中 i=1,2,……I,而数字化信号的各个采样y(i),覆盖的时间相当长(几秒钟),足以使从所关心的最深反射体来的主反射信号返回来。采样频率高到能将被观察的模拟地震信号中的有用最高频率考虑进去。在步骤56中运用已有技术中任何一个合适的自 相关技术对数字化的被观察地震信号y(i)进行自 相关,以产生自 相关信号R(yy,i),R(yy,i)也采用数字化形式并由以标号i表示的各个采样组成。然后再次使用已知技术求出这个自 相关函数的最大偏移滞后d,它是用使中心瓣与右边最靠近的瓣(例如,右边第一个正的或负的峰)分离的采样数i来求得的。图2说明了自 相关函数42-0的这个滞后d值,在这个例子中d是从起始点到最大正峰或负峰的距离。总之,虽然不是必须但最好使用到最大峰值(以绝对值表示)的距离,不管它是正峰还是负峰。然而,至少在某些迭代中有可能要使用到第一正峰的距离,或使用到第一负峰的距离。在步骤60中要做一次检查,以确定a项是否是趋于零或无穷大,这种情况当震源子波不是δ函数而它的自 相关是δ的函数时可能发生。如果在步骤60的检查确定二项算子的非零非1项a并不趋于零或无穷大,则该过程转到步骤62,以确定根据a是其最后一项的算子导出的Varimax量V(A)是否大于或等于选定的乘子B剩以根据其中a是第一项的算子导出的Varimax量V(A′)的积,例如乘子B可以是1.05。一个较方便的Varimax量V是V=〔y(i)〕4/〔y(i)22,式中y(i)是被现有的和所有前面的算子累积修正的被观察信号。当从步骤62获得肯定结果时,则使用已知的褶积技术用二项算子A去对现有的被观察信号进行褶积运算以便在步骤64导出新的被观察信号。当在步骤62获得否定结果时,二项算子A′(其中a是最后一项)为了相同目的而在步骤66中应用。假如步骤60中的检查表明:a项趋近于零或无穷大,则该过程转到步骤68以求出a项,a项在用做上述二项式算子的非零项和非1项时,使诸如V这样的Varimax项最大化(或者至少大大提高它),即把通过使用这样的算 子所导出的修正的被观察信号的熵减至最小(或至少是大大减小),或使该结果最大地偏离高斯分布。步骤68的结果用的步骤70中作为二项算子A″的一部分,A″由现有的被观察信号褶积以导出新的被观察信号。如上所述在步骤64,66或70中导出的新的被观察信号在步骤27中检查,看它的熵是否明显地减小。例如,做一个检查来确定新的被观察信号的Varimax量(诸如V)是否比它前面紧接的被观察信号的Varimax值加上阈值t之和还要大。如果结果是肯定的,表明有明显的改善,该过程回到步骤54,而且这时这里所用的被观察地震信号是由步骤64,66和70中适合的一个导出的。如果步骤72指出由最新导出的二项算子对熵的减小并无明显的改善,那么现在的被观察信号于步骤74被贮存。这个信号可以用在步骤75中作为消褶积的被观察地震信号。一个典型的应用是通过使用上面讨论的已知技术制作所关心的地下岩层的地震图。这样的地震图得益于通过使用前面叙述的本发明处理过程而获得的改进的消褶积结果。最终得到的地震图可以在一个屏上显示出来,或用绘图仪绘在纸上,或者用其它方法贮存和(或)使它能用在地下资源的勘探和(或)开发上。
按照本发明,对存在于被观察信号中的几乎全通的子波的识别和处理表示在图8中。子波x几乎是一个全通子波,其中它的自    相关R几乎是一个尖峰脉冲。正如本技术中已知的那样,在被观察信号中有这样一个子波存在几乎不或者根本就不对被观察信号的振幅产生影响,但它会引起明显的相移。一个做为本发明的一部分所揭示的用来证实被观察信号不含有全通子波的实用方法是检查子波(被观察信号)与它的立方函数之间的相关。该相关在图8中表示为G(xx)。如果没有全通子波存在,则相关G(xx)象自    相关一样是一个尖峰形。图8中的例子情况不是这样的。相关G(xx)不是象x的自    相关那样的尖峰形,这表明有几乎全通的子波x存在。甚至在由全通子波与另一个不是全通的波进行褶 积产生复杂得多的函数的情况下,该检查也能从中探测全通子波的存在。
图9说明了在甚至带有含有复数因式分解的子波的被观察信号的情况下,本发明的技术的有效性。图9中的三项子波A具有复数因式分解。已有技术中的维纳(Wiener)(或最小二乘)技术被采用来获得子波A的反演B。表示在C中并根据本发明使用二项算子系列(它们仅为实数)导出的EGD反演给出了反演C。
图10说明使用的被观察信号具有某些高斯分布特征时本发明技术的有效性。在每种情况中的剩余子波是由设计的子波与一个100点随机系列的褶积的结果,然后按照本发明求出和应用二项算子系列。图10给出的P-值表明被观察信号的广义高斯参数(设计子波与100点随机系列褶积的结果)。如图中所见,本发明至少在被观察信号所具有的P-值大约比1.5的情况下能给出稳定的消褶积。
上面结合图6所述的典型处理过程可用图1所示意表示的存储和处理系统36来完成,该系统可由适当规模和结构的通用数字计算机及完成上述功能所必须的适当的外部设备组成。一个适用的数字计算机的例子可以从数字设备公司得到,型号为VAX,它可以按编制的程序(如用Fortran一类适用的高级语言描述图10的各步骤并运用现成的子程序)完成上述的处理过程。换一个方案,上面结合图10讨论的某些或所有步骤可以用硬线连接和(或)固线连接电路来完成。
应该指出,为了画图方便,图3-5中位于(c1)和(c2)线上的各算子都表示成沿时间轴(横轴)有一定的宽度,然而实际上它们只代表振幅,没有宽度。在图8-10中,用了另一种表示方法去表示无宽度的振幅(沿纵轴)。还应该指出,在特定情况下有可能滞后(如图2-5中的d值所示)如此之小,以至于二项算子的各项不被零分开。在这种情况下,算子可具有〔1,a〕或〔a,1〕的形式,人们会明白,这样的算子只是算子〔1,0,……a〕和〔a,0,……1〕的特例,而且该特例也包括在本 发明的范围内。
对本发明系统和处理过程应用于地震信号的有效性的概念说明在图7中给出。正如前面指出的,如果地震数据有占据图上标记的矩形区域的趋势,则已有技术中最小熵和维纳消褶积技术应用的数据范围比本发明的熵引导消褶积技术应用的数据范围要小得多。

Claims (17)

1、一种地震勘探方法,其特征在于:
导出一个被观察地震信号,该信号与地震源子波和地下岩层的相互作用有关;
用一系列算子对被观察信号进行累积修正,而这些算子是根据最小二乘方法和根据它们对减少熵的有效性来加以选择的;
使用累积修正的被观察信号作为对上述地下岩层的地震性质的估计。
2、根据权利要求1所述的方法,其特征在于:在反射地震勘探中,被观察地震信号是根据一个或多个地音探测器的输出导出的。
3、根据权利要求1所述的方法,其特征在于:被观察地震信号是从垂直地震剖面图中的一个或多个接收机的输出导出的。
4、根据权利要求1所述的方法,其特征在于:被观察地震信号不是最小相位或者离高斯分布不很远。
5、根据权利要求1所述的方法,其特征在于:被观察地震信号是混合相位。
6、根据权利要求1所述的方法,其特征在于:算子有被零项分隔的两个非零项。
7、根据权利要求6所述的方法,其特征在于:每个算子具有〔1,0,0,……a〕或〔a,0,0,……1〕的形式。
8、根据权利要求7所述的方法,其特征在于:二项算子的a项是作为被观察信号的自相关函数被导出的,而被观察信号应用任何前面算子而被修正。
9、根据权利要求8所述的方法,其特征在于:算子的非零项之间的间隔是这样导出的:即根据被观察信号的自相关的中心和选定旁瓣间的距离导出的,而这个被观察信号由于应用了任何前面的算子而被修正。
10、根据权利要求9所述的方法,其特征在于:每个算子的〔1,0,0,……a〕和〔a,0,0,……1〕形式的选择是这样做出的:即根据当把它们应用于被观察信号时,看哪个形式造成熵的减小幅度更大来决定,而被观察信号由于任何前面算子的应用被累积修正。
11、根据权利要求10所述的方法,其特征在于:熵的减小是用Varimax量来衡量的。
12、根据权利要求11所述的方法,其特征在于:Varimax量取V=〔y(i)〕4/〔y(i)22的形式,式中y(i)是现有的累积修正被观察信号。
13、一种地球物理勘探方法,其特征在于:
导出一个被观察信号,它可能不是相位最小或者是远离高斯分布,它与地震源波和地下岩层之间的相互作用有关;
运用一系列算子修正被观察信号,每个算子的选择则根据最小二乘方法以及根据把它应用于被观察信号所能达到的熵减小程度来决定,被观察信号由于前面各个算子的应用被累积修正;
使用累积修正的被观察信号作为上述地下岩层特征的非时变序列的估计。
14、根据权利要求13所述的方法,其特征在于:被观察信号是地震信号。
15、根据权利要求14所述的方法,其特征在于:被观察信号是反射地震勘探信号。
16、根据权利要求14所述的方法,其特征在于:被观察信号垂直地震剖面信号。
17、根据权利要求13所述的方法,其特征在于:当应用新的算子不能使熵的减小量新的算子超过预定值时,上述修正就停止。
CN85108418A 1984-12-24 1985-11-23 由熵引导的地震信号的消褶积 Expired CN1007555B (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US684,811 1984-12-24
US06/684,811 US4688198A (en) 1984-12-24 1984-12-24 Entropy guided deconvolution of seismic signals

Publications (2)

Publication Number Publication Date
CN85108418A CN85108418A (zh) 1986-07-16
CN1007555B true CN1007555B (zh) 1990-04-11

Family

ID=24749678

Family Applications (1)

Application Number Title Priority Date Filing Date
CN85108418A Expired CN1007555B (zh) 1984-12-24 1985-11-23 由熵引导的地震信号的消褶积

Country Status (7)

Country Link
US (1) US4688198A (zh)
EP (1) EP0186569B1 (zh)
CN (1) CN1007555B (zh)
AT (1) ATE71738T1 (zh)
DE (1) DE3585220D1 (zh)
EG (1) EG19383A (zh)
OA (1) OA08152A (zh)

Families Citing this family (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2616920B1 (fr) * 1987-06-19 1989-10-13 Schlumberger Prospection Inversion d'un profil sismique vertical en minimisant une fonction du type entropie
US4922362A (en) * 1988-03-04 1990-05-01 Schlumberger Technology Corporation Methods for deconvolution of unknown source signatures from unknown waveform data
FR2632734B1 (fr) * 1988-06-10 1990-09-14 Schlumberger Prospection Procede d'etablissement d'un modele stratigraphique du sous-sol a partir d'un profil d'impedance acoustique et d'une section sismique
IL97345A0 (en) * 1991-02-24 1992-05-25 Univ Ramot Method and apparatus for blind deconvolution
WO1994003857A1 (en) * 1992-08-10 1994-02-17 Advanced Logic Research, Inc. Computer interface for concurrently performing plural seeks on plural disk drives
US5454106A (en) * 1993-05-17 1995-09-26 International Business Machines Corporation Database retrieval system using natural language for presenting understood components of an ambiguous query on a user interface
US5400299A (en) * 1993-08-20 1995-03-21 Exxon Production Research Company Seismic vibrator signature deconvolution
US5648937A (en) * 1995-01-18 1997-07-15 Atlantic Richfield Company Method and apparatus for correlating geological structure horizons from velocity data to well observations
US5638338A (en) * 1995-10-16 1997-06-10 Peterson; Fred M. Seismic processing apparatus and method
WO1998011455A1 (en) 1996-09-13 1998-03-19 Pgs Tensor, Inc. Method for time lapse reservoir monitoring
WO2002097472A2 (en) 2001-05-25 2002-12-05 Exxonmobil Upstream Research Company Multiple suppression for ocean bottom seismic data
MXPA05010458A (es) * 2003-04-01 2006-03-21 Exxonmobil Upstream Res Co Fuente vibratoria de alta frecuencia conformada.
CN101813786B (zh) * 2010-04-02 2012-05-23 中国石油集团西北地质研究所 子波处理二步法反褶积方法
RU2632249C1 (ru) 2013-11-08 2017-10-03 Шлюмбергер Текнолоджи Б.В. Определение режима течения для адаптации модели потока
WO2015070008A1 (en) * 2013-11-08 2015-05-14 Schlumberger Canada Limited Spectral analysis with spectrum deconvolution
CN106772573B (zh) * 2016-11-16 2018-06-26 电子科技大学 基于最大相关熵的地震子波信号提取方法
CN106896407B (zh) * 2017-03-28 2018-07-13 吉林大学 一种基于近似负熵的微地震信号初至拾取方法
CN113296157B (zh) * 2020-02-21 2023-05-26 中国石油天然气集团有限公司 一种利用广义高斯分布进行储层的预测方法及装置

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR1495419A (fr) * 1966-02-02 1967-09-22 Mobil Oil Corp Traitement de données géophysiques
US4025721A (en) * 1976-05-04 1977-05-24 Biocommunications Research Corporation Method of and means for adaptively filtering near-stationary noise from speech
US4335373A (en) * 1980-11-07 1982-06-15 Fairchild Camera & Instrument Corp. Method for analyzing a digital-to-analog converter with a nonideal analog-to-digital converter
US4458267A (en) * 1981-10-06 1984-07-03 Analogic Corporation Digital x-ray system

Also Published As

Publication number Publication date
EG19383A (en) 1995-02-28
EP0186569A2 (en) 1986-07-02
DE3585220D1 (de) 1992-02-27
US4688198A (en) 1987-08-18
ATE71738T1 (de) 1992-02-15
EP0186569A3 (en) 1989-02-22
OA08152A (en) 1987-03-31
EP0186569B1 (en) 1992-01-15
CN85108418A (zh) 1986-07-16

Similar Documents

Publication Publication Date Title
AU612552B2 (en) Model-based depth processing of seismic data
Coppens First arrival picking on common‐offset trace collections for automatic estimation of static corrections
Bacon et al. 3-D seismic interpretation
EP1660914B1 (en) Method for attenuating water layer multiples in a seismic survey
CN1007555B (zh) 由熵引导的地震信号的消褶积
US4293933A (en) Well logging apparatus and method: synthetic logs and synthetic seismograms with extrapolated reflector dip from log measurements
US9081107B2 (en) Shot scheduling limits for seismic acquisition with simultaneous source shooting
US5012453A (en) Inverse vertical seismic profiling while drilling
EP0117280B1 (en) Vertical seismic profiling
US6188964B1 (en) Method for using global optimization to the estimation of surface-consistent residual statics
RU2282877C2 (ru) Способ корректировки сейсмических данных при морской сейсмической разведке
AU687590B2 (en) Method of processing seismic data having multiple reflection noise
EP2548052B1 (en) System and method of 3d salt flank vsp imaging with transmitted waves
White et al. Variations in oceanic upper crustal structure in a small area of the north-eastern Atlantic
Pritchett Acquiring better seismic data
US20120269035A1 (en) Evaluating Prospects from P-Wave Seismic Data Using S-Wave Vertical Shear Profile Data
Pasasa et al. Prestack Kirchhoff depth migration of shallow seismic data
US3629798A (en) Method and system for refraction seismic exploration
Cassell Vertical seismic profiles-an introduction
Bennett 3-D seismic refraction for deep exploration targets
US3530430A (en) Method of and apparatus for stacking electrical seismic traces
Talagapu 2D and 3D land seismic data acquisition and seismic data processing
Ziolkowski Simplified wavelet estimation using source-signature measurements
Gaiser et al. VSP fundamentals that improve CDP data interpretation
Willis et al. Characterization of Scattered Waves from Fractures by Estimating the Transfer Function Between Reflected Events Above and Below Each Interval

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C13 Decision
GR02 Examined patent application
C14 Grant of patent or utility model
GR01 Patent grant
C19 Lapse of patent right due to non-payment of the annual fee
CF01 Termination of patent right due to non-payment of annual fee