基于Lorenz規(guī)范條件下磁矢勢(shì)和標(biāo)勢(shì)耦合方程的頻率域電磁法三維正演

陳輝,殷長春,鄧居智

1 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春　130026 2 東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南昌　330013

?

基于Lorenz規(guī)范條件下磁矢勢(shì)和標(biāo)勢(shì)耦合方程的頻率域電磁法三維正演

陳輝1,2,殷長春1*,鄧居智2

1 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春130026 2 東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南昌330013

為了克服空氣層和地表耦合以及避免一次場計(jì)算，開發(fā)適合不同類型場源、不同應(yīng)用范圍的頻率域三維正演模擬統(tǒng)一平臺(tái)，本文從麥克斯韋基本方程出發(fā)，推導(dǎo)基于Lorenz規(guī)范條件的磁矢勢(shì)和標(biāo)勢(shì)耦合方程；通過將不同類型場源分解成一系列短導(dǎo)線(電性)源組合，采用交錯(cuò)網(wǎng)格采樣和有限體積技術(shù)對(duì)方程進(jìn)行離散得到對(duì)稱大型稀疏線性方程組，并采用Jacobi迭代預(yù)處理QMR(Quasi-Minimum-Residual，擬最小殘差)算法進(jìn)行求解，我們成功實(shí)現(xiàn)不同類型場源、不同應(yīng)用范圍的頻率域電磁法三維正演模擬.通過層狀模型下大地電磁法以及有限長接地導(dǎo)線和大回線磁性源激發(fā)下的電磁場響應(yīng)模擬，并與一維解析解對(duì)比驗(yàn)證算法的有效性.進(jìn)而，我們利用該算法平臺(tái)的模擬結(jié)果對(duì)典型地電模型在不同場源激發(fā)下頻率域電磁法響應(yīng)特征進(jìn)行對(duì)比分析.本文算法研究及實(shí)現(xiàn)為建立頻率域電磁法三維正反演統(tǒng)一框架打下基礎(chǔ).

頻率域電磁法；有限體積法；正演模擬；Lorenz規(guī)范；耦合勢(shì)方程

1　引言

頻率域電磁法(Frequency-Domain Electromagnetic,FDEM)通過觀測和分析不同頻率的人工場源或天然場源激發(fā)地下介質(zhì)產(chǎn)生的二次場或總電磁場分布規(guī)律以探明地下電性分布特征.該方法種類繁多，按照?qǐng)鲈葱再|(zhì)(Nabighian，1988)可分為平面波場的大地電磁法(MT)、音頻大地電磁法(AMT)、甚低頻電磁法(VLF)；接地導(dǎo)線激發(fā)的可控源電磁法(包括地面可控源音頻大地電磁CSAMT和海洋可控源電磁法MCSEM)；不接地回線激法的電磁剖面法(包括地面電磁和航空電磁法AEM)等.這些方法在深部地球結(jié)構(gòu)探測 (Selway,2014)、地?zé)峥碧?Muoz,2014)、礦產(chǎn)勘查(Smith,2014)、油氣勘探(Strack,2014)以及壞境與工程勘探(Sheard et al.,2005)等各領(lǐng)域起到關(guān)鍵作用.

