大地電磁三維正演聚集多重網(wǎng)格算法

陳 輝，尹 敏，殷長(zhǎng)春，鄧居智

1.吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026 2.東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南昌 330013

0 引言

經(jīng)過(guò)幾十年的發(fā)展，大地電磁三維正演已基本趨于成熟[1-3]。其基本步驟一般為：先對(duì)已知電阻率或電導(dǎo)率模型的整個(gè)求解區(qū)域進(jìn)行結(jié)構(gòu)或非結(jié)構(gòu)性離散，然后利用有限差分法(FD)[4-7]、有限體積法(FV)[8-10]、有限單元法(FE)[9-15]、積分方程法(IE)[16-18]等數(shù)值方法對(duì)雙旋度磁場(chǎng)或電場(chǎng)方程或矢量和標(biāo)量勢(shì)耦合方程進(jìn)行離散，并施加邊界條件得到龐大的稀疏復(fù)數(shù)線性方程組。對(duì)于該線性方程組，通常先采用預(yù)處理方法降低矩陣的條件數(shù)(比如不完全LU分解[19]、不完全Cholesky分解[20]等)，然后采用Krylov子空間迭代方法進(jìn)行求解，如雙共軛梯度法(BICG)[21]、準(zhǔn)殘量最小化法(QMR)[22]、最小殘量法(MRM)[23]等。并利用散度校正技術(shù)能夠減小迭代次數(shù)，加快正演速度，提高計(jì)算精度，特別在低頻時(shí)可取得明顯效果[24-25]。另外，隨著矩陣排序技術(shù)和計(jì)算機(jī)性能的提高，直接求解技術(shù)也逐步應(yīng)用于線性方程組求解。2009年，Streich[26]利用MUMPS數(shù)學(xué)求解庫(kù)在并行機(jī)計(jì)算機(jī)上直接求解該線性方程組，實(shí)現(xiàn)電磁法三維正演計(jì)算。2016年，Puzyrev等[27]對(duì)各種直接求解數(shù)學(xué)庫(kù)在地球物理電磁正演模擬中的應(yīng)用效果和計(jì)算效率進(jìn)行系統(tǒng)分析。但是，直接求解算法對(duì)計(jì)算機(jī)內(nèi)存及計(jì)算機(jī)性能要求較高，傳統(tǒng)迭代算法隨著模型尺度的增大，迭代次數(shù)和求解時(shí)間急劇增加，均難以滿足大規(guī)模大地電磁三維正演問(wèn)題需求。多重網(wǎng)格(MG)法目前已發(fā)展成為大型稀疏線性方程組迭代求解的最有效方法之一，在地球物理電磁問(wèn)題中得到應(yīng)用[28-30]。該方法可分為幾何多重網(wǎng)格(GMG)法和代數(shù)多重網(wǎng)格(AMG)法。AMG法是在GMG法原理上建立起來(lái)的一種求解矩陣方程的迭代方法，它不依賴于待求解問(wèn)題的幾何和物理屬性，僅利用離散后的系數(shù)矩陣信息即可對(duì)線性方程組進(jìn)行求解[31]。為了改善AMG算法的穩(wěn)定性和適用性，2010年，Yvan Notay[32]提出了利用矩陣圖論思想的光滑聚集構(gòu)建方法以實(shí)現(xiàn)AMG的粗化(簡(jiǎn)稱AGMG)，對(duì)算法的基本原理和計(jì)算效率進(jìn)行了詳細(xì)闡述，并在對(duì)流擴(kuò)散方程和橢圓方程求解應(yīng)用中實(shí)現(xiàn)了并行計(jì)算。

本文引入新型代數(shù)多重網(wǎng)格算法——聚集多重網(wǎng)格(AGMG)算法求解大地電磁三維正演問(wèn)題，通過(guò)典型地電模型計(jì)算驗(yàn)證方法的有效性，通過(guò)對(duì)不同尺度和不同地電模型的數(shù)值模擬分析AGMG及AGMG預(yù)處理技術(shù)的計(jì)算效率，為大地電磁法三維大規(guī)模、高精度和快速正演模擬提供技術(shù)支撐。

1 大地電磁三維有限體積正演

由大地電磁法理論可知，在忽略位移電流的條件下，麥克斯韋方程組的積分形式為

(1)

