自適應核時頻分析方法的改進及應用

岳 龍， 劉懷山， 徐秀剛

(中國海洋大學 海底科學與探測技術教育部重點實驗室，青島 266100)

自適應核時頻分析方法的改進及應用

岳 龍， 劉懷山， 徐秀剛

(中國海洋大學 海底科學與探測技術教育部重點實驗室，青島 266100)

提高地震信號的時頻分辨率是時頻分析研究的重點，傳統(tǒng)的線性時頻分析方法(短時傅里葉變換、小波變換、S變換)時頻分辨率不高，而雙線性時頻分析方法(如魏格納分布)時頻分辨率高，但存在嚴重的交叉項干擾。自適應核時頻分析方法，屬于雙線性時頻分布，計算過程中采用隨時間變化的最優(yōu)核函數(shù)，在保證高分辨率的前提下，壓制信號模糊域的交叉項干擾，提高了時頻譜的可讀性。對比了四種時頻分析方法對于線性調頻信號的分析效果，并將聚焦性更好的自適應核方法應用于實際資料，并利用改進的自適應核時頻分析方法在瞬時屬性提取和譜分解識別斷層的應用中進行了初步試驗，計算效率明顯提高且取得良好的效果。

自適應核時頻分析； 瞬時屬性提??； 譜分解； 斷層識別

0 引言

時頻分析是從時間域和頻率域兩方面對信號進行描述，獲取信號在時-頻域的分布特征。傳統(tǒng)的時頻分析方法主要有：短時傅里葉變換、小波變換、S變換、魏格納分布和Cohen類時-頻分析[1]等。短時傅里葉變換、小波變換和S變換[2-4]屬于線性時頻分析方法，短時傅里葉變換受固定窗函數(shù)的影響，時間分辨率和頻率分辨率相互制約[5]。S變換結合了傅里葉變換和小波變換，使信號的時頻譜與傅里葉譜相關聯(lián)，但是S變換的基函數(shù)固定，在應用方面受到一定地限制[6]。魏格納分布[7]屬于雙線性時頻分布，具有很高的時頻分辨率，但是其存在嚴重的交叉項干擾，影響了信號有效時頻譜地識別。Cohen給出了雙線性時頻分布一般的表示形式，在變換公式中加入了核函數(shù)來壓制交叉項干擾的存在。但是不同的Cohen類時頻分布核函數(shù)單一，對某一類型的信號有效，不具有普適性。針對上述問題，Baraniuk等[8]提出了徑向高斯核時頻分布，其核函數(shù)能自適應不同類型信號，但是對于同一類型信號，其信號特征隨時間變化的情況不能夠很好地解決?；趶较蚋咚购说淖赃m應最優(yōu)核時頻分析方法(AOK)[9]解決了徑向高斯核時頻分布存在的問題，其核函數(shù)不僅可以適應不同類型的信號，而且可以隨信號的變化而變化，具有局部自適應的能力。

我們對幾種常用的時頻分析方法進行對比分析，并探討了自適應核時頻分析方法，在瞬時屬性提取和譜分解識別斷層的應用效果。

1 方法原理

1.1 雙線性時頻分布

目前已知的雙線性時頻分析方法，都可以歸結為Cohen類時頻分布[1]。首先給出信號模糊函數(shù)的定義:

(1)

其中：s*(t)是s(t)的共軛函數(shù);θ、τ是信號模糊域變量。

利用信號的模糊函數(shù)，得到信號s(t)的Cohen類雙線性時頻分布如式(2)。

e-jθt-jτωdθdτ

(2)

其中：P(t,ω) 是Cohen類時頻譜；Φ(θ,τ)是信號s(t)的核函數(shù)。

從式(2)可以看出，信號的雙線性時頻分布，就是對信號的模糊函數(shù)加權再進行二維傅里葉變換得到的。通過設定不同形態(tài)的核函數(shù)Φ(θ,τ)，使信號模糊函數(shù)的自項能量也就是反映信號特征的有效能量通過，壓制交叉項的干擾。

當設定核函數(shù)Φ(θ,τ)=1時，Cohen類時頻分布就變成了魏格納時頻分布，由于沒有對信號的交叉項干擾進行壓制，因此，魏格納時頻分布存在嚴重的交叉項干擾。

