王健, 陳浩, 王秀明

1 中國科學院聲學研究所聲場與聲信息國家重點實驗室， 北京　100190 2 中國科學院大學， 北京　100049

?

用于固體礦床多分量感應測井響應模擬的矢量有限元法

王健1,2, 陳浩1, 王秀明1

1 中國科學院聲學研究所聲場與聲信息國家重點實驗室， 北京100190 2 中國科學院大學， 北京100049

摘要本文開發(fā)了基于非結(jié)構化四面體網(wǎng)格的三維矢量有限元法，實現(xiàn)了固體礦床井眼中多分量感應測井響應的數(shù)值模擬，并分析了多分量感應測井儀器在復雜礦床模型中的響應特征.本文通過采用幾何因子背景場，有效地避免了源的奇異性問題；同時，在井眼邊界采用非均質(zhì)網(wǎng)格并用Gauss-Legendre積分計算四面體單元的等效電導率.利用LU分解求解線性方程組，實現(xiàn)了一次網(wǎng)格劃分多點的數(shù)值計算，提高了計算效率，從而實現(xiàn)快速連續(xù)的多分量感應測井模擬.非結(jié)構化的四面體網(wǎng)格確保了該方法可以模擬實際問題中所能遇到的復雜的礦體模型.基于水平三層分層和徑向分層模型，驗證了算法在各向同性和各向異性介質(zhì)中的可靠性.我們還以三個不同的礦床模型為例，研究了多分量感應測井儀的不同分量的探測特性，結(jié)果表明，結(jié)合九個分量的信息，可以探測礦體的深度，也可以識別礦體的方位和走向，為精確地描述礦體的三維分布特征打下了基礎.

關鍵詞矢量有限元； 多分量感應測井； 固體礦床

Due to the characteristic of solid deposit, such as wide variety, complex distribution, we use the vector finite element method to simulate the response of multicomponent induction instrument, which has an advantage for medium with complex structure and boundary based on unstructured mesh. Moreover, the method could avoid defects of nodal finite element method. Specifically, geometry factor background field is used to eliminate the singularities of source efficiently. Meanwhile, heterogeneous grids are introduced at the boundary of the borehole and the equivalent conductivities of the tetrahedron elements are calculated by using Gauss-Legendre integral. To take advantage of the LU decomposition to solve the linear system of equations, we could obtain multipoint numerical results for each meshing. It greatly improves the efficiency. So as to ensure a rapid and successive simulation of the multicomponent induction logging.

It has been tested and verified that this method is valid and reliable in the isotropic and anisotropic medium based on the horizontal three-layered and radial model. In the dipping layered formation, all nine components are studied. TheXXandYYresponses are more sensitive to the dipping angle than theZZwith sharp dipping angle; the cross-componentsXZandZXare influenced by border and anisotropy. The apparent conductivity curves ofXZandZX-components show horns in the opposite direction. While in the spherical and cylindrical ore body model, we show azimuthal sensitivity of all nine components of apparent conductivity. Circle lines for response ofZZ-component confirm that the lack of azimuthal sensitivity of the conventionalZZ.The remaining eight components do show azimuthal sensitivity, but four of them (XX,XY,YY,YX) have the period of 180°, while the other four (XZ,YZ,ZX,ZY) have the period of 360°.It is quite obvious that having any combination of measurement from the first group will not allow a resolution between left and right because of the 180°symmetry. The pairs of signals from the second group (ZX,ZY) or (XZ,YZ) uniquely define the azimuthal direction of the ore body.

We develop a three dimensional vector finite element program for multicomponent induction logging based on unstructured tetrahedral mesh in frequency domain, then we use it to simulate the response of multicomponent instrument in solid deposit. The numerical results show that it is able to identify anisotropy of ore body as same as in oil and gas reservoir. TheZZ-component indicates the vertical conductivity, and theXXandYY-components indicate the horizontal conductivity. The cross-componentsXZ,ZX,YZ,ZYhave the period of 360°, therefore provide azimuthal information by which you could laid a foundation to accurately describe the distribution of ore deposit.

1引言

