李建慧, Colin G. Farquharson, 胡祥云, 曾思紅

1 中國(guó)地質(zhì)大學(xué)地球物理與空間信息學(xué)院，地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢　430074 2 紐芬蘭紀(jì)念大學(xué)地球科學(xué)系，圣約翰斯， 加拿大 3 中國(guó)礦業(yè)大學(xué)深部巖土力學(xué)與地下工程國(guó)家重點(diǎn)實(shí)驗(yàn)室， 江蘇徐州　221116

基于電場(chǎng)總場(chǎng)矢量有限元法的接地長(zhǎng)導(dǎo)線(xiàn)源三維正演

李建慧1,2,3, Colin G. Farquharson2, 胡祥云1*, 曾思紅1

1 中國(guó)地質(zhì)大學(xué)地球物理與空間信息學(xué)院，地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢430074 2 紐芬蘭紀(jì)念大學(xué)地球科學(xué)系，圣約翰斯， 加拿大 3 中國(guó)礦業(yè)大學(xué)深部巖土力學(xué)與地下工程國(guó)家重點(diǎn)實(shí)驗(yàn)室， 江蘇徐州221116

摘要電磁場(chǎng)數(shù)值模擬的背景場(chǎng)/異常場(chǎng)算法是三維正演的有效策略之一，優(yōu)點(diǎn)為采用解析法計(jì)算電磁場(chǎng)背景場(chǎng)代替場(chǎng)源項(xiàng)、克服了場(chǎng)源奇異性，缺點(diǎn)為不適用于發(fā)射源布置于起伏地表或背景模型復(fù)雜的情形.總場(chǎng)算法是直接對(duì)電磁場(chǎng)總場(chǎng)開(kāi)展數(shù)值模擬，其難點(diǎn)是有效加載場(chǎng)源、保證近區(qū)與過(guò)渡區(qū)數(shù)值解精度.本文以水平電偶源形式分段加載接地長(zhǎng)導(dǎo)線(xiàn)源，并以電場(chǎng)總場(chǎng)Helmholtz方程為矢量有限元法控制方程，實(shí)現(xiàn)了基于非結(jié)構(gòu)化四面體網(wǎng)格剖分的接地長(zhǎng)導(dǎo)線(xiàn)源頻率域電磁法三維正演.通過(guò)與均勻全空間中水平電偶源產(chǎn)生的電場(chǎng)解析解對(duì)比，驗(yàn)證了本文算法的正確性，并分析了四面體外接圓半徑與其最短棱邊的最大比值和四面體二面角最小值對(duì)數(shù)值解精度的影響規(guī)律.通過(guò)與塊狀高導(dǎo)體地電模型的積分方程法、有限體積法和基于磁矢量勢(shì)Helmholtz方程的有限元法數(shù)值解對(duì)比，進(jìn)一步驗(yàn)證了本文算法正確性，同時(shí)說(shuō)明了非結(jié)構(gòu)化四面體網(wǎng)格能夠更加精細(xì)地剖分電性異常體，利于獲得精確數(shù)值解.

關(guān)鍵詞接地長(zhǎng)導(dǎo)線(xiàn)源； 三維正演； 矢量有限元法； 電場(chǎng)總場(chǎng)； 非結(jié)構(gòu)化網(wǎng)格

The numerical accuracy of FE method depends on meshing quality to some extent, which is evaluated by the maximum allowable radius-edge ratio and the minimum allowable dihedral angle for the unstructured tetrahedral elements generated by the TetGen software. Taking the electric field excited by a HED within a 0.01 S·m-1whole-space as examples, the results show that as the maximum allowable radius-edge ratio decreases or the minimum allowable dihedral angle increases, the numerical accuracy will be improved. By this model, we also validated the algorithm presented here.

For a conductive prism buried in a homogeneous half-space with a 100m-long grounded wire, the electric field calculated by our algorithm was compared with those calculated by the integral equation method based on secondary electric field, the finite volume method based on total electric field and the FE method based on magnetic vector potential. The results show that these four numerical solutions coincide well with each other, and the behavior of the electric field well indicates the conductive prism. Then, we applied the algorithm presented here to compute the amplitude and phase of the electric field for the model of disc-shape hydrocarbon buried in seabed. Through this model, the validity and ability of modelling the electromagnetic field for irregular bodies were tested simultaneously. Finally, an inclined fault, contact zone and metalliferous vein are always approximated by an inclined plate in three-dimensional modelling. For the inclined plate model with different conductivities, the electric field excited by a 1000m-long grounded wire source was calculated. The electric field for the inclined plate of 0.0333 S·m-1has stronger anomaly responses than those for the inclined plate of 0.0167 S·m-1.

In the following study, this algorithm will be applied to three-dimensional modelling for a long grounded wire source laid on a rugged surface.

1引言

電磁法勘探中，電性源裝置已成功應(yīng)用于金屬礦產(chǎn)勘探(Hu et al.，2013)、環(huán)境水文調(diào)查與監(jiān)測(cè)(Fu et al.，2013；Grayver et al.，2014)、海底油氣藏勘探等領(lǐng)域(Zhdanov et al.，2014).為了更加精細(xì)地反映地電結(jié)構(gòu)、進(jìn)一步改善勘探效果，開(kāi)發(fā)精度高、速度快、能夠模擬復(fù)雜地電模型的三維正演技術(shù)已成為必由之路.

電性源裝置具體包括電偶源和接地長(zhǎng)導(dǎo)線(xiàn)源.除一維解析法外，電性源的正演按照維度還可分為二維、2.5維和三維.二維正演問(wèn)題是指地電模型為二維、接地長(zhǎng)導(dǎo)線(xiàn)源沿地電模型走向無(wú)限延伸的情況，即二維模型、二維場(chǎng)源問(wèn)題.閻述和陳明生(1999a，1999b)發(fā)展和完善了無(wú)限長(zhǎng)接地導(dǎo)線(xiàn)的有限單元法二維正演；陳小斌和胡文寶(2002)利用有限元法二維正演分析了頻率域電磁測(cè)深近區(qū)觀測(cè)的可行性；王若等(2006)開(kāi)展有限元法正演時(shí)，在計(jì)算區(qū)域外邊界加載了一階吸收邊界條件，并將偽δ函數(shù)引入場(chǎng)源項(xiàng)，提高了近區(qū)數(shù)值解精度.

