劉福平， 王安玲， 劉華群， 楊長春

1 北京印刷學院， 北京　102600　2 中國科學院地質(zhì)與地球物理研究所， 北京　100029

?

過套管電阻率測井的地層電阻率快速反演算法

劉福平1, 2， 王安玲1， 劉華群1， 楊長春2

1 北京印刷學院， 北京1026002 中國科學院地質(zhì)與地球物理研究所， 北京100029

摘要過套管電阻率測井是通過測量套管壁電勢實現(xiàn)測量地層的視電阻率，基于傳輸線方程理論，針對層狀地層，給出了套管壁電勢、電流對地層橫向電阻導數(shù)的微分方程(稱Jacobi矩陣微分方程)及邊界條件；利用Jacobi矩陣微分方程邊值問題導出了過套管電阻率測井反演地層參數(shù)的Jacobi矩陣系數(shù)的解析表示, 利用Marquardt方法實現(xiàn)了過套管測井的地層電阻率反演；通過計算對Jacobi矩陣的特性進行了探討，并獲得了較快的計算速度(因為Jacobi矩陣是用解析解表示的),反演結果與地層模型取得了較好的逼近.本文實現(xiàn)了過套管電阻率測井地層參數(shù)的Jacobi系數(shù)矩陣的快速計算及地層電阻率反演，為進一步開展電阻率測井數(shù)據(jù)處理提供了理論依據(jù)和快速反演算法.關鍵詞過套管電阻率測井； Jacobi矩陣； 傳輸線方法； 測井響應； 反演計算

1引言

過套管電阻率測井是油氣井在開發(fā)過程中油氣藏動態(tài)監(jiān)測和剩余油地層重新評價的重要測井方法之一(Askey et al., 2002; Kaufman,1990; Kaufman and Wightman,1993; 尹軍強等，1998；高杰等,2008; 劉福平等,2007; Vail, 1993; 謝樹棋等,1999)，但自從20世紀30年代過套管電阻率測井被提出后的幾十年間一直未取得突破性進展(Vail,1989,1991; Vail et al.,1993,1995; Schenkel and Morrision,1994,1990;Xie et al.,1999; 劉福平等,2013)．直到20世紀90年代初Kaufman發(fā)表了基于傳輸線方程的套管井電阻率測井近似理論模型和測量理論后(Benimeli et al., 2002; Kaufman,1990; Kaufman and Wightman,1993; 劉福平等,2013)，在過套管電阻率測井理論上和在測井儀器的研制上才取得了實質(zhì)性進展(Aulia et al., 2001; Zhou et al.,2002；高杰等,2008; 劉福平等, 2007)．然而到目前為止，尚未涉及到過套管電阻率測井的參數(shù)反演問題.由于金屬套管的電阻率與地層的電阻率存在巨大的差別(有的相差1010)，致使有限元、有限差分法離散方程的矩陣系數(shù)大小存在巨大差別，用這兩種方法實現(xiàn)套管電阻率測井的正演計算仍然存在困難(目前還沒見到成功的算例)，傳輸線方法已成為過套管電阻率測井成熟的正演計算方法(高杰等,2008; Fondyga et al.,2004；劉福平等,2013)，因此過套管電阻率測井參數(shù)反演的正演計算利用傳輸線方法目前應為一種好的選擇(目前過套管電阻率測井的地層參數(shù)反演及Jacobi矩陣的計算尚未見到算法報道).要實現(xiàn)地層參數(shù)的反演關鍵是如何實現(xiàn)Jacobi矩陣的計算，反演計算時間的長短和精度取決于Jacobi矩陣的計算方法(魏寶君和張庚驥,2002; 丁繼才等,2007; 張美根等，2003; 刑光龍等，2007；劉福平和楊長春,2003; Liu et al.,2012)，因此針對層狀地層開展Jacobi矩陣算法研究具有重要意義.本文研究了套管壁電勢、電流對地層橫向電阻導數(shù)的Jacobi矩陣計算方法，導出了過套管電阻率測井反演地層參數(shù)Jacobi矩陣的解析表示，并探討了Jacobi矩陣的特性，為進一步實現(xiàn)電阻率測井資料反演和測井數(shù)據(jù)解釋都將具有實際意義.

2場分布的正演計算

設金屬套管單位長度的電阻為Rc，金屬套管單位長度所對應的地層的橫向電阻為T(漏電電阻)，則傳輸線方程可寫為(高杰等,2008; 劉福平等,2007)

(1)

(2)

(3)

其中

(4)