式中：假定電磁場(chǎng)隨時(shí)間變化的因子為e-iωt；l為環(huán)路線積分步長(zhǎng)；S為面積分單元；μ為磁導(dǎo)率；σ為電導(dǎo)率；ω為角頻率；E為電場(chǎng)；H為磁場(chǎng)。

在笛卡爾坐標(biāo)系中采用正六面體網(wǎng)格剖分，每個(gè)單元采用Yee氏交錯(cuò)網(wǎng)格進(jìn)行采樣(圖1)。電場(chǎng)分量定義在六面體每個(gè)邊的中點(diǎn)；磁場(chǎng)分量定義為六面體的面中心，方向與面的法向一致。對(duì)方程組(1)第一式和第二式采用有限體積進(jìn)行離散，并消去磁場(chǎng)得到電場(chǎng)的線性方程組：

Ae=b。

(2)

式中：A為大型復(fù)稀疏矩陣；e為待求的電場(chǎng)分量；b為與邊界條件有關(guān)的右端項(xiàng)。邊界條件采用第一類的Dirichlet邊界條件。為了提高方程組(2)的求解速度和精度，通常在迭代求解過(guò)程中施加散度校正技術(shù)。

圖1 三維有限差分模擬網(wǎng)格Fig.1 Three-dimension finite-difference grid

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

對(duì)線性方程組(2)的求解是大地電磁三維正演的關(guān)鍵之一。本文采用AGMG算法進(jìn)行求解。該算法是一種新型的代數(shù)多重網(wǎng)格算法，通過(guò)基于矩陣圖論思想的光滑聚集實(shí)現(xiàn)矩陣粗化，提高傳統(tǒng)AMG算法的穩(wěn)定性。該方法可分為聚集粗化以及迭代求解2個(gè)階段。

2.1 聚集粗化

