王健, 陳浩, 王秀明

1 中國科學(xué)院聲學(xué)研究所聲場與聲信息國家重點(diǎn)實(shí)驗(yàn)室， 北京　100190 2 中國科學(xué)院大學(xué)， 北京　100049

?

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

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

1 中國科學(xué)院聲學(xué)研究所聲場與聲信息國家重點(diǎn)實(shí)驗(yàn)室， 北京100190 2 中國科學(xué)院大學(xué)， 北京100049

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

關(guān)鍵詞矢量有限元； 多分量感應(yīng)測井； 固體礦床

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

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

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

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

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

2基本原理

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

(1)

(2)

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

E=Eb+Es,

(3)

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

(4)

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

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

(5)

(6)

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

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

表1　數(shù)學(xué)符號(hào)

Ax=L Ux=b，

(7)

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

(8)

其中y=L-1b.

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

3精度驗(yàn)證

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

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

(9)

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

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

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

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

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

4實(shí)例分析

4.1傾斜各向異性地層模型

(10)

(11)

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

圖5　多分量感應(yīng)測井儀器的基本結(jié)構(gòu)等效圖

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

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

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

4.2球形礦體模型

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

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

圖9　不同地層傾斜角度的視電導(dǎo)率曲線

圖11　多分量感應(yīng)測井在球形礦體模型中的響應(yīng)

(12)

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

圖10　球形礦體模型

4.3傾斜圓柱形礦體模型

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

5結(jié)論

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

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

圖12　球形礦體方位角計(jì)算

圖13　傾斜圓柱型礦體模型

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

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

圖A1　四面體單元

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

邊端點(diǎn)1端點(diǎn)21ab2ac3ad4bc5bd6cd

(A1)

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

(A2)

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

(A3)

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

(A4)

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

(A5)

(A6)

(A7)

其有離散形式如下：

(A8)

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

(A9)

有關(guān)后處理方法更詳細(xì)的討論參見文獻(xiàn)(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.

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

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

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

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

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

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

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

(本文編輯胡素芳)

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

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

doi:10.6038/cjg20160130 中圖分類號(hào)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. 用于固體礦床多分量感應(yīng)測井響應(yīng)模擬的矢量有限元法.地球物理學(xué)報(bào)，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.