頻率域電磁法三維正、反演已成為目前電磁勘探研究的熱點(diǎn)和難點(diǎn)問題.其中，MT三維正反演技術(shù)已基本成熟，并開始走向?qū)嵱没?Siripunvaraporn，2012；Miensopust et al.,2013).Avdeev(2005)，B?rner(2010)和Newman(2014)對(duì)電磁法三維正反演技術(shù)發(fā)展和挑戰(zhàn)進(jìn)行綜述，指出在正演算法中積分方程法(Wannamaker,1991;Song et al.,2002;Farquharson et al.,2006;Zhdanov et al.,2006;Avdeev and Knizhnik,2009;Bakr and Mannseth,2009)通常利用Green函數(shù)求解關(guān)于二次電場的Fredholm積分式，方法已基本成熟，但對(duì)復(fù)雜模型的正演精度不高；有限單元法(Yoshimura and Oshiman,2002;Mukherjee and Everett,2011;Pardo et al.,2011;Yahya et al.,2012;Schwarzbach and Haber,2013；Cai et al.,2014)通常利用變分原理或加權(quán)余量法求解關(guān)于電場、磁場或耦合勢(shì)的微分方程，能夠模擬復(fù)雜地形和地下結(jié)構(gòu)，但計(jì)算量巨大、求解速度慢；有限差分或有限體積法既能模擬復(fù)雜模型，又具有較快計(jì)算速度，已成為三維電磁正反演中主流算法.Mackie等(1993)使用有限差分法實(shí)現(xiàn)了三維大地電磁正演模擬并將計(jì)算結(jié)果與積分方程法進(jìn)行了對(duì)比，驗(yàn)證了方法的正確性；Smith(1996a,1996b)對(duì)交錯(cuò)采樣有限差分法的原理、誤差分析以及加快計(jì)算速度的雙共軛梯度法進(jìn)行了詳細(xì)闡述；譚捍東等(2003),Tan等(2006)系統(tǒng)論述了消去電場分量得到關(guān)于磁場的大地電磁三維交錯(cuò)采樣有限差分?jǐn)?shù)值模擬算法，并對(duì)實(shí)現(xiàn)過程中交錯(cuò)網(wǎng)格剖分、積分公式離散化、邊界條件、方程組求解、三維張量阻抗的計(jì)算等內(nèi)容做出研究，并實(shí)現(xiàn)了并行計(jì)算；陳輝等(2011)在此基礎(chǔ)上提出利用磁場散度實(shí)校正方法(RRCM)提高三維正演精度和加快正演速度；沈金松(2003)利用交錯(cuò)網(wǎng)格有限差分法消去磁場分量得到關(guān)于電場的線性方程組，實(shí)現(xiàn)了三維頻率域電磁響應(yīng)正演模擬，并且在文中提出在迭代過程中施加了散度校正的方法；Weiss和Constable(2006)利用交錯(cuò)網(wǎng)格有限差分求解二次電場的雙旋度方程實(shí)現(xiàn)海洋電磁法三維數(shù)值模擬；鄧居智等(2011)利用交錯(cuò)網(wǎng)格有限差分求解二次磁場的雙旋度方程實(shí)現(xiàn)可控源音頻大地電磁模擬；殷長春等(2014)利用交錯(cuò)網(wǎng)格有限差分求解二次電場的雙旋度方程實(shí)現(xiàn)海洋電磁法各向異性三維正演模擬；Haber等(2000)針對(duì)基于電或磁場的雙旋度方程在空氣電阻率無窮大引起非橢圓方程和地面強(qiáng)烈耦合作用造成求解不準(zhǔn)問題，提出利用有限體積法求解Coulomb規(guī)范條件下矢勢(shì)和標(biāo)勢(shì)耦合方程，以實(shí)現(xiàn)頻率域電磁法三維正演，其離散線性方程組是正定的，但非對(duì)稱.Hou等(2006)和張燁等(2012)利用有限體積法求解Coulomb規(guī)范條件下二次矢勢(shì)和標(biāo)勢(shì)耦合方程，實(shí)現(xiàn)感應(yīng)測井三維各向異性模擬；周建美等(2014)利用該方法實(shí)現(xiàn)海洋電磁法三維各向異性正演模擬.綜上所述，為了解決人工場源源項(xiàng)處理和場源畸變問題，通常采用數(shù)值算法求解二次電場和磁場雙旋度方程；考慮到空氣層和地下介質(zhì)耦合問題，人們近年提出求解Coulomb規(guī)范條件下二次場標(biāo)勢(shì)和矢勢(shì)耦合方程，但該算法具有網(wǎng)格節(jié)點(diǎn)數(shù)巨大，求解一次場非常耗時(shí)及Coulomb規(guī)范條件產(chǎn)生的耦合方程非對(duì)稱等缺點(diǎn).

本文從麥克斯韋方程組出發(fā)，在Lorenz規(guī)范條件下建立磁矢勢(shì)和標(biāo)勢(shì)耦合天然對(duì)稱方程，采用交錯(cuò)采樣網(wǎng)格和有限體積法離散磁矢勢(shì)和標(biāo)勢(shì)方程，并將不同場源類型分解成一系列短導(dǎo)線(電性)源組合，對(duì)離散線性方程組采用Jacobi迭代的預(yù)處理和擬最小殘差法(QMR)求解，從而實(shí)現(xiàn)任意場源激發(fā)下頻率域三維電磁正演.最后，我們通過與各種不同類型一維理論解進(jìn)行對(duì)比驗(yàn)證了該算法的有效性，并進(jìn)一步對(duì)比分析不同場源激發(fā)下典型地電模型的電磁響應(yīng)特征.

2　磁矢勢(shì)和標(biāo)勢(shì)耦合方程及邊界條件

在電磁勘探中，對(duì)大多數(shù)地下巖石我們可以忽略位移電流.假定時(shí)間變化因子為eiω t，則準(zhǔn)靜態(tài)條件下麥克斯韋方程式可寫為

