井地-高密度電法聯(lián)合數(shù)值模擬研究

侯智超， 戴前偉， 尹俊濤,3

(1. 河南省交通規(guī)劃勘察設(shè)計(jì)研究院有限責(zé)任公司， 鄭州 450052;2.中南大學(xué)有色金屬成礦預(yù)測教育部重點(diǎn)實(shí)驗(yàn)室， 長沙 410083; 3.中南大學(xué) 地球科學(xué)與信息物理學(xué)院， 長沙 410083)

0 引言

井地電法是通過在井中供入大功率的電流作為激發(fā)源，在地表布置一定數(shù)量的電極進(jìn)行接收的勘探方式，是以地下介質(zhì)的電阻率、極化率等差異為物理基礎(chǔ)。井地電法在1958年由前蘇聯(lián)人[1]率先提出，并且首次用于圈定煤層水平邊界；Alfano[2]計(jì)算了當(dāng)在地下供電時(shí)，三層介質(zhì)模型的地表電位解析解的分布，與常規(guī)電測深方法相比有更好的優(yōu)勢；Daniels[3]研究了球體和n層層狀介質(zhì)視電阻率的異常響應(yīng)，其結(jié)果表明在井-井，井-地不同裝置下并不存在“最好”的方法，其方法的選擇取決于目標(biāo)體的埋深、井井之間的間距等；在國內(nèi)譚河清等[4]利用井地電位技術(shù)在孤東八區(qū)進(jìn)行了剩余油分布的研究，其結(jié)果表明井地電位測量技術(shù)與其他監(jiān)測的結(jié)果吻合很好，顯示了井地電法技術(shù)在剩余油分布等方面良好的應(yīng)用前景；王志剛等[5]進(jìn)行了室內(nèi)水槽的三維物理模型實(shí)驗(yàn)，研究物理場的分布特征和能有效圈定油藏邊界和不同儲層分布范圍的最佳供電激發(fā)點(diǎn)的選擇依據(jù)；屈有恒等[6]推導(dǎo)了傾斜線狀電流源的理論電位，并且利用有限差分法計(jì)算了在垂直和傾斜線源情況下地表視電阻率分布的變化；湯井田等[7]根據(jù)井地有限元數(shù)值模擬的計(jì)算結(jié)果，提出了采用井地電阻率法歧離率確定高阻油氣藏邊界，討論了不同深度油氣藏和油氣藏偏離鉆井位置時(shí)歧離率的變化規(guī)律，結(jié)果表明歧離率可以很好地確定油氣藏的異常并且對其橫向邊界圈定有很好的效果；戴前偉等[8]基于有限單元法對垂直線源井地電位法中的地下動態(tài)導(dǎo)體進(jìn)行了三維數(shù)值模擬，其結(jié)果表明垂直線源井地電位對探測地下動態(tài)導(dǎo)體可行有效的；陳德鵬[9]利用有限元實(shí)現(xiàn)了任意線源井地電法三維正演模擬并分析了其異常規(guī)律；潘紀(jì)順等[10]提出地面/井地/井間超高密度泛裝置，為電阻率成像新方法的應(yīng)用提供了新的思路和前景；高密度電法以直流電阻率法為基礎(chǔ)，通過向地下發(fā)送電流形成穩(wěn)定的電流場，在地表上布置多個(gè)電極進(jìn)行電位測量，然后通過轉(zhuǎn)換裝置對布設(shè)電極的斷面進(jìn)行自動探測和記錄的方法[11]；在進(jìn)行野外測量時(shí)，井地電法一般是將源放于井中[12]，地表觀測，使其有更大的探測深度，但是對于其橫向的分辨率不夠精確，因此作者采取井地電法和高密度電法進(jìn)行聯(lián)合正演模擬，將高密度電法和井地電法聯(lián)合起來對地下目標(biāo)體進(jìn)行探測，能夠使我們得到更多的地下場的信息，有助于在生產(chǎn)中實(shí)際問題的解決，具有很好的生產(chǎn)實(shí)踐意義。

1 井地-高密度電法觀測方式