考慮n+1層地層后有

(5)

由于在n+1層中z可以取無限遠，為保證電流I(z)有限，應取An+1=0， 設電源在坐標原點，則I1(0)=I0，其中I0為由電源流出而流向上半個空間的電流.利用電流源條件和方程(5)可得Bn+1、A1、B1，再利用公式(3)可得任意一層中的解.

3過套管電阻率測井的Jacobi矩陣

(6)

矩陣形式為

(7)

其中

(8)

(9)

(8)式即為過套管電阻率測井觀測數(shù)據(jù)對地層模型參數(shù)的導數(shù)(稱Jacobi矩陣或敏感系數(shù)).

4Jacobi矩陣的計算方法

為了利用解析結果計算Jacobi矩陣，設觀測數(shù)據(jù)也有Nz個(這樣可以得到計算Jacobi矩陣的解析結果，利用這個解析解可以計算任意一個觀測數(shù)據(jù)的Jacobi矩陣)，與網(wǎng)格個數(shù)相同.將方程(1)對Tj求導得第p層中電勢、電流對Tj的導數(shù)方程

(10)

其中

(11)

在縱向多層狀地層的界面電流、電勢是連續(xù)的(方程同時對地層參數(shù)求導等式仍成立)，因此由于某一網(wǎng)格ΔTj的微小變化，在相鄰兩層中所引起的變化ΔU、ΔI(電荷守恒)也應該相等，所以Gpj、Dpj在地層界面是連續(xù)的(Zhangetal.,2003)，即Gp j(dp)=G(p+1)j(dp)，Dp j(dp)=D(p+1)j(dp).

(12)

(13)

(14)

(15)

由(13)(14)及(15)式得

(16)

實際上當觀測數(shù)據(jù)少于網(wǎng)格個數(shù)時，因為p只是取整數(shù)，當p≠j時(第p、j層不在同一地層)，方程(10)為齊次方程，其通解為

(17)

由(13)式得

(18)

當p=j時(即第p、j層為同一地層)，觀測點和方程(10)為非齊次方程，其特解

(19)

非齊次方程解為

(20)

同理

(21)

其中Ip、Up為方程(1)的解.若縱向有Nz層地層，利用Gpj、Dpj在地層界面連續(xù)可得兩相鄰地層均為齊次方程時的系數(shù)遞推關系

(22)

其中

(23)

在方程(10)中每選一個j則有Nz個微分方程，只有當p=j時，微分方程為非齊次，在該層中由上界面連續(xù)邊界條件(12)及(20)、(21)式得

(24)

=Ap+1ηp+1eαp+1dp-Bp+1ηp+1e-αp+1dp,

(25)

寫成矩陣形式有

(26)

其中

(27)

在第p層中由下界面連續(xù)邊界條件(12)及(20)、(21)式得

(28)

(29)

由(22)式得

(30)

(31)

將(30)、(31)式代入(26)、(28)式得

(32)

由電流源(16)式得

(33)

在(32)、(33)式中只有A1、B1、An+1、Bn+1為未知常數(shù)，考慮到在z→∞時Gpj有限，則An+1=0，利用(32)、(33)式可計算A1、B1、Bn+1，再利用 (26)(28)及(22)式可確定全部待定系數(shù)，從而實現(xiàn)Jacobi矩陣的解析計算.

5過套管電阻率測井的地層電阻率反演方法

設反演目標函數(shù)為

(34)

其中Ud為過套管電阻率測井電勢實際測量結果，U為地層模型電勢響應(由地層模型計算的電勢分布).設f(m)=U-Ud，則f(m)為地層模型電勢響應與過套管電阻率測井電勢實際測量結果所構成的參差函數(shù)，若設mk為m的第k次反演迭代近似.使反演目標函數(shù)極小得迭代公式

(35)

(36)

(37)

其中I為n階單位矩陣，αk為一個正實常數(shù)，根據(jù)計算精度要求通過試算適當選取.

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

6.1Jacobi矩陣特性分析