(1a)

(1b)

(1c)

(1d)

式中，E為電場，H為磁場，Je為電流密度，μ為磁導(dǎo)率，ε為介電常數(shù)，σ為電導(dǎo)率.

(2)

假設(shè)介質(zhì)磁導(dǎo)率μ和介電常數(shù)ε均為常數(shù)，將式(2)代入式(1b)可得

(3)

為了保證磁矢勢(shì)A和標(biāo)勢(shì)φ唯一及方程的對(duì)稱性，我們引入洛倫茲(Lorenz)規(guī)范條件：

(4)

其中c=iωμ.將式(4)代入式(3)可得

(5)

對(duì)式(3)兩邊取散度,并代入Lorenz規(guī)范條件(4)可得

(7)

可以看出，該耦合方程具有天然對(duì)稱性，對(duì)離散后形成的線性方程組求解精度和穩(wěn)定性有一定提升.由于地下介質(zhì)為有損介質(zhì)，電磁場隨傳播距離增加呈指數(shù)衰減，在正演模擬中可選擇足夠大的求解區(qū)域Ω，在求解區(qū)域Ω外的電磁場值非常小.因此可以采用簡單的截?cái)噙吔鐥l件.與其對(duì)應(yīng)的Lorenz規(guī)范條件下磁矢勢(shì)和標(biāo)勢(shì)的邊界條件可表示為

(8)

3　頻率域三維有限體積電磁正演模擬

為了求解磁矢勢(shì)和標(biāo)勢(shì)耦合方程(7)，我們將系統(tǒng)闡述有限體積法離散過程、場源處理及線性方程組合成求解等.

3.1區(qū)域離散

首先，將求解區(qū)域Ω采用一系列平行于x、y、z坐標(biāo)軸的三組平面剖分成若干個(gè)小的六面體網(wǎng)格單元.假設(shè)沿x軸方向被剖分成Nx段，含有l(wèi)個(gè)節(jié)點(diǎn)(l=Nx+1)；節(jié)點(diǎn)編號(hào)沿x軸方向遞增i=1,2,…,l，網(wǎng)格間距為dxi(i=1,…,Nx)；沿y軸方向被剖分成Ny段，含有m個(gè)節(jié)點(diǎn)(m=Ny+1)，節(jié)點(diǎn)編號(hào)沿y軸方向遞增，j=1,2,…,m，網(wǎng)格間距為dyj(j=1,…,Ny)；沿z軸方向被剖分成Nz段，含有n個(gè)節(jié)點(diǎn)(n=Nz+1)，節(jié)點(diǎn)編號(hào)沿z軸方向遞增，k=1,2,…,n，網(wǎng)格間距為dzk(k=1,…,Nz);共計(jì)有Nx×Ny×Nz個(gè)長方體剖分單元.

圖1　磁矢勢(shì)和標(biāo)勢(shì)交錯(cuò)采樣網(wǎng)格Fig.1　Staggered grids of magnetic vectors and scalar potentials

3.2方程離散

為實(shí)現(xiàn)有限體積法離散，需要定義磁矢勢(shì)三個(gè)分量和標(biāo)勢(shì)體積單元大小.我們以Ax分量為例.參見圖2a，體積單元沿x方向邊長為dxi，沿y和z方向的邊長分別為(dyj-1+dyj)/2和(dzk-1+dzk)/2.因此，該單元體積為

Vi+1/2,j,k=dxi(dzk-1+dzk)(dyj-1+dyj)/4.

(9)

按照等效電阻率原理，該體積單元電導(dǎo)率可以表示為

(10)

圖2　磁矢勢(shì)和標(biāo)勢(shì)體積單元示意圖Fig.2　Sketch showing cells of magnetic vectors and scalar potentials

同理可得Ay,Az分量體積單元表達(dá)式.對(duì)于標(biāo)勢(shì)φ體積單元(參見圖2b)，沿x，y和z方向的邊長分別為(dxi-1+dxi)/2，(dyj-1+dyj)/2，(dzk-1+dzk)/2.由此，單元體積為

(11)

而其體積單元等效電導(dǎo)率為

(14)

(15)

(16)

對(duì)于磁矢勢(shì)和標(biāo)勢(shì)耦合方程組(2)第二式，在某個(gè)標(biāo)勢(shì)φ體積單元上積分可得

(18)

仿照前面離散過程，(18)式中左端四項(xiàng)的離散公式為

(19)

(20)

(21)

(22)

3.3場源設(shè)置

3.4離散線性方程組合成及求解

