多重网格法二维Helmholtz方程解算及其在电磁法正演模拟中的应用

2017-10-23杨振威赵秋芳杨双安

石油地球物理勘探 2017年1期

关键词：网格法计算精度剖分

杨振威冯磊赵宁赵秋芳杨双安

（①河南理工大学资源环境学院，河南焦作454000；②中原经济区煤层（页岩）气河南省协同创新中心，河南焦作454000；③河南理工大学计算机学院，河南焦作454000）

·非地震·

多重网格法二维Helmholtz方程解算及其在电磁法正演模拟中的应用

杨振威①②冯磊①②赵宁*③赵秋芳①②杨双安①②

为了提高Helmholtz方程数值计算效率和精度，研究了多重网格算法，并对比研究了该算法与共轭梯度法、预处理共轭梯度法和超松弛法求解二维Helmholtz方程的计算精度和收敛速度，网格剖分采用可实现网格自动细化的Delaunay三角网格算法。研究结果表明：多重网格法在计算时间和迭代收敛效率方面具有较大优势，但其迭代计算误差大于其他算法，这或许与不规则网格剖分导致网格层间插值、限制算子扩大了计算误差有关。最后，初步研究了基于多重网格算法的大地电磁二维正演模拟响应。

多重网格 Helmholtz方程大地电磁共轭梯度

1 引言

大地电磁三维正、反演是计算地球物理研究的热点和难点，其面临的一个重要问题是计算量大、效率低，且数值解的精度不高。大地电磁正演模拟基于Helmholtz方程，开展快速求解Helmholtz方程的数值算法研究，对于提高正、反演计算速度具有重要的意义［1-3］；多重网格法（MG）在不同尺度的网格层上迭代求解线性方程组，具有收敛速度快、计算效率高的特点，对于计算地球物理正演问题具有一定的研究价值［4-7］。

Mulder［8］采用多重网格算法研究了时谐电磁场的扩散规律；Plaks等［9］研究了多重网格算法在地球物理研究中的开边界问题；Zaslavsky等［10］应用自适应多重网格技术进行分子生物学模拟；Tang等［11］研究了基于自适应多重网格有限元法的三维直流电正演模拟，证明该算法相对于解析解的误差小于1%；柳建新等［12］和 Mitsuhata等［13］应用广义傅里叶谱分析了基于多重网格算法的一维大地电磁正演计算的收敛性。

上述研究在MG法及其计算地球物理领域的应用取得了一定成果，但也存在不足之处，即未详尽地研究MG法的计算效率及收敛性。本文在研究了多重网格法求解二维Helmholtz方程的迭代次数、运算时间和误差的基础上，初步研究了该算法在大地电磁测深正演模拟中的应用。

2 多重网格算法

MG法可以分为代数多重网格法和几何多重网格法。多重网格法源自迭代求解线性方程时，误差分量的傅里叶分量可以分为频率高、变化快的高频分量和频率低、变化慢的光滑分量，而高频分量和光滑分量是相对于网格尺度而言的，在不同尺度的网格层上迭代求解方程，将在细网格层上平滑后的误差准确地投射到粗网格层上，可以达到快速收敛的效果，从而减小计算量，提高计算效率。近年来，该方法在流体力学、结构力学等涉及大型稀疏矩阵计算的领域得到广泛应用［14-17］。

MG法的迭代计算包含以下过程：

①预处理，在细网格D h上进行若干次松弛迭代，即A hu h=f h；

②残差估计，rh=f h-A hu h；

④粗网格方程求解，

⑥循环迭代，直至最大迭代次数或达到误差限。

上述计算过程中，D为求解区域，h为细网格层的网格尺寸，A为系数矩阵，A h为细网格层系数矩阵，f h为细网格层上矩阵方程的右端项，R为从粗网格层至细网格层的限制算子，I为从细网格层至粗网格层的插值算子。