下面取四層狀地層模型為算例，電極供電電流I0=6 A，電源位于坐標原點，套管電阻率為ρc=2.0×10-7Ωm, 套管半徑a=0.1 m，套管厚度Δa=0.01 m,四層狀地層界面位于dp=15,18,22 m，取觀測數(shù)據(jù)點的坐標為z0=10,17,20,24 m(表1及圖1均使用上述統(tǒng)一參數(shù)).表1的3種情況的差別僅是地層的電阻率不同，表的第2，3，4，5列分別為在z0=10,17,20,24 m測量點電勢對第1，2，3，4層地層橫向電阻的導數(shù)，表1中Case 1四層狀地層的電阻率比較接近(ρb=1.1, 1.0,1.1,1.0 Ωm),表1中Case 3四層狀地層的電阻率差距較大(ρb=20, 1,20,1 Ωm).表1顯示地層電阻率對電勢敏感系數(shù)(電勢對地層橫向電阻的導數(shù))影響較大.圖1為電勢對地層橫向電阻的導數(shù)隨z坐標的變化曲線，由于金屬套管的高導電性的影響(電勢沿套管的變化緩慢，接近線性)，使得電勢敏感系數(shù)隨z坐標的變化基本呈線性關系，曲線1、2、3、4分別為測量點位置變化時(z坐標是變化的)，測量電勢對第1，2，3，4層地層橫向電阻的導數(shù)，圖1中曲線1，2是重合的.在該算例中第4層是一個半無界的空間層，電勢敏感系數(shù)與其他層的差別較大.

表1　電勢對地層橫向電阻的導數(shù)

cond(E)=2.372855356002384，4.141677504127008，3.873081055865031

6.2數(shù)值反演算例

為考察Jacobi矩陣的可靠性及反演方法的可行性,下面以三層狀和五層狀地層實現(xiàn)了過套管地層電阻率反演,除地層參數(shù)外,其余計算條件與前面圖1算例相同.圖2為三層狀地層模型,其初值為15 Ωm,反演的最大相對誤差為1.3%,反演結果與地層模型取得很好的逼近.圖3為五層狀地層模型, 其初值分為2段,第1,2層為8 Ωm,第3,4,5層為13 Ωm,反演結果見表2,其最大相對誤差為2.2%,反演結果與地層模型也取得很好的逼近,驗證了方法的可靠性與可行性.算例說明該方法是可以用于過套管電阻率的反演.

圖1　電勢對地層橫向電阻的導數(shù)曲線Fig.1　Curves of derivations of potential with respect to transverse resistance of formation

地層序號12345模型參數(shù)Rs-mod/Ωm5.0010.0020.0010.0020.00反演初值Rs-st/Ωm8.008.0013.0013.0013.00反演結果Rs-inv/Ωm5.1110.2319.829.7819.79

7結論

本文針對層狀地層，導出了過套管電阻率測井套管壁電勢、電流對地層橫向電阻導數(shù)(Jacobi系數(shù)矩陣)的解析表示，實現(xiàn)了地層參數(shù)反演中Jacobi系數(shù)矩陣的計算，選擇Marquardt反演方法實現(xiàn)了過套管測井的地層電阻率反演.由于Jacobi系數(shù)矩陣計算的解析表示，可獲得較快的計算速度.本文通過算例對Jacobi系數(shù)矩陣的特性進行了探討，結果表明由于金屬套管的高導電性影響使得電勢敏感系數(shù)隨z坐標的變化接近線性關系，且變化緩慢，地層電阻率的差異對電勢敏感系數(shù)有較大影響；本算法是基于傳輸線方程實現(xiàn)的，同樣克服了由于金屬套管的電阻率與地層電阻率存在巨大差別給有限元、有限差分法模擬過金屬套管問題所造成的困難.算例還考察了Jacobi矩陣特性，發(fā)現(xiàn)由Jacobi矩陣組成的反演方程屬病態(tài)方程,針對反演迭代矩陣通過試驗發(fā)現(xiàn)，利用每次反演迭代矩陣每行代數(shù)和的模作為Marquardt反演方法中每行的迭代校正因子,較好地改善了反演迭代系數(shù)矩陣條件數(shù).本文提出了地層參數(shù)反演中Jacobi系數(shù)矩陣的計算方法及地層電阻率反演算法，并用解析的形式實現(xiàn)了Jacobi系數(shù)矩陣的計算，不僅計算精度高而且計算速度快(推導過程沒有新的近似)，反演算例表明,反演結果能較好地逼近地層模型,為進一步的過套管電阻率測井數(shù)據(jù)處理提供了理論依據(jù)和方法.

圖2　地層電阻率反演結果(三層狀地層)Fig.2　Inversion result of formation resistance (3-layer model)

圖3　地層電阻率反演結果(五層狀地層)Fig.3　Inversion result of formation resistance (5-layer model)

References

Askey S, Farag S, Logan J, et al. 2002. Cased Hole Resistivity Measurements Optimize Management of Mature Waterflood in Indonesia. ∥ SPWLA 43rd Annual Logging Symposium. Oiso, Japan: Society of Petrophysicists and Well-Log Analysts.

