電導(dǎo)率各向異性的海洋電磁三維有限單元法正演

蔡紅柱， 熊彬， Michael Zhdanov

1 College of Mines & Earth Sciences, University of Utah, Salt lake city, UT 84112, USA 2 桂林理工大學(xué)地球科學(xué)學(xué)院， 桂林 541004

?

電導(dǎo)率各向異性的海洋電磁三維有限單元法正演

蔡紅柱1， 熊彬2*， Michael Zhdanov1

1 College of Mines & Earth Sciences, University of Utah, Salt lake city, UT 84112, USA 2 桂林理工大學(xué)地球科學(xué)學(xué)院， 桂林 541004

本文提出了一種基于非結(jié)構(gòu)化網(wǎng)格的海洋電磁有限單元正演算法.為了回避場源奇異性，文中選用二次場算法，將背景電阻率設(shè)置為水平層狀且各向異性，場源在水平層狀各向異性介質(zhì)中所激發(fā)的一次場通過漢克爾積分得到.基于Coulomb規(guī)范得到二次矢量位和標(biāo)量位所滿足的Maxwell方程組，通過Galerkin加權(quán)余量法形成大型稀疏有限元方程，采用不完全LU分解(ILU)預(yù)條件因子的quasi-minimum residual(QMR)迭代解法對有限元方程進(jìn)行求解得到二次矢量位和標(biāo)量位；進(jìn)而，利用滑動平均方法得到二次矢量位和標(biāo)量位在空間的導(dǎo)數(shù)，由此得到二次電磁場；通過一維模型對算法的可靠性進(jìn)行驗(yàn)證，與此同時(shí)，針對實(shí)際復(fù)雜海洋電磁模型，比較有限元模擬結(jié)果與積分方程模擬結(jié)果，進(jìn)一步驗(yàn)證算法精度.若干計(jì)算結(jié)果均表明，文中算法具有良好的通用性，適用于井中電磁、航空電磁，環(huán)境地球物理等非均勻且各向異性介質(zhì)中的電磁感應(yīng)基礎(chǔ)研究.

海洋電磁; 非結(jié)構(gòu)化網(wǎng)格; Coulomb規(guī)范; 三維正演; 有限單元法; 滑動平均

1 引言

陸地可控源電磁法在油氣和礦產(chǎn)資源勘探領(lǐng)域早已得到廣泛應(yīng)用并取得了良好的實(shí)際效果(Ward and Hohmann, 1988；Zhdanov et al., 2006). 近年來，海洋電磁法(MCSEM)也被引入油氣勘探領(lǐng)域以降低勘探風(fēng)險(xiǎn)(Srnka et al.,2006；Constable and Srnka, 2007；Um and Alumbaugh, 2007；Andréis and MacGregor, 2008；Zhdanov, 2010).在海洋環(huán)境中，受海底地形起伏和沉積因素所造成的電阻率不均勻分布和各向異性的影響，地電模型往往會非常復(fù)雜.在這種情況下，三維正演對于能否正確地解釋可控源電磁數(shù)據(jù)顯得至關(guān)重要(da Silva et al.,2012).目前常用的電磁三維正演方法包括有限單元法(FE)、有限差分(FD)和積分方程法(IE).積分方程法只需要對異常體區(qū)域進(jìn)行剖分，因而對于簡單模型中有若干局部異常體非常有效.同時(shí)，積分方程法不需要額外的邊界條件.對于復(fù)雜幾何形狀的異常體，譬如，地形起伏、積分方程法只能用階梯狀網(wǎng)格來近似，且該方法所形成的矩陣為滿陣，當(dāng)異常體區(qū)域較大時(shí)，其模擬效率會明顯降低(Zhdanov et al.,2006).有限差分法相對而言較為簡單，通過有限差分法形成的總體系數(shù)矩陣為稀疏矩陣，其求解效率很高(Newman and Alumbaugh, 1995)，可解決積分方程中模擬大型異常區(qū)域的困難.類似于積分方程法，有限差分法也只能用規(guī)則網(wǎng)格對模擬區(qū)域進(jìn)行剖分，對于復(fù)雜的幾何形狀，有限差分法只能用階梯狀網(wǎng)格來近似.相比基于結(jié)構(gòu)化網(wǎng)格的有限差分和積分方程法，基于非結(jié)構(gòu)化網(wǎng)格的有限單元法更適合于復(fù)雜電阻率模型下的可控源電磁模擬，其形成的總體系數(shù)矩陣的稀疏性和有限差分法相當(dāng)，求解線性方程組的效率很高(Cai et al.,2014).有限單元法求解Maxwell方程可直接用矢量有限單元法形成關(guān)于電磁場本身的大型方程(Cai et al.,2014)，也可以通過解關(guān)于位的方程來間接求取電場和磁場.相對而言，矢量有限元更為復(fù)雜，且形成的有限單元矩陣條件數(shù)更大，用迭代法求解收斂較慢(Jin, 2002).

