二維直流電阻率傾斜各向異性自適應有限元正演

嚴 波,韓 波

(1.湖南文理學院國際學院,湖南常德415000;2.中國海洋大學“海底科學與探測技術”教育部重點實驗室,山東青島266100)

二維直流電阻率傾斜各向異性自適應有限元正演

嚴 波1,韓 波2

(1.湖南文理學院國際學院,湖南常德415000;2.中國海洋大學“海底科學與探測技術”教育部重點實驗室,山東青島266100)

提出了一種二維直流電阻率傾斜各向異性自適應有限元數(shù)值模擬算法,利用對偶加權后驗誤差估計因子指導網(wǎng)格自動細化過程,并在易于模擬起伏地形和傾斜界面等復雜構造的非結構三角網(wǎng)格上實現(xiàn)了這種算法。模擬了水平兩層電阻率傾斜各向異性模型的直流視電阻率,并將其與解析解以及電阻率各向同性模型的直流視電阻率進行了比較,結果表明,除了電源點附近測點的有限元數(shù)值解誤差較大外,其它測點處相對誤差都在1%以內(nèi),且電阻率各向異性對視電阻率曲線的峰值和形態(tài)都有明顯影響。二維電流密度矢量圖能夠很好地解釋視電阻率曲線峰值和形態(tài)隨電阻率各向異性傾角變化而變化的現(xiàn)象。

自適應有限元;非結構三角單元;后驗誤差估計因子;傾斜各向異性

在石油物探中,電法勘探占有重要地位。特別是在區(qū)域普查中,它可用于研究區(qū)域構造,確定沉積盆地基底起伏情況,圈定含油氣遠景區(qū)等[1-2]。而在進行電法勘探時,發(fā)現(xiàn)地下巖石的導電性常常表現(xiàn)出非常明顯的宏觀各向異性特性,例如板巖和頁巖,它們在層理方向與垂直層理方向的電導率之比達到了1.2～5.0[3-5]。在解釋這些區(qū)域采集到的地電數(shù)據(jù)時,如果忽略了電阻率各向異性的影響,往往會得到不精確甚至錯誤的地電構造信息[6-7],因此,研究電阻率各向異性地電斷面的正演算法具有重要的理論和實際意義。

對于直流電阻率各向異性數(shù)值模擬算法的研究,國內(nèi)外學者已經(jīng)取得了一定的研究成果。WAIT[8],LI等[9]和PERVAGO等[10]推導出了一維層狀任意各向異性介質(zhì)點電源電位的解析解,并分析了各向異性對層狀介質(zhì)直流視電阻率的影響;徐世浙[3]實現(xiàn)了點源二維垂向各向異性地電斷面的直流電場有限元算法,并分析了各向異性對電測深曲線的影響;VERNER等[11]基于規(guī)則四邊形網(wǎng)格離散微分方程,采用狄利克雷邊界條件,實現(xiàn)了二維直流電阻率任意各向異性有限差分數(shù)值模擬算法;BIBBY[12]針對軸對稱模型實現(xiàn)了直流電阻率有限元算法;ZHOU等[13]利用高斯積分網(wǎng)格實現(xiàn)了二維和三維直流電阻率傾斜各向異性正演算法;LI等[4]實現(xiàn)了三維直流電阻率任意各向異性有限元正演算法。以上這些二維、三維數(shù)值解法均基于規(guī)則的結構網(wǎng)格,很難準確模擬地形起伏、傾斜界面等復雜的地質(zhì)結構,為此,WANG等[14]利用基于非結構四面體網(wǎng)格的有限元法,實現(xiàn)了三維直流電阻率任意各向異性正演。非結構網(wǎng)格不僅在模擬復雜地電結構方面相對于結構網(wǎng)格有明顯的優(yōu)勢,還特別適合模擬多尺度構造,比如在大尺度的區(qū)域構造中存在很多小尺度的不均勻體。