Aulia K, Poernomo B, Richmond W, et al. 2001. Resistivity behind casing.OilfieldReview, 13(1): 1-25.

Benimeli D, Levesque C, Rouault G, et al. 2002. A new technique for faster resistivity measurements in cased holes. ∥ SPWLA 43rd Annual Logging Symposium. Oiso, Japan: Society of Petrophysicists and Well-Log Analysts.

Ding J C, Chang X, Liu Y K, et al. 2007. Layer by layer waveform inversion of seismic reflection data.ChineseJ.Geophys. (in Chinese), 50(2): 574-580.

Fondyga A, Sherba G. 2004. Cased hole density, neutron, sonic and resistivity logging in South America, case histories, lessons learned, and a methodology for including the technology in formation evaluation programs.∥ SPWLA 45th Annual Logging Symposium. Noordwijk, Netherlands: Society of Petrophysicists and Well-Log Analysts.Gao J, Liu F P, Bao D Z, et al. 2008. Responses simulation of through-casing resistivity logging in heterogeneous-casing wells.ChineseJ.Geophys. (in Chinese), 51(4): 1255-1261.Kaufman A A. 1990. The electrical field in a borehole with a casing.Geophysics, 55(1): 29-38.

Kaufman A A, Wightman W E. 1993. A transmission-line model for electrical logging through casing.Geophysics, 58(12): 1739-1747.Liu F P, Yang C C. 2003. The numerical calculation of sensitivity coefficients of permeability field.ChineseJ.Geophys. (in Chinese), 46(1): 131-137.Liu F P, Gao J, Sun B D, et al. 2007. Forward calculation of the resistivity logging response through casing by transmission line equation for multi-layer formations.ChineseJ.Geophys. (in Chinese), 50(6): 1905-1913.Liu F P, Meng X J, Xiao J Q, et al. 2012. Applying accurate gradients of seismic wave reflection coefficients (SWRC) to the inversion of seismic wave velocities.ScienceChinaEarthSciences, 55(12): 1953-1960.

Liu F P, Chen X A, Zhao B C, et al. 2013. Instability analysis of through casing resistivity logging responses in higher resistance formation.WellLoggingTechnology(in Chinese), 37(1): 85-89.Schenkel C J, Morrision H F. 1994. Electrical resistivity measurement through metal casing.Geophysics, 59(7): 1072-1082.Schenkel C J, Morrision H F. 1990. Effects of well casing on potential field measurements using downhole current sources.GeophysicalProspecting, 38(6): 663-686.

Vail W B. 1989. Methods and apparatus for measurement of electronic properties of geological formations through borehole casing. U. S. Patent No. 4882542. November 21.

Vail W B. 1991. Methods and apparatus for measurement of the resistivity of geological formations from within cased well in presence of acoustic and magnetic energy sources. U. S. Patent No. 5043669. August 27.

Vail W B. 1993. Measuring resistivity changes from within a first cased well to monitor fluids injected into oil bearing geological formations from a second cased well while passing electrical current between the two cased wells. U. S. Patent No. 5187440. February 16.Vail W B, Momii S T, Woodhouse R, et al. 1993. Formation resistivity measurements through metal casing.∥ SPWLA 34th Annual Logging Symposium. Calgary, Alberta, Canada: Society of Petrophysicists and Well-Log Analysts.

Vail W B, Momii S T, Haines H, et al. 1995. Formation resistivity measurements through metal casing at the MWX-2 well in Rifle, Colorado.∥ SPWLA 36th Annual Logging Symposium. Paris, France: Society of Petrophysicists and Well-Log Analysts.Wei B J, Zhang G J. 2002. Forward modeling and inversion of 3-D cross-hole electromagnetic fields.ChineseJ.Geophys. (in Chinese), 45(5): 735-743. Xie S Q, Chu Z T, Li K P, et al. 1999. On methodology of resistivity logging through casing.WellLoggingTechnology(in Chinese), 23(5): 338-343.

Xing G L, Yang S D, Li S G. 2007. On the fast algorithm and characteristic analysis of Jacobi matrix in the inversion of electromagnetic wave logs.ChineseJ.Geophys. (in Chinese), 50(2): 642-650.

Yin J Q, Feng Q N, Zhu L D, et al. 1998. Feasibility study on formation resistivity measurement through metal casing.WellLoggingTechnology(in Chinese), 22(5): 376-379.