2.5維正演問(wèn)題是指地電模型為二維、電性源為電偶源或接地長(zhǎng)導(dǎo)線(xiàn)源的情況，即二維模型、三維場(chǎng)源問(wèn)題.開(kāi)展2.5維正演時(shí)，首先利用傅里葉變換將與地質(zhì)構(gòu)造走向一致的電磁場(chǎng)分量轉(zhuǎn)換至波數(shù)域，再利用有限元法等數(shù)值方法求取一系列波數(shù)的電磁場(chǎng)響應(yīng)；最后利用傅里葉逆變換將波數(shù)域電磁場(chǎng)響應(yīng)轉(zhuǎn)換至空間域，進(jìn)而獲得二維模型的三維電磁場(chǎng)響應(yīng).基于電磁場(chǎng)數(shù)值模擬的背景場(chǎng)/異常場(chǎng)算法，Unsworth等(1993)利用有限元法實(shí)現(xiàn)了水平電偶源的頻率域電磁法2.5維正演問(wèn)題，并有效克服了場(chǎng)源的奇異性問(wèn)題；底青云等(2004)將該算法應(yīng)用于復(fù)雜介質(zhì)的正演計(jì)算，為反演解釋提供了良好依據(jù).Li和Key(2007)基于背景場(chǎng)/異常場(chǎng)算法，利用自適應(yīng)有限元法實(shí)現(xiàn)了海洋可控源電磁法的2.5維正演問(wèn)題，場(chǎng)源為水平電偶源.2.5維正演問(wèn)題的另一種策略是采用電磁場(chǎng)數(shù)值模擬的總場(chǎng)算法，如何克服場(chǎng)源奇異性對(duì)數(shù)值解的不利影響是該策略的關(guān)鍵問(wèn)題.Mitsuhata(2000)以偽δ函數(shù)代替δ函數(shù)處理場(chǎng)源項(xiàng)、克服了場(chǎng)源奇異性，實(shí)現(xiàn)了基于有限元法的水平電偶源頻率域電磁法2.5維正演；Mitsuhata等(2002)延續(xù)上述思路，開(kāi)展了接地長(zhǎng)導(dǎo)線(xiàn)源頻率域電磁法的2.5維正演和反演研究.

三維正演問(wèn)題是指地電模型為三維、電性源為電偶源或接地長(zhǎng)導(dǎo)線(xiàn)源的情況，即三維模型、三維場(chǎng)源問(wèn)題.基于電磁場(chǎng)數(shù)值模擬的背景場(chǎng)/異常場(chǎng)算法，國(guó)內(nèi)外許多學(xué)者開(kāi)展了電性源頻率域電磁法的三維正演問(wèn)題，F(xiàn)arquharson 和Oldenburg(2002)采用了積分方程法；Weiss 和Constable(2006)、韓波等(2015)采用了有限體積法；Sasaki 和Meju(2009)、Streich(2009)采用了交錯(cuò)網(wǎng)格有限差分法；Schwarzbach等(2011)、李勇等(2015)采用了矢量有限元法.采用電磁場(chǎng)數(shù)值模擬總場(chǎng)算法是電性源頻率域電磁法三維正演的另一有效策略.閻述和陳明生(2000)對(duì)利用矢量有限元法開(kāi)展電偶源頻率域電磁法三維正演進(jìn)行了探索研究.張繼鋒等(2009)基于電場(chǎng)矢量波動(dòng)方程，利用節(jié)點(diǎn)有限元法實(shí)現(xiàn)了水平電偶源的頻率域電磁法三維正演，以偽δ函數(shù)代替δ函數(shù)處理場(chǎng)源項(xiàng)，改善了近區(qū)數(shù)值解精度；并以水平電偶源的一維解析解作為計(jì)算區(qū)域外邊界條件.徐志鋒和吳小平(2010)基于磁矢量勢(shì)Helmholtz方程，采用穩(wěn)定化節(jié)點(diǎn)有限元法實(shí)現(xiàn)了水平電偶源的頻率域電磁法三維正演，也以水平電偶源的一維解析解作為計(jì)算區(qū)域外邊界條件.Ansari和Farquharson(2014)基于磁矢量勢(shì)Helmholtz方程，采用基于非結(jié)構(gòu)化四面體的矢量有限元法實(shí)現(xiàn)了水平電偶源和接地長(zhǎng)導(dǎo)線(xiàn)源的頻率域電磁法三維正演，將接地長(zhǎng)導(dǎo)線(xiàn)看作有限多個(gè)水平電偶源的組合，消除了場(chǎng)源奇異性對(duì)近區(qū)數(shù)值解的不利影響.Jahandari 和Farquharson(2014)基于電場(chǎng)總場(chǎng)Helmholtz方程，采用基于非結(jié)構(gòu)化四面體和Vorono?網(wǎng)格的有限體積法實(shí)現(xiàn)了水平電偶源和接地長(zhǎng)導(dǎo)線(xiàn)源的頻率域電磁法三維正演，也將接地長(zhǎng)導(dǎo)線(xiàn)看作有限多個(gè)水平電偶源的組合.

以接地長(zhǎng)導(dǎo)線(xiàn)源作為發(fā)射裝置的電磁法有廣域電磁法(何繼善，2010)、短偏移距瞬變電磁法(薛國(guó)強(qiáng)等，2014)、海洋可控源電磁法(韓波等，2015)，這些方法的共同點(diǎn)為測(cè)點(diǎn)不僅位于遠(yuǎn)區(qū)，而且還分布于過(guò)渡區(qū)和近區(qū).如果仍將發(fā)射源看作水平電偶源開(kāi)展數(shù)值模擬，將不利于獲得真實(shí)的電磁場(chǎng)響應(yīng).本研究將接地長(zhǎng)導(dǎo)線(xiàn)源看作有限多個(gè)水平電偶源的組合，基于電場(chǎng)總場(chǎng)Helmholtz方程，采用基于非結(jié)構(gòu)化四面體網(wǎng)格剖分的矢量有限元法開(kāi)展頻率域電磁法三維正演研究.矢量有限元法將自由度定義在棱邊上，自動(dòng)滿(mǎn)足電場(chǎng)切向分量連續(xù)、法向分量不連續(xù)的連續(xù)條件，因此可避免節(jié)點(diǎn)有限元法中強(qiáng)加法向分量連續(xù)而產(chǎn)生的偽解.目前，該方法已廣泛應(yīng)用于電磁法勘探的三維正演研究中(Schwarzbach et al.，2011；Li et al.,2011；Ansari and Farquharson，2014；李勇等，2015).非結(jié)構(gòu)化四面體網(wǎng)格剖分不僅使得本文算法能夠模擬更加復(fù)雜的地電模型，而且能夠通過(guò)采取局部加密的方法，細(xì)化場(chǎng)源、測(cè)點(diǎn)和異常體所在區(qū)域網(wǎng)格剖分，進(jìn)而提高數(shù)值解精度.

