海洋可控源電磁三維非結(jié)構(gòu)矢量有限元數(shù)值模擬

楊軍， 劉穎， 吳小平*

1 中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院 地震與地球內(nèi)部物理實驗室， 合肥 230026 2 中國海洋大學(xué)海洋地球科學(xué)學(xué)院， 青島 266100

?

海洋可控源電磁三維非結(jié)構(gòu)矢量有限元數(shù)值模擬

楊軍1， 劉穎2， 吳小平1*

1 中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院 地震與地球內(nèi)部物理實驗室， 合肥 230026 2 中國海洋大學(xué)海洋地球科學(xué)學(xué)院， 青島 266100

本文實現(xiàn)了海洋可控源電磁三維矢量有限元數(shù)值模擬.由于采用非結(jié)構(gòu)四面體單元進行三維網(wǎng)格剖分，該方法可以模擬復(fù)雜電性異常體和海底地形.一維模型的數(shù)值模擬結(jié)果表明，電場實、虛部均與解析解吻合得相當(dāng)好，計算誤差基本小于1%.二維模型的計算結(jié)果與已有的二維自適應(yīng)非結(jié)構(gòu)有限元模擬結(jié)果吻合很好.帶地形模型的數(shù)值模擬結(jié)果顯示，海底地形對電場影響很大，有可能掩蓋海底油氣藏產(chǎn)生的異常.

海洋可控源電磁; 三維正演; 矢量有限元; 非結(jié)構(gòu)網(wǎng)格

1 引言

電磁法是近地表地球物理勘探中的一種重要方法，已被廣泛應(yīng)用于環(huán)境監(jiān)測、工程勘察、固體礦產(chǎn)和油氣勘探等領(lǐng)域.進入21世紀，能源的緊缺使海洋油氣勘探、開發(fā)成了各國關(guān)注的焦點.海洋地震勘探可以獲得高分辨率的地下結(jié)構(gòu)，指示潛在的海底油氣圈閉構(gòu)造，但地震剖面無法直接區(qū)分圈閉構(gòu)造含油還是含水(da Silva et al., 2012).海洋可控源電磁(CSEM)勘探通過人工發(fā)射源在海底激發(fā)和測量低頻電磁波可以獲得海底地層的電阻率信息，從而識別高阻油氣藏，降低勘探風(fēng)險，提高鉆井成功率.實測資料表明，有明顯海洋CSEM異常的遠景區(qū)，鉆井成功率達70%，反之則鉆井成功率將至35%以下(Hesthammer et al., 2010).近年來，ExxonMobile、Shell等國際石油公司和EMGS等海洋電磁勘探技術(shù)服務(wù)公司已在世界各海域開展了大量海洋電磁勘探工作(何展翔和余剛, 2008).海洋CSEM目前的數(shù)據(jù)處理主要還是通過某一頻率實測電磁場的振幅隨偏移距變化曲線(MVO)和相位隨偏移距變化曲線(PVO)做定性解釋(Constable and Srnka, 2007).要充分挖掘海洋CSEM電磁觀測數(shù)據(jù)中包含的海底物性信息，對復(fù)雜海底構(gòu)造的電性參數(shù)做定量分析，海洋CSEM三維反演是必不可少的，而要實現(xiàn)反演，首先要開發(fā)可靠的三維正演模擬算法.

現(xiàn)有的地球電磁三維正演模擬算法主要基于積分方程法 (Hohmann, 1975; Wannamaker, 1991; Xiong, 1992; Avdeev et al., 2002; Zhdanov et al., 2006; Avdeev and Knizhnik, 2009; 殷長春和樸化榮, 1994; 陳桂波等, 2009; Ren et al., 2013a)、有限差分法(Newman and Alumbaugh, 1995; Weiss and Constable, 2006; Maa?, 2007; Streich, 2009; 岳建華等, 2008)、有限體積法(Haber and Ascher, 2001; 楊波等, 2012)、邊界單元法(Ren et al., 2013a)和有限單元法(Pridmore et al., 1981; Badea et al., 2001; da Silva et al., 2012; Puzyrev et al., 2013; Sommer et al., 2013; Weiss, 2013; Koldan et al., 2014; 付長民等, 2009; 徐志鋒和吳小平, 2010; 李勇等, 2015)等方法，Avdeev(2005)和B?rner(2010)對近年來地球電磁場數(shù)值模擬方法發(fā)展做了較為全面的綜述，眾多的數(shù)值計算方法在不同的物理問題中各有優(yōu)劣，不能一概而論.本文以下主要討論有限單元法在電磁三維數(shù)值模擬中的應(yīng)用.

