多重網(wǎng)格法二維Helmholtz方程解算及其在電磁法正演模擬中的應(yīng)用

楊振威 馮 磊 趙 寧 趙秋芳 楊雙安

（①河南理工大學(xué)資源環(huán)境學(xué)院，河南焦作454000；②中原經(jīng)濟(jì)區(qū)煤層（頁巖）氣河南省協(xié)同創(chuàng)新中心，河南焦作454000；③河南理工大學(xué)計(jì)算機(jī)學(xué)院，河南焦作454000）

·非地震·

多重網(wǎng)格法二維Helmholtz方程解算及其在電磁法正演模擬中的應(yīng)用

楊振威①②馮 磊①②趙 寧*③趙秋芳①②楊雙安①②

（①河南理工大學(xué)資源環(huán)境學(xué)院，河南焦作454000；②中原經(jīng)濟(jì)區(qū)煤層（頁巖）氣河南省協(xié)同創(chuàng)新中心，河南焦作454000；③河南理工大學(xué)計(jì)算機(jī)學(xué)院，河南焦作454000）

為了提高Helmholtz方程數(shù)值計(jì)算效率和精度，研究了多重網(wǎng)格算法，并對(duì)比研究了該算法與共軛梯度法、預(yù)處理共軛梯度法和超松弛法求解二維Helmholtz方程的計(jì)算精度和收斂速度，網(wǎng)格剖分采用可實(shí)現(xiàn)網(wǎng)格自動(dòng)細(xì)化的Delaunay三角網(wǎng)格算法。研究結(jié)果表明：多重網(wǎng)格法在計(jì)算時(shí)間和迭代收斂效率方面具有較大優(yōu)勢(shì)，但其迭代計(jì)算誤差大于其他算法，這或許與不規(guī)則網(wǎng)格剖分導(dǎo)致網(wǎng)格層間插值、限制算子擴(kuò)大了計(jì)算誤差有關(guān)。最后，初步研究了基于多重網(wǎng)格算法的大地電磁二維正演模擬響應(yīng)。

多重網(wǎng)格 Helmholtz方程 大地電磁 共軛梯度

1 引言

大地電磁三維正、反演是計(jì)算地球物理研究的熱點(diǎn)和難點(diǎn)，其面臨的一個(gè)重要問題是計(jì)算量大、效率低，且數(shù)值解的精度不高。大地電磁正演模擬基于Helmholtz方程，開展快速求解Helmholtz方程的數(shù)值算法研究，對(duì)于提高正、反演計(jì)算速度具有重要的意義［1-3］；多重網(wǎng)格法（MG）在不同尺度的網(wǎng)格層上迭代求解線性方程組，具有收斂速度快、計(jì)算效率高的特點(diǎn)，對(duì)于計(jì)算地球物理正演問題具有一定的研究價(jià)值［4-7］。

Mulder［8］采用多重網(wǎng)格算法研究了時(shí)諧電磁場的擴(kuò)散規(guī)律；Plaks等［9］研究了多重網(wǎng)格算法在地球物理研究中的開邊界問題；Zaslavsky等［10］應(yīng)用自適應(yīng)多重網(wǎng)格技術(shù)進(jìn)行分子生物學(xué)模擬；Tang等［11］研究了基于自適應(yīng)多重網(wǎng)格有限元法的三維直流電正演模擬，證明該算法相對(duì)于解析解的誤差小于1%；柳建新等［12］和 Mitsuhata等［13］應(yīng)用廣義傅里葉譜分析了基于多重網(wǎng)格算法的一維大地電磁正演計(jì)算的收斂性。

上述研究在MG法及其計(jì)算地球物理領(lǐng)域的應(yīng)用取得了一定成果，但也存在不足之處，即未詳盡地研究MG法的計(jì)算效率及收斂性。本文在研究了多重網(wǎng)格法求解二維Helmholtz方程的迭代次數(shù)、運(yùn)算時(shí)間和誤差的基礎(chǔ)上，初步研究了該算法在大地電磁測(cè)深正演模擬中的應(yīng)用。

