沈楠,段友祥,孫岐峰,李娜

(中國石油大學(華東)計算機科學與技術學院,山東青島266558)

0 引 言

伽馬成像測井解釋是對地層伽馬測井數(shù)據(jù)進行成像,通過對圖像的分析獲得地層構造和拾取地層傾角,進而更好地指導地質導向和地層評價工作。地層分界線在伽馬成像測井圖中呈現(xiàn)正弦形狀,拾取伽馬圖像中具有正弦構造的地層輪廓是測井解釋工作的重要內容[1-8]。在傳統(tǒng)的伽馬成像測井解釋中,對地層輪廓的拾取往往依賴于人工解釋或者解釋軟件輔助[9]。國內外代表性的測井軟件,如國內Forward[10]、Logview[11]、LEAD[12]系統(tǒng)以及國外的GeoFrame、Express、DPP系統(tǒng)[13]等,僅能提供測井圖像的繪制工具與基本的圖像處理工具,而在識別地質現(xiàn)象時,通常都是使用人工交互識別的方法。但是隨著測井數(shù)據(jù)量的增加和對解釋精度的要求不斷提高,人工圖像解釋工作量越來越大,難度增加,解釋效率很低。因此,中國學者開始對自動識別地質現(xiàn)象的方法展開研究,陸敬安等[14]利用霍夫變換提取出測井圖像中正弦曲線;張程恩[15]利用蟻群算法對測井圖像進行分割;張曉峰等[16]對測井圖像進行垂向和橫向二維小波分解,識別測井圖像中的裂縫信息;劉英明等[17]利用灰度直方圖實現(xiàn)測井圖像中正弦構造的拾取,并通過遺傳算法對其進行擬合。上述方法中,有些方法僅能從測井圖像中識別出地層輪廓的正弦構造,并沒有進行輪廓提取和傾角計算;有些方法雖然可以將地層輪廓從伽馬圖像中進行分割并提取出來,但是分割精度不夠,影響后續(xù)的計算結果。

本文結合圖像處理技術,提出了一種基于卷積神經(jīng)網(wǎng)絡UNet3+[18]的伽馬測井自動解釋方法,在像素級別上完成了地層輪廓的自動拾取,并對拾取輪廓進行細化處理。利用目標提取算法[19-21],確定地層變化位置以及每條輪廓的擬合區(qū)域,在區(qū)域內進行正弦擬合,最終計算出地層傾角和方位角。該方法消除了人為因素的不確定性,不需要具體的解釋公式以及選取大量的參數(shù),顯著提高了伽馬成像測井解釋的效率和準確率。

1 基于UNet3+的地層輪廓提取

1.1 基于UNet3+的地層輪廓識別模型

地層輪廓自動拾取是利用神經(jīng)網(wǎng)絡將地層輪廓從伽馬圖像中分割出來,完成此工作需要建立一個基于神經(jīng)網(wǎng)絡的圖像語義分割模型。根據(jù)圖像語義分割的處理方式,基于神經(jīng)網(wǎng)絡的圖像語義分割方法可以分為2大類:基于區(qū)域分類的圖像語義分割方法(Image Semantic Segmentation Based on the Regional Classification,ISSBRC)和基于像素分類的圖像語義分割方法(Image Semantic Segmentation Based on the Pixel Classification,ISSBPC)[22-24]。地層傾角是伽馬圖像解釋中需要獲得的重要參數(shù),而擬合地層輪廓的像素點是傾角計算過程中的關鍵一步。為了提高傾角計算結果的準確度和精度,本文采用基于像素點分類的圖像語義分割方法,對伽馬測井圖像實現(xiàn)像素級分割,提取出地層輪廓像素點。

圖1 UNet3+基本模型結構