電磁場有限單元模擬可直接由麥克斯韋方程組求解電磁場(Schwarzbach et al., 2011; da Silva et al., 2012; Um et al., 2013; Grayver and Bürg, 2014)，也可以先求解庫倫規(guī)范或洛倫茲規(guī)范下的電磁勢，再求電磁場(Badea et al., 2001; Mukherjee and Everett, 2011; Puzyrev et al., 2013; Weiss, 2013; Koldan et al., 2014; 徐志鋒和吳小平, 2010).求解電磁勢避免了地球介質(zhì)中電場的法向不連續(xù)問題，可以直接利用節(jié)點有限元方法進行求解，但由于節(jié)點有限元的插值函數(shù)不滿足散度為零的條件，會產(chǎn)生偽解.在變分中添加罰項可以在一定程度上緩解偽解問題，但并不能完全消除，而且可能降低精度(王長清, 2005)，先求解電磁勢再通過差分格式求解電磁場時顯然也會降低數(shù)值模擬的精度(Badea et al., 2001).另一種思路是直接利用麥克斯韋方程組求解電場和磁場，若采用常規(guī)的節(jié)點有限元進行求解，會遇到三個問題：(1)電場不滿足法向連續(xù)條件(徐志鋒和吳小平, 2010)；(2)電磁場散度為零的條件得不到保證；(3)方程中旋度算子的處理和邊界條件的施加較為麻煩(Jin, 2002).

為了解決以上問題，通常采用矢量有限元(或稱

邊有限元)方法進行求解.McMahon(1952)最早開始采用棱邊單元插值形函數(shù)，但直到Nédélec(1980)發(fā)表關(guān)于矢量有限元的理論文章，并由Bossavit和Vérité(1983)成功運用于渦旋電流的計算，這種類型的插值形函數(shù)才得以受到關(guān)注并廣泛應(yīng)用于計算電磁學(xué)中.矢量有限元通過將待求物理量置于離散單元的邊上，只要求該物理量切向連續(xù)，解決了電場法向不連續(xù)問題.同時，矢量有限元的形函數(shù)自動滿足散度為零的條件，不需要在變分問題中強加散度校正項，因此矢量有限元方法在地球電磁數(shù)值模擬中備受關(guān)注.Schwarzbach等(2011)在有限元庫FEMSTER的基礎(chǔ)上發(fā)展了基于二次場的海洋可控源電磁正演.Ren等(2013b, 2014b)實現(xiàn)了三維平面電磁波場數(shù)值模擬的自適應(yīng)網(wǎng)格矢量有限元計算和域分解矢量有限元并行計算.Um等(2013)采用預(yù)條件迭代解法的矢量有限元方法研究了電磁場在地球介質(zhì)中的擴散. Ren等(2014a)在三維平面電磁波場正演中采用了混合邊界元-有限元方法，并利用快速多極子算法實現(xiàn)了加速.Grayver和Bürg(2014)通過將復(fù)數(shù)形式的雙旋度電場方程分解為兩個實型方程實現(xiàn)了基于六面體網(wǎng)格的有限元方法穩(wěn)定求解.近年來國內(nèi)地球電磁數(shù)值模擬中也有部分矢量有限元方面的研究.王燁(2008)研究了基于矢量有限元的高頻大地電磁三維數(shù)值模擬，孫向陽等(2008)利用矢量有限元模擬了隨鉆測井儀在傾斜各向異性地層中的電磁響應(yīng)，劉長生等(2010)實現(xiàn)了自適應(yīng)矢量有限元的三維大地電磁正演模擬，Li等(2011)完成了中心回線激發(fā)的瞬變電磁場的矢量有限元三維正演模擬.在海洋CSEM三維數(shù)值模擬方面，采用較多的仍是有限差分法和節(jié)點有限元方法，基于矢量有限元的三維數(shù)值模擬還未見相關(guān)文章發(fā)表.

本文從麥克斯韋方程組出發(fā)，導(dǎo)出關(guān)于電場的雙旋度方程，利用廣義變分原理求得對應(yīng)的泛函，再利用矢量有限元方法在非結(jié)構(gòu)四面體網(wǎng)格中進行求解，實現(xiàn)了非均勻各向同性介質(zhì)中電場的求解.程序中電偶極子源的姿態(tài)可以任意設(shè)定.與一維解析解及二維有限元計算結(jié)果的對比顯示了我們的算法具有較高的計算精度.同時，采用非結(jié)構(gòu)四面體單元網(wǎng)格進行空間離散，可以模擬復(fù)雜的三維海底地形和電性結(jié)構(gòu)，便于研究海底地形對實際海洋CSEM勘探的影響.

2 基本原理

2.1 邊值問題

時間域中的電磁場可由以下麥克斯韋方程組描述，

(1a)

(1b)

(1c)

(1d)

式中E和H分別為電場和磁場強度，σ為介質(zhì)電導(dǎo)率，ε為介電常數(shù)(電容率)，μ為磁導(dǎo)率，Js是外加電流源.通常我們假設(shè)地球介質(zhì)中不存在自由電荷，即ρ=0.式(1a)—(1d)中只有兩式是獨立的，對于頻率域問題，我們以(1a)、(1b)式為控制方程，取時諧因子eiω t，可推導(dǎo)出關(guān)于電場強度E的一個雙旋度方程(NewmanandAlumbaugh, 1995):