Zhang M G, Wang M Y, Li X F, et al. 2003. Full wavefield inversion of anisotropic elastic parameters in the time domain.ChineseJ.Geophys. (in Chinese), 46(1): 94-100.

Zhou Q, Julander D, Penley L. 2002. Experiences with casedhole resistivity logging for reservoir monitoring. ∥ SPWLA 43rd Annual Logging Symposium. Oiso, Japan: Society of Petrophysicists and Well-Log Analysts.

附中文參考文獻

丁繼才, 常旭, 劉伊克等. 2007. 反射地震數(shù)據(jù)的逐層波形反演. 地球物理學報, 50(2): 574-580.

高杰, 劉福平, 包德洲等. 2008. 非均勻套管井中的過套管電阻率測井響應. 地球物理學報, 51(4): 1255-1261.

劉福平, 楊長春. 2003. 滲透率場敏感系數(shù)的數(shù)值計算. 地球物理學報, 46(1): 131-137.

劉福平, 高杰, 孫寶佃等. 2007. 實際井眼條件下過套管電阻率測井響應的傳輸線方程正演算法. 地球物理學報, 50(6): 1905-1913.

劉福平, 陳小安, 趙寶成等. 2013. 過套管電阻率測井高電阻率層段不穩(wěn)定性分析. 測井技術, 37(1): 85-89.

魏寶君, 張庚驥. 2002. 三維井間電磁場的正反演計算. 地球物理學報, 45(5): 735-743.

謝樹棋, 儲昭坦, 李克沛等. 1999. 套管井電阻率測井方法研究. 測井技術, 23(5): 338-343.

刑光龍, 楊善德, 李曙光. 2007. 電磁波測井資料反演中Jacobi矩陣的快速算法及其特性分析. 地球物理學報, 50(2): 642-650.

尹軍強, 馮啟寧, 朱龍德等. 1998. 通過金屬套管測量地層電阻率的可行性研究. 測井技術, 22(5): 376-379.

張美根, 王妙月, 李小凡等. 2003. 時間域全波場各向異性彈性參數(shù)反演. 地球物理學報, 46(1): 94-100.

(本文編輯何燕)

基金項目北京市自然科學基金重點項目B類(KZ201510015015)；北京市自然科學基金(4142016)；北京市教委項目 (KM201510015009，KM201410015006，PXM2016_014223_000025)資助.

作者簡介劉福平，男，1960年生，教授，1994年于石油大學勘探系獲碩士學位，2002年于中國科學院地質(zhì)與地球物理研究所獲博士學位，2004年于中國科學院地質(zhì)與地球物理研究所博士后出站，主要從事電磁波傳播理論、數(shù)值計算及油藏數(shù)值模擬方面的研究和教學工作. E-mail: fupingliu@bigc.edu.cn, fupingliu60@sina.com.cn

doi:10.6038/cjg20160634 中圖分類號P631

收稿日期2014-12-29，2015-12-05收修定稿

A fast inversion algorithm of formation resistivity for resistivity logging through casing

LIU Fu-Ping1,2， WANG An-Ling1， LIU Hua-Qun1， YANG Chang-Chun2

1BeijingInstituteofGraphicCommunication,Beijing102600,China2InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China

AbstractThe resistivity logging through casing is that the potential distribution on the metal casing is measured to realize the measurement of formation resistivity. For this method, this study derived the differential equations of Jacobi matrix from the transmission line equation and its boundary conditions for a multi-layered formation, which are the derivative equations of the potential and current with respect to the transverse resistance of formation. With the differential equations we have given the analytic formula of the Jacobi matrix and realized the inversion of formation resistivity (Marquardt inversion method). With computing examples, the characteristics of Jacobi matrix were discussed, and the fast computing speed was obtained, in which the results of inverse are in excellent agreement with the model of formation. The realization of formation resistivity inversion and its fast computation of Jacobi matrix provide a theoretical basis and a fast inversion algorithm for the further development of the resistivity logging and data processing.

KeywordsResistivity logging through casing; Jacobi matrix; Transmission line approach; Logging response; Inversion calculation

劉福平， 王安玲， 劉華群等. 2016. 過套管電阻率測井的地層電阻率快速反演算法. 地球物理學報，59(6):2326-2332,doi:10.6038/cjg20160634.

Liu F P, Wang A L, Liu H Q, et al. 2016. A fast inversion algorithm of formation resistivity for resistivity logging through casing.ChineseJ.Geophys. (in Chinese),59(6):2326-2332,doi:10.6038/cjg20160634.