王永斐， 柳建新,3， 郭榮文,3*， 劉嶸， 李健，陳杭， 楊剛強(qiáng)

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院， 長(zhǎng)沙 410083 2 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室， 長(zhǎng)沙 410083 3 中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室， 長(zhǎng)沙 410083

0 引言

大地電磁(MT)法作為電磁(EM)法(盧杰和李予國(guó)，2019；楊海燕等，2019；底青云等，2020；何展翔等，2020)的一種，廣泛應(yīng)用于解決各種工程和科學(xué)問題，例如油氣勘探(夏訓(xùn)銀等，2012；He et al.，2015)、地?zé)嵴{(diào)查(Newman et al.，2008)、礦產(chǎn)勘探(高才坤等，2009)和深部電性結(jié)構(gòu)的研究(金勝等，2012；劉營(yíng)等，2013；王剛等，2017；崔騰發(fā)等，2020)等.它的成功應(yīng)用離不開MT數(shù)據(jù)的解釋，其中重要一環(huán)是MT數(shù)據(jù)的反演解釋.數(shù)據(jù)的反演(Avdeev, 2005; Yang et al., 2014, 2016; 周思杰和黃清華，2018；彭榮華等，2019；阮帥等，2020)是建立在正演問題的求解基礎(chǔ)上，正演求解的效率很大程度上決定了反演的效率.然而MT勘探范圍較大，比如我國(guó)的SinoProbe項(xiàng)目等(Dong et al., 2013；楊文采等，2020)，精細(xì)反演解釋需要大量模型單元，在現(xiàn)有正演算法的計(jì)算效率下，要開展全區(qū)域數(shù)據(jù)的精細(xì)反演解釋，即便是采用目前最先進(jìn)的高性能計(jì)算機(jī)也具有挑戰(zhàn)性(Zhdanov et al., 2011).因此快速M(fèi)T的正演研究一直是MT勘探的一個(gè)研究熱點(diǎn)(李焱等，2012；殷長(zhǎng)春等，2017；秦策等，2020).

由于地球深部結(jié)構(gòu)比較復(fù)雜，MT正演結(jié)果往往無法解析獲得，需要通過數(shù)值近似.常用的數(shù)值模擬方法有有限元法(FEM)(Ren et al., 2013, 2014；蔡紅柱等，2015；湯文武等，2018；顧觀文和李桐林，2020；惠哲劍等，2020；齊彥福等，2020)、積分方程法(IEM)(陳桂波等，2009；任政勇等，2017；Liu et al., 2018；Wang et al., 2019)和有限差分法(FDM)(Egbert and Kelbert, 2012；李焱等，2012；董浩等，2014；Dong and Egbert, 2019；孫懷鳳等，2021).基于電磁場(chǎng)雙旋度方程的交錯(cuò)網(wǎng)格FDM，由于其實(shí)現(xiàn)簡(jiǎn)單和靈活，已廣泛應(yīng)用于3D MT數(shù)據(jù)的反演(董浩等，2014；Kelbert et al., 2014).

MT正演的數(shù)值求解往往涉及的計(jì)算區(qū)域比較大，其正演方程的離散得到的是大型稀疏矩陣，通常很難用直接求解方法進(jìn)行求解.一般采用的是基于Krylov子空間的迭代算法進(jìn)行求解，如穩(wěn)定雙共軛梯度法(BiCGstab)(Van der Vorst, 1992)、擬最小殘差法(QMR)(Freund, 1993)等.由于離散的Maxwell方程存在豐富的零空間，常用的迭代方法很難收斂，且當(dāng)頻率趨向于零時(shí)離散方程無法保證電流散度處處為零.為了提高正演效率，通常需要采用預(yù)條件技術(shù)，如超松弛預(yù)條件(SSOR)(Chen et al., 2002)、分塊不完全LU分解預(yù)條件(block ILU)(柳建新等，2009)和高斯預(yù)條件(GS)(Adams et al., 2003)技術(shù)，以及散度校正(Dong and Egbert，2019)來改善迭代算法的收斂性.但是隨著網(wǎng)格數(shù)量的增加和頻率的降低，這些方法的收斂性同樣會(huì)變差(Koldan et al., 2014).

多重網(wǎng)格法(MG)(Oosterlee, 1997；魯晶津等，2009；Horesh and Haber, 2011；柳建新等，2011; Pan et al., 2017)將細(xì)網(wǎng)格上的大型方程求解問題轉(zhuǎn)換到較粗網(wǎng)格上的更小的方程組求解問題，有效的消除高頻誤差和低頻誤差，從而達(dá)到快速求解方程的目的.且其計(jì)算量隨著網(wǎng)格數(shù)的增加呈線性增加，因此廣泛應(yīng)用于求解大型橢圓型微分方程問題.它包含代數(shù)多重網(wǎng)格法(AMG)(魯晶津，2010; Koldan et al., 2014；Chen et al., 2017)和幾何多重網(wǎng)格法(GMG)(J?nsth?vel, 2006；Mulder, 2006, 2008; Li et al., 2016a, b; Guo et al., 2022).由于GMG構(gòu)造各算子簡(jiǎn)單、高效，該方法常用于構(gòu)建高效的大型方程組求解器或預(yù)條件技術(shù)(Aruliah and Ascher, 2002；Koldan et al., 2014；Hassan et al., 2019).但是隨著方程組各向異性的增強(qiáng)(由模型網(wǎng)格拉伸的加大等引起)，GMG的收斂效率會(huì)有所降低.Mulder (2006)表明將GMG與基于Krylov子空間的迭代求解方法相結(jié)合可以改善其收斂性.

但是如果采用傳統(tǒng)的平滑算法如GS等，當(dāng)頻率比較低時(shí)，即便采用散度校正，GMG也無法有效的把高頻誤差平滑掉(Li et al., 2020).這時(shí)可以將場(chǎng)的方程轉(zhuǎn)化為矢量勢(shì)和標(biāo)量勢(shì)方程(Bochev et al., 2008)或采用分塊GS平滑算法，即通過一次計(jì)算一個(gè)局部小型方程組(比如一個(gè)點(diǎn)所連接的六個(gè)分量)來改善GMG的收斂性(J?nsth?vel, 2006).但是這種逐點(diǎn)迭代涉及大量的迭代過程，通常計(jì)算比較耗時(shí).針對(duì)該問題，Li等(2020)提出四色分塊GS平滑算法，避免了逐點(diǎn)迭代，顯著提高了分塊GS平滑算法效率.如果該方法能用來作為GMG的平滑算法，將會(huì)大大提高GMG的計(jì)算效率.

本文結(jié)合四色分塊GS為平滑算法GMG和BiCGstab，提出一種高效的MT正演方法.該方法將四色分塊GS為平滑算法GMG作為BiCGstab的預(yù)條件技術(shù)，以最大程度地利用這兩種方法的優(yōu)點(diǎn).為了驗(yàn)證算法的正確性和高效性，本文基于三個(gè)合成模型對(duì)比了本文提出的正演方法和基于BiCGstab的傳統(tǒng)預(yù)條件技術(shù)(GS、block ILU和SSOR)的迭代次數(shù)和收斂時(shí)間.

1 交錯(cuò)網(wǎng)格有限差分法

由于MT所采用的頻率比較低，因此在大地中可以忽略位移電流的作用.假設(shè)時(shí)諧因子為eiω t并且使用國(guó)際標(biāo)準(zhǔn)單位.在頻率域中可以將電磁場(chǎng)的麥克斯韋方程組進(jìn)行簡(jiǎn)化，得到關(guān)于電場(chǎng)的二階偏微分方程：

(1)

式中E表示電場(chǎng)矢量，σ和μ分別表示電導(dǎo)率和磁導(dǎo)率，ω為角頻率.結(jié)合邊界條件，就可以對(duì)(1)式進(jìn)行數(shù)值和/或解析求解，本文采用的邊界條件為狄利克雷第一類邊界條件.

圖1 交錯(cuò)有限差分離散示意圖：電場(chǎng)定義在單元的棱邊， 磁場(chǎng)定義在單元面中心Fig.1 Staggered finite difference discretization with electric fields on cell edges and magnetic fields placed on cell faces

我們采用交錯(cuò)網(wǎng)格FDM對(duì)方程(1)進(jìn)行離散化.如圖1所示，分別將電場(chǎng)和磁場(chǎng)分量定義在單元棱邊和單元面中心.采用矢量形式表示，則方程(1)離散形式可以表示為：

[C*C+diag(iωμσe)]e=0,

(2)

式中：e表示離散電場(chǎng)矢量，diag為對(duì)角矩陣，σe為單元棱邊的電導(dǎo)率(通過體積平均圍著棱邊的四個(gè)單元的電導(dǎo)率得到).C表示離散旋度算子，它將單元的棱邊上的矢量映射到單元的面上,C*是C的伴隨矩陣，表示將單元面上的矢量映射回到單元棱邊上(Egbert and Kelbert, 2012; Kelbert et al., 2014; Cherevatova et al., 2018).

將方程(2)簡(jiǎn)化為如下矩陣形式的線性方程：

Ae=b,

(3)

式中A表示方程的系數(shù)矩陣，b表示右端源項(xiàng).方程(3)可采用直接法或迭代法進(jìn)行求解，一旦求解得到e，磁場(chǎng)矢量可以通過以下關(guān)系得到：

h=(-iωμ)-1Ce.

(4)

對(duì)于MT正演，其計(jì)算區(qū)域往往比較大，方程(3)中系數(shù)矩陣通常是大型稀疏的，采用直接求解法求解該方程需要耗費(fèi)大量的計(jì)算時(shí)間和內(nèi)存.因此，針對(duì)這類問題，往往采用迭代法進(jìn)行求解(如本文中使用的BiCGstab).但是由于雙旋度算子存在豐富的零空間，迭代法在求解線性方程(3)的時(shí)候收斂速度較慢，為了提高計(jì)算效率，通常需要采用預(yù)條件技術(shù).同時(shí)隨著ω的減小，數(shù)值解不能充分地模擬電導(dǎo)率不連續(xù)處的電荷積累，這也會(huì)造成迭代法的收斂速度變慢，這時(shí)需要通過強(qiáng)加散度校正來提高收斂速度，其具體形式如下：

(5)

2 預(yù)條件技術(shù)

通常情況下，構(gòu)建一個(gè)矩陣P，稱之為預(yù)條件矩陣，將對(duì)于方程(3)的求解轉(zhuǎn)化為求解以下方程：

P-1Ae=P-1b.

(6)

通過這種轉(zhuǎn)換，可以有效降低系數(shù)矩陣A的條件數(shù)，相比于方程(3)，方程(6)更容易求解.在理想情況下預(yù)條件矩陣P取為A，式(6)直接可以得到電場(chǎng)的解.但是這時(shí)式(6)的預(yù)條件求解等效于式(3)的求解，失去了使用預(yù)條件技術(shù)的意義.實(shí)際應(yīng)用中，我們通常采用更容易求解的P來近似A，比如常用的有SSOR、block ILU和 GS等.

多重網(wǎng)格法(MG)將最細(xì)網(wǎng)格上的殘差逐步投影到更粗網(wǎng)格上，在最粗網(wǎng)格上花更少的代價(jià)進(jìn)行方程求解，然后將所求結(jié)果逐步插值到最細(xì)網(wǎng)格上對(duì)解進(jìn)行校正，從而達(dá)到快速求解的目的.它主要包含有以下三個(gè)過程，首先在較細(xì)網(wǎng)格上進(jìn)行少數(shù)幾次光滑迭代，消除高頻誤差；然后將低頻誤差投影到較粗網(wǎng)格上，細(xì)網(wǎng)格上的部分低頻誤差在粗網(wǎng)格上成為高頻成分，通過數(shù)次迭代可以有效剔除；重復(fù)以上過程直到最粗網(wǎng)格，這時(shí)最粗網(wǎng)格上未知數(shù)通常比較少，可以采用直接解法進(jìn)行快速精確求解；最后從最粗網(wǎng)格將解逐級(jí)插值到細(xì)網(wǎng)格上以校正細(xì)網(wǎng)格的解，在每一次插值都需要進(jìn)行數(shù)次平滑迭代以消除插值帶來的高頻誤差，直到最細(xì)網(wǎng)格.

圖2 多重網(wǎng)格法迭代過程示意圖Fig.2 Diagram of iterative process of multigrid method

多重網(wǎng)格有多種循環(huán)方式(例如：V-循環(huán)，W-循環(huán)和F-循環(huán))，根據(jù)求解問題的復(fù)雜程度可以選擇不同的循環(huán)方式，其中V-循環(huán)是最簡(jiǎn)單的循環(huán)方式，后二者都是由V-循環(huán)嵌套而成.最為簡(jiǎn)單的V-循環(huán)二重網(wǎng)格法有以下5個(gè)步驟(h和2h分別代表細(xì)網(wǎng)格和粗網(wǎng)格的步長(zhǎng))：

(5)后光滑：在插值過程中往往會(huì)引入新的誤差，因此還需要對(duì)方程Aheh=bh進(jìn)行n2次光滑迭代以消除新的誤差.

(7)

(8)

圖3 基于棱邊體積的加權(quán)限制算子的二維示意圖 (J?nsth?vel，2006)Fig.3 A two-dimensional schematic of the weighting restriction operator based on edge volumes (J?nsth?vel, 2006)

粗網(wǎng)格操作算子A2h產(chǎn)生采用文獻(xiàn)中的模塊思想(Egbert and Kelbert 2012; Kelbert et al., 2014; Cherevatova et al., 2018)，很容易產(chǎn)生粗網(wǎng)格上的方程(3)中的系數(shù)矩陣.平滑算法采用四色分塊GS算法(Li et al., 2020)，該平滑算法將所有節(jié)點(diǎn)分成四種顏色且相同顏色的每一節(jié)點(diǎn)相連接的分量組成的系統(tǒng)完全解耦(具有高度并行性)，可以同時(shí)進(jìn)行更新.四色分塊GS算法在顏色間進(jìn)行循環(huán)，當(dāng)一種顏色更新完成后，相應(yīng)的解可以作為下一種顏色更新的邊界條件.一旦獲得限制算子、插值算子、粗網(wǎng)格操作算子及平滑算法，我們就可以按照前面的流程進(jìn)行多重網(wǎng)格計(jì)算.

3 數(shù)值計(jì)算

在這一節(jié)中，我們測(cè)試了GMG 預(yù)條件技術(shù)的正確性和效率.首先，設(shè)計(jì)了一個(gè)層狀模型，通過比較數(shù)值解和解析解來驗(yàn)證該算法的正確性.然后，用一個(gè)經(jīng)典的雙阻模型和一個(gè)國(guó)際標(biāo)準(zhǔn)模型從迭代次數(shù)和計(jì)算時(shí)間兩方面測(cè)試了算法的效率.我們的算法與傳統(tǒng)的預(yù)條件技術(shù)的BiCGstab求解器進(jìn)行比較，包括GS、SSOR和block ILU.為了加速收斂，對(duì)于傳統(tǒng)預(yù)條件技術(shù)的 BiCGstab求解器，我們采用每40次迭代使用一次散度校正，而GMG預(yù)條件技術(shù)不需要散度校正.所有測(cè)試在AMD (R) 7-3700x 3.59 Hz的計(jì)算機(jī)上進(jìn)行，當(dāng)?shù)较鄬?duì)殘差小于等于10-10或達(dá)到1000次，迭代終止.