本文基于完全非結構三角網(wǎng)格,實現(xiàn)了二維直流電阻率傾斜各向異性自適應有限元數(shù)值模擬算法。算法利用對偶加權后驗誤差估計因子來控制非結構三角網(wǎng)格自動細化過程[15-17],網(wǎng)格剖分過程自動完成,無需人工干預,不僅可以很好地模擬起伏地形和復雜的地質(zhì)構造,還可以大大減小人為網(wǎng)格剖分所導致的誤差,增加了模擬的靈活性,此外,混合邊界條件的使用能進一步提高數(shù)值解的精度。對多個理論模型進行了正演計算,在驗證算法有效性的同時,深入分析了傾斜各向異性對直流電阻率測深的影響。

1 方法理論

1.1 控制方程和邊界條件

在傾斜各向異性介質(zhì)中,電導率σ為張量,建立如圖1所示的坐標系。

在x′y′z′坐標系中(x′與x平行,y′與y軸、z′與z軸的夾角為β),假設沿x′,y′,z′方向上的電導率分別為σx,σy,σz,則電導率張量σ′可以表示為:

(1)

圖1 傾斜各向異性介質(zhì)三維空間坐標系

(2)

當點電源I位于地表,地下介質(zhì)為電阻率各向異性時,電位u滿足的微分方程為[3]:

(3)

式中:δ(A)表示以A為中心的狄拉克函數(shù)。

(4)

(4)式是一個三維微分方程,這里選擇用傅里葉變換將它轉(zhuǎn)換為二維問題。將(4)式關于x方向進行傅里葉變換,得到:

(5)

式中:U為波數(shù)域的電位值。

或者:

(6)

假定邊界?！尬挥跓o窮遠處,使得電阻率各向異性不均勻體在該邊界上的影響近似為零,這時?！奚系碾娢豢山瓢袋c電源的電位確定。假設Γ∞周圍介質(zhì)是均勻各向異性的,且各向異性主軸與坐標軸重合,這時點電源的電位可表示為:

(7)

其中,λy=σx/σy,λz=σx/σz,c為比例系數(shù)。對(7)式進行傅里葉變換:

(8)

式中:K0為第二類零階修正貝塞爾函數(shù)。基于(8)式對y求偏導數(shù),得到:

(9)

(10)

同理,基于(8)式對z求偏導數(shù),并代入比例系數(shù)c,有:

(11)

空氣不導電,地表Γs電流密度J的法向分量為零,故有:

(12)

對(12)式進行傅里葉變換,得到:

(13)

(6)式、(10)式、(11)式和(13)式組成了電阻率傾斜各向異性二維點電源直流電阻率波數(shù)域的邊值問題。

1.2 有限元近似

將邊值問題轉(zhuǎn)化為有限元方程的方法主要有兩種,即加權余量法和變分法[18],前者利用虛功原理,后者利用最小位能原理,兩種方法是等價的。本文利用加權余量法推導電阻率傾斜各向異性二維點電源直流電阻率邊值問題的有限元方程,將等式(6)兩邊乘以轉(zhuǎn)換電位的任意小量δU,并在模型區(qū)域Ω內(nèi)積分,可得:

(14)

(15)

將邊界條件(10)式、(11)式、(13)式代入(15)式,可得到:

(16)

其中,θ為?！薜耐夥ㄏ蚍较蚺cy軸的夾角。將模型區(qū)域Ω離散成若干個三角單元,則(16)式的積分等于所有三角單元的積分總和,形成如下線性方程組:

(17)

式中:U是所有網(wǎng)格節(jié)點處電位值構成的列向量;P=(0,…,0,I/2,0,…,0)T,只有電源點對應的節(jié)點處其值為I/2;系數(shù)矩陣K對稱且高度稀疏。為了節(jié)省內(nèi)存空間,采用按行壓縮存儲方式,并利用直接求解器SuperLU[19]求解線性方程組(17)。對于若干個波數(shù)k,解方程(17)即可求得波數(shù)域的電位U,再對U進行反傅里葉變換[20],即可求得空間域所有網(wǎng)格節(jié)點處的電位u,再利用三角單元中3個節(jié)點處的電位值進行線性插值可以求得接收點處的電位值。根據(jù)不同的測量裝置[21],由接收點處的電位值就可以計算出視電阻率值。

