劉旭躍， 周 巍， 張 兵， 黃 駿

(中國石化石油 物探技術(shù)研究院，南京 211103)

一種基于圖像學(xué)的地震層位自動追蹤方法

劉旭躍， 周 巍， 張 兵， 黃 駿

(中國石化石油 物探技術(shù)研究院，南京 211103)

構(gòu)造復(fù)雜地區(qū)地震數(shù)據(jù)受信噪比低、噪聲干擾等影響，層位追蹤結(jié)果精度不高，為解決該問題，提出了一種適用于低信噪比，同相軸不穩(wěn)定的基于圖像學(xué)的地震層位自動追蹤方法。該方法首先采用圖像學(xué)的全變分方法對地震數(shù)據(jù)進行保護邊緣的降噪處理，然后采用基于相似系數(shù)的相干圖像處理技術(shù)對地震剖面進行計算獲得相干剖面，同時提取出構(gòu)造的傾角，接著對同相軸進行自動追蹤，程序中采用鏈?zhǔn)浇Y(jié)構(gòu)存儲數(shù)據(jù)，便于檢通過實際數(shù)據(jù)測試，該方法能有效去除噪聲，增強同相軸連續(xù)性，提高自動追蹤效果。

全變分; 相干圖像處理; 同相軸傾角; 層位自動追蹤

0 引言

地震野外獲得的原始資料，經(jīng)過室內(nèi)處理后，得到可供解釋的地震剖面和其他成果圖件。對這些圖件進行分析研究，可以推斷出地震剖面上各反射層所相當(dāng)?shù)牡刭|(zhì)層位，從而了解地下地質(zhì)的情況。其中層位地識別和追蹤是地震層位解釋的重要環(huán)節(jié)，它對后續(xù)的處理和解釋工作有很大影響，是地震層位反射同相軸空間分布地直觀體現(xiàn)[1]。

目前常規(guī)應(yīng)用的層位追蹤方法主要包含：①人工拾取層位；②計算機自動追蹤[3]。人工拾取層位是在二維地震剖面上，利用波形相似性對地層連續(xù)反射同相軸進行手動追蹤得到層位線[4]。它的工作量大，效率低，準(zhǔn)確性差。隨著計算機和勘探技術(shù)的發(fā)展，層位自動追蹤技術(shù)得到了快速發(fā)展，它通過計算機計算，自動拾取出同相軸的位置。隨著自動追蹤效率的不斷提高，逐漸取代了人工拾取層位。國、內(nèi)外對層位追蹤的方法主要分為兩類:①基于波形特征的自動追蹤；②基于相關(guān)的自動追蹤[5]。波形特征追蹤方法只尋找搜索時窗內(nèi)特征點的像素波形結(jié)構(gòu)形態(tài)，不考慮道與道之間的相關(guān)性，準(zhǔn)確度不高，(邊緣檢測法[6])?；谙嚓P(guān)的自動追蹤是根據(jù)在搜索時窗內(nèi)，經(jīng)過隨機過程的互相關(guān)分析，計算相鄰地震道的互相關(guān)來反映同相軸的不連續(xù)性，又稱為相干圖像處理技術(shù)，主要應(yīng)用于描述地震數(shù)據(jù)的空間連續(xù)性，檢測相干性變差位置處的異常地質(zhì)體(如斷層、河道等)。Bahorich, M.等[7]提出了第一代相干性計算方法即互相關(guān)法，其缺點是利用數(shù)據(jù)較少，分辨率不高且對噪音的壓制不好；Marfurt等[8]提出了第二代相干性計算方法，相似系數(shù)法，利用了多道的信息進行相干性的計算，提高了相干剖面的分辨率，而且對噪音的壓制有著明顯的增強；Gersztenkorn等[9]提出了第三代相干性計算方法即本征結(jié)構(gòu)法，該方法在穩(wěn)定性方面以及噪音壓制方面均優(yōu)于前兩代相干性算法，但其計算量也遠遠大于前兩代算法。此外，還有其他一些算法，如局部結(jié)構(gòu)熵、高階統(tǒng)計量等。利用相干體技術(shù)檢測同相軸、提取同相軸傾角的原理是利用同相軸的橫向連續(xù)性較好，相干性較強；而地質(zhì)異常體、噪音的橫向連續(xù)性較差，相干性較弱。同相軸在相干計算后的相干剖面上更為明顯。同時由于在計算相干性的過程中對傾角進行掃描，使用與界面傾角相一致的傾角進行相干性計算時得到的相干性最大，因此可以獲得界面的傾角，并且能夠用于同相軸追蹤等后續(xù)工作。

