基于同步擠壓改進(jìn)短時傅立葉變換的分頻螞蟻追蹤在斷裂識別中的應(yīng)用

黃葦,周捷,高利君,王勝利,嚴(yán)海滔

(1.中國石油化工股份有限公司 西北油田分公司勘探開發(fā)研究院，新疆 烏魯木齊 830011； 2.成都理工大學(xué) 地球物理學(xué)院，四川 成都 610059)

0 引言

近些年由于油氣勘探力度不斷增大，淺層構(gòu)造油氣藏幾乎已完成開采；油氣勘探工作者逐漸將淺層構(gòu)造油氣藏的目光轉(zhuǎn)向深層油氣藏。由于大地濾波、吸收衰減等因素使得地震波能量在向下傳播過程中發(fā)生嚴(yán)重衰減，因此，對深層油氣藏的探尋變得愈發(fā)艱難。油氣從生成層轉(zhuǎn)移至儲集層時需要通道，所以，儲層發(fā)育的裂縫扮演著不可或缺的角色；它將影響著儲層的滲透性以及連通性，直接決定著油氣的產(chǎn)量；因此，對于微小斷層的刻畫或者裂縫發(fā)育帶檢測已然成為油氣儲層預(yù)測中非常重要而又迫切的問題。

裂縫、斷層檢測方法經(jīng)過多年的發(fā)展研究，自1995年Bahorich M[1]提出地震相干數(shù)據(jù)方法后，Amoco公司將地震相干技術(shù)應(yīng)用于地震資料斷裂系統(tǒng)檢測，取得了很好的效果。1998年Marfurt[2-4]基于相似性算法，提出了比第一代相干算法更能描繪斷層的第二代相干體算法；1999年Gersztenkorn等基于本征結(jié)構(gòu)以及協(xié)方差矩陣，提出了第三代相干算法，從而求取多道地震數(shù)據(jù)彼此的相關(guān)性；第三代相干算法相比第一代和第二代相干算法，其提高了地震數(shù)據(jù)的縱向、橫向分辨率，且將該算法推廣到三維地震數(shù)據(jù)體中，不需要層位約束，但是對于許多微小斷裂以及次生裂縫無法進(jìn)行精細(xì)的刻畫。王西文等[5]把小波變換和相干算法相融合，該算法比傳統(tǒng)相干算法有更高的分辨率，但刻畫小裂縫的效果依舊不理想。1999年Dorigo等[6]提出了基于蟻群算法的螞蟻追蹤技術(shù)，該算法模仿螞蟻覓食遇到障礙所釋放的信息素來追蹤斷層以及裂縫。目前，螞蟻追蹤技術(shù)是公認(rèn)的一種能有效刻畫斷層以及裂縫的技術(shù)。但在實際應(yīng)用中該技術(shù)仍存在些許精度上的缺陷，如許多高陡構(gòu)造，由大斷層而衍生出來的微小裂縫、斷層等，螞蟻追蹤技術(shù)將無法準(zhǔn)確地雕刻。

Partyka等[7]于1999年提出了將全頻帶的地震數(shù)據(jù)體分解為各個頻率的單頻數(shù)據(jù)體的頻譜分解技術(shù)，以此反映出全頻帶地震數(shù)據(jù)中無法凸顯出的地質(zhì)信息。2009年Zeng等[8]利用頻譜分解技術(shù)得到不同頻率的地震數(shù)據(jù)體，結(jié)果顯示單頻地震數(shù)據(jù)體比全頻帶數(shù)據(jù)體對地下異常體的區(qū)域范圍能有更清晰的勾勒，并且能刻畫更多細(xì)節(jié)。以上學(xué)者都為頻率域的資料解釋做了極大的貢獻(xiàn)。

Daubechies等[9]將小波變換后的頻譜值在時頻域內(nèi)進(jìn)行擠壓以及重新排列，提出了將時頻聚焦能力提高的同步擠壓算法。2017年Yu G等[10]將短時傅立葉變換窗函數(shù)延展，使其結(jié)果比短時傅立葉變換的分辨率更高，但是，它仍會因為海森堡原理的制約以及交叉項干擾使時頻分析精度受到很大限制。