將(14)—(17)式代入方程組(7)可以建立磁矢勢(shì)Ax分量和標(biāo)勢(shì)φ的離散表達(dá)式(見圖4a)，同理將(19)—(22)式代入方程組(7)可建立標(biāo)勢(shì)φ和磁矢勢(shì)Ax、Ay、Az分量的離散表達(dá)式(見圖4b).將方程組(7)離散后Ax、Ay、Az、φ四個(gè)等式，并且結(jié)合邊界條件(8)式可以合成總體線性方程組：

(23)

式中Cx，Cy，Cz，Cφ及Dx，Dy，Dz為系數(shù)項(xiàng)，而rx，ry，rz，rφ為源項(xiàng)，待求未知數(shù)個(gè)數(shù)為

N=(nx-1)×(ny-2)×(nz-2)

+(nx-2)×(ny-1)×(nz-2)

+(nx-2)×(ny-2)×(nz-1)

+(nx-2)×(ny-2)×(nz-2).

(24)

對(duì)于線性方程組求解目前主要采用直接求解和迭代求解.迭代求解通常在Krylov子空間內(nèi)進(jìn)行.首先采用不完全LU分解、不完全Cholesky分解等預(yù)處理方法降低系數(shù)矩陣的條件數(shù)，然后采用共軛梯度法(CG),雙共軛梯度法(BiCG)，廣義最小殘量法(GMRES),擬最小殘差法(QMR)等迭代方法進(jìn)行求解.該算法需要的內(nèi)存小，計(jì)算速度快，適合求解場源較少的正反演問題.直接求解方法一般直接對(duì)系數(shù)矩陣進(jìn)行LU分解，然后利用各種成熟的求解庫MUMPS、PARDISO等進(jìn)行求解.該方法所需內(nèi)存大、計(jì)算速度慢，但特別適合多源問題.本文離散線性方程組(23)的系數(shù)矩陣為大型對(duì)稱稀疏矩陣，由于僅考慮單一發(fā)射源問題，我們選用Jacobi迭代的QMR算法進(jìn)行求解.

圖3　頻率域電磁法不同類型場源類型Fig.3　Transmitting sources of different types for FDEM

圖4　標(biāo)勢(shì)和磁矢勢(shì)離散關(guān)系圖Fig.4　Template for discretizing magnetic vectors and scalar potentials

4　數(shù)值模擬與結(jié)果分析

為了驗(yàn)證本文的頻率域三維電磁正演算法的可行性和有效性，我們利用本文算法結(jié)果和一維正演結(jié)果進(jìn)行對(duì)比，研究不同場源激發(fā)條件下頻率域視電阻率和相位響應(yīng)特征.

4.1算法驗(yàn)證

選取水平層狀三層H模型(如圖5a).第一層厚度500 m，電阻率100 Ωm；第二層厚度1000 m，電阻率10 Ωm；第三層電阻率為1000 Ωm.求解區(qū)域大小為10000 m×10000 m×10000 m，網(wǎng)格剖分為100×100×100，分別向外延拓16層，延拓指數(shù)為1.5.

首先進(jìn)行平面波場X極化和Y極化模式下的電磁場三維模擬，然后計(jì)算大地電磁阻抗視電阻率和相位.計(jì)算頻率為1～1000 Hz.圖5b為水平層狀介質(zhì)三維和一維正演結(jié)果響應(yīng)對(duì)比.由圖可見三維模擬的視電阻率和相位與一維結(jié)果完全吻合，視電阻率相對(duì)誤差最大0.3%，相位最大相對(duì)誤差2.9%，說明本文算法對(duì)大地電磁法三維模擬具有較高精度.

為驗(yàn)證本文算法的廣譜有效性，我們進(jìn)一步對(duì)有限長接地線源的CSAMT進(jìn)行三維正演模擬.網(wǎng)格參數(shù)設(shè)置和上面模型一致.假設(shè)接地長導(dǎo)線發(fā)射源沿x方向置于地面，長度為2000 m，中心點(diǎn)相對(duì)投影(0,0,0)，發(fā)射頻率為10 Hz.接收點(diǎn)設(shè)在位于發(fā)射源中垂線的正向y軸上，點(diǎn)距間隔為100 m.圖6為層狀介質(zhì)CSAMT三維和一維模擬結(jié)果對(duì)比.從圖可以看出，Ex分量的實(shí)部和虛部三維模擬結(jié)果幾乎與一維結(jié)果重合，最大相對(duì)誤差不超過3%；Hz分量實(shí)部和虛部也基本與一維正演結(jié)果吻合，除了Hz分量實(shí)部受場源的畸變影響在發(fā)射源附近兩個(gè)測點(diǎn)相對(duì)誤差達(dá)到10%，其他測點(diǎn)相對(duì)誤差均小于5%.考慮到CSAMT通常在遠(yuǎn)區(qū)或過渡區(qū)進(jìn)行觀測，在場源附近發(fā)生的畸變并不會(huì)影響數(shù)據(jù)解釋精度.

