劉福平, 王安玲, 劉華群, 楊長春
1 北京印刷學(xué)院, 北京 102600 2 中國科學(xué)院地質(zhì)與地球物理研究所, 北京 100029
?
過套管電阻率測(cè)井的地層電阻率快速反演算法
劉福平1, 2, 王安玲1, 劉華群1, 楊長春2
1 北京印刷學(xué)院, 北京1026002 中國科學(xué)院地質(zhì)與地球物理研究所, 北京100029
摘要過套管電阻率測(cè)井是通過測(cè)量套管壁電勢(shì)實(shí)現(xiàn)測(cè)量地層的視電阻率,基于傳輸線方程理論,針對(duì)層狀地層,給出了套管壁電勢(shì)、電流對(duì)地層橫向電阻導(dǎo)數(shù)的微分方程(稱Jacobi矩陣微分方程)及邊界條件;利用Jacobi矩陣微分方程邊值問題導(dǎo)出了過套管電阻率測(cè)井反演地層參數(shù)的Jacobi矩陣系數(shù)的解析表示, 利用Marquardt方法實(shí)現(xiàn)了過套管測(cè)井的地層電阻率反演;通過計(jì)算對(duì)Jacobi矩陣的特性進(jìn)行了探討,并獲得了較快的計(jì)算速度(因?yàn)镴acobi矩陣是用解析解表示的),反演結(jié)果與地層模型取得了較好的逼近.本文實(shí)現(xiàn)了過套管電阻率測(cè)井地層參數(shù)的Jacobi系數(shù)矩陣的快速計(jì)算及地層電阻率反演,為進(jìn)一步開展電阻率測(cè)井?dāng)?shù)據(jù)處理提供了理論依據(jù)和快速反演算法.關(guān)鍵詞過套管電阻率測(cè)井; Jacobi矩陣; 傳輸線方法; 測(cè)井響應(yīng); 反演計(jì)算
1引言
過套管電阻率測(cè)井是油氣井在開發(fā)過程中油氣藏動(dòng)態(tài)監(jiān)測(cè)和剩余油地層重新評(píng)價(jià)的重要測(cè)井方法之一(Askey et al., 2002; Kaufman,1990; Kaufman and Wightman,1993; 尹軍強(qiáng)等,1998;高杰等,2008; 劉福平等,2007; Vail, 1993; 謝樹棋等,1999),但自從20世紀(jì)30年代過套管電阻率測(cè)井被提出后的幾十年間一直未取得突破性進(jìn)展(Vail,1989,1991; Vail et al.,1993,1995; Schenkel and Morrision,1994,1990;Xie et al.,1999; 劉福平等,2013).直到20世紀(jì)90年代初Kaufman發(fā)表了基于傳輸線方程的套管井電阻率測(cè)井近似理論模型和測(cè)量理論后(Benimeli et al., 2002; Kaufman,1990; Kaufman and Wightman,1993; 劉福平等,2013),在過套管電阻率測(cè)井理論上和在測(cè)井儀器的研制上才取得了實(shí)質(zhì)性進(jìn)展(Aulia et al., 2001; Zhou et al.,2002;高杰等,2008; 劉福平等, 2007).然而到目前為止,尚未涉及到過套管電阻率測(cè)井的參數(shù)反演問題.由于金屬套管的電阻率與地層的電阻率存在巨大的差別(有的相差1010),致使有限元、有限差分法離散方程的矩陣系數(shù)大小存在巨大差別,用這兩種方法實(shí)現(xiàn)套管電阻率測(cè)井的正演計(jì)算仍然存在困難(目前還沒見到成功的算例),傳輸線方法已成為過套管電阻率測(cè)井成熟的正演計(jì)算方法(高杰等,2008; Fondyga et al.,2004;劉福平等,2013),因此過套管電阻率測(cè)井參數(shù)反演的正演計(jì)算利用傳輸線方法目前應(yīng)為一種好的選擇(目前過套管電阻率測(cè)井的地層參數(shù)反演及Jacobi矩陣的計(jì)算尚未見到算法報(bào)道).要實(shí)現(xiàn)地層參數(shù)的反演關(guān)鍵是如何實(shí)現(xiàn)Jacobi矩陣的計(jì)算,反演計(jì)算時(shí)間的長短和精度取決于Jacobi矩陣的計(jì)算方法(魏寶君和張庚驥,2002; 丁繼才等,2007; 張美根等,2003; 刑光龍等,2007;劉福平和楊長春,2003; Liu et al.,2012),因此針對(duì)層狀地層開展Jacobi矩陣算法研究具有重要意義.本文研究了套管壁電勢(shì)、電流對(duì)地層橫向電阻導(dǎo)數(shù)的Jacobi矩陣計(jì)算方法,導(dǎo)出了過套管電阻率測(cè)井反演地層參數(shù)Jacobi矩陣的解析表示,并探討了Jacobi矩陣的特性,為進(jìn)一步實(shí)現(xiàn)電阻率測(cè)井資料反演和測(cè)井?dāng)?shù)據(jù)解釋都將具有實(shí)際意義.
2場(chǎng)分布的正演計(jì)算
設(shè)金屬套管單位長度的電阻為Rc,金屬套管單位長度所對(duì)應(yīng)的地層的橫向電阻為T(漏電電阻),則傳輸線方程可寫為(高杰等,2008; 劉福平等,2007)
(1)
(2)
(3)
其中
(4)
考慮n+1層地層后有
(5)
由于在n+1層中z可以取無限遠(yuǎn),為保證電流I(z)有限,應(yīng)取An+1=0, 設(shè)電源在坐標(biāo)原點(diǎn),則I1(0)=I0,其中I0為由電源流出而流向上半個(gè)空間的電流.利用電流源條件和方程(5)可得Bn+1、A1、B1,再利用公式(3)可得任意一層中的解.
3過套管電阻率測(cè)井的Jacobi矩陣
(6)
矩陣形式為
(7)
其中
(8)
(9)
(8)式即為過套管電阻率測(cè)井觀測(cè)數(shù)據(jù)對(duì)地層模型參數(shù)的導(dǎo)數(shù)(稱Jacobi矩陣或敏感系數(shù)).
4Jacobi矩陣的計(jì)算方法
為了利用解析結(jié)果計(jì)算Jacobi矩陣,設(shè)觀測(cè)數(shù)據(jù)也有Nz個(gè)(這樣可以得到計(jì)算Jacobi矩陣的解析結(jié)果,利用這個(gè)解析解可以計(jì)算任意一個(gè)觀測(cè)數(shù)據(jù)的Jacobi矩陣),與網(wǎng)格個(gè)數(shù)相同.將方程(1)對(duì)Tj求導(dǎo)得第p層中電勢(shì)、電流對(duì)Tj的導(dǎo)數(shù)方程
(10)
其中
(11)
在縱向多層狀地層的界面電流、電勢(shì)是連續(xù)的(方程同時(shí)對(duì)地層參數(shù)求導(dǎo)等式仍成立),因此由于某一網(wǎng)格ΔTj的微小變化,在相鄰兩層中所引起的變化ΔU、ΔI(電荷守恒)也應(yīng)該相等,所以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í)際上當(dāng)觀測(cè)數(shù)據(jù)少于網(wǎng)格個(gè)數(shù)時(shí),因?yàn)閜只是取整數(shù),當(dāng)p≠j時(shí)(第p、j層不在同一地層),方程(10)為齊次方程,其通解為
(17)
由(13)式得
(18)
當(dāng)p=j時(shí)(即第p、j層為同一地層),觀測(cè)點(diǎn)和方程(10)為非齊次方程,其特解
(19)
非齊次方程解為
(20)
同理
(21)
其中Ip、Up為方程(1)的解.若縱向有Nz層地層,利用Gpj、Dpj在地層界面連續(xù)可得兩相鄰地層均為齊次方程時(shí)的系數(shù)遞推關(guān)系
(22)
其中
(23)
在方程(10)中每選一個(gè)j則有Nz個(gè)微分方程,只有當(dāng)p=j時(shí),微分方程為非齊次,在該層中由上界面連續(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→∞時(shí)Gpj有限,則An+1=0,利用(32)、(33)式可計(jì)算A1、B1、Bn+1,再利用 (26)(28)及(22)式可確定全部待定系數(shù),從而實(shí)現(xiàn)Jacobi矩陣的解析計(jì)算.
5過套管電阻率測(cè)井的地層電阻率反演方法
設(shè)反演目標(biāo)函數(shù)為
(34)
其中Ud為過套管電阻率測(cè)井電勢(shì)實(shí)際測(cè)量結(jié)果,U為地層模型電勢(shì)響應(yīng)(由地層模型計(jì)算的電勢(shì)分布).設(shè)f(m)=U-Ud,則f(m)為地層模型電勢(shì)響應(yīng)與過套管電阻率測(cè)井電勢(shì)實(shí)際測(cè)量結(jié)果所構(gòu)成的參差函數(shù),若設(shè)mk為m的第k次反演迭代近似.使反演目標(biāo)函數(shù)極小得迭代公式
(35)
(36)
(37)
其中I為n階單位矩陣,αk為一個(gè)正實(shí)常數(shù),根據(jù)計(jì)算精度要求通過試算適當(dāng)選取.
6數(shù)值算例及分析
6.1Jacobi矩陣特性分析
下面取四層狀地層模型為算例,電極供電電流I0=6 A,電源位于坐標(biāo)原點(diǎn),套管電阻率為ρc=2.0×10-7Ωm, 套管半徑a=0.1 m,套管厚度Δa=0.01 m,四層狀地層界面位于dp=15,18,22 m,取觀測(cè)數(shù)據(jù)點(diǎn)的坐標(biāo)為z0=10,17,20,24 m(表1及圖1均使用上述統(tǒng)一參數(shù)).表1的3種情況的差別僅是地層的電阻率不同,表的第2,3,4,5列分別為在z0=10,17,20,24 m測(cè)量點(diǎn)電勢(shì)對(duì)第1,2,3,4層地層橫向電阻的導(dǎo)數(shù),表1中Case 1四層狀地層的電阻率比較接近(ρb=1.1, 1.0,1.1,1.0 Ωm),表1中Case 3四層狀地層的電阻率差距較大(ρb=20, 1,20,1 Ωm).表1顯示地層電阻率對(duì)電勢(shì)敏感系數(shù)(電勢(shì)對(duì)地層橫向電阻的導(dǎo)數(shù))影響較大.圖1為電勢(shì)對(duì)地層橫向電阻的導(dǎo)數(shù)隨z坐標(biāo)的變化曲線,由于金屬套管的高導(dǎo)電性的影響(電勢(shì)沿套管的變化緩慢,接近線性),使得電勢(shì)敏感系數(shù)隨z坐標(biāo)的變化基本呈線性關(guān)系,曲線1、2、3、4分別為測(cè)量點(diǎn)位置變化時(shí)(z坐標(biāo)是變化的),測(cè)量電勢(shì)對(duì)第1,2,3,4層地層橫向電阻的導(dǎo)數(shù),圖1中曲線1,2是重合的.在該算例中第4層是一個(gè)半無界的空間層,電勢(shì)敏感系數(shù)與其他層的差別較大.
表1 電勢(shì)對(duì)地層橫向電阻的導(dǎo)數(shù)
cond(E)=2.372855356002384,4.141677504127008,3.873081055865031
6.2數(shù)值反演算例
為考察Jacobi矩陣的可靠性及反演方法的可行性,下面以三層狀和五層狀地層實(shí)現(xiàn)了過套管地層電阻率反演,除地層參數(shù)外,其余計(jì)算條件與前面圖1算例相同.圖2為三層狀地層模型,其初值為15 Ωm,反演的最大相對(duì)誤差為1.3%,反演結(jié)果與地層模型取得很好的逼近.圖3為五層狀地層模型, 其初值分為2段,第1,2層為8 Ωm,第3,4,5層為13 Ωm,反演結(jié)果見表2,其最大相對(duì)誤差為2.2%,反演結(jié)果與地層模型也取得很好的逼近,驗(yàn)證了方法的可靠性與可行性.算例說明該方法是可以用于過套管電阻率的反演.
圖1 電勢(shì)對(duì)地層橫向電阻的導(dǎo)數(shù)曲線Fig.1 Curves of derivations of potential with respect to transverse resistance of formation
地層序號(hào)12345模型參數(shù)Rs-mod/Ωm5.0010.0020.0010.0020.00反演初值Rs-st/Ωm8.008.0013.0013.0013.00反演結(jié)果Rs-inv/Ωm5.1110.2319.829.7819.79
7結(jié)論
本文針對(duì)層狀地層,導(dǎo)出了過套管電阻率測(cè)井套管壁電勢(shì)、電流對(duì)地層橫向電阻導(dǎo)數(shù)(Jacobi系數(shù)矩陣)的解析表示,實(shí)現(xiàn)了地層參數(shù)反演中Jacobi系數(shù)矩陣的計(jì)算,選擇Marquardt反演方法實(shí)現(xiàn)了過套管測(cè)井的地層電阻率反演.由于Jacobi系數(shù)矩陣計(jì)算的解析表示,可獲得較快的計(jì)算速度.本文通過算例對(duì)Jacobi系數(shù)矩陣的特性進(jìn)行了探討,結(jié)果表明由于金屬套管的高導(dǎo)電性影響使得電勢(shì)敏感系數(shù)隨z坐標(biāo)的變化接近線性關(guān)系,且變化緩慢,地層電阻率的差異對(duì)電勢(shì)敏感系數(shù)有較大影響;本算法是基于傳輸線方程實(shí)現(xiàn)的,同樣克服了由于金屬套管的電阻率與地層電阻率存在巨大差別給有限元、有限差分法模擬過金屬套管問題所造成的困難.算例還考察了Jacobi矩陣特性,發(fā)現(xiàn)由Jacobi矩陣組成的反演方程屬病態(tài)方程,針對(duì)反演迭代矩陣通過試驗(yàn)發(fā)現(xiàn),利用每次反演迭代矩陣每行代數(shù)和的模作為Marquardt反演方法中每行的迭代校正因子,較好地改善了反演迭代系數(shù)矩陣條件數(shù).本文提出了地層參數(shù)反演中Jacobi系數(shù)矩陣的計(jì)算方法及地層電阻率反演算法,并用解析的形式實(shí)現(xiàn)了Jacobi系數(shù)矩陣的計(jì)算,不僅計(jì)算精度高而且計(jì)算速度快(推導(dǎo)過程沒有新的近似),反演算例表明,反演結(jié)果能較好地逼近地層模型,為進(jìn)一步的過套管電阻率測(cè)井?dāng)?shù)據(jù)處理提供了理論依據(jù)和方法.
圖2 地層電阻率反演結(jié)果(三層狀地層)Fig.2 Inversion result of formation resistance (3-layer model)
圖3 地層電阻率反演結(jié)果(五層狀地層)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.
附中文參考文獻(xiàn)
丁繼才, 常旭, 劉伊克等. 2007. 反射地震數(shù)據(jù)的逐層波形反演. 地球物理學(xué)報(bào), 50(2): 574-580.
高杰, 劉福平, 包德洲等. 2008. 非均勻套管井中的過套管電阻率測(cè)井響應(yīng). 地球物理學(xué)報(bào), 51(4): 1255-1261.
劉福平, 楊長春. 2003. 滲透率場(chǎng)敏感系數(shù)的數(shù)值計(jì)算. 地球物理學(xué)報(bào), 46(1): 131-137.
劉福平, 高杰, 孫寶佃等. 2007. 實(shí)際井眼條件下過套管電阻率測(cè)井響應(yīng)的傳輸線方程正演算法. 地球物理學(xué)報(bào), 50(6): 1905-1913.
劉福平, 陳小安, 趙寶成等. 2013. 過套管電阻率測(cè)井高電阻率層段不穩(wěn)定性分析. 測(cè)井技術(shù), 37(1): 85-89.
魏寶君, 張庚驥. 2002. 三維井間電磁場(chǎng)的正反演計(jì)算. 地球物理學(xué)報(bào), 45(5): 735-743.
謝樹棋, 儲(chǔ)昭坦, 李克沛等. 1999. 套管井電阻率測(cè)井方法研究. 測(cè)井技術(shù), 23(5): 338-343.
刑光龍, 楊善德, 李曙光. 2007. 電磁波測(cè)井資料反演中Jacobi矩陣的快速算法及其特性分析. 地球物理學(xué)報(bào), 50(2): 642-650.
尹軍強(qiáng), 馮啟寧, 朱龍德等. 1998. 通過金屬套管測(cè)量地層電阻率的可行性研究. 測(cè)井技術(shù), 22(5): 376-379.
張美根, 王妙月, 李小凡等. 2003. 時(shí)間域全波場(chǎng)各向異性彈性參數(shù)反演. 地球物理學(xué)報(bào), 46(1): 94-100.
(本文編輯何燕)
基金項(xiàng)目北京市自然科學(xué)基金重點(diǎn)項(xiàng)目B類(KZ201510015015);北京市自然科學(xué)基金(4142016);北京市教委項(xiàng)目 (KM201510015009,KM201410015006,PXM2016_014223_000025)資助.
作者簡介劉福平,男,1960年生,教授,1994年于石油大學(xué)勘探系獲碩士學(xué)位,2002年于中國科學(xué)院地質(zhì)與地球物理研究所獲博士學(xué)位,2004年于中國科學(xué)院地質(zhì)與地球物理研究所博士后出站,主要從事電磁波傳播理論、數(shù)值計(jì)算及油藏?cái)?shù)值模擬方面的研究和教學(xué)工作. E-mail: fupingliu@bigc.edu.cn, fupingliu60@sina.com.cn
doi:10.6038/cjg20160634 中圖分類號(hào)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. 過套管電阻率測(cè)井的地層電阻率快速反演算法. 地球物理學(xué)報(bào),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.