王磊，劉英明，王才志，范宜仁，巫振觀

（1. 中國石油大學（華東）地球科學與技術學院，山東青島 266580；2. 海洋國家實驗室海洋礦產(chǎn)資源評價與探測技術功能實驗室，山東青島 266071；3. 中國石油大學CNPC測井重點實驗室，山東青島 266580；4. 中國石油勘探開發(fā)研究院，北京 100083）

0 引言

隨鉆電磁波測井是水平井地質導向的常用工具[1-2]，提供的相位差和幅度比電阻率曲線常出現(xiàn)分離和“犄角”等復雜響應，導致其定量乃至定性解釋困難。同時，隨鉆時提供的測量曲線數(shù)量有限且缺乏方位識別能力，也進一步增加了地層界面預測的難度[3]。

目前，水平井隨鉆電磁波測井資料的處理主要基于一維（1D）水平層狀地層模型[4-5]和梯度算法來進行的，國外學者已經(jīng)實現(xiàn)了對水平井資料的井斜校正和各向異性地層電阻率反演[5-6]。然而，該資料處理方法需依賴其他測井資料提供的界面位置信息，故僅適用于鉆后評價。同時，梯度算法初值依賴性強，導致解釋結果可能不唯一。為實時確定井周界面位置，水平井導向或儲集層評價中普遍采用交互式解釋方法，即不斷人工調整界面位置并進行正演，以獲得模擬數(shù)據(jù)與實測數(shù)據(jù)一致性高的地層模型[7]。通過引入導眼井/鄰井測井資料和已知地質信息等約束，交互式解釋方法有效地降低了水平井測井解釋的不確定性。但是該方法是基于試錯法的手動尋優(yōu)，理論上其工作量大且處理精度有限，如何提升現(xiàn)有解釋方法的效率和地層界面的預測精度是目前隨鉆測井亟待解決的關鍵問題。

快速 1D正演是水平井隨鉆電磁波測井資料處理的核心，其難點之一在于核函數(shù)的高效積分[8-12]。水平井測井資料處理過程中，積分核函數(shù)振蕩性強、收斂速度慢，直接積分方法不再適用。目前業(yè)界普遍采用快速Hankel變換（FHT）方法正演隨鉆電磁波測井響應[13]。然而，該方法仍存在諸多局限：發(fā)射與接收縱向位置相同時，存在積分奇點[14]；積分精度難以控制；不同接收點的場需進行重復推導。因此，探索新的積分方法，避免陣列線圈響應的重復計算，是提高隨鉆電磁波測井正演速度和精度的必然要求。

本文首先研究了水平井隨鉆電磁波測井正演模型構建方法；然后，提出層界面最優(yōu)匹配技術和索末菲高效直接積分新方法，實現(xiàn)了隨鉆電磁波測井正演加速；然后，將交互式模型調整策略和梯度尋優(yōu)算法相結合，形成基于交互式反演的水平井隨鉆電磁波測井地層界面預測方法；最后，將正反演方法應用于現(xiàn)場資料處理，取得了良好的應用效果。

1 隨鉆電磁波測井正演模型構建

隨鉆電磁波測井采用單發(fā)雙收的軸向天線結構，可同時提供幅度比和相位差視電阻率曲線。該類測井儀器的源距通常小于1.5 m，能夠探測數(shù)米內的地層。一般而言，井周附近各地層近似認為相互平行[15-17]。此時，正演模型中應考慮井眼、儀器、地層和儀器的相對夾角θ等，即圖1a所示的3D正演模型。正演模擬時，該模型可簡化為一維。這主要基于兩點考慮：①隨鉆電磁波測井響應通常受井眼影響較小[4]；②井眼影響可經(jīng)相關校正進行消除。