2矢量有限單元法

2.1控制方程

在準(zhǔn)靜態(tài)、采用正諧時(shí)eiω t的條件下，麥克斯韋方程組中的電場(chǎng)E和磁場(chǎng)H旋度方程為：

(1)

(2)

式(1)中，B為磁感應(yīng)強(qiáng)度，μ為磁導(dǎo)率，B=μH.本研究不考慮磁導(dǎo)率變化對(duì)電磁場(chǎng)的影響，故取真空中磁導(dǎo)率μ0.ω為角頻率，其值等于2πf，f為頻率.式(2)中，σ為地電模型電導(dǎo)率，Js為場(chǎng)源電流密度.將式(2)代入式(1)，得頻率域電場(chǎng)總場(chǎng)Helmholtz方程：

(3)

當(dāng)場(chǎng)源為沿x方向分布的電偶源時(shí)，Js為(WardandHohmann，1988)：

(4)

其中I與ds分別表示電流強(qiáng)度與x方向長(zhǎng)度，δ表示狄拉克函數(shù).

2.2單元分析

以式(3)為控制方程，應(yīng)用Galerkin法推導(dǎo)有限元方程.設(shè)余量r為

(5)

令r在計(jì)算區(qū)域Ω滿(mǎn)足

(6)

(7)經(jīng)矢量恒等變換，(7)式改寫(xiě)為：

(8)

本研究中，將計(jì)算區(qū)域剖分為有限多個(gè)四面體單元e.每個(gè)四面體e中，點(diǎn)(x,y,z)處的電場(chǎng)矢量場(chǎng)可由各條棱邊的電場(chǎng)和矢量基函數(shù)表示，有

(9)

J為單元e的棱邊數(shù)，四面體單元的棱邊數(shù)為6，j值范圍為1至6.矢量基函數(shù)Nj可寫(xiě)作(Jin，2002)：

(10)

其中j1和j2為四面體頂點(diǎn)編號(hào)，其值范圍為1至4.lj為某一棱邊長(zhǎng)度，Nj1與Nj2為該棱邊對(duì)應(yīng)兩節(jié)點(diǎn)的標(biāo)量基函數(shù)，棱邊編號(hào)與節(jié)點(diǎn)編號(hào)如圖1所示.將(9)式代入(8)式，有

-∫Ωiωμ0N·JsdΩ,

(11)

圖1　四面體單元中，節(jié)點(diǎn)與棱邊編號(hào)Fig.1　The numbering diagram of the edges and nodes in a tetrahedral element

式中，矢量基函數(shù)的旋度可表示為(Jin，2002)：

(12)

關(guān)于式(11)左端項(xiàng)的推導(dǎo)過(guò)程詳見(jiàn)Jin(2002)的著作，而其右端項(xiàng)僅存在于場(chǎng)源分布單元.

(13)

2.3求解方程組

與迭代法相比，直接法求解大型方程組具有精度高、穩(wěn)定性好的優(yōu)點(diǎn)，但其內(nèi)存需求大、求解效率低.近年來(lái)計(jì)算機(jī)硬件和并行計(jì)算的快速發(fā)展，直接法效率已有明顯提升，并已大量應(yīng)用于電磁法三維數(shù)值模擬(Streich，2009；Schwarzbach et al.，2011；Jahandari and Farquharson，2014).本研究將采用直接法開(kāi)源軟件MUMPS(Amestoy et al.，2006)求解大型線(xiàn)性方程組；由于矢量有限元法生成的系數(shù)矩陣具有對(duì)稱(chēng)性，因此只需存儲(chǔ)和輸入下三角矩陣.獲得各棱邊電場(chǎng)值后，即可通過(guò)(9)式計(jì)算單元內(nèi)任意點(diǎn)的電場(chǎng)矢量場(chǎng).

3網(wǎng)格剖分與算法驗(yàn)證

本研究采用的開(kāi)源代碼或程序如下：求解方程組采用MUMPS 5.0.0(Amestoy et al.，2006)、非結(jié)構(gòu)化四面體網(wǎng)格剖分采用TetGen 1.5.1-beta1(Si，2015)，而非結(jié)構(gòu)化四面體網(wǎng)格顯示采用ParaView 4.3.1(Ayachit，2015).另外，計(jì)算環(huán)境為：采用Ubuntu 14.04系統(tǒng)和GNU Fortran編譯器，硬件為Inter I5處理器、16 G內(nèi)存臺(tái)式機(jī).

有限元法數(shù)值解精度一定程度上取決于網(wǎng)格剖分質(zhì)量.在TetGen軟件中，非結(jié)構(gòu)化四面體形狀和尺度由四面體外接圓半徑與其最短棱邊比值和四面體二面角(如圖1所示)角度控制(Si，2013，2015).Ansari和Farquharson(2014)采用外接圓半徑與四面體最短棱邊最大比值為1.414、二面角最小值為16°的網(wǎng)格剖分方案，使用矢量有限元法取得了理想的數(shù)值結(jié)果；而Jahandari和Farquharson(2014)將四面體外接圓半徑與其最短棱邊的最大比值設(shè)置在1.1～1.2之間，并未限制四面體二面角角度，采用有限體積法也取得了高精度的數(shù)值解.

本研究將以置于0.01 S·m-1均勻全空間中的水平電偶源為例，分析四面體外接圓半徑與最短棱邊比值和二面角角度對(duì)有限元法數(shù)值解精度的影響.為了避免計(jì)算區(qū)域尺度不足引起的數(shù)值解誤差，本例計(jì)算區(qū)域尺度為400 km ×400 km ×400 km.表1為網(wǎng)格剖分方案技術(shù)參數(shù)，四面體外接圓半徑與其最短棱邊的最大比值(R)分別為1.2、1.3和1.4，二面角最小值(DA)分別為未限制、12°、14°、16°和18°.圖2為R值等于1.4、DA值未限制條件下的非結(jié)構(gòu)化網(wǎng)格剖分示意圖，其中圖2b為圖2a中白色線(xiàn)條所圍區(qū)域的放大圖.非結(jié)構(gòu)化四面體網(wǎng)格剖分時(shí)，采用局部加密策略對(duì)電偶源和測(cè)點(diǎn)分布區(qū)域進(jìn)行加密.加密方案具體如下：將電偶源或測(cè)點(diǎn)作為一個(gè)正四面體中心，該四面體的四個(gè)頂點(diǎn)即為加密點(diǎn).本例中，正四面體棱邊長(zhǎng)度為1 m.