本文采用UNet3+的基本結構,即編碼器-解碼器結構,網(wǎng)絡左部分為編碼器部分,由卷積層和池化層構成,池化層選用最大池化運算降低網(wǎng)絡的規(guī)模,其作用是利用卷積、最大池化運算對圖像進行下采樣,將圖像特征進行提取,生成不同尺度的特征圖;右部分為解碼器部分,由卷積層、池化層和反卷積層構成,其作用是采用全尺度的跳躍連接并利用卷積、最大池化、反卷積運算將不同尺寸特征圖進行融合,最后利用濾波器對其進行卷積,輸出融合后的特征圖,從而逐步將特征圖恢復至原始的尺寸,實現(xiàn)對像素點的分類。深度監(jiān)督分布在每一個解碼器中,將每個解碼器的輸出與輸出標準進行比較計算損失,從全尺度的聚合特征圖中學習層次表達。

在本文基于UNet3+的地層輪廓拾取模型中,每個卷積層后采用ReLU激活函數(shù),在輸入輸出之間構建非線性映射關系,其計算如式(1)所示。ReLU函數(shù)是卷積神經(jīng)網(wǎng)絡中的主流激活函數(shù),相對于Sigmoid函數(shù)和Tanh函數(shù),計算比較簡單,只需要一個合理的閾值就可以輸出激勵值。而且,ReLU函數(shù)具有良好的可導性,可以克服梯度消失的問題,其導數(shù)如式(2)所示。

ReLU(x)=max(0,x)

(1)

(2)

式中,x為模型輸入。

(3)

在右半邊的解碼器結構中,每層解碼器加入深度監(jiān)督,通過將每個解碼塊的輸出與標準輸出進行比較,計算網(wǎng)絡損失,對邊界進行精準分割,損失函數(shù)的計算公式為

ms-ssim=

(4)

1.2 地層輪廓的提取

使用UNet3+對伽馬圖像進行分割,生成地層輪廓概率圖。地層輪廓概率圖中顯示每個點屬于或不屬于地層輪廓的概率,利用Argmax函數(shù)處理概率值,將二維張量轉化為{0,1}的單個輸出,表示屬于/不屬于地層輪廓。則該二值圖像就是提取出的地層輪廓圖像。

1.2.1數(shù)據(jù)集

數(shù)據(jù)集采用中國某油田X井伽馬成像測井圖,并由專業(yè)人員標注地層。按照地層標注位置裁剪伽馬圖像,生成子塊圖像數(shù)據(jù)集,子塊圖像中含有一個或多個地層輪廓。鉆孔與地層分界面相交,在伽馬圖像中呈現(xiàn)正弦構造,鉆孔穿過地層的方式(由上至下或由下至上)的不同使正弦構造呈現(xiàn)出不同的凹凸特征。由于UNet3+基于像素對伽馬圖像進行分割,所以在訓練模型時,為了保證分割精度,需考慮不同的測井軟件在生成伽馬圖像時對色標的選擇。基于此,為了提高模型的特征提取能力與泛化能力,通過翻轉、拉伸、填充對數(shù)據(jù)集進行擴充,最終得到400張256×256像素的原始伽馬樣本圖像和400張256×256像素的專家標注的輪廓標簽圖像。選用90%的數(shù)據(jù)集,即360張原始伽馬樣本圖像和360張輪廓標簽圖像作為訓練數(shù)據(jù)集;選用10%的數(shù)據(jù)集,即40張原始伽馬樣本圖像和40張輪廓標簽圖像作為測試數(shù)據(jù)集,部分數(shù)據(jù)集見圖2。

圖2 部分樣本數(shù)據(jù)集示例圖

1.2.2實驗

為了驗證本文基于UNet3+的地層輪廓識別模型對伽馬圖像的分割性能,采用相同的伽馬訓練集對FCN中圖像分割效果最好的FCN-8s、UNet、UNet++網(wǎng)絡模型和本文模型進行訓練,并在測試集和實際伽馬圖像上,從可視化效果和定量評價指標2方面對比上述4種網(wǎng)絡模型的輪廓識別結果。

將360組訓練數(shù)據(jù),分別輸入FCN-8s、UNet、UNet++、UNet3+中進行訓練與學習,得到預測模型。實驗采用Keras和Tensorflow深度學習框架搭建模型,訓練過程中,FCN-8s網(wǎng)絡學習率初始值設置為0.001,批次為2,運行迭代次數(shù)為500,損失函數(shù)采用交叉熵損失函數(shù)為

(5)