2 基于后驗誤差估計的網(wǎng)格自動細化

本文中用到的基于后驗誤差估計的網(wǎng)格自動細化的基本理論在文獻[15]中已經(jīng)有詳細論述,自適應網(wǎng)格細化有限元正演流程如圖2所示。首先,輸入初始網(wǎng)格和模型參數(shù)數(shù)據(jù),并計算初始網(wǎng)格的有限元近似解;然后,利用有限元近似解計算初始網(wǎng)格中每個三角單元的局部誤差估計值和相應對偶問題的有限元解;最后,計算每個三角單元的后驗誤差估計因子值,對后驗誤差估計因子值比較大的三角單元進行網(wǎng)格細化;再次計算網(wǎng)格細化后的有限元解,并且重復迭代上述過程,直到有限元數(shù)值解收斂為止。

圖2 二維直流電阻率自適應有限元算法流程

3 模型計算

下面的模型計算中,設定所有模型的走向方向均為x方向,采用二維網(wǎng)格剖分器triangle[22]剖分網(wǎng)格,每次網(wǎng)格細化的三角單元數(shù)占三角單元總數(shù)的20%,網(wǎng)格細化的終止條件是相鄰兩次有限元解的相對誤差小于1%。所有計算均在3.20GHz CPU,4G RAM的個人計算機上完成。

3.1 模型一

圖3 兩層電阻率各向同性模型

圖4 兩層電阻率各向異性模型

采用二極裝置,設電源點位于坐標原點。接收電極沿著y軸正方向布置,極距從1m逐漸增加到50m,共有20個接收電極。為了選擇適合本模型的離散化波數(shù),我們采用阮百堯等[23]提出的適合單個點電源情形的波數(shù)和加權系數(shù)進行計算,如表1所示。模擬網(wǎng)格大小為300m×800m,經(jīng)過11次迭代后,網(wǎng)格剖分細化終止,有限元數(shù)值解收斂。具體的網(wǎng)格剖分信息和計算結果如圖5和圖6所示。

表1 單點電源情形下所用波數(shù)及加權系數(shù)

注:ki代表波數(shù);gi代表加權系數(shù)。

從圖5b可以看出,電源點處和接收點處的網(wǎng)格加密程度要明顯高于其它區(qū)域,實現(xiàn)了網(wǎng)格局部加密的目的,這不僅可以使得有限元數(shù)值解快速收斂于精確解,還能減少計算資源要求,提高計算效率。從圖6 可以看出,兩層電阻率各向同性和各向異性模型的二維直流電阻率數(shù)值解都與解析解[8]吻合得很好,除電源點附近數(shù)值解的相對誤差較大外,其它接收點處的相對誤差都在1%以內(nèi)。

3.2 模型二

假定在ρ0=5Ω·m的低阻均勻介質(zhì)中埋藏有一個電阻率傾斜各向異性正方體,并假定其在走向x方向上無限延伸,y方向和z方向上的邊長都為5m,其電阻率為ρxx=ρzz=100Ω·m,ρyy=5Ω·m,不均勻體距離水平地面的距離為0.5m,如圖7所示。

圖5 兩層電阻率各向異性模型初始網(wǎng)格(a)與最終網(wǎng)格局部放大結果(b)

圖6 解析解與自適應有限元結果(a)及其相對誤差(b)

采用單極裝置,設電源點位于坐標原點,接收電極沿著y軸正負兩個方向布設,其距坐標原點的距離為:0.45,0.90,1.35,1.80,2.25,2.70,3.15,3.60,4.05,4.50,4.95,5.40m。我們依然采用阮百堯等[23]提出的波數(shù)和加權系數(shù)(表1),自適應有限元模擬網(wǎng)

