CN112904430A - 非线性直接叠前地震泊松阻抗反演的计算机实现的方法 - Google Patents

Abstract

本文公开了一种非线性直接叠前地震泊松阻抗反演的计算机实现的方法，以及实现该方法的系统，用于使用非线性直接叠前地震泊松阻抗反演来计算弹性属性的最终模型。在勘测区域的各个入射点处的各个入射角下得到用户输入和地球模型数据。然后计算各种模型以用于岩性识别和流体鉴定，并用于初步地震勘探和储层表征。因此，由于它们对多孔介质的储层有局限性，埋藏深度、压实度和上覆压力的变化会导致这些模型需要进一步优化。通过系统计算机来执行使用非线性直接叠前地震模型的进一步优化，产生了弹性属性的最终模型。这一模型可应用于岩性识别和流体鉴定，识别油气勘探中的潜在目标，并估算非常规页岩气应用中的甜点。

Description

非线性直接叠前地震泊松阻抗反演的计算机实现的方法

技术领域

本公开总体上涉及地震勘探和储层表征领域中的用于岩性识别和流体鉴定的计算机实现的方法和系统。

背景技术

1.总论

本领域的普通技术人员已经知道，岩性识别和流体鉴定在地震勘探和储层表征中起着重要的作用。岩性识别和流体鉴定的定义是使用叠前地震数据的加权叠加方法来估算的，最早由Smith,G.和Gidlow,P.在他们的研究“Weighted stacking for rock propertyestimation and detection of gas:Geophysical Prospecting(1987)”中提出。后来，Goodway，B.、T.Chen和J.Downton在SEG第67届国际年会上介绍了“Improved AVO fluiddetection and lithology discrimination using Lamépetrophysical parameters:“λρ”,“μρ”&“λ/μfluid stack”from P and S inversions”，发现拉梅参数对烃饱和度来说比弹性参数更敏感，并提出了用于岩性识别和流体鉴定的LMR方法。之后发现，该方法会取决于埋藏深度、压实度和上覆压力而发生变化，因此对于具有多孔介质的实际储层而言有一定的限制。在2006年，Quakenbush,M.、B.Shang和C.Tuttle引入了泊松阻抗的概念，将其用于岩性识别和流体鉴定。实际上，在他们的被认为较新颖的研究论文“Poisson impedance:The Leading Edge,v.25,no.2,pp.128-138”中，首次根据区分岩性和流体含量的所述方法开发了一种敏感性工具。这开启了开发计算机实现的方法以获取流体属性、执行流体鉴定但不包括烃预测的时代。

泊松阻抗属性随后定义为声阻抗和剪切阻抗之间的加权差异，其通过轴旋转进行优化，以产生更好的岩性和流体含量的分辨率(例如参见Sharma,R.和S.Chopra，2013，Poisson impedance inversion for characterization of sandstone reservoirs，SEG第83届国际年会，扩展摘要，第2549-2553页)。考虑到这一研究，然后通过结合泊松比的识别特性和堆积密度(两者都是用于储层划分的参数)而以地震方式推导泊松阻抗。结果，在过去的十年中，泊松阻抗已成为本领域的技术人员的非常有利的工具，可用于识别新的勘探前景并表征弹性储层的属性，在描述页岩中的含烃和盐砂岩方面显示出许多成功的应用。

另一方面，根据地震振幅以及振幅-偏移关系来进行烃预测仍然是一项艰巨的任务。一种方法是使用地震反射来将其与地下岩石特性紧密相关。然而，地震数据中最强的振幅-偏移关系(AVO)通常是由岩石中的烃饱和引起的。使用叠前地震反演来提取地震数据的关于地下弹性参数方面的信息的进步极大地帮助了以最小的误差来表征岩相和预测储层属性，从而减少了世界上某些盆地中的干井数量和钻井风险(例如参见Russel,B.，2014，Prestack seismic amplitude analysis:an integrated overview:Interpretation,v.2,no.2,SC19-SC36)。这种叠前地震反演模型已经常规地应用于岩性预测和流体探测，以识别油气勘探中的潜在目标。最近，它已被广泛用于估算非常规页岩气应用中的甜点。

叠前地震反演模型的本质在于这样的事实：当流体饱和度改变时，岩石的剪切模量不会改变。当流体饱和度改变时，体积模量则发生显著变化。由于岩石骨架的剪切模量等于具有孔隙流体的岩石的剪切模量，因此可以将剪切阻抗视为主要连接到岩石骨架的地震属性，而声阻抗则由孔隙流体和岩石骨架两者主导。结果，泊松阻抗可以最佳地从剪切模量和声阻抗中消除岩石骨架的影响，从而可以产生更好的流体含量的分辨率。通过叠前地震数据，传统上是根据可以直接从地震数据中求反的P波速度、S波速度和密度而间接地计算出泊松阻抗(见Goodway，同上)。但是，参数估算的间接方法会因间接计算而产生更多的不确定性(例如参见Zong,Z.、X.Yin和G.Wu，2013，Elastic impedance parameterizationand inversion with Young’s modulus and Poisson’s ratio，Geophysics，v.78,no.6,N35-N42)。在2006年，Wang,B.、X.Yin和F.Zhang在他们的研究论文“Laméparametersinversion based on elastic impedance and its application:Applied Geophysics,v.3,pp.120-123”中将从叠前地震数据中估计拉梅参数的直接方法与使用反演P波和S波来估计拉梅参数的间接方法进行了比较，但是得出的结论是，它仍然会造成更大的偏差。

因此，由于所有这些试图减少累积误差以提供准确而合理的结果的尝试均告失败，直接和间接地从叠前地震反演中提取流体因子近年来受到了广泛的关注，特别是考虑到解决上述大型科学和工程问题所需要的具有高性能计算系统的系统计算机的更广泛的可用性和可得性。

2.叠前地震反演的基本依据

Zoeppritz方程(Zoeppritz,K.，1919，Erdbebenwellen VII,VII B,

und DurchgangseismischerWellendurch

Nachrichten von der

Gesellschaft der Wissenschaftenzu

Mathematisch-physikalischeKlasse,pp.66-84.)是叠前地震反演的基础，其描述了在给定入射角下的界面反射率与界面两侧的介质弹性属性之间的非线性关系。但是，Zoeppritz方程的固有的复杂性和高度的非线性使它们无法求解大规模的多层和多维地质模型。在实践中，叠前地震反演方法基于Zoeppritz方程的一阶线性近似，并基于如下假设：两个均匀的半空间在弹性界面处焊接，并且跨边界只发生弹性参数的相对小的变化。由于它们在数学形式和计算效率方面的简单性，在过去的几十年中已经推导了一系列的Zoeppritz方程的线性近似，并用于各种计算机实现的方法和系统(参阅Bortfeld,R.，1961，Approximations to thereflection and transmission coefficients of plane longitudinal and transversewaves:Geophysical Prospecting,v.9,no.4,pp.485-502；Aki Richards方程：Aki,K.和P.Richards，1980，Quantitative seismology:Theory and methods:W.H.Freeman andCo.；Shuey,R.，1985，A simplification of the Zoeppritz equations:Geophysics,v.50,no.4,pp.609-614；以及Fatti,J.、G.Smith、P.Vail、P.Strauss和P.Levitt，1994，Detection of gas in sandstone reservoirs using AVO analysis:A 3-D seismiccase history using the geostack technique,Geophysics,v.59,no.9,pp.1362-1376)。

