隨鉆方位電磁波測井響應快速正演方法與地質導向應用

岳喜洲, 劉天淋, 李國玉, 聶在平, 馬明學， 孫向陽

1 電子科技大學電子科學與工程學院， 成都 611731 2 中國海油田服務股份有限公司油田技術研究院, 北京 101149

0 引言

隨鉆方位電磁波測井技術是大斜度井/水平井實時地質導向的關鍵，在海上和非常規油氣藏勘探開發中發揮著重要作用(Wu et al., 2020a,b; Pardo et al., 2015; Fang et al., 2008; 王磊等, 2018).相比于傳統無方位的隨鉆電阻率儀器，隨鉆方位電磁波測井可實現對地層界面、電阻率等信息實時反演，從而將被動式導向升級為主動式導向，極大地提高了儲層鉆遇率(Wang et al., 2019; Wilson et al., 2019;王磊等, 2020; 胡旭飛等, 2018).如何實現隨鉆方位電磁波測井快速、精確正演，是其實時反演的關鍵.各向異性介質電磁場正演方法可分為三類：一維解析法、二維半解析半數值法、三維數值法.其中，三維數值法可適用于復雜地層條件，可對儀器結構進行精細模擬(Teixeira, 2008)，但其計算量大，復雜度高，不適用于實時地質導向；二維半解析半數值法計算速度快，可應用于直井地層侵入等模型下的儀器響應研究(Wang et al., 2009)，不適用于大斜度井/水平井實時地質導向；一維解析法針對徑向或縱向分層地質模型，對儀器響應解析求解，具有速度快、精度高等優點(肖加奇等, 2013)，是當前隨鉆電磁波測井實時反演的基礎.國際各大油服公司的隨鉆方位電磁波測井儀均搭載了相應的一維正演算法，如斯倫貝謝PeriScope(Li et al., 2005)、阿特拉斯Azitrak(Wang et al., 2007)、哈里伯頓ADR(Bittar et al., 2007)等.

近年來，各向異性介質電磁場解析計算得到快速發展，介質模型由單軸各向異性發展到雙軸各向異性，甚至全張量各向異性(L?seth and Ursin, 2007; Yuan et al., 2010; He et al., 2016).在數理意義上，雙軸各向異性、全張量各向異性介質模型更具有普適性(鄧少貴等, 2020；Hu et al., 2018)，單軸各向異性介質模型可由其退化而得.但對于隨鉆電磁波測井尺度而言，特別是海上隨鉆測井，絕大部分地層表現為各向同性，或單軸各向異性.雙軸各向異性、全張量各向異性介質模型電磁場計算通常需要雙重積分，計算復雜度較高，難以適用于實時反演；而單軸各向異性介質模型電磁場通常為含貝塞爾函數的一重積分，高斯積分、快速漢克爾變換等積分方法均可實現其快速計算(莫平華, 2007; Guptasarma and Singh, 1997; Kong, 2007).因此，目前在隨鉆電磁波測井領域，單軸各向異性地層模型及其快速算法更具代表性與適用性.

針對平面分層各向異性介質麥克斯韋方程組一維求解方法，可分為鏡像源法(Wang and Liu, 2014)、并矢Green函數法(魏寶君等, 2007)、TE波與TM波分離法以及Hertz勢函數法(Zhong et al., 2008；陳赟等，2021)等.當利用邊界條件進一步確定各層中的場系數時，又可分為廣義反射法和系數遞推法(姚東華等, 2010)等.雖然物理意義上場的傳播是收斂的，但在解的推導過程中，中間變量常出現指數增大項，受計算機位數限制，易出現數值上溢或下溢情況.為解決數值計算中的場溢出問題，Zhong等(2008)對系數遞推過程中的傳遞矩陣進行變形，壓制了因層厚等原因造成的數值溢出，取得了一定成效.Hong等(2014)引入振幅消除了TE波、TM波在傳播過程中的指數增大項，使得算法適用于任意多層介質模型，壓制了積分上限等原因造成的數值溢出.在早期工程計算中，為避免因層厚導致數值溢出，有學者將厚層分解為多個電阻率一致的薄層，但該方法將引入“假界面”，從而降低計算效率.因此，基于層狀各向異性介質電磁場全張量一維正演方法，消除工程計算中的數值溢出，可實現隨鉆方位電磁波測井響應收斂、快速正演模擬.

