蔣旭東,曹俊興,胡江濤

(成都理工大學(xué)地球物理學(xué)院,四川成都610059)

精準(zhǔn)的層位信息是地震資料解釋的基礎(chǔ),是儲(chǔ)層預(yù)測(cè)的重要依據(jù)。層位信息被廣泛應(yīng)用于各種地震處理方法中,如速度分析和層析成像等,但人工層位追蹤需要巨大的人力成本,且非常耗時(shí),影響了地震勘探的效率。因此,開發(fā)一種穩(wěn)定且高效的自動(dòng)層位追蹤技術(shù)對(duì)地震資料的處理和解釋非常重要。

目前自動(dòng)層位追蹤的方法主要有三大類。第一類是基于相關(guān)的自動(dòng)追蹤技術(shù),主要利用第一代或第三代相干算法[1-2]。該方法通過道與道之間的相關(guān)性進(jìn)行同相軸追蹤,其抗噪性和穩(wěn)定性較好,但計(jì)算量大,而且當(dāng)相鄰層位存在相似波形時(shí),容易出現(xiàn)串層現(xiàn)象,降低了結(jié)果的可靠性。第二類是基于人工神經(jīng)網(wǎng)絡(luò)的自動(dòng)追蹤技術(shù),由ALBERTS等[3]提出,由于其在層位追蹤中取得了較好的效果,受到廣泛關(guān)注[4-5],但該技術(shù)需要大量的樣本進(jìn)行學(xué)習(xí)和訓(xùn)練,計(jì)算量較大;同時(shí),不同工區(qū)的地震數(shù)據(jù)具有不同的特征,導(dǎo)致工區(qū)間不能共享學(xué)習(xí)。第三類是基于圖像的自動(dòng)追蹤技術(shù)。HALE等[6]和NAEINI等[7]將結(jié)構(gòu)張量應(yīng)用于斷層及層位識(shí)別,利用一系列的方向?yàn)V波器,提取特征方向,進(jìn)而獲得層位發(fā)育的主方向。該方法計(jì)算效率高,但當(dāng)層位形態(tài)較復(fù)雜時(shí),會(huì)影響追蹤結(jié)果的精度和穩(wěn)定性。

本文利用希爾伯特變換對(duì)原始地震數(shù)據(jù)進(jìn)行預(yù)處理,降低假頻和噪聲的影響,增強(qiáng)其信噪比和連續(xù)性,進(jìn)而提高結(jié)構(gòu)張量算法估計(jì)角度信息的穩(wěn)定性,然后對(duì)角度信息進(jìn)行預(yù)處理,進(jìn)一步改善角度信息的可靠性。在角度信息和一系列約束條件下,沿層位發(fā)育方向進(jìn)行雙方向的層位匹配追蹤,獲得優(yōu)選結(jié)果。最后用理論沉積模型和實(shí)際數(shù)據(jù)進(jìn)行了測(cè)試并與常規(guī)方法結(jié)果進(jìn)行了對(duì)比,以驗(yàn)證方法的穩(wěn)定性和適應(yīng)性。

1 方法原理

1.1 基于希爾伯特變換的數(shù)據(jù)預(yù)處理

地震資料在采集過程中,往往存在噪聲干擾,進(jìn)而導(dǎo)致同相軸連續(xù)性弱和目的層位信噪比低等;同時(shí),由于道間距較大等原因,會(huì)出現(xiàn)空間假頻現(xiàn)象。這些因素會(huì)降低結(jié)構(gòu)張量估算角度信息的穩(wěn)定性。為此,利用希爾伯特變換進(jìn)行預(yù)處理,求取信號(hào)的包絡(luò),提高數(shù)據(jù)的信噪比和同相軸連續(xù)性,降低假頻和噪聲的影響,再利用包絡(luò)數(shù)據(jù)進(jìn)行角度信息估計(jì)。

(1)

其深度域包絡(luò)可表示為:

(2)

