在英特尔® 集成众核架构(英特尔® MIC 架构)上优化地震成像处理

简介

地震成像的主要目标是详细地了解地下参数,以便定位特定目标。所有地震采集工作都从传播和反射(还有其它地震波现象)到地面的地震波入手。地震波反射是由于不同地层之间的强烈参数对比而产生的,如密度和压缩波速度。这些反射波将被地球表面或专用井内的水听器(海洋地震)或地音探听器(陆上地震)记录下来。

震源信号的初始频率及其相关波长在很大程度上决定了可能的探测深度,因为地面首先会显著降低最高频率,然后再减小可见对象,这是因为空间分辨率与频率组成密切相关。这便使地震成像有了多种用途,比如在土木工程行业中,通过使用从几厘米到几米的空间采样,可将目标定位到前几米到几百米之间。地震成像运用最广泛的莫过于石油与天然气行业,一个几米到 25 米的空间采样最多可将目标定位到地面以下 10 至 12 公里的位置。

从更大范围而言,波的传播被广泛用于地震学中,包括全球对地成像和本地地震危害分析。在地震学中,震源信号是由地震产生的,而不是由气枪(海洋地震勘探所用)或振动卡车(陆上地震勘探所用)产生的。由于不知道源点位置,因此可以分析许多其他方面的差异,这与波发射的初始时刻不同,这一时刻的源点位置是已知的。此外震源信号也复杂得多,这是因为地震的裂缝是沿着一个断层断裂的,而地震勘测有一个源反射点。

地震科学和地震勘测之间最重要的区别之一就是所记录的地震数据量。地震学家可以利用全球多个地震检波器网络所提供的某次特定地震的一些记录(震波图)。在另一方面,石油行业进行了专门的勘测(比如说针对 1,000 km² 的面积)。在勘测中,一艘或多艘震波勘测船带有 8 10 公里长的地震检波器拖缆,且每隔 25 就配备一台水听器。随后,气枪将会每隔 10 到 20 秒的时间和 25 50 的距离发出地震信号,并且每台水听器都会以 2 毫秒的典型抽样时间来记录反射波。

这样可以轻松、详细地计算数据集的大小,并生成 PB 级数据,从而在进行进一步的地质解释、油藏描述或油藏模拟之前进入知名的“地震数据处理序列”(图 1),此时可获得地震图像。

图1 地震数据处理序列概要。实地探测后的第一步就是使用几个信号处理过滤器来处理观测到的数据,以消除所有多余的噪音。第二步包括针对地震波的传播位置定义初始模型(速度、密度和各向异性),然后使用叠前资料和初始模型执行地震成像。

地震数据处理

地震数据必须要在一个所谓的预处理序列中进行分类和清理。

这种对实地捕获数据的清理在很大程度上取决于我们模拟地震波传播以及将我们透过标准地震波反射所看到的每个现象模拟或“解释”为多个反射波、杂波和折射波的能力。

在 20 世纪 80 年代,地震成像被描述并重新拟定为一个反问题,涉及到特定地下物理描述(或模型)的观测数据(实地记录)和计算数据(在计算机上模拟)之间成本函数的最小化。

也就是说,预处理必须要消除计算机模拟无法轻易重现的实地记录的所有地震波现象。

地震成像或所谓的“偏移”(在我们试图将表面记录的事件偏移回实际深度时)极度依赖于波传播技术。

偏移技术的分类是输入数据的一个函数(叠后或叠前),偏移算子的应用领域(时间或深度)以及用于波模拟的波动方程的近似。

为了独立提取每条地震道(例如在水听器记录中)的所有信息,可以将输入数据视为 “叠前”,这样对复杂的地质区域最为直观。叠后数据,即“堆叠”的输入数据,拥有数据集小、信噪比高的显著优势,但它破坏了每台独立水听器所记录的所有特定信息。

就应用领域而言,时间偏移通常会导致一种快速方法,由于在传播介质上执行的强逼近,这种方法展现复杂结构的能力较差。如前所述,从反问题方面进行重新拟定的深度偏移将使用完整的地下物理和异构描述并将取得较好的结果。

