CN105474047A - 用于有限差分建模的统一应用框架 - Google Patents

Abstract

有限差分波场建模框架通过使用以数据流图形式的平台无关的中间表示使物理建模与硬件-软件优化的任务分离。在至少一些实施例中，波场模拟方法包括指定将要在空间中的每个点处应用的运算的内核以及将内核表示为平台无关的数据流图。对于多个实现平台中的每个平台，该方法还包括：从平台无关的图中导出用于实现平台的定制图；将定制图转化为用于实现平台的配置信息；将配置信息提供给实现平台，以获得波场模拟器；以及采用波场模拟器以利用指定的输入信号或输出信号将运算的所述内核反复应用到空间中的每个点。

Description

用于有限差分建模的统一应用框架

背景技术

科学家和工程师经常采用地球物理勘测用于勘探、考古研究和工程项目。地球物理勘测可以提供关于地下结构的信息，包括地层边界、岩石类型以及流体储层的存在或不存在。这些信息可以极大地帮助寻找水、地热储层和矿藏(诸如油气和矿石)。石油公司尤其经常投入大量的地震和电磁勘测以选择用于勘探油井的地点。

地球物理勘测可以在陆地上或水中进行。如在图1中所指示的，靠近感兴趣的区域104的能量源102生成波106，其传播进感兴趣的区域并且从诸如岩床边界的内部特征反射。最终，反射波108到达表面112上或附近或者钻孔中的接收器阵列110。记录系统114捕获接收到的信号用于存储和处理。利用多个不同的源位置并且任选地利用不同的接收器位置重复该过程。虽然存在各种方法用于将接收到的波信号转换为地下结构的图像，但是，最流行的这些技术采用有限差分波场建模过程，该过程使用离散时间步长和波函数导数的快速近似及时地向前或向后传播波。利用足够小的时间步长和足够精确的近似，有限差分波场建模过程提供了高效、高保真的解决方案。但是，有限差分模拟过程强加了大量的计算需求。

已经做出了很多努力来满足有限差分波场建模过程的计算需求。已实现了许多组合的硬件-软件解决方案来试图最大化运行时间效率。除了利用由不同的硬件系统提供的不同的并行策略之外，软件通常被定制，以结合不同的地球模型、不同的波函数近似以及不同的传播策略。这些因素导致严重的代码库分裂(fragmentation)，使得难以维护和重新使用应用。在许多情况下，由于略有不同的底层算法或硬件实现，可能甚至难以比较不同应用的结果。

附图说明

相应地，本文公开了采用用于有限差分建模的统一应用框架的系统和方法。

图1是地震勘测的示意图。

图2是具有有限差分波场模拟器的地震勘测系统的框图。

图3示出了说明性的数字化接收器信号迹线。

图4示出了表示感兴趣的地下区域的数据空间。

图5是说明性的地震成像方法的流程图。

图6是适于实现有限差分波场建模的说明性的硬件平台的框图。

图7说明了代码库分裂的特性形式。

图8说明了用于有限差分建模的统一应用框架。

图9示出了适于在基于数据流的中间表示中使用的说明性的节点类型的层级结构。

图10示出了说明性的基于数据流的中间表示。

但是，应该理解，在下面的附图和具体实施方式中给出的具体实施例并不限制本公开内容。相反，它们为普通技术人员提供了辨别涵盖在所附权利要求的范围内的可替代形式、等价物和其它修改的基础。

具体实施方式

以下的公开内容呈现了某些波场模拟系统和方法，其适于包括在各种硬件平台上地震成像的许多用途。从有限差分内核的平台无关的数据流图中，用于各种硬件平台的定制软件是可导出的。

图2提供了用于所公开的波场模拟系统和方法的说明性的上下文。地震接收器阵列110将地震波转换为放大的和数字化的电信号。(说明性的信号在图3中示出)。记录系统114经由总线202或其它通信路径收集数字化信号并且将数字化信号存储在信息存储介质上用于以后处理。通常，每个数字化信号与接收器位置和炮点位置204以及系统设计者认为值得的其它此类信息相关联。记录系统114可以执行一些初始化处理，以过滤和/或压缩数据，并且至少在一些情况下，执行质量控制。

