王欣，陳浩，王秀明，曲敏

(中國科學(xué)院聲學(xué)研究所，北京100190)

1 引言

油氣儲(chǔ)層評價(jià)中的電法測井，主要是通過對井下電法測井儀器測得的視電阻率資料進(jìn)行處理，反演出地層的真電阻率或者接近于真實(shí)的巖石電阻率，配合聲波測井和核物理測井，可以最終評價(jià)儲(chǔ)層中流體的性質(zhì)和定量評價(jià)孔隙中含油飽和度，它是油氣儲(chǔ)層測井評價(jià)中重要的手段。側(cè)向測井是電法測井中常用的一種。1934年，Schlumberger兄弟首次以“監(jiān)督電極”的方式提出了側(cè)向測井的思想［1］。隨后，Doll在1951年提出并完善了側(cè)向測井的概念，發(fā)展了三側(cè)向和七側(cè)向［2］?；谇叭斯ぷ鳎琒uau等人在1972年提出并發(fā)展了雙側(cè)向(九側(cè)向)測井技術(shù)［3］，如圖1所示。為了準(zhǔn)確得到地層的真電阻率，就需要高精度、高靈敏度的電法測井儀器，而要設(shè)計(jì)出這樣的電法測井儀器，就要了解雙側(cè)向測井儀器對實(shí)際測井的視電阻率響應(yīng)特征。由于該儀器結(jié)構(gòu)較為復(fù)雜，即便對最簡單的地層模型，其響應(yīng)也沒有解析解。

圖1 雙側(cè)向測井模型

因此只能采用數(shù)值的方法來模擬實(shí)際的儀器響應(yīng)，進(jìn)而優(yōu)化設(shè)計(jì)相關(guān)電法測井儀器的結(jié)構(gòu)，如儀器的長短、電極系的幾何尺寸等等。地層電參數(shù)差異較大、側(cè)向測井儀器自身結(jié)構(gòu)復(fù)雜的特點(diǎn)，使擅長復(fù)雜模型數(shù)值計(jì)算的有限元方法從一開始就成為人們研究側(cè)向測井問題的最主要手段。但由于有限元在自由剖分網(wǎng)格的同時(shí)，會(huì)形成一個(gè)龐大的稀疏矩陣，計(jì)算量非常巨大。為了解決大規(guī)模矩陣運(yùn)算問題，1970年，Irons提出了一種用于有限元分析的波前法(Frontal Solver)［4］。1971 年，Coggon 首次將有限元方法引入電法測井中［5］。1975年，李大潛等人對采用有限元素法進(jìn)行視電阻率測井?dāng)?shù)值計(jì)算的問題進(jìn)行了比較系統(tǒng)的考察［6］。1982年，Anderson用有限元法模擬計(jì)算了感應(yīng)測井響應(yīng)［7］。1986年，張庚驥首次將“波前法”引入電測井問題的有限元模擬中，大大提高了有限元素法在電阻率測井響應(yīng)數(shù)值模擬中的計(jì)算效率［8］。1989年，Khokhar利用有限元法對側(cè)向測井在薄地層中的響應(yīng)進(jìn)行了模擬計(jì)算［9］。1993年，Lovell發(fā)展了一種基于頻率相關(guān)的電位偏微分方程模擬方法［10］。1998年，Chen等人通過有限元的數(shù)值模擬方法檢驗(yàn)了一種陣列側(cè)向測井模型的分辨率和精度［11］。1999年，Liu等人研究了一種鹽水泥漿動(dòng)態(tài)侵入的雙側(cè)向測井響應(yīng)的模型并應(yīng)用于中國渤海灣油田［12］。2002年，劉振華、胡啟利用有限元法對陣列側(cè)向電極系在不同地層模型下的測井響應(yīng)進(jìn)行了數(shù)值模擬，討論了井眼效應(yīng)、侵入特性和層厚影響等問題，對電極系儀器常數(shù)的計(jì)算也進(jìn)行了研究［13］，Galli等人介紹了一種由Schlumberger和ENI-Agip共同開發(fā)的陣列側(cè)向二維電阻率正反演技術(shù)在挪威近海一處碎屑巖儲(chǔ)集層的應(yīng)用情況［14］。2006 年，Michael等人基于 LRI(Localized Rapid Inversion)新方法提出了一種針對陣列側(cè)向測井儀器的實(shí)時(shí)解釋技術(shù)［15］。2007年，Cozzolino等人研究了軸對稱地層模型中雙側(cè)向儀的合成聚焦和模擬［16］，Yang等人用二階有限元方法解釋了中東碳酸鹽巖儲(chǔ)層的雙側(cè)向測井?dāng)?shù)據(jù)［17］，朱軍等通過對儀器電極系的改進(jìn)對高分辨率雙側(cè)向測井響應(yīng)進(jìn)行了數(shù)值模擬分析［18］，范曉敏等對雙側(cè)向測井曲線形狀的影響因素進(jìn)行了研究［19］。2008年，仵杰等人利用二維有限元法對陣列側(cè)向測井儀器的響應(yīng)進(jìn)行了數(shù)值模擬，分析了不同地層模型下目的層厚度和侵入深度變化對測井響應(yīng)的影響［20］。2009年，范宜仁等利用有限元方法，在二維非均勻介質(zhì)模型下對陣列側(cè)向測井響應(yīng)進(jìn)行了數(shù)值模擬，研究了不同層厚、侵入帶、井眼條件下陣列側(cè)向測井響應(yīng)，與雙側(cè)向測井響應(yīng)進(jìn)行了對比［21］。

