在地震勘探中, 基于射线理论的地震初至波走时层析方法用于反演近地表速度模型, 对静校正量计算及偏移速度模型的建立至关重要。该算法适用于地形起伏较大、近地表结构复杂的地区, 但是实现过程较复杂, 计算量大、计算效率不高^[1-3]。

经过多年的研究, 许多学者提出了多种优化方法以提高近地表层析反演的计算效率, 其中并行计算是解决层析反演存储量大和计算量大的有效途径之一。罗省贤等^[4]使用消息传递接口(MPI)技术将求解大型稀疏矩阵方程的LSQR串行算法并行化, 在地震表层速度模型层析反演中得到了有效应用; 张建中等^[5]应用MPI技术在Linux环境下实现了Fresnel层析成像的并行计算, 具有较高的并行度和加速比; 魏星等^[6]应用高性能计算系统进行三维正演模拟, 对模型进行切割, 将计算任务划分至多节点, 通过MPI通信协同完成并行计算; 蔡志成等^[7]基于MPI并行变网格方法进行大规模三维正演数值模拟, 该方法能够减少内存消耗, 提高效率。然而这种基于共享存储、采用MPI并行的高性能计算方法存在一些不足: 共享存储的物理结构决定了所有进程都向主进程传递规约数据体, 当并行进程增多时, 极易形成网络I/O竞争, 不能通过扩大集群规模提高计算效率^[8]; 尽管MPI提供了并行计算标准及消息传递接口, 但是程序员需要通盘考虑程序的容错性、可扩展性、负载平衡等问题^[9], 使得程序设计及实现效果差异性较大。

近年来, 随着大数据技术的不断发展, Spark技术逐渐应用于地球物理的高性能计算领域。蒋楠^[10]基于Spark大数据处理框架进行了逆时偏移成像研究, 在保证成像精度的前提下实现计算效率的提升; 朱丽萍等^[11]进行了Spark框架下地震属性处理方法研究, 解决了单机环境下地震属性计算效率低的问题。为了进一步提高近地表层析反演的计算效率, 增强该技术的生产应用能力, 本文在大数据技术研究的基础上, 提出了一种基于Spark技术的近地表层析反演实现方法, 该方法充分利用Spark弹性分布式数据集(resilient distributed datasets, RDD)的技术特征, 构建基于RDD的数据组织及数据处理逻辑, 将三维反演数据体扁平化为深度方向的一维反演数据, 作为分布式数据规约单元, 有效解决了高性能环境下基于常规MPI并行实现产生的网络I/O问题, 提高了并行计算效率。

近地表层析反演的基本原理是将地下模型分解成一系列网格化的基本速度单元体, 通常依据地震初至波射线追踪进行正演, 以实际走时与正演走时之差迭代反演所有单元体的速度。主要步骤如下: ①在完成初至拾取的基础上, 建立慢度模型; ②根据观测系统参数, 正演计算所有炮的射线路径和旅行时; ③根据实际走时与正演走时之差以及射线路径, 利用同时迭代重建技术(simultaneous iterative reconstruction technique, SIRT)^[12]计算各基本速度单元的慢度修正量; ④利用慢度修正量修正慢度模型; ⑤如果走时残差无法达到给定标准, 则重新正演, 否则输出速度模型。

上述步骤中射线追踪正演最为耗时, 其计算时间与速度模型体、单炮道数以及三维网格单元体的不同面上射线穿过的节点密度等都有关系。在某些工区, 单炮计算时间多达30s, 甚至更多, 当炮集数据达到1.0×10⁵炮时, 串行运行至少需要833h才能完成一次迭代。事实上这些影响因素在不同炮内互不影响, 可以采用并行计算解决这些问题, 且理论上并行节点越多, 计算效率越高。

