基于變窗函數(shù)譜重排的譜分解技術

崔震，曹思遠，劉偉，馬媛媛（中國石油大學a.油氣資源與探測國家重點實驗室；b.CNPC物探重點實驗室，北京102249）

基于變窗函數(shù)譜重排的譜分解技術

崔震，曹思遠，劉偉，馬媛媛
（中國石油大學a.油氣資源與探測國家重點實驗室；b.CNPC物探重點實驗室，北京102249）

譜分解是將地震信號由時間域變換到時間—頻率域，然后再對其局部時頻屬性和瞬時頻率特性進行分析的技術。常規(guī)譜分解方法受Heisenberg準則限制，不能同時兼顧時間分辨率和頻率分辨率，時頻譜重排有效規(guī)避了Heisenberg準則的限制，同時提高了時頻譜的聚集性，且對噪音等干擾聚焦。變窗函數(shù)譜重排算法是在常規(guī)時頻譜重排的基礎上，利用Hermite窗的正交性和可遞推性，得到多個變化窗口的譜重排結果，然后利用最小二乘估計原理，對多個譜重排剖面進行加權近似組合。合成記錄和實際例子均表明該方法在保證了時頻譜高聚焦性的同時，分散了噪音干擾，能準確歸位時頻能量的真實位置，對低頻陰影有很好的顯示效果。

譜分解；Hermite窗函數(shù)；變窗口；譜重排；加權組合

譜分解技術在地震數(shù)據(jù)處理中應用十分廣泛。文獻［1］發(fā)現(xiàn)在油氣儲集層下方存在低頻伴影現(xiàn)象；文獻［2］總結了產生低頻陰影區(qū)的至少10個原因；文獻［3］將短時傅里葉變換引入油氣勘探領域，對頻率的橫向變化進行定量描述，產生了最初的地震譜分解技術。隨著時頻分析技術的不斷發(fā)展，譜分解也向高精度、高保真度方向發(fā)展，目前已經廣泛應用于油氣檢測、河道檢測、斷裂解釋、AVO分析等方面[4-6]。

時頻分析技術是譜分解的基本手段，文獻［7］提出以高斯函數(shù)為線函數(shù)構建基函數(shù)，將一維時間信號映射為二維時頻變信號。文獻［8］提出短時傅里葉變換（Short Time Fourier Transform，STFT），使用固定時窗分析局部時間的頻率成分，但是無法兼顧時間和頻率的分辨率。伽柏變換是使用高斯窗的一種特殊的短時傅里葉變換，由于高斯窗的特性，使得它具有短時傅里葉變換的所有性質，并且在一定程度上兼顧了時間和頻率分辨率[9]。文獻［10］提出的連續(xù)小波變換對于非平穩(wěn)信號的刻畫有了很大的改善，其窗函數(shù)的寬度隨頻率變化，在高頻時使用窄窗，在低頻時使用寬窗，這樣就可以用不同的分辨率來分析信號，相對于短時傅里葉變換，小波變換具有更好的能量聚集性和多分辨率特征。文獻［11］S變換時頻分析方法，其基于滑動高斯窗函數(shù)，同時具有傅里葉變換和小波變換的優(yōu)點。

文獻［12］提出了時頻譜重排的思想，主要應用在非平穩(wěn)信號的分析中。文獻［13］重新解釋了譜重排的實際含義，把譜重排看成一個時頻分布的反平滑過程。時頻重排譜對信號具有更好的局部刻畫能力，能夠把原始譜中的每一個點的時頻能量集中到信號的瞬時頻率和延遲時的位置。前人的研究主要集中在信號分析領域方面，在地球物理領域的研究相對較少。時頻譜重排能夠生成聚焦的時頻分布，同時也對隨機噪音能量聚焦，文獻［14］提出多窗口重排方法，對于信號的抖動具有更好的穩(wěn)定性。本文將變窗口思想與時頻譜重排相結合，選取正交Hermite函數(shù)作為窗函數(shù)，得到了時頻分辨率較高的時頻分布，同時壓制了噪音的能量。

1　方法原理

1.1時頻重排原理