Badea等(2001)提出了一種基于Coulomb位的有限單元法井中電磁模擬算法.基于Badea所提出的方法，并考慮各向同性的全空間海水背景、各向異性的沉積巖和復(fù)雜地形因素，Puzyrev等(2013)實(shí)現(xiàn)了基于Coulomb電位的海洋電磁三維有限元并行模擬.取海水電導(dǎo)率為背景電阻率固然可簡化問題，一次場也可由解析公式得到(Puzyrev et al.,2013)，然而，更為真實(shí)的背景電阻率模型為層狀介質(zhì).在Badea和Puzeyrev等工作的基礎(chǔ)上，本文提出了一種以水平層狀各向異性電阻率為背景的海洋電磁有限單元法三維正演，計(jì)算中采用猶他大學(xué)電磁正反演協(xié)會(CEMI)的格林庫函數(shù)來計(jì)算一次場，利用COMSOL Multiphysics形成優(yōu)化的非結(jié)構(gòu)化四面體網(wǎng)格，對有限單元方程的總體系數(shù)矩陣做高效的不完全LU分解，并將結(jié)果作為QMR迭代解法的預(yù)條件因子加快收斂速度.通過對不完全LU分解預(yù)條件因子和雅克比預(yù)條件因子加以比較，我們發(fā)現(xiàn)QMR法的收斂效率得到了顯著提高.

在文章的模型研究中，通過對比有限單元法模擬結(jié)果和解析解以及積分方程解，驗(yàn)證了算法的可靠性；同時(shí)，算例中也充分考慮了各向異性和地形起伏對電磁場響應(yīng)的影響.

2 二次電位算法

在地球物理電磁勘探領(lǐng)域，通常采用低頻的電磁信號，忽略位移電流，采用正弦諧變時(shí)間因子e-iω t，Maxwell方程可寫成如下形式：

(1)

(2)

(3)

在海洋環(huán)境中由于沉積因素的影響，通常只考慮橫向各向異性：σx=σy=σh為水平電導(dǎo)率，σz=σv為垂直電導(dǎo)率.

根據(jù)Coulomb位的定義，電磁場E和H可以寫成如下形式：

(4)

(5)

(6)

以及

).

(7)

通過總場算法解方程(6)和(7)一般需要用三維脈沖函數(shù)來近似偶極源中的電流，這種近似的做法通常會引入額外的計(jì)算誤差.選用二次場算法，激發(fā)源中的電流項(xiàng)可以從式(6)和(7)中消除，假定總的電導(dǎo)率張量為背景電導(dǎo)率張量和異常電導(dǎo)率張量之和:

(8)

Coulomb矢量位和標(biāo)量位可類似地寫成一次位和異常位之和：

A=Ap+As，ψ=ψp+ψs.

(9)

將式(9)代入(6)和(7)中可得到如下方程組：

(11)

通過對比式(6)和(7)，方程(10)和(11)的右端項(xiàng)中已不再含有電流項(xiàng)，取而代之的為一次位.因此，要解方程(10)和(11)，須首先計(jì)算層狀各向異性電導(dǎo)率介質(zhì)中的一次矢量位和標(biāo)量位，一次矢量位和標(biāo)量位可通過漢克爾變換求取.注意到公式(10)和(11)中需要求取一次標(biāo)量位的梯度，考慮到如下關(guān)系式

(12)

公式(10)和(11)可寫成如下形式

(13)

(14)

通過公式 (12) 用一次電場代替一次矢量位和標(biāo)量位的梯度之和可避免求取一次標(biāo)量位所引入的數(shù)值誤差.

3 有限單元法分析

將式(13)與(14)沿x,y,z軸依次展開，利用Galerkin法將公式兩端加權(quán)(Jin, 2002)，同時(shí)考慮矢量恒等式和散度定理，可得到如下體積分方程組：

(15)

(16)

(17)

(18)

式中，(u,v)Ω=?ΩuvdΩ，N為有限單元法的形函數(shù).用四面體單元離散整個(gè)計(jì)算區(qū)域，在各個(gè)單元中，電導(dǎo)率為常數(shù)，二次矢量位As及電標(biāo)量位ψs呈線性變化:

(19)

(20)

一次場的各個(gè)分量也可寫成類似公式(19)和(20)的形式.將式(19)和(20)代入到(15)、(16)、(17)和(18)中可得到離散化的線性方程組：

Ku=b,

(21)

假設(shè)有限單元網(wǎng)格有L個(gè)內(nèi)部節(jié)點(diǎn)，K則為一個(gè)4L×4L的方陣(Puzyrev et al.,2013)：

(22)

式中I33為3×3的單位矩陣.方程(21)的右端項(xiàng)可以類似地寫成如下形式：

(23)

為了得到線性方程(21)的唯一解，需要加入合適的邊界條件，文中采用了均勻Dirichlet邊界條件(吳建平等，2004)

(As,ψs)?！?0,0)，

(24)