為驗證上述正演模型降維方法的可行性，建立兩層地層模型，其中儀器以80°的相對夾角穿過地層。當模型中有/無井眼存在時，基于3D/1D算法計算的視電阻率曲線如圖2a、2b所示。模擬時，井眼半徑為0.101 6 m，井內分別填充油基鉆井液（電阻率為1 000 Ω·m）、淡水鉆井液（電阻率為1.0 Ω·m）和鹽水鉆井液（電阻率為0.1 Ω·m）。由圖可知，油基鉆井液填充的鉆井的隨鉆電磁波測井響應和無井眼1D測井響應一致，此時降維方法可行。相比而言，使用水基鉆井液的鉆井的視電阻率值與1D響應差異明顯，且鉆井液電阻率越低差異越大。此時，建立考慮對應井眼情況的刻度圖版，對隨鉆電磁波測井響應進行井眼校正，校正后的結果見圖 2c、2d?？梢钥闯?，井眼校正后的測井響應與1D模擬結果基本重合，這表明模型降維正演方法同樣適用于水基鉆井液填充的鉆井。

圖2 正演模擬校正前后幅度比電阻率與校正前后相位差電阻率

2 隨鉆電磁波測井快速正演方法

2.1 1D層狀各向異性介質偶極子源電磁場遞推算法

在圖1b所示的1D介質中，源激發(fā)的電磁場可通過遞推算法獲得。對隨鉆電磁波測井，接收點的磁場H可表示為水平磁偶極子（HMD）和垂向磁偶極子（VMD）激發(fā)的磁場的疊加，即：

在 1D介質中，VMD只激發(fā)橫電（TE）波，而HMD既產(chǎn)生TE波又激發(fā)橫磁（TM）波。同時，磁偶極子源產(chǎn)生的兩種波相互分離。假定軸向發(fā)射線圈處為地層坐標系原點，當發(fā)射和接收線圈分別位于m和n層（n＞m）時，（1）式可展開為：

為簡化起見，本文僅給出 TM波的傳播系數(shù)表達式，TE波傳播系數(shù)的推導過程可見文獻[9]。對水平磁偶極子激發(fā)的TM波，F(xiàn)TM,h表達式如下：

邊界處，反射系數(shù)和透射系數(shù)的公式可寫為：

2.2 積分核函數(shù)特性

為計算接收線圈的磁場，還需實現(xiàn)（2）式所示的索末菲積分。儀器相對傾角分別為0°，30°，60°，80°和 90°時，圖 3給出了兩層模型中核函數(shù)隨積分變量kρ的變化關系。其中，儀器源距和頻率分別為0.5 m和2 MHz，且發(fā)射線圈固定于地層界面處?？梢钥闯?，隨著傾角的增加，積分核函數(shù)的振蕩性加劇、收斂性降低。當儀器與地層相互平行時，核函數(shù)不收斂，導致總場的計算困難。

總場的振蕩性是由散射場和一次場共同引起，且后者具有解析表達式。圖3b給出了譜域散射場核函數(shù)（hs）的收斂性。與總場相比，散射場的振蕩性變弱，且在水平井時仍收斂。利用相同的地層模型，圖 4進一步給出了水平井散射場積分譜與儀器距地層界面距離（D）的關系圖。隨D的不斷增大，核函數(shù)急劇衰減、收斂性迅速增強。一般而言，相對傾角小于 60°或D大于0.1 m時，積分上限截斷至100，即可保證磁場的計算誤差在0.1%以內。

圖3 不同井斜角度條件下核函數(shù)隨積分變量的變化關系

圖4 散射場的收斂性與儀器距地層界面距離的關系

2.3 散射場直接積分新方法

理論上，散射場可分解為多次反射/透射波的疊加，且多次波的貢獻隨反射/透射次數(shù)的增加而迅速減小。利用該性質，筆者提出一種基于一次反射/透射波扣除的積分新方法。該方法分為 3步：①譜域散射場中扣除一次反射/透射場的貢獻；②被扣除項解析計算；③對剩余譜域核函數(shù)直接積分。散射場Hs來自兩部分，即0階和1階貝塞爾函數(shù)積分。將一次反射/透射場分離，變?yōu)椋?/p>

