复电阻率法二维正演并行计算研究

0 引言

复电阻率法(CR)又称为谱激发极化法(SIP)是以岩矿石电阻率的频谱特性差异为基础，通过观测大地的视复电阻率频谱或激发极化场的衰减曲线，来寻找极化目标体的一种电法勘探方法[1]。相对其他电法勘探分支方法，复电阻率法有以下优点：①该方法能得到多个电性参数的信息，结合这些参数的对比解释可提供更丰富的地电信息；②采用轻便的采集装备，使其在复杂地形地区开展勘探工作更有利。目前，复电阻率法已成功应用在固体矿产勘探[2-4]、水文地质[5]、油气资源[6] 、监测环境[7]、工程地质[8]等方面。

目前计算机技术发展迅速，将MPI等并行计算技术引入到地球物理数值模拟中是当前研究热点。当网格剖分单元较多，多发射源、多频率的情况下，完成复电阻率法二维有限元数值模拟需要的计算时间较长，因此有必要开发复电阻率法二维正演并行算法来提高其计算效率。

1 MPI简介

MPI(Message Passing Interface)是一个“库”而非一门语言，其具有C++和FORTRAN90的函数调用接口[9]。笔者采用Fortran和MPICH2相结合的开发工具编写并行程序。

MPI程序设计包含两种并行模式：①主从模式；②对等模式[10](图1)。

1)主从模式，也就是在所有进程中定义其中的一个作为主进程，除主进程外的所有进程都当作子进程。通过主进程可以进行广播数据或发送数据等一系列操作，同时也可以将分配好的任务分发给子进程。等到所有的子进程都完成计算以后，主进程通过收集操作从每个子进程得到任务结果，然后再看是否需要下一步的操作。

2)对等模式，就是所实现的MPI程序中各个进程间的计算量和功能大体一致。

图1 并行模式分类
Fig.1 Parallel mode classification
(a)对等模式；(b)主从模式

2 复电阻率法二维正演

2.1 二维数值模拟

二维地电模型如图2所示，其中构造走向方向为y轴，假设其电性参数仅仅在x-z平面是变化的而沿y轴方向不变。

图2 二维地电模型示意图
Fig.2 2D geoelectric model

假设时间因子是eiωt，忽略位移电流的影响，二次场满足的麦克斯韦方程组为[11]

(1)

式中：Es为电场二次场；Hs为磁场二次场；Ep为电场一次场；为角频率；μ0为真空中磁导率；σ0、σ分别为模型和背景电导率，σa=σ-σ0是异常电导率。

傅里叶变换

(2)

沿走向方向对式(1)作傅里叶变换，整理后得到波数域中满足的麦克斯韦方程组

(3)

(4)

其中:“^”为波数域的值；ky为波数值；在忽略位移电流的条件下，

利用四节点的等参单元法进行数值模拟，将式(3)和式(4)分别运用伽里金法、格林公式，且正演使用第一类边界条件(边界网格上的电磁场值为零)，可以得到本文二维数值模拟满足的有限单元法计算公式[12-16]：

(6)

其中:Ne为剖分单元的总数量；表示在第e个剖分单元中第i个节点的插值形函数。

2.2 Cole-Cole模型

为同时考虑激电效应和电磁效应，将大地的实际电阻率由Cole-Cole模型定义得到的复电阻率替换。Pelton等[17]通过对大量岩石、矿石标本的一些研究，总结出就均匀岩石、矿石而言，它们的复电阻率频谱可以用Cole-Cole模型来表示。

(7)

其中:ρ(iω)是复电阻率；m是极化率；ρ0是零频电阻率；τ是时间常数； c是频率相关系数。

2.3 有限单元计算公式

引入Cole-Cole后，通过单元分析和系数矩阵存储，最终正演计算形成的矩阵方程可写成如下形式

(8)

其中；和分别是网格节点处待求的波数域二次电磁场值；S是与网格节点处波数域一次场值有关的向量。k1、k2、k3、k4、s1和s2可表示如下:

k1(i,j)=

(9)

(10)

(11)

k4(i,j)=

(12)

(13)

(14)

2.4 矩阵的压缩存储

正演矩阵方程(8)中；左端项系数矩阵K是稀疏且对称的，全部存储K矩阵需要的内存较大。由于矩阵K每行的非零元素个数有限，充分利用这个特点，可以有效地节约内存空间。结合有限单元法的特点，单元分析时每个节点仅与其周围的几个节点参与运算。以笔者采用的4节点等参单元为例，最终形成的矩阵中每行最多只有18个非零元素，在网格剖分后，可根据有限元单元分析特点得到非零元素的行列号。采用CSR(按行压缩存储)方法，通过三个一维向量就可完成矩阵的存储，便于使用共轭梯度法来求解方程组。

2.5 求解方程组

正演形成的方程组采用林绍忠[18]介绍的SSOR—PCG(对称逐步超松弛预处理共轭梯度法)方法解方程。

解形如Ax=b的方程时，SSOR法的预处理矩阵M为：

M=(2-ω)-1(D/ω+L)(D/ω)-1(D/ω+L)T

(15)

其中:D为A的对角阵;L为A的严格下三角矩阵;ω为松弛因子，本文中将松弛因子取为“1”，初始x为零向量。那么迭代格式可以简化为