在準(zhǔn)靜態(tài)、采用正諧時(shí)eiω t的條件下，沿x方向布置的電偶源在均勻全空間中激發(fā)的電場(chǎng)Ex分量解析解為(Ward and Hohmann，1988)：

表1　基于非結(jié)構(gòu)化四面體的網(wǎng)格剖分方案

+(k2r2-ikr-1)],

(14)

式中，I為電流強(qiáng)度，ds為電偶源長(zhǎng)度，σ為電導(dǎo)率，波數(shù)k由k2=-iμ0σ ω定義，ω為角頻率，r為收發(fā)距.本例中，坐標(biāo)原點(diǎn)布置于電偶源幾何中心，26個(gè)測(cè)點(diǎn)y和z坐標(biāo)均為0 m，x坐標(biāo)范圍為50 m 至10000 m，即x坐標(biāo)即為收發(fā)距.

圖3為電場(chǎng)Ex分量的有限元法數(shù)值解及其與解析解的相對(duì)誤差曲線(xiàn)，圖中黑色虛線(xiàn)表示3%的誤差限.如圖3a所示：采用3種網(wǎng)格剖分方案，有限元法計(jì)算的數(shù)值解與解析解吻合良好；僅收發(fā)距為9000 m和10000 m，且R值為1.4的數(shù)值解精度略差.圖3b至3f展現(xiàn)了電場(chǎng)實(shí)部和虛部相對(duì)誤差隨二面角最小值變化的規(guī)律.不難發(fā)現(xiàn)實(shí)部相對(duì)誤差均大于虛部相對(duì)誤差，因此本文采用實(shí)部相對(duì)誤差作為分析指標(biāo).綜合分析圖3b、3c和3d，我們發(fā)現(xiàn)：隨著收發(fā)距的增大，實(shí)部與虛部相對(duì)誤差總體呈現(xiàn)增大趨勢(shì).當(dāng)R值為1.4時(shí)，僅有收發(fā)距小于1000 m的少數(shù)測(cè)點(diǎn)實(shí)部相對(duì)誤差小于3%；當(dāng)R值為1.3時(shí)，收發(fā)距小于6000 m的大多數(shù)測(cè)點(diǎn)實(shí)部相對(duì)誤差小于3%；當(dāng)R值為1.2時(shí)，僅有收發(fā)距大于6000 m的少數(shù)測(cè)點(diǎn)實(shí)部相對(duì)誤差大于3%.圖3e和3f中，僅有收發(fā)距大于7000 m的少數(shù)實(shí)部相對(duì)誤差大于3%.如圖3b至3f所示：DA值固定時(shí)，電場(chǎng)實(shí)部與虛部相對(duì)誤差隨著R值減小而呈現(xiàn)出遞減的變化規(guī)律.另外如表1所示：DA值為18°的網(wǎng)格剖分過(guò)密，單元數(shù)與棱邊數(shù)激增.比如R值為1.2、DA值為18°的網(wǎng)格棱邊數(shù)高達(dá)74095813.這一網(wǎng)格剖分規(guī)模遠(yuǎn)大于本研究所采用硬件的允許計(jì)算規(guī)模(經(jīng)測(cè)試約為1800000條棱邊)，因此并未計(jì)算該網(wǎng)格剖分的有限元法數(shù)值解.

結(jié)合表1和圖3綜合分析，隨著R值的減小或DA值的增大，非結(jié)構(gòu)化網(wǎng)格剖分的四面體和棱邊數(shù)目呈遞增趨勢(shì)、剖分網(wǎng)格逐漸精細(xì)，有限元法數(shù)值解精度也隨之提高.綜合考慮計(jì)算效率和精度，本文后續(xù)研究均采用R值為1.3和DA值為16°的網(wǎng)格剖分參數(shù).

4算例

4.1塊狀高導(dǎo)體模型

如圖4所示：一個(gè)電導(dǎo)率為0.2 S·m-1的塊狀體埋置于電導(dǎo)率為0.02 S·m-1的均勻半空間中，其三維尺度(依次為x、y、z方向)為120 m×200 m×400 m，中心位于(1000,0,300) m；空氣電導(dǎo)率設(shè)置為10-8S·m-1.測(cè)線(xiàn)布設(shè)于地表，起始位置為(400,0,0)m，終止位置為(1600,0,0)m.發(fā)射源是長(zhǎng)度為100 m的接地長(zhǎng)導(dǎo)線(xiàn)，其中心坐標(biāo)為(50,0,0)m.發(fā)射電流為1 A、發(fā)射頻率為3 Hz.本例計(jì)算區(qū)域?yàn)?0 km×50 km×50 km，共剖分325207個(gè)四面體、378508條棱邊.為了改善數(shù)值解精度，采用局部加密技術(shù)對(duì)接地長(zhǎng)導(dǎo)線(xiàn)和測(cè)點(diǎn)所在區(qū)域網(wǎng)格進(jìn)行了細(xì)分，并限制塊狀高導(dǎo)體所在區(qū)域的剖分單元最大體積為5000 m3(如圖5所示).另外，本例計(jì)算時(shí)間約為600 s.