(2)

其中μ為相對磁導(dǎo)率，k2=iωμ0(σ+iωε)，μ0=4π×10-7為真空中的磁導(dǎo)率.

為了簡化變分問題，通常采用理想導(dǎo)體(PEC)邊界條件(Schwarzbachetal., 2011;Grayveretal., 2013)，即n×E=0，為了達到較理想的精度，需要取非常大的計算域.這里我們采用吸收邊界條件，通過構(gòu)造一種m階微分算子，將邊界上電磁場的前m項抵消或湮滅，本文采用一階算子，即索末菲輻射算子，或稱為索末菲輻射邊界條件(王長清, 2005):

(3)

式(2)—(3)就組成了關(guān)于電場E的邊值問題.該邊值問題適用于非均勻電阻率、介電常數(shù)和磁導(dǎo)率的各向同性介質(zhì).

有限元方法求解邊值問題的思路是將邊值問題轉(zhuǎn)換成對應(yīng)的泛函極值問題，通過將計算域離散成有限個單元，由基函數(shù)插值得到一組有限元方程，再進行求解.由廣義變分原理，可以得到與邊值問題(2)—(3)對應(yīng)的泛函極值問題(Jin, 2002):

(4)

2.2 矢量有限元求解

本文采用非結(jié)構(gòu)四面體單元網(wǎng)格進行空間離散，選用Nédélec型(也稱Whitney型)矢量形函數(shù)(或稱基函數(shù))(Nédélec, 1980; Jin, 2002):

(5)

其中l(wèi)i為第i條邊的長度，Li1、Li2分別為第i條邊的i1和i2節(jié)點的節(jié)點形函數(shù).二維三角形單元和三維四面體單元中邊的編號規(guī)則如圖1所示，則離散單元內(nèi)任一點的電場值可由以下插值得到，

(6)

式中n=6，為單元上邊的條數(shù)，Ei為待求系數(shù).

圖1 四面體單元中邊的編號規(guī)則Fig.1 Edge definition for tetrahedral elements

(7)

其中

E6×1=(E1,E2,E3,E4,E5,E6)T,

(8)

(9)

(10)

(11)

(12)

這里我們考慮電偶極子源，即

Js=Idlδ(r-rs),

(13)

rs為電偶極子源所在的空間位置.

對(7)式取一階變分即得有限元方程組

KE=P,

(14)

其中K=K1+K2+K3.若令K3=0，則轉(zhuǎn)化為PEC邊界條件下的電場邊值問題的求解.對于線性方程組(14)，我們利用PARDISO求解器(SchenkandG?rtner, 2004)進行求解，系數(shù)矩陣采用CSR格式進行壓縮存儲.

3 數(shù)值模擬結(jié)果

3.1 一維層狀模型

為了驗證有限元算法的正確性并分析計算精度，我們設(shè)計了M1模型(空氣-海水-沉積層模型)，其幾何與物理參數(shù)如圖2a所示.由于偶極子源產(chǎn)生的場具有一定的對稱性，我們將有限元算法的計算域設(shè)計為圓柱體狀.相比常規(guī)的六面體，在相同半徑下，采用圓柱體可以減小計算域，從而減少網(wǎng)格節(jié)點，提高求解效率.偶極子源置于距海底50 m高處，測線置于海底，沿x方向布置，測點間距為50 m.圖2b是利用開源網(wǎng)格生成程序Gmsh(Geuzaine and Remacle, 2009)生成的網(wǎng)格，生成網(wǎng)格時在源和測線附近進行了局部加密處理，最終M1模型被離散成700704個四面體，共118191個節(jié)點，821044條邊.計算中采用了三種裝置：(1)inline裝置，即電偶極子源與測點共線，且偶極子方向與測線方向一致；(2)broadside裝置，即電偶極子源與測點共線，且偶極子方向與測線方向垂直；(3)電偶極子源與測點共線，且偶極子方向與測線方向成45°夾角.我們分別計算了發(fā)射頻率為0.01 Hz、0.1 Hz和1.0 Hz的偶極子源產(chǎn)生的電場.分配8個主頻為2.4 GHz的處理器做并行計算時，單個頻點的計算時間約為2.0 min.圖3給出了不同裝置、不同發(fā)射頻率的電場數(shù)值解與解析解的對比.

在M1模型的基礎(chǔ)上，我們設(shè)計了如圖4所示的二維模型M2.相比M1，M2多了一個電阻率為10 Ωm的基底層，基底層與沉積層之間是一個起伏界面.左側(cè)界面距離海底1.5 km，右側(cè)界面距離海底1.0 km，中部存在一個寬度為1.5 km的地壘狀突起.這個二維模型中電偶極子源產(chǎn)生的電磁場無法得到精確的解析解，我們選擇了將矢量有限元的數(shù)值模擬結(jié)果與國際上知名的海洋電磁二維自適應(yīng)非結(jié)構(gòu)有限單元程序MARE2DCSEM(Li and Key, 2007)的計算結(jié)果進行對比，以驗證我們程序的正確