一旦方程組(21)加上合適的邊界條件，就可以利用適當(dāng)?shù)臄?shù)值方法對其進(jìn)行求解.求解方程(21)的方法總體上可分為直接解法和迭代解法(FreundandNatchtigal, 1991).直接解法的優(yōu)點(diǎn)是對矩陣K進(jìn)行分解，求近似逆和方程右端項(xiàng)無關(guān)，分解的結(jié)果可重復(fù)用于其他右端項(xiàng)，方程右端項(xiàng)與一次場以及場源相關(guān)，因此，這種解法的優(yōu)點(diǎn)是對于有多個(gè)激發(fā)源的正演問題只需將矩陣K做一次分解.直接解法的缺點(diǎn)是對內(nèi)存要求較高且求解速度較慢.基于Krylov子空間的迭代解法對內(nèi)存需求較低，每次迭代僅僅需要進(jìn)行矩陣和數(shù)組相乘等基本運(yùn)算，文中采用QMR方法對方程進(jìn)行求解.相比于傳統(tǒng)的雅克比預(yù)條件因子，這里采用了不完全LU分解的預(yù)條件因子(Um and Newman, 2013).

4 滑動平均求位的導(dǎo)數(shù)

上述有限單元法求解得到的是關(guān)于電磁場的二次位(As,ψs)，為了進(jìn)一步得到物理意義上的電磁場(Es,Hs)各分量，必須借助有效的數(shù)值微分方法來獲取:

(25)

(26)

在規(guī)則化網(wǎng)格中，常規(guī)做法就是在節(jié)點(diǎn)周圍的單元中分別計(jì)算場的導(dǎo)數(shù)，再通過面積加權(quán)平均或體積加權(quán)平均來獲得節(jié)點(diǎn)上場的導(dǎo)數(shù)；或者，直接利用節(jié)點(diǎn)附近沿某一坐標(biāo)軸方向若干點(diǎn)上的場值，通過等距差商的辦法來獲取節(jié)點(diǎn)上場的偏微分值(徐世浙, 1994).

對于完全非結(jié)構(gòu)化網(wǎng)格上的二次位，利用傳統(tǒng)方法得到導(dǎo)數(shù)的誤差很大.通常采用滑動平均的方法來求取二次位的導(dǎo)數(shù)，設(shè)單元中位的線性描述為

fi=axi+byi+czi+d，

(27)

顯然，式中系數(shù)a,b,c即為所需的位沿各坐標(biāo)軸的導(dǎo)數(shù).利用點(diǎn)(xi,yi,zi)附近若干點(diǎn)的場值，形成若干個(gè)類似于(27)的方程，進(jìn)而對這些方程做最小二乘擬合得到系數(shù)a,b,c，即為f在(x,y,z)方向的導(dǎo)數(shù).不同于傳統(tǒng)的最小二乘擬合，這里將這些方程進(jìn)行加權(quán)，離點(diǎn)(xi,yi,zi)距離較近的點(diǎn)將被賦予更高的權(quán)重.文中選取高斯函數(shù)作為加權(quán)函數(shù)，通過加權(quán)最小二乘擬合可得到

(a,b,c)T=(XTWX)-1XTWF，

(28)

其中

X=(x;y;z;I)，

(29)

x=[x1,x2,…,xN]T，

(30)

y=[y1,y2,…,yN]T，

(31)

z=[z1,z2,…,zN]T，

(32)

(33)

h=max(d).

(34)

一旦得到二次位的導(dǎo)數(shù)，就可將結(jié)果代入式(25)和(26)求取二次電場和磁場.

5 數(shù)值算例

為了驗(yàn)證文中所提方法的可靠性和有效性，這里將以四例典型地電模型從不同角度來進(jìn)行分析和闡述.5.1 模型1：具有準(zhǔn)解析解的一維模型

背景電導(dǎo)率為均勻半空間模型，上半空間空氣電導(dǎo)率為10-6S·m-1，下半空間背景電導(dǎo)率0.1S·m-1.從地下400 m到500 m深處為一水平方向無限延伸的高阻體，其電導(dǎo)率為0.01 S·m-1.沿x方向的水平電偶極源位于(0,500,100)m處，其頻率為1 Hz，偶極矩為100 Am.一次場所對應(yīng)的電導(dǎo)率為均勻半空間模型，相對應(yīng)的二次場為水平無限延伸的高阻層所產(chǎn)生的，觀測點(diǎn)位于地表y=0，z=0 處.

對于這個(gè)簡單的層狀模型，可以通過快速漢克爾變換精確地得到近似于解析解的二次場(Ward and Hohmann, 1988；Guptasarma and Singh, 1997).在有限單元法數(shù)值模擬中，將模擬區(qū)域Ω={-4 km≤x,y<4 km, -2 km≤z≤4 km}剖分為非結(jié)構(gòu)化的四面體網(wǎng)格(z軸向下為正).圖1為非結(jié)構(gòu)化四面體網(wǎng)格在y=0處的斷面圖，網(wǎng)格在激發(fā)源和觀測點(diǎn)處有所加密，網(wǎng)格大小可根據(jù)電磁波的趨膚深度來確定，空氣中的網(wǎng)格較為稀疏.這個(gè)模型的單元數(shù)為 1116520，內(nèi)部單元節(jié)點(diǎn)數(shù)為191207，方程自由度為764828，剛度矩陣維數(shù)為764828×764828.圖2為剛度矩陣稀疏特征圖，不難發(fā)現(xiàn)雖然矩陣維數(shù)很大，但是絕大部分的元素值為0，在采取稀疏存儲后剛度矩陣只需要650 M左右的存儲空間.