图1a为多重网格迭代过程示意图，G代表松驰迭代算法，N代表迭代次数，一般N=3～5，H代表最粗网格层网格尺寸，uih、G ih、eih和f ih（i=0，1，2，…）分别是定义在第ih网格层上的迭代解、迭代矩阵、迭代误差和右端项。经过若干次迭代，获得满足精度要求的迭代近似解。图1b是层间网格解及误差传递示意图，图中红黑圆点表示细网格层的迭代解，红黑圆圈表示粗网格层上的迭代解，该解可由细网格解插值生成。

图1 MG法迭代计算示意图

网格剖分采用Delaunay三角剖分算法，初次网格剖分后再进行5次左右的网格自动加密处理（如图2所示），进而提高解的精度，但同时计算量和计算时间呈指数级增长。

图2 Delaunay三角网格剖分示意图

3 二维Helmholtz方程的多重网格法

在研究区域D及边界∂D上的二维Helmholtz方程边值问题

上述边值问题与下列变分问题等价

式中：F（u）表示u的泛函；c为未知系数。

网格剖分是进行有限元解算的关键步骤，网格剖分密度和质量对有限元法计算精度、计算效率和收敛性起着十分重要的作用。本次网格剖分基于Delaunay三角剖分算法［13］，为了提高解的精度，在7个网格层上迭代求解，网格数量逐层倍增，最底层网格数量约2万个。

在区域D上对Helmholtz方程三角网格离散化，c2=0.3，D=［-1，1］2，分别采用 MG法、共轭梯度法（CG）、预处理共轭梯度法（PCG）及超松弛法（SOR）求解离散化后的有限元线性方程组AU=f，U表示方程组的解，f表示右端项。

数值计算采用Matlab软件编程并在PC机（CPU2.2GHz，内存2GB）双精度型下进行，在 MG法计算过程中，采用不规则三角网格V循环、七层网格套迭代。当最精细网格层上的误差残量的L2范数小于控制收敛准则ε（ε取1.0×10-4）时，计算即为收敛。图3为基于MG法和SOR法的二维Helmholtz方程解的曲面图。

图3 Helmholtz方程解曲面图

4 计算结果分析

4.1 误差分析

为了估计多重网格法的计算精度，采用欧几里得范数对迭代解进行误差分析。误差分析公式如下

式中：el为第l次迭代后的相对误差，即本次迭代近似解与上次迭代近似解的欧氏范数；α为网格节点序号。

为了深入分析MG法的计算精度和收敛性，对比研究了该算法与CG法、PCG法、SOR法的相对误差，见表1。由表可知，无论网格单元个数多与少，该算法在计算精度方面均不具有优势，或者说其计算误差大于其他三种算法，这或许与其在多个网格层进行不规则三角网格剖分有关，具体原因有待进一步研究。

4.2 收敛性分析

图4是网格数分别为11425和45377个时，MG和PCG法迭代的收敛曲线。从图中不难看出，MG法的收敛速度明显优于PCG法，说明MG法的收敛速度与网格尺度无关，且初始收敛速度明显快于PCG法。随着迭代次数的增加，误差减小，这两种算法的收敛速度差距逐步缩小。

表1 不同方法计算相对误差对比

图4 MG法（a）和PCG法（b）收敛曲线对比图

4.3 计算时间分析

计算时间是衡量算法计算效率的重要参考指标。在设置相同误差精度的条件下，采用Matlab中计算程序运行时间的函数（tic／toc），分别得到 MG法、CG法、SOR法和PCG法的计算时间如表2所示。从表中可以看出，当网格单元个数较少时，MG法在计算时间上优势并不明显；当网格个数较多时（如722177），MG法的计算时间分别是SOR法的近1／10、CG法的约1／2、PCG法的约1／6，具有明显优势。

表2 不同算法计算时间对比

5 基于多重网格法的大地电磁二维正演模拟