基于MPI的近地表层析反演常规并行程序设计通常基于共享存储, 并行程序在结构上采用主从模式, 主进程承担组织分配、数据收集、数据规约、模型平滑等工作, 从进程承担计算工作。在通信方面, 采用标准通信模式, 调用MPI消息传递库函数MPI_SEND()和MPI_RECV()分别实现对数据的发送和接收。

图 1为基于MPI技术的近地表层析反演并行实现流程图。中间为主进程, 左、右为从进程。其中, 蓝色箭头表示主进程使用MPI_SEND()函数分发各从进程的计算任务, 即有效炮集野外文件号(ffid)列表; 绿色箭头表示从进程使用MPI_RECV()函数传递其计算的中间数据给主进程, 主进程收集全部通信信息后, 求解慢度变化体, 更新速度模型体; 粉色箭头表示当走时残差不满足条件时, 需要进行迭代, 主进程通过MPI_SEND()函数传递给所有从进程更新后的速度模型和继续迭代的指令, 进入下一次迭代; 黄色箭头表示当走时残差满足条件或者达到迭代次数上限时, 主进程通过MPI_SEND()函数传递给所有从进程结束工作的指令, 主进程输出速度模型和射线密度模型, 程序结束。

图 1(Fig. 1)

图 1 基于MPI技术的近地表层析反演并行实现流程

当所有进程并行完成各自的计算任务后, 中间结果以网络通信方式传递给主进程, 生成慢度变化体和走时残差。当所有从进程的反演数据体传递给主进程时, 如绿色箭头线所示, 传输速率将受限于网络传输量、网络带宽、主进程所在节点的性能上限等。具体表现为: ①由于共享存储的I/O带宽有限, 当层析反演数据规模不变时, 使用更多的计算节点虽然能够有效缩短计算时间, 但I/O量也随之增长, 当节点规模持续增大时, 必将导致I/O竞争, 甚至堵塞; ②高性能计算环境中, 每个计算节点一般采用相同的物理配置, 因此每次迭代, 各个计算节点完成反演体计算的时间基本相同, 传输过程中几乎所有从进程同时向主进程发送规约请求, 此时程序从计算密集阶段进入I/O密集阶段, 会表现出明显的突发I/O特性, 短时间内产生大量的I/O请求, 甚至超过存储系统的最大处理能力, 造成严重的I/O竞争, 降低了I/O效率。

Spark最初是在2009年由UC Berkeley AMP实验室在Hadoop MapReduce基础上提出的高效大数据并行计算技术。目前Spark已经成为Apache基金会的一个顶级开源工程软件, 并成为互联网领域广为使用的基于内存的大数据并行计算系统^[13]。

大数据计算环境具有分布式存储、统一资源管理和并行计算框架等特征。结合地球物理反问题特点, 基于Spark技术的大数据环境可以利用下述5个特征提高算法计算效率。

1) Spark的核心概念是RDD, 它能统一使用集群上各节点内存, 充分利用集群总体内存存储数据;

2) Spark程序的核心是描述RDD分区的数据变换及RDD的生成逻辑链条有向无环图(directed acyclic graph, DAG), 尤其适合迭代型高性能计算的算法实现;

3) Spark与Hadoop无缝集成, 充分利用Hadoop分布式资源调度机制, 使得计算最大程度地向数据迁移, 这样能有效减缓网络I/O压力, 提高数据I/O效率;

4) Spark具有高度的容错机制, 若出现个别节点异常, 可根据RDD血统在其它节点对丢失的数据进行恢复, 有利于提高算法实现的稳定性;

5) Spark基于Scala语言开发, Scala是一种基于Java虚拟机(java virtual machine, JVM)纯面向对象语言, 完全兼容Java, 对于计算密集型算法可采用Java本地接口(java native interface, JNI)集成相应的C++算法, 提高计算效率。

本文算法具体实现过程如图 2所示, 其中, ffid_i表示第i炮的野外文件号, 通过多次正反演迭代逼近实际地下速度。相对于基于共享存储、常规MPI并行的实现方式, 在大数据环境下基于分布式存储, 采用Spark并行技术的程序设计过程包含持久化炮属性信息、复杂核心算法实现、设计规约逻辑和基于Spark Shuffle实现规约逻辑几部分。文中所使用的数据结构定义如下。