Aki和Richards提出了最著名的针对Zoeppritz方程的工业标准近似，它们将P波反射率表示为一系列的针对40度以下的入射角和跨边界的弹性属性变化很小的速度和密度对比。Aki和Richards方程是针对Zoeppritz方程的重要线性近似，这导致了其他近似方程的发展，这些近似方程也被业界广泛应用，例如由Shuey和Fatti提出的(同上)。但是，这些方程式是基于对Zoeppritz方程式的一阶近似推导的，并且仅对跨界面的相对较小的弹性属性对比以及界面处的相对较小的入射角有效。在跨过勘测区域的反射边界的弹性和非弹性参数存在较大的相对变化时，Zoeppritz方程的非线性近似将适用于弹性参数的精确确定。通过推导高阶近似或将校正项添加到线性近似中，许多研究人员提高了大角度的非线性近似的准确性(参见Wang,Y.，1999，Approximations to the Zoeppritz equationsand their use in AVO analysis，Geophysics，v.64，no.6，pp.1920-1927；Stovas，A.和B.Ursin，2001，Second order approximations of the reflection and transmissioncoefficients between two viscoelastic isotropic media，Journal of SeismicExploration，v.9，pp.223-233；Zhi，L.、S.Chen和X.Li，2016，Amplitude variation withangle inversion using the exact Zoeppritz equations-Theory and methodology，Geophysics，v.81，no.2，N1-N15)。相反，其他人则得出结论，即使在小角度的范围内，地震反射AVO响应的非线性也很重要，主要是在引起反射的弹性属性具有足够大的对比的情况下。

3.具有泊松阻抗的二次AVO近似

考虑到入射在两种焊接的半无限各向同性均质介质的水平弹性界面上的稳态平面P波，Aki和Richards所给出的P波反射系数可以表示为：(1)入射角的函数；(2)透射角；以及(3)六个弹性参数，包括一个勘测区域上的上、下介质的速度和密度。对于通过介质传播的平面波，借助于斯涅尔定律，水平慢度或射线参数将保留在原波和模式转换波中。因此，P波反射系数的确切公式可以用透射射线参数p表示为有理函数的形式(Mallick，S.，1993，Asimple approximation to the P-wave reflection coefficient and its implicationin the inversion of amplitude variation with offset data，Geophysics，58，第4期，544-552)：

其中，

A＝(ρ₂q_α1-ρ₁q_α2)(ρ₂q_β1+ρ₁q_β2)，

B＝4(Δμ)²q_α1q_α2q_β1q_β2-4Δμ(p₂q_α1q_β1+ρ₁q_α2q_β2)-(Δρ)²

C＝4(Δμ)²(q_α1q_β1-q_α2q_β2)+4ΔμΔρ，

D＝4(Δμ)²，

E＝(ρ₂q_α1+ρ₁q_α2)(ρ₂q_β1+ρ₁q_β2)，

F＝4(Δμ)²q_α1q_α2q_β1q_β2-4Δμ(ρ₂q_α1q_β1+ρ₁q_α2q_β2)+(Δρ)²

G＝4(Δμ)²(q_α1q_β1+q_α2q_β2)-4ΔμΔρ

其中，ρ₁和ρ₂是跨界面的密度；q_α1、q_α2、q_β1和q_β2是P波和S波的垂直慢度，可以作为射线参数p、P波速度α₁和α₂以及S波速度β₁和β₂的函数来读取。下标数字1和2代表上层和下层。Δρ＝ρ₂-ρ₁定义了密度的对比，而

是剪切模量的对比。

通过在透射射线参数设置为零的点处将泰勒展开式应用于方程(1)，并且省略任何高阶项，可以相对于射线参数p以二次形式表示方程(1)的近似

其中，θ₁和θ₂分别是P波跨界面的入射角和透射角。方程(3)可以理解为流体-流体反射系数R_f，它是在两种介质中的相应横波速度都设为零时两种介质之间的反射系数。

最后，在方程(1)的各种替换和迭代之后，本发明中的泊松阻抗的Zoeppritz方程的二次四项近似计算如下：

方程(4)是在本发明中关于在界面处的P波速度、泊松阻抗和密度的三个参数中的相对对比的所计算的二次近似。

4.非线性AVO泊松阻抗反演

二次AVO近似方程(4)提供了P波反射系数与三个弹性参数(P波速度反射率、泊松阻抗反射率和密度反射率)之间的非线性关系，这三个参数对于在单个或多个入射点处估算岩性至关重要，不论偏倚因素如何。可以使用高斯-牛顿方案来实现所述非线性AVO反演。然而，大多数AVO反演方法是一个不适定的问题，因为某些方法需要适当的平滑或滤波技术以使解决方案唯一且稳定。因此，为了使计算机实现的方法和系统产生准确而有效的结果，有必要将这种非线性问题作为归一化的反演方法加以解决。在此，采用阻尼高斯-牛顿法求解非线性AVO方程，考虑叠前地震数据，并利用约束和分析趋势的组合。例如，由于地震数据的带限性质，将从叠加速度和测井曲线中计算出的低频模型合并到目标函数中，以提高反演的稳定性和准确性(Xia，F.、S.Zhao、W.Ding和Y.Chen，2019年，Constraints guidedbasis pursuit prestack AVA inversion，SEG第89届国际年会，扩展摘要，第585-589页)。此外，可以在目标函数中包含横向约束，该目标函数用于增强相邻迹线之间的空间相关性(Hamid，H.和A.Pidlisecky，2015年，Multitrace impedance inversion with lateralconstraints，Geophysics，第80卷，第6期，M101-M111)。

为了反演P波速度、泊松阻抗和密度，方程(4)自动地重新排列为：

其中，

方程(6)是Zoeppritz方程的线性近似，等效于关于泊松阻抗重新排列Aki-Richards方程。给定一个具有M个角度的叠前角度道集，并且每个迹线包含N个样本，则方程(6)可以重写为矩阵矢量形式：

