基于徑向剖分的井地三維有限元電阻率與極化率模擬

黃俊革, 李林杰, 農觀海,2, 閆懷超, 高文利, 林振洲

(1.上海應用技術學院 城市建設與安全工程學院，上海 201418;

2.桂林理工大學 地球科學學院，桂林 541004;

3.中國地質科學院 地球物理地球化學勘查研究所，河北 廊坊 065000)

?

基于徑向剖分的井地三維有限元電阻率與極化率模擬

黃俊革1, 李林杰1, 農觀海1,2, 閆懷超1, 高文利3, 林振洲3

(1.上海應用技術學院 城市建設與安全工程學院，上海 201418;

2.桂林理工大學 地球科學學院，桂林 541004;

3.中國地質科學院 地球物理地球化學勘查研究所，河北 廊坊 065000)

[摘要]采用徑向剖分的有限元法對井地三維電阻率法進行數值模擬，以適應不同觀測裝置對數值模擬的要求。以井口為中心，將區域在環鉆孔方向以45°夾角均勻剖分為八等分，再在直徑方向上等距離剖分。徑向剖分形成的單元呈環形分布，更符合井中電場的分布規律。與這種剖分相對應，設計了以井口為中心布置4條放射狀測線(“米”字形)的井地觀測裝置，解決了井旁地質體位置未知的情況下難以確定觀測方向的問題。通過對典型的地質模型進行有限元正演模擬計算，結果表明徑向剖分數值模擬方法是可行的、準確的；以這種井地觀測方式探測位置未知的地質體是行之有效的。

[關鍵詞]井地；徑向剖面；有限元；電阻率； 極化率

傳統的井中電阻率測深大都采用單方向或向量激發的方式通過平行測線進行觀測，在探測井旁盲礦、油氣藏邊界探查等領域得到廣泛的應用[1]。基于勘查需求和應用領域的拓展，推動了井中電阻率測深的理論研究，如Wang等、Avdeev 等分別用有限差分和積分方程法對電阻率測井進行了三維模擬[2,3]；王志剛等對井地直流電法三維數值模擬及異常規律研究[4]。劉雪軍等應用正演模擬方法研究供電電極在井中針對儲集目標供電激發時，地下與地面電磁場的異常特征[1]；戴前偉等利用雙頻激電井地電位技術研究剩余油分布[5]；毛立峰等對井-地交流電法的矢量有限元三維正演問題進行了研究[6]；何展翔等在圈定油氣藏、研究儲層分布、油井的注水和注漿的動態監測方面均取得了一定的效果[7]；李長偉等研究了不同測量方式( 井-地，地-井，井-井) 下點源場井中電法的三維有限元數值模擬[8,9]；呂玉增等研究了地井五方位測量方法及正反演[10,11]；岳建華等研究了井-地三維電阻率成像技術及實現方法[12]。這些方法中，大都以地質體所在方位為主要探測方向，以井口為中心布置平行測網。在這樣的探測方式中，地質體位置與探測異常的幅值有很大的關系，在探測之前必須確定正確的觀測方向，才可以有效突出地質體異常，減少干擾異常；但在地下地質體位置未知的情況下，很難預先確定測線的布置方向。為此我們設計了以井口為中心布置4條放射狀測線的井地徑向觀測裝置，井中供電，地面測量相鄰2個電極間的電位差，這樣可以有效解決井旁盲礦體位置未知的情況下難以確定觀測方向的問題，提高對井旁地質體的勘探能力，在探測地下水流向的工作中這種電極布設方式被廣泛采用。針對這種觀測方式，井地電法的數值模擬進行區域離散時，難以將其離散成六面體單元，且網格疏密不易控制，在剖分單元的3個方向上長度不能同時放大，造成網格剖分困難[8]。本文針對井中電場分布的特點和井地徑向觀測裝置，設計了井地電場基于徑向剖分的三維電阻率法有限元數值模擬方法，對目標區域進行統一的四面體剖分，以適應徑向剖分多種網格單元并存的情況。通過對典型模型模擬計算并分析其異常特征, 結果表明該方法效果良好。