自從Doll(Doll,1949)提出了幾何因子理論以來，感應測井在油氣儲層勘探和開發(fā)中得到了廣泛的應用(Moran and Kunz,1962).這種方法通過測量井下的載流線圈在井附近介質(zhì)中產(chǎn)生的感應電流所形成的二次場與周圍介質(zhì)電導率的密切關系來獲得介質(zhì)的電導率信息.在礦床勘探中，應用感應測井方法識別電性異常和磁性異常礦床的想法早已有之(Broding et al.,1952;Cheryauka and Sato,2002).然而，包括感應測井方法在內(nèi)的測井方法在礦床勘探中的應用迄今遠不如在石油勘探中普及.礦床測井發(fā)展比較滯后的主要原因是長期以來礦床的勘探和開采主要集中在淺部，鉆井和取樣成本比較低，因而對礦床的精細評價主要是對取樣的礦體進行直接的測量，而非石油勘探中的間接的地球物理測量(Baltosse and Lawrence,1970).但是隨著礦床的埋深增加，地面地球物理方法的分辨率和精度降低，鉆井和取樣成本增加，人們開始逐步重視測井方法在深部礦床勘探中的作用.目前一些井中地球物理方法已經(jīng)在深部找礦中取得了良好的效果.例如：加拿大的薩德伯里銅、鎳礦區(qū)通過深部鉆孔中瞬變電磁方法組合，從20世紀80年代中期以來相繼發(fā)現(xiàn)了一批深部銅、鎳硫化物礦床(滕吉文等，2007).雖然礦床測井與石油測井有一定的相似性，但前者也有獨特之處：除了礦床種類繁多以外，固體的礦床分布復雜，形態(tài)多樣，空間分布有限，由此礦體的巖石物理模型一般都為復雜的三維模型，而油氣多為層狀的二維模型，因此對礦體的空間分布的評價除了要給出礦體的深度信息，還要求給出礦體的方位和走向信息.

多分量感應測井儀器由正交的三個發(fā)射線圈和正交的三個接收線圈系陣列組成，能夠測量XX，YY，XY，YX，XZ，YZ，ZX，ZY和ZZ分量(Kriegshauser et al.,2000;Zhdanov et al.,2001).其中分量的第一個字母代表發(fā)射線圈的方向，第二個字母代表接收線圈的方向.對這些分量的分析處理能夠識別儲層或礦體的各向異性電導率、裂縫、斷裂、礦體邊界、方位和走向等信息.因此該方法在未來的礦床評價中會有獨特的優(yōu)勢，但要從九個分量中正確獲取這些信息遠比只從ZZ分量進行電導率評價時困難的多(Rabinovich et al.,2006;Davydycheva,2010;肖加奇等,2013)，因此需要開展大量的理論研究和數(shù)值模擬分析，充分了解和認識不同分量的探測特性，并進一步提出反演或評價的方法.目前有關多分量感應測井的研究大多是針對具有平面分布特征的油氣儲層(Wang et al.,2008;汪宏年等,2008；魏寶君等,2009;Davydycheva et al.,2009)，而對其在具有復雜結(jié)構分布的固體礦床勘探中的研究則較少.

在電磁場正演模擬中，最常用的數(shù)值方法是有限差分法(Finite difference method:FDM)和有限元法(Finite element method:FEM)，其他數(shù)值方法諸如積分方程法(Avdeevet al.,2002;Fang et al.,2006)和有限體積方法(Weiss,2013)也有一定的應用.有限差分法廣泛用于解決2-D和3-D時域和頻域問題(Newman and Alumbaugh,1995;Weiss and Newman,2002).基于結(jié)構化網(wǎng)格的有限差分方法比有限元方法更容易實現(xiàn)，但它對網(wǎng)格進行局部優(yōu)化十分困難，任何網(wǎng)格尺寸的改變對于整個問題的計算量有很大的影響(Puzyrev et al.,2013).而且由于復雜的模型邊界不能夠用矩形網(wǎng)格精確的近似從而對結(jié)果的精度造成極大的影響.因此如果目標模型是復雜的幾何體，有限元法由于支持非結(jié)構化的網(wǎng)格而顯得更有優(yōu)勢.其所產(chǎn)生的網(wǎng)格能夠準確地模擬復雜邊界和高電導率對比度的情況，這種情況僅需在介質(zhì)的分界面、源和接收器處進行精細的網(wǎng)格劃分.

