基于GPU的航空瞬变电磁法响应模拟

沈 磊， 董春晖,2

(1.山东交通职业学院 公路与建筑学院，潍坊261206； 2.大连海事大学 道路与桥梁工程研究所，大连 116026)

基于GPU的航空瞬变电磁法响应模拟

沈 磊1， 董春晖1,2

(1.山东交通职业学院 公路与建筑学院，潍坊261206； 2.大连海事大学 道路与桥梁工程研究所，大连 116026)

针对航空瞬变电磁法正演算法结构与GPU的并行性运算相结合，推导出了航空瞬变电磁法一维正演公式，在2007年NVIDIA公司推出的CUDA平台上实现了一维正演算法在GPU上的调用。为验证CUDA C语言程序的正确性，设计了几个模型，并与C代码模拟的响应做了对比，结果表明，基于GPU的CUDA C代码所编的程序在航空瞬变电磁法的模拟中拥有明显优势，速度提高非常明显，在复杂的地球物理学模型计算和反演中，拥有巨大的发展潜力。

GPU； CUDA； 航空瞬变电磁法； 一维正演

0 引言

随着地球物理勘探工作的深度和广度不断提高，数据处理结果的精度要求也越来越高。巨大的数据量和较高的数据处理精度对计算机处理技术提出了越来越高的要求。提高计算精度和效率，吸引了更多研究者地注意。近年来，由于计算机图形处理器(Graphic Processing Unit，GPU)成本和功耗比较低且并行处理能力强大的特点，使它在通用计算方面得到了极大地发展[1]。

2007年，NVIDIA公司发布了一种用于GPU通用运算的C扩展语言，被称为CUDA(Compute Unified Device Architecture)平台。CUDA是一种不需要借助图形学API(Application Programming Interface，应用程序接口)，就可以使用类C语言进行GPU通用计算的开发环境和软件体系。

航空瞬变电磁法是一种使用直升机或固定翼飞机作为机载工具，利用物质电导率(或电阻率)的差异来探测地球电磁响应的非常实用的地球物理探测方法，由于它拥有快速、高效的扫描广大区域的能力，且空间密度大大高于典型的地面电磁法，所以很受欢迎。

笔者编写了基于GPU的航空瞬变电磁法一维正演模拟程序，并且证明了它是高效、可行的。

1 CUDA C程序

CUDA C与C语言非常类似[2]，仅仅是增加了少量说明符。在用CUDA C编写的程序中， CPU为“主机(host)”,GPU为“装置(device)”。在使用GPU时，我们编写了能在GPU上运行的数据转换代码和函数定义。首先分配GPU上的内存，数据要在GPU上成功运算要传送到GPU。GPU上的计算通过一个定义了特殊说明“__global__”的函数传送。这个函数叫做“kernel”函数，然后当“主机”调用它的时候，GPU执行“kernel”函数。

CPU与GPU协同工作，各司其职。CPU负责逻辑性强的事务处理和串行计算(如数据准备和设备初始化)；GPU则专注于执行高线程化的并行处理。CPU与GPU拥有相互独立的存储器地址空间。一旦确定了程序中的并行部分，我们就可以考虑把这部分计算工作交给GPU来执行。本程序模型[3]如图1所示。

图1 CUDA程序模型Fig.1 CUDA programming model

2 航空瞬变电磁一维正演模型

2.1 一维正演公式推导

正演模型的收发系统为中心回线装置。一维地电模型[4]如图2所示。

图2 一维地电模型Fig.2 One dimensional geoelectric model

根据参考文献[5]，推导出的回线中心处磁场的垂直分量表达式为：

(1)

式中：I为源电流；a是发射线圈的半径；ρ是接收线圈半径；rTE是水平层状介质的反射系数；λ是波长；u0是自由空间的磁导率；ui为大地的磁导率，通常我们取大地的磁导率为自由空间磁导率，ui=u0=4π×10-7H/m；h是发射线圈的高度；z为接收线圈距地面的距离,在航空电磁法中，z=-h,模型共有n层;hi是每一层的厚度;σi是地层的电导率。

