劉一灼, 汪勇*, 桂志先
(1.長江大學地球物理與石油資源學院, 武漢 430100; 2.油氣資源與勘探技術教育部重點實驗室(長江大學), 武漢 430100)
傅里葉變換是一種線性變換,針對實際地震波信號中更多的非線性問題,傅里葉變換明顯作用有限。在這一問題的研究中,1946年Gabor首先提出了短時傅里葉變換(Short-time Fourier transform,STFT),可以在局部范圍內(nèi)對非平穩(wěn)信號進行分析,經(jīng)過后續(xù)研究表明該方法在實際測試中存在分辨率整體不高的問題[1-3]。在此基礎之上,1986年由Mallat提出的連續(xù)小波變換(continuous wavelet transform,CWT),其本質(zhì)是一種可調(diào)節(jié)窗口的傅里葉變換,在信號仿真過程中初端分辨率得到有效提高,但末端分辨率發(fā)散的問題依然存在[4-6]。而近些年提出的Chirplet變換(Chirplet transform,CT)是對短時傅里葉變換和連續(xù)小波變換的延申說明,由于其核函數(shù)不僅具有短時傅里葉變化和連續(xù)小波變換的時移、伸縮、時間和頻率上的4種線性調(diào)頻變換,還具備對核函數(shù)調(diào)制的特性,因此Chirplet變換對待分析信號在時頻域分析中分辨率提升明顯[7]。國外學者Zhu等[8]利用高分辨率勘探技術對Chirplet變換中的物性參數(shù)進行了對比分析,得出了不同取值下的Chirplet變換時頻譜特征;Pukhova等[9]則是對Chirplet變換的波形上弦波和下弦波的展布特征進行了分析,得出了Chirplet變換在處理復雜波形的優(yōu)勢所在,國內(nèi)學者田亞軍等[10]、聶偉東等[11]對STFT、CWT以及CT幾種時頻分析方法在沉積旋回識別上的差異情況進行了分析,頻率隨時間的變化關系得到準確定量描述:頻率隨時間依次降低,代表地層層間厚度依次變厚,反之則代表地層厚度依次變薄。邱劍鋒等[12]則是通過修改Chirplet變換中的核函數(shù)控制因子的取值,提高時間域的頻次以此來增強Chirplet變換的整體分辨率。
在前人研究的基礎之上,依據(jù)短時傅里葉變換、連續(xù)小波變換以及Chirplet變換的原理,現(xiàn)分別利用單分量和多分量仿真信號對3種變換進行測試,分析對比不同時頻分析方法下的待分析信號頻譜特征,同時再利用多分量仿真信號對原始Chirplet、所改進的2次Chirplet以及1.5次Chirplet進行測試,以此來突出在復雜時域波形中Chirplet變換的優(yōu)勢所在。隨后利用改進后的不同頻次的Chirplet變換對經(jīng)過雷克子波與正負反射系數(shù)褶積后的混合旋回仿真信號進行測試,分析不同頻次下的CT刻畫沉積旋回薄互層準確位置的能力。同時針對塔里木地區(qū)工區(qū)一段特殊地震道下的地震信號進行分析,該道地震信號分布清晰且有一定規(guī)律,在該段地震信號中包含四組尺度不同的儲層位,利用傳統(tǒng)線性時頻方法STFT、CWT以及原始CT和所改進后的變頻次CT分別對該地震信號進行頻譜分析,突出改進后的修改控制因子的Chirplet變換在應對實際地震勘探中的時頻譜能量團集中度以及壓制頻散效果的優(yōu)勢,因此在后續(xù)地震信號分析中該方法應該受到重視與發(fā)展,以期為實際儲層地震勘探提供理論依據(jù)支撐。
短時傅里葉變換由Gabor最先提出,針對傳統(tǒng)傅里葉變換對非平穩(wěn)信號局部信號區(qū)域時頻分布無法進行分析的缺陷,短時傅里葉變換得以對其進行處理。由于所選擇的信號相較于完整信號長度較短,因此該變換被稱為短時傅里葉變換(STFT)。短時傅里葉變換的基本原理是:首先通過窗函數(shù)來截取給定平穩(wěn)信號;再利用傅里葉變換處理窗內(nèi)的信號,這部分的目標是準確定位測試局部信號時間內(nèi)的頻率變化;最后追蹤信號定向時移傳函數(shù),即可準確獲得這部分局域頻率隨時間的變化關系。簡言之,變換過程首先對給定函數(shù)與窗函數(shù)進行褶積處理,對處理后的多元函數(shù)進行傅里葉變換,隨即時移窗函數(shù)得到所需變換后的結果,即是短時傅里葉變換整個處理過程[13]。
已知信號z(t),其時域短時傅里葉變換為
(1)
已知信號z(t),其頻域短時傅里葉變換為
(2)
式中:z(t)為源信號;g(t)為窗函數(shù);t為時間;f為頻率。
利用式(1)、式(2)可推導得出短時傅里葉變換的逆變換為
(3)
通過逆變換分析,短時傅里葉變換主要是利用參數(shù)τ來對信號進行描述時間、頻率的變化關系,即通過小段信號來對整體信號進行描述。短時傅里葉變換的時窗控制范圍為[u-η,u+η],其中參數(shù)η用來控制時頻分辨率大小,因此短時傅里葉變換受限于固有的窗函數(shù)被選定,本身變換的矩形時頻原子無法進行定量變化。
連續(xù)小波變換1986年由Mallat所提出,與短時傅里葉變換的全區(qū)域變換所不同的是,該變換是在時頻域局部限定區(qū)域進行的,因此,可以針對所研究的局域信號進行這部分信息的提取。連續(xù)小波變換的基本原理是:在傅里葉變換中,利用小波基函數(shù)的平移、伸縮等變化關系替代傳統(tǒng)伸縮函數(shù),替代后的時頻變換就是連續(xù)小波變換。簡單來說,連續(xù)小波變換實質(zhì)是在傳統(tǒng)傅里葉變換下其本身窗口可調(diào)節(jié),但由于其小波窗內(nèi)的信號與短時傅里葉讓變換要求一致,必須是平穩(wěn)信號,因而仍然存在傳統(tǒng)傅里葉變換以及短時傅里葉讓變換對信號的分析限制[14]。
已知信號z(t),其連續(xù)小波變換為
(4)
通過式(4)可知,連續(xù)小波變換主要是通過時移參數(shù)τ、尺度參數(shù)a來對整個變換進行控制和調(diào)節(jié)的,利用這兩個主要參數(shù),得到小波族函數(shù)為
(5)
通過小波族函數(shù)進行分析,該變換主要是由尺度因子a來控制基本小波的大小,參數(shù)a增大,時域帶寬變小,頻率增大。由于該變換時間隨頻率變化而變化,此類特性在應對二維信號分析中效果好于短時傅里葉變換,但該變換沒有確切的限定值域限制時頻域的增加或減少。
Chirp信號最早由Winkler于1962年提出,該信號由于是非平穩(wěn)信號,因此又被命名為線性調(diào)頻信號。而以線性調(diào)頻信號為研究中心的Chirplet變換是一種新的時頻分析方法,該變換本質(zhì)上來講是對短時傅里葉變換與連續(xù)小波變換兩種變換的擬合,所以說可以看作是兩種變換的推廣[3,9]。由于Chirplet變換中所研究的對象是利用Chirp信號函數(shù)族與待分析信號的褶積計算得出,在褶積運算中同時運用了時間、頻率以及單位尺度變化,這也再次證實了Chirplet變換是在短時傅里葉變換與連續(xù)小波變換的基礎之上的一種線性變換的說法。Chirplet變換的基本原理是:利用對窗函數(shù)g(t)做時間上的平移和頻率調(diào)制、對基本小波作平移和伸縮以及矩形時頻原子在斜方向的拉伸與旋轉(zhuǎn)變化這5種變化關系來對給定待分析信號進行的一種線性調(diào)頻變換方式[12]。
給定一個函數(shù)s(t),Chirplet變換表達式為
(6)
式(6)中:s(τ)為給定函數(shù);h(t)為族函數(shù);參量t、f、a、c、d分別為時間、頻率、伸縮量、時域上的線性調(diào)頻、頻域上的線形調(diào)頻這5種時頻變化關系。
通過對該變換的原理進行分析,該變換主要是通過將待分析信號投影到族函數(shù)中實現(xiàn)的,而族函數(shù)一般是通過修正原始窗函數(shù)得到的,窗函數(shù)一般都是選擇壓制頻散效果較好的高斯窗函數(shù),其表達式為
(7)
整個線性調(diào)頻變換過程主要包括在時間上的線性調(diào)頻以及在頻率上的線性調(diào)頻兩部分,時間上的線性調(diào)頻主要是利用線性調(diào)頻信號chirp1與高斯窗函數(shù)進行乘積得到,表示為
(8)
而頻率上的線性調(diào)頻主要是利用另一組線性調(diào)頻信號chirp2與高斯窗函數(shù)進行乘積得到,表示為
(9)
因此,Chirplet變換的線性調(diào)頻過程即核函數(shù)表示可以利用式(8)與式(9)的乘積后引入分析信號得到,即
(10)
Chirplet變換對于解析信號的局部特性充分考慮,在重建非線性瞬時頻率中去除誤差效果很好。
通過后續(xù)相關仿真結果可以得出,參數(shù)d在整個時頻變換中是不進行相關計算的,因而Chirplet變換的線性調(diào)頻表達式可簡化為
(11)
由于Chirplet變換的主要變化就是在線性調(diào)頻信號chirp也就是核函數(shù)中進行的,因此通過修改核函數(shù)增加核函數(shù)中的參數(shù)因子可以對整個變換進行有效改進
(12)
(13)
為了分析不同變換方法在對給定不同類型的待分析信號下的時頻譜分辨能力以及壓制頻散效果,依次仿真了3類信號進行測試對比。
信號1:一次單分量仿真線性調(diào)頻信號,s(t)=sin(2πut2)。
信號2:二次單分量信號,s(t)=sin(πut3)。
信號3:多分量線性調(diào)頻信號,x(t)=cos[2π×(2t-t2+t1.5)]+cos[2π(3t-t2+3t1.5)]。
其中u=300,輸入其他參量及參數(shù)選取為:采樣頻率fs=1 000 Hz、采樣時間t=1 s、采樣時間間隔dt=0.001 s。3種信號的時頻域展布如圖1所示。
圖1 3種信號時域波形Fig.1 Three kinds of signal time-domain waveforms
對3種信號分別做短時傅里葉變換(STFT)、連續(xù)小波變換(CWT)以及Chirplet變換(CT)的時頻分布對比如圖2~圖4所示。
圖2 信號1的3種變換時頻對比圖Fig.2 Comparison of three transformations of signal 1 in time frequency
圖3 信號2的3種變換時頻對比圖Fig.3 Comparison of three transformations of signal 2 in time frequency
圖4 信號3的3種變換時頻對比圖Fig.4 Comparison of three transformations of signal 3 in time frequency
從對給定單分量信號進行仿真測試結果來看,短時傅里葉變換(STFT)由于受到固定窗函數(shù)的制約,無法對本身時頻原子進行有效變換,因此該變換下的待分析信號在時域分析中分辨率整體較低。連續(xù)小波變換(CWT)在時頻分析中有高頻部分頻散嚴重這類情況,其有效分辨率也不是太高。而Chirplet變換(CT)更能對單分量信號的時間與頻率相對變化關系進行準確描述,相比于短時傅里葉變換(STFT)以及連續(xù)小波變換(CWT),Chirplet變換(CT)在時域分析下能量分辨率最高。
對給定二次單分量信號進行短時傅里葉變換(STFT)、連續(xù)小波變換(CWT)、Chirplet變換(CT)可以發(fā)現(xiàn),產(chǎn)生的差異變化基本與一次單分量信號相一致。由于知道地層層間走向通常是彎曲而非平穩(wěn)的,因此選擇地震信號要盡可能貼近于實際地質(zhì)構造,二次分量線性調(diào)頻信號相較于一次信號能更好地模擬實際地震信號時頻域分布變化關系,二次信號時頻變化帶有角度彎曲這一特性,更符合地震信號在實際地下介質(zhì)中的傳播,更有利于進行此類數(shù)值模擬的研究。
從圖4可以看出,在符合對單分量信號分析的特性同時,在對給定多分量信號分析中可以更加清楚發(fā)現(xiàn),短時傅里葉變換下的時頻域頻散現(xiàn)象最為嚴重,連續(xù)小波變換其次,Chirplet變換壓制頻散能力相對最好,相對而言能更好地區(qū)分兩個合成分量。所對應的Chirplet變換時頻域分辨率最佳,連續(xù)小波變換次之,短時傅里葉變換分辨率最低,也就是說CT變換相較于其他兩種時頻譜方法更能對兩組頻率相近的信號進行有效的分辨。
為了更直觀反映出修改核函數(shù)后的二次Chirplet變換、1.5次Chirplet變換與原始Chirplet變換的差異,因此對給定多分量線性調(diào)頻信號分別做1次CT、2次CT和1.5次CT,其時頻譜特征如圖5所示。
圖5 信號3的不同頻次的Chirplet變換時頻對比圖Fig.5 Comparison of Chirplet transform time frequency for signal 3 at different frequencies
從圖5可以看出,在對給定多分量信號分析中可以清楚發(fā)現(xiàn),3種不同頻次的Chirplet變換在對兩個分頻信號進行分析中,原始Chirplet變換下的時頻域頻散現(xiàn)象最為嚴重,二次Chirplet變換其次,1.5次Chirplet變換壓制頻散能力相對最好,能更好地區(qū)分兩個合成分量。所對應的1.5次CT變換時頻域分辨率最佳,二次CT次之,原始CT分辨率最低,也就是說1.5次CT變換相較于原始CT變換以及二次CT變換更能對兩組頻率相近的信號進行有效的分辨[10-11]。
對沉積旋回特征分析有助于更精準定位層序地層中的薄互層層間厚度變化的研究。沉積旋回信號的頻率隨著地層層間的厚度時刻發(fā)生變化,通常表現(xiàn)為:薄互層變厚,沉積旋回信號頻率逐漸變小,反之薄互層變薄,沉積旋回信號頻率隨即變大。換言之,對沉積旋回進行準確分析是確定層序地層層間厚度的關鍵前提,而準確掌握地層的層間厚度變化對尋找油氣藏意義顯著。
通常采用傳統(tǒng)子波與反射系數(shù)序列進行褶積,疊加合成混合型旋回信號。其中參量取值為:采樣頻率fs=1 000 Hz,采樣時間t=0.7 s,采樣時間間隔ds=0.07 s,這里選擇主頻f=50 Hz的雷克子波,反射系數(shù)序列選擇正負相反的反射序列進行褶積處理,其中反射系數(shù)周期T=20個長度單元,每段間隔為1個長度單元,雷克子波、反射系數(shù)序列以及褶積模型如圖6所示。
圖6 構建褶積模型圖Fig.6 Construct the convolution model
對如圖6所示得到的合成旋回信號記錄,依次進行原始Chirplet變換、2次Chirplet變換、1.5次Chirplet變換,所得不同頻次的CT時頻譜如圖7所示。
圖7 混合旋回信號CT變換時頻分布Fig.7 Time-frequency distribution of the original CT transformation of the mixed rotary signal
從以上分析可以看出,由于不同頻次下的菱形時頻原子表征物性參量數(shù)量級各不相同,這里1.5次Chirplet變換物性參數(shù)最為豐富,因此1.5次Chirplet變換相較于原始Chirplet變換和2次Chirplet變換,能量集中度最高,分辨率最強,同時沉積旋回信號旁瓣頻散情況最低,更能準確觀測出信號的頻率隨時間的變化趨勢。通過對核函數(shù)中的曲率等物性參數(shù)在一定范圍內(nèi)進行調(diào)節(jié)可以對地層旋回變化趨勢進行準確刻畫,有利于地震勘探中對層間屬性中的每一段薄層進行準確預測實際分布位置,為后續(xù)尋找油氣提供幫助[15]。
為了表明Chirplet變換及其改進的變頻次的Chirplet變換在對沉積旋回儲層檢測的有效性,對塔里木1個實際工區(qū)的地震信號進行處理分析,在反復對比了該地區(qū)各道地震數(shù)據(jù)的情況后,某一地震道信號分布簡單且對比度明顯,該道地震道范圍為8 000 m,整體振幅1 m,主要包含4個厚度大小不同的儲層,儲層厚度大致依次為800、400、200、100 m,該道實際地震道信號的時域波形如圖8所示。
圖8 工區(qū)地震道時域波形Fig.8 Time domain waveform of seismic channel in the work area
對該地震道時域波形依次做短時傅里葉變換、連續(xù)小波變換、原始Chirplet變換、2次Chirplet變換以及1.5次Chirplet變換,其不同方法的時頻譜特征如圖9所示。
通過對比在實際工區(qū)下的地震道信號做時頻分析可以得出:短時傅里葉變換時頻譜分辨率最低,能量保持相對最不集中,頻散最為嚴重;連續(xù)小波變換時頻譜可以對大尺度儲層進行刻畫,小尺度儲層分辨率變低,出現(xiàn)頻散現(xiàn)象和交叉項的干擾;原始Chirplet變換相較于兩種線性變換分辨率提高明顯,但對尺寸給臨界刻畫不夠清晰,會出現(xiàn)對儲層位置的誤判;二次Chirplet變換在原始變換基礎上,能量得到更加集中,儲層分界面刻畫更加清晰,只是在應對小尺度儲層分辨率不是很高,會出現(xiàn)一定的頻散問題;而1.5次Chirplet變換其分辨率最高,能量更加集中,無明顯頻散問題,儲層邊界刻畫清晰,對于小尺度儲層同樣可以清晰識別。
(1)分別對3種時頻變換短時傅里葉變換、連續(xù)小波變換、Chirplet變換的基本原理進行了對照總結,并利用給定單分量、多分量信號仿真對3種變換進行控制變量分析,結果表明:Chirplet變換由于其特定的時頻原子物性豐富度的優(yōu)勢,該變換相較于短時傅里葉變換、連續(xù)小波變換在時頻分布下壓制頻散效果最佳,同時分辨率相對最高。
(2)分別對不同頻次的Chirplet變換調(diào)頻實際過程進行分析,對各自調(diào)頻的核函數(shù)物性參數(shù)的差異進行對比描述,并再次利用給定多分量信號仿真對不同頻次的Chirplet變換進行控制變量分析,結果表明:1.5次Chirplet變換由于對其調(diào)頻核函數(shù)引入了曲率參量,時頻分布下頻率隨時間的變化程度得到明顯增強,在對多分量信號進行測試下,1.5次CT變換更能對頻率相近的兩組信號進行有利準確地劃分,在不同頻次下的時頻分布關系1.5次CT變換分辨率最為出色。
(3)在實際應用中,利用雷克子波與反射系數(shù)序列進行褶積,得到混合沉積旋回信號,在前人僅針對短時傅里葉變換、連續(xù)小波變換、Chirplet變換下的混合旋回信號時頻分布特征對比分析下,對3種不同頻次的Chirplet變換進行定量對比分析,結果表明:1.5次Chirplet變換在分析沉積旋回信號具有比傳統(tǒng)CT變換、2次CT變換分辨率高的優(yōu)勢,同時在旁瓣壓制頻散能力最強,時頻域信號信息豐富度最全面。在實際對地震道信號進行時頻譜分析中,改進的Chirplet變換分辨率更高,能量更為集中,壓制頻散效果明顯,對儲層尤其是小尺度尺寸給刻畫清晰,因此1.5次CT變換可以更好地利用時頻分布關系來確定薄層厚度變化趨勢,為后續(xù)尋找油氣層提供理論支撐。