圖2 (a) M1模型簡圖； (b) M1的離散網(wǎng)格Fig.2 (a) Sketch of the M1 model; (b) discretization of the M1 model

圖3 不同裝置、不同發(fā)射頻率電場數(shù)值解與解析解對比. 圖中實線和虛線分別代表電場解析解的實虛部，圓圈和三角形代表電場三維數(shù)值解的實虛部，點虛線為相對誤差曲線.Fig.3 Comparisons between the numerical and analytical solutions. Solid and dashed lines represent the real and imaginary parts of the analytical electrical fields, respectively. Circles and triangles stand for the real and imaginary parts of the numerical solutions. Dash-dotted lines are the relative errors.

圖4 (a)M2模型簡圖；(b)MARE2DCSEM采用的二維網(wǎng)格，由上到下分別為空氣層(紅色)、海水層(綠色)、沉積層(天藍色)和基底(深藍色)；(c)矢量有限元正演采用的三維網(wǎng)格Fig.4 (a) Sketch of the M2 model; (b) The triangular meshes applied in 2D FE algorithm MARE2DCSEM.Layers from top to bottom are air (red), sea (green), sediment (sky blue) and basement (deep blue); (c) Tetrahedral meshes applied in our 3D simulation

圖5 三維矢量有限元與二維自適應(yīng)有限元(MARE2DCSEM)計算結(jié)果對比.圖中實線和虛線分別代表電場解析解的實虛部，圓圈和三角形代表電場三維數(shù)值解的實虛部Fig.5 Comparisons between the 3D vector FE and MARE2DCSEM solutions.Solid and dashed lines represent the real and imaginary parts of the analytical electrical fields, respectively. Circles and triangles stand for the real and imaginary parts of the numerical solutions

圖6 (a) M3系列模型簡圖； (b) M32模型的網(wǎng)格剖分，由上到下分別為空氣層、海水層和沉積層，沉積層中紫色扁橢球體為油藏模型； (c)M34模型的網(wǎng)格剖分Fig.6 (a) Sketchs of the M3 model series; (b) Meshes of the M32 model.Layers from top to bottom are air (yellow), sea (blue) and sediment (green) and reservoir (slim purple ellipsoid); (c) Meshes of the M34 model

圖7 Inline裝置和broadside裝置的電場MVO曲線Fig.7 MVO curves with inline and broadside geometries

圖8 (a) 帶三維海底地形的海底油藏模型，圖中淺藍色平面為海平面，橙色界面為海底，紫色界面為沉積層與基地交界面；(b) 海底地形及電偶極子源和測線的位置；(c) 海底地形網(wǎng)格；(d) 沉積層與基底界面地形網(wǎng)格Fig.8 (a) A reservoir model with 3D submarine and interface topographies. Interfaces from top to bottom are sea level (light blue), submarine (orange) and sediment-basement interface (purple); (b) Submarine topography and the survey geometry; (c) Meshes of the submarine topography; (d) Meshes of the interface topography between the sediment and basement

圖9 Inline和broadside裝置對M41和M42模型的MVO響應(yīng)曲線和測線所在剖面的地形Fig.9 MVO responses to Model M41 and M42 with inline and broadside geometries. Topography along the survey line is also shown in the figures

性.如圖5所示，無論是inline裝置還是broadside裝置，我們的計算結(jié)果都和MARE2DCSEM的吻合得很好.3.3 海底地形對電場的影響分析

實際海底地形條件遠比前面討論的模型復(fù)雜，為此有必要分析海底地形對海洋CSEM勘探實測數(shù)據(jù)的影響.為此，我們設(shè)計了如圖6所示的M3系列模型，其中M31為無油藏的水平層狀模型，海水深度1000 m，M32在M31基礎(chǔ)上加入了油藏模型，M33在M31基礎(chǔ)上加入了一個半橢球狀的凸起地形，M34則在M33的基礎(chǔ)上添加了油藏模型.位于沉積層內(nèi)的油藏利用一個高阻的扁平橢球體模擬，中心坐標為(4000，0，-2000)，長半軸a=2000 m，沿x方向，b=1000 m，沿y方向，c=100 m，沿z方向.M33和M34模型中的半橢球形凸起地形位于油藏正上方，長半軸a=2000 m，沿y方向，b=1000 m，沿x方向，c=200 m，沿z方向.電偶極子源置于(-1000,0，-950)處，根據(jù)油藏埋深，我們選擇了0.75 Hz和1.0 Hz兩個發(fā)射頻率，分別計算了inline和broadside裝置的電場.位于海底的測線沿x軸布置，范圍在x=1000 m到x=7000 m之間，點距100 m.圖7所示為數(shù)值模擬結(jié)果，從圖7(a—d)M31和M32曲線可以看到，inline模式下水平層狀模型中油藏引起了明顯MVO異常，而broadside模式下卻未能獲得明顯的異常，這主要是因為油藏的長半軸沿x方向，因此y方向的異常要小于x方向.M31和M33曲線的差異表明地形對電場的影響很大，M33和M34曲線的對比表明地形的存在會削弱甚至掩蓋油藏產(chǎn)生的MVO異常，因此當(dāng)采用MVO曲線對潛在油藏做定性分析時，必須考慮海底地形的影響.

