彭榮華， 胡祥云， 韓波， 蔡建超

1 中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院， 武漢　430074 2 不列顛哥倫比亞大學(xué)地球、海洋與大氣科學(xué)學(xué)院， 溫哥華 V6T 1Z4， 加拿大 3 中國(guó)海洋大學(xué)海洋地球科學(xué)學(xué)院， 青島　266100

?

基于擬態(tài)有限體積法的頻率域可控源三維正演計(jì)算

彭榮華1,2， 胡祥云1*， 韓波3， 蔡建超1

1 中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院， 武漢430074 2 不列顛哥倫比亞大學(xué)地球、海洋與大氣科學(xué)學(xué)院， 溫哥華 V6T 1Z4， 加拿大 3 中國(guó)海洋大學(xué)海洋地球科學(xué)學(xué)院， 青島266100

大規(guī)模地球物理電磁數(shù)據(jù)的定量解釋需要發(fā)展高效、穩(wěn)定的三維正反演算法.本文通過求解離散化的三維電場(chǎng)矢量Helmholtz方程，實(shí)現(xiàn)了基于有限體積法的頻率域可控源電磁(CSEM)三維正演算法.為模擬具有強(qiáng)電性差異的三維電性介質(zhì)，該算法采用擬態(tài)有限體積法(MFV)對(duì)Maxwell方程組進(jìn)行離散化；另外，為獲得穩(wěn)定、高精度的正演數(shù)值結(jié)果，采用直接矩陣分解技術(shù)來求解離散所得到的大型稀疏線性方程組.對(duì)于具有多個(gè)發(fā)射源的CSEM測(cè)量來說，一次矩陣分解結(jié)果能夠用于同頻率下所有場(chǎng)源的正演計(jì)算.為降低場(chǎng)源奇異性及邊界條件對(duì)數(shù)值精度的影響，采用虛擬場(chǎng)源校正技術(shù)，避免了散射場(chǎng)公式中在構(gòu)建場(chǎng)源項(xiàng)時(shí)所需的大量時(shí)間.對(duì)于具有多個(gè)頻率的CSEM的模擬計(jì)算，采用分頻并行策略來加快三維正演計(jì)算.最后，通過與一維層狀模型及三維模型的數(shù)值結(jié)果的對(duì)比驗(yàn)證了本文所開發(fā)的正演算法對(duì)頻率域CSEM模擬計(jì)算的準(zhǔn)確性及有效性，表明該正演算法能夠有效應(yīng)用于三維介質(zhì)的數(shù)值計(jì)算.另外，對(duì)于多頻率CSEM的并行測(cè)試結(jié)果表明基于分頻并行策略的并行計(jì)算能夠顯著地降低正演計(jì)算時(shí)間.

可控源電磁法； 有限體積法； 虛擬場(chǎng)源校正技術(shù)； 三維正演； 直接分解法

An important step in this 3D modeling scheme is to solve the large linear equation system resulting from MFV discretization. In order to obtain stable and accurate numerical solutions, the linear equation system is solved by a direct-matrix solver, namely MUMPS, which is more robust than commonly used iterative solvers (e.g, Krylov subspace iterative techniques) for numerically difficult cases. The algorithm is very suitable for multi-source CSEM modeling by separating the solving process into a single expensive matrix factorization and relatively inexpensive forward backward substitutions for many right-hand sides.

For total-field formulation, dense gridding in the vicinity of source points is usually required to mitigate source singularities, and extra padding cells at the boundaries are needed to meet boundary conditions. Both of them can quickly increase the size of the linear equation system to be solved, making it computationally expensive for direct solutions. To avoid substantial increase of the size of the discretized model, a source correction technique is applied to reduce source singularities and boundary effects. In addition, considering the independence of computation for different frequencies, a parallel forward modeling scheme based on frequency partition is implemented to speed up the simulation of multi-frequency CSEM problems.

Numerical experiments have been carried out to evaluate the performance of our forward modeling algorithm. Numerical solutions using this algorithm show good agreement with quasi-analytic solutions for 1D layered models, and numerical errors are significantly reduced with source correction in the vicinity of source points. In addition, comparison of simulated data generated by our algorithm to published 3D data for a typical marine 3D model validates our algorithm. Besides, statistical results demonstrate that parallel computing based on simple frequency partition approach can achieve nearly linear speedup due to the independence of computation between frequencies for multi-frequency CSEM simulation.

1　引言

頻率域可控源電磁法(controlled-source electromagnetic method，CSEM)由于具有勘探深度大，分辨力高且野外抗干擾能力強(qiáng)等優(yōu)點(diǎn)，經(jīng)過30多年快速發(fā)展，目前已被廣泛應(yīng)用于包括陸地、海洋及航空等領(lǐng)域內(nèi)的油氣、礦產(chǎn)資源勘探和近地表地質(zhì)勘查.隨著電磁儀器、數(shù)據(jù)采集技術(shù)以及解釋方法的快速發(fā)展，可控源電磁法在復(fù)雜地質(zhì)環(huán)境下的三維勘探成為當(dāng)前電磁勘探的重要趨勢(shì)及研究熱點(diǎn).在電磁數(shù)據(jù)的處理解釋中，反演是必不可少的步驟，而正演計(jì)算是反演的核心.因此，電磁數(shù)據(jù)的定量解釋首先需要發(fā)展有效的正演算法.從數(shù)值模擬的角度來說，目前常用于三維電磁正演的數(shù)值方法主要包括：積分方程法(Avdeev et al., 2002；Hursán and Zhdanov, 2002；Zhdanov et al., 2006；Avdeev and Knizhnik, 2009)、有限差分法(Mackie et al.,1994；Newman and Alumbaugh, 1995；Smith, 1996；Weiss and Newman, 2002；沈金松, 2003； Streich, 2009)、有限體積法(Haber and Ascher, 2001；Weiss and Constable, 2006； 楊波等, 2012； Jahandari and Farquharson, 2014；韓波等, 2015a)和有限單元法(Badea et al., 2001； 徐志鋒和吳小平, 2010；Schwarzbach et al., 2011；Puzyrev et al., 2013； Ansari and Farquharson, 2014；李勇等, 2015；楊軍等, 2015).另外，其他方法如偽譜法(Huang et al., 2010; Liu et al., 2013)、Lanczos譜分解法(Knizhnerman et al.,2009)也被廣泛用于三維電磁正演問題的求解.有關(guān)電磁法三維正演方面進(jìn)展的綜述可以詳見(Avdeev, 2005；B?rner, 2010).總體來說，目前三維正演算法研究的重點(diǎn)在于提高精度和計(jì)算效率.