由于存在多种组合,因此产生了多种用于处理各种区域(包括非常简单的区域和非常复杂的区域)的方法,如叠后时间偏移和叠前深度偏移。

就数值法而言,偏移技术与用于解决正演问题的方法有关。您可以基于渐进论来使用波动方程积分解,由此可以引出一种称为叠前深度克希霍夫偏移法或射线束偏移法的偏移技术。

其它实施方法均基于微偏分波动方程的数值积分,如用于这项工作的单向传播近似。这可以引出波场延拓偏移法,也就是波动方程偏移 (WEM)

对于最复杂的探测和传播区域,处理全波场的最佳方法是仍使用双向微分方程,这可以引出知名的逆时偏移 (RTM) 和全波形反演 (FWI) 算法。

分步傅里叶叠前深度偏移 (PsDM)

在这项工作中,我们会考虑一种在频域中进行波传播的单向近似表达方法。这种方法称之为分步傅里叶 PsDM (SSF PsDM) 算法。这项技术极为高效,特别适合轻度的横向速度变化,是中石化 iCluster* 地震成像系统的关键模块之一,曾在国内石油工业领域内应用并普及开来。

我们将 SSF PsDM 算法移植到 2010 5 月推出的英特尔® 至强融核协处理器中。

3D SSF PsDM 是一种基于共炮的偏移方法,这意味着与每次波发射相对应的每个炮集以及每个接收器所记录的相关反射波可以进行独立处理。这些波是从源位置传播的,并从每台接收器的位置传播回的,但由于单向近似,因此仅以一个方向传播,这在上文中已经提到过。两个波场的互关联是在相同的深度进行的,以获取图像。这一过程将重复进行,一直到最大深度。最终,每个偏移的虚炮点将堆叠在相同的成像点上,以形成最终图像。

图2 分步傅里叶 PsDM 流程的卡通表现方法。波是从源位置传播的,并从接收器位置传播回的。通过互关联两个波场,成像条件将被运用到各个深度切片中。

对于特定炮点,使用图 2 中的成像滤波方程可获得分图像,其中 Pimg 是炮集的偏移基值,Ps* (Xs,W) 是源波场,Ps* (Xs,W) 是记录的向下传播的波场,Xr, Xs 是源和接收器的空间坐标。由于传播是在频率空间域中进行的,我们将 nw 表示为所涉及的频率数。最终图像是通过叠加每个炮集基值获得的。

在英特尔® 至强融核协处理器上实施

基准代码和工作负载简介

图 3 显示了采用伪代码的 SSF PsDM 算法。请注意,SSF PsDM 内核使用了由一个深度域循环、一个频域循环和两个空间循环组成的四嵌套循环。我们的应用程序采用分布式并行(通过英特尔® MPI 库实施)和主从模式,在多个英特尔® 至强融核协处理器上同时处理多个炮集数据集。在英特尔® 至强融核协处理器或单个节点的 CPU 上,我们使用 OpenMP* 来并行化处理每个炮集的计算,这两种架构的代码是统一的。


// SSF PsDM kernel
!$omp parallel
IZ loop: Depth domain loop

!$omp do
   IW  loop: Frequency domain loop
      2D DFT forward for source wavefield
      2D DFT forward for recorded wave field
      phase shift calculation
      2D DFT backward for source wavefield
      2D DFT backward for recorded wave field
      Time-shift calculation, heavy sin/cos computation
   End of IW loop

   Imaging: wavefield imaging and stack up 
End of IZ loop 

图3 SSF PsDM 内核伪代码

在这项研究中,我们用了一个总大小为 73 GB 的否定结构工作负载,该负载有 4,164 个炮点,6 秒的记录时间以及 2 毫米的采样间隔。由于这些炮点具有相同的特征(就接收器数量和记录长度而言),因此我们在下面的性能评估中只计算了几个炮点。

我们使用了基准代码英特尔® 数学核心函数库中的离散傅里叶变换 (DFT) 例程。如图 3 所示,单次迭代中有四个 2D DFT 调用,导致实际工作负载中有大量的 (IZ*IW*4) DFT 调用。  例如在否定结构工作负载中,每个炮集约调用 2D DFT 1,250,000 次,这是一个计算密集型和高度并行化的任务,在这样的任务中,卸载和传输数据的时间可以忽略不计。