3.1 算法驗(yàn)證

首先設(shè)計(jì)了一個(gè)一維層狀模型，該模型在電導(dǎo)率為200 Ωm的均勻半空間上方有一個(gè)電導(dǎo)率為100 Ωm的薄層，層厚2590 m.采用非均勻網(wǎng)格對(duì)64×64×128個(gè)單元進(jìn)行離散，最小的網(wǎng)格單元為30 m×30 m×10 m，研究區(qū)域?yàn)?115 km×7115 km×2372 km.對(duì)于周期選取，我們?cè)?.1～1000 s的范圍內(nèi)采用對(duì)數(shù)間隔選取20個(gè)計(jì)算周期.分別使用解析法和數(shù)值法(求解器為GMG-BiCGstab)對(duì)模型進(jìn)行計(jì)算得到視電阻率和相位曲線(圖4)，并求出相應(yīng)的相對(duì)誤差.結(jié)果顯示：視電阻率和相位的最大相對(duì)誤差分別小于0.6%和0.4%(圖4c)，表明兩種方法得出的結(jié)果具有很好的一致性，驗(yàn)證了本文提出的算法的正確性.

3.2 雙塊狀電阻率模型

為了進(jìn)一步說明算法的準(zhǔn)確性和高效性，在本小節(jié)我們?cè)O(shè)計(jì)了一個(gè)雙塊狀電阻率模型.如圖5所示：異常體的電阻率分別為10Ωm和1000 Ωm，尺寸大小均為20 km×20 km×20 km，埋深均為10 km，模型的背景電阻率為100 Ωm，其中最小單元的大小為2 km×2 km×2 km，網(wǎng)格數(shù)量為128×128×128(不包括空氣層)，正演模擬區(qū)域?yàn)?81 km×281 km×311 km.我們選擇1、10、100、1000 s四個(gè)計(jì)算周期進(jìn)行正演.