CSEM正演中常采用散射場(chǎng)方法來避免場(chǎng)源的奇異性，然而構(gòu)建散射場(chǎng)的場(chǎng)源項(xiàng)需要計(jì)算幾乎每個(gè)離散網(wǎng)格采樣點(diǎn)處的一次電場(chǎng)或磁場(chǎng).對(duì)于三維問題，離散網(wǎng)格中的電場(chǎng)或磁場(chǎng)采樣點(diǎn)個(gè)數(shù)大約是網(wǎng)格數(shù)的3倍，很容易達(dá)到幾百萬甚至上千萬.因此構(gòu)建散射場(chǎng)場(chǎng)源往往需要大量的時(shí)間，隨著發(fā)射場(chǎng)源的增加，甚至有可能超過解線性方程需要的時(shí)間(韓波等, 2015b).相較之下，總場(chǎng)方法直接對(duì)物理場(chǎng)源離散化，不存在這個(gè)問題.但總場(chǎng)方法中一般需要在物理場(chǎng)源附近加密網(wǎng)格來削弱場(chǎng)源奇異性的影響，大量網(wǎng)格數(shù)同樣會(huì)導(dǎo)致大量的計(jì)算時(shí)間.除了加密網(wǎng)格以外，還可采用場(chǎng)源校正技術(shù)等手段來降低場(chǎng)源奇異性的影響.

正演計(jì)算中的一個(gè)關(guān)鍵步驟是求解離散化后產(chǎn)生的大型線性方程組.在早期受限于計(jì)算機(jī)內(nèi)存，三維正演中大型線性方程組幾乎都采用對(duì)內(nèi)存需求極低的Krylov子空間迭代法.隨著大型稀疏矩陣分解算法的不斷優(yōu)化(Amestoy et al., 2001，2006； Li, 2005； Schenk and G?rtner, 2006)以及計(jì)算機(jī)硬件技術(shù)的快速發(fā)展，利用直接解法求解三維正演的大型線性方程組成為可能.與迭代解法相比，對(duì)于同一線性方程組，直接解法通常需要消耗更多的計(jì)算內(nèi)存和更長(zhǎng)的計(jì)算時(shí)間.盡管如此，直接解法相對(duì)于迭代解法有兩大顯著優(yōu)點(diǎn)：一是求解時(shí)間和精度基本不受系數(shù)矩陣的條件數(shù)影響，因此網(wǎng)格的剖分方式或介質(zhì)的電性差異對(duì)正演算法影響很小，使得求解精度更高，正演算法更穩(wěn)定；另一個(gè)是直接解法的主要開銷集中在系數(shù)矩陣的分解階段，對(duì)于單一頻率多個(gè)場(chǎng)源位置的問題，只需進(jìn)行一次矩陣分解加上多次計(jì)算量極小的前代-回代過程即可，而不像迭代解法對(duì)每個(gè)頻率每個(gè)場(chǎng)源都要單獨(dú)求解，因此直接解法特別適合具有多個(gè)場(chǎng)源的CSEM三維模擬(Oldenburg et al.，2013).正因?yàn)檫@些優(yōu)點(diǎn)，基于矩陣分解的直接解法最近幾年開始應(yīng)用于CSEM三維正演計(jì)算(Streich, 2009；Chung et al., 2014；楊軍等, 2015).

本文采用擬態(tài)有限體積法(Mimetic Finite Volume, MFV)，通過求解離散化的三維電場(chǎng)矢量Helmholtz方程，實(shí)現(xiàn)了頻率域CSEM三維正演算法.為獲得穩(wěn)定且高精度的數(shù)值結(jié)果，采用直接矩陣分解法來求解離散所得到的大型稀疏線性方程組；為消弱場(chǎng)源奇異性及邊界條件的影響，采用場(chǎng)源校正技術(shù)來提高求解精度.對(duì)于多頻率正演問題，采用分頻并行策略來降低正演計(jì)算時(shí)間；最后通過多個(gè)理論合成模型正演計(jì)算測(cè)試，驗(yàn)證了本文所開發(fā)的正演算法對(duì)頻率域CSEM模擬計(jì)算的準(zhǔn)確性及有效性.

2　方法理論

2.1控制方程

在CSEM所采用的頻率范圍內(nèi)，位移電流的影響可以忽略不計(jì).假定時(shí)諧因子取為eiwt，則頻率域Maxwell方程組可表示為

(1)

(2)

其中，E為電場(chǎng)強(qiáng)度(V·m-1)，B為磁感應(yīng)強(qiáng)度(T)，ω為角頻率(rad/s).σ為介質(zhì)電導(dǎo)率(S·m-1)，對(duì)于各向同性介質(zhì)，σ為標(biāo)量.μ0為真空中磁導(dǎo)率(H·m-1). Js為外部激發(fā)場(chǎng)源的電流密度(A·m-2).對(duì)(1)式取旋度并將其帶入(2)式，消去磁感應(yīng)強(qiáng)度B，得到關(guān)于電場(chǎng)E的二階矢量Helmholtz方程為

(3)

2.2擬態(tài)有限體積法

