基于MPI瞬变电磁测深一维反演并行算法探究

柳建新，刘鹏茂，刘 颖，童孝忠

(中南大学地球科学与信息物理学院，长沙 410083)

基于MPI瞬变电磁测深一维反演并行算法探究

柳建新，刘鹏茂，刘 颖，童孝忠

(中南大学地球科学与信息物理学院，长沙 410083)

在瞬变电磁测深反演中运用并行技术可以减少计算时间，提高反演的运算效率。MPI(MessagePassingInterface)是目前最重要的并行编程工具，它具有移植性好、功能强大、效率高等多种优点。这里基于在Windows系统下使用FORTRAN和MPICH2相结合的开发工具，编写瞬变电磁并行算法程序，对瞬变电磁一维采用直接反演法，通过理论模型对该算法进行试算，计算结果证明了该算法的正确性、高效性和稳定性。

瞬变电磁;直接反演;并行计算;MPI

0 前言

二十世纪以来，计算机技术的广泛应用，极大地促进了科学的发展、社会的进步，推动了计算机的不断更新和完善。随着计算机在硬件方面的快速发展，已经大大提高了计算机微处理器的处理速度。在矿产地球物理勘探中，信息技术是勘探的支撑，而对高性能计算的需求和依赖，是勘探信息技术的前进方向。因此，促使了地球物理领域并行计算技术的飞速发展。

并行计算可以加快速度，在更短的时间内解决相同的问题，或者在相同的时间内解决更多更复杂的问题;也可以节省投入，以较低的成本完成同量的任务。因此，并行计算技术将会成为地球物理数据处理的主流，在国内，并行处理各种算法研究取得了一定的成果，陈金窗、戴光明等［17］在微机网络上实现了2.5维CSAMT正演的并行计算;谭捍东等［18］实现了三维模型的大地电磁正演并行计算，随后又实现了三维快速松弛反演的并行计算，证明了此方法的有效性和稳定性;刘羽［12］实现了二维大地电磁Occam反演的并行计算;陈露军、王绪本等［19］以电磁场理论应用为背景，把并行思想和电磁勘探中的电磁计算问题结合起来，研究了在曙光TC－4000微机集群上，基于的分布式并行编程环境下的求解三维电磁场正演模型的问题;王月英［19］实现了三维波动方程的有限元并行正演模拟，采用有限元并行算法，克服了单机对数据容量和计算速度的限制，可以较串行算法模拟更广泛区域的三维空间;这些都大大地推进了国内并行计算在电磁领域的发展。在国外，Newman和 Alumbaugh［20］在1995年运用并行算法来进行电磁成像;1997年美国洛伦茨贝克利国家实验室的谢干权等［18］人把并行计算技术运用到全局积分局部微分并行域分解的新算法中进行电磁场的反演，并取得了一定的效果;2000年，Zyserman和 Santos［19］实现了大地电磁三维有限元的并行计算。

A.G.Nekut［7］在 1987 年发表了一篇摘要，主要介绍了瞬变电磁一种简单的一维直接反演方法。这种方法是把瞬变电磁法测深数据转换为电导率，测深剖面图的方法是基于对衰变斡旋电流的近似图像表示方法，它可以为任何时域电磁测深数据提供一维反演。作者是在Windows系统上使用FORTRAN和MPICH2相结合的开发工具，编写了瞬变电磁测深的一维直接反演并行程序，通过理论模型对该算法进行试算，计算结果证明了该算法的正确性、可行性和稳定性。可以预知，当计算量比较大时，特别是对于瞬变电磁二维、三维的计算时，这种并行技术优越性可以得到更好的体现。

1 MPI简介

MPI(MessagePassingInterface)是为了统一不同厂商的消息传递API(ApplicationInterface)，有来自高性能计算领域的专家和MPP厂商的代表组成的委员会制定的工业标准，是消息传递并行程序设计的标准之一。它能为高性能并行计算提供一个方便的环境，具有移植性好、功能强大、效率高等多种优点。MPI是一个库，而不是一门语言。按照并行语言的分类，可以把FORTRAN+MPI或C+MPI，看作是一种在原来串行语言基础之上扩展后得到的并行语言，MPI库可以被FORTRAN77/C/FORTRAN90、C++调用。

MPI并行程序设计平台由标准消息传递函数及相关辅助函数构成，多个进程通过调用这些函数，进行通信。MPI程序采用的是SPMD执行模式，即一个程序同时启动多份，形成多个独立进程，在不同的处理机上运行。每个进程开始执行时。将获得一个唯一的序号。