圖3和圖4分別為二次電場與磁場各個(gè)分量的有限單元數(shù)值解和通過漢克爾變換得到的準(zhǔn)解析解的比較.其中二次電場的x，y，z分量的有限單元法解與解析解的歸一化誤差分別為0.6%, 0.45%, 0.8%.二次磁場的x，y，z分量的有限單元法解與解析解的歸一化誤差分別為0.9%, 0.7%, 1.1%.如前所述，在用QMR算法解有限元方程時(shí)，采用了不完全LU分解的結(jié)果作為預(yù)條件因子.圖5為采用LU預(yù)條件因子的QMR收斂曲線和采用傳統(tǒng)雅克比預(yù)條件因子的QMR收斂曲線的比較.容易發(fā)現(xiàn)，采用雅克比預(yù)條件因子的QMR需要1000次迭代收斂到歸一化余量為10-8， 而采用不完全LU分解預(yù)條件因子的QMR僅僅需要300次迭代歸一化余量就可以收斂到10-10.故而，通過改進(jìn)預(yù)條件因子，迭代解法的收斂特性得到了顯著提高，同時(shí)，不完全LU分解所需要的計(jì)算資源很少，甚至可以忽略不計(jì).

圖1 模型1的非結(jié)構(gòu)化網(wǎng)格在y=0處的斷面示意Fig.1 The vertical section of the unstructured mesh for Model 1 at y=0

圖2 模型1剛度矩陣的稀疏特征Fig.2 The sparsity pattern of the stiffness matrix for Model 1

圖3 模型1電場的有限單元解和解析解比較.子圖a, b, c分別為二次電場的x, y, z分量Fig.3 A compassion of the finite element and analytical solution for the electric field for Model 1. Panel a, b and c show the x, y and z components of the secondary electric field

圖4 模型1磁場的有限單元解和解析解比較.子圖a, b, c分別為二次磁場的x, y, z分量Fig.4 A compassion of the finite element and analytical solution for the magnetic field for Model 1. Panel a, b and c show the x, y and z components of the secondary magnetic field

圖5 模型1的QMR迭代解收斂曲線.紅色實(shí)線為采用雅克比預(yù)條件因子的收斂曲線，藍(lán)色星號曲線為采用不完全LU分解的預(yù)條件因子的收斂曲線Fig.5 QMR convergence plot for Model 1. The solid red curve represents for QMR convergence with Jacobian preconditioner while the blue star indicates the QMR convergence with ILU preconditioner

5.2 模型2：無地形起伏的油藏模型

為簡易起見和便于同積分方程結(jié)果進(jìn)行比較，這里將忽略海底地形.海水深度為1 km，其各向同性電導(dǎo)率為3.3 S·m-1，沉積巖水平方向電導(dǎo)率為1 S·m-1，垂直方向電導(dǎo)率為0.8 S·m-1，油藏所在區(qū)域的水平方向電導(dǎo)率為0.01 S·m-1，垂直方向電導(dǎo)率為0.001 S·m-1.圖 6為模型電導(dǎo)率在x=0和y=0處的切片圖，背景電導(dǎo)率呈層狀的空氣—海水—沉積巖三層分布.有限單元法的計(jì)算模擬區(qū)域設(shè)為Ω={-5 km≤x,y<5 km, -1 km≤z≤4 km}，油藏所在區(qū)域?yàn)棣竌={-1 km≤x,y<1 km, -2 km≤z≤2.1 km}.限于篇幅，略去網(wǎng)格剖分信息.激發(fā)源為一位于(-3000,0,900)m處，沿x方向布設(shè)的水平電偶極源，偶極矩為100 Am，激發(fā)頻率為1 Hz，測線位于y=0，z=1000 m處.

圖7 模型2電導(dǎo)率各向同性情況下的FE解(紅色圓圈)和IE解(藍(lán)色實(shí)線).子圖a，b分別為電場的實(shí)部與虛部Fig.7 A comparison of FE solution (red circle) and IE solution (solid blue) for Model 2 with isotropic conductivity. Panel a and b are the real and imaginary parts of the electric field

圖8 模型2在電導(dǎo)率各向同性情況下y=0處的背景場和總場Fig.8 A comparison of total and background field at y=0 for Model 2 with isotropic conductivity

首先考慮各向同性模型，沉積巖和油藏區(qū)域的電導(dǎo)率都為其所對應(yīng)的水平方向電導(dǎo)率值.之后，再行考慮沉積巖和油藏電導(dǎo)率呈各向異性的情況.對于這樣相對比較簡單的油藏模型，可通過積分方程法直接計(jì)算出二次場(Zhdanov et al.,2006).圖7為模型2在電導(dǎo)率各向同性情況下，二次電場Ex的有限單元數(shù)值解(實(shí)部與虛部)和積分方程法解(實(shí)部與虛部)的比較.Ex分量的有限單元法解與解析解的歸一化誤差為0.5%.圖8為各向同性情況下，Ex的背景場和總場在y=0處的比較，油藏電導(dǎo)率所引起的總場畸變在圖8中清晰可見.