因此，根據(jù)同步擠壓算法原理，通過對改進(jìn)短時傅立葉變換后的時頻譜擠壓重排來發(fā)展同步擠壓改進(jìn)短時傅立葉變換，理論信號表明，該方法具有比傳統(tǒng)時頻分析方法更高的時頻分辨率。結(jié)合之前學(xué)者的研究成果[11-19]，本文利用同步擠壓改進(jìn)短時傅立葉變換在時頻聚焦上的優(yōu)勢，對某工區(qū)的三維地震數(shù)據(jù)體進(jìn)行頻譜分解從而獲取各個頻率的單頻數(shù)據(jù)體；再應(yīng)用螞蟻追蹤技術(shù)對單頻數(shù)據(jù)體進(jìn)行計算。其結(jié)果表示該方法可以對微小裂縫進(jìn)行細(xì)膩的勾勒，并證明了其在實際應(yīng)用中的準(zhǔn)確性及實用性。

1 方法原理

傅立葉變換可以算是時頻分析領(lǐng)域的鼻祖，但其無法對非平穩(wěn)信號進(jìn)行很好地刻畫；而現(xiàn)實中的信號幾乎為非平穩(wěn)信號，尤其是地震信號，傅立葉變換早已無法達(dá)到人們對精度的要求。短時傅立葉變換通過增加窗函數(shù)去截斷非平穩(wěn)信號，并將整個非平穩(wěn)信號分割成數(shù)個小段，這時便可把窗函數(shù)里的每段信號看作平穩(wěn)信號，再對每個小段信號分別做傅立葉變換便得到了二維時頻譜圖。

假設(shè)信號x(t)∈L2(R)，則信號的數(shù)學(xué)表達(dá)式為：

x(t)=Asin(ct+φ) ，

(1)

式中：A為幅值,cm;c為角頻率,Hz;t為時間，s；φ為初始相位角；L2(R)是實數(shù)域中的平方可積函數(shù)空間；x(t)表示信號。

信號x(t)的短時傅立葉變換(Sstft)表達(dá)式為：

式中：f為頻率,Hz；τ為時窗長度，常數(shù)；i表示虛數(shù)；g(t-τ)為隨τ的變化在時間軸上運動的窗函數(shù)。但短時傅立葉變換的窗函數(shù)固定且唯一，這意味著信號截斷的位置以及時窗長度對計算結(jié)果的影響特別大。

將式(2)離散化得到離散傅立葉變換(Sstftx)表達(dá)式為：

(3)

1.1 改進(jìn)短時傅立葉變換

由于短時傅立葉變換時窗無法改變這一缺點限制了該方法提高時頻聚焦性，想要克服這一缺陷必須要對時窗函數(shù)進(jìn)行改進(jìn)。由式(2)，分別求出信號和窗函數(shù)的復(fù)共軛，則可以得到：

(4)

(5)

令t-r=r，則有

(6)

(7)

原始Sstft的表達(dá)式為式(2)，則改進(jìn)短時傅立葉變換(Gstft)的表達(dá)式為：

1.2 同步擠壓改進(jìn)短時傅立葉變換

根據(jù)改進(jìn)短時傅立葉變換的公式，首先，計算Gstft對時間t導(dǎo)數(shù)為：

(9)

則可以得到信號的瞬時頻率表達(dá)式：

(10)

根據(jù)式(10)得到的瞬時頻率，可以得到同步擠壓改進(jìn)短時傅立葉變換(Sgstft)值的表達(dá)式為：

式中：η為同步擠壓變換后的時頻譜頻率值，Hz。

同步擠壓改進(jìn)短時傅立葉變換則是對短時傅立葉變換改進(jìn)后的時頻譜值在時頻域內(nèi)進(jìn)行擠壓以及重排并將其置于瞬時頻率處，進(jìn)而提高信號的時頻聚焦性。

1.3 螞蟻追蹤技術(shù)

螞蟻追蹤技術(shù)是在蟻群算法的基礎(chǔ)上發(fā)展起來的，模仿螞蟻覓食所釋放的信息素指導(dǎo)其他螞蟻進(jìn)行追蹤，對斷裂信息有著非常強的識別能力。但是對于許多伴生褶皺以及次生裂縫，基于原始數(shù)據(jù)進(jìn)行的螞蟻追蹤無法進(jìn)行有效的刻畫。為避免破壞實際地震資料的縱橫向分辨率，所以實際地震資料處理對時頻分辨率的要求較高。傳統(tǒng)時頻分析方法，如短時傅立葉變換，S變換，WVD等，受自身窗函數(shù)的制約，時頻分辨率無法達(dá)到最佳，甚至還會產(chǎn)生許多虛假信息；尤其是WVD，由于WVD方法在計算時會產(chǎn)生意想不到的交叉項，導(dǎo)致其結(jié)果存在很大的誤差。而同步擠壓算法利用壓制海森堡不確定性以及交叉項影響上的優(yōu)勢使得結(jié)果具有很高的分辨率，這一點在實際應(yīng)用中被很好地驗證。