本文基于赫茲勢函數法求解各向異性介質中的麥克斯韋方程組，不同于Zhong等(2008)，對解中的指數增大項進行了改造，并消除系數遞推矩陣中的溢出因子，使得遞推過程完全收斂.然后，將所得的無溢出各向異性介質全張量電磁場快速正演算法應用到隨鉆方位電磁波測井領域，對國產隨鉆方位電磁波測電阻率井儀DWPR信號響應進行分析，原始信號正演分析表明，與國外同類儀器相比，DWPR特有的雙斜線圈系進一步提升了探邊距離，達到了國際領先水平.最后以物理實驗與地質導向實例驗證算法的可靠性與高效性，以及儀器的穩定性與實用性.

1 一維各向異性介質正演解析算法

1.1 層狀各向異性介質中麥克斯韋方程組的解

各向異性層狀介質中若只存在時諧磁偶極子源時(圖1)，麥克斯韋方程組可寫為(時諧因子取e-iω t)：

(1)

(2)

(3)

式中，σh指介質水平方向電導率，σv指介質垂直方向電導率.

引入赫茲勢函數Π，與電磁場E、H關系為：

(4)

(5)

以及洛倫茲規范條件：

(6)

代入式(1)、(2)，得到Π的波動方程后對其進行求解，即可得無限厚介質中勢函數Π通解，層狀各向異性介質中勢函數Π表達式可寫為通解與特解之和.

當磁偶極子源沿z方向發射時，Πx、Πy均為0，Πz可寫為：

×kρJ0(kρρ)dkρ，

(7)

×kρJ0(kρρ)dkρ，

(8)

-ξv,jP′je-λjξv,jz+ξv,jQ′jeλjξv,jz

+S′je-ξh,jz+T′jeξh,jz]J1(kρρ)dkρ，

(9)

圖1 層狀各向異性介質模型Fig.1 Layered anisotropic formation model

1.2 系數收斂遞推法

式(6)、(7)、(8)中，由于z可為任意實數，積分變量kρ∈(0,∞)，故z>0時eξh,jz、eλjξv,jz項為指數增大項，z<0時e-ξh,jz、e-λjξv,jz項為指數增大項，在計算過程中將會存在嚴重的數值溢出問題.因此，考慮到場的收斂性與衰減性，在上述表達式中引入地層邊界位置.磁偶極子源沿z向發射時，將式(7)改寫為：

×kρJ0(kρρ)dkρ，

(10)

磁偶極子源沿x向發射時，式(8)、(9)改寫為：

+Qjeλjξv,j(z-zj))kρJ0(kρρ)dkρ，

(11)

-ξv,jPje-λjξv,j(z-zj-1)+ξv,jQjeλjξv,j(z-zj)

+Sje-ξh,j(z-zj-1)+Tjeξh,j(z-zj)]J1(kρρ)dkρ，

(12)

式中，zj-1、zj分別表示第j-1、j個界面.Aj、Pj、Sj可認為是Hertz勢在第j-1個界面處向上傳播的波幅，e-ξh,j(z-zj-1)、e-λjξv,j(z-zj-1)則代表其傳播的衰減模式；Bj、Qj、Tj可認為是Hertz勢在第j個界面處向下傳播的波幅，eξh,j(z-zj)、eλjξv,j(z-zj)則代表其傳播的衰減模式.

由于第1層為向下延伸的半無限厚層，不存在下界面，故A1、P1、S1均為0；類似的，第N層為向上延伸的半無限厚層，不存在上界面，故BN、QN、TN均為0.相應的電磁場表達式可以此類推.

在第j個界面處，電磁場切向分量連續：Hx,j=Hx,j+1，Hy,j=Hy,j+1，Ex,j=Ex,j+1，Ex,j=Ex,j+1.首先考慮無源區，對于沿z向發射的磁偶極子源：

(13)

可得：

ξh,iAie-ξh,i(zi-zi-1)-ξh,iBieξh,i(zi-zi)=ξh,i+1Ai+1e-ξh,i+1(zi-zi)-ξh,i+1Bi+1eξh,i+1(zi-zi+1)，

(14)