1井地徑向電阻率、極化率觀測方式

一般來說，為探測地下地質體，獲得最大觀測異常，應首先確定地下目標體的走向等信息，再確定觀測主剖面、布置測網。井地徑向三維電阻率、極化率測量方法使用的是單極—偶極裝置,其基本思想是在井下以滾動形式進行點電源供電，地面上以井口為中心布置4條放射狀(“米”字形)測線(圖1)，測量測線上相鄰電極之間的電位差，得到地下目標體視電阻率、視極化率信息。這種測量方式，無須預先確定觀測主方位，便可得到地下地質體分布信息。

圖1　井地徑向剖面三維電阻率觀測方式示意圖Fig.1　Sketch of borehole-ground radial section 3-D resistivity survey method

2區域離散與單元分析

以有限單元法對區域進行三維模擬計算時，大都將觀測區域均勻剖分成六面體，遠離井口的大規模地質體需在3個方向(X、Y、Z)的剖分中占用很多剖分單元，剖分的適應性不佳。為了提高井地電法模擬的單元適應能力，我們采用以井口為中心，在直徑方向(R方向)上、深度方向上(Z方向)等距離剖分，在環鉆孔方向以45°夾角均勻剖分為八等分(圖2)。采用這樣的網格單元剖分方法對區域進行剖分，一方面使單元體積與單元到井孔的距離呈正比，即單元距離井口越遠，單元體積越大，同樣的剖分區域占用的網格數量大幅減少；另一方面，由于區域剖分單元以井口為中心環形分布，單元的一個邊與電場等位線基本重合，與地井電場分布的實際情況更加吻合，提高了該單元節點之間的電位插值精度，從而可以提高正演的精度和速度。

圖2　三維空間剖分示意圖Fig.2　Sketch of 3-D spatial division(A)目標區網格剖分俯視圖;(B)目標區網格剖分側視圖

將區域進行剖分時，與井口相鄰的單元被剖分成三棱柱，其他單元被剖分成水平截面為梯形的六面體(圖3)。為了將網格剖分和算法統一，我們將每個六面體單元首先剖分成2個三棱柱單元, 再對每個三棱柱單元按照右手螺旋規則剖分成4個四面體，單元剖分規則和節點編號如圖3所示。這樣就可以實現在不增加節點數又能得到盡量多的剖分單元,從而達到提高有限單元計算精度的目的。

下面我們以節點編號為1, 2, 3, 5的四面體單元說明單元節點電位求解的過程。假設1, 2, 3, 5節點坐標分別為(x1,y1,z1), (x2,y2,z2), (x3,y3,z3), (x4,y4,z4);p(x,y,z)為四面體內部的一個點(圖4)，p點的位置可寫為以下4個表達式[13]

Li是x,y,z的線性函數，例如

其中

圖3　模型單元剖分示意圖Fig.3　Sketch of the model elements division(A)三棱柱體單元剖分示意圖; (B)六面體單元剖分示意圖

圖4　四面體單元示意圖Fig.4　Sketch of tetrahedron elements

是與四面體頂點坐標有關的常數。類似地，可得到L2,L3,L4。

單元中電位采用線性插值，u=α1x+α2y+α3z+α4,根據4個頂點的坐標和函數值，可確定4個系數α1, …,α4,可寫成

其中形函數Ni=Li(i=1,…,4)。

為提高計算精度,采用三維電場的異常電位法求解[13],單元的電位列向量ue擴展成全體節點的電位列向量u，可以得到以下方程組

Ku=-K′u0

(1)

其中K和K′為系數矩陣，u為待求異常電位向量，u0為正常電位向量。井-地觀測方式均勻半空間中的正常電位表達式為

(2)

式中的r′為鏡像電源到觀測點之間的距離。由于觀測點位于地表，顯然有r=r′，解線性方程組(1)，可得各結點的異常電位u，與預先計算的各節點正常電位u0相加，便可得到節點的總電位[14]。

