随着电磁勘探和计算机技术的发展，电磁法多维正反演已基本成熟，正逐步推向实用化(赵国泽等，2007；Kelbert et al., 2008).然而，勘探精细化和三维探测急需发展，大规模和海量数据的电磁正反演计算成为当前难点和热点问题之一.正演是反演的基础，如何实现大规模网格的快速高精度电磁法正演模拟是问题的关键.电磁法正演通常利用有限差分法(FD)、有限单元(FE)、积分方程(IE)等方法对基本方程进行离散，结合边界条件形成大型稀疏线性方程组，然后利用迭代法或直接法进行求解，从而得到电磁场的解.电磁法正演的实质是大型稀疏线性方程组的存储和求解，随着剖分网格的增加，线性方程组的所需内存和计算时间成指数增加，因而对于大规模剖分问题的求解难于在PC机上实现.为提高多维电磁法正演效率，地球物理学者引入了大型稀疏线性方程组的多种快速求解算法，如预处理共轭梯度法(CG)(Smith，1996)、预处理双共轭梯度法(BICG)(Lin et al., 2011)、多重网格法(MG)(Mulder，2008；邓居智等，2011)、最小残差法(QMR)(Siripunvaraporn et al., 2002)等.这些算法虽然能够大大加快电磁法正演速度，但是所需内存空间并未减少甚至增加，因而仍难于实现大规模多维问题的计算.区域分解算法(DDM)通过将大规模问题分解为若干个小问题进行求解(Toselli and Widlund, 2005)，缩小计算规模，极大节省内存空间，已成为多个科学工程领域解决大规模工程数值模拟的有力手段，如：电磁工程(厉天威等，2007；吕志清和安翔，2009；Shang and He, 2009；盛亦军等，2010；胡俊等，2011)、热传导(顾彩梅，2006；田敏和羊丹平，2008；修航，2014；金鑫禹，2015)和结构力学(刘佩和姚谦峰，2008；瞿叶高等, 2012, 2013；周忠双，2013)等领域.为此，引起地球物理学者极大兴趣，他们在多维正演模拟中引入区域分解算法.对于地球物理电磁正演问题，Zyserman等(1999)及Zyserman和Santos(2000)采用非重叠型区域分解算法分别实现电磁法2D和3D正演；Xiong(1999)采用自适应Schwarz重叠型区域分解算法实现可控源大地电磁法三维正演；Rung-Arunwan和Siripunvaraporn(2010)利用改进的分级区域分解算法实现二维MT的正演.

有鉴于此，本文以TM模式下大地电磁二维正演为例，将Schwarz重叠型区域分解算法与有限差分算法结合，实现基于Schwarz重叠型区域分解的大地电磁二维正演模拟，同时探讨子域数目、子域分解方式及子域重叠区域大小三个因素对计算效率的影响，为大规模问题的求解提供参考和依据.

在二维地电结构情况下，假设x轴为走向方向，y轴为与走向垂直的方向，z轴垂直向下.在准静态条件下，Rao和Babu(2006)给出磁场极化(TM模式)下大地电磁法二维正演模拟的控制方程和边界条件，可表示为

(1)

其中，k²=-iσωμ，ω为角频率，σ为电导率, μ为磁导率，i为虚数单位.

首先，将求解区域采用规则网格剖分成若干个小的矩形网格；然后对于任意一个节点H_x采用5点有限差分格式对等式(1) 进行离散；结合给出的边界条件，能够合成待解的线性方程组，可表示为

(2)

其中，U是大型复稀疏对称矩阵.X表示待解的磁场H_x分量，G为右端项，与边界条件和场源相关.对于线性方程组(2) 通常采用直接分解法或Krylovz子空间迭代法进行求解得到磁场H_x分量的值，然后利用法拉第定律求得电场E_y分量，最后通过卡尼亚电阻率定义公式求得MT响应.

区域分解法分为重叠型(ODDM)和非重叠型(NDDM)两种(Toselli and Widlund, 2005).Schwarz交替法是重叠型区域分解算法(ODDM)中的一种，于1870年由德国学者H.A.Schwarz提出，经过长期推广和应用已成为重叠型区域分解法的经典理论.本文将Schwarz重叠区域分解算法与有限差分算法结合实现大地电磁二维正演模拟.