矩陣粗化是多重網(wǎng)格中最重要的部分。在幾何多重網(wǎng)格中利用粗細(xì)網(wǎng)格之間的幾何關(guān)系進(jìn)行粗化。代數(shù)多重網(wǎng)格脫離了網(wǎng)格幾何關(guān)系，可通過(guò)代數(shù)中的矩陣運(yùn)算實(shí)現(xiàn)粗化。假設(shè)系數(shù)矩陣A，變量數(shù)為n，我們可以構(gòu)建n×nc的延拓矩陣P，其中nc為粗化矩陣變量數(shù)。該矩陣能夠?qū)⒓?xì)網(wǎng)格上的變量集變換到粗網(wǎng)格變量集，實(shí)現(xiàn)矩陣粗化。粗化后系數(shù)矩陣Ac的變量數(shù)為nc(nc

Ac=PTAP。

(3)

同時(shí)，我們定義延拓矩陣的轉(zhuǎn)置為限制矩陣。

在AMG聚集粗化過(guò)程中，定義聚集變量集Gi(i=1，...，nc)，它在幾何意義上可認(rèn)為是網(wǎng)格粗化后的網(wǎng)格節(jié)點(diǎn)；而在代數(shù)上Gi是矩陣粗化過(guò)程中的聚集變量集。Gi是系數(shù)矩陣A的互不相交子集，其大小等于粗化矩陣變量數(shù)nc。由此，延拓矩陣P定義為

(4)

式中：i表示矩陣A的變量序號(hào)；j表示聚集變量集Gj的序號(hào)。由(4)式可以看出延拓矩陣P為布爾矩陣，每行最多只有一個(gè)非零元素。將式(4)代入式(3)可得粗化矩陣Ac表達(dá)式為

(5)

式中，i、j分別表示粗化后的系數(shù)矩陣Ac行號(hào)和列號(hào)。

在AMG算法中，粗網(wǎng)格矩陣Ac的選取直接影響算法的穩(wěn)定性和運(yùn)行效率，因此延拓矩陣P的構(gòu)建至關(guān)重要。本文采用Notay等[32]提出的成對(duì)聚集算法構(gòu)建延拓矩陣P。該算法通過(guò)矩陣圖論的正則覆蓋，沿連通最好的方向進(jìn)行聚集。對(duì)于任意一個(gè)系數(shù)矩陣A中的虛擬網(wǎng)格節(jié)點(diǎn)，通過(guò)尋找最強(qiáng)連通點(diǎn)進(jìn)行配對(duì)并成對(duì)聚集。為此，首先定義虛擬節(jié)點(diǎn)i強(qiáng)連通點(diǎn)集為

(6)

式中，β為強(qiáng)弱連接閾值，需要根據(jù)系數(shù)矩陣A的特性進(jìn)行選擇。對(duì)于橢圓方程問(wèn)題，β一般取0.25[32]。對(duì)于電磁問(wèn)題，β一般取0.4～0.6。由此，聚集變量集Gi可寫為

(7)

實(shí)際計(jì)算時(shí)，利用式(7)構(gòu)建延拓矩陣P，再通過(guò)式(5)即可實(shí)現(xiàn)矩陣粗化。

2.2 迭代求解

在多重網(wǎng)格算法中，需要一定的套迭代技術(shù)將各層網(wǎng)格或矩陣連接起來(lái)。常用的套迭代技術(shù)循環(huán)方式有V循環(huán)、W循環(huán)、FMV循環(huán)。在實(shí)際應(yīng)用中采用較多的是V或W循環(huán)。W循環(huán)能保持收斂因子不隨網(wǎng)格而變化，具有Robust性，但耗費(fèi)計(jì)算時(shí)間長(zhǎng)；當(dāng)網(wǎng)格層不多時(shí)，V循環(huán)具有與W循環(huán)相同的性質(zhì)且計(jì)算量小。因此，本文采用V循環(huán)套迭代技術(shù)。

圖2給出V循環(huán)套迭代技術(shù)流程圖。V循環(huán)從細(xì)網(wǎng)格出發(fā)，逐漸粗化到最粗網(wǎng)格，最后又返回得到細(xì)網(wǎng)格上的解。圖2中:PT表示把細(xì)網(wǎng)格上的殘余誤差限制到粗網(wǎng)格上的算子,稱為限制算子;P表示把粗網(wǎng)格上的結(jié)果插值到細(xì)網(wǎng)格上的算子，稱為插值算子;u為線性方程組迭代解;A表示某一層網(wǎng)格上的系數(shù)矩陣，而f表示針對(duì)該層網(wǎng)格線性方程組的右端項(xiàng);Gv(A,f)表示光滑迭代;γ為粗網(wǎng)格修正量；k為迭代次數(shù)。由多重網(wǎng)格法的循環(huán)流程可以看出，迭代高頻誤差通過(guò)前光滑消去，而低頻誤差通過(guò)矩陣粗化轉(zhuǎn)化為高頻誤差，在粗化矩陣上采用前光滑迭代消除；在最粗層上矩陣已經(jīng)足夠小，可采用直接求解算法進(jìn)行求解。最后，通過(guò)限制算子逐層返回并進(jìn)行后光滑處理消除殘余誤差，實(shí)現(xiàn)線性方程組的快速求解。對(duì)于光滑迭代Gv(A,f)，采用Gauss- Seidel迭代法，迭代次數(shù)通常為1～3次，既能保持AGMG算法具有良好數(shù)值穩(wěn)定性，又能很好地消除迭代過(guò)程中出現(xiàn)的高頻誤差。

圖2 AGMG算法的V循環(huán)示意圖Fig.2 V-cycle scheme in AGMG method

前人研究[33]表明，當(dāng)系數(shù)矩陣A嚴(yán)重病態(tài)時(shí)，AMG算法中利用傳統(tǒng)套迭代技術(shù)(V循環(huán)、W循環(huán)及FMV循環(huán))難以取得好的效果，解容易發(fā)散。為此，通常將AMG算法與共軛梯度(CG)類算法相結(jié)合，將AMG算法用作共軛梯度類預(yù)處理算子，有效提高了傳統(tǒng)AMG算法的計(jì)算速度和穩(wěn)定性。本文設(shè)計(jì)2種求解策略求解線性方程組：傳統(tǒng)的V循環(huán)AGMG算法(簡(jiǎn)稱AGMG)及將AGMG作為傳統(tǒng)Krylov子空間迭代算法的預(yù)處理算子。我們選用共軛梯度算法(CG)和廣義共軛殘差法(GCR)作為Krylov子空間迭代算法[34-35]，分別簡(jiǎn)稱為AGMG-CG和AGMG-GCR。

3 模型算例

為檢驗(yàn)本文實(shí)施3種聚集多重網(wǎng)格算法AGMG、AGMG-CG、AGMG-GCR的準(zhǔn)確性和計(jì)算效率，我們參照2008年MT 3D Inversion Workshop提出的Dublin Test Model 1(DTM1)地電模型[36]，設(shè)計(jì)一個(gè)典型的地電模型進(jìn)行三維數(shù)值模擬，進(jìn)而分別從計(jì)算速度、迭代次數(shù)、誤差衰減特性及計(jì)算時(shí)間等方面進(jìn)行分析，并與Kelbert等[37]開發(fā)的大地電磁三維正反演代碼ModEM進(jìn)行對(duì)比。本文所有算法均在Microsoft Visual Studio 2015開發(fā)平臺(tái)上運(yùn)行，采用Inter Visual Fortran 2016編譯，程序運(yùn)行環(huán)境為Interl(R) Core(TM) i7-6700HQ CPU @ 2.60 GHz四核八線程、16 G內(nèi)存、64位win10系統(tǒng)。

3.1 精度驗(yàn)證

為了驗(yàn)證本文算法精度，設(shè)計(jì)一個(gè)含有3個(gè)異常體的復(fù)雜地電模型。圍巖電阻率為100 Ω·m；異常體ρ1為一個(gè)4 000 m×500 m×1 600 m的六面體，頂部埋深為500 m，電阻率為10 Ω·m；異常體ρ2為一個(gè)1 600 m×2 400 m×500 m的六面體，頂部埋深為2 100 m，電阻率為1 Ω·m；異常體ρ3為一個(gè)1 600 m×2 400 m×2 900 m的六面體，頂部埋深為2 100 m，電阻率為1 000 Ω·m。圖3給出該模型在yz方向的剖面圖和xy方向的平面圖。在22.6 km×30.3 km×20.6 km的求解區(qū)域內(nèi)，我們采用不等距進(jìn)行網(wǎng)格剖分，剖分網(wǎng)格數(shù)為108×114×78。計(jì)算頻率設(shè)為0.01 Hz，在求解區(qū)域正中心x方向上觀測(cè)40個(gè)測(cè)點(diǎn)，間距為400 m。在完成x和y極化模式下的大地電磁場(chǎng)三維正演模擬后計(jì)算得到大地電磁阻抗。

圖4為復(fù)雜地電模型不同求解算法的阻抗Zxy

a.yz剖面圖；b.xy平面圖。圖3 復(fù)雜地電模型示意圖Fig.3 Sketch map of complex geoelectrical model

圖4 不同求解算法的大地電磁三維正演結(jié)果對(duì)比Fig.4 Comparison of 3D MT numerical solutions for various iterative methods

和Zyx三維模擬結(jié)果對(duì)比。數(shù)值模擬分別采用本文實(shí)施的3種AGMG求解策略和Kelbert等[37]開發(fā)的大地電磁三維正演模擬程序ModEM求解，其中ModEM采用的迭代求解算法為QMR。由圖4可見(jiàn)，本文采用的3種不同AGMG求解算法(AGMG、AGMG-CG和AGMG-GCR)與ModEM計(jì)算的阻抗Zxy和Zyx實(shí)、虛部吻合很好，驗(yàn)證了算法的準(zhǔn)確性。

3.2 AGMG算法效率分析

為了分析3種不同AGMG算法的計(jì)算效率，對(duì)圖3中的模型設(shè)計(jì)4種不同規(guī)模的網(wǎng)格，分別為36×38×26、72×76×52、108×114×78、144×152×104，然后分別采用AGMG、AGMG-CG、AGMG-GCR和ModEM(QMR)對(duì)該模型進(jìn)行x和y極化模式大地電磁三維正演模擬。求解區(qū)域?yàn)?2.6 km×30.3 km×20.6 km，計(jì)算頻率為0.01 Hz，求解相對(duì)誤差設(shè)為10-8。