3.4 三維復(fù)雜海底地形下的油藏模型

為了分析海底地形對海洋CSEM勘探的影響，我們設(shè)計了如圖8a所示的帶海底地形的油藏模型.海水深度在931～1479 m之間，海底存在一個走向沿y軸的海嶺，最大高差約為550 m，測線垂直海嶺走向布置，長度為6 km(圖8b).沉積層和基底層的界面也存在地形起伏，其深度在2509～4088 m之間.位于沉積層內(nèi)的油藏同樣采用一個高阻的扁平橢球體模擬，其幾何和物理參數(shù)與M3系列模型中的一樣.為了提高數(shù)值模擬精度，在生成非結(jié)構(gòu)四面體網(wǎng)格時對兩個起伏界面的所有高程控制點均進行了加密處理，測線所在剖面也進行了網(wǎng)格加密，詳見圖8c、8d.為了分析地形的影響和MVO曲線對海底油藏的響應(yīng)，我們分別對不存在海底油藏(M41)和存在海底油藏(M42)兩種情況進行了數(shù)值模擬.我們選取了0.25 Hz、0.50 Hz和0.70 Hz三個發(fā)射頻率，計算結(jié)果如圖9所示.

4 結(jié)論與討論

本文實現(xiàn)了基于非結(jié)構(gòu)網(wǎng)格的海洋CSEM三維矢量有限元正演模擬.與常規(guī)的節(jié)點有限元方法相比，矢量有限元避免了電場法向不連續(xù)帶來的問題，也消除了節(jié)點有限元插值基函數(shù)散度不為零帶來的偽解問題，可以直接從電場的雙旋度矢量方程出發(fā)求解電場.通過與一維模型的解析解及已有二維有限元計算結(jié)果對比，驗證了我們的算法具有較高的計算精度.由于采用非結(jié)構(gòu)四面體單元網(wǎng)格，我們可以模擬復(fù)雜海底地形和異常體結(jié)構(gòu)的海洋CSEM響應(yīng)，計算結(jié)果表明海底地形對電場影響很大，有可能掩蓋深部油藏產(chǎn)生的MVO異常.

海洋可控源電磁數(shù)值模擬的最終目的是服務(wù)于海洋地質(zhì)調(diào)查和油氣勘探，而實際的海洋地質(zhì)環(huán)境非常復(fù)雜，要求采用的數(shù)值計算方法在模型網(wǎng)格和算法實現(xiàn)上都要有較大的靈活性.有限差分方法的優(yōu)點是簡單、容易實現(xiàn)，但因其要求結(jié)構(gòu)化的正交網(wǎng)格，模擬復(fù)雜地形和地質(zhì)構(gòu)造存在較大困難.而且正交網(wǎng)格不易實施局部加密處理，若對局部區(qū)域進行加密，易導(dǎo)致邊界單元質(zhì)量較差，影響計算精度.為了實現(xiàn)正交網(wǎng)格的局部加密，國際上也提出了一些新的方法，其中一個是存在懸掛點正交網(wǎng)格(Bangerth et al., 2007; Bangerth and Kayser-Herold, 2009)，但其在算法實現(xiàn)上還是較為復(fù)雜.近來在地震波場數(shù)值模擬中為了應(yīng)對復(fù)雜地形和地質(zhì)構(gòu)造發(fā)展了基于曲邊網(wǎng)格的有限差分算法(Zhang et al., 2012; Li et al., 2015)，但地球電磁場正演模擬中還未見相關(guān)研究.積分方程方法自然滿足電磁場的邊界條件，但實際勘探中地質(zhì)異常體構(gòu)造的復(fù)雜性和分布的不確定性，以及求解背景模型基本解的困難很大程度上限制了其在電磁三維數(shù)值模擬和反演中的應(yīng)用.常規(guī)的節(jié)點有限元方法可以通過求解電磁勢的方式實現(xiàn)復(fù)雜地質(zhì)模型中電磁場的正演計算，但在非結(jié)構(gòu)網(wǎng)格中，需要先利用插值方法獲得差分格點的電磁勢，再通過差分格式得到電磁場分量，前后兩步處理都會帶來計算精度的損失.采用基于非結(jié)構(gòu)四面體網(wǎng)格的矢量有限元方法進行海洋CSEM正演模擬的一大優(yōu)勢是可以方便地模擬復(fù)雜的三維結(jié)構(gòu)和海底地形，并直接獲得電磁場.進一步提高矢量有限元計算精度需采用高階基函數(shù)，但是Nédélec型高階基函數(shù)不能向下兼容，即不同階次的基函數(shù)的形式不同，因此采用高階基函數(shù)的海洋電磁三維矢量有限元模擬還面臨相當(dāng)多的困難，也是下一階段的研究方向.