圖5　三層水平層狀模型及大地電磁響應(yīng)結(jié)果對(duì)比Fig.5　Comparison of 1D and 3D MT responses for a 3-layer model

圖6　接地長導(dǎo)線水平層狀模型CSAMT一維和三維電磁響應(yīng)對(duì)比Fig.6　Comparison of 1D and 3D CSAMT responses for the 3-layer model of Fig.4

圖7　不接地回線三層水平層狀模型一維和三維電磁響應(yīng)對(duì)比Fig.7　Comparison of 1D and 3D EM responses for a wire-loop transmitter over the 3-layer model in Fig.4

最后對(duì)大回線源激發(fā)條件下的三維電磁場進(jìn)行正演模擬.網(wǎng)格參數(shù)設(shè)置和上面模型一致.設(shè)置發(fā)射源為3000 m×3000 m方形回線，置于地面，中心點(diǎn)相對(duì)投影(0,0,0)，發(fā)射頻率為10 Hz.接收點(diǎn)位于發(fā)射源中垂線的正向y軸上，點(diǎn)距100 m.圖6為層狀介質(zhì)大回線激發(fā)下的電磁場三維和一維模擬結(jié)果對(duì)比.從圖中可以看出，Ex分量的實(shí)部和虛部三維模擬和一維模擬結(jié)果重合，最大相對(duì)誤差不超過3%；Hz分量實(shí)部和虛部也基本與一維正演結(jié)果吻合，除了Hz分量實(shí)部受場源的畸變影響在靠近發(fā)射源1.5 km附近相對(duì)誤差可達(dá)10%，但其余各測點(diǎn)相對(duì)誤差均小于5%.因此，本文算法對(duì)大回線激發(fā)下的電磁場三維模擬具有較高精度.

4.2典型地電模型

通過水平三層模型不同類型場源的數(shù)值模擬結(jié)果與解析解的對(duì)比，證明了本文算法適合各類場源的頻域電磁法三維正演，并且具有較高精度.下面，我們研究不同類型場源復(fù)雜模型的頻率域電磁響應(yīng)特征，特別是進(jìn)行視電阻率和相位響應(yīng)結(jié)果的對(duì)比.如圖8，我們?cè)O(shè)計(jì)一個(gè)低阻異常體模型大小為800 m×800 m×500 m，頂部埋深為500 m，異常體電阻率為10 Ωm，圍巖電阻率為100 Ωm.將10000 m×10000 m×10000 m計(jì)算區(qū)域剖分為100×100×100個(gè)單元，沿x、y、z三個(gè)方向的擴(kuò)展層數(shù)為16，空氣擴(kuò)展層為12.將異常體剖分成8×8×5個(gè)單元,計(jì)算頻率為1～1000 Hz.人工源頻率域電磁法采用赤道觀測裝置，有線長接地導(dǎo)線長度為2000 m，磁偶源為600 m×600 m回線，收發(fā)距都為6 km.

圖8　低阻異常模型Fig.8　Low-resistivity anomaly model in a homogeneous half-space

圖9為不同場源激發(fā)下低阻模型中心點(diǎn)視電阻率和相位結(jié)果對(duì)比.從圖可以看出，在100 Hz以上四種不同類型場源(XY模式、YX模式、接地長導(dǎo)線、大回線磁偶源)的高頻段視電阻率和相位曲線基本重合；隨著頻率降低XY模式、YX模式大地電磁視電阻率響應(yīng)受三維異常體的影響出現(xiàn)分離，而相位曲線基本重合；其中XY模式受三維異常體影響要大于YX模式.對(duì)于接地長導(dǎo)線、大回線源，隨頻率降低觀測場由遠(yuǎn)區(qū)逐漸進(jìn)入近區(qū)，電阻率呈現(xiàn)快速上升，相位快速下降趨勢(shì).其中大回線磁偶源進(jìn)入近區(qū)頻率要高于接地長導(dǎo)線源，而且幅度大于長導(dǎo)線源.因此，實(shí)際電磁勘探工作中，為便于在遠(yuǎn)區(qū)觀測通常采用接地長導(dǎo)線源作為發(fā)射源.