本文采用米字型網(wǎng)格，合理分配網(wǎng)格疏密程度和節(jié)點(diǎn)編號，在保證精度的前提下提高了程序的運(yùn)行速度，通過對多種模型參數(shù)進(jìn)行改變，模擬和分析了各種因素對雙側(cè)向測井響應(yīng)的影響，并首次探討了三種近似處理后的侵入模型下測井響應(yīng)的特征，為儀器的優(yōu)化設(shè)計(jì)及電阻率反演作基礎(chǔ)，為測井資料解釋方面提供參考和指導(dǎo)。

2 雙側(cè)向測井問題的模型建立

雙側(cè)向電極系共包括九個(gè)電極:中間七個(gè)體積較小的環(huán)狀電極，兩端各一個(gè)體積較大的柱狀電極。如圖2所示，A0是主電極，M、M'和N、N'是兩對監(jiān)督電極，A1、A'1是一對屏蔽電極。

圖2 電極系模型

當(dāng)A2、A'2作為一對屏蔽電極時(shí)，電極系為深側(cè)向電極系;當(dāng)A2、A'2作為一對回路電極時(shí)，電極系為淺側(cè)向電極系。以主電極為中心，這四對電極對稱地排列在兩端，每對電極之間用導(dǎo)線連接，因此M和M'電極具有相同的電位，N 和 N'，A1和 A'1，A2和A'2電極也具有相同的特點(diǎn)。所要解決的問題是，在這種雙側(cè)向測井儀器結(jié)構(gòu)下，其視電阻率響應(yīng)是如何受儀器結(jié)構(gòu)、地層參數(shù)和井眼參數(shù)影響的。數(shù)學(xué)形式可以歸結(jié)為:求出一個(gè)連續(xù)而且適當(dāng)光滑的位函數(shù)U，它滿足以下各條件:

(a)在求解區(qū)Ω內(nèi)需要滿足式(1):

(b)方程采用柱坐標(biāo)，其中R為介質(zhì)的電阻率。在恒流電極表面上，需要滿足式(2):

其中Rm為泥漿電阻率，IA為恒流電極所供電流，是一已知量，n為邊界的外法線。

(c)在絕緣邊界，需要滿足式(3):