在電磁勘探中，有限元方法雖然已經(jīng)得到了廣泛應用(Chang and Anderson,1984;Key and Ovall,2011;da Silva et al.,2012)，但也還存在二個不足.(1)傳統(tǒng)的節(jié)點有限元方法不能夠正確處理介質(zhì)分界面處電場強度法向方向的不連續(xù)性.(2)有散度不等于0的非物理解存在(Jin,2002).對于上述問題，一種解決的辦法是用規(guī)范的磁矢勢方程(Biro and Preis,1989;Badea et al.,2001)而不是直接的電場強度E或者磁場強度H方程.另一種方法是選用矢量有限元方法(Nédélec,1980;孫向陽等,2008;Li et al.,2011;Mukherjee and Everett,2011)，它的自由度是沿單元邊的切向分量.矢量有限元在單元內(nèi)強加了散度為0的條件，并允許界面處的電場強度法向分量的不連續(xù)性.

本文采用基于非結(jié)構化四面體網(wǎng)格的三維矢量有限元方法對多分量感應測井儀器的響應進行數(shù)值模擬.首先，給出了矢量有限元在各向異性電導率介質(zhì)中的原理與方法.其次，通過與解析結(jié)果的對比驗證了方法和代碼的正確性并分析了井眼網(wǎng)格剖分方法對數(shù)值精度的影響.然后，模擬了多分量感應測井儀器在傾斜各向異性地層中的響應，分析了傾角和各向異性對九個分量的影響，著重分析了交叉的XZ和ZX分量的曲線特征.最后對距井壁一定距離的球形礦體和傾斜的圓柱狀礦體計算了儀器相對于礦體旋轉(zhuǎn)時不同分量所具有的探測特性.本文嘗試利用數(shù)值方法對多分量感應測井儀器在復雜模型下的響應進行數(shù)值模擬，以便對測量得到的視電導率進行有效的校正和反演，從而推動該儀器在固體礦床測井評價中的應用.

2基本原理

在各向異性的導電介質(zhì)中，電場E和磁場H所滿足的微分Maxwell方程為

(1)

(2)

求解方程(2)相對于求解方程(1)的優(yōu)點是其有較少的變量，從而有限元離散得到的自由度較小(da Silva et al.,2012).由于感應測井的頻率比較低，位移電流可以被忽略不計.為了提高計算效率和避免直接用網(wǎng)格離散源所帶來的困難，本文將總場分解成背景場和散射場(Newman and Alumbaugh,1995；王健等,2015).這樣(2)式中的源可以由一系列已知的背景場Eb替代.定義散射場Es和異常電導率張量如下：

E=Eb+Es,

(3)

將式(3)代入式(2)中，可得到關于散射場滿足的方程為

(4)

在柱坐標系下，有限尺寸的載流線圈源在電導率為σb的各向同性均勻介質(zhì)中形成的背景場為

×J1(λac)J1(λρ)λdλ，

(5)

(6)

一般來說，求解域的范圍取4～5個趨膚深度就可滿足要求，趨膚深度由模型中的最小電導率決定.由于四面體網(wǎng)格能更好地處理復雜邊界，本文采用開源軟件NETGEN(Sch?berl，1997)生成非結(jié)構化的四面體網(wǎng)格進行空間離散.

通過有限元分析(附錄A)將計算得到的單元矩陣Ke，F(xiàn)e和Ge分別組成矩陣A和右端項b，其中矩陣A是復對稱稀疏矩陣，

表1　數(shù)學符號

Ax=L Ux=b，

(7)

(7)式中L和U是下和上三角矩陣.一旦L和U已經(jīng)計算得到，那么對于任意的右端項b，可以用下式求得x，

(8)

其中y=L-1b.

求解時LU分解所花費的時間要遠高于計算式(8)所需要的時間.在感應測井問題中A矩陣只與網(wǎng)格有關，b矩陣與源相關.因此當儀器沿井軸移動時，如果只進行一次網(wǎng)格剖分，則可以只進行一次LU分解并利用(8)式實現(xiàn)多個測量點響應的快速計算.但是值得注意的是：為了滿足邊界條件并保證數(shù)值精度，需要網(wǎng)格范圍取的較大，并且需要對儀器軌跡上的網(wǎng)格進行加密.考慮到LU分解對計算機的內(nèi)存需求較大，我們設定一個自由度的閾值，當自由度超過閾值時，測量點就被分成Ns組，每組進行一次網(wǎng)格剖分.這樣既盡可能減少網(wǎng)格剖分次數(shù)并利用LU求解的特點，同時又節(jié)省了內(nèi)存.本文使用多線程并行求解器UMFPACK(Davis,2004)求解線性方程組.