并行化和卸载 SSF PsDM 内核

第一代英特尔® 至强融核协处理器采用 22 纳米制程。该协处理器中包含一个 PCIe 卡,共有 61 个内核,时钟频率为 1.1 GHz。英特尔® 集成众核架构是一个基于 x86 的众核处理器架构,采用小型的按次序内核,独一无二地融合了如今通用 CPU 架构的全面可编程性与现代 GPU 架构的计算吞吐量和内存带宽能力。

使用英特尔® 至强融核协处理器的开发人员将可以获得以下重要优势:能够使用现有的标准编程工具和方法,如 CC++ Fortran,以及相应的并行化技术,如 OpenMP*MPI 和英特尔® 线程构建模块(英特尔® TBB)。

英特尔® Composer XE 2013 包含语言扩展,能够对节进行卸载,以便在英特尔 MIC 架构上运行。面向英特尔® 至强融核协处理器卸载计算的语言扩展(包括相关代码和数据)由 C/C++ 编译指示和 Fortran 指令组成。

为英特尔® 至强融核协处理器编写的相同源代码可以在标准的英特尔® 至强® 处理器上编译和运行。熟悉的编程模式消除了开发人员培训障碍,使得他们能够集中精力解决问题,而不是进行软件工程设计。

例如,标准 OpenMP* 可用来并行化处理图 3 中的内核,在图 3 中我们可以看到,OpenMP* 并行化是通过内部 IW 循环进行的。由于包含数据依赖关系,因此外部 IZ 循环没有进行并行化处理。

接下来,我们使用 Fortran 指令将整个 OpenMP* 区域(见图 4)卸载到英特尔® 至强融核协处理器上。(图 4 所示的 SSF_PSDM_MIC_WRAPPER() 函数是一个包装程序,包含图 3 中的并行内核)。


!dec$ offload target(mic)
     +     in (…)   // data needs to be copied in
     +     inout (…)  // data needs to copied in and out 

        CALL SSF_PSDM_MIC_WRAPPER (…)

图4 将内核卸载至英特尔® MIC 架构。

多协处理器支持

在真实的工作负载中,使用多个英特尔至强融核协处理器可以进一步提升性能。我们已经在主从模式下,使用英特尔 MPI 库针对炮点级并行化实现了 ICluster SSF PsDM 代码。因此,凭借英特尔® Composer XE 2013 中包含的多协处理器支持,可轻松实现使用两个或多个英特尔® 至强融核协处理器。在图 5 中,我们将展示如何使用简单表达式来卸载到特定协处理器:(myid-1)%M + 1,其中 M 是协处理器的总数,myid MPI 进程数。

此外,为了开发使用 CPU 和英特尔® 至强融核协处理器资源的异构代码,我们创建了一个布尔变量 use_mic,以控制特定 MPI 进程中调用的内核应该卸载还是应该在 CPU 上。


!dec$ offload if (use_mic) target(mic: myid)
     +     in (…)   // data needs to be copied in 
     +     inout (…)  // data needs to copied in and out 
        CALL SSF_PSDM_MIC_WRAPPER (…)

图5 SSF PsDM 代码的多协处理器支持。

在下面的测试中,我们将设置一台配有两个英特尔® 至强融核协处理器、基于一个双插槽英特尔® 至强® 处理器 E5-2680 的服务器。我们实施了一个完全混合的异构模型,以便 SSF PsDM 代码能够在 CPU 和协处理器上同时运行。具体而言,我们创建了 4 个 MPI 进程:1 个主进程和 3 个相关的工作进程。其中一个工作进程被绑定到 2 个 CPU 插槽上,另外两个工作进程分别被卸载到 2 个协处理器上。我们为 CPU 上运行的 MPI 进程创建了 12 条线程,并在每个协处理器上创建了 93 条线程。

英特尔® 至强融核协处理器包含多达 61 个内核,每个内核最多可支持 4 个硬件线程。在初步研究中,我们发现每个内核使用 3 条线程可实现最佳性能。

线程关联

