专利名称：形成地下区域中代表一物理量分布的模型的方法
技术领域：
本发明涉及从不均匀介质区探测得到的数据形成该区中物理量分布的有代表性的模型的方法，所述模型(至少部分)不存在数据中可能含有的相关噪声。
该方法运用如地下区域中的声阻抗的定量法。
背景技术：
该方法在于寻找在几乎所有科技领域业已开发用来调整实验测量的模型。这种方法称为下面的种种名称对反演问题解法的参数估计最小2乘方法。对这一方法在球科学范围内的较好表示，可参考——Tarantola，A.《Inverse Problem TheoryMethed for Data Fitting andModel Parameter Estimation》，Elsevier，Amsterdam，1987。
应该指出，术语“最小2乘方”是指数据空间中用来定量在模型的响应(它是由以前选定的模拟算符的模型的映像)与数据之间差异的平均数的平方，一种必须使其最小以解决问题的成本函数。用平均数的平方来定义成本函数只是实用上的方便，从根本上说并不是本质的。此外，许多作者由于种种理由采用成本函数的另一种定义，但这种定义仍基于数据空间中平均数或半平均数的应用。最后，我们有较大自由度选用数据空间中的平均数(或半平均数)(我们决定被迫使用欧几里德平均数)。在含有噪声的数据的场合，解决大体上决定于这一阶段作出的选择。对有关这个问题的较多进层可参考下列出版物-—Tarantola，A.《Inverse Problem TheoryMethed 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=1JBjC(m,β1,β2,...βJ)=||F(m)+Σj=1JTjβj-d||D2(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标出的样值，-由这些传感器实施的测量值的物理特性在我们的实验中，传感器测量波传播导致的垂直位移，-记录时间在我们的实验中，传感器测量在时间区
中的振动状态，这些数据每Δt＝1ms被采样1次，从而获得从0至N标出的样值。
我们称在深度Xmin+iΔx和时间nΔt上记录的数据样值为din。当然有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)∂2y∂t2-∂∂x(σ(x)∂y∂x)=0]]>在]xmin，+∞[×
中，Neumann边界条件-σ(xmin)∂y∂x(xmin,t)=h(t)on
]]>以及初始条件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有关的(非线性)算符。
与各相关的噪声有关的模拟算符的说明本节中我们说明用于模拟各相关噪声的步骤。各相关噪声由假设已知或至少近似知道的对其相关方向特性来表征通过部分(cjx(x，t)，cjt(x，t))(在域[xmin，xmax]×
的任一点上需被确定)，我们可以等效的方法通过代表与cj(x，t)有关的场线的相关线束来规定相关方向。按照申请人提交的专利EP-354,112(US-4,972,383)中所述的技术可以得到这种相关线，从采样噪声的一些特性通过传统的内插和/或外插步骤，信息被扩展至整个域[xmin，xmax]×
。相关方向的说明就是说明了以关系cj(x,t)=cjx(x,t)/cjt(x,t)]]>与相关矢量关联的噪声cj(x，t)的传播速度分布。假设在波传播过程中波幅度值不变(不同情况在第2实施例中提出)。如上所述相关噪声的模拟基于下述的观察波沿相关线传播而无幅值改变是下述传输方程的解▿→b(x,t).c→(x,t)=0]]>对此，我们引入-噪声-发生函数的空间Bj，它将是支持函数βj(t)在[tjmin，tjmax]
内的空间，我们用0在整个区间
中扩展，-噪声模拟算符Tj
Tjβj(t)∈Bj→bj(x，t)∈D传输方程的解▿→bj(x,t).c→j(x,t)=0]]>在[xmin，xmax]×
中以及满足初始条件 应该指出，相关线的束的几何形状(或说传播速度分布)不能是任意的几何形状相关线不能相交，否则传输方程的解将不是正确的组。同样理由，相关线不能与被规定的初始条件的边缘交切两次。后一种考虑可导致(为模拟对应于上升波的相关噪声)规定初条件不再在x＝xmin而在x＝xmax，或甚至中间值上，即使它意味着将Cauchy问题分解成两个与位于该中间值边缘任一侧上的域有关的问题。
我们进一步假设相关线的几何形状为(考虑更一般情况，而后面的结果取不同形式)cj(x，t)＝ξj(x)×τ(t)ou encore:cjx(x,t)=ξj(x)etcjt(x,t)=1τ(t)]]>最后我们假设相关线的束在x＝xmin区间[tjmin，tjmax]相交，而不在t＝T区间[xmin，xmax]相交。
事实上，用数字方案来解传输方程。如可用传统的有限差分对中方案。因此，我们选择离散化区间Δxj和Δtj(它们不必等于Δx和Δt，但我们假设为它们的约数，即使可考虑更一般的情况)，并引入一格栅，其节点是坐标(xmin+i′Δxj′，n′Δtj′)，i′和n′∈N的点。传统的有限差分对中方案(以后说明)使能从初始条件来逐步计算各格栅节点上函数bj(x，t)所取的不同的值(为简化符号我们将略去与相关噪声有关的下标j)。
ξi′+122Δt′[bi′+1n′+1+bi′n′+1-bi′+1n′-bi′n′]+]]>12Δx′τn′+12[bi′+1n′+1+bi′+1n′-bi′n′+1-bi′n′]=0]]>上式中，量αin′表示在坐标点(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=(∫xminxmaxdx∫0Tdt1τ(t)u2(x,t))12]]>我们也可使用半平均数，从下节可见这是较好的选择||u||D=(ΔxΔtΣi=0I-1Σn=0N-11τn+12(uin+1+uin)2)12]]>式中u为数据空间的任意矢量(i和n为分别表示空间和时间中样值数目的索引)。
此外，其他选择也是可能的。
在噪声-发生函数的各空间Bj中各噪声模拟算符Tj为之在噪声-发生函数空间和数据空间(有§5中定义的半平均数)之间建立等比例关系(或近似等比例关系)的平均数的说明本节中我们说明在与各相关噪声有关联的噪声-发生函数空间中用来定义平均数的步骤。这个步骤对各噪声都相同，我们一般地加以说明并在公式中略去与所述噪声有关的下标j。
首先选有离散平均数的各噪声-发生函数空间B，离散平均数就是指由下式在连续基础上定义的平均数的近似(通过求积公式)||β||B=(X∫0Tdt1τ(t)β2(t))12]]>在这种情况下，利用在§4提出的假设，我们有bi(x，t)＝0，x∈[xmin，xmax]，线性的算符Tj实现Bj和D之间的等比例关系。然而，由于Tjβj的计算是以数字方式实现的，算符T的等比例性质并不是理想的但只是近似而已，当离散化间隔Δx，Δt，Δxj′，Δtj′小时，近似更好。仅有近似上的等比例的事实将导致有关在后面第8节提出的算法实现的较低性能。
较好的选择是有半平均数的各噪声-发生函数空间Bj||β||B=(ΔxΔtIΣn=0N-11τn+12(βn+1+βn)2)12]]>在这种情况下，算符Tj严格地在Bj和D之间实现等比例，至少在uiN＝0Ai＝0，I-1时如此。从下式离散的能量等式得到结果IΣn=0N-11τn+1/2[uon+1+uon)2]=Σi=0I-1Σn=0N-11τn+1/2[(uin+1+uin)2]]]>+Σi=0I-1Σk=0iΔxξk+1/2Δt[(uk+1N+ukN)2]-IΔx(uoN)2ξ-1/2Δt]]>在Δx′＝Δx、Δt′＝Δt并定义ξ-1/2＝ξ1/2的情况下等式有效。
成本函数的说明我们通过公式定义成本函数如下C(m,β1,β2,...βJ)=||F(m)+Σj=1JTjβj-d||D2]]>寻求模型和使成本函数最小的噪声-发生函数，通过采取噪声模拟算符的等比例特性的显性中隐性优点的算法来实现这一搜索(见ξ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)以及等比例性是理想的，则 的决定只要求估计该二次式的梯度。如存在几种噪声的叠加(J＞1)，可考虑各种算法以打取算符Tj的等比例特性的优点。可通过例如块驰豫法(Gauss-Seidel混淆)来决定β~1(m),β~2(m),···,β~J(m),]]>该法将块与一组未知数 相关联，迭弛法只要求(仍在理想的等比例情情况下)估计与所述块有关的二次式的梯度。对于使得难以接近各种噪声类型的相关方向的问题，可考虑用预处理的共轭梯度法，该预处理矩阵具有对称块(弛豫混淆)Gauss-Seidel型式。因此，预处理再一次是非常容易实现的，与块相关的问题的解只要求估计与所述块有的二次式的梯度。
图6至11示出了用该方法的三类实验得到的结果。
在既有相关噪声又有随机噪声的叠加的含噪声数据(图4)的反演的情况下，图6给出所得到的阻抗分布，图7示出反演剩余。
在既有相关噪声又有随机噪声的叠加的含噪声数据(图5)的反演情况下，图8给出所得到的阻抗分布，图9示出反演的剩余。
在既有相关的噪声又有随机噪声的叠加的含噪声数据(图4)的反演，其中流体波的传播特性并不精确知晓的情况下，自然要寻求包括两相关噪声的叠加的相关噪声模型，所述两相关噪声具有不准确的传播速度，但通过它构建真实的流体波的传播速度。图10验出所得的阻抗特性，图11示出反演剩余。
通过线性反演由多次反射所混杂的多低消数据来确定模型的方法适用的第2实施例本应用不同于第1例的主要在于考虑到噪声幅值沿相关方向变化时说明本方法提供的可能性。另一个差别在于相关噪声的传播性质，多次反射的传播没有色散。
线性化反演是地球物理工作者所用的一种先进的方法，来获得地下不均匀性的定量模型，这些信息对油层特性是明显地宝贵的。这里给出的数据是表面地震检测数据。例如如下的刊物说明了这一方法-Bourgeois，A.et al.(1989)《The Linearized Seismic Inverse Probleman 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的轨迹，
-传感器测量从波传播产生的压力�。�-记录时间在我们的实验中，传感器测量在时间区
的振动状态，这些数据在每△t＝1ms时被采样。得到编号从0至N的样值。
-我们称对抵消xmin+i△x和在时间n△t时记录的数据样值为din。当然有T＝N△t，xmax＝xmin+I△x。
按照实验协议的数据的获得从合成数据实现我们的实验这些数据通过下面2维波动方程的数字解(有限差分方法)得到1σc∂2P∂t2-▿σc▿→P=δ(x,z-zs)w(t)]]>在 中边界条件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的分量。
与各相关噪声有关的模拟算符的说明本节中我们说明用于模拟各相关噪声的步骤。各相关噪声由假设已知或至少近似知道的对其相关方向的特性来表征。通过分量(cjx(x，t)，cjt(x，t))(在域[xmin，xmax]×
的任一点上需被确定)的相关矢量 的场规定这些相关方向。我们可以等效的方法通过代表与 有关的场线的相关线束来规定相关方向。按照申请人提交的专利EP-354,112(US-4,972,383)中所述的技术可以得到这种相关线，从采集噪声的一些特性通过传统的内插和/或外插步骤，如申请人提交的上述专利Ep-354,112所述，信息被扩展至整个域[xmin，xmax]×
。对于要消除上述的多次反射的第2应用例，我们采集这一多次反射(如幅值峰值)，它能使我们确定作为抵消点的函数的到达时间改变，通过将所采集线的简单垂直移动确定相关线的连续区(在这些条件下，矢量 不决定于t)。相关方向的说明就是说明了以关系cj(x,t)=cjx(x,t)/cjt(x,t)]]>与相关矢量关联的噪声cj(x，t)的传播速度分布。
与上述第1应用例不一样，这里噪声幅值沿相关线改变。假设我们能根据数据的测量或据理论上的考虑来估计这些幅值改变，则这一测量确定了函数g(x)。
因而通过两个算符的合成；噪声得以模拟传输算符(如第1应用例中)；幅值调制，即以函数g(x)乘以传输方程的解。
如果我们再一次利用第1应用例的方案(显然如使用相同的假设和符号)，则我们引入-噪声发生函数的空间Bj，它将是支持函数βj(t)在[ttjmin，tjmax]
中的空间，我们用0在整个区间
中扩展；噪声模拟算符TjTj∶βj(t)∈Bj→bj(x，t)∈D传输方程的解▿→bj(x,t).c→j(x,t)=0]]>在[xmin，xmax]×
并满足初始条件
事实上，用数字方案解传输方程。如利用传统的有限差分对中方案是可能的。所以我们选择离散化间隔Δxj′和Δtj′(不必等于Δx和Δt，但我们假定为这些量的约数，尽管可考虑更一般的情况)，并引入一格栅，其节点是坐标(xmin+i′Δxj，n′Δt′j)，i′和n′∈N的点。下面说明的传统有限差分对中方案使可根据初始条件逐步计算函数bj(x，t)在各格栅节点所取的不同的值(为简化符号，略去与所述相关噪声有关的下标j)。
ξi′+122Δt′[bi′+1n′+1+bi′n′+1-bi′+1n′-bi′n′]+]]>12Δx′τn′+12[bi′+1n′+1+bi′+1n′-bi′n′+1-bi′n′]=0]]>上式中，量αi′n′代表在坐标(xmin+i′Δxj′+n′Δtj′)处量α(α是一般项)的估值。
这里，必须选择离散化间隔Δxj′和Δtj′相对于波长为足够小，因为我们要求给出传输方程的函数解的精确近似的数字方案(我们要求不色散)。
噪声模拟步骤的最后一步是选择值bi′n′它属于两个格栅共有的节点并以g(iΔx)与它们相乘从而我们得到矢量(数据空间的)Tj(βj)的各分量。
数据空间中平均数或半平均数的说明首先取任意离散的平均数作为数据空间的平均数，这意味着求由下式在连续基础上确定的平均数的近似值(通过求积公式)||u||D=(∫xminxmaxdx∫0Tdt1τ(t)u2(x,t))12]]>我们也可使用半平均数，从下节可见这是较好的选择||u||D=(ΔxΔtΣi=0I-1Σn=0N-11τn+12(uin+1+uin)2)12]]>式中u为数据空间的任意矢量(i和n为分别表示空间和时间中样值数目的索引)。
此外，其他选择也是可能的。
在噪声-发生函数的各空间Bj中各噪声模拟算在Tj为之在噪声-发生函数空间和数据空间(有§5中定义的半平均数)之间建立等比例关系(或近似等比例关系)的平均数的说明。
本节中我们说明在与各相关噪声有关联的噪声-发生函数空间中用来定义平均数的步骤。
首先选有离散平均数的各噪声-发生函数空间Bj，离散平均数就是指由下式在连续基础上定义的平均数的近似(通过求积公式)(这里我们再一次略去下标以简化符号)xj′||β||B=(∫xminxmaxdxg2(x))(∫0Tdt1τ(t)β2(t))]]>在这种情况下，利用bj(x，T)＝0，x∈[xmin，xmax]的事实，线性算符Tj实现Bj和D之间的等比例关系。然而由于Tjβj的计算是以数字方式实现的，算符Tj的等比例性质并不是理想的但只是近似而已，当离散化间隔Δx，Δt，Δxj′，Δtj′小时，近似是很好的。公有近似上的等比例的事实将导致有关第8节提出的算法实现的较低性能。
较好的选择是半平均数的各噪声一发生函数空间Bj(再一次略去下标j以简化符号)||β||B=(∫xminxmaxdx g2(x))(∫0Tdt1τ(t)β2(t))]]>在这种情况下，算符Tj严格地在Bj和D之间实现等比例，至少如果传输方程的离散化方案的解对n＝N同样为零时是这样。
成本函数的说明我们通过公式定义成本函数如下C(δm,β1,β2,...βJ)=||FNL(mo)+F(δm)+Σj=IJTjβj-d||D2]]>
应该指出，上述函数中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)。
以上说明了模拟以下分布的物理参数是声阻抗的实施例。显然本方法在它的最一般定义上可用来寻求任何不均匀介质中被相关的噪声影响的物理量的分布。
权利要求
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=0I-1Σn=0N-11τn+12(uin+1+uin)2)12]]>以及为噪声发生函数空间所选的半均数是||β||B=(ΔxΔtIΣn=0N-11τn+12(βn+1+βn)2)12]]>
8.如权利要求1所述的方法，其特征在于，寻求与以前所选的基准模型中有关的阻抗扰动的分布和所述介质区中速度的分布，影响数据的相关噪声是由于多重的反射引起的，其随抵消而变化的动能与幅值业已估计，所测的数据用地震检测表面传感器检拾，传感器的位置、地震激励模式、记录时间及时间采样点均被确定，模拟算符通过基准模型附近使波动方程线性化加以确定。
9.如权利要求8所述的方法，其特征在于，所述定量该差异的成本函数是数据空间中该差异的半平均数的平方。
10.如权利要求8或9所述的方法，其特征在于，通过块驰豫法来消除对应于各相关噪声发生函数的未知数，得到模型和噪声发生函数的调整，该驰豫法在共轭梯度算法的迭代中实施用于计算该模型。
全文摘要
一种从不均匀区域(如地下区域)探测得到的数据来形成该区中至少一个物理量分布的有代表性的模型的方法，所述模型不存在数据中可能含有的相关噪声。该方法包括下列步骤获取给定有关该区特性的信息的数据，说明模拟算符，它将模型的响应与模型相关联(合成数据)，通过对噪声发生函数(n、g、f)施加特定模拟算符模拟各相关噪声，说明数据空间中的半平均数，说明各n、g、f空间中的平均数，各噪声模拟算符对平均数构成等比例关系，以及借助采取这些等比例关系特性优点的算法寻求模型和使成本函数最小化的n、g、f，成本函数通过半平均数来定量在测量数据与模型的响应和与n、g、f有关联的噪声的叠加之间的差分，应用例如寻求地下区域中声阻抗、传播速度、磁导率等的分布。
文档编号G01V1/36GK1503006SQ0315302
公开日2004年6月9日 申请日期2003年8月5日 优先权日2002年8月5日
发明者P·莱利, F·雷纳德, L·佩莱, F·德尔普拉－詹诺德, P 莱利, 傻, 绽 詹诺德 申请人:法国石油研究所
专业数控铣床产品供应商
20年生产经验，实力工厂，精工品质，质量有保障