3視極化率的計算

(3)

故

ηs=1-ρs2/ρs1

(4)

其中:ΔU為極化場電位差; ΔU1為一次場電位差; ΔU2為二次場電位差; ηs表示視極化率。

4數值算例與結果分析

4.1井旁低阻高極化模型

圖5為一個直立井模型，井深20m。模型中井旁有一個棱柱地質體(位于南北方向測線第5斷面下部，地質體的尺寸以及形狀見圖6)。其電阻率為5Ω·m，極化率為20%；圍巖電阻率為100Ω·m，極化率為5%；地質體中心與井孔間的距離4m，距離地面6m。地面上以井口為中心布置4條放射狀測線(圖1)，供電點依次由淺及深(-1～-13m)供電，測量地面相鄰兩個電極的電位差，測線上2個測量電極之間的距離為1m。視電阻率、視極化率斷面圖成圖規則是，以供電點深度為Z坐標，觀測點中心坐標為X、Y坐標，共形成南北(5-1)、南西-北東(6-2)、西東(7-3)、西北-東南(8-4)方向4個斷面。

圖5　井旁低阻高極化模型示意圖Fig.5　Schematic diagram of model 1

圖6　異常體三維模型Fig.6　3-D model anomalous body

從圖7(A)、(B)的低阻棱柱視電阻率南-北、東北-西南方向剖面可見，鉆井兩側附近出現一對高幅值、直立脈沖形異常，鉆井靠異常體一側的脈沖異常為低阻，與地質體電性對應,另一側為高阻。低阻脈沖異常頂部深度約為6m，呈下窄上寬的形態，低阻異常上部的寬度與地質體寬度基本一致。從圖7(C)、(D)的視電阻率西-東、北西-南東方向剖面的異常特征圖可見，除在井孔的兩側附近出現一對直立脈沖形異常外，脈沖異常的頂部出現一個扁平馬鞍形低阻異常，與低阻地質體在斷面上的投影位置和電性基本對應。總體來說，各方向斷面均可見低阻體異常，且異常位置與地質體位置或地質體在斷面上的投影位置相對應。

圖7　低阻高極化地質體井地徑向剖分數值模擬視電阻率斷面Fig.7　Apparent resistivity sections on low resistivity & high polarizability geological body borehole-ground radially split numerical simulation

從視極化率剖面圖(圖8)可以看到，井孔兩側同樣出現脈沖形高、低視極化率異常，高極化異常位于地質體一側或與地質體投影位置相對應，最大異常達到9.8%，呈上寬下窄之勢，異常外邊緣與地質體外邊緣對應良好。在各個方向的斷面中均見極化率異常，位置對應良好，異常值大小與觀測斷面和地質體間距離相關。

為了檢驗計算結果是否準確，計算精度是否滿足要求，利用較為成熟的井地三極測深裝置六面體剖分數值模擬解與井地徑向剖分的數值模擬結果進行對比。如圖9(A)所示，低阻立方體棱長4m，電阻率5Ω·m，位于南北方向測線的下方，中心距離井孔4m，背景電阻率100Ω·m，供電點位于井中-1～-13m，觀測點極距1m。對比圖9(B)所示的視電阻率斷面和圖7(A)，異常形態基本相似。從表1的異常幅值的差異來看，均小于5%，兩種模擬結果有差異主要是由于地質體體積差異所致。

4.2井旁高阻高極化模型

井旁高阻高極化模型的異常體位置、形態與井旁低阻高極化模型完全相同，只在電性上不同。異常地質體的電阻率為1 000Ω·m，極化率20%；圍巖電阻率100Ω·m，極化率5%。

從模型的視電阻率斷面(圖10)來看，井孔深部兩側有2個縱向、性質相反的條帶狀異常，高阻異常位于南側，與鉆孔南側地質體的性質相對應，異常頂部深度與地質體頂部基本吻合，但異常的位置與地質體對應不佳。