圖7 電阻率傾斜各向異性模型

格大小為300m×600m。分別計算電阻率各向異性傾角β=30°,45°,60°,75°情形下的視電阻率測深曲線,結果如圖8所示。從圖8可以看出,當電阻率各向異性傾角不同時,視電阻率曲線異常的峰值會發(fā)生變化,隨著電阻率各向異性傾角的增大,視電阻率曲線的峰值會變小,可以用地下電流密度的分布特征來解釋這種現(xiàn)象。在二維傾斜各向異性介質(zhì)中,電流密度是電場的橫向分量和縱向分量的線性組合[24]。這時,σyx=σyz=0,yz平面內(nèi)y,z方向上電流密度分別為:

式中:σyz,σyy,σzz分別表示測點電導率σ在不同方向

圖8 不同電阻率各向異性傾角時的視電阻率曲線a β=30°; b β=45°; c β=60°; d β=75°

將接收電極置于地面以下,并按照(18)式和(19)式分別計算電流密度y方向和z方向上的分量。

圖9顯示了當電阻率各向異性傾角β=30°,45°,60°,75°時,電流密度矢量在yz平面內(nèi)的分布情況。因為電流密度在遠離電源點的地方衰減很快,電流密度的分布與r2(r表示測點與點電源I之間的距離)成反比,所以,在遠離電源點的區(qū)域電流密度值比較小。電流密度主要沿著低阻介質(zhì)流動,其大致流動方向和低阻地層傾斜的方向一致,角度約等于傾角,所以,隨著低阻地層傾角的增大,地面測得的電位差就會變小,計算出來的視電阻率值也就會隨之變小。

3.3 模型三

假設在ρ0=100Ω·m的地下均勻介質(zhì)中埋藏有一個低阻傾斜各向異性長方體,它在x方向上無限延伸,y方向和z方向上的邊長分別為6m和2m,其電阻率為ρxx=ρyy=4Ω·m,ρzz=25Ω·m,長方體距離水平地面的距離為2m,如圖10所示。利用PIDLISECKY等[25]的最優(yōu)化波數(shù)計算程序,獲得了適用于極距范圍0.25～220.00m,雙點電源情形下的7個波數(shù)和加權系數(shù)(表2)。

表2 雙點電源情形下所用波數(shù)和加權系數(shù)

注:ki代表波數(shù);gi代表加權系數(shù)。

圖10 二維電阻率傾斜各向異性模型

圖11顯示了當間距因子n=5時偶極裝置不同各向異性傾角情形的視電阻率曲線,各向異性長方體的主軸電阻率ρxx,ρyy,ρzz為4,4,25Ω·m,各向異性傾角β=0,30°,60°,90°,電源電極和測量電極距為AB=MN=2m,測量范圍為-20～20m。從圖11 可以看出,各向異性傾角β對視電阻率的影響很大。盡管二維模型是對稱的,但當β=30°,60°時,視電阻率曲線卻不對稱,視電阻率曲線的極小值偏離了中心位置并偏向了左邊,這時的視電阻率曲線與向右傾斜的電阻率各向同性長方形板體的視電阻率曲線非常相似。

圖11 視電阻率曲線(偶極裝置)

4 結論

本文基于完全非結構三角網(wǎng)格,實現(xiàn)了二維直流電阻率傾斜各向異性自適應有限元數(shù)值模擬算法。利用對偶加權后驗誤差估計因子控制非結構三角形網(wǎng)格自動細化過程,網(wǎng)格剖分過程自動完成,無需人工干預,不僅可以很好地模擬起伏地形和復雜的地質(zhì)構造,還可以大大減小人為網(wǎng)格剖分所導致的誤差。

電阻率各向異性傾角不僅會影響視電阻率曲線的峰值大小,還會影響視電阻率曲線的對稱性,從而對異常體的形態(tài)和分布位置的解釋都產(chǎn)生影響,所以,在地質(zhì)解釋時考慮地下介質(zhì)電阻率各向異性的影響很有必要。利用電流密度矢量圖可以方便地解釋二維傾斜各向異性模型視電阻率曲線的特征。