(d)在無窮遠(yuǎn)邊界，需要滿足式(4)。

雙側(cè)向測井問題的定解形式直接求解起來非常復(fù)雜，可轉(zhuǎn)化為泛函極值的形式，其表述為在X中使泛函J(U)取極值(最小值)的函數(shù)U。

其中IE和UE為電極的電流和電位。如圖3所示，為了保持對稱性，采用米字型網(wǎng)格?？拷酆驮瓲畹貙犹幍木W(wǎng)格劃分得較密，以保證數(shù)值模擬的精確度，圍巖及靠近無窮遠(yuǎn)邊界處的網(wǎng)格劃分得較稀疏，以提高程序計(jì)算效率。

圖3 網(wǎng)格剖分

模型的幾何尺寸比例盡量保持和實(shí)際井下情況一致。編號的合理與否將直接影響程序的運(yùn)行效率，采用的具體實(shí)施過程是:從左至右、自下而上進(jìn)行編號。把單元電導(dǎo)陣Ge安裝到總體電導(dǎo)陣G上。其中，單元電導(dǎo)陣:

其中

{i，j，k}={1，2，3}或{3，1，2}或{2，3，1}。測井時(shí)電流調(diào)節(jié)要保持式(7)成立，

令主電極電流為1，代入式(8)，

求解式(8)，即可得到求解域的電場分布，求解過程采用的是波前法。

主電極電流為1，當(dāng)K已知時(shí)，將求得的監(jiān)督電極電位UM代入式(9)，即可得到最終所要求解的視電阻率Ra。因此，在求解視電阻率值時(shí)，不必求解出整個(gè)電場，只要知道UM即可，這樣可以進(jìn)一步改進(jìn)線性方程組的求解過程，采用改進(jìn)后的波前法，編號時(shí)將所在節(jié)點(diǎn)最后編號，這就可以在邊安裝邊消元的過程完成后，直接得到UM，從而節(jié)省大量回代過程，大大提高了程序的計(jì)算效率。淺雙側(cè)向的情況只要將屏蔽電極A2、A'2改成回流電極即可，具體實(shí)施過程這里不再贅述。

3 雙側(cè)向測井響應(yīng)的數(shù)值模擬

本文通過正演模擬，研究了多種模型下井眼、侵入帶、層厚對雙側(cè)向測井響應(yīng)的影響，探討了不同侵入模式下的測井響應(yīng)特征。模型采用的雙側(cè)向儀器直徑為0.089m，全長為9.36m，電極系尺寸為(單位:m):

計(jì)算得到的深側(cè)向儀器常數(shù)為0.739，淺側(cè)向儀器常數(shù)為1.485。與可查文獻(xiàn)中相同模型的K值0.733(深側(cè)向)和 1.505(淺側(cè)向)基本相同［8］，說明程序正確。

(a)井徑的影響

如圖4a建立模型，取井內(nèi)泥漿電阻率為1Ω·m，無侵入，原狀地層厚度為3m，原狀地層電阻率為50Ω·m，圍巖電阻率為1Ω·m。井徑分別取0.1m、0.2m、0.3m、0.5m，對比深、淺雙側(cè)向響應(yīng)。模擬結(jié)果中，目標(biāo)層淺側(cè)向的響應(yīng)大于相同條件下深側(cè)向的響應(yīng)。

圖4 (LLS為淺側(cè)向，LLD為深側(cè)向)

由圖4b可知，隨著井徑的增大，深、淺側(cè)向響應(yīng)都減小。曲線形狀漸趨平緩，雙峰現(xiàn)象逐漸消失，低阻地層電阻率模擬值與真實(shí)值基本趨于相同。

(b)井內(nèi)泥漿電阻率的影響

