孫歧峰，李娜，段友祥，李洪強，唐海全

（1.中國石油大學(xué)（華東）計算機科學(xué)與技術(shù)學(xué)院，山東青島 266580；2.中國石化勝利石油工程有限公司鉆井工藝研究院，山東東營 257000；3.中國石化勝利石油工程有限公司測控技術(shù)研究院，山東東營 257000）

0 引言

隨鉆測井與地層評價技術(shù)通過實時獲取鉆頭附近的軌跡參數(shù)、地質(zhì)參數(shù)和工程參數(shù)來實時評價所鉆遇地層的情況，進而精確控制鉆具運動軌跡以命中最佳地質(zhì)目標(biāo)[1]，是高精度地質(zhì)導(dǎo)向鉆井的關(guān)鍵支撐技術(shù)。隨鉆測井技術(shù)在分辨地層變化及確定巖性[2-5]、評價泥砂/砂泥巖含量[6]、預(yù)測高壓地層[7-8]等方面有明顯效果，能夠提高鉆井工程的效益[9]。由于實時數(shù)據(jù)傳輸?shù)男畔⒎浅Ｓ邢耷医忉岆y度大[10]，將人工智能與隨鉆測井解釋相結(jié)合[11-13]以提高實時數(shù)據(jù)處理的準(zhǔn)確性和效率有重要的意義。

針對隨鉆伽馬測井解釋過程中面臨的問題，結(jié)合長短期記憶神經(jīng)網(wǎng)絡(luò)（LSTM）及小波變換等[14-16]，本文提出一種實時地層傾角解釋方法：首先提取井的層位特征并計算動態(tài)閾值，初步判定地層邊界位置，然后訓(xùn)練分類器模型，篩選真實的地層邊界，最后計算地層相對傾角及方位角。對方法進行詳細闡述并開展現(xiàn)場應(yīng)用。

1 隨鉆方位伽馬數(shù)據(jù)解釋流程

對隨鉆儀器實時傳輸上來的鉆井液脈沖信號進行解碼得到井下各種傳感器的測量數(shù)據(jù)，進行坐標(biāo)轉(zhuǎn)換、時深轉(zhuǎn)換及校正，然后對數(shù)據(jù)進行插值得到方位伽馬的成像結(jié)果。在此基礎(chǔ)上，進行隨鉆方位伽馬數(shù)據(jù)的解釋，流程如下。

①地層邊界初步劃分。首先，采用Lowess（局部加權(quán)回歸散點平滑法）濾波方法對方位伽馬曲線進行過濾，消除與地層變化無關(guān)的隨機噪聲[17]引起的測井曲線中細微尖銳的變化，保留地層變化時的趨勢性特征。然后，基于小波變換模極大值計算測井曲線上的極大值點，用Lipschitz指數(shù)（Lip）表征測井曲線上極大值點的奇異性以驗證各極大值點的真?zhèn)危瑢崿F(xiàn)測井曲線層位特征的初步提取。

②動態(tài)閾值計算。根據(jù)伽馬圖像上不同層序的特征與范圍，以熵和類間差為衡量標(biāo)準(zhǔn)計算所有層序的局部閾值。

③地層識別。依據(jù)各層序的局部閾值，計算各層序的輪廓點集合。將輪廓點集合作為輸入數(shù)據(jù)，通過LSTM分類模型篩選真假地層邊界輪廓。

④傾角計算。用非線性最小二乘法將各輪廓點集合擬合到正弦曲線，得到正弦曲線的振幅、初相位、角速度等關(guān)鍵信息。由傾角計算公式計算得到地層的相對傾角，根據(jù)井斜參數(shù)得到地層真傾角及方位角。

2 關(guān)鍵技術(shù)

2.1 基于小波變換的層位特征提取與動態(tài)閾值計算

小波變換作為一種處理信號的重要手段，在層序劃分[15]、圖像增強[18]等方面有廣泛的應(yīng)用。在信號奇異值計算中，通常使用卷積運算的連續(xù)小波變換。

2.1.1 基于小波變換的奇異檢測原理

設(shè)曲線函數(shù)為 f（t），Wf（s,t）為 f（t）在對應(yīng)母波函數(shù) ψs（t）上的分解。s為尺度因子，在實際應(yīng)用中，通常取值為 s=2j（j∈Z）[19]，且尺度的選擇必須與具體信號結(jié)合，通常選定j=4或j=5[20]，Wf（s,t）即為f（t）的卷積型二進制小波變換。