圖10為低阻異常體模型不同類型場源激發(fā)條件下視電阻率和相位擬斷面圖.由圖可以看出，當(dāng)頻率高于100 Hz時(shí)，即處于遠(yuǎn)區(qū)條件下，CSAMT、回線源赤道裝置的視電阻率響應(yīng)基本相同，在異常體上方近地表位置有一個(gè)不太明顯的小范圍高阻區(qū)，下方為低阻異常區(qū)；低阻異常區(qū)的分布范圍與模型中低阻異常體的實(shí)際寬度基本相同；相位特征也具有類似特征，呈現(xiàn)高值異常.隨著頻率降低，XY模式視電阻率呈現(xiàn)低阻視電阻率異常和高值相位異常；YX模式呈現(xiàn)低阻加兩個(gè)對(duì)稱高阻視電阻率異常和高值相位外加兩個(gè)對(duì)稱低值異常；然而接地長導(dǎo)線和大回線源的赤道裝置受近場影響，呈現(xiàn)高阻假異常和低值相位假異常.其中大回線源進(jìn)入近場頻率要高于接地長導(dǎo)線，而且異常響應(yīng)范圍和幅值都要小.因此開發(fā)不同場源的頻率域電磁三維反演算法可克服三維異常體和場源引起的虛假異常，有助于提高數(shù)據(jù)解釋質(zhì)量.

圖9　不同場源激發(fā)下低阻模型中心點(diǎn)視電阻率和相位結(jié)果對(duì)比Fig.9　Comparison of apparent resistivities and phases at central point in low-resistivity model for different transmitting sources

圖10　不同類型場源低阻異常模型的視電阻率和相位響應(yīng)擬斷面圖橫坐標(biāo)為算術(shù)坐標(biāo)，縱坐標(biāo)頻率為對(duì)數(shù)坐標(biāo).Fig.10　Pseudo-cross section of apparent resistivity and phase for the model in Fig.8 for different transmitting sources

5　結(jié)論

本文從麥克斯韋方程組出發(fā)，在Lorenz規(guī)范條件下建立磁矢勢(shì)和標(biāo)勢(shì)耦合對(duì)稱方程.采用交錯(cuò)網(wǎng)格和有限體積法離散磁矢勢(shì)和標(biāo)勢(shì)耦合方程，并將不同場源類型分解成一系列短導(dǎo)線電性源的組合；對(duì)離散線性方程組采用Jacobi迭代的預(yù)處理擬最小殘差法(QMR)進(jìn)行求解，我們成功實(shí)現(xiàn)任意場源下頻率域三維電磁正演模擬.通過對(duì)水平層狀模型不同類型場源的三維電磁模擬可以發(fā)現(xiàn)，本文開發(fā)的算法能夠?qū)崿F(xiàn)不同類型場源各種不同應(yīng)用范圍頻率域三維電磁正演模擬，避免耗時(shí)的一次場計(jì)算，計(jì)算精度較高.由于采用特殊的源處理技術(shù)，使得頻率域三維電磁正反演統(tǒng)一數(shù)據(jù)處理和解釋平臺(tái)成為可能.對(duì)不同場源激發(fā)下典型地電模型模擬發(fā)現(xiàn)，不同類型場源的視電阻率和相位特征和規(guī)律不同，在實(shí)際勘探中需要考慮場源的選擇和設(shè)置，同時(shí)在電磁數(shù)據(jù)解釋中應(yīng)考慮場源的畸變效應(yīng).

References

Avdeev D B.2005.Three-dimensional electromagnetic modelling and inversion from theory to application.Surveys in Geophysics,26(6):767-799.Avdeev D,Knizhnik S.2009.3D integral equation modeling with a linear dependence on dimensions.Geophysics,74(5):F89-F94.

Bakr S A,Mannseth T.2009.Feasibility of simplified integral equation modeling of low-frequency marine CSEM with a resistive target.Geophysics,74(5):F107-F117.

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

Cai H Z,Xiong B,Han M R,et al.2014.3D controlled-source electromagnetic modeling in anisotropic medium using edge-based finite element method.Computers &Geosciences,73:164-176.

Chen H,Deng J Z,Tan H D,et al.2011.Study on divergence correction method in three-dimensional magnetotelluric modeling with staggered-grid finite difference method.Chinese Journal of Geophysics (in Chinese),54(6):1649-1659,doi:10.3969/j.issn.0001-5733.2011.06.025.