3精度驗證

為了驗證算法的精度，本文采用具有解析結(jié)果的三層模型.將開發(fā)的有限元程序?qū)﹄p線圈儀器計算的結(jié)果與解析結(jié)果進行對比，三層地層的電導率分別為1 S·m-1，0.002 S·m-1，20 S·m-1，中間層厚2 m.發(fā)射源電流強度等于1 A，頻率為20 kHz.發(fā)射線圈半徑和源距分別為0.03 m和1 m.模型中感應電流在電導率最低的介質(zhì)中的趨膚深度δ≈80 m.因而將模型區(qū)域設定為一個半徑和高分別為398 m和896 m的圓柱.最終網(wǎng)格的總節(jié)點數(shù)為64212，四面體單元總數(shù)為360128，相應的自由度為429156.模擬的結(jié)果如圖1所示.有限元模擬的結(jié)果與解析結(jié)果吻合較好：最大誤差小于2%.模型中電導率對比度高達104，這說明本文的算法有能力處理實際問題中出現(xiàn)的高電導率對比度的情況.文中所使用的計算機配置為4核3.1 GHz的i5-2400 CPU，內(nèi)存12G，每個深度點平均計算耗時2.8 s.

圖1　水平三層地層模型的有限元法與解析結(jié)果對比

(9)

圖2是采用均質(zhì)網(wǎng)格和非均質(zhì)網(wǎng)格進行模型離散的網(wǎng)格截面圖.圖2a由于采用均質(zhì)網(wǎng)格，井眼內(nèi)的網(wǎng)格尺寸受到井壁的限制而不能擴大，在距離源較遠的區(qū)域產(chǎn)生了大量不必要的小尺寸單元;圖2b中由于采用非均質(zhì)網(wǎng)格，網(wǎng)格尺寸能夠迅速增大，減少了總的單元數(shù)量.同時，這種方法也能夠保證結(jié)果的精度.當井內(nèi)外介質(zhì)電導率對比度增大時，一般可以通過增加積分點來提高精度.圖3中我們將采用非均勻介質(zhì)網(wǎng)格的有限元計算的結(jié)果與解析結(jié)果在1D徑向分層地層上進行對比.其中泥漿電導率為0.005S·m-1，井外地層電導率是0.01～50S·m-1的對數(shù)等間隔插值，儀器結(jié)構與上例相同.程序采用45點Gauss-Legendre積分計算非均勻質(zhì)單元的平均電導率.從圖中可以看到：模擬的結(jié)果與解析結(jié)果吻合很好，證實了這種處理井眼的方法在井內(nèi)外電導率對比度高達104時仍然可取得令人滿意的精度.

為了驗證本文的有限元算法和程序在各向異性介質(zhì)中的正確性，我們比較了用有限元方法和格林函數(shù)法計算的多分量儀器在三層各向異性地層(SunandNie,2008)模型中的響應.地層的參數(shù)詳見圖4a.圖4b的結(jié)果表明，由有限元計算的水平電導率和垂直電導率與格林函數(shù)法計算的十分吻合.

圖3　徑向分層地層模型的有限元法與解析結(jié)果對比

圖2　感應測井網(wǎng)格切面圖. (a) 均質(zhì)網(wǎng)格； (b) 非均質(zhì)網(wǎng)格

圖4　各向異性模型驗證. (a)三層各向異性模型; (b)結(jié)果對比

4實例分析

4.1傾斜各向異性地層模型

(10)

(11)

其中T為轉(zhuǎn)動矩陣，θ為繞Y軸的地層傾角，φ為繞Z軸的地層方位角.

圖5　多分量感應測井儀器的基本結(jié)構等效圖

圖6　傾斜各向異性地層模型

圖7　均勻各向異性傾斜地層的XZ分量和ZX分量