记录系统114经由互联网或一些其它通信机制206将地震勘测数据提供给具有足够的计算资源的数据处理中心，以用于成像处理。数据处理中心包括充当有限差分波场模拟器208的一个或多个计算机。虽然这模拟器208(一个或多个)可以是被适当编程的通用计算机，但是能够通过使用专用硬件(诸如超级计算机、矢量处理单元、图形处理单元、现场可编程门阵列(FPGA)、数据流处理器以及它们的组合)来获得更高的效率。结果产生的波场被合并以形成地下结构的地图，其可以被在监视器210上显示并且被存储在持久性存储装置中用于以后使用。

为了模拟波场和构建图像，数据处理中心将感兴趣的区域建模为被划分成元素网格402的二维或三维空间。图4示出了三维空间的示例。在元素内的所有点的属性被假定是均匀的并且因此可以被表示为单个点。模拟器逐个元素地跟踪相关的属性(例如，压力、声速、各向异性)。

图5是可以由数据处理中心中的计算机执行的说明性的地震成像方法的流程图。在方框502中，计算机获得勘测数据，其包括数字化信号和相关联的源以及接收器位置信息。也可以在这个框中获得地下的初始模型，例如，以将估计的声速指定为位置的函数。通常，采用均匀速度模型作为起始点并且在成像过程中逐步细化。

在方框504中，计算机通过将源信号(或合适的替代物)从源位置传播到空间中来模拟由源生成的波场。也就是说，源生成按照波动方程从源位置向外传播的压力波、剪切波和/或其它类型的波。有限差分模型模拟这种传播，以确定源波场，即，模型空间中的每个点处作为时间的函数的压力或加速度。

在方框506中，类似的操作被应用到接收器信号。认识到，所记录的压力或加速度到达接收器位置，这些压力或加速度被根据波动方程在时间上向后传播到模型空间中。尽管具有反转的时间索引，但是可以使用与之前相同的有限差分传播模型来确定接收器波场。(随着时间反转，接收器波形被像源波形一样对待)。因为地球中的反射器将源波场转换为接收器波场，这两个波场在反射点处匹配。相应地，在框508中，源和接收器波场作为时间和位置的函数被分析，以识别匹配的那些区域并且因此指示反射点。针对每个地震炮点这一过程被重复并且匹配信息被相加(“堆叠”)，以对整个感兴趣的区域成像。

在方框510中，图像被分析，以识别估计速度中的不精确性。由于这种不精确性在图像数据中产生特定图案，因此它们使得速度模型能够被细化。方框504-510可以被迭代地重复，以逐步细化和改进速度模型(以及得到的图像)。一旦速度模型收敛，在方框512中计算机可以显示图像。然后可以基于显示的图像采取行动，包括转向决定、着陆决定和完成决定。

一个特定的注意点是，在达到收敛之前可能需要多次重复成像过程，并且成像过程中的每次迭代需要传播方程针对每次放炮以及针对每个时间步长在模型空间中的每个点处的解。因此，有限差分传播过程可能需要数千亿次计算。幸运的是，该过程包括显著的并行度，并且已经开发了多种硬件平台来增加这一过程的效率。在传播过程中这种效率增益的益处不限于成像。例如，也可以采用有限差分建模来为模拟目的生成合成炮点道集数据。

图6是适于实现有限差分波场建模的说明性的硬件平台的框图。一个或多个多核处理器602基于存储在存储器606中和/或存储在经由I/O桥608可访问的持久性存储装置中的软件协调系统的运算。总线桥604经由总线607将处理器602(一个或多个)和存储器606耦接到系统的其余部分。各种外围设备可以耦接到总线，以使得可被处理器602(一个或多个)访问。图示的平台包括一组具有可配置处理器和用于存储模型数据的大存储器的加速器卡610A-610D。加速器卡还可以包括可以提供有限差分传播过程的流水线硬件实现的现场可编程门阵列(FPGA)。处理器602(一个或多个)可以配置加速器卡以高效地实现波场迁移，并且(例如，在模型空间的不同部分上或者在不同的炮点上)并行地运行它们。

