CN1503006A - 形成地下区域中代表一物理量分布的模型的方法 - Google Patents

形成地下区域中代表一物理量分布的模型的方法 Download PDF

Info

Publication number
CN1503006A
CN1503006A CNA031530249A CN03153024A CN1503006A CN 1503006 A CN1503006 A CN 1503006A CN A031530249 A CNA031530249 A CN A031530249A CN 03153024 A CN03153024 A CN 03153024A CN 1503006 A CN1503006 A CN 1503006A
Authority
CN
China
Prior art keywords
noise
data
model
function
space
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.)
Granted
Application number
CNA031530249A
Other languages
English (en)
Other versions
CN1333267C (zh
Inventor
P・莱利
P·莱利
F·雷纳德
普拉-詹诺德
L·佩莱
F·德尔普拉-詹诺德
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.)
IFP Energies Nouvelles IFPEN
Original Assignee
IFP Energies Nouvelles IFPEN
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 IFP Energies Nouvelles IFPEN filed Critical IFP Energies Nouvelles IFPEN
Publication of CN1503006A publication Critical patent/CN1503006A/zh
Application granted granted Critical
Publication of CN1333267C publication Critical patent/CN1333267C/zh
Anticipated expiration legal-status Critical
Expired - Fee Related 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
    • 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/34Displaying seismic recordings or visualisation of seismic data or attributes
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/30Noise handling
    • G01V2210/32Noise reduction

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

一种从不均匀区域(如地下区域)探测得到的数据来形成该区中至少一个物理量分布的有代表性的模型的方法,所述模型不存在数据中可能含有的相关噪声。该方法包括下列步骤:获取给定有关该区特性的信息的数据,说明模拟算符,它将模型的响应与模型相关联(合成数据),通过对噪声发生函数(n、g、f)施加特定模拟算符模拟各相关噪声,说明数据空间中的半平均数,说明各n、g、f空间中的平均数,各噪声模拟算符对平均数构成等比例关系,以及借助采取这些等比例关系特性优点的算法寻求模型和使成本函数最小化的n、g、f,成本函数通过半平均数来定量在测量数据与模型的响应和与n、g、f有关联的噪声的叠加之间的差分,应用:例如寻求地下区域中声阻抗、传播速度、磁导率等的分布。

Description