1.2 自適應核時頻分布

自適應核時頻分布采用短時模糊函數(shù)和隨時間變化的核函數(shù)[10]，即給信號加一個短時窗函數(shù)，在窗函數(shù)內部計算信號的模糊函數(shù)，同時在該時窗內計算最優(yōu)的核函數(shù)，以達到核函數(shù)隨信號的變化而變化的目的。

首先給出信號短時模糊函數(shù)的定義:

(3)

自適應核時頻分布的定義，如式(4)。

Φopt(t;θ,τ)e-jθt-jτωdθdτ

(4)

其中：Φopt(t;θ,τ)是給定窗函數(shù)內的最優(yōu)核函數(shù)；ω是頻率變量；Φopt(t;θ,τ)是通過求解式(5)最優(yōu)化問題得到的。

(5)

最優(yōu)化問題的約束條件為

(6)

(7)

1.3 自適應核實現(xiàn)算法

最優(yōu)核函數(shù)的快速計算需要在極坐標下完成[10]，計算過程需要用到極坐標下的短時模糊函數(shù)，其離散表達式為

bΔrsinqΔΨ

(8)

其中：Δr和ΔΨ是極坐標半徑和角度變量的步長；b=0、1、…、B-1;q=0、1、…、Q-1。

核函數(shù)極坐標下離散表達式為式(9)。

(9)

其中：σq是控制核函數(shù)的形狀的擴展因子。

σq=σ(qΔΨ)q=0,1,…，Q-1

(10)

計算最優(yōu)核函數(shù)的過程，即計算最優(yōu)擴展因子的過程，利用最速上升法來實現(xiàn)。計算步驟如式(11)和式(12)所示。

(11)

(12)

(13)

(14)

將式(13)和式(14)計算的結果代入式(11)，計算出最優(yōu)的核函數(shù)和短時模糊函數(shù)之后，利用式(4)的離散格式就可以計算自適應核時頻譜了。具體步驟如下：

1)在直角坐標系中計算信號的短時模糊函數(shù)A(n;θ,τ)。

2)將A(n;θ,τ)轉換到極坐標下表示得到Ap(n;b,q)。

3)利用式(11)～式(12)的迭代算法計算最優(yōu)的短時核函數(shù)的擴展因子σ。

4)計算最優(yōu)的核函數(shù)Φp(n;b,q)。

5)將Φp(n;b,q)轉換到直角坐標，并與短時模糊函數(shù)A(n;θ,τ)相乘，然后做二維傅里葉變換得到自適應核時頻譜。

1.4 自適應核參數(shù)

在實際應用過程中，自適應核時頻分析有兩個關鍵參數(shù)：①時窗寬度；②核函數(shù)體積。利用自適應核時頻分析完成不同的分析任務時，時窗寬度的設定有一個大致的準則，進行瞬時屬性提取、譜分解時，選擇比較小的時窗寬度，增大信號時頻譜的時間分辨率。經(jīng)過前期大量的實際地震資料的試驗和應用，得出時窗寬度在二分之一主周期和一個主周期(主頻對應的周期)之間效果較好。

在自適應核時頻分析算法中，核函數(shù)體積的大小沒有明確規(guī)定，只是給出了核函數(shù)體積范圍一般在1～5之間。根據(jù)大量地震資料的試驗和應用，得到在進行瞬時屬性提取、譜分解時，核函數(shù)體積對結果影響不是很大。考慮到計算過程盡量保持各頻率分量的衰減比不變[12]，選擇比較小的核函數(shù)體積。

1.5 改進的自適應核時頻分析

自適應核時頻分析方法在計算時頻譜的過程中，在每個樣點所在的數(shù)據(jù)段上，都需要計算一個最優(yōu)的核函數(shù)，因此計算量比較大。特別是在譜分解及其他屬性提取的計算過程中，需要對每個地震道都需要進行最優(yōu)核函數(shù)的計算，降低了運算效率。由于地下界面的連續(xù)性，地震剖面相鄰地震道之間具有很強的相關性。因此，在利用自適應核來進行譜分解和屬性提取的過程中，可以考慮利用相鄰地震道之間的相關性，來提高自適應核時頻方法的計算效率。