其它硬件平台可以提供迁移过程的不同实现，例如，采用图形处理器单元(GPU)、诸如XeonPhi^TM协处理器的协处理器、或其它矢量处理器来实现单指令多数据(SIMD)计算，或者依赖通用处理器用于测试新的算法。许多“通用”处理器现在包括对SIMD指令(例如，流SIMD扩展(SSE)和高级矢量扩展(AVX))的支持，其中每个都需要其自己的制造商特定的优化。作为这种情况的自然结果，趋向于发生代码分裂。图7表示由于为五种不同的平台的各种定制已被分裂为超过20个软件包的七种不同的迁移过程(例如，用于不同的近似，不同的优化和/或诸如各向同性、垂直横向各向同性(VTI)、倾斜横向各向同性(TTI)或完全各向异性的不同的地层模型的过程)。维护这种断裂的代码库是不可行的。

与此对照，图8说明了用于有限差分建模的统一应用框架，其使得单个有限差分建模应用802能够是平台无关的并且仍然能够被转化为高效的、平台特定的程序816、826、836、846。该框架使得应用开发人员能够指定传播“内核”，其描述将要在模型中的每个点处在每个时间步长执行的有限差分计算。开发人员可以采用具有有用函数和关键字的领域特定的语言用于这一目的。前端模块804将这个内核转换为按照有向、无环图特征化数据流的中间表示(“IR”)。一组平台特定的后端模块810、820、830、840中的一个执行针对中间表示的一次或多次优化通过(optimizationpasses)以获得细化的或定制的表示812、822、832、842，其是已针对给定平台优化的数据流图。相应的转化器(translator)模块814、824、834、844将定制表示转化为平台特定的程序816、826、836、846，其可以采取配置信息、定制的源代码、对象库和/或完全链接的可执行二进制文件的形式。将用于给定平台的后端和转化器模块合在一起可以被称为“构建器”模块，并且至少在一些实施例中，它们被合并为单个可执行模块。主机应用可以与平台特定的程序集成、链接或结合。例如，主机程序可以将平台特定的程序下载到加速板上，以将它们配置为用于传播过程的主要实现组件。

以这种方式，框架使得能够独立于用于选择的硬件平台的优化过程执行传播内核的规范和优化。其的使用从可能必须由应用开发人员解决的许多场景(例如，物理模型复杂性、不同算法选择)以及在异构硬件架构的设计中出现的许多开发中产生，以提供比同构或通用架构更好的计算效率。通常，这些技术集由不同的人(例如，地球物理学家相对计算机科学家)持有，并且硬件平台具有比通过良好维护的软件可获得的短得多的生命周期(2-3年)。平台无关的、领域特定的语言提供了物理建模与硬件-软件优化问题的分离，使得开发人员通过存在统一代码库的帮助能够封装和扩展其传播核心策略，所述统一代码库使得能够将一个内核用于每种底层物理的近似。优化专家可以独立于硬件开发周期而专注于促进良好定义的结构/范式。可以在用于代码库中的每个内核的硬件平台之间进行直接比较。

中间表示将内核数据流特征化为具有表示运算的结果的节点和连接节点以表示节点之间的依赖性的有向边的有向无环图。进入节点的边表示输入到那个节点，并且离开节点的边表示从那个节点输出。没有输入边的节点是表示从外部源提取的值的源节点。没有输出边的节点是表示被发送到外部目的地的值的汇聚节点(sinknodes)。

如在图9中所指示的，中间表示可以采用各种节点类型。各种类型在这里被呈现为(部分)面向对象的类层级结构。抽象类型以圆角示出，在中间表示图中可采用的具体类型以方角示出。例如，抽象有限差分节点(FDNode)类型在该图的顶部示出，有四个抽象类型从它垂挂：一元运算节点(NodeUnaryOp)、二元运算节点(NodeBinaryOp)、源节点和汇聚节点。一元节点类型的示例包括三种具体类型：否定、查找和卷积。否定是不言自明的。查找是利用输入数据作为对表的索引的表查找。卷积对输入数据应用卷积(过滤器)模板。二元运算节点既包括算术运算(加、减、乘等)又包括逻辑运算(与、或、异或等)。源节点类型包括常数值、“海绵”模板(对于衰减来自模型空间边界的伪影(artifacts)有用)、模型空间属性(例如，作为模型空间中位置的函数的声速)以及输入波场(例如，在模型空间中每个位置处的当前压力)。可以使用输入波场和输出波场节点来指定解压缩/压缩技术和用于存储时间步长之间的波场信息的存储器缓冲区。