线程关联当然有用,以便尽可能减少两个处理器之间的迁移,从而提高高速缓存效率。在 CPU 端,我们使用KMP_AFFINITY=compact OpenMP* 线程绑定到相同插槽上,以便在双插槽的英特尔至强处理器服务器上启动 2 个进程。在英特尔® 至强融核协处理器端,我们使用“平衡”模式(见图 6)(而非“紧凑”模式)来确保所有线程均匀分布于所有协处理器内核上,只有英特尔® 集成众核架构上的 OpenMP* 运行时才能识别该模式。

实际上,“分散”模式也能够保持良好运行,因为与每个计算线程相关的数据是独立的。否则,如果有相邻线程在相同的数据上工作,那么在“平衡”模式下的性能更佳,因为线程可以共享高速缓存。或者,KMP_AFFINITY=proclist=[…] 可用于显式设置线程号。有关面向 OpenMP* 线程的处理器关联性的更多详细信息,请参阅 OpenMP* 实施文档 [3]


call kmp_set_defaults("KMP_LIBRARY=balanced")

图6 英特尔 MIC 架构上的线程关联性

性能评估和特性

性能对比

如表 1 所示,单个英特尔® 至强融核协处理器在双插槽的 16 内核英特尔至强处理器服务器(基准代码)上实现了 2.2 倍的性能提升。通过利用完全异构模式(无疑是同时利用 CPU 和协处理器的最佳模式,带有通过英特尔® 至强处理器 E5-26802 CPU)和两个英特尔® 至强融核协处理器配置的两个节点),我们实现了比基于英特尔至强处理器的服务器基准性能高 10.3 倍的性能。

表1 使用英特尔® 至强® 处理器和英特尔® MIC 架构的性能结果

否定结构工作负载

英特尔® 至强® 处理器 E5-2680 2.7Ghz 基准性能(2 CPU– 16 线程

英特尔至强处理器 E5-2680 2.7Ghz 优化性能(2 CPU– 16 线程

1 个英特尔® 至强融核协处理器卸载性能 – 180 线程

两个节点(4 个英特尔至强处理器 E5-2680 + 4 个英特尔® 至强融核协处理器)

花费的时间

2802

1342

1273

272

加速

1

2.08

2.2

10.3

注:在异构测试中,我们使用了统一的优化代码。

我们还评估了英特尔® 至强融核协处理器上的内核扩展情况(我们在每个内核上锁定了三条线程),如图 7 所示。整体可扩展性非常好 我们在这个工作负载的 60 个内核中实现了 42 倍的性能提升。

图7 英特尔® 至强融核协处理器上的线程可扩展性

性能分析

为了更好地了解性能特征,我们使用了 VTune™ Amplifier XE 2013 来分析英特尔至强融核协处理器端的代码。英特尔® 至强处理器 E5-2680 和英特尔® 至强融核协处理器的热点如下面的表 2 所示:

表2 性能配置文件常用函数

模块

函数

英特尔® 至强® 处理器 E5-2680 2.7Ghz 基准

英特尔® 至强融核™ 协处理器

备注

英特尔MKL

DFT

31.4%

37.3%

大小为 512x256 2D C2C

sin/cos/vccis

11.7%

21.5%

VML 的矢量化实施

PsDM

phase_shift

5.3%

11.7%

未在任一平台上执行矢量化处理。大量分支。

ssf_test

2.9%

1.7%

复杂运算,矢量化

imaging

32.1%

3.2%

英特尔 MIC 架构上的手工编码优化

libiomp

 

6.8%

16.2%

线程开销

总计

 

90.2%

91.6%

 

首先,为了使用英特尔® 至强融核协处理器上的 512 SIMD 引擎并获得最佳性能,我们通过以下三种方式对热循环中的大部分计算进行了矢量化处理:1) 通过调用英特尔® MKL 函数的优化实施(例如英特尔® MKL DFTsin/cos/vccis 和英特尔® MKL VML),2) 通过使用英特尔® Composer XE 2013 自动矢量化,在 ssf_test 函数中,我们使用“fno-alias”选项将“没有失真且成功地对函数进行了进一步的矢量化处理”信息告诉运行中的英特尔® Composer XE 20133) 通过手工编码优化,例如面向波场成像和叠加的成像功能。我们重新构建了循环,以便对本来包含许多随机内存读取和许多复杂数据类型运算的内部循环进行矢量化处理。通过分析,我们发现 32.1% 的计算时间用在了基准代码上,优化之后这一时间被显著缩短。图 8 显示了原始版本的基准代码。


