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

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

?

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

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

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

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

可控源電磁法； 有限體積法； 虛擬場源校正技術(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)前電磁勘探的重要趨勢及研究熱點(diǎn).在電磁數(shù)據(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)在于提高精度和計算效率.

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

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

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

2　方法理論

2.1控制方程

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

(1)

(2)

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

(3)

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

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

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

(4)

(5)

其中W和F分別為與電場E和磁感應(yīng)強(qiáng)度B屬于相同Hilbert空間的任意測試函數(shù)(Hyman and Shashkov, 1999a).對(5)式左端項進(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)都能用來對(4)式和(8)式的積分進(jìn)行離散化.本文采用基于Yee網(wǎng)格(Yee, 1966)的交錯網(wǎng)格.圖1展示了交錯網(wǎng)格電磁場采樣方式及電場分量和磁場分量各自的控制體積.其中電場定義在單元棱邊的中心，磁感應(yīng)強(qiáng)度定義在單元面的中心，介質(zhì)電導(dǎo)率和磁導(dǎo)率定義在單元中心.

(10)

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

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

(11)

其中

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

(12)

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

圖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　網(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虛擬場源校正技術(shù)

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

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

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

(13)

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

(14)

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

(15)

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

2.4大型線性方程組的求解

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

2.5多頻率并行計算

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

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

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

3　理論模型計算

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

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

為了驗證本文所開發(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ī)上完成.該小型并行集群系統(tǒng)配置有8個計算節(jié)點(diǎn)，每個計算節(jié)點(diǎn)含有2個4核Intel○RXeon○RE5410型處理器，主頻為2.33 GHz，每個計算節(jié)點(diǎn)配置有32G內(nèi)存.操作系統(tǒng)為CentOS 5.5.層狀模型的擬解析解均通過Dipole1D程序(Key，2009)求得.

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

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

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

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

3.2三維模型的對比測試

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

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

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

3.3并行計算效率分析

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

程序運(yùn)行模式計算節(jié)點(diǎn)數(shù)每個節(jié)點(diǎn)分配的頻點(diǎn)個數(shù)模型剖分網(wǎng)格大小程序運(yùn)行時間(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) 本文采用直接解法來求解離散化的電場矢量Helmholtz方程，實現(xiàn)了基于有限體積法的頻率域CSEM三維正演計算.對多個理論模型的正演測試結(jié)果表明本文的算法對于典型海洋和陸地CSEM測量是有效的.

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

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

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

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

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

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é)報, 58(8): 2812-2826, doi: 10.6038/cjg20150816.

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

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

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

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

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

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

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

彭榮華，男，1988年生，博士研究生，研究方向為電磁法三維正演與反演模擬.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)有限體積法的頻率域可控源三維正演計算. 地球物理學(xué)報，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.