段友祥,閆亞男,孫歧峰,李洪強,2

(1.中國石油大學(華東)計算機與通信工程學院,山東 青島 266580;2.中國石化集團勝利鉆井工藝研究院,山東 東營 257000)

0 引 言

隨鉆方位伽馬測井以常規(guī)自然伽馬測井為基礎(chǔ),利用儀器不同方向上的多個伽馬探測器進行測量,使伽馬數(shù)據(jù)具有方位特性,可以根據(jù)方位特性反演得到地層構(gòu)造情況,從而獲得儲層的邊界、儲層厚度和地層傾角等相關(guān)信息。在地質(zhì)導向中可根據(jù)隨鉆方位伽馬數(shù)據(jù)指導鉆頭準確入層,及時獲得地層傾角,幫助分析地層結(jié)構(gòu),使井眼軌跡更準確地保持在目的層中,提高油氣采收率。

地層傾角計算的2個必要參數(shù)分別為同一地層上不同2點之間的高度差和水平距離。目前,計算地層傾角主要有2種方案。①利用鄰井信息,根據(jù)不同井鉆遇同一地層時,2個鉆點參數(shù)計算地層傾角。楊錦舟[1]應(yīng)用鄰井自然伽馬建立先驗?zāi)Ｐ徒Y(jié)合人工經(jīng)驗進行了地質(zhì)導向鉆井;王謙等[2]根據(jù)鄰井數(shù)據(jù)建立模型完成地層傾角計算和儲層預測。但是,應(yīng)用鄰井信息計算傾角易受資源限制,鄰井資料完整性影響著地層信息的準確性。②利用隨鉆方位伽馬數(shù)據(jù),根據(jù)方位數(shù)據(jù)計算地層傾角。侯爽[3]通過判斷數(shù)據(jù)增量識別地層,但是不同的增量閾值會造成較大的效果差異;袁超[4]根據(jù)自然伽馬圖像進行了地層傾角及厚度計算,但需人工拾取計算的相關(guān)參數(shù);McKinny K[5],Jiaxin Wang等[6]根據(jù)穩(wěn)定區(qū)域均值識別邊界并計算地層傾角,但是在實際鉆井中數(shù)據(jù)很難穩(wěn)定不變;Uzoh等[7]模擬不同地層條件下的方位伽馬圖像和曲線模型為人工識別服務(wù);唐海全、肖紅兵等[8]利用方位伽馬曲線走勢判斷鉆遇地層,但是無法計算地層傾角;Tyurin E等[9]在三維空間上計算地層真傾角,但仍需人工拾取地層特征。利用隨鉆方位伽馬數(shù)據(jù)是計算地層傾角較為準確快速的方式,但是如何提取地層特征是目前難點之一。識別地層可以采用人工拾取、人機交互、計算機自動識別等方式。人工方式受經(jīng)驗影響較大,且速度慢;效率低;人機交互方式是人工拾取地層關(guān)鍵點,再由計算機計算,這種方式效率較高,但仍難以擺脫經(jīng)驗影響;計算機自動識別可以避免人為造成的影響和誤差,而且效率更高,有利于快速作出正確的綜合決策,提高開發(fā)效益,是目前的主要研究方向。

針對上述問題,本文以隨鉆方位伽馬數(shù)據(jù)為基礎(chǔ),結(jié)合計算機圖像處理等相關(guān)技術(shù),對隨鉆方位伽馬圖像進行圖像分割;通過地層邊界曲線獲取地層信息并完成地層傾角計算,以達到指導地質(zhì)導向鉆井的目的。該方法不僅實現(xiàn)了隨鉆方位伽馬圖像中的地層自動識別與傾角計算,而且滿足了地質(zhì)導向中的實際應(yīng)用需求。

1 算法基本思想及流程

本文主要研究地層傾角自動識別,對隨鉆方位伽馬數(shù)據(jù)進行歸一化,并對歸一化結(jié)果構(gòu)造色階直方圖,綜合考慮類內(nèi)聚性和類間離散度,引入綜合指標選取分割閾值;根據(jù)所選閾值對所有數(shù)據(jù)點進行閾值標記,自適應(yīng)伽馬數(shù)據(jù)測距,確定跟蹤范圍,按自定義跟蹤順序進行輪廓跟蹤,獲取地層邊界;將哈夫變換正弦擬合與最小二乘擬合相結(jié)合,確定地層正弦公式;根據(jù)地層正弦參數(shù)計算地層傾角。算法流程見圖1。

圖1 算法流程圖

2 伽馬圖像地層提取

2.1 綜合指標閾值圖像分割