圖7給出了在均勻各向異性地層中θ=60°和-60°的XZ和ZX分量，當不存在界面時，XZ分量和ZX分量相等，且θ=60°時有負的幅度，而θ=-60°有正的幅度.圖8a—8c分別給出了傾角θ=0°，60°和-60°時9個分量的視電導率曲線.θ=0°時只有共面的XX，YY和共軸的ZZ分量不為0，其余6個分量均為0，XX分量和YY分量相等.在圍巖的各向同性地層中，XX和YY分量小于ZZ分量說明共面的XX，YY分量受趨膚效應影響更大.當θ=60°和-60°時，各向異性目標層的XZ分量和ZX分量不等于0，XX分量和YY分量開始分離，并和ZZ分量接近.圖8d給出了θ=60°時且目標層電導率為0.25 S·m-1的各向同性地層的9個分量的視電導率曲線.與相同角度的各向異性地層比較，我們可以看到：XZ分量和ZX分量僅在界面附近存在非0的極化角且其方向和各向異性地層相同，這說明極化角的方向與是否存在各向異性無關，而只與地層的傾斜角度以及邊界兩側(cè)介質(zhì)電導率的大小有關.當?shù)貙诱齼A斜(0°<θ<90°)時，儀器從低阻地層進入高阻地層時XZ分量的極化角方向為正.當儀器從高阻地層進入低阻地層時，極化角方向為負且其幅度要小于正向幅度.當?shù)貙迂搩A斜(-90°<θ<0°)時，儀器從低阻地層進入高阻地層XZ分量的極化角方向為負.當儀器從高阻地層進入低阻地層時，極化角方向為正且其幅度要小于負向幅度.ZX曲線的規(guī)律恰好和XZ曲線的相反.

4.2球形礦體模型

圖10是本文計算的球形礦體模型，其中各向同性的球形礦體半徑為2 m，球心距離井軸2.6 m，礦體的電導率為5 S·m-1，周圍地層電導率為0.05 S·m-1.儀器X方向發(fā)射線圈的法向方向與礦體的相對角度φ從0°變化到360°.儀器固定在Z軸上且沒有穿過礦體，其線圈系中點平行于礦體球心.圖11給出了儀器在此模型下響應的9個分量，在極坐標下清晰地展示了各個分量所具有的周期性.當?shù)V體的方位從0°到360°變化時，常規(guī)儀器得到的ZZ分量視電導率曲線呈圓形，因此其不具備方位識別的能力，而且當儀器沒有穿過礦體時，感應電流串聯(lián)地通過地層和礦體，因此其值接近于在均勻地層中的響應，不能反映出礦體的存在.共面的XX分量和YY分量具有一定的方位敏感性，具有180°的周期.XY分量和YX分量也具有180°周期，其與XX分量和YY分量的區(qū)別是它們的幅度存在正負變化，在45°和225°取得極大值，在135°和315°取得極小值，在0°，90°，180°，270°不存在XY和YX分量.顯然，由于具有180°的對稱性，無論如何組合上述4個分量，都無法確定礦體的方位.如圖所示，XZ，YZ，ZX，ZY這四個交叉分量具有360°周期，且他們的極大值或極小值總是對應著礦體的方位或者相反的方位.XZ和ZX與YZ和ZY有90°相位差，因此利用(XZ，YZ)或者是(ZX，ZY)可以唯一確定礦體與井的相對位置.如果儀器是旋轉(zhuǎn)測量的，那么僅僅依靠ZX或者XZ分量的信息就能充分地判斷礦體的方位.

由于多分量感應測井響應的復雜性，對其測量得到的視電導率曲線的直觀解釋十分困難.一種利用分量XZ-ZX和XZ+ZX以分離各向異性和邊界影響的方法被提出(Omeragic et al.,2006;Davydycheva,2010)，以減少多分量感應測井曲線解釋的復雜性.本文將該方法用于計算球形礦體的方位角中，計算方法如下式所示：

圖9　不同地層傾斜角度的視電導率曲線

圖11　多分量感應測井在球形礦體模型中的響應

(12)

如圖12所示，縱坐標是通過式(12)計算得到的方位角，橫坐標是實際的方位角.計算的結(jié)果和實際的方位角吻合得很好，說明利用交叉分量計算礦體的方位角是可行的.盡管反正切變換的角度不能超過90°，但是根據(jù)(YZ-ZY)和(XZ-ZX)的符號可以很容易確定其所在的象限.

圖10　球形礦體模型

4.3傾斜圓柱形礦體模型