為了更好地驗(yàn)證GMG預(yù)條件算法的正確性，在周期為100 s，我們對(duì)比了GMG預(yù)條件技術(shù)和傳統(tǒng)(例如：SSOR)預(yù)條件技術(shù)在兩種極化模式下的視電阻率和相位正演結(jié)果.如圖6所示：視電阻率的最大差異小于2×10-6，相位的最大差異小于1.5×10-6，說明兩者結(jié)果吻合得很好，也進(jìn)一步驗(yàn)證本文所提出的GMG預(yù)條件算法的正確性.

圖4 層狀模型GMG-BiCGstab與解析解的(a)視電阻率 和(b)相位對(duì)比及(c)相對(duì)誤差Fig.4 Comparison of the GMG-BiCGstab solution with analytic solution based on the layered model in terms of (a) apparent resistivity and (b) phase. (c) are the relative errors

圖7所示的是在不同周期下基于GMG預(yù)條件技術(shù)的BiCGstab在有、無散度校正(分別寫為GMG-DC和GMG)時(shí)的迭代次數(shù)和計(jì)算時(shí)間對(duì)比(每進(jìn)行2次迭代計(jì)算執(zhí)行1次散度校正).結(jié)果顯示有無散度校正對(duì)迭代次數(shù)無任何影響；散度校正需要消耗額外的計(jì)算時(shí)間，加散度校正反而會(huì)使計(jì)算時(shí)間顯著增加.由于GMG預(yù)條件技術(shù)采用了四色分塊GS平滑算法，它滿足局部電流散度為零，因此無需額外強(qiáng)加散度校正算法.

