基于COMSOL的井地電阻率正演研究

馮碩 ，劉得軍 ，張穎穎 ，秦民君 ，兀鳳娜

（1.中國(guó)石油大學(xué)（北京）地球物理與信息工程學(xué)院，北京 102249；2.中國(guó)石油集團(tuán)測(cè)井有限公司生產(chǎn)測(cè)井中心，陜西 西安 710201；3.中國(guó)石油長(zhǎng)慶油田分公司第一采油廠，陜西 延安 716000）

0 引言

目前我國(guó)大部分油田已進(jìn)入高含水開發(fā)期，隨著開發(fā)的不斷深入，為實(shí)現(xiàn)開發(fā)后期油田的穩(wěn)產(chǎn)挖潛，對(duì)剩余油的開發(fā)與研究提出了更高的要求［1-15］。在諸多地球物理勘探方法中，井地電阻率測(cè)量方法具有對(duì)地層影響小、測(cè)量成本低等優(yōu)勢(shì)，在剩余油分布研究中得到了廣泛的應(yīng)用。

井地電阻率測(cè)量方法，是利用鉆孔套管、井液作為線電極供電，通過井口周圍安置的陣列接收電極測(cè)量電位差，然后進(jìn)行數(shù)據(jù)處理，同時(shí)綜合已知地質(zhì)數(shù)據(jù)分析結(jié)果，得到地下局部電阻率的勘探方法［16］。其主要應(yīng)用于圈定油氣藏邊界、油田開發(fā)注水檢測(cè)以及壓裂裂縫成像等方面［17］。為了實(shí)現(xiàn)井地電阻率方法的精確定量解釋，必然需要進(jìn)行快速高效的正演算法研究。

井地電阻率法的正演仿真，主要是通過已知空間的電阻率分布，利用數(shù)學(xué)物理方法計(jì)算得到地表的電位分布。目前解決井地電阻率測(cè)量的正演方法主要有α中心法、有限差分法、有限元法、積分方程法等。α中心法［3］簡(jiǎn)單快捷但精度不高，現(xiàn)在已經(jīng)很少使用。有限差分法［2］的數(shù)學(xué)表達(dá)式簡(jiǎn)單直觀，發(fā)展較早且應(yīng)用較廣，但邊界刻畫能力弱，只使用矩形網(wǎng)格，無(wú)法對(duì)復(fù)雜地質(zhì)構(gòu)造進(jìn)行準(zhǔn)確地模擬。有限元法［4］將求解區(qū)域劃分為離散網(wǎng)格，將每個(gè)網(wǎng)格看作獨(dú)立的單元進(jìn)行求解，這樣可以得到若干個(gè)有限單元的解，這些解的集合構(gòu)成整體函數(shù)的解。積分方程法［5］對(duì)于復(fù)雜地電模型中差異邊界面無(wú)法進(jìn)行有效離散化，尤其是模型中存在多個(gè)非閉合異常體時(shí)，因此該方法現(xiàn)在只是處于試驗(yàn)階段。

水永秀樹等［4］利用有限差分法實(shí)現(xiàn)了垂直線元三維井地電位正演，隨后Tetsuo Aono等［7］討論了垂直線電流三維成像理論。徐凱軍、劉地淵等［8-9］對(duì)垂直線元、任意線元三維地電模型進(jìn)行了正演研究。屈有恒等［10］實(shí)現(xiàn)了傾斜線元的三維電場(chǎng)正演。劉海飛等［11］實(shí)現(xiàn)了連續(xù)介質(zhì)三維線元井地電位的正演計(jì)算。何展翔、王志剛等［12］研究了利用積分方程法對(duì)井地電阻法的三維正演模擬。

本文首先利用COMSOL軟件進(jìn)行建模，然后利用COMSOL可靈活處理各類邊界條件特點(diǎn)，施加混合邊界條件，最后用有限元方法，計(jì)算高、低阻異常體，分析模型變化對(duì)地表電位分布的影響。

1 理論和方法