從視極化率斷面來看，井孔的深部兩側出現了2組高低極化條帶狀異常，高極化異常位于高極化地質體一側。除東西方向的視極化率斷面(圖11-C)與地質體的位置對應不佳外，其余斷面的高極化與地質體的位置對應良好。這種情況，與東西方向的觀測斷面與地質體距離較大有關，而這一特征，更有利于地下地質體的定位。

圖8　低阻高極化地質體井地徑向剖分數值模擬視極化率斷面Fig.8　Apparent polarizability sections on low resistivity high polarizability geological & body borehole-ground radially split numerical simulation

圖9　六面體模型及視電阻率斷面圖Fig.9　Schematic diagram of cube geological model and apparent resistivity section

觀測點位視電阻率/Ω·mXYZ徑向剖分六面體剖分相對誤差/%觀測點位視電阻率/Ω·mXYZ徑向剖分六面體剖分相對誤差/%0-18.5-3102.4103.5-1.10-8.5-389.486.63.20-17.5-3102.1103.3-1.20-7.5-387.084.03.50-16.5-3101.6102.9-1.30-6.5-385.483.81.90-15.5-3101.0102.4-1.40-5.5-385.586.4-1.00-14.5-3100.2101.7-1.50-4.5-387.791.0-3.70-13.5-399.2100.6-1.40-3.5-391.395.7-4.50-12.5-397.999.0-1.10-2.5-395.399.6-4.60-11.5-396.396.8-0.50-1.5-399.1102.4-1.30-10.5-394.393.90.50-0.5-3103.7106.0-2.20-9.5-392.090.31.800.5-399.495.14.4

4.3井孔穿過低阻高極化體模型

如圖10(A)所示，模型中井孔穿過八棱柱低阻體(低阻地質體的尺寸以及形狀參見圖10-B)。八棱柱上表面距離地面4 m，縱向延伸4 m，低阻

圖10　高阻高極化地質體井地徑向剖分數值模擬視電阻率斷面Fig.10　Apparent resistivity sections of high resistivity and polarizability geological body borehole-ground radially split numerical simulation

八棱柱電阻率5 Ω·m，極化率為20%；圍巖電阻率為100 Ω·m，極化率為5%。這種模型主要針對地下水流速測定時，井孔中注以鹽水等低電阻率模型而設計，主要勘察低阻高極化率物質被注入時視電阻率、視極化率的異常。

由于模型呈軸對稱，各個方向的視電阻率、視極化率斷面完全相同，因此僅呈現南北方向斷面即可(圖11-A、B)。圖11(A)的視電阻率斷面圖顯示，低阻馬鞍狀異常在寬度、深度方向的延展范圍與低阻體形狀基本一致，異常形態在水平方向上以井孔為中心對稱分布。圖11(B)所示的視極化率異常，在井孔周圍出現一個圓錐形高低極化異常，上部中心為低極化異常，最大幅值-22%，

與異常體中心位置對應；兩翼為10%左右的中極化異常，寬度與地質體吻合，呈對稱分布；下部為以井孔為中心的高極化異常，最高異常幅值達到30%。 通過視電阻率和視極化率異常的位置與性質，可以推斷地質體空間分布和電性。

4.4井孔穿過非對稱形式的低阻高極化體模型

圖12所示模型，為測定地下水流速時，井中注入低阻高極化率介質一段時間后，形成的低阻高極化率模型。上表面距離地面為4 m，縱向延伸6 m，低阻八棱柱電阻率10 Ω·m，極化率為10%；圍巖電阻率為100 Ω·m，極化率為5%。

從圖13和圖11的視電阻率斷面圖對比來看，圖13(A)的視電阻率異常在井孔南北兩側呈明顯的非對稱現象，井孔的南側(-5～0 m)、視深度-4～-9 m之間存在明顯的低阻異常，深部在井孔的南北兩側出現條帶狀異常。圖13(B)的視極化率異常同樣出現不對稱現象，井孔的南側(-5～0 m)、視深度-4～-9 m之間存在明顯的高極化率異常，深部在井孔的南北兩側出現條帶狀異常。從視電阻率和視極化率的斷面圖中可以很好地判斷水流的方向。