在准静态条件下：

(2)

(3)

(4)

(5)

μ0=λ

(6)

(7)

由频谱分析理论可以知道，时间域电磁场与频率域电磁场通过傅立叶变换相互关联[6]。

(8)

或

(9)

为了减少横向的不均匀性，实际工作中我们一般测量的是测量磁场的垂直分量的时间导数Hz。

(10)

或

(11)

现在我们得到公式：

(12)

又因为

(13)

(14)

显然，时间域响应是含有贝塞尔函数的双重无穷积分。含有贝塞尔函数的积分为式(15)。

(15)

其中:Jν是ν阶第一类贝塞尔函数;实数ν>-1。则可以进行快速汉克尔变换，成为如下形式[4]。

(16)

式中：ω是快速汉克尔变换系数。本设计中所使用的汉克尔变换系数是根据参考文献[5]求取的，Δ=ln10/10。

笔者发现FHT(快速汉克尔变换)后的双重无穷积分非常适用于并行计算。为此，编写了一个CUDA C语言程序，为了提高数值计算的稳定性，tanh(uihi)使e-uihi用分解运算。使用式(17)变换直接分离实部与虚部。

(17)

2.2 CUDA C程序验证

为验证CUDA C语言程序的正确性和精度，笔者计算了几种典型地电模型的瞬变电磁响应。

由式(15)可知，发射电流强度为一个常数倍数，所以本设计不在考虑发射电流的影响，取发射电流强度I=1A。

2.2.1 不同发射线圈半径的瞬变电磁响应

在均匀半空间模型的条件下，地层电导率设为0.1 s/m，飞行高度为30 m，发射线圈半径不同的条件下垂直方向航空瞬变电磁Hz响应的比较(图3)。取双对数坐标，线圈半径分别为10 m、50 m、100 m，采样时间为0.100 000 0 E-07 s到0.158 500 0 E-01 s取对数平均。

由图3可见，在地电模型相同的情况下，发射线圈半径越小，瞬变电磁响应的衰减越快，在实验所用的10 m～50 m内，由于迅速衰减，在观测时间内，线圈半径越大响应值越大，但在线圈半径为100 m时，与50 m的响应有交叉，可以推断，线圈半径越大，其早期响应值越小。在晚期响应中，曲线显然是趋近于平行的。

2.2.2 不同电导率的瞬变电磁响应特征

依旧采用均匀半空间模型，飞行高度为30 m，发射半径为50 m，模型电导率分别取10 S/m、 0.1 S/m、0.01 S/m。其垂直方向航空瞬变电磁响应如图4所示。采样时间为0.100 000 0E-05 s到0.158 500 0E-01 s取对数平均。

由图4可见，地下电导率越小，进入晚期的双对数坐标下的负斜率直线越早，在早期观测中的响应值越大，反之亦然。晚期曲线趋近于平行。

图3 不同发射线圈半径瞬变电磁响应特征曲线Fig.3 Airborne transient electromagnetic response characteristic curves for different transmitting coil radius

图4 不同电导率瞬变电磁响应特征曲线Fig.4 Airborne transient electromagnetic response characteristic curves for different conductivity

2.2.3 均匀半空间下不同飞行高度的瞬变电磁响应特征

在均匀半空间条件下，发射线圈半径为50 m，飞行高度分别取10 m、30 m、100 m，为了使趋势更加明显，我们增加了采样时间，由0.100 000 0E-04 s到0.158 500 0E+01 s取对数进行采样，垂直方向的航空瞬变电磁响应曲线如图5。

图5 不同飞行高度的瞬变电磁响应特征曲线图Fig.5 Airborne transient electromagnetic response characteristic curves for different flight altitude