图 2(Fig. 2)

图 2 近地表层析反演算法流程

速度模型体网格点数量N_V=N_IL×N_XL×N_D, 其中, N_IL、N_XL、N_D分别表示整个工区主测线、联络线和深度方向的网格单元个数; 第i炮的有效速度模型网格点数量N_Vⁱ=N_ILⁱ×N_XLⁱ×N_D, 其中, N_ILⁱ、N_XLⁱ分别表示第i炮在主测线、联络线方向的网格单元个数; 第i炮的炮属性信息SA_i, 包含初至信息{f_ik}、炮点S_i的坐标信息{S_xi, S_yi, S_zi}; 与炮点对应的检波点D_ik的坐标信息为{D_xik, D_yik, D_zik}, 其中, k∈{1, 2, …, n_i}, n_i表示第i炮的检波点个数; 反演数据体IV=U_j=1^N_Vg_j, 第i炮有效反演数据体VIVⁱ=U^N_Vi_j=1·g_j, 深度方向的一维反演数据指对应主测线和联络线平面上搜索网格点(IL_m, XL_n)的深度方向的一维数据VIVⁱ_{(ILm, XLn)}=U_j=1^N_Dg_jN_D, 其中, IL_m∈[0, N_IL], XL_n∈[0, N_XL], g_j={t_j, w_j, rn_j}; 旅行时时差变化体ΔT=U_j=1^N_Vt_j, 第i炮有效旅行时时差变化体ΔVTⁱ=U^N_Vi_j=1t_j; 旅程权重体ΔW=U_j=1^N_Vw_j, 第i炮有效旅程权重体ΔWⁱ=U^N_Vi_j=1w_j; 射线密度体RD=U_j=1^N_Vrn_j, 第i炮有效射线密度体VRDⁱ=U^N_Vi_j=1rn_j, 其中, t_j=∑_k=1^rnj(a_kjΔt_k/l_k), w_j=∑_k=1^rn_ja_kj, Δt_k是经过第j个速度单元第k条射线的实际走时与正演走时的差, l_k是第k条射线的总长度, a_kj是第k条射线在第j个速度单元的长度, rn_j是经过第j个单元的射线总数, j∈[1, N_V]。

根据层析反演算法, 每次迭代都要使用炮属性信息。在单机中或是基于MPI并行方法中, 通常实时计算获取该数据信息, 但用时较长。本文利用分布式内存管理技术存储中间数据的特征, 计算一次, 然后持久化到分布式内存中。采用的持久化级别为MEMORY_AND_DISK, 即当内存不足时, 存放到磁盘, 因此能充分利用内存空间, 减少计算时间。例如, 某测试数据共3317炮, 每炮3200道, 需要持久化的炮属性信息数据量为493.6MB, 若在具有13个计算节点的集群上进行如下资源配置: 每个计算节点启动2个执行(Executor)进程, 单个进程的内存为5GB, 一个执行进程中同时运行两个任务(Task)线程。在这种条件下, 一次迭代计算炮属性信息时间为102s, 持久化炮属性信息时间则为30s, 因此, 迭代次数越多, 节省的计算用时就越多。

层析反演算法中, 炮检关系匹配、射线追踪、速度模型平滑等算法较为耗时。本文以C++函数实现, 程序设计时采用JNI接口调用, 并通过序列化和反序列化进行C++和Java空间内的数据传输和交换, 充分利用C++程序执行效率高的优势, 提高层析反演方法的运行效率。

根据前文定义的数据结构, 第i炮有效反演数据体比反演数据体小很多。例如A工区处理面积为427.20km², 采集范围21.36km×20.00km, 单炮覆盖面积8.0km×6.0km, 近地表模型反演深度2.5km, 网格距分别为50, 50, 10m, 则反演数据体网格点数量N_V=429×401×251, 第i炮有效反演数据体网格点数量N_Vⁱ=161×121×251, 占比为11.2%。