圖5 雙塊狀電阻率模型示意圖Fig.5 Diagram of two-block resistivity model

圖6 雙塊狀電阻率模型在周期為100 s時(shí)，不同求解器(GMG-BiCGstab與SSOR-BiCGstab) 在中心測(cè)線的視電阻率和相位的相對(duì)差異Fig.6 For the two-block resistivity model, the relative differences of apparent resistivity and phase along a middle survey line using different solvers (GMG-BiCGstab and SSOR-BiCGstab) at 100 s

圖7 GMG-BiCGstab法有、無散度校正的數(shù)值表現(xiàn) (a)(c) x-極化模式； (b)(d) y-極化模式； (a)(b) 迭代次數(shù)； (c)(d) 計(jì)算時(shí)間.Fig.7 Numerical comparison of the GMG-BiCGstab method with and without divergence correction (a)(c) x-polarization mode; (b)(d) y-polarization mode; (a)(b) Iteration numbers; (c)(d) Compute time.

接下來將對(duì)比GMG、GS、block ILU、SSOR四種不同的預(yù)條件技術(shù)的數(shù)值表現(xiàn)，它們的迭代過程和計(jì)算時(shí)間對(duì)比結(jié)果分別如圖8—9所示.圖8結(jié)果顯示，對(duì)所有周期GMG預(yù)條件技術(shù)的迭代次數(shù)比其他預(yù)條件技術(shù)的迭代次數(shù)大大減少，GMG預(yù)條件技術(shù)基本在7次以內(nèi)就能收斂.在傳統(tǒng)預(yù)條件技術(shù)中，SSOR效果最好.GMG預(yù)條件技術(shù)迭代次數(shù)基本和相對(duì)殘差的對(duì)數(shù)成線性關(guān)系，而傳統(tǒng)預(yù)條件技術(shù)隨著相對(duì)殘差的減少收斂會(huì)明顯變差.計(jì)算時(shí)間對(duì)比如圖9所示，對(duì)所有的周期，GMG預(yù)條件技術(shù)的計(jì)算時(shí)間遠(yuǎn)遠(yuǎn)少于其他預(yù)條件技術(shù).相比其他方法，在x-極化模式下GMG預(yù)條件技術(shù)的計(jì)算時(shí)間減少了73%～90%，在y-極化模式下的計(jì)算時(shí)間減少了76%～88%.