致謝 感謝中國海洋大學(xué)李予國教授課題組提供二維自適應(yīng)非結(jié)構(gòu)有限單元計算結(jié)果和兩位匿名評審專家非常好的建議.

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.

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

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.

Bangerth W, Hartmann R, Kanschat G. 2007. Deal. II—a general-purpose object-oriented finite element library. ACM Transactions on Mathematical Software (TOMS), 33(4): Article No. 24. Bangerth W, Kayser-Herold O. 2009. Data structures and requirements for hp finite element software. ACM Transactions on Mathematical Software (TOMS), 36(1): Article No. 4. B?rner R U. 2010. Numerical modelling in geo-electromagnetics: advances and challenges.SurveysinGeophysics, 31(2): 225-245.

Bossavit A, Vérité J. 1983. The "TRIFOU" Code: Solving the 3-D eddy-currents problem by using H as state variable.IEEETransactionsonMagnetics, 19(6): 2465-2470.

Chen G B, Wang H N, Yao J J , et al. 2009. Modeling of electromagnetic responses of 3-D electrical anomalous body in a layered anisotropic earth using integral equations.ChineseJournalofGeophysics(in Chinese), 52(8): 2174-2181, doi: 10.3969/j.issn.0001-5733.2009.08.028.

Constable S, Srnka L J. 2007. An introduction to marine controlled-source electromagnetic methods for hydrocarbon exploration.Geophysics, 72(2): WA3-WA12.

da Silva N V, Morgan J V, MacGregor L, et al. 2012. A finite element multifrontal method for 3D CSEM modeling in the frequency domain.Geophysics, 77(2): E101-E115.

Fu C M, Di Q Y, Wang M Y. 2009. 3D numeric simulation of marine controlled source electromagnetics (MCSEM).OilGeophysicalProspecting(in Chinese), 44(3): 358-363.

Geuzaine C, Remacle J F. 2009. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities.InternationalJournalforNumericalMethodsinEngineering, 79(11): 1309-1331.

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.

Grayver A V, Bürg M. 2014. Robust and scalable 3-D geo-electromagnetic modelling approach using the finite element method.GeophysicalJournalInternational, 198(1): 110-125.

Haber E, Ascher U M. 2001. Fast finite volume simulation of 3D electromagnetic problems with highly discontinuous coefficients.SIAMJournalonScientificComputing, 22(6): 1943-1961.He Z X, Yu G. 2008. Marine EM survey technology and its new advances.ProgressinExplorationGeophysics(in Chinese), 31(1): 2-9.

Hesthammer J, Stefatos A, Boulaenko M, et al. 2010. CSEM performance in light of well results.TheLeadingEdge, 29(1): 34-41.

Hohmann G W. 1975. Three-dimensional induced polarization and electromagnetic modeling.Geophysics, 40(2): 309-324.

Jin J M. 2002. The Finite Element Method in Electromagnetics. 2nd ed. New York: Wiley-IEEE Press.

Koldan J, Puzyrev V, de la Puente J, et al. 2014. Algebraic multigrid preconditioning within parallel finite-element solvers for 3-D electromagnetic modelling problems in geophysics.GeophysicalJournalInternational, 197(3): 1442-1458.

Li H, Zhang W, Zhang Z G, et al. 2015. Elastic wave finite-difference simulation using discontinuous curvilinear grid with non-uniform time step: two-dimensional case.GeophysicalJournalInternational, 202(1): 102-118.

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 G, Key K. 2007. 2D marine controlled-source electromagnetic modeling: Part 1—An adaptive finite-element algorithm.Geophysics, 72(2): WA51-WA62.

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.ChineseJournalofGeophysics(in Chinese), 58(3): 1072-1087, doi: 10.6038/cjg20150331.

Liu C S, Tang J T, Ren Z Y, et al. 2010. Three-dimension magnetotellurics modeling by adaptive edge finite-element using unstructured meshes.JournalofCentralSouthUniversity(ScienceandTechnology) (in Chinese), 41(5): 1855-1859.

Maa? F A. 2007. Fast finite-difference time-domain modeling for marine-subsurface electromagnetic problems.Geophysics, 72(2): A19-A23.

McMahon J. 1952. Lower bounds for the electrostatic capacity of a cube. Proceedings of the Royal Irish Academy.SectionA:MathematicalandPhysicalSciences, 55: 133-167. Mukherjee S, Everett M E. 2011. 3D controlled-source electromagnetic edge-based finite element modeling of conductive and permeable heterogeneities.Geophysics, 76(4): F215-F226.Nédélec J C. 1980. Mixed finite elements in R3.NumerischeMathematik, 35(3): 315-341.

Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences.GeophysicalProspecting, 43(8): 1021-1042.