2 多重網(wǎng)格算法

MG法可以分為代數(shù)多重網(wǎng)格法和幾何多重網(wǎng)格法。多重網(wǎng)格法源自迭代求解線性方程時(shí)，誤差分量的傅里葉分量可以分為頻率高、變化快的高頻分量和頻率低、變化慢的光滑分量，而高頻分量和光滑分量是相對(duì)于網(wǎng)格尺度而言的，在不同尺度的網(wǎng)格層上迭代求解方程，將在細(xì)網(wǎng)格層上平滑后的誤差準(zhǔn)確地投射到粗網(wǎng)格層上，可以達(dá)到快速收斂的效果，從而減小計(jì)算量，提高計(jì)算效率。近年來，該方法在流體力學(xué)、結(jié)構(gòu)力學(xué)等涉及大型稀疏矩陣計(jì)算的領(lǐng)域得到廣泛應(yīng)用［14-17］。

MG法的迭代計(jì)算包含以下過程：

①預(yù)處理，在細(xì)網(wǎng)格D h上進(jìn)行若干次松弛迭代，即A hu h=f h；

②殘差估計(jì)，rh=f h-A hu h；

④粗網(wǎng)格方程求解，

⑥循環(huán)迭代，直至最大迭代次數(shù)或達(dá)到誤差限。

上述計(jì)算過程中，D為求解區(qū)域，h為細(xì)網(wǎng)格層的網(wǎng)格尺寸，A為系數(shù)矩陣，A h為細(xì)網(wǎng)格層系數(shù)矩陣，f h為細(xì)網(wǎng)格層上矩陣方程的右端項(xiàng)，R為從粗網(wǎng)格層至細(xì)網(wǎng)格層的限制算子，I為從細(xì)網(wǎng)格層至粗網(wǎng)格層的插值算子。

圖1a為多重網(wǎng)格迭代過程示意圖，G代表松馳迭代算法，N代表迭代次數(shù)，一般N=3～5，H代表最粗網(wǎng)格層網(wǎng)格尺寸，uih、G ih、eih和f ih（i=0，1，2，…）分別是定義在第ih網(wǎng)格層上的迭代解、迭代矩陣、迭代誤差和右端項(xiàng)。經(jīng)過若干次迭代，獲得滿足精度要求的迭代近似解。圖1b是層間網(wǎng)格解及誤差傳遞示意圖，圖中紅黑圓點(diǎn)表示細(xì)網(wǎng)格層的迭代解，紅黑圓圈表示粗網(wǎng)格層上的迭代解，該解可由細(xì)網(wǎng)格解插值生成。

圖1 MG法迭代計(jì)算示意圖

網(wǎng)格剖分采用Delaunay三角剖分算法，初次網(wǎng)格剖分后再進(jìn)行5次左右的網(wǎng)格自動(dòng)加密處理（如圖2所示），進(jìn)而提高解的精度，但同時(shí)計(jì)算量和計(jì)算時(shí)間呈指數(shù)級(jí)增長。

圖2 Delaunay三角網(wǎng)格剖分示意圖

3 二維Helmholtz方程的多重網(wǎng)格法

在研究區(qū)域D及邊界?D上的二維Helmholtz方程邊值問題

上述邊值問題與下列變分問題等價(jià)

式中：F（u）表示u的泛函；c為未知系數(shù)。

網(wǎng)格剖分是進(jìn)行有限元解算的關(guān)鍵步驟，網(wǎng)格剖分密度和質(zhì)量對(duì)有限元法計(jì)算精度、計(jì)算效率和收斂性起著十分重要的作用。本次網(wǎng)格剖分基于Delaunay三角剖分算法［13］，為了提高解的精度，在7個(gè)網(wǎng)格層上迭代求解，網(wǎng)格數(shù)量逐層倍增，最底層網(wǎng)格數(shù)量約2萬個(gè)。

