韓波， 胡祥云， 黃一凡， 彭榮華， 李建慧， 蔡建超

中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 430074

?

基于并行化直接解法的頻率域可控源電磁三維正演

韓波， 胡祥云*， 黃一凡， 彭榮華， 李建慧， 蔡建超

中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 430074

電磁法的三維數(shù)值模擬是一個(gè)對(duì)數(shù)值算法和計(jì)算機(jī)硬件要求都非常高的問題.對(duì)常用的微分類方法如有限單元法和有限差分法而言，求解最后所得的大型線性方程組是至關(guān)重要的一步，直接影響到正演算法的實(shí)用性.如何高效、穩(wěn)定且準(zhǔn)確地解線性方程長(zhǎng)期以來(lái)一直是被探討的問題.本文實(shí)現(xiàn)了基于線性系統(tǒng)直接求解技術(shù)的頻率域可控源電磁(CSEM)三維正演.使用交錯(cuò)網(wǎng)格有限體積法(FV)來(lái)離散化關(guān)于二次電場(chǎng)的Helmholtz方程；使用直接解法取代傳統(tǒng)的迭代解法來(lái)求解離散線性系統(tǒng)，即對(duì)系統(tǒng)矩陣進(jìn)行完全LU分解，具體通過調(diào)用大規(guī)模并行矩陣直接求解器(MUMPS)來(lái)實(shí)現(xiàn).基于理論模型做了一系列數(shù)值實(shí)驗(yàn)，首先證明了直接解法的高精度和穩(wěn)定性，并考察了其內(nèi)存需求、計(jì)算時(shí)間和并行可伸縮性等主要計(jì)算性能， 最后檢驗(yàn)了所開發(fā)的算法快速模擬多場(chǎng)源CSEM問題的能力以及對(duì)常規(guī)海洋和陸地CSEM模擬的有效性.

可控源電磁法； 三維模擬； LU分解； 直接解法

1 引言

頻率域可控源電磁法(controlled-source electromagnetic methods, CSEM)通過觀測(cè)人工場(chǎng)源在地層、空氣和海水等介質(zhì)中激發(fā)的電磁場(chǎng)來(lái)獲取介質(zhì)的電性分布，長(zhǎng)期以來(lái)一直是金屬礦勘探中最重要的地球物理方法(底青云和王若，2008；Hu et al.，2013)，并在最近十多年逐漸成為海洋油氣勘探中一種不可或缺的地球物理技術(shù)(Constable，2010).為了通過CSEM觀測(cè)數(shù)據(jù)盡可能真實(shí)地還原出介質(zhì)的電性結(jié)構(gòu)以提升勘探效果，實(shí)施三維勘探并對(duì)數(shù)據(jù)進(jìn)行三維模擬是必需的.正演模擬是反演模擬的基礎(chǔ).CSEM中最常用的三維正演方法包括積分方程法(IE)(Avdeev and Knizhnik，2009；王若等，2009；陳桂波等，2009；Zaslavsky et al.，2011)，有限單元法(FE)(徐志鋒和吳小平，2010；Schwarzbach et al.，2011；Li et al.，2011；da Silva et al.，2012；Puzyrev et al.，2013)，有限差分法(FD)(Newman and Alumbaugh，1995；沈金松，2003；Streich，2009；鄧居智，2011)以及與FD密切相關(guān)的有限體積法(FV)(Weiss and Constable，2006；楊波等，2012；Jahandari and Farquharson，2014；韓波等，2015).對(duì)這幾種方法的評(píng)述可以在很多文獻(xiàn)中找到(Avdeev，2005；吳桂桔等，2010； B?rner，2010).

使用FE、FD和FV法對(duì)Maxwell方程離散化后一般會(huì)得到大型、稀疏線性系統(tǒng)方程，在三維情況下待求解的未知數(shù)很容易達(dá)到上百萬(wàn)個(gè)，要高效地求解這樣的方程并不是一個(gè)容易的問題.普遍采取的求解手段是利用Krylov子空間迭代技術(shù)(Saad，2003)來(lái)逐步逼近方程的真實(shí)解.這一類方法不必顯式存儲(chǔ)系數(shù)矩陣的元素，而只需存儲(chǔ)矩陣與向量的乘積(Weiss，2001)，因此對(duì)計(jì)算機(jī)內(nèi)存的需求很小，完全可以在現(xiàn)代的PC機(jī)上執(zhí)行，并且計(jì)算速度很快(在收斂性好的情況下).然而，迭代解法的收斂性能很大程度上取決于系數(shù)矩陣譜的性質(zhì).如果譜的性質(zhì)不佳，就需要大量的迭代來(lái)收斂到滿足精度要求的解，但也可能會(huì)無(wú)法收斂甚至發(fā)散，求解失敗.此時(shí)需要對(duì)原始方程進(jìn)行"預(yù)優(yōu)"處理來(lái)改善系數(shù)矩陣的譜的性質(zhì)，再應(yīng)用迭代解法才有可能快速收斂.在電磁法數(shù)值模擬中，空氣、地層和海水之間巨大的電導(dǎo)率差異，以及不同區(qū)域的離散網(wǎng)格尺寸的巨大差異等往往導(dǎo)致最后的線性系統(tǒng)方程譜性質(zhì)極差(比如具有很大的譜條件數(shù))，必須要預(yù)優(yōu).預(yù)優(yōu)因子的選取也是一個(gè)問題，好的預(yù)優(yōu)因子往往計(jì)算起來(lái)較困難，并且沒有一種預(yù)優(yōu)因子能保證迭代解法對(duì)于任意復(fù)雜的模型都能穩(wěn)定收斂.

基于矩陣分解的直接解法是求解線性方程的另一個(gè)途徑.與迭代解法相比，直接解法分解矩陣要占用大量的內(nèi)存，且計(jì)算時(shí)間一般也更長(zhǎng).尤其是對(duì)于系統(tǒng)矩陣極其龐大的3D正演，直接解法的這兩個(gè)缺點(diǎn)變得非常突出，令普通PC機(jī)難以承受，因此在以往基本不被考慮.然而直接解法相對(duì)于迭代解法也有顯著的優(yōu)點(diǎn)，首先是求解精度明顯更高，且計(jì)算時(shí)間主要取決于系數(shù)矩陣的大小，與其條件數(shù)關(guān)系不大；此外，求解多個(gè)系數(shù)矩陣相同而右端項(xiàng)不同的方程組，只需要進(jìn)行一次矩陣分解，因此計(jì)算量與解一個(gè)方程組相當(dāng)(Grayver，2013)，這一點(diǎn)不僅對(duì)于多場(chǎng)源的CSEM模擬非常有利(Oldenburg et al.，2008)，而且有利于電磁反演中靈敏度的快速構(gòu)建(韓波等，2012).近年來(lái)，直接解法在2D/2.5D電磁模擬中已有應(yīng)用(Abubakar et al.，2008；Streich et al.，2011；Luo et al.，2014).隨著數(shù)值算法的發(fā)展，在科學(xué)計(jì)算領(lǐng)域出現(xiàn)了許多性能優(yōu)良且免費(fèi)的矩陣直接求解器軟件包，其中包括針對(duì)大規(guī)模稀疏矩陣作了高度優(yōu)化的并行求解器，如MUMPS(Amestoy et al.，2006)和PARDISO(Schenk and G?rtner，2004)，加上計(jì)算機(jī)性能的快速提升，使得在小型服務(wù)器、工作站和PC機(jī)群上利用直接解法求解3D 電磁問題成為可能.Streich(2009)、da Silva等(2012)以及Jahandari和Farquharson(2014)分別使用FD、FE和FV法進(jìn)行CSEM三維正演計(jì)算，他們均利用MUMPS作為線性求解器；Schwarzbach等(2011)則在其高階FE算法中調(diào)用PARDISO解線性方程.

