CN102906728B - 在仿真期间用于检查指示的方法和系统 - Google Patents

在仿真期间用于检查指示的方法和系统 Download PDF

Info

Publication number
CN102906728B
CN102906728B CN201180024861.8A CN201180024861A CN102906728B CN 102906728 B CN102906728 B CN 102906728B CN 201180024861 A CN201180024861 A CN 201180024861A CN 102906728 B CN102906728 B CN 102906728B
Authority
CN
China
Prior art keywords
checkpoint
association
storage
data
time
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.)
Expired - Fee Related
Application number
CN201180024861.8A
Other languages
English (en)
Other versions
CN102906728A (zh
Inventor
L·谭
J·E·安德森
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.)
ExxonMobil Upstream Research Co
Original Assignee
Exxon Production Research Co
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 Exxon Production Research Co filed Critical Exxon Production Research Co
Publication of CN102906728A publication Critical patent/CN102906728A/zh
Application granted granted Critical
Publication of CN102906728B publication Critical patent/CN102906728B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V99/00—Subject matter not provided for in other groups of this subclass
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/001—Acoustic presence detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/003—Seismic data acquisition in general, e.g. survey design
    • G01V1/005—Seismic data acquisition in general, e.g. survey design with exploration systems emitting special signals, e.g. frequency swept signals, pulse sequences or slip sweep arrangements
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/301—Analysis for determining seismic cross-sections or geostructures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/67—Wave propagation modeling
    • G01V2210/675—Wave equation; Green's functions
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/67—Wave propagation modeling
    • G01V2210/679—Reverse-time modeling or coalescence modelling, i.e. starting from receivers

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Retry When Errors Occur (AREA)
  • Processing Or Creating Images (AREA)

Abstract

本发明涉及用于例如在迁移（326）或反演地震数据中的使正向（328）和反向（308）传播的波互关联（316）中更有效的检查指示策略的系统和方法。检查指示策略包括在存储器中在已检查指示时间步存储正向仿真数据，其中已存储数据足以在该时间步实行互关联，但不重启正向仿真。在其它检查点，足以重启仿真的更大量的数据可以存储在存储器中（314）。披露了用于发现优化组合，即为给定量的计算机存储器（1004），并为在优化时间步定位检查点（306，1214，1310）使两种类型的检查点的计算时间最小化的优化组合的方法。优化检查指示策略也使快速（1402）对慢速（1404）存储的使用优化化。

Description