! baseline code 
!$omp parallel do private(ix,iy,iw)
do iy = 1, NY
 do ix = 1, NX
  do iw = 1, NW
    rimag(ix,iy)=rimag(ix,iy)+dat_r(ix,iy,iw)*CONJG(dat_s(ix,iy,iw))
  enddo
 enddo
enddo 

图8 基准成像代码

由于频繁的高速缓存缺失,原始代码的性能很差。在这种情况下,优化策略是非常简单的 交换循环,以确保内部循环(IX 循环)的内存访问具有连续性。这还可以让英特尔® Composer XE 2013 基于编译报告对循环进行自动矢量化处理。图 9 显示了优化代码。


! optimized code 
!$omp parallel do private(ix,iy,iw)
do iy = 1, NY
 do iw = 1, NW
  do ix = 1, NX
    rimag(ix,iy)=rimag(ix,iy)+dat_r(ix,iy,iw)*CONJG(dat_s(ix,iy,iw))
  enddo
 enddo
enddo

图9 优化成像代码

将 IW 循环作为最外侧的循环并对其进行并行化处理并不会取得更好的结果,因为在这种情况下,我们还需要添加一个原子操作来保护对 rimag(ix,iy) 的调用,这样代码的质量将远远低于图 9 中的代码。

其次,正如我们在表 2 中所看到的,英特尔 MKL 中的二维复杂对复杂 (C2C) DFT sin/cos/vccis 函数是这一应用中最热门的函数,占到英特尔® 至强融核协处理器计算时间的 58.8%phase_shift 函数占计算时间的 11.7%,由于大量分支和间接内存访问而没有在任一平台上进行矢量化处理。

模块名称

60 条线程

120 条线程

180 条线程

240 条线程

Libiomp5.so

12.1%

12.2%

16.2%

26.36%

表3 英特尔 MIC 架构上的线程开销

最后我们还发现了一个随着线程数增加而产生的负载平衡问题,如表 4 所示。根据 VTune™ Amplifier XE 2013 分析结果,这一问题与 __kmp_wait_sleep() 函数有关。具体而言,这个问题是由于工作共享结构端的隐式障碍引起的,因为不同的线程号无法均匀地分配我们并行化处理的频域循环的运行次数。不平衡部分的比例将随着线程数量的增加而上升。调度策略,比如 OpenMP* 中的动态调度或引导调度,在其它情况下可能会有一些好处。但对于这种情况,实验表明静态调度策略能够实现最佳性能,原因在于动态调度和引导调度所涉及的开销。

总结

在这份报告中,我们评估了第一代英特尔至强融核协处理器上的分步傅里叶 PsDM 性能,通过一台基于双插槽英特尔® 至强处理器 E5-2680 的服务器(主频为 2.7Ghz)实现了 2.2 倍的性能加速。我们还对英特尔® 至强处理器和英特尔® 至强融核协处理器上的多协处理器编程和异构计算模型进行了评估,在配有两个英特尔® 至强融核协处理器、基于双插槽至强处理器 E5-2680 的服务器的两个节点上实现了比基准性能高 10.3 倍的性能。在这种情况下,我们只添加了 8 行 Fortran 代码,以卸载 PsDM 内核并管理英特尔® 至强融核协处理器的数据流。更重要的是,英特尔® 至强融核协处理器上执行的优化(如交换 imaging() 函数中的嵌套循环并从英特尔® MKL 库中调用 sin/cos/vccis 例程)也适用于 CPUCPU 上的 PsDM 性能也提高了 2.08 倍。

与此同时,我们还进行了广泛的性能分析,包括可扩展性和性能分析。分析结果显示,内核数方面实现了近线性的性能扩展(在每个内核有 3 条线程的 60 个内核上实现了 42 倍的性能加速)。

参考

