起伏地表?xiàng)l件下的井中激電井地觀測正演模擬研究

王 智,吳愛平,李 剛,2

(1.長江大學(xué)電子信息學(xué)院,湖北荊州434023;2.中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院,湖北武漢430074)

井中激發(fā)極化(井中激電)時(shí)將供電裝置/測量裝置/供電裝置和測量裝置放入井中使其靠近異常體,以提高激發(fā)極化強(qiáng)度或激發(fā)極化響應(yīng)強(qiáng)度。相較井井觀測和地井觀測,井地觀測具有觀測數(shù)據(jù)量大,可精確反演的優(yōu)點(diǎn)。目前主要采用井地電阻率法來圈定油氣藏的邊界[1],張?zhí)靷怺2-3]分別在井中無套管與井中有套管的情況下進(jìn)行物理水槽實(shí)驗(yàn)以確定油氣藏邊界,并證明了利用井地電阻率法確定異常體邊界的可行性;張?zhí)靷惖萚4-6]利用非無窮遠(yuǎn)處三極剖面法確定油氣藏邊界,并在新疆等地成功發(fā)現(xiàn)了小塊油氣藏;湯井田等[7]利用井地電阻率法得到的歧離率確定高電阻率油藏邊界。

激發(fā)極化的正、反演以直流電阻率法正、反演為基礎(chǔ)。直流電阻率正演主要采用有限差分法[8-15]、有限單元法[16-22]和邊界單元法[23-25]。有限差分法受限于曲面邊界的處理而有限單元法的優(yōu)勢在于處理曲面邊界;邊界單元法適用于純地形因素的起伏地表模擬,具有速度快、精度高與節(jié)約內(nèi)存等優(yōu)點(diǎn),但無法實(shí)現(xiàn)復(fù)雜模型的模擬;有限單元法作為一種主要的直流電阻率正演方法,處理復(fù)雜的模型和邊界時(shí)有較強(qiáng)的靈活性和適應(yīng)性,可以自動滿足其內(nèi)部邊界條件。在實(shí)際勘探中,地形的起伏會造成高電阻率或低電阻率干擾,并對反演結(jié)果造成不可預(yù)計(jì)的影響,有限單元法適合模擬復(fù)雜物性和起伏地表模型。直流電阻率正演存在電源奇異性的問題,即源點(diǎn)對應(yīng)節(jié)點(diǎn)處的電位無窮大,引起模擬誤差增大,目前解決該問題的方法主要有異常電位法[14,26]和源點(diǎn)附近局部加密網(wǎng)格[19,20]兩種方法。在源點(diǎn)附近加密網(wǎng)格雖然提高了局部的模擬精度,但增加了計(jì)算量,此外網(wǎng)格結(jié)構(gòu)的任意性也破壞了有限單元法固有的收斂結(jié)構(gòu),在提高局部精度同時(shí)降低了整體的收斂速度[27]。LOWRY等[26]提出了異常去除技術(shù),該技術(shù)將總電位分為一次場電位(含奇異值部分)和二次場電位(非奇異值部分),一次場電位如平坦地形下的均勻半空間、水平層狀介質(zhì)或者垂直接觸面等電位,可用簡單模型的背景場電位方程解析解來替代,將該解析解代入原總電位滿足的方程,可消除方程右端電源項(xiàng)狄拉克函數(shù)引起的奇異性。起伏地表?xiàng)l件下,簡單模型的背景場電位方程無解析解,因此無法直接使用異常電位法。BLOME[25]利用邊界單元法對背景場電位方程單獨(dú)求解;PENZ[28]重新定義了異常電位法,背景場電位仍然為平坦地形下均勻半空間電位方程的解析解,將起伏地表與地下不均勻體電位方程作為異常場電位方程,將平坦地形下均勻半空間電位方程的解析解作為背景場方程的解,采用廣義有限差分法對異常電位的新邊界條件進(jìn)行求解。

