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

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

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

?

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

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

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

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

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

1 引言

頻率域可控源電磁法(controlled-source electromagnetic methods, CSEM)通過觀測人工場源在地層、空氣和海水等介質(zhì)中激發(fā)的電磁場來獲取介質(zhì)的電性分布，長期以來一直是金屬礦勘探中最重要的地球物理方法(底青云和王若，2008；Hu et al.，2013)，并在最近十多年逐漸成為海洋油氣勘探中一種不可或缺的地球物理技術(shù)(Constable，2010).為了通過CSEM觀測數(shù)據(jù)盡可能真實地還原出介質(zhì)的電性結(jié)構(gòu)以提升勘探效果，實施三維勘探并對數(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).對這幾種方法的評述可以在很多文獻(xiàn)中找到(Avdeev，2005；吳桂桔等，2010； B?rner，2010).

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

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

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

2 方法理論

2.1 控制方程

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

(1)

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

(2)

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

(3)

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

(4)

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

2.2 有限體積法

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

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

(5)

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

(6)

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

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

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

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

Ae=b,

(7)

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

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

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

2.4 計算流程

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

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

3 直接解法的性能測試

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

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

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

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

表1 對海洋1D模型，頻率與空氣層電導(dǎo)率在較大范圍內(nèi)變化時直接解法的精度和矩陣分解時間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)存需求與計算時間

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

表2 對陸地1D模型，頻率與空氣層電導(dǎo)率在較大范圍內(nèi)變化時直接解法的精度和矩陣分解時間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 對海洋1D模型，對于不同的模型網(wǎng)格剖分?jǐn)?shù)直接解法的(a)內(nèi)存需求與(b)解方程時間Fig.5 (a) Memory requirement and (b) solving time of the direct solver for the marine 1D model for various grid sizes

3.3 并行可伸縮性

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

4 三維正演程序的有效性檢驗

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

4.1 與解析解的對比

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

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

圖6 對海洋層狀模型，直接解法解方程的時間與內(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 對海洋層狀模型，網(wǎng)格為86×50×57時水平(a, b)電場和(c, d)磁場的SFV數(shù)值解(圓圈)與解析解(線條)對比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 對于海洋層狀模型，在3種不同網(wǎng)格剖分下SFV解相對于解析解的誤差Fig.8 Relative errors of SFV solutions compared to analytic solutions for the marine 1D model for three different grid sizes

圖9 對陸地層狀模型，坐標(biāo)為(8000, 0)的測點的水平(a, b)電場、(c, d)磁場和(e, f)卡尼亞視電阻率的SFV解(圓圈)與解析解(線條)對比，以及(g, h)視電阻率的SFV解相對于解析解的誤差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

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

(8)

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

4.2 多場源問題的模擬

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

4.3 海洋三維模型

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

分別考慮海水無窮深、不含空氣層和海水1000 m深、包含空氣層這兩種情況.參數(shù)設(shè)置為：場源沿x方向，位于高阻體x負(fù)向邊界中點的正上方，海底上方50 m，坐標(biāo)為(0, 0, -50 m)；發(fā)射頻率為1 Hz；接收器沿x方向布設(shè)、分布于海底y=0的測線上(inline觀測).對高阻體截面邊長(設(shè)為a)為1 km、2 km和5 km的三種情況進(jìn)行了正演.所得到的水平電場Ex分量如圖12所示，其中還給出了不含高阻體的半空間(a=0)和高阻體在水平方向無限延伸(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ū)分開來，振幅曲線在對數(shù)坐標(biāo)系中的斜率趨于常數(shù)，說明振幅趨于按指數(shù)衰減.a=1 km與半空間模型的響應(yīng)幾乎重合，說明該尺寸的異常體難以分辨；a=2 km與半空間的響應(yīng)雖然在大收發(fā)距時的衰減速度保持一致，但二者的絕對差值較為明顯，說明該尺寸的異常體容易分辨；a=5 km與半空間的響應(yīng)差別非常明顯，且在收發(fā)距大于高阻體邊長時與高阻體無限延伸(a=∞)的響應(yīng)區(qū)分開來，衰減速度逐漸趨于半空間模型的衰減速度.根據(jù)以上分析，可以再次確認(rèn)Constable和Weiss

圖10 對于海洋層狀模型，SFV正演總體計算時間及各主要階段消耗的時間隨右端項數(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é)論：對于此類水平薄板狀高阻體，只有當(dāng)其水平尺寸至少為埋深的2倍時，才可以利用CSEM將其分辨出來.

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

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

圖12 三維儲層模型的儲層截面邊長變化時的Ex響應(yīng)(a, b)為不考慮空氣層、假設(shè)海水無限深的情況，(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 陸地三維模型

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

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

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

5 結(jié)論與討論

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

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

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

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

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

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

圖14 陸地三維模型在x=6000 m的測線上的卡尼亞視電阻率響應(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的測線上的卡尼亞視電阻率響應(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時的卡尼亞視電阻率響應(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é)報, 52(8): 2174-2181, doi: 10.3969/j.issn.0001-5733.2009.08.028.

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

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

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

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

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

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

王若, 底青云, 王妙月等. 2009. 用積分方程法研究源與勘探區(qū)之間的三維體對CSAMT觀測曲線的影響. 地球物理學(xué)報, 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é)報, 53(8): 1931-1939, doi: 10.3969/j.issn.0001-5733.2010.08.019.

楊波, 徐義賢, 何展翔等. 2012. 考慮海底地形的三維頻率域可控源電磁響應(yīng)有限體積法模擬. 地球物理學(xué)報, 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

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

韓波，男，1987年生，博士研究生，研究方向為電磁法的正演與反演模擬.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é)報，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.