在實(shí)際的海洋環(huán)境中，由于沉積因素影響，電導(dǎo)率通常情況下會呈橫向各向異性，水平電導(dǎo)率往往大于垂向電導(dǎo)率.圖9為沉積層和油藏電導(dǎo)率均呈各向異性下的有限單元法模擬結(jié)果和積分方程模擬結(jié)果的對比，兩者依然非常吻合.Ex分量的有限單元法解與解析解的歸一化誤差為2.3%.圖10則為背景電場Ex和總場在y=0處的比較.通過對比圖7和圖9，容易看到由電導(dǎo)率各向異性所引起二次場的明顯畸變.因此，在實(shí)際的海洋電磁數(shù)據(jù)處理、反演以及解釋中，需要考慮電導(dǎo)率各向異性這一重要地電特征.

圖9 模型2在電導(dǎo)率各向異性情況下的FE解(紅色圓圈)和IE解(藍(lán)色實(shí)線).子圖a，b分別為電場的實(shí)部與虛部Fig.9 A comparison of FE solution (red circle) and IE solution (solid blue) for Model 2 with anisotropic conductivity. Panel a and b are the real and imaginary parts of the electric field

圖10 模型2在電導(dǎo)率各向異性情況下y=0處的背景場和總場Fig.10 A comparison of total and background field at y=0 for Model 2 with anisotropic conductivity

5.3 模型 3：梯形山模型

在海洋電磁勘探中，復(fù)雜的海底地形往往能引起很大的電磁場畸變.Sasaki(2011)提出了通過有限差分模擬地形和進(jìn)行地形校正的算法(Sasaki, 2011).為了驗(yàn)證算法對于地形模擬的可靠性，這里將首先考慮較為簡單的梯形山模型.模型中海水深度為1000 m，電導(dǎo)率為3.3 S·m-1，沉積層和油藏電導(dǎo)率均為各向同性，電導(dǎo)率值分別為1 S·m-1和0.05 S·m-1，背景介質(zhì)設(shè)為層狀的空氣-海水-沉積物三層結(jié)構(gòu).為減小計(jì)算量，文中考慮較小的油藏模型位于Ωa={-500 m≤x,y<500 m, 2000 m≤z≤2100 m}，有限單元法的模擬計(jì)算區(qū)域?yàn)棣?{-5 km≤x,y<5 km, -1 km≤z≤4 km}，海底地形起伏位于Ωt={-500 m≤x,y<500 m}范圍內(nèi)，在{-100 m≤x,y<100 m}范圍內(nèi)，海水深度為900 m，然后線性增加到1000 m.激發(fā)源為一位于(-1500,0,900)處、沿x方向布設(shè)的水平電偶極源，偶極矩為100 Am，激發(fā)頻率為1 Hz.

圖11為梯形山示意圖，其中藍(lán)色部分代表較淺海水，紅色部分代表較深海水.圖12為梯形山模型在y=0處的斷面圖.圖13為模型非結(jié)構(gòu)化網(wǎng)格在y=0處的斷面圖，從圖中可以看出網(wǎng)格在異常區(qū)域以及海水沉積層分界面處均做了加密.

有限單元法模擬的一個(gè)很大的優(yōu)點(diǎn)是可以用非結(jié)構(gòu)化的網(wǎng)格精確地模擬地形.然而，積分方程只能使用結(jié)構(gòu)化、相對比較規(guī)則的網(wǎng)格，通過積分方程法模擬地形需要用規(guī)則化的階梯狀模型來近似地形.圖14為梯形山模型的階梯狀近似，對于這種近似的地形模擬，在測點(diǎn)離地形相對較遠(yuǎn)時(shí)(例如航空電磁)模擬結(jié)果相對光滑可靠；測點(diǎn)離地形越近，由階梯狀近似所引起的離散誤差越大.

為了比較由純地形因素引起的二次場，文中預(yù)設(shè)油藏區(qū)的電導(dǎo)率和沉積層電導(dǎo)率相同，另外，為最

圖11 梯形山三維示意圖Fig.11 A 3D view of the trapezoidal hill model

圖12 梯形山模型在y=0處的斷面圖Fig.12 The vertical section of the trapezoidal model at y=0

圖13 梯形山模型非結(jié)構(gòu)化網(wǎng)格在y=0處的斷面圖Fig.13 The vertical section of the unstructured mesh for the trapezoidal model at y=0

圖14 梯形山地形的階梯狀近似Fig.14 A staircase approximation for the trapezoidal hill model