Deng J Z,Tan H D,Chen H,et al.2011.CSAMT 3D modeling using staggered-grid finite difference method.Progress in Geophysics (in Chinese),26(6):2026-2032,doi:10.3969/j.issn.1004-2903.2011.06.017.Farquharson C G,Duckworth K,Oldenburg D W.2006.Comparison of integral equation and physical scale modeling of the electromagnetic responses of models with large conductivity contrasts.Geophysics,71(4):G169-G177.Haber E,Ascher U M,Aruliah D A,et al.2000.Fast simulation of 3D electromagnetic problems using potentials.Journal of Computational Physics,163(1):150-171.

Hou J S,Mallan R K,Torres-Verdín C.2006.Finite-difference simulation of borehole EM measurements in 3D anisotropic media using coupled scalar-vector potentials.Geophysics,71(5):G225-G233.Mackie R L,Madden T R,Wannamaker P E.1993.Three-dimensional magnetotelluric modeling using difference equations—Theory and comparisons to integral equation solutions.Geophysics,58(2):215-226.Miensopust M P,Queralt P,Jones A G.2013.Magnetotelluric 3-D inversion—a review of two successful workshops on forward and inversion code testing and comparison.Geophysical Journal International,193(3):1216-1238.Mukherjee S,Everett M.2011.3D controlled-source electromagnetic edge-based finite element modeling of conductive and permeable heterogeneities.Geophysics,76(4):F215-F226.Muoz G.2014.Exploring for geothermal resources with electromagnetic methods.Surveys in Geophysics,35(1):101-122.Nabighian M N.1988.Electromagnetic Methods in Applied Geophysics.Tulsa:Society of Exploration Geophysicists.Newman G A.2014.A review of high-performance computational strategies for modeling and imaging of electromagnetic induction data.Surveys in Geophysics,35(1):85-100.

Pardo D,Nam M J,Torres-Verdín C,et al.2011.Simulation of marine controlled source electromagnetic measurements using a parallel Fourier hp-finite element method.Computational Geosciences,15(1):53-67.Schwarzbach C,Haber E.2013.Finite element based inversion for time-harmonic electromagnetic problems.Geophysical Journal International,193(2):615-634.

Selway K.2014.On the causes of electrical conductivity anomalies in tectonically stable lithosphere.Surveys in Geophysics,35(1):219-257.

Sheard S N,Ritchie T J,Christopherson K R,et al.2005.Mining,environmental,petroleum,and engineering industry applications of electromagnetic techniques in geophysics.Surveys in Geophysics,26(5):653-669.

Shen J S.2003.Modelling of 3-D electromagnetic responses in frequency domain by using staggered-grid finite difference method.Chinese Journal of Geophysics (in Chinese),46(2):281-288.Siripunvaraporn W.2012.Three-dimensional magnetotelluric inversion:an introductory guide for developers and users.Surveys in Geophysics,33(1):5-27.

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

Smith J T.1996b.Conservative modeling of 3-D electromagnetic fields,Part II:Biconjugate gradient solution and an accelerator.Geophysics,61(5):1319-1324.Smith R.2014.Electromagnetic Induction Methods in Mining Geophysics from 2008 to 2012.Surveys in Geophysics,35(1):123-156.Song Y,Kim H J,Lee K H.2002.An integral equation representation of wide-band electromagnetic scattering by thin sheets.Geophysics,67(3):746-754.Strack K M.2014.Future directions of electromagnetic methods for hydrocarbon applications.Surveys in Geophysics,35(1):157-177.

Tan H D,Tong T,Lin C H.2006.The parallel 3D magnetotelluric forward modeling algorithm.Applied Geophysics,3(4):197-202.

Tan H D,Yu Q F,Booker J,et al.2003.Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method.Chinese Journal of Geophysics (in Chinese),46(5):705-711.Wannamaker P E.1991.Advances in three-dimensional magnetotelluric modeling using integral equations.Geophysics,56(11):1716-1728.

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.

Yahya N,Akhtar M N,Nasir N,et al.2012.Guided and direct wave evaluation of controlled source electromagnetic survey using finite element method.Journal of Electromagnetic Analysis and Applications,4(3):135-146.

Yin C C,Ben F,Liu Y H,et al.2014.MCSEM 3D modeling for arbitrarily anisotropic media.Chinese Journal of Geophysics (in Chinese),57(12):4110-4122,doi:10.6038/cjg20141222.Yoshimura R,Oshiman N.2002.Edge-based finite element approach to the simulation of geoelectromagnetic induction in a 3-D sphere.Geophysical Research Letters,29(3):9-1-9-4.Zhang Y,Wang H N,Tao H G,et al.2012.Finite volume algorithm to simulate 3D responses of multi-component induction tools in inhomogeneous anisotropic formation based on coupled scalar-vector potentials.Chinese Journal of Geophysics (in Chinese),55(6):2141-2152,doi:10.6038/j.issn.0001-5733.2012.06.036.

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.