分離后的積分項變多，且一次反射/透射場的積分項無解析解。問題變?yōu)槿绾螌崿F(xiàn)一次反射/透射波的高效計算。

當kρ趨近于無窮大時，kn,z≈kρ，且廣義反射系數(shù)約等于狹義反射系數(shù)。此時，對TE和TM波的狹義反射系數(shù)進行泰勒展開，即：

基于（13）式和（14）式，可進一步將（11）式變?yōu)椋?/p>

上式中，第2項為一次反射/透射波的逼近。此時，只需對第 1項進行直接積分即可。利用同樣的思路，可處理一階貝塞爾函數(shù)積分，此處不再贅述。

為驗證新方法的有效性，建立圖5a所示的地層模型。當儀器與地層的夾角為90°和89°時，圖5b對比了新方法處理前后散射場的收斂性。由圖5b可知，將一次透射/反射波高階逼近且扣除后，剩余散射場的核函數(shù)隨kρ的增加而急劇衰減。采用新方法后，積分核函數(shù)收斂性可提高3個數(shù)量級以上。此時，積分區(qū)間取[0，50]，即可滿足工程精度要求。表1進一步對比了新的直接積分方法和傳統(tǒng)FHT方法的計算效率及精度。當計算誤差為萬分之一時，F(xiàn)HT方法至少需要260個濾波點[18]，而新方法只需40個高斯勒讓德積分（GLQ）點。本文給出的計算方法精度更高，且速度是傳統(tǒng)FHT方法的6倍以上。同時，新方法的速度基本不受接收線圈數(shù)量的影響，更適用于陣列化隨鉆電磁波測井儀器的模擬。

圖5 地層模型（a）及其新方法處理前后散射場收斂性（b）

表1 索末菲積分方法對比（θ=89°）

2.4 模型邊界自適應匹配方法

模型的地層界面數(shù)是制約隨鉆電磁波測井正演速度的關鍵之一，且兩者基本成線性反比關系[15]。地質導向模型通常包含眾多地層，故正演中必須進行適當截取。目前多采用固定閾值的方法截斷模型邊界，如正演中只考慮距儀器10 m以內的地層界面。由于不同電阻率和頻率條件下儀器探測范圍不一，該邊界截斷方法可能會導致正演精度不夠或計算耗時等問題。為此，本節(jié)提出了一種模型邊界自適應匹配方法。

在均勻地層中，電磁場的衰減通常用趨膚深度δ表示：

1D介質中，源激發(fā)的電磁場在界面處的衰減率s是趨膚深度、層厚和透射系數(shù)的函數(shù)。假定源在m層，儀器上部l層界面處s可近似為：

當s小于1%時，電磁波基本衰減完畢，更遠處地層對測井響應影響極小，故正演中取s=1%為模型截斷邊界。為驗證上述方法的可行性，建立圖6a所示各向異性俄克拉荷馬模型[9]。儀器自上而下穿過地層時，采用不同邊界截斷方法計算的幅度比電阻率響應和相對誤差見圖6?？梢钥闯?，當截斷距離固定為8 m時，傳統(tǒng)方法的模擬結果與精確解一致。將截斷距離減小至4 m時，傳統(tǒng)方法的模擬精度急劇下降。特別是在橫向深度25～170 m處，傳統(tǒng)方法的相對誤差常高于5%。相比而言，采用自適應邊界匹配方法的計算結果相對誤差可控制在0.5%以內，這表明新方法準確、可靠。

圖6 不同截斷方法幅度比電阻率響應與計算精度對比（儀器源距和頻率分別為1.016 m和400 kHz）

為進一步說明新方法的優(yōu)勢，圖 7對比了兩種邊界截斷方法所需的正演模型層數(shù)。與傳統(tǒng)方法相比，新方法能夠對鉆遇地層進行自適應截斷，正演模型層數(shù)急劇減少。采用新方法后，隨鉆電磁波測井模擬速度可提升一倍以上。采用Intel（R）i7-8700 CPU計算時，本文給出的正演方法每秒可計算約16 000個測井點，能夠滿足水平井隨鉆電磁波測井實時交互反演的需要。