3.3 DTM1模型

本小節(jié)我們將GMG-BiCGstab法應(yīng)用于求解一個(gè)更復(fù)雜的Dublin模型1(DTM1)(Miensopust et al., 2013).如圖10所示，在均勻半空間(100 Ωm)中，存在三個(gè)埋深、尺寸、電阻率均不同的異常體，異常體的參數(shù)如表1所示.我們將模型均勻剖分，單元的尺寸均為2.5 km×2.5 km×2.5 km，網(wǎng)格數(shù)量為128×128×128(不包括空氣層)，總的區(qū)域?yàn)?20 km×320 km×320 km.與3.2節(jié)一樣，我們采用的周期為1 s、10 s、100 s和1000 s.

圖8 雙塊狀電阻率模型使用不同預(yù)條件技術(shù)的BiCGStab在不同周期的收斂過程 (a)(b)(c)(d) x-極化模式； (e)(f)(g)(h) y-極化模式.Fig.8 The convergence process of BiCGStab using different preconditioners for different periods based on the two-block resistivity model (a)(b)(c)(d) x-polarization mode; (e)(f)(g)(h) y-polarization mode.

圖9 雙塊狀電阻率模型使用不同預(yù)條件技術(shù)的 BiCGStab在不同周期段內(nèi)的計(jì)算時(shí)間對(duì)比圖 (a) x-極化模式； (b) y-極化模式.Fig.9 Comparison of computer time of BiCGStab using different preconditioners for different periods based on the two-block resistivity model (a) x-polarization mode; (b) y-polarization mode.