本文使用交錯(cuò)網(wǎng)格FV法來(lái)進(jìn)行頻率域CSEM的三維正演計(jì)算，嘗試?yán)么笠?guī)模并行矩陣求解器MUMPS解線性系統(tǒng)方程.通過大量的理論模型正演測(cè)試來(lái)考察MUMPS的求解精度、穩(wěn)定性、內(nèi)存需求、計(jì)算時(shí)間以及并行可伸縮性等計(jì)算性能，并驗(yàn)證所開發(fā)的算法和程序的有效性和廣泛適用性.

2 方法理論

2.1 控制方程

對(duì)于一般CSEM勘探所使用的頻率范圍和大部分已知的傳播介質(zhì)的電阻率變化范圍，位移電流的影響可以忽略，則在電流型場(chǎng)源激發(fā)下關(guān)于電場(chǎng)的二階頻率域微分方程為(取時(shí)諧因子為eiω t)

(1)

其中，E為電場(chǎng)強(qiáng)度(V/m)，σ為電導(dǎo)率(S/m)，ω為角頻率，μ0為真空中的磁導(dǎo)率(H/m)，J為外加場(chǎng)源的電流密度(A·m-2).J為沖激函數(shù)，在場(chǎng)源點(diǎn)具有奇異性(無(wú)窮大)(WeissandConstable，2006)，使得數(shù)值求解方程(1)有一定的困難.為此，可假設(shè)有一個(gè)參考模型σp，在場(chǎng)源J的激發(fā)下同樣滿足上述形式的方程：

(2)

式(1)與式(2)相減，可得

(3)

若將實(shí)際模型的場(chǎng)看做總場(chǎng)，參考模型的場(chǎng)看做一次場(chǎng)，令二次場(chǎng)Es=E-Ep，式(3)可重寫為

(4)

其中Js=(σ-σp)Ep，可看做激發(fā)二次場(chǎng)的電流密度.方程(4)便是關(guān)于二次電場(chǎng)的二階控制方程，也是本文使用FV法直接求解的方程.與總場(chǎng)方程(1)相比，由于Js不是奇異函數(shù)，方程(4)求解起來(lái)更加容易.對(duì)于一次場(chǎng)，一般將參考模型選為所需計(jì)算模型的一維“背景模型”，可以通過解析法求得；一次場(chǎng)與二次場(chǎng)相加便得到總場(chǎng).

2.2 有限體積法

圖1 (a)電磁場(chǎng)的交錯(cuò)采樣和(b)采樣點(diǎn)處的控制體積Fig.1 (a) Sampling positions of EM fields on a staggered grid and (b) the control volume for a sampling site

求解形如式(4)的Helmholtz方程，一種簡(jiǎn)單而有效的方法是交錯(cuò)網(wǎng)格有限差分法(Staggered Finite-Difference, SFD)，其將研究區(qū)域離散化為一系列規(guī)則長(zhǎng)方體單元，然后在長(zhǎng)方體單元上交錯(cuò)地對(duì)電場(chǎng)和磁場(chǎng)采樣，如圖1a所示.該方法在電磁法3D正演中被廣泛使用，其基本理論比較簡(jiǎn)單，可以在很多文獻(xiàn)中找到(e.g.，Newman and Alumbaugh，1995)，這里不再贅述.使用SFD法對(duì)(4)進(jìn)行離散化直接得到的線性方程系數(shù)矩陣是不對(duì)稱的，為了使其對(duì)稱以便于方程的求解，每個(gè)采樣點(diǎn)的線性方程兩端需要同時(shí)乘以一個(gè)對(duì)應(yīng)于該采樣點(diǎn)的“尺度因子”V，V具有某種體積的意義，即SFD線性方程實(shí)際上對(duì)應(yīng)的微分方程為

(5)

對(duì)SFD法稍加改動(dòng)，便可以實(shí)現(xiàn)交錯(cuò)網(wǎng)格有限體積法(SFV)，即在每個(gè)采樣點(diǎn)處定義一個(gè)包圍該采樣點(diǎn)的控制體積，由共享該采樣點(diǎn)所在棱邊的相鄰4個(gè)網(wǎng)格單元體積的1/4組成.以x分量為例，如圖1b所示.對(duì)式(4)在某個(gè)體積單元Ω內(nèi)求積分：

(6)

然后對(duì)式(6)求解.事實(shí)上，式(5)中的尺度因子剛好是SFV中的控制體積值，因此式(5)與式(6)在物理意義上的區(qū)別在于前者用某個(gè)離散點(diǎn)處的變量值代表包圍該點(diǎn)的一片區(qū)域內(nèi)每一處的變量值，后者則考慮到變量在區(qū)域內(nèi)的變化，對(duì)區(qū)域進(jìn)行體積分.式(6)左端第一項(xiàng)可以完全展開(楊波等，2012)，并且與(5)左端第一項(xiàng)具有完全相同的離散形式.對(duì)于(6)的左端第二項(xiàng)，通常認(rèn)為二次電場(chǎng)在控制體積范圍內(nèi)的變化不大，因此Es可提取到積分號(hào)以外；采樣點(diǎn)處的電導(dǎo)率可與SFD中一樣，使用采樣點(diǎn)所在棱邊的相鄰4個(gè)單元電導(dǎo)率的體積加權(quán)平均值，也不必求積分.經(jīng)過這樣的簡(jiǎn)化，SFV的左端項(xiàng)就與SFD完全相同，因此二者的線性方程系數(shù)矩陣相同，唯一不同的是右端項(xiàng).由于CSEM的電磁場(chǎng)可能在很小的空間范圍內(nèi)發(fā)生較大變化(如場(chǎng)源附近)，因此SFV比SFD更能準(zhǔn)確地代表電磁場(chǎng)的分布.

CSEM的場(chǎng)源是局部的，因此可讓剖分域的邊界離場(chǎng)源足夠遠(yuǎn)，從而運(yùn)用齊次的Dirichlet邊界條件，即令邊界上的切向電場(chǎng)為0；SFV法右端項(xiàng)的一次場(chǎng)體積積分可使用數(shù)值積分實(shí)現(xiàn)，例如Gauss-Legendre積分(DavisandRabinowitz，2007)；一次場(chǎng)的計(jì)算為一維正演問題，使用基于Hankel變換的準(zhǔn)解析法(Key，2009).

2.3 大型線性系統(tǒng)方程的求解