1.1 線電源場(chǎng)基本分析方法

在井地電阻率法中，電流主要通過鉆孔套管和井液傳到地下，在套管直徑遠(yuǎn)小于套管長(zhǎng)度的情況下，可以將套管簡(jiǎn)化為有限長(zhǎng)線電流源來(lái)處理［13］。

井地電阻率分布規(guī)律滿足以下拉普拉斯形式的微分方程：

式中：σ為線電源所處地層模型電阻率，Ω·m；Φ為模型中各節(jié)點(diǎn)電位，V；f為電流源，A。

線電源的電流強(qiáng)度為I時(shí)，f為

式中：δ為狄拉克函數(shù)，在源點(diǎn)處為無(wú)窮大，其余各點(diǎn)為0；r，rs分別為觀測(cè)點(diǎn)與線電源的位置半徑，m。

對(duì)于均勻半空間介質(zhì)，線電流源產(chǎn)生的電位可以由點(diǎn)電源激發(fā)的電位積分得到：

式中：l2，l1分別為線電流源的頂、底埋深，m；zp為測(cè)點(diǎn)埋深，m；rp為測(cè)點(diǎn)到線電源的徑向距離，m。

地層模型中，電阻率為非均勻的情況下，利用解析法無(wú)法求出模型中的電位分布，這時(shí)需要使用數(shù)值算法進(jìn)行求解。利用有限元進(jìn)行數(shù)值計(jì)算時(shí)，需要確定邊界條件進(jìn)而求解。

在計(jì)算區(qū)域Ω中，井地電阻率法的邊界條件主要包括2種：

1）地面邊界Γs，由于認(rèn)為空氣是絕緣體，選擇使用諾依曼邊界條件 σ·?Ф/?n=0。

2）無(wú)限遠(yuǎn)處邊界面?！?，為保持電位的物理特性，消除電性突變，選擇使用混合邊界條件?Ф/?n-λФ=0，Ф為均勻介質(zhì)情況下利用解析公式得到的電位值。

井地電阻率模型中電位滿足：

1.2 COMSOL正演準(zhǔn)確性驗(yàn)證

COMSOL是一款以MATLAB為內(nèi)核的仿真軟件［14］，相比大多數(shù)仿真軟件，利用COMSOL進(jìn)行電法仿真的主要優(yōu)勢(shì)在于能夠自由定義各類邊界條件，得到更為準(zhǔn)確的正演結(jié)果；正確利用COMSOL中的physics選項(xiàng)可以有效簡(jiǎn)化對(duì)邊界條件的處理。

地表邊界面Γs，可通過選擇COMSOL中的Electric insulation邊界進(jìn)行設(shè)定。對(duì)無(wú)限遠(yuǎn)處邊界面?！?，需要選擇 Distributed resistance 邊界條件 n·J=σ（V-Vref）/d,并將混合邊界條件變形為n·J=-λσФ，然后把Vref選項(xiàng)賦值0，d選項(xiàng)賦值-1/λ，即完成混合邊界條件的處理。

采用驗(yàn)證模型對(duì)該方法的精確性進(jìn)行驗(yàn)證。驗(yàn)證模型空間為1 000 m×1 000 m×1 000 m的正方體，電阻率為10 Ω·m，線電流長(zhǎng)度為100 m，供電電流為1 A，線電源通過地表供電，供電點(diǎn)坐標(biāo)為（0，0，0）。分別對(duì)無(wú)限遠(yuǎn)邊界施加第1類、第2類及混合邊界條件，經(jīng)過正演計(jì)算，得到 5，10，15，20，30，40 m 處的電位值（見表1）。由表1可知，在觀測(cè)范圍區(qū)域內(nèi)，最大絕對(duì)差值為0.001 1 V，最大相對(duì)誤差小于0.25%，這樣的計(jì)算精度是令人滿意的［15］。

表1 不同邊界條件電位值計(jì)算結(jié)果

2 算例分析