如圖5a建立模型，取井徑為0.2 m，侵入帶厚度為0.3m，侵入帶電阻率為2Ω·m，原狀地層厚度為3m，原狀地層電阻率為50Ω·m，圍巖電阻率為1Ω·m。井內(nèi)泥漿電阻率分別取0.1Ω·m、1Ω·m、10Ω·m，對比深、淺雙側(cè)向響應(yīng)。模擬結(jié)果中，目標(biāo)層淺側(cè)向的響應(yīng)大于相同條件下深側(cè)向的響應(yīng)。

圖5 (LLS為淺側(cè)向，LLD為深側(cè)向)

由圖5b可知，隨著泥漿電阻率的增大，深側(cè)向響應(yīng)減小，對應(yīng)的淺側(cè)向響應(yīng)增大。說明淺側(cè)向受泥漿電阻率的影響較大，變化趨勢和泥漿電阻率的變化趨勢保持一致。當(dāng)泥漿電阻率增大時(shí)，深側(cè)向受侵入帶的影響變大，故其模擬值呈下降趨勢。

(c)侵入帶厚度的影響

如圖6a建立模型，取井徑為0.2 m，井內(nèi)泥漿電阻率為1Ω·m，侵入帶電阻率為2Ω·m，原狀地層厚度為3m，原狀地層電阻率為50Ω·m，圍巖電阻率為1Ω·m。侵入帶厚度分別取0m、0.1m、0.3m、0.7m、1.1m、1.5m，對比深、淺雙側(cè)向響應(yīng)。

圖6(LLS為淺側(cè)向;LLD為深側(cè)向)

由圖6b可知，隨著侵入帶厚度的增大，深、淺側(cè)向響應(yīng)都減小。當(dāng)侵入很少時(shí)，雙峰現(xiàn)象比較明顯;隨著侵入加深，曲線形狀漸趨平緩。當(dāng)侵入厚度小于1.5m時(shí)，淺側(cè)向的響應(yīng)大于相同條件下深側(cè)向的響應(yīng);當(dāng)侵入厚度大于1.5m時(shí)，淺側(cè)向的響應(yīng)小于相同條件下深側(cè)向的響應(yīng)。

(d)侵入帶電阻率的影響

如圖7a建立模型，取井徑為0.2 m，井內(nèi)泥漿電阻率為1Ω·m，侵入帶厚度為0.3m，原狀地層厚度為3m，原狀地層電阻率為50Ω·m，圍巖電阻率為1Ω·m。侵入帶電阻率分別取1Ω·m、10Ω·m、30Ω·m、50Ω·m，對比深、淺雙側(cè)向響應(yīng)。模擬結(jié)果中，目標(biāo)層淺側(cè)向的響應(yīng)大于相同條件下深側(cè)向的響應(yīng)。

圖7(LLS為淺側(cè)向，LLD為深側(cè)向)

由圖7b可知，隨著侵入帶電阻率的增大，深、淺側(cè)向響應(yīng)都增大，曲線的雙峰變得比較明顯，在地層交界處出現(xiàn)兩個(gè)突變，低阻地層出現(xiàn)極小值，高阻地層出現(xiàn)極大值。

(e)層厚影響

如圖8a建立模型，取井徑為0.2 m，井內(nèi)泥漿電阻率為1Ω·m，無侵入，原狀地層電阻率為50Ω·m。原狀地層厚度分別取0m、0.2m、0.4m、0.6m、0.8m、1m、2m、3m、4m、5m、6m、7m、8m、9m，圍巖電阻率分別取1Ω·m、3Ω·m、5Ω·m、10Ω·m，對比深、淺雙側(cè)向響應(yīng)。

圖8(LLS為淺側(cè)向，LLD為深側(cè)向)

由圖8b可知，當(dāng)圍巖電阻率相對原狀地層比較低時(shí)，其影響很小;當(dāng)原狀地層厚度為0.8m左右時(shí)，深、淺側(cè)向響應(yīng)出現(xiàn)一個(gè)極大值;當(dāng)原狀地層厚度小于1m時(shí)，深、淺側(cè)向響應(yīng)都很小;當(dāng)原狀地層厚度大于1m時(shí)，深、淺側(cè)向視電阻率值比較能反映原狀地層真實(shí)電阻率值。