大限度地減少積分方程中由階梯狀模型近似地形所帶來的誤差，這里將測線放在相對遠(yuǎn)離地形的水平位置z=500 m處.圖15為有限單元法和積分方程法模擬結(jié)果的比較，兩種方法模擬的結(jié)果在遠(yuǎn)離地形的地方的歸一化誤差約為0.8%.實(shí)際的海洋電磁勘探中，接收器一般位于海底，圖16為將測點(diǎn)移到海底后的有限單元法和積分方程模擬結(jié)果的比較，兩種方法的模擬結(jié)果的歸一化誤差約為15%，誤差很可能是由于積分方程中用階梯狀模型近似地形所引起.

接下來，文中同時(shí)考慮地形和油藏對二次場的影響，設(shè)油藏電導(dǎo)率為0.05 S·m-1.圖17為考慮油藏異常電導(dǎo)率之后有限單元法和積分方程法計(jì)算結(jié)果的比較.兩種方法所得到的結(jié)果的歸一化誤差大約為15%.對比圖16和圖17可以得出結(jié)論，相對于油藏，地形對二次場的影響更大.

5.4 模型 4：復(fù)雜海底地形的油藏模型

模型4的各項(xiàng)參數(shù)和模型2中各向異性模型相

圖15 無油藏的梯形山模型的有限單元解(紅色圓圈)和積分方程解(藍(lán)色實(shí)線)在水平剖面y=0， z=500 m處的比較.子圖a，b分別為電場的實(shí)部與虛部Fig.15 A comparison of FE solution (red circle) and IE solution (solid blue) at y=0 and z=500 m for the trapezoidal hill model without reservoir. Panel a and b are the real and imaginary parts of the electric field

圖16 無油藏的梯形山模型的有限單元解(紅色圓圈)和積分方程解(藍(lán)色實(shí)線)沿地形在y=0處剖面的比較.子圖a，b分別為電場的實(shí)部與虛部Fig.16 A comparison of FE solution (red circle) and IE solution (solid blue) at y=0 along the bathymetry for the trapezoidal hill model without reservoir. Panel a and b are the real and imaginary parts of the electric field

圖17 含油藏的梯形山模型的有限單元解(紅色圓圈)和積分方程解(藍(lán)色實(shí)線)沿地形在y=0處剖面的比較.子圖a，b分別為電場的實(shí)部與虛部Fig.17 A comparison of FE solution (red circle) and IE solution (solid blue) at y=0 along the bathymetry for the trapezoidal hill model with reservoir. Panel a and b are the real and imaginary parts of the electric field

圖20 模型4二次場用一次場歸一化后在y=0處的比較圖.紅色圓圈為帶地形的歸一化場，黑色星號為不含地形的歸一化場Fig.20 Secondary field normalized by the background field at y=0. The red circle represents the normalized field with seafloor bathymetry while the black star indicates the normalized field without the seafloor bathymetry

同，相對于模型2的海水-沉積層水平分界面，模型4更為復(fù)雜.圖18為海底地形模型的三維示意圖，海水深度從750 m到1000 m不等，與模型2類似，這里考慮地電背景為層狀的空氣-海水-沉積物三層模型，地形區(qū)域和油藏都被視為異常區(qū).在無油藏情況下總場定義為背景場與由地形引起的二次場的總和，在有油藏和地形的情況下，總場為背景場與由地形和油藏引起的二次場總和.

圖19為Ex的總場和背景場，其中黑色虛線為不含地形影響的總場.通過對比模型4的總場和不含地形的總場曲線，容易看到由地形所產(chǎn)生的電磁場畸變明顯.將二次場振幅用一次場振幅進(jìn)行歸一化后可以得到圖20所示的歸一化場剖面圖，其中紅色圓圈為帶地形的二次場歸一化圖，黑色星號為不帶地形的二次場歸一化圖，從中也能發(fā)現(xiàn)地形對電磁場的影響.

7 結(jié)論

本文提出了一種基于完全非結(jié)構(gòu)化網(wǎng)格的海洋電磁有限單元三維正演算法.相比于傳統(tǒng)的結(jié)構(gòu)化網(wǎng)格，非結(jié)構(gòu)化網(wǎng)格可以大大地減少計(jì)算量，同時(shí)可以模擬復(fù)雜的地電模型，譬如，海底地形.文中求解基于Coulomb規(guī)范的矢量位和標(biāo)量位Maxwell方程組，利用不完全LU分解的QMR迭代算法有效提高收斂速度，采用二次場算法回避總場算法的場源奇異性，選用更符合實(shí)際情況的層狀地電模型作為背景介質(zhì)，討論了地形起伏及電導(dǎo)率各向異性的海洋電磁三維響應(yīng)規(guī)律和特征.通過對若干典型地電模型對比正演數(shù)值解與準(zhǔn)解析解和積分方程解，驗(yàn)證了文中所提算法的可靠性，該算法適用于復(fù)雜三維地形和各向異性電導(dǎo)率介質(zhì)的電磁模擬.

圖18 復(fù)雜海底地形模型三維示意圖Fig.18 A 3D view of the model with complex seafloor bathymetry