基于此，在三維地震斷裂解釋中，本文將螞蟻追蹤技術(shù)與同步擠壓改進(jìn)短時傅立葉變換在時頻分析中的優(yōu)勢相結(jié)合，利用不同頻率信息對地震數(shù)據(jù)斷裂系統(tǒng)進(jìn)行雕刻，主要操作步驟為(圖1)：

1)應(yīng)用同步擠壓改進(jìn)短時傅立葉變換對三維地震數(shù)據(jù)體做頻譜分解，分別獲得低頻、中頻、高頻的數(shù)據(jù)體。

2)對單頻數(shù)據(jù)體開展邊緣檢測、邊緣增強等工作，如：提取單頻數(shù)據(jù)體的方差體、混沌體、相干體、曲率屬性等；再對這些屬性體平滑增強，從而使數(shù)據(jù)體中存在的不連續(xù)處暴露出來。

3)再將螞蟻追蹤技術(shù)應(yīng)用于上述屬性體，得到單頻的螞蟻屬性體。

圖1 基于同步擠壓改進(jìn)短時傅立葉變換分頻螞蟻追蹤流程框架Fig.1 Improved ant tracing process framework of short-time Fourier transform based on synchronous extrusion

2 模型數(shù)據(jù)測試

為了驗證該方法的時頻分辨能力以及抗噪性，合成由x1、x2、x3、x4疊加的線性調(diào)頻信號：

(12)

應(yīng)用短時傅立葉變換、改進(jìn)短時傅立葉變換、同步擠壓改進(jìn)短時傅立葉變換分別對信號x進(jìn)行時頻變換，對比分析其時頻聚焦能力。

對比圖2中三幅時頻譜圖可以得出：改進(jìn)短時傅立葉變換的時頻分辨率較短時傅立葉變換得到了一定的提升，但海森堡不確定性以及交叉項干擾仍然是阻礙分辨率提高的重要因素；同步擠壓改進(jìn)短時傅立葉變換的時頻譜圖中，其結(jié)果不僅能量匯聚度強、時頻分辨率高，并且解決了短時傅立葉變換時窗不變的缺陷。在短時傅立葉變換以及改進(jìn)短時傅立葉變換的基礎(chǔ)上，計算得到的同步擠壓改進(jìn)短時傅立葉變換值是對改進(jìn)短時傅立葉變換的再處理過程，由于改進(jìn)短時傅立葉變換分辨率的局限性，使得信號在每一個時間點所對應(yīng)的瞬時頻率不止一個，因此產(chǎn)生了虛假帶寬的問題。同步擠壓改進(jìn)短時傅立葉變換在改進(jìn)短時傅立葉變換的基礎(chǔ)上，對虛假帶寬進(jìn)行擠壓，重新排列于信號的瞬時頻率附近，使得信號的時頻分辨率得到了極大的提升。

a—原始合成信號;b—短時傅立葉變換的時頻譜;c—改進(jìn)短時傅立葉變換的時頻譜;d—同步擠壓改進(jìn)短時傅立葉變換的時頻譜a—original synthesized signal;b—time frequency spectrum of short time Fourier transform;c—improved time frequency spectrum of short time Fourier transform;d—synchronous to improve time frequency spectrum of short time Fourier transform圖2 合成信號及利用不同分析方法得到的時頻譜Fig.2 Synthesized signal and time spectrum obtained by different analysis methods

對信號x進(jìn)行隨機(jī)加噪，從而驗證本文方法的抗噪性(圖3)。通過觀察時頻譜信息可以知道，在對信號添加隨機(jī)噪聲后，短時傅立葉變換、改進(jìn)短時傅立葉變換、同步擠壓改進(jìn)短時傅立葉變換的時頻能量都受到噪聲干擾，但是同步擠壓改進(jìn)短時傅立葉變換的有效信息保留得較好，整個時頻譜圖依舊比較干凈，受噪聲影響很小。因此，同步擠壓改進(jìn)短時傅立葉變換算法具有一定的抗噪能力。