[1] http://software.intel.com/mic-developer
[2]
http://www.openmp.org/
[3] 2011
北京 IDF课程 STFS009通过英特尔® 多核处理器和英特尔® 集成众核架构实现统一的异构编程方法
[4]Virieux J., Calandra H., Plessix R.E
面向地球物理成像的光谱、伪谱、有限差分和有限元建模技术回顾。地球物理勘探法,2011, 59, 794–813

系统配置

·         英特尔 MIC 架构软件开发平台(Shady Cove

·         英特尔测试级软件(英特尔 MIC 编译器、驱动程序等)

硬件规格

2 个英特尔® 至强融核协处理器

ES2 Si1.1Ghz),8 GB GDDR55.5GT/s

主机:英特尔® 至强® E5u

一个配有 2 个英特尔至强处理器 X26802.6 Ghz8 内核,20 MB 三级高速缓存)和 48 GB DDR3(主频为 1333 MHz)的双插槽平台。操作系统:RHEL 6.3

英特尔 MIC 架构软件堆栰

英特尔 MIC 架构内核驱动程序,2.1.3653-8 版本(测试版)
Flash
映像/uOS 2.1.02.0314/2.6.34.11-g4af9302

英特尔® Composer XE 2013 Beta VTune™ Amplifier XE 2013 Beta. Mpich1.2.7

声明

本文件中包含关于英特尔产品的信息。本文件不构成对任何知识产权的授权,包括明示的、暗示的,也无论是基于禁止反言的原则或其他。除英特尔产品销售的条款和条件规定的责任外,英特尔不承担任何其他责任。英特尔在此作出免责声明:本文件不构成英特尔关于其产品的使用和/或销售的任何明示或暗示的保证,包括不就其产品的(i)对某一特定用途的适用性、(ii)适销性以及(iii)对任何专利、版权或其他知识产权的侵害的承担任何责任或作出任何担保。

“关键业务应用”是指当英特尔产品发生故障时,可能会直接或间接地造成人员伤害或死亡的应用。如果您针对此类关键业务应用购买或使用英特尔产品,您应当对英特尔进行赔偿,保证因使用此类关键业务应用而造成的产品责任、人员伤害或死亡索赔中直接或间接发生的所有索赔成本、损坏、费用以及合理的律师费不会对英特尔及其子公司、分包商和分支机构,以及相关的董事、管理人员和员工造成损害,无论英特尔及其分包商在英特尔产品或其任何部件的设计、制造或警示环节是否出现疏忽大意的情况。

英特尔有权随时更改产品的规格和描述而毋需发出通知。设计者不应信赖任何英特产品所不具有的特性,设计者亦不应信赖任何标有“保留权利”或“未定义”说明或特性描述。对此,英特尔保留将来对其进行定义的权利,同时,英特尔不应为因其日后更改该等说明或特性描述而产生的冲突和不相容承担任何责任。此处提供的信息可随时改变而毋需通知。请勿根据本文件提供的信息完成一项产品设计。

本文件所描述的产品可能包含使其与宣称的规格不符的设计缺陷或失误。这些缺陷或失误已收录于勘误表中,可索取获得。

在发出订单之前,请联系当地的英特尔营业部或分销商以获取最新的产品规格。索取本文件中或英特尔的其他材料中提的、包含订单号的文件的复印件,可拨打1-800-548-4725,或登陆:http://www.intel.com/design/literature.htm

英特尔、Intel 标识、Vtune、融核、Phi、至强和 Xeon 是英特尔公司在美国和其他国家(地区)的商标。
* 其他的名称和品牌可能是其他所有者的资产
英特尔公司© 2012 年版权所有。所有权保留。

在性能检测过程中涉及的软件及其性能只有在英特尔微处理器的架构下方能得到优化。诸如SYSmark和MobileMark等测试均系基于特定计算机系统、硬件、软件、操作系统及功能, 上述任何要素的变动都有可能导致测试结果的变化。请参考其他信息及性能测试(包括结合其他产品使用时的运行性能)以对目标产品进行全面评估。

英特尔处理器标号不是性能的指标。处理器标号仅用于区分同属一个系列的处理器的特性,而不能够用于区分不同系列的处理器。详情敬请登陆 了解有关英特尔® 处理器号的信息

有关编译器优化的更完整信息,请参阅优化通知