Helmholtz方程經(jīng)SFD或SFV離散化后得到線性系統(tǒng)方程：

Ae=b,

(7)

其中A為大型、對(duì)稱、稀疏復(fù)矩陣(僅對(duì)角元為復(fù)數(shù))，其行(列)數(shù)也即待求的未知數(shù)個(gè)數(shù)約為3Nx·Ny·Nz(Nx、Ny和Nz分別為三個(gè)方向的網(wǎng)格數(shù))，每行最多只有13個(gè)非零元.作為一個(gè)簡(jiǎn)單示例，圖2中顯示了一個(gè)網(wǎng)格數(shù)為5×6×7的模型的SFD/SFV離散線性方程系數(shù)矩陣的稀疏結(jié)構(gòu).

本文使用基于矩陣三角分解的直接解法求解方程(7)，具體通過一個(gè)開源的軟件包——多波前大規(guī)模并行求解器(MultifrontalMassivelyParallelSolver，MUMPS)(Amestoyetal.，2006)來(lái)實(shí)現(xiàn).MUMPS專門針對(duì)大型稀疏矩陣，并且還能利用矩陣的對(duì)稱性，就內(nèi)存需求和運(yùn)行時(shí)間作了高度優(yōu)化，因此非常適合解方程(7)；此外MUMPS還能通過消息傳遞接口(MPI)實(shí)現(xiàn)并行化求解.MUMPS解方程包括3個(gè)階段：①矩陣的分析與預(yù)處理；②矩陣的LU(LDLT)分解；③三角方程的求解.其中矩陣的分解占用絕大部分計(jì)算時(shí)間和內(nèi)存.

圖2 網(wǎng)格數(shù)為5×6×7的模型的SFD/SFV離散線性方程系數(shù)矩陣A的稀疏結(jié)構(gòu)Fig.2 Sparsity pattern of A arising from the SFD/SFV discretization on a staggered-grid with size 5×6×7

對(duì)照控制方程(4)或(6)來(lái)看方程(7)，在模型剖分不變的情況下，改變頻率，則系數(shù)矩陣A與右端項(xiàng)b都會(huì)改變，需要重新分解A并求解；若模型剖分和頻率均固定，只改變發(fā)射源的方位或位置，則一次場(chǎng)會(huì)改變，導(dǎo)致b發(fā)生改變，但A不會(huì)受到影響，因此無(wú)需重新進(jìn)行矩陣分解，而可以利用原有的分解因子，對(duì)新的右端項(xiàng)重新進(jìn)行解三角方程就夠了.這樣就大大節(jié)省了計(jì)算量.在實(shí)際CSEM勘探中(Oldenburgetal.，2008；Constable，2010；Grayver，2013)，為了更全面地獲取地下信息，常需要頻繁改變發(fā)射源的極化方向(方位)和空間位置，此時(shí)利用直接解法可以快速計(jì)算多方位、多空間位置發(fā)射源的電磁響應(yīng).

2.4 計(jì)算流程

根據(jù)前文的描述，本文的SFV算法程序的計(jì)算流程如圖3所示，從中可以清晰地看出MUMPS的各個(gè)階段需要執(zhí)行的次數(shù)與頻率數(shù)、場(chǎng)源數(shù)之間的關(guān)系.

圖3 本文SFV三維正演算法的計(jì)算流程Fig.3 Workflow of the presented SFV modeling scheme

3 直接解法的性能測(cè)試

為了測(cè)試MUMPS求解大型線性系統(tǒng)的性能，基于兩個(gè)較典型的層狀地電模型(圖4)進(jìn)行了一系列正演測(cè)試.圖4a為海洋CSEM中常用的模擬高阻油氣薄層的模型(Constable and Weiss，2006)，圖4b為一個(gè)陸地三層H型模型.本節(jié)的測(cè)試主要與模型參數(shù)和頻率有關(guān)，而與一些具體的收發(fā)設(shè)置如測(cè)點(diǎn)設(shè)置無(wú)關(guān)，因此將不提及這些設(shè)置.本文的所有計(jì)算均在中科曙光“天闊W580I-G10”工作站上完成，該工作站包含兩顆Intel XEON E5-2603型CPU，每顆CPU為4核心4線程，主頻1.8 GHz.為了應(yīng)對(duì)矩陣分解的巨大內(nèi)存需求，配置了128 GB的內(nèi)存.操作系統(tǒng)為Ubuntu 12.04.

3.1 解方程的精度與穩(wěn)定性

從基本原理上看，直接解法未作任何近似，因此其求解結(jié)果在理論上是精確的，但實(shí)際執(zhí)行時(shí)由于計(jì)算機(jī)的字長(zhǎng)有限，對(duì)浮點(diǎn)數(shù)運(yùn)算會(huì)有舍入誤差，導(dǎo)致直接解法也不可避免地存在誤差，但這種誤差與迭代解法中由于對(duì)數(shù)學(xué)問題本身的近似帶來(lái)的誤差有本質(zhì)區(qū)別.

圖4 用來(lái)測(cè)試直接解法性能的(a)海洋與(b)陸地一維層狀模型Fig.4 The (a) marine layered 1D model and (b) land layered 1D model used to test the performance of the direct solver

方程(7)中，系數(shù)矩陣A的數(shù)值特性與頻率、空氣電導(dǎo)率等密切相關(guān)，因此方程的求解也可能會(huì)受這些因素的影響，尤其是迭代解法往往對(duì)這些因素非常敏感.為此，對(duì)以上兩個(gè)1D模型正演時(shí)，大范圍改變頻率和空氣電導(dǎo)率來(lái)測(cè)試直接解法的穩(wěn)定性.令空氣電導(dǎo)率從10-5S/m變化到10-11S/m，頻率在實(shí)際勘探常用的范圍內(nèi)變化，表1和表2給出了直接解法解方程的殘差和矩陣分解耗時(shí).這里殘差的定義式為‖Ae-b‖，相對(duì)殘差為‖Ae-b‖/‖b‖；兩個(gè)模型的網(wǎng)格剖分?jǐn)?shù)分別為68×40×43和64×42×46；使用了8個(gè)線程并行計(jì)算.可以看出，空氣電導(dǎo)率的變化對(duì)相對(duì)殘差的影響非常??；頻率由高到低變化，相對(duì)殘差會(huì)有明顯的增加，但最大的相對(duì)殘差也僅僅在10-11數(shù)量級(jí)，而迭代解法收斂的相對(duì)殘差門檻值一般設(shè)在10-5～10-8之間(Grayver，2013).在矩陣分解耗時(shí)方面，頻率和空氣電導(dǎo)率變化對(duì)海洋1D模型的影響很小，對(duì)陸地模型的影響稍大，但時(shí)間的波動(dòng)相對(duì)于整體耗時(shí)并不大.由此可見，在實(shí)際勘探常用的頻率范圍內(nèi)，對(duì)于包含很大電性差異的模型，直接解法都具有相當(dāng)高的精度，計(jì)算速度也很穩(wěn)定.

