0 引言

大地电磁测深法是天然电磁场探测方法，由于它具有探测深度大、不受高阻屏蔽的影响、对低阻层反应灵敏等优点，该方法始终为人们所关注和研究(刘国栋和陈乐寿，1984；石应骏等，1985；万汉平，2010).

近年来，随着计算机技术的迅速发展，大地电磁测深二维正、反演算法的研究已经趋于成熟；二维反演方法已成为大地电磁测深生产与科学研究中主要的反演计算方法.但是目前大地电磁二维程序在计算较细网格多频率模型时需要计算的时间较长.我们迫切需要发展一种更高效的算法来提高计算效率.

在本项研究中，我们实现了MT二维正演和反演的MPI并行算法，本文将具体阐述大地电磁二维正反演MPI并行算法的实现流程，并结合算例进行讨论，给出程序的正确性验证和效率分析.

1 大地电磁二维正演和反演 1.1 二维正演

按照二维大地电磁问题的惯例，把整个地球模拟为半空间导体(z≥0) 作为二维地质模型, 其上方为绝缘空气层.电磁场源为在地表或地表上方z=－h处的平面电流层.

(1) 问题的数学物理模型

描述电磁场分布的麦克斯韦方程组的微分形式为：

(1)

(2)

选取笛卡儿直角坐标系：坐标原点在地面，z轴垂直向下，x轴平行于模型走向方向.对于给定的地球模型物理参数，麦克斯韦方程组可以分解成为横电模式(TE)和横磁模式(TM)两种极化方式(Madden，1972；考夫曼和凯勒，1987；Wannamaker et al., 1987；陈乐寿等，1989).

(2) 有限差分法做二维大地电磁数值模拟

利用“传输面相似对比”原理，从而两种极化方式下的麦克斯韦方程组可以写成如下的统一形式，这里V和I分别是电压和电流；Z和Y分别是单位长度上的分布阻抗和分布导纳参数.公式为：

(3)

(4)

(5)

大地电磁二维正演问题实质上是解偏微分方程，要得到偏微分方程的解，结合涉及问题的边界条件(周熙襄和钟本善，1986；纳比吉安，1992；徐世浙，1994；谢飞，2006；曾国，2008).当网格参数和边界条件给定后，根据麦克斯韦方程组(3~5) 就可以对每一个节点写出电流连续性方程，公式为

(6)

其中，i=1, …, M为水平网格数，j=1, …, N为竖直网格数，S_ij为外加线电流源，仅在TE极化方式情况下在网格的顶部z=－h处存在，Y_ij为所论节点对地的集总导纳，Z_connecting为相邻节点间的集总阻抗.

对每个节点上的电磁场分量按顺序排列，可将方程(6) 写成标准矩阵形式，公式为

(7)

其中K为系数矩阵，v在不同的极化方式下分别对应为电场E_x或磁场H_x，s为方程右端向量.

1.2 非线性共轭梯度二维反演

反演问题可以写为：

(8)

其中，d为数据向量，m为模型向量，e为误差向量，F为正演模拟函数.

数据向量d=[d¹, d², …, d^N]^T，dⁱ是某一特定极化方式(TM或TE)，某一观测点上，某一频率(ω)视电阻率ρ_app的振幅对数值或相位值.模型向量m=[m¹, m², …, m^M]^T是模型参数向量.为了和正演数值模拟方案一致，取M为模型块的个数，取m^j为其中某一块模型的电阻率对数值logρ.

目标函数定义为：

(9)

其中，正则化参数λ为正数；正定矩阵V为误差向量e的方差；取矩阵L为简单的二次差分拉普拉斯算子，使得当模型网格是均匀时，Lm可以近似为logρ的拉普拉斯算子.这目标函数具有数据拟合差最小和模型最光滑的双重约束.目标函数的第一项为数据的拟合差，目标函数的第二项为模型的光滑度.

J为目标函数F的雅可比矩阵，非线性共轭梯度反演算法不需要实际计算雅可比矩阵J(Smith and Booker, 1991；Mackie and Madden, 1993；Rodi and Mackie, 2001；胡祖志等，2005)，而只需要计算雅可比矩阵J或它的转置J^T和一任意向量x相作用的结果(吴小平和徐果明，2000；王家映，2002；刘小军，2007).下面，先给出非线性共轭梯度迭代反演的计算流程.

从流程图中可以看出，整个反演过程的主要计算量(除了常规的第一次正演以外)在于：计算目标函数的梯度g=－2J^TV^－1e+2λL^TLm和计算f=Jp.只要计算出了J^TV^－1e(程序中用Jty表示)和Jp的值，α=－p^Tg/2(f^TV^－1f+λp^TL^TLp); 就可以得到步长α，更新模型m=m+αp(谭捍东，2000；林昌洪，2004)，继续重复上述计算直至达到结束条件.

2 大地电磁二维正反演MPI并行算法 2.1 MPI并行计算简介