圖11　高阻高極化地質體井地徑向剖分數值模擬視極化率斷面Fig.11　Apparent polarizability sections of high resistivity & high polarizability geological body borehole-ground radially split numerical simulation

圖12　模型示意圖Fig.12　Schematic diagram of models

圖13　南-北方向視電阻率、視極化率斷面Fig.13　Apparent resistivity and polarizability sections in south-north direction(A)視電阻率斷面圖; (B)視極化率斷面圖

圖14　模型示意圖Fig.14　Schematic diagram of a model

圖15　南-北方向視電阻率、視極化率斷面Fig.15　Apparent resistivity and polarizability sections in south-north direction(A)視電阻率斷面圖; (B)視極化率斷面圖

5結 論

本文給出一種適合井地電阻率、極化率觀測的三維徑向剖分的有限單元法數值模擬技術，通過與井地三極測深裝置六面體剖分數值模擬結果進行比較，證明井地徑向剖分數值模擬方法計算結果正確，滿足異常分析的要求。通過對幾個典型的地質模型的正演模擬計算, 分析其視電阻率和視極化率異常的分布規律，結果表明井地徑向觀測系統在探測位置未知的地質體是可行、有效的，尤其在地下水流向觀測方面，徑向剖分的數值模擬方式與觀測裝置的吻合性更好。

[參考文獻]

[1] 劉雪軍,王家映,何展翔,等.研究油氣儲集目標的井中地面電磁新技術[J].勘探地球物理進展,2006,29(2):98-101.

Liu X J, Wang J Y, He Z X,etal. Study of hydrocarbon accumulation by borehole-ground EM method[J]. Progress in Exploration Geophysics, 2006, 29(2): 98-101. (In Chinese)

[2] Wang T, Signorelli. Finite difference modeling of electromagnetic tool response for logging while drilling[J]. Geophysics, 2004, 69(1): 152-160.

[3] Avdeev D B, Kuvshinov A V, Pankratov O V,etal. Three dimensional induction logging problems, part 1: Anintegral equation solution and model comparisons[J]. Geophysics, 2002, 67(2): 413-426.

[4] 王志剛,何展翔,劉昱.井地直流電法三維數值模擬及異常規律研究[J].工程地球物理學報,2006,3(2):87-92.

Wang Z G, He Z X, Li Y. Research of three-dimensional modeling and anomalous rule on borehole-ground dc method[J]. Chinese Journal of Engineering Geophysics, 2006, 3(2): 87-92. (In Chinese)

[5] 戴前偉,陳德鵬,劉海飛,等.雙頻激電井地電位技術研究剩余油分布[J].地球物理學進展,2009,24(3):959-964.

Dai Q W, Chen D P, Liu H F,etal. Research of the residual oil distribution with dual frequency induced polarization and the borehole-to-surface potential method[J]. Progress in Exploration Geophysics, 2006, 29(2): 98-101. (In Chinese)

[6] 毛立峰,王緒本,何展翔.矢量有限元法在井-地交流電法三維正演問題中的應用[C]//中國地球物理學會第22屆年會論文集.成都：四川科學技術出版社,2006:703-704.

Mao L F, Wang X B, He Z X. Application of vector finite-element method in 3D modeling on borehole-ground AC method[C]//22th Annual Conference of Chinese Geophysical Society. Chengdu: Sichuan Science and Technology Press, 2006: 703-704. (In Chinese)

[7] 何展翔,劉雪軍,裘尉庭,等.大功率井-地電法油藏邊界預測技術及效果[J].石油勘探與開發,2004,1(5):74-76.

He Z X, Liu X J, Qiu W T,etal. High-power surface-borehole electrical method in predicting reservoir boundary and its application[J]. Petroleum Exploration and Development, 2004, 31(5): 74-76. (In Chinese)