在區(qū)域D上對(duì)Helmholtz方程三角網(wǎng)格離散化，c2=0.3，D=［-1，1］2，分別采用 MG法、共軛梯度法（CG）、預(yù)處理共軛梯度法（PCG）及超松弛法（SOR）求解離散化后的有限元線性方程組AU=f，U表示方程組的解，f表示右端項(xiàng)。

數(shù)值計(jì)算采用Matlab軟件編程并在PC機(jī)（CPU2.2GHz，內(nèi)存2GB）雙精度型下進(jìn)行，在 MG法計(jì)算過程中，采用不規(guī)則三角網(wǎng)格V循環(huán)、七層網(wǎng)格套迭代。當(dāng)最精細(xì)網(wǎng)格層上的誤差殘量的L2范數(shù)小于控制收斂準(zhǔn)則ε（ε取1.0×10-4）時(shí)，計(jì)算即為收斂。圖3為基于MG法和SOR法的二維Helmholtz方程解的曲面圖。

圖3 Helmholtz方程解曲面圖

4 計(jì)算結(jié)果分析

4.1 誤差分析

為了估計(jì)多重網(wǎng)格法的計(jì)算精度，采用歐幾里得范數(shù)對(duì)迭代解進(jìn)行誤差分析。誤差分析公式如下

式中：el為第l次迭代后的相對(duì)誤差，即本次迭代近似解與上次迭代近似解的歐氏范數(shù)；α為網(wǎng)格節(jié)點(diǎn)序號(hào)。

為了深入分析MG法的計(jì)算精度和收斂性，對(duì)比研究了該算法與CG法、PCG法、SOR法的相對(duì)誤差，見表1。由表可知，無論網(wǎng)格單元個(gè)數(shù)多與少，該算法在計(jì)算精度方面均不具有優(yōu)勢(shì)，或者說其計(jì)算誤差大于其他三種算法，這或許與其在多個(gè)網(wǎng)格層進(jìn)行不規(guī)則三角網(wǎng)格剖分有關(guān)，具體原因有待進(jìn)一步研究。

4.2 收斂性分析

圖4是網(wǎng)格數(shù)分別為11425和45377個(gè)時(shí)，MG和PCG法迭代的收斂曲線。從圖中不難看出，MG法的收斂速度明顯優(yōu)于PCG法，說明MG法的收斂速度與網(wǎng)格尺度無關(guān)，且初始收斂速度明顯快于PCG法。隨著迭代次數(shù)的增加，誤差減小，這兩種算法的收斂速度差距逐步縮小。

表1 不同方法計(jì)算相對(duì)誤差對(duì)比

圖4 MG法（a）和PCG法（b）收斂曲線對(duì)比圖

4.3 計(jì)算時(shí)間分析

計(jì)算時(shí)間是衡量算法計(jì)算效率的重要參考指標(biāo)。在設(shè)置相同誤差精度的條件下，采用Matlab中計(jì)算程序運(yùn)行時(shí)間的函數(shù)（tic／toc），分別得到 MG法、CG法、SOR法和PCG法的計(jì)算時(shí)間如表2所示。從表中可以看出，當(dāng)網(wǎng)格單元個(gè)數(shù)較少時(shí)，MG法在計(jì)算時(shí)間上優(yōu)勢(shì)并不明顯；當(dāng)網(wǎng)格個(gè)數(shù)較多時(shí)（如722177），MG法的計(jì)算時(shí)間分別是SOR法的近1／10、CG法的約1／2、PCG法的約1／6，具有明顯優(yōu)勢(shì)。

表2 不同算法計(jì)算時(shí)間對(duì)比

5 基于多重網(wǎng)格法的大地電磁二維正演模擬

5.1 邊值問題

為了檢驗(yàn)MG法在大地電磁二維正演模擬中的計(jì)算特性，初步研究了均勻半空間條件下二維地電模型的大地電磁正演響應(yīng)。大地電磁測(cè)深法二維邊值問題可以表示為

