嚴(yán) 波,韓 波

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

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

嚴(yán) 波1,韓 波2

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

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

自適應(yīng)有限元;非結(jié)構(gòu)三角單元;后驗(yàn)誤差估計(jì)因子;傾斜各向異性

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

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

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

1 方法理論

1.1 控制方程和邊界條件

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

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

(1)

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

(2)

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

(3)

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

(4)

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

(5)

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

或者:

(6)

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

(7)

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

(8)

式中:K0為第二類零階修正貝塞爾函數(shù)?；?8)式對y求偏導(dǎo)數(shù),得到:

(9)

(10)

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

(11)

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

(12)

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

(13)

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

1.2 有限元近似

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

(14)

(15)

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

(16)

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

(17)

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

2 基于后驗(yàn)誤差估計(jì)的網(wǎng)格自動(dòng)細(xì)化

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

圖2 二維直流電阻率自適應(yīng)有限元算法流程

3 模型計(jì)算

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

3.1 模型一

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

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

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

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

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

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

3.2 模型二

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

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

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

采用單極裝置,設(shè)電源點(diǎn)位于坐標(biāo)原點(diǎn),接收電極沿著y軸正負(fù)兩個(gè)方向布設(shè),其距坐標(biāo)原點(diǎn)的距離為: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ù)和加權(quán)系數(shù)(表1),自適應(yīng)有限元模擬網(wǎng)

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

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

式中:σyz,σyy,σzz分別表示測點(diǎn)電導(dǎo)率σ在不同方向

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

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

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

3.3 模型三

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

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

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

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

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

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

4 結(jié)論

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

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

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

[2] 馬永生,張建寧,趙培榮,等.物探技術(shù)需求分析及攻關(guān)方向思考——以中國石化油氣勘探為例[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].北京:科學(xué)出版社,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] 嚴(yán)波,劉穎,葉益信.基于對偶加權(quán)后驗(yàn)誤差估計(jì)的2.5維直流電阻率自適應(yīng)有限元正演[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] 阮百堯,徐世浙.電導(dǎo)率分塊線性變化二維地電斷面電阻率測深有限元數(shù)值模擬[J].地球科學(xué),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

嚴(yán)波,韓波.二維直流電阻率傾斜各向異性自適應(yīng)有限元正演[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。

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

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

國家自然科學(xué)基金青年基金項(xiàng)目(41604063)資助。

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