圖19 模型4正演結(jié)果的總場和背景場在y=0處的比較圖.藍(lán)色實(shí)線為背景場，紅色點(diǎn)線為含地形效應(yīng)的總場， 黑色虛線為不含地形效應(yīng)的總場Fig.19 A comparison of the total field and background at y=0 for Model 4. The solid blue line represents the background field, the dotted red lines indicates the total field with seafloor bathymetry while the dashed black line is the total field without the bathymetry

致謝 作者感謝美國猶他大學(xué)電磁法正反演科研組(CEMI)提供本人優(yōu)良的科研環(huán)境.

Andréis D, MacGregor L. 2008. Controlled-source electromagnetic sounding in shallow water: Principles and applications.Geophysics, 73(1): F21-F32, doi: 10.1190/1.2815721. 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, doi: 10.1190/1.1444968.

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, doi: 10.1016/j.cageo.2014.09.008.

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

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, doi: 10.1190/geo2010-0398.1.

Srnka L J, Carazzone J J, Ephron M S, et al. 2006. Remote reservoir resistivity mapping.TheLeadingEdge, 25(8): 972-975, doi: 10.1190/1.2335169.

Freund R W, Nachtigal N M. 1991. QMR: a quasi-minimal residual method for non-Hermitian linear systems.Numer.Math., 60(1): 315-339, doi: 10.1007/BF01385726.

Guptasarma D, Singh B. 1997. New digital linear filters for Hankel J0and J1transforms.Geophys.Prospect., 45(5): 745-762, doi: 10.1046/j.1365-2478.1997.500292.x.

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

Newman G A, Alumbaugh D L. 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences.GeophysicalProspecting, 43(8): 1021-1042, doi: 10.1111/j.1365-2478.1995.tb00294.x.

Puzyrev V, Koldan J, de la Puente J, et al. 2013. A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling.Geophys.J.Int., 193(2): 678-693, doi: 10.1093/gji/ggt027.Sasaki Y. 2011. Bathymetric effects and corrections in marine CSEM data.Geophysics, 76(3): F139-146, doi: 10.1190/1.3552705.

Um E S, Alumbaugh D L. 2007. On the physics of the marine controlled-source electromagnetic method.Geophysics, 72(2): WA13-WA26, doi: 10.1190/1.2432482.

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.Geophys.J.Int., 193(3): 1460-1473, doi: 10.1093/gji/ggt071. Ward S H, Hohmann G W. 1988. Electromagnetic theory for geophysical applications // Nabighian M N ed. Electromagnetic Methods in Applied Geophysics—Theory. SEG, 130-311.

Wu J P, Wang Z H, Li X M. 2004. Efficient Solution and Parallel Computation for Sparse Linear Systems (in Chinese). Changsha: Hunan Science and Technology Press. Xu S Z. 1994. FEM in Geophysics (in Chinese). Beijing: Science Press.

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, doi: 10.1190/1.2358403.

Zhdanov M S. 2010. Electromagnetic geophysics: Notes from the past and the road ahead.Geophysics, 75(5): 75A49-75A66, doi: 10.1190/1.3483901.

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

吳建平, 王正華, 李曉梅. 2004. 稀疏線性方程組的高效求解與并行計(jì)算. 長沙: 湖南科學(xué)技術(shù)出版社.

徐世浙. 1994. 地球物理中的有限單元法. 北京: 科學(xué)出版社.

(本文編輯 胡素芳)

Three-dimensional marine controlled-source electromagnetic modelling in anisotropic medium using finite element method

CAI Hong-Zhu1, XIONG Bin2*, Michael Zhdanov1

1CollegeofMines&EarthSciences,UniversityofUtah,Saltlakecity，UT84112,USA2CollegeofEarthSciences,GuilinUniversityofTechnology,Guilin541004,China