假設(shè) t0是曲線函數(shù) f（t）在尺度 s0下的局部極值點，則有：

如果t0某一鄰域內(nèi)的任意點t都滿足如（2）式所示的條件，則稱 t0為小波變換的模極大值點，稱Wf（s0,t0）為小波變換的模極大值。

基于小波變換模極大值計算曲線函數(shù)f（t）的奇異值，通過Lip指數(shù)表征f（t）在某點上的奇異性。若曲線函數(shù)f（t）在某點處連續(xù)可微，導(dǎo)數(shù)有界且不連續(xù)，則該點的Lip指數(shù)為1，該點無奇異性；若f（t）在該點處不連續(xù)但有界，則該點的Lip指數(shù)為0，該點有奇異性。f（t）在某點處的尖銳程度隨Lip指數(shù)的減小而增大，通過計算Lip指數(shù)可以得知f（t）上某點的奇異程度，輔助校驗曲線極大值點的突變特性。

2.1.2 層位特征提取與閾值計算

基于小波變換的層位特征提取與動態(tài)閾值計算方法流程為：首先對隨鉆伽馬測井曲線進行濾波處理，然后依據(jù)小波變換模極大值及 Lip指數(shù)對伽馬曲線進行層序劃分，最后根據(jù)動態(tài)劃分結(jié)果計算局部閾值。

2.1.2.1 測井曲線噪聲影響及濾波處理

對伽馬曲線的平滑要突出不同層巖性的強烈變化，削弱同層內(nèi)巖性的細微變化，避免產(chǎn)生由伽馬測井儀器的放射性測量統(tǒng)計漲落、數(shù)據(jù)傳輸時的外界干擾、同層巖性伽馬波動等因素[17]導(dǎo)致的曲線小直徑波峰現(xiàn)象，使供層位提取使用的測井曲線突變位置更加清晰。

首先采用Lowess濾波方法對方位伽馬曲線進行過濾，其基本思路參見文獻[21]。以隨鉆伽馬平均曲線為例，各濾波方法對測井曲線平滑處理的結(jié)果如圖 1所示。圖 1表明，與直線型形態(tài)學(xué)濾波、半圓型形態(tài)學(xué)濾波相比，Lowess濾波不會丟失重要的曲線突變特征；與Savitzky-Golay濾波相比，Lowess濾波保留的重要突變特征的位置更吻合原曲線形態(tài)。

圖1 濾波效果對比

2.1.2.2 測井曲線層序劃分

基于Lowess濾波后的平均伽馬曲線，采用小波變換模極大值計算平均伽馬曲線的奇異點，使用Lip指數(shù)檢測曲線上點的奇異性。平均伽馬曲線上的突變點易被檢測為奇異點，而突變的位置往往代表可能存在的地層變化，如圖 2所示。因此，可以基于小波變換模極大值定位地層可能存在變化的深度，作為計算局部閾值和地層邊界輪廓點的依據(jù)。

圖2 曲線與地層變化關(guān)系示意圖

結(jié)合測井圖像特點及地層分布的規(guī)律，對于分層點 a0，若在[a0-C，a0+C]內(nèi)無其他分層點，則 a0的局部閾值計算范圍為[a0-C，a0+C]；若在[a0-C，a0]內(nèi)有其他分層點（設(shè)距離最近的分層點為a1），在[a0，a0+C]內(nèi)無其他分層點，則 a0的局部閾值計算范圍為[(a0-a1)ρ，a0+C]；若在[a0-C，a0]內(nèi)無其他分層點，在[a0，a0+C]內(nèi)有其他分層點（設(shè)距離最近的分層點為a2），則 a0的局部閾值計算范圍為[a0-C，(a2-a0)ρ]；若在[a0-C，a0]和[a0，a0+C]內(nèi)均有其他分層點（設(shè)兩個區(qū)間內(nèi)距離最近的分層點分別為 a1和 a2），則 a0的局部閾值計算范圍為[(a0-a1)ρ，(a2-a0)ρ]。其中，C=δ(D+2H)/λ，為常量，決定分層點a0在查找上下有無其他分層點時的最大范圍。δ和ρ分別為與上下無分層點、上下有分層點時的劃分范圍相關(guān)的參數(shù)，根據(jù)相對傾角在伽馬成像圖中展開的正弦曲線振幅規(guī)律，本文設(shè)δ和ρ的值分別為20和0.85。H為伽馬圖像的成像深度，依據(jù)文獻[22]研究成果計算后作為常數(shù)值處理，本文設(shè)為0.035 m。