表1 DTM1 模型中異常體在三個(gè)方向的大小Table 1 The anomalies size in three directions for DTM1

針對(duì)該模型，在周期為100 s，我們對(duì)比了GMG-BiCGstab的正演結(jié)果和MTNet(http:∥www.mtnet.info/workshops/mt3dinv/2008_Dublin/dublin_forward.html)提供的正演結(jié)果，MTNet結(jié)果由不同的算法(Mackie et al., 1994; Siripunvaraporn et al., 2002; Farquharson et al., 2002; Nam et al., 2007)得到.結(jié)果顯示，我們的結(jié)果與MTNet的結(jié)果基本吻合，進(jìn)一步驗(yàn)證本文提出的算法的正確性.

GMG預(yù)條件技術(shù)和三種傳統(tǒng)預(yù)條件技術(shù)(GS、block ILU和SSOR)的迭代收斂過程如圖12所示，結(jié)果顯示收斂過程對(duì)比與前一模型的對(duì)比結(jié)果相似.對(duì)所有周期，GMG預(yù)條件技術(shù)迭代次數(shù)與相對(duì)殘差的對(duì)數(shù)成線性關(guān)系，其他預(yù)條件技術(shù)隨著相對(duì)殘差的減小收斂明顯變差.GMG預(yù)條件技術(shù)收斂需要的迭代次數(shù)遠(yuǎn)遠(yuǎn)小于其他三種傳統(tǒng)的預(yù)條件技術(shù).不同預(yù)條件技術(shù)收斂需要的計(jì)算時(shí)間如圖13所示，對(duì)于所有周期，GMG預(yù)條件技術(shù)的耗時(shí)遠(yuǎn)遠(yuǎn)小于其他預(yù)條件技術(shù).其中，在周期為1000 s時(shí)，GMG預(yù)條件技術(shù)在x-極化模式下的計(jì)算時(shí)間相對(duì)GS減小84%(188 s vs.1245 s)，相對(duì)block ILU減小83%(188 s vs.1131 s)，相對(duì)SSOR減小79%(188 s vs.926 s).y-極化模式下計(jì)算時(shí)間相對(duì)GS減小84%(175 s vs.1119 s)，相對(duì)block ILU減小83%(175 s vs.1089 s)，相對(duì)SSOR減小78%(175 s vs.825 s).

圖10 DTM1模型示意圖 綠色、藍(lán)色和紅色分別表示電阻率為10 Ωm、1 Ωm和10000 Ωm的異常體.背景電阻率為100 Ωm.Fig.10 2D view of the DTM1 model Green,blue and red represent 10 Ωm, 1 Ωm and 10000 Ωm, respectively with background conductivity of 100 Ωm.