Pridmore D F, Hohmann G W, Ward S H, et al. 1981. An investigation of finite-element modeling for electrical and electromagnetic data in three dimensions.Geophysics, 46(7): 1009-1024.

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.Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2013a. Boundary element solutions for broad-band 3-D geo-electromagnetic problems accelerated by an adaptive multilevel fast multipole method.GeophysicalJournalInternational, 192(2): 473-499.

Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2013b. A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modelling.GeophysicalJournalInternational, 194(2): 700-718.

Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2014a. A hybrid boundary element-finite element approach to modeling plane wave 3D electromagnetic induction responses in the Earth.JournalofComputationalPhysics, 258: 705-717.

Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2014b. A finite-element-based domain-decomposition approach for plane wave 3D electromagnetic modeling.Geophysics, 79(6): E255-E268.

Schenk O, G?rtner K. 2004. Solving unsymmetric sparse systems of linear equations with PARDISO.FutureGenerationComputerSystems, 20(3): 475-487.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.Sommer M, H?lz S, Moorkamp M, et al. 2013. GPU parallelization of a three dimensional marine CSEM code.Computers&Geosciences, 58: 91-99.

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.

Sun X Y, Nie Z P, Zhao Y W, et al. 2008. The electromagnetic modeling of logging-while-drilling tool in tilted anisotropic formations using vector finite element method.ChineseJournalofGeophysics(in Chinese), 51(5): 1600-1607, doi: 10.3321/j.issn:0001-5733.2008.05.036.

Um E S, Commer M, Newman G A. 2013. Efficient pre-conditioned iterative solution strategies for the electromagnetic diffusion in the Earth: finite-element frequency-domain approach.GeophysicalJournalInternational, 193(3): 1460-1473. Wang C Q. 2005. Advanced Computational Electromagnetics (in Chinese). Beijing: Peking University Press.

Wang Y. 2008. A study of 3D high frequency magnetotellurics modeling by edge-based finite element method [Ph. D. thesis] (in Chinese). Changsha: Central South University. 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. Weiss C J. 2013. Project APhiD: A Lorenz-gaugedA-Φdecomposition for parallelized computation of ultra-broadband electromagnetic induction in a fully heterogeneous Earth.Computers&Geosciences, 58: 40-52.