識別地層就是要對伽馬圖像進行圖像分割,即將伽馬數(shù)據(jù)像素點分為不同的類,既要提高類內(nèi)的相似度又要降低類間的相似度。本文首先對給定的隨鉆方位伽馬數(shù)據(jù)進行歸一化,再根據(jù)成像結(jié)果進行閾值分割。傳統(tǒng)閾值分割采用圖像灰度構(gòu)造直方圖,本文將伽馬值映射為0～127級色階,再轉(zhuǎn)換為圖像,但是在伽馬數(shù)據(jù)→色階→RGB值→灰度值的轉(zhuǎn)化過程中,會丟失灰度值與伽馬數(shù)據(jù)的正相關(guān)關(guān)系,所以本文將伽馬數(shù)據(jù)映射后的色階值作為數(shù)據(jù)輸入,構(gòu)造色階直方圖。

本文引用汪啟偉[13]提出的綜合指標,綜合考慮衡量類內(nèi)內(nèi)聚性的熵和衡量類間離散度的類間差,選取閾值。熵值由式(1)計算,類間差由式(2)計算,綜合指標由式(3)計算得到。

Hent(T)=

(1)

Hotsu(T)=pb(T)po(T)[mb(T)-mo(T)]2

(2)

(3)

式中,T為當前閾值;pb(T)、po(T)為2個類別所占總樣本的比例;mb(T)、mo(T)為2個類別的類均值;went和wotsu分別為熵和類間差在計算綜合指標時的權(quán)重,計算所有可選閾值下的綜合指標,選取綜合指標最大的值作為圖像分割閾值。

2.2 自適應(yīng)擴充鄰域輪廓跟蹤

圖像分割閾值點即為表現(xiàn)伽馬射線活動對比度平面上的點,確定閾值后,進行輪廓跟蹤。輪廓跟蹤通過順序找出邊緣點獲取地層邊界。具體過程如下

(1)列優(yōu)先遍歷圖像中的所有數(shù)據(jù)點,對數(shù)據(jù)點(x,y),如果f(x,y)=value(value為閾值),則標記該點為閾值點;否則如果滿足f(x,y)與f(x+1,y)位于閾值點兩側(cè),則選取(x,y)與(x+1,y)中與閾值差距小的點標記為閾值點。

(2)建立二維鏈表BoundaryList作為邊界鏈表。

(3)按照從上到下,從左到右的順序找到一個未標記邊的閾值點A,建立一個一維鏈表BList。

(4)若A不是噪聲點或邊界端點,在A點周圍一定有一個閾值點,一般來講應(yīng)該按照左、左上、上、右上、右、右下、下、左下的順序跟蹤下一個閾值點。但是以閾值點作為邊界有可能出現(xiàn)如圖2所示情況,圖2中紅色點為閾值點,P1與P2同屬一邊界,但按照上述跟蹤順序,無法跟蹤到P2。當出現(xiàn)P2、P3、P4的情況,從P2開始,按照上述順序,會先跟蹤到P4,導致丟失P3,但是針對實際情況,期望跟蹤點依次為P2、P3、P4。根據(jù)隨鉆方位伽馬數(shù)據(jù)的特性,只要邊界連續(xù),那么在圖3中灰色區(qū)域內(nèi)一定有閾值點。2點在縱向上相差行數(shù)由測距決定,可設(shè)定最小值,當測距較小時,2點縱向相差行數(shù)可能較大,根據(jù)實際情況,選取L內(nèi)含有的數(shù)據(jù)行數(shù)為檢測范圍。

圖2 輪廓跟蹤圖示 圖3 伽馬數(shù)據(jù)輪廓區(qū)域示例

本文分析了伽馬數(shù)據(jù)中閾值點的分布屬性,創(chuàng)建了適用于方位伽馬數(shù)據(jù)的自適應(yīng)輪廓跟蹤算法,根據(jù)測距自動收縮或擴張跟蹤范圍。以先右后左、順時針、近距離優(yōu)先為原則建立了如圖4所示跟蹤順序。

圖4 自定義輪廓跟蹤順序

參數(shù)N由式(4)計算

N=4(L/step)+12

(4)

式中,Step為伽馬數(shù)據(jù)測距;L為單邊檢測范圍,本文中取0.2。按上述順序跟蹤下一個邊界點,標記該點所屬邊,將該點加入到邊界鏈表BList。

(5)以最新點為A點,重復步驟(4),直至跟蹤不到閾值點,將單邊鏈表BList加入到總邊界鏈表BoundaryList。

(6)跟蹤一條邊界結(jié)束,查找下一條邊界,重復步驟(3)到步驟(5)。

3 地層傾角計算

3.1 改進正弦擬合

經(jīng)過輪廓跟蹤獲取到多條地層邊界后,要對地層邊界進行曲線擬合,確定每個地層平面在三維空間中的位置。鉆孔與表現(xiàn)伽馬射線活動對比度的平面相交,其形態(tài)見圖5。在三維空間內(nèi)投影在地層平面上顯示為橢圓,在井壁的平面展開圖上,觀察到特征正弦響應(yīng),所以要將地層展開面上的地層邊界線擬合為正弦曲線。

