王 智,吳愛平,李 剛,2

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

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

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

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

1 井中激電正演模擬算法

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

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

(1)

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

(2)

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

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

(5)

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

1.2 改進的異常電位法

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

(6)

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

(7)

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

(8)

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

(9)

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

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

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

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

圖1 三角形單元

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

(10)

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

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

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

(11)

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

(12)

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

(13)

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

(14)

i=1,2,3

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

(15)

i,j=1,2,3

(16)

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

(17)

(18)

由文獻[33]可得:

(19)

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

(20)

式中:u2為二次場電位。

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

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

2.1 算法驗證

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

取變換函數(shù):

(21)

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

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

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

(22)

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

(23)

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

(24)

式中:I為供電電流。

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

圖3 二維純山谷地形示意

2.2 計算實例

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

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

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

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

圖4 二維山谷模型

圖5 二維山脊模型

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

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

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

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

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

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

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

3 結(jié)論

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