圖5為x極化模式下不同迭代方法的誤差衰減曲線對(duì)比圖。其中，AGMG算法在網(wǎng)格大小為36×38×26時(shí)聚集粗化層數(shù)為6層，在72×76×52時(shí)聚集粗化層數(shù)為7層，在108×114×78時(shí)聚集粗化層數(shù)為8層，在144×152×104時(shí)聚集粗化層數(shù)為10層。從圖5c中可以看出：傳統(tǒng)QMR迭代算法在迭代早期呈現(xiàn)快速衰減，迭代誤差衰減曲線光滑性較差，同時(shí)隨著迭代次數(shù)增加衰減速度變慢；AGMG算法起始呈現(xiàn)快速衰減，隨著迭代次數(shù)增加衰減速度變緩，但總體比QMR下降速度快，而且迭代誤差衰減曲線較光滑；AGMG預(yù)處理的Krylov子空間迭代算法(AGMG-CG、AGMG-GCR)呈現(xiàn)近似線性下降特征，隨著迭代次數(shù)增加衰減速度基本保持不變；另外，AGMG-GCR相對(duì)于AGMG-CG衰減速度更快，迭代誤差曲線更光滑。綜合對(duì)比不同網(wǎng)格尺度(圖5a—d)的迭代衰減曲線可以看出，QMR和AGMG衰減速度隨網(wǎng)格數(shù)增大衰減速度逐漸變慢，迭代次數(shù)隨網(wǎng)格數(shù)增大急劇增加，而AGMG-GCR和AGMG-CG隨網(wǎng)格數(shù)增大均保持近似線性快速下降，迭代次數(shù)隨網(wǎng)格數(shù)增大增加較為緩慢，同時(shí)AGMG-GCR迭代速度優(yōu)于AGMG-CG，迭代誤差曲線更加光滑。