Zhou J M,Zhang Y,Wang H N,et al.2014.Efficient simulation of three-dimensional marine controlled-source electromagnetic response in anisotropic formation by means of coupled potential finite volume method.Acta Physica Sinica (in Chinese),63(15):159101.

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

陳輝,鄧居智,譚捍東等.2011.大地電磁三維交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬中的散度校正方法研究.地球物理學(xué)報(bào),54(6):1649-1659,doi:10.3969/j.issn.0001-5733.2011.06.025.

鄧居智,譚捍東,陳輝等.2011.CSAMT三維交錯(cuò)采樣有限差分?jǐn)?shù)值模擬.地球物理學(xué)進(jìn)展.26(6):2026-2032,doi:10.3969/j.issn.1004-2903.2011.06.017.

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

譚捍東,余欽范,Booker J等.2003.大地電磁法三維交錯(cuò)采樣有限差分?jǐn)?shù)值模擬.地球物理學(xué)報(bào),46(5):705-711.

殷長春,賁放,劉云鶴等.2014.三維任意各向異性介質(zhì)中海洋可控源電磁法正演研究.地球物理學(xué)報(bào),57(12):4110-4122,doi:10.6038/cjg20141222.

張燁,汪宏年,陶宏根等.2012.基于耦合標(biāo)勢(shì)與矢勢(shì)的有限體積法模擬非均勻各向異性地層中多分量感應(yīng)測井三維響應(yīng).地球物理學(xué)報(bào),55(6):2141-2152,doi:10.6038/j.issn.0001-5733.2012.06.036.

周建美,張燁,汪宏年等.2014.耦合勢(shì)有限體積法高效模擬各向異性地層中海洋可控源的三維電磁響應(yīng).物理學(xué)報(bào),63(15):159101.

(本文編輯何燕)

A finite-volume solution to 3D frequency-domain electromagnetic modelling using Lorenz-gauged magnetic vector and scalar potentials

CHEN Hui1,2,YIN Chang-Chun1*,DENG Ju-Zhi2

1 Geo-Exploration Science and Technology Institute,Jilin University,Changchun 130026,China 2 Key Laboratory of Radioactive Geology and Exploration Technology Fundamental Science for National Defense， East China Institute of Technology，Nanchang 330013,China

To eliminate the coupling between air and the earth and to avoid calculation of primary field,we develop a universal algorithm on forward modelling of the frequency-domain electromagnetic (EM)method for different source types and applications.We first present a magnetic vector and a scalar potential with Lorenz gauge based on Maxwell equation,then divide different sources into a lot of short wires.Next,we use staggered grids and the finite volume method to discrete the potential equations and implement the quasi-minimum-residual (QMR)iteration with Jacob to solve the large sparse and symmetrical linear equations system.Finally,we accomplish frequency-domain EM modelling for different source types and application areas.Through analyses and comparison of 1D and 3D MT,CSAMT and loop source responses,we demonstrate the efficiency and accuracy of the algorithm proposed in this study.Based on this,we analyze EM responses for different source types.The algorithm and numerical results presented in this paper build a framework for 3D frequency-domain EM forward modelling and inversion.KeywordsFDEM;Finite-volume method;Forward modelling;Lorenz-gauged;Coupled potentials

陳輝,殷長春,鄧居智.2016.基于Lorenz規(guī)范條件下磁矢勢(shì)和標(biāo)勢(shì)耦合方程的頻率域電磁法三維正演.地球物理學(xué)報(bào)，59(8):3087-3097,

10.6038/cjg20160831.

Chen H,Yin C C,Deng J Z.2016.A finite-volume solution to 3D frequency-domain electromagnetic modelling using Lorenz-gauged magnetic vector and scalar potentials.Chinese J.Geophys.(in Chinese),59(8):3087-3097,doi:10.6038/cjg20160831.

國家青年基金項(xiàng)目(41404057),國家自然科學(xué)基金項(xiàng)目(41164003),國家重大科研裝備研究項(xiàng)目(ZDYZ2012-1-03和20130523MTEM05)聯(lián)合資助.

陳輝，男，1985年8月生，博士生，主要從事電磁勘查技術(shù)正反演研究.E-mail:schoolhui@163.com

殷長春，男，1965年生，教授，國家“千人計(jì)劃”特聘專家，主要從事電磁勘探理論，特別是航空和海洋電磁方面的研究.

E-mail:yinchangchun@jlu.edu.cn

10.6038/cjg20160831

P631

2015-01-21，2015-10-08收修定稿