根據(jù)分層點及其局部閾值計算范圍，將伽馬成像數(shù)據(jù)劃分為不同的層序，然后計算每一層序的局部閾值。

2.1.2.3 局部閾值計算

本文使用熵和類間差作為指標(biāo)衡量閾值[13]，該指標(biāo)綜合考慮了類內(nèi)聚和類間離散的特征，局部閾值的計算方法不再贅述。

2.2 基于LSTM的地層識別模型

通過層位特征提取計算的層序受曲線形態(tài)、濾波或算法效果等因素的影響勢必會導(dǎo)致提取到非真實地層，因此需要再對各層序進行篩選。LSTM模型的記憶門結(jié)構(gòu)使得該模型較適合解決具有序列特征的問題[23]，在多種場景下有廣泛的應(yīng)用[24-25]，而地層的輪廓數(shù)據(jù)具有序列的特征，因此采用LSTM模型設(shè)計地層識別分類器。

2.2.1 地層識別分類器設(shè)計

井筒與地層邊界相交，反映在隨鉆伽馬成像圖上為一條具有正弦特征的曲線，對輪廓點集合的篩選可以抽象為LSTM網(wǎng)絡(luò)的正弦特征二分類問題。根據(jù)計算的輪廓點集合有限點的特征及LSTM設(shè)計原則，構(gòu)建地層識別分類器模型，如圖3所示。

圖3 基于LSTM的地層識別分類器模型

訓(xùn)練集由兩部分組成：第 1部分為人工測井解釋并標(biāo)注的樣本數(shù)據(jù)；第 2部分為計算機根據(jù)地層邊界特征生成的樣本數(shù)據(jù)。規(guī)定正樣本是具有正弦特征的邊界點集合，負(fù)樣本是不具有正弦特征的點集合。訓(xùn)練集總共包含300組正樣本和300組負(fù)樣本，每組的數(shù)據(jù)點數(shù)量為 N（N∈[50，100]），N的范圍與伽馬成像預(yù)處理時的方位伽馬插值數(shù)量相關(guān)。從后續(xù)的分類結(jié)果可見，樣本數(shù)量設(shè)置為600組即可達到預(yù)期精度。設(shè)計正樣本訓(xùn)練數(shù)據(jù)時，應(yīng)結(jié)合正弦曲線表達式及地層邊緣的特征，在此基礎(chǔ)上產(chǎn)生的訓(xùn)練數(shù)據(jù)集示例如圖4所示。模型最終輸出結(jié)果為1，0標(biāo)簽，分別表示數(shù)據(jù)點集合正弦特征的有、無。

圖4 建模產(chǎn)生的訓(xùn)練數(shù)據(jù)集示例

LSTM 模型隱藏層神經(jīng)元的選擇對訓(xùn)練結(jié)果有較大的影響，通常依據(jù)經(jīng)驗公式確定神經(jīng)元個數(shù)[26]：

采用交叉熵函數(shù)作為損失函數(shù)計算分類器結(jié)果與樣本標(biāo)簽的誤差，通過誤差反向傳播進行訓(xùn)練，在交叉熵平穩(wěn)變化后，用 Adam梯度優(yōu)化算法繼續(xù)調(diào)整，與其他梯度優(yōu)化器相比，Adam優(yōu)化器在LSTM模型的實際訓(xùn)練中效果更優(yōu)[27]。

2.2.2 模型驗證與分析

結(jié)合（3）式設(shè)定LSTM模型中不同的隱藏層神經(jīng)元節(jié)點數(shù)，對比不同神經(jīng)元節(jié)點數(shù)的交叉熵值，選取最小誤差的節(jié)點數(shù)作為地層分類器的參數(shù)。本文不同神經(jīng)元節(jié)點數(shù)的交叉熵值如圖5所示。

圖5 不同神經(jīng)元節(jié)點數(shù)的交叉熵

由圖5可見，節(jié)點數(shù)在12附近時訓(xùn)練效果最優(yōu)，在4，6，8處效果較差。因此，本文選用12作為隱藏層的神經(jīng)元節(jié)點數(shù)，在該節(jié)點數(shù)下的地層識別分類器模型訓(xùn)練結(jié)果如圖6所示。該模型訓(xùn)練在epoch=80附近收斂，epoch=120附近時平均損失值下降至1.2×10-4左右，模型測試集的平均損失值約為2.3×10-4，說明該模型能夠正確地對有正弦特征的輪廓點數(shù)據(jù)和無正弦特征的輪廓點數(shù)據(jù)進行分類。