目前對井中激發(fā)極化法的應(yīng)用以及正、反演研究非常少[29-31],大多數(shù)的電阻率法與激發(fā)極化法的正、反演研究都是基于平坦地形的。但實(shí)際應(yīng)用中地形的影響不可忽略,因此本文基于PENZ[28]的研究,將改進(jìn)的異常電位法引入起伏地表?xiàng)l件下井中激電的正演模擬。基于二維線源場,采用加權(quán)余量法得到有限單元積分方程,并引入直接解法求解器求解最終的線性方程組。將幾個(gè)典型的起伏地表模型試算結(jié)果與采用傳統(tǒng)的總電位法得到的結(jié)果進(jìn)行比較,驗(yàn)證了本方法的正確性與有效性。

1 井中激電正演模擬算法

1.1 直流電阻率法二維線源邊值問題

在直角坐標(biāo)系中,假設(shè)Z軸垂直向下,源點(diǎn)A的坐標(biāo)為(xA,zA),則在二維線源場中總電位滿足以下方程:

(1)

式中:δ是狄拉克函數(shù);σ為地層模型的電導(dǎo)率;Ω為整個(gè)計(jì)算區(qū)域;in表示在區(qū)域內(nèi)。求解式(1)還需要相應(yīng)的邊界條件,邊界分為地表邊界Γs與無窮遠(yuǎn)邊界Γ∞。因地表無電流流向空氣中,故滿足紐曼邊界條件:

(2)

式中:n為邊界Γs的外法向方向;?！逓闊o窮遠(yuǎn)邊界;on表示在邊界上。COGGON[32]采用了以下兩種邊界條件:

試算結(jié)果表明,采用(3)式計(jì)算得到的電位偏小,采用(4)式計(jì)算得到的電位偏大,因此DEY等[10]引入了以下混合邊界條件:

(5)

式中:n為邊界?！薜耐夥ㄏ蚍较?r為源點(diǎn)到無窮遠(yuǎn)邊界?！薜木嚯x;r為源點(diǎn)到無窮遠(yuǎn)邊界?！薜挠邢蚓嚯x;(r,n)為r,n之間的夾角。(1)式,(2)式和(5)式共同稱為總電位的邊值問題。

1.2 改進(jìn)的異常電位法

(1)式右端項(xiàng)包含的狄拉克函數(shù)會造成源點(diǎn)處的奇異值,從而引起源點(diǎn)附近的模擬誤差增大。目前有兩種解決該問題的方法,一種為加密源點(diǎn)附近的網(wǎng)格[20];另一種為文獻(xiàn)[26]提出的異常消除技術(shù)。該技術(shù)將總電位u分離為一次場電位up(奇異值部分)與二次場電位us(非奇異值部分),一次場電位為平坦地形下的均勻半空間、層狀空間或者垂直接觸面的電位等,可使用簡單模型的背景場電位方程解析解來替代,總電位定義為:

(6)

異常電位的邊值問題如下[14,18,33]:

(7)

式中:σ為介質(zhì)電導(dǎo)率;σ0為背景場電導(dǎo)率[14]。當(dāng)源點(diǎn)在地面時(shí),平坦地形條件下均勻半空間的解析解為up=I/(πσ0)ln(1/r0),其中r0為源點(diǎn)到觀測點(diǎn)的距離。(7)式不包含點(diǎn)源項(xiàng),且在起伏地表?xiàng)l件下背景場電位方程無解析解,通常采用邊界單元法[23-25]來進(jìn)行數(shù)值求解,增加了計(jì)算量。在起伏地表?xiàng)l件下,本文采用PENZ等[28]提出的方法,將總電位分為一次場電位和二次場電位,一次場電位仍然按平坦地形條件下的正常電位計(jì)算up=1/(2πσ0r0),修改自然邊界條件,定義新的自然邊界條件為:

(8)

由(8)式可知,一次場電位滿足非齊次的邊界條件,理論上說明存在流向空氣的電流,將這部分的電位定義為“虛擬電位”,因?yàn)閷?shí)際中不存在流向空氣的電流[34]。我們認(rèn)為二次場電位us是由起伏地表與地下不均勻體產(chǎn)生的異常電位。由(5)式,(6)式和(7)式可得無窮遠(yuǎn)邊界條件為[34]:

(9)

(7)式、(8)式和(9)式組成了起伏地表?xiàng)l件下基于線源的異常電位法的邊值問題,通常采用有限差分法、有限單元法等方法求解該類邊值問題,本文采用有限單元法并使用非結(jié)構(gòu)化網(wǎng)格對地電模型及起伏地表進(jìn)行剖分近似。