Xiong Z H. 1992. Electromagnetic modeling of 3-D structures by the method of system iteration using integral equations.Geophysics, 57(12): 1556-1561.Xu Z F, Wu X P. 2010. Controlled source electromagnetic 3-D modeling in frequency domain by finite element method.ChineseJournalGeophysics(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.ChineseJournalofGeophysics(in Chinese), 55(4): 1390-1399, doi: 10.6038/j.issn.0001-5733.2012.04.035.

Yin C C, Piao H R. 1994. Three-dimensional electromagnetic modelling and its application in frequency-domain electromagnetic sounding methods.ActaGeophysicaSinica(in Chinese), 37(6): 820-827.

Yue J H, Yang H Y, Hu B. 2008. 3D finite difference time domain numerical simulation for TEM in-mine.ProgressinGeophysics

(in Chinese), 22(6): 1904-1909, doi: 10.3969/j.issn.1004-2903.2007.06.036.

Zhang W, Zhang Z G, Chen X F. 2012. Three-dimensional elastic wave numerical modelling in the presence of surface topography by a collocated-grid finite-difference method on curvilinear grids.GeophysicalJournalInternational, 190(1): 358-378.

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.

附中文參考文獻

陳桂波, 汪宏年, 姚敬金等. 2009. 用積分方程法模擬各向異性地層中三維電性異常體的電磁響應(yīng). 地球物理學(xué)報, 52(8): 2174-2181, doi: 10.3969/j.issn.0001-5733.2009.08.028.

付長民, 底青云, 王妙月. 2009. 海洋可控源電磁法三維數(shù)值模擬. 石油地球物理勘探, 44(3): 358-363.

何展翔, 余剛. 2008. 海洋電磁勘探技術(shù)及新進展. 勘探地球物理進展, 31(1): 2-9.

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

劉長生, 湯井田, 任政勇等. 2010. 基于非結(jié)構(gòu)化網(wǎng)格的三維大地電磁自適應(yīng)矢量有限元模擬. 中南大學(xué)學(xué)報 (自然科學(xué)版), 41(5): 1855-1859.

孫向陽, 聶在平, 趙延文等. 2008. 用矢量有限元方法模擬隨鉆測井儀在傾斜各向異性地層中的電磁響應(yīng). 地球物理學(xué)報, 51(5): 1600-1607, doi: 10.3321/j.issn:0001-5733.2008.05.036.

王長清. 2005. 現(xiàn)代計算電磁學(xué)基礎(chǔ). 北京: 北京大學(xué)出版社.

王燁. 2008. 基于矢量有限元的高頻大地電磁法三維數(shù)值模擬 [博士論文]. 長沙: 中南大學(xué).

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

殷長春, 樸化榮. 1994. 三維電磁模擬技術(shù)及其在頻率測探法中應(yīng)用. 地球物理學(xué)報, 37(6): 820-827.

岳建華, 楊海燕, 胡搏. 2008. 礦井瞬變電磁法三維時域有限差分數(shù)值模擬. 地球物理學(xué)進展, 22(6): 1904-1909, doi: 10.3969/j.issn.1004-2903.2007.06.036.

(本文編輯 胡素芳)

3D simulation of marine CSEM using vector finite element method on unstructured grids

YANG Jun1, LIU Ying2, WU Xiao-Ping1*

1LaboratoryofSeismologyandPhysicsofEarth′sInterior,SchoolofEarthandSpaceSciences,UniversityofScienceandTechnologyofChina,Hefei230026,China2CollegeofMarineGeosciences,OceanUniversityofChina,Qingdao266100,China

Nodal finite element (FE) has become a popular numerical method in computational geo-electromagnetics since 1970 s. However, the zero-divergence property of the electromagnetic field and the discontinuity of the normal component of the electrical field harass the elegant implementation of nodal FE. The former brings in the spurious solution in the simulation while the latter makes it impossible for regular nodal FE to solve electromagnetic field directly. This paper focuses on 3D forward modeling of marine controlled-source electromagnetics (CSEM) with an electrical dipole source, using vector finite element method on unstructured tetrahedral grids which can simulate complicated geometry and submarine topography without spurious solution.The curl-curl electrical field full-wave equation, which is derived from Maxwell′s equations, is used as the governing equation in our marine CSEM simulation. Somerfeld radiation operator is chosen as the boundary condition on the truncated boundary. By applying the generalized variational principle to the boundary value problem of the electrical field, the equivalent variational problem is obtained. Unstructured tetrahedral grids are employed for spatial discretization in the 3D calculation domain using a non-commercial mesh generator Gmsh. Nédélec type vector basis functions are adopted for the interpolation of the electrical field within each tetrahedral element. These basis functions are divergence-free intrinsically and they guarantee the continuity of the tangential field component when crossing the interface between two different media. Unlike the nodal FE, the unknowns are assigned to the edges of the elements rather than the nodes. Substituting the electrical field by the interpolation, a system of linear equations is constructed. We solve the system of linear equations with a direct solver.A 1D stratified air-sea-sediment model is designed for the validation of our 3D vector FE algorithm. We output the numerical solution of the electrical field of 0.01 Hz, 0.1 Hz and 1.0 Hz with three measurement configurations: inline, broadside and a 45-degree angle between the survey line and the dipole source direction. The comparison between the analytical and numerical solution shows that the relative error are less than 1% for most cases when the observation points are a few hundred meters away from the source. The CSEM response of a 2D air-sea-sediment-basement model with a horst type interface between the sediment and basement layers is then compared with the numerical solution of the well-known 2D FE algorithm MARE2DCSEM and it shows high consistency between the 2D and 3D numerical solutions. The influence of submarine topography to CSEM response is also studied. The simulation indicates that submarine topography has significant impact on the electrical field in marine CSEM. The existence of submarine topography may cover the anomalous signal induced by oil reservoir. Finally, the CSEM response to a 3D reservoir model covered by complicated submarine topography and fluctuant stratum interface is simulated and the results show that the complicated topography makes unsmooth changes to the MVO (Magnitude Versus Offset) curve.We develop a 3D marine CSEM simulation algorithm using vector FE on unstructured meshes. The main advantages of vector FE we use in marine CSEM over regular nodal FE are: (1) no divergence correction is needed because of the vector basis functions are divergence-free; (2) the electrical field is solved directly without worrying about the discontinuity of the normal component of the electrical field. On the other hand, unstructured tetrahedral grids can easily handle complicated submarine topography and refined reservoir structures. The comparisons with 1D analytical solution and 2D numerical solution show high accuracy of our algorithm. The simulations also show that submarine topography may be an important contributor to the uncertainties of CSEM MVO anomalies.

Marine controlled-source electromagnetic； 3D simulation； Vector finite element method； Unstructured grids

國家自然科學(xué)基金(41374076，41130420),國家863計劃(2014AA06A610、2012AA09A201),天然氣水合物資源勘查與試采工程國家專項(GZH201400305)和安徽省國土資源科技項目(2013-05)聯(lián)合資助.

楊軍，男，1987年生，博士研究生，主要從事地球物理數(shù)值模擬與反演成像研究.E-mail:youngjun@mail.ustc.edu.cn

*通訊作者吳小平，男，1967年生，教授，博士生導(dǎo)師，研究方向為電磁測深和地球內(nèi)部物理. E-mail:wxp@ustc.edu.cn

10.6038/cjg20150817.

10.6038/cjg20150817

P631

2014-11-07，2015-06-13收修定稿

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

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.