圖13是本文計算的傾斜礦體模型，其中各向同性的圓柱型礦體的長度為6 m，半徑1.5 m，傾斜角度θ=45°，礦體電導率為1 S·m-1，周圍地層電導率為0.05 S·m-1.儀器X方向發(fā)射線圈的法向方向與礦體的相對角度φ從0°變化到360°.儀器沿平分礦體的Z軸正方向移動.圖14給出了儀器在此模型下測量得到的9個分量，其中ZZ分量提供了礦體的深度信息，在中間位置出現(xiàn)了高視電導率層，厚度約為3 m，反映出圓柱型礦體半徑的大小.但ZZ分量從0°到360°沒有任何變化，因此無法判斷礦體走向等信息.XX和YY分量也反映了礦體的存在，但與ZZ分量相比所反映的視厚度更小，在礦體長度的邊緣位置呈現(xiàn)出極化角.XY和YX分量的幅值只在礦體的邊緣才出現(xiàn)明顯的周期正負變化.這是因為當儀器處于礦體中時，礦體的走向?qū)Y和YX分量影響不明顯.XZ，YZ，ZX，ZY具有360°周期，它們的極大值或極小值總是對應著礦體的走向或者正交方向，因此利用(XZ，YZ)或者是(ZX，ZY)可以唯一確定礦體的走向.我們注意到：XZ和ZX分量在礦體的上和下是不對稱的，同時彼此是反對稱的，因此可以依靠XZ或者ZX分量確定某一深度點儀器是在礦體之上或者之下，以及礦體的傾斜方向.上述二個例子說明利用多分量感應測井儀器的交叉分量可以獲得礦體的方位和走向等信息.當然利用ZX，XZ，ZY，YZ分量判斷礦體的方位和走向的前提條件是儀器具有合適的探測深度和分辨率.

5結(jié)論

本文在頻率域內(nèi)實現(xiàn)了基于非結(jié)構化四面體網(wǎng)格的三維矢量有限元方法的多分量感應測井響應模擬算法，并模擬了多分量感應測井儀器在礦床測井模型中的響應.為了提高計算效率，算法在井眼邊界采用非均質(zhì)網(wǎng)格并利用Gauss-Legendre積分計算四面體單元的平均電導率.同時結(jié)合感應測井模型和LU分解方法的特點，實現(xiàn)了一次網(wǎng)格劃分多深度點的數(shù)值計算.通過三個算例，驗證了算法在各向同性介質(zhì)和各向異性介質(zhì)中有很好的計算精度和速度.

數(shù)值結(jié)果表明，多分量感應測井儀器能夠識別礦體的各向異性.在水平地層情況下，其ZZ分量能夠給出水平電導率，XX和YY分量給出垂直電導率.對于儀器探測深度范圍內(nèi)的有限尺寸礦體模型，由于XZ，YZ，ZX，ZY分量具有360°的方位敏感性，因而能夠提供礦體的方位和走向等信息.

圖12　球形礦體方位角計算

圖13　傾斜圓柱型礦體模型

附錄A矢量基函數(shù)和有限元分析

四面體單元內(nèi)的節(jié)點和棱邊的關系如圖A1和表A1所示.對于一個四個頂點分別為(a，b，c，d)的四面體單元，其形函數(shù)φa(ra)有如下定義(SilvesterandFerrari,1996)：

圖A1　四面體單元

圖14　多分量感應測井儀器在傾斜圓柱形礦體(θ=45°)模型中的響應

邊端點1端點21ab2ac3ad4bc5bd6cd

(A1)

(A1)式中r表示四面體中任一點的位置，其他頂點的形函數(shù)有相似的定義.基于矢量有限元，散射場Es可表示為

(A2)

其中ai是待確定的系數(shù)，αi(r)是矢量基函數(shù)并有如下形式：

(A3)

參數(shù)(a，b)是網(wǎng)格中第i條邊的始末端點.應用伽遼金方法(ZienkiewiczandTaylor,2005)，并令αj為試函數(shù)，則得到方程(4)的弱解形式：

(A4)

方程(A4)可進一步展開：

(A5)

(A6)

(A7)

其有離散形式如下：

(A8)

在得到了網(wǎng)格邊上的散射場值后，需要將其轉(zhuǎn)換為節(jié)點上的散射場值.如果端點為(a，b)的棱邊屬于總計P個四面體，并且端點a被D條邊共用，則節(jié)點a上的散射場可表示為

(A9)

有關后處理方法更詳細的討論參見文獻(MukherjeeandEverett,2011).

References

Avdeev D B, Kuvshinov A V, Pankratov O V, et al. 2002. Three-dimensional induction logging problems, Part I: An integral equation solution and model comparisons.Geophysics, 67(2): 413-426.