(16)

2.6 正演正确性验证

如图3所示，在背景电阻率为100 Ω·m的地下介质中存在电阻率为10 Ω·m的二维棱柱体，其顶界面埋深为120 m，长为200 m，厚为100 m。发射源沿x轴水平方向放置，且位于x = y = z = 0处，从x = 50 m 到 x = 1 000 m放置25个接收点，频率f = 8 Hz、16 Hz。将计算结果与Kerry Key的2DCSEM程序[19]计算结果进行对比，如图4所示。2DCSEM计算结果用红色圆圈表示，本文程序计算结果用蓝色点表示，电场分量实部曲线及Ex虚部曲线基本吻合，验证了该程序的正确性。

图3 二维棱柱体模型
Fig.3 Two-dimensional prism model

3 二维正演并行算法总体设计

在对开发的复电阻率法二维正演串行程序深入分析后，可以发现正演中，其大部分计算时间都花费在频率循环部分求解电磁场值上面，每个发射源对应频率的方程组的求解等都是相互独立的且互不影响的，并行性较好，因此我们采用主从模式，基于多发射源和多频率来实现正演并行计算，并对其并行效率进行分析。并行的基本思路是：①将MPI环境初始化MPI_INIT( )；②主进程读入所需的参数文件；③主进程广播数据到各子进程MPI_Bcast( )；④各进程分别计算自己的任务；⑤主进程收集和整理所有进程的结果MPI_GATHERV()；⑥结束MPI并行环境MPI_FINALIZE( ) [19-21]。

图4 二维棱柱体模型计算得到Ex实虚部对比曲线
Fig.4 The Ex response comparison between the 2D finite element results and the results from 1DCSEM
(a)f=8.0 Hz; (b)f=8.0 Hz; (c)f=16.0 Hz; (d)f=16.0 Hz

3.1 模型

如图5所示，在背景电阻率为100 Ω·m的地下介质中存在一个长200 m厚100 m的异常体，其顶部埋深为160 m，其中ρ0=10 Ω·m，m=0.5，c=0.5，τ=30 s。采用的装置为偶极-偶极，沿x方向水平放置电偶极子，xz方向网格剖分为121×64。以下正演并行程序均在18个节点的集群上运算，其计算环境均为：Linux操作系统，CPU-Intel(R) Xeon(R) CPU E5-2620 v2 @ 2.1GHz，内存128G，开发语言为Fortran，编译器为ifort，并行环境为MPICH2。

图5 低阻高级化模型
Fig.5 Low resistivity and high polarization model

各个发射源对应频率的计算是独立的，基于发射源对应的频率来进行并行计算的思路，其流程如图6所示。

3.1.1 效率分析

采用图5所示的模型，21个发射源(其位置从X = -1 000 m到X = 1 000 m之间每隔100 m一个，共21个源)，从X= -1 000 m到X= 1 000 m之间每隔50 m布置一个测点，每个排列除发射源点外共40个接收点。5个频率(f= 0.1 Hz, 1.0 Hz, 8.0 Hz, 32.0 Hz, 128.0 Hz)，用不同节点数来进行并行计算，将其计算时间和串行程序计算时间统计如表1所示。

图6 正演并行流程图
Fig.6 Forward parallel flow chart

表1 并行计算时间统计表

Tab.1 Parallel running time

从表1可以看出，在每个节点只有一个进程参与计算的前提下，当参与并行计算的节点数增加时，加速比增大。当有7个节点参与并行计算时，其计算速度约为串行程序计算速度的5.847倍，很大程度地降低了复电祖率法二维正演所需时间。但是，从表2中并行效率这一列看出，并不是参与并行的节点数越多，并行效率越高。当参与并行的节点数从2增加到3时，并行效率增加，而当参与并行的节点数从3增加到4(或节点数从5增加到6)时，并行效率反而降低，这是由于节点等待或节点间的数据交换等通信耗费一些时间。

3.1.2 正确性验证

采用多发射源和多频率进行正演计算时，笔者仅展示其中一部分数据的对比结果。如图7所示，图7中Tx、Rx分别代表发射源和接收点的位置，串行计算结果用红色圆圈表示，并行计算结果用蓝色点表示，计算结果完全吻合。

4 结论

当发射源、频率较多时，完成复电阻率法二维正演所需时间较长，我们通过引入MPI并行计算，实现了复电阻率法二维正演并行算法，极大地提高了正演的计算效率。通过理论模型的试算结果表明，该并行算法是可靠的，稳定的，较好地解决了复电阻率法二维数值模拟的计算速度问题，为其反演并行计算打下良好的基础。对其并行计算时间分析可知，当每个节点的进程数都为“1”时，随着节点数的增多，其并行加速比逐渐增大，但其并行效率并不与参与计算的节点数成正比。

图7 正演并行结果与串行结果对比图
Fig.7 The comparison between forward modeling parallel and serial results
(a)Tx=-1 000 m Rx=-300 m; (b)Tx=-1 000 m Rx=-300 m;(c)Tx=-1 000 m Rx=0 m; (d)Tx=-1 000 m Rx=0 m;(e)Tx=-1 000 m Rx=300 m; (f)Tx=-1 000 m Rx=300 m;(g)Tx=-500 m Rx=100 m; (h)Tx=-500 m Rx=100 m