其中

是各角度处的线性反射率矢量

以及

A_i＝diag[A_i，1 A_i，2 … A_i，N]^T (15)

B_i＝diag[B_i，1 B_i，2 … B_i，N]^T (16)

C_i＝diag[C_i，1 C_i，2 … C_i，N]^T (17)

其中，m_α、m_PI和m_ρ分别是P波速度反射率、泊松阻抗反射率和密度反射率的矢量。

方程(7)是Zoeppritz方程的近似的非线性反射率，包含一个作为三个弹性参数m_α、m_pI和m_ρ的函数的二次项。类似地，方程(7)可重写为矩阵矢量形式，如下黑体所示：

其中

是各角度处的线性反射率矢量

D_i＝diag[D_i，1 D_i，2 … D_i，N]^T (23)

从反演的角度来看，泊松阻抗对背景的扰动将从叠前地震数据中反演。方程(9)中的系数B对于将整体数据失配归因于泊松阻抗起着重要作用。这样，使方程(9)等于零，并且发现对于所有入射角而言，旋转优化因子c都存在奇异性，例如表示为：

c＝α/β (27)

为了简单起见，将方程(5)重写为非线性正向问题，得出：

r＝f(m) (28)

其中，r是P波反射系数，f是包含矩阵A、B、C、D和Ψ的非线性算子，m是包含要进行反演的m_α、m_PI和m_ρ的未知的反射率对比。

有噪声的地震道d可以写为与反射系统r和误差向量e卷积的与给定角度相关的小波W的总和：

d＝Wr+e (29)

合并方程(28)和(29)，得到：

d＝Wf(m)+e (30)

然后，采用高斯-牛顿法来求解方程(30)的非线性AVO反演。假设使用一阶泰勒级数展开将模型更新m表示为初始参考模型m₀和模型扰动Δm的组合，则非线性AVO方程(30)的线性化可重写为：

d＝Wf(m₀)+WFΔm+e (31)

其中，f(m₀)是在参考模型m₀处评估的P波反射系数，F代表弗雷歇导数，即数据d关于模型参数m的偏导数的矩阵，如下所示：

高斯-牛顿法始于对m₀之类的模型的初始猜测，然后进行如下迭代：(1)将数据残差线性化，以预测预测值与当前模型更新附近的正确值之间的距离；以及(2)使用从来自(1)中的数据残差来计算线性最小二乘范数，直到在范数小于或等于用户定义的容差值时实现反演过程中的收敛。因此，高斯-牛顿法将目标函数迭代地近似为模型扰动的二次函数：

其中，

和

分别是数据误差和模型参数的协方差矩阵。

为了平滑或增强求解方程(33)的稳定性，添加了其他约束条件，以促进反演过程的快速收敛，从而将低频约束条件作为惩罚函数纳入目标函数中(Xia，F.等，同上)，然后在反演模型中添加结构约束以改善横向连续性(Hamid，H.和A.Pidlisecky，2015，Multitraceimpedance inversion with lateral constraints，Geophysics，80，no.6，M101-M111)。最小化目标函数(33)并结合模型约束，产生最终的目标函数：

其中，

是低频趋势模型，可以通过测井或叠加速度容易地估算。

方程(34)的最后三项定义了模型空间上的平滑函数。在此，逆协方差矩阵

与数据拟合的加权算符乘以其伴随关系有关。运算符L将模型空间从垂直方向转换为水平方向；协方差逆矩阵

与粗化算子有关，例如，水平方向的二阶导数乘以其伴随。

在高斯-牛顿方法下，为了最小化二次目标函数(34)，就数据d-Wf(m₀)中的剩余误差而言，模型Δm的扰动是线性最小二乘问题。因此，线性最小二乘问题可以通过公认的共轭梯度方法解决(Fletcher，R.和C.Reeves，1964，Function minimization by conjugategradients，The Computer Journal，第7卷，第149-154页)。共轭梯度法中的梯度是目标函数的梯度，通过取目标函数相对于模型m^T的导数来确定，它用于描述函数最迅速降低的方向，表示为：

其中，r是数据残差，定义为r＝d-Wf(m₀)。

对于地质环境中包含鲜明的界面对比的勘测区域，基本追踪法(Chen，S.、D.Donoho和M.Saunders，2001，Atomic decomposition by basis pursuit，Society forIndustrial and Applied Mathematics Review，第43卷，第1期，129-159)或混合l¹/l²范数方法(Li，Y.、Y.Zhang和J.Claerbout，2012年，Hyperbolic estimation of sparse modelsfrom erratic data，Geophysics，第77卷，第1期，第V1-V9页)确保了反射系数的稀疏性和模型参数的块状结构，这比解决l²范数优化问题的平滑模型更可取。

5.井震连接方法

井震连接方法在本领域中也通常称为连井方法，或者就称为井震方法。无论如何称呼，本领域的普通技术人员都应理解，它们在地震解释的标定，在阻抗和AVO反演中使用地震振幅以及最终推论到风险分析中的基础都是至关重要的。它们提供了以下手段：(1)正确识别要拾取的地层，以及(2)改进或更新用于叠前地震反演的与角度相关的角度或小波。

为了计算连井方法，将合成地震图与区域地震道相匹配，并将井的特征与地震数据相关。构建综合模型中的主要概念是卷积模型。这表示地震反射信号是振幅和极性不同但形状相同的干扰反射脉冲的序列。这种脉冲形状是地震小波或角度小波，形式上是目标深度处单位强度的隔离反射器返回的反射波形。因为反射边界在深度上的间隔比反射脉冲的长度要细得多，所以干扰程度通常很严重，只有最强的反射器或反射复合体才能在反射信号中突出。

用于执行从井震方法的几种地球物理软件包在本领域中是已知的。此类产品的示例包括可从HampsonRussell获得的WellTie，它利用了可协同工作的其他软件应用程序，如也可从Hampson Russel获得的Evaluate Wavelets和Well Log Editor。根据这些已知的计算机实现的技术，最佳的角度相关性小波被确定为在与测井反射率卷积时所得到的合成波，其以最小二乘法的方式与输入地震最匹配。然后，通过评估候选角度小波的对输入道执行地震反演的能力，进一步审查该候选角度小波。可以使用各种质量控制来评估其性能。

不管使用哪种软件包，在连井时都将更新角度相关性小波，以与测井结果进行交互，从而根据井眼生成地震的综合表示。这是通过以下步骤开始的：(1)在井位处获得测得的地震图；(2)从测井中获取并计算时深关系；(3)实时计算反射系数；(4)通过将反射系数级数与地震角相关性小波卷积来制作合成地震图；(5)将合成生成的地震图与井位处的测得地震图进行关联。