(f)侵入模型的三種近似處理

如圖9a建立模型，取井徑為0.2 m，井內(nèi)泥漿電阻率為0.03Ω·m，圍巖電阻率為1Ω·m，侵入地層電阻率分布情況如圖9b所示。

圖9

如圖9c建立三種侵入模型，對比三種侵入模型下的響應(yīng)特征可以看出(圖9d、圖9e)，深側(cè)向測井以終點(diǎn)電阻率為侵入帶電阻率更接近于原狀地層的真實(shí)電阻率，而淺側(cè)向測井以起點(diǎn)電阻率為侵入帶電阻率更能反映侵入帶的真實(shí)電阻率。對比深淺側(cè)向測井響應(yīng)特征可以看到，淺側(cè)向的縱向分辨率大概為3m左右，而深側(cè)向的縱向分辨率要大于常規(guī)測井解釋中所認(rèn)為的3m。如圖9a計(jì)算7層36m深的地層模型，計(jì)算步長為0.2m，共計(jì)算181個(gè)視電阻率值點(diǎn)，在主頻1.81GHz、內(nèi)存2.00GB的單核臺式電腦上的運(yùn)行耗時(shí)大約為4.05s，速度比較快，滿足了生產(chǎn)應(yīng)用中對測井資料精細(xì)解釋的要求，為電阻率反演工作的有效進(jìn)行提供了可能。

4 結(jié)論

采用米字型網(wǎng)格，對稱模型下雙側(cè)向測井響應(yīng)的數(shù)值模擬結(jié)果也保持對稱。模擬7層36m深的地層模型，計(jì)算步長為0.2m，共計(jì)算181個(gè)視電阻率值點(diǎn)，較不分開設(shè)置縱向圍巖截?cái)噙吔绾涂v向侵入帶截?cái)噙吔?，分開設(shè)置縱向圍巖截?cái)噙吔绾涂v向侵入帶截?cái)噙吔?，模擬速度提高了10倍以上。

隨著井徑的增大，曲線形狀漸趨平緩，深、淺側(cè)向響應(yīng)都減小。隨著泥漿電阻率的增大，深側(cè)向響應(yīng)減小，對應(yīng)的淺側(cè)向響應(yīng)增大。隨著侵入帶厚度的增大，深、淺側(cè)向響應(yīng)都減小。當(dāng)侵入厚度小于1.5m時(shí)，淺側(cè)向的響應(yīng)大于相同條件下深側(cè)向的響應(yīng);當(dāng)侵入厚度大于1.5m時(shí)，淺側(cè)向的響應(yīng)小于相同條件下深側(cè)向的響應(yīng)。隨著侵入帶電阻率的增大，深、淺側(cè)向響應(yīng)都增大。當(dāng)圍巖電阻率相對原狀地層比較低時(shí)，其影響很小;當(dāng)原狀地層厚度為0.8m左右時(shí)，深、淺側(cè)向響應(yīng)出現(xiàn)一個(gè)極大值;當(dāng)原狀地層厚度小于1m時(shí)，深、淺側(cè)向響應(yīng)都很小;當(dāng)原狀地層厚度大于1m時(shí)，深、淺側(cè)向視電阻率值比較能反映原狀地層真實(shí)電阻率值。

深側(cè)向測井以終點(diǎn)電阻率為侵入帶電阻率時(shí)，計(jì)算得到的視電阻率更接近原狀地層的真實(shí)電阻率;淺側(cè)向測井以起點(diǎn)電阻率為侵入帶電阻率時(shí)，計(jì)算得到的視電阻率更能反映侵入帶的真實(shí)電阻率。淺側(cè)向的縱向分辨率大概為3m左右，深側(cè)向的縱向分辨率要大于常規(guī)測井解釋中所認(rèn)為的3m。