圖5 地層與井筒相交示意圖

井壁展開面上的正弦曲線表達式見式(5)

f(x)=Asin (ωx+φ)+y0

(5)

式中,A為振幅;ω為角速度;φ為初相位;y0為偏距,振幅是傾角計算中的關(guān)鍵參數(shù)。如圖5所示,井壁展開面上是一個完整的正弦周期,即周期T是井壁展開面的寬度,角速度可由式(6)計算

(6)

在正弦公式中,只需要擬合3個參數(shù):振幅A,初相位φ,偏距y0。

最小二乘擬合、迭代擬合[10]、基于哈夫變換的擬合[11-12]是測井數(shù)據(jù)處理中常用的3種正弦擬合算法,它們各有自己的優(yōu)勢和不足。最小二乘算法計算精度高,但是當范圍過大時,計算量過大,計算速度慢;迭代擬合算法與最小二乘算法相比速度大幅提高,在單參數(shù)擬合時收斂效果較好,但在三參數(shù)擬合時難以保證收斂于最優(yōu)值;哈夫變換算法速度最快,適合大量點中選擇少數(shù)點擬合,是一個拾取和擬合相結(jié)合的過程,但是不適用于本文應(yīng)用場景。為此,本文在研究分析3種擬合算法優(yōu)缺點的基礎(chǔ)上,根據(jù)伽馬數(shù)據(jù)本身的特點,提出了一種適用于伽馬數(shù)據(jù)的綜合擬合算法。該算法過程如下:

(1)邊界點排序。對輪廓跟蹤中得到的二維鏈表BoundaryList,遍歷每1條一維鏈表BList,對鏈表中的數(shù)據(jù)以x為關(guān)鍵字排序。

(2)確定初始偏移量。遍歷鏈表BList中的點,對一個點(xa,ya)在BList中查詢滿足xb=xa+T/2的點(xb,yb),偏移y0根據(jù)式(7)計算

(7)

式(7)的原理是正弦本身的對稱性。如圖6所示,正弦上2點S1、S2橫坐標相差半周期,2點中點Q一定在直線y=y0上。對所有計算結(jié)果,建立投票機制選擇y0。

圖6 正弦基線計算圖示

(3)確定最小二乘參數(shù)范圍。遍歷鏈表BList中每一個點,得到y(tǒng)的最大值ymax、最小值ymin,以及y取最大值時對應(yīng)的橫坐標xmax,正弦公式中偏移量的遍歷范圍是[y0-(ymax-ymin)/6,y0+(ymax-ymin)/6],振幅的遍歷范圍是[0,(ymax-ymin)×2/3];由于正弦曲線本身的性質(zhì),初相位的遍歷范圍由式(8)確定

(8)

(4)計算殘差。對3個參數(shù)的可取值自由組合,對每種組合(yi,Aj,φk)計算殘差,殘差Errorijk由式(9)計算

(9)

(5)確定參數(shù)。比較所有殘差,殘差最小的三參數(shù)取值即為該條邊界的參數(shù)值。

(6)計算下一條邊界。取BoundaryList中的下一條邊界,重復步驟(2)到步驟(5),直至所有邊界計算完畢。

3.2 傾角計算

地層傾角是地質(zhì)導向中指導鉆頭行進的關(guān)鍵參數(shù),根據(jù)正弦公式即可計算地層傾角。如圖5所示,平面邊界與鉆孔相交的相對傾角可以用式(10)確定

(10)

式中,D為鉆孔直徑,計算中可認為等于鉆頭尺寸;H為正弦曲線在井壁展開面上的展布高度,等于正弦曲線振幅的2倍;根據(jù)Jiaxin Wang等[6]論文中的定義,DI為伽馬圖像的成像深度。Jiaxin Wang詳細研究了不同地層傾角、地層密度、鉆孔尺寸等環(huán)境對DI的影響,本文采用該思想,利用式(10),根據(jù)上文所得平面正弦參數(shù)計算地層傾角。

4 實驗與分析

在隨鉆測井中,采用泥漿脈沖遙測方式將隨鉆方位伽馬測井數(shù)據(jù)實時傳輸至地面,并對傳輸數(shù)據(jù)進行實時預處理與伽馬成像。在實際鉆井中,地質(zhì)情況復雜,地層構(gòu)造信息不精確。基于人工構(gòu)造的地層模型,對測井響應(yīng)進行正演模擬,可達到實測數(shù)據(jù)效果,而且可以提供準確的地層界面信息作對比參考。本文分析了前人研究的各種正演算法,最終采用邵才瑞[14]、Zhen Qin等[15]所用積分思想進行了隨鉆方位伽馬數(shù)據(jù)模擬。