通过MPI实现了数据的广播，发送，接收和数据的同步. MPI也支持多种数据类型，包括复数.虽然看起来MPI的程序是一套程序，但是可以通过进程ID号进行区分，各进程可分配不同的计算任务.

MPI程序中，所有的变量(包括全局变量和局部变量)、函数，只要没有明确的区分，那么每个进程都享有这个变量和函数，这些变量对各进程来说名字是相同的，但是装载的信息很可能不相同，各进程若想了解其他进程的数据须进行通信.同时，对于alloctable的数据类型来说，每个进程可以为指针分配大小不一的内存空间.特别是MPI应用在频率域的时候，每个进程分配的任务很可能不完全相同，那么就需要灵活的分配空间.

2.2 MT的MPI并行反演算法

在MT二维正、反演计算中，频率循环部分求解电磁场值消耗了主要计算时间，占据了全部求解时间的90%以上，是进行并行化的重点环节.二维反演根据不同的频率求解Jty和Jp.对于给定的不同频率值，方程组求出的磁场值是相互独立的，根据这一特点我们可以将串行算法按频率进行划分，将每个频率对应的计算任务分配到不同进程上同时进行计算(都志辉，2001；张武生等，2009).

程序使用的9个频率分别是freq(1~9) 100、31.6、10、3.16、1、0.316、0.1、0.0316、0.01 Hz.每个频率的计算所需的时间随着频率的变小递增.该算法采用主从并行模式，分为主进程和子进程，0进程为主进程，负责维护全局的数据和任务的分配、参数的广播、计算结果的回收、整合、广播以及计算结果文件的输出，子进程从主进程处接受参数，执行分配任务的计算，并将结果发回主进程.由于主进程用于进程间交互的计算量不大，为了充分利用计算资源，主进程同时也参与高频部分的计算.二维正反演并行计算的基本思路是：

(1) 利用函数MPI_INIT()初始化并行环境，主进程读入所有频率值，模型网格剖分和观测数据，采用MPI_Bcast将其广播给所有子进程，其中模型全部使用背景电阻率作为初始模型.

(2) 主从进程按照分配的频率(表 1)，同时各自独立的计算，通过第一次正演求解得到各节点的磁场值，观测点的视电阻率ρa(程序中为modrestm_mpi)，通过第一次“拟正演”求解得到Jty(程序中为Jty_mpi)，所有进程按频率计算完后，主进程使用MPI_Gatherv收集所有进程的视电阻率，其中Jty_all记录了9个频率所产生的Jty值，还要对其进行求和；收集到的视电阻率值顺序为1、7、2、3、6、4、8、5和9，还需对视电阻率值进行排序.主进程完成上述操作后，MPI_Bcast将视电阻率值，广播给所有子进程.

(3) 计算目标函数Ψ，梯度函数g.

(4) 主进程可得到拟合差rms，若rms < 1.5，循环中止转6，否则程序进行第2次拟正演求解，转5.

(5) 主从进程按照分配的频率同时各自独立的计算Jp(程序中为Jp_all)，待所有进程完成求解Jp后，MPI_gatherv()，主进程收集所有进程的Jp，Jp记录了根据9个频率计算出的Jp值，顺序为1, 7, 2, 3, 6, 4, 8, 5, 9.还需对Jp进行排序.MPI_Bcast()将Jp广播给所有子进程.得到Jp后我们可以计算模型更新步长α，对模型更新m=m+αp，程序转2进行正演求解.

(6) 达到循环中止条件后，主进程将视电阻率和相位输出，MPI_FINALIZE()结束并行环境.图 2为大地电磁二维反演并行计算的流程图.

2.3 并行算法的开发和运行环境

并行程序的开发环境见表 2，编译和运行指令:Mpirun-np N./nlcg，其中N为进程数.

3 算例与结果讨论 3.1 低阻模型

假设一个低阻异常棱柱体，其顶面埋深1500 m，大小为2500 m×5000 m.背景电阻率为100 Ω·m，异常体电阻率为10 Ω·m，如图 5~9所示.模型网格剖分为400×400，数据采集点为地表150、170、180、190、200、201、211、221、231和251处.使用从低到高为0.01 Hz到100 Hz的9个频率(100、31.6、10、3.16、1、0.316、0.1、0.0316、0.01 Hz).各网格单元沿y方向的剖分间隔(单位：m)：

100000 50000 20000 10000 5000 2000 1000 500 200 100 100 100 100 100 100 100 100 100 100 100

…

100 100 100 100 100 100 100 100 100 100 100 200 500 1000 2000 5000 10000 20000 50000

100000

各网格单元沿z方向的剖分间隔(单位：m)：

100 100 100 100 100 100 100 100 100 100

100 100 100 100 100 100 100 100 100 100

…

100 200 500 1000 2000 5000 10000 20000 50000

100000

通过并行计算得到的TM模式200号点正演演效果对比图(视电阻率和相位)如下：图 4中左侧的图为视电阻率随频率变化图；右侧为相位随频率变化图.横轴都为频率的导数取对数，mpi计算的结果用o表示，串行计算的结果用+表示，计算结果完全重合.