圖6 模型收斂曲線

在實際地層識別中，需要先計算各層序的局部閾值，得到各層序內(nèi)的輪廓點集合：假設(shè)某層序的閾值為th1，遍歷該層序中的所有數(shù)據(jù)點，若f(xi,yi)=th1，則標(biāo)記(xi,yi)為輪廓點；若 f(xi,yi)-th1與 f(xi+1,yi)-th1結(jié)果符號相反，則其中結(jié)果較小的數(shù)據(jù)點也將標(biāo)記為輪廓點。依此類推，計算所有層序的輪廓點集合并保存至邊界集合，邊界集合即由各輪廓點集合組成。運用地層識別分類器模型對邊界集合中的各輪廓點集合進行正弦性檢查，去除無正弦特征的輪廓點集合，使邊界集合能反映真實的地層邊界。

通過上述方法對某深度上的18組輪廓點數(shù)據(jù)進行篩選，結(jié)果如表1所示。結(jié)果表明，訓(xùn)練的基于LSTM的地層識別分類器具備篩選真、假地層的能力。

表1 輪廓點集合識別結(jié)果

2.3 傾角計算

通過擬合邊界集合中經(jīng)過正弦性檢查后的所有輪廓點集合，可獲得地層邊界輪廓對應(yīng)的正弦信息。擬合方法為：結(jié)合非線性最小二乘法，通過重復(fù)迭代尋找函數(shù)的局部最小值，從而使目標(biāo)函數(shù)取得最優(yōu)解。根據(jù)擬合后得到的正弦信息（見（4）式）及文獻[17]中的方法計算地層相對傾角和方位角，公式分別如（5）式、（6）式所示，具體推導(dǎo)過程不再贅述。

3 方法驗證與應(yīng)用

為驗證本文方法的準(zhǔn)確性和可行性，從方位伽馬數(shù)據(jù)解釋和隨鉆實時數(shù)據(jù)處理兩方面進行應(yīng)用及分析。

3.1 方位伽馬數(shù)據(jù)解釋

計算解釋的數(shù)據(jù)選用中國某油田X井井場8道方位伽馬數(shù)據(jù)（2 100～2 380 m井段），該井方位伽馬數(shù)據(jù)前期處理完善、伽馬成像結(jié)果清晰。通過對真實的井?dāng)?shù)據(jù)進行相關(guān)地層解釋，驗證本文方法的準(zhǔn)確性，包括算法的地層識別和傾角解釋能力。

本文方法的分層結(jié)果、擬合結(jié)果、傾角識別結(jié)果及與全局閾值法[13]的對比如圖 7所示?？梢钥闯?，本文方法有較好的分層及輪廓擬合效果，可以準(zhǔn)確地描述地層的變化、刻畫地層邊界輪廓曲線；與全局閾值法相比，本文方法能解釋出更多的地層傾角。局部井段的地層傾角本文方法解釋結(jié)果與人工解釋結(jié)果對比如表 2所示，可見本文方法解釋結(jié)果與人工解釋結(jié)果相當(dāng)。全部井段長度為280 m，測點22 400個，人工解釋層數(shù)為50，本文方法解釋層數(shù)為46，本文方法在較長井段、較多測點數(shù)的情況下，層數(shù)識別率為92%，傾角解釋準(zhǔn)確度達到96.6%，證明本文方法準(zhǔn)確性較好。

表2 采用本文方法及人工方法對真實方位伽馬數(shù)據(jù)進行解釋后得到的地層傾角解釋結(jié)果對比

圖7 本文算法與全局閾值法計算結(jié)果對比

3.2 隨鉆實時數(shù)據(jù)處理

Y井是中國東部油田某采油廠完善采油井網(wǎng)中的一口井，是為挖潛小層高部位剩余油、進一步提高儲量動用程度、改善開發(fā)效果而鉆探的一口開發(fā)水平井。地質(zhì)導(dǎo)向服務(wù)采用了中國自主研制的近鉆頭地質(zhì)導(dǎo)向系統(tǒng)+隨鉆方位伽馬+一體化隨鉆測量解釋系統(tǒng)組成的新型地質(zhì)導(dǎo)向系統(tǒng)。地質(zhì)導(dǎo)向技術(shù)服務(wù)井段從A靶前開始，從井深2 184 m鉆進到井深3 075 m完鉆，累計進尺891 m，井下工作82 h。系統(tǒng)不僅實現(xiàn)常規(guī)MWD（隨鉆測量）井斜、方位、工具面、溫度等軌跡姿態(tài)測量參數(shù)上傳，同時上傳的近鉆頭井斜、兩道伽馬（上、下）測量曲線也為實時地質(zhì)導(dǎo)向提供了及時準(zhǔn)確的決策依據(jù)。