[1] 王家映.石油電法勘探[M].北京:石油工業(yè)出版社,1992:1-4 WANG J Y.Electrical prospecting for petroleum[M].Beijing:Petroleum Industry Press,1992:1-4

[2] 馬永生,張建寧,趙培榮,等.物探技術需求分析及攻關方向思考——以中國石化油氣勘探為例[J].石油物探,2016,55(1):1-9 MA Y S,ZHANG J N,ZHAO P R,et al.Requirement analysis and research direction for the geophysical prospecting technology of SINOPEC[J].Geophysical Prospecting for Petroleum,2016,55(1):1-9

[3] 徐世浙.地球物理中的有限單元法[M].北京:科學出版社,1994:131-178 XU S Z.The finite element method in geophysics[M].Beijing:Science Press,1994:131-178

[4] LI Y G,SPITZER K.Finite element resistivity modelling for three-dimensional structures with arbitrary anisotropy[J].Physics of the Earth & Planetary Interiors,2005,150(1):15-27

[5] 唐杰,方兵,孫成禹,等.基于各向異性流體替換的裂隙介質(zhì)波傳播特征研究[J].石油物探,2015,54(1):1-8 TANG J,FANG B,SUN C Y,et al.Study of seismic wave propagation characteristics based on anisotropic fluid substitution in fractured medium[J].Geophysical Prospecting for Petroleum,2015,54(1):1-8

[6] ASTEN M W.The influence of electrical anisotropy on mise-a-la-masse surveys[J].Geophysical Prospecting,1974,22(2):238-245

[7] LI Y G,DAI S K.Finite element modelling of marine controlled-source electromagnetic responses in two-dimensional dipping anisotropic conductivity structures[J].Geophysical Journal International,2011,185(2):622-636

[8] WAIT J R.Current flow into a 3-dimensionally anisotropic conductor[J].Radio Science,1990,25(5):689-694

[9] LI P,UREN N F.Analytical solution for the point source potential in an anisotropic 3-D half-space I:two-horizental-layer case[J].Mathematical & Computer Modelling,1997,26(5):9-27

[10] PERVAGO E,MOUSATOV A,SHEVNIN V.Analytical solution for the electric potential in arbitrary anisotropic layered media applying the set of Hankel transforms of integer order[J].Geophysical Prospecting,2006,54(5):651-661

[11] VERNER T,PEK J.Numerical modelling of direct currents in 2-D anisotropic structures[C].Deutsche Geophysikalische Gesellschaft Expanded Abstracts,1998:228-237

[12] BIBBY H M.Direct current resistivity modeling for axially symmetric bodies using finite element method[J].Geophysics,1978,43(3):550-562

[13] ZHOU B,GREENHALGH M,GREENHALGH S A.2.5-D/3-D resistivity modelling in anisotropic media using Gaussian quadrature grids [J].Geophysical Journal International,2009,176(1):63-80

[14] WANG W,WU X P,SPITZER K.Three-dimensional DC anisotropic resistivity modelling using finite elements on unstructured grids[J].Geophysical Journal International,2013,193(2):734-746

[15] 嚴波,劉穎,葉益信.基于對偶加權后驗誤差估計的2.5維直流電阻率自適應有限元正演[J].物探與化探,2014,38(1):145-150 YAN B,LIU Y,YE Y X.2.5D direct current resistivity adaptive finite-element numerical modeling based on dual weighted posteriori error estimation[J].Geophysical and Geochemical Exploration,2014,38(1):145-150

[16] LI Y G,KERRY K.2D marine controlled-source electromagnetic modeling:part 1—an adaptive finite-element algorithm[J].Geophysics,2007,72(2):WA51-WA62

[17] LI Y G,JOSEF P.Adaptive finite element modelling of two-dimensional magnetotelluric fields in general anisotropic media[J].Geophysical Journal International,2008,175(3):942-954