通过并行计算得到的TM模式反演效果图 5.色棒的取值为视电阻率取对数，背景视电阻率为100 Ω·m，取对数为2，异常体视电阻率为10 Ω·m，取对数为1，通过图 5可看出程序较好的恢复出了异常体的位置.

3.2 程序正确性的验证和效率分析

由于本文的并行计算完全是在原有串行算法上实现的，只是根据各频率的计算是独立的特点，将不同频率分配给各进程并行处理，对算法本身未做大的修改，因此，并行计算结果应该和串行计算结果一致.本研究对两种计算方式得到的正演结果进行了比对，图 4比对结果显示两种计算方式正演结果一致.图 5是TM模式反演出来的电阻率断面图，从图中可以看出比较好的恢复出了低阻体的空间位置和异常值的大小，从而验证了并行程序是正确的.

由于各频率的反演计算相互独立，因此将多个频率的计算任务分配给多个进程并行处理，从而提高了反演的计算效率.这里要讨论大地电磁非线性共轭梯度二维反演MPI并行程序的计算效率，由于采取不同的进程数来分配计算任务，这里用并行加速比和并行效率来量化并行算法的效率.并行加速比＝单进程串行算法的运算时间/多个进程并行算法的运算时间.加速比越高，说明并行算法的加速效果好，加速比低，说明并行算法的加速效果差.并行效率＝并行加速比/N，N为并行计算的进程个数.并行效率高，说明并行算法中每个进程的利用率较高，反之，说明算法中每个进程的利用率较低.通过统计，我们得到表 3，表中数据显示了400×400网格计算时，每个正演求解所用的时间.通常情况下高频的计算时间短，低频的计算时间较长，时间差距能达到15倍.因此在进行任务分配时，可以通过频率的高低频搭配使得各进程需要执行的任务更加均衡，从而提高计算效率.

下面通过三个数组来控制每个进程需要计算的频率数据.下面分别介绍了开辟不同进程个数所获得的三个数组.

parr(1:4)=(/2, 3, 2, 2/) //程序开4个进程，进程的计算频率个数分别为2, 3, 2, 2;

pstart(1:4)=(/1, 3, 6, 8/) //每个进程从数组p的1, 3, 6, 8开始取频率；

p(1:9)=(/1, 7, 2, 3, 6, 4, 8, 5, 9/) //第1个进程要计算的频率为1, 7，第2个进程要计算的频率为2, 3, 6；第3个进程要计算的频率为4和8，第4个进程要计算的频率为5, 9；

parr(1:2)=(/4, 5/) //程序开2个进程；

pstart(1:2)=(/1, 5/)；

p(1:9)=(/2, 4, 6, 8, 1, 3, 5, 7, 9/)；

parr(1:6)=(/1, 2, 2, 2, 1, 1/) //程序开6个进程；

pstart(1:6)=(/1, 2, 4, 6, 8, 9/)；

p(1:9)=(/4, 1, 5, 2, 6, 3, 7, 8, 9/)；

parr(1:8)=(/1, 2, 1, 1, 1, 1, 1, 1/) //程序开8个进程；

pstart(1:8)=(/1, 2, 4, 5, 6, 7, 8, 9/)；

p(1:9)=(/1, 2, 3, 4, 5, 6, 7, 8, 9/).

对表 4进行分析，当使用两个进程计算时效率并行效率最高，当进程数为6和8时，加速比变化较小，但是并行效率反而下降，这是因为当节点数增加时，进程通信会占用一定的时间，同时在按频率进行运算时，高频计算时，求解方程需要的时间短，计算低频时，求解方程需要的时间长，详细情况见表 3.所以计算高频的进程需要等待计算低频的进程，例如当进程数为8时，除了2号进程外，其他进程只需执行1个频率的计算任务，而第8个进程需要计算最低频部分，求解方程需要较长的时间，因此其他7个进程将运行结果发送给1进程后，只能等待第8进程，等待的时间从1 s到20 s不等.导致其他进程处于等待状态，所以并行效率降低，加速比没有达到预期.当进程数为2时，由于各进程分配的频率比较均衡，需要运算的时间相当，CPU的利用率较高，并行效率也高，达到81.5%.从表 4中可以看出并行计算的效果还是比较明显的，当采用8进程进行计算时，计算时间由3376 s缩短为1303 s.

4 结论

当大地电磁网格剖分较细时，多频率的反演计算时比较耗时，采用MPI并行计算应用在正反演求解中是提高计算速度的关键，并行计算有效解决二维大地电磁较细网格计算耗时的问题，本文在工作站计算平台上，以MPI为并行编程环境，实现了大地电磁二维反演并行算法.通过对不同频率数据的计算和统计，发现高频数据计算花费的时候少，低频数据计算花的时间多，在对任务进行划分的时候可以通过高低频搭配来提高进程的计算效率.通过模型试算和分析，结果证明了该算法的正确性和高效性，很好的解决了大网格情况下大地电磁非线性共轭梯度二维正、反演计算速度问题.

致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！