為了提高復(fù)雜地區(qū)低信噪比資料層位自動追蹤的精確度，這里提出了一種基于圖像學(xué)的地震層位自動追蹤方法。圖像學(xué)去噪的方法很多，常見有小波變換[9]、全變分方法等。我們先應(yīng)用圖像學(xué)的全變分去噪技術(shù)壓制噪聲，然后基于相似系數(shù)的相干性算法對傾角進行掃描，把最大相干值對應(yīng)的傾角做為界面傾角，最后進行同相軸追蹤。

1 方法原理

我們主要應(yīng)用了圖像學(xué)中的全變分(TV)去噪技術(shù)和基于相似系數(shù)的相干成像技術(shù)。

1.1 全變分去噪技術(shù)

全變分去噪方法可以在保護邊緣信息的同時，降低圖像的噪聲?？紤]一個加噪的模型，令μ為原始的清晰圖像，μ0為被噪聲污染的圖像，即：

μ0(x,y)=μ(x,y)+n(x,y)

(1)

其中：n是具有零均值；方差為σ2的隨機噪聲；如果Ω表示圖像的定義域，像素點(x，y)∈Ω。

設(shè)Ω是Rn中的有界開子集，μ為局部可積函數(shù)，則其全變分定義為式(2)。

(2)

定義的全變分去噪能量泛函為式(3)。

(3)

通常有噪聲的全變分比沒有噪聲的全變分明顯大，最小化全變分可以消除噪聲，因此基于全變分的降噪可以歸結(jié)為如下最小化問題:

(4)

滿足約束條件：

(5)

(6)

如果一味地對圖像的全變分進行極小化，那么代表細節(jié)與紋理的許多自身特征也會被一并抹掉。為此，ROF模型用于圖像降噪的出發(fā)點是最小化“能量泛函”，即

(7)

它的前一項要求輸出μ的全變分盡可能小，稱為平滑項；后一項則要求μ與μ0盡可能相近，稱為數(shù)據(jù)保真項，它主要起保留圖像特征和降低圖像失真度的作用。參數(shù)λ用來平衡這兩個相互沖突的要求。

由梯度下降法，可以得到TV平滑模型：

(8)

在全變分中引入一個小的正數(shù)β，對它進行正則化，即用

(9)

(10)

最終的正則化模型為式(11)。

(11)

求解 公式(7)的偏微分方程，利用邊界約束條件推出：

(12)

(13)

=>

λ(u-u0)-Δu=0

(14)

離散迭代格式為：

(15)

正則化參數(shù)λ的取值非常關(guān)鍵，算法步驟如下：

1)讀入帶有噪聲的數(shù)據(jù)u0。

2)初始化參數(shù)：n=0，Δt=0.25，u0=u0，divp=0。

3)當(dāng)n小于最大迭代次數(shù)時，循環(huán)執(zhí)行如下操作：n=n+1，根據(jù)離散迭代公式計算下一步的un,計算擴散項的值。

4)結(jié)束迭代，最后一次為去噪的結(jié)果。

1.2 相似系數(shù)的相干性計算

Chopra, S.等[10]提出基于相似系數(shù)的相干性計算，可以利用與中心道相鄰的多道信息進行計算，如圖1所示。