井地-高密度電法的測量分為兩種：①將供電電極放入井下慢慢移動，在地表利用測量電極測量地表電位分布情況，然后利用地表電位求取視電阻率(A為供電電極，M為測量電極，如圖1所示)；②將電流源放在地表，其原理與高密度電法的測量方式相同(A為供電電極，M為測量電極 ，如圖2所示)。

圖1 井地電法測量方式示意圖Fig．1 The borehole-ground electrical method measurement diagram

圖2 高密度電法測量方式示意圖Fig．2 The high-density resistivity method measurement diagram

兩種測量方式主要是因?yàn)樵次恢玫牟煌鼈兊臏y量電極都是在地表進(jìn)行電位測量，因此在實(shí)際生產(chǎn)實(shí)踐中可以將兩種測量方式合二為一，通過源在地表和井中移動，分別采集數(shù)據(jù)。我們將采集的數(shù)據(jù)分別進(jìn)行處理，當(dāng)源在地表時(shí)，可以將采集的數(shù)據(jù)按照高密度電法數(shù)據(jù)進(jìn)行處理；而當(dāng)源在井中時(shí)，將其采集的數(shù)據(jù)作為井地電法數(shù)據(jù)進(jìn)行處理。

2 2.5維井地-高密度電法的變分問題

井地電法和高密度電法都是點(diǎn)源的變分問題，其不同之處僅為點(diǎn)源的位置和最后記錄點(diǎn)的選取，所以利用徐世浙[13]所述的方法，點(diǎn)源三維變分問題對z方向進(jìn)行傅里葉變換可得2.5維邊值問題

(1)

式(1)與下列變分問題等價(jià)

(2)

其中U是傅式變換的總電位；k為波數(shù)；σ為介質(zhì)的電導(dǎo)率；I為電流；Γs是地面邊界；?！奘菬o窮遠(yuǎn)邊界;Ω為Γs和?！藿M成的封閉區(qū)域；n為電流密度的法線方向；r是通過坐標(biāo)原點(diǎn)、與走向垂直的截面內(nèi)，Γ∞邊界上的點(diǎn)至坐標(biāo)遠(yuǎn)點(diǎn)的距離；K0是第二類零階修正貝塞爾函數(shù)；K1是第二類一階修正貝塞爾函數(shù)；cos(r,n)是矢徑r與外法線方向n的夾角余弦。

3 有限單元法

有限單元法的求解步驟：

(1)網(wǎng)格剖分。首先將求解的二維區(qū)域剖分成矩形單元，如圖1所示

圖3 網(wǎng)格剖分及節(jié)點(diǎn)編號示意圖Fig．3 Sketch map of mesh division and node number

(2)線性插值。圖4(a)是母單元示意圖，圖4(b)是子單元示意圖，兩個(gè)單元間的坐標(biāo)變換關(guān)系為：

x=x0+(a/2)·ξ，z=z0+(b/2)·η

(3)

其中x0、z0是子單元中點(diǎn)的坐標(biāo);a、b是子單元的兩個(gè)邊長，微分關(guān)系為：

dx=(a/2)·dξ，dz=(b/2)·dη，dxdy=(ab/4)·dξdη

(4)

單元形函數(shù)為雙線性插值形函數(shù)，可寫為

(5)

其中ξi、ηi是點(diǎn)i(i=1,2,3,4)的坐標(biāo)，形函數(shù)的分量可寫為

N1=(1-ξ)(1-η)/4

N2=(1+ξ)(1-η)/4

N3=(1+ξ)(1+η)/4

N4=(1-ξ)(1+η)/4

(6)

圖4 母單元、子單元示意圖Fig．4 Sketch map of parent element and sub-element(a)母單元；(b)子單元

(3) 單元分析及總體合成。式(2)中第一項(xiàng)積分為

(7)

其中K1e=(k1ij)；k1ij=k1ji；ue=(ui)T；i、,j=1、2、3、4。

式(2)中第二項(xiàng)積分為

(8)

式(2)中第四項(xiàng)積分為

(9)

根據(jù)δ函數(shù)的積分性質(zhì)，式(2)中第三項(xiàng)積分為

(10)

式(10)說明只有點(diǎn)電源的UA對該項(xiàng)積分才有貢獻(xiàn)。