μi(Aie-ξh,i(zi-zi-1)+Bieξh,i(zi-zi))=μi+1(Ai+1e-ξh,i+1(zi-zi)+Bi+1eξh,i+1(zi-zi+1))，

(15)

以hj表示第j層層厚，即hj=zj-zj-1.令：

(16)

(17)

可得系數遞推式：

(18)

同樣的，對于沿x向發射的磁偶極子源：

μi(-ξv,iPie-λiξv,ihi+ξv,iQi)=μi+1(-ξv,i+1Pi+1+ξv,i+1Qi+1e-λi+1ξv,i+1hi+1)，

(19)

(20)

μiSie-ξh,ihi+μiTi=μi+1Si+1+μi+1Ti+1e-ξh,i+1hi+1，

(21)

ξh,iSie-ξh,ihi-ξh,iTi=ξh,i+1Si+1-ξh,i+1Ti+1e-ξh,i+1hi+1，

(22)

令：

(23)

(24)

可得系數遞推式：

(25)

(26)

值得注意的是，式(18)、(25)、(26)中的溢出因子eξh,jhj、eλjξv,jhj等項被單獨提取，以方便觀察、處理.

然后考慮有源區.假設源位于第m層，在界面z=zm處，朝z向、x向發射的磁偶極子源，分別有：

(27)

(28)

(29)

在界面z=zm-1處，有：

(30)

(31)

(32)

結合系數遞推式(12)、(13)、(14)，可得：

(33)

(34)

(35)

1.3 算法驗證

當儀器坐標系與地層坐標系不重合時，已有大量文獻證明地層坐標系下的場與儀器坐標系下的場可由旋轉矩陣R進行轉換，如Zhang等(2004)、鄧少貴等(2021)，這里不再贅述.

本文采用Guptasarma和Singh(1997)所述120/140點快速Hankel變換計算0階/1階貝塞爾函數廣義積分，以實現單軸各向異性介質中電磁場快速計算.本文所述算法與格林函數法(魏寶君等, 2007)、截斷溢出法在Oklahoma地層模型(圖2a)下進行對比，計算結果如圖2b、c所示，其中儀器頻率為2 MHz，源距為101.6 cm，相對井斜角為60°.由圖2可知，截斷溢出法在層厚較大和電阻率變化較大的層位出現明顯誤差，其原因為核函數在該情況下容易溢出，導致積分被截斷.而本文所述方法與格林函數法結果一致，從而證明本文所述方法的正確性.

2 隨鉆方位電磁波測井響應

近年來，隨鉆方位電磁波測井技術得到廣泛應用，極大地提高了大斜度井、水平井儲層鉆遇率，為油氣資源勘探與開發提供了強有力的支撐.國際三大油服公司10余年前便相繼推出隨鉆方位電磁波測井儀，占有國際絕大部分隨鉆探邊市場；國內中海油服于2019年推出了自主研發的商用隨鉆方位電磁波測井儀DWPR，在渤海油氣田等取得良好應用；DWPR儀器結構異于斯倫貝謝、貝克休斯與哈里伯頓的隨鉆探邊儀器，本文結合DWPR儀器結構，從正演角度闡述儀器響應特征及算法應用情況.

圖2 Oklahoma地層模型及Hzz分量響應Fig.2 The Oklahoma formation model and the responses of Hzz

2.1 DWPR原始信號響應

如圖3所示，DWPR儀器中包含軸向發射-軸向接收、軸向發射-傾斜接收、傾斜發射-軸向接收、傾斜發射-傾斜接收四種線圈系結構.在井下測量時，儀器旋轉一周進行采樣，軸向發射-軸向接收線圈系(常規線圈系)接收信號為場zz分量信號，軸向發射-傾斜接收(單斜線圈系)與傾斜發射-軸向接收磁場信號類似，均可用一階三角函數表示(馬明學等, 2018)：

(36)

圖3 DWPR儀器結構Fig.3 The structure of DWPR tool

式中，φ表示儀器自旋角度.傾斜發射-傾斜接收線圈系(雙斜線圈系)磁場信號可用二階三角函數表示：

(37)