傳統(tǒng)時頻分析方法是對每一道地震數(shù)據(jù)進(jìn)行處理，因此改進(jìn)短時傅立葉變換在縱向分辨率上得到了很大的提升，然而斷裂尺度等信息是和橫向分辨率相關(guān)的。同步擠壓改進(jìn)短時傅立葉變換在提升信號縱向分辨率的同時，開展了多道間的擠壓，即對多道地震記錄同時做擠壓以便兼顧地震數(shù)據(jù)的橫向分辨率。

通過對比分析圖4原始地震剖面及單頻剖面可以發(fā)現(xiàn)，對單道地震記錄循環(huán)做同步擠壓短時傅立葉變換所得到的單頻剖面，其縱向分辨率相較于改進(jìn)短時傅立葉變換得到了很大的提升，但是其橫向分辨率卻存在較大的弊端；因此開展了多道間的同步擠壓短時傅立葉變換，對比單道的同步擠壓短時傅立葉變換，多道的同步擠壓短時傅立葉變換無論其縱向分辨率還是橫向分辨率都得到了極大的提升。

本文所使用的數(shù)據(jù)屬于某工區(qū)三維地震數(shù)據(jù)體。頻譜分解技術(shù)首先需要確定頻率值的大小，頻率值直接決定著單頻剖面中的有效信息能否被準(zhǔn)確的描繪出來；因此，頻率值的選取在整個流程中占據(jù)十分重要的位置。本文首先對實際地震數(shù)據(jù)的井旁道做時頻變換處理，通過確定該地震資料的優(yōu)勢頻帶范圍，從而確定頻率值的大小，以突出所需要的地震數(shù)據(jù)有效信息。

首先對原始地震資料分別使用短時傅立葉變換以及同步擠壓改進(jìn)短時傅立葉變換提取井旁道做頻譜分析得到圖5頻譜。分析圖5中同步擠壓改進(jìn)短時傅立葉變換時頻能量分布可知，該地震數(shù)據(jù)體的優(yōu)勢頻帶范圍位于25～60 Hz。因此，選取25、55 Hz作為單頻數(shù)據(jù)體的頻率值。

a—原始加噪信號;b—短時傅立葉變換的時頻譜;c—改進(jìn)短時傅立葉變換的時頻譜;d—同步擠壓改進(jìn)短時傅立葉變換的時頻譜a—original noise signal;b—time frequency spectrum of short time Fourier transform;c—improved time frequency spectrum of short time Fourier transform;d—synchronous to improve time frequency spectrum of short time Fourier transform圖3 加噪信號及利用不同分析方法得到的時頻譜Fig.3 Noise-added signal and time spectrum obtained by different analysis methods

a—原始地震記錄;b—改進(jìn)短時傅立葉變換單頻剖面;c—單道循環(huán)擠壓改進(jìn)短時傅立葉變換單頻頻譜;d—多道循環(huán)同步擠壓改進(jìn)短時傅立葉變換單頻剖面a—original seismic signal;b—improved short time Fourier transform single frequency profile;c—single channel cycle squeeze to improve short time Fourier transform single frequency profile;d—multi channel cycle squeeze to improve short time Fourier transform single frequency profile圖4 原始地震記錄及單頻剖面Fig.4 Original seismic record and single frequency profile

a—井旁單道地震記錄;b—短時傅立葉變換得到的時頻譜;c—同步擠壓改進(jìn)短時傅立葉變換得到的時頻譜a—single channel seismic record near the well;b—time frequency obtained by short time Fourier transform;c—synchronous squeezing to improve the time spectrum obtained by short time Fourier transform圖5 井旁地震道及利用不同分析方法得到的時頻譜Fig.5 Seismic trajectory by the well and the time spectrum obtained by different analysis methods

3 應(yīng)用實例