基本思想是，以奇數(shù)道為一個組合(例如3道、5道)，然后計算中間地震道的時頻譜，得到該道每個樣點所在數(shù)據(jù)段最優(yōu)的核函數(shù)。在計算相鄰道時頻譜時，直接使用中間道計算得到的最優(yōu)的核函數(shù)。這樣可以節(jié)省最優(yōu)核函數(shù)的一些計算過程，提高了計算效率。特別是數(shù)據(jù)量比較大時，計算效果比較明顯。

選取某二維地震剖面，并分別以3道和5道作為一個組合，采用改進前、后的自適應核時頻分析方法計算單個組合的時頻譜，并比較改進前、后的計算時間。計算參數(shù)分別為：時窗寬度為64個樣點，核函數(shù)體積為“1”，最優(yōu)核函數(shù)迭代次數(shù)分別為2和3，程序語言matlab，計算時間統(tǒng)計如表1所示。

表1 改進前后自適應核時頻分析計算時間對比Tab.1 Calculation time comparison of adaptivekernel time frequency analysis

從表1可以看出，隨著組合內道數(shù)的增加，計算效率提高的越多。因此當?shù)貙悠鸱兓容^小時，可以適當增大同一個組合內的道數(shù)，當?shù)貙悠鸱容^大時，適當減小組合內道數(shù)。如果增加最優(yōu)核迭代次數(shù)，改進的自適應核時頻分析加速效果更加明顯，一般情況下由于相鄰樣點最優(yōu)核函數(shù)變化不大，因此迭代次數(shù)一般為2次～3次即可。

2 自適應核時頻分析的應用

2.1 幾種時頻分析方法對比

時頻分析方法的時頻譜聚焦性能可以通過線性調頻信號進行檢驗，為此設計一個線性調頻信號如式(15)所示。該信號是兩個線性調頻信號的和，兩個信號的頻率變化范圍分別是0.5Hz～40Hz、0.5Hz～80Hz，隨著時間的增大信號的瞬時頻率逐漸增大，每個時間點的瞬時振幅是相同的。該信號長度為2s，采樣間隔為1ms，圖1為該線性調頻信號。

y= sin(2π(0.5+10*t)*t)+

sin(2π(0.5+20*t)*t)

(15)

其中y是兩個線性調頻信號的和。

圖2為采用不同的時頻分析方法對線性調頻信號進行時頻分析的結果。由圖2可以看出：魏格納分布(圖2(a))的兩條時頻譜能量線清晰，時頻譜聚焦性好，但是存在嚴重的交叉項干擾，影響了有效時頻譜的識別；高斯小波變換(圖2(b))的時頻譜不存在交叉項干擾，隨著信號頻率的增大，信號的時頻譜能量減小了，沒有能夠真實地反映信號的瞬時振幅特征；S變換(圖2(c))能較好地反映信號的時頻譜的變化規(guī)律，隨著頻率的變化信號的時頻譜能量保持一致，但在高頻端，S變換的時頻譜的寬度明顯增大，時頻譜的聚焦度明顯降低；自適應核時頻分布(圖2(d))采用最優(yōu)的核函數(shù)壓制交叉項干擾，同時保持高時頻譜聚焦度，兩條時頻譜能量線清晰，能量變化反映了線性調頻信號的特征。無論從聚焦度還是交叉項干擾的控制，自適應核時頻分布都是理想的對非平穩(wěn)信號進行分析的工具。

圖1 線性調頻信號Fig.1 The linear FM signal

圖2 不同時頻分析方法的效果對比Fig.2 Result of different time-frequency analysis methods(a)魏格納分布；(b)高斯小波變換；(c)S變換；(d)自適應核

2.2 瞬時屬性參數(shù)提取