圖7 不同邊界截斷方法的正演模型層數(shù)對比

3 水平井隨鉆電磁波測井交互式反演方法

3.1 隨鉆電磁波測井交互反演方法

為確定井周地層界面位置，本節(jié)給出了一種水平井隨鉆電磁波測井交互式反演方法，具體流程見圖8。該方法的核心為鄰井建模、反演模型人工調整和地層界面梯度尋優(yōu)。鄰井建模是指基于導眼井/鄰井信息，提取地層電阻率與巖性序列，以建立參考導向/解釋模型。將參考模型與鄰近測井點的先驗約束（已知的反演結果）相結合，可以確定初始反演模型。然后，利用梯度算法和1D快速正演程序，對地層界面尋優(yōu)。若反演結果不滿足精度誤差，則不斷手工調整地層界面位置進行梯度尋優(yōu)，直至模擬與實測結果相吻合。

圖8 水平井隨鉆電磁波測井反演流程圖

3.2 邊界反演算法與初值選取

隨鉆電磁波測井反演可轉換為求實測資料Ra與模擬響應S的最小二乘問題，反演目標函數(shù)滿足下式：

更新方式如下：

為獲取目標函數(shù)的最小值，采用正則化 Gauss-Newton方法求取。在第t次迭代時，儀器附近的地層界面位置可用下式確定：

采用正則化梯度方法時，反演結果的精度常取決于反演初值的精度。為準確預測地層界面位置，本文采用交互式多初值反演策略，該策略實現(xiàn)方式如下：①利用鄰近測井點界面信息作初值；②基于初始模型反演結果，結合 Ciflog軟件，可視化手工調整地層并建立新的初始反演模型；③重復步驟②，直至反演結果滿足精度誤差。交互式調整反演初值時，應當遵循以下原則：①視電阻率值遠高于地層序列電阻率值時，則將井眼距地層上或下界面的距離變小；②儀器在高電阻率地層時，若重構響應小于實測值，則適當縮小儀器距地層界面的距離d，反之則增大d。一般而言，對模型進行3～5次交互式調整和梯度迭代，即可獲取準確的地層界面位置。利用1D快速正演算法，該交互式反演方法每秒可處理幾十至幾百個測井點，滿足了隨鉆電磁波測井資料實時處理的需要。

4 應用實例

將交互式正反演方法應用于吉林油田某水平井，該井靶層為粉砂巖且層內夾薄泥巖層。該井水平段地層呈水平展布，地層結構沿橫向變化較為緩慢，井眼與地層界面的相對夾角為60°～133°。根據(jù)鄰井資料提取地層序列，并建立初始地質模型（見圖9）。在橫向距離214.6 m處，自然伽馬曲線急劇下降，而視電阻率值迅速增大至 20 Ω·m，故判斷儀器由此入靶。入靶后，通過電阻率曲線難以直接判斷儀器進出層情況。此時，參考地震剖面地層走向，然后根據(jù)自然伽馬曲線確定井軌跡進/出地層界面的位置。對比第1道內的兩條曲線可知，初步調整后的地質模型的自然伽馬值與實測結果基本吻合。然而，基于該模型正演后的隨鉆電磁波測井響應與實測值差距較大，無法滿足精細解釋的需要。

進一步利用隨鉆電磁波測井反演方法，結合具體儀器模式（斯倫貝謝公司的MCR儀器），以準確獲取井周地層界面位置。經(jīng)逐點交互式反演后，最終確定的水平井解釋模型見圖10。與圖9相比，綜合解釋后的模型更為平滑，且實測相位電阻率曲線與模擬重構響應一致性高，證明了解釋模型的準確性和可靠性。解釋結果表明本井儲集層鉆遇率高，故對全井段進行壓裂。該水平井放噴初期日產(chǎn)油86.5 m3，為純油層。

圖9 基于隨鉆自然伽馬測井資料的水平井測井解釋結果（不同顏色代表不同地層）