［1］ Schlumberger C，Schlumberger M，Leonardon E G.Electrical coring:a method of determining bottom-hole data by electrical measurements［J］. Petroleum Transactions，AIME，1934(110):237-272.

［2］ Doll H G.The Laterolog:A New Resistivity Logging Method With Electrodes Using An Automatic Focusing System［J］.Petroleum Transactions，AIME，1951(192):305-316.

［3］ Suau J，Grimaldi P，Poupon A，Souhaite P.The Dual Laterolog-Rxo Tool［J］.Fall Meeting of the Society of Petroleum Engineers of AIME，1972:1-44.

［4］ Irons B M.A Frontal Solution Scheme for Finite Element Analysis.Numer［J］.Meth Eng，1970(2):5-32.

［5］ Coggon J H.Electromagnetic and electrical modeling by the finite element method［J］.J Geophysics，1971，36(1):132-155.

［6］ 李大潛，鄭宋穆，譚永基.有限元素法在電法測井中的應(yīng)用［M］.北京:石油工業(yè)出版社，1980.

［7］ Anderson B，Chang S K.Synthetic induction logs by the finite element method［J］.23th Ann，The Log Analyst，1982:17-26.

［8］ 張庚驥.電法測井(下)［M］.北京:石油工業(yè)出版社，1986.

［9］ Khokhar R W，Johnson WM J.A deep laterolog for ultrathin formation evaluation［J］.30th Ann，Society of Professional Well Log Analysts Logging Symposium，1989:1-10.

［10］ Lovell JR.Finite element method in resistivity logging［D］.Delft:Delft University of Technology，1993.

［11］ Chen Y H，Chew W C，Zhang G J.A novel array laterolog method［J］.39th Ann，The Log Analyst，1998:22-33.

［12］ Liu Z H，Oyang J，Zhang JH.Dynamic duallaterolog responses:model and field applications in the Bohai Gulf of China［J］.J Petroleum Science and Engineering，1999，23(1):1-11.

［13］ 劉振華，胡啟.陣列側(cè)向測井響應(yīng)的計(jì)算及其特征［J］.西安石油學(xué)院學(xué)報(bào)(自然科學(xué)版)，2002，17(1):53-57.

［14］ Galli M T，Gonfalini M，Mele M.Resistivity Modeling of Array Laterolog Tools:An Application in an Offshore Norway Clastic Reservoir［J］.SPE Annual Technical Conference and Exhibition，2002.

［15］ Frenkel M A.Real-Time Interpretation Technology for New Multi-Laterolog Array Logging Tool［J］.SPE Russian Oil and Gas Technical Conference and Exhibition，2006.

［16］ Cozzolino K，Silva J C.Synthetic focusing and simulation of dual laterolog tool inaxisymmetric subsurface models［J］.J Applied Geophysics，2007，61(2):102-110.

［17］ Yang W，Carlos T V，Ridvan A.Interpretation of frequency-dependent dual-laterolog measurements acquiredin middle-east carbonate reservoirs using second-order finite element method［J］.SPWLA 48th Annual Logging Symposium，2007，1-13.

［18］ 朱軍，馮琳偉.高分辨率雙側(cè)向測井響應(yīng)數(shù)值模擬分析［J］.石油地球物理勘探，2007，42(4):457-462.

［19］ 范曉敏.雙側(cè)向測井曲線形狀與地層侵入關(guān)系研究［J］.地球物理學(xué)進(jìn)展，2007，22(1):142-146.

［20］ 仵杰，謝尉尉，解茜草.陣列側(cè)向測井儀器的正演響應(yīng)分析［J］.西安石油大學(xué)學(xué)報(bào)(自然科學(xué)版)，2008，23(1):73-80.

［21］ 范宜仁，蔣建亮，鄧少貴.高分辨率陣列側(cè)向測井響應(yīng)數(shù)值模擬［J］.測井技術(shù)，2009，33(4):333-336.