包絡(luò)數(shù)據(jù)突出了數(shù)據(jù)的低頻信息,可增強(qiáng)數(shù)據(jù)的連續(xù)性,降低噪聲及空間假頻對(duì)結(jié)果的影響,為角度信息的估計(jì)提供更好的基礎(chǔ)數(shù)據(jù)。

1.2 結(jié)構(gòu)張量算法

在角度信息提取中,常用的方法有結(jié)構(gòu)張量算法、局部傾斜疊加和平面波分解法[8]等。就計(jì)算速度而言,結(jié)構(gòu)張量算法最優(yōu);就分辨率而言,結(jié)構(gòu)張量算法雖然略低于局部傾斜疊加法,但結(jié)構(gòu)張量算法求取的是一個(gè)鄰域內(nèi)角度的平均值,所以結(jié)果更光滑連續(xù),更利于同相軸穩(wěn)定追蹤[9]。所以,我們選用結(jié)構(gòu)張量算法進(jìn)行角度信息提取,同時(shí),為了提高其穩(wěn)定性,使用二維高斯函數(shù)進(jìn)行迭代濾波處理,并對(duì)其結(jié)果進(jìn)行平滑處理,以提高穩(wěn)健性。最后采用WEICKER[10]提出的算子進(jìn)行角度求取。

二維地震包絡(luò)剖面可表示為一個(gè)M×N(M為地震道深度采樣點(diǎn),N為地震道數(shù)目)的矩陣I,其某一元素的梯度向量為:

(3)

i=1,2,…,Mj=1,2,…,N

式中:Ix(i,j)為剖面在x方向的梯度;Iz(i,j)為剖面在z方向的梯度。為了避免噪聲等影響,VLIET等[11]提出了結(jié)構(gòu)張量算法:

TS(i,j)=G(i,j)×GT(i,j)=

(4)

式中:TS(i,j)為點(diǎn)的結(jié)構(gòu)張量。利用(4)式能消除局部平均的影響,同時(shí)能增強(qiáng)同一方向的梯度矢量。進(jìn)一步利用二維高斯函數(shù)((5)式)對(duì)結(jié)構(gòu)張量進(jìn)行濾波:

(5)

式中:Kp=1/(2πρ2)·exp[-(x2+y2)/2ρ2]為二維高斯函數(shù);ρ代表分析結(jié)構(gòu)的鄰域大小;Iσ=KσI為原始數(shù)據(jù)I經(jīng)過高斯函數(shù)Kσ平滑后的結(jié)果;σ為標(biāo)準(zhǔn)方差。求解(6)式半正定矩陣TS1的最大特征向量即可獲得方向信息:

(6)

其中,λ表示特征值向量;E為單位矩陣;x表示特征向量。特征值可解析表示為:

(7)

定義結(jié)構(gòu)一致性因子α和最大特征值向量的斜邊距g為:

進(jìn)一步可得層位方向信息為[10]:

式中:v,w為方向向量的兩個(gè)標(biāo)量;ε為平滑系數(shù);θ為該層位點(diǎn)處地層的發(fā)育角度。平滑處理能夠提高角度值求取的穩(wěn)定性,進(jìn)一步降低層位突變和噪聲的影響。

1.3 層位追蹤

HALE等[12]利用結(jié)構(gòu)張量提取層位的主發(fā)育方向,進(jìn)而實(shí)現(xiàn)結(jié)構(gòu)導(dǎo)向?yàn)V波。基于角度信息,沿層位發(fā)育方向,結(jié)合一系列約束條件進(jìn)行層位追蹤。如圖1 所示,首先在層位追蹤的目的層處選定一個(gè)層位較清晰且同相軸較連續(xù)的點(diǎn)作為層位追蹤的種子點(diǎn),然后對(duì)該種子點(diǎn)進(jìn)行自動(dòng)修正,使其為準(zhǔn)確的波峰點(diǎn)或波谷點(diǎn)。利用角度信息可確定該層位的相鄰空間點(diǎn)位置,由于層位的橫向空間位置x可由數(shù)據(jù)道頭獲得,因此層位追蹤僅需獲取層位的深度變化:

(12)

式中:z為層位的深度值;Δx表示兩道之間的間隔;θ為該層位點(diǎn)處地層的發(fā)育角度。為了提高相鄰層位點(diǎn)的提取精度,在上述計(jì)算的深度位置處沿深度方向取一定窗口,優(yōu)選層位結(jié)果。該深度時(shí)窗內(nèi)的數(shù)據(jù)可表示為以下數(shù)據(jù)集合:

(13)

其中,I(i,j+1)為深度zj+1位置點(diǎn)對(duì)應(yīng)的包絡(luò)值,K為窗口長度。由于角度指向性強(qiáng),所以K值可相對(duì)較小(即小窗口數(shù)據(jù)),它能有效縮短計(jì)算時(shí)間,進(jìn)而提高效率。在追蹤層位時(shí),將上一層位點(diǎn)對(duì)應(yīng)的包絡(luò)值I(i,j)與窗口內(nèi)地震數(shù)據(jù)集合進(jìn)行一一求差,選取差值最小點(diǎn)作為最終層位點(diǎn)。同時(shí),在層位追蹤中,通常沿波形的波谷或波峰(對(duì)應(yīng)于數(shù)據(jù)中的極值點(diǎn))進(jìn)行追蹤,并進(jìn)行點(diǎn)位自動(dòng)修正,使其為準(zhǔn)確的波峰或者波谷點(diǎn)。在追蹤波峰或波谷的一些特殊情況下,還可加上額外的對(duì)比準(zhǔn)則,即將I(i,j)所對(duì)應(yīng)的梯度信息G(i,j)與窗口內(nèi)數(shù)據(jù)集的梯度信息進(jìn)行一一對(duì)比,選取包絡(luò)值和梯度值差異最小的點(diǎn)作為最終層位點(diǎn)。層位自動(dòng)追蹤流程如圖2所示。

圖2 層位自動(dòng)追蹤流程

2 數(shù)值試驗(yàn)

2.1 模型數(shù)據(jù)測(cè)試

首先采用3層傾斜模型進(jìn)行測(cè)試,該模型深度域剖面如圖3a所示。3層的地層的傾斜度分別為10°,15°和20°;道間距為30m,水平距離為10000m;深度采樣間隔為10m,深度為4000m。道間距較大時(shí)會(huì)產(chǎn)生空間假頻,所以采用該模型來測(cè)試方法對(duì)假頻的適應(yīng)能力。對(duì)該模型原始數(shù)據(jù)進(jìn)行頻譜分析,結(jié)果如圖3b所示,可以看出,該數(shù)據(jù)的空間假頻明顯。對(duì)原始數(shù)據(jù)進(jìn)行角度提取,結(jié)果如圖3c所示,可以明顯看出,在傾斜層處角度提取結(jié)果不穩(wěn)定。

采用本文方法對(duì)原始數(shù)據(jù)進(jìn)行處理,結(jié)果如圖3d 所示,可以看出,預(yù)處理后的數(shù)據(jù)連續(xù)性明顯增強(qiáng),從其頻譜(圖3e)可看出假頻也被有效地消除。采用本文方法對(duì)預(yù)處理后的數(shù)據(jù)進(jìn)行角度提取,結(jié)果如圖3f所示。對(duì)比圖3c和圖3f可看出,采用本文方法提取的角度結(jié)果更可靠。

再采用經(jīng)典的深度域理論模型進(jìn)行測(cè)試。該模型包含多種典型的地質(zhì)特征,如圖4a所示。由于理論模型數(shù)據(jù)不包含噪聲,為接近實(shí)際情況,對(duì)該模型數(shù)據(jù)添加信噪比為0.3的高斯隨機(jī)噪聲,如圖4b所示。經(jīng)過預(yù)處理后,數(shù)據(jù)的方向性和連續(xù)性均有所加強(qiáng)(圖4c)。從圖4中明顯可以看出,噪聲導(dǎo)致同相軸不連續(xù),經(jīng)過預(yù)處理后的數(shù)據(jù)連續(xù)性明顯增強(qiáng),預(yù)處理后的數(shù)據(jù)與原始數(shù)據(jù)也有很好的一致性。