時頻譜重排是在原有頻譜的基礎上進行能量重新分配，提高譜圖的時頻聚集性，使得時頻定位準確度提高。由于伽柏譜圖在一定程度上兼顧時間分辨率和頻率分辨率，可以為重排提供較好的基礎。伽柏時窗內的時頻振幅譜可以定義為地震信號振幅譜和高斯窗函數(shù)的振幅譜的乘積：

高斯窗h(t)=exp(-t2/2σ)/ 2πσ，固定σ后，根據(jù)Heisenberg測不準準則，窗口大小固定為，則由（1）式和（2）式可以看出，伽柏譜即是用的時頻窗，分別在時間和頻率方向對真實點譜的平滑。

文獻［15］根據(jù)時頻譜窗函數(shù)與真實點譜的關系，提出譜重排算法，其目的是將伽柏譜中每一個時頻點處的能量進行移動，使其聚集在二維高斯窗內，譜能量的重心處即為真實譜坐標點，從而規(guī)避了Heisenberg測不準準則，提高了時頻分辨率。由文獻［15］提出的譜能量重排算法可以得到：

任意點(t',f')的值是所有重排到這一點的點(t,f)值的和：

1.2多窗口統(tǒng)計方法

在計算中，時頻分析實際上是對無限長平穩(wěn)信號序列的截斷，這種截斷效應往往使譜圖分辨率降低。針對這一現(xiàn)象，文獻［14］在發(fā)現(xiàn)多個正交窗口的組合可以削弱這種截斷效應的基礎上提出了多窗口方法，利用多個正交窗口{hk(t)，k∈N}獲得獨立的近似功率譜估計，綜合統(tǒng)計這些近似估計，最終得到平穩(wěn)信號功率譜估計，即

文獻［16］則將其推廣到時頻域：

1.3改進的時頻重排

由于譜重排是譜圖平滑的逆過程，相比于普通的譜圖，譜重排對平穩(wěn)譜的截斷效果更加明顯。從（6）式的分析可以看出，多窗口加權可以壓制畸變，重新得到穩(wěn)定的光滑譜，與原譜圖相比具有更高的聚焦性。

1.3.1Hermite窗函數(shù)

為了得到平穩(wěn)信號譜估計，要求多窗口的選取為正交函數(shù)，其時頻分布近似橢圓或圓形分布，由于Hermite函數(shù)是正交函數(shù)，其時頻聚集性好，能夠較好地刻畫大多數(shù)非平穩(wěn)信號[17]，而且在計算過程中可以通過遞推方式減少計算量，因此選擇Hermite函數(shù)作為窗口函數(shù)進行時頻譜重排。

k階Hermite函數(shù)表示為

由（3）式可以看出，計算過程中需要對th（t）和dh（t）/dt進行多次迭代，為此，對dh（t）/dt化簡得到遞推關系式以減少計算量：

1.3.2變窗統(tǒng)計組合

利用不同階數(shù)Hermite窗口函數(shù)按照（3）式進行時頻譜重排，可以得到多個窗口聚焦時頻譜，即變窗口譜重排（Multiple Window Reassigned Spectrum，MWRS），對這些譜進行統(tǒng)計組合，即

估算系數(shù)dk利用最小二乘法計算：

離散化得到矩陣形式：

對（12）式中dk求導，令

求解得

這樣得到加權估計系數(shù)dk，對P的最小二乘估計為

其物理意義是利用反演的方法，尋找多窗口譜圖重排的最優(yōu)組合。

2　模型測試

圖1a中的信號由50Hz雷克子波、50Hz90°相移雷克子波、20Hz 180°相移雷克子波和2個時間相近的30Hz雷克子波構成。分別對該信號進行伽柏變換，伽柏譜重排（圖1c）和變窗口譜重排，可以看出，與伽柏譜（圖1b）相比，伽柏重排譜（圖1c）時頻聚集性有較大提升，時頻分辨率都有所提高，然而仔細查看不難發(fā)現(xiàn)，伽柏重排譜中產生了多個亮點（譜放大區(qū)域），這是譜重排截斷效應引起的，而且對于距離較近的子波其時頻譜會出現(xiàn)相互干涉，無法準確指示頻率的真實位置（紅色區(qū)域）。變窗口譜重排（圖1d）在保證較高時頻聚集性的同時克服了多個極值的問題，能準確指示頻率的真實位置，且能很好地區(qū)分距離相近的子波。