對(duì)于上邊界，u=1；對(duì)于側(cè)邊界，對(duì)于底邊界，

5.2 限制算子

MG法的關(guān)鍵是殘量和近似解在網(wǎng)格層間進(jìn)行精確傳遞，即合理設(shè)計(jì)粗—細(xì)網(wǎng)格層間插值和限制算子。

圖5 限制算子設(shè)計(jì)（面積率）示意圖

限制算子的作用是將細(xì)網(wǎng)格層經(jīng)過松弛迭代后的殘差投射至較粗網(wǎng)格層，粗網(wǎng)格層上每一個(gè)網(wǎng)格節(jié)點(diǎn)及其周圍，有9個(gè)細(xì)網(wǎng)格層節(jié)點(diǎn)與其對(duì)應(yīng)（圖5），殘量由細(xì)網(wǎng)格層投射至粗網(wǎng)格層過程中，細(xì)網(wǎng)格層上節(jié)點(diǎn)對(duì)粗網(wǎng)格對(duì)應(yīng)節(jié)點(diǎn)的“貢獻(xiàn)”不同，當(dāng)前，限制算子較多的采用完全加權(quán)算子，其設(shè)計(jì)思想一般采用面積率，即將某一粗網(wǎng)格進(jìn)行16等分，設(shè)a為每個(gè)粗網(wǎng)格面積的1／16，則其上每一細(xì)網(wǎng)格節(jié)點(diǎn)對(duì)粗網(wǎng)格的“貢獻(xiàn)”可以歸納為：a0=4a，a1=a2=a3=a4=2a，a5=a6=a7=a8=a，a0、a1、a2、a3、a4、

a5、a6、a7、a8分別表示粗網(wǎng)格周圍對(duì)應(yīng)的細(xì)網(wǎng)格點(diǎn)，如圖5所示，基于此，完全加權(quán)算子簡化為

式中下標(biāo)i，j代表網(wǎng)格節(jié)點(diǎn)編號(hào)。式（6）用矩陣的形式表示為

5.3 插值算子

插值算子將解的修正量從粗網(wǎng)格傳遞至細(xì)網(wǎng)格，算子構(gòu)建原則是與粗網(wǎng)格相對(duì)應(yīng)的細(xì)網(wǎng)格點(diǎn)值保持不變，其余細(xì)網(wǎng)格點(diǎn)值由其相鄰粗網(wǎng)格點(diǎn)值進(jìn)行算術(shù)平均得到，即

式中：v為網(wǎng)格節(jié)點(diǎn)解的近似量；上標(biāo)h、2h分別表示細(xì)、粗網(wǎng)格層?；谑剑?），可得 MG法的粗、細(xì)網(wǎng)格層間插值算子

5.4 模型計(jì)算

設(shè)計(jì)模型電阻率為10Ω·m，向下極限深度為1.0×105m，考慮到Gauss-Seidel等一般迭代算法在計(jì)算過程中存在收斂性不穩(wěn)定的缺陷，本次大地電磁二維正演模擬采用了收斂性較穩(wěn)定的雙共軛梯度法（Bicgstab），多重（二重）網(wǎng)格算法的細(xì)網(wǎng)格松弛迭代算法也采用Bicgstab算法。正演模型采用二維均勻半空間模型，橫、縱向網(wǎng)格數(shù)均為20，為提高計(jì)算精度，設(shè)計(jì)單元網(wǎng)格長為50，寬為10，最大迭代次數(shù)為10000，誤差限為1.0×10-20。分別采用Bicgstab算法和多重網(wǎng)格算法（細(xì)網(wǎng)格松弛迭代采用Bicgstab算法）對(duì)模型進(jìn)行計(jì)算，模型響應(yīng)如圖6所示。由圖可知，Bicgstab算法在計(jì)算精度方面較多重網(wǎng)格算法有一定優(yōu)勢(shì)。

圖6 大地電磁二維正演模型響應(yīng)示意圖