UNet、UNet++、UNet3+使用相同的訓練參數(shù),學習率設置為0.000 1,批次為4,迭代次數(shù)為400。UNet、UNet++網(wǎng)絡的損失函數(shù)采用Loss函數(shù)[見式(5)],UNet3+的損失公式采用式(4)。

為考量模型的有效性,即檢測輪廓的靈敏度與對輪廓的分割精度,從2個方面進行分析評價:①對提取的地層輪廓的可視化效果進行評價,即評價地層輪廓的準確率、連續(xù)性和平滑性;②定量指標評價,即準確率(Precision)、召回率(Recall)[29]、F1[30]、IOU[31]的評價,F1是對準確率和召回率綜合考量的指標,IOU是真實值與預測值的交集除以像素的真實值和預測值的并集。因為模型是對伽馬圖像實現(xiàn)像素級的分割,所以評價指標均是基于像素點數(shù)量進行計算。

在40組驗證集上對各模型進行輪廓提取實驗,視覺結果如圖3,定量指標結果見表1。由圖3可見,4種網(wǎng)絡模型都能從伽馬圖像中分割出地層輪廓。由表1可見,UNet++、UNet3+的Precision、IOU值高于FCN-8s、UNet,說明UNet++、UNet3+的分割精度高于FCN-8s、UNet模型的分割精度。

圖3 部分測試集分割結果

表1 測試集-模型定量評價結果

利用各模型對實際伽馬數(shù)據(jù)圖像進行輪廓提取實驗,即用訓練好的4種網(wǎng)絡模型FCN-8s、UNet、UNet++、UNet3+對實際伽馬圖像進行分割,得到地層輪廓,并在伽馬圖像中用綠色標簽顯示,各模型提取輪廓的視覺結果見圖4,定量評價結果見表2。

圖4 輪廓分割效果對比圖

由圖4可以看出,FCN-8s網(wǎng)絡僅能檢測出部分地層輪廓;UNet網(wǎng)絡雖然能檢測出地層輪廓,但是在一些位置上出現(xiàn)了誤檢,而且檢測出的地層輪廓比較粗糙且不完整;UNet++網(wǎng)絡和UNet3+網(wǎng)絡都可以全面且正確地檢測地層輪廓,UNet3+網(wǎng)絡對細節(jié)更加敏感,分割精度最好。由表2可見,相比其他模型,UNet3+模型在4個評價指標上都達到了最好,能夠更精準地分割地層輪廓。

表2 實際伽馬圖像-模型定量評價結果

2 地層輪廓細化

2.1 輪廓細化的基本方法

2.1.1基于非極大值抑制的輪廓細化

為了使地層輪廓更好地呈現(xiàn)出正弦構造、提高傾角計算的精度,需要對地層輪廓進行細化操作。將伽馬圖像輸入到UNet3+中,得到地層輪廓概率圖,該圖顯示每個像素點屬于地層輪廓的概率,利用其對地層輪廓進行細化[32]。計算輪廓概率圖中每一點概率的梯度方向,在梯度方向上使用非極大值抑制(Non-Maximum Suppression,NMS)方法得到概率極大值點的集合[33-34]。

2.1.2地層輪廓概率梯度方向的計算

使用概率圖的二階導數(shù)求出輪廓中每個像素點的概率方向,概率圖在X、Y方向上的二階導數(shù)如式(6)、式(7)所示。

Gx=P(x-1,y)+P(x+1,y)-2P(x,y)

(6)

Gy=P(x,y-1)+P(x,y+1)-2P(x,y)

(7)

式中,Gx為概率圖在X方向上的二階導數(shù);Gy為概率圖在Y方向上的二階導數(shù);(x,y)為圖中像素點的坐標;P(x,y)為點(x,y)的概率。則點(x,y)的梯度方向的計算公式見式(8)。

θ=arctan(Gy/Gx)

(8)

式中,θ為點(x,y)的概率梯度方向,(°)。

2.1.3概率極大值像素點的獲取

對于任意一個像素,比較該像素點和其梯度方向正負方向像素點的地層輪廓概率梯度強度。若該像素點地層輪廓概率值最大,則認為該像素點是邊緣點;否則認為該像素點屬于后景,則刪除該像素點。最終得到地層輪廓概率極大值點的像素點的集合,進而得到細化后的地層輪廓。