圖2是在圖1信號的基礎上加入高斯噪音之后進行時頻分析的結果，不難看出，伽柏譜重排對一些隨機噪音能量進行了聚焦，變窗譜重排在保證較高的時頻分辨率的同時，對噪音能量也有很好的壓制作用。圖3是在實際地震資料中抽取單道數(shù)據(jù)的時頻分析結果，變窗口重排譜具有更高的時頻分辨率，且對噪聲也有較好的壓制，低頻陰影現(xiàn)象突出明顯。

圖1　合成地震記錄時頻分析結果

圖2　含噪合成地震記錄時頻分析結果

3　實例

由于大地濾波和吸收衰減作用，地震波頻譜呈現(xiàn)非平穩(wěn)性，低頻陰影現(xiàn)象就是非平穩(wěn)性之一。低頻陰影是在含油氣層下方出現(xiàn)的低頻強能量區(qū)域，該區(qū)域高頻能量則比較弱[18]。這一特性常被用來探測含氣儲集層，在低頻剖面上含氣儲集層下方出現(xiàn)高亮能量團，隨著譜分解剖面頻率增加，亮度削弱直到消失。

圖4為新疆某地區(qū)地震疊后剖面，聲波時差測井曲線表明，在該位置速度保持低位，為高儲含氣層，與鉆井結果相吻合。利用常規(guī)伽柏變換和變窗譜重排分別對整個剖面進行譜分解，分別截取伽柏譜數(shù)據(jù)體和變窗重排譜數(shù)據(jù)體的20 Hz頻率切片如圖5a和圖5b所示。由圖5可見，由于伽柏譜時頻分辨率過低，能量團上下延時接近0.1 s，影響了目的層的精確定位，而變窗重排譜時頻分辨率較高，能夠準確刻畫目的層位置，并且可以清晰判別能量展布。

圖3　單道實際記錄時頻分析結果

圖4　實際含氣地震剖面（紅色區(qū)域為含氣儲集層）

圖5　頻率切片對比

圖5c和圖5d分別為截取伽柏譜數(shù)據(jù)體和變窗重排譜數(shù)據(jù)體的35Hz頻率切片。從圖5可以看出，隨著頻率增加，能量逐漸減弱，變窗重排譜相較伽柏譜時頻分辨率更高，能夠更加準確刻畫目的層位置。此外，變窗重排譜對低頻陰影現(xiàn)象的敏感度要明顯高于伽柏譜，是一種實用的地震記錄譜分解方法。

4　結論

本文在常規(guī)伽柏時頻分析方法的基礎上，討論了伽柏譜重排的優(yōu)點，同時指出了其對于信號能量收斂聚焦的不穩(wěn)定性，結合多窗口壓制不穩(wěn)定震蕩的優(yōu)點，提出了變窗口譜重排法，該方法利用多階Hermite函數(shù)作為窗函數(shù)進行時頻譜重排，通過最小二乘原理對其加權組合得到穩(wěn)定的高分辨率時頻譜。通過對合成數(shù)據(jù)和實際資料處理，可以看出變窗口譜重排相對于常規(guī)的伽柏變換和伽柏譜重排具有更高的時頻分辨率，同時對噪聲也有一定的壓制效果。

符號注釋

［1］TanerM T，Koehler F，SheriffR E.Complex seismic trace analysis［J］.Geophysics，1979，44（6）：1 041-1 063.

［2］Ebrom D.The low?frequency gas shadow on seismic sections［J］.The Leading Edge，2004，23（8）：772.

［3］Partyka G，Gridley J，Lopez J.Interpretational applications ofspec?tral decomposition in reservoir characterization［J］.The Leading Edge，1999，18（3）:353-360.

［4］MarfurtK J，Kirlin R L.Narrow?band spectralanalysisand thin?bed tuning［J］.Geophysics，2001，66（4）：1 274-1 283.

［5］Prasad R，DawwasM，Tanoli S.Detecting channelsands using spec?traldecomposition on 3?D seismic data:a case study［R］.EAGEDe?tective Stories Behind Prospect Generation Workshop?Challenges and theWay Forward，2009.

［6］Wei X D.Interpretational applications ofspectral decomposition in identifyingminor faults［R］.72nd EAGE Conference&Exhibition. Barcelona，2010.