基于以上理論基礎(chǔ)以及模型分析可以知道，同步擠壓短時傅立葉變換相較于傳統(tǒng)的時頻分析算法具有很高的時頻聚焦能力，可以壓制信號時頻變換過程中產(chǎn)生的虛假帶寬問題，將信號的真實瞬時頻率凸顯出來。因此應(yīng)用同步擠壓短時傅立葉變換所得到的單頻數(shù)據(jù)體可以刻畫出傳統(tǒng)時頻分析算法所無法描繪的斷裂系統(tǒng)信息。基于上述知識，本文將同步擠壓改進(jìn)短時傅立葉變換應(yīng)用于實際地震資料，分析得出同步擠壓改進(jìn)短時傅立葉變換在處理實際地震數(shù)據(jù)時的準(zhǔn)確性。分別對三維地震數(shù)據(jù)體提取25、55 Hz的單頻數(shù)據(jù)體，再對單頻數(shù)據(jù)體做邊緣增強，計算得到單頻螞蟻體。

圖6為原始振幅沿層切片，圖7為全頻帶沿層螞蟻體切片，圖8、圖9分別為25、55 Hz單頻螞蟻體沿層切片。

在單頻螞蟻體中，許多微小斷裂信息被很好地雕刻出來，25 Hz單頻螞蟻體與全頻帶螞蟻體都對于大斷裂信息進(jìn)行了詳細(xì)的描繪，而55 Hz相對高頻的螞蟻體對許多微小斷裂信息做了很清晰的刻畫。

全頻帶螞蟻體可以對大斷層、中斷層進(jìn)行了非常清晰地描繪；單頻螞蟻體中，隨著頻率升高，剖面中一些細(xì)微的斷裂信息被凸顯出來；但是隨著頻率的增加，噪聲信息也會被放大，因此，在高頻剖面中會出現(xiàn)一些采集腳印等干擾信息；但是單頻螞蟻體技術(shù)能夠識別出全頻帶螞蟻體無法刻畫的小斷裂，具有很高的實用價值。

結(jié)合上述基于同步擠壓短時傅立葉變換所得到的中、高頻螞蟻體對微小斷裂信息有著較強的刻畫能力，低頻螞蟻體可以對大、中斷層進(jìn)行詳細(xì)雕刻。

圖6 振幅沿層切片F(xiàn)ig.6 Amplitude slice along the layer

圖7 全頻帶沿層螞蟻體切片F(xiàn)ig.7 Full-band ant slice along the layer

圖8 25 Hz沿層螞蟻體切片F(xiàn)ig.8 25 Hz ant slice along the layer

圖9 55 Hz沿層螞蟻體切片F(xiàn)ig.9 55 Hz ant slice along the layer

4 結(jié)論及建議

進(jìn)行分頻屬性檢測，時頻分析是核心，傳統(tǒng)的時頻分析方法在縱橫分辨率上已無法達(dá)到人們對精度的要求，尤其是縱向分辨率；而同步擠壓變換算法可以很好地提高地震數(shù)據(jù)的縱向分辨率；因此，在同步擠壓變換算法上進(jìn)行分頻屬性提取，會達(dá)到比傳統(tǒng)時頻分析算法更好的效果。

本文應(yīng)用同步擠壓改進(jìn)短時傅立葉變換對地震數(shù)據(jù)進(jìn)行頻譜分解，有效地壓制噪聲的影響，對后續(xù)的螞蟻體計算提供了高信噪比數(shù)據(jù)體；并且高頻螞蟻數(shù)據(jù)體可以刻畫出全頻帶螞蟻體無法刻畫的微小斷裂以及伴生褶皺等信息；尤其在裂縫發(fā)育帶，該方法可以較為準(zhǔn)確地雕刻其分布區(qū)域以及特征；數(shù)據(jù)體融合重構(gòu)技術(shù)能夠匯集低頻、中頻、高頻等裂縫信息，對斷層可以更進(jìn)一步的清晰描繪?；谕綌D壓改進(jìn)短時傅立葉變換的分頻螞蟻追蹤方法能夠精確地描繪出區(qū)域斷裂系統(tǒng)信息，為后續(xù)儲層產(chǎn)量預(yù)測以及井位的確定提供現(xiàn)實依據(jù)。后期為了得到更加精細(xì)的斷裂特征，可以參照RGB技術(shù)的思想，將低頻、中頻、高頻螞蟻數(shù)據(jù)體進(jìn)行融合重構(gòu)，得到一個新的螞蟻數(shù)據(jù)體；該螞蟻數(shù)據(jù)體將同時具備低頻的信息，也將具備中頻、高頻的裂縫信息；它可以很清晰地描繪出斷層的發(fā)育走向、微小裂縫以及伴生褶皺等綜合信息。