由图5可见，不同飞行高度的早期为平行的斜线，但高度越大，早期响应值越小，经过中间的渐变过程后，晚期趋近于同一条负斜率直线，这与使用相同的地电模型是相对应的。

2.2.4 不同参数地电模型的瞬变电磁响应特征。

图6是D型地电模型利用CUDA C语言程序的计算结果。地层介质的导电率分别是0.01 S/m，10 S/m，第一层的厚度是100 m，发射线圈半径为50 m，飞行高度30 m。

图6 D型地电模型航空瞬变电磁响应特征曲线图Fig.6 Airborne transient electromagnetic response characteristic curves for Type D geoelectric model

图7是三层地电模型利用CUDA C语言程序的计算结果。K型地电模型的地层电导率分别为0.1 S/m、0.01 S/m、10 S/m。Q模型的地层电导率分别为0.01 S/m、0.1 S/m、10 S/m。H模型的地层电导率为0.1 S/m、10 S/m、0.01 S/m。地层厚度为100 m,发射线圈半径为50 m,飞行高度30 m。

图7 三层地电模型航空瞬变电磁响应特征曲线图Fig.7 Airborne transient electromagnetic response characteristic curves for three layers geoelectric model

2.3 CUDA C程序与C程序效果对比

为了比较运算性能，我们又用相同的运算方法，使用C语言程序运算了一次。使用当前主流计算机平台，CPU型号AMD Athlon II X2 240；GPU型号Nvidia GeForce 9500 GT，我们设计了以下三个实验。

2.3.1 运算精度对比

在确保GPU上的高性能运算是精确的同时，我们构建一个均匀半空间模型。它的导电率为0.1 S/m.发射线圈半径为50 m，飞行高度为30 m。图8为由C语言源代码双精度浮点数计算值与CUDA C的单精度浮点数计算值比较。

从图8中，我们可以看到两程序数据的结果是一致的，CUDA C程序的精确度是可以信赖的。

图8 均匀半空间航空瞬变电磁响应数值模拟曲线GPU与CPU对比图Fig.8 Comparison chart the airborne transient electromagnetic response characteristic on GPU and CPU

首先定义相对误差，

GPUS为均匀半空间模型在使用GPU运算的响应，CPUS为均匀半空间模型使用CPU运算的响应，则相对误差曲线如图9所示。

图9 均匀半空间航空瞬变电磁响应数值模拟GPU与CPU运算相对误差曲线Fig.9 The relative error curves of GPU and CPU in the numerical simulation of airborne transient electromagnetic response in homogeneous half-space model

可见，在均匀半空间比较计算比较稳定的情况下，计算误差在10-6左右，且相对比较稳定，因此计算精度满足要求。

2.3.2 运算速度对比

(1)不同地电模型的运行时间。在正演运算中，地电模型层数与运算时间直接相关。仅通过改变地电模型层数，我们得到了两程序执行时间如表1所示。Z表示地电模型层数，GT表示GPU执行时间，CT表示CPU执行时间，图10所示为GPU和CPU代码执行时间的对比结果。事实上，结论非常明显，随着地电模型层数的增加，GPU执行时间相对CPU进步一缩短了。

图10 GPU与CPU代码执行时间对比图Fig.10 Execution time comparison chart about CPU program and GPU program

地层数ZGPU耗时GT/msCPU耗时CT/msCT/GT278132817.036125418733.5011171782845.78162191151552.58212651507856.90

