雙側向測井響應的有限元法模擬

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

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

1 引言

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

圖1 雙側向測井模型

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

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

2 雙側向測井問題的模型建立

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

圖2 電極系模型

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

(a)在求解區Ω內需要滿足式(1):

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

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

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

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

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

其中IE和UE為電極的電流和電位。如圖3所示，為了保持對稱性，采用米字型網格。靠近井眼和原狀地層處的網格劃分得較密，以保證數值模擬的精確度，圍巖及靠近無窮遠邊界處的網格劃分得較稀疏，以提高程序計算效率。

圖3 網格剖分

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

其中

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

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

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

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

3 雙側向測井響應的數值模擬

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

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

(a)井徑的影響

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

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

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

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

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

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

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

(c)侵入帶厚度的影響

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

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

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

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

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

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

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

(e)層厚影響

如圖8a建立模型，取井徑為0.2 m，井內泥漿電阻率為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，對比深、淺雙側向響應。

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

由圖8b可知，當圍巖電阻率相對原狀地層比較低時，其影響很小;當原狀地層厚度為0.8m左右時，深、淺側向響應出現一個極大值;當原狀地層厚度小于1m時，深、淺側向響應都很小;當原狀地層厚度大于1m時，深、淺側向視電阻率值比較能反映原狀地層真實電阻率值。

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

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

圖9

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

4 結論

采用米字型網格，對稱模型下雙側向測井響應的數值模擬結果也保持對稱。模擬7層36m深的地層模型，計算步長為0.2m，共計算181個視電阻率值點，較不分開設置縱向圍巖截斷邊界和縱向侵入帶截斷邊界，分開設置縱向圍巖截斷邊界和縱向侵入帶截斷邊界，模擬速度提高了10倍以上。

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

深側向測井以終點電阻率為侵入帶電阻率時，計算得到的視電阻率更接近原狀地層的真實電阻率;淺側向測井以起點電阻率為侵入帶電阻率時，計算得到的視電阻率更能反映侵入帶的真實電阻率。淺側向的縱向分辨率大概為3m左右，深側向的縱向分辨率要大于常規測井解釋中所認為的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］ 李大潛，鄭宋穆，譚永基.有限元素法在電法測井中的應用［M］.北京:石油工業出版社，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］.北京:石油工業出版社，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］ 劉振華，胡啟.陣列側向測井響應的計算及其特征［J］.西安石油學院學報(自然科學版)，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］ 朱軍，馮琳偉.高分辨率雙側向測井響應數值模擬分析［J］.石油地球物理勘探，2007，42(4):457-462.

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

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

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