圖11 DTM1模型在周期為100 s時(shí)，中心測(cè)線的響應(yīng)對(duì)比 (a)(c) x-極化模式； (b)(d) y-極化模式； (a)(b) 視電阻率； (c)(d) 相位.Fig.11 Based on the DTM1 model, the response comparison along a middle survey line at period of 100 s (a)(c) x-polarization mode; (b) (d) y-polarization mode; (a)(b) Apparent resistivity; (c) (d) Phase.

圖12 DTM1模型使用不同預(yù)條件技術(shù)的BiCGStab在不同周期的收斂過程 (a)(b)(c)(d) x-極化模式； (e)(f)(g)(h) y-極化模式.Fig.12 The convergence process of BiCGStab using different preconditioners for different periods based on the DTM1 model (a)(b)(c)(d) x-polarization mode; (e) (f)(g)(h) y-polarization mode.

圖13 DTM1模型使用不同預(yù)條件技術(shù)的BiCGStab 在不同周期的計(jì)算時(shí)間對(duì)比圖 (a) x-極化模式； (b) y-極化模式.Fig.13 Comparison of computer time of BiCGStab using different preconditioners for different periods based on the DTM1 model (a) x-polarization mode; (b) y-polarization mode.

4 結(jié)論

3D MT 在礦產(chǎn)資源勘探和科學(xué)研究等領(lǐng)域中發(fā)揮著越來越重要的作用，但是對(duì)于覆蓋范圍大的MT測(cè)量數(shù)據(jù)，進(jìn)行全區(qū)的3D反演仍然非常具有挑戰(zhàn)性，需要提高大規(guī)模3D MT正演模擬的計(jì)算效率.為此，本文提出了一種高效的GMG預(yù)條件技術(shù)用于3D MT FDM正演模擬.GMG將細(xì)網(wǎng)格上的大型稀疏方程求解通過若干次簡(jiǎn)單平滑迭代逐次轉(zhuǎn)化到更粗網(wǎng)格上以更小的代價(jià)進(jìn)行高效求解，特別適合大型稀疏矩陣的求解.該預(yù)條件技術(shù)的平滑算法采用四色分塊GS法，滿足局部電流散度為零，避免了散度校正的使用.通過設(shè)置一個(gè)簡(jiǎn)單層狀模型驗(yàn)證了本算法的正確性，然后設(shè)置一個(gè)雙塊體電阻率模型和一個(gè)國(guó)際標(biāo)準(zhǔn)DTM1模型，基于BiCGStab，通過與傳統(tǒng)的預(yù)條件技術(shù)(GS、block ILU、SSOR)進(jìn)行對(duì)比測(cè)試了本文所提出算法的數(shù)值表現(xiàn).

結(jié)果顯示GMG預(yù)條件技術(shù)在所有的測(cè)試實(shí)例中具有比較穩(wěn)定的數(shù)值表現(xiàn)，相對(duì)殘差的對(duì)數(shù)隨著迭代次數(shù)的增加呈線性衰減，而傳統(tǒng)預(yù)條件技術(shù)隨著相對(duì)殘差的對(duì)數(shù)減小收斂明顯變慢.傳統(tǒng)預(yù)條件技術(shù)的計(jì)算效率往往隨著周期變長(zhǎng)明顯變差，而GMG預(yù)條件技術(shù)的收斂速度不隨周期的變長(zhǎng)而變慢.對(duì)本文所有實(shí)例，GMG預(yù)條件技術(shù)所耗費(fèi)的計(jì)算時(shí)間明顯小于傳統(tǒng)預(yù)條件，比傳統(tǒng)預(yù)條件技術(shù)減少最少73%以上.以上結(jié)果顯示本文提出的GMG預(yù)條件技術(shù)具有很好的數(shù)值表現(xiàn)，適合應(yīng)用于大型3D MT正演計(jì)算.

致謝作者感謝美國(guó)俄勒岡州立大學(xué)Gary D. Egbrt教授對(duì)我們的論文工作提供了很多有意義的討論和建議，同時(shí)也向三位匿名審稿專家以及本文的編輯提出的寶貴修改意見表示衷心的感謝.