常规MPI并行方法以反演数据体为基本规约单元, 如1.2节介绍, 每个从进程在计算任务结束后将反演数据体通过网络传输给主进程, 反演数据体包含大量的无效数据, 传输过程中, 过多的无效数据加剧了网络I/O竞争, 降低了网络I/O性能。因此只传输有效数据是解决网络I/O瓶颈的关键。

本文以深度方向的一维反演数据作为规约基本单元, 且只对有效数据进行规约。在Spark并行计算中使用flatMap算子高效地将三维反演数据体扁平化为深度方向的一维反演数据集合。

文中规约逻辑的实现在执行过程中分为两个阶段, 同一节点上所有以主测线和联络线平面上搜索网格点(IL_m, XL_n)为键(Key)的规约基本单元进行规约; 不同节点上同一搜索网格点规约后的深度方向的一维反演数据再进行规约。两个阶段的特点如下。

1) 避免同一节点内的各进程出现网络I/O。

如图 3所示, 红色范围表示第i炮的覆盖范围, 绿色范围表示第j炮的覆盖范围, 两炮中有一部分重合, 定义红色点标识搜索网格点为(IL_m, XL_n), 假设某一计算节点启动两个执行进程, 一个进程负责计算第i炮的反演体, 另一个进程负责计算第j炮的反演体, 通过扁平化操作后, 两炮中以(IL_m, XL_n)为Key的规约基本单元数据不同, 规约后的值称为对应搜索网格点(IL_m, XL_n)的规约深度方向的一维反演数据。

图 3(Fig. 3)

图 3 集群某一节点上深度方向的一维反演数据的规约过程

Spark技术提供两种不同的算子完成上述功能, 分别为groupByKey和reduceByKey, 本文选择reduceByKey执行规约操作, 这是由于在规约阶段前, reduceByKey会在RDD存在的每个节点上对每一个Key做一次聚合, 即输出一个中间结果, 规约之后, reduceByKey算子内的匿名函数被调用, 此时会将RDD所有分区中间结果做规约, 从而产生一个最终的结果。groupByKey算子只是简单地把搜索网格点为(IL_m, XL_n)的键值对放在一起, 并不做任何额外的操作, 因此所有的键值对都被规约, 使得网络传输数据量增多, 网络I/O量增大。

2) 减少主节点规约I/O量。

利用Spark Shuffle减少主节点数据规约量。程序设计中不是所有的节点都同时链接主节点, 而是不同网格位置的反演数据在不同的节点间规约, 如图 4所示。所有节点以粉色点标识位置作为Key的反演数据, 都规约到节点1上, 所有节点以黑色点标识位置作为Key的反演数据, 都规约到节点2上, 同理, 所有节点以紫色点标识位置作为Key的反演数据, 都规约到节点3上。这种规约方式能够有效减少瞬间所有节点都与主节点链接的压力, 降低网络I/O堵塞发生的概率。

图 4(Fig. 4)

图 4 基于Spark技术的近地表层析反演实现在集群间进行数据规约过程

根据前文对层析反演算法数据流的分析和并行规约逻辑的阐述, 文中共设计8类不同的RDD, 如表 1所示。有效炮集列表RDD描述工区内有效炮集野外文件号列表集合; 炮属性RDD描述炮属性信息集合; 反演体RDD描述所有炮的反演数据体集合; 一维反演数据RDD是反演数据体RDD中有效数据的扁平化形式, 即一系列的基本规约单元集合; 一维规约反演数据RDD是对所有一维反演数据RDD的规约, 采用2.2节中的方式进行规约; 旅行时残差RDD描述了旅行时残差值‖ΔT‖_iter; 慢度变化体RDD描述了根据一维反演数据计算得到的慢度变化体; 射线密度体RDD描述迭代中的射线密度模型。