地震波在傳播過程中，地震子波受到地層吸收、波前擴散和地層反射等影響，其頻率、振幅和相位不斷發(fā)生變化，并攜帶了大量的地下信息[10]，因此對地震記錄瞬時屬性的分析可以間接地反映地下地層的某些特性。利用地震瞬時屬性的多屬性聯(lián)合解釋技術[13]在實際應用(如儲層預測[14]、地震相分析[15]、薄層厚度反演[16]以及頻譜成像[17])中發(fā)揮了重要作用。利用自適應核時頻分析方法可以提取地震記錄的瞬時屬性參數(shù)，并且相對于常用的小波變換有更高的時頻譜分辨率。為了檢驗自適應核時頻分析方法的實用性，選擇某油田二維剖面進行分析，得到該剖面的瞬時屬性剖面。該數(shù)據(jù)采樣間隔為2ms，采樣點為2 001，共有221道。選取時間段為1 000ms～2 500ms、145道～185道的部分地震數(shù)據(jù)進行瞬時屬性提取，圖3即為選取的二維剖面的局部。

從圖4和圖5可以看出：自適應核得到的瞬時屬性和小波變換得到的結果形態(tài)相近，但分辨率要高于小波變換；圖4(b)瞬時振幅剖面能量聚焦，相鄰層的振幅能量界限更加明顯；圖5(b)瞬時頻率的變化更清晰，層間關系明確。因此，自適應核時頻分析可以用來計算地震記錄的瞬時屬性，而且提取的屬性分辨率相對于常規(guī)的時頻分析方法要高。

圖3 地震記錄Fig.3 Seismic record

圖4 瞬時振幅剖面Fig.4 The instantaneous amplitude profile(a)小波變換瞬時振幅；(b)自適應核瞬時振幅

圖5 瞬時頻率剖面Fig.5 The instantaneous frequency profile(a)小波變換瞬時頻率；(b)自適應核瞬時頻率

2.3 譜分解識別斷層

譜分解技術是利用時頻分析將地震信號變換到時頻域，提取某一頻率的振幅信息作為分頻剖面，然后利用不同的分頻剖面，解釋目的層的地質特征[18-19]。在實際應用中，可以利用譜分解分析確定斷層的空間形態(tài)[20]。為了檢驗自適應核方法在斷層識別方面的應用，利用自適應核時頻分析對某二維剖面，提取其多個單一頻率的振幅剖面，分析斷層的形態(tài)。該數(shù)據(jù)采樣間隔為2ms，采樣點為3 000，共有593道。選取1道～151道的部分地震數(shù)據(jù)進行譜分解，圖6(a)即為選取的二維剖面的局部。

從圖6中的標注位置可以看出，340ms～390ms的斷層，在原始地震剖面上很模糊，但是在不同頻率的譜分解剖面上，斷層的斷點變清晰，尤其是50Hz剖面，斷層的空間范圍和斷面位置很明確。

因此，對于斷距比較小的斷層，可以通過提取高于地震資料主頻的頻率切片來分析斷層的形態(tài)特征。由于地震波存在吸收衰減，深部地層高頻能量很弱，單純依靠高頻切片識別斷層有一定的限制。因此可以利用自適應核譜分解技術得到不同頻率的譜分解剖面，利用多個譜分解剖面綜合識別地震剖面上不同深度的斷層的空間形態(tài)。

圖6 原始地震剖面及不同頻率的譜分解剖面Fig.6 Seismic record and spectral decomposition profiles of different frequency(a) 地震剖面；(b)10Hz；(c)30Hz；(d)50Hz

3 結束語

這里論述了自適應核時頻分析方法的基本原理和計算的實現(xiàn)，并簡單探討了自適應核時頻分析方法在瞬時屬性提取和譜分解識別斷層等方面的應用，得到以下結論。

1) 自適應核時頻分析方法屬于雙性時頻分布，相對于小波變換和S變換等線性時頻分布，具有更高的時頻分辨率；相對于魏格納分布等Cohen類雙線性時頻分布，能夠很好地壓制交叉項干擾，精確地刻畫非平穩(wěn)信號的時頻譜。

2)自適應核時頻分析方法在屬性提取方面，提取的瞬時屬性具有比較高的聚焦度，層間界限明顯；譜分解剖面識別斷層，可以識別常規(guī)地震剖面不易識別的斷層，高頻切片在小斷距斷層的識別方面效果明顯。

3)在地震剖面相鄰道相似性的基礎上，提出了改進的自適應核時頻分析方法，提高了計算效率。

[1]COHENL.TimeFrequencyAnalysis:TheoryandApplication[M].EnglewoodCliffs,NJ:Prentice-Hall,1995.