以隨鉆方位伽馬模擬數(shù)據(jù)為基礎(chǔ),實驗環(huán)境采用Visual Studio 2013和Open Scene Graph,實現(xiàn)了隨鉆方位伽馬圖像的地層自動識別與傾角計算。采用滑動窗口分析伽馬圖像,以最新更新數(shù)據(jù)為滑動窗口結(jié)束位置,使用本文算法處理滑動窗口內(nèi)的圖像數(shù)據(jù)?；瑒哟翱诘拇伴L和滑動步長可根據(jù)實際情況設(shè)定,窗長設(shè)定過長,會導致數(shù)據(jù)量過多,降低分析效率;窗長設(shè)定過短,可能造成大傾角地層的傾角信息不完整,測得傾角存在誤差;滑動步長設(shè)定過長可能會導致發(fā)現(xiàn)地層較晚,錯過最佳控制鉆頭的時機;滑動步長設(shè)定過短會導致頻繁分析伽馬圖像,降低分析效率,嚴重可能導致延誤地層發(fā)現(xiàn)時機。設(shè)定更新井軌跡上一個采樣點八方位數(shù)據(jù)的時間為一個單位時間,伽馬數(shù)據(jù)的測距為0.1 m,八方位伽馬數(shù)據(jù)橫向內(nèi)插為64個點數(shù)據(jù),表1為該條件下不同窗長和滑動步長下的單位時間處理數(shù)據(jù)點數(shù)的比較。

表1 不同滑動窗口處理數(shù)據(jù)量比較

正演模型中地層模型為砂巖和泥巖,地層模型相關(guān)參數(shù)見表2。K每秒平均發(fā)射光子數(shù)3.4/gs;Th每秒平均發(fā)射光子數(shù)1×104/gs;U每秒平均發(fā)射光子數(shù)2.8×104/gs[16-17]。

表2 正演模型相關(guān)地層參數(shù)

模擬4組數(shù)據(jù)作為實驗數(shù)據(jù),圖7、圖8、圖9、圖10分別為4組數(shù)據(jù)實驗結(jié)果,表3為4組數(shù)據(jù)的識別結(jié)果與實際角度的對比及誤差。

表3 本文算法地層傾角識別準確性

圖7 第1組數(shù)據(jù)識別結(jié)果

圖8 第2組數(shù)據(jù)識別結(jié)果

圖9 第3組數(shù)據(jù)識別結(jié)果

圖10 第4組數(shù)據(jù)識別結(jié)果

實驗組別擬合點數(shù)角度/(°)擬合算法耗時/s殘差第1組7284.29迭代最小二乘擬合0.04213.582哈夫變換0.0040.1011傳統(tǒng)最小二乘擬合25.7940.0606本文算法14.9580.0606第2組6578.69迭代最小二乘擬合0.0112.9477哈夫變換0.0040.0238傳統(tǒng)最小二乘擬合5.8290.0101本文算法3.2810.0101第3組6563.43迭代最小二乘擬合0.0030.5174哈夫變換0.0030.0116傳統(tǒng)最小二乘擬合1.2260.0072本文算法0.6150.0072第4組6526.57迭代最小二乘擬合0.0010.3158哈夫變換0.0030.0167傳統(tǒng)最小二乘擬合0.1160.0027本文算法0.0620.0027

本文所用正弦擬合算法,與迭代擬合和哈夫變換擬合相比殘差更小;與傳統(tǒng)最小二乘算法相比,耗時更少,表4為4種算法耗時及殘差對比。

實驗證明,本文所用圖像分割、輪廓跟蹤、正弦擬合算法取得了較好效果,可自動識別隨鉆方位伽馬圖像中的地層并計算地層傾角,且精確度較高。

5 結(jié)論

(1)針對利用隨鉆方位伽馬數(shù)據(jù)獲取地層傾角的問題,研究分析了各種計算方式的局限性,根據(jù)相關(guān)研究的成果與原理啟發(fā),結(jié)合計算機圖像處理技術(shù),建立了利用隨鉆伽馬測井數(shù)據(jù)進行地層自動識別與傾角計算的方法,較好地解決了計算機自動進行地層傾角識別的問題。

(2)該方法采用綜合指標計算分割閾值;提出了根據(jù)隨鉆方位伽馬數(shù)據(jù)測距自定義跟蹤順序和擴充鄰域進行輪廓跟蹤;結(jié)合哈夫變換擬合思想改進擬合算法,縮小最小二乘計算范圍,提高擬合速度,擬合地層分界線并計算地層傾角。

(3)本文方法針對不同的地層角度和地層構(gòu)造均可有效識別地層并計算傾角,并且具有較高的正確率和運行效率。