E. 最終迭代相對(duì)誤差；NI. 迭代次數(shù)。a.36×38×26；b.72×76×52；c.108×114×78；d.144×152×104。圖5 x極化模式不同求解算法的大地電磁三維正演誤差衰減曲線Fig.5 Convergence behaviors of AGMG, AGMG-CG, AGMG-GCR and QMR for x-polarization 3D MT modelling with different grid sizes

圖6為y極化模式下不同迭代方法的誤差衰減曲線對(duì)比圖，不同網(wǎng)格的聚集粗化層數(shù)與x極化模式相同。從圖6中可以看出，4種方法的誤差衰減曲線特征與x極化模式基本一致。AGMG-GCR相對(duì)于其他3種迭代算法，誤差呈現(xiàn)近似線性下降，迭代次數(shù)隨網(wǎng)格數(shù)增加較為緩慢，誤差迭代曲線更光滑。因此，AGMG算法與傳統(tǒng)Krylov子空間迭代算法結(jié)合，不僅能夠改善AGMG的穩(wěn)定性，而且能夠加快收斂速度，而AGMG-GCR是一種更為快速穩(wěn)定的迭代求解算法。

為了進(jìn)一步對(duì)比AGMG算法的計(jì)算效率，本文對(duì)不同網(wǎng)格、不同迭代算法的迭代結(jié)果進(jìn)行統(tǒng)計(jì)，結(jié)果見(jiàn)表1。從表1可以看出，4種迭代算法在不同網(wǎng)格數(shù)下均能滿足設(shè)定的精度。在小尺度網(wǎng)格36×38×26時(shí)：QMR算法的迭代次數(shù)較多，迭代時(shí)間和正演模擬時(shí)間最少；AGMG-GCR迭代次數(shù)最少，迭代時(shí)間和正演模擬時(shí)間和QMR相當(dāng)。隨著網(wǎng)格尺度增加，4種迭代算法的求解時(shí)間和正演模擬時(shí)間急劇增加，這是由于網(wǎng)格尺度增加導(dǎo)致計(jì)算量增加的結(jié)果。但AGMG-GCR和AGMG-CG迭代次數(shù)增加緩慢，迭代時(shí)間和正演模擬時(shí)間增加的幅度遠(yuǎn)小于QMR和AGMG算法。在大尺度網(wǎng)格108×114×78時(shí)：QMR算法的迭代次數(shù)和求解時(shí)間均最大，而AGMG-GCR迭代次數(shù)和求解時(shí)間最小，相對(duì)于傳統(tǒng)QMR算法加速了十幾倍。綜上所述，將AGMG算法作為傳統(tǒng)Krylov子空間迭代法的預(yù)處理算子，在網(wǎng)格尺度較大時(shí)具有迭代次數(shù)少、計(jì)算速度快等優(yōu)點(diǎn)，同時(shí)迭代誤差呈近似線性下降，具有很好的穩(wěn)定性。隨著網(wǎng)格尺度規(guī)模增加，AGMG預(yù)處理算法迭代次數(shù)增加緩慢，加速效率變得越來(lái)越明顯。特別地，AGMG-GCR是一種快速、穩(wěn)定、高效的求解算法，為大規(guī)模問(wèn)題的大地電磁三維模擬提供技術(shù)保障。

a.36×38×26；b.72×76×52；c.108×114×78；d.144×152×104。圖6 y極化模式不同求解算法的大地電磁三維正演的誤差衰減曲線Fig.6 Convergence behaviors of AGMG, AGMG-CG, AGMG-GCR and QMR for y-polarization 3D MT modelling with different grid sizes