圖10 水平井隨鉆電磁波測井交互式反演結果（不同顏色代表不同地層）

5 結論

針對水平井中地層界面位置預測的難題，本文給出了一種隨鉆電磁波測井實時正反演方法，并將其應用于實際資料處理，得到以下認識。

利用散射場一次反射/透射波的逼近與扣除方法，索末菲積分的收斂性可提高 3個數(shù)量級。相比于傳統(tǒng)方法，新積分方法適用于任意傾角情況，其精度高、速度快，且更適用于陣列電測井儀器響應的正演。

模型邊界自適應截斷方法大大減少了正演模型的層數(shù)，平衡了傳統(tǒng)方法計算速度與精度之間的矛盾。將新的積分方法和邊界處理方法相結合后，本文研發(fā)的隨鉆電磁波測井正演算法每秒可計算16 000個測井點以上，滿足了實時正反演的需求。

水平井隨鉆電磁波測井交互式反演方法充分利用了鄰井信息和典型電磁波測井響應規(guī)律等先驗信息和人為經(jīng)驗，能夠實時、準確提取井周地層界面位置，減少了地層解釋的多解性和不確定性。

符號注釋：

An——第n個地層中波場的幅度，無因次；C(m)——反演目標函數(shù)；d——地層層厚，m；D——儀器距離地層界面的垂直距離，m；FTE,h，F(xiàn)TE,v——水平和垂直方向的源在接收點所在層激發(fā) TE波的傳播系數(shù)；FTM,h——水平方向的源在接收點所在層激發(fā)TM波的傳播系數(shù)；h——磁場的積分核函數(shù)，A/m；hs——散射磁場的積分核函數(shù)，A/m；H——磁場強度，A/m；Hpq——p方向的源激發(fā)的q方向的磁場分量，A/m；Hs——磁場的散射場強度，A/m；——與 0階貝塞爾函數(shù)相關的散射磁場強度，A/m；——與1階貝塞爾函數(shù)相關的散射磁場強度，A/m；I——單位矩陣；J——雅克比矩陣；J0——0階貝塞爾函數(shù)；J1——1階貝塞爾函數(shù)；k——波數(shù)，無因次；kρ——積分變量；kn,h,z，kn,v,z——第n層地層水平和垂直方向的波數(shù)；m，n，l——地層序號；m——待反演的參數(shù)矢量；N——模型中地層總數(shù)；mref——已知參考模型；O——取高階無窮小函數(shù)；Ra——實測視電阻率數(shù)據(jù)，Ω·m；RP33H——頻率2 MHz、源距0.838 2 m（33 in）的相位差電阻率，Ω·m；，——zn處水平和垂向源激發(fā)的上行TE波的狹義反射系數(shù)；，——zn處上行和下行TM波的狹義反射系數(shù)；，——zn處上行和下行 TM波的廣義反射系數(shù)；，——zn處水平和垂向源激發(fā)的上行TE波的廣義反射系數(shù)；r——接收點距發(fā)射點的距離，m；s——源激發(fā)的電磁場在界面處的衰減率；S(m)——正演響應，Ω·m；t——迭代次數(shù)；——zn處上行 TM波的狹義透射系數(shù)；z——垂向坐標，m；zn——第n個地層界面的垂向坐標，m；σn,h，σn,v——第n層地層水平和垂向電導率，S/m；ρ——徑向距離，m；χ——正則化系數(shù)；χ0——初始正則化系數(shù)；δ——趨膚深度，m；θ——儀器軸與地層界面法向的相對夾角，（°）；ω——角頻率，rad/s；μ——自由空間的磁導率，H/m；λ——地層電阻率各向異性系數(shù)；ξ1，ζ1——和積分中上行波的幅度項；ξ2，ζ2——和積分中下行波的幅度項；，——一次反射/透射波的幅度；——譜域上行波扣除項系數(shù)；——譜域下行波扣除項系數(shù)；——求向量的2范數(shù)。