在仿真期间用于检查指示的方法和系统
关联申请的交叉参考
本申请要求在2010年5月19日提交的标题为“METHOD AND SYSTEM FORCHECKPOINTING DURING SIMULATIONS”的美国临时专利申请61/346,298的权益，其全部内容包含在此作为参考。
技术领域
本发明一般涉及地球物理勘探的领域，并更特别涉及地震数据处理。特殊地，本发明涉及用于存储检查点，从而改善例如在地震数据的迁移中以逆时(time reverse)顺序访问仿真数据的计算机仿真的效率的方法和系统。
背景技术
本节意图介绍可以与本技术的示范实施例关联的本领域各方面。该讨论据信帮助提供框架从而促进本技术的特别方面的更优理解。因此，应理解本节应据此阅读并且不必需作为现有技术的承认来阅读。
许多参数估计、反演和成像方法计算应用从初始状态在时间上正向外推的步骤序列的正向仿真。关联的反演或成像方法应用从最终时间状态时间反向步进的伴随算子。该方法可以然后通过将正向仿真与在相同时间步的伴随仿真互关联，生成关于系统的信息。生成的信息可以然后用来改善用来拟合可用数据的仿真参数以及用于其它目的。例如，该信息可以用来生成目标函数的梯度或Hessian关联改变从而生成仿真参数的改变。在另一例子中，该信息可以用来创造图像，如在逆时深度迁移(RTM)中完成。
这些技术为互关联计算需要以逆时顺序访问的时间正向(forward-in-time)仿真。然而，为所有时间步保存正向仿真需要的信息的大小经常超过可用存储器存储。
例如，叠前逆时深度迁移(RTM)普遍用来在实行烃勘探时用地震数据将地下地球结构成像。在RTM期间，在地震试验中记录的接收器波场时间反向传播，并与震源波场的时间正向仿真互关联。这意味着时间正向源仿真必须以逆时顺序访问。直接途径是计算并保存在全部地下位置的正向传播震源波场的全部时间步。由于涉及的海量数据，因此这可以是不实际的。例如，对于标准测试项目，勘探地球物理学家协会(SEG)先进建模(SEAM)项目，这可以意味着为总数约3.36拍它字节的数据保存84吉字节地下体积的40，000个时间步。不仅需要的存储量是有问题的，而且移动该量的数据需要的访问时间也可以是有问题的。
为减少需要的存储，在此称为全状态检查点的正向仿真波场状态可以定义为包括执行互关联或从给定时间步重启正向仿真需要的全部信息。全状态检查点可以然后在较小数目的仔细挑选时间步保存，因此正向仿真可以从已保存检查点重启，并向希望的时间步传播(如与在前述仔细挑选时间步保存的数据对照，术语检查点也可以用来指代该前述仔细挑选时间步)。因此，可以按需要从已保存检查点开始以逆时顺序在时间步重计算仿真。只要不再需要检查点存储器，那么其再使用于新检查点。折衷是执行更多计算从而最小化输入/输出和存储需求。只要涉及的数据大小比可用的存储器存储大得多，那么该技术可以是有用的。
Griewank和Walther提出给定其中保存检查指示(checkpointed)正向仿真波场状态的可用特定存储器量的挑选应检查指示的时间步的“优化”方式。见于例如Griewank，A.，“Achieving logarithmic growth of temporal and spatial complexity in reverseautomatic differentiation”，1 Optimization Methods and Software 35-54(1992)；Griewank，A.，Evaluating Derivatives：Principles and Techniques of AlgorithmicDifferentiation，Society for Industrial and Applied Mathematics(Philadelphia，PA，2000)；Griewank，A.和A.Walther，“Algorithm 799：An implementation ofcheckpointing for the reverse or adjoint mode of computationaldifferentiation”，26 ACM Transactions on Mathematical Software 19-45(2000)。
在SEAM例子中，可以假设在正向仿真中具有40，000个时间步。使用存储或存储器100在任何一个时间保存100个震源波场状态缓冲检查点，正向仿真可以使用114，747个正向仿真时间步，使用Griewank检查指示方案逆序访问。因此，以实行正向仿真的计算2.87倍的计算为代价，节省约400倍的存储空间。对于RTM应用，在慢速磁盘和快速计算的情况下，从检查点为正向仿真重计算时间步可以经常比实行访问已存储正向仿真时间步需要的输入/输出更快，即使充足的磁盘空间可用于存储正向仿真时间步中的全部。
Symes应用Griewank优化检查指示作为逆时迁移的有效实施策略。Symes，W.W.，2007，Reverse time migration with optimal checkpointing，72Geophysics(No.5)，P.SM213-SM221。由于使用波场检查指示的正向仿真的时间反转策略总是稳定的，因此检查指示法可以对包括波场衰减的物理现象的RTM是特别重要的(例如，使用P-波和S-波成像)。比较地，保存边界值和最终波场状态，并实行震源波场仿真的逆时外推的可替换RTM实施策略可以是不稳定的。进一步地，当包括衰减作为正向仿真中物理现象的部分时，在该技术中的不稳定性使其是不合适的应用策略。
用于实行正向仿真的时间反转的检查指示法具有比仅用于RTM多得多的一般应用。应用是非常一般的且与任何时间步进正向仿真器关联。这些包括储层仿真器，流体流动、热传递、地质盆地建模，以及地震全波场反演(FWI)。见于例如Tarantola，A.，1984，Inversion of seismic reflection data in the acoustic approximation：49Geophysics 1259-1266；Tarantola，A.，Inverse Problem Theory：Method for DataFitting and Model Parameter Estimation，Elsevier 125-133，144-258(1987)；Plessix，R.E.，“A review of the adjoint-state method for computing the gradientof a functional with geophysical applications”，167 Geophysical JournalInternational 495-503(2006)。Krebs，J R.，J.E.Anderson，D.Hinkley，R.Neelamani，S.Lee，A.Baumstein，M.D.Lacasse，“Fast Full-Wavefield Seismic Inversion UsingEncoded Sources”，74 Geophysics P.WCC177-WCC188(2009)。该技术也可以应用于自动差分的方法。见于例如Griewank，A.，Juedes，D.和Srinivasan，“ADOL-C，a package for theautomatic differentiation of algorithms written in C/C++”，Preprint MCS-180-1190，Mathematics and Computer Science Division，Argonne National Laboratory，Argonne，Illinois(1990)；Griewank，A.，Evaluating Derivatives：Principles andTechniques of Algorithmic Differentiation，Society for Industrial and AppliedMathematics(Philadelphia，PA，2000)。
发明内容
本技术的示范实施例提供用于降低计算机仿真的计算成本的方法。该方法包括执行检查指示策略，其中该检查指示策略包括在时间步存储关联检查点(在此有时称为“关联缓冲”)，其中该关联检查点包含用来使由计算机仿真生成的正向传播值关联的数据，并且其中计算机仿真不可以从存储在关联检查点中的值重启。检查指示策略也包括为关联检查点分配存储空间，并在多个时间步中的每个运行计算机仿真。在多个时间步中的每个，源自测量数据的反向传播值与源自计算机仿真的正向传播值关联。正向传播值存储在关联检查点中。
该方法也可以包括存储全状态检查点，并从该全状态检查点重启仿真。进一步地，该方法可以包括为存储检查点确定优化位置，其中该优化位置是在图形处理器单元(GPU)存储器、随机访问存储器(RAM)、RAM盘、磁盘驱动器或其任何组合中分配的存储空间。优化位置的确定至少部分基于存储空间的访问速度。用于关联检查点的存储空间可以小于用于全状态检查点的存储空间。优化化检查指示策略可以包括确定与存储关联检查点关联的计算成本。可以至少部分基于多个存储单元中的每个的访问速度，以及全状态检查点需要的存储空间对关联检查点需要的存储空间的比率将该计算成本最小化。包括时间步和在时间步中的每个存储的检查点类型的表格可以生成。
确定计算成本可以包括计算可以存储在快速存储单元中的关联检查点的最大数目，并计算可以存储在慢速存储单元中的关联检查点的最大数目。如果任何关联检查点可以存储在快速存储中，则计算如果仅使用快速存储的扫描优化数目并添加到列表。如果仅使用快速存储的扫描最小数目可以计算并添加到列表。如果任何关联检查点可以存储在慢速存储中，则如果使用全部类型存储的扫描最小数目可以计算并添加到列表。如果使用快速存储中的全部并在慢速存储中存储至少一个关联检查点的扫描优化数目可以计算，并如果小于列表上的最后值则添加到列表。最接近列表上的每个值的整数值可以添加到列表。可以计算列表上每个值的成本并且可以返回最小成本。
该方法可以包括从地震数据生成地下区域的图像。该方法也可以包括将储层数据与储层仿真历史匹配。
另一示范实施例提供用于比较收集数据与仿真数据的方法。该方法包括跨栅格反向传播收集数据并运行计算机仿真，从而跨栅格生成正向传播数据。数个关联缓冲可以在计算机仿真期间存储，其中该仿真不可以从关联缓冲中的任何重启。存储在关联缓冲中的数据可以和在栅格中每个点的反向传播数据比较。
该方法可以包括在仿真期间存储全状态检查点(在此也称为“全状态缓冲”)，其中仿真可以从该全状态检查点重启。该方法也可以包括为存储全状态检查点计算成本，并至少部分基于该成本为存储全状态检查点确定时间步。正向传播数据可以与下个时间步中的后向传播数据关联。
该方法可以包括在快速存储媒体中存储多个关联检查点，并在慢速存储媒体中存储至少一个关联检查点。进一步地，该方法可以包括存储至少一个全状态检查点，其中仿真可以从该全状态检查点重启。该方法可以包括在快速存储媒体中存储至少一个全状态检查点。
本技术的另一示范实施例提供用于使仿真数据与收集数据关联的系统。该系统包括处理器和存储系统，其中该存储系统包括表示测量数据和经配置跨栅格时间反向传播测量数据的传播算法。该存储系统也包括经配置跨栅格生成时间正向仿真数据的计算机仿真。该系统进一步包括存储器，其中该存储器包括代码从而引导处理器跨栅格时间反向传播测量数据，用源自计算机仿真的时间正向仿真数据填充栅格，在计算机仿真期间在时间步存储关联检查点，并使反向传播数据与存储在关联检查点中的仿真数据关联。
该处理器可以包括单核、多核、图形处理单元或其任何组合。关联检查点可以存储在存储器中。该存储器可以包括随机访问存储器(RAM)、RAM盘、图形处理单元存储器或其任何组合。该存储器可以包括经配置在时间步存储全状态检查点的代码，其中计算机仿真可以从该全状态检查点重启。
本技术的另一示范实施例提供用于迁移地震数据从而生成地震图像的方法。该方法包括跨地震栅格时间反向传播测量地震数据、跨地震栅格时间正向传播仿真激励脉冲，以及存储至少一个关联检查点，其中该关联检查点仅包含用于感兴趣互关联的仿真激励脉冲。使用源自关联检查点将测量地震数据的强度与在地震栅格中每个点的仿真激励脉冲互关联。
该方法可以包括在时间步存储至少一个全状态检查点，其中该至少一个全状态检查点包含从该时间步重启传播需要的全部数据。可以使用全双向波动方程执行仿真激励脉冲时间正向传播。可以基于基尔霍夫、波束或单向波动方程的迁移技术执行仿真激励脉冲时间正向传播。
附图说明
本技术的优点通过参考以下详细描述和附图更优理解，其中：
图1是感兴趣区域的图示，根据本技术的示范实施例图解在陆地环境中的单地震试验；
图2是根据本技术的示范实施例图解逆时迁移(RTM)的操作的示意图；
图3A是根据本技术的示范实施例示出用于执行逆时互关联例如RTM的方法的过程流程图；
图3B是根据本技术的示范实施例为逆时互关联示出波场的正向传播的过程流程图；
图4是在根据本技术的示范实施例图解在时间指数k使用全状态检查点的最小计算量中有用的图示；
图5A-B是在根据本技术的示范实施例图解时间步中有用的图示，全状态或关联检查点可以在该时间步存储；
图6A-C是在根据本技术的示范实施例图解具有55个时间步的简化仿真的Griewank和优化检查指示方案的操作中有用的图示；
图7A-C是根据本技术的示范实施例将标准关联的检查指示与在伴随计算中过去状态影响与正向计算中当前状态的互关联时的检查指示比较的图示，
图8A-E是在根据本技术的示范实施例图解数个状况中有用的图示，其中快速和慢速存储器的量变化，同时假设每个类型的检查点需要的存储空间相同；
图9A-E是在根据本技术的示范实施例图解数个状况中有用的图示，其中快速和慢速存储器的量变化，其中关联检查点需要的存储空间是二并且全状态检查点需要的存储空间是三；
图10是根据本技术的示范实施例的用于实施使用关联检查点的仿真的方法的过程流程图；
图11是根据本技术的示范实施例的用于在没有任何检查点的情况下计算最小成本的方法的过程流程图；
图12是根据本技术的示范实施例的计算用于使n个状态互关联的最小成本和优化位置的方法的过程流程图；
图13是根据本技术的示范实施例示出怎样预计算与在优化k位置取得的检查点关联的最小成本的表格的过程流程图；
图14是根据本技术的示范实施例的用于计算最小成本和优化检查指示策略，以便在起始检查点在StartType存储中的状况下使n个状态与nf_used快速存储器状态缓冲和ns_used慢速存储器状态缓冲关联的过程流程图；
图15是根据本技术的示范实施例示出在S/C＝12，并且Rs＝Ws＝Rc＝Wc＝0时使用四个缓冲的266个时间步的互关联的第一全状态检查点(k)的最小成本和位置之间关系的图表；
图16是根据本技术的示范实施例示出在S/C＝12，Rs＝Ws＝12，并且Rc＝Wc＝1时使用四个缓冲的266个时间步的互关联的第一全状态检查点的最小成本和位置之间关系的图表；
图17是根据本技术的示范实施例示出带有快速和慢速存储的优化检查指示策略超过仅有快速存储的优化检查指示策略和超过Griewank检查指示策略的性能改善的图表；
图18是根据本技术的示范实施例的可以用来实施本技术的示范集群计算系统的框图。
具体实施方式
在以下详细描述章节中，本技术的特定实施例连同例子实施例描述。然而，就以下描述对本技术的特别实施例或特别使用特定来说，其意图仅用于示范目的并简单提供示范实施例的描述。因此，本技术不限于在下面描述的特定实施例，但相反这样的技术包括落入附加权利要求的真实精神和保护范围内的全部替换、修改和等效。
起初，并为容易参考，阐述在本申请中使用的某些术语及其在该上下文中使用的意义。就在此使用的术语不在下面定义来说，应给予其如在至少一个印刷出版物或已提交专利中反映的关联领域技术人员给予的最广泛定义。进一步地，由于服务于相同或相似目的的全部等效、同义词、新发展和术语或技术认为在本权利要求的保护范围内，因此本技术不受下面示出术语的用法限制。
“算法”表达其正常意义，并无限制指代导致离散“值”的任何系列的可重复步骤。例如，算法可以包括任何数目的用户指定、预设、自动确定和/或工业或系统可接受的数据元素之间的任何数学、统计、位置和/或关系的计算。在若干实施例中，各种算法可以在主题数据上关于先前定义数据评估样本执行，以便产生单个的有意义数据值。
“计算机可读媒体”或“非暂时性的计算机可读媒体”如在此使用，指代参与提供指令到处理器以便执行的任何非暂时性存储和/或传输媒体。这样的媒体可以包括但不限于非易失性媒体和易失性媒体。非易失性媒体包括例如NVRAM或磁或光盘。易失性媒体包括动态存储器例如主存储器。计算机可读媒体的共同形式包括例如软盘、软磁盘、硬盘、硬盘阵列、磁带，或任何其它磁媒体、磁光媒体、CD-ROM、全息媒体、任何其它光媒体、RAM、PROM、EPROM、FLASH-EPROM、固态媒体如存储器卡、任何其它存储器芯片或盒式磁带，或计算机可以从其读取数据或指令的任何其它有形媒体。
如在此使用，“坐标系”是带有空间坐标(x，y，z)和时间坐标t的直角(笛卡尔)坐标域。地震数据通常在坐标系中收集和存储。空间坐标x和y可以表示正交水平坐标方向，例如其中获取数据的勘查的纵线方向和交叉线方向。空间坐标z通常表示垂直坐标方向，例如正向下测量的深度。
“互关联”是测量两个时间序列的数字有多少相互类似的过程。互关联可以是线性的或非线性的，可以考虑交互信息的方面(例如基于熵的交互信息)，或可以由认为对用户有用的任何其它技术执行。如在此使用，在时间正向投射的激励脉冲和时间反向投射的地震道之间的互关联用来迁移接收能量到地下中的反射点，鉴别位置(见于迁移)。
如在此使用，“显示”或“正在显示”包括导致显示物理物体的图形表示的直接行动，以及促进显示物理物体的图形表示的任何间接行动。间接行动包括提供用户通过其能够影响显示器的网站、超链接到这样的网站，或与执行这样的直接或间接行动的实体协作或合作。因此，第一方可单独操作或与第三方卖主协作，从而使信息能够在显示器设备上生成。显示器设备可以包括适合显示参考图像的任何设备，无限制例如虚拟现实显示器、3-d显示器、CRT监视器、LCD监视器、等离子设备、平板设备或打印机。显示器设备可以包括已用任何常规软件校准的设备，该常规软件意图用于评估、校正和/或改善显示效果。
“示范”在此专门用来意谓“用作例子、实例或说明”。在此描述为“示范”的任何实施例不解释为超过其它实施例优选或有利。
“地层”指代岩体或充分特别并连续以使其可以例如由地震技术绘图的其它地下固体。地层可以是主要是一个类型的岩体或类型组合的岩体。地层可以含有一个或更多含烃区。注意术语地层、储层和层段(interval)可以互换使用，但一般用来表示逐渐较小的地下区域、地带或体积。更特定地，地层一般是最大的地下区域，储层一般是在地层内的区域，并且一般是含烃地带(含有石油、气体、重油或其任何组合的地层、储层或层段)，并且层段一般指代储层的子区域或部分。含烃地带可以由较低渗透度的地带例如泥岩、页岩或页岩状(高度致密的)砂子从其它含烃地带分离。在一个或更多实施例中，含烃地带除砂子、泥土或其它多孔固体之外还包括重油。
“地震检波器”指代对地震信号灵敏的有源元件和支撑该有源元件的承载体(或结构)。有源元件通常包含压电元件，但也可以包括光学元件、显微机械加工的机电传感器元件等。
“烃”一般定义为主要由碳和氢原子形成的分子，例如石油和天然气。烃也可包括其它元素，例如但不限于卤素、金属元素、氮、氧和/或硫。烃可以通过穿透含烃地层的阱从烃储层产生。从烃储层获得的烃可以包括但不限于干酪根、沥青、焦性沥青、沥青烯、石油或其组合。烃可以在地球内位于矿物基质(mineral matrices)内或邻近矿物基质。基质可以包括但不限于沉积岩、砂子、沉积石英岩、碳酸盐、硅藻土和其它多孔介质。
“水听器”指代在水体环境中对压力(声音)波灵敏的有源元件。水听器也包括支撑该有源元件的承载体(或结构)。有源元件通常包含压电元件，但也可以包括其它声检测元件。
“阻抗”是地震波速和密度的乘积。阻抗通常在不同岩层之间变化，例如分界面的相反侧面具有不同阻抗。两个类型的阻抗术语Ip和Is一般定义，其中Ip是P-波阻抗，也称为声阻抗，并且Is是S-波阻抗，也称为剪切阻抗(shear impedance)。分界面的反射系数一般取决于这些阻抗中的对比度和分界面任意一侧上的岩石密度。特定地，在地质层的这些性质中的对比度在分离两个层的边界影响反射系数。基于记录的地震反射数据确定地下区域的阻抗和/或密度结构的一个地球物理过程是地震反演。
“反演”或“地震反演”是人尝试发现在公差内重产生测量地震响应的地下性质的模型，并满足地质和地球物理约束的过程。具有大量众所周知的地震反演方法。这些众所周知的方法落入两个类别，迭代反演和非迭代反演中的一个。非迭代反演通过假设一些简单背景模型并基于输入数据更新该模型来实现。比较地，迭代反演使用已更新模型作为反演的另一步骤的输入(见于迁移)。
“基尔霍夫迁移”是依靠信号的统计相长干涉和噪声的相消干涉的反演的逆反向散射法。其为两步操作，首先从每个深度点向上投射或射线追踪到地表，并建立到地表位置的潜在射线路径的走时表。然后在基于如由走时表定义的样本源和接收器位置的时间为每个周围轨迹将样本求和。基尔霍夫迁移可以基于单向波动方程，并可以相比全双向波动方程法，例如在此讨论的逆时迁移(RTM)较不准确。
“水体环境”指代任何离岸位置。离岸位置可以是在浅水中或在深水中。水体环境可以是海洋体、海湾、大湖、河口、海或水道。
“迁移”通常在地震成像的数据处理阶段期间执行，以便准确定位地下地震反射体。因为可变地震波速和倾斜反射体在未迁移的地震图像中导致地震反射从而在不正确位置出现，所以出现地震迁移的需要。地震迁移是其中地震反射移动或“迁移”到其真实地下位置的反演操作。具有许多不同的地震迁移技术。这些迁移技术中的一些在数据轨迹的同中点(CMP)堆叠之后应用。如在本领域中已知，CMP堆叠是数据处理进程，其中数个地震数据轨迹具有相同的震源-反射体终点，但不同的偏移求和从而形成为所讨论中点将零偏移数据轨迹仿真的堆叠数据轨迹。这样的“叠后”迁移可以例如通过沿衍射曲线(称为“基尔霍夫”迁移)积分、通过数值有限差分或波场的相移向下延拓，或通过在频率-波数或其它域中的等效操作来完成。
其它地震迁移技术可以在地震数据轨迹的堆叠之前应用。这些“叠前”迁移技术应用于个别非零偏移数据轨迹，并且迁移结果然后堆叠从而形成最终图像。叠前迁移通常产生比叠后迁移更好的图像。然而，叠前迁移一般比叠后迁移昂贵得多(即，计算需求高)。因此，叠前迁移的使用通常限于其中叠后迁移不提供可接受结果的状况，例如其中反射体陡峭倾斜。在一些情况下，反射体倾斜可以超过90度。如在地震勘探领域中已知，可能使用源自地震“转弯射线(turning rays)”的数据将这些“倾覆的”反射体成像。
具有两种一般类型的叠前迁移，叠前时间迁移和叠前深度迁移。在地震成像中需要用来描述在地下中地震波传播速度的背景地震波传播速度模型。在其中地下地震波速仅在垂直方向上变化的区域中，使用的地震成像法是叠前时间迁移(PSTM)。在其中地下地震波传播速度在垂直和横向(或水平)方向上变化的情况下，叠前深度迁移(PSDM)需要用来给出准确结果。
“模式转换地震波”是具有与震源地震波不同模式的反射P-波或S-波。地震波具有两个模式。在第一模式中，地震波是在地震波的传播方向上延伸的P-波(或压缩波)。在第二模式中，地震波是一般垂直于地震波的传播方向定向的S-波(或剪切波)。由震源(例如，爆炸，气枪等)产生的震源地震波引入可以具有地下反射体(例如，到烃层、水层或感兴趣的另一层的分界面)的地下结构。地震波从地下反射体反射，其中反射地震波可以由地震传感器(或复数个地震传感器)(例如，地震检波器或其它地震传感器)测量。在许多情况下，模式转换可以紧接着源自地下反射体的反射发生。例如，震源P-波可以反射为S-波，或可替换地，震源S-波可以反射为P-波。具有不同于震源地震波的模式的反射地震波称为模式-转换地震波。根据一些实施例的地震勘查技术为模式-转换地震波估计吸收参数，这允许地震勘查技术考虑与模式-转换地震波关联的吸收效应。考虑反射模式-转换地震波的吸收效应，以及反射纯(或单模)地震波的吸收效应，允许地震勘查技术比常规地震勘查技术(其通常仅计算纯波吸收参数)更准确。
“多分量勘查”指代使用记录在接收器上入射的地震能量的两个或更多分量的多分量接收器的勘查。例如，3分量(3-C)地震接收器含有三个正交地震检波器，并因此可以在接收器记录质点运动的x、y和z分量(质点运动可以是质点位移、质点速度或质点加速度，以至在原理上可以是质点位移的高阶导数)。在水体地震勘查中，4分量(4-C)地震接收器可以可替换使用。4-C接收器除了三个正交传感器之外还含有压力传感器例如水听器，并因此除了质点运动的x、y和z分量之外还可以记录水柱的压力(其为标量)。
“单向波动方程”是近似描述波场的向下或向上传播的空间方向，但不同时描述该两者的发展方程。相似方程存在从而涉及各种类型的电磁数据。
如在此使用，术语“优化的”、“正在使优化化”、“使优化化”、“优化性”和“优化化”(以及语言学上涉及的单词和短语)不意图在需要本发明发现优化解从而做出决策的意义上限制。尽管数学优化解可以实际上达到全部数学可用的可能性的优化，但优化化路线、方法、模型和过程的现实实施例可以向这样的目标工作而不曾实际实现优选。因此，具有本披露的益处的本领域技术人员认识到这些术语在本描述范围的背景下是更一般的。该术语可以描述向解工作，其可以是优化的可用解、优选解或在约束的范围内供应特定益处的解；进一步地，该术语可以指代改善或精炼当前解。该术语也可以为目标函数搜索高点或最大值，或处理从而减少补偿函数。
“像素”意思是在图像或仿真栅格中的最小可寻址点。例如，在地震成像中，在地震栅格中的一个计算单元可以表示在图像数据中的一个像素。
“叠后”指代在个别传感器记录已求和或堆叠之后的处理。爹后属性包括例如反射强度、瞬时频率、反射非均匀性、声阻抗、速度、倾角、深度和方位角。
“PP”和“PS”表示在地震勘查期间发生的模式转换事件。由震源发射的声能可以主要是压力波(或P-波)。当声能在分界面经历反射时，其可以也经历部分模式转换成剪切波(S-波)。结果，在接收器取得的地震波场可以因此含有P-波和S-波。由于起因于P-波到达的事件涉及作为P-波发射并在接收器作为P-波记录的声能，因此该事件可以称为PP事件。由于起因于S-波到达的事件涉及作为P-波发射，但紧接着反射经历模式转换成S-波并因此在接收器上作为S-波记录的声能，因此该事件可以称为PS事件。PP事件在所获取地震数据的垂直分量中更主要发生，而PS事件在所获取地震数据的水平分量中更主要出现。
“P-波”和“S-波”是在地下发生的地震波的类型。P-波是作为交替压缩和膨胀波用垂直于波前的质点运动传播的纵向压缩波。P-波接收器响应关于地球表面的垂直质点运动。S-波或剪切波平行于波前极化，并为同位介质(isotopic media)分类为SH波和SV波。在本披露的上下文中，对于SH波，质点运动在横向于剖面线的平面中是水平的。，SV波的质点运动在垂直于SH质点运动并平行于剖面线的垂直平面中。因为流体没有剪切强度，所以剪切波不可以在流体中传播。因为一些介质是非均匀的，所以其对于S-波是双折射的。即，在传输通过介质并用不同于纯SH或SV波的极化方向传输期间，声能分裂进入由不同传播速度表征为快速(Sf)和慢速(Ss)波的普通和特别射线路径。
P-波探测包含地震勘查的主要设备。然而，特殊研究，例如使用在此描述的技术，可以允许由于应力和断裂的所选择岩层的非均匀特性的另外探测。这些研究可以通过在单独勘查中组合S-波技术和P-波技术来执行。
“接收器”是通常放置在地球表面上或刚好在下面的阵列或栅格状图案中，用来从岩层检测振动的反射的设备。在众多位置的正在到达的反射波的振幅和到达时间的测量值允许岩层的映射，并提供关于岩层(例如层)的厚度和组成的信息。接收器可以包括地震检波器、水听器、振动检测器、加速度计，或能够准确测量反射波的振幅的任何其它检测器。三维接收器可以例如沿x、y和z轴的每个具有加速度计，从而确定全部三个方向上的运动。四维接收器可以将三维接收器与振幅检测器例如水听器组合，从而确定P-波的强度。三维或四维接收器都可以用来检测从地下传播的P-波和S-波。
“储层”或“储层地层”定义为包括砂岩、石灰岩、白垩、煤和一些类型页岩的产油层(pay zone)(例如，烃生产层)。产油层可以在厚度上从小于一英尺(0.3048m)到数百英尺(数百m)变化。储层地层的渗透度为生产提供潜力。
“储层性质”和“储层性质值”定义为表示含有储层流体的岩石的物理属性的量。术语“储层性质”如在本申请中使用，包括可测量和描述性的属性。可测量储层性质值的例子包括对p-波的阻抗、对s-波的阻抗、孔隙度、渗透度、含水饱和度和裂缝密度。描述性储层性质值的例子包括外观、岩石学(例如砂岩或碳酸盐)，以及沉积环境(EOD)。储层性质可以构成计算单元的储层构架从而生成储层模型。
“岩石物理模型”将岩石地层的岩石物理的和生产关联的性质与岩石地层的体积弹性性质相联系。岩石物理的和生产关联的性质可以包括但不限于孔隙度、孔隙几何形状、页岩或泥土的孔隙连通性体积、估计的上覆岩层应力(overburden stress)或关联数据、孔隙压力、流体类型和含量、含泥量、矿物学、温度以及不均匀性，并且体积弹性性质的例子可以包括但不限于P-阻抗和S-阻抗。岩石物理学模型可以提供可以为地震勘查用作速度模型的值。
“地震衰减”是在地震波远离震源传递时由于微观摩擦力和源自薄层的散射的地震波中振幅或能量的频率关联的减少。其经常以地震质量因数Q描述。地震衰减受流体饱和度、含泥量和薄成层影响。如果充分衰减发生，那么由地震数据的增加的空间采样(例如，测量P和S-波)提供的改善可以是最小的。地震质量因数Q估计对地震成像、处理和储层表征应用是有价值的。这样应用的例子包括振幅和相位补偿、小波处理、获取设计，以及岩石学/流体识别。此外，不同于许多其它地震属性，地震衰减可以直接涉及渗透度(例如，经喷射流机制)。在本技术的示范实施例中，岩石物理学模型可以与成像中的近来进步组合，从而确定在估计地震质量因数和储层参数之间的连结。
“地震数据获取”意思是地震波的感测、处理和记录。
“地震数据轨迹”是所选择地震属性(例如，地震振幅或声阻抗)在单独x-y(映射)位置的记录。地震轨迹可以表示为单元的堆叠，或由振幅在每个z(或t)数据点为所讨论的x-y位置反映属性值的连续曲线(称为“轨迹”)表示。对于三维或四维接收器，为每个维度收集轨迹。因此，清晰理解在下面描述的技术应用到在维度的全部中收集的轨迹。然而，为解释简单，在此的描述集中在每个接收器的单独轨迹上。
“地震数据”指代含有关于使用地震方法获得的关于矿区地下结构中的点的信息的多维矩阵。在地震数据中的至少三个维度可以表示空间-时间中每个点的位置，例如x、y和t，其中x和y表示在地表上x、y栅格中点的位置，并且t表示反射弹性波到达地表的时间。取决于地质层的性质，t可以一般表示在地球表面下面的点的深度。多位矩阵含有表示在每个位置的信号值的至少一个其它维度。例如，该维度可以表示可以连接到矿区中特别点的任何地震数据，包括地震数据，尤其例如偏移距道集(offset gathers)、偏移堆叠、角道集、角堆叠Ip、Is、Vs、ρ和μ。
“地震反射勘查”是最普遍类型的地震勘查，通过在地球表面发动震波并在多个地表位置监视源自下伏地下地层的该扰动的反射来执行。这些反射从其中具有地球的声阻抗变化的一般是邻近地层之间分界面的区域发生。用来监视反射的设备称为地震检波器或水听器。由每个地震检波器记录的信号表示为由该地震检波器检测到反射振幅的时间的函数。为良好近似，由每个地震检波器检测到的反射从位于经过震源和地震检波器之间中点的垂直线上的每个反射面上的点发生。因此，对于每个地震扰动(“发射”)，每个地震检波器记录表示在地球表面上已知点垂直下方的地层特征的信号(“轨迹”)。
“结构化栅格”是其中每个单元都可以在二维中由指数(i，j)，或在三维中由指数(i，j，k)寻址的栅格。结构化栅格的全部单元具有相似形状和相同数目的顶点、边缘和面。这样，栅格的拓扑结构(即，在单元、面、边缘和顶点之间的边界和邻接关系)完全通过指数化来定义(例如，单元(i，j)邻近单元(i+1，j))。最普遍使用的结构化栅格是每个单元在二维中具有四个边缘或在三维中具有六个面的笛卡尔栅格或径向栅格。结构化栅格通常与地震数据体积一起使用。
概述
在本技术的示范实施例中，实施在存储空间用来保存包括在给定时间步重启正向仿真器需要的全部信息的完整正向仿真状态(全状态检查点)对仅用于实行关联计算(关联检查点)之间区别的优化检查指示方案。进一步地，本技术为多个类型的存储在访问速度和存储容量差异之间区别，从而优化化已保存的全状态和关联检查点的混合。对于在集群计算环境中的地震逆时深度迁移(RTM)和全波场反演(FWI)应用，与本发明关联的超过传统Griewank检查指示策略的计算加速可以范围从百分之二十到超过百分之三百。然而，在此描述的本技术不限于RTM，或实际上任何地震成像应用，但可以在需要使以逆时顺序呈现的第一数据流与以正时顺序呈现的第二数据流关联的任何数目的应用中使用。
Griewank检查指示策略不包括保存关联缓冲。Griewank仅保存全状态检查点。为精确，Griewank为关联并为重启仿真使用全状态检查点。在此披露的本技术通过利用执行关联需要的存储可以在访问速度和存储大小上低于存储完整正向仿真状态需要的存储，来获得相对于Griewank检查指示策略的优点。例如，完整地震数据记录可以具有记录的九个分量；三个表示在三维的每个中的质点速度；以及六个应力张量场，关于P-波和S-波的信息可以从该六个应力张量场提取。在地震正向仿真中，全状态检查点包括重启正向仿真时间步的计算需要的全部参数和变量，这可以需要保存许多波场分量的密集栅格加上实施边界条件需要的另外缓冲。见于例如Marcinkovich，C.和K.Olsen，“On the implementation Ofperfectly matched layers in a three-dimensional fourth-order velocity-stressfinite difference scheme”，108 JOURNAL OF GEOPHYSICAL RESEARCH(No.B5)2276(2003)。
相反，关联波场(关联检查点)可以仅在全模型体积的子集中的较小栅格上需要，并可以仅包括波场分量的子集，并可以需要较低的数值精度。例如，如果在RTM中获得P-P图像是特定应用的唯一目标，则可以仅需要使一个波场(即，波场的压力分量)关联，即使九个分量用于计算。作为另一例子，可以为计算希望双精度，但互关联仅需要单精度。
由于计算资源变得更复杂，例如为特别应用优化化，因此本技术也提供计算在保存全状态检查点和关联检查点之间的优化平衡。在示范实施例中，本技术在存储器大小和速度之间区别，并因此，超过Griewank技术的节省可以在包括技术例如图形处理器单元(GPU)和固态硬盘驱动器(SSD)以及传统随机访问存储器(RAM)和磁盘存储的层次环境中更显著。因此，使用本技术开发的应用软件在变化的计算机环境中可以是可移植的和灵活的。
已推导的优化检查指示策略一般应用于时间步进方法。因此，其可以用于如实行时域全波形反演以及逆时深度迁移需要的伴随状态梯度计算。然而，如上面提到，RTM关于单震源激励或“试验”更明确解释，如在图1中图解。
图1是感兴趣区域100的图示，根据本技术的示范实施例图解在陆地环境中的单地震试验。本领域技术人员理解巨大数目的地震试验在感兴趣区域100上执行，从而使地下结构有效成像。在图1中示出的感兴趣区域100中，盐丘102向地表104突出。盐丘102的侧面可以是垂直的以至仰拱，例如侧面106，使得难以由许多技术成像。
盐丘102可以捕集在冠岩层110下面的油砂或多孔岩石中，位于盐丘102上面的烃沉积物，例如石油108。其它沉积物例如石油或天然气112可以位于沿盐丘102的侧面例如侧面106的岩层114中。更其它的烃例如沉积物116可以位于更远离盐丘102的岩层118中。
地震试验可以通过沿感兴趣区域100的地表104放置一个或更多震源120与一个或更多接收器122来执行。在实施例中，接收器122的每个可以是能够检测P-波和S-波的三维或四维接收器。在触发时，震源120创造脉冲，该脉冲可以传导进入地下124作为波126。尽管脉冲波126在全部方向上传播，但图1中的图解被简化，例如仅示出脉冲波朝向盐丘102的传播。脉冲波126的能量的部分可以没有从任何表面反射，但代替地直接行进通过地下124到达接收器122，例如作为直达波128。
当脉冲波126在反射体，例如带有不同阻抗的层如盐丘102上冲击时，可以从波126创造反射波130。反射波130可以部分或完全模式转换，例如冲击P-波可以具有S-波分量。如果反射波130具有S-波分量，那么S-波可以由含有充实浓度流体例如沉积物116的岩石衰减。该衰减可以称为Q衰减。S-波的Q衰减的有效确定可以对增强地震成像的准确度、允许层例如支撑流体例如烃的层例如岩层118、含水层等的定位。较简单的成像技术经常仅为成像使用P-波反射。然而，在本技术的示范实施例中，P-波和S-波都用于成像，允许Q衰减的确定和量化。
在RTM中双向波动计算的使用允许确定其它效应。例如，脉冲波126的一部分可以不在盐丘102的表面反射，但可以代替地折射进入盐丘102。折射波132可以在由盐丘102的另一表面朝向接收器122折射返回之前行进通过盐丘102作为波134。较简单的成像技术，例如基尔霍夫迁移或基于单向波动方程的迁移不能捕捉该各向异性。因此，这些技术可以丢失地下124的分量，例如陡峭反射体(侧面106)、折射地下层(盐丘102)等。
为地震成像计算的目的，感兴趣区域102可以分为计算单元的三维栅格。在地震计算中，栅格一般是规则栅格，允向计算单元分配表示位置的坐标。由于可以使用非结构化计算网格，因此本技术不限于规则栅格。然而，由于每个计算单元可以含有位置信息，因此非结构化计算网格可以添加显著的开销。每个计算单元都可以表示或含有具有关联性质，例如速度、阻抗等的地震数据体积。进一步地，地震波前例如P-波或S-波的速度和方向可以在每个地震数据体积中使用该地震数据体积的性质计算。实际上，迁移过程“移动”震源和接收器，从而将反射与其中该反射发源的计算单元关联。本技术可以用来生成地下的图像。通过沿震源和接收器检查感兴趣区域100的剖面，RTM过程可以更明显可见，如关于图2进一步详细讨论。
图2是根据本技术的示范实施例图解逆时迁移(RTM)的操作的示意图200。计算单元的每个可以在示意图200中图解为计算单元202。计算单元202可以不是地震数据在其上收集或建模的最小大小，但可以表示例如已组合从而简化计算复杂度，或例如已分离从而改善图像准确度的较小体积的群。计算单元202的大小在下面更详细讨论。
在基本形式中，RTM分离处理每个地震试验，并通过将用来收集数据的物理过程仿真来迁移数据。如在图2中所示，算法利用正在成像的区域的波传播速度的速度模型206。速度模型206一般包含在计算单元202的每个中的岩石性质，例如密度、阻抗等，该岩石性质可以由波传播算法用来计算波的速度和方向。
通常是有限差分或伪谱的波传播算法可以用来从震源208的位置传播脉冲进入地下。如在此讨论，这可以使用在比计算单元202更精细的时标上的计算单元执行。该传播可以提供可以从每个地下反射体例如盐丘210反射的能量的近似作为时间的函数，并在接收器212记录。由于震源波场214主要向下行进，因此其可以标记为D(xS，x，t)，其中xS是表示震源位置的在三维空间中的向量，x是表示图像点216位置的在三维空间中的向量，以及t是表示从地震试验开始的时间的值(即，t0表示脉冲波的初始生成)。
正当震源波场时间正向传播从而使反射体的照明仿真时，在接收器212记录的数据轨迹218可以时间反向传播从而使能量从迄今未知的反射点到接收器212的运动仿真。后向传播可以由使轨迹作为波场220映射进入地下的伴随算子执行。由于接收器波场220一般向上行进，因此其可以标记为U(xS，xR，x，t)，其中xS是表示震源位置的在三维空间中的向量，xR是表示接收器位置的在三维空间中的向量，x是表示图像点216位置的在三维空间中的向量，以及t是表示从地震试验开始的时间的值(tmax表示数据收集试验的结束)。
在RTM中使用的成像原理是由接收器212同时记录的能量必须已从位置例如图像点216反射，以使震源波场214非零，并且向后传播的接收器波场220在相同传播时间t非零。RTM将已由倒转能量照明的位置成像的能力将RTM从其它基于波动方程的迁移算法区别开。为在给定地下位置x为在xS的给定发射计算图像，在已接收波场同样非零时的时间t必须具有在源自震源的能量可以从x反射并到达接收器时非零的量。试验的总时间是tmax，因此事件时间t范围从0到tmax。这在数学上可以表达为震源波场214和接收器波场220的互关联。
在方程1中，I(x)是在地下位置x，例如图像点216的迁移图像。在方程1中的加权函数表明震源波场214(D)沿时间t正向传播，并且接收器波场220(U)可以从最终时间tmax反向传播(因此，名称是逆时迁移)。两个波场214和220的乘积(例如，互关联)可以在每个图像点216(x)求和。这为每个震源位置xS和每个接收器位置xR完成。产生的图像可以为全部震源和接收器位置求和从而提供最终图像。
成像条件的计算需要震源波场214和接收器波场220的值在相同时间t可用。最直接途径是为感兴趣的全部时间t(例如t0到tmax)计算并存储波场中的一个(例如震源波场214)。其它波场(例如接收器波场220)可以然后计算并与震源波场214相乘。然而，对于大多数问题，由于存储震源波场214需要大量数据存储，因此该途径可以是不实际的。因此，RTM的实际应用需要人周期性重计算震源波场214。
其中x是位置向量(x，y，z)的图像I(x)表示地下的反射率作为位置的函数。反射率在其中阻抗可观改变的位置是巨大的。然而，如从方程1可见，RTM算法在其中具有在震源和接收器波场214和220之间的任何关联的任何点产生非零图像值。这可以在除对应于真实反射的位置例如图像点216之外的位置发生，导致将作为在真实反射体上叠加的低频噪声出现的伪迹成像。用于移除这些伪迹的数个方法，尤其低阻滤波、拉普拉斯滤波的应用等可用于该目的。
如先前提到，RTM使用全波动方程使能量传播通过地下，并因此其是计算上昂贵的迁移算法。由于成本有效的计算集群的出现，因此其中个别发射分离迁移的叠前RTM仅在近几年变得计算上可行。通常，数千个多CPU节点的集群可以链接在一起从而形成计算集群。有待迁移的每个发射可以分配到集群中的一组节点，并且围绕该发射的区域的合适模型(速度、密度等)可以读入与已分配节点关联的存储器。可以然后使用在上面描述的波动方程的正向和反向解使该发射迁移，并且结果可以在磁盘上存储，从而与其它迁移发射的结果求和。在本技术的示范实施例中，优化化方案用来计算在计算机存储器和计算开销之间的折衷。如下面讨论，该方案可以存储全状态检查点(关联可以从其执行，或传播操作可以从其重新开始)、关联检查点(仅包括在特别时间的互关联计算需要的数据)，或其组合。检查指示的使用使得RTM的实际实施可行。进一步地，检查指示对考虑非弹性传播效应(“Q”)或计算许多波场分量例如S-波场和P-波场的实施例甚至更重要。
计算模型可以分为计算单元202或较小体积204。计算单元的间距使得每最小波长具有若干栅格块。对于通常的模型，间距可以是约25米。因此，由于感兴趣区域的模型具有约几十千米的尺寸，因此在每个维度上栅格单元的数目可以是几百到多于一千。因此，计算单元202的总数可以是例如大于十亿。计算可以通过在整个栅格上运行有限差分模板来执行，以便计算在每个计算单元202将波动方程建模需要的波场导数。因此，可以在每个时间步执行数十亿个浮点运算。
需要的栅格间距Δx，例如计算单元202或较小体积204的大小可以使用在方程2中示出的公式计算。
在方程2中，vmin是最小速度，fmax是最大频率，并且N是每最小波长的点数目，其可以是约5-10个点。该间距或更小间距可以避免数值弥散引入传播波场。
时间步的总数也可以是非常巨大的。地震数据在约十秒的总时间记录。该数据通常用2-4毫秒的简单间隔存储，意味着每个轨迹可以包括约2500-5000个样本。然而，RTM波动方程解必须在由方程3中示出的Courant-Friedrichs-Levy(CFL)稳定条件给出的精细得多的步骤上计算。
在方程3中，k是取决于使用的维数和精确算法的一阶常数。例如，对于二维方程k可以是或对于三维方程k可以是进一步地，如果方程在时间上是更高阶的，那么k可以更大，或如果方程在时间上是更低阶的，那么k可以更小。对于通常条件，Δt～0.5msec，并因此，十秒的总传播时间需要20000个时间步。这表明每个昂贵的时间步必须计算非常多的数目。此外，使用方程2和3，可以注意到总计算成本与f4 max成正比，因此将有待迁移的最大频率加倍使计算成本提高16倍。由于在频带宽度较大时地震图像的解释能力通常较优，因此成本对频率的依赖性趋向于迅速提高获得高质量图像需要的时间。
图3A是根据本技术的示范实施例示出用于执行逆时互关联例如RTM的方法的过程流程图。在该实施例中，该过程可以在方框302用逆时迁移(RTM)关联图像的初始化来开始，例如通过设定图像体积中的全部点为零。在方框304，伴随波场可以初始化，例如通过设定伴随波场在最大记录时间匹配接收器波场。在方框306，可以计算存储检查点缓冲和/或全状态缓冲的优化时间位置。优化时间位置也可以在运行仿真之前计算，并存储在稍后仿真中使用的表格中。
伴随波场可以然后在第一时间步j计算，如在方框308示出。第一时间步j可以在仿真中对应于最终时间。在方框310，做出关于正向关联缓冲是否在时间步j存在的确定。如果关联缓冲存在，无论作为独立关联检查点或作为全状态检查点的一部分，过程流程都进展到方框312。如果在方框310确定在特别时间步j关联缓冲不存在，那么流程可以进展到在图3B中示出的步骤，从而计算与在该时间步实行关联所必需的正向传播震源激励波场关联的必需信息。
在方框312，从在存储单元(例如存储器或硬盘驱动器)中含有的已存储正向仿真状态314的数据库访问在关联缓冲中的波场。在方框316，正向仿真波场状态与伴随波场在时间步j互关联，其中结果求和进入关联图像内的关联体积。关联图像可以然后由公式生成，其中是“互关联”运算，并且ci是全震源激励波场状态si的关联信息子集，其表达为ci＝p(si)。关联算子可以表示比仅是一些应用的简单数学互关联的算子更一般的算子。重要特征是该算子需要在相同时间步j的正向仿真状态的子集ci和伴随状态作为输入。
在方框318，伴随波场可以时间反向向后传播一步，即j＝j–1。然后(n+1)个伴随状态的另一序列可以用终点和根据生成的在的n个过程调用生成，例如源自已记录轨迹的向后传播。在方框320，j可以减量，并且在方框322，可以做出关于伴随传播是否完成的确定。如否，那么过程流程返回到方框308从而在下个时间步继续。
如果在方框322伴随传播完成，那么过程流程进展到方框324，在方框324各种后处理和输出步骤可以执行。对于用来生成地震图像的实施例，这样的后处理可以包括例如从迁移点形成子表面偏移距道集。偏移距道集可以转换成反射角道集从而形成最终图像。进一步地，可以实施滤波计算从而减少噪声，例如拉普拉斯算法，如上面讨论。过程在方框326结束，其中图像的数据表示输出或存储。该数据表示可以包括结果的子集的图像。例如，在地震成像中，输出可以仅包括反射角道集。在其它实施例中，生成数据表示可以包括形成完整图像的其它步骤。例如，在地震成像中，可以执行静噪和堆叠操作从而生成完整地震图像。然而，本领域技术人员认可本技术不限于地震成像，并可以用于将正向仿真数据与以逆时顺序投影的数据互关联的任何数目的其它仿真。
图3B是根据本技术的示范实施例为逆时互关联示出波场的正向传播的过程流程图。如果在方框310中没有发现关联缓冲，那么可以执行在图3B中的步骤。在方框328，正向仿真的波动方程可以跨栅格传播。可以通过取得用s0开始的(n+1)个状态si的序列和根据si+1←fi(si)生成的在si的n个过程调用fi生成，例如正向激励脉冲的传播。
如在上面的方程中表明，接收器波场可以时间反向传播，并在每个计算单元(x)与震源波场的时间正向仿真互关联。这意味着时间正向震源仿真必须以逆时顺序访问。由于存储器和访问时间限制，为互关联存储时间点中正向的全部可以是不可行的，因此在方框330，正向仿真波场314可以在各种时间步存储，如在此讨论。可以认为这些正向仿真波场是为与逆向数据投影的互关联含有关联缓冲的检查点。数据库314中检查点中的每个可以是可以用来重启传播计算的全状态检查点，或仅是用于在特别时间步确定互关联的较小的关联检查点。在方框332，震源波场可以时间正向传播一步。在方框336，做出关于正向传播是否完成的确定。如否，那么流程进展到方框328从而继续计算。如是，那么流程返回到方框312(图3A)，互关联在方框312重新开始。其中保存全状态检查点对仅保存关联检查点的时间步的确定(例如，在方框306的检查指示过程)在下面讨论。
在本技术的示范实施例中，数个假设被做出从而改善检查指示过程，并在以下讨论中使用。特定地，该每个全状态检查点si需要大小为S的相同量的存储器。进一步地，该每个关联检查点ci需要大小为C的相同量的存储器。另外，该每个过程调用fi需要相同量的计算量l，而且没有或ci＝p(si)的计算量。在fi或的过程期间另外存储器需求不认为是存储器约束的部分。
Griewank检查指示策略假设保存全状态检查点(变量S)需要的存储器等于保存为成像步骤实行互关联必需的信息需要的存储器(变量C)。例如，如果执行基于速度-应力传播函数的逆时迁移，那么可以仅具有使压力关联的需要，代替使质点速度和应力的全部分量关联的需要。这可以使S/C的比容易地提高到4以至高达9。如果使用优选匹配层(PML)边界条件，那么S/C比可以接近地球模型的边界提高两到三倍。因此，波的传播可以在精细栅格上发生，但成像可以在较粗糙栅格上执行。对于其中关联也可以在每个空间维度中每个其它样本上完成的情况，产生的S/C比可以是八。如果在互关联步骤中的一些准确度可以通过存储压缩的关联检查点来满足，那么S/C比可以提高两倍。这里提到的S/C比分量相乘从而产生该算法需要的S/C比。因此，实际上在S和C之间可以具有充分的不同。假定逆时深度迁移非常计算上昂贵，那么本技术可以通过利用S和C之间的不同提供计算开销的充分减少。
在没有存储器约束情况下的优化策略
如果没有存储器约束，那么优化策略(在最少计算量的意义上)可以详述为：
对于i＝0，n-1
ci＝p(si)，存储ci
si+1←fi(si)
对于i＝n-1，1
在该策略中，总存储器消耗是(n–1)C，并且总正向计算量是n(假设每个伴随步骤需要与正向步骤相同的工作量，那么总正向计算量是2n-1)。可以注意到由于状态0通常对应于逆时迁移应用的零波场，因此忽略在状态零的关联。然而，在此描述的技术可以容易修改从而适应非零初始波场。
在以下讨论中，符号J(n，M)可以用来表示在存储器约束M下用于使n个状态(状态1到状态n)关联的计算量。在此情况下，可以假设伴随状态n是已知的，并且可以总是在计算期间通过将一切重设到零或从存储器加载状态0(其不是M的部分)，不用计算量重获状态0。在计算之后，使用在伴随状态1(代替在伴随状态0)的伴随状态。
进一步地，符号Jno_cp(n，M)在此用来表示在不使用任何全状态检查点的情况下(即仅存储关联检查点)，在存储器约束M下用于使n个状态关联的计算量。在其中M巨大并且具有充足存储空间从而在合适缓冲中保存关联检查点中的全部的情况下，这可以简单表达为对于M≥(n-1)C，Jno_cp(n，M)＝n。
有存储器约束的优化策略
通常，在计算的一个扫描nC＝1+M/C中，关联缓冲可以保存在存储器中，即nC是可能给出的存储器M和关联缓冲大小C的数目。这假设fi的函数调用不在存储器上放置约束，或源自正向仿真的状态缓冲结果可以用来执行关联。因此，在保存关联缓冲但没有检查指示任何正向仿真状态时需要的计算可以由在方程4中示出的公式计算。
对于M<(n-1)C，nC＝1+M/C，t＝n/nC，r＝n–tnC
在方程4中，t＝使用全部nC关联缓冲实行的扫描的数目，并且r＝剩余＝在其中仅使用关联缓冲的一部分或不使用该关联缓冲的最后一次实行的关联的数目。注意如果r等于零，那么该最后的部分次不需要完成。应理解在方程4和其它相似方程中，使用整数运算代替浮点运算。
图4是根据本技术的示范实施例图解在时间指数k使用全状态检查点的最小计算量的示意图。如果在此采用的优化化目标是在时间步k使用一个检查点402优化化J(n，M)，那么对于在时间步指数k和高于时间步指数k的关联步骤404，k个正向步可以用来为时间正向得到检查点位置加上J(n-k，M-S)，并得到时间步指数k+1到n需要的伴随步进。由于在时间步指数k的状态缓冲的检查点402，因此可用的存储器减少S。指数0到k–1的正向和伴随时间步进的计算量406是J(k-1，M)，如在方程5中示出。
Jadd_one_cp(n,M,k)＝k+J(n-k,M-S)+J(k-1,M) 方程5
进一步地，方程6允许在没有使用的检查点或在一些状态k的至少一个检查点(状态k在状态0和状态n之间)的限制下最小化J(n，M)。
如在上面公式中可见，存储器约束大小总是以S为单位变化，并且在存储器等于nSS的情况下具有整数数目的检查点状态缓冲nS。与nC个关联缓冲关联的存储器等于(nC-1)C。与两种检查点状态缓冲和关联缓冲的组合关联的总存储器在方程7中示出。
nSS+(nC-1)C≤M 方程7
在示范实施例中，上面方程迭代求解从而为希望的问题大小和存储器覆盖范围寻找数值解。
图5A-B是根据本技术的示范实施例图解时间步的示意图，全状态或关联检查点可以在该时间步存储。例如，如在图5A中示出，关联检查点502可以在位于时间步k的全状态检查点504之前存储。然而，对于在全状态检查点之前使用关联检查点的任何情况，具有在至少同样有效的全状态检查点之前不使用关联检查点的对应情况。在图5A中，在位于k的全状态检查点504之前到达nC个关联检查点508中的第一个的成本506可以表达为J＝J(k-1–nC，M)。关联检查点508的成本可以表达为J＝nC+2。最终，从全状态检查点504到达结束时间步512的成本510可以表示为J＝J(n-k，M-S–nCC)。基于这些，如果优化情况在位于一些时间步k的全状态检查点504之前存储nC个缓冲，那么总计算量可以在方程8中给出。
∑J＝k+nC+2+J(k-1–nC，M)+J(n-k，M-S–nCC) 方程8
在图5B中，到达在k–nC的全状态检查点516的成本514可以表达为J＝J(k-1-nC，M)。从全状态检查点516到达在结束时间步522之前存储的nC个关联检查点中的第一个的成本可以表示为J＝J(n-k，M-S)。因此，在其中全状态检查点516之前没有关联检查点存储的情况下，总计算量在方程9中示出。
∑J＝k+nC+2+J(k-1–nC，M)+J(n-k，M-S)。 方程9
如通过在方程8和9中的公式可见，后面的情况至少与前面的情况(其在全状态检查点之前存储关联检查点)一样良好。因此，在示范实施例中，仅需要考虑后面的情况。
优化检查指示对Griewank检查指示
全状态检查点对关联检查点的存储比(S/C)可以经常在10左右，因此，这提供方便的比率从而用于将在此描述的优化检查指示的比率与Griewank检查指示比较。在下面表1中，n提供正向和伴随计算需要的时间步的数目。标记为“理想的”的列提供在其中存储器充足并且成本可能最小的环境中使得正向和伴随仿真能够关联的时间步进的计算成本。标记为“Griewank检查指示”的列假设nS＝M/S个状态缓冲可用，并且不使用关联检查点。标记为“优化检查指示”的列示出使用全状态检查点和关联检查点获得的结果。标记为“比率”的列通过Griewank计算结果除以优化检查指示结果的比率来计算，并提供与仅使用Griewank优化检查指示策略比较的用该新方法可能的加速。
可以从表1中的结果做出两个观察结论。第一，通过利用状态缓冲和关联缓冲之间的不同，可以减少计算量。进一步地，该优化检查指示策略的实际应用需要一些预计算从而寻找优化策略。然而，预计算量与应用例如逆时迁移和梯度计算中的总计算量比较是可忽略的。
进一步地，可以注意到应用于3-D逆时迁移的普通Linux集群应用软件的通常参数可以对应于具有4000个时间步、10个状态缓冲并且M＝100的直线，其中比率(速度提高)约为1.28。如果公司具有专门为逆时反演工作的$2000万的计算机系统，那么该算法改变的经济价值在增加的生产中约为$560万。
表1：优化检查指示对Griewank检查指示的成本
图6A-C是在根据本技术的示范实施例图解具有55个时间步的简化仿真的Griewank和优化检查指示方案的操作的示意图。对于例子的每个，步数(55)和存储空间的量(20)是相同的。图6A图解对于n＝55，M＝20，S＝10的情况具有210的成本J(即，没有关联检查点)。图6B图解对于n＝55，M＝20，C＝1，S＝10的情况，带有J＝102的优化检查策略。图6C图解对于n＝55，M＝20，C＝2，S＝10的情况，带有J＝148的优化检查策略(中间)。在图6A中示出的Griewank检查指示算法中，仅保存全状态(由水平线表明)，允许传播算法从每个已保存状态完全重启。在操作中，算法执行正向仿真的多个迭代。例如，在第一轮602中，在时间步35(606)和在时间步50(604)的激励脉冲的正向传播期间保存全状态。这完全利用可用存储器M＝20。算法然后将正向传播时间步前进到时间步55，并执行与源自接收器的反向传播波场的互关联。然后反向传播波场时间反向延伸一个时间步到时间步54。正向传播波场在时间步50(64)恢复并前进到时间步54，并且再次执行互关联。该进程为时间步53、52、51重复，并且在时间步50(604)使用全状态检查点自身执行关联。
用于时间步50的存储器然后释放，并且互关联与从时间步35(606)运行的正向仿真一起进行。已释放存储器用于从时间步35(606)存储检查点，例如在时间步45(608)的检查点。在时间步45(608)的检查点用来在时间步49到46运行正向仿真，其中在每个时间步，正向传播波场与源自接收器的反向传播波场互关联。在时间步45(608)使用全状态检查点完成互关联之后，存储器释放，并用于在时间步41(610)的另一检查点。该进程重复直到已为全部时间步执行互关联。可以注意到在检查点之间的间距在图6A中是不均匀的：在第一次扫描中，例如检查点在50、45、41、38和36，而不是例如在50、45、40和36。这是因为使用优化化算法(在图6A的情况下，Griewank优化化)确定在附图中示出的数字为优化，而检查点位置的其它可能选择可以或可以不是优化的。尽管利用的存储器的量减少，但可见计算的量充分增加。由于存储器的量和存储器的访问时间经常有限，因此该进程可以充分改善运行时间，使全RTM可行。然而，在Griewank检查指示中使用的计算数和存储器的量对于复杂地震成像计算例如RTM可以仍然是有问题的。
在本技术的示范实施例中，保存两个不同类型，仿真状态可以从其重启(或互关联可以从其执行)的全状态检查点和仅保存互关联需要的波场数据的较小检查点。这可以为计算和存储器提供显著益处。例如，如果假设用来保存关联检查点的存储器是用来保存全状态检查点的存储器大小的1/10，那么可以导致在图6B中示出的优化存储器利用。在该条件下，可以需要少得多的全仿真轮。如在图6C中使用，较长直线表明全状态检查点，而较短直线表明关联检查点。在图6B中的第一轮期间，全仿真轮612从零点614开始。在全仿真轮612期间，关联检查点(或关联缓冲)在从时间步35(616)到时间步54的每个点保存，因此耗尽带有20个关联缓冲的可用存储器M＝20。由于在正向仿真数据自身中含有关联信息，因此不需要在时间步55存储该信息。可以然后在这些点中的每个执行与源自接收器的反向传播波场的互关联。一旦为时间步55到35完成互关联，那么从零时间步执行全仿真的第二轮620。在全仿真的第二轮620期间，在从时间步14(622)到时间步33(624)的每个时间步保存关联检查点，再次耗尽M＝20。在这些时间步中的每个使用在该时间步的关联检查点重复互关联，并然后为全仿真的第三轮626重复该进程。如可见，对于相等资源，优化检查指示策略提供计算时间的充分节省，其中与Griewank检查指示进程的209个正向波传播步相比，优化只是策略仅执行101个正向波传播步。
在图6B的概述中，可以通过尝试满足方程7的nS和nC的全部组合，并选择使用方程6为J给出最小值的一个(或若干个中的一个)来确定的对于C＝1的本发明的优化策略证明通过完全不使用全状态缓冲，而是使用完全用于关联缓冲的可用存储器M＝20被满足。如接下来在图6C的讨论中所见，对于C的更大值(C＝2)的优化策略包括同时使用一个全状态缓冲和五个关联缓冲。
即使用于关联检查点的存储空间大小更接近全状态检查点的存储空间，但优化检查指示进程仍可以提供充分的节省。例如，图6C示出以与6A和6B中相同的存储器量，但用全状态检查点和关联检查点之间大小的5比1的关系运行的优化检查指示。该进程与图6B中一般相同，但全状态检查点628在第一全状态仿真轮630期间保存。在第一轮全状态仿真期间，关联检查点在时间步49(632)和时间步54(634)之间保存。互关联可以然后在时间步55，而且用时间步54和49之间的值的每个执行。该全仿真轮从时间步31(628)重复三次或更多，如由参考号636表明，并且在每个时间步运行互关联。存储器是空闲的，并且三个最终全仿真轮638从零时间步执行，其中互关联为已保存关联检查点中的每个执行。如从图6C可见，对于该状况正向波传播步仅为147，这仍显著小于Griewank检查指示的209个正向波传播步。
优化检查指示技术可以为RTM提供超过较早的Griewank检查指示策略百分之10到超过百分之400之间的加速，取决于给定RTM应用的计算机硬件环境和问题大小。进一步地，在此讨论的检查指示策略不限于直接关联，但可以用于可以关于一些FWI梯度计算和一些自动微分公式需要的其中当前状态取决于过去状态的关联。
图7A-C是根据本技术的示范实施例将标准关联的检查指示与在伴随计算中过去状态影响与正向计算中当前状态的互关联时的检查指示比较的示意图。图7A图解标准关联条件，而图7B图解除例如在运行其中过去条件影响当前条件的储层仿真时的移位关联条件之外的标准关联条件。图7C图解优化检查指示策略对具有12个时间步、大小为6的存储器，并且其中C＝2和S＝3情况的应用。在图7A-C中，向上箭头702表示正向仿真，并且向下箭头704表示伴随仿真，例如其中在先前时间步计算的值影响在稍后时间步计算的值。水平双向箭头706表示标准关联条件，并且对角双向箭头708表示移位关联条件。进一步地，粗水平线710表示大小为S的全状态检查点，并且虚线712表示大小为C的关联检查点。
对于许多应用(例如，逆时迁移或梯度计算算法中的一些)，仅需要使对应状态关联(即，如在图7A中示出)。然而，在一些情况下(例如对于在储层仿真中使用的一些梯度计算算法)，我们不仅需要使对应状态关联，但也需要使邻近状态关联(即，如在图7B中示出)。在本技术的示范实施例中，修改优化检查指示策略从而为后面的情况工作，如在图7C中表明。
以在图7C中描述的方式实施优化检查指示策略提供最小存储器使用。然而，使用三种类型的关联，即，使伴随波场与正向波场关联、使伴随波场与存储中的快照缓冲(snapshot buffer)关联，以及使伴随波场与存储中的关联缓冲关联。这可以通过将互关联问题视作逆向链问题来更优理解。本技术的实施例可以提供将链式计算逆向的更有效方式，例如逆向链，该链式计算使源自正向仿真的结果与源自反向或伴随仿真的结果关联。其可以应用于如在图7A中图解的标准互关联。其也可以应用于如在图7B中图解并在图7C中更详细图解的已修改互关联。其它互关联条件可以从本技术受益，只要实施基于使用逆向链。
优化检查指示到带有不同访问速度的二级存储的延伸
除了全状态检查点或关联检查点的存储之外，本技术的示范实施例也可以基于可用存储的速度使检查点的存储优化化。在本计算架构中，普遍具有带有不同访问速度的多级存储。例如，计算可以在使用GPU存储器(最快存储)的图形处理单元(GPU)上执行，或在CPU存储器、固态硬盘存储器和磁盘存储器(最慢存储)中执行。为简化该概念的解释，可以假设二级存储，例如CPU存储器和磁盘。
在确定检查点的优化分布中，可以做出两个假设。第一，可以假设大小MF的快速存储需要分别写入全状态检查点和关联检查点的的成本，以及分别读取全状态检查点和关联检查点的 的成本。进一步地，可以假设慢速存储(大小MS)需要分别写入全状态检查点和关联检查点的的成本，以及分别读取全状态检查点和关联检查点的的成本。另外，可以假设将任何状态发展到下个状态的计算量是一，全状态检查点的大小是S，并且关联检查点的大小是C。
与前面相似，成本函数可以如在方程10中示出来定义。在方程左侧上的成本函数是在右侧上列出的三个选择的最小值。
在方程10中，*可以用F或S代替，从而表明起始检查点是否在快速或慢速存储中。因此，对于其中第一状态在快速存储中的情况(即，如果第一状态需要恢复以便重计算，那么该信息在快速存储中)，JF(n，MF，MS)是使n个状态与MF快速存储和MS慢速存储互关联的成本。项JS(n，MF，MS)是相似的，除了第一状态在慢速存储中之外。JF no_cp(n，MF，MS)和JS no_cp(n，MF，MS)遵循带有不允许全状态检查点并因此仅存储关联检查点的约束的相同定义。另外，可以注意因为尽管全状态检查点在存储中，但如果目的不是重启仿真那么仅读取全状态检查点的关联缓冲部分，所以或用于第二和第三行中的第四项，代替或然而，这可以取决于恢复关联缓冲或已存储全状态检查点的实际实施成本。
成本函数J*no_cp(n，MF，MS)的计算比在一个存储的情况下使用Jno_cp(n，M)的显式形式更复杂。然而，其可以使用相似的递归技术执行，如在方程11中示出。
在方程11中，i和j是我们要使用的关联缓冲的数目。通过不表达n-1-i-j≥0的事实来简化该公式。此外，对于n-1-i-j＝0的情况，在最小值方程中的第二行不再计算。如在下面进一步讨论，可以明确估算J*no_cp(n，MF，MS)，代替依靠上面的递归公式。
确定在具有两个不同访问速度的存储中优化检查点的结果在图8A-图8E和图9A-图9E中示出。例如在图8A-图8E和9A-图9E中示出，为简化解释，假设访问快速存储器的成本为零(即， )。然而，在此描述的技术为具有到快速存储的非零访问成本的任何实施工作。进一步地，长水平线表示全状态检查点，短水平线表示关联检查点，并且在末端具有菱形的任何直线存储在慢速存储器中。
图8A-E是根据本技术的示范实施例图解数个状况的示意图，其中快速和慢速存储器的量变化，同时每个类型的检查点需要的存储空间保持相同。在图8A-E中，假设慢速存储器的速度为 假设每种类型的检查点使用的存储空间的量为C＝S＝2(匹配Griewank假设)。在图8A中，存储器的量是MF＝3并且MS＝0。在8B中，MF＝2并且MS＝3。在8C中，MF＝2并且MS＝12。在8D中，MF＝14并且MS＝0。最终，在8E中，MF＝3并且MS＝20。如在图8A-E中可见，存在的存储的量和类型可以对优化化的结果具有充分效果。例如，在图8C和8D中，总存储的量相同，但存储速度的不同导致所获得优化检查指示策略的不同。
一般地，访问存储的成本(即，I/O的时间)基于信息的大小是不同的。例如，如在此讨论，在S充分大于C的情况下，用于全状态检查点的存储可以充分大于用于关联检查点的存储。
图9A-E是根据本技术的示范实施例图解数个状况的示意图，其中快速和慢速存储器的量变化，其中用于关联检查点的存储空间是二并且用于全状态检查点的存储空间是三。因此，在图9A-E中， C＝2并且S＝3。在图9A中，存储器的量是MF＝3并且MS＝0。在9B中，MF＝2并且MS＝3。在9C中，MF＝2并且MS＝12。在9D中，MF＝14并且MS＝0。最终，在9E中，MF＝3并且MS＝20。可以注意在图9A-图9E的每个中使用的存储器的量对应于在图8A到E的每个中使用的存储器的量，示出关于全状态检查点对关联检查点的存储器利用不同时发生的策略的改变。
图8E和9E的比较图解使用慢速存储的后果可以完全改变策略，即使总存储对于两个情况显然足够大。因此，用检查指示非常大的磁盘来执行RTM可以慢于向并行计算机上小得多的数目的分布式RAM缓冲相同地使用检查指示。
表2提供在此披露的优化检查点技术与Griewank检查点策略之间求解时间(包括正向和伴随计算)的比较。在该例子中，1000个状态可以互关联，C＝1Gb并且S＝10Gb。进一步地，示出存储的不同组合，并且假设访问成本为和
硬件存储 检查指示策略 求解时间
20G存储器 Griewank 3.66小时
20G存储器 优化 2.75小时
200G磁盘 Griewank 3.47小时
200G磁盘 优化 1.89小时
20G存储器和200G磁盘 优化 1.07小时
200G存储器 Griewank 1.04小时
200G存储器 优化 47分钟
无限存储器 33分钟
表2：在优化检查指示策略和Griewank策略之间的求解成本比较，其中C＝1Gb，S＝10Gb，并且n＝1000。
成本的显式估算
代替使用递归公式，可以明确估算J*nocp(n，MF，MS)。为达到该解，考虑J*nocp(n，m，MF，MS)，其为在*存储(F或S)和MF、MS的存储器约束中，通过不使用全状态检查点的m个正向扫描使n个时间步与状态0互关联的的互关联的成本。因此，在方程12中示出的进程可以用来以存储全状态检查点来确定成本。
最小值可以然后从在方程13中示出的公式确定，如在方程13中示出。
为最小化J*nocp(n，m，MF，MS)，可以考虑三个情况。第一，可以完全不使用存储，其中m＝n。第二，可以仅使用快速存储，这在方程14中提供公式。
在第三情况下，快速存储和慢速存储都使用，并且在此情况下，通过使用尽可能多的快速存储来改善速度(否则不需要使用慢速存储)。该情况在方程15中提供公式。
基于这些情况，可以围绕在方程16中示出的5个临界值搜索m。
方程16
可以用来得到优化m和对应J*nocp(n，MF，MS)。在应用这些上面公式时应考虑合适边界。例如，如果不在和n之间，那么不需要考虑它。
优化检查指示策略的有效计算
尽管J*no_cp(n，MF，MS)可以明显确定，但递归算法仍用来寻找J*(n，MF，MS)。使用递归公式，使用nF个快速存储缓冲(nF＝MF/S)和nS个慢速存储缓冲(nS＝MS/S)使n个状态互关联的成本可以用计算量～n2×nF×nS预确定。这里符号“～”用来表示“大约”。该预制表数据需要大小为～n2×nF×nS的存储。一旦具有预制成本数据，那么仅需要工作量～n2从而为两种类型存储的给定容量约束计算优化二级存储检查指示策略。最计算便宜的时间步是用工作量～n2×nF×nS预计算成本表J*(n，MF，MS)(用于使高达n个状态与存储约束MF、MS关联)。尽管仅是一次计算量，但具有使该时间步有效的价值。在示范实施例中，这可以使用k的另外制表(第一全状态检查点位置)来执行。检查指示策略的实际确定采取～n的工作量代替～n2。在上面讨论的技术和公式可以用来为检查点确定优化数目和位置，如关于在图10-14中示出的方法所讨论。
图10是根据本技术的示范实施例的用于实施使用关联检查点的仿真的方法的过程流程图。过程在方框1002用优化检查指示策略的计算来开始。该计算可以使用上面讨论的公式执行。在方框1004，分配需要的存储。这可以受可用存储、全状态检查点和关联检查点的大小等限制。在方框1006，实施检查指示策略和关联计算，例如，循环通过一系列检查指示指令。
例如，检查指示策略可以用指令的排序列表的形式表达。下面表3提供可以用于实施检查指示策略的指令集的例子。在该例子中，“指令号”与在遇到该“指令号”时完成的一系列操作关联。可以具有定义这样的指令集的多种方式，但具有一个版本例如下面给出的版本足以使本技术成为可能。
表3：用于实施检查指示策略的示范指令
实施检查指示进程的指令的使用例子在下面呈现。在此情况下，关联具有266个时间步，其中S＝12，C＝1，MF＝10，MS＝40，慢速存储的访问成本Rs＝Ws＝2且Rc＝Wc＝0.6，并且快速存储的访问成本是零。
(0 74 S)，(0 71 S)，(0 66 S)，(2 55 10F 4S)，(3 S)，(2 40 10F 4S)，(3 S)，(225 10F 3S)，(3 S)，(2 11 10F OS)，(1 S S)，(2 65 10F 16S)，(3 S)，(2 38 10F 16S)，(3S)，(2 11 10F 0S)，(1 S S)，(2 70 10F 28S)，(3S，2 31 10F 9S)，(3 S)，(11 10F 0S)，(1S F)，(2 73 10F 40S)，(3 F)，(2 22 10F OS)，(3 F)，(2 11 10F OS)，
附录A和B是可以用来生成优化检查指示策略，例如通过上面指令列表，仅给定上面指定的参数，即S、C等来实施的优化检查指示策略的计算机程序的两个例子的源代码。附录A和B的程序使用本发明的检查指示策略；然而，在Griewank检查指示优化的情况下，该程序选择Griewank检查指示。
图11是根据本技术的示范实施例的用于在没有任何全状态检查点的情况下计算最小成本的方法的过程流程图。该方法在方框1102用空列表的初始化从而存储扫描的潜在优化数目来开始。在方框1104，该列表的第一个值初始化从而等于n，在正向仿真中时间步的数目。
在方框1106，通过公式和来计算可以保存在快速存储和慢速存储中的关联缓冲的最大数目。在方框1108，做出关于至少一个快速关联检查点缓冲是否可用，即的确定。如果没有快速检查点缓冲可用，那么过程流程进展到方框1118。
如果至少一个快速关联缓冲可用，那么在方框1110，使用公式计算在仅使用快速关联缓冲的情况下扫描的最小数目。在方框1112，使用公式计算在仅使用快速关联缓冲的情况下扫描的优化数目。在方框1114，如果并且仅如果mF<mF opt<n，那么mF opt的值添加到列表。在方框1116，mF添加到列表。
在方框1118，做出关于至少一个慢速关联检查点缓冲是否可用，即的确定。如果没有慢速关联检查点缓冲可用，那么过程流程进展到方框1128。
如果至少一个慢速关联缓冲可用，那么在方框1120，使用公式计算在使用全部类型的关联缓冲的情况下需要的扫描最小数目。在方框1122，使用公式计算在使用快速关联缓冲的全部和该至少一个关联缓冲的情况下扫描的优化数目。在方框1124，如果并且仅如果mF&S<<列表中的最后值，那么添加到列表。在方框1126，mF&S添加到列表。
在方框1128，列表重生成从而在相同范围内包括最近整数。例如，如果原列表是(17，10，1)，那么新列表应是(17，16，11，10，9，2，1)。这可以帮助解决实数到整数的舍入。
在方框1130，使用以下公式计算新列表中每个潜在优化扫描数目的对应成本。
其中mF&S用来在使用前述公式之前定义MF和MS。在方框1132，返回最小成本。
图12是根据本技术的示范实施例的计算用于使n个状态互关联的最小成本和优化位置的方法的过程流程图。该方法假设nf_used快速存储器检查点状态缓冲和ns_used慢速存储器检查点状态缓冲已经使用，而且起始全状态检查点存储在StartType存储中，并且第一全状态检查点存储在KType存储中(斜体术语在本发明方法的C和C++实施中使用，并在下面公式中再次出现)。StartType表明与正向仿真从其发动的全状态检查点关联的存储类型，快速或慢速。KType表明与在时间指数k存储的全状态检查点关联的存储类型，快速或慢速。“StartType”可以取两个值，对于快速存储器取F或对于慢速存储器取S。“nf_used”是当前使用的快速存储器检查点状态缓冲的数目。用于状态或关联信息的新检查点的当前可用快速存储器不包括用于当前使用的快速存储器检查点状态缓冲的空间。同样，“ns_used”是当前使用的慢速存储器检查点状态缓冲的数目，并且当前可用慢速存储器不包括用于当前使用的慢速存储器检查点状态缓冲的空间。
在方框1202，做出关于状况是否可能的确定。如果该状况不可能，那么过程流程进展到方框1204，其中返回巨大成本(例如1e10)从而鉴别其是不可能的。
在方框1206，做出关于是否n＝2的确定。如是，那么优化第一全状态检查点位置必须是1，并且对应成本在方框1208返回。
在方框1210，从预计算表获得先前优化k。计算在表中的该值以便使n-1个状态与相同存储器约束以及相同起始全状态检查点和第一全状态检查点存储类型需求互关联。
在方框1212，为离开先前优化k的某个距离范围例如S/C+1内的每个k，计算放置在位置k的第一全状态检查点的对应成本。如果第一全状态检查点需要在快速存储中(即KTYPE＝F)，那么成本由下面公式给出：
如果第一全状态检查点不需要在快速存储中，那么成本由下面公式给出：
在上面方程中，以下公式：
已用来表明存储以S为单位减少。在方框1214，预计算表用第一检查点的已计算最小成本和优化位置更新。
图13是根据本技术的示范实施例示出怎样预计算与在优化k位置取得的全状态检查点关联的最小成本的表格的过程流程图。该方法在方框1302用分配数组从而存储表格开始。例如，两个五维数组：J_with_cp[2][2][n][nfMax][nsMax]和BestK[2][2][n][nfMax][nsMax]可以分配从而保存第一检查点的最小成本和优化位置的轨迹。作为分配的部分，在数组中的全部值初始化为-1。在两个数组中，第一维表示起始检查点的存储类型。第二维表示第一检查点的存储类型。第三维表示有待关联的状态的数目。第四维表示在快速存储中检查点的最大数目。第五维表示在慢速存储中检查点的最大数目。
在方框1304，追踪计数器以便从2到n循环通过I，其中I表示有待关联的状态的数目。在方框1306，追踪另一计数器以便从0到nfMax循环通过J，其中nfMax表示已使用的快速存储器状态缓冲的数目。在方框1308，追踪第三计数器以便从0到nsMax循环通过K，其中nsMax表示已使用的慢速存储器状态缓冲的数目。
在方框1310，计算第一检查点的最小成本和优化位置，以便使I个状态和已经用于4个情况的J个快速存储器全状态检查点缓冲与K个慢速存储器全状态检查点缓冲互关联：1)用在快速存储中的全状态检查点开始，为在慢速存储中的第一全状态检查点计算成本。2)用在快速存储中的全状态检查点开始，为在慢速存储中的第一全状态检查点计算成本。3)用在慢速存储中的全状态检查点开始，为在快速存储中的第一全状态检查点计算成本。4)最终，用在慢速存储中的全状态检查点开始，为在慢速存储中的第一全状态检查点计算成本。获得的结果用来为J_with_cp和BestK更新对应表格条目。方框1304的这四个步骤中的每个可以使用图12的方法计算。成本计算1212需要指定存储类型的全状态存储的成本、在指定时间步为给定全状态检查点完成多少次恢复、以及指定存储类型的全状态恢复的成本的认识。在nS≥2时，上面概述的过程可以递归应用从而优化定位另外的全状态检查点。
图14是根据本技术的示范实施例的用于计算最小成本和优化检查指示策略，以便在起始检查点在StartType存储中的状况下使n个状态与nf_used快速存储器状态缓冲和ns_used慢速存储器状态缓冲关联的过程流程图。该方法在方框1402开始，其中与快速存储中第一全状态检查点关联的成本从表格J_with_cp获得，并且全状态检查点的对应优化位置从表格BestK获得。在方框1404，与慢速存储中第一全状态检查点关联的成本从表格J_with_cp获得，并且全状态检查点的对应优化位置从表格BestK获得。由图13的方法生成的表格可以在方框1402和方框1404使用。在方框1406，确定在没有使用任何检查点和对应优化扫描数目的情况下的成本(见于图11)。在方框1408，根据哪个成本(源自步骤1、2或3)最小，优化检查指示策略递归输出。
在检查指示计算中最小成本的鉴别
图15是根据本技术的示范实施例示出在S/C＝12，并且Rs＝Ws＝Rc＝Wc＝0时使用四个缓冲的266个时间步的互关联的第一全状态检查点(k)的最小成本和位置之间关系的图表。在S/C＝1，并且Rs＝Ws＝Rc＝Wc＝0的特殊情况下，优化第一检查点的位置具有显式形式(如由Griewank检查指示策略确定)。然而，一旦具有S>C，那么优化第一检查点的位置变得复杂，即使在Rs＝Ws＝Rc＝Wc＝0仍然保持的简单情况下。在图15中示出的例子中，对于第一检查点的选择具有三个位置最小值。当Rs、Ws、Rc、Wc不再是0时，行为变得更非线性，如在图16中示出。
图16是根据本技术的示范实施例示出在S/C＝12，Rs＝Ws＝12，并且Rc＝Wc＝1时使用四个缓冲的266个时间步的互关联的第一全状态检查点的最小成本和位置之间关系的图表。在该例子中，假设从正向仿真状态零恢复不需要计算量来生成图表。如从图16中的图表可见，可以使最小值1602的鉴别复杂。源自用递归公式的众多试验的一个重要观察是当(1)保持可用存储相同(即MF和MS相同)和(2)起始全状态检查点类型和第一全状态检查点类型相同(即在快速或慢速存储中)时，在增加n时第二局部最小值一般变为全局最小值。另一重要观察是在两个局部最小值之间的距离约在离开～S/C+1的距离内。由于这两个观察，可以在小数目的可能性之间搜索第一检查点位置(k)。这导致与～n2×nF×nS比较仅用～n×nF×nS的计算量预计算成本表的算法。
仅接近预期检查点位置搜索用策略质量的一些牺牲显著减少策略计算成本～n倍。然而，本技术可以提供带有离开从全搜索获得的优化指示策略小于1％的另外成本的检查指示策略。
对于非常大的n，用于cost和BestK表格的存储器和表格计算时间可以是巨大的。用于非常大的n的实际实施策略可以是为nmax_table的合理大小计算优化表格，并然后切换到Griewank全状态检查指示策略，从而如果n>nmax_table那么寻找BestK，并且如果n<nmax_table那么使用表格。可以在不需要表格的情况下估计BestK的Griewank值。通常这产生仅比用全尺寸表格寻找的优化值昂贵几个百分比的检查指示策略。
从优化检查指示策略获得的性能提升
图17是根据本技术的示范实施例示出优化检查指示策略超过使用Griewank进程的检查指示的性能改善的图表。在性能和可用存储器之间的关系通常具有在图17中示出的形状。在图表17中，顶实线1702是为使用二级存储，例如快速存储和慢速存储计算的优化性能。中实线1704是为使用单级存储计算的优化性能。下虚线1704是为使用Griewank检查指示计算的优化性能。如在图表中示出，由于需要全状态检查点的较少重计算，因此性能随着更多的可用存储(通常在添加更多计算节点时)初始提高。
然而，在可用存储超过某点(例如，点A 1708或点F 1710)之后，放大损失变得显著，并且性能下降。放大损失是向问题添加更多计算单元的结果，其中由于通过使用多个单元导致的增加开销，例如在单元之间或单元和存储系统之间的通信，因此每个计算单元的效率降低。放大损失可以是由于缺少应用的理想可放大性。进一步地，源自重计算减少的益处可以变得临界。源自优化检查指示策略的性能提升取决于可用存储。例如，如果充足存储可用，那么优化检查指示策略超过Griewank检查指示的性能提升可以是在点A 1708和点F1710之间的比率。如果存储受约束(通常是该情况)，那么性能提升可以是在点A 1712和点D1714之间的比率。
应用
本技术的优化检查指示策略可以用来为任何数目的计算密集的大容量存储过程，例如使用表格驱动的软件实施获得增强效率。一般地，从本技术受益的过程包括以逆时顺序访问从而与反向时间步进计算关联的时间步进正向仿真。
例如，在示范实施例中，优化检查指示进程在地震成像应用中使用。在这些实施例中，优化检查指示进程可以为叠前发射-记录RTM、叠前仿真震源RTM、叠前发射-记录FWI梯度、叠前发射-记录FWI Hessian、叠前仿真震源FWI梯度、叠前仿真震源Hessian以及相似地震计算提供改善的计算效率。在这些应用的每个中的计算机仿真器可以基于地震波动方程，例如各向异性粘-弹物理过程。已记录数据可以包括质点速度、压力、应力、加速度，或其任何组合。用于逆时优化检查指示策略的一些表格可以预计算并保存。在这些应用的每个中，正向时间步进地震仿真可以应用范围从恒定密度声仿真、可变密度声仿真、粘性-各向异性声仿真、各向同性弹性仿真、各向异性弹性仿真直到全各向异性-粘-弹仿真的不同水平的物理过程。正向仿真与从测量数据生成的伴随仿真状态的关联可以用来生成地震图像，从而尤其纠正速度模型、估计刚度系数、估计震源子波，以及估计衰减质量因数。
可以注意关联时间步仅需要震源激励波场和接收器时间反向传播波场为全部时间步同时互关联并求和。因此，传播或仿真可以用逆向时间步顺序执行，或用正向时间步顺序执行，并且结果相同。在上面讨论的RTM例子具有以在逆时顺序访问的正时顺序自然可访问的震源激励波场，同时接收器波场时间后向行进。可替换地，优化检查指示策略也可以应用于实施，其中接收器波场在时间后向传播期间检查指示，并然后以正时顺序访问，从而匹配时间正向传播的震源激励波场。两个实施都由在此描述的方法覆盖。
除这些应用之外，本技术可以用于任何数目的其它应用。例如，在实施例中，优化检查指示策略可以用于烃储层仿真的历史匹配。在该应用中，从储层收集的数据可以尤其包括从阱产生的流体体积和类型、底孔温度以及地层岩石记录。计算机仿真可以基于流体流动的达西定律的多相表达式。关联可以用来调整参数尤其例如孔隙度、渗透度、流体饱和度和压力分布。关联也可以用来为储层并为储层管理任务生成生产预测。这样的任务可以尤其包括为新阱确定位置，并确定是否使生产阱转换成注入阱。
在另一实施例中，本技术可以用于盆地建模。在该应用中，地层岩石记录可以是测量参数中的一个。计算机仿真基于盆地的构造变形和形成该盆地的沉积过程的模型。源自关联的输出可以帮助确定沉积速率和可用容纳空间作为位置和地质时间的函数。进一步地，输出可以包括随地质时间推移的盆地历史形状，以及目标储层和源岩的埋藏深度历史和热历史。
本技术的另一实施例可以在确定重油储层，例如在加拿大冷湖中的沥青储层中的热传递中使用。在该应用中，测量数据可以包括温度测量值、岩石类型和烃类型测量值、热膨胀、与烃生产关联的沉降等。计算机仿真可以尤其基于热传递的扩散方程，与岩石和重油刚度方程的方程耦合作为温度的函数，与重油关联系耦合作为温度和压力的函数。关联的结果可以用来尤其确定在覆盖层中的导热率和地下温度，以及储层层段和烃成分的相位变化。
本技术不限于烃系统。例如，本技术的实施例可以用来确定通过或跨过物体的热传递。例如，在航天器的大气层再入中，本技术可以用来将在航天器表面的温度与从正向有限元分析生成的仿真结果比较。本技术也可以使其它不同应用如气候仿真、核武器爆炸的仿真等受益。
集群计算系统的例子
图18是根据本技术的示范实施例的可以用来实施本技术的示范集群计算系统的框图。图解的集群计算系统具有四个计算单元1802，其每个都可以为仿真成像的部分执行计算。然而，本领域技术人员认可由于可以选择任何数目的计算配置，因此本技术不限于该配置。例如，较小仿真如地震成像计算可以在单个计算单元1802例如工作站上运行，而巨大仿真可以在具有10个、100个、1000个以至更多计算单元1802的集群计算系统上运行。
集群计算系统可以从一个或更多客户端1804经由网络1806访问，例如通过高速网络接口1808访问。网络1806可以包括局域网(LAN)、广域网(WAN)、互联网或其任何组合。客户端系统1804中的每个可以为操作代码和程序的存储具有非暂时性计算机可读存储器1810，包括随机访问存储器(RAM)和只读存储器(ROM)。操作代码和程序可以包括用来实施在此讨论的方法中的全部或任何部分的代码，例如关于图3A、图3B、图10、图11A、图11B和图12-图14讨论。进一步地，非暂时性计算机可读媒体可以保存全状态检查点、关联检查点和仿真结果，例如地下空间的数据表示。客户端系统1804也可以具有其它非暂时性计算机可读媒体，例如存储系统1812。存储系统1812可以包括一个或更多硬盘驱动器、一个或更多光盘驱动器、一个或更多闪存驱动器、这些单元的任何组合，或任何其它合适存储设备。存储系统1812可以用于检查点、代码、模型、数据和用于实施在此描述的方法的其它信息的存储。例如，数据存储系统可以为优化检查指示策略保存检查点。
高速网络接口1808可以在集群计算系统中耦合到一条或更多通信总线，例如通信总线1814。通信总线1814可以用来将指令和数据从高速网络接口1808通信到集群存储系统1816，并通信到集群计算系统中计算单元1802的每个。通信总线1814也可以用于在计算单元1802和存储阵列1816之间通信。除了通信总线1814之外，高速总线1818可以存在从而提高计算单元1802和/或集群存储系统1816之间的通信速率。
集群存储系统1816可以具有一个或更多有形的、计算机可读的媒体设备例如存储阵列1820，以便存储检查点、数据、视觉表示、结果、代码或其它信息，例如关于图3A、图3B、图10、图11A、图11B和图12-图14的方法的实施和源自该方法的结果的信息。存储阵列1820可以包括硬盘驱动器、光盘驱动器、闪存驱动器、全息存储阵列或任何其它合适设备的任何组合。
计算单元1802中的每个可以具有处理器1822和关联的本地有形计算机可读媒体，例如存储器1824和存储1826。处理器1822中的每个可以具有多核单元，例如多核CPU或GPU。存储器1824可以包括ROM和/或RAM，该ROM和/或RAM用来存储代码，例如用来引导处理器1822实施关于图3A、图3B、图10、图11A、图11B和图12-图14讨论的方法的代码。存储1826可以包括一个或更多硬盘驱动器、一个或更多光盘驱动器、一个或更多闪存驱动器，或其任何组合。存储1826可以用来为检查点、中间结果、数据、图像、或包括用来实施图3A、图3B、图10、图11A、图11B和图12-图14的方法的与操作关联的代码提供存储。
本技术不限于在图18中图解的架构或单元配置。例如，任何合适的基于处理器的设备可以用来实施本技术的实施例的全部或一部分，无限制包括个人计算机、网络个人计算机、膝上计算机、计算机工作站、GPU、移动设备，以及带有(或没有)共享存储器的多处理器服务器或工作站。此外，实施例可以在专用集成电路(ASIC)或超大规模集成(VLSI)电路上实施。事实上，本领域技术人员可以利用能够根据实施例执行逻辑操作的任何数目的合适结构。
例子
在背景节中讨论的勘探地球物理学家协会先进建模(SEAM)项目例子可以为不同计算机配置使用Griewank“优化”检查指示和使用本发明工作。
保存正向仿真的全部时间步的直接的即完全没有检查指示的方法可以意味着保存84吉字节的40，000个时间步，或约3.36拍它字节的数据。在SEAM数据上实行RTM的例子表明在RAM的大小与应用问题的大小比较非常有限时，本发明向RAM提供超过Griewank检查指示的最大优点。
这些SEAM RTM应用例子为“快速存储器”使用RAM，并为“慢速存储器”使用磁盘，并假设从磁盘读取重启状态缓冲或写入重启状态缓冲到磁盘是五个成本单位，从磁盘读取关联缓冲或写入关联缓冲到磁盘的成本是两个单位，并且计算正向时间步的成本是一个单位，以及实行时间后向步的成本是一个单位。使用RAM读取或写入重启状态缓冲或关联缓冲的成本取为零。同样，假设S/C比为12。“正向成本”是以逆时顺序提供正向仿真信息需要的成本单位的数目，并包括正向时间步进的成本与实行存储和恢复操作的成本。“总成本”是实行正向仿真逆时步进和时间后向步进从而实施RTM需要的成本单位的数目。
以下信息为每个例子制表。使用的算法需要为每个类型的存储器并为每个类型的操作计算存储和恢复操作的数目。该信息用来计算正向仿真逆时的成本。每当C大于或等于S时，使用的算法提供Griewank解。在此情况下，由可以仅通过存储实行关联需要的信息来完成的Griewank完成的全状态检查点由算法报告，如同该全状态检查点是关联信息的存储。
标记为A、B和C的三个例子比较(1)向RAM使用Griewank检查指示的逆时策略的计算机成本，(2)仅向RAM使用检查指示的本发明，以及(3)向RAM和磁盘使用检查指示的本发明。这些例子假设n等于40，000个时间步，以及可以在RAM中支持100个重启状态缓冲并在磁盘中支持100个重启状态缓冲的计算机环境。本发明可以提供超过其中全部Griewank检查点进入RAM的Griewank检查指示策略1.30x的加速。
例子A：在存储器中的Griewank
S＝84G
快速存储器＝8316G(99个状态缓冲)
正向成本＝114750
总成本＝154646
nreadc[F]＝39248 nwritec[F]＝38540 nreads[F]＝751 nwrites[F]＝708nreadc[S]＝0 nwritec[S]＝0 nreads[S]＝0 nwrites[S]＝0
例子B：在存储器中的优化
S＝84G
C＝7G
快速存储器＝8316G(用作状态或关联缓冲)
正向成本＝79262
总成本＝119261
nreadc[F]＝39957 nwritec[F]＝39916 nreads[F]＝42 nwrites[F]＝41nreadc[S]＝0 nwritec[S]＝0 nreads[S]＝0 nwrites[S]＝0
关联缓冲的使用为正向仿真逆时提供1.44x加速，并为实行RTM提供超过使用“在存储器中的Griewank”1.29x的总加速。
例子C：在存储器和磁盘中的优化
S＝84G
C＝7G
快速存储器＝8316G(用作状态或关联缓冲)
慢速存储器＝8400G(用作状态或关联缓冲)
正向成本＝79163
总成本＝119162
nreadc[F]＝39957 nwritec[F]＝39934 nreads[F]＝1 nwrites[F]＝0 nreadc[S]＝32 nwritec[S]＝0 nreads[S]＝32 nwrites[S]＝32
尽管可用的磁盘空间可以保存100个重启状态缓冲或1200个关联缓冲，但本发明的算法选择仅使用32个重启状态缓冲。向磁盘实行输入和输出操作的高成本阻碍算法使用可用磁盘空间中的全部。在此情况下加速的大部分来自S/C比，并且增加的加速仅从具有可用磁盘获得。例子C为实行RTM提供超过使用“在存储器中的Griewank”1.30x的总加速。
接下来的三个例子假设n等于40，000个时间步，以及可以在RAM中支持20个重启状态缓冲并在磁盘中支持100个重启状态缓冲的计算机环境。注意可用RAM小于上面例子中的可用RAM，并且与本发明关联的优点更大。在例子G中应用的本发明可以提供超过例子D中其中全部Griewank检查点进入RAM的Griewank检查指示策略1.67x的加速，并提供超过例子E中其中全部Griewank检查点进入磁盘的Griewank检查指示策略的甚至更高的加速。
例子D：在存储器中的Griewank
S＝84G
快速存储器＝1680G(20个重启状态缓冲)
正向成本＝185055
总成本＝225054
nreadc[F]＝33275 nwritec[F]＝27439 nreads[F]＝6724 nwrites[F]＝5836nreadc[S]＝0 nwritec[S]＝0 nreads[S]＝0 nwrites[S]＝0
例子E：在磁盘上的Griewank
S＝84G
慢速存储器＝8400G(100个重启状态缓冲)
正向成本＝513685
总成本＝553684
在磁盘中的Griewank：S＝84G，C＝84G，n＝40，000，M[F]＝0，M[S]＝0(100个缓冲)
正向成本＝513685＝114750+(39248+38540+751+708)*5(源自从先前情况添加I/O成本的计算)
注意由于对磁盘的输入和输出操作，因此在该例子中对磁盘的Griewank检查指示比对RAM的Griewank检查指示更昂贵。
例子F：在存储器和磁盘中的优化
S＝84G
C＝7G
快速存储器＝1680G(用作状态或关联缓冲)
慢速存储器＝8400G(用作状态或关联缓冲)
正向成本＝94655
总成本＝134654
nreadc[F]＝39710 nwritec[F]＝39692 nreads[F]＝33 nwrites[F]＝18nreadc[S]＝100 nwritec[S]＝0 nreads[S]＝156 nwrites[S]＝100该解决方案保存100个重启状态指示点到磁盘，并从磁盘实行重启检查点信息的156个恢复。另外18个重启状态检查点在RAM中完成。总共33个恢复操作在RAM中的重启状态检查点上完成。没有状态缓冲保存到磁盘。
接下来的两个例子假设可以仅在RAM中支持10个重启状态缓冲并在磁盘上支持50个重启状态缓冲的计算机环境。在该环境中，本发明可以提供超过其中全部检查点在RAM中保存的Griewank检查指示策略1.97x的加速。
例子G：在存储器中的Griewank
S＝84G
快速存储器＝840G(10个重启状态缓冲)
正向成本＝269620
总成本＝309619
nreadc[F]＝24586 nwritec[F]＝14547 nreads[F]＝15413 nwrites[F]＝10039nreadc[S]＝0 nwritec[S]＝0 nreads[S]＝0 nwrites[S]＝0
例子H：在存储器和磁盘中的优化
S＝84G
C＝7G
快速存储器＝840G(用作状态或关联缓冲)
慢速存储器＝4200G(用作状态或关联缓冲)
正向成本＝116736
总成本＝156735
nreadc[F]＝39356 nwritec[F]＝39327 nreads[F]＝44 nwrites[F]＝29nreadc[S]＝266 nwritec[S]＝0 nreads[S]＝333 nwrites[S]＝266例子H的正向传播激励波场的逆时的成本具有超过例子G的2.3x的加速。例子H的全RTM应用具有超过例子G的1.97x的加速。
接下来的两个例子假设可以仅在RAM中支持5个重启状态缓冲并在磁盘上支持25个重启状态缓冲的计算机环境。在该环境中，本发明可以提供超过其中全部检查点在RAM中保存的Griewank检查指示策略2.9x的加速。
例子J：在存储器中的Griewank
S＝84G
快速存储器＝420G(5个重启状态缓冲)
正向成本＝483735
总成本＝523734
nreadc[F]＝12113 nwritec[F]＝3224 nreads[F]＝27886 nwrites[F]＝8889nreadc[S]＝0 nwritec[S]＝0 nreads[S]＝0 nwrites[S]＝0
例子K：在存储器和磁盘中的优化
S＝84G
C＝7G
快速存储器＝420G(用作状态或关联缓冲)
慢速存储器＝2100G(用作状态或关联缓冲)
正向成本＝138318
总成本＝178317
nreadc[F]＝38915 nwritec[F]＝38830 nreads[F]＝176 nwrites[F]＝85nreadc[S]＝350 nwritec[S]＝0 nreads[S]＝558 nwrites[S]＝350
接下来的两个例子假设可以仅在RAM中支持2个重启状态缓冲并在磁盘上支持10个重启状态缓冲的计算机环境。在该环境中，本发明可以提供超过其中全部检查点在RAM中保存的Griewank检查指示策略5.7x的加速。
例子L：在存储器中的Griewank
S＝84G
快速存储器＝168G(2个重启状态缓冲)
正向成本＝1804685
总成本＝1844684
nreadc[F]＝1914 nwritec[F]＝62 nreads[F]＝38085 nwrites[F]＝1852nreadc[S]＝0 nwritec[S]＝0 nreads[S]＝0 nwrites[S]＝0
例子M：在存储器和磁盘中的优化
S＝84G
C＝7G
快速存储器＝168G(用作状态或关联缓冲)
慢速存储器＝840G(用作状态或关联缓冲)
正向成本＝280124
总成本＝320123
nreadc[F]＝37505 nwritec[F]＝37295 nreads[F]＝757 nwrites[F]＝210nreadc[S]＝461 nwritec[S]＝0 nreads[S]＝1276 nwrites[S]＝461
接下来的两个例子假设可以仅在RAM中支持1个重启状态缓冲并在磁盘上支持5个重启状态缓冲的计算机环境。在该环境中，本发明可以提供超过其中全部检查点在RAM中保存的Griewank检查指示策略20x的加速。
例子N：在存储器中的Griewank
S＝84G
快速存储器＝84G(1个重启状态缓冲)
正向成本＝7502798
总成本＝7542797
nreadc[F]＝282 nwritec[F]＝1 nreads[F]＝39717 nwrites[F]＝281 nreadc[S]＝0 nwritec[S]＝0 nreads[S]＝0 nwrites[S]＝0
例子O：在存储器和磁盘中的优化
S＝84G
C＝7G
快速存储器＝84G(用作状态或关联缓冲)
慢速存储器＝420G(用作状态或关联缓冲)
正向成本＝336109
总成本＝376108
nreadc[F]＝35758 nwritec[F]＝35758 nreads[F]＝9 nwrites[F]＝0 nreadc[S]＝1262 nwritec[S]＝0 nreads[S]＝2970 nwrites[S]＝1262
最终的两个例子假设可以仅在RAM中支持0/5个重启状态缓冲并在磁盘上支持2.5个重启状态缓冲的计算机环境。在该环境中，本发明可以提供超过其中全部检查点在RAM中保存的Griewank检查指示策略893x的加速。
例子P：在存储器中的Griewank
S＝84G
快速存储器＝42G(0个重启状态缓冲)
正向成本＝8000200
总成本＝8000600
nreadc[F]＝0 nwritec[F]＝0 nreads[F]＝39999 nwrites[F]＝0 nreadc[S]＝0 nwritec[S]＝0 nreads[S]＝0 nwrites[S]＝0
例子Q：在存储器和磁盘中的优化
S＝84G
C＝7G
快速存储器＝42G(用作关联缓冲)
慢速存储器＝210G(用作状态或关联缓冲)
正向成本＝855154
总成本＝895153
nreadc[F]＝21768 nwritec[F]＝21768 nreads[F]＝27 nwrites[F]＝0 nreadc[S]＝14604 nwritec[S]＝14254 nreads[S]＝3600 nwrites[S]＝350因此如果小于重启状态缓冲的存储足够含有关联缓冲，那么本发明的算法能够使用该存储。
尽管本技术可以易受各种修改和可替换形式，但在上面讨论的示范实施例仅作为例子示出。然而，应再次理解本技术不希望限于在此披露的特别实施例。当然，本技术包括落入附加权利要求的真实精神和保护范围内的全部替换、修改和等效。
附录A
在下面示出的源代码是本发明的优化检查指示策略的未试验实施，其中仅对“状态缓冲”概念的“关联缓冲”概念的添加用于表明该设想。
附录B
该源代码提供包括二级存储优化检查指示策略的本发明方法的例子实施。