圖3 模型及頻譜結(jié)果對(duì)比a 原始模型深度域剖面;b 原始剖面頻譜;c 原始剖面角度提取結(jié)果;d 預(yù)處理后的模型深度域剖面;e 經(jīng)過變換后的剖面頻譜;f 預(yù)處理后剖面的角度提取結(jié)果

為進(jìn)一步驗(yàn)證預(yù)處理的效果,對(duì)圖4所示的3種數(shù)據(jù)分別進(jìn)行角度提取,結(jié)果如圖5所示。圖5a為對(duì)原始無噪數(shù)據(jù)進(jìn)行角度提取的結(jié)果,圖5b為對(duì)加噪后的模型數(shù)據(jù)進(jìn)行角度提取的結(jié)果。相比于圖5a,噪聲的存在使得獲取的角度有較大的誤差,并且存在不穩(wěn)定現(xiàn)象。圖5c為對(duì)預(yù)處理后的模型數(shù)據(jù)進(jìn)行角度提取的結(jié)果,該結(jié)果與圖5a的結(jié)果具有高度一致性。經(jīng)過預(yù)處理之后,能夠獲得與無噪數(shù)據(jù)類似的角度估計(jì)結(jié)果,說明該預(yù)處理方法有較強(qiáng)的抗噪性,能夠提高角度提取穩(wěn)定性。

圖4 預(yù)處理前、后的理論沉積模型數(shù)據(jù)a 理論模型數(shù)據(jù);b 添加噪聲后模型數(shù)據(jù);c 預(yù)處理后的模型數(shù)據(jù)

圖5 采用不同數(shù)據(jù)的角度提取結(jié)果a 原始數(shù)據(jù);b 加噪處理后的數(shù)據(jù);c 希爾伯特預(yù)處理后角度提取數(shù)據(jù)

在該理論模型上選取兩種地質(zhì)特征的同相軸進(jìn)行層位自動(dòng)追蹤,結(jié)果如圖6所示,圖中,紅色五角星為在兩種地質(zhì)特征的同相軸上選取的種子點(diǎn),綠線為傾斜層位的追蹤結(jié)果,藍(lán)線為褶皺層位的追蹤結(jié)果。從圖6中可以看出,層位追蹤結(jié)果可靠,特別是在有斷層的地方,能夠準(zhǔn)確地識(shí)別出斷層的走向和位置。

2.2 實(shí)際數(shù)據(jù)應(yīng)用

選取龍門山前陸盆地深層海相碳酸鹽巖深度域疊后數(shù)據(jù)進(jìn)行測(cè)試。龍門山前陸盆地海相地層主要為碳酸鹽巖,從震旦系到中三疊統(tǒng),其中缺失志留系、泥盆系和石炭系海相地層。圖7為該區(qū)原始地震剖面,顯示該區(qū)地震數(shù)據(jù)信噪比較低,斷層較多。在層位追蹤時(shí),容易出現(xiàn)串層等現(xiàn)象,追蹤難度大。經(jīng)過預(yù)處理后,對(duì)該剖面進(jìn)行角度提取,結(jié)果如圖8所示,可見,提取的角度信息與構(gòu)造角度吻合度較好,斷層在角度剖面上均有顯示。

圖6 理論沉積模型層位提取結(jié)果

選取層位1和層位2作為目的追蹤層位。采用本文方法的層位追蹤結(jié)果如圖9所示,圖中黃色五角星為層位追蹤的種子點(diǎn)。由于層位1的中間部分有較大的破碎帶,所以其種子點(diǎn)分別選在破裂帶的兩邊。從圖9可以看出,層位追蹤結(jié)果與實(shí)際層位吻合度高。