表 1

表 1 基于Spark技术的近地表层析反演RDD数据操作及特点 RDD 操作 函数目的 键 值 数据体量 有效炮集列表 parallelize 分发任务 炮索引号 ffid号 长整形数组 炮属性 map 组织炮结构 ffid号 炮属性结构 炮检位置信息和初至信息 反演体 map 反演函数 ffid号 反演体 速度模型体的3m倍, m表示炮数目 一维反演数据 flatMap 扁平化反演体 (ILm, XLn) 一维反演数据 所有反演数据体的有效部分, 可根据单炮覆盖范围占比计算 一维规约反演数据 reduceByKey 规约 (ILm, XLn) 规约后一维反演数据 速度模型体的3倍 旅行时残差 reduce 合并 旅行时残差 double类型值 慢度变化体 collect 收集 慢度变化体 速度模型体的1倍 射线密度体 collect 收集 射线密度体 速度模型体的1倍

表 1 基于Spark技术的近地表层析反演RDD数据操作及特点

当所有RDD设计完成后, 按照图 5所示的有向无环图实现基于Spark技术的层析反演计算, 该图随着迭代次数的增多而变长, 其依赖关系明确, 当某一块数据丢失时, 可以根据上述依赖关系重新生成所需数据, 其中绿色模块表示持久化的中间数据。

图 5(Fig. 5)

图 5 基于Spark技术的近地表层析反演一次迭代过程的有向无环图

图 6给出了大数据环境下基于Spark技术的近地表层析反演并行实现流程, 其在集群上的实现过程主要包含以下步骤:

图 6(Fig. 6)

图 6 大数据环境下基于Spark技术的近地表层析反演并行实现流程

1) 主节点负责进行炮检关系匹配、初始速度模型平滑、工区内ffid列表收集等任务, 并向其它执行节点广播初始速度模型, 由于节点内的不同进程可以共享数据, 因此向各执行节点广播一次即可;

2) Spark调度器根据有向无环图在各个执行节点分配计算任务, 主节点通过调用SparkContext的parallelize方法, 创建有效炮集列表RDD, 并根据初至数据所在位置将任务分发给各个节点的执行进程;

3) 每个执行进程首先整理并持久化炮属性信息, 然后以单炮为单位进行射线追踪, 得到反演数据体RDD;

4) Spark调度器在每个执行进程中将反演数据体RDD扁平化为一维反演数据RDD, 该RDD的每个分区是深度方向的一维反演数据;

5) Spark调度器在各个节点对一维反演数据进行规约, 形成一维规约反演数据RDD;

6) Spark调度器在各个节点执行reduce操作, 生成旅行时残差值, reduce后的一维规约反演数据仍旧分布在各个节点;

7) 主节点对一维规约反演数据执行collect操作, 将各个执行节点上的规约数据收集到主节点, 计算出慢度变化体和射线密度体, 并根据慢度变化体求出速度模型体, 对其进行平滑, 根据旅行时残差值确定是否进入下一次迭代。

为验证本文基于Spark技术的近地表层析反演程序设计应用效果, 选取国内西南地区A工区实际资料进行测试。工区地表高程变化比较大, 有些地表基岩出露, 静校正问题比较严重。该数据共48286炮, 所需参数如2.2.3节所述。测试集群包含125个节点, 每个节点配有2块16核CPU, 型号为Intel^Ⓡ Xeon^Ⓡ CPU E%-2682 v4@2.50GHz, 单个节点内存256GB; Spark作业并行进程数目为2000, 执行进程的内存配置为16GB, 每个节点上启动16个进程。

图 7给出了采用相同的输入数据和参数得到的基于Spark技术的近地表层析反演速度模型与基于MPI技术的近地表层析反演速度模型差值。两条剖面为第201条主测线和第198条联络线, 整体误差范围位于[-0.001, 0.001], 局部误差值较小且具有随机性, 表明反演结果一致。

图 7(Fig. 7)