Table1Error,runningtimeinsecondsanditerationnumbersofAGMG,AGMG-CG,AGMG-GCRandQMRfor3DMTmodellingwithdifferentgridsizes

網(wǎng)格方法x極化y極化E/10-9NIts/sta/sE/10-9NIts/sta/s144×152×104AGMG10.006892704.26866.310.007903149.17954.7AGMG-CG14.00132638.91513.59.60154727.61726.9AGMG-GCR9.4088432.31044.09.6095459.71071.3QMR9.788695389.011207.79.9810176290.013069.4108×114×78AGMG9.90414920.41961.99.904671021.22164.3AGMG-CG10.0099284.6563.49.60102287.7562.1AGMG-GCR12.0066186.4397.89.4070190.1400.1QMR9.925321331.92898.69.957521803.73914.172×76×52AGMG10.00186234.4363.69.90242305.2486.6AGMG-CG9.806694.5140.79.5072102.1146.2AGMG-GCR10.004463.2107.110.005072.3116.5QMR9.47161121.2273.79.63239181.3391.136×38×26AGMG12.007719.331.31000.0010624.437.2AGMG-CG7.105312.520.19.305412.919.3AGMG-GCR9.70328.312.716.00338.514.4QMR9.60586.210.59.60778.313.3

注：ts. 迭代求解時(shí)間；ta. 包含散度校正的正演求解時(shí)間。

4 結(jié)論與建議

本文從準(zhǔn)靜態(tài)麥克斯韋方程出發(fā)，利用Yee氏交錯(cuò)網(wǎng)格有限體積法進(jìn)行離散，并采用第一類Dirichlet邊界條件得到關(guān)于待求解的電場(chǎng)線性方程組，然后引入新型多重網(wǎng)格算法——聚集多重網(wǎng)格算法進(jìn)行求解，實(shí)施3種不同的多重網(wǎng)格求解算法(AGMG、AGMG-CG、AGMG-GCR)。通過(guò)典型三維地電模型的正演模擬，并與已有的ModEM正反演模擬軟件對(duì)比驗(yàn)證了本文算法的準(zhǔn)確性。同時(shí)，從數(shù)值模擬結(jié)果可得出如下結(jié)論：

1)在大地電磁法三維正演中，將傳統(tǒng)AGMG算法作為Krylov子空間的迭代算法預(yù)處理算子，能夠極大改善AGMG和Krylov子空間迭代算法的穩(wěn)定性，加快誤差衰減速度、減少迭代次數(shù)，有效提升正演求解效率。

2)AGMG-GCR算法具有隨網(wǎng)格尺度增加迭代次數(shù)緩慢增加、誤差呈近似線性下降的特征，是一種高效、快速、穩(wěn)定的迭代求解算法，特別適合當(dāng)前“精細(xì)化”和“大尺度”電磁正反演問(wèn)題的求解。

致謝：本文驗(yàn)證程序?yàn)镵elbert等開發(fā)的大地電磁三維正反演代碼ModEM，聚集代數(shù)多重網(wǎng)格理論是Yvan Notay等有關(guān)多重網(wǎng)格法的研究成果。感謝Anna等和Yvan Notay提供公開源代碼及幫助的相關(guān)專家。

[1] Newman G A. A Review of High-Performance Com-putational Strategies for Modeling and Imaging of Electromagnetic Induction Data[J]. Surveys in Geophysics, 2014, 35(1): 85-100.

[2] Smith R. Electromagnetic Induction Methods in Mi-ning Geophysics from 2008 to 2012[J]. Surveys in Geophysics, 2014, 35(1): 123-156.

[3] Siripunvaraporn W. Three-Dimensional Magnetotellu-ric Inversion: An Introductory Guide for Developers and Users[J]. Surveys in Geophysics, 2012, 33(1): 5-27.

[4] 譚捍東，余欽范，Booker John，等. 大地電磁法三維交錯(cuò)采樣有限差分?jǐn)?shù)值模擬[J]. 地球物理學(xué)報(bào)，2003，46(5)：706-711.

Tan Handong, Yu Qinfan, Booker J, et al. Magnetotelluric Three-Dimensional Modelling Using the Staggered-Grid Fnite Difference Method[J]. Chinese Journal of Geophysics, 2003, 46(5): 706-711.

[5] 沈金松. 用交錯(cuò)網(wǎng)格有限差分法計(jì)算三維頻率域電磁響應(yīng)[J]. 地球物理學(xué)報(bào)，2003，46(2)：281-289.