為求解方程(3)，本文采用擬態(tài)有限體積法(MFV)來對(duì)連續(xù)的Maxwell方程組進(jìn)行離散化.MFV方法的最大特點(diǎn)是能夠?qū)B續(xù)的微分算符進(jìn)行精確的離散模擬，并確保離散化后的矢量場(chǎng)仍然滿足其對(duì)應(yīng)的連續(xù)形式的矢量性質(zhì)及物理特性(Hyman and Shashkov, 1999a, b).如離散化的電磁場(chǎng)的能量守恒能夠自動(dòng)滿足，有效避免了偽解的產(chǎn)生.對(duì)于正交規(guī)則網(wǎng)格及各向同性介質(zhì)來說，該離散化方法與常用的交錯(cuò)網(wǎng)格有限差分法(Staggered Finite-Difference, SFD)相一致(Smith, 1996).盡管最近也有將非正交交錯(cuò)網(wǎng)格有限差分法應(yīng)用于電磁數(shù)值模擬的研究(邱稚鵬等，2013)，但MFV適用范圍更廣，其不僅對(duì)于非正交網(wǎng)格同樣適用(Hyman and Shashkov, 1997)，而且對(duì)于具有高度電性差異及各向異性介質(zhì)，MFV能夠自然地得到對(duì)稱的離散化形式(Hyman et al., 1997；Haber and Ruthotto, 2014)，這對(duì)于復(fù)雜三維介質(zhì)的模擬十分有效.

為獲得Maxwell方程組的離散形式，首先考慮方程(1)和(2)的弱解形式為(Hyman and Shashkov, 1999a)：

(4)

(5)

其中W和F分別為與電場(chǎng)E和磁感應(yīng)強(qiáng)度B屬于相同Hilbert空間的任意測(cè)試函數(shù)(Hyman and Shashkov, 1999a).對(duì)(5)式左端項(xiàng)進(jìn)行分部積分，可得公式為

(6)

(7)

將(6—7)式代入到(5)式中，得到公式為

(8)

由于MFV是在離散單元的控制體積內(nèi)進(jìn)行積分，因此正交結(jié)構(gòu)網(wǎng)格(Weiss and Constable, 2006；韓波等, 2015a)、半結(jié)構(gòu)網(wǎng)格(Haber and Heldmann, 2007)和非結(jié)構(gòu)網(wǎng)格(Jahandari and Farquharson, 2014)都能用來對(duì)(4)式和(8)式的積分進(jìn)行離散化.本文采用基于Yee網(wǎng)格(Yee, 1966)的交錯(cuò)網(wǎng)格.圖1展示了交錯(cuò)網(wǎng)格電磁場(chǎng)采樣方式及電場(chǎng)分量和磁場(chǎng)分量各自的控制體積.其中電場(chǎng)定義在單元棱邊的中心，磁感應(yīng)強(qiáng)度定義在單元面的中心，介質(zhì)電導(dǎo)率和磁導(dǎo)率定義在單元中心.

(10)

由于F和W為任意的網(wǎng)格測(cè)試函數(shù)，將(9)式代入(10)式并進(jìn)行簡(jiǎn)化可得(3)式的離散形式為

(CurlTMfμCurl+iωMe σ)E=S，

(11)

其中

(11)式可簡(jiǎn)化成大型線性方程組為

(12)

其中系數(shù)矩陣A為稀疏、正定、對(duì)稱復(fù)矩陣.圖2給出了網(wǎng)格數(shù)為5×5×5的模型的MFV離散所得到的系數(shù)矩陣A的稀疏結(jié)構(gòu).

圖1　電磁場(chǎng)分量交錯(cuò)采樣方式及其控制體積(a) 電場(chǎng)x分量采樣位置及其控制體積； (b) 磁場(chǎng)x分量采樣位置及其控制體積.Fig.1　Staggered discretization of MFV(a) Integration volume for x-component of electric field； (b) Integration volume for x-component of magnetic field.

圖2　網(wǎng)格數(shù)為5×5×5的模型的MFV離散線性方程系數(shù)矩陣的稀疏結(jié)構(gòu)Fig.2　Sparse structure of coefficient matrix A resulted from MFV discretization on a 5×5×5 staggered grid

2.3虛擬場(chǎng)源校正技術(shù)

在可控源電磁法中，由于外加場(chǎng)源的存在，待求解的電磁場(chǎng)在場(chǎng)源附近會(huì)發(fā)生急劇變化，場(chǎng)源的奇異性會(huì)極大降低數(shù)值計(jì)算結(jié)果的精度.目前主要有兩種處理方法來降低場(chǎng)源奇異性的影響. 第一種方法是采用散射場(chǎng)公式(Newman and Alumbaugh, 1995；Streich, 2009；韓波等, 2015a)，將待求解的電磁場(chǎng)分解為由參考模型(一般為均勻空間或一維層狀模型)在外加場(chǎng)源激發(fā)下所產(chǎn)生的一次場(chǎng)和三維模型所產(chǎn)生的二次場(chǎng)的疊加.由于外加場(chǎng)源已包含在一次場(chǎng)的計(jì)算中，二次場(chǎng)的場(chǎng)源不再具有奇異性.另一種是基于電磁場(chǎng)總場(chǎng)公式，通過對(duì)場(chǎng)源點(diǎn)附近進(jìn)行網(wǎng)格加密的方法來減少場(chǎng)源奇異性的影響(Plessix et al., 2007).

在散射場(chǎng)公式中，為避免計(jì)算二次場(chǎng)時(shí)出現(xiàn)場(chǎng)源的奇異性，必須確保選擇的參考模型的電阻率與場(chǎng)源點(diǎn)附近的電阻率一致.對(duì)于海洋CSEM測(cè)量來說，由于發(fā)射場(chǎng)源一般位于海水中，一般容易滿足該條件.但對(duì)于陸地CSEM測(cè)量來說，通常CSEM勘探會(huì)在較大測(cè)區(qū)內(nèi)布設(shè)多個(gè)發(fā)射位置(Streich et al., 2011)，當(dāng)場(chǎng)源處的電阻率發(fā)生變化時(shí)，需要求解不同參考模型的一次電磁場(chǎng)響應(yīng)，從而增加計(jì)算量.另外，在散射場(chǎng)公式中，每個(gè)場(chǎng)源處都需要計(jì)算三維空間內(nèi)所有電磁場(chǎng)采樣點(diǎn)處的一次電磁場(chǎng)來構(gòu)建場(chǎng)源項(xiàng)或進(jìn)行總場(chǎng)合成.隨著發(fā)射場(chǎng)源位置的增加，場(chǎng)源項(xiàng)的構(gòu)建所消耗的時(shí)間會(huì)快速增加，甚至?xí)^系數(shù)矩陣分解所消耗的時(shí)間(韓波等, 2015b).與之對(duì)比，在總場(chǎng)公式中，網(wǎng)格加密的策略雖然能夠降低場(chǎng)源奇異性的影響，但與此同時(shí)會(huì)顯著地增加待求解的未知量個(gè)數(shù)，從而增大計(jì)算開銷，降低求解效率.另外，為滿足邊界條件，無論是采用散射場(chǎng)公式還是總場(chǎng)公式，通常都需要將計(jì)算區(qū)域的外邊界設(shè)置得足夠遠(yuǎn)，這使得待求解的方程組會(huì)非常巨大.