发明内容

通常，勘探和储层表征是在勘测区域上针对其土壤和流体潜在特性进行的。取决于在勘测区域中发现的特性，可以建立一个或多个烃(即，石油和天然气)储层。一旦在勘测区域内建立了这些储层或井，就可以使用计算机实现的方法来自动地解释地震结果。这些计算机实施的方法中的许多方法都使用反演模型，其包括一种分析方法，其中将与来自地球地下层之间的反射界面处的声能反射相对应的深度域信号转换为代表地层的物理属性的一条或多条道。因此，与之相反，本发明的一个目的是在地震勘探和储层表征的领域中使用非线性直接叠前模型来提供一种用于岩性识别和流体鉴定的计算机实现的方法。

本发明的基于反射率的计算机实现的反演方法和系统是在时域中在地震数据的振幅取决于反射系数的假设下开发的，这意味着在反演之前地震数据已经完全通过先进的叠前偏移技术而针对波传播效应(如传输损耗、层间多次散射和模式转换)进行了校正。

因此，本发明的一个目的是将其实现为自动计算机方法，用于估计泊松阻抗，以避免在具有至少一个烃井或储层位置的勘测区域中的因间接计算所带来的不确定性、传播和累积误差。

在本发明的一个实施例中，该计算机实现的方法检索在时域中表示的经放大的测井数据，包括在勘测区域的各个入射点上的纵波速度、横波速度和密度数据。该方法还从同一勘测区域中检索角度图像道集的集合，其经过预处理以保留各种入射角的信号幅度信息；并且还检索包含地层信息的地球模型数据，以及一段时间内的地震速度数据。然后，该方法通过插值来生成各种低频模型，这些模型基于检索到的地球模型数据和经放大的测井数据分组为集合，并且还生成勘测区域上的角度小波的集合。然后，该方法对角度图像道集的集合和经放大的测井数据应用井震连接方法，以便更新角度小波的集合。使用取决于所检索到的数据点的计算，该方法通过将角度图像道集的集合中的第一角度图像道集和低频模型的集合中的相应的第一低频模型与角度小波的集合进行组合，以计算非线性直接叠前地震反演模型。该计算机实现的方法自动地重复针对勘测区域的选择和计算的步骤，直到角度图像道集的集合中的所有角度图像道集与低频模型的集合中的所有相应的低频模型都已经组合到角度小波的集合中为止。作为质量控制方法，该计算机实现的方法验证在自动重复步骤期间是否通过将角度图像道集的集合中的最终角度图像道集和低频模型集合中的相应的最终低频模型与角度小波的集合相结合而执行了最终的非线性直接叠前地震反演模型。最后，通过计算弹性反射率的指数积分，该方法从所述非线性直接叠前地震反演模型的结果中产生勘测区域的弹性属性的最终模型。

由于前述原因，需要开发一种计算机实现的方法和系统，其在使用泊松阻抗的同时计算非线性直接叠前地震反演模型，以避免因使用间接计算的其他方法和系统导致的不确定性、传播和累积误差。所提出的实施例考虑了在预临界入射角处的反射响应，并且不受弱对比界面或数据的偏置部分的假设的影响。实际上，通过将其与经典的Aki-Richards方程(也参见图9)以及与精确的Zoeppritz方程进行比较，对所提出的计算机实现的方法的准确性进行了检验，发现所获得的反演结果产生了针对泊松阻抗的更稳定和更精确的估计。此外，测试结果表明，使用了非线性直接叠前反演的计算机实现的方法和系统的本发明对于其他的实际应用也具有很大的潜力。

附图说明

通过结合附图并考虑以下详细描述，可以容易地理解本发明的教导。

图1是示出了根据本公开的一个实施例的具有地震源的多个入射点的勘测区域的平面示意图；

图2是根据本公开的一些实施例的说明性环境的剖视示意图，该环境具有地震源的入射点、地震检波点、井位、井眼、多个透射射线，以及多个入射角；

图3是示出了根据本公开的一个实施例的计算机实现的方法的流程图，该方法采用了非线性直接叠前地震泊松阻抗反演；

图4是示出了根据本公开的一个实施例的用于生成角度图像道集的集合的方法的流程图；

图5是示出了根据本公开的一个实施例的用于生成低频模型的集合的操作的流程图；

图6是示出了根据本公开的一个实施例的用于生成角度小波的集合的操作的流程图；

图7以图形方式显示了根据本公开的一个实施例的由用于生成角度图像道集的集合的步骤所产生的结果(图7A)、由用于生成低频模型的集合的步骤所产生的结果(图7B)，以及由用于生成角度相关性小波的集合的步骤所产生的结果生(图7C)；

图8以图形方式显示了根据本公开的一个实施例的井震方法的测试结果，用于比较井位处的合成地震数据与测量地震数据；

图9以图形方式显示了图9A所示的根据本公开的一个实施例的使用本发明的计算机实现的方法和系统的反演泊松阻抗(PI)和图9B所示的使用间接的Aki-Richards近似的反演泊松阻抗之间的比较；

图10以图形方式显示了根据本公开的一个实施例的由本发明的计算机实现的方法和系统所产生的弹性属性的最终模型，其中图10A示出了P波速度，图10B示出了泊松阻抗，图10C示出了密度；

图11是根据本公开的一个实施例的用于执行该计算机实现的方法的数字高性能计算系统的电气图，以框图的形式显示；和

图12是根据本公开的一个实施例的流程图，示出了由高性能数字计算系统的系统计算机所执行的用于生成非线性直接叠前地震反演模型的子例程。

具体实施方式

现在将详细参考本公开的几个实施例，其示例在附图中示出。应注意的是，在适用之处，附图中使用相同或相似的附图标记来相同或相似的功能。这些附图仅出于说明的目的描绘了本公开的实施例。本领域的技术人员将从以下描述中容易地认识到，在不脱离本公开所描述的原理的情况下，可以采用本文所示的结构、系统和方法的替代实施例。

图1示出了地震勘测区域101，其可适用本发明的优选实施例。重要的是应注意，图1的勘测区域位于陆地区域102。本领域的普通技术人员可以认识到，地震勘测区域可产生局部地质的详细图像，以确定可能的烃(油气)储层的位置和大小，以及井位103。在这些勘测区域中，声波在冲击过程中在各个入射点104处从地下岩层中弹起，而反射回地面的波被地震数据记录传感器105捕获，并由数据传输系统1101从所述传感器105中进行无线传输1102，然后存储起来以供以后处理，并由图11的数字高性能计算系统进行分析。

在图2中，勘测区域中的一部分地层的截面图被示为201，其显示了不同类型的地层102、203、204，其包括本发明中的0地震勘测数据。特别地，本领域的普通技术人员将很快认识到，本示例示出了共中心点道集，其中地震数据道通过表面几何形状分类，以近似地层中的单个反射点。在此示例中，来自多个炮点和检波点的数据可以组合为单个图像道集，也可以根据要进行的分析的类型而单独地使用。尽管本示例可以示出平面反射器和相应的图像道集类别，但是可以使用本领域中已知的其他类型或类别的图像道集，并且其选择可以取决于各种地层状况或事件的存在。

如图2所示，来自多个入射点104的地震能量将从不同地层之间的界面处反射。这些反射将由多个地震数据记录传感器105(各传感器放在相互不同的位置偏移210处)和井103捕获。因为所有入射点104和所有地震数据记录传感器105都位于不同的偏移210处，地震勘探数据或道(在本领域中也被称为道集)将以不同的入射角208来记录。入射点104在地层中产生向下的透射射线205，它们通过其向上的透射反射而被记录传感器105捕获。在该示例中，井位103被示为具有连接到井眼209的现有钻井，沿井眼209可使用本领域中已知的技术获得多个测量值。该井眼209用于获得测井数据，其中包括P波速度、S波速度和密度等。在勘测区域内还布置了图2中未示出的其他传感器，以便还捕获翻译工具和本领域的普通技术人员执行各种地球物理分析所需的地层数据信息。在本示例中，将从实地记录中对这些道集进行分类，以便检查振幅、信噪比、偏移、频率组成、相位和其他地震属性对于入射角208、偏移量测量210、方位角和对于数据处理和成像来说重要的为本领域的普通技术员所知的其它地质属性的依赖性。

然而，为了执行弹性反射率分析，在本领域中通常有用的是利用地震数据的冗余，并产生具有比物理空间的两个坐标更大的尺寸的图像。在本发明中，弹性属性的最终模型314是通过计算主要在时域中表示的各种数据而制成的三维立方体。当创建弹性属性的最终模型314时，没有包括其他的尺寸，但在其他的情况下，这些尺寸可以包括绝对偏移量和方位角。在本发明的优选实施例中，由于弹性属性的最终模型314是根据三维坐标创建的，所以该立方体将最终包括主测线方向(in-line)、横测线方向(cross-line)和时间输入变量。

参照图3，其示出了根据本发明的优选实施例的总体流程图。检索阶段301启动该过程，并且包括从勘测区域中获得并准备数据和信息。特别地，从勘测区域的地震数据中检索出四种不同类型的输入：在时域中表示的放大后的测井数据302，角度图像道集的集合310，地层信息304，以及地震速度数据305。测井数据可以直接从井位获得，并远程传输到数据库1103中以进行进一步的编译，如图11中所设想的那样。它也可以由系统计算机1105从数据库1103中检索。对于本发明来说，该测井数据包括勘测区域上的P波速度、S波速度和密度。此外，角度图像道集的集合310(进一步由图7A示例)包括多个图像道集，其对于在不同入射角208处的在时域中表示的勘测区域101中的入射点104的数量来说是恒定的。尽管如此，由于该地震勘探数据402太原始、噪声太多或来自多个入射点104，因此它有待进一步处理。角度图像道集的集合的这种进一步细化是在图4中示出的本发明的计算机实现方法中的子例程过程303。该过程303检索401勘测区域101的地震勘测数据402和地震速度数据305，所有这些都将被存储在数据库(例如存储设备1104)中，并由系统计算机1105从中检索以进行进一步处理。该子例程使用图11的高性能计算执行叠前时间迁移方法403，该方法解决了具有不同叠加速度和数据的冲突倾角的问题。然后，生成404偏移图像道集，进行叠加405，并使用最初为炮检偏移而引入的径向道映射转换回时域406。针对勘测区域101上的各个入射点104生成310角度图像道集的最终集合。图7A示出了本发明中的包括角度图像道集的集合的一个例子。此外，图11的高性能计算系统设想了将角度图像道集数据的这些集合310存储在数据库(例如存储设备1104)中，并由系统计算机1105从中检索。

继续图3的计算机实现的方法，尤其是所检索的地层信息数据304可以被远程传输到数据库1103上，以进行进一步的编译，例如图11中所设想的那样。此外，数据304也可以由系统计算机1105从数据库1103中检索。出于本发明的目的，所检索的地层信息数据304通常表示为地球中的层、表面和界面，及其在地震数据中的表现形式。在其最纯粹的形式中，数据304是无法实现的，因为本领域的技术人员很少能确切知道地层在地下的位置。然而，地层信息数据304包括三维平面(x，y，z)中的离散样本的矩阵，其可以存储在ASCII文件中。因此，地层是可以在地图上明确绘制的东西，并且可以像光栅图像一样处理。一些软件工具甚至调用属性映射地层，进一步模糊了其定义，并因此将其创建为软件的人工产物。

本发明收集的另一组数据是地震速度305。地震速度(纵波速度和横波速度)305是材料的基本性质，该性质随外部(应力，温度)和内部(流体饱和度，裂纹密度)的条件的变化而变化。监测这些外部或内部条件的变化是地球物理研究的目标，例如地震预测(通过应力变化监测)和油藏开采(通过流体饱和度监测)。

一旦已经检索到所有数据变量(302、310、304和305)，该计算机实现的方法便使用系统计算机1105继续执行两个子例程。

使用系统计算机1105的第一子例程涉及通过基于地球模型数据和勘测区域101的各个入射点104上的放大后的测井数据的插值而生成低频模型307的集合。然后，第一子例程307输出低频模型的集合309(如图7B更好地显示)，并将其存储在设备1104上。在此之后执行子例程306，其包括通过井震连接方法605对角度图像道集312的集合以及放大后的测井数据302进行计算，生成来自勘测区域101的多个入射角208的角度相关性小波308的集合(由图7C示例性示出)。

图5进一步描述了用于生成低频模型的集合的第一子例程307。通过执行针对在时域302中表示的经放大的测井数据的另一检索操作500，并且与经解释的地层信息数据304和地震速度数据305相结合，本发明生成在勘测区域101的多个入射点204上的低频模型的集合(图7B进一步示出)，其使用插值处理506，也基于302和305。然而，由于如本发明中所描述的那样从地质传感器206中检索到的地震数据是带限的，所以检索到的信息缺乏对于相对-绝对属性转换而言必不可少的低频。这些未滤波的或有噪声的低频模型的集合可以被认为是数据的平均趋势，而较高的频率是与该平均趋势的偏差。因此，本发明对302和305进行滤波或平滑，以生成低频模型的集合(图7B)，它是本发明的关键步骤，因为它们用于解释地震数据下的丢失信息。此外，本发明解决了如何在不同的入射点205之间对测井进行插值的问题。

因为解释器通常不确定如何获得如图7B所示的低频模型的集合，并且常常对于结果不满意，所以本发明公开了一种与本发明的优选实施例兼容并且能产生更好的结果的方法。为了获得所述更好的结果，本发明公开了生成低频模型的集合的步骤，其进一步包括：通过在一组低通滤波器上运行所检索的数据，从而对所检索的地震速度数据305进行数字化平滑501，并且对经放大的测井数据302进行数字化平滑502。该滤波计算幅度谱以满足频域中的规格，从而切出高频，并显著降低来自所检索的经放大的测井数据302以及所检索的地震速度305中的噪声。一旦将地震速度305平滑到503，就使用所检索的经解释的地层信息304并遵循Dix变换反演方法(Durbaum，1954；Dix，1955)将其转换为内部速度数据504。所述方法主要依赖于通过众所周知的下述公式根据拾取的均方根(rms)或堆叠的速度以及相应的行进时间来计算层速度：

其中，

是时间间隔Δt_i＝t_i-t_i-1上的局部内部层速度，其中i是从1开始到第n个值的迭代器。参数

和

是层的顶部和底部界面处的均方根速度。瞬时速度

假定为分段常数，

并且在界面处具有不连续性。

然后将转换后的内部速度504与所解释的地层信息304和平滑的放大后的测井数据505相结合，以执行插值步骤506。由于数据即使经过平滑处理后输入数据仍具有噪声并需要进一步完善的事实，因此需要执行步骤506。这样，本发明的计算机实现的方法使用在本领域中称为cokriging的插值方法。该方法的特殊特征以及为什么在本发明中具有相关用途，是因为它可以利用不同性质的数据来对特定属性进行建模、优化和插值。一般而言，本发明的计算机实现的方法将一种算法应用于编码为可接受s组变量的插值程序(与典型的cokriging算法中的仅2组变量相反)，产生了遵循以下方程的线性系统：

其中，C是所有已知变量对的协方差(或其估计)矩阵，C₀是未知变量和所有其他变量之间的成对协方差的矢量。其中，μ是所有拉格朗日乘子μ₁到μ_s的矢量。E是矩阵I₁到I_s的矢量。每个矩阵I_i(i∈{1…s})的大小为(第i个变量集中的点数)x s。I_i的第i列中的所有元素均为1，所有其他项均为零。其中，T是所有系数的矢量，I₀是方程右侧C₀下所有元素的的列矢量，大小为s x 1。

一旦本发明的计算机实现的方法已按照上述算法对304、504和505的所有变量进行插值，它将生成低频模型的集合309，包括各种单独的低频模型，其数量与角度图像道集的数量相同。

类似地，子例程306涉及使用系统计算机1105从勘测区域101上的多个入射角208处生成角度小波的集合(在图7C中示出)的步骤，并且由图6表示。完全执行该子例程所需的步骤还包括检索601角度图像道集的集合310，以及经放大的测井数据302。然后，从检索到的角度图像道集312中提取角度小波的初始集合602，其用作针对最终的角度相关性小波应看起来像什么的初始估计。在本发明中使用统计小波提取程序来计算提取，该程序仅使用地震道来提取角度相关性小波604。然而，为了提取每个角度小波604，使用了角度图像道集的集合310中的对应道。然后，提取一个分析窗口，并使所提取的分析窗口的开始和结尾逐渐变细，其长度等于10个样本中的较小者，或分析窗口的1/4。在数据窗口中计算自相关过程，其中，该自相关的长度等于所需波长的1/2。计算自相关的振幅谱，并计算自相关的平方根。执行该操作以近似小波的振幅谱。此后，添加期望的相位(例如0°，5°，10°，15°等)，并且通过1105在602处计算快速傅立叶逆变换，以产生或提取初始角度小波604。提取过程602验证是否之前还产生了其他小波，并且生成初始小波。如果创建了其他小波，则从相同角度图像道集312的其他道中在先前生成的小波旁边生成新的小波。提取过程继续进行，直到已经分析了来自角度图像道集312的所有角度为止。一旦所有角度小波已经被统计性地提取，就生成604初始角度小波的提取集合。

根据本发明的优选实施例，子例程306中的下一步是使用方程(38)的二次近似来计算角度相关性反射率603，其将θ_i表示的入射角208相对于勘测区域101的以下参数进行结合：(a)放大的测井数据302；(b)P波速度α；(c)S波速度β；(d)泊松阻抗PI；以及(e)密度ρ。

接下来，通过对所提取的用于勘测区域的角度小波的初始集合604执行井震方法605(在图8中图示)，生成角度相关性小波的集合308，并将其存储在设备1104上。所述角度相关性小波的集合在图7A中进一步示例。

根据本发明的优选实施例，一旦生成角度图像道集的集合310、角度小波的集合308和低频模型的集合309，则接下来执行过程311。该过程使用系统计算机1105将三个不同的数据集合组合为子例程312的前身，并将其存储在设备1104中。

然后，通过由用户经由用户输入设备(如键盘1108和鼠标1109)的组合来输入1201用户定义的容差值1202，从而启动子例程312(图12所示)；本领域的普通技术人员可以理解，该容差值小于1(即10^-6)。用户还输入包括变量m₀的初始估计模型1203，本领域的普通技术人员也将其视为初始模型预测。然后，从用户计算机系统发送这两个输入，以存储在数据库设备1104中，以供系统计算机1105在步骤1204之后进行处理。过程1204使用初始估计模型1203和用户定义的容差值1202来计算一组初始数据残差。在地球物理学术语中，数据残差是由测量数据和预测数据之间的差得出的。在此过程中，使用以下表达式计算所计算出的初始数据残差，这意味着本领域的技术人员将很容易将其识别为从其他已知数据残差表达式中得出的术语，用于计算观测数据和预测数据之间的数据残差：

r＝d-Wf(m₀) (39)

此后，本发明的计算机实现的方法开始准备系统计算机1105运行循环系统，以获得最终的非线性直接叠前地震反演313。但是，用户必须首先根据值m₀到m_i在个人计算机中设置初始估计模型1205。这是过渡计算步骤，向系统计算机1105发送消息：如在系统计算机1154中所定义的那样，模型m_i将被重复，直到满足特定要求为止。根据本发明的优选实施例，表达式m_i中的i是存储在设备1104中的迭代器，被系统计算机1105用作指针，用于抽象不同结构(即多个线性对象)的元素的地址。结果，i允许数据结构与本发明的算法进行接口，而本发明的算法并不知道它们所要进行操作的结构。此外，所述迭代器由系统计算机1105使用，以基于其当前可用的迭代器来生成新的迭代器，并且在将其通过各种所需的过程(参见1206、1207、1208、1209和1210)之后返回到模型变量m的序列中的下一个迭代器。

当所有值都已经输入到系统计算机1105，以及已经计算出的初始数据残差1204时，使用下述表达式来计算1206一组模型梯度：

图11的数字计算系统使用系统计算机1105，针对Zoeppritz方程的近似的线性和非线性反射率计算一组数据残差变化1207，还包括以下表达式：

Δr＝WFg(m_i) (41)

根据本发明的优选实施例，更新输入估计模型1208的过程由系统计算机1105执行，其计算从任何先前模型(即m₀)到当前模型(m_i)的变化，并将其添加到当前模型。然后生成更新的输入模型m_i+1，并且迭代器的值加1。

接下来执行使用表达式r＝d-Wf(m_i+1)来计算数据残差1209的步骤。在完成步骤1209之后，系统计算机1105为所计算的数据残差1209生成矩阵范数1210。然后，系统计算机1105使用范数计算的数据残差1209与所输入的用户定义的容差值1202来执行一系列迭代1211，直到确认所计算出的数据残差1209小于容差值1202。一旦系统计算机1209确认了这种不等关系，就结束迭代过程1211，并以最终的更新模型的集合的形式生成1212反演结果，并将其和在子例程312期间生成的整个数据集合输出到存储设备1106。然后，通过选择导出功能，用户可以使用个人计算机1107检索数据。

在子程序312结束后，图3所示的本发明的计算机实现的方法执行迭代313，其中它重复将低频模型的集合309与角度相关性小波的集合308相组合的步骤，直到来自角度图像道集的集合310的所有道集以及来自低频模型的集合309中的所有相应的低频模型都已经被组合到角度小波的集合308中为止。然后，系统计算机1105验证重复步骤是否执行了最终的非线性直接叠前地震反演模型，换句话说，不再有更多的最终角度图像道集(来自角度图像道集的集合)与相应的最终低频模型(来自低频模型集合)组合到角度小波的集合中。如果系统计算机1105验证313已经组合了所有模型的集合，则生成弹性属性的最终模型314，如图10所示。如果系统计算机1105不能验证313将要产生弹性属性的最终模型314，它将再次启动组合所有模型和数据集合的步骤以及子例程312，直到能够产生弹性属性的最终模型313为止，然后它将由系统计算机1105输出到存储设备1106上。然后，个人计算机1107的用户可以检索并显示弹性属性的最终模型314，以进行进一步分析或建议。

对于任何前述实施例，计算机系统1107可以单独地或彼此组合地包括以下任何一个部件。系统计算机1105执行的功能包括输入、检索、转换、组合、生成、产生、校准、执行、提取、过滤和输出的功能；一系列不同的值，包括测井数据、地层信息、地震速度、P波速度、泊松阻抗、密度、反射率，角度图像道集、低频模型、角度相关性小波，以及更复杂的功能，例如基于与包含至少一个井位和穿透地下地层的井眼的勘测区域相关的信息和从中获得的信息来执行本发明的非线性直接叠前地震反演模型。

参照图7，它以图形的形式示出了作为前述子例程(303、306和307)的输出而生成的不同的数据集(310、308和309)。同样，图8以图形的形式示出了通过使用角度相关性小波和所计算出的角度相关性反射率来将合成地震图与勘测地震数据进行比较而执行井连方法的过程。

参照图9，其以图形方式显示了使用本发明的计算机实现的方法(采用非线性直接叠前地震泊松阻抗反演，图9A示出)和另一计算机实现的方法(使用间接Aki-Richards近似，图9B示出)之间的比较。在两种图形表示中，线901代表真实模型，线902代表从各计算机实现的方法中获得的结果。显然，图9A示出了比图9B更准确的结果，由线901和线902之间的差异产生的间隙显示。

在图10中示出了计算机实现的弹性属性的最终模型315的图形表示。当根据它们的各个表示来划分参数时，图10A示出了P波速度，图10B示出了泊松阻抗，图10C示出了密度。线1001代表真实参数，线1002代表初始参数，而线1003代表在由本发明的弹性属性的最终模型315中产生的参数。正如所期望的那样，由本发明产生的参数与速度、泊松阻抗和密度的真实值吻合得很好。

如本文所使用的那样，术语“计算”包括各种各样的动作，包括计算、确定、处理、推导、勘测、查找(例如，在表、数据库或另一数据结构中查找)、确定，等等。它还可以包括接收(例如接收信息)、访问(例如访问存储器中的数据)，等等。而且，“计算”可以包括解析、选择、建立，等等。

根据本发明的优选实施例，仅作为示例实施例详细介绍了某些硬件和软件描述，但这并不限制所公开实施例的实现的结构。例如，尽管已经描述了图12所示的计算机系统的许多内部和外部部件，本领域的普通技术人员将理解，这样的部件及其互连是众所周知的。另外，所公开的发明的某些方面可以体现为使用一个或多个处理单元1105中执行的软件。该技术的程序方面通常被视为可执行形式的“产品”或“制品”，和/或在一种类型的机器可读介质上承载或体现的关联数据。有形的非暂时性“存储”类型的介质和设备包括用于计算机、进程等的任何或所有存储器或其他存储元件，或其相关模块，例如各种半导体存储器、磁带驱动器、磁盘驱动器、光盘或磁盘，以及可以随时为软件编程提供存储的器件。

另外，附图中的流程图和框图示出了根据本公开的各种实施例的系统、方法和计算机程序产品的可能实现的架构、功能和操作。还应当注意，在一些替代实施方式中，方框中指出的功能可以不按附图中指出的顺序发生。例如，取决于所涉及的功能，实际上可以基本上同时执行连续示出的两个方框，或者有时可以以相反的顺序执行这些方框。还应注意，框图和/或流程图图示的每个方框以及框图和/或流程图图示中的方框的组合可以由执行指定硬件功能或动作的基于专用硬件的系统来实现，或由专用硬件和计算机指令的组合来实现。

尽管在前面的说明书中已经针对本公开的某些优选实施例描述了本公开，并且出于说明性的目的已经阐述了许多细节，但是本发明不应不当地限于已经出于说明性目的而阐述的前述内容。相反，对于本领域的技术人员来说明显的是，在不背离如所附权利要求限定的本发明的真实范围的前提下，很多种修改和替代实施例将是显而易见的。另外，应当理解，在本文的任何一个实施例中示出或描述的结构特征或方法步骤也可以在其他实施例中使用。

符号表

Claims (17)