Shen Jinsong. Modelling of 3-D Eectromagnetic Responses in Frequency Domain by Using Straggered-Grid Finite Difference Method[J]. Chinese Journal of Geophysics, 2003, 46(2): 281-289.

[6] Mackie R L, Madden T R, Wannamaker P E. Three-Dimensional Magnetotelluric Modeling Using Difference Equations: Theory and Comparisons to Integral Equation Solutions[J]. Geophysics, 1993, 58(2): 215-226.

[7] 李焱, 胡祥云, 楊文采, 等. 大地電磁三維交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬的并行計(jì)算研究[J]. 地球物理學(xué)報(bào)，2012，55(12)：4036-4043.

Li Yan, Hu Xiangyun, Yang Wencai, et al. A Study on Parallel Computation for 3D Magnetotelluric Modeling Using the Staggered-Grid Finite Difference Method[J]. Chinese Journal of Geophysics, 2012, 55(12): 4036-4043.

[8] Haber E, Ascher U M. Fast Finite Volume Simulation of 3D Electromagnetic Problems with Highly Discontinuous Coefficients[J]. SIAM Journal on Scientific Computing, 2000, 22(6): 1943-1961.

[9]Haber E, Ruthotto L. A Multiscale Finite Volume Method for Maxwell’s Equations at Low Frequencies[J]. Geophysical Journal International, 2014, 199(2): 1268-1277.

[10] 陳輝，殷長(zhǎng)春，鄧居智. 基于Lorenz規(guī)范條件下磁矢勢(shì)和標(biāo)勢(shì)耦合方程的頻率域電磁法三維正演[J]. 地球物理學(xué)報(bào)，2016，59(8)：3087-3097.

Chen Hui, Yin Changchun, Deng Juzhi. A Finite-Volume Solution to 3D Frequency-Domain Electromagnetic Modelling Using Lorenz-Gauged Magnetic Vector and Scalar Potentials[J]. Chinese Journal of Geophysics, 2016, 59(8): 3087-3097.

[11] Ren Z, Kalscheuer T, Greenhalgh S, et al. A Finite-Element-Based Domain-Decomposition Approach for Plane Wave 3D Electromagnetic Modeling[J]. Geophysics, 2014, 79(6): E255-E268.

[12] 黃臨平，戴世坤. 復(fù)雜條件下3D電磁場(chǎng)有限元計(jì)算方法[J]. 地球科學(xué)：中國(guó)地質(zhì)大學(xué)學(xué)報(bào)，2002，27(6)：775-779.

Huang Linping, Dai Shikun. Finite Element Calculation Method of 3D Electromagnetic Field under Complex Condition[J]. Earth Science: Journal of China University of Geoscineces, 2002, 27(6): 775-779.

[13] Mitsuhata Y, Uchida T. 3D Magnetotelluric Mode-ling Using the T-Omega Finite-Element Method[J]. Geophysics, 2004, 69(1): 108-119.

[14] 李俊杰，嚴(yán)家斌，皇祥宇. 無(wú)單元Galerkin法大地電磁三維正演模擬[J]. 地質(zhì)與勘探，2015，51(5)：946-952.

Li Junjie, Yan Jiabin, Huang Xiangyu. Three-Dimensional Forward Modeling of Magnetotellurics Using the Element-Free Galerkin Method[J]. Geology and Exploration, 2015, 51(5): 946-952.