表1 對(duì)海洋1D模型，頻率與空氣層電導(dǎo)率在較大范圍內(nèi)變化時(shí)直接解法的精度和矩陣分解時(shí)間Table 1 The accuracy and factorization time of the direct solver for the marine 1D model for a wide range of frequencies & air conductivities

3.2 內(nèi)存需求與計(jì)算時(shí)間

再來(lái)考察直接解法的內(nèi)存需求與計(jì)算時(shí)間對(duì)模型剖分網(wǎng)格數(shù)的敏感程度.對(duì)海洋1D模型，發(fā)射頻率為1 Hz，空氣電導(dǎo)率取10-11S/m(本文從此處開始，空氣電導(dǎo)率值均取該值).設(shè)計(jì)了從小到大21種網(wǎng)格，其中最小的網(wǎng)格數(shù)為7×8×8，最大的網(wǎng)格數(shù)為98×80×57.圖5顯示了直接解法的內(nèi)存需求與計(jì)算時(shí)間(MUMPS的三個(gè)階段之和)隨未知數(shù)個(gè)數(shù)(約為3Nx·Ny·Nz)的變化.在調(diào)用MUMPS時(shí)執(zhí)行了“Out-Of-Core”選項(xiàng)，即利用硬盤來(lái)存儲(chǔ)矩陣的分解因子，這樣能節(jié)約很大一部分內(nèi)存.對(duì)于較大的網(wǎng)格數(shù)，使用了多個(gè)進(jìn)程并行計(jì)算(進(jìn)程的增加會(huì)導(dǎo)致內(nèi)存需求的增加，見下一小節(jié)).可以看到，內(nèi)存需求與計(jì)算時(shí)間均隨模型網(wǎng)格數(shù)的增加而快速增加.當(dāng)網(wǎng)格數(shù)為78×50×57時(shí)(對(duì)應(yīng)的未知數(shù)個(gè)數(shù)約為650000)，內(nèi)存需求在16GB左右，而當(dāng)網(wǎng)格數(shù)達(dá)到98×80×57時(shí)，內(nèi)存需求則達(dá)到了接近120GB.雖然無(wú)法給出網(wǎng)格數(shù)與內(nèi)存需求之間的具體關(guān)系式，但不難看出整體上內(nèi)存需求增加的速度要遠(yuǎn)快于網(wǎng)格數(shù)增加的速度.

表2 對(duì)陸地1D模型，頻率與空氣層電導(dǎo)率在較大范圍內(nèi)變化時(shí)直接解法的精度和矩陣分解時(shí)間Table 2 The accuracy and factorization time of the direct solver for the land 1D model for a wide range of frequencies & air conductivities

圖5 對(duì)海洋1D模型，對(duì)于不同的模型網(wǎng)格剖分?jǐn)?shù)直接解法的(a)內(nèi)存需求與(b)解方程時(shí)間Fig.5 (a) Memory requirement and (b) solving time of the direct solver for the marine 1D model for various grid sizes

3.3 并行可伸縮性

在并行計(jì)算中，可伸縮性(Scalability)是指系統(tǒng)通過改變可用計(jì)算資源和調(diào)度方式來(lái)動(dòng)態(tài)調(diào)整自身計(jì)算性能的能力.為了測(cè)試MUMPS的可伸縮性，依次使用1～8個(gè)進(jìn)程對(duì)海洋1D模型進(jìn)行正演計(jì)算，頻率為1 Hz，使用3種網(wǎng)格剖分，統(tǒng)計(jì)了解方程的時(shí)間、相對(duì)于單進(jìn)程運(yùn)算的加速因子、總內(nèi)存消耗和平均每個(gè)進(jìn)程的內(nèi)存消耗等，如圖6所示.可以看到，計(jì)算時(shí)間隨進(jìn)程數(shù)的增加穩(wěn)定下降，當(dāng)進(jìn)程數(shù)達(dá)到8時(shí)，計(jì)算時(shí)間是單進(jìn)程的1/4左右.但從加速因子來(lái)看，還遠(yuǎn)不能達(dá)到理論上最理想的并行(加速因子等于進(jìn)程數(shù))，這是由矩陣分解的可并行程度有限所決定的.對(duì)于同一種網(wǎng)格，總體內(nèi)存消耗隨進(jìn)程數(shù)的增加呈上升趨勢(shì)(個(gè)別情況除外)，8進(jìn)程運(yùn)行的總內(nèi)存消耗是單進(jìn)程的4倍左右，這是由于進(jìn)程之間需要共享很多變量.因此在試圖增加進(jìn)程數(shù)來(lái)提升計(jì)算速度時(shí)，還要考慮到進(jìn)程增多帶來(lái)的內(nèi)存需求增大.盡管如此，內(nèi)存消耗平均到每個(gè)進(jìn)程上，還是會(huì)隨進(jìn)程數(shù)增加而下降的，這一點(diǎn)對(duì)機(jī)群系統(tǒng)比較有利.

4 三維正演程序的有效性檢驗(yàn)

接下來(lái)通過具體正演結(jié)果的比較與分析來(lái)檢驗(yàn)所開發(fā)的SFV三維正演程序的有效性.

4.1 與解析解的對(duì)比

首先，還是基于第3節(jié)中的兩個(gè)層狀模型，比較SFV數(shù)值解與解析解.所有的解析解均使用開源的Dipole1D程序(Key，2009)獲得.

對(duì)于海洋層狀模型，觀測(cè)系統(tǒng)設(shè)置為：發(fā)射電偶極沿x方向，位于海底上方50 m，取偶極中心在x-y平面內(nèi)的投影為坐標(biāo)原點(diǎn)，發(fā)射頻率為1 Hz；接收器位于海底表面，沿x軸方向布設(shè)(即inline觀測(cè)，參考圖4a).圖7給出了網(wǎng)格數(shù)為86×50×57時(shí)Ex分量和Hy分量的振幅-收發(fā)距(MVO)曲線與相位-收發(fā)距(PVO)曲線的SFV解與解析解對(duì)比，圖8則為在3種不同網(wǎng)格剖分下SFV解相對(duì)于解析解的誤差.可以看出，SFV解的誤差整體上隨收發(fā)距的增加而增大；網(wǎng)格為54×32×27時(shí)誤差非常明顯，Hy分量振幅的最大相對(duì)誤差可達(dá)20%以上；當(dāng)網(wǎng)格加密為68×40×43時(shí)，數(shù)值解的精度得到了極大的提高，振幅最大相對(duì)誤差在8%左右，相位最大相差6°左右；進(jìn)一步細(xì)化網(wǎng)格至86×50×57時(shí)，精度繼續(xù)得以提升，但并不顯著，這也許跟網(wǎng)格加密的方式有一定關(guān)系(不在本文討論范圍以內(nèi)).

圖6 對(duì)海洋層狀模型，直接解法解方程的時(shí)間與內(nèi)存需求隨所使用的進(jìn)程數(shù)的改變Fig.6 Solving time & memory requirement of the direct solver vary as the number of computing processors changes for the marine layered model