2)相同地电模型的执行时间。在航空瞬变电磁法一维正演模拟中，GPU程序和CPU程序执行时间都非常短。执行时间短会造成测量误差，为了减少观测误差，我们做了以下改进。在函数主体语句前加一个代码：for(int dd=0;dd

表2 实验2执行时间数据记录

3 结论

笔者提出了基于GPU高性能运算的CUDA C代码编程的航空瞬变电磁法正演模拟法，通过模型实验，我们确认基于GPU的CUDA C代码所编的程序相比较C语言编程更有优势，同时，数据量越大，GPU耗时明显越少，总而言之，GPU的并行计算能力使它在电磁法数值计算方面应用前景广阔。

[1] 丁科，谭营.GPU通用计算极其在计算智能领域的应用[J].智能系统学报，2015，10(1)：1-11. DING K，TAN Y． A review on general purpose computing on GPU and its applications in computational intelligence[J]．CAAI Transactions on Intelligent Systems，2015，10( 1) : 1-11．(In Chinese)

[2] TAKEHIKO NAWATA，REIJI SUDA. APTCC:Auto Parallelizing Translator From C To CUDA[J].Procedia Computer Science,2011(4):352-361.

[3] 张舒，褚艳利.GPU高性能运算之CUDA[M].北京：中国水利水电出版社，2009. ZHANG S,CHU Y L. CUDA GPU high-performance computing[M].Beijing:China Water & Power Press,2009.(In Chinese)

[4] 朱凯光，林君，刘长胜，等.频率域航空电磁法一维正演与探测深度[J].地球物理学进展，2008，6(23)：1943-1946. ZHU K G，LIN J，LIU C S，et al. One-dimensional forwand and prospecting depth for airborne frequency domain electromagnetic method [J]. Progress in Geophysics，2008，6(23)：1943-1946. (In Chinese)

[5] 米萨克·N·纳比吉安.勘察地球物理电磁法第一卷理论[M].赵经祥，译.北京：地质出版社，199. Misac N Nabighian.Explore geophysics-electromagnetic method[M].Translated by ZHAO jingxiang.Beijing:Geological Publishing House,1992.(In Chinese)

[6] 刘继东，林长佑，杨长福.瞬变电磁资料反演解释方法[M].兰州：兰州大学出版社，2001 LIU J D，LIN Y, YANG C F. Transient electromagnetic information inversion method[M]. Lanzhou: Lanzhou University Press, 2001. (In Chinese)

[7] 牛之琏.时间域电磁法[M].长沙：中南大学出版社，2007. NIU Z L. Time domain electromagnetic method[M]. Changsha: Central South University Press, 2007. (In Chinese)

[8] 朴化荣.电磁测深法原理[M].北京：地质出版社，1990. PIAO H R. Principle of electromagnetic soundingmethod [M].Beijing:Geological publishing, 1990. (In Chinese)

第39卷 第2期2017年3月物探化探计算技术COMPUTING TECHNIQUES FOR GEOPHYSICAL AND GEOCHEMICAL EXPLORATIONVol.39 No.2Mar. 2017

Airborne TEM response simulation based on GPU

SHEN Lei1，DONG Chunhui1,2

(1. School of highway and architecture, Shandong Transport Vocational College, Weifang 261206,China； 2. Institute of Highway and Bridge Engineering,Dalian Maritime University，Dalian 116026,China)

An efficient algorithm for time-domain airborne EM forward simulation is presented. It is based on CUDA (compute unified device architecture) that was released by NVIDIA in 2007. We derivate the time domain airborne EM formula and then write the code with GPU language accordingly. The modeled results demonstrate that the parallel processing capabilities of graphics processors (GPUs). This significant improvement in performance is gained compared to the corresponding CPU-based solver, while maintaining the numerical accuracy. To the best of the authors' knowledge, this is the fastest time-domain solver for modeling the AEM simulation. It would have potential for complex geophysics model calculation and inversion.

GPU; CUDA; airbome transient electromagnetic method; one-dimensional forward

2016-03-18 改回日期：2016-03-30

沈磊(1978-)，男，讲师，主要从事计算机软件应用、工程岩土方面的教学与研究工作，E-mail：shenlei116@163.com。

1001-1749(2017)02-0170-06

P 631.2

A

10.3969/j.issn.1001-1749.2017.02.03