[2] 馬見青,李慶春.S變換和TT變換聯(lián)合壓制地震面波[J].物探化探計算技術,2011,33(3):252-25.MAJQ,LIQC.JointStransformationandTT-transformationmethodofsurfacewavesuppression[J].ComputingTechniquesforGeophysicalandGeochemicalExploration, 2011,33(3):252-257.(InChinese)

[3] 段俊,張白林,潘樹林. 廣義S變換去除面波技術及應用[J].物探化探計算技術,2011, 33(3):280-285.DUANJ,ZHANGBL,PANSL.GeneralizedS-transformanditsapplicationtoremovesurfacewave[J].ComputingTechniquesforGeophysicalandGeochemicalExploration, 2011, 33(3):280-285.(InChinese)

[4]HLAWATSCHF,BOUDREAUX-BARTELSGF.Linearandquadratictime-frequencysignalrepresentation[J].IEEESignalProcessingMag, 1992, 9(2): 21-67.

[5] 劉小龍,王華,趙淑娥,等. 自適應最優(yōu)核時頻分布在地震儲層預測中的應用[J]. 中南大學學報(自然科學版), 2012,48(8): 3114-3120.LIUXL,WANGH,ZHAOSE,etal.Applicationofadaptiveoptimalkerneltime-frequencyrepresentationinreservoirprediction[J].JournalofCentralSouthUniversity(ScienceandTechnology), 2012, 48(8): 3114-3120.(InChinese)

[6] 劉喜武,年靜波,黃文松. 利用廣義S變換提取地震旋回的方法[J]. 石油物探, 2006,45(2): 129-133.LIUXW,NIANJB,HUANGWS.SeismiccyclesextractionusinggeneralizedS-transform[J].GeophysicalProspectingforPetroleum, 2006, 45(2):129-133.(InChinese)

[7] 劉喜武,劉婉瑩,劉洪,等.地震信號廣義時頻分析及其數(shù)值實現(xiàn)[J].物探化探計算技術,2007, 29(5):386-390.LIUXW,LIUWY,LIUH,etal.Generalizedseismicsignaltime-frequencyanalysisandnumericalalgorithms[J].ComputingTechniquesforGeophysicalandGeochemicalExploration, 2007, 29(5):386-390.(InChinese)

[8]BARANIUKRG,JONESDL.Signal-dependenttime-frequencyanalysisusingaradiallyGaussiankernel[J].SignalProcessing,1993,32(3):263-284.

[9]JONESDL,BARANIUKRG.Anadaptiveoptimal-kerneltime-frequencyrepresentation[J].IEEETrans.onSignalProcessing,1995,43 (10):2361-2371.

[10]馬成英. 時頻屬性提取方法及應用[D].青島：中國海洋大學,2012.MACY.Time-frequencyattributesextractionmethodsandapplication[D].Qingdao:OceanUniversityofChina, 2012.(InChinese)

[11]BARANIUKRG,JONESDL.Asignal-dependenttime-frequencyrepresentation:optimalkerneldesign[J].IEEETrans.onSignalProcessing, 1993, 41 (4):1589-1602.

[12]張宇輝,金國彬,李天云. 基于自適應最優(yōu)核時頻分布理論的間諧波分析新方法[J]. 中國電機工程學報, 2006, 26(18): 84-89.ZHANGYH,JINGB,LITY.Anovelapproachtointerharmonicsanalysisbasedonadaptiveoptimalkerneltime-frequencydistribution[J].ProceedingsoftheCSEE, 2006, 26(18):84-89. (InChinese)

[13]周懷來,李錄明. 廣義S變換在地震信號特征信息提取中的應用[J]. 新疆石油地質,2008,29(6):758-760.ZHOUHL,LILM.ApplicationofgeneralizedS-transformtoextractionofseismicsignalcharacteristicinformation[J].XinjiangPetroleumGeology, 2008,29(6):758-760.(InChinese)

[14]劉其,趙永勤,竇慧媛,等. 利用地震波主參數(shù)預測巖性和油氣[J]. 石油物探, 2003, 42(1): 77-81.LIUQ,ZHAOYQ,DOUHY,etal.Predictionoflithologyandhydrocarbonaccumulationwithseismicprincipalparameters[J].GeophysicalProspectingforPetroleum, 2003, 42(1):77-81.(InChinese)