[8] 李長偉,熊彬,呂玉增.電法測井的三維有限元模擬[J].物探與化探，2012,36(4)：585-590.

Li C W, Xiong B, Lyu Y Z. Three-dimensional finite element modeling of electrical well logging[J]. Geophysical and Geochemical Exploration, 2012, 36(4): 585-590. (In Chinese)

[9] 李長偉,阮百堯,呂玉增,等.點源場井-地電位測量三維有限元模擬[J].地球物理學報,2010,53(3):729-736．

Li C W, Ruan B Y, Lyu Y Z,etal. Three-dimensional FEM modeling of point source borehole-ground electrical potential measurements[J]. Chinese Journal of Geophysics, 2010, 53(3): 729-736．(In Chinese)

[10] 呂玉增，阮百堯.地-井五方位IP異常特征與解釋[C]//第九屆中國國際地球電磁學術討論會論文集.桂林：桂林工學院,2009:177-180.

Lyu Y Z, Ruan B Y. Anomalous characteristics and interpretation of surface-borehole 5-direction IP[C]//9th China International Geo-electromagnetic Workshop. Guilin: Guilin University of Technology, 2009: 177-180. (In Chinese)

[11] 呂玉增,阮百堯,彭蘇萍.地-井方位激電觀測異常特征研究[J].地球物理學進展,2012,27(1)：201-216.

Lyu Y Z, Ruan B Y, Peng S P. A study on anomaly of surface-borehole direction induced polarization survey[J]. Progress in Geophysics, 2012, 27(1): 201-216. (In Chinese)

[12] 岳建華,劉志新.井-地三維電阻率成像技術[J].地球物理學進展，2005,20(2):407-411.

Yue J H, Liu Z X. Three dimension resistivity tomography of mine-ground[J]. Progress In Geophysics, 2005, 20(2): 407-411. (In Chinese)

[13] 徐世浙.地球物理中的有限單元法[M].北京:科學出版社,1994.

Xu S Z. The Finite Element Methods in Geophysics[M]. Beijing: Science Press, 1994. (In Chinese)

[14] 黃俊革，阮百堯，王家林，等.鉆井-地表電極聯合電阻率觀測裝置的異常特征研究[J].地球物理學報, 2009,52(5):1348-1362.

Huang J G, Ruan B Y, Wand J L,etal. A study on anomaly of borehole-to-ground joint resistivity surveying system [J]. Chinese Journal of Geophysics, 2009, 52(5):1348-1362. (In Chinese)

3-D FEM modeling of borehole-ground resisitivity and

polarizability based on radial subdivision

HUANG Jun-Ge1, LI Lin-jie1, NONG Guan-hai1,2,

YAN Huai-chao1, GAO Wen-Li3, LIN Zhen-Zhou3

1.SchoolofConstructionandSafetyEngineering,ShanghaiInstituteofTechnology,Shanghai201418,China;

2.CollegeofEarthScience,GuilinUniversityofTechnology,Guilin541004,China;

3.InstituteofGeophysicalandGeochemicalExploration,CAGS,Langfang, 065000,China

Abstract:In order to fit different detecting installations, this paper makes 3-D modeling of the borehole-ground resisitivity and polarizability based on the radial division. Taking the well mouth as the centre, the region is equably divided into 8 sectors by 4 radial lines on the ground with the angle of 45° between two lines. Each sector is divided by a regular distance on radial line. The radial subdivision cells present a circular form, which is more coincidental with the distribution of an electric field. The paper designs the borehole-ground observation installation by setting 4 radial survey lines taking the well mouth as the center, thus, the surveying direction can be set under the unknown geologic bodies. The calculating of typical models shows that the subdivision method is correct and feasible. The surveying method by arraying around wells to detect the unknown position of geologic bodies is effective.

Key words:borehole-ground; radial subdivision; FEM; resisitivity; polarizability

[文獻標志碼][分類號] P631.3 A