1.3 有限單元-加權(quán)余量法

有限單元法求解變分問題,必須先得到對應(yīng)的有限單元積分方程,將滿足邊界條件的微分方程轉(zhuǎn)換為對應(yīng)的積分方程,通常采用兩種方法:里茲(Ritz)法和加權(quán)余量法。里茲法利用變分原理將邊界條件轉(zhuǎn)化為等價(jià)泛函的極值問題;加權(quán)余量法引入權(quán)函數(shù)將微分方程轉(zhuǎn)化為弱形式的積分方程,當(dāng)基函數(shù)為權(quán)函數(shù)時(shí)稱為伽遼金(Galerkin)法[33]。變分問題的推導(dǎo)較復(fù)雜,依據(jù)標(biāo)準(zhǔn)的變分原理只能處理自伴正定,且邊界條件為第一類或第三類齊次的微分算子方程,大量的偏微分方程均無對應(yīng)的泛函。通常采用較為簡單的Galerkin方法求解變分問題,前人的研究表明,如果邊值問題存在相應(yīng)的泛函,采用Galerkin法會得到與采用Ritz法相同弱形式的積分方程。

將區(qū)域Ω剖分成一系列的三角形單元e,節(jié)點(diǎn)i,j,m為任意3個(gè)節(jié)點(diǎn),如圖1所示。

圖1 三角形單元

假設(shè)單元內(nèi)的電位是線性分布的,則單元中的電位u可以表示為:

(10)

式中:ui,uj,um分別為i,j,m這3個(gè)節(jié)點(diǎn)上的電位;Ni,Nj,Nm分別為3個(gè)節(jié)點(diǎn)對應(yīng)的形函數(shù):

其中,中間變量ai,aj,am,Δ分別為:

根據(jù)加權(quán)余量法原理,(7)式的余量為:

(11)

將形函數(shù)Ni作為權(quán)函數(shù),得到三角形單元e的加權(quán)余量積分表示并令其等于0:

(12)

將(11)式代入(12)式可得:

(13)

式中:σ′為異常電導(dǎo)率,σ′=σ-σ0。(13)式可表示為:

(14)

i=1,2,3

根據(jù)散度定理,將(10)式代入(14)式可得:

(15)

i,j=1,2,3

(16)

(16)式可表示為以下矩陣形式:

(17)

(18)

由文獻(xiàn)[33]可得:

(19)

當(dāng)?shù)叵陆橘|(zhì)無極化特性時(shí),正演模擬得到的一次場電位為u1;當(dāng)?shù)叵码娮杪蕿棣训慕橘|(zhì)存在極化特性時(shí),正演模擬得到極化場總電位為u,兩次正演模擬結(jié)束之后,得到視極化率如下(基于二級裝置):

(20)

式中:u2為二次場電位。

2 數(shù)值算例及分析

網(wǎng)格剖分是有限單元法重要的環(huán)節(jié),因?yàn)榉墙Y(jié)構(gòu)化網(wǎng)格可以局部加密網(wǎng)格,因而其在精確模擬各種復(fù)雜模型時(shí)有著天然的優(yōu)勢。對于非結(jié)構(gòu)化網(wǎng)格,我們利用IntelMKL庫提供的Pardiso求解器求解線性方程組。為了檢驗(yàn)改進(jìn)的異常電位法的效果,首先利用二維線源條件下特殊地形存在解析解這一優(yōu)勢,來驗(yàn)證該方法的正確性;然后選擇兩個(gè)理論模型:山谷與山脊模型檢驗(yàn)該方法對起伏地表的計(jì)算效果;最后分析地形對二次場電位的影響規(guī)律,并與總電位法得到的結(jié)果進(jìn)行對比分析。

2.1 算法驗(yàn)證

在二維線源的情況下,經(jīng)保角變換可得到特殊地形下的解析解,保角變換的原理如下:

取變換函數(shù):

(21)

式中:z為Z平面里的任意一點(diǎn)的坐標(biāo),ω為W平面里任意一點(diǎn)的坐標(biāo),d為如圖2a所示的邊界DE長度。