為降低場(chǎng)源奇異性及邊界條件對(duì)數(shù)值精度的影響，本文采用虛擬場(chǎng)源校正技術(shù)(Pidlisecky et al., 2007)來對(duì)總場(chǎng)控制方程(11)式中的源項(xiàng)進(jìn)行校正，避免了在場(chǎng)源點(diǎn)附近進(jìn)行網(wǎng)格局部加密所造成的計(jì)算量的增大.Pidlisecky等(2007)、Pidlisecky和Knight(2008)將該校正技術(shù)運(yùn)用到直流電阻率的正演計(jì)算中，取得了很好的效果.為得到校正后的場(chǎng)源項(xiàng)，考慮存在解析解的均勻半空間模型σ0，利用上節(jié)發(fā)展的有限體積法對(duì)其離散化可得公式為

(13)

其中系數(shù)矩陣A(σ0)和右端項(xiàng)S0均由(11)式確定.假設(shè)E0是模型σ0在外加場(chǎng)源Js激勵(lì)下的網(wǎng)格采樣點(diǎn)處的電場(chǎng)響應(yīng)，可通過解析求解很容易得到(Ward and Hohmann, 1988)，由于離散化誤差及邊界條件的影響，待求解的電場(chǎng)分量E與真實(shí)的電場(chǎng)響應(yīng)E0會(huì)有一定差異.此時(shí)，我們可以得到場(chǎng)源校正量為

(14)

將(14)式所得到的場(chǎng)源校正量添加到(12)式右端，得公式為

(15)

通過求解(15)式便可以得到待求解的電場(chǎng)值.

2.4大型線性方程組的求解

由于離散后所得到的系數(shù)矩陣為稀疏、正定、對(duì)稱復(fù)矩陣，本文通過調(diào)用基于多波前分解算法的MUMPS (Amestoy et al., 2001，2006)線性運(yùn)算庫(kù)對(duì)其進(jìn)行LDLH分解，然后計(jì)算所有發(fā)射場(chǎng)源的CSEM電磁響應(yīng).有關(guān)MUMPS的運(yùn)算效率及與迭代法的對(duì)比可參考Streich(2009)和Oldenburg等(2013)，另外關(guān)于MUMPS的求解精度、穩(wěn)定性及內(nèi)存需求等性能的詳細(xì)分析，可參考韓波等(2015b).

2.5多頻率并行計(jì)算

為增大勘探深度及提高數(shù)據(jù)空間分辨率， CSEM測(cè)量中通常會(huì)采集多個(gè)頻率的數(shù)據(jù).對(duì)于三維正演來說，為降低計(jì)算時(shí)間，通常需考慮對(duì)其進(jìn)行并行計(jì)算.一種思路是利用并行矩陣分解算法對(duì)多個(gè)頻率依次進(jìn)行求解.由于矩陣并行分解算法涉及到不同進(jìn)程之間大量的通信需求，并行效率不甚理想；并且，隨著并行進(jìn)程數(shù)的增加，總的內(nèi)存消耗也會(huì)迅速增加(Pardo et al., 2012；韓波等, 2015b).對(duì)于頻率域CSEM來說，考慮到不同頻率之間的正演計(jì)算是相互獨(dú)立的特征，在計(jì)算資源允許的條件下，可以對(duì)多個(gè)頻率同時(shí)進(jìn)行并行矩陣分解計(jì)算.Grayver等(2013)的測(cè)試結(jié)果表明，分頻并行策略比采用矩陣并行分解算法依次對(duì)多個(gè)頻率進(jìn)行分解計(jì)算的效率要高很多.因此本文采用分頻并行策略，將所有頻點(diǎn)盡可能均勻地分配到所有進(jìn)程中進(jìn)行并行計(jì)算.

并行程序采用主從并行模式.其中主進(jìn)程作為控制節(jié)點(diǎn)，負(fù)責(zé)模型網(wǎng)格、收發(fā)參數(shù)的發(fā)送以及計(jì)算任務(wù)的分配，并對(duì)各子進(jìn)程的計(jì)算結(jié)果進(jìn)行回收以及結(jié)果的輸出，不參與計(jì)算；各子進(jìn)程負(fù)責(zé)從主進(jìn)程處接受計(jì)算任務(wù)并將計(jì)算結(jié)果發(fā)回給主進(jìn)程.在整個(gè)并行計(jì)算過程中，各子進(jìn)程只與主進(jìn)程進(jìn)行通信.由于本文所使用的并行計(jì)算平臺(tái)計(jì)算資源(可調(diào)用的計(jì)算節(jié)點(diǎn)數(shù)和計(jì)算內(nèi)存)的限制，在子進(jìn)程中，對(duì)于每個(gè)頻點(diǎn)的計(jì)算，并未采用矩陣并行分解算法.當(dāng)分配有多個(gè)頻點(diǎn)時(shí)，需依次對(duì)多個(gè)頻率進(jìn)行分解計(jì)算.整個(gè)并行算法的計(jì)算流程如圖3所示.

圖3　CSEM 正演并行計(jì)算流程圖Fig.3　Flowchart of 3D parallel CSEM forward modeling

3　理論模型計(jì)算

本節(jié)通過對(duì)多個(gè)理論模型的正演數(shù)值結(jié)果的對(duì)比分析來檢驗(yàn)本文所開發(fā)的正演算法對(duì)頻率域CSEM模擬計(jì)算的準(zhǔn)確性及有效性.

3.1一維層狀模型的對(duì)比測(cè)試