[18] ZIENKIEWICZ O C,TAYLOR R L.The finite-element method[M].5thed.Oxford:Butterworth-Heinemann,2000:54-97

[19] DEY A,Morrison H F.Resistivity modelling for arbitrarily shaped three-dimensional structures[J].Geophysics,1979,44(4):753-780

[20] XU S Z,DUAN B C,ZHANG D H.Selection of the wave numbers k using optimization method for the inverse Fourier transform in 2.5D electrical modeling[J].Geophysical Prospecting,2000,48(5):789-796

[21] 李金銘.地電場與電法勘探[M].北京:地質(zhì)出版社,2005:137-145 LI J M.Geo-electric field and electric exploration[M].Beijing:Geological Publishing House,2005:137-145

[22] SHEWCHUK J R.Delaunay refinement mesh generation[D].USA:Carnegie Mellon University,1997

[23] 阮百堯,徐世浙.電導率分塊線性變化二維地電斷面電阻率測深有限元數(shù)值模擬[J].地球科學,1998,23(3):303-307 RUAN B Y,XU S Z.FEM for modeling resistivity sounding on 2-D geo-electric model with line variation of conductivity within each block[J].Earth Science,1998,23(3):303-307

[24] YIN C,MAURER H W.Electromagnetic induction in layered earth with arbitrary anisotropy [J].Geophysics,2001,66(5):1405-1415

[25] PIDLISECKY A,KNIGHT R.FW2_5D:a MATLAB 2.5-D electrical resistivity modeling code[J].Computers & Geosciences,2008,34(12):1645-1654

(編輯：陳 杰)

Adaptivefiniteelementmodelingofdirectcurrentresistivityin2-Dtiltedanisotropystructures

YAN Bo1,HAN Bo2

(1.Internationalcollege,HunanUniversityofArtsandScience,Changde415000,China;2.KeyLabofSubmarineGeosciencesandProspectingTechniquesofMinistryofEducation,OceanUniversityofChina,Qingdao266100,China)

In this paper,we propose an adaptive finite element (FE) numerical simulation algorithm for direct current (DC) resistivity in two-dimensional tilted anisotropy structures.This algorithm uses a dual weighting posteriori error estimation factor to guide the mesh automatic refining process,and it is implemented on unstructured triangular elements that facilitate simulation of complex structures characterized by undulating topography and tilted interfaces.We simulated the DC apparent resistivity of a two-layer resistivity tilted anisotropy model and compared it with the analytic solution and the DC apparent resistivity of the resistivity isotropy model.The results showed that the relative error of the finite element numerical solution in most measuring points was within 1%,the relative error of the measuring points near the power points was large,and the resistivity anisotropy had a significant effect on the peak and shape of the apparent resistivity curves.The phenomenon of the peak and the shape of the resistivity curves varying with the change in resistivity anisotropy dip angles can be explained using a two-dimensional current density vector diagram.

adaptive finite element,unstructured triangular element,posteriori error estimation factor,tilted anisotropy

P631

：A

1000-1441(2017)05-0766-08DOI:10.3969/j.issn.1000-1441.2017.05.017

嚴波,韓波.二維直流電阻率傾斜各向異性自適應有限元正演[J].石油物探,2017,56(5):773

YAN Bo,HAN Bo.Adaptive finite element modeling of direct current resistivity in 2D tilted anisotropy structures

[J].Geophysical Prospecting for Petroleum,2017,56(5):773

2015-12-16;改回日期：2017-06-27。

嚴波(1986—),男,博士,主要從事二維直流電阻率任意各向異性自適應有限元正演算法研究。

韓波(1987—),男,博士,主要從事海洋可控源電磁與大地電磁三維數(shù)據(jù)空間法聯(lián)合反演研究工作。

國家自然科學基金青年基金項目(41604063)資助。

This research is financially supported by the National Natural Science Foundation of China (Grant No.41604063).