In recent years, we have observed a strong interest in applying controlled-source electromagnetic (CSEM) method in hydrocarbon exploration to reduce the uncertainty resulting from seismic method, especially in the marine environment. Comparing to seismic method, CSEM method can be potentially used to distinguish the oil-bearing and water bearing reservoir. In practical application, the subsurface conductivity is characterized by three-dimensional inhomogeneity and anisotropy. In the meantime, the CSEM data can also be distorted significantly by the seafloor bathymetry in marine environment. To correctly interpret the subsurface structure with the observed CSEM data, it is desirable to develop full three-dimensional forward modeling algorithm which can address the conductivity anisotropy and the variation of seafloor bathymetry.Generally, there are three main numerical methods for the modeling of CSEM data: integral equation, finite difference, and finite element methods. Among these three methods, the finite element method with unstructured tetrahedral mesh is more suitable for simulating complex geological structures such as seafloor bathymetry. Comparing to the integral equation method, the finite element stiffness matrix is much sparser and as such the method can be used for large scale modeling of CSEM data. We adopt the node-based finite element method in this paper. For CSEM modeling, the accuracy is affected significantly by the representation of source. In order to avoid the source singularity problem caused by the total field formulation, we use the secondary field formulation and consider that the background conductivity is formed by a series of infinite horizontal layers and each with a constant conductivity value. We consider that the background conductivity is characterized by transverse anisotropy which is typically the case for marine CSEM. For the secondary field formulation, one needs to compute the background field on the mesh. We use the fast Hankel transform method to compute the background field. We adopt the classic node-based finite element method for the forward modeling. In order to avoid the spurious solution which is usually encountered in solving the Maxwell's equation directly using node-based finite element, we solve the secondary scalar and vector electromagnetic potential based on the Coulomb gauge. We re-formulated the original Maxwell's equation for the secondary scalar and vector potential. The computed background field from fast Hankel transform is converted to the primary scalar and vector potential. After applying the Galerkin finite element analysis, one can obtain a large sparse system of finite element equations. We applied the effective Dirichlet boundary condition and consider that both the scalar and vector secondary electromagnetic potential vanish on the boundary of the study domain. We solve the sparse finite element system of equations using quasi-minimum residual (QMR) method with an incomplete LU decomposition (ILU) preconditioner. Once the secondary scalar and vector potential is solved, we can compute the secondary electric and magnetic field by taking the differential of the secondary potential. The general numerical differential method can introduce significant noise. In order to solve this problem, we use a weighted moving least square method to calculate the secondary electric and magnetic field from the computed secondary potential.To validate the effectiveness and correctness of the developed algorithm, we have implemented several model studies to compare the electromagnetic field computed from our method with the 1D analytical solution and integral equation solution. We first consider a 1D model with air and half-space as the background. We consider that there exists an infinite thin layer embedded in the half-space background. The anomalous field caused by this infinite thin layer can be computed either with finite element or fast Hankel transform. Our numerical result shows that the solution from finite element method is practically the same as the 1D analytical solution obtained from fast Hankel transform either for secondary electric or magnetic field. We have also done a comparison study for the QMR solver with Jacobian and ILU preconditioner.The result demonstrates that the convergence can be speeded up significantly by adopting the ILU preconditioner. We also consider a 3D marine reservoir model with flat bathymetry. For this model the secondary field can be effectively computed by integral equation method. As such, we have compared our results with integral equation solution for this model. In the model, we also considered the conductivity anisotropy for both reservoir and the marine sediment.Numerical result shows that the finite element solution for the problem with conductivity anisotropy is practically the same as integral equation solution. Following this, we consider a marine reservoir model with simple seafloor bathymetry represented by a trapezoidal hill. For integral equation method, we simulate this trapezoidal hill model by a staircase approximation. The numerical results show that the bathymetry effect simulated from finite element method is very close to the integral equation method especially when the receivers are far away from the bathymetry. The difference between the finite element and integral equation method can be attributed to the staircase approximation of the bathymetry for integral equation modeling. Finally, we consider a complex seafloor bathymetry model with reservoir and conductivity anisotropy. Our numerical modeling result using the finite element method have demonstrated that the observed CSEM data can be significantly distorted by the seafloor bathymetry in practical application and this effect should be addressed by robust forward modeling algorithm for the correct interpretation of subsurface geo-electric structures.We have implemented a three-dimensional numerical modeling method for marine CSEM data using node-based finite element method based on a fully unstructured tetrahedral mesh. Comparing to the conventional structured mesh, the computation load can be significantly reduced by adopting the unstructured mesh. In the meantime, the unstructured mesh is capable of simulating complex geo-electric structures such as seafloor bathymetry. We have formulated our problem for the scalar and vector potential instead of the electric and magnetic field to address the spurious solution. In order to avoid the source singularity problem for CSEM modeling, we adopt a secondary field formulation. The sparse finite element system of equations is solved using a quasi-minimum residual method with incomplete LU decomposition result as the preconditioner to speed up the convergence. The developed algorithm was tested for several synthetic and realistic models. The numerical results have demonstrated the effectiveness of the algorithm for the modeling of complex geo-electric structures with conductivity anisotropy.

Marine CSEM; Unstructured mesh; Coulomb gauge; 3D modeling; Finite element; Moving least-squares interpolation

國家自然科學(xué)基金項(xiàng)目(40974077、41164004)，廣西自然科學(xué)基金項(xiàng)目(2011GXNSFA018003、2013GXNSFAA019277)聯(lián)合資助.

蔡紅柱，男，1989年出生，2009年畢業(yè)于中南大學(xué)，猶他大學(xué)在讀博士. 主要從事電磁場三維正反演、位場以及其全張量的三維正反演與快速成像方面的研究. E-mail:caihongzhu@hotmail.com

*通訊作者熊彬，男，1974年生，博士，教授，主要從事電磁場理論及反演成像方面的教學(xué)與科研工作.E-mail:xiongbin@msn.com

10.6038/cjg20150818.

10.6038/cjg20150818

P631

2014-07-02，2014-12-09收修定稿

蔡紅柱，熊彬，Michael Zhdanov. 2015. 電導(dǎo)率各向異性的海洋電磁三維有限單元法正演.地球物理學(xué)報(bào)，58(8):2839-2850,

Cai H Z, Xiong B, Michael Zhdanov. 2015. Three-dimensional marine controlled-source electromagnetic modelling in anisotropic medium using finite element method.ChineseJ.Geophys. (in Chinese),58(8):2839-2850,doi:10.6038/cjg20150818.