图 7 基于Spark技术计算的近地表速度模型与基于MPI计算的近地表速度模型差值

以该工区为例, 对比分析采用MPI并行技术和Spark并行技术随着计算节点增多时程序的计算性能。令主测线方向、联络线方向和深度方向的网格距分别为50, 50, 10m, 3个方向的网格数目分别为429, 401, 251, 在这种网格规模下, 速度模型大小为328.7MB。基于MPI并行计算和基于Spark并行计算的层析反演程序分别在5节点、20节点、40节点、60节点、80节点和100节点集群规模下运行, 迭代次数为30次, 每个节点内存为256GB, 配有2块16核CPU。表 2给出了在MPI并行框架和Spark并行框架下该工区不同计算节点的网络I/O量和程序运行总时间。图 8给出了不同并行框架下的运行时间和并行加速比。

表 2

表 2 不同并行框架下集群规模与网络I/O量及程序运行时间 节点数 进程数 单进程计 算炮数 单炮平均计 算时间/s MPI网络 I/O量/GB Spark网络 I/O量/GB MPI程序运 行时间/h Spark程序 运行时间/h 5 100 483 40 98 0.81 164.0 162.0 20 500 97 40 478 2.30 36.9 33.1 40 1000 49 40 956 3.69 24.6 17.9 60 1500 33 40 1432 4.11 23.2 14.6 80 2000 25 40 1910 4.19 28.1 11.5 100 2500 20 40 2392 4.74 45.6 9.2

表 2 不同并行框架下集群规模与网络I/O量及程序运行时间

图 8(Fig. 8)

图 8 不同并行框架下的运行时间(a)和并行加速比(b)

由表 2和图 8可以看出:

1) 随着节点规模的不断扩大, 两种并行实现的网络I/O量都会增加, 但Spark并行实现的网络I/O量只有常规MPI并行实现的1/500~1/120, 且节点规模越多, 网络I/O量占比越小, 因此, Spark并行框架下通过网络传输的数据量相对于MPI并行具有明显优势;

2) 随着节点规模的不断扩大, 尤其是从5节点增加至40节点时, 两种并行实现程序的运行时间都在大幅减小, 当节点规模继续扩大时(达到80左右), MPI并行实现程序由于出现了网络I/O风暴导致程序运行时间延长, 而Spark并行实现程序随着并行节点的增多, 程序运行时间不断减少, 当节点规模达到100时, Spark并行程序运行时间只有MPI并行程序运行时间的1/4;

3) 随着节点规模的不断扩大, 基于Spark技术的层析反演和基于MPI的层析反演实现的并行加速比都逐步增大, MPI并行的加速比先增大后降低, Spark并行的加速比呈现类线性增长的趋势, 总体上, 节点规模越多, Spark并行计算相对于常规MPI并行计算的运行效率越高。

大数据环境下的Spark技术在地震勘探中的应用尚处于起步阶段, 虽然在地震数据属性提取、地震偏移成像等方面已经有一些研究, 但还没有像互联网领域的应用那么成功。本文针对常规MPI并行程序设计在近地表速度模型层析反演方面存在的不足, 采用Spark技术进行并行实现。

通过与基于共享存储、常规MPI并行计算层析反演实现方式的对比, 基于Spark技术的并行实现, 以一维深度方向的反演数据为基本单位的分布式规约效果明显, 均衡了网络I/O量, 有效解决了随着计算节点数的增大, 常规MPI高性能计算环境下近地表层析反演并行计算导致的网络I/O瓶颈问题, 实现了并行加速比的类线性增长。

近地表速度模型层析反演是一类迭代型密集计算问题, 本文基于Spark大数据技术, 构建了基于RDD的数据组织及数据处理逻辑, 形成了多次迭代的DAG, 该并行实现思路和实现框架具有较好的稳定性和对大规模集群的适应性, 为地震资料处理领域迭代型密集计算的实现提供了一种可借鉴的解决方案。