若某像素點的地層輪廓概率值為0.85,地層輪廓概率的梯度方向為45°,梯度正方向的地層輪廓概率值為0.63,負方向的地層輪廓概率值為0.74。該點的地層輪廓概率值為0.85,大于該點梯度正方向的地層輪廓概率值0.63和梯度負方向的地層輪廓概率值0.74,該像素點是其地層輪廓概率梯度強度的極大值點,保留該點。

2.2 輪廓細化的具體步驟

以圖5(b)中紅色框內的部分輪廓[放大后見圖5(e)]為例,對細化步驟進行描述。

首先,利用地層輪廓概率圖的二階導數(shù),求出輪廓中每個像素點的地層輪廓概率方向θ。其次,由于1個像素點周圍有8個像素點,將每個像素點的地層輪廓概率方向θ按相近程度離散為8個方向,按上下左右和45°方向,劃分為8個區(qū)域。例如,某像素點的梯度方向θ為56°,在[45°,90°]區(qū)域內,計算θ和45°、90°的差值,按照相近原則,該點地層輪廓概率梯度方向θ為56°,其被近似為45°。圖5(e)對應的地層輪廓概率圖見圖5(f),其中每個像素點的數(shù)值代表了該點屬于地層輪廓的概率,紅色箭頭表明該像素點的地層輪廓概率梯度方向。

圖5 輪廓細化示意圖

最后,對求出梯度方向的地層輪廓概率圖采用非極大值抑制方法,保留每個像素點上地層輪廓概率梯度強度的極大值,刪除其他值,得到地層輪廓概率極大值點的像素點的集合。圖5(g)顯示對地層輪廓概率圖進行非極大值抑制的結果,對其進行可視化,可以得到細化后的地層輪廓[見圖5(h)]。

3 地層傾角計算

對地層輪廓提取細化后,根據(jù)地層輪廓的位置信息計算擬合區(qū)域,通過擬合區(qū)域內地層輪廓點,獲得地層邊界輪廓對應的正弦信息并計算地層傾角。

3.1 輪廓擬合區(qū)域劃分

伽馬圖像的分割結果通常不只有1個地層輪廓,而地層傾角是對1個目標輪廓進行擬合得到,所以需要將伽馬圖像中已被識別的m個目標輪廓分別提取出來,并計算其對應的擬合區(qū)域。圖6為擬合區(qū)域示例圖,其中圖6(a)顯示了正確的劃分方法,含有2個地層輪廓,正確劃分為2個不重疊、不交叉且完全包含各自輪廓點的紅色、藍色矩形框區(qū)域。圖6(b)、(c)、(d)為錯誤劃分方法,圖6(b)中一個正弦擬合區(qū)域包含多條地層輪廓;圖6(c)中不同地層輪廓的擬合區(qū)域有交叉現(xiàn)象;圖6(d)中紅色矩形框區(qū)域不能完全包含一個輪廓。

圖6 擬合區(qū)域示例圖

為解決以上問題,實現(xiàn)輪廓擬合區(qū)域的精準劃分,采用Selective Search算法,將m個目標地層輪廓分別固定在m個矩形目標候選框中,之后的傾角計算在矩形目標候選框中進行。Selective Search算法[35]是將初始的分割區(qū)域作為輸入,使用相似度算法計算這些分割區(qū)域之間的相似性,相似性標準主要有顏色、紋理、大小等。最后根據(jù)區(qū)域相似度對子區(qū)域不斷地進行迭代合并,每次迭代過程中對這些合并的分割區(qū)域做外切矩形。

輪廓擬合區(qū)域劃分的具體步驟:①劃分出分割區(qū)域,將目標輪廓從背景圖中分離出來,從原像素圖中分離到與原圖一樣大小的黑色背景圖中;②采用Selective Search算法提取目標,產(chǎn)生目標候選框;③篩選目標候選框,經(jīng)過UNet3+分割后的地層輪廓圖像仍含有少量雜波,通過計算候選框面積并根據(jù)預先設置的閾值,去除面積較小的雜波,保留真實的地層輪廓信息。