為了驗(yàn)證本文所開發(fā)的CSEM三維正演算法的準(zhǔn)確性，首先考慮具有準(zhǔn)解析解的一維層狀模型.圖4a是海洋CSEM中用于模擬高阻油氣層的標(biāo)準(zhǔn)模型(Constable and Weiss，2006).為模擬強(qiáng)電性差異，油氣層電阻率設(shè)為1000 Ωm，厚度為100 m，頂部距海底1000 m；海底沉積層電阻率為1 Ωm，海水層厚度為1000 m，電阻率為0.33 Ωm.圖4b為典型陸地一維層狀模型.電阻率為100 Ωm沉積層中包含一層電阻率為10 Ωm的低阻目標(biāo)層，該目標(biāo)層頂部埋深200 m，厚度為400 m.所有模型的空氣層電阻率均取109Ωm.本文所有計(jì)算均在一個(gè)小型并行機(jī)上完成.該小型并行集群系統(tǒng)配置有8個(gè)計(jì)算節(jié)點(diǎn)，每個(gè)計(jì)算節(jié)點(diǎn)含有2個(gè)4核Intel○RXeon○RE5410型處理器，主頻為2.33 GHz，每個(gè)計(jì)算節(jié)點(diǎn)配置有32G內(nèi)存.操作系統(tǒng)為CentOS 5.5.層狀模型的擬解析解均通過Dipole1D程序(Key，2009)求得.

對(duì)于海底層狀模型(圖4a)，其觀測(cè)系統(tǒng)設(shè)置為：電偶極子位于海底上方100 m，發(fā)射頻率取0.25 Hz和1 Hz，電流強(qiáng)度為1A.接收器位于海底表面，布設(shè)方向與偶極方向一致(inline觀測(cè)方式).模型剖分網(wǎng)格數(shù)為132×52×52.圖5給出了不同頻率下Ex分量和By分量的振幅-收發(fā)距(Magnitude Versus Offset, MVO)曲線與相位-收發(fā)距(Phase Versus Offset, PVO)曲線的MFV三維數(shù)值解與解析解的對(duì)比.從圖中可以看出，不同頻率的電磁場(chǎng)分量的MFV數(shù)值解與解析解均高度一致.圖6則給出了MFV解相對(duì)于解析解的誤差情況.當(dāng)未進(jìn)行場(chǎng)源校正時(shí)，可以看到在發(fā)射偶極附近(收發(fā)距小于1000 m)，電磁場(chǎng)的各個(gè)分量均具有較大的誤差，振幅最大誤差可達(dá)10%左右，相位最大誤差可達(dá)2°，這主要是由于場(chǎng)源的奇異性造成的.需要說明的是，為降低待求解的線性方程組的尺寸，減少直接求解的計(jì)算量，本文并未在場(chǎng)源點(diǎn)附近進(jìn)行精細(xì)的網(wǎng)格加密來降低場(chǎng)源奇異性的影響.當(dāng)采用場(chǎng)源校正技術(shù)后，在場(chǎng)源附近，場(chǎng)源奇異性對(duì)于數(shù)值精度的影響會(huì)迅速降低，振幅最大誤差降為5%以內(nèi)，相位最大誤差降為1.5°以內(nèi).另外，隨著收發(fā)距的增加，場(chǎng)源奇異性的影響迅速降低，誤差也隨之迅速減小，電磁場(chǎng)分量的振幅誤差降為2%以內(nèi)，相位誤差在2°以內(nèi).對(duì)于海洋CSEM來說，近場(chǎng)區(qū)內(nèi)(收發(fā)距小于1000 m)一次場(chǎng)占據(jù)主導(dǎo)地位，一般并不包含海底目標(biāo)層信息(Zhdanov et al., 2014).另外，從數(shù)值計(jì)算結(jié)果對(duì)比來看，場(chǎng)源奇異性對(duì)于遠(yuǎn)場(chǎng)區(qū)的計(jì)算精度影響有限，因此采用總場(chǎng)公式對(duì)于海洋CSEM的正演計(jì)算是完全可行的.

圖4　海洋(a)及陸地(b)一維層狀模型Fig.4　Sketches of (a) marine-based and (b) land-based 1D layered models

圖5　海洋層狀模型CSEM數(shù)值解與擬解析解對(duì)比(a)(c)為電磁場(chǎng)分量的振幅，(b)(d)為對(duì)應(yīng)電磁場(chǎng)的相位.Fig.5　Accuracy comparison between MFV numerical solutions and quasi-analytic solutions for marine 1D layered model(a)(c)are amplitudes of field components,(b)(d)are corresponding phases.

圖6　海洋層狀模型CSEM數(shù)值解相對(duì)于擬解析解的誤差(a)(c)為電磁場(chǎng)分量振幅的相對(duì)誤差，(b)(d)為對(duì)應(yīng)電磁場(chǎng)的相位差異.Fig.6　Errors between MFV numerical solutions and quasi-analytic solutions for marine 1D layered model(a)(c)indicate relative errors in amplitude,(b)(d)are corresponding phase difference.

圖7　層狀模型CSEM數(shù)值解相對(duì)于擬解析解的誤差(a)(c)為電磁場(chǎng)分量的振幅，(b)(d)為對(duì)應(yīng)電磁場(chǎng)的相位.Fig.7　Accuracy comparison between MFV numerical solutions and quasi-analytic solutions for land 1D layered model(a)(c)are amplitudes of field components,(b)(d)are corresponding phases.