Badea E A, Everett M E, Newman G A, et al.2001. Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials.Geophysics, 66(3): 786-799.

Baltosser R W, Lawrence H W. 1970.Application of well logging techniques in metallic mineral mining.Geophysics, 35(1): 143-152.

Biro O, Preis K. 1989. On the use of the magnetic vector potential in the finite-element analysis of three-dimensional eddy currents.IEEETransactionsonMagnetics,25(4): 3145-3159.Broding R A, Zimmerman C W, Somers E V, et al. 1952. Magnetic well logging.Geophysics, 17(1): 1-26.

Chang S K, Anderson B. 1984. Simulation of induction logging by the finite-element method.Geophysics, 49(11): 1943-1958.

Cheryauka A B, Sato M. 2002. Directional induction logging for evaluating layered magnetic formations.Geophysics, 67(2): 427-437.

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.Davis T A. 2004. Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method.ACMTransactionsonMathematicalSoftware(TOMS), 30(2): 196-199.Davydycheva S, Homan D, Minerbo G. 2009.Triaxial induction tool with electrode sleeve: FD modeling in 3D geometries.JournalofAppliedGeophysics, 67(1): 98-108.

Davydycheva S. 2010. Separation of azimuthal effects for new-generation resistivity logging tools—Part I.Geophysics, 75(1):E31-E40.

Doll H G. 1949. Introduction to induction logging and application to logging of wells drilled with oil base mud.JournalofPetroleumTechnology, 1(6): 148-162.

Druskin V L, Knizhnerman L A, Lee P. 1999. New spectral Lanczos decomposition method for induction modeling in arbitrary 3-D geometry.Geophysics, 64(3): 701-706. Fang S, Gao G, Torres-Verdín C. 2006. Efficient 3D electromagnetic modelling in the presence of anisotropic conductive media, using integral equations.ExplorationGeophysics, 37(3): 239-244.

Guptasarma D, Singh B. 1997. New digital linear filters for Hankel J0 and J1 transforms.GeophysicalProspecting, 45(5): 745-762.

Jin J M. 2002. The Finite Element Method in Electromagnetics. New York: Wiley.

Kennedy D, Peksen E, Zhdanov M. 2001. Foundations of tensor induction well-logging.Petrophysics, 42(6):588-610.

Key K, Ovall J. 2011. A parallel goal-oriented adaptive finite element method for 2.5-D electromagnetic modelling.GeophysicalJournalInternational, 186(1): 137-154.Kriegshauser B, Fanini O, Forgang S, et al. 2000. A new multicomponent induction logging tool to resolve anisotropic formations. ∥ SPWLA 41st Annual Logging Symposium. Dallas,Texas.

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.Moran J H, Kunz K S. 1962. Basic theory of induction logging and application to study of two-coil sondes.Geophysics, 27(6): 829-858.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 in3.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.

Omeragic D, Dumont A, Esmersoy C, et al. 2006. Sensitivities of directional electromagnetic measurements for well placement and formation evaluation while drilling. ∥ 76th SEG Annual Meeting. New Orleans, LA.

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.Rabinovich M, Tabarovsky L, Corley B, et al. 2006. Processing multi-component induction data for formation dips and anisotropy.Petrophysics, 47(6): 506-526.

Rathod H T, Venkatesudu B, Nagaraja K V, et al. 2007. Gauss Legendre-Gauss Jacobi quadrature rules over a tetrahedral region.AppliedMathematicsandComputation, 190(1): 186-194.

Sch?berl J. 1997. NETGEN An advancing front 2D/3D-mesh generator based on abstract rules.ComputingandVisualizationinScience, 1(1): 41-52.Silvester P P, Ferrari R L.1996. Finite Elements for Electrical Engineers. 3rd ed. Cambridge:Cambridgeuniversity Press.

Sun X Y, Nie Z P. 2008. Vector finite element analysis of multicomponent induction response in anisotropic formations.ProgressinElectromagneticsResearch, 81:21-39.

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.ChineseJ.Geophys. (in Chinese), 51(5): 1600-1607.

Teng J W, Yang L Q, Yao J Q, et al. 2007. Deep discover ore、exploration and exploitation for metal mineral resources and its deep dynamical process of formation.ProgressinGeophysics(in Chinese), 22(2):317-334.