5.1 边值问题

为了检验MG法在大地电磁二维正演模拟中的计算特性，初步研究了均匀半空间条件下二维地电模型的大地电磁正演响应。大地电磁测深法二维边值问题可以表示为

对于上边界，u=1；对于侧边界，对于底边界，

5.2 限制算子

MG法的关键是残量和近似解在网格层间进行精确传递，即合理设计粗—细网格层间插值和限制算子。

图5 限制算子设计（面积率）示意图

限制算子的作用是将细网格层经过松弛迭代后的残差投射至较粗网格层，粗网格层上每一个网格节点及其周围，有9个细网格层节点与其对应（图5），残量由细网格层投射至粗网格层过程中，细网格层上节点对粗网格对应节点的“贡献”不同，当前，限制算子较多的采用完全加权算子，其设计思想一般采用面积率，即将某一粗网格进行16等分，设a为每个粗网格面积的1／16，则其上每一细网格节点对粗网格的“贡献”可以归纳为：a0=4a，a1=a2=a3=a4=2a，a5=a6=a7=a8=a，a0、a1、a2、a3、a4、

a5、a6、a7、a8分别表示粗网格周围对应的细网格点，如图5所示，基于此，完全加权算子简化为

式中下标i，j代表网格节点编号。式（6）用矩阵的形式表示为

5.3 插值算子

插值算子将解的修正量从粗网格传递至细网格，算子构建原则是与粗网格相对应的细网格点值保持不变，其余细网格点值由其相邻粗网格点值进行算术平均得到，即

式中：v为网格节点解的近似量；上标h、2h分别表示细、粗网格层。基于式（7），可得 MG法的粗、细网格层间插值算子

5.4 模型计算

设计模型电阻率为10Ω·m，向下极限深度为1.0×105m，考虑到Gauss-Seidel等一般迭代算法在计算过程中存在收敛性不稳定的缺陷，本次大地电磁二维正演模拟采用了收敛性较稳定的双共轭梯度法（Bicgstab），多重（二重）网格算法的细网格松弛迭代算法也采用Bicgstab算法。正演模型采用二维均匀半空间模型，横、纵向网格数均为20，为提高计算精度，设计单元网格长为50，宽为10，最大迭代次数为10000，误差限为1.0×10-20。分别采用Bicgstab算法和多重网格算法（细网格松弛迭代采用Bicgstab算法）对模型进行计算，模型响应如图6所示。由图可知，Bicgstab算法在计算精度方面较多重网格算法有一定优势。

图6 大地电磁二维正演模型响应示意图

Bicgstab算法的计算时间为8.003234s，MG算法的计算时间为4.01562s，由此可见，MG算法在计算时间方面有明显优势。

6 结论和认识

本文采用 MG、SOR、PCG法，研究了二维Helmholtz方程的Dirichlet边值问题的计算精度和收敛速度，并在此基础上，初步研究了均匀半空间条件下基于MG法的大地电磁二维正演模拟响应，取得了如下认识。

（1）MG法的迭代收敛速度与网格尺度无关，且最初几次迭代的收敛速度明显优于其他算法。

（2）MG法在网格单元个数较少时，相同误差限条件下的计算时间优势并不明显；在网格单元个数较多时，计算时间最多只是其他算法的一半；

（3）本次MG法的网格剖分采用不规则三角网格算法，导致采用网格层间插值和限制算子放大了计算误差，进而使该算法的计算误差较其他三种算法大，如何降低误差有待进一步研究。

［1］刘小军，王家林，陈冰等.二维大地电磁数据的聚焦反演算法探讨.石油地球物理勘探，2007，42（3）：338-342.Liu Xiaojun，Wang Jialin，Chen Bing et al.Discussion on focus inversion algorithm of 2-D MT data.OGP，2007，42（3）：338-342.