將Z平面(圖2a)的邊界CDEFG及其下部平面變換成W平面的實(shí)軸及其下部平面(圖2b),由(21)式可得:

圖2 Z平面變換到W平面示意a Z平面; b W平面

(22)

將z=x+iy,w=u+iv代入式(22)并展開得:

(23)

式中:x的正負(fù)與u相同;y的正負(fù)與v相同。對于給定的d和w(u,iv),可根據(jù)(23)式計(jì)算w(u,iv)的對應(yīng)點(diǎn)z在Z平面內(nèi)的坐標(biāo)(x,iy);對于W平面內(nèi)的直線R′S′,可根據(jù)(23)式得到Z平面內(nèi)對應(yīng)的曲線RS。W平面內(nèi)(相當(dāng)于平坦地形下)的電位為:

(24)

式中:I為供電電流。

圖3為二維純山谷地形,厚度為20m,電流強(qiáng)度為1A,電導(dǎo)率為1.0S/m,五角星為源位置。表1為改進(jìn)的異常電位法求得的總電位與解析解、總電位法的結(jié)果對比。該異常電位法得到的結(jié)果和解析解的誤差小,同時(shí)提高了源點(diǎn)附近的總電位精度。

圖3 二維純山谷地形示意

2.2 計(jì)算實(shí)例

模型1:圖4,圖5分別為山谷與山脊模型,五角星均為源位置,均勻介質(zhì),各項(xiàng)參數(shù)如下:電阻率為1Ω·m,三極測量裝置中測量電極M,N變化范圍為0～200m,測量電極距2m,山谷模型網(wǎng)格剖分的計(jì)算區(qū)域?yàn)?00m×100m,擴(kuò)展邊界區(qū)域?yàn)?0000m,三角形網(wǎng)格單元數(shù)為21151,節(jié)點(diǎn)數(shù)為10864;山脊模型網(wǎng)格剖分的計(jì)算區(qū)域?yàn)?00m×100m,擴(kuò)展邊界長度為10000m,三角形網(wǎng)格單元數(shù)為21261,節(jié)點(diǎn)數(shù)為10919。

從圖6與圖7可以看出采用總電位法得到的山谷模型視電阻率曲線出現(xiàn)了高阻異常,山脊模型的視電阻率曲線出現(xiàn)了低阻異常,這與前人研究得出的地形與視電阻率呈鏡像對稱的規(guī)律相符[35];在源點(diǎn)附近,采用總電位法得到的視電阻率出現(xiàn)了畸變,這也驗(yàn)證了改進(jìn)的異常電位法的正確性,與總電位法得到的結(jié)果相比,改進(jìn)的異常電位法得到的視電阻率精度更高。

模型2:該山谷模型的各項(xiàng)參數(shù)如圖8所示,圍巖電阻率ρ0=100Ω·m,極化率η0=0.01,低阻異常體電阻率ρ1=5Ω·m,極化率η1=0.2,異常體寬度20m,厚度10m,埋深10m,網(wǎng)格剖分的計(jì)算區(qū)域?yàn)?00m×100m,擴(kuò)展邊界長度為10000m,生成三角形網(wǎng)格單元數(shù)為16169,節(jié)點(diǎn)數(shù)為8253。五角星為源在井中的位置,采用井地二極觀測方式,視電阻率計(jì)算公式為:ρs=πu/ln(1/r)。

表1 采用改進(jìn)的異常電位法求得的總電位與解析解、總電位法求得的結(jié)果對比(源點(diǎn)x=0)

圖4 二維山谷模型

圖5 二維山脊模型

圖6 山谷模型的視電阻率曲線

圖7 山脊模型的視電阻率曲線

圖8 二維山谷模型(含低阻高極化異常體)