對(duì)于陸地層狀模型(圖4b)，觀測(cè)系統(tǒng)采用可控源音頻大地電磁法(CSAMT)(底青云和王若，2008)的觀測(cè)方式.其中發(fā)射源為1000 m的接地導(dǎo)線，發(fā)射頻點(diǎn)數(shù)為14個(gè)，頻率在0.25～4096 Hz之間呈對(duì)數(shù)等間隔分布；測(cè)線平行于發(fā)射源布設(shè).模型剖分網(wǎng)格數(shù)為99×56×40.圖7給出了與場(chǎng)源中心距離為6 km的測(cè)點(diǎn)上所有頻率的Ey分量和Bx分量的MFV數(shù)值解與解析解對(duì)比.從圖中可以看出，所有頻點(diǎn)的數(shù)值解與解析解均非常吻合.圖8則給出了數(shù)值解相對(duì)于解析解的誤差情況.從圖中可以看出，當(dāng)未進(jìn)行場(chǎng)源校正時(shí)，大多數(shù)頻點(diǎn)的電磁場(chǎng)分量的振幅誤差為6%～8%；與之對(duì)比，采用場(chǎng)源校正技術(shù)后，電磁場(chǎng)分量的振幅誤差降為5%以內(nèi).另外，在進(jìn)行場(chǎng)源校正后，電磁場(chǎng)分量的相位誤差略有增加，但仍保持在2°以內(nèi).在野外CSAMT測(cè)量中，由于各種人文噪聲的影響，視電阻率誤差和相位誤差一般會(huì)比較大(Hu et al., 2013)，經(jīng)過場(chǎng)源校正后，采用總場(chǎng)公式來進(jìn)行CSAMT的正演模擬完全能夠滿足實(shí)測(cè)精度要求.

3.2三維模型的對(duì)比測(cè)試

在1D層狀模型測(cè)試的基礎(chǔ)上，本節(jié)對(duì)海洋3D模型進(jìn)行了正演測(cè)試，并與已經(jīng)發(fā)表的采用散射場(chǎng)公式的正演算法(韓波等, 2015b)所得到的結(jié)果進(jìn)行了對(duì)比驗(yàn)證.

圖9為典型海洋三維油氣藏模型.油氣藏模型尺寸為4 km×4 km×0.1 km，電阻率為100 Ωm，其頂部埋深距海底1000 m.均勻海底沉積層電阻率為1 Ωm，海水層厚度為1 km，電阻率為0.33 Ωm.觀測(cè)系統(tǒng)設(shè)置為：發(fā)射偶極子沿x方向，位于海底上方100 m處，坐標(biāo)為(0, 0, -100 m).接收站沿y=0的測(cè)線布設(shè)于海底，測(cè)線方向與偶極方向一致.發(fā)射頻率為1 Hz.

圖10給出了對(duì)于海洋三維模型，采用本文算法所得到的數(shù)值結(jié)果與散射場(chǎng)算法的正演結(jié)果對(duì)比圖.從圖中可以看出，無論是電場(chǎng)分量(圖10a,b)還是磁場(chǎng)分量(圖10c,d)，利用本文算法所得到的正演結(jié)果均與散射場(chǎng)算法的結(jié)果吻合的非常好，盡管在大偏移距處(>5 km)，兩者的結(jié)果有略微的偏差.

3.3并行計(jì)算效率分析

為了檢驗(yàn)本文所開發(fā)CSEM三維正演算法的并行效率，對(duì)陸地1D層狀模型的多頻率正演響應(yīng)進(jìn)行了并行計(jì)算測(cè)試，按對(duì)數(shù)等間隔在0.25～4096 Hz之間取12個(gè)頻率.表1給出了使用不同數(shù)量的計(jì)算節(jié)點(diǎn)時(shí)的運(yùn)行時(shí)間和并行加速比.可以看出，隨著計(jì)算節(jié)點(diǎn)數(shù)的增加，總的計(jì)算時(shí)間迅速降低.由于各個(gè)頻率之間的計(jì)算是完全獨(dú)立的，主進(jìn)程與各個(gè)子進(jìn)程只涉及到非常少的通信(計(jì)算參數(shù)的分配及結(jié)果回收)，從并行加速比來看，并行效率非常接近于理想情況(加速比等于計(jì)算進(jìn)程數(shù)).

程序運(yùn)行模式計(jì)算節(jié)點(diǎn)數(shù)每個(gè)節(jié)點(diǎn)分配的頻點(diǎn)個(gè)數(shù)模型剖分網(wǎng)格大小程序運(yùn)行時(shí)間(s)并行加速比串行模式11299×56×405315.4無并行模式11299×56×405357.20.992699×56×402688.31.984399×56×401372.63.876299×56×40935.95.68

4　結(jié)論

(1) 本文采用直接解法來求解離散化的電場(chǎng)矢量Helmholtz方程，實(shí)現(xiàn)了基于有限體積法的頻率域CSEM三維正演計(jì)算.對(duì)多個(gè)理論模型的正演測(cè)試結(jié)果表明本文的算法對(duì)于典型海洋和陸地CSEM測(cè)量是有效的.

(2) 由于采用直接解法來求解線性方程組，一次矩陣分解的結(jié)果可以用于不同發(fā)射源位置的計(jì)算，特別適合具有多發(fā)射源的CSEM測(cè)量；此外，在矩陣分解中，系數(shù)矩陣的病態(tài)程度對(duì)于矩陣分解結(jié)果影響很小，使得基于矩陣分解的直接解法對(duì)于包含強(qiáng)電性差異的地電模型都可以獲得穩(wěn)定且高精度的解.

圖9　海洋三維油氣藏模型平面

圖10　海洋三維模型數(shù)值結(jié)果對(duì)比(a)(c)電磁場(chǎng)分量振幅； (b)(d)電磁場(chǎng)分量相位.其中方框?yàn)楸疚恼菟玫降慕Y(jié)果，圓圈為韓波等(2015a)計(jì)算的結(jié)果.Fig.10　Comparison of data obtained from our modeling scheme to data resulted from Han et al. (2015a) for 3D marine model(a) Ex amplitudes； (b) Ex phases； (c) By amplitudes； (d) By phases.

(3) 采用了總場(chǎng)公式，場(chǎng)源的模擬采用直接離散化的方式.與基于散射場(chǎng)公式的正演算法相比，場(chǎng)源項(xiàng)的構(gòu)建并不需要多次求解均勻模型或1D模型的一次電磁場(chǎng)，發(fā)射源位置的增加只引起計(jì)算量的少量增加，同樣有利于多發(fā)射源的CSEM測(cè)量.與此同時(shí)，為避免對(duì)場(chǎng)源點(diǎn)附近進(jìn)行網(wǎng)格局部加密所引起的計(jì)算量的增加，采用虛擬場(chǎng)源校正技術(shù)來降低場(chǎng)源奇異性及邊界條件對(duì)數(shù)值精度的影響，取得了較好的效果.