使用該井?dāng)?shù)據(jù)驗證本文方法在實時鉆井生產(chǎn)中插值成像、地層邊界解釋效果。在該場景下本文方法的處理結(jié)果如圖 8所示?？梢钥闯?，本文方法基于數(shù)據(jù)插值的方法得到了實時伽馬數(shù)據(jù)的清晰成像；基于地層劃分及局部閾值計算的方法得到了符合實際情形的地層劃分；基于LSTM分類器模型得到了準(zhǔn)確的地層邊界，特別是圖8中井深2 440～2 450 m位置，本文方法通過分類器模型篩查出了多條無正弦特征的輪廓點集合，并結(jié)合快速正弦擬合的方法，將輪廓點擬合到了正確的位置。由表 3可見，本文方法在真實井實時數(shù)據(jù)中仍能夠保持可接受的誤差，與人工解釋的結(jié)果大致相同，該段耗時約為0.27 s，比人工解釋效率高。通過對圖8與表3所示結(jié)果的分析，證明了本文方法在實際生產(chǎn)中的效果。

表3 采用本文方法及人工方法對隨鉆實時數(shù)據(jù)進行解釋后得到的地層傾角解釋結(jié)果對比

圖8 井場隨鉆實時數(shù)據(jù)的解釋結(jié)果

以上計算結(jié)果與分析表明，本文方法在自動地層識別方面有較高的準(zhǔn)確率且計算速度較快，具備處理實時隨鉆測井?dāng)?shù)據(jù)的能力。實際應(yīng)用中，利用近鉆頭井斜測量確保軌跡準(zhǔn)確入靶并在靶體中延伸，借助近鉆頭上下方位伽馬及電阻率測量數(shù)據(jù)進行處理分析，保證了目的層鉆遇率。

4 結(jié)論

相較于全局閾值法和固定閾值法，基于小波變換的動態(tài)閾值法使地層界面解釋的準(zhǔn)確率顯著提升，且在地層巖性變化穩(wěn)定時減少了無用計算，并適用于地質(zhì)構(gòu)造變化較快的情況。

運用LSTM設(shè)計了基于隨鉆伽馬數(shù)據(jù)特征的地層邊界分類器模型，建立了地層邊界特征數(shù)據(jù)集，在訓(xùn)練過程中調(diào)整神經(jīng)元節(jié)點數(shù)及學(xué)習(xí)率等參數(shù)避免導(dǎo)致欠/過擬合的問題。分析對比發(fā)現(xiàn)，該分類器模型能準(zhǔn)確識別出具有正弦特征的輪廓，從而篩選出真實的地層邊界。

以人工智能為基礎(chǔ)的隨鉆測井解釋方法能夠為鉆井工程技術(shù)人員的決策提供輔助，為實時地質(zhì)導(dǎo)向應(yīng)用提供技術(shù)支持。

符號注釋：

a0，a1，a2——分層點，即層位提取步驟中的小波變換模極大值點位置；A——振幅；D——鉆孔直徑，m；f(·)——曲線函數(shù)；H——伽馬圖像的成像深度，m；i——曲線上點的序號；j——與尺度因子相關(guān)的系數(shù)；K——常數(shù)，K∈[0,10]；m，n——輸入層和輸出層的節(jié)點數(shù)；M——神經(jīng)元個數(shù)；N——每組樣本中數(shù)據(jù)點數(shù)量；s——尺度因子；th1——某層序的局部閾值；t——曲線 f(t)上的某點；t0——曲線函數(shù) f(t)在尺度 s0下的局部極值點；Wf(s,t)——f(t)在對應(yīng)母波函數(shù)ψs(t)上的分解；x——直角坐標(biāo)系x軸坐標(biāo)，m；xi，yi——曲線上某點的坐標(biāo)，m；y0——偏距，m；Z——整數(shù)集；α，β——地層相對傾角和方位角，（°）；δ——與上下無分層點時的劃分范圍相關(guān)的參數(shù)；λ——測井儀器測點的間距，m；ρ——與上下有分層點時的劃分范圍相關(guān)的參數(shù)；φ——初相；ψs(t)——母波函數(shù)；ω——角速度。