圖7 原始地震剖面

圖8 角度提取結(jié)果

圖9 層位追蹤結(jié)果

為了進(jìn)一步檢驗(yàn)層位追蹤的結(jié)果,將其與兩種常用的層位自動(dòng)追蹤方法的結(jié)果進(jìn)行了對(duì)比,結(jié)果如圖10所示?；谙喔傻淖詣?dòng)追蹤方法利用相干算法在設(shè)定時(shí)窗范圍中提取相關(guān)值最大的點(diǎn)作為層位點(diǎn);基于波形相似度的自動(dòng)追蹤方法通過衡量相鄰地震道在時(shí)窗范圍內(nèi)的相似度,同時(shí)加入斷層等構(gòu)造信息進(jìn)行追蹤。從圖10可以看出,3種方法的結(jié)果基本相同。在使用方便程度方面,本文方法僅需一個(gè)種子點(diǎn)進(jìn)行追蹤,即可實(shí)現(xiàn)同相軸的雙向追蹤;基于相干的自動(dòng)追蹤方法需給定追蹤的方向和多個(gè)種子點(diǎn);基于波形相似的自動(dòng)追蹤方法也僅需一個(gè)種子點(diǎn),但當(dāng)波形相似度變?nèi)鯐r(shí)會(huì)追蹤停止,在一些較大的斷層區(qū)域,追蹤結(jié)果的連續(xù)性弱。

圖11為追蹤結(jié)果的局部放大顯示。從圖11a中可以看出,層位1和層位2連續(xù)性好且信噪比高,采用這3種方法都取得了較好的結(jié)果。圖11b中層位1發(fā)育一些小斷層,采用基于相干的方法時(shí)追蹤結(jié)果出現(xiàn)串層現(xiàn)象,采用基于相似度的方法時(shí)追蹤結(jié)果未能在小斷層處斷開層位,而采用本文方法能夠很好地與斷層信息相吻合;層位2處的結(jié)果也顯示采用本文方法的追蹤結(jié)果更貼合層位的走勢(shì)。圖11c為采用這3種方法在破碎帶區(qū)域的追蹤結(jié)果,可以看出,采用基于相干和基于相似度的兩種方法的追蹤結(jié)果出現(xiàn)明顯的不貼合層位現(xiàn)象,而采用本文方法的追蹤結(jié)果在一定時(shí)窗內(nèi)優(yōu)選層位信息,能夠較好地貼合連續(xù)性強(qiáng)的層位信息,為后續(xù)進(jìn)一步修正層位信息提供依據(jù)。

圖10 不同方法的自動(dòng)追蹤結(jié)果

圖11 不同方法自動(dòng)追蹤結(jié)果的局部放大顯示a 區(qū)域A;b 區(qū)域B;c 區(qū)域C

3 結(jié)論

常用自動(dòng)層位追蹤算法通常存在穩(wěn)定性不高且適應(yīng)性不強(qiáng)的問題。本文結(jié)合結(jié)構(gòu)張量算法,提出了一種自動(dòng)層位追蹤方法。對(duì)地震剖面利用希爾伯特變換提取包絡(luò)信息,改善了結(jié)構(gòu)張量方法求取方向信息的穩(wěn)定性;對(duì)結(jié)果進(jìn)行經(jīng)過高斯濾波和平滑等處理,進(jìn)一步提高角度信息的穩(wěn)定度;最后在角度信息和同相軸特征的約束下實(shí)現(xiàn)層位自動(dòng)追蹤。理論沉積模型及實(shí)際二維地震剖面數(shù)據(jù)的測(cè)試表明該方法結(jié)果相對(duì)穩(wěn)定且與同相軸吻合度高。與基于相干和基于相似度的兩種常用方法的對(duì)比結(jié)果表明,本文方法在層位自動(dòng)追蹤中具有較高的穩(wěn)定性。