正演計(jì)算井地電阻率模型為一個(gè)正方體，為1 000 m×1 000 m×1 000 m，背景電阻率為 10 Ω·m。 異常模型的高阻體電阻率為100 Ω·m，大小為200 m×200 m×100 m，中心點(diǎn)位置在（110，110，-80）。 線電流長(zhǎng)度為200 m，供電電流為5 A，線電源通過地表供電，供電點(diǎn)坐標(biāo)為（0，0，0）。

通過COMSOL計(jì)算，得到模型中每個(gè)有限元節(jié)點(diǎn)位置、地表面上的電位值（見圖1）。

圖1 高阻異常體模型電位分布

由圖1可以看出，在含有高阻異常體的情況下（異常體位置如虛線所示），0.062 63 V等勢(shì)線出現(xiàn)明顯等位圓壓縮現(xiàn)象，并在壓縮曲線上有一個(gè)小突起，其方向指向高阻異常體方向。電位分布情況與井地模型中高阻物質(zhì)對(duì)電流排斥作用相符合。

對(duì)模型中異常體的電阻率進(jìn)行修改，改為1 Ω·m，電阻率與背景電阻率的比例保持不變，電位分布情況見圖2。

圖2 低阻異常體模型電位分布

由圖2可以看出，在有低阻異常體的情況下，0.087 25 V等勢(shì)線明顯向低阻體位置拉伸。電位分布情況與井地模型中低阻物質(zhì)對(duì)電流吸引作用相符。

模型中2個(gè)異常體，電阻率分別為1，100 Ω·m，大小均為200 m×200 m×100 m，其中心點(diǎn)位置分別在（-110，-110，-80）、（110，110，-80），電位分布見圖 3。

圖3 異常體水平排列時(shí)模型電位分布

由圖3可以看出，當(dāng)高、低阻異常體位于同一水平面時(shí)，0.074 31 V等勢(shì)線在高阻異常電阻區(qū)域出現(xiàn)壓縮現(xiàn)象，并依然保留一個(gè)沿高阻異常體對(duì)角線方向的小突起，在低阻異常電阻區(qū)域明顯被拉大。仿真結(jié)果表明，井地電阻率法能夠大致分辨出同一高度面上，不同電阻率異常體的位置。

模型中設(shè)置2個(gè)不同高度面異常體，電阻率分別為 1，100 Ω·m，大小均為 200 m×200 m×100 m。將 2個(gè)異常體上下排列，異常體中心點(diǎn)位置分別放置在（-110，-110，-180）、（-110，-110，-80），電位分布見圖 4。

由圖4可以看出，模型中同時(shí)存在高阻異常體與低阻異常體，且當(dāng)?shù)妥梵w置于高阻體上方時(shí)，0.06300 V等勢(shì)線表現(xiàn)出低阻異常體拉伸現(xiàn)象；與此同時(shí)，在異常體存在方向上，等勢(shì)線也存在一個(gè)小突起，說明電位分布圖依然保留有高阻體存在時(shí)的特征。仿真結(jié)果表明，當(dāng)高阻異常體在下，低阻異常體在上時(shí)，根據(jù)地表電位能夠大致判斷地下電阻分布情況。

圖4 異常體上下排列時(shí)模型電位分布

將異常體調(diào)換位置，當(dāng)高阻異常體置于低阻異常體上方時(shí)，電位分布見圖5。

圖5 高阻異常體置于低阻異常體上方模型電位分布

由圖5可知，高、低電阻體安置位置的互換，地表電位分布已經(jīng)幾乎無(wú)法判斷是否有異常體的存在，等勢(shì)線只存在微弱的不對(duì)稱現(xiàn)象，無(wú)法分析地下異常體電阻率分布情況。仿真結(jié)果表明，對(duì)井地電阻率法的研究應(yīng)當(dāng)考慮上、下不同電阻率分布的地層的影響，在進(jìn)行反演研究時(shí)也需要參照其他前期地質(zhì)資料，對(duì)主要反演區(qū)域的上、下地層進(jìn)行考慮。

