崔震，曹思遠(yuǎn)，劉偉，馬媛媛（中國(guó)石油大學(xué)a.油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室；b.CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京102249）

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

崔震，曹思遠(yuǎn)，劉偉，馬媛媛
（中國(guó)石油大學(xué)a.油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室；b.CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京102249）

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

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

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

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

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

1　方法原理

1.1時(shí)頻重排原理

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

高斯窗h(t)=exp(-t2/2σ)/ 2πσ，固定σ后，根據(jù)Heisenberg測(cè)不準(zhǔn)準(zhǔn)則，窗口大小固定為，則由（1）式和（2）式可以看出，伽柏譜即是用的時(shí)頻窗，分別在時(shí)間和頻率方向?qū)φ鎸?shí)點(diǎn)譜的平滑。

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

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

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

在計(jì)算中，時(shí)頻分析實(shí)際上是對(duì)無限長(zhǎng)平穩(wěn)信號(hào)序列的截?cái)?，這種截?cái)嘈?yīng)往往使譜圖分辨率降低。針對(duì)這一現(xiàn)象，文獻(xiàn)［14］在發(fā)現(xiàn)多個(gè)正交窗口的組合可以削弱這種截?cái)嘈?yīng)的基礎(chǔ)上提出了多窗口方法，利用多個(gè)正交窗口{hk(t)，k∈N}獲得獨(dú)立的近似功率譜估計(jì)，綜合統(tǒng)計(jì)這些近似估計(jì)，最終得到平穩(wěn)信號(hào)功率譜估計(jì)，即

文獻(xiàn)［16］則將其推廣到時(shí)頻域：

1.3改進(jìn)的時(shí)頻重排

由于譜重排是譜圖平滑的逆過程，相比于普通的譜圖，譜重排對(duì)平穩(wěn)譜的截?cái)嘈Ч用黠@。從（6）式的分析可以看出，多窗口加權(quán)可以壓制畸變，重新得到穩(wěn)定的光滑譜，與原譜圖相比具有更高的聚焦性。

1.3.1Hermite窗函數(shù)

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

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

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

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

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

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

離散化得到矩陣形式：

對(duì)（12）式中dk求導(dǎo)，令

求解得

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

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

2　模型測(cè)試

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

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

圖1　合成地震記錄時(shí)頻分析結(jié)果

圖2　含噪合成地震記錄時(shí)頻分析結(jié)果

3　實(shí)例

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

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

圖3　單道實(shí)際記錄時(shí)頻分析結(jié)果

圖4　實(shí)際含氣地震剖面（紅色區(qū)域?yàn)楹瑲鈨?chǔ)集層）

圖5　頻率切片對(duì)比

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

4　結(jié)論

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

符號(hào)注釋

［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］周家雄，張國(guó)棟，尚帥.基于重排Gabor變換的高分辨率譜分解［J］.世界地質(zhì)，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變換在地震信號(hào)特征信息提取中的應(yīng)用［J］.新疆石油地質(zhì)，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］.武漢：華中師范大學(xué)出版社，2005.

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

［18］張波，王真理，周水生.譜分解在含氣檢測(cè)中的應(yīng)用［J］.地球物理學(xué)進(jìn)展，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

國(guó)家科技重大專項(xiàng)（2011ZX05024-001-01）

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