采用改進(jìn)的異常電位法對該模型進(jìn)行正演模擬,得到了線源在不同深度時(shí)的視電阻率曲線,為了便于比較起伏地表?xiàng)l件下異常體的視電阻率響應(yīng)特征,還給出了純山谷地形的視電阻率曲線和包含低阻異常體的平坦地形視電阻率曲線(圖9a,圖9b)。從圖9a可以看出,山谷地形造成了高阻異常的畸變,隨著源深度的增加,山谷地形造成的高阻異常逐漸變小;圖9b 為包含低阻異常體的平坦地形(異常體的各項(xiàng)參數(shù)與模型2一致)隨著源深度的增加視電阻率的變化規(guī)律,源在異常體上方時(shí),視電阻率曲線表現(xiàn)出高阻特征,源在異常體頂端時(shí)(源在地下20m處),視電阻率最大;源在異常體下方時(shí),視電阻率曲線表現(xiàn)出低阻特征,當(dāng)源在異常體底部時(shí)(源在地下30m處),視電阻率的幅值達(dá)到最小,隨著源深度增加,視電阻率曲線趨于平緩;從圖9c可以看出,地形因素造成的高阻異常與低阻異常體產(chǎn)生的異常疊加,掩蓋了原有低阻異常體的視電阻率變化規(guī)律,造成了高阻異常的假象,這會產(chǎn)生錯誤的地震解釋成果;當(dāng)源位于地下20m處時(shí),分別采用總電位法和改進(jìn)的異常電位法對模型2進(jìn)行正演模擬,為了便于對比分析,又采用總電位法對純山谷地形與含低阻異常體的平坦地形(異常體的各項(xiàng)參數(shù)與模型2一致)進(jìn)行正演模擬,共獲得了4條視電阻率曲線,如圖9d所示,可以看出起伏地表?xiàng)l件下,總電位法與異常電位法的結(jié)果幾乎一致,再次驗(yàn)證了異常電位法的正確性與有效性,純山谷地形導(dǎo)致的視電阻率曲線異常與含有低阻異常體的平坦地形產(chǎn)生的異常相似,這也意味著起伏地表對視電阻率影響非常大。

圖9 視電阻率曲線a 純山谷地形的視電阻率曲線; b 含低阻異常體的平坦地形視電阻率曲線; c 模型2的視電阻率曲線; d 不同地形條件下采用不同方法得到的視電阻率曲線(源在地下20m處)

從圖10a的視極化率曲線可以看出,源在低阻高極化異常體(異常體的各項(xiàng)參數(shù)與模型2一致)上方時(shí),視極化率低于圍巖的極化率,為低值異常;源在該異常體頂端時(shí),視極化率變化輻值達(dá)到最大;源在異常體下方時(shí),視極化率曲線變化趨勢相反。圖10b為模型2的視極化率曲線,視極化率的變化規(guī)律幾乎和圖10a一致,僅值不同。圖10c與圖10d分別為含低阻高極化異常體的平坦地形與模型2的視極化率斷面圖,從圖10c中可以看到,兩個(gè)異常的邊界幾乎對應(yīng)了低阻高極化異常體的上、下邊界(白色虛線為異常體位置);圖10d中兩個(gè)異常邊界與圖10c中的異常邊界幾乎一樣,這也說明了山谷地形未導(dǎo)致激電假異常,僅對視極化率的值產(chǎn)生影響。

圖10 視極化率曲線a 含低阻高極化異常體的平坦地形視極化率曲線; b 模型2的視極化率曲線; c 含低阻高極化異常體的平坦地形視極化率斷面; d 模型2的視極化率斷面

3 結(jié)論

針對起伏地表?xiàng)l件下傳統(tǒng)的異常電位法中背景場電位方程無解析解這一問題,從二維線源出發(fā),采用有限單元法,對起伏地表?xiàng)l件下的二維地電模型進(jìn)行正演模擬,在前人的研究基礎(chǔ)之上,重新定義異常電位法,將起伏地表與地下不均勻體視作異常場,采用加權(quán)余量法對新的邊界條件進(jìn)行推導(dǎo),得到了對應(yīng)的有限單元積分方程,引入直接法求解器,求解最終的線性方程組得到總電位。將二維線源條件下山谷模型的解析解以及總電位法得到的結(jié)果與改進(jìn)的異常電位法得到的結(jié)果對比,驗(yàn)證了該方法的正確性與有效性。含低阻高極化異常體的山谷模型正演模擬結(jié)果表明,山谷地形改變了原有的低阻高極化異常體變化規(guī)律,造成了高阻假異常,但未導(dǎo)致激電假異常,僅影響視極化率的值。