在Schwarz重叠区域分解算法中，将待求解区域Ω沿Z和Y方向分别分解N和M重叠子域，子域总数为N×M，子域用Ω_{i, j}(i=1, …, N; j=1, …M)表示，子域边界记为 .子域Ω_{i, j}的边界与待求解区域边界重叠部分记为Γ_{i, j}(Γ_{i, j}=Ω_{i, j}∩Γ)，子域Ω_{i, j}其他边界采用虚拟边界，记为Γ′_{i, j}=∩Ω_{r, k}.图 1给出等距剖分情况下分解为3×3个重叠子域示意图，子域与子域重叠两个节点，图中“○”表示子域非重叠节点，“●”表示子域重叠节点，实线框表示重叠子域，阴影表示子域重叠区域.由图可见，对求解区域分别沿y方向(横向)和z方向(纵向)分解成q和p个重叠子域，总个数为p×q个；子域重叠大小通常沿y方向和z方向保持一致，用重叠节点数恒量，分别用m_y和m_z表示.因此，Schwarz交替法求解的收敛速度与重叠域的大小有直接关系，在实际使用中优化子域重叠大小和子域数目.

图 1

Figure 1

图 1 大地电磁二维正演3×3重叠区域分解示意图 (○：子域非重叠节点，●：子域重叠节点，实线框：重叠子域，阴影：子域重叠区域.) Figure 1 Schematic diagram about overlapping domain decomposition method of MT 2-D forward and subdivided into 3×3 (○: Non overlapping node，●: overlapping node，Solid line box: sub-domain，Shadow: Overlapping region of sub-domain.)

在大地电磁二维正演模拟中，根据Schwarz交替法的思想，任意一个重叠子域Ω_{i, j}满足的控制方程和边界条件可表示为

(3)

其中，H′_x为邻域内虚拟边界Γ′ _{i, j}磁场H_x分量.对于控制方程(3) 可以采用与全域方程(1) 一样的离散方式，合成重叠子域待解线性方程为

(4)

其中，U_{i, j}是重叠子域Ω_{i, j}的复稀疏系数矩阵.X_{i, j}表示重叠子域Ω_{i, j}内待解磁场H_x分量，G_{i, j}为右端项，与虚拟边界条件和邻域相关.对比式(4) 和式(2) 可以发现，全域待求解方程组能够转化为若干子域线性方程组进行求解，因而可缩小求解规模.

在大地电磁法二维正演求解过程中，本文利用Schwarz交替法通过子域重叠部分实现参数传递然后通过全域迭代实现求解.首先给定所有节点初值H_x⁰，这样所有子域控制方程(3) 的边界条件就可确定，按顺序独立求解每个重叠子域的线性方程组(4)，把上一个子域的解作为下一个子域的新初值，然后再求解子域的定解问题式(4)，如此反复迭代直到相邻两次解的误差达到某个精度δ，具体流程见图 2.公式为

(5)

图 2

Figure 2

图 2 MT二维正演的重叠型区域分解算法流程图 Figure 2 Flowchart based on overlapping domain

在实际迭代求解中，只要地电模型和重叠区域分解确定后，对于每个子域的稀疏矩阵U_{i, j}保持不变，只需重新计算右端项G_{i, j}.由于重叠子域稀疏矩阵U_{i, j}相对全域系数矩阵U规模缩小很多，在迭代过程只需进行一次分解就能够通过回代实现快速迭代求解.因此本文重叠子域求解算法选用LU直接分解算法.

为了验证和对比本文提出的基于Schwarz重叠区域分解的大地电磁二维正演算法的可行性和有效性，本文算法结果和传统有限差分结果进行对比，同时研究重叠子域大小、子域数的影响规律.本文所有数值模拟均在笔记本电脑(内存4.0 GB, 处理器Intel(R)Core(TM)i5，主频1.70 GHz)上运行，开发环境为Intel Visual Fortran 2013和Visual Studio 2010.

为了验证本文算法的准确性，设计低阻异常模型.低阻异常体大小为26 km×11 km，电阻率为100 Ω·m，顶部埋深为4 km，围岩电阻率为500 Ω·m.整个求解区域大小为60 km×20 km采用等距剖分，剖分间隔500 m，剖分网格数为：120×40，计算频率为0.37 Hz.在重叠区域分解算法中，将求解区域剖分2×1个子域，重叠区域大小m_y=3，m_z=0，子域求解算法采用LU直接求解法，全域迭代精度设为10E-8.图 3为低阻异常体模型大地电磁响应结果对比，(a)为视电阻率，(b)为阻抗相位.由图可知，采用重叠型区域分解算法的计算结果与传统有限差分法的计算结果基本重合，因此验证本文提出算法的准确性和可行性.

图 3

Figure 3

图 3 低阻异常体模型大地电磁响应结果对比和算法验证 (a)视电阻率；(b)相位. Figure 3 Comparison of electromagnetic response of low resistivity model and algorithm verification

为了进一步分析Schwarz重叠区域分解算法的优势和计算效率，本文对子域数目、重叠区域大小和不同子域求解算法进行研究.