［7］GaborD.Theory ofcommunication.Part 1:the analysis ofinforma?tion［J］.Journal of the Institution ofElectrical Engineers-PartⅢ：Radio and Communication Engineering，1946，93（26）：429-441.

［8］Potter R K，Kopp G A，Kopp H G.Visible speech［M］.New York: DoverPublications Inc.，1947.

［9］周家雄，張國棟，尚帥.基于重排Gabor變換的高分辨率譜分解［J］.世界地質，2013，32（1）：153-157.

Zhou Jiaxiong，Zhang Guodong，Shang Shuai.High?resolution spec?trum decomposition based on reassigned Gabor transform［J］.World Geology，2013，32（1）：153-157.

［10］Morlet J，Arens G，F(xiàn)ourgeau E，etal.Wave propagation and sam?pling theory—PartⅠ：complex signaland scattering inmultilayered media［J］.Geophysics，1982，47（2）：203-221.

［11］周懷來，李錄明.廣義S變換在地震信號特征信息提取中的應用［J］.新疆石油地質，2008，29（6）：758-760.

Zhou Huailai，Li Luming.Application ofgeneralized S?transform to extraction ofseismic signal characteristic information［J］.Xinjiang Petroleum Geology，2008，29（6）：758-760.

［12］Kodera K，De Villedary C，Gendrin R.A new method for the nu?merical analysis ofnon?stationary signals［J］.Physics of the Earth and Planetary Interiors，1976，12（2）：142-150.

［13］Auger F，F(xiàn)landrin P.Improving the readability of time?frequency and time?scale representationsby the reassignmentmethod［J］.Sig?nal Processing，1995，43（5）：1 068-1 089.

［14］Thomson D J.Spectrum estimation and harmonic analysis［J］.Pro?ceedings ofthe IEEE，1982，70（9）：1 055-1 096.

［15］Flandrin P，Auger F，Chassande?Mottin E.Applications in time?fre?quency signalprocessing［M］.Arizona：CRCPress，2003.

［16］Bayram M.Multiplewindow time?frequency analysis［D］.Houston：Rice University，1996.

［17］丁夏畦，丁毅.Hermite展開與廣義函數(shù)［M］.武漢：華中師范大學出版社，2005.

Ding Xiagui，Ding Yi.Hermite expansion and generalized functions［M］.Wuhan:Huazhong NormalUniversity Press，2005.

［18］張波，王真理，周水生.譜分解在含氣檢測中的應用［J］.地球物理學進展，2010，25（4）：1 360-1 364.

Zhang Bo，Wang Zhenli，Zhou Shuisheng.Application ofspectrum decomposition to gas detection［J］.Progress in Geophysics，2010，25（4）：1 360-1 364.

SpectralDecomposition Technology Based on M ultipleW indowsReassigned Spectrum

CUIZhen,CAOSiyuan,LIUWei,MA Yuanyuan

(China University ofPetroleum,a.State Key Laboratory ofPetroleum Resource and Prospecting; b.CNPCKey Lab ofGeophysics Exploration,Beijing 102249,China)

Spectral decomposition decomposes a time domain signal to the time?frequency domain,which can characterize the local time?frequency properties and instantaneous frequency.The time?frequency resolution of conventional spectral decomposition methods is not high due to Heisenberg uncertainty principle.The reassigned spectrum method avoids the Heisenberg uncertainty principle and the time?frequency resolution is increased.Besides,italso generates concentrated energy for random noise.Themultiple windows reassigned spec?trum is based on reassigned spectrum and gets reassigned spectrums using orthogonality and recurrence.At last,the reassigned spectrums are distributed by theweighted value based on least square estimation.Synthetic data and field data show that the proposed approach can gethigh time?frequency resolution and disperse the random noise;meanwhile the realposition oftime?frequency energy can be located and indicates the low?frequency shadow related with hydrocarbon effectively.

spectraldecomposition;Hermitewindow function;multiplewindows;reassigned spectrum;weighted value combination

P631.445

A

1001-3873（2015）05-0602-05

10.7657/XJPG20150520

2014-12-20

2015-05-22

國家科技重大專項（2011ZX05024-001-01）

崔震（1991-），男，河北滄州人，碩士研究生，地震資料處理，（Tel）18810906151（E-mail）761683429@qq.com.