將式(7)、式(8)、式(9)、式(10)帶入到式(2)中得

(11)

將單元向量Ue擴(kuò)展成ND×1的向量U，Ke擴(kuò)展成ND×ND的矩陣K得：

KU=P

(12)

其中Ke=K1e+K2e+K3e；

(4) 方程組求解和視電阻率的計(jì)算。利用LDLT法求解線性方程組式(12)可得各節(jié)點(diǎn)傅氏電位U，然后通過反傅氏變換即可得到三維空間中的電位u。作者采用井地電法和高密度電法均為二極裝置[15]進(jìn)行觀測，其視電阻率ρS計(jì)算公式為

(13)

其中K為裝置系數(shù)；Δu為電位差。

4 模型算例及分析

在上述理論的基礎(chǔ)上，編制了相應(yīng)的程序，為了驗(yàn)證程序的正確性和有效性，首先對一地下均勻半空間模型進(jìn)行正演模擬計(jì)算，并與解析解進(jìn)行了對比分析，然后對一系列的模型進(jìn)行了正演模擬計(jì)算。

4.1 模型1 地下均勻半空間模型

為了驗(yàn)證正演算法的可行性，作者選取均勻半空間進(jìn)行計(jì)算。地下介質(zhì)的電阻率為100 Ω·m，電流大小為1 A，計(jì)算區(qū)域?yàn)?00 m×70 m，采用雙線性插值，網(wǎng)格剖分為176×40。利用高密度正演模擬程序求出各節(jié)點(diǎn)的電位u，然后通過將電位代入視電阻率公式，求出視電阻率與地下真實(shí)電阻率進(jìn)行對比。

圖5為高密度電法得到的地下視電阻率的二維斷面圖，將視電阻率與地下真實(shí)電阻率進(jìn)行對比分析， 經(jīng)過計(jì)算得出的視電阻率最大值為 99.79 Ω·m，最小值為95.59 Ω·m，而且最小值基本存在于邊緣區(qū)域，其最小相差為0.21 Ω·m，其誤差為0.21%，最大相差為4.41 Ω·m，其誤差為4.41%，確定了正演算法的有效性。

圖6為井地電法得到的視電阻率曲線圖，將視電阻率與地下真實(shí)電阻率進(jìn)行對比分析，經(jīng)過計(jì)算得出的視電阻率最大值為99.91 Ω·m，最小值為99.01 Ω·m，而且最小值基本存在于點(diǎn)源點(diǎn)附近，其最小相差為0.09 Ω·m，其誤差為0.09%，最大相差為0.99 Ω·m，其誤差為0.99%，確定了正演算法的正確性。

4.2 模型2高阻異常體模型

模型2如圖7所示為單個(gè)對稱的高阻體示意圖，背景介質(zhì)的電阻率大小為100 Ω·m，高阻體電阻率為1000 Ω·m，其頂部埋深為4 m，長為16 m，寬為8 m，與井中心對稱。計(jì)算區(qū)域的大小為100 m×40 m，點(diǎn)電流源的供電大小為1 A。

首先利用井地電法，對上述模型進(jìn)行有限元正演模擬，得到圖8井地電法正演模擬的視電阻率擬斷面圖和圖9井地電法正演模擬的視電阻率曲線圖。從圖9中可以看出，在高阻體的上方，視電阻率曲線表現(xiàn)出相對高阻的異常響應(yīng)特征，當(dāng)隨著點(diǎn)電流源的向下移動，視電阻率響應(yīng)逐漸增大，且幅度逐漸增大；當(dāng)達(dá)到高阻體的上表面時(shí)，即當(dāng)點(diǎn)電源在h=-4 m，其響應(yīng)視電阻率達(dá)到最大，其異常極大值為ρs=186.56 Ω·m，隨著點(diǎn)電源繼續(xù)向下移動，即點(diǎn)電源在高阻體內(nèi)部時(shí)，其視電阻率曲線逐漸減小，但仍然呈現(xiàn)相對高阻的響應(yīng)特性，當(dāng)點(diǎn)電源向下移動至高阻體的下表面，即h=-12 m，其響應(yīng)視電阻率達(dá)到最小，其異常極小值為ρs=68.01 Ω·m，之后隨著點(diǎn)電源深度的繼續(xù)增加，其響應(yīng)視電阻率曲線逐漸增大，慢慢趨向于背景介質(zhì)的電阻率。從圖8中可以看出，高阻體的上下表面正好對應(yīng)兩個(gè)異常中心深度的的位置。從圖9中還可以看到，在高阻體的上方，視電阻率曲線圖的曲線中間有一個(gè)拐點(diǎn)，其原因?yàn)楫?dāng)點(diǎn)電源離地表較近時(shí)，其地表測量電位受點(diǎn)電源的影響，而隨著點(diǎn)電源深度的增大，其影響很小。