[15]鄒文,陳愛萍,賀振華,等. 基于S變換的地震相分析技術[J]. 石油物探,2006,45(1):48-51.ZOUW,CHENAP,HEZH,etal.SeismicfaciesanalysisbasedonS-transform[J].GeophysicalProspectingforPetroleum, 2006, 45(1):48-51.(InChinese)

[16]張顯文,韓立國,王宇,等. 地震信號譜分解匹配追蹤快速算法及其應用[J]. 石油物探,2010,49(1):1-6.ZHANGXW,HANLG,WANGY,etal.Seismicspectraldecompositionfastmatchingpursuitalgorithmanditsapplication[J].GeophysicalProspectingforPetroleum, 2010, 49(1):1-6.(InChinese)

[17]劉蘭鋒,曹思遠,夏竹,等. 廣義S變換及其在儲集層頻譜成像中的應用[J]. 新疆石油地質,2009,30(5):629-632.LIULF,CAOSY,XIAZ,etal.GeneralizedS-Transformanditsapplicationtoreservoirimagingbyenergyspectrum[J].XinjiangPetroleumGeology, 2009,30(5):629-632.(InChinese)

[18]尚帥,韓立國,胡瑋,等. 壓縮小波變換地震譜分解方法應用研究[J]. 石油物探,2015,54(1):51-55.SHANGS,HANLG,HUW,etal.Appliedresearchofsynchrosqueezingwavelettransforminseismicspectraldecompositionmethod[J].GeophysicalProspectingforPetroleum, 2015, 54(1):51-55.(InChinese)

[19]胡水清,韓大匡,劉文嶺,等. 應用譜分解技術預測河道砂體[J]. 新疆石油地質,2008,29(3):370-372.HUSQ,HANDK,LIUWL,etal.Applicationofspectrumdecompositiontechniquetopredictionofriverchannelsands[J].XinjiangPetroleumGeology, 2008,29(3):370-372.(InChinese)

[20]張延慶,魏小東,王亞楠,等. 譜分解技術在QL油田小斷層識別與解釋中的應用[J]. 石油地球物理勘探, 2006, 41(5):584-586,591.ZHANGYQ,WEIXD,WANGYN,etal.ApplicationofspectrumfactorizationinrecognitionandinterpretationofminorfaultsinQLoilfield[J].OilGeophysicalProspecting, 2006, 41(5): 584-586,591.(InChinese)

Improvement of adaptive kernel time-frequency analysis and its application

YUE Long, LIU Huaishan, XU Xiugang

(Key Lab of Submarine Geosciences and Prospecting Techniques, Ministry of Education,Ocean University of China, Qingdao 266100,China)

Improving the seismic signal time-frequency resolution is the focus of time-frequency analysis research. The time-frequency resolution of the traditional linear time-frequency analysis methods, including short-time Fourier transform, wavelet transform and S transform is poor. Bilinear time-frequency analysis method such as Wigner distribution has a high time-frequency resolution but it is contaminated by serious cross-terms. Adaptive kernel time-frequency analysis belongs to bilinear time-frequency analysis. Under the premise of ensuring high resolution, it adopts the optimal kernel function changing with time to suppress cross-terms interference and improves the readability of time-frequency spectrum. This article compares the time-frequency spectrum of linear frequency modulated signal using for different methods. The adaptive kernel method has a better time-frequency concentration and been applied to real data and several tests have been carried out in instantaneous attributes extraction, fault identification using spectral factorization and sedimentary cycle analysis. The calculation effectiveness is improved and the results are satisfactory.

adaptive kernel time-frequency analysis; instantaneous attributes extraction; spectral factorization; fault identification

2015-12-30 改回日期：2016-03-01

國家自然科學基金(41176077，41230318，41204088);國家863 項目(2013AA092501);國土資源部海洋油氣資源與環(huán)境地質重點實驗室開放基金(MRE201303)

岳龍(1988-)，男，博士，主要從事地震勘探資料處理方法研究工作，E-mail: 957115928@qq.com。

1001-1749(2017)01-0044-08

P 631.4

A

10.3969/j.issn.1001-1749.2017.01.07