Claims (28)

1.一种用于降低计算机仿真的计算成本的方法，包含：
执行检查指示策略，其中所述检查指示策略包含：
选择用于存储关联检查点的时间步，其中所述关联检查点包含用于关联由所述计算机仿真生成的正向传播值的数据，并且其中所述计算机仿真不能从存储在所述关联检查点中的值重启；
其中所述选择用于存储关联检查点的时间步以在(i)和(ii)之间平衡：(i)根据存储关联检查点而不是全状态检查点减小计算机存储需求，而没有能力重启所述计算机仿真，(ii)根据存储全状态检查点减小计算时间，所述计算机仿真能够根据所述全状态检查点重启；
为所述关联检查点分配存储空间；以及
在多个时间步中的每个处运行所述计算机仿真，其中在所述多个时间步中的每个处：
源自测量数据的反向传播的仿真值与源自所述计算机仿真的所述正向传播值关联，其中所述正向传播值存储在所述关联检查点中。
2.根据权利要求1所述的方法，进一步包含：
存储全状态检查点；以及
从所述全状态检查点重启所述仿真。
3.根据权利要求1所述的方法，进一步包含确定存储全状态检查点或关联检查点的优化位置，其中所述优化位置是在图形处理器单元存储器即GPU存储器、随机访问存储器即RAM、RAM盘、磁盘驱动器或其任何组合中分配的存储空间，并且其中确定所述优化位置至少部分基于所述存储空间的访问速度。
4.根据权利要求1所述的方法，其中用于所述关联检查点的存储空间小于用于全状态检查点的存储空间。
5.根据权利要求1所述的方法，进一步包含通过多个步骤优化所述检查指示策略，所述多个步骤包含：
确定与存储所述关联检查点相关的计算成本；
至少部分基于多个存储单元中的每个的访问速度和全状态检查点需要的存储空间对所述关联检查点需要的存储空间的比率，将所述计算成本最小化；以及
生成包含所述多个时间步和在所述多个时间步中的每个处有待被存储的检查点类型的表格。
6.根据权利要求5所述的方法，其中确定所述计算成本包含：
计算能够存储在快速存储单元中的关联检查点的最大数目；
计算能够存储在慢速存储单元中的关联检查点的最大数目；
如果任何关联检查点能够存储在快速存储中：
计算如果仅使用快速存储的扫描的优化数目并添加到列表；
计算如果仅使用快速存储的扫描的最小数目并添加到所述列表；
如果任何关联检查点能够存储在慢速存储中：
计算如果使用全部类型的存储的扫描的最小数目并添加到所述列表；
计算如果使用全部的所述快速存储并在慢速存储中存储至少一个关联检查点的扫描的优化数目，并且如果小于列表上的最后值则添加到列表；
将最接近所述列表上的每个值的整数值添加到所述列表；
计算所述列表上每个值的所述成本；以及
返回最小成本。
7.根据权利要求1所述的方法，进一步包含从地震数据生成地下区域的图像。
8.根据权利要求1所述的方法，进一步包含将储层数据与储层仿真历史匹配。
9.根据权利要求1所述的方法，其中使用全双向波动方程执行时间正向传播仿真激励脉冲。
10.根据权利要求1所述的方法，其中使用基尔霍夫、波束或单向波动方程的迁移技术执行时间正向传播仿真激励脉冲。
11.根据权利要求1所述的方法，其中所述正向传播数据与下个时间步中的反向传播数据关联。
12.根据权利要求1所述的方法，进一步包含：
在快速存储媒体中存储多个关联检查点；以及
在慢速存储媒体中存储至少一个关联检查点。
13.根据权利要求12所述的方法，进一步包含存储至少一个全状态检查点，其中所述仿真能够从所述全状态检查点重启。
14.根据权利要求13所述的方法，进一步包含在快速存储媒体中存储所述至少一个全状态检查点。
15.根据权利要求1所述的方法，其中所述关联检查点由于以下项中的至少一项包含少于重启所述计算机仿真所需要的数据：
(a)用于所述计算机仿真的计算与用于所述关联的计算相比，利用更高的精确度执行；
(b)用于所述计算机仿真的计算需要边界条件信息，而用于所述关联不需要所述边界条件信息；
(c)不是所有被仿真的数据分量都被关联；以及
(d)所述计算机仿真比用于所述关联的所述计算使用更密集的栅格。
16.一种用于比较收集数据与仿真数据的方法，包含：
跨栅格反向传播所述收集数据；
运行计算机仿真从而跨所述栅格生成正向传播数据，其中：
在所述计算机仿真期间存储多个关联缓冲；以及
所述仿真不能从所述多个关联缓冲中的任何重启；以及
使存储在所述关联缓冲中的所述数据与在所述栅格中每个点处的所述反向传播数据关联；
其中至少部分通过平衡(i)和(ii)来选择所述计算机仿真中的时间步，所述关联缓冲在所述时间步处被存储：(i)根据存储关联缓冲而不是全状态缓冲减小计算机存储需求，而没有能力重启所述计算机仿真，(ii)根据存储全状态缓冲减小计算时间，所述计算机仿真能够根据所述全状态缓冲重启。
17.根据权利要求16所述的方法，包含在所述仿真期间存储全状态检查点，其中所述仿真能够从所述全状态检查点重启。
18.根据权利要求16所述的方法，进一步包含：
为存储全状态检查点计算成本；以及
至少部分基于所述成本确定存储全状态检查点的时间步。
19.一种用于使仿真数据与收集数据关联的系统，包含：
处理器；
存储系统，所述存储系统包含表示以下内容的数据结构：
测量数据；
传播算法，所述传播算法经配置跨栅格时间反向传播所述测量数据；
计算机仿真，所述计算机仿真经配置跨所述栅格生成时间正向仿真数据；以及
存储器，其中所述存储器包含代码从而引导所述处理器：
跨所述栅格时间反向传播所述测量数据；
用源自所述计算机仿真的时间正向仿真数据填充所述栅格；
在所述计算机仿真期间在选择时间步处存储关联检查点，所述关联检查点包含足以实行关联但不足以重启仿真的信息；其中选择用于存储关联检查点的时间步以在(i)和(ii)之间平衡：(i)根据存储关联检查点而不是全状态检查点减小计算机存储需求，而没有能力重启所述计算机仿真，(ii)根据存储全状态检查点减小计算时间，所述计算机仿真能够根据所述全状态检查点重启；以及
使所述时间反向传播数据与存储在所述关联检查点中的所述仿真数据关联。
20.根据权利要求19所述的系统，其中所述处理器包括单核、多核、图形处理单元或其任何组合。
21.根据权利要求19所述的系统，其中所述关联检查点存储在所述存储器中。
22.根据权利要求19所述的系统，其中所述存储器包含随机访问存储器即RAM、RAM盘、图形处理单元存储器或其任何组合。
23.根据权利要求19所述的系统，其中所述存储器包含代码从而在时间步处存储全状态检查点，其中所述计算机仿真能够从所述全状态检查点重启。
24.一种用于解决成像或反演地震数据的问题的计算机实施方法，其中时间步进正向地震仿真必须以逆时顺序访问从而与时间反向步进伴随计算关联，包含：
确定并实施检查指示策略，所述检查指示策略通过选择大小为C的nC个存储器缓冲和大小为S的nS个存储器缓冲的组合，针对包括数据存储的固定量的可用计算机存储器M减少计算时间，所述大小为C的nC个存储器缓冲按大小从正向仿真存储在已检查指示的时间步处执行关联需要的数据，而不是存储从所述时间步重启所述仿真需要的全部数据，并且所述大小为S的nS个存储器缓冲按大小存储从已检查的指示时间步重启所述正向仿真需要的全部数据，其中S>C，nC C+nS S≤M，nC≥1并且nS≥0；
其中至少部分通过平衡(i)和(ii)来选择所述已检查指示的时间步：(i)根据存储大小为C的缓冲而不是存储大小为S的缓冲减小计算机存储需求，而没有能力重启所述计算机仿真，(ii)根据存储大小为S的缓冲减小计算时间，所述计算机仿真能够根据所述大小为S的缓冲重启；
使用源自已选择的存储器缓冲的已保存数据在所述计算机上执行所述成像或反演。
25.根据权利要求24所述的方法，其中所述检查指示策略至少部分基于存储器资源的不同访问速度。
26.根据权利要求24所述的方法，其中所述检查指示策略包括为所述正向仿真中预选总数的时间步优化哪些时间步选为重启检查点。
27.根据权利要求26所述的方法，其中所述优化哪些时间步选为重启检查点基于使计算成本函数最小化。
28.根据权利要求24所述的方法，其中C取决于在所述成像或反演中使用哪个类型的关联。
CN201180024861.8A 2010-05-19 2011-03-14 在仿真期间用于检查指示的方法和系统 Expired - Fee Related CN102906728B (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US34629810P 2010-05-19 2010-05-19
US61/346,298 2010-05-19
PCT/US2011/028350 WO2011146161A1 (en) 2010-05-19 2011-03-14 Method and system for checkpointing during simulations

Publications (2)

Publication Number Publication Date
CN102906728A CN102906728A (zh) 2013-01-30
CN102906728B true CN102906728B (zh) 2017-04-26

Family

ID=44973196

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201180024861.8A Expired - Fee Related CN102906728B (zh) 2010-05-19 2011-03-14 在仿真期间用于检查指示的方法和系统

Country Status (11)

Country Link
US (1) US8756042B2 (zh)
EP (1) EP2572292A4 (zh)
KR (1) KR101931935B1 (zh)
CN (1) CN102906728B (zh)
AU (1) AU2011256799B2 (zh)
BR (1) BR112012025845A2 (zh)
CA (1) CA2796631C (zh)
MY (1) MY175034A (zh)
RU (1) RU2573269C2 (zh)
SG (1) SG184808A1 (zh)
WO (1) WO2011146161A1 (zh)

Families Citing this family (71)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
ES2341065B1 (es) * 2007-10-25 2011-06-13 Airbus España, S.L. Metodos y sistemas para mejorar mallas usadas en dinamica de fluidos computacional.
CA2771865C (en) * 2009-10-23 2016-04-05 Exxonmobil Upstream Research Company Method for optimization with gradient information
US8694299B2 (en) 2010-05-07 2014-04-08 Exxonmobil Upstream Research Company Artifact reduction in iterative inversion of geophysical data
SG185364A1 (en) * 2010-06-02 2012-12-28 Exxonmobil Upstream Res Co Efficient computation of wave equation migration angle gathers
US9366776B2 (en) * 2010-11-30 2016-06-14 Schlumberger Technology Corporation Integrated formation modeling systems and methods
US9046626B2 (en) * 2011-01-10 2015-06-02 Westerngeco L.L.C. Performing reverse time imaging of multicomponent acoustic and seismic data
US9291733B2 (en) * 2011-01-31 2016-03-22 Cggveritas Services Sa Device and method for determining S-wave attenuation in near-surface condition
EP2691795A4 (en) 2011-03-30 2015-12-09 CONVERGENCE SPEED OF COMPLETE WAVELENGTH INVERSION USING SPECTRAL SHAPING
WO2012160430A2 (en) * 2011-05-24 2012-11-29 Geco Technology B.V. Data acquisition
US9176930B2 (en) * 2011-11-29 2015-11-03 Exxonmobil Upstream Research Company Methods for approximating hessian times vector operation in full wavefield inversion
RU2612896C2 (ru) 2012-03-08 2017-03-13 Эксонмобил Апстрим Рисерч Компани Ортогональное кодирование источника и приемника
US20130311149A1 (en) * 2012-05-17 2013-11-21 Yaxun Tang Tomographically Enhanced Full Wavefield Inversion
US9665604B2 (en) * 2012-07-31 2017-05-30 Schlumberger Technology Corporation Modeling and manipulation of seismic reference datum (SRD) in a collaborative petro-technical application environment
US9429912B2 (en) * 2012-08-17 2016-08-30 Microsoft Technology Licensing, Llc Mixed reality holographic object development
CN102866423B (zh) * 2012-09-13 2015-07-22 浪潮(北京)电子信息产业有限公司 地震叠前时间偏移的处理方法和系统
CA2883169C (en) * 2012-09-28 2021-06-15 Exxonmobil Upstream Research Company Fault removal in geological models
US9857488B2 (en) 2012-11-20 2018-01-02 International Business Machines Corporation Efficient wavefield compression in seismic imaging
CA2892041C (en) 2012-11-28 2018-02-27 Exxonmobil Upstream Research Company Reflection seismic data q tomography
US10591638B2 (en) * 2013-03-06 2020-03-17 Exxonmobil Upstream Research Company Inversion of geophysical data on computer system having parallel processors
US9262560B2 (en) * 2013-03-13 2016-02-16 Saudi Arabian Oil Company Automatic recovery of reservoir simulation runs from processing system failures
WO2014149923A1 (en) * 2013-03-15 2014-09-25 Chevron U.S.A. Inc. Beam inversion by monte carlo back projection
US10088588B2 (en) * 2013-04-03 2018-10-02 Cgg Services Sas Device and method for stable least-squares reverse time migration
US9702993B2 (en) 2013-05-24 2017-07-11 Exxonmobil Upstream Research Company Multi-parameter inversion through offset dependent elastic 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
CN104376026B (zh) * 2013-08-18 2018-04-13 复旦大学 基于网格和多维树混合结构的表格查找方法
CA2913496C (en) 2013-08-23 2018-08-14 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
US10267937B2 (en) 2014-04-17 2019-04-23 Saudi Arabian Oil Company Generating subterranean imaging data based on vertical seismic profile data and ocean bottom sensor data
US9562983B2 (en) * 2014-04-17 2017-02-07 Saudi Arabian Oil Company Generating subterranean imaging data based on vertical seismic profile data
WO2015171215A1 (en) 2014-05-09 2015-11-12 Exxonmobil Upstream Research Company Efficient line search methods for multi-parameter full wavefield inversion
US10185046B2 (en) 2014-06-09 2019-01-22 Exxonmobil Upstream Research Company Method for temporal dispersion correction for seismic simulation, RTM and FWI
US9689999B2 (en) * 2014-06-13 2017-06-27 Pgs Geophysical As Seismic imaging using higher-order reflections
CA2947410A1 (en) 2014-06-17 2015-12-30 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
US9551210B2 (en) * 2014-08-15 2017-01-24 Carbo Ceramics Inc. Systems and methods for removal of electromagnetic dispersion and attenuation for imaging of proppant in an induced fracture
WO2016054008A1 (en) * 2014-10-02 2016-04-07 Conocophillips Company Pulsed marine source
US10386511B2 (en) 2014-10-03 2019-08-20 Exxonmobil Upstream Research Company Seismic survey design using full wavefield inversion
MY182815A (en) 2014-10-20 2021-02-05 Exxonmobil Upstream Res Co Velocity tomography using property scans
CA2963416A1 (en) * 2014-10-31 2016-05-06 Exxonmobil Upstream Research Company Handling domain discontinuity in a subsurface grid model with the help of grid optimization techniques
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
CA2972028C (en) 2015-02-13 2019-08-13 Exxonmobil Upstream Research Company Efficient and stable absorbing boundary condition in finite-difference calculations
US10670750B2 (en) 2015-02-17 2020-06-02 Exxonmobil Upstream Research Company Multistage full wavefield inversion process that generates a multiple free data set
US10317551B2 (en) * 2015-06-01 2019-06-11 Pgs Geophysical As Using seabed sensors and sea-surface reflections for structural imaging of a subsurface location in a geological formation
AU2016270000B2 (en) 2015-06-04 2019-05-16 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
KR102020759B1 (ko) 2015-10-02 2019-09-11 엑손모빌 업스트림 리서치 캄파니 Q-보상된 전 파동장 반전
CN108139498B (zh) 2015-10-15 2019-12-03 埃克森美孚上游研究公司 具有振幅保持的fwi模型域角度叠加
CN105354385B (zh) * 2015-11-13 2019-10-01 电子科技大学 一种半导体工艺仿真中的网格排序实现方法
EP3394642B1 (en) * 2015-12-22 2022-02-16 Shell Internationale Research Maatschappij B.V. Method and system for generating a seismic gather
WO2017188873A1 (en) * 2016-04-29 2017-11-02 Telefonaktiebolaget Lm Ericsson (Publ) Encoding and decoding using a polar code
US10768324B2 (en) 2016-05-19 2020-09-08 Exxonmobil Upstream Research Company Method to predict pore pressure and seal integrity using full wavefield inversion
CA3031422A1 (en) * 2016-07-22 2018-01-25 Schlumberger Canada Limited Modeling of oil and gas fields for appraisal and early development
WO2018053051A1 (en) * 2016-09-13 2018-03-22 Purdue Research Foundation System and method for divide-and conquer checkpointing
US10345466B2 (en) * 2017-07-25 2019-07-09 Advanced Geophysical Technology Inc. Memory efficient Q-RTM computer method and apparatus for imaging seismic data
US11644593B2 (en) 2018-10-11 2023-05-09 Landmark Graphics Corporation Calibrating time-lapse seismic images for production operations
CA3107201C (en) * 2018-10-11 2024-04-09 Landmark Graphics Corporation Calibrating time-lapse seismic images for production operations
US11073629B2 (en) * 2018-10-16 2021-07-27 Halliburton Energy Services, Inc. Method to improve DAS channel location accuracy using global inversion
US11080141B2 (en) 2019-01-22 2021-08-03 International Business Machines Corporation Automatic restarting and reconfiguration of physics-based models in event of model failure
US11360224B2 (en) 2019-05-03 2022-06-14 Exxonmobil Upstream Research Company Inversion, migration, and imaging related to isotropic wave-mode- independent attenuation
US11500116B2 (en) * 2019-05-15 2022-11-15 Saudi Arabian Oil Company Identifying characteristics of a subterranean region using vector-based wavefield separation of seismic data from the subterranean region
US11080856B1 (en) * 2019-06-18 2021-08-03 Euram Geo-Focus Technologies Corporation Methods for digital imaging of living tissue
US11346967B2 (en) * 2019-09-10 2022-05-31 Advanced Geophysical Technology Inc. Systems and methods for providing amplitude versus offset compliant migration gathers
US20220236439A1 (en) * 2021-01-23 2022-07-28 Manzar Fawad Rock physics model for shale volume estimation in subsurface reservoirs
KR102424206B1 (ko) * 2021-03-23 2022-07-22 국방과학연구소 시뮬레이션 모델기반 국방체계 설계를 위한 역방향 시뮬레이션 엔진 서버 및 그 역방향 시뮬레이션 방법
US20220413176A1 (en) * 2021-06-28 2022-12-29 Halliburton Energy Services, Inc. Annulus Velocity Independent Time Domain Structure Imaging In Cased Holes Using Multi-Offset Secondary Flexural Wave Data
US11867857B2 (en) * 2021-07-13 2024-01-09 Saudi Arabian Oil Company Method and system for updating a seismic velocity model
CN113706603B (zh) * 2021-10-28 2022-02-22 中国科学院地质与地球物理研究所 页岩有机质孔隙连通性分类表征的方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5504678A (en) * 1994-06-03 1996-04-02 Exxon Production Research Company Method for seismic data processing using depth slice decomposition
CN101017204A (zh) * 2006-02-09 2007-08-15 Pgs地球物理公司 三维双向声波动方程叠前成像系统和方法
CN101681394A (zh) * 2006-09-28 2010-03-24 埃克森美孚上游研究公司 来自并发地球物理源的数据的迭代反演

Family Cites Families (135)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US3812457A (en) 1969-11-17 1974-05-21 Shell Oil Co Seismic exploration method
US3864667A (en) 1970-09-11 1975-02-04 Continental Oil Co Apparatus for surface wave parameter determination
US3984805A (en) 1973-10-18 1976-10-05 Daniel Silverman Parallel operation of seismic vibrators without phase control
US4168485A (en) 1974-08-12 1979-09-18 Continental Oil Company Simultaneous use of pseudo-random control signals in vibrational exploration methods
US4545039A (en) 1982-09-09 1985-10-01 Western Geophysical Co. Of America Methods for seismic exploration
US4675851A (en) 1982-09-09 1987-06-23 Western Geophysical Co. Method for seismic exploration
US4575830A (en) 1982-10-15 1986-03-11 Schlumberger Technology Corporation Indirect shearwave determination
US4594662A (en) 1982-11-12 1986-06-10 Schlumberger Technology Corporation Diffraction tomography systems and methods with fixed detector arrays
US4562540A (en) 1982-11-12 1985-12-31 Schlumberger Technology Corporation Diffraction tomography system and methods
JPS606032A (ja) 1983-06-22 1985-01-12 Honda Motor Co Ltd 内燃エンジンの作動状態制御方法
US4924390A (en) 1985-03-04 1990-05-08 Conoco, Inc. Method for determination of earth stratum elastic parameters using seismic energy
US4715020A (en) 1986-10-29 1987-12-22 Western Atlas International, Inc. Simultaneous performance of multiple seismic vibratory surveys
FR2589587B1 (fr) 1985-10-30 1988-02-05 Inst Francais Du Petrole Procede de prospection sismique marine utilisant un signal vibratoire code et dispositif pour sa mise en oeuvre
US4707812A (en) 1985-12-09 1987-11-17 Atlantic Richfield Company Method of suppressing vibration seismic signal correlation noise
US4823326A (en) 1986-07-21 1989-04-18 The Standard Oil Company Seismic data acquisition technique having superposed signals
US4686654A (en) 1986-07-31 1987-08-11 Western Geophysical Company Of America Method for generating orthogonal sweep signals
US4766574A (en) 1987-03-31 1988-08-23 Amoco Corporation Method for depth imaging multicomponent seismic data
US4953657A (en) 1987-11-30 1990-09-04 Halliburton Geophysical Services, Inc. Time delay source coding
US4969129A (en) 1989-09-20 1990-11-06 Texaco Inc. Coding seismic sources
US4982374A (en) 1989-10-23 1991-01-01 Halliburton Geophysical Services, Inc. Method of source coding and harmonic cancellation for vibrational geophysical survey sources
GB9011836D0 (en) 1990-05-25 1990-07-18 Mason Iain M Seismic surveying
US5469062A (en) 1994-03-11 1995-11-21 Baker Hughes, Inc. Multiple depths and frequencies for simultaneous inversion of electromagnetic borehole measurements
GB2293010B (en) 1994-07-07 1998-12-09 Geco As Method of processing seismic data
US5583825A (en) 1994-09-02 1996-12-10 Exxon Production Research Company Method for deriving reservoir lithology and fluid content from pre-stack inversion of seismic data
US5924049A (en) 1995-04-18 1999-07-13 Western Atlas International, Inc. Methods for acquiring and processing seismic data
DE69625978T2 (de) 1995-04-18 2003-11-13 Western Atlas Int Inc Gleichmässige unteroberflächendeckung bei anwesenheit von steilen neigungen
US5721710A (en) 1995-09-29 1998-02-24 Atlantic Richfield Company High fidelity vibratory source seismic method with source separation
US5719821A (en) 1995-09-29 1998-02-17 Atlantic Richfield Company Method and apparatus for source separation of seismic vibratory signals
US5822269A (en) 1995-11-13 1998-10-13 Mobil Oil Corporation Method for separation of a plurality of vibratory seismic energy source signals
US5790473A (en) 1995-11-13 1998-08-04 Mobil Oil Corporation High fidelity vibratory source seismic method for use in vertical seismic profile data gathering with a plurality of vibratory seismic energy sources
US5715213A (en) 1995-11-13 1998-02-03 Mobil Oil Corporation High fidelity vibratory source seismic method using a plurality of vibrator sources
US5838634A (en) 1996-04-04 1998-11-17 Exxon Production Research Company Method of generating 3-D geologic models incorporating geologic and geophysical constraints
US5798982A (en) 1996-04-29 1998-08-25 The Trustees Of Columbia University In The City Of New York Method for inverting reflection trace data from 3-D and 4-D seismic surveys and identifying subsurface fluid and pathways in and among hydrocarbon reservoirs based on impedance models
GB9612471D0 (en) 1996-06-14 1996-08-14 Geco As Method and apparatus for multiple seismic vibratory surveys
US5878372A (en) 1997-03-04 1999-03-02 Western Atlas International, Inc. Method for simultaneous inversion processing of well log data using a plurality of earth models
US5999489A (en) 1997-03-21 1999-12-07 Tomoseis Inc. High vertical resolution crosswell seismic imaging
US6014342A (en) 1997-03-21 2000-01-11 Tomo Seis, Inc. Method of evaluating a subsurface region using gather sensitive data discrimination
US5920838A (en) 1997-06-02 1999-07-06 Carnegie Mellon University Reading and pronunciation tutor
FR2765692B1 (fr) 1997-07-04 1999-09-10 Inst Francais Du Petrole Methode pour modeliser en 3d l'impedance d'un milieu heterogene
GB2329043B (en) 1997-09-05 2000-04-26 Geco As Method of determining the response caused by model alterations in seismic simulations
US5999488A (en) 1998-04-27 1999-12-07 Phillips Petroleum Company Method and apparatus for migration by finite differences
US6219621B1 (en) 1998-06-30 2001-04-17 Exxonmobil Upstream Research Co. Sparse hyperbolic inversion of seismic data
US6388947B1 (en) 1998-09-14 2002-05-14 Tomoseis, Inc. Multi-crosswell profile 3D imaging and method
US6574564B2 (en) 1998-10-01 2003-06-03 Institut Francais Du Petrole 3D prestack seismic data migration method
FR2784195B1 (fr) 1998-10-01 2000-11-17 Inst Francais Du Petrole Methode pour realiser en 3d avant sommation, une migration de donnees sismiques
US6225803B1 (en) 1998-10-29 2001-05-01 Baker Hughes Incorporated NMR log processing using wavelet filter and iterative inversion
US6021094A (en) 1998-12-03 2000-02-01 Sandia Corporation Method of migrating seismic records
US6754588B2 (en) 1999-01-29 2004-06-22 Platte River Associates, Inc. Method of predicting three-dimensional stratigraphy using inverse optimization techniques
AU3229900A (en) 1999-02-12 2000-08-29 Prange, Michael Uncertainty constrained subsurface modeling
US6058073A (en) 1999-03-30 2000-05-02 Atlantic Richfield Company Elastic impedance estimation for inversion of far offset seismic sections
FR2792419B1 (fr) 1999-04-16 2001-09-07 Inst Francais Du Petrole Methode pour obtenir un modele optimal d'une caracteristique physique dans un milieu heterogene, tel que le sous-sol
GB9927395D0 (en) 1999-05-19 2000-01-19 Schlumberger Holdings Improved seismic data acquisition method
US6327537B1 (en) 1999-07-19 2001-12-04 Luc T. Ikelle Multi-shooting approach to seismic modeling and acquisition
FR2798197B1 (fr) 1999-09-02 2001-10-05 Inst Francais Du Petrole Methode pour former un modele d'une formation geologique, contraint par des donnees dynamiques et statiques
EP2296013B1 (en) 1999-10-22 2016-03-30 CGG Services (NL) B.V. Method of estimating elastic and compositional parameters from seismic and echo-acoustic data
US6480790B1 (en) 1999-10-29 2002-11-12 Exxonmobil Upstream Research Company Process for constructing three-dimensional geologic models having adjustable geologic interfaces
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
CN1188711C (zh) 2000-01-21 2005-02-09 施鲁博格控股有限公司 用于地震波场分离的系统和方法
CN1188710C (zh) 2000-01-21 2005-02-09 施鲁博格控股有限公司 估算地震介质特性的系统和方法
US6826486B1 (en) 2000-02-11 2004-11-30 Schlumberger Technology Corporation Methods and apparatus for predicting pore and fracture pressures of a subsurface formation
FR2805051B1 (fr) 2000-02-14 2002-12-06 Geophysique Cie Gle Methode de surveillance sismique d'une zone souterraine par utilisation simultanee de plusieurs sources vibrosismiques
GB2359363B (en) 2000-02-15 2002-04-03 Geco Prakla Processing simultaneous vibratory seismic data
US6687659B1 (en) 2000-03-24 2004-02-03 Conocophillips Company Method and apparatus for absorbing boundary conditions in numerical finite-difference acoustic applications
US6317695B1 (en) 2000-03-30 2001-11-13 Nutec Sciences, Inc. Seismic data processing method
CA2426160A1 (en) 2000-10-17 2002-04-25 David Lee Nyland Method of using cascaded sweeps for source coding and harmonic cancellation
AU2002239619A1 (en) 2000-12-08 2002-06-18 Peter J. Ortoleva Methods for modeling multi-dimensional domains using information theory to resolve gaps in data and in theories
FR2818753B1 (fr) 2000-12-21 2003-03-21 Inst Francais Du Petrole Methode et dispositif de prospection sismique par emission simultanee de signaux sismisques obtenus en codant un signal par des sequences pseudo aleatoires
FR2821677B1 (fr) 2001-03-05 2004-04-30 Geophysique Cie Gle Perfectionnements aux procedes d'inversion tomographique d'evenements pointes sur les donnees sismiques migrees
US6751558B2 (en) 2001-03-13 2004-06-15 Conoco Inc. Method and process for prediction of subsurface fluid and rock pressures in the earth
US6927698B2 (en) 2001-08-27 2005-08-09 Larry G. Stolarczyk Shuttle-in receiver for radio-imaging underground geologic structures
US6545944B2 (en) 2001-05-30 2003-04-08 Westerngeco L.L.C. Method for acquiring and processing of data from two or more simultaneously fired sources
US6882958B2 (en) 2001-06-28 2005-04-19 National Instruments Corporation System and method for curve fitting using randomized techniques
GB2379013B (en) 2001-08-07 2005-04-20 Abb Offshore Systems Ltd Microseismic signal processing
US6593746B2 (en) 2001-08-27 2003-07-15 Larry G. Stolarczyk Method and system for radio-imaging underground geologic structures
US7672824B2 (en) 2001-12-10 2010-03-02 Westerngeco L.L.C. Method for shallow water flow detection
FR2833384B1 (fr) * 2001-12-10 2004-04-02 Tsurf Procede, dispositif et produit programme de modelisation tridimensionnelle d'un volume geologique
US7069149B2 (en) 2001-12-14 2006-06-27 Chevron U.S.A. Inc. Process for interpreting faults from a fault-enhanced 3-dimensional seismic attribute volume
US7330799B2 (en) 2001-12-21 2008-02-12 Société de commercialisation des produits de la recherche appliquée-Socpra Sciences et Génie s.e.c. Method and algorithm for using surface waves
US6842701B2 (en) 2002-02-25 2005-01-11 Westerngeco L.L.C. Method of noise removal for cascaded sweep data
GB2387226C (en) 2002-04-06 2008-05-12 Westerngeco Ltd A method of seismic surveying
FR2839368B1 (fr) 2002-05-06 2004-10-01 Total Fina Elf S A Methode de decimation de traces sismiques pilotee par le trajet sismique
US6832159B2 (en) 2002-07-11 2004-12-14 Schlumberger Technology Corporation Intelligent diagnosis of environmental influence on well logs with model-based inversion
FR2843202B1 (fr) 2002-08-05 2004-09-10 Inst Francais Du Petrole 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
AU2003279870A1 (en) 2002-10-04 2004-05-04 Paradigm Geophysical Corporation Method and system for limited frequency seismic imaging
GB2396448B (en) 2002-12-21 2005-03-02 Schlumberger Holdings System and method for representing and processing and modeling subterranean surfaces
US6735527B1 (en) 2003-02-26 2004-05-11 Landmark Graphics Corporation 3-D prestack/poststack multiple prediction
US6999880B2 (en) 2003-03-18 2006-02-14 The Regents Of The University Of California Source-independent full waveform inversion of seismic data
US7184367B2 (en) 2003-03-27 2007-02-27 Exxonmobil Upstream Research Company Method to convert seismic traces into petrophysical property logs
WO2004095073A2 (en) 2003-04-01 2004-11-04 Exxonmobil Upstream Research Company Shaped high frequency vibratory source
US7072767B2 (en) 2003-04-01 2006-07-04 Conocophillips Company Simultaneous inversion for source wavelet and AVO parameters from prestack seismic data
NO322089B1 (no) 2003-04-09 2006-08-14 Norsar V Daglig Leder Fremgangsmate for simulering av lokale prestakk dypmigrerte seismiske bilder
GB2400438B (en) 2003-04-11 2005-06-01 Westerngeco Ltd Determination of waveguide parameters
US6970397B2 (en) 2003-07-09 2005-11-29 Gas Technology Institute Determination of fluid properties of earth formations using stochastic inversion
US6882938B2 (en) 2003-07-30 2005-04-19 Pgs Americas, Inc. Method for separating seismic signals from two or more distinct sources
US6944546B2 (en) 2003-10-01 2005-09-13 Halliburton Energy Services, Inc. Method and apparatus for inversion processing of well logging data in a selected pattern space
US6901333B2 (en) 2003-10-27 2005-05-31 Fugro N.V. Method and device for the generation and application of anisotropic elastic parameters
US7046581B2 (en) 2003-12-01 2006-05-16 Shell Oil Company Well-to-well tomography
US20050128874A1 (en) 2003-12-15 2005-06-16 Chevron U.S.A. Inc. Methods for acquiring and processing seismic data from quasi-simultaneously activated translating energy sources
FR2872584B1 (fr) 2004-06-30 2006-08-11 Inst Francais Du Petrole Methode pour simuler le depot sedimentaire dans un bassin respectant les epaisseurs des sequences sedimentaires
US7646924B2 (en) 2004-08-09 2010-01-12 David Leigh Donoho Method and apparatus for compressed sensing
US7480206B2 (en) 2004-09-13 2009-01-20 Chevron U.S.A. Inc. Methods for earth modeling and seismic imaging using interactive and selective updating
GB2422433B (en) 2004-12-21 2008-03-19 Sondex Wireline Ltd Method and apparatus for determining the permeability of earth formations
US7373251B2 (en) 2004-12-22 2008-05-13 Marathon Oil Company Method for predicting quantitative values of a rock or fluid property in a reservoir using seismic data
US7230879B2 (en) 2005-02-12 2007-06-12 Chevron U.S.A. Inc. Method and apparatus for true relative amplitude correction of seismic data for normal moveout stretch effects
WO2006090374A2 (en) 2005-02-22 2006-08-31 Paradigm Geophysical Ltd. Multiple suppression in angle domain time and depth migration
US7840625B2 (en) 2005-04-07 2010-11-23 California Institute Of Technology Methods for performing fast discrete curvelet transforms of data
US7271747B2 (en) 2005-05-10 2007-09-18 Rice University Method and apparatus for distributed compressed sensing
US7405997B2 (en) 2005-08-11 2008-07-29 Conocophillips Company Method of accounting for wavelet stretch in seismic data
RU2440604C2 (ru) 2005-10-18 2012-01-20 Синвент Ас Визуализация данных отклика геологической среды с использованием потоковых процессоров
US7373252B2 (en) 2005-11-04 2008-05-13 Western Geco L.L.C. 3D pre-stack full waveform inversion
FR2895091B1 (fr) 2005-12-21 2008-02-22 Inst Francais Du Petrole Methode pour mettre a jour un modele geologique par des donnees sismiques
GB2436626B (en) 2006-03-28 2008-08-06 Westerngeco Seismic Holdings Method of evaluating the interaction between a wavefield and a solid body
US7620534B2 (en) 2006-04-28 2009-11-17 Saudi Aramco Sound enabling computerized system for real time reservoir model calibration using field surveillance data
US20070274155A1 (en) 2006-05-25 2007-11-29 Ikelle Luc T Coding and Decoding: Seismic Data Modeling, Acquisition and Processing
US7725266B2 (en) 2006-05-31 2010-05-25 Bp Corporation North America Inc. System and method for 3D frequency domain waveform inversion based on 3D time-domain forward modeling
US7599798B2 (en) 2006-09-11 2009-10-06 Westerngeco L.L.C. Migrating composite seismic response data to produce a representation of a seismic volume
ATE543109T1 (de) 2007-01-20 2012-02-15 Spectraseis Ag Zeitumkehr-reservoir-lokalisierung
JP2009063942A (ja) 2007-09-10 2009-03-26 Sumitomo Electric Ind Ltd 遠赤外線カメラ用レンズ、レンズユニット及び撮像装置
US20090070042A1 (en) 2007-09-11 2009-03-12 Richard Birchwood Joint inversion of borehole acoustic radial profiles for in situ stresses as well as third-order nonlinear dynamic moduli, linear dynamic elastic moduli, and static elastic moduli in an isotropically stressed reference state
US20090083006A1 (en) 2007-09-20 2009-03-26 Randall Mackie Methods and apparatus for three-dimensional inversion of electromagnetic data
US20090164186A1 (en) 2007-12-20 2009-06-25 Bhp Billiton Innovation Pty Ltd. Method for determining improved estimates of properties of a model
US8577660B2 (en) 2008-01-23 2013-11-05 Schlumberger Technology Corporation Three-dimensional mechanical earth modeling
EP2105765A1 (en) 2008-03-28 2009-09-30 Schlumberger Holdings Limited Simultaneous inversion of induction data for dielectric permittivity and electric conductivity
US8494777B2 (en) 2008-04-09 2013-07-23 Schlumberger Technology Corporation Continuous microseismic mapping for real-time 3D event detection and location
US8345510B2 (en) 2008-06-02 2013-01-01 Pgs Geophysical As Method for aquiring and processing marine seismic data to extract and constructively use the up-going and down-going wave-fields emitted by the source(s)
WO2010019957A1 (en) 2008-08-15 2010-02-18 Bp Corporation North America Inc. Method for separating independent simultaneous sources
US20100054082A1 (en) 2008-08-29 2010-03-04 Acceleware Corp. Reverse-time depth migration with reduced memory requirements
US8127154B2 (en) 2008-10-02 2012-02-28 International Business Machines Corporation Total cost based checkpoint selection
US8296069B2 (en) 2008-10-06 2012-10-23 Bp Corporation North America Inc. Pseudo-analytical method for the solution of wave equations
US7616523B1 (en) 2008-10-22 2009-11-10 Pgs Geophysical As Method for combining pressure and motion seismic signals from streamers where sensors are not at a common depth
US9213119B2 (en) 2008-10-29 2015-12-15 Conocophillips Company Marine seismic acquisition
US20100118651A1 (en) * 2008-11-10 2010-05-13 Chevron U.S.A. Inc. Method for generation of images related to a subsurface region of interest
US20100142316A1 (en) 2008-12-07 2010-06-10 Henk Keers Using waveform inversion to determine properties of a subsurface medium
US8223587B2 (en) * 2010-03-29 2012-07-17 Exxonmobil Upstream Research Company Full wavefield inversion using time varying filters
US8437998B2 (en) * 2010-09-27 2013-05-07 Exxonmobil Upstream Research Company Hybrid method for full waveform inversion using simultaneous and sequential source method

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5504678A (en) * 1994-06-03 1996-04-02 Exxon Production Research Company Method for seismic data processing using depth slice decomposition
CN101017204A (zh) * 2006-02-09 2007-08-15 Pgs地球物理公司 三维双向声波动方程叠前成像系统和方法
CN101681394A (zh) * 2006-09-28 2010-03-24 埃克森美孚上游研究公司 来自并发地球物理源的数据的迭代反演

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Reverse time migration with optimal checkpointing;William W. Symes;《GEOPHYSICS》;20071031;第72卷(第5期);摘要,第2-4页 *

Also Published As

Publication number Publication date
AU2011256799A1 (en) 2012-11-29
US8756042B2 (en) 2014-06-17
BR112012025845A2 (pt) 2016-06-28
CA2796631A1 (en) 2011-11-24
WO2011146161A1 (en) 2011-11-24
CA2796631C (en) 2017-03-28
KR101931935B1 (ko) 2019-03-20
CN102906728A (zh) 2013-01-30
EP2572292A1 (en) 2013-03-27
RU2012154894A (ru) 2014-06-27
MY175034A (en) 2020-06-03
SG184808A1 (en) 2012-11-29
KR20130105792A (ko) 2013-09-26
EP2572292A4 (en) 2017-11-29
RU2573269C2 (ru) 2016-01-20
US20110288831A1 (en) 2011-11-24
AU2011256799B2 (en) 2016-02-25

Similar Documents

Publication Publication Date Title
CN102906728B (zh) 在仿真期间用于检查指示的方法和系统
Regone et al. Geologic model building in SEAM Phase II—Land seismic challenges
Zeng et al. An improved vacuum formulation for 2D finite-difference modeling of Rayleigh waves including surface topography and internal discontinuities
Xu et al. Numerical investigation of implementation of air-earth boundary by acoustic-elastic boundary approach
US20140129479A1 (en) Method to aid in the exploration, mine design, evaluation and/or extraction of metalliferous mineral and/or diamond deposits
US20100054082A1 (en) Reverse-time depth migration with reduced memory requirements
US8659974B2 (en) System and method of 3D salt flank VSP imaging with transmitted waves
CN103293552A (zh) 一种叠前地震资料的反演方法及系统
CN103926619A (zh) 一种三维vsp数据的逆时偏移方法
Lumley The Role of Geophysics in Carbon Capture
Rusmanugroho et al. 3D, 9C seismic modeling and inversion of Weyburn Field data
Huynh et al. Near-Surface Seismic Arrival Time Picking with Transfer and Semi-Supervised Learning
Hakso Analysis of microseismic multiplets in a hydraulic fracturing operation
Zheng et al. Improved numerical solution of anisotropic poroelastic wave equation in microseismicity: Graphic process unit acceleration and moment tensor implementation
Saleh et al. Uncertainty Quantification Through the Assimilation of CO2 Plume Size from 4D Seismic Survey
Rodgers Computational Approach for Improving Three-Dimensional Sub-Surface Earth Structure for Regional Earthquake Hazard Simulations in the San Francisco Bay Area
Levin SEG Technical Program Expanded Abstracts 2010
Gao Waveform tomography and its application at a ground water contamination site
Bilgi Investigation of Seismic Surveys and Enhancement of Seismic Images
Popovici et al. SEG Technical Program Expanded Abstracts 2017
Hardage SEG Technical Program Expanded Abstracts 2007
Alumbaugh et al. SEG Technical Program Expanded Abstracts 2018
Fomel Testing SEG bibliography from 16001 to 17001
Bartel et al. Modeling of time lapse seismic reflection data for CO2 sequestration West Pearl Queen Field.
Fomel Testing SEG bibliography from 19001 to 20000

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20170426

Termination date: 20210314

CF01 Termination of patent right due to non-payment of annual fee