［2］王青平，白武明，王洪亮.多重网格在二维泊松方程有限元分析中的应用.地球物理学进展，2010，25（4）：1467-1474.Wang Qingping，Bai Wuming，Wang Hongliang.Multigrid finite element analysis for 2D modeling of Poisson equation.Progress in Geophysics，2010，25（4）：1467-1474.

［3］何继善，李帝铨，戴世坤.广域电磁法在湘西北页岩气探测中的应用.石油地球物理勘探，2014，49（5）：1006-1012.He Jishan，Li Diquan，Dai Shikun.Shale gas detection with wide field electromagnetic method in Northwestern Hunan.OGP，2014，49（5）：1006-1012.

［4］胡祖志，陈英，何展翔等.大地电磁并行模拟退火约束反演及应用.石油地球物理勘探，2010，45（4）：597-601.Hu Zuzhi，Chen Ying，He Zhanxiang et al.MT parallel simulated annealing constrained inversion and its application.OGP，2010，45（4）：597-601.

［5］鲁晶津，吴小平，Klaus Spitzer.三维泊松方程数值模拟的多重网格方法.地球物理学进展，2009，24（1）：154-158.Lu Jingjin，Wu Xiaoping，Klaus Spitzer.Multigrid method for 3D modeling of Poisson equation.Progress in Geophysics，2009，24（1）：154-158.

［6］韩波，胡祥云，何展翔等.大地电磁反演方法的数学分类.石油地球物理勘探，2012，47（1）：177-187.Han Bo，Hu Xiangyun，He Zhanxiang et al.Mathematical classification of magnetotelluric inversion methods.OGP，2012，47（1）：177-187.

［7］ Airaksinen T，Heikkola E，Pennanen A et al.An algebraic multigrid based shifted-Laplacian preconditioner for the Helmholtz equation.Journal of Computational Physics，2007，226（1）：1196-1210.

［8］ Mulder W A.A multigrid solver for 3D electromagnetic diffusion.Geophysical Prospecting，2006，54（5）：633-649.

［9］ Plaks A，Tsukerman I，Painchaud S et al.Multigrid methods for open boundary problems in geophysics.IEEE Transactions on Magnetics，2000，36（4）：633-638.

［10］ Zaslavsky L Y，Schlick T.An adaptive multigrid technique for evaluating long-range forces in biomolecular simulations.Applied Mathematics and Computation，1998，97（2-3）：237-250.

［11］ Tang J T，Wang F Y，Ren Z Y et al.3-D direct current resistivity forward modeling by adaptive multigrid finite element method.Journal of Central South University，2010，17（3）：587-592.

［12］柳建新，郭荣文，童孝忠.基于多重网格法的MT正演模型边界截取.中南大学学报（自然科学版），2011，42（11）：3429-3437.Liu Jianxin，Guo Rongwen，Tong Xiaozhong.Boundary truncation of magnetotelluric modeling based on multigrid method.Journal of Central South University（Science and Technology），2011，42（11）：3429-3437.

［13］ Mitsuhata Y J，Uchida T.3D magnetotelluric modeling using the T-Ωfinite element method.Geophysics，2004，69（1）：108-119.

［14］ Moucha R，Bailey R C.An accurate and robust multigrid algorithm for 2D forward resistivity modeling.Geophysical Prospecting，2004，52（3）：197-212.

［15］ William L B.A Multigrid Tutorial.Society for Industrial and Applied Mathematics，2000.

［16］ Liu Yuanqing，Yuan Jiansheng.A finite element domain decomposition combined with algebraic multigrid method for large-scale electromagnetic field computation.IEEE Transactions on Magnetics，2006，42（4）：655-658.

［17］高晓峰.2D-Delaunay三角网格的数据结构与遍历.天津理工大学学报，2006，22（2）：66-70.Gao Xiaofeng.Data-structure and traverse of 2DDelaunay triangulation.Journal of Tianjin University of Technology，2006，22（2）：66-70.