a)重叠子域分解方案

为了研究重叠子域数目和分解组合的影响规律，设计不同子域分解方案对低阻异常体模型进行数值模拟，模型参数和算法验证模型一致.模型求解区域剖分网格数为120×40，并分解成p×q个重叠子域.其中p和q分别为Y和Z方向重叠子域数，从1到10和1到5变化，子域重叠节点数均为3，子域求解选用LU直接求解法.对于传统有限差分LU直接求解法，剖分网格大小为M_y×M_z总内存近似计算公式为(M_y－1)(M_z－1)(3M_z－1)×16，16为双精度复数所占bites数(Tawat et al., 2010)，所需内存约为8.4 Mbites，CPU运行时间T为7.094 s.由于本文采用不存储重叠子域分解矩阵，因此对于重叠型区域分解算法的内存计算公式可近似表示为：((m_y×(p－1)+M_y)/p－1)((m_z×(q－1)+M_z)/q－1)(3×(m_z×(q－1)+M_z)/q－1)×16.可以看出，重叠区域分解算法的内存使用量与子域组合、子域数目及子域重叠大小有关.表 1为不同重叠子域分解方案的CPU计算时间T_ODDM和内存使用量统计情况.为了清晰描述重叠区域分解算的优势，本文定义相对内存百分比(U_ODDM-U_FD)/U_FD×100%及相对时间百分比(T_ODDM-T_FD)/T_FD×100%，其中FD表示传统FD算法，ODDM表示重叠区域分解算法.图 4为重叠区域分解算法相对内存百分比和相对时间百分比等值线图，横轴为Y方向子域数目p，纵轴为Z方向子域数目q.由表 1和图 4可知，内存使用量随着子域数目增大而减小；CPU计算时间不随子域数目增大而严格减小，这是由于重叠子域缩小到一定规模后直接分解时间很小而全域迭代次数将增加引起的.在重叠子域分解为6×3时CPU计算时间最小为0.391 s，仅为传统FD算法CPU计算时间的5.5%，同时内存使用量也仅为2.5%.因此，在实际求解过程中应该选择合适的子域数目和剖分组合.

表 1 (Table 1)

表 1 不同子域组合CPU计算时间TODDM及内存使用量 Table 1 CPU time of TODDM and memory usage of different sub-domain combination

表 1 不同子域组合CPU计算时间TODDM及内存使用量 Table 1 CPU time of TODDM and memory usage of different sub-domain combination

图 4

Figure 4

图 4 不同子域数目的大地电磁二维正演模拟效率和内存使用量对比 (a)内存使用量百分比;(b) 计算时间百分比. Figure 4 Contour map of the relative percent about CPU time and memory usage

b)重叠区域大小

重叠型区域分解算法中，子域重叠规模大小是影响算法内存和计算时间的重要参数.为研究重叠区域大小对计算效率的影响，设计不同子域重叠规模大小对低阻模型进行数值模拟.模型参数与算法验证模型一致，求解区域大小和剖分网格数与上个模型一致.采用三种不同子域分解组合，分别为2×1，3×2，2×2，子域重叠规模大小从2到60变化.图6为不同子域重叠规模大小的CPU计算时间与内存使用量相对百分比.由图 5可知，不同子域组合的趋势基本一致，CPU计算时间随子域重叠规模大小呈线性增加；内存使用量相对量也随子域重叠规模大小呈线性增加.这是由于子域采用LU直接分解算法，对于子域间参数传递只有虚拟边界节点进行.因此，在实际求解过程中，子域重叠规模大小最佳选择为2.

图 5

Figure 5

图 5 不同子域重叠规模大小的内存和CPU计算时间变化规律 (a)CPU计算时间；(b)与全域求解法相比的内存使用量相对比值. Figure 5 Variation about CPU time and relative memory usage of sub-domain overlap size

本文将大地电磁二维求解区域分解成若干重叠子域，通过Schwarz交替法构建子域虚拟边界，对于每个子域采用有限差分算法进行离散并采用LU直接分解算法进行独立求解，然后通过全域迭代实现大地电磁法二维正演模拟.通过典型低阻模型的数值模拟验证本文算法的准确性和可行性，并通过不同子域分解数目、子域分解组合及子域重叠规模的数值计算表明，重叠区域分解算法相对于传统有限差分算法，内存使用量和CPU计算时间都大大减少，在实际求解过程中子域重叠大小选择2为佳并选择合适子域数目和子域分解组合.因此，本文提出的重叠区域分解算法能够将大规模转化成若干小规模问题进行求解，为大规模大尺度电磁三维正演问题实现提供一种新的途径.本文程序设计存在不足，CPU计算时间结果仅供参考.