基于MPI+OpenMP的时间域航空电磁快速正演算法

doi:10.11720/wtyht.2020.1460

基于MPI+OpenMP的时间域航空电磁快速正演算法

任运通, 李貅, 齐彦福, 曹华科

长安大学地质工程与测绘学院,陕西西安 710054

Research on time domain airborne electromagnetic fast forward algorithm based on MPI+OpenMP

REN Yun-Tong, LI Xiu, QI Yan-Fu, CAO Hua-Ke

College of Geological Engineering and Geomatics, Chang'an University, Xi'an 710054, China

通讯作者: 李貅(1958-),男,教授、博导,主要从事瞬变电磁场的理论与应用方面的研究工作

责任编辑: 沈效群

收稿日期: 2019-09-24 修回日期: 2019-12-3 网络出版日期: 2020-04-20

基金资助:

国家重点研发计划项目. 2017YFC1502605
国家自然科学基金重点项目. 41830101

Received: 2019-09-24 Revised: 2019-12-3 Online: 2020-04-20

作者简介 About authors

任运通(1995-),男,长安大学硕士,研究方向为瞬变电磁探测。

摘要

时间域有限元算法已被广泛应用于航空电磁三维正演模拟当中,然而由于航空电磁测区面积大,且采样密集,造成正演计算量巨大,传统的串行算法已经无法满足计算效率要求,为此,开展了并行加速算法研究以解决计算效率不足的问题。基于航空电磁系统的影响范围有限,采用局部网格技术将计算任务划分成多个子网格,即每个发射源一套网格,各网格的正演计算相互独立,不存在数据依赖性,具有很好的可并行性;利用MPI技术对多个子网格正演任务进行分配,在各个进程上进行并行计算;针对每个正演子网格,在进行时间域有限元算法正演模拟过程中,采用OpenMP技术对单元矩阵进行并行计算。典型地电模型的数值模拟结果表明本文开发的MPI+OpenMP并行正演算法可以有效提高正演速度,最高加速比可达10倍。

关键词： 时间域航空电磁 ; MPI ; OpenMP ; 有限元

Abstract

The time domain finite element algorithm has been widely used in airborne electromagnetic three-dimensional forward modeling. However, due to the large airborne electromagnetic measurement area and dense sampling, the forward calculation is huge, and the traditional serial algorithm cannot meet the calculation efficiency requirements. In view of such a situation, the authors carried out parallel acceleration algorithm research to solve the problem of insufficient computational efficiency. Based on the limited range of influence of aviation electromagnetic system, the local grid technology is used to divide the computing task into multiple sub-grids, that is, a set of grids for each source, and the forward calculations of each grid are independent of each other, and hence no data dependency is existent and there is good parallelism. In this paper, MPI technology is used to allocate multiple sub-grid forward tasks, and parallel computing is performed on each process. For each forward subgrid, in the forward modeling of the time domain finite element algorithm, the parallel matrix is calculated by OpenMP technology. The numerical simulation results of the typical electrical model show that the MPI+OpenMP parallel forward algorithm developed in this paper can effectively improve the forward speed, and the maximum acceleration ratio can reach 10 times.

Keywords： time-domain airborne EM ; MPI ; OpenMP ; finite element