3 結(jié)論

1）仿真結(jié)果表明，井地電阻率法對(duì)低阻物質(zhì)較敏感，并不適用于石油領(lǐng)域勘探初期。注水油田進(jìn)入高含水開發(fā)期，地層將出現(xiàn)明顯的電阻率差異，井地電阻率法可以進(jìn)行剩余油氣分布測(cè)量、注水前沿檢測(cè)等項(xiàng)目。

2）實(shí)驗(yàn)表明，異常體縱向分布情況不同，就會(huì)得到不同的地表電位分布，在進(jìn)行反演研究時(shí)也要參照前期地質(zhì)資料。

3）用具體仿真實(shí)例進(jìn)行分析，驗(yàn)證了在使用混合邊界條件的情況下，運(yùn)用COMSOL軟件能夠得到更加精確的正演結(jié)果。

［1］付德奎，馮振雨，曲金明，等.剩余油分布現(xiàn)狀及展望［J］.斷塊油氣田，2007，14（2）：39-41.

［2］楊承林，于春娥，石淑燦.應(yīng)用數(shù)值模擬法研究復(fù)雜斷塊剩余油分布［J］.斷塊油氣田，2007，14（2）：36-38.

［3］Petrick Wm R，Sill Wm R，Wards S H.Three-dimensional resistivity inversionusingalphacenters［J］.Geophysics，1981，46（8）：1148-1162.

［4］水永秀樹，青野哲雄，牛島惠輔.流電電位法（充電法）三維模擬［J］.金旭，譯.物理探查，1991，44（4）：215-226.

［5］Yutaka Sasaki.3-D resistivity inversion using the finite-element method［J］.Geophysics，1995，59（1）：1839-1848.

［6］底青云，王妙月.積分法三維電阻成像［J］.地球物理學(xué)報(bào)，2001，44（6）：843-851.

［7］Aono T，Mizunaga H，Ushijima K.Direct imaging fractures by a 4-D geoelectrical technique at Sumikawa Geothermal Field［A］.Japan：Proc.of 27th Sanford Geothermal Workshop，2002.

［8］徐凱軍，李桐林.垂直有限線源三維地電場(chǎng)有限差分正演研究［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2006，36（1）：137-141.

［9］劉地淵，徐凱軍，趙廣茂，等.任意形狀線電流源三維地電場(chǎng)研究［J］.地球物理學(xué)進(jìn)展，2006，21（2）：395-399.

［10］屈有恒，張貴賓，晉風(fēng)明.傾斜線源的三維電場(chǎng)數(shù)值模擬研究［J］.物探化探計(jì)算技術(shù)，2007，29（5）：431-435.

［11］劉海飛，陳德鵬，戴前偉，等.連續(xù)電性介質(zhì)線源井-地電位三維有限元數(shù)值模擬［J］.桂林理工大學(xué)學(xué)報(bào)，2011，31（1）：28-37.

［12］王志剛，何展翔，劉昱.井地直流電法三維數(shù)值模擬及異常體規(guī)則研究［J］.工程地球物理學(xué)報(bào)，2006，3（2）：87-92.

［13］Kaufman A A，Wightman W E.A transmission-line model for electrical logging through casing［J］.Geophysics，1993，58（12）：1739-1747.

［14］王小龍，馮宏，田華光，等.基于COMSOL MULRIPHYSICS的直流電法正演模擬［J］.煤田地質(zhì)與勘探，2011，39（5）：76-80.

［15］湯井田，肖曉，杜華坤，等.ANSYS在直流電法正演中的應(yīng)用［J］.地球物理學(xué)進(jìn)展，2006，21（3）：987-992.

［16］高軍強(qiáng)，李書忠，李油建.激電法找油的原理及應(yīng)用［J］.斷塊油氣田，2003，10（1）：29-30.

［17］陳亮，牛艷平，張一偉.高含水期剩余油分布的多學(xué)科研究［J］.斷塊油氣田，1999，6（5）：29-33.