大多数节点类型独立于其它数据点对每个输入数据点运算。也就是说，通过其XYZ坐标i，j，k表示3D模型空间中的每个点，产生输出y的节点的函数f可以被如下表示：

y_i，j，k＝f(a_i，j，k，b_i，j，k，c_i，j，k，...)

其中a、b、c等是输入，并且y是输出。表示标准算术、逻辑和关系运算符的节点都属于这一类。但是，一些节点可以提供对其它数据点的依赖性。例如，NodeConvolve利用可以选自各种模具形状或用任意系数值指定的模板对波场执行1、2或3D空间卷积。NodeOffset执行空间移位，例如，

y_i，j，k＝a_{(i+dx)，(j+dy)，(k+dz)}

也可以提供高层次节点，诸如NodeSponge–利用运行时间可配置的海绵轮廓执行边界海绵擦拭(sponging)。NodeLookup执行表查找，即

y_i，j，k＝f(a_i，j，k)

其中a是整数，并且f是运行时间可配置的查找函数。NodeCounter计算结果为当前空间索引的指定函数，例如，

y_i，j，k＝j

NodeMux基于整数索引输入提供作为选定的输入中的一个的输出。NodeConstant和NodeInputScalar计算结果分别为编译时间和运行时间可配置的常数值。NodeInputWavefield、NodeOutputWavefield、NodeInputEarthmodel是用于封装用于输入流和输出流的多个压缩和缓冲技术的节点类型的示例。

可以通过平台特定的后端模块提供附加节点类型，以表示平台特定的概念或结构或者支持平台特定的优化过程。例如，DFE后端模块可能会引入NodeBuffer来表示硬件缓冲。中间表示还可以结合支持后端操作的节点结构。例如，节点类型可以提供有后端模块(一个或多个)采用来逐节点控制代码生成的平台特定的属性对象。这些属性对象的值将不被内核或前端模块设置。相反，在平台特定的后端模块中指定的规则将控制这些值。当然，平台无关的节点属性可以由内核来设置，例如，由卷积节点所采用的模具，或者由偏移节点所采用的移位。另外，节点可以逐边地指定将要在对应的输入边或输出边上传送的特定类型的数据。说明性的类型可以包括布尔型(Boolean)、整数型、浮点型、和表达式(用于表示子图，其结果意图要针对整个体积(诸如针对导出的地球模型参数)被预先计算)。

内核映射出在单个时间步长期间为模型空间中的所有点执行的运算。时间步长包括利用新的运行时间数据配置一个或多个源节点、为图中的每个节点重新计算输出、并且从一个或多个汇聚节点中收集新的运行时间数据。为所有时间步长重复内核的循环在外部结构中实现。

作为示例，考虑受由以下波动方程约束的二阶各向同性模型：

$\frac{{\∂}^{2}p}{\∂t^2}=v^2{\▿}^{2}p+s\left(t\right)$

转换为有限差分近似并且以具有适当的边界海绵擦拭的伪码表达，波动方程变为：

表面上，利用适当的模具的当前压力场p[x,y,z]的卷积被用于近似在模型空间中的每个点(x,y,z)处的拉普拉斯(Laplacian)。基于当前压力场p[x,y,z]、先前的压力场pp[x,y,z]、平方速度dvv[x,y,z]、拉普拉斯和源能量，在每个点处计算输出压力场pp_out。然后，边界海绵被应用到当前压力场和输出压力场以抑制边界伪影。

然后，应用开发人员可以利用包括在中间表示中可用的每个节点类型的表示的领域特定的语言起草源代码。也就是说，应用开发人员能够将其应用的传播部分表示为简洁的数据流描述，而不是表示为常规的过程代码。领域特定的语言是完整的，因为应用开发人员能够为每种可能的运算编码，而无需考虑硬件特定的细节。这种设计的结果是使物理建模与硬件-软件优化分离，如下面进一步讨论的。开发人员自己不必关心过程问题，诸如在问题领域上循环，因为内核的数据流描述被隐含地平行应用到模型空间中的每个点。我们注意到，得到的代码密切反映上述离散波动方程的伪码表达式。在这个示例中，源代码变为：

//来自先前两个时间步长的波场

FDVarp＝io.waveFieldInput("p",1.0,6/*halo*/)；