MPI的两种最基本的并行程序设计模式，即对等模式和主从模式。

(1)对等模式即参与运算的各进程地位相同，计算程序一致，只是处理的数据不同。

(2)主从模式则是在程序中设置主进程和从进程，主进程向备从进程发送数据，从进程进行计算。然后将计算结果返回主进程，主进程收到结果后发送结束标志，使从进程退出执行。

2 直接反演方法

2.1 方法原理

应用中心回线测深装置，在地表上放置半径为R的一个环形线圈中通有的电流为I，在t=0时刻突然断电，一个磁力计在环形线圈中心测量中心点垂直方向上的磁场H(t)，它的像源深度D(t)可以用H(t)的变化趋势来反映［17］，其关系如式(1)。

D=R{［H(t)/H0］－2/3－ 1} (1)式中 H0=I/2R。

依照此方法，通过与放置在距水平良导面正上方δ的电磁源进行类比，可得象源的深度D与探测深度(穿透深度)δ的关系式(2)。

对于电导率为δ的均匀大地介质，瞬变场的有效勘探深度可写为式(3)。

式(3)等价于频率域中趋肤深度的公式。

因此，只要给出了地下均匀介质的电导率σ，通过公式(1)、公式(2)、公式(3)，就可以得到H(t)的近似值。

为了校正像源场合瞬变场在均匀半空间条件下的误差，可以在穿透深度上加一个校正系数F。对于均匀半空间模型，校正系数可以用如下函数来表示:

或

F¯(δ/R)与δ/R的关系，我们可以通过均匀半空间中计算 δ/δ得到一组F¯(δ/R)和 δ/R，然后通过三次样条插值法，得到层状介质情况下的F¯。

因此，给定任何一个H(t)，通过公式(1)、公式(2)、公式(3)联立算出场的有效勘探深度，即

电导率的阶数n对应数据的n阶连续微分。

一阶视电阻率为:

二阶视电阻率为:

也可以间接表示为:

图1 (δ/R)与δ/R的关系图Fig.1 Thediagramof(δ/R)and δ/R

3 并行算法实现

该瞬变电磁一维并行算法程序中，正演方法采用数值滤波方法，系数为801。反演方法采用的是直接反演法，并把反演结果存储在 S0、S1、S2.dat文件中。通过对反演解释中微分算字符的处理，把线性坐标下的微分，转化为双曲坐标下的微分，取得的数值更加精确，反演精度S0－S2精度逐渐升高。实现的并行算法程序是在双核单机下模拟多进程进行，操作系统使用的是Windows系统，PC机配置是 GenuineIntel(R)CPUT1600，主频1.66GHz，1G 内存。

图2 (δ/R)与δ/R的关系图(σ=0.01，10－8＜t＜1s)Fig.2 Thediagramof)F(δ/R)and δ/R(σ =0.01，10－8＜t＜1s)

算例一:对电阻率为三层的理论模型进行并行计算，层数(NC)=3;每层的电导率:σ1=0.01、σ2=0.001、σ3=0.005;每层的深度:H(1)=100、H(2)=100;电流(AMP)=25;天线长度(RAD)=100。并行计算结束需要的时间如表1所示。

表1 理论模型1并行计算不同进程运行时间Tab． 1 Different processes running time of theoreticalmodels 1 parallel computing

该并行计算程序对直接反演方法没有做任何修改，反演结果如图3所示。

算例二:层数(NC)=3;每层的电导率:σ1=0.001、σ2=0.01、σ3=0.001;每层的深度:H(1)=100、H(2)=100;电流(AMP)=1;天线长度(RAD)=200。并行计算结束需要的时间如表2所示。

表2 理论模型2并行计算不同进程运行时间Tab． 2 Different processes running time of theoreticalmodels 2 parallel computing

从计算时间结果可以看出，应用并行计算可以加快速度，在更短的时间内解决相同的问题。反演结果如图4所示。

图3 模型1电导率—深度反演结果图Fig． 3 The inversion results map of electrical conductivity-depth of model1

图4 模型2电导率—深度反演结果图Fig． 4 The inversion results map of electrical conductivity-depth of model2

4 结论

(1)瞬变电磁一维正演计算，常采用数值滤波方法，即运用汉克尔变换计算频率域一维响应，然后通过余弦变换法从频率域转换到时间域中。数值滤波法具有精度高，时间带宽长，不同滤波系数之间无相关，因此特别适合并行计算。

(2)通过对均匀半空间模型的数值解与解析解对比，计算其相对误差，验证正演算法的正确性。