圖6為塊狀高導(dǎo)體的電場(chǎng)分量Ex響應(yīng)曲線(xiàn)，其中實(shí)部為正值、虛部為負(fù)值.為了進(jìn)一步驗(yàn)證本文算法的正確性，將本研究數(shù)值解與積分方程法(IE)(Farquharson and Oldenburg, 2002)、有限體積法(FV)(Jahandari and Farquharson, 2014)和基于磁矢量勢(shì)Helmholtz方程的有限元法(FE)(Ansari and Farquharson, 2014)數(shù)值解進(jìn)行了對(duì)比.由圖6可看出：除積分方程法外，其他三種方法的數(shù)值解實(shí)部與虛部曲線(xiàn)形態(tài)一致、吻合程度良好；與其他三種方法相比，積分方程法數(shù)值解在塊狀體分布區(qū)域(x值為900至1100 m)偏大.積分方程法采用矩形塊網(wǎng)格剖分方案，高導(dǎo)塊狀體沿x、y和z方向剖分的網(wǎng)格數(shù)目依次為8×8×8，被剖分為512個(gè)矩形塊單元；其他三種方法均采用了非結(jié)構(gòu)化四面體單元，比如本文算法采用TetGen軟件將塊狀體剖分為7311個(gè)四面體單元.因此，我們推斷積分方程法數(shù)值解與其他三種方法數(shù)值解存在上述差異的潛在原因?yàn)榉e分方程法求解時(shí)，塊狀高導(dǎo)體的剖分單元數(shù)目不足，影響了數(shù)值解精度.

圖4　塊狀高導(dǎo)體模型示意圖(a) x-z平面斷面圖(y=0 m)； (b) x-y平面俯視圖.Fig.4　The diagrams for the conductive prism(a) The section view (y=0 m) in the x-z plane; (b) The top view in the x-y plane.

圖5　塊狀高導(dǎo)體模型的非結(jié)構(gòu)化網(wǎng)格剖分?jǐn)嗝鎴D(y=0 m)Fig.5　The section view (y=0 m) of the unstructured meshing for the conductive prism

另外，這些數(shù)值方法計(jì)算的電場(chǎng)響應(yīng)對(duì)高導(dǎo)塊狀體反映顯著，電場(chǎng)實(shí)部和虛部曲線(xiàn)在x為900至1100 m范圍內(nèi)均呈現(xiàn)向下凹陷現(xiàn)象.這是由于地下電流易被高導(dǎo)體“吸引”，流經(jīng)地表的電流密度降低所致.

4.2海底油氣藏模型

海底油氣藏資源是近十年來(lái)電磁法勘探的熱門(mén)領(lǐng)域，本例模型與典型油氣藏模型(Weiss and Constable，2006)類(lèi)似.如圖7所示，該地電模型包括電導(dǎo)率為3.2 S·m-1的海水層、1 S·m-1的海底以及0.01 S·m-1的油氣藏.油氣藏分布呈扁平圓柱體狀，半徑為1000 m，厚度為100 m，頂部埋深950 m.長(zhǎng)度為100 m的接地長(zhǎng)導(dǎo)線(xiàn)源分布于海水層，距離海底高度為100 m，其中心坐標(biāo)為(-50,0,-100) m.發(fā)射電流為1 A，發(fā)射頻率為1 Hz.計(jì)算區(qū)域?yàn)?0 km×30 km×16 km，其中海水層厚度為1 km，該模型共剖分1342046個(gè)四面體以及1562622條棱邊.如圖8a所示，我們依然對(duì)接地長(zhǎng)導(dǎo)線(xiàn)、測(cè)點(diǎn)和電性異常體所在區(qū)域進(jìn)行了局部加密，而圖8b說(shuō)明了采用TetGen軟件能夠精細(xì)地剖分扁平圓柱體狀油氣藏資源，同時(shí)也避免了采用結(jié)構(gòu)化網(wǎng)格易導(dǎo)致計(jì)算區(qū)域邊界處出現(xiàn)扁平單元和狹長(zhǎng)單元的現(xiàn)象.另外，本例計(jì)算時(shí)間約為3000 s.

圖6　均勻半空間中塊狀高導(dǎo)體地電模型的電場(chǎng)總場(chǎng)Ex分量Fig.6　The Ex component of total electric field for the conductive prism in the homogeneous half-space

圖7　海底油氣藏模型示意圖(a) x-z平面斷面圖(y=0 m)； (b) x-y平面俯視圖.Fig.7　The diagrams for the hydrocarbon buried in seabed(a) The section view (y=0 m) in the x-z plane; (b) The top view in the x-y plane.

圖8　海底油氣藏模型的非結(jié)構(gòu)化網(wǎng)格剖分結(jié)果(a) x-z平面斷面圖(y=0 m)； (b) x-y平面切片圖(z=1000 m).Fig.8　The results of the unstructured meshing for the hydrocarbon buried in seabed (a) The section view (y=0 m) in the x-z plane; (b) The slice (z=1000 m) in the x-y plane.

圖9　海底油氣藏模型的電場(chǎng)總場(chǎng)Ex分量“FE”表示Ansari與Farquharson(2014)的有限元法數(shù)值解.(a) 實(shí)部與虛部；(b) 振幅；(c) 相位.Fig.9　The Ex component of total electric field for the hydrocarbon buried in seabed The symbol of “FE” denotes the solution of the finite element method (Ansari and Farquharson, 2014). (a) Real and imaginary parts; (b) Amplitude; (c) Phase.

圖9為海底油氣藏地電模型的電場(chǎng)分量Ex實(shí)部與虛部、振幅和相位曲線(xiàn).由圖9可知，電場(chǎng)實(shí)部與虛部數(shù)量級(jí)從10-5降低到10-15，并存在多個(gè)極小值，因此對(duì)數(shù)值方法精度和穩(wěn)定性要求極高.如圖9a所示，本文算法計(jì)算的電場(chǎng)實(shí)部和虛部與基于磁矢量勢(shì)Helmholtz方程的有限元法(Ansari and Farquharson, 2014)計(jì)算結(jié)果吻合程度良好.兩種數(shù)值算法的主要區(qū)別為極小值的數(shù)值不同，即本研究計(jì)算結(jié)果偏大，這種現(xiàn)象由采樣不足引起.本研究測(cè)點(diǎn)間距為50 m，而Ansari與Farquharson(2014)的測(cè)點(diǎn)間距為5 m，其更能反映真實(shí)場(chǎng)值.頻率域海洋可控源電磁法中，通常結(jié)合航行數(shù)據(jù)，將實(shí)測(cè)信號(hào)轉(zhuǎn)換為振幅-收發(fā)距曲線(xiàn)與相位-收發(fā)距曲線(xiàn)(李予國(guó)和段雙敏，2014).因此，本研究將電場(chǎng)實(shí)部與虛部轉(zhuǎn)換為了電場(chǎng)振幅與相位曲線(xiàn).由圖9b和9c可知兩種數(shù)值方法的電場(chǎng)振幅與相位完全吻合，可見(jiàn)采樣間隔的這種差異對(duì)電場(chǎng)振幅與相位影響微小.