[15] 嚴(yán)家斌，皇祥宇. 大地電磁三維矢量有限元正演[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2016，46(5)：1538-1549.

Yan Jiabin, Huang Xiangyu. 3D Forward Modeling of Magnetotelluric Field by Vector Finite Element Method[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(5): 1538-1549.

[16]Kruglyakov M, Geraskin A, Kuvshinov A. Novel Accurate and Scalable 3-D MT Forward Solver Based on a Contracting Integral Equation Method[J]. Computers & Geosciences, 2016, 96: 208-217.

[17]Wannamaker P E. Advances in Three-Dimensional Magnetotelluric Modeling Using Integral Equations[J]. Geophysics, 1991, 56(11): 1716-1728.

[18] 王書明，李德山，胡浩. 三維/三維構(gòu)造下大地電磁相位張量數(shù)值模擬[J]. 地球物理學(xué)報(bào)，2013，56(5)：1745-1752.

Wang Shuming, Li Deshan, Hu Hao. Numerical Modeling of Magnetotelluric Phase Tensor in the Context of 3D/3D Formation[J]. Chinese Journal of Geophysics, 2013, 56(5): 1745-1752.

[19]de Groot-Hedlin C. Finite-Difference Modeling of Magnetotelluric Felds: Error Estimates for Uniform and Nonuniform Grids[J]. Geophysics, 2006, 71(3): G225-G233.

[20] Han N, Nam M J, Kim H J, et al. A Comparison of Accuracy and Computation Time of Three-Dimensional Magnetotelluric Modelling Algorithms[J]. Journal of Geophysics and Engineering, 2009, 6(2): 136.

[21] Smith J T. Conservative Mmodeling of 3-D Electro-magnetic Fields: Part II: Biconjugate Gradient Solution and an Accelerator[J]. Geophysics, 1996, 61(5): 1319-1324.

[22] Siripunvaraporn W, Egbert G, Lenbury Y. Nume-rical Accuracy of Magnetotelluric Modelling: A Comparison of Finite Difference Approximations[J]. Earth Planets Space, 2002, 54: 721-725.

[23]Mackie R L, Madden T R. Conjugate Direction Relaxation Solutions for 3-D Magnetotelluric Modeling[J]. Geophysics, 1993, 58(7): 1052-1057.

[24] Weiss C J, Newman G A. Electromagnetic Induction in a Generalized 3D Anisotropic Earth: Part 2: The LIN Preconditioner[J]. Geophysics, 2003, 68(3): 922-930.

[25] 陳輝，鄧居智，譚捍東，等. 大地電磁三維交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬中的散度校正方法研究[J]. 地球物理學(xué)報(bào)，2011，54(6)：1649-1659.

Chen Hui, Deng Juzhi, Tan Handong, et al. Study on Divergence Correction Method in Three-Dimensional Magnetotelluric Modeling with Staggered-Grid Finite Fifference Method[J]. Chinese Journal Geophysics, 2011, 54(6): 1649-1659.

[26]Streich R. 3D Finite-Difference Frequency-Domain Modeling of Controlled-Source Electromagnetic Data: Direct Solution and Optimization for High Accuracy[J]. Geophysics, 2009, 74(5): 95-105.

[27] Puzyrev V, Koric S, Wilkin S. Evaluation of Parallel Direct Sparse Linear Solvers in Electromagnetic Geophysical Problems[J]. Computers & Geosciences, 2016, 89: 79-87.

[28] Koldan J, Puzyrev V, de la Puente J, et al. Algebraic Multigrid Preconditioning Within Parallel Finite-Element Solvers for 3-D Electromagnetic Modelling Problems in Geophysics[J]. Geophysical Journal International, 2014, 197(3): 1442-1458.

[29] Mulder W A. Geophysical Modelling of 3D Electro-magnetic Diffusion with Multigrid[J]. Computing and Visualization in Science, 2008, 11(3): 129-138.

[30]Pan K, Tang J. 2.5-D and 3-D DC Resistivity Modelling Using an Extrapolation Cascadic Multigrid Method[J]. Geophysical Journal International, 2014, 197(3): 1459-1470.

[31]Trottenberg U, Clees T. Multigrid Software for Industrial Applications: From MG00 to SAMG[M]. Heidelberg: Springer, 2009: 423-436.

[32] Notay Y. An Aggregation-Based Algebraic Multigrid Method[J]. Electronic Transactions on Numerical Analysis, 2010, 37(6): 123-146.

[33] Henson V E, Yang U M. Boomer AMG: A Parallel Algebraic Multigrid Solver and Preconditioner[J]. Applied Numerical Mathematics, 2002, 41(1): 155-177.

[34] Saad Y. Iterative Methods for Sparse Linear Systems[M]. Philadelphia: SIAM, 2003.

[35] Pflaum C. A Multigrid Conjugate Gradient Method[J]. Applied Numerical Mathematics, 2008, 58: 1803-1817.

[36]MT 3D Inversion Workshop. Dublin Test Model 1 (DTM1)[EB/OL]. [2016-10-20] http://www.complete-mt-solutions.com/mtnet/workshops/mt3di-nv/ 2008_Dublin/Dublin/3dmodel.html.

[37] Kelbert A, Meqbel N, Egbert G D, et al. ModEM: A Modular System for Inversion of Electromagnetic Geophysical Data[J]. Computers & Geosciences, 2014, 66: 40-53.