FDVarpp＝io.waveFieldInput("pp",1.0,0)；

//地球模型(速度平方)

FDVardvv＝io.earthModelInput("dvv",1.0/4.0,0)；

//源函数

FDVarsource＝io.hostInput("source_in",1.0,0)；

//边界条件(海绵擦拭边界)

FDVarsponge＝boundaries.sponge(30)；

//计算拉普拉斯作为在所有轴上的卷积

FDVarlaplacian＝p.convolve(ConvolveAxes.XYZ,myStencil)；

//差分方程、源注入、以及海绵擦拭

FDVarresult＝p*2-pp*sponge+dvv*laplacian+source；

result＝result*sponge；

//输出

io.hostOutput("receiver",result)；

io.waveFieldOutput("pp_out",result)；

可注意到的差异包括变量类型的规范以及将类似Java对象的语法用于指定输入流和输出流的参数。出于优化原因，海绵擦拭运算相对于伪码以不同的顺序实现。最终，注意，从内核中省略循环结构。

图10以数据流形式示出了这个源代码的中间表示，即，作为有向图。它具有六个源节点902-912、八个操作节点914-928以及两个汇聚节点930-932。节点902是接受源流的主机-输入节点。节点904是接受用于模型空间中的每个点的平方速度属性的地球-模型输入节点。节点906是当前压力波场输入节点。节点908是常数值输入节点。节点910是先前压力波场输入节点。节点912是提供围绕模型空间的边界的衰减的海绵输入节点。

节点914产生来自节点908的常数与来自节点906的当前压力波场值的乘积。节点916产生来自节点910的先前压力波场值与来自节点912的海绵系数的乘积。节点918产生来自节点914和916的乘积的差。节点920是对当前压力波场运算以产生拉普拉斯的卷积节点。节点922产生拉普拉斯与来自节点904的平方速度属性的乘积。节点924产生来自节点918和922的值的和。节点926将该和与来自节点902的源输入值合并，以获得更新的压力波场值。节点928产生更新的压力波场值与海绵系数的乘积。

汇聚节点930收集节点928的输出用于存储为输出压力波场。汇聚节点932也收集那个输出用于供给到主机作为接收器信号。

该框架可以提供API机制，其使得能够逐节点地进行平台特定的优化，同时最小化对应用内核代码的污染。这种机制是具有可以被平台特定的后端模块应用的平台特定的“规则”的节点属性的设置。规则指定对数据流图(被称为组)的子集的属性-操作动作(被称为函子)。可以基于节点类型来定义组，以促进某些运算的优化。每个后端模块可以指定任意数量的定制函子；通常每个函子将设置一个或多个节点属性。例如，后端模块可以指定整数值节点属性foo，以及用于应用相关联的函子的方法setFoo(intfoo)。函子可以被合并，以创建较高层次的函子。例如，is(NodeConvolve.class)>>setFoo(3)使用预定义的“is”函子来创建新的函子，如果该节点表示卷积，则该新的函子将foo属性设置为3。这种API机制促进了计算机科学家和地球物理学家之间的分工，使得每个人都能够专注于其自己的能力。

该框架还可以提供用于中间表示的检查和操作的API机制。类似于那些在编译器图通过(graphpasses)期间所执行的，平台特定的后端模块可以采用API以将变换应用到数据流图。图通过通常被用于：图形检查(例如生成视觉表示)、优化(例如用左移位替代乘以2)、平台特定的操作(例如为DFE实现插入缓冲区)、以及最终目标呈现。该API机制使得定制图通过能够被定义并且被注册，以用于在执行中间表示的平台特定的优化中使用。

图通过可以为根据需要尽可能多的节点类型注册处理程序(handler)方法。然后，适当的处理程序可以适当地被后端模块自动地调用。缺省情况下，可以以拓扑顺序访问节点(意味着直到其全部上游节点都已被访问时节点才被访问)。但是，也可以为每个图通过指定其它的顺序，包括相反的顺序、按类型排序以及通过列表索引(用于任意处理)。

一些框架实施例还可以提供用于相对于给定物理内核评估多个硬件目标的API机制，从而促进未来硬件平台的评估和相关联的购买和管理决定。