圖7 對(duì)海洋層狀模型，網(wǎng)格為86×50×57時(shí)水平(a, b)電場(chǎng)和(c, d)磁場(chǎng)的SFV數(shù)值解(圓圈)與解析解(線條)對(duì)比Fig.7 Comparison of horizontal (a, b) electric fields and (c, d) magnetic fields obtained from the SFV method (circles) and the analytic solution (line) for the marine 1D model with grid size 86×50×57

圖8 對(duì)于海洋層狀模型，在3種不同網(wǎng)格剖分下SFV解相對(duì)于解析解的誤差Fig.8 Relative errors of SFV solutions compared to analytic solutions for the marine 1D model for three different grid sizes

圖9 對(duì)陸地層狀模型，坐標(biāo)為(8000, 0)的測(cè)點(diǎn)的水平(a, b)電場(chǎng)、(c, d)磁場(chǎng)和(e, f)卡尼亞視電阻率的SFV解(圓圈)與解析解(線條)對(duì)比，以及(g, h)視電阻率的SFV解相對(duì)于解析解的誤差Fig.9 Comparison of (a, b) electric fields, (c, d) magnetic fields and (e, f) Cagniard apparent resistivities of the site with coordinate (8000, 0) obtained from the SFV method (circles) and the analytic solution (line) for the land 1D model. (g, h) are the apparent resistivity errors of SFV solutions compared to analytic solutions

對(duì)于陸地層狀模型，采用地面可控源音頻大地電磁法(CSAMT)(底青云和王若，2008)的觀測(cè)方式，觀測(cè)系統(tǒng)設(shè)置為：發(fā)射源為沿y方向的1000 m的接地導(dǎo)線，取其中心為坐標(biāo)原點(diǎn)，發(fā)射0.125～8192 Hz之間按對(duì)數(shù)等間隔分布的17個(gè)頻率；測(cè)線平行于發(fā)射源布設(shè)，中點(diǎn)到場(chǎng)源中心的最大距離為8 km(參考圖4b).網(wǎng)格數(shù)為66×80×46.圖9給出了坐標(biāo)為(8000, 0)的測(cè)點(diǎn)上所有頻率的Ey分量、Hx分量和卡尼亞視電阻率ρyx的SFV解與解析解對(duì)比，以及ρyx的誤差.其中ρyx的計(jì)算式為

(8)

可以看出，數(shù)值解與解析解吻合較好.頻率位于64～256 Hz之間時(shí)誤差稍大，其余頻率的視電阻率相對(duì)誤差均在8%以內(nèi)，阻抗相位誤差均在±2°以內(nèi).

4.2 多場(chǎng)源問題的模擬

考慮海洋CSEM勘探中需要發(fā)射源在多個(gè)位置激發(fā)的情況，對(duì)上述海洋1D模型，令發(fā)射源沿x方向在-2000 m到8000 m之間移動(dòng)，其余觀測(cè)參數(shù)與4.1節(jié)中相同.如2.3節(jié)中描述的，頻率不變時(shí)，發(fā)射源每改變一次激發(fā)的位置(或方位)，就會(huì)給三維正演模擬中線性系統(tǒng)方程增加一個(gè)右端項(xiàng)，但不會(huì)增加系數(shù)矩陣.圖10為右端項(xiàng)數(shù)目由1增加到101時(shí)，SFV正演計(jì)算的總時(shí)間及各主要部分消耗的時(shí)間變化.可以看到，當(dāng)右端項(xiàng)比較少時(shí)，矩陣分解占用了絕大部分計(jì)算時(shí)間，而解三角方程(Solution)這一步消耗的時(shí)間幾乎可以忽略不計(jì)；隨著右端項(xiàng)的增多，由于SFV法每構(gòu)建一次右端項(xiàng)都要在控制體積內(nèi)進(jìn)行數(shù)值積分，構(gòu)建右端項(xiàng)所占用的時(shí)間比例越來(lái)越大；當(dāng)右端項(xiàng)增加到101個(gè)時(shí)，構(gòu)建右端項(xiàng)占用的時(shí)間比例已經(jīng)很顯著了，但解三角方程所占比重仍然很小，總計(jì)算時(shí)間也僅僅只比1個(gè)右端項(xiàng)時(shí)增加了約1倍.

4.3 海洋三維模型

若將前面的海洋層狀模型中高阻層在水平方向的邊界截?cái)啵瑒t得到更接近真實(shí)情況的三維油氣儲(chǔ)層模型，如圖11所示.已有不少其他作者研究過類似的模型(e.g.，Constable and Weiss，2006；Streich，2009)，本文的模型與這些作者略有差別，即高阻體在x-y平面內(nèi)的截面為正方形而非圓形.

分別考慮海水無(wú)窮深、不含空氣層和海水1000 m深、包含空氣層這兩種情況.參數(shù)設(shè)置為：場(chǎng)源沿x方向，位于高阻體x負(fù)向邊界中點(diǎn)的正上方，海底上方50 m，坐標(biāo)為(0, 0, -50 m)；發(fā)射頻率為1 Hz；接收器沿x方向布設(shè)、分布于海底y=0的測(cè)線上(inline觀測(cè)).對(duì)高阻體截面邊長(zhǎng)(設(shè)為a)為1 km、2 km和5 km的三種情況進(jìn)行了正演.所得到的水平電場(chǎng)Ex分量如圖12所示，其中還給出了不含高阻體的半空間(a=0)和高阻體在水平方向無(wú)限延伸(a=∞)這兩種極限情況的響應(yīng)(使用解析法獲得).

看圖12a和12b，Ex響應(yīng)的分段特征非常明顯：在小收發(fā)距內(nèi)，經(jīng)海水傳播的直達(dá)波占主導(dǎo)地位，振幅大約按1/r3衰減；隨著收發(fā)距的增大，地層波開始起主導(dǎo)作用，含高阻層(a>0)與不含高阻層(a=0)的響應(yīng)逐漸區(qū)分開來(lái)，振幅曲線在對(duì)數(shù)坐標(biāo)系中的斜率趨于常數(shù)，說(shuō)明振幅趨于按指數(shù)衰減.a=1 km與半空間模型的響應(yīng)幾乎重合，說(shuō)明該尺寸的異常體難以分辨；a=2 km與半空間的響應(yīng)雖然在大收發(fā)距時(shí)的衰減速度保持一致，但二者的絕對(duì)差值較為明顯，說(shuō)明該尺寸的異常體容易分辨；a=5 km與半空間的響應(yīng)差別非常明顯，且在收發(fā)距大于高阻體邊長(zhǎng)時(shí)與高阻體無(wú)限延伸(a=∞)的響應(yīng)區(qū)分開來(lái)，衰減速度逐漸趨于半空間模型的衰減速度.根據(jù)以上分析，可以再次確認(rèn)Constable和Weiss

圖10 對(duì)于海洋層狀模型，SFV正演總體計(jì)算時(shí)間及各主要階段消耗的時(shí)間隨右端項(xiàng)數(shù)目的變化Fig.10 The total running time and running time of the significant tasks in SFV modeling vs the numbers of right hand sides for the marine 1D model

(2006)的結(jié)論：對(duì)于此類水平薄板狀高阻體，只有當(dāng)其水平尺寸至少為埋深的2倍時(shí)，才可以利用CSEM將其分辨出來(lái).