(3)对瞬变电磁测深一维采用直接反演法，进行并行计算，通过理论模型对实现的并行程序进行试算，计算结果证明了该算法的正确性、高效性和稳定性。

(4)由于反演算法中)F(δ/R)与δ/R的关系是确定的，为了确定)F(δ/R)与δ/R的关系，我们通过均匀半空间中计算)δ/δ得到一组)F(δ/R)与δ/R，采用三次样条插值法，得到层状介质情况下任意δ/R值对应的)F(δ/R)。通过对多个均匀模型的计算，我们发现，σ=0.01的均匀半空间模型计算的)F(δ/R)与δ/R关系，可以用来作为其它任何模型的插值基数。

随着地球物理勘探向深度和广度方向发展，对勘探精度的要求也越来越高，这就会使需要处理的数据量越来越大，大量的数据量使得并行技术将会在地球物理勘探领域中起到越来越重要的地位和作用。并行计算技术可以加快速度，节省投入，在瞬变电磁反演资料处理技术中，对于其二维、三维的计算，并行技术将会得到更好的应用和发展。

［1］ 罗省贤，何大可.基于MPI的网络并行计算环境及应用［M］.成都:西南交通大学出版社，2001.

［2］ 李炎，胡祥云.基于MPI的一维大地电磁并行计算研究［J］.地球物理学进展，2010，25(5):1612.

［3］ 谢丽，胡文宝，陆辉.瞬变电磁响应计算的并行算法研究［J］.长江大学学报，2009(6):4.

［4］ 黄易，师学明.并行计算技术及其在勘探地球物理学中的现状和展望［J］.地球物理学进展，2010，25(2):642.

［5］ 荣莹，曹俊兴.基于MPI的机群并行计算系统平台构建［J］.物探化探计算技术，2005，27(4):89.

［6］ 牛之琏.时间域电磁法原理［M］.长沙:中南大学出版社，2007.

［7］ NUKET A G． Direct inversion of time． Domain in electromagneticdate． Geophysics， 1987，52( 10) : 1431．

［8］ 瞬变电磁测深法原理［M］.西安:西北工业大学出版社，1993.

［9］ JONES A G． On the equivalence of the“Bostick”and“Niblett”transformations in the magnetotelluric method［J］． Geophys， 1983， 53: 72．

［10］严良俊，胡文宝.中心回线瞬变电磁测深法快速电阻率成像方法及应用［J］.煤田地质与勘探，2002，30(6):58.

［11］郭文波，宋建平，韩俊明，等.中心回线瞬变电磁探测的一种快速解释方法［J］.物探与化探，2006，30(2):616.

［12］ MACNAE J C，LAMONTAGNE Y． Data processingtechniques for EM sounding and prospecting［J］． Geophys，1984， 29( 1) : 464．

［13］都志辉.高性能计算并行编程计算－MPI并行程序设计［M］.北京:清华大学出版社，2001.

［14］孙家旭，张林波，迟学斌，等.网络并行计算与分布式编程环境［M］.北京:科学出版社，1996.

［15］ MACNAE J C，LAMONTAGNE Y． Imaging quasi． layeredconductive structures by simple processing of transientelectromagnetic data． Geophysics， 1987，52: 45．

［16］石明娟，罗延钟.西方地面瞬变电磁法理论的发展现状［J］.国外地质勘探技术，1989(1):255.

［17］陈金窗，戴光明.微机网络并行计算及2.5维CSAMT正演的并行实现［J］.物探化探计算技术，1997，19(2):103.

［18］ TAN HANDONG，TONG YOU，LIN CHANGHONG． Heparallel 3D magnetotelluric forward modeling algotithm［J］． Applied Geophysics， 2006，3( 4) : 197．

［19］FABIO I，ZYSERMAN JUAN E，SANTOS． Parallel finiteelement algorithm with domain decomposition forthree － dimensional magnetotelluric modeling［J］． J． Appl．Geophys， 2000， 44: 337．

［20］ NEW MAN G A，ALUMBAUGH D L． Three-dimensionalmassiverly parallel electromagnetic inversion［J］．Theory，Geophys，J ． Int， 1997， 128: 345．

P631.3+25

A

1001—1749(2011)05—0491—05

湖南省高校科技创新团队(2008－244);湖南省重点实验室资助项目(2010TC2007);湖南省重大专项(2008FJ1006)

2011－01－16 改回日期:2011－06－17

柳建新(1962－)，男，博士，教授，现主要从事大地电磁理论与研究工作。