圖5 均勻半空間高密度正演模擬斷面圖Fig．5 The forward modeling of high density resistivity section in the homogeneous half space

圖6 均勻半空間井地電法正演模擬斷面圖Fig．6 The forward modeling of borehole-ground electrical section in the homogeneous half space

圖7 高阻異常體示意圖Fig．7 The high resistivity abnormality model

圖8 模型2點(diǎn)電源不同深度視電阻率曲線分布圖Fig．8 The curves of apparent resistivity with point current source in different depths of model 2

圖9 模型2視電阻率擬斷面圖Fig．9 The apparent resistivity pseudosection map of model 2

圖10為采用高密度電法對模型2計(jì)算得到的視電阻率斷面圖。從高密度電法正演模擬視電阻率斷面圖可以看出，在圖10中出現(xiàn)了一個(gè)基本對稱的相對高阻異常，其異常最大值為ρs=160.40Ω·m，在逐漸遠(yuǎn)離高阻體的位置，其響應(yīng)視電阻率逐漸趨向于其背景介質(zhì)的電阻率。在圖10中可以看出，模型頂部埋深為-4 m，但是從高密度正演模擬的視電阻率斷面圖中可以看出，其縱向位置，異常的位置頂部約在z= -2 m處，在垂直方向上移了約2 m，這主要是受異常高阻體的體積效應(yīng)的影響；而在橫向可以看出，高阻異常閉合圈是關(guān)于x= 0 m對稱的，其影響區(qū)域在x= -10m和x= 10m之間，與模型的位置一致。

圖10 模型2的高密度視電阻率斷面圖Fig．10 The high density apparent resistivity pseudosection map of model 2

從圖8和圖10中可以看出，井地電法和高密度電法均可以有效地對地下高阻異常體進(jìn)行正演模擬。井地電法正演模擬中存在兩個(gè)關(guān)于井對稱的異常，上面是一個(gè)相對高阻異常，下面是一個(gè)相對低阻異常，而高阻體的上下表面正好對應(yīng)兩個(gè)異常中心深度的位置；高密度電法正演模擬中出現(xiàn)一個(gè)明顯關(guān)于井對稱的相對高阻異常，但是高阻異常的縱向位置有一定的偏差。

4.3 模型3邊緣隱伏礦為高阻

圖11為邊緣高阻的模型示意圖，背景介質(zhì)的電阻率為100 Ω·m，高阻體的電阻率大小為 1000 Ω·m，高阻體頂部埋深為11 m，其長為10 m，寬為10 m，距井左側(cè)10 m，計(jì)算區(qū)域大小約為100 m×50 m，點(diǎn)電流源供電大小為1A。

圖11 邊緣高阻的模型示意圖Fig．11 The edge high resistivity abnormality model

圖12為井地電法正演模擬的邊緣高阻的擬斷面圖，從圖12中可以看出，在距井左側(cè)10 m左右有一個(gè)高阻閉合圈，而且閉合圈的中心在x= -15 m且z= -10 m的地方,而且在高阻閉合圈的下方，可以看出有一個(gè)類似低阻異常的閉合圈，低阻異常中心在x= -30 m且z= -22 m地方。這個(gè)從側(cè)面說明，利用井地電法模擬時(shí)，在相對高阻異常下面會出現(xiàn)一個(gè)相對低阻異常的正確性。依照前面的異常規(guī)律響應(yīng)，兩個(gè)異常中心的z方向近似對應(yīng)異常體的上下表面，從圖12中也可以得到近似的響應(yīng)，但是這時(shí)候并不能確定其異常體具體范圍。因此需要結(jié)合源在地面時(shí)高密度正演模擬得到的圖進(jìn)一步確定異常體的范圍。