再來(lái)看包含空氣層的情況.對(duì)比圖12c、d和a、b，可以看到，空氣波的影響非常明顯，在收發(fā)距超過5 km時(shí)，a=1 km、a=2 km和半空間模型的響應(yīng)中空氣波開始起主導(dǎo)作用，電場(chǎng)振幅衰減極其緩慢；而a=5 km的模型直到收發(fā)距達(dá)到8 km左右時(shí)才開受空氣波影響，a=∞的模型在所展示的收發(fā)距內(nèi)還未被空氣波影響.

圖11 海洋三維油氣儲(chǔ)層模型Fig.11 The marine 3D reservoir model

圖12 三維儲(chǔ)層模型的儲(chǔ)層截面邊長(zhǎng)變化時(shí)的Ex響應(yīng)(a, b)為不考慮空氣層、假設(shè)海水無(wú)限深的情況，(c, d) 為海水1000 m深、包含空氣層的情況.Fig.12 Ex responses for the 3D reservoir model with various horizontal edge lengths of the reservoir(a, b) is for the case air layer is not presented and the seawater is arbitrarily thick, while (c, d) is for the case the seawater is 1000 m thick and air layer is included.

4.4 陸地三維模型

最后對(duì)一個(gè)陸地三維模型進(jìn)行正演.如圖13所示，在100 Ωm的均勻半空間中植入一個(gè)400 m×500 m×400 m的低阻塊體，埋深為200 m，電阻率5 Ωm.使用兩個(gè)分別沿y方向和x方向的正交場(chǎng)源交替發(fā)射來(lái)模擬張量CSAMT觀測(cè)，發(fā)射導(dǎo)線長(zhǎng)度均為1000 m，發(fā)射0.125～8192 Hz之間按對(duì)數(shù)等間隔分布的17個(gè)頻率；發(fā)射源的中心為坐標(biāo)原點(diǎn)，到低阻體中心的水平距離為6000 m；網(wǎng)格數(shù)為62×44×53.

由于有兩個(gè)正交場(chǎng)源，因此可以使用類似于大地電磁法(MT)中的方式來(lái)計(jì)算張量阻抗，進(jìn)而計(jì)算卡尼亞視電阻率響應(yīng)，具體計(jì)算式見譚捍東等(2004).圖14和圖15分別為沿x方向和沿y方向穿過低阻體正上方的兩條測(cè)線上的響應(yīng)，圖16則為頻率為4 Hz的響應(yīng)切片，均顯示了XY和YX兩種模式的響應(yīng).整體上來(lái)看，這些對(duì)CSEM的電磁場(chǎng)用MT的方式處理得到的響應(yīng)與MT響應(yīng)有一定的相似度，異常體的空間輪廓在正演響應(yīng)中有良好的反映.XY和YX模式的響應(yīng)特征有所不同，XY響應(yīng)對(duì)x方向的電性邊界反映較好，YX響應(yīng)則對(duì)y方向的電性邊界反映較好；在與低阻體對(duì)應(yīng)的位置、對(duì)應(yīng)的頻率處，XY視電阻率比YX視電阻率在數(shù)值上更接近低阻體的電阻率值.在頻率低于2 Hz時(shí)，視電阻率值均超過了300 Ωm，甚至可達(dá)1000 Ωm，這是因?yàn)轭l率較低時(shí)，電磁波在背景地層中的趨膚深度較大，測(cè)點(diǎn)所處的位置遠(yuǎn)遠(yuǎn)不能滿足遠(yuǎn)區(qū)條件，此時(shí)卡尼亞視電阻率公式完全不適用.

圖13 陸地三維低阻體模型以及發(fā)射源的布置示意圖Fig.13 The land 3D conductive brick model and the source geometry

5 結(jié)論與討論

本文利用多波前大規(guī)模并行求解器軟件包MUMPS，實(shí)現(xiàn)了基于矩陣三角分解的頻率域可控源電磁三維有限體積正演.

針對(duì)MUMPS直接解法的計(jì)算性能做了一系列測(cè)試，結(jié)果表明：(1)直接解法非常穩(wěn)定，不易被系數(shù)矩陣的病態(tài)程度所影響，在實(shí)際CSEM勘探常用的頻率范圍內(nèi)，對(duì)于包含很大電性差異的模型，都能獲得高精度的解，計(jì)算速度也很穩(wěn)定；(2)直接解法的內(nèi)存需求與計(jì)算時(shí)間取決于模型網(wǎng)格數(shù)，且均隨模型網(wǎng)格數(shù)的增加而急劇增加，內(nèi)存需求很容易超出單臺(tái)普通PC機(jī)的物理內(nèi)存；(3)使用多進(jìn)程并行執(zhí)行MUMPS時(shí)，其計(jì)算速度能得到大幅度提升，也能有效地減少單節(jié)點(diǎn)的內(nèi)存消耗，但由于矩陣分解本身的可并行程度有限，MUMPS的并行可伸縮性也比較有限.

對(duì)理論模型正演結(jié)果的比較與分析表明：本文的SFV三維正演算法是有效的，適用于常規(guī)的陸地和海洋CSEM勘探，并且由于利用了直接解法，非常適合于CSEM中典型的多場(chǎng)源問題，場(chǎng)源的增加只會(huì)帶來(lái)極小的額外計(jì)算量.

就目前最普通、最可用的計(jì)算設(shè)備如個(gè)人電腦、工作站和小型服務(wù)器的平均性能而言，電磁法的三維數(shù)值模擬仍然是一個(gè)對(duì)計(jì)算資源需求很高的問題，使用直接解法時(shí)更是如此，因此模型剖分的網(wǎng)格數(shù)要嚴(yán)格控制(本文3.2節(jié)中對(duì)此給出了參考)，這無(wú)疑會(huì)影響計(jì)算結(jié)果的準(zhǔn)確性.本文的研究便受限于計(jì)算設(shè)備的性能，在并行計(jì)算的進(jìn)程數(shù)、地面CSAMT模擬中的收發(fā)距等方面受到了較大的限制.另外，在計(jì)算設(shè)備的性能足夠強(qiáng)勁的情況下，對(duì)多頻率的正演計(jì)算實(shí)施“多層次”并行化無(wú)疑能進(jìn)一步提高效率，即將不同頻率的正演任務(wù)分配到不同的節(jié)點(diǎn)組，然后每個(gè)節(jié)點(diǎn)組對(duì)單個(gè)頻率進(jìn)行并行正演計(jì)算.這些都是有待改進(jìn)和完善的地方.

隨著數(shù)值算法的發(fā)展以及計(jì)算機(jī)性能的快速提升，可以預(yù)見，直接解法由于其特殊的優(yōu)點(diǎn)，必將越來(lái)越多地被應(yīng)用于3D電磁模擬中.

致謝 感謝MUMPS軟件包和Dipole1D的開發(fā)者們，沒有他們的開源精神，本文的研究將無(wú)法完成；特別感謝美國(guó)俄勒岡州立大學(xué)(OSU)Adam Schultz教授的指導(dǎo).

圖14 陸地三維模型在x=6000 m的測(cè)線上的卡尼亞視電阻率響應(yīng)擬斷面圖Fig.14 The Cagniard apparent resistivity response pseudo-sections of the survey line at x=6000 m for the land 3D model

圖15 陸地三維模型在y=0 m的測(cè)線上的卡尼亞視電阻率響應(yīng)擬斷面圖Fig.15 The Cagniard apparent resistivity response pseudo-sections of the survey line at y=0 m for the land 3D model

圖16 陸地三維模型在頻率f=4 Hz時(shí)的卡尼亞視電阻率響應(yīng)頻率切片圖Fig.16 The Cagniard apparent resistivity response frequency slices at f=4 Hz for the land 3D model

Abubakar A, Habashy T M, Druskin V L, et al. 2008. 2.5D forward and inverse modeling for interpreting low-frequency electromagnetic measurements.Geophysics, 73(4): F165-F177.

Amestoy P R, Guermouche A, L′Excellent J, et al. 2006. Hybrid scheduling for the parallel solution of linear systems.ParallelComputing, 32(2): 136-156.

Avdeev D B. 2005. Three-dimensional electromagnetic modelling and inversion from theory to application.SurveysinGeophysics, 26(6): 767-799. Avdeev D, Knizhnik S. 2009. 3D integral equation modeling with a linear dependence on dimensions.Geophysics, 74(5): F89-F94. B?rner R U. 2010. Numerical modelling in Geo-Electromagnetics: advances and challenges.SurveysinGeophysics, 31(2): 225-245. Chen G B, Wang H N, Yao J J, et al. 2009. Modeling of electromagnetic responses of 3-D electrical anomalous body in a layered anisotropic earth using integral equations.ChineseJ.Geophys. (in Chinese), 52(8): 2174-2181, doi: 10.3969/j.issn.0001-5733.2009.08.028.

Constable S C, Weiss C J. 2006. Mapping thin resistors and hydrocarbons with marine EM methods: Insights from 1D modeling.Geophysics, 71(2): G43-G51.

Constable S C. 2010. Ten years of marine CSEM for hydrocarbon exploration.Geophysics, 75(5): 75A67-75A81.

da Silva N V, Morgan J V, MacGregor L, et al. 2012. A finite element multifrontal method for 3D CSEM modeling in the frequency domain.Geophysics, 77(2): E101-E115.

Davis P J, Rabinowitz P. 2007. Methods of Numerical Integration. 2nd ed. New York: Dover Publications.

Deng J Z. 2011. Studies of three dimensional stagged-grid finite difference for controlled source audio-frequency magnetotelluric numerical simulation [Ph. D. thesis] (in Chinese). Beijing: China University of Geosciences (Beijing).

Di Q Y, Wang R. 2008. CSAMT Forward Modeling and Inversion and Its Application (in Chinese). Beijing: Science Press. Grayver A V. 2013. Three-dimensional controlled-source electromagnetic inversion using modern computational concepts [Ph. D. thesis]. Berlin: Free University of Berlin.

Han B, Hu X Y, He Z X, et al. 2012. Mathematical classification of magnetotelluric inversion methods.OilGeophysicalProspecting(in Chinese), 47(1): 177-187.

Han B, Hu X Y, Schultz A, et al. 2015. Three-dimensional forward modeling of the marine controlled-source electromagnetic field with complex source geometries.ChineseJ.Geophys. (in Chinese), 58(3): 1059-1071, doi: 10.6038/cjg20150330.

Hu X Y, Peng R H, Wu G J, et al. 2013. Mineral Exploration using CSAMT data: Application to Longmen region metallogenic belt, Guangdong Province, China.Geophysics, 78(3): B111-B119. Jahandari H, Farquharson C G. 2014. A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids.Geophysics, 79(6): E287-E302. Key K. 2009. 1D inversion of multicomponent, multifrequency marine CSEM data: methodology and synthetic studies for resolving thin resistive layers.Geophysics, 74(2): F9-F20.

Li J H, Zhu Z Q, Liu S C, et al. 2011. 3D numerical simulation for the transient electromagnetic field excited by the central loop based on the vector finite-element method.JournalofGeophysicsandEngineering, 8(4): 560-567. Luo M, Li Y, Liu Y. 2014. 2.5D controlled-source electromagnetic modeling by an adaptive finite-element method with an arbitrarily oriented finite length dipole source. 22nd EM Induction Workshop. Abstracts. Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences.GeophysicalProspecting, 43(8): 1021-1042.

Oldenburg D W, Haber E, Shekhtman R. 2008. Forward modelling and inversion of multi-source TEM data. ∥ 78th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 559-563. Puzyrev V, Koldan J, de la Puente J, et al. 2013. A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling.Geophys.J.Int., 193(2): 678-693. Saad Y. 2003. Iterative Methods for Sparse Linear Systems. 2nd ed. Philadelphia: SIAM.

Schenk O, G?rtner K. 2004. Solving unsymmetric sparse systems of linear equations with PARDISO.FutureGenerationComputerSystems, 20(3): 475-487.

Schwarzbach C, B?rner R U, Spitzer K. 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics—a marine CSEM example.Geophys.J.Int., 187(1): 63-74. Shen J S. 2003. Modeling of 3-D electromagnetic responses in frequency domain by using staggered-grid finite difference method.ChineseJ.Geophys. (in Chinese), 46(2): 281-289.

Streich R. 2009. 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy.Geophysics, 74(5): F95-F105.

Streich R, Becken M, Ritter O. 2011. 2.5D controlled-source EM modeling with general 3D source geometries.Geophysics, 76(6): F387-F393.

Tan H D, Wei W B, Deng M, et al. 2004. General use formula in MT tensor impedance.OilGeophysicalProspecting(in Chinese), 39(1): 113-116. Wang R, Di Q Y, Wang M Y, et al. 2009. Research on the effect of 3D body between transmitter and receivers on CSAMT response using Integral Equation method.ChineseJ.Geophys. (in Chinese), 52(6): 1573-15821, doi: 10.3969/j.issn.0001-5733.2009.06.019.

Weiss C J. 2001. A matrix-free approach to solving the fully 3D electromagnetic induction problem. ∥71th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1451-1454. Weiss C J, Constable S C. 2006. Mapping thin resistors and hydrocarbons with marine EM methods, Part II — Modeling and analysis in 3D.Geophysics, 71(6): G321-G332.

Wu G J, Hu X Y, Liu H. 2010. Progress in CSAMT three-dimensional forward numerical simulation.ProgressinGeophys. (in Chinese), 25(5): 1795-1801, doi: 10.3969/j.issn.1004-2903.2010.05.037.

Xu Z F, Wu X P. 2010. Controlled source electromagnetic 3-D modeling in frequency domain by finite element method.ChineseJ.Geophys. (in Chinese), 53(8): 1931-1939, doi: 10.3969/j.issn.0001-5733.2010.08.019.

Yang B, Xu Y X, He Z X, et al. 2012. 3D frequency-domain modeling of marine controlled source electromagnetic responses with topography using finite volume method.ChineseJ.Geophys. (in Chinese), 55(4): 1390-1399, doi: 10.6038/j.issn.0001-5733.2012.04.035.

Zaslavsky M, Druskin V, Davydycheva S, et al. 2011. Hybrid finite-difference integral equation solver for 3D frequency domain anisotropic electromagnetic problems.Geophysics, 76(2): F123-F137.

附中文參考文獻(xiàn)

陳桂波, 汪宏年, 姚敬金等. 2009. 用積分方程法模擬各向異性地層中三維電性異常體的電磁響應(yīng). 地球物理學(xué)報(bào), 52(8): 2174-2181, doi: 10.3969/j.issn.0001-5733.2009.08.028.

鄧居智. 2011. 可控源音頻大地電磁法三維交錯(cuò)采樣有限差分?jǐn)?shù)值模擬研究[博士論文]. 北京: 中國(guó)地質(zhì)大學(xué)(北京).

底青云, 王若. 2008. 可控源音頻大地電磁數(shù)據(jù)正反演及方法應(yīng)用. 北京: 科學(xué)出版社.

韓波, 胡祥云, 何展翔等. 2012. 大地電磁反演方法的數(shù)學(xué)分類. 石油地球物理勘探, 47(1): 177-187.

韓波, 胡祥云, Schultz A等. 2015. 復(fù)雜場(chǎng)源形態(tài)的海洋可控源電磁三維正演. 地球物理學(xué)報(bào), 58(3): 1059-1071, doi: 10.6038/cjg20150330.

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

譚捍東, 魏文博, 鄧明等. 2004. 大地電磁法張量阻抗通用計(jì)算公式. 石油地球物理勘探, 39(1): 113-116.

王若, 底青云, 王妙月等. 2009. 用積分方程法研究源與勘探區(qū)之間的三維體對(duì)CSAMT觀測(cè)曲線的影響. 地球物理學(xué)報(bào), 52(6): 1573-1582, doi: 10.3969/j.issn.0001-5733.2009.06.019.

吳桂桔, 胡祥云, 劉慧. 2010. CSAMT三維正演數(shù)值模擬研究進(jìn)展. 地球物理進(jìn)展, 25(5): 1795-1801, doi: 10.3969/j.issn.1004-2903.2010.05.037.

徐志鋒, 吳小平. 2010. 可控源電磁三維頻率域有限元模擬. 地球物理學(xué)報(bào), 53(8): 1931-1939, doi: 10.3969/j.issn.0001-5733.2010.08.019.

楊波, 徐義賢, 何展翔等. 2012. 考慮海底地形的三維頻率域可控源電磁響應(yīng)有限體積法模擬. 地球物理學(xué)報(bào), 55(4): 1390-1399, doi: 10.6038/j.issn.0001-5733.2012.04.035.

(本文編輯 胡素芳)

3-D frequency-domain CSEM modeling using a parallel direct solver

HAN Bo， HU Xiang-Yun*， HUANG Yi-Fan， PENG Rong-Hua， LI Jian-Hui， CAI Jian-Chao

HubeiSubsurfaceMulti-scaleImagingKeyLaboratory,InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China

Three-dimensional modeling of electromagnetic data is a computationally demanding problem. For frequently-used numerical techniques such as finite-element and finite-difference methods, solving the large linear systems arising from the discretization of Maxwell's equations is a key step which has a major impact on the applicability of the solution, and it has always been a research topic to solve the linear equations efficiently, robustly and accurately. A 3D modeling scheme based on direct solutions of the linear system is presented for frequency domain controlled-source electromagnetic (CSEM) surveys. The Helmholtz equation in terms of secondary electric fields is discretized using a finite-volume (FV) method over a staggered grid. Taking advantage of recent developments in numerical algorithms and the availability of computational resources, the resulting linear system of FV equations is solved directly using the massively parallel solver, namely MUMPS, instead of the most commonly used linear solvers, i.e. Krylov subspace iterative techniques. The direct solver carries out an LU (and possibly LDLT) decomposition of the system matrix and then computes solutions efficiently by applying forward and backward substitutions.To evaluate the computational performance of the direct solver, a series of numerical tests based on synthetic 1D models were conducted, and the results indicate that (1) Normalized residuals of solutions are almost independent of the conductivity value assigned to air layers but increase rapidly as the frequency value decreases. Nevertheless, the order of magnitude of the largest normalized residual is as small as 10-11. At the same time, although the matrix factorization time varies as either the air conductivity or the frequency changes, the variation is only a fraction of the total run time. (2) Both the execution time and required memory increase rapidly (more than linearly) with increasing grid sizes. (3) By executing MUMPS in parallel over multiple processors, not only the total run time but also the average memory used per processor can be reduced a lot. However, the total memory requirement increases with the number of processes. The scalability of MUMPS is limited. Additional numerical experiments considering specific survey settings were done to demonstrate the reliability and effectiveness of the code, and the results are as follows: (1) The FV numerical solutions show excellent agreement with semi-analytic solutions for the 1D models. (2) The computation time of a multitransmitter problem is comparable to that of a single-transmitter problem. (3) Reasonable modeling results of 3D models can be obtained for both typical land and marine survey scenarios.In summary, compared with iterative linear solvers, the direct solvers generally benefit CSEM modeling in three aspects. The first is they often provide more accurate solutions. The second is that direct solvers are much more stable for ill-conditioned linear systems, which are almost inevitable because of large electrical conductivity contrasts and/or non-uniform grid. The last is in multitransmitter problems, only a single matrix factorization is necessary, and multiple solutions can be achieved very easily by reusing the factors. The presented 3D CSEM modeling scheme, which employs the MUMPS direct solver, possesses all these advantages. In addition, solving linear systems can be executed in parallel to speed up the computation and to reduce the average memory used per node, although the parallel scalability of MUMPS is limited. In spite of the fact that matrix factorizations for large models can entail tremendous computational cost, it can be anticipated that direct solvers will be used more and more widely as the development of both numerical algorithms and computers.

Controlled-source electromagnetics; 3-D modeling; LU factorization; Direct solver

國(guó)家自然科學(xué)基金項(xiàng)目(41274077、41474055)和中國(guó)地質(zhì)調(diào)查局項(xiàng)目(12120113101800)聯(lián)合資助.

韓波，男，1987年生，博士研究生，研究方向?yàn)殡姶欧ǖ恼菖c反演模擬.E-mail：hanbo1735@163.com

*通訊作者胡祥云，男，1966年生，教授，主要從事電磁法的理論與應(yīng)用研究.E-mail：xyhu@cug.edu.cn

10.6038/cjg20150816.

10.6038/cjg20150816

P631

2015-02-12，2015-06-16收修定稿

韓波，胡祥云，黃一凡等. 2015. 基于并行化直接解法的頻率域可控源電磁三維正演.地球物理學(xué)報(bào)，58(8):2812-2826,

Han B, Hu X Y, Huang Y F, et al. 2015. 3-D frequency-domain CSEM modeling using a parallel direct solver.ChineseJ.Geophys. (in Chinese),58(8):2812-2826,doi:10.6038/cjg20150816.