本文所公开的框架采用基于适于指定有限差分建模内核的平台无关的中间表示阐述的领域特定的语言。该框架使物理建模与硬件-软件优化任务分离。它强加了用于代码库的基于物理的元结构并且产生非常便携且适用于后续优化的代码设计。因为开发人员可以将其注意力放到内核，代码库的维护变得容易。就期望平台特定的细化而言，它们可以被平台特定的后端模块和/或平台特定的转化器模块容纳，使得这些细化能够按照对那个平台的需要均匀地跨全部内核来应用。该框架还促进各种平台的性能之间的直接比较，因为各种平台中的每个平台都实现公共内核。

一旦以上公开内容被完全理解，对于本领域技术人员来说，许多其它的变化和修改将变得显而易见。以下权利要求意在被解释为涵盖所有此类变化和修改。

Claims (15)

1.一种波场模拟方法，所述波场模拟方法包括：

指定将要在空间中的每个点处应用的有限差分内核；

将内核表示为平台无关的有向无环图，其具有表示运算的结果的节点并且具有表示依赖性的边；以及

对于多个实现平台中的每个平台：

从平台无关的图中导出用于实现平台的定制图，其中所述导出包括执行基于图的优化；

将定制图转化为用于实现平台的配置信息；

将配置信息提供给实现平台，以获得波场模拟器；

采用波场模拟器以利用指定的输入信号或输出信号将运算的所述内核反复应用到空间中的每个点；以及

至少部分地基于由波场模拟器提供的波场生成图像。

2.如权利要求1所述的方法，其中所述信号包括地震炮点或接收到的地震波形。

3.如权利要求1所述的方法，其中所述配置信息包括用于处理器、协处理器或者图形处理器单元(GPU)的可执行软件或固件指令。

4.如权利要求1所述的方法，其中所述配置信息包括用于现场可编程门阵列(FPGA)、可编程逻辑器件(PLD)或专用集成电路(ASIC)的数据流规范。

5.如权利要求1所述的方法，其中所述基于图的优化包括将图运算拆分成组并且将平台特定的规则应用到所述组中的至少一个，以利用硬件特定的属性。

6.一种有限差分建模系统，所述有限差分建模系统包括：

具有用于多个实现平台中的每个平台的构建器模块的非临时性信息存储介质；以及

耦接到信息存储介质以执行构建器模块的一个或多个处理器，其中每个构建器模块致使所述一个或多个处理器：

获得有限差分内核的平台无关的数据流图；

利用针对平台无关的图的一次或多次优化通过从平台无关的图中导出平台定制的数据流图；以及

将平台定制的数据流图转化为配置信息，所述配置信息致使当容纳指定的换能器信号时，相关联的平台将有限差分内核应用到模型空间中的每个点。

7.如权利要求6所述的系统，其中所述信号包括地震炮点或接收到的地震波形。

8.如权利要求6所述的系统，还包括至少部分地基于由波场模拟器提供的波场显示图像的监视器。

9.如权利要求6所述的系统，其中所述配置信息包括用于处理器、协处理器或者图形处理器单元(GPU)的可执行软件或固件指令。

10.如权利要求6所述的系统，其中所述配置信息包括用于现场可编程门阵列(FPGA)、可编程逻辑器件(PLD)或专用集成电路(ASIC)的数据流规范。

11.如权利要求6所述的系统，其中作为所述优化的一部分，每个构建器模块处理图节点的组，以利用硬件特定的属性。

12.一种地震成像系统，所述地震成像系统包括：

多个可配置的实现平台，在被提供平台特定的配置信息之后，每个可配置的实现平台响应于指定的换能器信号模拟波场，所述波场是可合并的，以形成地下区域的地震图像；

计算机，其从有限差分内核的平台无关的数据流图中导出所述平台特定的配置信息、将所述平台特定的配置信息存储在非临时性信息存储介质中；以及

通用计算机，其合并模拟的波场，以获得地下图像。

13.如权利要求12所述的系统，其中所述平台包括至少一个通用多处理器计算机以及具有多个可配置硬件加速器的至少一个异构计算机。

14.如权利要求13所述的系统，其中所述加速器包括具有可配置数据流架构的现场可编程门阵列(FPGA)。

15.如权利要求13所述的系统，其中所述加速器包括图形处理器单元(GPU)、协处理器或数字信号处理器(DSP)。