(4) 三維電磁模擬對(duì)于計(jì)算資源要求很高，并且往往花費(fèi)較長(zhǎng)的計(jì)算時(shí)間.因此提高計(jì)算效率對(duì)于快速的電磁三維正反演至關(guān)重要.隨著個(gè)人工作站和小型服務(wù)器的普及，一個(gè)有效的策略是考慮并行計(jì)算.對(duì)于多頻率CSEM測(cè)量來說(如CSAMT)，基于分頻并行的策略能夠有效降低總的計(jì)算時(shí)間.另外，在計(jì)算資源允許的條件下，還可考慮對(duì)矩陣分解算法進(jìn)行并行，從而進(jìn)一步提高并行效率.

致謝感謝加拿大英屬哥倫比亞大學(xué)(UBC)Eldad Haber教授的指導(dǎo)及與GIF組內(nèi)同學(xué)進(jìn)行的有益討論，感謝MUMPS線性運(yùn)算庫(kù)的開發(fā)者們，感謝兩名匿名審稿人對(duì)本文提出的寶貴修改意見.

Amestoy P R, Duff I S, L′Excellent J Y, et al. 2001. A fully asynchronous multifrontal solver using distributed dynamic scheduling.SIAMJournalonMatrixAnalysisandApplications, 23(1): 15-41. Amestoy P R, Guermouche A, L′Excellent J Y, et al. 2006. Hybrid scheduling for the parallel solution of linear systems.ParallelComputing, 32(2): 136-156.

Ansari S, Farquharson C G. 2014. 3D finite-element forward modeling of electromagnetic data using vector and scalar potentials and unstructured grids.Geophysics, 79(4): E149-E165.

Avdeev D, Knizhnik S. 2009. 3D integral equation modeling with a linear dependence on dimensions.Geophysics, 74(5): F89-F94.

Avdeev D B. 2005. Three-dimensional electromagnetic modelling and inversion from theory to application.SurveysinGeophysics, 26(6): 767-799.Avdeev D B, Kuvshinov A V, Pankratov O V, et al. 2002. Three-dimensional induction logging problems, Part I: An integral equation solution and model comparisons.Geophysics, 67(2): 413-426.

Badea E A, Everett M E, Newman G A, et al. 2001. Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials.Geophysics, 66(3): 786-799.

B?rner R U. 2010. Numerical modelling in geo-electromagnetics: advances and challenges.SurveysinGeophysics, 31(2): 225-245.

Chung Y, Son J S, Lee T J, et al. 2014. Three-dimensional modelling of controlled-source electromagnetic surveys using an edge finite-element method with a direct solver.GeophysicalProspecting, 62(6): 1468-1483.

Di Q Y, Wang R. 2008. Controlled Source Audio-Frequency Magneto Tellurics (in Chinese). Beijing: Science Press.

Grayver A V, Streich R, Ritter O. 2013. Three-dimensional parallel distributed inversion of CSEM data using a direct forward solver.GeophysicalJournalInternational, 193(3): 1432-1446.

Haber E, Ascher U M. 2001. Fast finite volume simulation of 3D electromagnetic problems with highly discontinuous coefficients.SIAMJournalonScientificComputing, 22(6): 1943-1961. Haber E, Heldmann S. 2007. An octree multigrid method for quasi-static Maxwell′s equations with highly discontinuous coefficients.JournalofComputationalPhysics, 223(2): 783-796. Haber E, Ruthotto L. 2014. A multiscale finite volume method for Maxwell′s equations at low frequencies.GeophysicalJournalInternational, 199(2): 1268-1277.

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

Han B, Hu X Y, Schultz A, et al. 2015a. 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.

Huang Q H, Li Z H, Wang Y B. 2010. A parallel 3-D staggered grid pseudospectral time domain method for ground-penetrating radar wave simulation.JournalofGeophysicalResearch, 115(B12): B12101, doi: 10.1029/2010JB007711.

Hursán G, Zhdanov M S. 2002. Contraction integral equation method in three-dimensional electromagnetic modeling.RadioScience, 37(6): 1-1-1-13.

Hyman J, Shashkov M, Steinberg S. 1997. The numerical solution of diffusion problems in strongly heterogeneous non-isotropic materials.JournalofComputationalPhysics, 132(1): 130-148.

Hyman J M, Shashkov M. 1997. Adjoint operators for the natural discretizations of the divergence, gradient and curl on logically rectangular grids.AppliedNumericalMathematics, 25(4): 413-442.

Hyman J M, Shashkov M. 1999a. Mimetic discretizations for Maxwell′s equations.JournalofComputationalPhysics, 151(2): 881-909.

Hyman J M, Shashkov M. 1999b. The orthogonal decomposition theorems for mimetic finite difference methods.SIAMJournalonNumericalAnalysis, 36(3): 788-818.

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.

Knizhnerman L, Druskin V, Zaslavsky M. 2009. On optimal convergence rate of the rational Krylov subspace reduction for electromagnetic problems in unbounded domains.SIAMJournalonNumericalAnalysis, 47(2): 953-971. Li X S. 2005. An overview of SuperLU: Algorithms, implementation, and user interface.ACMTransactionsonMathematicalSoftware(TOMS), 31(3): 302-325.

Li Y, Wu X P, Lin P R. 2015. Three-dimensional controlled source electromagnetic finite element simulation using the secondary field for continuous variation of electrical conductivity within each block.ChineseJ.Geophys. (in Chinese), 58(3): 1072-1087, doi: 10.6038/cjg20150331.

Liu L B, Li Z H, Arcone S, et al. 2013. Radar wave scattering loss in a densely packed discrete random medium: Numerical modeling of a box-of-boulders experiment in the Mie regime.JournalofAppliedGeophysics, 99: 68-75, doi: 10.1016/j.jappgeo.2013.08.022.

Mackie R L, Smith J T, Madden T R. 1994. Three-dimensional electromagnetic modeling using finite difference equations: the magnetotelluric example.RadioScience, 29(4): 923-935.

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. 2013. Three dimensional inversion of multisource time domain electromagnetic data.Geophysics, 78(1): E47-E57.

Pardo D, Paszynski M, Collier N, et al. 2012. A survey on direct solvers for Galerkin methods.SEMAJournal, 57(1): 107-134.

Pidlisecky A, Haber E, Knight R. 2007. RESINVM3D: A 3D resistivity inversion package.Geophysics, 72(2): H1-H10.

Pidlisecky A, Knight R. 2008. FW2_5D: A MATLAB 2.5-D electrical resistivity modeling code.Computer&Geosciences, 34(12): 1645-1654.

Plessix R E, Darnet M, Mulder W A. 2007. An approach for 3D multisource, multifrequency CSEM modeling.Geophysics, 72(5): SM177-SM184.

Puzyrev V, Koldan J, de la Puente J, et al. 2013. A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling.GeophysicalJournalInternational, 193(2): 678-693. Qiu Z P, Li Z H, Li D Z, et al. 2013. Non-orthogonal-Grid-based three dimensional modeling of transient electromagnetic field with topography.ChineseJ.Geophys. (in Chinese), 56(12): 4245-4255, doi: 10.6038/cjg20131227.

Schenk O, G?rtner K. 2006. On fast factorization pivoting methods for sparse symmetric indefinite systems.ElectronicTransactionsonNumericalAnalysis, 23(1): 158-179.

Schwarzbach C, B?rner R U, Spitzer K. 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics—a marine CSEM example.GeophysicalJournalInternational, 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-288.

Smith J T. 1996. Conservative modeling of 3-D electromagnetic fields, Part I: Properties and error analysis.Geophysics, 61(5): 1308-1318.

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, Matzander U, et al. 2011. Strategies for land-based controlled-source electromagnetic surveying in high-noise regions.TheLeadingEdge, 30(10): 1174-1181.

Ward S H, Hohmann G W. 1988. Electromagnetic theory for geophysical applications.∥ Electromagnetic Methods in Applied Geophysics: Vol. I, Theory. Tulsa, OK: SEG. Weiss C J, Constable S. 2006. Mapping thin resistors and hydrocarbons with marine EM methods, Part II: Modeling and analysis in 3D.Geophysics, 71(6): G321-G332. Weiss C J, Newman G A. 2002. Electromagnetic induction in a fully 3-D anisotropic earth.Geophysics, 67(4): 1104-1114.

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.

Yang J, Liu Y, Wu X P. 2015. 3D simulation of marine CSEM using vector finite element method on unstructured grids.ChineseJ.Geophys. (in Chinese), 58(8): 2827-2838, doi: 10.6038/cjg20150817.

Yee K S. 1966. Numerical solution of initial boundary value problems involving Maxwell′s equations in isotropic media.IEEETransactionsonAntennasandPropagation, 14(3): 302-307.

Zhdanov M S, Lee S K, Yoshioka K. 2006. Integral equation method for 3D modeling of electromagnetic fields in complex structures with inhomogeneous background conductivity.Geophysics, 71(6): G333-G345.

Zhdanov M S, Masashi E, Cox L H, et al. 2014. Three-dimensional inversion of towed streamer electromagnetic data.GeophysicalProspecting, 62(3): 552-572.

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

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

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

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

李勇, 吳小平, 林品榮. 2015. 基于二次場(chǎng)電導(dǎo)率分塊連續(xù)變化的三維可控源電磁有限元數(shù)值模擬. 地球物理學(xué)報(bào), 58(3): 1072-1087, doi: 10.6038/cjg20150331.

邱稚鵬, 李展輝, 李墩柱等. 2013. 基于非正交網(wǎng)格的帶地形三維瞬變電磁場(chǎng)模擬. 地球物理學(xué)報(bào), 56(12): 4245-4255, doi: 10.6038/cjg20131227.

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

徐志鋒, 吳小平. 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.

楊軍, 劉穎, 吳小平. 2015. 海洋可控源電磁三維非結(jié)構(gòu)矢量有限元數(shù)值模擬. 地球物理學(xué)報(bào), 58(8): 2827-2838, doi: 10.6038/cjg20150817.

(本文編輯張正峰)

3D frequency-domain CSEM forward modeling based on the mimetic finite-volume method

PENG Rong-Hua1,2, HU Xiang-Yun1*, HAN Bo3, CAI Jian-Chao1

1InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China2DepartmentofEarth,OceanandAtmosphericSciences,UniversityofBritishColumbia,VancouverV6T1Z4,Canada3CollegeofMarineGeosciences,OceanUniversityofChina,Qingdao266100,China

Quantitative interpretation of large-scale controlled-source electromagnetic (CSEM) data in the frequency domain requires efficient and stable three-dimensional (3D) forward modeling and inversion codes. In this paper, we present a 3D forward modeling scheme for frequency-domain CSEM surveys based on the mimetic finite volume (MFV) method, which solves the Helmholtz equation for the total electric field.

Controlled-source electromagnetic； Mimetic finite volume method； 3D forward modeling； Source correction technique； Direct solver

10.6038/cjg20161036.

國(guó)家自然科學(xué)基金(41274077和41474055)、國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃項(xiàng)目(2013CB733200)、國(guó)家留學(xué)基金委(201406410020)和湖北省自然科學(xué)基金(2015CFA019)聯(lián)合資助.

彭榮華，男，1988年生，博士研究生，研究方向?yàn)殡姶欧ㄈS正演與反演模擬.E-mail:prhjiajie@163.com

胡祥云，男，1966年生，教授，博士生導(dǎo)師，主要從事電磁法的理論與應(yīng)用研究.E-mail：xyhu@cug.edu.cn

10.6038/cjg20161036

P631

2016-07-06，2016-09-12收修定稿

彭榮華， 胡祥云， 韓波等. 2016. 基于擬態(tài)有限體積法的頻率域可控源三維正演計(jì)算. 地球物理學(xué)報(bào)，59(10):3927-3939,

Peng R H, Hu X Y, Han B, et al. 2016. 3D frequency-domain CSEM forward modeling based on the mimetic finite-volume method.ChineseJ.Geophys. (in Chinese),59(10):3927-3939,doi:10.6038/cjg20161036.