基于擬態有限體積法的頻率域可控源三維正演計算

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

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

?

基于擬態有限體積法的頻率域可控源三維正演計算

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

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

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

可控源電磁法； 有限體積法； 虛擬場源校正技術； 三維正演； 直接分解法

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)由于具有勘探深度大，分辨力高且野外抗干擾能力強等優點，經過30多年快速發展，目前已被廣泛應用于包括陸地、海洋及航空等領域內的油氣、礦產資源勘探和近地表地質勘查.隨著電磁儀器、數據采集技術以及解釋方法的快速發展，可控源電磁法在復雜地質環境下的三維勘探成為當前電磁勘探的重要趨勢及研究熱點.在電磁數據的處理解釋中，反演是必不可少的步驟，而正演計算是反演的核心.因此，電磁數據的定量解釋首先需要發展有效的正演算法.從數值模擬的角度來說，目前常用于三維電磁正演的數值方法主要包括：積分方程法(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)也被廣泛用于三維電磁正演問題的求解.有關電磁法三維正演方面進展的綜述可以詳見(Avdeev, 2005；B?rner, 2010).總體來說，目前三維正演算法研究的重點在于提高精度和計算效率.

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

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

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

2　方法理論

2.1控制方程

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

(1)

(2)

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

(3)

2.2擬態有限體積法

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

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

(4)

(5)

其中W和F分別為與電場E和磁感應強度B屬于相同Hilbert空間的任意測試函數(Hyman and Shashkov, 1999a).對(5)式左端項進行分部積分，可得公式為

(6)

(7)

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

(8)

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

(10)

由于F和W為任意的網格測試函數，將(9)式代入(10)式并進行簡化可得(3)式的離散形式為

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

(11)

其中

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

(12)

其中系數矩陣A為稀疏、正定、對稱復矩陣.圖2給出了網格數為5×5×5的模型的MFV離散所得到的系數矩陣A的稀疏結構.

圖1　電磁場分量交錯采樣方式及其控制體積(a) 電場x分量采樣位置及其控制體積； (b) 磁場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　網格數為5×5×5的模型的MFV離散線性方程系數矩陣的稀疏結構Fig.2　Sparse structure of coefficient matrix A resulted from MFV discretization on a 5×5×5 staggered grid

2.3虛擬場源校正技術

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

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

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

(13)

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

(14)

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

(15)

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

2.4大型線性方程組的求解

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

2.5多頻率并行計算

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

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

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

3　理論模型計算

本節通過對多個理論模型的正演數值結果的對比分析來檢驗本文所開發的正演算法對頻率域CSEM模擬計算的準確性及有效性.

3.1一維層狀模型的對比測試

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

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

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

圖5　海洋層狀模型CSEM數值解與擬解析解對比(a)(c)為電磁場分量的振幅，(b)(d)為對應電磁場的相位.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數值解相對于擬解析解的誤差(a)(c)為電磁場分量振幅的相對誤差，(b)(d)為對應電磁場的相位差異.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數值解相對于擬解析解的誤差(a)(c)為電磁場分量的振幅，(b)(d)為對應電磁場的相位.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.

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

3.2三維模型的對比測試

在1D層狀模型測試的基礎上，本節對海洋3D模型進行了正演測試，并與已經發表的采用散射場公式的正演算法(韓波等, 2015b)所得到的結果進行了對比驗證.

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

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

3.3并行計算效率分析

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

程序運行模式計算節點數每個節點分配的頻點個數模型剖分網格大小程序運行時間(s)并行加速比串行模式11299×56×405315.4無并行模式11299×56×405357.20.992699×56×402688.31.984399×56×401372.63.876299×56×40935.95.68

4　結論

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

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

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

圖10　海洋三維模型數值結果對比(a)(c)電磁場分量振幅； (b)(d)電磁場分量相位.其中方框為本文正演所得到的結果，圓圈為韓波等(2015a)計算的結果.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) 采用了總場公式，場源的模擬采用直接離散化的方式.與基于散射場公式的正演算法相比，場源項的構建并不需要多次求解均勻模型或1D模型的一次電磁場，發射源位置的增加只引起計算量的少量增加，同樣有利于多發射源的CSEM測量.與此同時，為避免對場源點附近進行網格局部加密所引起的計算量的增加，采用虛擬場源校正技術來降低場源奇異性及邊界條件對數值精度的影響，取得了較好的效果.

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

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

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.

附中文參考文獻

底青云, 王若. 2008. 可控源音頻大地電磁數據正反演及方法應用. 北京: 科學出版社.

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

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

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

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

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

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

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

楊軍, 劉穎, 吳小平. 2015. 海洋可控源電磁三維非結構矢量有限元數值模擬. 地球物理學報, 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.

國家自然科學基金(41274077和41474055)、國家重點基礎研究發展計劃項目(2013CB733200)、國家留學基金委(201406410020)和湖北省自然科學基金(2015CFA019)聯合資助.

彭榮華，男，1988年生，博士研究生，研究方向為電磁法三維正演與反演模擬.E-mail:prhjiajie@163.com

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

10.6038/cjg20161036

P631

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

彭榮華， 胡祥云， 韓波等. 2016. 基于擬態有限體積法的頻率域可控源三維正演計算. 地球物理學報，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.