CN104755961B - 基于蒙特卡洛反向投影的波束反演 - Google Patents
基于蒙特卡洛反向投影的波束反演 Download PDFInfo
- Publication number
- CN104755961B CN104755961B CN201480002816.6A CN201480002816A CN104755961B CN 104755961 B CN104755961 B CN 104755961B CN 201480002816 A CN201480002816 A CN 201480002816A CN 104755961 B CN104755961 B CN 104755961B
- Authority
- CN
- China
- Prior art keywords
- data
- wave beam
- cross
- correlation
- time shift
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 claims abstract description 55
- 230000001427 coherent effect Effects 0.000 claims abstract description 10
- 238000005315 distribution function Methods 0.000 claims description 31
- 230000008859 change Effects 0.000 claims description 18
- 230000001186 cumulative effect Effects 0.000 claims description 15
- 238000005314 correlation function Methods 0.000 claims description 14
- 238000002922 simulated annealing Methods 0.000 claims description 8
- 238000005259 measurement Methods 0.000 claims description 6
- 238000013508 migration Methods 0.000 claims description 5
- 230000005012 migration Effects 0.000 claims description 5
- 238000004364 calculation method Methods 0.000 claims description 2
- 238000004590 computer program Methods 0.000 claims description 2
- 230000001052 transient effect Effects 0.000 claims 2
- 238000006073 displacement reaction Methods 0.000 description 10
- 230000006870 function Effects 0.000 description 7
- 230000008569 process Effects 0.000 description 7
- 238000003384 imaging method Methods 0.000 description 6
- 238000012545 processing Methods 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 4
- 238000004891 communication Methods 0.000 description 4
- 238000003325 tomography Methods 0.000 description 4
- 238000013459 approach Methods 0.000 description 3
- 238000004587 chromatography analysis Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000001816 cooling Methods 0.000 description 2
- 230000007423 decrease Effects 0.000 description 2
- 239000004744 fabric Substances 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 230000009191 jumping Effects 0.000 description 2
- 230000005055 memory storage Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000000926 separation method Methods 0.000 description 2
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000033228 biological regulation Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000005266 casting Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000011549 displacement method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000011221 initial treatment Methods 0.000 description 1
- 238000003780 insertion Methods 0.000 description 1
- 230000037431 insertion Effects 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 230000000873 masking effect Effects 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 150000003839 salts Chemical class 0.000 description 1
- 230000009329 sexual behaviour Effects 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 238000011282 treatment Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/51—Migration
- G01V2210/512—Pre-stack
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/66—Subsurface modeling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Remote Sensing (AREA)
- General Physics & Mathematics (AREA)
- Geology (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- Theoretical Computer Science (AREA)
- General Engineering & Computer Science (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Geophysics And Detection Of Objects (AREA)
- Radar Systems Or Details Thereof (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及一种用于最小化感兴趣的地下区域的地震图像中的伪像的方法和***,其中所述图像被确定为从记录的地震数据导出的数据波束集合以及至少部分地从与地下区域有关的速度模型导出的建模的波束集合。可能由周期跳跃和相干噪声产生的所述伪像导致了所述建模的波束集合与数据波束集合的未对齐。本发明利用蒙特卡洛反演技术以更新所述速度模型并且因此最小化所述地震图像中的所述伪像。
Description
技术领域
本发明一般涉及用于地震成像和地球建模的方法和***,并且更具体地涉及最小化周期跳跃(cycle-skipping)和与相干噪声的假对齐的影响的方法,所述影响降低先前公开的行程时间反射层析成像方法的性能。
背景技术
例如美国专利申请序列号12/606861(“所述′861申请”)中所描述的波束层析成像利用地震图像通过发现改善波束对齐的速度修正来改善对地震速度的估计。所述′861申请描述的波束层析成像方法包括确定速度模型修正的步骤,所述速度模型修正改善了由局部操纵记录的数据而形成的波束与由使用当前速度模型和当前图像对所述记录的波束建模而形成的波束之间的匹配。因此,速度模型更新可以通过最小化记录的数据与由前向建模计算形成的合成数据之间的差异来获得。
所述′861申请描述了所谓的最大相关移位法(MCS),所述MCS包括将数据和建模的波束互相关以确定数据和建模的波束对之间的定量匹配的步骤。所述波束对中的每一对均对应于穿过地球的单条射线路径。完整的射线路径具有始于震源的一段和终于接收器的另一段。所述波束是沿着它们关联的射线路径的行程时间的函数。成一对的两条波束在行程时间移位的范围内是互相关的。最大化每个波束对的互相关的时移将被获得并用作所述层析成像的数据值。为多条射线路径完成基于互相关的匹配,对当前应用通常是107的数量级。
接下来,沿每条射线路径的速度必须被改变以符合所述数据与建模的波束之间测得的时移。每条射线路径和对应的时移生成方程,所述方程将沿所述射线路径的所述速度改变与测得的行程时间移位线性联系起来。对于波束层析成像的MCS实施方法而言,这一大的稀疏方程组通过大的稀疏矩阵的迭代再加权的最小二乘反演来求解。用于MSC方法的输入数据是带有单个移位值的对齐移位,所述单个移位值代表每个建模的波束和数据波束对的恰当对齐。MCS实施方式的成功依赖于互相关提供波束对之间的未对齐的良好测量。不幸的是,使得所述互相关歪曲正确对齐的周期跳跃可降低所述测量的质量。
在某种程度上,周期跳跃可以通过驱使再加权的最小二乘解符合L1范数来容忍,这更好地处理了由周期跳跃导致的异常数据。然而,所述再加权过程仅能处理很小程度的周期跳跃,这意味着所述数据必须具有低级别的噪声。
图1(A)和1(B)示出了周期跳跃问题的一个例子。图1(A)和1(B)示出了多个建模的波束和数据波束对:深色轨迹是通过局部倾斜堆叠记录的数据而计算出的波束,而浅色轨迹是通过使用描述地下速度结构的地球模型和描述地下反射体的地震图像对计算建模而计算的。浅色和深色轨迹被成对示出,以使得每对波束都对应于穿过地球的一条射线路径。图1(A)示出了未对齐的轨迹,而图1(B)示出了对齐后的轨迹。注意,穿过图1(A)和1(B)中部的水平线是对应于射线轨迹到达时间的时间线。
参考图1(A),每对波束的两条轨迹均有由地球模型和图像中的误差引起的一定量的未对齐。大多数轨迹的恰当对齐是明确清晰的并且可通过寻找互相关最大值来容易地测量。在图1(B)中,通过对轨迹移位以获得浅色与深色轨迹之间的最大互相关来对齐所述轨迹。图1(B)顶部的线条(“移位”)示出了最大化所述互相关所需的时移。然而,注意图1(B)中部附近的中间偏移(mid offsets)。在这里,使用时移的快速变化来对齐轨迹。这些快速变化由周期跳跃引起,因此最大互相关对齐不同于基于波形的一个或多个周期的正确对齐。因此,只从数据中并不清楚互相关函数中的哪个局部最大值对应于正确对齐。由于噪声或数值近似,全局最大值可能并不正确。
此外,在数据包含由如图2(A)中的多个反射引起的强相干噪声的情形中,轨迹可被移位以形成图2(B)中心附近所示的明显对齐好的事件。然而,互相关最大值可能并不对应于恰当的轨迹对齐。虽然图2(B)中部附近看起来具有令人信服的事件的对齐,但图表顶部绘出的时移的极端不规则暗示存在由相干噪声引起的假对齐。
因此,先前公开的MCS方法存在的问题是:测量波束对的波束之间的时移的步骤独立于用于地球模型中的速度修正的反演的步骤。只通过寻找最大互相关而执行的对齐倾向于周期跳跃,因为没有区别基于移位是否对应于速度的合理改变,和基于此移位是否与其他波束对之间测得的移位一致。测得的移位是固定的并随后用在速度的线性反演中。所述线性反演没有对修正由周期跳跃引起的异常测量或对丢弃与相干噪声假对齐的移位做出规定。
因此,需要出现一种较少受到周期跳跃和与相干噪声的假对齐影响的行程时间反射层析成像方法。
发明内容
本发明提供了一种用于最小化感兴趣的地下区域的地震图像中的伪像(artifacts)的计算机实施的方法,所述图像基于从记录的地震数据导出的数据波束集合以及至少部分地从与地下区域有关的速度模型中导出的建模的波束集合。所述方法包括以下步骤:确定所述建模的波束集合与数据波束集合在预定时移范围内的互相关函数;使用所述互相关函数确定在所述预定时移范围内的概率分布函数;使用所述概率分布函数确定在所述预定时移范围内的累积分布函数;生成预定间隔内的随机数;从预定时移范围内选择时移值以便所述随机生成的数等于所述累积分布函数;并且使用所述时移值为所述速度模型计算速度修正。
因此,本发明是一种将诸如波束对齐中的周期跳跃和由于相干噪声引起的假对齐等的某些伪像的影响最小化的蒙特卡洛方法。所述方法组合了模拟退火蒙特卡洛技术和反向投影技术,所述模拟退火蒙特卡洛技术先前已被用于反射静力学(Rothman,1985,1986)和更一般的速度模型反演(Sambridge等,2002)。即,所述方法组合了移位值的蒙特卡洛概率性选择和反向投影速度反演技术的迭代。所述组合可以在无法只通过互相关实现正确的周期对齐并且相干噪声引起假对齐的情景中产生速度估计。
附图说明
参照附图中描述的本发明的具体实施例来详细描述本发明。所述图仅描述本发明的典型实施例,并且因此不被认为限制它的范围。
图1(A)和1(B)是示出对应于第一位置的建模的波束与地震数据轨迹的未对齐的图。
图2(A)和2(B)是示出对应于第二位置的建模的波束与地震数据轨迹的未对齐的图。
图3是根据本发明的实施例示出波束对齐方法的流程图。
图4是根据本发明的实施例示出蒙特卡洛反演方法的流程图。
图5(A)-(C)分别是以下的图形表示:(A)对应于方程(1)的互相关函数(时移(τ)的函数);(B)对应于方程(2)的概率分布函数;以及(C)对应于方程(3)的累积分布函数。
图6是累积分布函数,带有生成的示例随机数R以根据方程(4)寻找时移。
图7(A)-(B)分别是以下的图形表示:(A)慢度(slowness)的改变沿射线路径穿过的单元被平滑分布,以便所述行程时间改变等于τi;以及(B)单元内的行程时间改变作为穿过所述单元的所有射线的贡献的加权平均值来计算。
图8(A)和8(B)分别是以下的图形表示:(A)高温处对应于方程(2)的概率分布函数;以及(B)高温处对应于方程(3)的累积分布函数。
图9(A)和9(B)分别是以下的图形表示:(A)低温处对应于方程(2)的概率分布函数;以及(B)低温处对应于方程(3)的累积分布函数。
图10(A)和10(B)分别是以下的图形表示:(A)在大的Θ值的情况下对应于方程(2)的概率分布函数;以及(B)在小的Θ值的情况下对应于方程(2)的概率分布函数。
图11是根据本发明的***的示意图。
具体实施方式
本发明可以在要由计算机执行的计算机方法和***的一般背景中描述和实施。这种计算机可执行指令可以包括可被用以执行特定任务和处理抽象数据类型的程序、例程、对象、组件、数据结构、和计算机软件技术。本发明的软件实施方式可以以不同语言编码用于各种计算平台和环境中的应用。将要理解,本发明的范围和基本原理不限于任何具体的计算机软件技术。
而且,本领域技术人员会理解本发明可以使用硬件和软件配置的任何一种或组合来实践,包括但不限于,具有单个和/或多个处理器的计算机处理器***、手持设备、可编程消费电子、迷你计算机、大型计算机、超级计算机等。本发明还可以实践在分布式计算环境中,其中任务由通过一个或多个数据通信网络链接的服务器或其他处理设备执行。在分布式计算环境中,程序模块可以位于包括存储器存储设备的本地和远程计算机存储介质这两者中。
另外,诸如CD、预录盘或其他等价设备等的与计算机处理器一同使用的制成品,可以包括其上记录有用于使得所述计算机处理器便于本发明的实践和实施的计算机程序存储介质和程序装置。这种设备和制成品也落入到本发明的精神和范围内。
现在参考图,本发明的实施例会被描述。本发明可以多种方式被实施为,包括例如***(包括计算机处理***)、方法(包括计算机实施的方法)、装置、计算机可读介质、计算机程序产品、图形用户界面、门户网站、或有形地固定在计算机可读存储器中的数据结构。本发明的若干实施例在以下被讨论。所述附图仅说明了本发明的典型实施例并且因此不被认为限制本发明的范围和广度。
图3是根据本发明的计算机实施的波束对齐方法10的流程图,所述方法使用速度模型的蒙特卡洛更新。图3示出了类似美国专利申请序列号12/606,861中描述的MCS发明的方法,所述美国专利申请序列号12/606,861以其全部被并入在此。根据下面提供的旨在作为说明的方法10的操作的一个或多个实施例,所述方法10可用于改善地震速度模型并且生成与所述地球的地下区域有关的地震图像。在某些实施例中,方法10可以使用一个或多个未描述的附加操作,和/或不使用一个或多个已讨论的操作来完成。此外,图3和4中说明的方法10的操作顺序以及下面描述的顺序并不旨在作为限制。
在一种实施例中,方法10包括以下步骤:访问记录的地震数据的步骤12,以及变换所述记录的地震以形成数据波束集合14。所述数据波束集合可被存储于计算机存储介质中。通常,所述记录的地震数据样本已经经过初步处理以提高信噪比和调整这些数据用于后续图像处理。在某些实施例中,例如在Hill,N.R.,Gaussian Beam Migration,Geophysics,Volume55,pp.1416-28(1990)和Hill,N.R.,Prestack Gaussian BeamMigration,Geophsyics,Volume66,pp.1240-50(2001)中描述的波束形成变换可用于所述波束形成变换操作,然而,本领域的技术人员会理解,也可使用诸如由Sun,Y等的3-DPrestack Kirchhoff Beam Migration for Depth Imaging,Geophysics,Volume65,pp.1592-1603(2000)中描述的波束方法等的其他方法。
在步骤18处,初始地球模型16和数据波束集合14接着输入到迁移中,这形成地下区域的地震图像。所述初始地球模型16通常已经由地质学解释和地球物理分析构建。例如,它可以包含由现有地震图像的大量地质学解释构建的盐质量模型,并且所述地下速度模型可能已经被本领域技术人员已知的任何数量的层析成像方法确定。
在步骤18生成地震图像后,步骤19决定成像后是否跟着速度模型的更新。不进行更新的原因包括判定步骤18生成的图像是满意的,例如,由于所述图像与关于所述感兴趣的地下区域的已知的地质学信息或其他信息相符。不进行更新的另一原因会是前一次迭代未对速度模型有显著更新。
在步骤20处,例如,如美国专利申请序列号12/606,861中描述的由地震图像和地球模型导出的建模的波束(在下文中称为“模型波束”,“建模的波束”,“模型波束集合”),与对应的数据波束在所述建模的波束与数据波束之间的时移范围内互相关。在步骤20,用于第i个波束对的作为时移τ的函数的互相关函数Ei(τ)由如下的方程1定义:
其中第i条射线的建模的波束为并且数据波束为模型波束和数据波束可例如通过美国专利申请序列号12/606,861中描述的方法来计算。积分间隔[T1,T2]延续一个预定义的时间窗,所述时间窗以由射线路径行程时间计算出的时延为中心。所述延时通常被选择成是沿源射线段从源点到离接收器射线段最近的段上的点的行程时间,加上沿所述接收器射线段从接收器到离所述源射线段最近的接收器段上的点的行程时间。图5(A)示出了示例性互相关函数。(参见对于步骤140的以下讨论。)接着在步骤22中所述互相关可在预定时移范围内存储。
回来参考图3,在某些实施例中,一旦步骤24确定再也没有波束需要建模,本发明的方法继续到参考图4在下面描述的步骤28的反演操作。反演操作28试图寻找对速度模型的修正,所述修正会修正建模的波束与数据波束之间剩余的未对齐。每个建模的波束都与具有两段的射线路径相关联。一段代表源波场的传播(下文记为“源射线段”),并且另外一段代表接收器波场的传播(下文记为“接收器射线段”)。所述两条射线段不必相交。这两条射线段由它们关联的所述记录的地震数据的波束成分和所述地球模型速度的波束成分所确定;它们不依赖于反射体结构的任何模型。这些射线路径穿过速度模型修正网格中的单元。所有所述单元中的速度均被调整以对所述建模的波束与数据波束之间的测得的时移建模。此建模为在步骤22中存储的每个波束对生成方程。沿所述源射线和所述接收器射线的慢度修正的积分应当等于对齐所述建模的波束和数据波束所需的时移。慢度被定义为速度的倒数。来自步骤28的更新的速度模型(图3和4中的圆圈B)接着提供给步骤18的图像更新步骤。
注意,对于方法10的第一次迭代,步骤18至24以图3示出的顺序操作。对后续迭代,步骤18至24在下述的图4的步骤340与360之间操作。
图4示出了根据本发明的蒙特卡洛反演方法(图3的步骤28)的流程图。根据所述方法,测得的时移数据被反演用于速度模型修正。蒙特卡洛反演的输入数据是初始参数值(步骤120)以及流程图(图3)的步骤22处确定的在模型波束与数据波束对之间的互相关。最重要的参数是会在下面详述的温度T和时间间隔Θ。对齐移位的蒙特卡洛选择开始于图4的步骤160。在步骤140中,互相关函数Ei(τ)被检索。
下一步,在步骤160中,根据方程(2)确定概率分布函数pdf(τ):
对于τ∈[τ1,τ2],其中τ1和τ2是移位的预定极限。图5(B)示出了概率分布函数pdf(τ)的示例图。方程(2)示出的函数是在模拟退火的当前迭代期间选择时移值τ以对齐第i个轨迹对的非归一化概率。值是通过对沿射线路径的慢度改变进行积分计算出的时移。所述“慢度”包括对在反演过程的较早迭代期间已经累积的模型慢度改变的更新。参数Θ是当前时间间隔参数,T是当前数值温度参数。随着速度模型被更新,参数Θ和T将在下述的图4的步骤360中的迭代数值冷却过程期间下降。
仍然参考图4,在步骤180中,累积分布函数cdf(τ)通过方程(3)在间隔τ∈[τ1,τ2]内被计算:
图5(C)示出了累积分布函数cdf(τ)的示例略图。注意cdf(τ)的取值范围在0和1之间,且不随τ的增加而下降。在步骤200处,在区间[0,1]内产生均匀随机数R。τi的值在步骤220中确定以便
cdf(τi)=R 方程(4)
图6示出了为通过方程(4)寻找时移而生成的示例随机数R。τi的值用在步骤240中通过例如如美国专利申请序列号12/606,861中公开的反向投影方法来计算速度修正。图7(A)简略绘出了所述反向投影。注意,慢度的改变分布在被第i条射线穿过的速度模型单元内。值是从第i条射线计算出的第j个单元的慢度改变。所述分布被这样形成使得所述慢度改变沿射线路径平滑变化,并且使得所述改变会令总行程时间改变等于τi。解释性判定被经常作出以只允许慢度在所述速度模型中的受限区域中更新。例如,可以认为断层附近的速度是最不确定的。在这种情况下,所述速度模型可能被掩蔽以只允许在这些断层附近的更新。在这种情况下,只有当路径穿过模型的掩蔽区域时,才会允许沿着所述射线路径的慢度改变。
回去参考图4,如果步骤260处有更多波束对要处理,则接着选择下一对(步骤280),且所述方法返回步骤160计算下一个概率分布函数。然而,如果没有更多的波束对,则所述方法继续步骤320,其中为所述波束的所有波束累积的慢度改变计算加权平均值之和。通常来说,加权简单地是图7(B)中示出的单元中的射线路径的长度。第j个单元中慢度改变的加权平均值之和根据方程(5)确定:
其中lij是第i条射线在第j个单元中的长度。N的值是波束对的总数。慢度改变用以计算速度改变。然而,当单元中的速度值更新时,所述改变根据方程(6)在步骤340中急剧衰减:
Vj→Vj+α·ΔVj 方程(6)
在修正被添加到速度模型用于后续迭代之前,所述修正通常以范围0.05到0.1之间的因子α衰减。在每次迭代期间此急剧衰减仅允许对速度更新有小的改变,但这并非无效,因为在模拟退火冷却过程期间会使用多次迭代。
下一步,在步骤360处,数值温度T和时间Θ间隔都被降低。在一种实施例中,T和Θ的值分别被乘法因子AT和AΘ衰减,所述乘法因子AT和AΘ贯穿所述过程的迭代期间稍微不一致但保持不变。T的值在第一次迭代中通常很大,使得概率分布函数pdf(τ)初始并不很大程度上依赖于互相关Ei(τ),例如如图8(A)所示。由于T的这一大的初始值,累积分布函数cdf(τ)平滑增长(见图8(B)),这意味着τ∈[τ1,τ2]范围内的所有τ值都会以近似相等的概率被选择。因此,蒙特卡洛过程选择的τ的值会几乎均匀分布于区间[τ1,τ2]。
然而,例如如图9(A)所示,随着T在迭代过程中下降,具有大互相关的移位会比具有小互相关的移位更可能出现。一旦温度变小,累积分布函数cdf(τ)会具有明显的阶梯(见图9(B)),这些阶梯出现在pdf(τ)的峰处。在此情况中,满足方程(4)的τi的值最可能出现在概率分布函数的峰处。因此,在低温时,概率分布函数在互相关函数取得局部极大值时剧烈地达到峰值。这些峰导致累积分布函数出现明显的阶梯。蒙特卡洛过程可能选择互相关中的峰附近的τi。
温度T的初始值被选择使得概率分布函数不会出现例如如图8(A)所示的窄峰。选择的最终温度足够低使得存在有如图9(A)所示的明显的峰,但并不低到只存在一个主导峰。选择的时间间隔Θ的初始值比互相关函数的若干周期宽,产生如图10(A)所示的概率分布函数。Θ的最终值通常略小于互相关函数的典型周期,以便产生如图10(B)所示的概率分布函数。如果蒙特卡洛过程的迭代朝着速度模型的良好估计进行,则值应当在后面迭代中接近正确。接着随着Θ变小,方程(2)右边的第一项会抑制与当前建模的行程时间较远的大概率值,这假定了这些远离的峰是由周期跳跃或诸如多重反射等的与相干噪声的假对齐导致的。总之,Θ的初始值要选的足够大以避免对互相关的任何一个极大值有偏见。随着迭代进行,Θ的值下降,因为远离时移的峰可能是诸如周期跳跃或与相干噪声的假对齐等的伪像。给定乘法因子AT和AΘ的值使得T和Θ在预定次数的迭代内从预定初始值变为预定最终值。
例如如图3的步骤18所示,对于当前发明,在过程的外层循环中并不频繁进行图像的更新。然而,替代的,通过将成像步骤移至内层循环可以改善所述过程的性能。例如,步骤18至24可以从图3的外层循环中移动并作为内层循环***在图4中的步骤340与360之间。
替代的,当前基于行程时间的对齐方法可被一般化以包括波形匹配方法(Ribodetti等,2011)。
图11示意说明了其中可实施本文描述的各种实施例的实施方式的计算机网络84的示例。计算机网络84可以包括可被实施为任何传统个人计算机或服务器的数据处理***或计算机***88。然而,本领域技术人员会理解,本文描述的各种技术的实施方式可以实践在其他计算机***配置中,包括超文本传输协议(HTTP)服务器、手持设备、多处理器***、基于微处理器的或可编程的消费电子、网络PC、迷你计算机、Linux计算机、大型计算机等。
包括至少一个处理器的计算机***88可以与计算机存储介质通信,例如盘存储设备86或盘存储器设备96,所述设备也可以是外部硬盘存储设备。盘存储器86和96预期是传统硬盘驱动器,并且因此通过局域网或远程访问来实施。当然,尽管盘存储设备86和96被示出为分离设备,但是如所需,单个盘存储设备也可用以存储任意和全部程序指令、测量数据和结果。
在一种实施方式中,与感兴趣的地下区域有关的数据会作为计算机存储介质存储在盘存储设备96中。计算机***88可以从盘存储设备96中检索恰当的数据以根据对应于本文描述的各种技术的实施方式的程序指令来处理数据。所述程序指令可以用诸如C++、Java等的计算机编程语言写成。所述程序指令可以存储在诸如程序盘存储设备86等的计算机可读介质中。这种计算机可读介质可以包括计算机存储介质和通信介质。计算机存储介质可以包括以用于信息的存储的任何方法或技术实施的易失存储器和非易失存储器以及可移动的和不可移动的介质,所述信息是诸如计算机可读指令、地震数据、结构、程序模块或其他数据等。计算机存储介质还可包括RAM、ROM、可擦可编程只读存储器(EPROM)、电可擦可编程只读存储器(EEPROM)、闪存或其它固态存储器技术、CD-ROM、数字通用光盘(DVD)或其它光存储、磁带盒、磁带、磁盘存储或其它磁性存储设备、或其它任何可用以存储所需信息并可由计算***88访问的介质。通信介质可包含计算机可读指令、数据结构、程序模块或诸如载波或其它传输机构等的调制的数据信号中的其他数据,并可以包括任何信息传递介质。术语“调制的数据信号”可以意指信号具有一个或多个它的特性集或以这种方式改变以在所述信号中编码信息。作为示例且不是限制,通信介质可以包括诸如有线网络或直连连接等的有线介质,以及诸如声学、RF、红外和其他无线介质等的无线介质。以上任何介质的组合都可被包括在计算机可读介质的范围内。
在一种实施方式中,计算机***88可以包括诸如图形显示器90和键盘92等的图形用户界面(GUI)组件,所述图形用户界面组件可以包括定点设备(例如,鼠标、轨迹球等,未示出)以启动交互操作。所述GUI组件可用以显示数据和处理的数据结果,并且允许用户选择用于实施所述方法的多个方面的选项。计算机***88可以在盘存储器86中存储以上描述的方法的结果,用于后续使用和进一步分析。
计算机***88可以位于远离数据采集区域或处理设备(未绘出)的数据中心中。计算机***88可以与所述数据采集接收器(直接地或通过记录单元,未示出)通信,以接收指示感兴趣的地下区域的地球物理特性的信号。这些信号在传统的格式化和其它初始处理后,可由计算机***88存储在盘存储器96中作为数字数据用于以上述方式的后续检索和处理。虽然图11中显示盘存储器96直连连接到计算机***88,但也可预期盘存储器96通过局域网或远程访问的形式访问。而且,虽然盘存储设备86和96被示出为用于存储输入数据和分析结果的分离设备,但是盘存储设备86和96可以在单个盘驱动器(一起的或分离的)内实施,或以本领域技术人员参考本说明书会完全理解的任何其他传统方式。
尽管本发明已经在以上以替代实施例的方式描述,但是预期仍有其他替代、修改和应用将对已阅读了本公开的本领域技术人员变得清晰。例如,要理解,在可能的程度上本发明预期任何实施例的一个或多个特征都可与任何其他实施例的一个或多个特征组合。因此本公开旨在被认为是说明性的而不是限制,并且所附权利要求被解释为包括落入到本发明的真实精神和范围内的所有这种应用、替代、修改和实施例。
Claims (3)
1.一种用于最小化感兴趣的地下区域的地震图像中的伪像的计算机实施的方法,所述图像基于从记录的地震数据导出的数据波束集合以及至少部分地从与地下区域和所述地震图像有关的速度模型导出的建模的波束集合,所述方法包括:
确定所述建模的波束集合与数据波束集合在预定时移范围内的互相关函数;
执行包括反向投影和全局最小化的蒙特卡洛反演以更新所述速度模型;以及
基于更新的速度模型形成地下区域的地震图像,
其中所述蒙特卡洛反演包括:
使用所述互相关函数确定在所述预定时移范围内的概率分布函数;
使用所述概率分布函数确定在所述预定时移范围内的累积分布函数;
生成预定间隔内的随机数;
从所述预定时移范围内选择时移值以使得所述随机数等于所述累积分布函数;以及
使用所述时移值通过所述反向投影为所述速度模型计算速度修正,
其中所述全局最小化是模拟退火;以及
其中蒙特卡洛反演的温度参数和时间间隔参数在模拟退火的接连迭代内从初始值减小至预定最终值,由此减少了基于更新的速度模型形成的地震图像中由于周期跳跃或与相干噪声的假对齐而导致的伪像的存在。
2.根据权利要求1所述的方法,还包括:
计算慢度改变;
使用所述慢度改变来计算所述速度修正。
3.一种生成与地球的地下区域有关的地球模型和地震图像的计算机实施的方法,所述方法包括:
将数据波束集合存储在非暂态计算机存储介质中,所述数据波束集合从对所述地下区域的一部分进行采样的记录的地震数据生成;
使用包括一个或多个处理器的计算机***,所述一个或多个处理器被配置成与所述计算机存储介质通信并且执行被配置为执行操作的一个或多个计算机程序,所述操作包括:
迁移包括多个数据波束的所述数据波束集合和具有代表所述地下区域的初始速度模型的初始地球模型,以生成所述地下区域的地震图像;
将从所述地震图像和地球模型导出的建模的波束与所述数据波束集合内的对应数据波束互相关;
将对于测量的时移具有超过确定的阈值的互相关值的互相关的建模的波束与数据波束的波束对存储在非暂态计算机存储介质中,并且重复所述互相关和存储操作直到确定的数目的对应数据波束已经被建模;以及
执行所述测量的时移的包括反向投影和全局最小化的蒙特卡洛反演以生成具有更新的速度模型的更新的地球模型,其中所得到的更新的地球模型能够生成准确描述所述数据波束穿过所述地下区域的传播的更新的地震图像;
其中执行所述蒙特卡洛反演的步骤包括:
确定所述建模的波束集合与数据波束集合在预定时移范围内的互相关函数;
使用所述互相关函数确定在所述预定时移范围内的概率分布函数;
使用所述概率分布函数确定在所述预定时移范围内的累积分布函数;
生成预定间隔内的随机数;
从所述预定时移范围内选择时移值以使得所述随机数等于所述累积分布函数;以及
使用所述时移值通过所述反向投影为所述速度模型计算速度修正;
其中所述全局最小化是模拟退火;以及
其中蒙特卡洛反演的温度参数和时间间隔参数在模拟退火的接连迭代内从初始值减小至预定最终值,由此减少了通过产生的更新的速度模型生成的更新的地震图像中由于周期跳跃或与相干噪声的假对齐而导致的伪像的存在。
Applications Claiming Priority (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
US201361788951P | 2013-03-15 | 2013-03-15 | |
US61/788,951 | 2013-03-15 | ||
PCT/US2014/021545 WO2014149923A1 (en) | 2013-03-15 | 2014-03-07 | Beam inversion by monte carlo back projection |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104755961A CN104755961A (zh) | 2015-07-01 |
CN104755961B true CN104755961B (zh) | 2018-05-04 |
Family
ID=50349969
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201480002816.6A Active CN104755961B (zh) | 2013-03-15 | 2014-03-07 | 基于蒙特卡洛反向投影的波束反演 |
Country Status (6)
Country | Link |
---|---|
US (1) | US10360317B2 (zh) |
EP (1) | EP2972506B1 (zh) |
CN (1) | CN104755961B (zh) |
AU (1) | AU2014237711B2 (zh) |
CA (1) | CA2886798C (zh) |
WO (1) | WO2014149923A1 (zh) |
Families Citing this family (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US11175421B2 (en) * | 2014-01-10 | 2021-11-16 | Cgg Services Sas | Device and method for mitigating cycle-skipping in full waveform inversion |
US10365386B2 (en) | 2015-03-30 | 2019-07-30 | Chevron U.S.A. Inc. | System and method for salt surface updating via wavefield redatuming |
CN106154323B (zh) * | 2015-04-01 | 2018-08-17 | 中国石油化工股份有限公司 | 基于地震拓频处理的相控随机反演薄储层预测方法 |
US11182639B2 (en) | 2017-04-16 | 2021-11-23 | Facebook, Inc. | Systems and methods for provisioning content |
US10942286B2 (en) * | 2017-08-24 | 2021-03-09 | Chevron U.S.A. Inc. | System and method for image-domain full waveform inversion |
US11047999B2 (en) | 2019-01-10 | 2021-06-29 | Chevron U.S.A. Inc. | System and method for seismic imaging |
US11474267B2 (en) | 2020-06-11 | 2022-10-18 | China Petroleum & Chemical Corporation | Computer-implemented method and system employing compress-sensing model for migrating seismic-over-land cross-spreads |
CN111722273B (zh) * | 2020-06-12 | 2022-02-11 | 中国海洋大学 | 一种模拟退火虚反射压制的方法、海上地震勘探*** |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102508293A (zh) * | 2011-11-28 | 2012-06-20 | 中国石油大学(北京) | 一种叠前反演的薄层含油气性识别方法 |
CN102667529A (zh) * | 2009-10-27 | 2012-09-12 | 雪佛龙美国公司 | 使用射束层析成像进行地震成像和地层建模的方法和*** |
Family Cites Families (18)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5274605A (en) | 1992-06-26 | 1993-12-28 | Chevron Research And Technology Company | Depth migration method using Gaussian beams |
US5640368A (en) | 1993-07-26 | 1997-06-17 | Exxon Production Research Company | Migration velocity analysis using limited-aperture and monte carlo migration |
US5570321A (en) * | 1994-03-03 | 1996-10-29 | Atlantic Richfield Company | Seismic velocity model optimization method using simulated annearling to determine prestack travel-times |
US5508914A (en) * | 1994-09-28 | 1996-04-16 | Western Atlas International, Inc. | Method for calculating static corrections for seismic data |
US6754588B2 (en) * | 1999-01-29 | 2004-06-22 | Platte River Associates, Inc. | Method of predicting three-dimensional stratigraphy using inverse optimization techniques |
US6246963B1 (en) * | 1999-01-29 | 2001-06-12 | Timothy A. Cross | Method for predicting stratigraphy |
US7523024B2 (en) * | 2002-05-17 | 2009-04-21 | Schlumberger Technology Corporation | Modeling geologic objects in faulted formations |
US7480206B2 (en) * | 2004-09-13 | 2009-01-20 | Chevron U.S.A. Inc. | Methods for earth modeling and seismic imaging using interactive and selective updating |
AU2006235820B2 (en) * | 2005-11-04 | 2008-10-23 | Westerngeco Seismic Holdings Limited | 3D pre-stack full waveform inversion |
WO2010080366A1 (en) * | 2009-01-09 | 2010-07-15 | Exxonmobil Upstream Research Company | Hydrocarbon detection with passive seismic data |
US8729903B2 (en) * | 2009-11-09 | 2014-05-20 | Exxonmobil Upstream Research Company | Method for remote identification and characterization of hydrocarbon source rocks using seismic and electromagnetic geophysical data |
US8537638B2 (en) * | 2010-02-10 | 2013-09-17 | Exxonmobil Upstream Research Company | Methods for subsurface parameter estimation in full wavefield inversion and reverse-time migration |
US8756042B2 (en) * | 2010-05-19 | 2014-06-17 | Exxonmobile Upstream Research Company | Method and system for checkpointing during simulations |
AU2011337143B2 (en) * | 2010-12-01 | 2016-09-29 | Exxonmobil Upstream Research Company | Simultaneous source inversion for marine streamer data with cross-correlation objective function |
US8862405B2 (en) * | 2011-12-06 | 2014-10-14 | Schlumberger Technology Corporation | System and method for producing look-ahead profile measurements in a drilling operation |
CA2864807A1 (en) * | 2012-02-17 | 2013-08-22 | Andrei I. Davydychev | Inversion-based calibration of downhole electromagnetic tools |
EP3101225B1 (en) * | 2014-10-28 | 2022-06-29 | Services Pétroliers Schlumberger | Integrated interpretation of pressure and rate transients for production forecasting |
US10197704B2 (en) * | 2014-12-19 | 2019-02-05 | Baker Hughes, A Ge Company, Llc | Corrective scaling of interpreted fractures based on the microseismic detection range bias correction |
-
2014
- 2014-03-07 CA CA2886798A patent/CA2886798C/en active Active
- 2014-03-07 EP EP14712524.9A patent/EP2972506B1/en active Active
- 2014-03-07 CN CN201480002816.6A patent/CN104755961B/zh active Active
- 2014-03-07 US US14/200,308 patent/US10360317B2/en active Active
- 2014-03-07 WO PCT/US2014/021545 patent/WO2014149923A1/en active Application Filing
- 2014-03-07 AU AU2014237711A patent/AU2014237711B2/en active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102667529A (zh) * | 2009-10-27 | 2012-09-12 | 雪佛龙美国公司 | 使用射束层析成像进行地震成像和地层建模的方法和*** |
CN102508293A (zh) * | 2011-11-28 | 2012-06-20 | 中国石油大学(北京) | 一种叠前反演的薄层含油气性识别方法 |
Non-Patent Citations (3)
Title |
---|
地球物理资料非线性反演方法讲座(二)蒙特卡洛法;王家映;《工程地球物理学报》;20070430;第4卷(第2期);第81-85页 * |
基于模拟退火算法的地质统计学反演方法;孙思敏 等;《石油地球物理勘探》;20070228;第42卷(第1期);第38-43页 * |
基于量子蒙特卡罗的地球物理反演方法;魏超 等;《地球物理学报》;20080930;第51卷(第5期);第1494-1502页 * |
Also Published As
Publication number | Publication date |
---|---|
WO2014149923A1 (en) | 2014-09-25 |
CN104755961A (zh) | 2015-07-01 |
AU2014237711B2 (en) | 2017-06-01 |
EP2972506B1 (en) | 2020-04-22 |
US10360317B2 (en) | 2019-07-23 |
EP2972506A1 (en) | 2016-01-20 |
AU2014237711A1 (en) | 2015-03-26 |
US20140278299A1 (en) | 2014-09-18 |
CA2886798C (en) | 2022-08-30 |
CA2886798A1 (en) | 2014-09-25 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104755961B (zh) | 基于蒙特卡洛反向投影的波束反演 | |
US9075163B2 (en) | Interferometric seismic data processing | |
CN102667529B (zh) | 使用射束层析成像进行地震成像和地层建模的方法和*** | |
CN101688926B (zh) | 基于拉普拉斯域波形反演的地下结构成像装置和分析方法 | |
US8737165B2 (en) | Interferometric seismic data processing for a towed marine survey | |
CN108369289A (zh) | 使用全波场反演点扩展函数分析设计地球物理勘测的方法 | |
US20170115418A1 (en) | Seismic data processing using matching filter based cost function optimization | |
US9268776B2 (en) | Methods and apparatus for data collection | |
US9575205B2 (en) | Uncertainty-based frequency-selected inversion of electromagnetic geophysical data | |
US20110199858A1 (en) | Estimating internal multiples in seismic data | |
RU2631407C1 (ru) | Способ и устройство для обработки сейсмических сигналов | |
CN105474048B (zh) | 利用波束分解预测地震数据中的层间多次波 | |
CN103119472B (zh) | 利用同时和顺序源方法进行全波形反演的混合方法 | |
EP2166379A2 (en) | Interbed seismic multiple prediction | |
Herrmann et al. | Frugal full-waveform inversion: From theory to a practical algorithm | |
US20150066458A1 (en) | Providing an objective function based on variation in predicted data | |
CN111060961B (zh) | 基于多信息约束反演的品质因子确定方法、装置及*** | |
Liu et al. | Tomographic velocity model building of the near surface with velocity-inversion interfaces: A test using the Yilmaz model | |
EA030770B1 (ru) | Система и способ адаптивной сейсмической оптики | |
CN108508481B (zh) | 一种纵波转换波地震数据时间匹配的方法、装置及*** | |
CN108594300B (zh) | 地貌成像方法、装置及计算机存储介质 | |
Wu et al. | Shot repetition: An alternative seismic blending code in marine acquisition | |
CN111965705B (zh) | 地震单炮记录地质层位的标定方法、装置、设备及介质 | |
Dong et al. | Automatic migration velocity estimation for prestack time migration | |
Yildirim et al. | Machine learning-enabled traveltime inversion based on the horizontal source-location perturbation |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |