基于模型的时间保持不变的层析成像

摘要

使用时间保持不变的层析成像进行地震数据建模的系统和方法，包括存储表征初始地震数据模型的参数值的初始集。初始地震模型对应于至少两个或多个射线对。每个射线对可具有走时。通过对初始模型中的两个或更多射线对中的每个射线对改变参数值的初始集中的两个或更多个参数值，产生改变后的模型。改变所述两个或更多个参数值中的一个参数值而不改变其余参数值可对应于每个射线对的走时变化，而组合地改变所述两个或更多个参数值通常对应于每个射线对的走时没有净变化。

权利要求

1.一种使用时间保持不变的层析成像进行地震数据建模的方法，所述方法包括：存储表征初始地震数据模型的参数值的初始集，其中所述初始地震模型对应于至少多个射线对，每个射线对包括走时；以及通过对所述初始模型中的每个射线对改变所述参数值的初始集中的两个或更多个参数值，产生改变后的模型，其中改变一个参数值而不改变所述两个或更多个参数值中的其余参数值对应于每个射线对的走时变化，而组合地改变所述两个或更多个参数值对应于每个射线对的走时没有净变化。

2.权利要求1的方法，其中被改变的所述两个或更多个参数值中的至少一个表示所述初始模型的一个或更多个构造层位位置的改变。

3.权利要求1的方法，其中被改变的所述两个或更多个参数值中的至少一个表示所述初始模型的一个或更多个区域中介质性质的改变。

4.权利要求1的方法，其中使用所述时间保持不变的层析成像来在保持沿每个射线对的总走时不变的情况下计算用于模拟所述初始模型的不同扰动的层析成像系数。

5.权利要求1的方法，其中所述参数值的初始集以矩阵形式表示。

6.权利要求5的方法，其中产生改变后的模型包括对模型应用矢量。

7.权利要求1的方法，包括接受用户输入和响应用户输入来改变所述两个或更多个参数值。

8.权利要求1的方法，包括将所述改变后的模型存储在存储器中。

9.一种使用时间保持不变的层析成像进行地震数据建模的方法，所述方法包括：存储表征地震数据集的初始地震模型的参数值的初始集，其中所述初始模型具有大小，且其中所述初始地震模型对应于至少多个射线对，每个射线对包括走时；

改变所述参数值的初始集中的两个或更多个参数值，其中被改变的参数值中的一个对应于每个射线对的走时的改变，所述走时的改变与其余被改变的参数值所对应的改变大小相等符号相反，且其中所述两个或更多个参数值以互补方式被改变，使得在每个射线对的走时中没有净改变；以及产生至少是部分地通过所述两个或更多个改变后的参数值表示的地震数据集的改变后的模型，其中所述改变后的模型具有与所述初始地震模型的大小相同的大小。

10.权利要求9的方法，其中每个初始模型包括第一组维度且所述改变后的模型包括第二组维度，其中所述第一组维度和所述第二组维度的大小相同。

11.权利要求9的方法，其中被改变的所述两个或更多个参数值中的至少一个表示所述初始模型的一个或更多个构造层位位置的改变。

12.权利要求9的方法，其中被改变的所述两个或更多个参数值中的至少一个表示所述初始模型的一个或更多个区域中介质性质的改变。

13.权利要求9的方法，其中使用所述时间保持不变的层析成像来在保持沿每个射线对的总走时不变的情况下计算用于模拟所述初始模型的不同扰动的层析成像系数。

14.权利要求9的方法，包括将所述改变后的模型存储在存储器中。

15.一种用于为地震数据建模的方法，所述方法包括：

产生地震数据集的模型，其中所述模型对应于至少多个射线对，每个射线对包括走时；以及改变所述模型的第一参数值，造成每个射线对走时的第一非零变化；

改变所述模型的第二参数值，造成每个射线对走时的第二非零变化，其中每个射线对走时的所述第一非零变化和所述第二非零变化大小相等符号相反，从而组合地改变所述第一参数值和所述第二参数值没有造成沿每个射线对的走时的总体变化。

16.权利要求15的方法，其中改变所述模型的参数值没有改变所述模型的大小。

17.权利要求15的方法，其中被改变的所述第一参数值和所述第二参数值中的至少一个表示所述初始模型的一个或更多个构造层位位置的改变。

18.权利要求15的方法，其中被改变的所述第一参数值和所述第二参数值中的至少一个表示所述初始模型的一个或更多个区域中介质性质的改变。

19.权利要求15的方法，其中使用所述时间保持不变的层析成像来在保持沿每个射线对的总走时不变的情况下计算用于模拟所述初始模型的不同扰动的层析成像系数。

20.权利要求15的方法，包括将所述改变后的模型存储在存储器中。

21.一种使用时间保持不变的层析成像进行地震数据建模的系统，所述系统包括：存储器，用于存储表征初始地震数据模型的参数值的初始集，其中所述初始地震模型对应于至少多个射线对，每个射线对包括走时；以及处理器，用于通过对所述初始模型中的每个射线对改变所述参数值的初始集中的两个或更多个参数值，产生改变后的模型，其中当所述处理器改变所述两个或更多个参数值中的每个而不改变参数值中的其余参数值时，产生每个射线对的对应走时变化，且其中当所述处理器组合地改变所述两个或更多个参数值时，每个射线对的走时没有对应的净变化产生。

22.权利要求21的系统，其中被改变的第一参数值和第二参数值中的至少一个表示所述初始模型的一个或更多个构造层位位置的改变。

23.权利要求21的系统，其中被改变的第一参数值和第二参数值中的至少一个表示所述初始模型的一个或更多个区域中介质性质的改变。

24.权利要求21的系统，其中所述处理器使用时间保持不变的层析成像来在保持沿每个射线对的总走时不变的情况下计算用于模拟所述初始模型的不同扰动的层析成像系数。

25.权利要求21的系统，其中存储器存储所述改变后的模型。

说明书

技术领域

本发明涉及从地震数据导出的模型参数(如地下地质模型)的改变以便生成新的改变后的模型，例如通过改变模型参数。

背景技术

理想地，地震成像数据集可通过唯一的地下地质模型来表示。然而，通常由地震成像导出的地下地质模型不是唯一的；就是说，地下地质模型表示和模型参数值的许多组合可用来满足同一成像条件。通常，这种变化是由于在获取地震数据方面的限制造成的，如震源和接收器之间的炮检距小以及方位覆盖范围有限。对地下地质模型的非唯一性有贡献的其它因素例如是：有限的频带、噪声、影区以及所用的波传播的射线方法和理论的限制(例如，其中没有考虑某些因素，如各向异性(对方向的依赖性)或各向同性(不依赖于方向)、由衰减造成的频散、散射、高频损耗等)等。这样，导出准确的模型是基于从地震成像、测井和先验地质知识得到的信息所进行的解释过程。这一过程通常是要求很多的过程，通常需要大量计算能力和大量人力介入。
作为结果而导出的地下模型通常包括几百万个数据点，每个数据点包含关于结构的几何形状和介质性质的信息。所以，改变模型通常是很复杂和费力的任务。
需要一种准确和有效的手段来改变地下地质模型同时将计算障碍和对用户输入的需求最小化。

发明内容

本发明的实施例包括使用时间保持不变的层析成像进行地震数据建模的系统和方法，包括存储表征初始地震数据模型的参数值的初始集。初始地震模型可对应于至少多个射线对。每个射线对可具有走时。通过对初始模型中的多个射线对中的每个射线对改变参数值的初始集中的两个或更多个参数值，可产生改变后的模型。改变所述两个或更多个参数值中的一个参数值而不改变其余参数值可对应于每个射线对的走时变化，而组合地改变所述两个或更多个参数值通常对应于每个射线对的走时没有净变化。

附图说明

下述附图中显示了本发明的各个实施例，它们只是作为示例，不对本发明的范围构成限制，其中：
图1是根据本发明一个实施例的地震数据模型的示意图；
图2是根据本发明一个实施例的为形成新的地震数据模型而改变的初始地震数据模型的示意图；
图3是根据本发明一个实施例的被层位垂直位移代替的层位法向位移的示意图；
图4是可根据本发明一个实施例而使用的系统的示意图；以及
图5是根据本发明一个实施例的方法的流程图。
应该理解，为了图示简单和清楚，图中所示元素不一定精确地和成比例地绘制。例如，为了显示清楚，某些元素的尺度相对于其它元素的尺度被放大，或多个物理部件被包括在一个功能块或元素中。再有，在认为适当的地方，附图标记可在附图中重复，以指出相应的或相似的元素。而且，图中描绘的某些块可被组合为单个功能。

具体实施方式

在下文的描述中将描述本发明的各个方面。为了便于解释，将给出具体配置和细节，以便提供对本发明的透彻理解。然而，本领域技术人员显然知道也可以不在本文所述的具体配置和细节的情况下实施本发明。再有，为了使本发明清楚，公知的特点可能被略去或简化。
标准的地下模型可以是地震数据集的初始的、通用的或基本的表示，以用作该地震数据集的每个新的表示的构成块。标准的地下模型可(例如至少是部分地)由运动学特性或参数值的初始集来定义。该参数值的初始集通常包括被模型化的数据的初始估计的、基本的或固有的性质(例如地下位置、空间关系、介质性质等)。该参数值的初始集可设置成例如模板矩阵Q。
通过对标准模型进行改变或扰动来改变作为模板的、基本的或标准的模型，可以定义每个新的地震数据模型。造成这种改变的数据可设置成例如扰动矢量X。改变后的模型(例如在改变后的模型中的走时)可通过模板矩阵Q和扰动矢量X相乘得到。造成改变或扰动的数据或信息可存储在除扰动矢量以外的其它数据结构中。
与地震数据集的新模型是通过对数据的每个改变进行“从无到有(from scratch)”地构建的传统方法相比，本发明的实施例可以再使用先前产生的例如存储在模板矩阵Q中的标准模型数据。该模板矩阵可以是产生地震数据集的所有其它模型的构成块。例如，每个不同的扰动矢量X可应用于唯一的模板矩阵Q，以构建同一数据的模型的新变体。这些不同的模型可进行比较，以选出最佳的建模参数值。
初始模型可例如由射线追踪机制来处理，以至少产生相应的多个射线对。对应于初始模型的每个射线对可具有走时(例如，波沿着射线对路径从震源传播到接收器的时间)。由于地下模型是从地震数据集导出的，故假定沿模拟的反射射线对计算的走时是正确的。这样，各模型参数的改变必须保持这些走时不变。对模型参数集的任何改变都可以伴随有其它参数集的相应改变，以保持沿每个射线对的总走时不变。
可使用时间保持不变的层析成像通过初始模型来构建改变的模型。在时间保持不变的层析成像中，在初始模型和改变后的模型中，沿射线对的总走时通常是相同的。通常，对各个初始模型参数值的扰动，如水平位移或介质性质的改变，会造成沿所生成的改变后的模型的射线对的走时改变。然而，在时间保持不变的系统中，改变后的模型应没有这种走时改变。由此，层析成像机制可产生一个或更多个额外的扰动，造成沿改变后的模型的射线对的走时发生补偿性变化(例如，相消或大小相等符号相反的变化)。由此，对初始模型应用组合扰动可使得改变后的模型中的走时整体上不出现改变(例如，模型间的走时保持不变)。
1.标准地下模型
参考图1，它是根据本发明一个实施例的地震数据模型100的示意图。模型100可以是三维(3D)空间地下模型。可通过处理地震数据(例如未加工数据)来产生模型100。模型100可由对应于模型每个节点112的参数值集来定义。例如，参数值可包括用于TTI介质性质和层位110的位置的值，如轴向压缩波速VA和两个汤姆森(Thomsen)各向异性参数厄普西隆(ε)和德耳塔(δ)。
模型100可通过例如射线追踪机制进行处理。射线追踪机制可产生用于产生射线对的层析成像系数。这样，在模型100中的数据或由它导出的数据(例如层析成像系数)至少可对应于多个射线对。模型100中的数据可对应于射线对以外的其它数据，例如模型的结构布局(structurallayout)(例如大小或一组维度)、计算关系(例如时间保持不变)等。每个射线对可包括在反射点128相遇的入射射线124和反射射线126。每个射线对可从震源点130发出并在接收器点132被收集。每个射线对可具有走时，例如，波沿着射线对路径从震源点130传播到接收器点132的时间。
模型100可包括多个被界面或层位110(例如由反射点128形成的)分开的地层114。界面可由空间表面(例如层位110)确定，而该空间表面又可由节点112(例如，层位110与垂直网格线的交点)确定。可在层位节点112定义网格120，层位节点通常在横向被有规律地隔开，而在垂直方向可能被不规则地隔开。在另一实施例中，节点112可在垂直方向被有规律地隔开，而在横向被不规则地隔开。也可使用节点112的其它配置。
可在网格节点112定义地震数据，例如包括背景介质性质，如TTI(倾斜横向各向同性)的轴向压缩波速或各向同性波速。层位110的每一点可由相对于参考层位的局部法线122的方向来确定。法线122可由倾角和方位角表示。在每个网格节点可指定TTI介质对称轴的方向。另外，在每个网格节点112，还可指定其它参数，如汤姆森参数δ和ε。
对给定的地下模型100，影响系数或模板矩阵Q(可用除Q以外的其它名字)可表示未被扰动的或初始的模型100的性质，例如，具有定义参数值的初始集或特性的系数。可通过针对初始模型100中的震源点130和接收器点132的不同组合追踪射线来计算出参数值的初始集。通常，初始特性集包括模型100的固有特性，初始特性集可用其最基本的、通用的或标准的形式来描述模型100的物理约束。计算这些系数集可以作为模型100构建过程的一部分。模板矩阵Q可被存储为初始模型100的表示，并可在后面用于对初始模型100参数值的不同扰动的快速和准确的模拟/预测。地震数据集的初始模型100表示可用作产生同一地震数据集的各个新的模型变体的基础。
2.修改模型
参考图2，它是根据本发明一个实施例的为形成新的地震数据模型100a而改变的初始地震数据模型100的示意图。在一个实施例中，改变或扰动可包括例如反射层位110的位移。在图2中，反射点128可移动一个反射点位移142(例如，指向新位置反射点128a的矢量Δd)。由于层位110的表面可被认为是局部平坦的，反射点位移142可被分解成法向位移(例如沿指向层位110表面的法线122方向上的位移)和切向位移。通常，只有法向位移影响相应射线对的走时。
对于初始地下模型100的给定改变或扰动，可产生相应的扰动矢量X
(对于扰动矢量，可用其它名字)来表示该改变。扰动矢量X可应用于模板矩阵Q(例如通过右乘)，以产生对应于更新后模型100a的新走时。扰动矢量X的系数可包括模型参数的改变量。
在一个实施例中，可在精细网格(例如图1的网格120)上定义模板矩阵Q的初始模型100的参数值，而在稀疏网格上定义扰动矢量X的模型参数改变量。例如，精细网格的分辨率可为4毫秒，而稀疏网格的分辨率可为50毫秒。由于稀疏网格不如精细网格详细，对各个稀疏网格单元定义的模型参数值的改变通常在稀疏网格的单元内为分段常数，而在精细网格的单元内，值通常是连续变化的。可使用其它的网格细化或配置。例如，“稀疏”和“精细”可以是彼此相对而言的术语。
对于模型100，模板矩阵Q的层析成像系数可通过沿镜面射线对(例如入射射线124和反射射线126)的射线追踪来计算。针对一组张角(opening angle)和基本上所有的方位角，可选择不同的精细网格层位节点作为震源和接收器，来执行射线追踪。层析成像系数可设置在模板矩阵Q中，其中，矩阵的行数(例如线性方程)通常等于被追踪射线的数目。对于较大规模的3D模型，该数目可能太大，以致于得到的模板矩阵Q通常不能被有效地存储在某些盘上或某些计算机存储器中。本发明的实施例可包括将矩阵Q转换成相对较小的压缩矩阵Qij。每个压缩矩阵Qij的维度可等于例如在稀疏网格120中的节点个数。
本发明的实施例提供一种系统和方法，用于快速和准确地模拟/预测地下地质模型100介质性质或构造层位110位置的改变。
在一些实施例中，时间可保持不变。在时间保持不变的层析成像系统中，地下模型100可从地震成像中导出，其中假定沿镜面的入射射线124和反射射线126穿过模型100的走时是正确的(例如，与记录的地震数据完全一致)。此外，对任何扰动的模型100a，这些走时可保持不变。例如，沿射线124和126的走时可近似等于沿射线124a和126a的走时。
每个参数值的单个改变通常会造成沿模型100镜面射线走时的改变。然而，在时间保持不变的层析成像系统中可补偿参数值的改变，使得当参数值的改变被组合起来时，所造成的走时改变被抵消，结果没有沿镜面射线的走时净改变或总改变。
例如，不需要改变模型中所有点。本领域的技术人员可以理解，改变模型100可包括只改变模型100中相应于地震数据集的子集的那一部分(例如网格120的节点112)。在一些实施例中，模型的这一部分可只影响射线对的子集的走时。所以，可通过只改变初始模型的某些参数值的子集来改变模型。在一个实施例中，被改变的参数值可固有地只影响射线对的子集。
时间保持不变的层析成像可用于地球物理解释和生产任务，如不确定性分析、预测由于介质性质的不同表示造成的层位模型变化(例如从各向同性到各向异性，或从某种程度的各向异性到另一种程度的各向异性)。再有，时间保持不变的层析成像可用于确定由于在井中测量的标准层位(horizon marker)和模型层位位置之间发生错误联系或不一致所造成的各向异性参数值。如这里将进一步详细讨论的那样，这种不一致还可用于快速和准确的交互式模型100校正，例如，在地震数据收集或钻井过程中，由此得到最佳的地质导向(geosteering)方案。
3.目标
例如，本发明的实施例可包括：在给定介质性质扰动的情况下，对模型层位110位移的快速和准确的层析成像分析和预测。例如，本发明的实施例可包括：在给定层位位移和其它介质扰动的情况下，各向异性介质参数值的变化的快速和准确的模拟。
本发明的实施例提供层析成像工具，例如基于各向异性模型参数的小变化与其它模型参数值的扰动之间的线性关系的层析成像工具。参数的每个变化通常造成双向路径(例如入射射线124和反射射线126)上的走时出现残差或变化。然而，由于模型100满足成像条件，可以假定总走时得到保持(例如保持不变)。
地下地质模型100可包括由界面(例如地质层位110)分开的一组地层(例如地层114)。地层界面(层位110)可以是反射和/或折射层位。通常假定各向异性介质参数值的分布在每个地层114中是连续和平滑的，而在沿界面(例如层位110)的过渡区是不连续的。在一个实施例中，模型100的介质可以是倾斜横向各向同性(TTI)介质，例如，在每一点或在每个地质地层114具有给定的对称轴取向和给定的层位法线122的方向。也可以使用其它类型或介质，更一般地如横向各向同性(TI)介质。
在一个实施例中，模型100可由3D空间中的(例如两个)独立的方向来定义。一个独立的方向可以是介质(例如TTI)的对称轴的取向，它可由例如顶点(zenith)和方位角来描述。另一个独立方向可以是层位表面的法线，它可被定义在每个层位110的任何点(例如，但并非必需，在3D空间的任何点)。层位表面法线122的取向还可由例如顶点和方位角来描述。这样，在与层位110相交的每个节点112，法线112的取向110通常可不同于介质对称轴的取向。在3D空间中不对应于节点112的点处，通常只存在对称轴取向。
可变的模型100的参数可包括：例如界面或层位110的位置和TTI介质的性质，如轴向压缩波速VA，两个汤姆森各向异性参数ε和δ(例如，总共四个参数)。被扰动的模型参数可在层位节点112定义，例如按稀疏规则横向网格进行定义。在模型网格120的精细网格节点和稀疏网格节点之间可有3D空间内插。使用内插函数，如果给定其它三个参数变化，则可以在任何模型节点112处模拟任何参数的变化。这里将进一步详细解释这项技术。
4.主要假设
在一些实施例中，时间保持不变的层析成像可用于改变标准地下模型100。时间保持不变的层析成像可基于下列假设：
·通常保持沿镜面的入射射线124和反射射线126的走时不变。
·通常模型参数的扰动较小，例如在不同参数的扰动之间具有线性关系。
·通常每次针对一个参数集来进行预测/模拟(例如，假定其它参数是给定的)。
·通常入射射线124和反射射线126的轨迹是固定的。这样，介质性质的小扰动可造成走时变化，该走时变化是由于沿射线路径的射线速度的变化但通常不是由于射线轨迹的几何变化造成的。
其它实施例可使用不同的或其它的假设。虽然这里描述的实施例使用时间保持不变的层析成像，但本领域技术人员会理解这并不是必须的。本发明的其它实施例可使用这些假设中的一个或更多个假设不成立的其它层析成像方法。例如，当模型参数被扰动或改变，沿镜面的入射射线124和反射射线126的走时也改变。
5.在横向各向同性介质中的射线追踪
可沿反射表面(例如层面110)从一组节点112执行TI介质(例如3D TTI介质)中的射线追踪。例如，从每个反射点128，以针对基本上全部方位角的不同的张角来追踪一组镜面射线对(例如入射射线124和反射射线126)。沿入射射线124和反射射线126的轨迹，可计算走时、弧长和层析成像系数。数学上，射线追踪通常是对一组有初始条件(例如起始点位置、反射点128以及相速度的方向)的常微分方程(ODE)的数值解。
为简化射线追踪过程，在数值求解过程的每一步可转动参考系，从而使介质的对称轴变为垂直。这样，射线追踪问题可简化为追踪具有垂直对称轴的TI介质(VTI)。例如，可在VTI参考系内计算慢度分量px、py、pz和射线路径的增量Δx、Δy、Δz，然后返回全局参考系。可进行转动以实现各个新的VTI参考系(例如，当轴的取向不同于前一个反射点128的轴取向时)。射线追踪算法可使用介质性质的导数，这些导数通常是在VTI参考系中计算的。沿着射线轨迹，射线哈密顿算子G通常会消失。VTI哈密顿算子G可例如由下式定义：
$G(p_h,p_z,x,y,z)=\frac{K-{\mathrm{LV}}_P^2-V_P^{-2}}{2(1-f)}---\left(5.1\right)$
其中函数K和L可例如由下式定义：
$K=(1+f)·(p_h^2+p_z^2)+2ϵp_h^2$(5.2)
$L=f·{(p_h^2+p_z^2)}^2+2ϵp_h^2·({\mathrm{fp}}_h^2+p_z^2)-2\delta (1-f)p_h^2{p}_z^2$
如同这里描述的其它公式一样，也可使用其它公式。在国际单位制(SI)中，哈密顿算子G的单位可以是[G]＝s2/m2。参数pz和ph可以分别是射线慢度的垂直(例如轴向)和水平(例如垂直于轴向)分量，例如：
$p_h^2=p_x^2+p_y^2---\left(5.3\right)$
参数VP可以是轴向压缩波速，δ和ε可以分别是第一和第二汤姆森各向异性参数δ和ε。最后，参数f可以是轴向剪切波速VS的平方与轴向压缩波速VP的平方之比。例如，对于压缩波(例如，或准压缩波)，可假定下列常数参数：
$f=V_S^2/V_P^2\approx 1/4---\left(5.4\right)$
TI介质可由(例如五个)参数描述。然而，其中至少一个参数，即汤姆森γ(汤姆森gamma)，通常是与其它参数不耦合的。例如，汤姆森γ可描述在水平面中的纯剪切运动。参数f通常被假定为常数且已知。这样，未知参数通常例如为轴向压缩波速、汤姆森δ和ε。
6.更新慢度量值
由于哈密顿算子G通常沿射线(例如入射射线124和反射射线126)消失，所以例如在给定慢度方向的情况下，可在射线轨迹的每一点处更新慢度的量值p。假定θphs为慢度矢量的倾角(例如，或顶角(zenith angle))(例如，方位角通常对于VTI无关紧要)，慢度分量例如为：
ph＝p·sin α，  pz＝p·cos    (6.1)
求解关于哈密顿算子G的方程，则归一化相速度例如可以是：
$\frac{V_{\mathrm{phs}}^2}{V_P^2}=\frac{1+f}{2}+ϵ{\sin }^2{\theta }_{\mathrm{phs}}\pm \frac{1}{2}\sqrt{{(1-f+2ϵ{\sin }^2{\theta }_{\mathrm{phs}})}^2-2(ϵ-\delta )(1-f){\sin }^2{2\theta }_{\mathrm{phs}}}---\left(6.2\right)$
其中，当波在垂直平面极化时，如果研究的是准压缩波，则最右边一项通常为正值，而当研究的是准剪切波，它通常为负值。相速度的倒数是慢度的量值，而相速度的方向通常与慢度方向一致，例如有下列关系：
$p=V_{\mathrm{phs}}^{-1}\to p_h=\frac{\sin {\theta }_{\mathrm{phs}}}{V_{\mathrm{phs}}},$ $p_z=\frac{\cos {\theta }_{\mathrm{phs}}}{V_{\mathrm{phs}}}---\left(6.3\right)$
7.射线追踪方程组
通常可利用一组射线追踪方程来构建射线对轨迹。该方程组通常包括(例如6个)基本ODE，还包括一个或更多个附加补充ODE。基本ODE
方程可描述射线点位置的导数，例如下式：

和慢度分量的导数，例如下式：
$\frac{{\mathrm{dp}}_x}{\mathrm{d\sigma }}=-\frac{\partial G}{\partial x},$ $\frac{{\mathrm{dp}}_y}{\mathrm{d\sigma }}=-\frac{\partial G}{\partial y},$ $\frac{{\mathrm{dp}}_z}{\mathrm{d\sigma }}=-\frac{\partial G}{\partial z}---\left(7.2\right)$
其中σ可以是沿射线路径的独立参数，其SI单位是例如[σ]＝m2/s(例如，速度的平方×时间)。这样，dr/dσ的SI单位可以是[dr/dσ]＝s/m，其中r表示x、y、z，dp/dσ的SI单位可以是[dp/dσ]＝s2/m3。哈密顿算子G的笛卡尔偏导数通常与3D模型100体中的介质性质的变化有关，例如，按照下列关系：
$\frac{\partial G}{\partial x}=\frac{\partial G}{\partial V_P}·\frac{{\partial V}_P}{\partial x}+\frac{\partial G}{\partial \delta }·\frac{\partial \delta }{\partial x}+\frac{\partial G}{\partial ϵ}·\frac{\partial ϵ}{\partial x}$
$\frac{\partial G}{\partial y}=\frac{\partial G}{\partial V_P}·\frac{\partial V_P}{\partial y}+\frac{\partial G}{\partial \delta }·\frac{\partial \delta }{\partial y}+\frac{\partial G}{\partial ϵ}·\frac{\partial ϵ}{\partial y}---\left(7.3\right)$
$\frac{\partial G}{\partial z}=\frac{\partial G}{\partial V_P}·\frac{\partial V_P}{\partial z}+\frac{\partial G}{\partial \delta }·\frac{\partial \delta }{\partial z}+\frac{\partial G}{\partial ϵ}·\frac{\partial ϵ}{\partial z}$
哈密顿算子G对于介质性质的导数例如可以是：
$\frac{\partial G}{\partial V_P}=\frac{2}{(1-f)·V_P^3}-\frac{2ϵp_h^2}{(1-f)·V_P}-\frac{1+f}{1-f}·\frac{p_h^2+p_z^2}{V_P}$(7.4)
$\frac{\partial G}{\partial ϵ}=\frac{1-({\mathrm{fp}}_h^2+p_z^2)·V_P^2}{1-f}·p_h^2,$ $\frac{\partial G}{\partial \delta }=p_h^2{p}_z^2·V_P^2$
ODE附加方程组可包括：例如，弧长和/或走时(例如在题为“弧长和走时”的部分中进一步详细描述的)和/或层析成像系数(例如在题为“层析成像系数”部分中进一步详细描述的)。
8.弧长和走时
沿射线的弧长增量dL可例如由下式定义：
$\mathrm{dL}=\sqrt{{\mathrm{dx}}^2+{\mathrm{dy}}^2+{\mathrm{dz}}^2}---\left(8.1\right)$
这样，弧长导数dL/dσ可例如由下式定义：
$\frac{\mathrm{dL}}{\mathrm{d\sigma }}=\sqrt{{\left(\frac{\mathrm{dx}}{\mathrm{d\sigma }}\right)}^2+{\left(\frac{\mathrm{dy}}{\mathrm{d\sigma }}\right)}^2+{\left(\frac{\mathrm{dz}}{\mathrm{d\sigma }}\right)}^2}---\left(8.2\right)$
结合方程(8.2)和射线追踪方程(7.1)，例如得出：
$\frac{\mathrm{dL}}{\mathrm{d\sigma }}=\sqrt{{\left(\frac{\partial G}{\partial p_x}\right)}^2+{\left(\frac{\partial G}{\partial p_z}\right)}^2+{\left(\frac{\partial G}{\partial p_z}\right)}^2}=\sqrt{{\left(\frac{\partial G}{\partial p_h}\right)}^2+{\left(\frac{\partial G}{\partial p_z}\right)}^2}---\left(8.3\right)$
沿射线的走时导数可由链规则得到：
$\frac{\mathrm{dt}}{\mathrm{d\sigma }}=\frac{\partial t}{\partial x}·\frac{\mathrm{dx}}{\mathrm{d\sigma }}+\frac{\partial t}{\partial y}·\frac{\mathrm{dy}}{\mathrm{d\sigma }}+\frac{\partial t}{\partial z}·\frac{\mathrm{dz}}{\mathrm{d\sigma }}---\left(8.4\right)$
结合式(8.4)和式(7.1)，例如得出：
$\frac{\mathrm{dt}}{\mathrm{d\sigma }}=\frac{\partial G}{\partial p_x}·p_x+\frac{\partial G}{\partial p_y}·p_y+\frac{\partial G}{\partial p_z}·p_z---\left(8.5\right)$
弧长导数dL/dσ的SI单位是[dL/d σ]＝s/m，dt/dσ的SI单位是[dt/dσ]＝s2/m2。在一些实施例中，射线路径的弧长和走时可被“在线(on-the-fly)”积分，即在求解射线追踪方程组(例如基本ODE)期间积分。
9.走时的变化
TI介质性质(例如轴向压缩波速、ε、δ值)的小变化与生成的沿整个射线或双向射线对(例如入射射线124和反射射线126)的走时变化有关。走时变化遵循走时方程(8.5)，例如如下：
$\frac{\Delta \left(\mathrm{dt}\right)}{\mathrm{d\sigma }}=\frac{\partial G}{\partial p_x}·\Delta p_{x^{\prime }}+\frac{\partial G}{\partial p_y}·\Delta p_y+\frac{\partial G}{\partial p_z}·\Delta p_z$(9.1)
$+\Delta \left(\frac{\partial G}{\partial p_x}\right)·p_x+\Delta \left(\frac{\partial G}{\partial p_y}\right)·p_y+\Delta \left(\frac{\partial G}{\partial p_z}\right)·p_z$
注意，当参数σ是独立参数时，可给出任何任意的从属函数F(σ)，如走时、弧长等，例如下式：
$\Delta \left(\frac{\mathrm{dF}}{\mathrm{d\sigma }}\right)=\frac{\Delta \left(\mathrm{dF}\right)}{\mathrm{d\sigma }}---\left(9.2\right)$
用式(7.1)中的对应项代替式(9.1)中的右侧项或“尾端”项，例如给出：
$\frac{\Delta \left(\mathrm{dt}\right)}{\mathrm{d\sigma }}=\frac{\partial G}{\partial p_x}·\Delta p_x+\frac{\partial G}{\partial p_y}·\Delta p_y+\frac{\partial G}{\partial p_z}·\Delta p_z$(9.3)
$+\Delta \left(\frac{\mathrm{dx}}{\mathrm{d\sigma }}\right)·p_x+\Delta \left(\frac{\mathrm{dy}}{\mathrm{d\sigma }}\right)·p_y+\Delta \left(\frac{\mathrm{dz}}{\mathrm{d\sigma }}\right)·p_z$
当参数σ独立于射线路径位置x、y、z时，根据式(9.2)，得到例如：
$\frac{\Delta \left(\mathrm{dt}\right)}{\mathrm{d\sigma }}=\frac{\partial G}{\partial p_x}·\Delta p_x+\frac{\partial G}{\partial p_y}·\Delta p_y+\frac{\partial G}{\partial p_z}·\Delta p_z+\frac{\Delta \left(\mathrm{dx}\right)·p_x+\Delta \left(\mathrm{dy}\right)·p_y+\Delta \left(\mathrm{dz}\right)·p_z}{\mathrm{d\sigma }}---\left(9.4\right)$
当射线轨迹固定时，介质性质的小扰动通常不会造成可以影响走时的射线路径变化，所以射线路径变化可以被忽略，例如给出：
Δ(dx)＝Δ(dy)＝Δ(dz)＝0    (9.5)
结合式(9.4)和式(9.5)，例如给出：
$\frac{\Delta \left(\mathrm{dt}\right)}{\mathrm{d\sigma }}=\frac{\partial G}{\partial p_x}·\Delta p_x+\frac{\partial G}{\partial p_y}·\Delta p_y+\frac{\partial G}{\partial p_z}·\Delta p_z---\left(9.6\right)$
如这里描述的那样，哈密顿算子通常沿射线路径消失。当哈密顿算子为常数时，它的全变化(full variation)ΔG消失，例如下式：
G＝const    ΔG＝0    (9.7)
式(9.7)中的const表示常数。射线哈密顿算子G通常依赖于慢度分量和介质参数，因此例如有如下关系：
$\mathrm{\Delta G}=\frac{\partial G}{\partial p_x}·\Delta p_x+\frac{\partial G}{\partial p_y}·\Delta p_y+\frac{\partial G}{\partial p_z}·\Delta p_z+\frac{\partial G}{\partial V}·\mathrm{\Delta V}+\frac{\partial G}{\partial \delta }·\mathrm{\Delta \delta }+\frac{\partial G}{\partial ϵ}·\mathrm{\Delta ϵ}=0---\left(9.8\right)$
哈密顿变化之和ΔG的最左边三项可与走时变化有关。结合式(9.8)和式(9.6)，例如给出：
$\frac{\Delta \left({\mathrm{dt}}_{\mathrm{medium}}\right)}{\mathrm{d\sigma }}=\frac{d\left(\Delta t_{\mathrm{medium}}\right)}{\mathrm{d\sigma }}=\frac{\partial G}{\partial V_P}·\Delta V_P-\frac{\partial G}{\partial \delta }·\mathrm{\Delta \delta }-\frac{\partial G}{\partial ϵ}·\mathrm{\Delta ϵ}---\left(9.9\right)$
在等式(9.9)中的下标“medium(介质)”可表示：在本式中走时变化的导数通常是由于介质性质的小变化造成的，通常不考虑层位位移。
10.层析成像系数
层析成像系数可将层位位移及介质性质的小扰动与沿镜面射线(例如入射射线124和反射射线126)的走时误差联系起来。为得到由介质性质变化造成的走时残差，可沿射线路径对等式(9.9)进行积分，例如如下式所示：
$\Delta t_{\mathrm{medium}}=-\Delta V_P·\underset{\sigma }{\int }\frac{\partial G}{\partial V_P}\mathrm{d\sigma }-\mathrm{\Delta \delta }·\underset{\sigma }{\int }\frac{\partial G}{\partial \delta }\mathrm{d\sigma }-\mathrm{\Delta ϵ}·\underset{\sigma }{\int }\frac{\partial G}{\partial ϵ}\mathrm{d\sigma }---\left(10.1\right)$
层析成像系数AV、Aδ、Aε可以是走时对模型参数变化的导数，例如：
$A_V=\frac{\partial \left(\Delta t_{\mathrm{medium}}\right)}{\partial \Delta V_P},$ $A_{\delta }=\frac{\partial \left(\Delta t_{\mathrm{medium}}\right)}{\partial \mathrm{\Delta \delta }},$ $A_{ϵ}=\frac{\partial \left(\Delta t_{\mathrm{medium}}\right)}{\partial \mathrm{\Delta ϵ}}---\left(10.2\right)$
于是，例如：
$A_V=-\underset{\sigma }{\int }\frac{\partial G}{\partial V_P}\mathrm{d\sigma },$ $A_{\delta }=-\underset{\sigma }{\int }\frac{\partial G}{\partial \delta }\mathrm{d\sigma },$ $A_{ϵ}=-\underset{\sigma }{\int }\frac{\partial G}{\partial ϵ}\mathrm{d\sigma }---\left(10.3\right)$
其中可在各向同性介质中或各向异性介质中沿射线对轨迹进行积分。结合式(10.3)和(10.1)，例如给出：
Δtmedium＝AV·ΔVP+Aδ·Δδ+Aε·Δε     (10.4)
等式(10.4)对射线对的给定的具体轨迹成立。注意，等式(10.3)可排列成ODE方程组的形式，例如：
$\frac{{\mathrm{dA}}_V}{\mathrm{d\sigma }}=-\frac{\partial G}{\partial V_P},$ $\frac{d{A}_{\delta }}{\mathrm{d\sigma }}=-\frac{\partial G}{\partial \delta },$ $\frac{{\mathrm{dA}}_{ϵ}}{\mathrm{d\sigma }}=-\frac{\partial G}{\partial ϵ}---\left(10.5\right)$
方程组(10.5)表示定义模型100的射线追踪方程组的(例如三个)附加ODE。附ODE可在基本射线追踪方程组和/或计算走时和弧长的计算期间被积分。这样，可在例如射线追踪期间“在线(on-the-fly)”计算层析成像系数。对于模型参数值的每次改变，这些系数可用于预测其它模型参数值的相应改变，以保持沿模型100的射线对的总走时不变。然而，层析成像系数本身通常不因模型参数值的变化而变化。在题为“水平位移”的部分中将对考虑了水平位移变化的层析成像系数进一步详细描述。
11.层位位移
射线的折射点、反射点和层位位置中的每一个的变化通常会影响射线的走时。
再参考图2，其中层位从第一位置的层位110位移到第二位置的层位110a。模型100可包括由层位表面110分开的地质地层114(例如，如参考图1描述的那样)。图2中示出了双向射线路径，包括具有初始反射点128的初始反射射线126和初始入射射线124。例如，由于层位位移或其它扰动，初始反射点128可移动到改变的反射点128a，造成新的改变后的双向射线路径，该路径具有改变后的入射射线124a和改变后的反射射线126a。
首先讨论由反射点128位移造成的走时变化，然后讨论在该变化之后由折射点(例如震源点130和接收器点132)位移造成的走时变化。这里可将走时变化称作“走时残差”。数学上，首先对二维(2D)模型考虑层位位移，然后将针对2D模型得到的结果外推到一般3D模型。射线对走时Told可例如由下式定义：
$T_{\mathrm{old}}^2=\frac{{(x_S-x_I^{\mathrm{old}})}^2+{(z_S-z_I^{\mathrm{old}})}^2}{V_{\mathrm{in},\mathrm{Ray}}^2}+\frac{{(x_R-x_I^{\mathrm{old}})}^2+{(z_R-z_I^{\mathrm{old}})}^2}{V_{\mathrm{re},\mathrm{Ray}}^2}---\left(11.1\right)$
等式中的上标或下标“old”表示与旧的反射点有关。其中Vin，Ray和Vre，Ray可分别是入射射线124和在镜像或反射点128的反射射线126的群速度。震源点130和接收器点132可以分别是入射射线124和反射射线126的位于最近的上层位(例如折射层位)134的折射点。通常，只需要在反射点128的附近考虑反射点128的速度。
例如，在一个实施例中，为简化数学推导，可假定所考虑的地层中的速度是相同的，但这并不是必需的。这部分的结果对于任何非均一速度场都是有效的。例如，在其它实施例中，在地层中的速度可以是非均一的或分段均一的。
从初始反射点128到改变后的反射点128a的反射或镜像点的位移可包括两个对应矢量：矢量140，它可以是层位110的位移，以及矢量142，它可以是反射点沿层位110a的位移。可认为层位110a是局部平坦的。对于位移后的射线对124a和126a，走时Tnew可例如由下式定义：
$T_{\mathrm{new}}^2=\frac{{(x_S-x_I^{\mathrm{new}})}^2+{(z_S-z_I^{\mathrm{new}})}^2}{V_{\mathrm{in},\mathrm{Ray}}^2}+\frac{{(x_R-x_I^{\mathrm{new}})}^2+{(z_R-z_I^{\mathrm{new}})}^2}{V_{\mathrm{re},\mathrm{Ray}}^2}---\left(11.2\right)$
等式中的上标或下标“new”表示与新的反射点有关。其中“新的”改变后的反射点128a的位置可例如表示为：
$x_I^{\mathrm{new}}=x_I^{\mathrm{old}}+\Delta d_x+\mathrm{\Delta s}·\cos \alpha ,$ $z_I^{\mathrm{new}}=z_I^{\mathrm{old}}+\Delta d_z-\mathrm{\Delta s}·\sin \alpha ---\left(11.3\right)$
Δs可以是新反射点128a沿层位表面110a的位移矢量Δs142的带符号长度。注意，图2中的倾角α146是负的。图2中的带符号位移Δs也是负的。对于等于零的倾角146，正的位移方向通常与x轴148的正方向一致。折射点(例如震源点130和点132)可被认为是显著地远离初始的反射点128和改变后的反射点的128a。这样，旧的和新的入射射线124和124a的方向可被认为是近乎相同。类似地，两个反射射线126和126a的方向也可被认为是近乎相同。速度非均一性对走时变化的影响也可被认为是可以忽略的(例如，数学上是较高阶无限小)。走时残差通常主要是由射线路径长度的变化造成的。结合式(11.3)和式(11.2)，例如给出：
$T_{\mathrm{new}}^2=\frac{{(x_S-x_I^{\mathrm{old}}-\Delta d_x-\mathrm{\Delta s}·\cos \alpha )}^2+{(z_S-z_I^{\mathrm{old}}-\Delta d_z+\mathrm{\Delta s}·\sin \alpha )}^2}{V_{\mathrm{in},\mathrm{Ray}}^2}$(11.4)
$+\frac{{(x_R-x_I^{\mathrm{old}}-\Delta d_x-\mathrm{\Delta s}·\cos \alpha )}^2+{(z_R-z_I^{\mathrm{old}}-\Delta d_z+\mathrm{\Delta s}·\sin \alpha )}^2}{V_{\mathrm{re},\mathrm{Ray}}^2}$
沿层位的位移Δs通常是未知的，可利用斯奈尔(Snell)定律来计算，它使新的走时Tnew为最小，例如下式：
${\mathrm{dT}}_{\mathrm{new}}^2/\mathrm{d\Delta s}=0---\left(11.5\right)$
这导致
$\frac{\cos \alpha ·(x_S-x_I^{\mathrm{old}}-\Delta d_x-\mathrm{\Delta s}·\cos \alpha )-\sin \alpha ·(z_S-z_I^{\mathrm{old}}-\Delta d_z+\mathrm{\Delta s}·\sin \alpha )}{V_{\mathrm{in},\mathrm{Ray}}^2}+$(11.6)
$\frac{\cos \alpha ·(x_R-x_I^{\mathrm{old}}-\Delta d_x-\mathrm{\Delta s}·\cos \alpha )-\sin \alpha ·(z_R-z_I^{\mathrm{old}}-\Delta d_z+\mathrm{\Delta s}·\sin \alpha )}{V_{\mathrm{re},\mathrm{Ray}}^2}=0$
对于包括入射射线124和反射射线126的初始“旧”射线对或者背景射线对，斯奈尔定律也是成立的。在初始反射点128附近，点的位置可例如由下式给出：
$x_I^{\mathrm{old}}=a+s·\cos \alpha ,$ $z_I^{\mathrm{old}}=b-s·\sin \alpha ---\left(11.7\right)$
其中s是沿初始背景层位110(例如，也被描述为局部平表面)的标量参数。通过结合式(11.1)和(11.7)，可例如通过下式来定义初始射线对的走时Told：
$T_{\mathrm{old}}^2=\frac{{(x_S-a-s·\cos \alpha )}^2+{(z_S-b+s·\sin \alpha )}^2}{V_{\mathrm{in},\mathrm{Ray}}^2}$(11.8)
$+\frac{{(x_R-a-s·\cos \alpha )}^2+{(z_R-b+s·\sin \alpha )}^2}{V_{\mathrm{re},\mathrm{Ray}}^2}$
使走时最小则给出参数s，例如如下式所示：
${\mathrm{dT}}_{\mathrm{old}}^2/\mathrm{ds}=0---\left(11.9\right)$
这导致例如：
$\frac{\cos \alpha ·(x_S-a-s·\cos \alpha )-\sin \alpha ·(z_S-b+s·\sin \alpha )}{V_{\mathrm{in},\mathrm{Ray}}^2}+$(11.10)
$\frac{\cos \alpha ·(x_R-a-s·\cos \alpha )-\sin \alpha ·(z_R-b+s·\sin \alpha )}{V_{\mathrm{re},\mathrm{Ray}}^2}=0$
结合式(11.7)和式(11.10)，例如给出：
$\frac{\cos \alpha ·(x_S-x_I^{\mathrm{old}})-\sin \alpha ·(z_S-z_I^{\mathrm{old}})}{V_{\mathrm{in},\mathrm{Ray}}^2}+\frac{\cos \alpha ·(x_R-x_I^{\mathrm{old}})-\sin \alpha ·(z_R-z_I^{\mathrm{old}})}{V_{\mathrm{re},\mathrm{Ray}}^2}=0---\left(11.11\right)$
从等式(11.6)(例如在旧反射或镜像点128a处的斯奈尔定律)中减去等式(11.11)(例如在新反射或镜像点128b处的斯奈尔定律)，例如给出：
$\frac{\cos \alpha ·(\Delta d_x+\mathrm{\Delta s}·\cos \alpha )-\sin \alpha ·(\Delta d_z-\mathrm{\Delta s}·\sin \alpha )}{V_{\mathrm{in},\mathrm{Ray}}^2}+$(11.12)
$\frac{\cos \alpha ·(\Delta d_x+\mathrm{\Delta s}·\cos \alpha )-\sin \alpha ·(\Delta d_z-\mathrm{\Delta s}·\sin \alpha )}{V_{\mathrm{re},\mathrm{Ray}}^2}=0$
等式(11.12)可简化为例如下式：
$(\Delta d_x·\cos \alpha -\Delta d_z\sin \alpha +\mathrm{\Delta s})·(V_{\mathrm{in},\mathrm{Ray}}^{-2}+V_{\mathrm{re},\mathrm{Ray}}^{-2})=0---\left(11.13\right)$
未知参数s可被定义为例如：
Δs＝-Δdx·cosα+Δdzsinα         (11.14)
在向内法线122的分量可例如为：
nx＝sinα，nz＝cosα                (11.15)
反射点128的平面内位移Δs 142的分量可例如为：
Δsx＝Δscosα，Δsz＝-Δssinα     (11.16)
结合式(11.14)和(11.16)，例如给出：
Δsx＝-Δdx·cos2α+Δdz sinαcosα，Δsy＝+Δdx·sinαcosα-Δdzsin2α(11.17)
在等式(11.17)中定义的面内位移Δs 142可表示为矢量形式，例如：
$\overrightarrow{\mathrm{\Delta s}}=-\overrightarrow{n}\times \overrightarrow{\mathrm{\Delta d}}\times \overrightarrow{n}---\left(11.18\right)$
在镜像或反射点128的旧位置与反射或镜像点128a的新位置之间的总位移可例如为：
$\overrightarrow{\mathrm{\Delta I}}\equiv \overrightarrow{\mathrm{\Delta d}}+\overrightarrow{\mathrm{\Delta s}}=\overrightarrow{\mathrm{\Delta d}}-\overrightarrow{n}\times \overrightarrow{\mathrm{\Delta d}}\times \overrightarrow{n}---\left(11.19\right)$
其中双叉乘例如由下式给出：
$\overrightarrow{n}\times \overrightarrow{\mathrm{\Delta d}}\times \overrightarrow{n}\equiv \overrightarrow{n}\times (\overrightarrow{\mathrm{\Delta d}}\times \overrightarrow{n})=(\overrightarrow{n}\times \overrightarrow{\mathrm{\Delta d}})\times \overrightarrow{n}---\left(11.20\right)$
如果假定该射线对的两条射线，即入射射线124和反射射线126，到达镜像或反射点128(例如，分别从震源和接收器到达)，则由反射点位移造成的走时残差可例如为：
$\mathrm{\Delta t}=(\stackrel{\to \mathrm{in}}{p}+\stackrel{\to \mathrm{re}}{p})·\overrightarrow{\mathrm{\Delta I}}\equiv \overrightarrow{\mathrm{\Delta d}}+\overline{\mathrm{\Delta s}}=(\stackrel{\to \mathrm{in}}{p}+\stackrel{\to \mathrm{re}}{p})·(\overrightarrow{\mathrm{\Delta d}}-\overrightarrow{n}\times \overrightarrow{\mathrm{\Delta d}}\times \overrightarrow{n})---\left(11.21\right)$
对于3D模型，可使用局部参考系，其中两个矢量(例如反射面层位110的法线n122和层位位移Δd140)中的每一个，都在xz平面内(例如，由x轴148和z轴150形成的)。这意味着对上述2D模型导出的等式(11.19)
也对一般3D模型成立。引入符号来表示在镜像或反射点128的慢度和，例如：
$\overrightarrow{p}\equiv \stackrel{\to \mathrm{in}}{p}+\stackrel{\to \mathrm{re}}{p}---\left(11.22\right)$
等式(11.21)可例如表示为：
$\mathrm{\Delta t}=\overrightarrow{p}·\overrightarrow{\mathrm{\Delta d}}-\overrightarrow{p}·\overrightarrow{n}\times \overrightarrow{\mathrm{\Delta d}}\times \overrightarrow{n}=\overrightarrow{p}·\overrightarrow{\mathrm{\Delta d}}-\overrightarrow{p}·[\overrightarrow{n}\times (\overrightarrow{\mathrm{\Delta d}}\times \overrightarrow{n})]---\left(11.23\right)$
等式(11.23)的右侧“最右”项是混合乘积，可被重新表示为例如：
$\overrightarrow{p}·[\overrightarrow{n}\times (\overrightarrow{\mathrm{\Delta d}}\times \overrightarrow{n})]=(\overrightarrow{p}\times \overrightarrow{n})·(\overrightarrow{\mathrm{\Delta d}}\times \overrightarrow{n})---\left(11.24\right)$
例如对于任意各向异性介质的斯奈尔反射定律，可为：
$(\stackrel{\to \mathrm{in}}{p}+\stackrel{\to \mathrm{re}}{p})\times \overrightarrow{n}=\overrightarrow{p}\times \overrightarrow{n}=0---\left(11.25\right)$
由斯奈尔定律得到：镜像或反射点128沿反射层位的小位移通常不影响射线对(例如，包括入射射线124和反射射线126)的总走时。这样，对于走时计算，等式(11.18)可以不需要。等式(11.23)可简化为例如：
$\mathrm{\Delta t}=\overrightarrow{p}·\overrightarrow{\mathrm{\Delta d}}---\left(11.26\right)$
再有，层位位移矢量140可分解为两个分量，例如在层位110表面的法线122和在这一表面的切线。然而，切线位移通常没有层析成像意义，所以可认为是可忽略的。这样，可假定层位位移140只有法线分量，例如如下式所示：
$\overrightarrow{\mathrm{\Delta d}}=\mathrm{\Delta d}\overrightarrow{n}---\left(11.27\right)$
其中Δd是标量值。结合等式(11.26)和(11.27)可例如给出：
$\partial t/\partial d=\overrightarrow{p}·\overrightarrow{n}=(\stackrel{\to \mathrm{in}}{p}+\stackrel{\to \mathrm{re}}{p})·\overrightarrow{n}---\left(11.28\right)$
层位位移140的层析成像系数可以是法线122方向和射线对慢度和的标量积。对于向内法线122，两射线(例如入射射线124和反射射线126)可分别从震源点130和接收器点132出发，到达镜像或反射点128a。或者，这两条射线可从反射或镜像点出发，在这种情况下，法线122可以是“向外”指向的(例如，沿着与图2中所示方向相反的方向)。
考虑一个特定情况，其中反射点位移是垂直的。可能已注意到，平面(例如局部平面)层位110的任何位移可表示为垂直位移，例如：
Δdx＝0，Δdy＝0，Δdz＝Δz    (11.29)
参考图3，图中示意性示出层位法线方向位移Δd 140被层位垂直位移Δz 136代替。法线位移140和垂直位移136之间的关系可例如由下式给出：
Δd＝Δzcosα＝Δznz          (11.30)
结合等式(11.30)和等式(11.26)，例如给出：
$\mathrm{\Delta t}=\overrightarrow{p}·\overrightarrow{\mathrm{\Delta d}}=\overrightarrow{p}·{\overrightarrow{\mathrm{\Delta d}}}_n=p_n·\Delta d_n=\overrightarrow{p}·\overrightarrow{n}\Delta {\mathrm{zn}}_z---\left(11.31\right)$
展开等式(11.31)的标量积，例如给出：
$\mathrm{\Delta t}=(p_x{n}_x+p_y{n}_y+p_z{n}_z)\Delta {\mathrm{zn}}_z=(p_x{n}_x{n}_z+p_y{n}_y{n}_z+p_z{n}_z^2)\mathrm{\Delta z}---\left(11.32\right)$
如等式(11.25)中所示的斯奈尔定律也可被展开为例如：
$\overrightarrow{p}\times \overrightarrow{n}=0\to \begin{array}{c}{p}_y{n}_z=p_z{n}_y\\ p_z{n}_x=p_x{n}_z\\ p_x{n}_y=p_y{n}_x\end{array}---\left(11.33\right)$
由式(11.33)得到例如：
$p_x{n}_x{n}_z=p_z{n}_x^2,$ $p_y{n}_y{n}_z=p_z{n}_y^2---\left(11.34\right)$
结合等式(11.34)和(11.32)例如给出：
$\mathrm{\Delta t}=(p_z{n}_x^2+p_z{n}_y^2+p_z{n}_z^2)\mathrm{\Delta z}=p_z(n_x^2+n_y^2+n_z^2)\mathrm{\Delta z}=p_z\mathrm{\Delta z}---\left(11.35\right)$
这样，用于反射点128的垂直位移136的层析成像系数可例如为：
$\Delta t_{\mathrm{reflection}}/\mathrm{\Delta z}=p_z^{\mathrm{reflection}}=p_z^{\mathrm{in}}+p_z^{\mathrm{re}}---\left(11.36\right)$
等式中的“reflection”表示与反射有关。对于折射点(例如，震源点130和接收器点132)，射线对(例如，入射射线124和反射射线126)通常到达层位110的上侧，且从下侧发出。如对反射点128所示，到达射线的慢度应被认为是正的，因此，例如：
$\Delta t_{\mathrm{refraction}}/\mathrm{\Delta z}=p_z^{\mathrm{refraction}}=p_z^{\mathrm{upper}}-p_z^{\mathrm{lower}}---\left(11.37\right)$
等式中的“refraction”表示与折射有关，“upper”表示与上部有关而“lower”表示与下部有关。在等式(11.37)中的正号和/或负号的使用可例如由随后的例子验证。考虑“向下”传播的垂直或几乎垂直的射线(例如从折射层位134的震源点130和接收器点132到反射点128)。考虑在界面之上具有“慢”介质而在界面之下具有“快”介质的折射层位134。术语“慢”介质或“快”介质例如可指的是：射线或波传播穿过该介质的速度。在折射层位134中的正垂直位移(例如向下位移)可增大射线传播通过的慢区域，从而增大射线走时。在这种情况下，当走时trefraction增大时，它在等式(11.37)左侧的导数Δtrefraction/Δz可为正值。现在参考等式(11.37)的右侧，由于在本例中上部介质是“慢”介质而下部介质是“快”介质，对于向下传播的垂直射线，$p_z^{\mathrm{upper}}>p_z^{\mathrm{lower}}.$这样，等式(11.37)的右侧也可为正值。
12.管理方程
层析成像系数(例如AV、Aδ、Aε、As)可以是例如射线124和126的走时对介质性质、层位位移和/或其它改变的性质或模型扰动的导数。例如：
$A_{h,i,j,k}^V=\frac{\partial t_{h,i,j}}{\partial \Delta V_k},$ $A_{h,i,j,k}^{\delta }=\frac{\partial t_{h,i,j}}{\partial \Delta {\delta }_k}$(12.1)
$A_{h,i,j,k}^{ϵ}=\frac{\partial t_{h,i,j}}{\partial \Delta {ϵ}_k},$ $A_{h,i,j,k}^z=\frac{\partial t_{h,i,j}}{\partial \Delta z_k}$
等式(12.1)中的下标h、i、j、k可例如被定义如下(其它定义也可使用)：
·h是层位110的下标；
·i是为反射点128(例如在精细横向网格120上的射线对124和126的起点)的下标，它可定义在线上的和跨线的位置；
·j是张角的下标，它可定义张角的量值和方位；以及
·k是与给定系数关联的稀疏网格(例如，沿着它可定义模型扰动)上的模型节点的下标，它可定义层位、稀疏网格的跨线、稀疏网格在线上的位置。
前三个下标h、i、j可完全定义包括入射射线124和反射射线126的射线对。为简化各项，可引入全局下标l，使l＝{h，i，j}，其中l类似地完全定义该射线对。例如，在等式(12.1)中定义的四个矩阵可重写为：
$A_{l,k}^V=\frac{\partial t_l}{\partial \Delta V_k},$ $A_{l,k}^{\delta }=\frac{\partial t_l}{\partial \Delta {\delta }_k},$ $A_{l,k}^{ϵ}=\frac{\partial t_l}{\partial \Delta {ϵ}_k},$ $A_{l,k}^z=\frac{\partial t_l}{\partial \Delta z_k}---\left(12.2\right)$
矩阵中的每个元素与射线对(例如下标为l)和模型节点k关联，其中1≤l≤n，1≤k≤m，n是全部射线对的总数，m是模型输出节点数。由于通常射线对方位和反射点128的数量很大，n可以显著地大于m(例如，n＞＞m)。每个介质性质的改变(例如，层位的位移或另一个模型扰动)通常造成射线对走时的改变。然而，通过相应地改变在模型100每个节点的性质(例如，层位位移)，可以以补偿的方式(例如，提供走时的大小相等符号相反的改变)使每个射线对的总走时通常保持不变，例如如下式所示：
$\underset{k}{\Sigma }{A}_{m,k}^V·\Delta V_k+\underset{k}{\Sigma }{A}_{m,k}^{\delta }·{\mathrm{\Delta \delta }}_k+\underset{k}{\Sigma }{A}_{m,k}^{ϵ}·{\mathrm{\Delta ϵ}}_k+\underset{k}{\Sigma }{A}_{m,k}^z·\Delta z_k=0---\left(12.3\right)$
该等式描述了沿特定射线对路径(例如对特定下标m)的走时保持不变。如式(12.3)所示，对于改变后的参数值集ΔV，Δδ，Δε，Δz，走时的相应改变分别为走时的总变化等于这些贡献之和，对于时间保持不变的层析成像，所造成的走时残差消失。对于每个射线对，式(12.3)可写成矩阵形式，例如：
AV·ΔV+Aδ·Δδ+Aε·Δε+Az·Δz＝0   (12.4)
矩阵A可具有例如n×m的维度。由于n通常大于m(例如n＞＞m)，每个矩阵A可为矩形，该矩形的高度(例如n行)通常比长度(例如m列)大得多。
使用式(12.4)，如果给出其它三个参数值，可估计出每个未知参数，例如ΔV，Δδ，Δε，Δz。例如，可将等式(12.4)乘以相应的转置矩阵以便估计未知参数的值，例如如下式所示，：

例如，对等式(12.5)的矩阵乘积可引入下列符号：

等式(12.6)的矩阵Q通常为小方阵(例如两个维度相等且较小)。另外，矩阵Q通常为非负定矩阵。使用小的标准校正，这些矩阵可转换成正定矩阵。具有相同左右下标的四个矩阵QVV，Qδδ，Qεε，Qzz中的每一个通常是对称的(例如每个矩阵等于其转置矩阵)。另外，其它十二个矩阵Q的每对之间通常是对称的。例如，如果每个矩阵(例如QVδ)的左下标等于另一矩阵(例如QδV)的右下标且右下标等于所述另一矩阵的的左下标，则该矩阵相对于所述另一矩阵是对称的，例如如下式所示：
$Q_{\mathrm{V\delta }}=A_V^T·A_{\delta }={(A_{\delta }^T·A_V)}^T=Q_{\mathrm{\delta V}}^T---\left(12.7\right)$
等式(12.6)中的十二个矩阵可有例如下式所示的对称性：
$Q_{\mathrm{V\delta }}=Q_{\mathrm{\delta V}}^T,$ $Q_{\mathrm{Vϵ}}=Q_{\mathrm{ϵV}}^T,$ $Q_{\mathrm{Vz}}=Q_{\mathrm{zV}}^T$(12.8)
$Q_{\mathrm{\delta ϵ}}=Q_{\mathrm{ϵ\delta }}^T,$ $Q_{\mathrm{\delta z}}=Q_{\mathrm{z\delta }}^T,$ $Q_{\mathrm{ϵz}}=Q_{\mathrm{zϵ}}^T$
这样，只有十个矩阵Q需要存储：四个有重复下标的矩阵和六个有不同下标的矩阵。例如，利用式(12.8)中定义的关系，其它六个矩阵可通过各自的对称对导出。
结合等式(12.6)和(12.5)，例如给出：

矩阵Q是模板矩阵，通常由背景或标准模型100的性质完全定义。模板矩阵Q通常独立于模型参数值的输入和输出变化(例如扰动)，并通过改变模型参数值保持其不变。当矩阵Q定义大量数目n的射线时，添加额外射线通常只能稍微改变这些矩阵。
13.乔里斯基(Cholesky)分解
由模板矩阵Q定义的标准地下模型100可通过(例如右)乘以扰动矢量X来适应改变。当用扰动矩阵X乘以模板矩阵Q时，通常产生由改变后的参数值定义的新的被扰动的模型。
为模拟不同的模型，相同的(例如左侧)模板矩阵例如QVV，Qδδ，Qεε，Qzz(例如共同定义同一标准模型100)可乘以(例如右乘以)不同的右侧矢量。例如，用(例如三个)不同扰动矢量X1，X2，X3中的每一个乘以模板矩阵，可产生(例如三个)不同的新模型。通常应只改变扰动矢量(例如右侧)来产生模型改变。这样，可使用和存储单组模板矩阵作为用于快速和有效的建模算法的基础，以便产生同一数据的其它模型。
模板矩阵Q可被分解为下三角矩阵ST和上三角矩阵S的乘积。下三角矩阵ST是上三角矩阵S的转置，所以通常只存储上三角矩阵和下三角矩阵中的一个。模板矩阵Q可以是对称正定矩阵，例如由下式表示：
Q＝ST·S    (13.1)
可由矩阵方程来定义模板矩阵Q被扰动矢量X右乘，例如：
Q·X＝B     (13.2)
结合等式(13.1)和(13.2)，例如给出：
STS·X＝B   (13.3)
通过引入例如下列符号：
S·X≡Y     (13.4)
等式(13.3)中的S·X项可由Y替换(例如根据等式(13.4))，例如给出：
ST·Y＝B    (13.5)
在矩阵方程(13.1)被分成方程(13.4)和(13.5)之后，可进行两次向后替换：一次使用上三角矩阵S，一次使用下三角矩阵ST，例如如下式所示：
ST·Y＝B  求解Y    向后替换1
                                (13.6)
S·X＝Y   求解X    向后替换2
方程(13.1)可能求解得比较慢(例如具有大量的计算工作)，但通常对每类模拟只求解一次(例如，总共最多有四类模拟；模拟次数可基于在节点定义的参数个数)。方程(13.4)和(13.5)可被相对较快地解出(例如具有相对较小数量的计算)，且可顺序地求解或运算。例如，如果矩阵Q的维数是m且矩阵Q的带宽是b，那么解方程(13.1)通常需要大约mb2/3次浮点运算而解方程(13.4)和(13.5)中的每个通常需要大约mb次浮点运算。使用右应用顺序(right application sequence)，等式(13.6)可应用于方程(13.4)和(13.5)中的每个。
本领域的技术人员应该理解，本发明的实施例可使用除矩阵和/或矢量以外的其它数据排列以及除矩阵/矢量相乘以外的其它运算。具体而言，矩阵Q被矢量X右乘只是运算的举例，不意味着构成限制。例如，可使用其它方法求解方程系统。
扰动矢量X可包括用于保持沿每个射线对的走时不变的改变后的模型100参数值的残差值ΔVk，Δδk，Δεk，Δzk。k可以是模型网格120的稀疏网格间隔的下标。这四个参数值中的三个通常已知或被指定。等式(13.6)可用于确定扰动矢量X的未知部分。扰动矢量X的已知部分用于计算等式(13.6)的方程组的右侧。
14.从各向同性到各向异性
本发明的实施例提供了用于从标准的各向同性模型100构建各向异性模型100a的建模算法。将模型100的介质从各向同性介质改变为各向异性介质通常改变了沿模型射线对的走时。为了保持射线对的总走时，另一个参数，例如各向同性模型的层位110的位置，可以以补偿的方式位移以便提供大小相等符号相反或者相抵消的走时改变。
各向异性模型可由介质参数，诸如轴向压缩波速和两个汤姆森各向异性参数ε和δ来定义。初始时，可产生各向异性模型，该模型例如具有近似等于各向同性模型的压缩波速的轴向压缩波速，并具有都为零的汤姆森参数值。例如：
VP＝VIso，δ＝0，ε＝0      (14.1)
应该指出，即使ε和δ都为零值，哈密顿算子的导数和通常也不为零。由此，相应的层析成像系数通常是存在的。通常在张角大到30°之前，沿射线对的走时可保持不变。在其它实施例中，在张角大到其它(例如更大的)张角之前，走时可保持不变。
可以以给定的介质参数值，例如针对张角(例如可直到30°)，针对基本上所有方位角，来运行各向异性层析成像射线追踪模拟。在其它实施例中，例如，当走时针对其它(例如更大的)张角保持不变时，在张角大到所述其它(例如更大的)张角之前可运行射线追踪模拟。在又一个实施例中，可在超过使时间保持不变的张角的情况下运行模拟。
对于每个射线对(例如入射射线124和反射射线126)，可针对介质参数值和针对沿层位的震源点130和接收器点132和/或折射点128的垂直位移136，来计算层析成像系数。接着，可计算小的正定方阵Qzz以及其它(例如非对称)小方阵QzV，Qzδ，Qzε。方阵的维数可对应于模型节点112的个数。例如，维度为n×n的矩阵n＊n可对应于具有n个节点112的模型100。在每个节点112处可以存在多个(例如4个)参数值。这样，可以有表示模型100的多个(例如4n个)参数值。然而，一些参数值(例如3n个)可以是已知的。因此，可以计算剩余的(例如n个)未知参数值(例如，或参数值更新)。例如，模型100可具有：有n个参数值(例如模型的未知参数值的个数)的n×n维的n＊n矩阵。
可根据例如方程组(12.9)的最后一个方程来定义矩阵Qzz，它是：
Qzz·Δz＝Bz，Bz＝QzV·ΔV-Qzδ·Δδ-Qzε·Δε   (14.2)
可使用层析成像方法计算矩阵Q。矩阵Qzz可被分解成下三角矩阵和上三角矩阵的乘积，例如下式所示：
$Q_{\mathrm{zz}}=S_z^T·S_z---\left(14.3\right)$
在一个实施例中，可产生矩阵Q并将其分解一次。矩阵Q可被保存在例如盘上或在存储器中，以便为产生每个新的模型100a而重复访问和使用。
对于由标准模型100改变而成的每个新模型100a，可使用不同的参数值，造成残差ΔV，Δδ，Δε或输入变化的不同组合。可由式(14.2)计算矢量Bz。
结合方程(14.2)和(14.3)，例如给出：
$S_z^T{S}_z·\mathrm{\Delta z}=B_z---\left(14.4\right)$
通过对式(14.4)运行向后替换两次，例如：一次使用上三角矩阵Sz，一次使用下三角矩阵SzT，可得出垂直位移Δz 136(例如对应于将模型从各向同性模型改变成各向异性模型)。
15.利用井下研究预测汤姆森δ
本发明的实施例提供利用各向同性分析和井下信息对第一汤姆森参数δ的估计。各向同性分析给出层位节点zIso(例如反射点128和/或震源点130和接收器点132)的背景深度估计。使用井下信息，得到真实的或实验测量的深度zWell。然后，假定层位节点深度的深度残差为例如：
＝zWell-zIso        (15.1)
可假定汤姆森参数值的背景值例如为零。由方程组(12.9)的第二方程得到轴向速度残差和参数δ的残差之间的关系是：
Qδδ·Δδ＝Bδ，Bδ＝-QδV·ΔV-Qδε·Δε-Qδε·Δz  (15.2)
第二汤姆森参数ε的残差也可假定为例如零，Δε＝0。这样，式(15.2)的右侧简化为例如：
Bδ＝-QδV·ΔV-Qδz·Δz         (15.3)
按照不同模拟针对模型100的不同扰动来求解方程(15.2)(例如用不同的右侧扰动矢量去乘)，从而得到δ的残差Δδ。
16.应用
对于多个关键地球物理解释和生产任务，如不确定性分析、预测由于介质性质的不同表示造成的层位模型变化(例如从各向同性到各向异性，或从某种程度的各向异性到另一种程度的各向异性)以及确定由于模型层位位置和例如井中测量的标准层位之间的误关联造成的各向异性参数值等，时间保持不变的层析成像是特别有吸引力的解决方案。确定各向异性参数值还可用于在整个钻井过程进行快速和准确的交互式模型校正，以提供最佳的地质导向方案。本发明的实施例可具有其它好处和其它应用。
本发明的实施例可包括使用构造的层析成像矩阵(例如影响或模板矩阵Q)的特定对称结构的方法。这些矩阵的系数通常依赖于背景模型(如模型100)的性质，并且通常独立于模型扰动或模型参数残差。通过求解多个线性方程组(例如设置成模板矩阵Q)、利用不同的右侧矢量(例如设置成扰动矢量X)来确定标准模型，这类方法可以实现标准模型的改变。多个线性方程组中的每个方程可由模板矩阵Qij确定，每个模板矩阵是小的对称的正定方阵。在一些实施例中，模板矩阵可进行乔里斯基分解，例如通常只进行一次分解。然后可针对每个新的或改变后的模型，运行使用三角矩阵来求解方程组(例如，用参考方程(13.6)描述的向后替换)的处理。
本发明的实施例可实现面向任务的解释性询问。例如，本发明的实施例可用于解决以下问题：“由于汤姆森参数值和轴向压缩波速的给定扰动造成的层位节点的位置变化通常是怎样的？”；或“由于给定的层位位移和/或给定的轴向波速变化造成的汤姆森参数值δ的变化(例如汤姆森参数ε为零)通常是怎样的？”。可以假定造成的沿射线对的全双向路径的走时残差为零。
本发明的实施例可实现模型的改变或变换，例如，从各向同性介质到各向异性介质。例如，在各向同性介质(例如倾斜TI)的各向同性压缩波速与该介质的轴向压缩波速一致且例如汤姆森参数值δ和ε为零的情况下，各向同性介质可被认为是各向异性的。当前为零的汤姆森参数值δ和ε可被更新。本发明的实施例可实现模型层位的位置或介质参数值的不同的改变，并确定这些改变对其它参数的影响。具体而言，通过测量沿各向同性地震模型提取的层位与沿测井确定的标准层位之间的误关联，可确定汤姆森参数δ。
本发明的实施例可提供用于确定放大因子的不确定性分析。例如，在给定一个或更多个给定的模型节点的一个或更多个参数值中的小误差(例如，或者是不确定性)的情况下，则可确定在该模型的特定的一个或更多个节点或全部节点处的参数的不确定性。
本发明的实施例可实现地质导向。地质导向可包括调整井孔位置和方向以达到地质目标的处理。这些改变可基于在钻井期间收集的地质信息，并可在数据汇集过程中或“在线地”进行这种改变。该信息可包括沿井孔的介质性质的测量结果，如不同地质地层之间的过渡(例如标准层位)。使用时间保持不变的层析成像分析，可以更准确地进行这种对构造模型的“在线”调节以及井孔路径的相应调整。
17.用于改变地震数据模型的系统
参考图4，它示意性示出根据本发明一个实施例的系统。系统200可用于例如显示、转化或改变这里描述的地震数据模型。系统200可完成这里讨论的任何方法和/或其它运算或计算的具体实施。
系统200可包括发射器210、接收器220、计算系统230、显示器280以及输入装置290。发射器210(例如位于参考图1和图2描述的震源点130)可输出任何适当的信号，或产生入射信号。例如，可从多个位置中的每一个位置发射一系列声音或地震能量射线或波。接收器220(例如位于参考图1和图2描述的接收器点132)可接受对应于或关联于发射器210发送的入射信号的反射信号。在其它领域的成像情况中，例如医学成像，发射器210可输出能量，诸如超声波、磁能量、x-射线或其它适当的能量。
计算系统230可包括例如处理器240、存储器250和软件260。处理器240可处理数据，例如从接收器220接收的原始数据。存储器250可存储数据，例如原始的或处理后的地震数据，如表征数据集的通用标准模型的模板矩阵Q，通过该模板矩阵Q可产生数据的其它模型。存储器250可存储根据本发明实施例执行的操作的结果。例如可通过处理器240执行软件260来执行或计算根据本发明实施例执行的操作，例如产生模板矩阵、产生扰动矢量、诸如矩阵/矢量相乘的运算、求解前面提到的方程(5.1)至(15.3)中的任意一个、求解线性方程系统、二次方程、偏微分方程、常微分方程、射线追踪、转换模型等。其它单元或处理器也可执行根据本发明实施例的这些操作或其它操作。
显示器280可显示来自发射器210、接收器220、计算系统230或任何其它适当系统、装置或程序(例如成像程序或软件)或者发射器或接收器、追踪装置的数据。显示器280可显示根据本发明实施例执行的操作的结果。显示器280可包括一个或更多个输入，以显示来自多个数据源的数据。系统可包括多个显示器，例如每个显示器同时显示同一数据的不同模型，例如用于观察者的比较。显示器280可显示由数据产生的图像。
输入装置290可包括例如可由用户操作的鼠标或键盘，以录入用户的输入数据。也可使用其它输入装置。用户可将数据或命令输入到用户区域。
计算系统230可包括例如任何适当的处理系统、计算系统、计算装置、处理装置、计算机、处理器等，可通过使用硬件和/或软件的任何适当组合来实现。
处理器240可包括例如一个或更多个处理器、控制器或中央处理单元(“CPU”)。例如，软件260可部分地或全部地存储在存储器250中。例如，软件260可包括任何适当的软件，以便进行根据本发明实施例进行处理或成像。处理器240可至少是部分地根据软件260中的指令进行操作。
例如，系统200可以对目标地球物理区域进行成像，例如使用软件260和/或处理器240或其它部件，如专用的图像或信号处理器。18.用于改变地震数据模型的方法的流程图
参考图5，图5是根据本发明一个实施例的方法的流程图。
在操作300中，可产生参数值的初始集。该参数值的初始集可表示与地震数据集一致的初始模型(例如参考图1描述的模型100)。该地震数据集可表征地下或地球物理区域。可利用各种方位数据获取方法来获取该地震数据集。也可使用其它数据，如用于医学成像的数据。
可通过参数值集来定义模型。参数可以是模型的一种或一类性质，例如界面或层位110的位置和TTI介质的性质，如轴向压缩波速。参数值是相应模型参数的值(例如在节点112的参数值)。
例如，参数值的初始集可包括层位位置的值、汤姆森参数和轴向压缩波速。例如，参数值的初始集可包括由沿着模型射线对进行的射线追踪所产生的模型的层析成像系数。还可使用其它参数。
初始地震模型可用模拟同一地震数据集的各种扰动或模型的基础。初始地震模型可包括其每个都具有走时的多个射线对。初始模型可具有大小和具有大小的一组维度。
在一个实施例中，参数值的初始集可以以矩阵形式表示(例如这里讨论的模板矩阵Q)，或例如被表示成多个矩阵(如这里讨论的模板矩阵Qij)。
例如，可在参考图4描述的处理器240中产生参数值的初始集。参数值的初始集可存储在参考图4描述的存储器250中。
在操作310中，参数值的初始集可存储在存储器中，例如参考图4描述的存储器250，或存储在其它存储装置中，如盘驱动器。在每次将初始模型用于产生该地震数据集的新的改变后的模型时，可从存储器中重复地提取参数值的初始集。
在操作320中，可改变参数值的初始集中的两个或更多个参数值，从而对初始模型中的一些或全部射线对中的每个造成走时残差。
两个或更多个参数值可被改变，以响应用户对模型相应特性的调节，在操作350中将对此进行进一步的详细描述。
被改变的两个或更多个参数值中的每一个(例如单独地，或在不改变这两个或更多个参数值中的其余参数值的情况下)可对应于每个改变的射线对的走时变化。模型的两个或更多个参数值可被改变，分别造成每个射线对走时的第一和第二非零变化。被改变的参数值之一可对应于每个射线对的走时变化，该走时变化与其它两个或更多个被改变的参数值所对应的变化大小相等符号相反。
然而，两个或更多个参数值可以以互补的方式(例如组合地)被改变，对应于每个射线对的走时没有净变化。每个射线对走时的第一和第二非零变化(例如这里所描述的)可大小相等符号相反，因此组合地改变第一和第二参数值不造成沿每个射线对的走时的总改变。这样，可保持沿每个被改变的射线对的走时。
时间保持不变的层析成像机制可用于模拟初始模型的射线对中的一些或每个的不同扰动或改变。可计算层析成像系数以保持沿每个改变后的射线对的总走时不变。可产生两个或更多个参数值的变化，以响应例如经由参考图4描述的显示器280所提供的用户界面的用户输入。
为改变模型，通常用户不直接改变模型参数值。这一任务通常太费力。作为替选，例如用户可通过改变或将改变输入到用户界面区域中的定义中，来改变或操控模型的特征(features)或定义。用户界面区域可允许改变后的特征成为用户可见的模型表征或特性。用户界面区域可显示给用户(例如在计算机显示器280上)。在一个例子中，用户界面区域可被包括在模型本身的可视化显示中。可视化显示可以例如在监视器上呈现给用户，可视化显示可包括可修改的或者可改变的特征，如层位位置、反射点或其它空间模型对象，它们可由用户修改以使处理器改变模型。
例如。用户可输入命令、突出、选择、拖曳和/或点击，以将层位、反射点和/或其它空间区域从一个位置移动到另一个位置。在另一个例子中，用户界面区域可以是文本、复选项、框或值区域，以便输入用于改变模型特征的数据。在一个例子中，用户可对框进行复选以指出模型是各向同性还是各向异性和/或输入指出所希望的改变的模型的各向异性程度的数字。在另一个例子中，用户可输入轴向压缩波速改变1-2％。例如，用户可使用输入装置(例如输入装置290)，如鼠标或键盘。
这些对模型特征的改变通常是对模型100的子集或一部分局部地进行的。例如，用户可在模型的特定层或地质地层114改变轴向压缩波速1-2％。然而，通常这种局部改变可造成在整个模型中级联的走时改变，从而影响其它层。这样，对于模型一部分的每个局部改变，本发明的实施例提供一种机制，用于确定模型其它部分的相应全局改变，以保持沿整个模型的所有射线对的走时不变。
通过用户界面区域的这种局部改变可指示或触发处理器(例如处理器240)改变模型的地层的所有受全局影响的参数值，以定义新的模型。以这种方式，用户可通过在用户界面区域中输入对模型特性的改变来容易地修改模型，而不用考虑确定这些改变的剩余影响。
在一个实施例中，用户可输入对模型特性的两组不同的改变。作为响应，例如，可在用户界面上产生两个不同的模型(例如，表征同一地震数据集，但有不同的特征)以便用户比较。用户可从不同的模型中选择优选的模型。
在一个实施例中，例如，两个或更多个改变后的参数值可表示初始模型的一个或更多个结构层位的位置改变(例如，如参考图1描述的从层位110到层位110a的改变)。在另一实施例中，例如，两个或更多个改变后的参数值可表示初始模型的一个或更多个区域的介质性质的改变(例如，从各向同性到各向异性)。
在一个实施例中，两个或更多个参数值的改变可表示为矢量形式(例如这里描述的扰动矢量X)。当每个射线对的两个或更多个参数值由模板矩阵Q表示时，该矩阵可被扰动矢量X右乘。这种乘法运算可改变两个或更多个参数值。所述改变可被定制以采取例如根据用户输入的形式。
在操作330中，可产生由两个或更多个改变后的参数值表示的改变后的模型。由于每个模型由参数值集来定义，可通过改变参数值集将初始模型改变为改变后的模型。
改变后的模型可具有与初始地震模型的大小相同的大小。在一个实施例中，初始模型和改变后的模型中的每个可有一组维度，因此初始模型和改变后的模型维度的大小可以相同。例如，在一个实施例中，如果初始模型是n×m矩阵，则改变后的模型也是n×m矩阵。
对改变后的模型或改变后的模型本身进行定义的改变后的或最终的参数值集可被存储在例如存储器250中或其它存储装置中，如盘驱动器。
在操作340中，改变后的模型可显示给用户，例如在参考图4描述的显示器280上。
在操作350中，用户可再次调节特征，使处理器改变初始模型的两个或更多个参数值以产生第二改变模型。改变后的特征以及相应的改变后的参数值可不同于操作320中使用的那些。例如，可改变不同的特征，或以不同的方式或程度改变同一特征或改变为不同的值。再次调节可造成对另外两个或更多个参数值重复前述操作320-340。可模拟相应的不同模型扰动。可再次使用初始地震模型作为基础，以产生该地震数据集的第二改变模型。
可比较不同的模型，以选择最佳的建模参数值。用户可重复或“微调”这种再次调节，直至产生准确的模型。
尽管这里描述的实施例通常是指将初始模型改变为改变后的模型，但显然这种方法的实施例可应用于将任何第一模型改变为任何第二模型。例如，可重复地改变模型(例如，从第一模型到第二模型再到第三模型等)直至用户或自动控制机制确定该模型已足够准确为止(例如，根据已知的检验或测试达到预定的误差或变化范围)。在另一例子中，例如，通过选择用户界面区域中的复位按钮，改变后的模型可被改变回初始模型。
可使用其它操作或操作系列。
初始地震模型和改变后的模型均可具有一组维度。模型的这组维度可包括空间或坐标维度、参数维度和/或其它维度。例如，空间或坐标维度可由笛卡尔坐标x、y、z或极坐标定义。这组维度可对应于模型维度的度数(degree)或个数。例如，3D模型(例如3D立方体或棱柱)可具有三个空间维度，而2D模型(例如平面或3D模型的截面)可具有两个空间维度。在一个实施例中，改变后的模型可具有与初始模型相同的一组维度。用于初始模型和改变后的模型中的每个的这组维度可具有大小。模型在每个维度的大小可对应于被模型化的地下体在该维度中的大小。在一个实施例中，改变后的模型可具有维度组，其大小与初始模型的维度组相同。可使用其它维度、维度组和/或其大小。
本领域技术人员可以理解，本发明的实施例可应用于地震数据处理和建模中涉及的任何系统。本发明的实施例可用于在各种区域或领域的建模，如油气勘探和生产、用于环境研究的浅地球模型成像(例如使用由地震和/或探地雷达(GPR)方法收集的数据)、建筑工程(例如识别管道的位置)、建筑安全与防护(例如识别井孔和渠道)、医学成像(例如使用CT、MBI以及超声波装置)、无损材料探测、出于安全原因(例如国土安全)进行的对内部物体的探测、海中声纳、天线和雷达系统。
本发明的实施例可包括计算机可读介质，如存储器、盘驱动器或闪存，该计算机可读介质包括在被处理器或控制器执行时实现这里描述的方法的指令。
上文中提供的对本发明实施例的描述是为了进行说明和描述。并非旨在将本发明穷尽于或限制于所公开的具体形式。本领域的技术人员应理解，借助于上述内容，许多修改、改变、替代、变化和等效物是可能的。因此，应理解，所附权利要求旨在覆盖所有落入本发明真实精神内的修改和改变。