基于上述信號構成，建立7層地層模型，層厚均為5 m，以中頻(400 kHz)、長源距線圈系(208.28 cm、243.84 cm)為例，儀器以89°穿過地層，其單斜、雙斜線圈系原始磁場信號幅度、相位響應如圖4所示.圖4中Rh與Rv分別為地層水平、垂直電阻率.

這里假設儀器勻速鉆進，圖4中A82M、P82M表示208.28 cm(單斜)線圈系旋轉一周均勻采樣16個扇區的幅度、相位值在深度上的均勻展開，A96M、P96M表示243.84 cm(雙斜)線圈系旋轉一周均勻采樣16個扇區的幅度、相位值在深度上的均勻展開.由圖5可知，208.28 cm、243.84 cm原始相位、幅度信號在地層邊界處波形幅度明顯增大，并且儀器由高阻入低阻、低阻入高阻時，波形幅度各不相同，故取其信號波形幅度即可作為方位探邊地質信號.在不同電阻率的地層中，波形基值不同，表明其對地層電阻率敏感，故可將其用于反映地層電阻率.

2.2 DWPR地質信號響應

由式(28)、(29)可知，將儀器旋轉測量所得信號經一階、二階三角函數擬合，得到原始信號波形幅度，作為儀器方位探邊信號，如圖5所示.

由圖5可知，208.28 cm、243.84 cm地質信號在地層邊界處達到極值；且從高阻層入低阻層時，信號為正值，從低阻層入高阻層時信號為負值.即208.28 cm、243.84 cm地質信號準確反映了地層邊界位置及方位信息.

為對比單斜線圈系與雙斜線圈系地質信號響應特征，建立兩層地層模型，其電阻率為100 Ωm和1 Ωm，線圈系源距均為243.84 cm，頻率為100 kHz，以89°傾角由100 Ωm地層穿過1 Ωm地層，其幅度原始信號與幅度比地質信號響應如圖6所示.

由圖6a可知，儀器靠近界面時，單斜線圈系與雙斜線圈系幅度原始信號波形幅度均增大，在界面處達到最大值.將原始信號處理為地質信號，并以0.02 dB為門檻值，如圖6b所示，單斜線圈系幅度比地質信號探測距離為6.19 m，雙斜線圈系幅度比地質信號探測距離為6.8 m.由此可說明，雙斜線圈系具有比單斜線圈系更突出的探邊距離.

圖4 208.28 cm、243.84 cm線圈系信號響應Fig.4 The responses of receivers at 208.28 cm and 243.84 cm spacing

圖5 208.28 cm、243.84 cm線圈系地質信號響應Fig.5 The geo-signal of receivers at 208.28 cm and 243.84 cm spacing

圖6 單斜線圈系與雙斜線圈系響應特征對比Fig.6 The responses comparison between monoclinc coils and biclinc coils

圖7 55.88 cm、91.44 cm、152.4 cm線圈系電阻率信號響應Fig.7 The apparent resistivity of receivers at 55.88 cm, 91.44 cm and 152.4 cm spacing

2.3 DWPR電阻率信號響應

由式(28)可知，儀器旋轉均勻采樣，取均值后即為zz分量；將對應相位差、幅度比刻度，儀器可得55.88 cm、91.44 cm、152.4 cm三種尺寸、三種頻率的相位差、幅度比電阻率信號，如圖7所示.

由圖7可知，電阻率信號RA22H、RA36M、RA60L、RP22H、RP36H、RP60L(R代表電阻率；A、P分別代表幅度比、相位差；22、36、60代表源距；H、M、L分別代表高頻、中頻、低頻，即2 MHz、400 kHz、100 kHz)均反映了地層電性特征.同尺寸下，相位差電阻率比幅度比電阻率更易受地層邊界影響，其在界面處犄角更明顯.相較于相位差電阻率，幅度比電阻率探測深度更大，反映地層宏觀電阻率.而相位差電阻率分辨率更高，對地層電阻率反映更接近真值.

3 應用實例

將收斂遞推算法應用于隨鉆方位電磁波測井儀器電阻率刻度與地質信號驗證，以物理實驗角度證明算法及儀器響應的準確性；將算法應用于實時地質導向，則可證明算法的實用性與可靠性.

3.1 地質信號驗證