Bicgstab算法的計(jì)算時(shí)間為8.003234s，MG算法的計(jì)算時(shí)間為4.01562s，由此可見，MG算法在計(jì)算時(shí)間方面有明顯優(yōu)勢(shì)。

6 結(jié)論和認(rèn)識(shí)

本文采用 MG、SOR、PCG法，研究了二維Helmholtz方程的Dirichlet邊值問題的計(jì)算精度和收斂速度，并在此基礎(chǔ)上，初步研究了均勻半空間條件下基于MG法的大地電磁二維正演模擬響應(yīng)，取得了如下認(rèn)識(shí)。

（1）MG法的迭代收斂速度與網(wǎng)格尺度無關(guān)，且最初幾次迭代的收斂速度明顯優(yōu)于其他算法。

（2）MG法在網(wǎng)格單元個(gè)數(shù)較少時(shí)，相同誤差限條件下的計(jì)算時(shí)間優(yōu)勢(shì)并不明顯；在網(wǎng)格單元個(gè)數(shù)較多時(shí)，計(jì)算時(shí)間最多只是其他算法的一半；

（3）本次MG法的網(wǎng)格剖分采用不規(guī)則三角網(wǎng)格算法，導(dǎo)致采用網(wǎng)格層間插值和限制算子放大了計(jì)算誤差，進(jìn)而使該算法的計(jì)算誤差較其他三種算法大，如何降低誤差有待進(jìn)一步研究。

［1］ 劉小軍，王家林，陳冰等.二維大地電磁數(shù)據(jù)的聚焦反演算法探討.石油地球物理勘探，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］ 王青平，白武明，王洪亮.多重網(wǎng)格在二維泊松方程有限元分析中的應(yīng)用.地球物理學(xué)進(jìn)展，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］ 何繼善，李帝銓，戴世坤.廣域電磁法在湘西北頁巖氣探測(cè)中的應(yīng)用.石油地球物理勘探，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］ 胡祖志，陳英，何展翔等.大地電磁并行模擬退火約束反演及應(yīng)用.石油地球物理勘探，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.三維泊松方程數(shù)值模擬的多重網(wǎng)格方法.地球物理學(xué)進(jìn)展，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］ 韓波，胡祥云，何展翔等.大地電磁反演方法的數(shù)學(xué)分類.石油地球物理勘探，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］ 柳建新，郭榮文，童孝忠.基于多重網(wǎng)格法的MT正演模型邊界截取.中南大學(xué)學(xué)報(bào)（自然科學(xué)版），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三角網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)與遍歷.天津理工大學(xué)學(xué)報(bào)，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.

P631

A

10.13810／j.cnki.issn.1000-7210.2017.01.023

楊振威，馮磊，趙寧，趙秋芳，楊雙安.多重網(wǎng)格法二維Helmholtz方程解算及其在電磁法正演模擬中的應(yīng)用.石油地球物理勘探，2017，52（1）：167-172.

1000-7210（2017）01-0167-06

*河南省焦作市高新區(qū)世紀(jì)路2001號(hào)，454000。Email：zhaoning＠hpu.edu.cn

本文于2015年10月9日收到，最終修改稿于2016年11月4日收到。

本項(xiàng)研究受河南省教育廳重點(diǎn)基金項(xiàng)目（15A170008）和河南省博士后基金項(xiàng)目聯(lián)合資助。

（本文編輯：馮杏芝）

楊振威 博士，講師，1984年生；2007年本科畢業(yè)于河南理工大學(xué)獲地理信息系統(tǒng)專業(yè)學(xué)士學(xué)位，2010年畢業(yè)于中國礦業(yè)大學(xué)（北京），獲礦產(chǎn)普查與勘探專業(yè)碩士學(xué)位，2013年畢業(yè)于中國地質(zhì)科學(xué)院，獲地球探測(cè)與信息技術(shù)專業(yè)博士學(xué)位；現(xiàn)在河南理工大學(xué)資源環(huán)境學(xué)院從事地球信息科學(xué)與技術(shù)專業(yè)的教研。