3.2 傾角計算方法

首先在目標候選框內使用非線性最小二乘法的信賴區(qū)域法(Levenberg-Marquardt,LM)[36]對地層輪廓點進行擬合,獲得地層邊界輪廓對應的正弦信息;其次,根據(jù)擬合的正弦曲線,求出地層傾角。本文基于目標候選框的位置信息,提出了2種方法確定出地層輪廓點,從而利用LM算法對地層輪廓點進行正弦擬合。

方法一:基于UNet3+地層輪廓識別模型的分割結果直接確定地層輪廓點。對目標候選框內UNet3+分割伽馬圖像生成的目標輪廓點進行擬合,得出正弦曲線。

方法二:基于目標候選框的局部閾值法。①采用Selective Search算法處理黑色背景圖下的地層輪廓,生成目標候選框[見圖7(a)紅色矩形框];②根據(jù)目標候選框的位置在伽馬圖像中找到對應的擬合區(qū)域[見圖7(b)紅色矩形框];③在伽馬圖像的紅色矩形框內使用局部閾值法,求出地層邊界輪廓點,在圖7(c)中用綠色點表示;④對綠色地層輪廓點進行正弦擬合,正弦曲線見圖7(d)。

圖7 局部閾值法示意圖

通過上述2種方法得出地層邊界輪廓對應的正弦曲線,正弦曲線的標準公式為式(9);根據(jù)擬合出的正弦曲線,使用傾角計算公式計算出地層傾角和方位角,計算公式為式(10)、式(11)。

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

(9)

(10)

(11)

式中,φ為初相位,(°);A、ω、y0為常數(shù)(ω≠0);α為相對傾角,(°);D為鉆孔直徑,cm,通常假設其等于鉆頭尺寸,或者可以通過現(xiàn)場鉆孔卡尺測量確定;d為圖像的深度,cm;h為正弦曲線的高度,cm;β為方位角,(°)。

4 實驗結果對比與分析

為驗證本文方法在伽馬成像測井解釋中的有效性,采用中國某油田X井8道方位伽馬實際測井數(shù)據(jù)(2 270～2 350 m)進行實驗,并將本文方法與全局閾值法[8]、基于小波變換模極大值的局部閾值法[37]、專家解釋的結果進行對比。

本文利用UNet3+對地層輪廓進行拾取,由于模型輸入圖像的尺寸會影響地層輪廓的計算效率,最終將原始圖像(217×8 300)劃分為217×512的子塊圖像。將各子塊圖輸入到UNet3+中,得到子塊區(qū)域像素分割結果的概率值。

為了保證子塊圖像邊緣像素計算結果的準確率,相鄰2個子塊圖像間設有217×3的重疊區(qū)域,對于這些重疊區(qū)域,采用分割結果概率平均的策略對子塊進行拼接合成,最終得到大尺寸圖像的地層輪廓分割結果。拾取的輪廓呈現(xiàn)橫向趨勢,在計算輪廓的矩形擬合區(qū)域后,根據(jù)正弦曲線的構造特征以及原始圖像的大小,篩除寬度為1個像素或者面積小于100個像素的矩形擬合區(qū)域,并將篩選過后的矩形擬合區(qū)域寬度統(tǒng)一為217。解釋結果如圖8(a)所示。首先對伽馬圖像進行輪廓提取與細化,通過Selective Search算法計算出地層輪廓目標候選框,然后直接擬合目標候選框內UNet3+分割伽馬圖像生成的目標輪廓點集合,生成地層輪廓對應的正弦曲線。

圖8 本文方法測井解釋效果圖

圖8(b)為基于目標候選框的局部閾值法的伽馬測井解釋效果圖,與圖8(a)所用方法不同之處在于:在得出地層輪廓目標候選框后,根據(jù)目標候選框在伽馬圖像中的位置信息計算出伽馬圖像中的擬合區(qū)域,利用局部閾值法計算擬合區(qū)域內伽馬圖像的局部閾值,根據(jù)該局部閾值分割目標候選框內伽馬圖像得出地層輪廓點集合,最后對地層輪廓點集合進行擬合。

