地震逆时偏移成像计算是一种计算密集和数据密集的算法, 因此地震逆时偏移成像技术的推广应用在很大程度上依赖于高性能计算技术的发展[1-3]。近年来, 随着计算机技术的不断发展, 计算机的计算性能也随之有了较大提升, 与此同时, 多种并行计算框架的发展为高性能计算提供了强大的技术支持[4-6], 持续推动了地震逆时偏移成像技术的快速发展。FOLTINEK等[7]分析了基于GPU的逆时偏移在计算性能上的明显优势。ARAYA-POLO[8]提出了基于众核处理器的三维叠前逆时深度偏移成像技术。孔祥宁等[9]提出了利用多GPU联合的并行计算技术加速逆时偏移算法。李博等[10]提出了利用CPU/GPU作为数值计算中心并建立随机边界模型, 提升了逆时偏移的计算效率。张向阳等[11]基于GPU实现了逆时偏移的并行加速算法并与单程波动方程偏移进行了对比分析, 证明了基于GPU的逆时偏移的高效性。总的来看, 目前主要的效率优化方式采取的是集中存储框架下联合GPU的并行策略, 解决的是偏移算法中的计算密集问题。而随着地震勘探目标的日益复杂, 勘探技术向精细化发展, 单点高密度地震等技术的成功应用使得地震采集的数据量成倍增加, 常规计算框架的数据集中存储方式难以满足生产要求。此外, 地震数据处理周期不断缩短, 要求地震处理技术能够融合更多最新的计算机技术, 以更大程度提高地震处理技术的计算效率。分布式并行框架是一种可有效支持数据密集型应用程序的框架, 可进行海量数据的并行运算。Spark框架作为一种典型的分布式并行框架由伯克利大学开发, 属于该校的研究性项目, 由于其计算的优越性已经被应用于多个领域, 如聚类算法、文本分类以及数据监测和分析等[13-16]。
为了进一步提高逆时偏移成像技术的计算效率, 增强技术的生产应用能力, 本文以Spark框架为基础, 充分利用分布式文件存储系统海量数据管理以及Spark并行框架分布式内存计算高效等优势, 形成了一套计算高效的逆时偏移成像处理技术。
1 逆时偏移成像基本原理及并行计算逆时偏移成像是一种全波场成像技术, 而其算法核心则是正演问题, 对应的三维声波方程为:
| $ \frac{{{\partial ^2}\mathit{\boldsymbol{u}}}}{{\partial {x^2}}} + \frac{{{\partial ^2}\mathit{\boldsymbol{u}}}}{{\partial {y^2}}} + \frac{{{\partial ^2}\mathit{\boldsymbol{u}}}}{{\partial {z^2}}} = \frac{1}{{{\mathit{\boldsymbol{v}}^2}}}\frac{{{\partial ^2}\mathit{\boldsymbol{u}}}}{{\partial {t^2}}} $ | (1) |
式中:u为地表记录的压力波场; v为纵横向可变的介质速度。
成像条件的选择是地震偏移成像算法的关键之一, 可采用互相关算法进行成像, 即对同一时刻的接收点外推波场uR(x, y, z, t)和震源外推波场uS(x, y, z, t)进行互相关计算, 得到成像值, 如下式:
| $ I(x,y,z) = \int\limits_0^{{t_{{\rm{max}}}}} {{\mathit{\boldsymbol{u}}_{\rm{R}}}} (x,y,z,t){\mathit{\boldsymbol{u}}_{\rm{S}}}(x,y,z,t){\rm{d}}t $ | (2) |
式中, 被积函数uR(x, y, z, t)uS(x, y, z, t)表示t时刻对波场做一次成像运算, 积分说明成像空间I(x, y, z)中的像是多个时间步长成像的叠加。
逆时偏移成像的实现过程主要包括震源波场正向延拓, 接收点波场反向延拓以及互相关成像3个部分。计算时以单个炮集为基本单元, 每一个炮集的偏移计算相互独立, 各个进程无需进行数据交换。不同计算节点之间任务的起止时间互不影响, 而且同一个计算节点可以持续获取数据进行计算, 不必等待其它计算节点的任务同步, 这种多数据单程序的算法特征具有较强的并行性, 能在多种并行框架中实现。逆时偏移成像的并行计算流程如图 1所示, 将单炮数据进行任务分配, 各个进程完成各自作业后对成像结果进行规约, 得到最终的成像结果。
|
图 1 逆时偏移成像的并行计算流程 |
目前主要有3种高性能计算框架及并行编程接口可以实现逆时偏移算法的并行计算, 分别是MPI, MapReduce, Spark。下面对这3种实现方式进行对比分析。
MPI是一个信息传递应用程序接口, 其目标是高性能、大规模性和可移植性, 具有完备的异步通信功能以及较强的可扩展性, 是目前应用范围较为广泛的一种并行编程接口, 但其实现难度相对较大。MPI基于消息传递, 可以根据需求进行任务划分或数据划分, 一次作业的进程数由作业发布前通过一次调度完成, 若在集群上运行, 集群的各节点之间通过互联网络进行消息传递, 节点通信一般会产生较大的通信开销, 且各节点进程之间相互依存, 若一个进程任务失败, 则整个作业失败, 所有的计算将重新开始。此外, 作业的容错性较差。
MPI编程接口主要基于集中式存储, 底层缺乏分布式文件存储系统的支持, 因此难以满足海量数据的存储和读写。MPI的数据存储和数据处理是分开的, 每次计算开始时, MPI需要从数据存储节点读取需要处理的数据分配给各个计算节点, 然后进行数据处理。即使基于分布式文件存储系统进行MPI并行作业运算, 仍无法进行分布式文件读取, 而且在数据读写尤其针对TB级乃至PB级数据的数据密集型应用, 需要大量的数据交换, 此时网络通信时间将成为影响系统性能的重要因素, 导致并行性能大大降低。
基于MPI的偏移成像实现方式如图 2所示。
|
图 2 基于MPI框架的偏移成像实现方式 |
MapReduce是Google公司提出的一个软件架构, 用于大规模数据集的并行运算。MapReduce模型简单易用, 降低了软件开发人员在并行程序结构研发上的工作量, 为地震数据处理软件开发人员提供了新的算法开发思路。MapReduce实现过程主要分为两个阶段:Map(映射)阶段和Reduce(归纳)阶段。图 3给出了基于MapReduce框架的偏移成像实现方式, 其中Map阶段的主要任务是对分割好的数据进行并行处理, 处理后的结果传输给Reduce, 由Reduce对多个并发的Map进程的结果进行规约。每个运行作业中包含一个主节点和多个计算节点, 主节点负责所有的作业调度并对计算节点的作业情况进行监控, 各个计算进程间的通信纯粹是用文件去联系的, 即每个Map的处理结果都要以文件的形式写在分布式文件存储系统上, 然后Reduce阶段再读取Map阶段生成的数据, 这个特性使得MapReduce有着良好的作业容错性, 当某一级的某一个进程出错时, 主节点会重新调度这个进程到另外一个计算节点上重新运行, 但重新调度产生的开销会比较大。
|
图 3 基于MapReduce框架的偏移成像实现方式 |
Spark是一种专为大规模数据处理而设计的类MapReduce的通用并行框架, 它通过完善的分布式内存计算以及处理优化机制来提升作业的运行速度。与MapReduce不同的是, Spark处理的中间结果可以全部在内存中实现数据流转, 可根据需要将数据进行落盘。图 4给出了基于Spark框架的偏移成像实现方式, 在整个计算过程中不需要频繁读写数据, 以节省数据I/O带来的较大开销, 因而作业的计算性能得到大幅度的提升。
|
图 4 基于Spark框架的偏移成像实现方式 |
Spark的核心是建立一种弹性分布式数据集(RDD), 以实现数据在内存中的批计算, RDD是一种抽象概念, 代表了一种有容错机制的数据集, 其值存在于内存中, 并且只读不能修改, 但可以通过执行转换操作生成新的RDD, 即Spark算法实现过程中是通过RDD的转换等操作实现数据流转的。Spark作业执行过程中, 每个RDD的生成与转换都有记录并与对应的磁盘数据进行关联, 因此若某一阶段的RDD丢失, 则根据记录关系, 可以重新生成新的RDD, 以保证作业的高容错性。而记录RDD这种依赖关系的工具叫做有向无环图(DAG), Spark通过提前分析作业实现过程设计特定的DAG, 以设定作业运行所需要的全部操作。
3 基于Spark的逆时偏移成像技术基于Spark框架的逆时偏移成像技术, 其输入数据的存储以及读写基于Hadoop分布式文件存储系统(HDFS)。HDFS具有较为完备的错误检测和快速、自动的恢复机制, 实现了较高的容错性, 对于海量数据的管理具有重要作用。此外, HDFS是“块”级别的分布式系统, 即进行文件存储时会根据预设的数据块大小将文件进行分块, 并以键值(key-value)对的形式将各个数据块存储到不同节点上, 同时将键值对映射到内存中, 作为文件索引以便文件的快速定位和读取。不同于MPI的集中式存储方式, 基于HDFS的并行作业在计算时各节点优先读取存储在自己节点的数据进行处理, 避免了大量数据的网络传输, 实现“计算向存储的迁移”, 从而大大提高系统的计算效率, 这对处理TB级乃至PB级的海量数据具有很大优势。基于Spark框架的作业并行进程数是由存储数据的分块数目决定的, 为了保证作业运行的负载均衡以及并行处理, 本文主要采取了以下两项技术方案。
1) 为了充分利用集群计算资源并满足负载均衡, 根据需要处理的地震数据量和可用计算资源设定数据块大小(如128 MB)以生成合理大小的RDD。具体实现过程为, 若当前节点资源数为M, 总数据量大小为N, 那么一个数据块的大小可设定为n=N/(M×k), 其中k为并行粒度的合理动态调整, 以防止一个RDD数据量过大引起的并行粒度过粗。
2) 根据HDFS各计算节点读取数据时优先读取本节点数据的特点, 为了保证负载均衡, 采用数据均衡服务使得分块数据的存储满足均衡分布。具体实现过程为:首先HDFS的NameNode根据集群内每个DataNode磁盘使用情况生成一份数据分布分析报告, 并根据需要移动的数据分布情况计算具体数据块迁移路线图, 最后根据计算好的数据迁移路线图完成数据迁移, 以达到数据均衡标准。
以上两项技术方案的应用最大程度满足了作业的负载均衡, 输入数据的读取减少了网络通信开销, 极大地提升了数据读取效率。
由于Spark框架在算法过程中以输入RDD为并行单元进行偏移计算, 而RDD是以数据块的形式存储, 因此, 对数据块进行分割时, 需要对当前数据的道头先进行判断再进行分割, 以保证一个RDD内包含的炮集为整数。对于包含多个炮集的RDD数据在计算过程中各个炮集会按序进行计算并完成偏移结果的实时叠加, 中间结果无需落盘, 只在内存中流转, 因而无消息通信和数据等待时间, 降低了I/O耗时, 进一步提升了计算效率。每个输入的RDD(炮集数据)在偏移计算完成后会转换成新的RDD(偏移结果)并保存在内存之中, 新生成的RDD的键值对以主测线号-偏移剖面的形式存储, 每完成一个并行作业就会对新生成的偏移结果RDD按照既定的键值对进行叠加, 直至全部作业完成, 最后按照既定的键值对中的主测线号对数据进行排序生成最终的偏移结果并输出。
逆时偏移成像计算过程中另一个重要的环节是断点保护的实现, 由于本文方法基于Spark框架, 算法实现过程是以输入炮集RDD为并行单元进行计算而且数据均是在内存中流转, 因此在实现断点保护时以输入RDD的偏移结果为最小并行粒度将数据落盘, 同时记录当前已完成的RDD序号, 并将其作为断点续航的标志, 与MapReduce的中间结果全部落盘再规约的方式相比, Spark仅将已完成叠加后的偏移结果作为副本保存在HDFS上, 实际计算的数据体仍在内存中流转。作业续发时若启动断点续发, 作业首先会判断中间结果文件是否存在, 若存在就读取当前结果并跳过已经完成的RDD续发作业, 进行剩余RDD的偏移, 最终实现作业的断点保护。图 5给出了基于Spark框架的偏移成像算法实现流程图。
|
图 5 基于Spark框架的偏移成像算法实现流程 |
为了进一步说明本文方法的实现过程, 我们给出了主体算法的伪代码(如图 6所示), 同时也给出了断点保护的实现机制, 整个作业过程主要分为两个阶段(如图 7所示), 第一阶段主要是从分布式文件系统中读取地震炮集数据生成所需RDD即seisData, 并通过Spark的mapPartitions函数进行RDD转换和数据分区, 同时通过pair操作返回炮号-数据键值对, 根据返回的键值对进行单炮偏移并将每个单炮数据的偏移结果进行叠加存放于如图 6中所示的result变量中, 直至完成该RDD内的所有炮数据; 第二阶段主要是通过如图 7中所示的向量转换(Shuffle)操作对输入RDD的偏移结果按照键值重新转换得到新的成像结果RDD并进行偏移结果规约, 进一步对最终的偏移结果按键值排序并输出结果resData。
|
图 6 主体算法伪代码 |
|
图 7 有向无环图(DAG) |
为了检验基于Spark大数据处理框架逆时偏移成像技术的应用效果, 选取了国内西部某工区的一块实际资料进行了测试。该资料中浅层信噪比较高, 以串珠成像和小断裂成像为目标。工区总计施工36 516炮, 采集的数据体大小约为3.04 TB, 数据存储于HDFS文件系统, 并行读取, 使用的计算集群为125个Tesla P100的GPU节点, 每个节点两个GPU卡, 显存为16 GB, 共250个GPU卡, 每个节点的内存为256 GB, 每个节点配置24 TB的机械硬盘以及两个1 TB的固态硬盘, 配置集群HDFS总存储容量为2.6 PB。由此可见, HDFS充分利用了每个节点的存储资源, 无需大容量的集中存储就能实现PB乃至更高级别的数据存储。这里设定512 MB为数据块分割的基本单位, 共6 086个数据块, 即并行进程数为6 086, 平均每块卡处理约24块数据, 每块数据包含6炮数据。该设定的并行粒度合理, 节点间数据分配均衡, 因而可有效提升数据读取效率。成像结果参数为Nx=1 901, Ny=1 071, Nz=1 001, 即成像结果大小为7.6 GB。图 8给出了基于Spark框架下的两条线的偏移成像结果, 可以看出, 偏移成像结果具有较高的分辨率, 串珠收敛, 断层接触关系清晰。
|
图 8 偏移成像剖面 a line 390; b line 550 |
为检验基于Spark大数据处理框架的地震偏移成像技术的计算效率优势, 利用上述数据分别进行了3种并行计算实现方式的地震偏移成像计算, 由于3种方式的核心计算方法一致, 因此得到的偏移成像结果是一样的, 此外, 对应的单炮偏移时间效率相同。影响3种实现方式计算效率的主要因素是数据的I/O以及规约通信。其中, 利用基于MPI并行方式进行测试时, 由于存储和计算相互独立, 因而将数据集中存储在一套特定的存储节点上, 各计算节点通过网络通信进行输入炮集数据读取, 经统计, 网络传输速率均值稳定在150 MB/s。而基于MapReduce和Spark的并行方式均基于HDFS读写, 优先获取存储在本地磁盘的炮集数据, 而本地磁盘的读写速率一般稳定在400 MB/s, 即基于HDFS的数据读写效率是基于网络传输效率的2.6倍, 因此基于Spark框架的偏移算法在输入数据尤其针对海量数据时的效率要明显优于MPI。此外, 关于偏移结果数据规约, 此次数据共包含36 516炮数据, 对于MPI在此次测试集群中并行度设定为250, 总共250个进程, 平均每个进程处理146炮, 即在偏移过程中每个进程除了需要从管理输入的主进程0(图 2)进行146次数据读取之外, 在数据规约时还需和管理输出的主进程1进行146次成像结果的传输并落盘, 由此在实际计算过程产生了较大的数据通讯以及数据I/O开销。对于MapReduce, 在每完成一个Map偏移时, 其对应的偏移结果都会单独落盘, 进而在Reduce阶段对所有Map的结果进行叠加(图 3), 即基于MapReduce框架对应的偏移过程和规约过程完全分离, 该模式虽具有高度作业容错性但带来了巨大的数据I/O开销, 因而降低了作业效率。此次测试中, 炮集数据共分割为6 086个数据块, 即对应的Map数为6 086并对应产生6 086个成像结果文件, 每个文件大小为7.6 GB, 共产生45 TB中间结果, 在Reduce阶段, 对这些数据进行规约, 不仅对储存产生较大压力, 同时也造成巨大的I/O开销。而Spark框架基于分布式内存计算, 在偏移计算过程中数据在内存中流转直接完成叠加, 叠加耗时基本隐藏在偏移计算过程中, 没有额外的数据I/O开销, 计算效率的优势较为明显。数据测试时间的统计结果表明, 假设MPI的运行时间为1, 则基于MapReduce的运行时间为1.3, 而基于Spark框架的运行时间为0.9, 表 1给出了3种实现方式的数据读取速率和相对总运行时间。
|
表 1 3种并行实现方式运行时间对比 |
本次测试数据大小为3.04 TB, 这种规模的数据对集中存储方式的硬盘容量有较高的要求, 随着数据量的进一步提升, 集中存储的容量将难以满足存储要求。此外, 在实际计算过程中, 集中存储的数据频繁读取会大量使用网络通信从而严重影响数据读取效率, 此时基于Spark的分布式文件存储将体现出更大的优势。
5 结束语本文提出了一种基于Spark大数据处理框架的地震偏移成像技术, 并通过实际数据测试从数据I/O、规约通信等方面与MPI和MapReduce两种并行实现方式下的逆时偏移成像技术进行了对比论述, 验证了该技术的成像精度以及计算高效性。近年来, 随着地震勘探目标的逐渐精细化, 单点高密度地震技术不断发展, 由此产生了海量地震数据。本文方法充分利用了Spark框架分布式存储的特性以及分布式内存的优势, 实现了海量数据存储并有效提升了读写效率, 并且随着地震数据量的进一步增加, 分布式存储的优势将会得到更好的体现。本文方法尚未对逆时偏移的核心算法进行Spark框架的全面融合, 因此下一步可将分布式内存计算应用到逆时偏移核心算法如波场延拓过程中的波场存储等阶段, 从而进一步减少数据I/O, 提升算法效率。Spark作为一种专为海量数据处理而设计的高效通用计算框架, 可针对地震数据处理过程中的多个阶段进行优化, 基于Spark框架进行地震资料处理技术的开发, 对于提升海量数据的处理效率具有重要意义。
| [1] |
陈可洋. 地震波逆时偏移方法研究综述[J]. 勘探地球物理进展, 2010, 33(3): 153-159. CHEN K Y. Review of seismic reverse time migration methods[J]. Progress in Exploration Geophysics, 2010, 33(3): 153-159. |
| [2] |
YOON K, MARFURT K J, STARR W. Challenges in reverse-time migration[J]. Expand Abstracts of 74th Annual Internat SEG Mtg, 2004, 1057-1060. |
| [3] |
DUSSAUD E. Computational strategies for reverse time migration[J]. Expand Abstracts of 78th Annual Internat SEG Mtg, 2008, 2267-2271. |
| [4] |
孟祥宾, 隋志强, 王修银, 等. 地震处理多核异构并行计算通用框架研究[J]. 油气地球物理, 2014, 12(2): 11-16. MENG X B, SUI Z Q, WANG X Y, et al. The research of universal seismic processing frame on heterogeneous multi-core parallel[J]. Petroleum Geophysics, 2014, 12(2): 11-16. |
| [5] |
张凯, 秦勃, 刘其成. 基于GPU-Hadoop的并行计算框架研究与实现[J]. 计算机应用研究, 2014(8): 314-316, 322. ZHANG K, QIN B, LIU Q C. Study of parallel computing framework based on GPU-Hadoop[J]. Application Research of Computer, 2014(8): 314-316, 322. |
| [6] |
陈国良, 毛睿, 陆克中. 大数据并行计算框架[J]. 科学通报, 2015, 60(5/6): 566-569. CHEN G L, MAO R, LU K Z. Parallel computing framework for big data[J]. Chinese Science Bulletin, 2015, 60(5/6): 566-569. |
| [7] |
FOLTINEK D, EATON D, MAHOVSKY J, et al. Industrial scale reverse time migration on GPU hardware[J]. Expand Abstracts of 79th Annual Internat SEG Mtg, 2009, 2789-2793. |
| [8] |
ARAYA-POLO M, RUBIO F, CRUZ R D L, et al. 3D seismic imaging through reverse-time migration on homogeneous and heterogeneous multi-core processors[J]. Entific Programming, 2009, 17(1-2): 185-198. |
| [9] |
孔祥宁, 张慧宇, 刘守伟, 等. 海量地震数据叠前逆时偏移的多GPU联合并行计算策略[J]. 石油物探, 2013, 52(3): 288-293. KONG X N, ZHANG H Y, LIU S W, et al. Multi-GPUs parallel computational strategy of prestack reverse-time migration for mass seismic data[J]. Geophysical Prospecting for Petroleum, 2013, 52(3): 288-293. |
| [10] |
李博, 刘红伟, 刘国峰, 等. 地震叠前逆时偏移算法的CPU/GPU实施对策[J]. 地球物理学报, 2010, 53(12): 2938-2943. LI B, LIU H W, LIU G F, et al. Computational strategy of seismic pre-stack reverse time migration on CPU/GPU[J]. Chinese Journal of Geopyhsics, 2010, 53(12): 2938-2943. |
| [11] |
张向阳, 冯超敏, 赵书贵, 等. 一种基于GPU的逆时偏移并行算法[J]. 计算机应用与软件, 2013, 30(10): 304-307. ZHANG X Y, FENG C M, ZHAO S G, et al. A GPU-based parallel algorithm of reverse time migration[J]. Computer Applications and Software, 2013, 30(10): 304-307. |
| [12] |
崔艺馨, 陈晓东. Spark框架优化的大规模谱聚类并行算法[J]. 计算机应用, 2020, 40(1): 168-172. CUI Y X, CHEN X D. Spark framework based optimized large-scale spectral clustering parallel algorithm[J]. Journal of Computer Applications, 2020, 40(1): 168-172. |
| [13] |
于苹苹, 倪建成, 姚彬修, 等. 基于Spark框架的高效KNN中文文本分类算法[J]. 计算机应用, 2016, 36(12): 3292-3297. YU P P, NI J C, YAO B X, et al. Highly efficient Chinese text classification algorithm of KNN based on Spark framework[J]. Journal of Computer Applications, 2016, 36(12): 3292-3297. |
| [14] |
卞琛, 于炯, 修位蓉, 等. Spark框架并行度推断算法[J]. 电子科技大学学报, 2019, 48(4): 567-574. BIAN C, YU J, XIU W R, et al. Parallelism deduction algorithm for spark[J]. Journal of University of Electronic Science and Technology of China, 2019, 48(4): 567-574. |
| [15] |
陈天宇, 张龙信, 李肯立, 等. Spark框架中RDD缓存替换策略优化[J]. 小型微型计算机系统, 2019, 40(6): 1248-1253. CHEN T Y, ZHANG L X, LI K L, et al. Optimization of RDD cache replacement strategy optimization in spark framework[J]. Journal of Chinese Computer Systems, 2019, 40(6): 1248-1253. |
| [16] |
WHITE T. Hadoop:The definitive guide[M]. New York: O'Reilly Media, 2009: 1-250.
|