4.3傾斜板狀體模型

斷層破碎帶、接觸帶和金屬礦脈往往可由傾斜板狀體予以近似.如圖10所示：一頂部出露于地表的傾斜板狀體分布于電導(dǎo)率為0.01 S·m-1的均勻半空間，其三維尺度為200 m×3000 m×1500 m，電導(dǎo)率分別為0.0333 S·m-1和0.0167 S·m-1；空氣電導(dǎo)率設(shè)置為10-8S·m-1.板狀體沿走向方向(y方向)分布范圍為[1500,4500] m；沿x方向分布范圍： 頂部為[50,250] m，底部為[-250,-50] m.

圖10　傾斜板狀體模型示意圖(a) x-z平面斷面圖(y=3000 m)； (b) x-y平面俯視圖.Fig.10　The diagrams for the inclined plate(a) The section view (y=3000 m) in the x-z plane； (b) The top view in the x-y plane.

本例在地表(z=0 m)布置了y坐標(biāo)分別為1000、2000、3000、4000和5000 m的5條測(cè)線(xiàn)、共105個(gè)測(cè)點(diǎn)，其中每條測(cè)線(xiàn)長(zhǎng)度為2000 m、含21個(gè)測(cè)點(diǎn).布置于x方向的接地長(zhǎng)導(dǎo)線(xiàn)源長(zhǎng)度為1000 m，中心置于(0,0,0) m.發(fā)射電流為1 A、發(fā)射頻率為1 Hz和32 Hz.本例計(jì)算區(qū)域?yàn)?0 km×50 km×50 km，共剖分1518786個(gè)四面體、1766490條棱邊.該地電模型的非結(jié)構(gòu)化四面體網(wǎng)格剖分?jǐn)嗝鎴D(y=3000 m)如圖11所示.

圖12為傾斜板狀體模型的電場(chǎng)Ex分量振幅曲線(xiàn).由圖12所示，兩個(gè)頻率的電場(chǎng)振幅曲線(xiàn)隨測(cè)點(diǎn)x坐標(biāo)變化的形態(tài)一致；相比于均勻半空間模型，電場(chǎng)振幅對(duì)傾斜高導(dǎo)板狀體模型反映顯著，其中0.0333 S·m-1模型的電場(chǎng)振幅比0.0167 S·m-1情形的數(shù)值更大，尤其是x坐標(biāo)從-100 m和400 m的六個(gè)測(cè)點(diǎn).兩個(gè)頻率電場(chǎng)振幅的主要區(qū)別在于遠(yuǎn)離測(cè)線(xiàn)中心的測(cè)點(diǎn)對(duì)高導(dǎo)板狀體的反映能力.當(dāng)x坐標(biāo)為-1000～-200 m和500～1000 m時(shí)，傾斜板狀體1 Hz的電場(chǎng)振幅略大于均勻半空間相應(yīng)頻率的電場(chǎng)振幅；而32 Hz時(shí)，兩者電場(chǎng)振幅近似相等.

圖11　傾斜板狀體的非結(jié)構(gòu)化網(wǎng)格剖分?jǐn)嗝鎴D(y=3000 m)Fig.11　The section view (y=3000 m) of the unstructured meshing for the inclined plate

圖12　傾斜板狀體模型的電場(chǎng)總場(chǎng)Ex分量振幅(y=3 km)Fig.12　The amplitude of the Ex component for the inclined plate (y=3 km)

圖13　傾斜板狀體模型的電場(chǎng)總場(chǎng)Ex分量振幅(1 Hz)(a) 均勻半空間； (b) 0.0167 S·m-1的板狀體； (c) 0.0333 S·m-1的板狀體.Fig.13　The amplitude of the Ex component for the inclined plate (1 Hz)(a) The homogeneous half-space； (b) The plate of 0.0167 S·m-1； (c) The plate of 0.0333S·m-1.

圖13和14為傾斜板狀體模型的電場(chǎng)分量Ex振幅等值線(xiàn)圖，其共同點(diǎn)是電場(chǎng)振幅隨著y坐標(biāo)的增大而逐漸衰減.圖13中1 Hz的電場(chǎng)振幅從1.2×10-5V·m-1衰減至9.0×10-8V·m-1，而圖14中32 Hz的電場(chǎng)振幅從1.9×10-5V·m-1衰減至1.0×10-7V·m-1.這一現(xiàn)象反映了高頻電磁場(chǎng)隨著傳播距離的增大而衰減更快的客觀規(guī)律.均勻半空間模型中，電場(chǎng)振幅等值線(xiàn)關(guān)于y軸對(duì)稱(chēng)；傾斜板狀體模型中，電場(chǎng)振幅等值線(xiàn)并不對(duì)稱(chēng).如圖13b、13c、14b和14c所示：由于受高導(dǎo)板狀體影響，x方向[0,300]m、y方向[1800,4500]m區(qū)域內(nèi)的等值線(xiàn)向下凹陷；與0.0167 S·m-1的板狀體模型相比，0.0333 S·m-1板狀體模型的等值線(xiàn)呈現(xiàn)出更加顯著的向下凹陷現(xiàn)象.

5結(jié)論

本研究將接地長(zhǎng)導(dǎo)線(xiàn)源視為有限多個(gè)水平電偶源的組合，并以水平電偶源的形式實(shí)現(xiàn)了場(chǎng)源的分段加載，有效去除了電磁法三維正演中外加場(chǎng)源對(duì)近區(qū)數(shù)值解精度的不利影響.通過(guò)將本文算法數(shù)值解與均勻全空間中水平電偶源產(chǎn)生的電場(chǎng)解析解對(duì)比研究，分析得出非結(jié)構(gòu)化網(wǎng)格剖分中，四面體外接圓半徑與其最短棱邊最大比值的減小或者四面體二面角最小值的增大都能提高有限元法數(shù)值解精度.然后，將本文算法數(shù)值解與塊狀高導(dǎo)體地電模型的積分方程法、有限體積法和基于磁矢量勢(shì)Helmholtz方程的有限元法數(shù)值解對(duì)比分析，四種數(shù)值解吻合良好、對(duì)高導(dǎo)異常體反映顯著，進(jìn)一步驗(yàn)證了本文算法的正確性.本例中，基于矩形塊網(wǎng)格剖分的積分方程法在異常體區(qū)域的數(shù)值解比基于非結(jié)構(gòu)化四面體網(wǎng)格剖分的三種方法數(shù)值解偏大，說(shuō)明了非結(jié)構(gòu)化四面體網(wǎng)格能夠更加精細(xì)地剖分電性異常體，有利于獲得精確數(shù)值解.海底油氣藏和傾斜板狀體地電模型主要說(shuō)明了基于非結(jié)構(gòu)化四面體網(wǎng)格剖分的數(shù)值算法能夠精確刻畫(huà)復(fù)雜地質(zhì)體，避免了結(jié)構(gòu)化網(wǎng)格剖分復(fù)雜地質(zhì)體易導(dǎo)致的扁平單元和狹長(zhǎng)單元.后續(xù)研究中，我們將致力于起伏地表接地長(zhǎng)導(dǎo)線(xiàn)源電磁法的三維正演算法.致謝感謝紐芬蘭紀(jì)念大學(xué)地球科學(xué)系SeyedMasoud Ansari與Hormoz Jahandari博士在三維正演程序編寫(xiě)中的熱情幫助.感謝國(guó)家留學(xué)基金委全額資助本文第一作者在紐芬蘭紀(jì)念大學(xué)開(kāi)展為期1年的博士后研究.特別感謝兩位審稿專(zhuān)家的中肯建議，使得本文質(zhì)量顯著提升.