圖1 相似系數(shù)法參與計算的地震道Fig.1 The similarity coefficient method is involved in the calculation of seismic trace

在進行相干性計算時，首先定義一個分析窗口(橢圓窗或者矩形窗)，在分析窗之內(nèi)以分析點為中心包含J個地震道，如圖2所示。方位角Φ和長短軸a、b決定著分析窗的方向和大小，當(dāng)Φ0、Φa、Φ、a、b確定后分析窗內(nèi)包含的J的地震道就確定了相干時窗的大小一般根據(jù)地震反射波的視周期T而定，通常取T/2～T/3。當(dāng)計算的相干時窗小于T/2時，不能描述一個完整的波峰或波谷，噪音與層位的影響突出；當(dāng)計算時窗大于T/2時，包含了多個反射波同相軸，可能會降低分辨率。就參加運算的道數(shù)來說，從斷層成像清晰度和隨機噪聲壓制程度看，一般參與相干計算的道數(shù)越多，平均效應(yīng)越大，對斷層的分辨率越低，這時突出的主要是大斷層；相反，相干道數(shù)少，平均效應(yīng)小，就會提高對地層邊界、斷層、特別是突出了對小斷層的分辨率。

圖2 相似系數(shù)法的計算窗口選取[3]Fig.2 The calculation window selection of similarity coefficient method(a)Φ0 Inine方向的方位角;(b)Φa 分析窗長軸的方位角

(16)

其中：uj為分析窗內(nèi)第j道；xj、yj分別為該道距離中心道的距離；p、q分別表示x、y方向的傾角。圖2是當(dāng)p、q為0時的情況，當(dāng)p、q分別為p1、q1時，如圖3所示，p1、q1分別表示分析窗沿x、y方向的傾角。

Cs(t,p,q)=

(17)

用不同的傾角進行掃描可以獲得使Cs(t,p,q)最大的p和q，即可以獲得同相軸的傾角，然后可以進行同相軸追蹤等后續(xù)工作。

在二維情況下，相干性Cs可以表示為

(18)

基于相似系數(shù)的相干性計算利用了多道的信息，考慮同相軸整體的連續(xù)性，因此,分辨率優(yōu)于互相關(guān)算法，并且對噪音的壓制也較互相關(guān)算法好。同時在計算相干性的時候,可以根據(jù)文獻[11]計算得到同相軸上每一點的傾角，因此我們采用這種方法進行傾角的拾取與同相軸的追蹤。

2 自動追蹤技術(shù)流程及數(shù)據(jù)結(jié)構(gòu)

自動追蹤技術(shù)實現(xiàn)流程如圖3所示。

在相干計算后，得到的構(gòu)造傾角可以進行同相軸自動追蹤。通過計算機確定拾取層位起始點坐標(biāo)，自動計算出參考道的時間深度h，檢索出該點相干性最大的對應(yīng)傾角，并記錄。然后以下一道為參考道，將時延t賦給h,循環(huán)計算，直至計算到相干時窗的最后一道。實現(xiàn)的追蹤流程為圖4所示。

計算得到的一個同相軸的追蹤結(jié)果采用鏈?zhǔn)酱鎯Ψ绞奖４妫鋬?yōu)點是存儲量較小，相對順序存儲方式而言，它的插入和刪除元素速度較快，方便編輯追蹤結(jié)果。而且它沒有空間限制，存儲元素的個數(shù)無上限，基本只與內(nèi)存空間大小有關(guān)。給每一個鏈條分配一個編號用以區(qū)別剖面上其他同相軸追蹤結(jié)果，方便實現(xiàn)層位追蹤結(jié)果檢索的功能。Patel等[12]提出將識別層位轉(zhuǎn)化為網(wǎng)格模型的方法，分離錯誤層位連接。具體數(shù)據(jù)結(jié)構(gòu)如圖5所示，每條鏈表存儲一條同相軸，包含一個頭結(jié)點和若干子節(jié)點。頭結(jié)點中存儲同相軸的編號、節(jié)點個數(shù)等，每一個節(jié)點存儲每一個追蹤到的該同相軸上的點坐標(biāo)和局部傾角。最后所有鏈表的頭結(jié)點組成一個頭結(jié)點數(shù)組，便于檢索，數(shù)據(jù)結(jié)構(gòu)如圖5所示。