1.一种计算机实现的方法，该方法使用非线性直接叠前地震反演模型来估计泊松阻抗，并避免在具有至少一个井位的勘测区域中因间接计算产生的不确定性和累积误差，该方法包括：

检索在时域中表示的经放大的测井数据，包括所述勘测区域中的各个入射点处的纵波速度、横波速度和密度数据；

产生角度图像道集的集合，对其进行预处理以保留所述勘测区域的各个入射角的信号幅度信息；

检索包括针对所述勘测区域的地层信息数据的地球模型数据；

检索所述勘测区域的一段时间上的地震速度数据；

通过基于所述地球模型数据和所述勘测区域的各个入射点处的经放大的测井数据来进行插值，从而生成低频模型的集合；

通过对所述角度图像道集的集合和所述经放大的测井数据进行井震连接计算，生成来自于所述勘测区域的多个入射角的角度相关性小波的集合；

将所述角度图像道集的集合中的第一角度图像道集和所述低频模型的集合中的对应的第一低频模型与所述角度小波的集合相组合；

通过使用组合的第一角度图像道集、第一低频模型和角度小波的集合来执行非线性直接叠前地震反演模型；

针对所述勘测区域重复组合和执行的步骤，直到将所述角度图像道集的集合中的所有道集和所述低频模型的集合中的所有对应的低频模型都已经与所述角度小波的集合相组合为止；

验证重复步骤是否通过将所述角度图像道集的集合中的最终角度图像道集和所述低频模型的集合中的对应的最终低频模型与所述角度小波的集合相组合而执行了最终的非线性直接叠前地震反演模型；和

通过计算弹性反射率的指数积分，从所述非线性直接叠前地震反演模型的结果中产生所述勘测区域中的弹性属性的最终模型。

2.根据权利要求1所述的计算机实现的方法，其特征在于，生成所述角度图像道集的集合步骤包括以下步骤：

检索勘测区域的地震勘测数据；

检索勘测区域的地震速度数据；

使用针对勘测区域的检索到的地震勘测数据和检索到的地震速度数据来执行叠前时间偏移法；

产生偏移图像道集；

在不同的入射偏移处堆叠一段时间上所检索到的偏移图像道集；

将偏移图像道集转换为角度图像道集；和

产生角度图像道集的集合。

3.根据权利要求1所述的计算机实现的方法，其特征在于，生成所述低频模型的集合的步骤还包括以下步骤：

检索勘测区域的地震速度数据；

从勘测区域的地球模型数据中检索经解释的地层信息；

检索在时域中表示的经放大的测井数据，包括所述勘测区域中的各个入射点处的纵波速度、横波速度和密度数据；

对检索到的地震速度数据和检索到的经放大的测井数据进行平滑化；

使用检索到的经解释的地层信息来将平滑化的地震速度数据转换为层速度数据；

使用经解释的地层信息和平滑化的经放大的测井数据来对经转换的层速度数据进行插值；和

产生低频模型的集合。

4.根据权利要求1所述的计算机实现的方法，其特征在于，生成来自所述勘测区域的多个入射角的角度小波的集合的步骤包括以下步骤：

检索勘测区域的角度图像道集；

检索勘测区域的经放大的测井数据；

从检索到的角度图像道集中提取角度小波的初始集合；

从检索到的经放大的测井数据中计算在时域中表示的反射率；

对勘测区域的角度小波的初始集合执行井震方法；和

生成角度相关性小波的集合。

5.根据权利要求1所述的计算机实现的方法，其特征在于，产生弹性属性的最终模型的步骤包括以下步骤：

对所有非线性直接叠前地震反演模型的弹性反射率进行反演，包括表达式：d＝Wf(m₀)+e；

产生纵波速度、泊松阻抗和密度的体积模型；

将体积模型输出到具有存储器资源的计算机系统；和

产生弹性属性的最终模型。

6.根据权利要求1所述的计算机实现的方法，其特征在于，所述低频模型的集合在数量上等于所述角度图像道集的集合中的相同量的角度图像道集。

7.一种数字计算系统，用于执行具有至少一个井位的与多个入射点相关联的地震勘测区域的非线性直接叠前地震反演模型，以便从离多个入射点为多个偏移距离的位置处的多个地震数据记录器上获得地震勘测数据，所述数字计算系统包括：

一组存储资源，用于存储与具有至少一个井位的勘测区域相对应的数据，所述数据对应于：经放大的测井、地球模型、角度图像道集、低频模型、角度小波、非线性直接叠前地震反演模型，以及弹性属性的最终模型；

计算机系统输出设备；和

系统计算机，其耦合到存储器资源和计算机系统输出设备，用于执行生成非线性直接叠前地震反演模型的操作。

8.根据权利要求7所述的数字计算系统，其特征在于，所述系统计算机还被编程为执行以下操作：

输入用户定义的容差值；

输入等于用户输入值或零的初始估计模型；

使用所述初始估计模型计算一组初始数据残差；

建立等于输入估计模型的初始估计模型；

计算一组模型梯度；

计算一组数据残差变化；

使用所述一组模型梯度和所述一组数据残差来更新所述输入估计模型；

使用更新后的输入估计模型来计算一组数据残差；

根据所计算的一组数据残差来生成矩阵范数；

重复以下步骤：计算一组模型梯度，计算一组数据残差变化，更新输入估计模型，以及计算一组数据残差，直到生成的范数等于用户定义的容差值为止；和

生成最后一组数据残差。

9.根据权利要求8所述的数字计算系统，其特征在于，输入初始估计模型还包括表达式m₀。

10.根据权利要求8所述的数字计算系统，其特征在于，计算一组初始数据残差的步骤还包括表达式：r＝d-Wf(m₀)。

11.根据权利要求8所述的数字计算系统，其特征在于，将初始估计模型设置为等于输入估计模型还包括表达式m₀＝m_i。

12.根据权利要求8所述的数字计算系统，其特征在于，计算一组模型梯度的步骤还包括以下表达式：

13.根据权利要求8所述的数字计算系统，其特征在于，计算一组数据残差变化的步骤还包括表达式：Δr＝WFg(m_i)。

14.根据权利要求13所述的数字计算系统，其特征在于，所述项F还包括以下表达式：

15.根据权利要求8所述的数字计算系统，其特征在于，更新输入估计模型的步骤还包括表达式：m_i+1＝m_i+Δm。

16.根据权利要求8所述的数字计算系统，其特征在于，计算一组数据残差的步骤还包括以下表达式：r＝d-Wf(m_i+1)。

17.根据权利要求8所述的数字计算系统，其特征在于，每执行一次所述重复步骤，所述重复步骤就将迭代器i加1，直到所生成的范数小于或等于所述用户定义的容差值为止。

Images