由此可見,本文方法可準確分割伽馬圖像,拾取地層輪廓;能夠有效地對地層輪廓進行細化,準確計算地層輪廓的擬合區(qū)域,確定地層變化位置。

為進一步驗證本文方法在實際伽馬圖像中的有效性,采用全局閾值法和基于小波變換模極大值的局部閾值法處理相同的實際伽馬圖像。全局閾值法是利用最大類間差法和最大熵算法計算整個伽馬圖像的全局閾值,利用閾值分割技術找到地層輪廓的像素點,并對輪廓點進行正弦擬合,最后利用正弦曲線進行傾角計算。基于小波變換模極大值的局部閾值法是用小波變換模極大值的方法對整個伽馬圖像進行劃分,采用Lowess濾波處理伽馬曲線,確定地層變化位置,利用局部閾值法確定地層邊界的輪廓點集合,最后對輪廓點進行擬合,并通過擬合出的正弦曲線計算出地層傾角。

全局閾值法的計算結果如圖9(a)所示,因為該方法采用單一閾值,所以對圖像噪聲較敏感,對差異不明顯的圖像分割較差;而本文方法通過結合多尺度的低分辨率和高分辨率特征,自動學習合適的特征表示,對地層變化的靈敏度更強,可檢測全局閾值法檢測不到的地層位置,解決了圖像噪聲和分布不均勻的問題?；谛〔ㄗ儞Q模極大值的局部閾值法的計算結果如圖9(b)所示,圖9中藍色曲線為平均伽馬曲線,紅色曲線為使用Lowess濾波處理后的曲線,紅色橫線為地層識別結果。由圖9可以得出:本文方法與基于小波變換模極大值的局部閾值法均可以較好地找出地層變化位置,這是由于基于小波變換模極大值的局部閾值法采用了自適應閾值來提高分割精度,但是其分割結果過于依賴濾波器的選擇,且需要訓練分類模型對結果進行進一步地篩選。通過對比不同局部閾值法的解釋效果[見圖8(b)、圖9(b)],可以得出:本文提出的基于目標候選框的局部閾值法可以更加精確地求解地層邊界在圖像中的分割范圍,識別地層變化位置。

圖9 伽馬測井解釋效果對比圖

表3展示了本文基于UNet3+分割結果的傾角計算結果和基于目標候選框的局部閾值法的傾角計算結果,并與全局閾值法、基于小波變換模極大值的局部閾值法的傾角計算結果、專家解釋結果對比。

由表3可以得出,本文基于UNet3+分割結果的傾角計算結果和基于目標候選框的局部閾值法的傾角計算結果誤差最小,與人工解釋結果相當,可以準確且有效地對伽馬圖像進行解釋,同時大大減少了人工解釋的工作量。

表3 傾角計算結果對比表

5 結 論

(1)針對傳統(tǒng)的地層輪廓識別存在工作量大、識別精度不高、效率低等問題,本文通過建立基于UNet3+的地層輪廓識別模型在像素級別上實現(xiàn)了對伽馬圖像的分割,完成地層輪廓的自動識別。通過從視覺效果和定量指標方面與FCN-8s、UNet、UNet++的分割結果進行對比,最終得出本文基于UNet3+的地層輪廓識別模型能夠全面且正確地檢測出地層輪廓,對細節(jié)更加敏感,分割精度最好。

(2)基于地層輪廓模型的分割結果和地層輪廓擬合區(qū)域,本文提出了2種不同的地層輪廓點提取方法:①直接采用UNet3+的分割結果作為地層輪廓點集合;②利用局部閾值法計算地層輪廓擬合區(qū)域內的伽馬圖像的最優(yōu)局部閾值并對伽馬圖像進行分割得出地層輪廓點集合。在得出輪廓點集合后,結合LM擬合算法和傾角計算公式得出地層傾角。

(3)將本文方法應用于實際伽馬圖像并與全局閾值法、基于小波變換模極大值的局部閾值法進行對比,結果表明:本文方法對地層變化位置更敏感、識別精確度更高、傾角計算更準確,與專家解釋結果質量相當;本文方法能夠顯著提高伽馬圖像解釋的效率,為深度學習等新技術在測井解釋領域的應用進行了有益的探索。