圖14　傾斜板狀體模型的電場(chǎng)總場(chǎng)Ex分量振幅(32 Hz)(a) 均勻半空間； (b) 0.0167 S·m-1的板狀體； (c) 0.0333 S·m-1的板狀體.Fig.14　The amplitude of the Ex component for the inclined plate (32 Hz)(a) The homogeneous half-space； (b) The plate of 0.0167 S·m-1； (c) The plate of 0.0333 S·m-1.

References

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.

Ayachit U. 2015. The ParaView Guide: A Parallel Visualization Application. ParaView 4.3 ed. Carrboro, NC: Kitware, Inc.

Chen X B, Hu W B. 2002. Direct iterative finite element (DIFE) algorithm and its application to electromagnetic response modeling of line current source.ChineseJ.Geophys. (in Chinese), 45(1): 119-130, doi: 10.3321/j.issn:0001-5733.2002.01.015. Di Q Y, Unsworth M, Wang M Y. 2004. 2.5-D CSAMT modeling with the finite element method over 2-D complex earth media.ChineseJ.Geophys. (in Chinese), 47(4): 723-730, doi: 10.3321/j.issn:0001-5733.2004.04.026. Farquharson C G, Oldenburg D W. 2002. An integral equation solution to the geophysical electromagnetic forward-modelling problem. ∥ Zhdanov M S, Wannamaker P E eds. Three-Dimensional Electromagnetics. Amsterdam: Elsevier, Inc., 3-19. Fu C M, Di Q Y, An Z G. 2013. Application of the CSAMT method to groundwater exploration in a metropolitan environment.Geophysics, 78(5): B201-B209. Grayver A V, Streich R, Ritter O. 2014. 3D inversion and resolution analysis of land-based CSEM data from the Ketzin CO2storage formation.Geophysics, 79(2): E101-E114.

Han B, Hu X Y, Schultz A, et al. 2015. Three-dimensional forward modeling of the marine controlled-source electromagnetic field with complex source geometries.ChineseJ.Geophys. (in Chinese), 58(3): 1059-1071, doi: 10.6038/cjg20150330. He J S. 2010. Wide field electromagnetic sounding methods.JournalofCentralSouthUniversity(ScienceandTechnology) (in Chinese), 41(3): 1065-1072. Hu X Y, Peng R H, Wu G J, et al. 2013. Mineral Exploration using CSAMT data: Application to Longmen region metallogenic belt, Guangdong Province, China.Geophysics, 78(3): B111-B119. Jahandari H, Farquharson C G. 2014. A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids.Geophysics, 79(6): E287-E302. Jin J M. 2002. The Finite Element Method in Electromagnetics. 2nd ed. New York: John Wiley & Sons, Inc.

Li J H, Zhu Z Q, Liu S C, et al. 2011. 3D numerical simulation for the transient electromagnetic field excited by the central loop based on the vector finite-element method.JournalofGeophysicsandEngineering, 8(4): 560-567.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.3038/cjg20150331.

Li Y G, Duan S M. 2014. Data processing of marine controlled-source electromagnetic data.PeriodicalofOceanUniversityofChina(in Chinese), 44(10): 106-112.

Li Y G, Key K. 2007. 2D marine controlled-source electromagnetic modeling: Part 1 — An adaptive finite-element algorithm.Geophysics, 72(2): WA51-WA62.

Mitsuhata Y. 2000. 2D electromagnetic modeling by finite-element method with a dipole source and topography.Geophysics, 65(2): 465-475. Mitsuhata Y, Uchida T, Amano H. 2002. 2.5-D inversion of frequency-domain electromagnetic data generated by a grounded-wire source.Geophysics, 67(6): 1753-1768. Sasaki Y, Meju M A. 2009. Useful characteristics of shallow and deep marine CSEM responses inferred from 3D finite-difference modeling.Geophysics, 74(5): F67-F76.

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. Si H. 2013. TetGen, A quality tetrahedral mesh generator and 3D Delaunay triangulator, Version 1.5 User′s Manual. Technical Report 13, Weierstrass Institute for Applied Analysis and Stochastics.

Si H. 2015. TetGen, a Delaunay-based quality tetrahedral mesh generator.ACMTransactionsonMathematicalSoftware, 41(2), doi: 10.1145/2629697.

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.

Unsworth M J, Travis B J, Chave A D. 1993. Electromagnetic induction by a finite electric dipole source over a 2-D earth.Geophysics, 58(2): 198-214.

Wang R, Wang M Y, Di Q Y. 2006. Electromagnetic modeling due to line source in frequency domain using finite element method.ChineseJ.Geophys. (in Chinese), 49(6): 1858-1866, doi: 10.3321/j.issn:0001-5733.2006.06.035.

Ward S H, Hohmann G W. 1988. Electromagnetic theory for geophysical applications. ∥ Nabighian M N ed. Electromagnetic Methods in Applied Geophysics: Volume 1, Theory. Tulsa, OK: Society of Exploration Geophysicists, 167-183. 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. 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. Xue G Q, Yan S, Chen W Y. 2014. Research prospect to grounded-wire TEM with short-offset.ProgressinGeophysics(in Chinese), 29(1): 177-181, doi: 10.6038/pg20140124.