3 應(yīng)用實例

在基于linux平臺的中石化iCluster系統(tǒng)中開發(fā)了層位自動追蹤功能，進行實際數(shù)據(jù)測試。選用塔河資料為例，主測線范圍在450～1100，CDP范圍450～1500，工區(qū)內(nèi)地層較為清晰。

3.1 降噪效果分析

這里采用峰值信噪比(PSNR)[12]和EPI進行評價。PSNR值大時，降噪處理效果就好。EPI是邊緣細節(jié)保留情況的主要評價標(biāo)準(zhǔn)。

圖3 傾角為p1、q1時的計算窗口選取Fig.3 The calculation window selection of the inclination of p1 and q1(a)橢圓形分析窗的情況;(b)矩形分析窗的情況

圖4 技術(shù)實現(xiàn)流程Fig.4 Technical realization process

圖5 數(shù)據(jù)結(jié)構(gòu)設(shè)計Fig.5 Design of data structure

(19)

(20)

圖6是原始數(shù)據(jù)的部分剖面圖，圖7是采用全變分方法，λ=1.0時經(jīng)處理得到的部分剖面圖。處理后的地震剖面同相軸連續(xù)性更好，層間結(jié)構(gòu)也更清晰。

表1 降噪效果Tab.1 Noise reduction effect

圖6 原始的剖面圖Fig.6 The original profile

圖7 全變分處理后的剖面圖Fig.7 Profile of total variation after treatmene

3.2 層位追蹤結(jié)果

都采用基于相似系數(shù)的相干計算，分兩種情況進行對比：①是不經(jīng)過全變分降噪處理;②是經(jīng)過全變分降噪處理，圖7是不經(jīng)過全變分降噪處理的追蹤結(jié)果，綠色是追蹤軌跡。圖8是經(jīng)過全變分降噪處理的追蹤結(jié)果，綠色是追蹤軌跡。比較圖8、圖9發(fā)現(xiàn)，不用全變分處理的追蹤結(jié)果，跳動大，局部有異常，經(jīng)過全變分降噪處理后追蹤結(jié)果有改善，符合較好。

圖8 不經(jīng)過全變分處理的追蹤結(jié)果Fig.8 Tracking results without total variation

圖9 經(jīng)過全變分處理的追蹤結(jié)果Fig.9 The tracking results of the total variation

4 結(jié)論

全變分方法較好地實現(xiàn)了去噪和保護邊緣的統(tǒng)一，保護了同相軸。基于相似系數(shù)的相干算法分辨率高，計算后得到的同相軸傾角可用于后續(xù)的同相軸自動追蹤。該方法在實際資料測試中得到光滑連續(xù)的同相軸，與地質(zhì)層位相符，取得良好的效果。

[1] 謝輝.基于地震層位趨勢約束的小層系地質(zhì)層位隨機模擬[J].科技創(chuàng)新與應(yīng)用,2013.31:83-85.XH.Stochasticsimulationofgeologicalhorizonofsmallstratabasedonthetrendofseismichorizon[J].ScienceandtechnologyinnovationandApplication, 2013,31:83-85.(InChinese)

[2]HOYESJ,CHERETT．Areviewofglobalinterpretationmethodsforautomated3Dhorizonpicking[J].TheLeadingEdge，2011,30(1):38-47.