形成地下区域中代表一物理量分布的模型的方法
发明领域
本发明涉及从不均匀介质区探测得到的数据形成该区中物理量分布的有代表性的模型的方法,所述模型(至少部分)不存在数据中可能含有的相关噪声。
该方法运用如地下区域中的声阻抗的定量法。
背景技术
该方法在于寻找在几乎所有科技领域业已开发用来调整实验测量的模型。这种方法称为下面的种种名称:对反演问题解法的参数估计最小2乘方法。对这一方法在球科学范围内的较好表示,可参考:
——Tarantola,A.:《Inverse Problem Theory:Methed for Data Fitting andModel Parameter Estimation》,Elsevier,Amsterdam,1987。
应该指出,术语“最小2乘方”是指数据空间中用来定量在模型的响应(它是由以前选定的模拟算符的模型的映像)与数据之间差异的平均数的平方,一种必须使其最小以解决问题的成本函数。用平均数的平方来定义成本函数只是实用上的方便,从根本上说并不是本质的。此外,许多作者由于种种理由采用成本函数的另一种定义,但这种定义仍基于数据空间中平均数或半平均数的应用。最后,我们有较大自由度选用数据空间中的平均数(或半平均数)(我们决定被迫使用欧几里德平均数)。在含有噪声的数据的场合,解决大体上决定于这一阶段作出的选择。对有关这个问题的较多进层可参考下列出版物:
-—Tarantola,A.:《Inverse Problem Theory:Methed for Data Fittingand Model Parameter Estimation》,Elsevier,Amsterdam,1987;Renard andLailly,2001;Scales and Gersztenkorn,1988;Al-Chalabu,1992。
测量的结果常含误差。当将实验与模拟结果比较时,模拟噪声进一步加到这些测量误差上:模拟从来不是理想的,因此总是对应于现实的简化模式。此后,噪声将被描述为下述的叠加:
——非相关成分(如白噪声),
——相关成分,指测量样本上的噪声存在转移到某些相邻测量点上相同特性的噪声存在;模拟噪声一般属于这一类。
当数据含有相关噪声时,用解反演问题来评估的模型的质量从而受到严重影响。如已提及,没有理想的模拟算符。因此,整个人类社会的工作被卷入描述模型的参数(被相关噪声的存在所妨碍)的识别中。在这些人中,地震探测法从事者是在最有关的人群之中:事实上,他们的数据具有较差的或甚至很差的信噪比。这就是为什么相关噪声滤波技术是地震探测法数据处理软件的一个重要部分的理由。最传统的技术采用变换(如傅里叶变换),其中信号与噪声被置于不同的空间区域中,从而使信号与噪声分离。对地震探测数据中消除噪声的传统方法的一般表示,可参考Yilmaz的书:
——Yilmaz,0.1987:《Seismic data processing》,地球物理研究No.2,地球物理探测学会,Tulsa,1987。
然而,这些技术并不理想:他们预先假设使信号与噪声完全分离的变换业已找到。对于这一点,由于相关噪声难于从信号(也是被相关的)中分离且幅度大,故特别麻烦。因此经常不易处理下述的折衷:要末信号保留了但大的噪声仍在,要末噪声被消除但信号被失真。滤波技术可以在解反演问题之前实施,于是它们构成对数据施加预处理。于是反演问题解法的质量广泛地依赖于滤波器消除噪声而不使信号失真的能力。
Nemeth等人在下列书中
——Nemth T.et al.(1999)《Least-Square Migration of IncompleteReflection Data》,Gecphysics,64,208-221,引入的方法构成了对受大幅值相关噪声干扰的数据的反演的重要进步。这些作者提出通过解反演问题来消除噪声:通过线性算符T将相关噪声的空间定义为矢量空间B(噪声发生函数的空间)的映像空间,并通过模拟算符F(假设为线性)将信号定义为矢量空间M(模型空间)的映像空间(的信号?)后,他们发现F(M)中的信号与T(B)中的噪声的和尽可能地接近(在平均数L2的欧几里得平均数或连续基础上的意义上)所测的数据。它使得用传统的滤波技术难以从信号中消除(或至少实质上减少)大幅度值相关噪声(如表面波)这一限度内,这一技术构成一个重要的进步。然而,根据这些作者的方法,为了该传统的反演问题的解即相关的噪声分量不明显,为这些性能所付的代价是所需计算时间的数量级的增加。而且,用该方法得到的结果对由于算符T的定义引入的任何不精确性极为敏感,这是为该方法区分信号与噪声的高性能的不可避免的补偿。
发明概要
根据本发明的方法,能从以下区域探测得到的数据来估计该区域至少一个物理量分布的有代表性的模型,这个模型不存在数据中可能含有的相关噪声。该方法基本上包括下列步骤:
a)根据预定的实验协议,获取对给定有关该区的某物理特性的信息的测量值,
b)给出将构成模型的响应的合成数据与各物理量的模型相关联的模拟算符的详细说明,所述测量值和所述合成数据属于数据空间,
c)对标有从1至J的下标j的各相关噪声,选择模拟算符,所述算符将相关噪声与属于预定的噪声发生函数空间(Bj)的噪声发生函数相关联,
d)给出所述数据空间中平均数或半平均数的详细说明,
e)给出上述噪声发生函数空间中半平均数的详细说明,各噪声模拟算符对在该噪声发生函数空间和数据空间之间建立大体上的等比例的关系,
f)定义-成本函数,所述成本函数定量在一方为测量值和另一方为模型响应和与上述噪声发生函数有关的相关噪声的叠加之间的差异,以及
g)通过采取噪声模拟算符的等比例特性的优点的算法,使所述成本函数最小来调整所述模型和所述噪声发生函数。
按照第1实施模式,响应于局部地震激励,寻找作为介质中声阻抗的深度的函数的分布,影响所述数据的相关噪声是通过表征其传播参数每一个识别的管道波,所述测得的数据是通过适合于检测介质中质点位移的传感器得到的VSP数据,确定所述传感器的位置、记录时间以及时间采样点,并且所述的模拟算符使所述合成数据与声阻抗分布相关联作为在传播时间中估计深度的函数,并与测得的垂直应力相关联作为在第一传感器深度上时间的函数。
所选用来定量在一方为测量值与另一方为模型响应和与噪声发出函数关联的相关噪声的叠加之间的差异的成本函数是例如数据空间中这一差异的半平均数的平方。
按照一实施例,通过块驰豫法来消除对应于各相关的噪声发生函数的未知数,得到模型和噪声发生函数的调整,所述驰豫法在准牛顿算法的迭代中实施用于计算该模型。
按照一实施例,通过选择在传感器位置上和以前确定的时间采样点上的质点位移得到的值,并通过施加算符适当补偿球面发散和衰减效应,用所考虑模型的一维波方程的数字解,实现用模拟算符对模型映像的数字计算。
按照一实施例,所述数字噪声模拟算符是一使噪声适移方程离散化的有限差分对中数字方案,并且包括作为沿所述区的边缘初始条件的噪声发生函数使各相关噪声与在定时间间隔内由支持时间函数组成的噪声发生函数空间相关联。
为数据空间和噪声发生函数空间所选的各平均数或半平均数例子将在以后详述。
按照第2实施模式,寻求与以前所选的基准介质中有关的阻抗扰动的分布和所述介质区域速度的分布,影响数据的相关噪声是由于多重反射引起的,其随抵消而变化的动能与幅值业已估计,所测的数据用地震检测表面传感器检拾,传感器的位置、地震激励模式、记录时间以及时间采样点均被确定,模拟算符通过基准基质附近使波动方程线性化加以确定。
选择用来定时在一方的测量与另一方模型响应和与噪声发生函数有关联的相关噪声的叠加之间的差的成本函数是例如数据空间中该差的半平均数的平方。
按照一实施例,通过块驰豫法来消除对应于各相关噪声发生函数的未知数,得到模型和噪声-发生函数的调整,该块驰豫法在共轭梯度算法的失代中实施用于计算该模型。
附图说明
通过结合附图阅读下面给出的非限制性实施例的说明,根据本发明的方法的其他特点和优点将显而易见,附图中,
图1示出通过垂直地震探测法(VSP)得到的数据例,
图2示出声阻抗对深度的分布模型,
图3示出基于图2的阻抗模型得到的合成的VSP数据,
图4示出被单个相关噪声混杂的VSP数据例,
图5示出被两个相关噪声混杂的VSP数据例,
图6示出通过图4的含噪声数据的反演得到的作为深度的函数的声阻抗分布,
图7示出反演的剩余(图4的数据与图6的模型的地震响应之间的差),
图8示出通过图5的含噪声数据的反演得到的作为深度的函数的声阻抗分布,
图9示出对应的反演剩余(图5的数据与图8的模型的地震响应之间的差),
图10示出通过寻求两个具有不正确传播速度即不同于图4数据中出现的噪声的传播速度的相关噪声的叠加形式的相关噪声对图4的含噪声数据反演得到的作为深度的函数的阻抗分布。
图11示出对应的反演剩余(图4的数据与图10的模型的地震响应之间的差),
图12A、12B分别示出在传统的线性性的反演后得到的有多次反射的地震数据和模型的地震响应,
图13A、13B分别示出阻抗分布和传播速度的例子,图12A的地震数据对此构成地震响应,及
图14A、图14B分别示出用传统反演(14A)和用所提议方法反演(14B)得到的模型地震响应的比较,所提议的方法包括除去多次反射的检拾和幅值随抵消而变化的估计。
具体实施方法
问题的表述
我们拥有从一函数的采样测量得到的数据。此函数依变于几个变数(如空间或空间-时间变量)。
这些数据(称为d)对应于已执行的测量以获得有关模型m的信息。他们包含各种噪声:
-相关的噪声(或相关噪声的叠加),
-非相关的噪声。
问题是从数据d定量地确定模型m(或这个模型的函数)。
因此我们选择一将模型m与该模型的响应相关联的模拟算符F(线性或否)。这个算符实际上只不完美地模拟该真实的物理现象。这就是为什么相关的噪声出现在数据中的主要(但不是全部)原因:这些噪声相当于与模型有关的信号,但它们之间的关系显得太复杂以致不包括在由于各种原因使我们想保持相对简单的模拟算符F中。
在这种情况下,我们通过引入一或多个噪声-发生空间Bj(j从1到J)和称作Tj的有关的噪声模拟算符。为了从数据d寻找模型m,我们提出寻找这个模型作为对该最佳化问题的解。
min ( m , β 1 , β 2 , . . . β J ) ∈ M × ∏ j = 1 J B j C ( m , β 1 , β 2 , . . . β J ) = | | F ( m ) + Σ j = 1 J T j β j - d | | D 2 ( 1 )
其中‖‖D是数据空间中平均数(或可能为半平均数)。
上面的公式与Nemeth等人(1999)提出的方法的不同之处在于:
-模拟算符F不需要线性,
-我们可以以明显增加有关未知数(它使该最佳化问题的解特别复杂)数目和待执行的噪声模拟数目的计数次数为代价,考虑不同类型(以由下标j标明的不同算符模拟它们的意义上)的相关的噪声的叠加。
-我们保留选择平均数‖‖D或可能的半平均数的可能性以在数据空间适合我们的方便。
我们为‖‖D要作的选择将受到与解法的效率(即允许处理相关噪声叠加的情况的效率)有关的考虑的影响(下面说明)。这是一个重要的通道:除了有可能处理含有相关噪声与复杂特性的数据之外,我们可以这种方式来克服Nemeth等人的方法中遇到的对噪声模拟算符引起的不正确性高度敏感有关的困难,为克服这一困难,我们只得通过与近似的模拟算符Tj有关的噪声的叠加来寻找该噪声。
应该指出,涉及最佳化问题中的平均数‖‖D的平方的选择不是本质的,我们通过选择如平方函数那样是R+上增函数(或假如min变成max时的减函数)的替代函数并不改变该解。
解法
该方法在于适当地选择:
-数据空间中的平均数(或半均数)‖‖D
-各空间Bj中的平均数‖‖D
以获得对最佳化问题(1)的高性能算法解。“适当地选择”的表述意指注意到各算符Tj在起始空间和结束空间中分别对平均数‖‖j和‖‖D构成等比例关系(或近似等比例关系)。使各算符等比例的平均数的决定的主要思想是基于能量守恒型等式的存在,这是由偏微分方程的解所证明的(或者对用于使这些方程式离散化的某些数字方案的解来说存在离散的能量等式)。
如果我们能作这种选择,则通过消除与相关噪声发生函数有关的未知数有可能使成本函数最小化十分简单(确定Schur的补数)。可用合适的算法如块驰豫法(块与噪声-发生函数相关联)或具有对称块Gauss-Seidel先决条件的共轭梯度法等来实行这一消除,所有这些算法实施方法是本专业内的技术人员所熟知的,例如由下述著述提出的:
-Golub,G.H.et al.(1983)《Resolution numerique des grands systemeslineaires》(Eurolles)。
简化的使成本函数最小化可显著地减少数字解所需的时间。
适用于被管道波混杂了的VSP数据的1维反演的第1实施例
VSP(垂直地震检测分布)数据属于传统的较好勘察数据。它有多种应用,其中主要用途之一是用来在地震表面数据与记录测量之间建立联系。例如在
-Mace,D.and Lailly,P.(1984)A Solution of an Inverse problem with1 D Wave Equation applied to the Inversion of Vertical Seismic Profile;Proc of the 16 th Intern.Conf.In《Analysis and Optimisation of System》,Nice,2,309-323中说明的VSP数据的反演已为此而开发。在1维反演问题中,我们假设地面下模型是横向不变的并以平面波激励垂直传播。问题的未知数是作为深度(以垂直的传播时间并在第1个数据收器所在深度开始的整个范围内侧得的)的函数的声阻抗分布和该地震波激励模式(在第1个传感器深度上精确地表征边界条件的时间函数)。这是一个非线性反演问题,VSP数据与阻抗分布是非线性关系。
VSP数据经常被很大幅值的流体波的相关噪声所混杂。此波在侵入井体的混浆中传播。其传播速度相对于井周围的岩石的传播速度而言是低的。此外,此波的传播速度与频率有关(色散传播)。还有,此波在井底反射显著,而且产生另一种传播模式。典型的VSP数据在以下的图中给出:可以看出管道波的主要特征。其幅值和被混杂区的扩展使得难以利用含在被流体波所混杂的许多地震检测样值中的信号。
以下描述VSP数据的1维反演方法的实施。
实验协议
规定:
-接收器位于井(假设为垂直的)中的各深度:在我们的实验中,传感器覆盖以单程时间[Xmin=0.05s;Xmax=0.2s]测得的深度范围,每Δx=1ms(单程时间)处有一接收器;从而获得从0至I标出的样值,
-由这些传感器实施的测量值的物理特性:在我们的实验中,传感器测量波传播导致的垂直位移,
-记录时间:在我们的实验中,传感器测量在时间区[0,T=1.5s]中的振动状态,这些数据每Δt=1ms被采样1次,从而获得从0至N标出的样值。
我们称在深度Xmin+iΔx和时间nΔt上记录的数据样值为di n。当然有T=NΔt,Xmax=Xmin+IΔx。
按此实验协议的数据的获得
从合成数据即由数字模拟得到的数据进行我们的实验,分两步实现:
产生不含噪声的合成数据:用1维声学波动方程的数字解来实现这一产生,模型是图2所示的声阻抗分布(深度用单时间测定),激励模式(深度为0处的Neumann边界条件)是传统的Ricker子波。
用这一模拟得到的合成VSP数据如图3所示:垂直轴代表观察时间,从0至1.5s,水平轴代表各传感器的深度(单程传输时间),从0.05至0.2s。
数据上加噪声
我们将对如此计算得到的VSP数据加上噪声进行干扰,加上的噪声其一为幅值颇大的随机噪声(但在地震检测频率中过滤),其二为幅值很大的一或多种相关的噪声。各噪声以取决于其频率(色散传播)的低速度向下传播:这些相关噪声被假定为有代表性的流体波。它们已以恒定传播速度的传播方程的数字解(使用有限关分对中方案)所模拟。在有限差分方案中利用大的离散化区间得到波的色散。在几种相关噪声叠加的情况下,与各噪声有关的传播速度是不同的。图4、5表示由单个相关的噪声混杂(图4)的和由两个相关的噪声叠加(图5)的VSP数据。
模拟算符的说明
通过解波动方程来实现模拟:
σ ( x ) ∂ 2 y ∂ t 2 - ∂ ∂ x ( σ ( x ) ∂ y ∂ x ) = 0 在]xmin,+∞[×[0,T]中,
Neumann边界条件:
- σ ( x min ) ∂ y ∂ x ( x min , t ) = h ( t ) on [ 0 , T ]
以及初始条件:
y ( x , 0 ) = ∂ y ∂ t ( x , 0 ) = 0 在[xmin,+∞]上。
模拟算符的定义首先要求使§1中所述实验协议定形的观察算符的定义。因此观察算符是与函数y(x,t)、上述系统的解、样值y(xi,tn),xi=xmin+iΔx和tn=nΔt,i=0,…,I,n=0,…,N有关的算符。
因而我们称为F(m)的模拟算符是与函数对m=(σ(x),h(t))样值,y(xi,tn),xi=xmin+iΔx和tn=nΔt,i=0,…I,N=0,…,N有关的(非线性)算符。
与各相关的噪声有关的模拟算符的说明
本节中我们说明用于模拟各相关噪声的步骤。各相关噪声由假设已知或至少近似知道的对其相关方向特性来表征:通过部分(cj x(x,t),cj t(x,t))(在域[xmin,xmax]×[0,T]的任一点上需被确定),我们可以等效的方法通过代表与cj(x,t)有关的场线的相关线束来规定相关方向。按照申请人提交的专利EP-354,112(US-4,972,383)中所述的技术可以得到这种相关线,从采样噪声的一些特性通过传统的内插和/或外插步骤,信息被扩展至整个域[xmin,xmax]×[0,T]。相关方向的说明就是说明了以关系 c j ( x , t ) = c j x ( x , t ) / c j t ( x , t ) 与相关矢量关联的噪声cj(x,t)的传播速度分布。假设在波传播过程中波幅度值不变(不同情况在第2实施例中提出)。如上所述相关噪声的模拟基于下述的观察:波沿相关线传播而无幅值改变是下述传输方程的解:
▿ → b ( x , t ) . c → ( x , t ) = 0
对此,我们引入:
-噪声-发生函数的空间Bj,它将是支持函数βj(t)在[tj min,tj max][0,T]内的空间,我们用0在整个区间[0,T]中扩展,
-噪声模拟算符Tj
Tj:βj(t)∈Bj→bj(x,t)∈D
传输方程的解
▿ → b j ( x , t ) . c → j ( x , t ) = 0 在[xmin,xmax]×[0T]中
以及满足初始条件:
Figure A0315302400132
应该指出,相关线的束的几何形状(或说传播速度分布)不能是任意的几何形状:相关线不能相交,否则传输方程的解将不是正确的组。同样理由,相关线不能与被规定的初始条件的边缘交切两次。后一种考虑可导致(为模拟对应于上升波的相关噪声)规定初条件不再在x=xmin而在x=xmax,或甚至中间值上,即使它意味着将Cauchy问题分解成两个与位于该中间值边缘任一侧上的域有关的问题。
我们进一步假设相关线的几何形状为(考虑更一般情况,而后面的结果取不同形式):
              cj(x,t)=ξj(x)×τ(t)
ou encore : c j x ( x , t ) = ξ j ( x ) et c j t ( x , t ) = 1 τ ( t )
最后我们假设相关线的束在x=xmin区间[tj min,tj max]相交,而不在t=T区间[xmin,xmax]相交。
事实上,用数字方案来解传输方程。如可用传统的有限差分对中方案。因此,我们选择离散化区间Δxj和Δtj(它们不必等于Δx和Δt,但我们假设为它们的约数,即使可考虑更一般的情况),并引入一格栅,其节点是坐标(xmin+i′Δxj′,n′Δtj′),i′和n′∈N的点。传统的有限差分对中方案(以后说明)使能从初始条件来逐步计算各格栅节点上函数bj(x,t)所取的不同的值(为简化符号我们将略去与相关噪声有关的下标j)。
ξ i ′ + 1 2 2 Δ t ′ [ b i ′ + 1 n ′ + 1 + b i ′ n ′ + 1 - b i ′ + 1 n ′ - b i ′ n ′ ] +
1 2 Δ x ′ τ n ′ + 1 2 [ b i ′ + 1 n ′ + 1 + b i ′ + 1 n ′ - b i ′ n ′ + 1 - b i ′ n ′ ] = 0
上式中,量αi n′表示在坐标点(xmin+i’Δxj’,n’Δtj’)上量α(α是一般项)的估计。
当然,如我们要想该数字方案给出传输方程函数解的精确近似,则使用相对于波长是充分小的离散化区间Δxj′和Δtj′。
然而使用较大离散化区间(但仍然小于波长)允许模拟更复杂的、波传播速度依赖于即色散传播的传播现象,如管道波的模拟。选择适当的离散化区间用于模拟特定的管道波模式可以这样做:
检查数据后确定这种模式的传播速度(不仅决定于空间和时间,也决定于频率);对此使用下列刊物所载的层析技术:
-Ernst,F.et al.(2000);Tomography of Dispersive Mdeia;J.AcousticSoc.Am.,108,105-116。
知道与数字方案有关的色散传播特性后如下述刊物所述从数字色散关系的分析中得到特性:
-Alford,R.M.et al.(1974),Accuracy of Finite Difference Modelingof the Acoustic Wave Equation;Geophysics,6,834-842。
数据空间的平均数或半平均数的说明
首先取任意离散的平均数作为数据空间的平均数,这意味着求由下式在连续基础上确定的平均数的近似值(通过求积公式):
| | u | | D = ( ∫ x min x max dx ∫ 0 T dt 1 τ ( t ) u 2 ( x , t ) ) 1 2
我们也可使用半平均数,从下节可见这是较好的选择:
| | u | | D = ( ΔxΔt Σ i = 0 I - 1 Σ n = 0 N - 1 1 τ n + 1 2 ( u i n + 1 + u i n ) 2 ) 1 2
式中u为数据空间的任意矢量(i和n为分别表示空间和时间中样值数目的索引)。
此外,其他选择也是可能的。
在噪声-发生函数的各空间Bj中各噪声模拟算符Tj为之在噪声-发生函数空间和数据空间(有§5中定义的半平均数)之间建立等比例关系(或近似等比例关系)的平均数的说明
本节中我们说明在与各相关噪声有关联的噪声-发生函数空间中用来定义平均数的步骤。这个步骤对各噪声都相同,我们一般地加以说明并在公式中略去与所述噪声有关的下标j。
首先选有离散平均数的各噪声-发生函数空间B,离散平均数就是指由下式在连续基础上定义的平均数的近似(通过求积公式):
| | β | | B = ( X ∫ 0 T dt 1 τ ( t ) β 2 ( t ) ) 1 2
在这种情况下,利用在§4提出的假设,我们有bi(x,t)=0,x∈[xmin,xmax],线性的算符Tj实现Bj和D之间的等比例关系。然而,由于Tjβj的计算是以数字方式实现的,算符T的等比例性质并不是理想的但只是近似而已,当离散化间隔Δx,Δt,Δxj′,Δtj′小时,近似更好。仅有近似上的等比例的事实将导致有关在后面第8节提出的算法实现的较低性能。
较好的选择是有半平均数的各噪声-发生函数空间Bj
| | β | | B = ( ΔxΔtI Σ n = 0 N - 1 1 τ n + 1 2 ( β n + 1 + β n ) 2 ) 1 2
在这种情况下,算符Tj严格地在Bj和D之间实现等比例,至少在ui N=0Ai=0,I-1时如此。从下式离散的能量等式得到结果:
I Σ n = 0 N - 1 1 τ n + 1 / 2 [ u o n + 1 + u o n ) 2 ] = Σ i = 0 I - 1 Σ n = 0 N - 1 1 τ n + 1 / 2 [ ( u i n + 1 + u i n ) 2 ]
+ Σ i = 0 I - 1 Σ k = 0 i Δx ξ k + 1 / 2 Δt [ ( u k + 1 N + u k N ) 2 ] - IΔx ( u o N ) 2 ξ - 1 / 2 Δt
在Δx′=Δx、Δt′=Δt并定义ξ-1/2=ξ1/2的情况下等式有效。
成本函数的说明
我们通过公式定义成本函数如下:
C ( m , β 1 , β 2 , . . . β J ) = | | F ( m ) + Σ j = 1 J T j β j - d | | D 2
寻求模型和使成本函数最小的噪声-发生函数,通过采取噪声模拟算符的等比例特性的显性中隐性优点的算法来实现这一搜索(见ξ6)。
该算法利用Schur补数的特性,一种与算符Tj的等比例特性有关联的特性。可通过实现使C(m,β1,β2,…βJ)最小化相当于使 C ′ ( m ) = C ( m , β ~ 1 ( m ) , β ~ 2 ( m ) , · · · , β ~ J ( m ) 最小化(其中对给定的 m , β ~ 1 ( m ) , β ~ 2 ( m ) , · · · , β ~ J ( m ) 使C(m,β1,β2,…βJ)最小化)来达到。应该注意,与这一二次式有关的Hession是Schur的补数。C′(m)可通过任何最佳化方法如准牛顿法的BFGS版本来最小化。显然,如果分别选中了选择(2)和(3)来确定‖‖D和‖‖Bj由于算符Tj的等比例特性, ( β ~ 1 ( m ) , β ~ 2 ( m ) , · · · , β ~ J ( m ) ) 的决定非常容易。如果只存在一种类型的噪声(J=1)以及等比例性是理想的,则
Figure A0315302400162
的决定只要求估计该二次式的梯度。如存在几种噪声的叠加(J>1),可考虑各种算法以打取算符Tj的等比例特性的优点。可通过例如块驰豫法(Gauss-Seidel混淆)来决定 β ~ 1 ( m ) , β ~ 2 ( m ) , · · · , β ~ J ( m ) , 该法将块与一组未知数
Figure A0315302400164
相关联,迭弛法只要求(仍在理想的等比例情情况下)估计与所述块有关的二次式的梯度。对于使得难以接近各种噪声类型的相关方向的问题,可考虑用预处理的共轭梯度法,该预处理矩阵具有对称块(弛豫混淆)Gauss-Seidel型式。因此,预处理再一次是非常容易实现的,与块相关的问题的解只要求估计与所述块有的二次式的梯度。
图6至11示出了用该方法的三类实验得到的结果。
在既有相关噪声又有随机噪声的叠加的含噪声数据(图4)的反演的情况下,图6给出所得到的阻抗分布,图7示出反演剩余。
在既有相关噪声又有随机噪声的叠加的含噪声数据(图5)的反演情况下,图8给出所得到的阻抗分布,图9示出反演的剩余。
在既有相关的噪声又有随机噪声的叠加的含噪声数据(图4)的反演,其中流体波的传播特性并不精确知晓的情况下,自然要寻求包括两相关噪声的叠加的相关噪声模型,所述两相关噪声具有不准确的传播速度,但通过它构建真实的流体波的传播速度。图10验出所得的阻抗特性,图11示出反演剩余。
通过线性反演由多次反射所混杂的多低消数据来确定模型的方法适用的 第2实施例
本应用不同于第1例的主要在于考虑到噪声幅值沿相关方向变化时说明本方法提供的可能性。另一个差别在于相关噪声的传播性质,多次反射的传播没有色散。
线性化反演是地球物理工作者所用的一种先进的方法,来获得地下不均匀性的定量模型,这些信息对油层特性是明显地宝贵的。这里给出的数据是表面地震检测数据。例如如下的刊物说明了这一方法:
-Bourgeois,A.et al.(1989)《The Linearized Seismic Inverse Problem:an Attractive Method for a Sharp Descriptipn of Strategic Traps》:59th Anm.Mtg.Soc.Expl.Geoph.Expanded Abstract,pp.793-976。
此外,需要有一参考模型,包括(在声学模型方面)速度分布和声阻抗分布:这是一个在其周围波动方程中可以实现线性化(Bornr的近似)的模型。为了说明方便,我们取1维情况,但由于实际应用主要涉及3维地震检测方法,所以这种假设较少意义。在1维范围内,所有地震检测信息包含在单个射击点内。此外,参考模型也是1维的:只与深度有关。问题的未知数是作为深度的函数的声阻抗的扰动与速度分布。另一方面,与第1实施例不同,假设地震激励模式(以及明显的地震子波)是知道的。这是一个通过线性化使成为线性的反演问题。
表面地震数据经常受到幅值大的多次反射的混杂。在1维范围内用扩展来表征的运动学特征通常十分不同于初次反射的:这两种波通常穿过具有不同传播速度的地质层。这就是为什么线性的反演对衰减多次反射是一种有效技术的原因:主反射的运动学(由Born近似模拟的事件)完全由参考模型内速度分布所确定,所以如果多次反射的扩散不同于主要反射的,则模拟的事件不能够调整多次反射。然而,由于零抵消点的固定性,所以在实际上存在部分调整,以及线性反演结果出现被包括多次反射的相关噪声的存在所混杂。图12A的地震检测数据是已考虑了多次反射(只有3次主反射,相继在500、1000、1650ms的时刻到达零抵消,第2次到达具有特高的振幅)的模拟。
在传统的线性化的反演之后得到的模型的地震响应中(见图12B),在1500ms附近到达零抵消的高振幅只是微弱的衰减,并且模型被多次反射所混杂。于是自然地试用Nemeth方法提供的可能性以地更好地区别多次的和主要的反射,显然具有我们的方法提供的改进。
以下说明该方法对表面的震检测数据器的1维线性化反演的实施。
实验协议
规定:
-地震激励模式:在我们的实验中,地震源是用函数w(t)(地震子波)时间调制的源点,
-接收器所在的各位置:在我们的实验中,安排传感器在深度10m并覆盖抵消区间[xmin=0,xmax=1500m],接收器每Δx=100m;得到编号从0至I的轨迹,
-传感器测量从波传播产生的压力场,
-记录时间:在我们的实验中,传感器测量在时间区[0,T=1800ms]的振动状态,这些数据在每△t=1ms时被采样。得到编号从0至N的样值。
-我们称对抵消xmin+i△x和在时间n△t时记录的数据样值为di n。当然有T=N△t,xmax=xmin+I△x。
按照实验协议的数据的获得
从合成数据实现我们的实验:这些数据通过下面2维波动方程的数字解(有限差分方法)得到:
1 σc ∂ 2 P ∂ t 2 - ▿ σ c ▿ → P = δ ( x , z - z s ) w ( t )
Figure A0315302400182
边界条件:
    P(z=0)=0
初始条件
P ( t = 0 ) = ∂ P ∂ t ( t = 0 ) = 0
方程中:
-x,z,t分别为横向坐标、深度和时间,
-σ和c分别为声阻抗和传播速度分布(深度的函数)。
所选的子波是中心频率为30Hz的常用Ricker子波,阻抗σ和速度c分布在图13A、13B给出。所以该模型的地震检测响应包括图12A包括的数据。将地震数据与函数对(σ(z),c(z))相关联的算符称为FNL
模拟算符的说明
我们选择由函数对(σo(z),co(z))描述的参考模型。在以后提出的实验中,我们已选择co(z)=c(z)和σo(z)=cte(=1)。我们选在点(σo(z),co(z))上的算符FNL的Jacobian算符作为模拟算符。事实上,由于已选Co(z)=c(z),故唯一未知数是δσ(z)=σ(z)-σo(z)以及唯一对应的Jacobian分量是重要的。
因此我作为F(δm)的模拟算符是(线性)算符,它将样值δy(xi,tn)与函数对δm=(δσ(z0,δc(z)相关联,样值δy(xi,tn)是矢量F(δm),xi=xmin+i△x和tn=n△t,i=0,…,I;N=0,…,N的分量。
与各相关噪声有关的模拟算符的说明
本节中我们说明用于模拟各相关噪声的步骤。各相关噪声由假设已知或至少近似知道的对其相关方向的特性来表征。通过分量(cj x(x,t),cj t(x,t))(在域[xmin,xmax]×[0,T]的任一点上需被确定)的相关矢量
Figure A0315302400191
的场规定这些相关方向。我们可以等效的方法通过代表与 有关的场线的相关线束来规定相关方向。按照申请人提交的专利EP-354,112(US-4,972,383)中所述的技术可以得到这种相关线,从采集噪声的一些特性通过传统的内插和/或外插步骤,如申请人提交的上述专利Ep-354,112所述,信息被扩展至整个域[xmin,xmax]×[0,T]。对于要消除上述的多次反射的第2应用例,我们采集这一多次反射(如幅值峰值),它能使我们确定作为抵消点的函数的到达时间改变,通过将所采集线的简单垂直移动确定相关线的连续区(在这些条件下,矢量 不决定于t)。相关方向的说明就是说明了以关系 c j ( x , t ) = c j x ( x , t ) / c j t ( x , t ) 与相关矢量关联的噪声cj(x,t)的传播速度分布。
与上述第1应用例不一样,这里噪声幅值沿相关线改变。假设我们能根据数据的测量或据理论上的考虑来估计这些幅值改变,则这一测量确定了函数g(x)。
因而通过两个算符的合成;噪声得以模拟:
传输算符(如第1应用例中);
幅值调制,即以函数g(x)乘以传输方程的解。
如果我们再一次利用第1应用例的方案(显然如使用相同的假设和符号),则我们引入:
-噪声发生函数的空间Bj,它将是支持函数βj(t)在[ttj min,tj max][0,T]中的空间,我们用0在整个区间[0,T]中扩展;
噪声模拟算符Tj:
      Tj∶βj(t)∈Bj→bj(x,t)∈D
传输方程的解
▿ → b j ( x , t ) . c → j ( x , t ) = 0 在[xmin,xmax]×[0,T]
并满足初始条件:
事实上,用数字方案解传输方程。如利用传统的有限差分对中方案是可能的。所以我们选择离散化间隔Δxj′和Δtj′(不必等于Δx和Δt,但我们假定为这些量的约数,尽管可考虑更一般的情况),并引入一格栅,其节点是坐标(xmin+i′Δxj,n′Δt′j),i′和n′∈N的点。下面说明的传统有限差分对中方案使可根据初始条件逐步计算函数bj(x,t)在各格栅节点所取的不同的值(为简化符号,略去与所述相关噪声有关的下标j)。
ξ i ′ + 1 2 2 Δ t ′ [ b i ′ + 1 n ′ + 1 + b i ′ n ′ + 1 - b i ′ + 1 n ′ - b i ′ n ′ ] +
1 2 Δ x ′ τ n ′ + 1 2 [ b i ′ + 1 n ′ + 1 + b i ′ + 1 n ′ - b i ′ n ′ + 1 - b i ′ n ′ ] = 0
上式中,量αi′ n′代表在坐标(xmin+i′Δxj′+n′Δtj′)处量α(α是一般项)的估值。
这里,必须选择离散化间隔Δxj′和Δtj′相对于波长为足够小,因为我们要求给出传输方程的函数解的精确近似的数字方案(我们要求不色散)。
噪声模拟步骤的最后一步是选择值bi′ n′它属于两个格栅共有的节点并以g(iΔx)与它们相乘:从而我们得到矢量(数据空间的)Tjj)的各分量。
数据空间中平均数或半平均数的说明
首先取任意离散的平均数作为数据空间的平均数,这意味着求由下式在连续基础上确定的平均数的近似值(通过求积公式):
| | u | | D = ( ∫ x min x max dx ∫ 0 T dt 1 τ ( t ) u 2 ( x , t ) ) 1 2
我们也可使用半平均数,从下节可见这是较好的选择:
| | u | | D = ( ΔxΔt Σ i = 0 I - 1 Σ n = 0 N - 1 1 τ n + 1 2 ( u i n + 1 + u i n ) 2 ) 1 2
式中u为数据空间的任意矢量(i和n为分别表示空间和时间中样值数目的索引)。
此外,其他选择也是可能的。
在噪声-发生函数的各空间Bj中各噪声模拟算在Tj为之在噪声-发生函数空间和数据空间(有§5中定义的半平均数)之间建立等比例关系(或近似等比例关系)的平均数的说明。
本节中我们说明在与各相关噪声有关联的噪声-发生函数空间中用来定义平均数的步骤。
首先选有离散平均数的各噪声-发生函数空间Bj,离散平均数就是指由下式在连续基础上定义的平均数的近似(通过求积公式)(这里我们再一次略去下标以简化符号):
xj
| | β | | B = ( ∫ x min x max dx g 2 ( x ) ) ( ∫ 0 T dt 1 τ ( t ) β 2 ( t ) )
在这种情况下,利用bj(x,T)=0,x∈[xmin,xmax]的事实,线性算符Tj实现Bj和D之间的等比例关系。然而由于Tjβj的计算是以数字方式实现的,算符Tj的等比例性质并不是理想的但只是近似而已,当离散化间隔Δx,Δt,Δxj′,Δtj′小时,近似是很好的。公有近似上的等比例的事实将导致有关第8节提出的算法实现的较低性能。
较好的选择是半平均数的各噪声一发生函数空间Bj(再一次略去下标j以简化符号):
| | β | | B = ( ∫ x min x max dx g 2 ( x ) ) ( ∫ 0 T dt 1 τ ( t ) β 2 ( t ) )
在这种情况下,算符Tj严格地在Bj和D之间实现等比例,至少如果传输方程的离散化方案的解对n=N同样为零时是这样。
成本函数的说明
我们通过公式定义成本函数如下:
C ( δm , β 1 , β 2 , . . . β J ) = | | F NL ( m o ) + F ( δm ) + Σ j = I J T j β j - d | | D 2
应该指出,上述函数中FNL(mo)是恒定矢量。
寻求模型和使成本函数最小的噪声-发生函数,通过采取噪声模拟算符的等比例特性的显性或隐性优点的算法来实现这一搜索(见§6)
该算法利用Schur补数的特性,一种与算符Tj的等比例特性有关联的特性。可通过实现使C(δm,β1,β2,…βJ)最小化相当于使 C ′ ( δm ) = C ( δm , β ~ 1 ( δm ) , β ~ 2 ( δm ) , · · · , β ~ J ( δm ) 最小化(其中对给定的 δm , β ~ 1 ( δm ) , β ~ 2 ( δm ) , · · · , β ~ J ( δm ) 使C(δm,β1,β2,…βJ)最小化)来达到。应该注意,与这一二次式有关的Hession是Schur的补数。C′(δm)可通过任何最佳化方法如准牛顿法的BFGS版本来最小化;成本函数C′(δm)是二次式,这里用共轭梯度法特别合适。
显然,如果分别选中了选择(4)和(5)来确定‖‖D和‖‖Bj,则由于算符Tj的等比例特性, ( β ~ 1 ( m ) , β ~ 2 ( m ) , · · · , β ~ J ( m ) ) 的决定非常容易。如果只存在一类噪声(J=1)(这是我们说明的情况),并如果等比例性是理想的(这也是由于作出选择的情况),则 的决定只要求估计该二次式的梯度。如存在几种噪声的叠加(J>1),可考虑各种算法以采取算符Tj的等比例特性的优点,如第1应用例中说明的那样。
图14B示出当选择函数g(x)为很简单形式(仿射函数)时实施该方法得到的(解模型的地震响应)结果。可以看出对传统的反演结果的改进(图14A)。
以上说明了模拟以下分布的物理参数是声阻抗的实施例。显然本方法在它的最一般定义上可用来寻求任何不均匀介质中被相关的噪声影响的物理量的分布。

Claims (10)

1.一种从不均匀介质区域探测得到的数据来估计该区中代表至少一个物理量分布的模型的方法,所述模型不存在数据中可能含有的相关噪声,其特征在于,它包含下列步骤:
a)根据预定的实验协议,获取对给定有关该区的某物理特性的信息的测量值,
b)给出将构成模型的响应的合成数据与各物理量的模型相关联的模拟算符的详细说明,所述测量值和所述合成数据属于数据空间,
c)对标有从1至J的下标j的各相关噪声,选择模拟算符,所述算符将相关噪声与属于预定的噪声发生函数空间(Bj)的噪声发生函数相关联,
d)给出所述数据空间中平均数或半平均数的详细说明,
e)给出上述噪声发生函数空间中半平均数的详细说明,各噪声模拟算符对在该噪声发生函数空间和数据空间之间建立大体上的等比例的关系,
f)定义-成本函数,所述成本函数定量在一方为测量值和另一方为模型响应和与上述噪声发生函数有关的相关噪声的叠加之间的差异,以及
g)通过采取噪声模拟算符的等比例特性的优点的算法,使所述成本函数最小来调整所述模型和所述噪声发生函数。
2.如权利要求1所述的方法,其特征在于,响应于局部地震激励,寻找作为介质中声阻抗的深度的函数的分布,影响所述数据的相关噪声是通过表征其传播参数每一个识别的管道波,所述测得的数据是通过适合于检测介质中质点位移的传感器得到的VSP数据,确定所述传感器的位置、记录时间以及时间采样点,并且所述的模拟算符使所述合成数据与声阻抗分布相关联作为在传播时间中估计深度的函数,并与测得的垂直应力相关联作为在第一传感器深度上时间的函数。
3.如权利要求2所述的方法,其特征在于,定量所述差异的所述成本函数是所述数据空间中这个差异的半平均数的平方。
4.如权利要求2或3所述的方法,其特征在于,通过块驰豫法来消除对应于各相关的噪声发生函数的未知数,得到模型和噪声发生函数的调整,所述驰豫法在准牛顿算法的迭代中实施用于计算该模型。
5.如权利要求2至4中任一项所述的方法,其特征在于,通过选择在传感器位置上和以前确定的时间采样点上的质点位移得到的值,并通过施加算符适当补偿球面发散和衰减效应,用所考虑模型的一维波方程的数字解,实现用模拟算符对模型映像的数字计算。
6.如权利要求2至5中任一项所述的方法,其特征在于,所述数字噪声模拟算符是一使噪声迁移方程离散化的有限差分中心数字方案,并且包括作为沿观察区边缘初始条件的噪声发生函数属于在给定时间间隔内由支持时间函数组成的空间(Bj)。
7.如权利要求2至6中任一项所述的方法,其特征在于,为数据空间所选的半平均数是:
| | u | | D = ( ΔxΔt Σ i = 0 I - 1 Σ n = 0 N - 1 1 τ n + 1 2 ( u i n + 1 + u i n ) 2 ) 1 2
以及为噪声发生函数空间所选的半均数是:
| | β | | B = ( ΔxΔtI Σ n = 0 N - 1 1 τ n + 1 2 ( β n + 1 + β n ) 2 ) 1 2
8.如权利要求1所述的方法,其特征在于,寻求与以前所选的基准模型中有关的阻抗扰动的分布和所述介质区中速度的分布,影响数据的相关噪声是由于多重的反射引起的,其随抵消而变化的动能与幅值业已估计,所测的数据用地震检测表面传感器检拾,传感器的位置、地震激励模式、记录时间及时间采样点均被确定,模拟算符通过基准模型附近使波动方程线性化加以确定。
9.如权利要求8所述的方法,其特征在于,所述定量该差异的成本函数是数据空间中该差异的半平均数的平方。
10.如权利要求8或9所述的方法,其特征在于,通过块驰豫法来消除对应于各相关噪声发生函数的未知数,得到模型和噪声发生函数的调整,该驰豫法在共轭梯度算法的迭代中实施用于计算该模型。
CNB031530249A 2002-08-05 2003-08-05 形成地下区域中代表一物理量分布的模型的方法 Expired - Fee Related CN1333267C (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
FR02/09.919 2002-08-05
FR02/09919 2002-08-05
FR0209919A FR2843202B1 (fr) 2002-08-05 2002-08-05 Methode pour former un modele representatif de la distribution d'une grandeur physique dans une zone souterraine, affranchi de l'effet de bruits correles entachant des donnees d'exploration

Publications (2)

Publication Number Publication Date
CN1503006A true CN1503006A (zh) 2004-06-09
CN1333267C CN1333267C (zh) 2007-08-22

Family

ID=27839441

Family Applications (1)

Application Number Title Priority Date Filing Date
CNB031530249A Expired - Fee Related CN1333267C (zh) 2002-08-05 2003-08-05 形成地下区域中代表一物理量分布的模型的方法

Country Status (5)

Country Link
US (1) US7092823B2 (zh)
CN (1) CN1333267C (zh)
FR (1) FR2843202B1 (zh)
GB (1) GB2391665B (zh)
NO (1) NO328835B1 (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102112894A (zh) * 2008-08-11 2011-06-29 埃克森美孚上游研究公司 用地震表面波的波形评估土壤性质
CN102037379B (zh) * 2008-03-31 2014-07-23 斯塔特伊石油公司 匹配第一和第二地震反射数据集的反射时移的方法

Families Citing this family (53)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070104028A1 (en) * 2005-11-04 2007-05-10 Dirk-Jan Van Manen Construction and removal of scattered ground roll using interferometric methods
CA2664352C (en) * 2006-09-28 2011-09-27 Exxonmobil Upstream Research Company Iterative inversion of data from simultaneous geophysical sources
US8504342B2 (en) * 2007-02-02 2013-08-06 Exxonmobil Upstream Research Company Modeling and designing of well drilling system that accounts for vibrations
US9279899B2 (en) * 2007-07-18 2016-03-08 Westerngeco L.L.C. System and technique to estimate physical propagation parameters associated with a seismic survey
CA2703588C (en) 2007-12-12 2015-12-01 Exxonmobil Upstream Research Company Method and apparatus for evaluating submarine formations
US8812282B2 (en) 2008-03-21 2014-08-19 Exxonmobil Upstream Research Company Efficient method for inversion of geophysical data
FR2933499B1 (fr) * 2008-07-03 2010-08-20 Inst Francais Du Petrole Methode d'inversion conjointe de donnees sismiques representees sur des echelles de temps differentes
US8451687B2 (en) * 2009-02-06 2013-05-28 Westerngeco L.L.C. Imaging with vector measurements
FR2945869B1 (fr) 2009-05-20 2011-05-20 Inst Francais Du Petrole Methode pour imager une zone cible du sous-sol a partir de donnees de type walkaway
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
US8223587B2 (en) * 2010-03-29 2012-07-17 Exxonmobil Upstream Research Company Full wavefield inversion using time varying filters
US8694299B2 (en) 2010-05-07 2014-04-08 Exxonmobil Upstream Research Company Artifact reduction in iterative inversion of geophysical data
US8756042B2 (en) 2010-05-19 2014-06-17 Exxonmobile Upstream Research Company Method and system for checkpointing during simulations
US8767508B2 (en) 2010-08-18 2014-07-01 Exxonmobil Upstream Research Company Using seismic P and S arrivals to determine shallow velocity structure
US9063882B1 (en) * 2010-09-09 2015-06-23 Sas Ip, Inc. Matrix preconditioners for simulations of physical fields
EP2622457A4 (en) 2010-09-27 2018-02-21 Exxonmobil Upstream Research Company Simultaneous source encoding and source separation as a practical solution for full wavefield inversion
US8437998B2 (en) 2010-09-27 2013-05-07 Exxonmobil Upstream Research Company Hybrid method for full waveform inversion using simultaneous and sequential source method
CN102043168B (zh) * 2010-10-15 2012-11-07 中国石油化工股份有限公司 一种对数字信号进行仿真加噪的处理方法
CN102096093B (zh) * 2010-11-29 2013-04-10 杨本才 一种利用微震点作为震源计算矿区地震波传播速度的方法
KR101797451B1 (ko) 2010-12-01 2017-11-14 엑손모빌 업스트림 리서치 캄파니 상호상관 목적 함수를 통한 해양 스트리머 데이터에 대한 동시 소스 반전
KR101931488B1 (ko) 2011-03-30 2018-12-24 엑손모빌 업스트림 리서치 캄파니 스펙트럼 성형을 이용하는 전 파동장 반전의 수렴 레이트
WO2012134609A1 (en) 2011-03-31 2012-10-04 Exxonmobil Upstream Research Company Method of wavelet estimation and multiple prediction in full wavefield inversion
US9140812B2 (en) 2011-09-02 2015-09-22 Exxonmobil Upstream Research Company Using projection onto convex sets to constrain full-wavefield inversion
US9176930B2 (en) 2011-11-29 2015-11-03 Exxonmobil Upstream Research Company Methods for approximating hessian times vector operation in full wavefield inversion
CA2861863A1 (en) 2012-03-08 2013-09-12 Exxonmobil Upstream Research Company Orthogonal source and receiver encoding
CN102681015B (zh) * 2012-05-30 2014-06-18 中国地质大学(北京) 矿区地层结构划分方法
CN103675905B (zh) * 2012-09-14 2016-10-05 中国科学院地质与地球物理研究所 一种基于优化系数的地震波场模拟方法及装置
MY178811A (en) 2012-11-28 2020-10-20 Exxonmobil Upstream Res Co Reflection seismic data q tomography
US10088588B2 (en) * 2013-04-03 2018-10-02 Cgg Services Sas Device and method for stable least-squares reverse time migration
CN105308479B (zh) 2013-05-24 2017-09-26 埃克森美孚上游研究公司 通过与偏移距相关的弹性fwi的多参数反演
US10459117B2 (en) 2013-06-03 2019-10-29 Exxonmobil Upstream Research Company Extended subspace method for cross-talk mitigation in multi-parameter inversion
US9702998B2 (en) 2013-07-08 2017-07-11 Exxonmobil Upstream Research Company Full-wavefield inversion of primaries and multiples in marine environment
WO2015026451A2 (en) 2013-08-23 2015-02-26 Exxonmobil Upstream Research Company Simultaneous sourcing during both seismic acquisition and seismic inversion
US10036818B2 (en) 2013-09-06 2018-07-31 Exxonmobil Upstream Research Company Accelerating full wavefield inversion with nonstationary point-spread functions
US9910189B2 (en) 2014-04-09 2018-03-06 Exxonmobil Upstream Research Company Method for fast line search in frequency domain FWI
MX2016013366A (es) 2014-05-09 2017-01-26 Exxonmobil Upstream Res Co Metodos de busqueda de linea eficientes para la inversion de campo de ondas completo de multi-parametros.
US10185046B2 (en) 2014-06-09 2019-01-22 Exxonmobil Upstream Research Company Method for temporal dispersion correction for seismic simulation, RTM and FWI
EP3158367A1 (en) 2014-06-17 2017-04-26 Exxonmobil Upstream Research Company Fast viscoacoustic and viscoelastic full-wavefield inversion
US10838092B2 (en) 2014-07-24 2020-11-17 Exxonmobil Upstream Research Company Estimating multiple subsurface parameters by cascaded inversion of wavefield components
US10422899B2 (en) 2014-07-30 2019-09-24 Exxonmobil Upstream Research Company Harmonic encoding for FWI
US10386511B2 (en) 2014-10-03 2019-08-20 Exxonmobil Upstream Research Company Seismic survey design using full wavefield inversion
WO2016064462A1 (en) 2014-10-20 2016-04-28 Exxonmobil Upstream Research Company Velocity tomography using property scans
WO2016099747A1 (en) 2014-12-18 2016-06-23 Exxonmobil Upstream Research Company Scalable scheduling of parallel iterative seismic jobs
US10520618B2 (en) 2015-02-04 2019-12-31 ExxohnMobil Upstream Research Company Poynting vector minimal reflection boundary conditions
SG11201704620WA (en) 2015-02-13 2017-09-28 Exxonmobil Upstream Res Co Efficient and stable absorbing boundary condition in finite-difference calculations
CA2972033C (en) 2015-02-17 2019-07-23 Exxonmobil Upstream Research Company Multistage full wavefield inversion process that generates a multiple free data set
WO2016168280A1 (en) * 2015-04-14 2016-10-20 Schlumberger Technology Corporation Generating an accurate model of noise and subtracting it from seismic data
EP3304133A1 (en) 2015-06-04 2018-04-11 Exxonmobil Upstream Research Company Method for generating multiple free seismic images
US10838093B2 (en) 2015-07-02 2020-11-17 Exxonmobil Upstream Research Company Krylov-space-based quasi-newton preconditioner for full-wavefield inversion
US10310113B2 (en) 2015-10-02 2019-06-04 Exxonmobil Upstream Research Company Q-compensated full wavefield inversion
CN108139498B (zh) 2015-10-15 2019-12-03 埃克森美孚上游研究公司 具有振幅保持的fwi模型域角度叠加
US10768324B2 (en) 2016-05-19 2020-09-08 Exxonmobil Upstream Research Company Method to predict pore pressure and seal integrity using full wavefield inversion
US11030780B2 (en) * 2018-03-26 2021-06-08 The Board Of Trustees Of The Leland Stanford Junior University Ultrasound speckle reduction and image reconstruction using deep learning techniques

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4809239A (en) * 1987-07-14 1989-02-28 Schlumberger Technology Corporation Method for evaluating parameters related to the elastic properties of subsurface earth formations
US4839522A (en) 1987-07-29 1989-06-13 American Screen Printing Company Reflective method and apparatus for curing ink
IT1263156B (it) * 1993-02-05 1996-08-01 Agip Spa Procedimento e dispositivo di rilevamento di segnali sismici per ottenere profili sismici verticali durante le operazioni di perforazione
US5715213A (en) * 1995-11-13 1998-02-03 Mobil Oil Corporation High fidelity vibratory source seismic method using a plurality of vibrator sources
US6131071A (en) * 1996-12-06 2000-10-10 Bp Amoco Corporation Spectral decomposition for seismic interpretation
US5835883A (en) * 1997-01-31 1998-11-10 Phillips Petroleum Company Method for determining distribution of reservoir permeability, porosity and pseudo relative permeability
FR2763702B1 (fr) * 1997-05-23 1999-07-09 Elf Exploration Prod Methode d'elaboration de cartes de risques de positionnement d'un puits dans un milieu
US6263284B1 (en) * 1999-04-22 2001-07-17 Bp Corporation North America Inc. Selection of seismic modes through amplitude characteristics
FR2800473B1 (fr) * 1999-10-29 2001-11-30 Inst Francais Du Petrole Methode pour modeliser en 2d ou 3d un milieu heterogene tel que le sous-sol decrit par plusieurs parametres physiques

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102037379B (zh) * 2008-03-31 2014-07-23 斯塔特伊石油公司 匹配第一和第二地震反射数据集的反射时移的方法
CN102112894A (zh) * 2008-08-11 2011-06-29 埃克森美孚上游研究公司 用地震表面波的波形评估土壤性质
CN102112894B (zh) * 2008-08-11 2015-03-25 埃克森美孚上游研究公司 用地震表面波的波形评估土壤性质

Also Published As

Publication number Publication date
US7092823B2 (en) 2006-08-15
GB0318295D0 (en) 2003-09-10
NO328835B1 (no) 2010-05-25
GB2391665B (en) 2005-09-28
GB2391665A (en) 2004-02-11
FR2843202A1 (fr) 2004-02-06
NO20033459D0 (no) 2003-08-04
US20050075791A1 (en) 2005-04-07
CN1333267C (zh) 2007-08-22
FR2843202B1 (fr) 2004-09-10

Similar Documents

Publication Publication Date Title
CN1503006A (zh) 形成地下区域中代表一物理量分布的模型的方法
CN1254696C (zh) 处理地震数据的方法和设备
CN1188710C (zh) 估算地震介质特性的系统和方法
CN1148585C (zh) 反射剪切波地震方法
Liu et al. Migration velocity analysis: Theory and an iterative algorithm
CN1203324C (zh) 自适应地震噪声和干扰衰减方法
Lebedev et al. Mapping the Moho with seismic surface waves: A review, resolution analysis, and recommended inversion strategies
Héloïse et al. Site effect assessment using KiK-net data: Part 1. A simple correction procedure for surface/downhole spectral ratios
CN1188711C (zh) 用于地震波场分离的系统和方法
CN101052896A (zh) 改进的地震检波器校准技术
CN1278132C (zh) 海洋延时地震勘测的方法
Toldi Velocity analysis without picking
CN1496485A (zh) 齐发震动地震数据的处理
CA2792176A1 (en) A process for characterising the evolution of a reservoir
US20120053839A1 (en) Method of detecting or monitoring a subsurface hydrocarbon reservoir-sized structure
Al Atik et al. A methodology for the development of 1D reference VS profiles compatible with ground‐motion prediction equations: Application to NGA‐West2 GMPEs
CN1275052C (zh) 处理地震数据
CN100344993C (zh) 叠加前地震数据的向下外推方法
CN1299127C (zh) 地震观测系统优化设计的层状介质双聚焦方法及其应用
Kamath et al. Full-waveform inversion of multicomponent data for horizontally layered VTI media
Ikeda et al. Two-station continuous wavelet transform cross-coherence analysis for surface-wave tomography using active-source seismic data
Lin et al. BEHAVIOR OF APPARENT DISPERSION CURVE AND ITS IMPLICATION TO MASW TESTING.
Price et al. Surface microseismic imaging in the presence of high-velocity lithologic layers
Liu et al. Post-critical SsPmp and its applications to Virtual Deep Seismic Sounding (VDSS)–2: 1-D imaging of the crust/mantle and joint constraints with receiver functions
Gao et al. Acquisition and processing pitfall with clipped traces in surface-wave analysis

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
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