Yan S, Chen M S. 1999a. Line source frequency electromagnetic sounding two-dimensional forward modeling I.CoalGeology&Exploration(in Chinese), 27(5): 60-62.

Yan S, Chen M S. 1999b. Line source frequency electromagnetic sounding two-dimensional forward modeling II.CoalGeology&Exploration(in Chinese), 27(6): 56-59. Yan S, Chen M S. 2000. Finite element solution of three-dimensional geoelectric models in frequency electromagnetic sounding excited by a horizontal electric dipole.CoalGeology&Exploration(in Chinese), 28(3): 50-56. Zhang J F, Tang J T, Yu Y, et al. 2009. Three dimensional controlled source electromagnetic numerical simulation based on electric field vector wave equation using finite element method.ChineseJ.Geophys. (in Chinese), 52(12): 3132-3141, doi: 10.3969/j.issn.0001-5733.2009.12.023. Zhdanov M S, Endo M, Cox L H, et al. 2014. Three-dimensional inversion of towed streamer electromagnetic data.GeophysicalProspecting, 62(3): 552-572.

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

陳小斌, 胡文寶. 2002. 有限元直接迭代算法及其在線(xiàn)源頻率域電磁響應(yīng)計(jì)算中的應(yīng)用. 地球物理學(xué)報(bào), 45(1): 119-130, doi: 10.3321/j.issn:0001-5733.2002.01.015.

底青云, Unsworth M, 王妙月. 2004. 復(fù)雜介質(zhì)有限元法2.5維可控源音頻大地電磁法數(shù)值模擬. 地球物理學(xué)報(bào), 47(4): 723-730, doi: 10.3321/j.issn:0001-5733.2004.04.026.

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

何繼善. 2010. 廣域電磁測(cè)深法研究. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 41(3): 1065-1072.

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

李予國(guó), 段雙敏. 2014. 海洋可控源電磁數(shù)據(jù)預(yù)處理方法研究. 中國(guó)海洋大學(xué)學(xué)報(bào), 44(10): 106-112.

王若, 王妙月, 底青云. 2006. 頻率域線(xiàn)源大地電磁法有限元正演模擬. 地球物理學(xué)報(bào), 49(6): 1858-1866, doi: 10.3321/j.issn:0001-5733.2006.06.035.

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

2010.08.019.

薛國(guó)強(qiáng), 閆述, 陳衛(wèi)營(yíng). 2014. 接地源短偏移瞬變電磁法研究展望. 地球物理學(xué)進(jìn)展, 29(1): 177-181, doi: 10.6038/pg20140124.

閻述, 陳明生. 1999a. 線(xiàn)源頻率電磁測(cè)深二維正演(一). 煤田地質(zhì)與勘探, 27(5): 60-62.

閻述, 陳明生. 1999b. 線(xiàn)源頻率電磁測(cè)深二維正演(二). 煤田地質(zhì)與勘探, 27(6): 56-59.

閻述, 陳明生. 2000. 電偶源頻率電磁測(cè)深三維地電模型有限元正演. 煤田地質(zhì)與勘探, 28(3): 50-56.

張繼鋒, 湯井田, 喻言等. 2009. 基于電場(chǎng)矢量波動(dòng)方程的三維可控源電磁法有限單元法數(shù)值模擬. 地球物理學(xué)報(bào), 52(12): 3132-3141, doi: 10.3969/j.issn.0001-5733.2009.12.023.

(本文編輯何燕)

A vector finite element solver of three-dimensional modelling for a long grounded wire source based on total electric field

LI Jian-Hui1,2,3, Colin G. Farquharson2, HU Xiang-Yun1*, ZENG Si-Hong1

1HubeiSubsurfaceMulti-scaleImagingKeyLaboratory,InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China2DepartmentofEarthSciences,MemorialUniversityofNewfoundland,St.John′s,NL,A1B3X5,Canada3StateKeyLaboratoryforGeomechanicsandDeepUndergroundEngineering,ChinaUniversityofMining&Technology,JiangsuXuzhou221116,China

AbstractDue to the easiness of removing source singularity, the primary/secondary field algorithm is one of popular algorithms for three-dimensional modelling in geophysical electromagnetic methods. However, it is cumbersome to implement this algorithm if a transmitting source is laid on a rugged surface of the Earth or background models are complex. The total field algorithm is an alternative scheme to numerically simulate the electromagnetics, and its difficulty is how to enforce the source item in the Helmholtz equation of electric field or its magnetic vector potential, especially observation points are located near the transmitting source. In controlled-source electromagnetic methods, transmitting sources including a long grounded wire can be viewed as a combination of many horizontal electric dipoles (HEDs). In our three-dimensional scheme, the source item in the Helmholtz equation of total electric field, the governing equation for vector finite element (FE) method, could be dealt with in the form of HED.

KeywordsLong grounded wire； Three-dimensional modelling； Vector finite element； Total electric field； Unstructured tetrahedrons

基金項(xiàng)目國(guó)家自然科學(xué)基金項(xiàng)目(41274077，41474055，41504088)、中國(guó)礦業(yè)大學(xué)深部巖土力學(xué)與地下工程國(guó)家重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金(SKLGDUEK1312)和同濟(jì)大學(xué)海洋地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金(MGK1405)聯(lián)合資助.

作者簡(jiǎn)介李建慧, 男, 1982年生, 講師, 從事電磁法數(shù)值計(jì)算與資料處理研究. E-mail: ljh_geophysics@163.com *通訊作者胡祥云, 男, 1966年生, 教授, 博士生導(dǎo)師, 從事電磁法理論及其應(yīng)用方面研究. E-mail: xyhu@163.com

doi:10.6038/cjg20160432 中圖分類(lèi)號(hào)P631

收稿日期2015-04-24，2015-08-21收修定稿

李建慧, Farquharson C G, 胡祥云等. 2016. 基于電場(chǎng)總場(chǎng)矢量有限元法的接地長(zhǎng)導(dǎo)線(xiàn)源三維正演.地球物理學(xué)報(bào)，59(4):1521-1534,doi:10.6038/cjg20160432.

Li J H, Farquharson C G, Hu X Y, et al. 2016. A vector finite element solver of three-dimensional modelling for a long grounded wire source based on total electric field.ChineseJ.Geophys. (in Chinese),59(4):1521-1534,doi:10.6038/cjg20160432.