地質信號驗證實驗實照如圖8所示，儀器水平放置于玻璃鋼架，地面以水泥膠結厚0.5 cm的鋁板，每轉動儀器22.5°后靜止采樣20 min，測量一周后可得地質信號.調整儀器距地面高度后，可得不同高度下的地質信號，如圖9所示，為儀器距地面1.16 m、0.76 m、0.36 m的儀器響應與正演響應對比圖.

圖8 地質信號驗證實驗Fig.8 Geo-signal verification test

由圖9可知，儀器從距鋁板2.5 m處往鋁板移動時地質信號增大，在距鋁板約0.63 m時達到最大值，繼續往鋁板移動時地質信號減小.實測信號與正演響應信號對應良好，只在高度為0.76 m的位置信號略有誤差.該次實驗證明了算法正演響應與儀器地質信號實測響應的正確性.

3.2 電阻率信號驗證

為驗證隨鉆方位電磁波測井電阻率信號，將儀器置于直徑4 m、高9.5 m的大型NaCl溶液罐中(圖10所示)，為減小外界環境影響，罐體外圍采用玻璃鋼架.通過調整NaCl濃度可改變溶液電阻率，從而模擬不同導電能力地層.圖11對比了儀器實測電阻率信號與正演電阻率刻度圖版，可以看出二者吻合，實現了正演計算和儀器性能的相互驗證.

3.3 隨鉆地質導向實例

進一步將正演算法應用于隨鉆方位電磁波測井資料處理，圖12為渤海某井DWPR儀器實時反演結果.圖中，第一道為伽馬曲線，第二道為深度道，第三道為電阻率曲線，第四道為伽馬成像，第五道為井眼軌跡及反演窗簾圖.

圖9 208.28 cm、243.84 cm方位信號驗證Fig.9 The results of 208.28 cm and 243.84 cm geo-signal verification

圖10 電阻率信號驗證環境Fig.10 The environment of apparent resistivity verification

根據鄰井對比，該井目的層下部存在水淹，在水平段導向過程中，為確保后期投產效果，控制在儲層中上部鉆進，保證油柱高度.根據電阻率反演剖面判斷，儲層橫向變化較大，層界面沿井軌跡方向呈現“凹”字形，先下傾后上傾，局部傾角達到2.5°，與地震剖面地層整體下傾的趨勢有所區別.測深2218 m(位置4)，垂深1505.63 m，井斜88.70°，根據測井顯示，電阻率已升高至10 Ωm以上，目的層物性變好，但探邊顯示儲層頂界面繼續下傾，目前軌跡已緊貼頂界面，為防止鉆出儲層，降斜至87.5°鉆進.測深2258 m(位置5)，垂深1506.95 m，井斜87.41°，據探邊顯示，距頂界面約2 m，頂界面下傾趨勢變緩，儲層橫向展布變化，地層變厚，下界面(水淹層界面)超出探邊探測范圍，判斷形成錐狀水淹，水淹層界面非均勻推進，增斜至89°鉆進.實鉆過程中的認識調整井斜角度，尋找優質儲層，至測深2448 m(位置8)處成功完鉆.實時導向過程中，儀器電阻率信號、方位信號響應穩定，搭載本快速正演算法的反演剖面清晰、穩定，體現了儀器響應及算法的正確性和高效性.

4 結論

本文基于赫茲勢函數法求解水平層狀各向異性介質中的電磁場，通過引入地層界面，提取并消除溢出因子，可解決各向異性介質全場張量響應解析解中的數值溢出問題.介紹并分析了國產DWPR儀器設計原理及其地層邊界探測性能.真實作業環境中DWPR儀器地質導向結果證明了基于本文算法的資料快速處理軟件具有穩定性和高效性.

圖11 55.88 cm、152.4 cm電阻率信號驗證Fig.11 The results of 55.88 cm and 152.4 cm apparent resistivity verification

圖12 渤海某井隨鉆地質導向實時反演結果Fig.12 The real-time geo-steering inversion in Bohai base

附錄A

令：

(A1)

(A2)

可得：

(A3)

(A4)

(A5)

(A6)

(A7)

(A8)

AN、B1、PN、Q1、SN和T1表達式已知后，第i層中的系數Ai、Bi、Pi、Qi、Si和Ti表達式可根據系數遞推式推導，從而可獲得各層中場的表達式.