PDF (5724KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

任运通, 李貅, 齐彦福, 曹华科. 基于MPI+OpenMP的时间域航空电磁快速正演算法. 物探与化探[J], 2020, 44(2): 290-299 doi:10.11720/wtyht.2020.1460

REN Yun-Tong, LI Xiu, QI Yan-Fu, CAO Hua-Ke. Research on time domain airborne electromagnetic fast forward algorithm based on MPI+OpenMP. Geophysical and Geochemical Exploration[J], 2020, 44(2): 290-299 doi:10.11720/wtyht.2020.1460

0 引言

时间域航空电磁法采用机载移动平台,通过线圈发射大功率一次场电磁波信号,利用接收线圈接收到感应的二次磁场,实现对地下电性结构及能源矿产分布的精细探测。该方法具有分辨率高、探测速度快、通行性极好等诸多优点,已被广泛应用于油气资源、环境工程和金属矿产等领域^{[1,2,3,4,5,6]}。然而,传统的电阻率成像和一维反演技术已经无法满足三维问题解释的要求,开展三维反演解释技术研究是当前的热点问题,作为反演的基础,时间域航空电磁三维正演算法受到越来越多的关注。

目前,时间域航空电磁三维有限元正演方法主要分为间接法和直接法两种。其中的间接法是基于频率域麦克斯韦方程组,首先计算频率域响应,再通过时频转换方法转换到时间域^[7,8,9,10]。由于间接算法存在稳定性差等诸多问题,人们普遍将研究方向转移到了时间域直接算法,直接求解时间域电磁响应。2010 年,Um等^[11]采用非结构化网格剖分策略,对电性源的海洋三维瞬变电磁场进行了时间域矢量有限元的三维正演模拟计算。2016年,李贺等^[12]人研究了直接时间域矢量有限元瞬变电磁三维正演模拟方法。2017年,齐彦福等人^[13,14]基于航空电磁footprint的局部网格策略,利用非结构矢量有限元方法,对复杂介质情况下的时间域三维航空电磁响应进行了进一步模拟研究。然而,由于航空电磁法采用同时移动发射源和测点的方式进行观测,且采样密集,每次移动测点均需求解一次正演方程,产生巨大的计算量,传统的串行计算方法无法满足未来三维反演的效率要求。考虑到并行计算可以有效提高程序的运算效率,解决相同数量的问题所需执行的时间更短,开展了基于MPI+OpenMP并行技术的时间域航空电磁快速正演算法研究。

目前广为使用的并行加速技术有两种,分别是MPI(Message Passing Interface)和OpenMP(Open Multi-Processing)并行技术。MPI和OpenMP均需要与计算机语言结合使用,因两者较好的通信性和可移植性而受到广泛地应用。在尺度上,MPI多为基于多台处理机上的跨节点并行方法,OpenMP主要通过一些指令集对现有的Fortran或C/C++等程序进行扩展。并行技术在提高计算效率上具有显著优点,将该方法应用到地球物理数值模拟和数据处理中具有良好的前景。早在1997年,Newman和Alumbaugh^[15]对三维大地电磁正演方法以及灵敏度矩阵进行了深入研究,成功将并行加速技术应用到了电磁反演当中,实现了多个节点的并行加速计算,极大推动了并行计算技术在地球物理领域的发展。国内对于并行计算技术应用到地球物理的研究,起步虽然稍晚于国外,发展却并不逊色。2005年荣莹等^[16]实现了基于MPI的处理机集群并行计算系统平台的构建,在现有的硬件条件下提高计算力; 2006年,谭捍东^[17]等结合MPI的优越性,通过频点并行的方式成功实现了大地电磁三维正演的并行计算,得到了很高的加速比,验证了并行算法的稳定性以及高效性;2011年,李小康^[18]研究了频率域航空电磁法有限单元二维正演的并行计算;2015年,陈辉^[19]研究了基于MPI的航空瞬变电磁一维正反演,提高了航空电磁一维模拟技术的计算效率。目前,对于并行加速技术在电磁方法数值模拟中的应用主要集中在大地电磁法,而航空电磁正演的并行算法研究依然停留在一维正反演和二维的频率域正演阶段。

本文将并行技术与时间域有限元算法相结合,基于三维航空电磁局部网格之间相互独立的特性^[20],通过采用MPI+OpenMP的并行加速策略,在保证正演结果可靠性的基础上,极大提高了时间域航空电磁有限元三维正演的计算效率。首先通过与串行有限元法和有限体积法进行对比检验本文并行程序的可靠性,然后进行计算效率分析,讨论并行加速比与硬件条件的关系,最后将该并行算法应用于复杂起伏地表模型,模拟其航空电磁响应。

1 方法理论

1.1 时间域航空电磁有限元正演理论

时间域麦克斯韦方程组可以表示为^[21]

(1)

\nabla \times e (r, t) = - \frac{\partial b (r, t)}{\partial t}

(2)

\nabla \times \frac{1}{μ} b (r, t) = j (r, t) + \frac{\partial d (r, t)}{\partial t}

(3)

\nabla \cdot d (r, t) = q

(4)

\nabla \cdot b (r, t) = 0

其中,t为时间,r为位置矢量,e(r,t)、j(r,t)、d(r,t)及b(r,t)分别表示r处在t时刻电场强度、电流密度、电位移矢量以及磁感应强度,q和μ分别表示累积电荷和磁导率。本文假定介质的磁导率与自由空间磁导率相同,即μ=μ₀=4π×10^-7 H/m,介电常数与自由空间介电常数ε₀相同,即ε=ε₀=8.85×10^-12 F/m。电磁场之间具有如下本构关系:

(5)

j (r, t) = σe (r, t) + j_{s} (r, t),

(6)

d (r, t) = εe (r, t),

(7)

b (r, t) = μh (r, t),

其中:j_s(r,t)表示外部施加的源电流密度,h(r,t)和σ分别表示磁场强度和电导率。通过消去磁场,可以获得电场扩散方程

(8)

\frac{1}{μ} \nabla \times \nabla \times e (r, t) + σ \frac{\partial e (r, t)}{\partial t} + \frac{\partial j_{s} (r, t)}{\partial t} = 0 。

本文采用非结构矢量有限元方法进行空间离散,四面体单元中任意位置的电场可以表示为:

(9)

e^{k} (r, t) = \overset{6}{\sum_{i = 1}} e_{i}^{k} (t) n_{i}^{k} (r),

其中:e^k(r,t)是第k个单元第i条棱边上的电场值; $n_{i}^{k}$ (r)是矢量插值基函数^[22],其表达形式为:

(10)

\begin{matrix} n_{i}^{k} (r) = (L_{i_{1}}^{k} \nabla L_{i_{2}}^{k} - L_{i_{2}}^{k} \nabla L_{i_{1}}^{k}) {l_{i}^{k}}_{} \end{matrix},

$l_{i}^{k}$ 是第k个单元第i条棱边上的长度, $L_{i_{1}}^{k}$ 和 $L_{i_{2}}^{k}$ 是第i条棱边上的两个节点i₁和i₂的标量插值基函数^[23]。通过伽辽金方法可以获取有限元控制方程:

(11)

M \frac{d e (t)}{d t} + Se (t) + J = 0,

其中:M、S分别为质量和刚度矩阵,J为电流源项,针对每个单元,可由下面三式给出

(12)

M_{i, j}^{k} = \underset{V^{k}}{\int \int \int} σ^{k} n_{i}^{k} (r) \cdot n_{j}^{k} (r) d V,

(13)

S_{i, j}^{k} = \frac{1}{μ} \underset{V^{k}}{\int \int \int} \nabla \times n_{i}^{k} (r) \cdot \nabla \times n_{j}^{k} (r) d V,

(14)

J_{i}^{k} = \underset{V^{k}}{\int \int \int} n_{i}^{k} (r) \cdot \frac{\partial j_{s} (r, t)}{\partial t} d V,

其中V是单元的体积。由于各个单元相互独立,因此M^k、S^k和J^k可以采用OpenMP并行技术加速计算。针对发射源,将发射线圈分解为若干段导线,每段近似为一个电偶极子^[24]每个电偶极子可以表示为:

(15)

j_{s} (r, t) = δ (r - r_{s}) vI (t) d l 。

上式中:I(t)表示t时刻的电流强度,dl表示电偶极子长度,v表示电流方向,δ表示脉冲函数。然后,利用二阶后推欧拉公式对有限元方程(11)进行时间域离散:

(16)

(3 A + 2 Δ tB) e^{i + 2} (t) = A [4 e^{i + 1} (t) - e^{i} (t)] - 2 Δ t S^{i + 2},

其中:Δt为时间步长,eⁱ(t)表示第i时刻的电场值。式(16)亦可简写为:

(17)

Fe = P,

其中:F表示大型稀疏系数矩阵,P是右端项。当发射阶跃电流波形时,在0时刻之前发射线圈中供恒定电流,在全空间产生稳定的磁场,根据楞次定律可知空间中任意位置的电场均为0,故电场的初始条件为:

(18)

e (r, 0) = 0 。

采用狄利克雷边界条件,根据电磁波在空间中的几何指数衰减规律,可以通过加大计算区域的扩边范围,使其满足

(19)

e |_{Γ} = 0,

其中Γ表示计算区域的外边界。最后,采用直接求解器Pardiso对线性方程组进行求解,即可获得全空间的电场值,再利用法拉第电磁感应定律(式(1))计算磁场响应。

1.2 局部网格正演计算技术

大量的数值实验表明三维数值模拟结果的精度受到网格质量的影响非常大,因此需要对导电大地和发射源进行精细的网格剖分来获得高精度正演结果。传统的时间域航空电磁三维正演模拟算法采用全局网格加密方式,即通过一套网格对所有测点的电磁响应进行模拟。然而,由于航空电磁法采样密集且观测剖面长,导致正演网格单元数量巨大。如果采用全局网格进行正演模拟,将产生巨大的计算量,严重降低计算效率。考虑到航空电磁系统影响范围有限(如图1所示),本文采用局部网格加密技术,针对每个测点或者相邻的几个观测点分别设计独立的网格,并分别在局部网格上进行正演模拟,以此在保证结果可靠性的同时提高计算效率。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 航空系统影响范围示意

Fig.1 Moving footprint for time-domian airborne EM system

局部网格技术根据航空电磁系统的影响范围设计合理的局部精细网格。图2展示了全局网格和局部网格,其中图2a表示全局网格,图2b、c、d则分别表示针对测线上不同测点所设计的相互独立的局部网格。从图中可知全局网格对测线上所有测点、测线下方地空分界面以及异常体部分均进行了相对细致的网格剖分,因此造成网格数量较大,导致正演模拟的计算量剧增,计算效率降低。而局部网格仅仅对异常体、当前测点位置以及其下方的地空分界面进行局部加密,对单个测点进行单独求解,因此所使用的网格数量大大降低,计算量减小。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 全局网格与局部网格对比

(a为全局网格;b、c、d为局部网格)

Fig.2 Global grid vs. local grid

((a) is a Global grid; (b) (c) (d) is a local grid)

2 并行加速

2.1 基于OpenMP的并行加速

OpenMP采用分叉—合并(Fork-Join)执行模式。一个OpenMP程序开始于一个单一的线程,该进程又称作主线程,通过并行指令对程序的并行区间进行定义,在这个区间中程序块由多个线程自动分配任务及并行执行,线程数由主线程决定,在一个大型程序中可以嵌套多个OpenMP并行区间,其中并行区间由一组$OMP指令进行开启(!OMP PARALLEL)和关闭(!OMP END PARALLEL)。图3是OpenMP并行模式的示意图。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 OpenMP并行模式示意

Fig.3 OpenMP parallel mode schematic

2.2 基于MPI的并行加速

MPI需要用头文件use mpi进行声明,然后在主进程(即节点)中进行初始化,开启MPI并行环境,获取进程个数及编号,再将多个任务合理地分配到各个进程并传输对应的数据,然后分别进行计算,最后汇集到主进程并关闭并行环境,如图4所示。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 MPI并行模式示意

Fig.4 MPI parallel mode schematic

2.3 并行程序设计及实现

由于航空电磁系统在每个观测点处的灵敏区域相对于总计算区域来说远远不及,因此仅仅在观测点的一定位置范围内对网格进行加密即可满足计算的精度要求,针对每个观测点分别单独设计独立的网格来提高计算效率。

时间域航空电磁三维正演各个测点计算相互独立,不存在数据依赖关系,具有非常好的并行性,在三维正演的基础上实现了基于MPI的多测点并行计算,在独立网格内实现了基于OpenMP的单元矩阵并行计算。程序采用主从模式,整个程序的基本结构及运行流程由主进程把控,主要负责读取模型参数、数据传输以及对并行任务的分配,子进程负责接收来自主进程传递的消息、对分配的任务进行计算以及输出本进程的计算结果。为了充分利用计算资源提高计算效率,主进程在完成任务分配和数据传输的工作后也参与到测点的并行计算当中。

具体实现的流程如图5所示,主进程读取模型的整套网格参数,随后将其传递给其他子进程,之后各进程根据分配的网格参数计算模型的响应结果,为了防止数据传输时发生通信冲突,采用各进程分别输出各自的计算结果。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 航空瞬变电磁三维正演并行计算流程

Fig.5 Airborne transient electromagnetic three-dimensional Forward modeling parallel computing flow chart

3 模型算例

3.1 精度验证

为检验本文所开发的时间域航空电磁多粒度并行正演方法的正确性,设计了如图6所示的水平板状体模型,发射波形为阶跃波,发射电流强度为435 A,通过全局网格和局部网格两种剖分方式,对异常体上方y=0 m剖面上-400~400 m范围内的32个测点,分别进行正演模拟,且将正演结果与有限体积方法的正演结果进行比对,进行数值精度验证。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 水平板状体模型示意

Fig.6 Schematic diagram of the horizontal plate model

使用的集群服务器有16个型号为Intel(R) Xeon(R) CPU E5-2609 v4 @1.70 GHz的逻辑CPU,每个CPU拥有2G内存和4个线程。采用时间域有限元全局网格法、基于多粒度并行加速策略的时间域有限元局部网格法以及有限体积方法分别计算的电磁响应结果如图7所示,可以清楚地看出,上述3种不同方法所计算的结果吻合得非常好,且均对模型中的地下板状体有明显反映。这一结果有效验证了本文并行算法的精度。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 全局网格串行计算和局部网格并行计算电磁响应结果对比

Fig.7 Comparison of electromagnetic response results between global grid serial computing and local grid parallel computing

3.2 并行性能分析

为了分析基于MPI+OpenMP并行加速策略的时间域航空电磁三维正演程序的并行加速效果,采用上述水平块状体模型进行测算,根据控制变量法,控制计算测点数为32和线程数为4保持不变,依次对不同数量的进程进行比较,通过加速比和并行效率两个参数来评估本文设计的三维正演并行程序。

令原有的串行程序在集群服务器单进程上的运行时间为T_s,经过多粒度并行优化后,运行所需时间为T_p,其中P表示开启的进程数,则加速比S_p可表示为:

(20)

S_{p} = T_{s} / T_{p},

并行效率可定义为:

(21)

e_{p} = S_{p} / P 。

测试结果见表1。在计算测点数一定时,随着进程数的增加,计算时间呈非线性下降趋势,加速比逐渐增大,相反并行效率逐渐降低;当开启16个进程时,采用多粒度并行加速后执行耗费的时间仅仅是串行程序耗时的1/10,同时并行效率达到最低。

表1 基于MPI+OpenMP并行加速策略的行正演计算统计

Table 1 Line forward calculation statistics based on MPI+OpenMP parallel acceleration strategy

程序类型	进程数	计算时间/s	加速比	并行效率
串行	1	2136.4303	1
	1	2139.7464	0.998	99.85%
	3	745.2973	2.867	95.55%
并行	5	499.5326	4.277	85.54%
	7	372.3317	5.738	81.97%
	12	266.8348	8.007	66.72%
	16	210.9645	10.127	63.29%

新窗口打开| 下载CSV

由图8可以发现以下特点:并行化后的正演程序加速比与开启的进程数并不是整数倍关系,而是呈非线性增长,且并行效率整体呈降低趋势。这是因为在集群服务器上开启并行环境后,会占用一定内存,且进程间进行数据传递需要占用时间,而且任务量与开启的进程数不匹配时,计算任务少的进程会率先完成计算量,但该进程并不会结束,而是要等待未完成任务的进程,从而导致时间的损耗。所以开启P个进程,加速比S_p并不能达到P,且随着进程数的增加S_p与P的差值越来越大。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 不同进程数的加速比(a)和并行效率(b)

Fig.8 Acceleration ratio (a) and parallel efficiency (b) for different node numbers

3.3 复杂模型计算

考虑到地形起伏不平对航空电磁正演结果的影响较为严重,采取四面体网格进行加密剖分后的计算量较大,从而导致计算时间较长,因此利用多粒度并行加速技术进行加速计算。导入地形文件数据模拟起伏地表,并设计地下埋藏有块状及倾斜板状良导体(图9)。设高阻围岩的电导率为0.01 S/m,块状体和板状体的电导率为1 S/m;倾斜版状体顶部埋深60 m,垂直深度250 m,块状体顶部埋深50 m,边长200 m。采用中心回线装置,飞行高度30 m,发射线框半径15 m,线圈匝数为1,采用阶跃波激发。全区设置11条测线,每条测线布设171个测点,共计1 881个测点,每个网格单元数约90 000个,正演过程共需进行9 350次矩阵分解,327 250次回代。若采用传统串行计算时,总耗时约1 968.623 min(32.8 h),采用并行加速开启16个进程后总计算时间约194.336 min(3.24 h),消耗内存约5.2 G,效率提高了10.13倍,如表2所示。

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 复杂模型计算

Fig.9 Calculation of complex models

表2 计算情况统计

Table 2 Calculation statistics table

	网格单元数	消耗内存	分解次数		回代次数		计算总耗时/h
串行并行	约90 000	约5 350M 5.2G	9 350	327 250		32.8 3.24

新窗口打开| 下载CSV

图10是主剖面(x=0测线)的多测道图,可以看出在早期表现出地形的响应,到了晚期,地形影响逐渐变小,表现出异常体的响应。图11呈现出整个测区在四个时间点时的航空电磁响应结果,从图中可以明显看出在时间早期,仅存在地形的响应,由于地形凸起,电磁响应呈现出相对低异常,在x=0,y=-500附近,相对低异常达到极大值,由于块状异常体处于凸起地形下方,相对于倾斜板埋藏较深,故随着时间延长,倾斜板的异常响应率先显示出来,随后出现块状体异常响应,且由于地下埋藏有异常体,导致地形的影响减弱,到晚期时,电磁波穿过异常体, 仅剩地形的电磁响应。

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 复杂模型主剖面(x=0测线)多测道曲线

Fig.10 Multiple trajectory map of complex model main section (x=0 line)

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 复杂地形响应曲面

Fig.11 Complex terrain anomaly response surface map

4 结论

利用MPI实现了三维正演模拟中的多测点并行计算和实现了独立网格内单元矩阵计算的OpenMP并行加速,最终实现了基于MPI+OpenMP并行加速策略对时间域航空电磁三维正演方法的并行化,并得出并行化后的正演程序,其加速比随着开启的进程数增加呈非线性增长。

由于“并行开销”的存在,开启的进程数越多,并行效率整体呈现出下降趋势,即使只开启一个进程,效率也比串行效率低。地形起伏模型的数值模拟体现了本文并行优化后的正演方法的高效性,正确高效的三维正演为三维反演提供了可能。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Pemberton R

Airborne electromagnetics in review

[J]. Geophysics, 1962,27(5):691-713.