圖13為采用高密度電法對模型3計(jì)算得到的視電阻率斷面圖，從上圖13可以看出，在距離井左側(cè)10 m的位置，可以看出有一個(gè)高阻異常閉合圈，在縱向上從圖13中可以看出，高阻異常閉合圈的頂部是在z= -5 m左右，而高阻模型頂部埋深在z= -10 m，其高阻異常在垂向方向上移了約5 m左右；在橫向可以看出，這個(gè)高阻異常閉合圈是關(guān)于x= -15 m左右對稱的，其影響區(qū)域約在x= -10 m和x= -20 m之間，這與模型的橫向位置是一致的。

從圖12和圖13得出：通過井地電法可以圈定邊緣高阻體的縱向上下表面的深度位置在z= -10 m到z= -22 m，從高密度電法可以近似圈定高阻體橫向的范圍約在x= -10 m和x= -20 m之間且關(guān)于x= -15 m對稱，因此結(jié)合兩種方法，可以近似確定邊緣高阻異常體范圍。

圖12 模型3點(diǎn)電源不同深度視電阻率曲線分布圖Fig．12 The curves of apparent resistivity with point current source in different depths of model 3

圖13 模型3的高密度視電阻率斷面圖Fig．13 The high density apparent resistivity pseudosection map of model 3

5 結(jié)論

作者采用有限元法實(shí)現(xiàn)了井地-高密度2.5維視電阻率的正演模擬，給出了井地-高密度電法2.5維有限元正演模擬的詳細(xì)解法；在井地電法正演模擬的基礎(chǔ)上，對目標(biāo)體進(jìn)行高密度電法的正演計(jì)算，從而結(jié)合兩者對目標(biāo)體進(jìn)行正演模擬，對地下目標(biāo)體的范圍可以更好地確定。

參考文獻(xiàn)：

[1] 黃仲良.石油重·磁·電法勘探[M]. 東營:石油大學(xué)出版社,1999．

[2] ALFANO L. Geoelectric prospecting with underground electrodes[J]. Geophysics Prospecting,1961,10:290-303.

[3] DANIELS J J. Three-dimensional resistivity and induced polarization modeling using buried electrodes[J]. Geophysics,1977,42:1006-1019.

[4] 譚河清,沈金松,周超,等. 井地電位成像技術(shù)及其在孤東八區(qū)剩余油分布研究中的應(yīng)用[J]. 石油大學(xué)學(xué)報(bào)：自然科學(xué)版, 2004,28(2):31-37.

[5] 王志剛,何展翔,魏文博,等. 井地電法三維物理模型試驗(yàn)[J]. 石油地球物理勘探, 2005,40(05):594-597.

[6] 屈有恒,張貴賓,晉風(fēng)明. 傾斜線源的三維電場數(shù)值模擬研究[J]. 物探化探計(jì)算技術(shù), 2007,29(5):431-435,370.

[7] 湯井田, 張繼鋒, 馮兵. 井地電阻率法歧離率確定高阻油氣藏邊界[J].地球物理學(xué)報(bào), 2007, (03):926-931.

[8] 戴前偉,陳德鵬,熊健奎,等. 垂直線源井-地電位探測地下動態(tài)導(dǎo)體的三維有限元數(shù)值模擬[J]. 工程地球物理學(xué)報(bào), 2008,5(6):643-647.

[9] 陳德鵬. 任意線源井地電法三維有限元正演研究[D]. 長沙:中南大學(xué),2009.

[10] 潘紀(jì)順,葛為中,折京平. 地面/井地/井間超高密度電阻率成像技術(shù)[J]. 華北水利水電學(xué)院報(bào), 2010,31(2):74-78.

[11] 王鵬飛.超高密度激電數(shù)據(jù)采集與正反演解釋方法研究[D]. 長沙:中南大學(xué),2012.

[12] 李天成. 電阻率成像技術(shù)的二維三維正反演研究[D]. 北京: 中國地質(zhì)大學(xué), 2008.

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