[3] 彭文,熊曉軍,韓小俊.基于高階統(tǒng)計量的層位自動追蹤方法[J].新疆石油地質(zhì)，2006,6(27):743-745.PENGW,XIONGXJ,HANXJ.Horizonautomatictrackingmethodbasedonhigherorderstatistics[J].XinjiangPetroleumGeology,2006,6(27):743-745.(InChinese)

[4] 鄭公營,曾婷婷.地震層位自動追蹤技術(shù)研究[J].物探化探計算技術(shù),2013,35(6):711-716.ZHENGGY,ZENGTT.Researchonautomatictrackingtechnologyofseismichorizon[J].Geophysicalprospectingandgeochemicalexploration, 2013,35(6):711-716.(InChinese)

[5] 李紅星，劉財，陶春輝．圖像邊緣檢測方法在地震剖面同相軸自動檢測中的應(yīng)用研究[J].地球物理學(xué)進展，2007，22(5)：1607-1610.LIHX,LIUC,TAOCH.Applicationofimageedgedetectionmethodintheautomaticdetectionofseismicprofileinphase[J].ProgressinGeophysics，2007，22(5)：1607-1610.(InChinese)

[6]BAHORICHMS,FARMERSL.3Dseismicdiscontinuityforfaultsandstratigraphicfeatures[J].TheLeadingEdge,1995,14(10):1053-1058.

[7]MARFURTKJ,KIRLINR,FARMERSL,etal.3DSeismicAttributesUsingaSemblance-basedCoherencyAlgorithm[J].Geophysics,1998,63(4):1150-1165.

[8]GERSZTENKORNADAM,MarfurtKJ.Eigenstructure-basedCoherenceComputationsasanAidto3DStructuralandStratigraphicMapping[J].Geophysics,1999,64(5): 1468-1479.

[9]GALIANA-MERINOJJ，ROSA-HERRANZJL，RosaCintasS,eta1.SeismicWaveTool：continuousanddiscretewaveletanalysisandfilteringformultichannelseismicdata[J].ComputerPhysicsCommunications,2013,184(1):162-171.

[10]CHOPRAS,MARFURTKJ.Seismicattributeexpressionofdifferentialcompaction[J].TheLeadingEdge, 2012, 31(31): 1418-1422.

[11]MARFURTKJ．Robustestimatesof3Dreflectordipandazimuth[J].Geophysics, 2006,71(4)：29-40.

[12]PATELD，BRUCKNERS，VIOLAI,eta1.Seismicvolumevisualizationforhorizonextraction[C].PacificVisualizationSymposium(PacificVis),2010IEEE.IEEE,Taipei,2010:73-80.

[13]公成敏. 全變分原理在地震數(shù)據(jù)去噪中的應(yīng)用[J].計算機與數(shù)字工程,2014,297(42):1271-1274.GONGCM.Thevariationalprincipleanditsapplicationindenoisingofseismicdata[J].Compuuteranddigitalengineering,2014,297(42):1271-1274.(InChinese)

An automatic tracking method for seismic horizons based on image theory

LIU Xuyue, ZHOU Wei, ZHANG Bing, HUANG Jun

(Sinopec Geophysical Research Institute,Nanjing 211103,China)

Seismic data in complex area is affected by low signal-to-noise ratio, noise interference, horizon tracking results are not high, this paper proposes a method to solve the problem. The seismic horizon automatic tracking method based on image theory, the whole variation method to calculate the seismic profile, and then the coherent image processing technology to calculate the seismic profile. And then the phase axis automatic tracking. It is proved that the method can achieve good results using the actual data test.

total variation; coherent image processing; the angle of the same phase; automatic tracing of horizons

2015-12-28 改回日期：2016-03-01

中石化科技部項目(P14152)

劉旭躍(1985-)，男，工程師，現(xiàn)從事物探軟件工作， E-mail：Liuxy@sinopec.com。

1001-1749(2017)01-0064-07

P 631.4

A

10.3969/j.issn.1001-1749.2017.01.10