Wang H N,Tao H G, Yao J J, et al. 2008. Study on the response of a multicomponent induction logging tool in deviated and layered

anisotropic formations by using numerical mode matching method.ChineseJ.Geophys. (in Chinese), 51(5):1591-1599.Wang H M, Wu P, Rosthal R, et al. 2008. Modeling and understanding the triaxial induction logging in borehole environment with dip anisotropic formation. ∥ 2008 SEG Annual Meeting, Society of Exploration Geophysicists.309-313

Wang J, Chen H, Wang X M, et al. 2015. Research on selection method of background field for finite element simulation of induction logging.ChineseJ.Geophys. (in Chinese),58(6): 2177-2187.

Wei B J, Wang T T, Wang Y. 2009. Computing the response of multi-component induction logging in layered anisotropic formation by the recursive matrix method for magnetic-current-source dyadic Green′s function.ChineseJ.Geophys. (in Chinese),52(11): 2920-2928.Weiss C J, Newman G A. 2002. Electromagnetic induction in a fully 3-D anisotropic earth.Geophysics, 67(4): 1104-1114.Weiss C J. 2013. Project APhiD: A Lorenz-gauged A-Φ decomposition for parallelized computation of ultra-broadband electromagnetic induction in a fully heterogeneous Earth.Computers&Geosciences, 58: 40-52.

Xiao J Q, Zhang G Y, Hong D C, et al .2013.Fast forward modeling and data processing of 3D induction logging tool in layered anisotropic formation.ChineseJ.Geophys. (in Chinese),56(2): 696-706.Zienkiewicz O C, Taylor R L. 2005. The Finite Element Method for Solid and Structural Mechanics. Oxford: Butterworth-H.

附中文參考文獻

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

滕吉文，楊立強，姚敬全等. 2007. 金屬礦產(chǎn)資源的深部找礦，勘探與成礦的深層動力過程.地球物理學進展，22(2):317-334.

汪宏年, 陶宏根, 姚敬金等. 2008. 用模式匹配算法研究層狀各向異性傾斜地層中多分量感應測井響應. 地球物理學報, 51(5):1591-1599.

王健, 陳浩, 王秀明等. 2015. 有限元感應測井模擬的背景場選擇方法研究. 地球物理學報, 58(6):2177-2187.

魏寶君, 王甜甜,王穎. 2009. 用磁流源并矢 Green 函數(shù)的遞推矩陣方法計算層狀各向異性地層中多分量感應測井響應. 地球物理學報, 52(11):2920-2928.

肖加奇, 張國艷, 洪德成等. 2013. 層狀各向異性地層中三維感應測井響應快速計算及資料處理. 地球物理學報, 56(2):696-706.

(本文編輯胡素芳)

基金項目國家重大科研裝備研制項目“深部資源探測核心裝備研發(fā)”(ZDYZ2012-1-07)資助.

作者簡介王健，男，1987年生，博士生，主要從事電磁波場數(shù)值模擬研究. E-mail: wangjianshinian@163.com

doi:10.6038/cjg20160130 中圖分類號P631

收稿日期2015-02-16，2015-10-22收修定稿

Response modeling of multi-component induction logging tool in the mineral logging using vector finite element

WANG Jian1,2， CHEN Hao1, WANG Xiu-Ming1

1StateKeyLaboratoryofAcoustics,InstituteofAcoustics,ChineseAcademyofSciences,Beijing100190，China2UniversityofChineseAcademyofSciences,Beijing100049，China

AbstractWith the development of economy, the shallow resources have been exploited excessively, the deep-seated resources become the investigation targets for the geophysics. Borehole geophysics become more important due to the reduction of resolution and accuracy of the surface methods, and the increase of the cost for drilling and sampling. So how to realize fine evaluation of solid deposit by logging is a crucial problem. This article′s main objective is to determine the relative position of solid mineral deposit and well based on the numerical simulation of multicomponent induction logging method.

KeywordsEdge finite element；Multi-component induction logging；Solid mineral deposit

王健, 陳浩, 王秀明. 2016. 用于固體礦床多分量感應測井響應模擬的矢量有限元法.地球物理學報，59(1):355-367,doi:10.6038/cjg20160130.

Wang J, Chen H, Wang X M. 2016. Response modeling of multi-component induction logging tool in the mineral logging using vector finite element.ChineseJ.Geophys. (in Chinese),59(1):355-367,doi:10.6038/cjg20160130.