一種基于匹配追蹤的Wigner高階譜地震信號譜分解方法
【技術領域】
[0001] 本發(fā)明屬于非平穩(wěn)信號時頻分析及地震信號處理領域，具體涉及一種基于匹配追 蹤的Wigner高階譜地震信號譜分解方法。
【背景技術】
[0002] 譜分解是一種地震信號解釋技術，它將地震數(shù)據(jù)分解到時頻域，由此揭示出時頻 域中包含有用的含油氣信息。多項研究都證明了該方法在儲層解釋和預測方面的實用性， 例如進行儲層厚度估計，地層解釋以及流體識別等。傳統(tǒng)的譜分解方法，如短時傅里葉變 換，它通過引入窗函數(shù)來計算信號的時頻分布，故而其時頻分辨率被相應的窗函數(shù)所限制， 不能滿足高精度地震勘探對精細儲層預測的要求。Ville (1948)將Wigner分布（1932)引 入了信號處理領域，提出著名的 Wigner-Ville 分布（Wigner-Ville Distribution，WVD)， 它在時頻平面上直接定義能量密度，不需要受分辨率的限制。業(yè)已普遍承認，沒有任何一種 時頻聯(lián)合分布的時頻分辨率能出其右，然而在實際應用中，該變換的二次性產生了干擾，即 WVD的交叉項問題。
[0003] 在眾多譜分解方法中，由Mallat和Zhang(1993)提出的匹配追蹤（Matching pursuit，MP)算法具有很高的時頻分辨率，得到了廣泛的應用。MP是一種靈活的自適應信 號分解方法，它結合信號的先驗信息構造合適的冗余字典，將信號在該字典上分解，自適應 地得到最能夠匹配信號實際結構的分解表達式。同時，由于字典中的原子具有很好的時頻 聚集性，將其與信號的WVD相結合，一方面可以利用WVD的最佳時頻分辨率提高MP時頻譜 的分辨率，另一方面可以解決WVD的交叉項問題。
[0004] Chakraboty和0kaya(1995)首先將MP算法引入地震信號譜分解中，取得了良好的 效果，自此，MP被廣泛地應用于地震信號低頻陰影檢測，時頻屬性提取等。但MP是一個密集 計算過程，計算效率低，并且其結果依賴于字典的構建。Liu等（2004, 2005)引入Morlet小 波和Ricker子波來構建字典，并將復數(shù)道分析加入地震信號的MP分解中，提出利用Morlet 小波的動態(tài)MP，提高了 MP的效率。Wang (2007)給出了完整的基于復數(shù)域動態(tài)最優(yōu)搜索的MP 分解方法，并導出一種快速計算表達式。此后，該方法被廣泛應用于地震信號時頻分析中， 并得到了一定程度的發(fā)展和改進。張繁昌等（2010,2013)提出了利用動態(tài)子波庫的雙參數(shù) 和單參數(shù)掃描算法，使算法效率得到大幅度提尚。黃桿東等（2012)在Morlet小波中引入 能量衰減因子，改進了算法靈活性和重構精度。Zhao等（2012)提出將字典擴展為Ricker 子波、Morlet小波和多相位地震子波的集合，表明多成分字典更能夠有效地反映地震信號 中包含的信息。
[0005] 隨著信號處理技術的發(fā)展，信號高階統(tǒng)計量具有一些重要的性質，使得人們對信 號高階統(tǒng)計量的研究逐漸重視起來，許多學者從高階統(tǒng)計量的角度來討論WVD，即高階統(tǒng)計 量與WVD相結合的一種產物--Wigner高階譜時頻分布。Gerr (1988)首先提出三階Wigner 分布。Swami (1992)對時變高階譜的定義、性質及應用作了進一步研究，取得了較大的進展。 Wigner高階譜是WVD在高階譜域的延伸，既繼承了時頻分布同時反映信號頻譜成分和時間 之間關系的優(yōu)點，還引入了高階譜對高斯噪聲的良好抑制能力，同時可以獲得時頻聚集性 更好的時頻譜。然而，Wigner高階譜與WVD-樣同樣存在交叉干擾項的問題。并且隨著階 數(shù)的增加，計算量大大提升，因此Wigner高階譜的計算主要是Wigner雙譜和Wigner三譜。

【發(fā)明內容】

[0006] 本發(fā)明提供了一種基于匹配追蹤的Wigner高階譜地震信號譜分解方法，旨在去 除Wigner高階矩譜的交叉項，并獲得高時頻聚集性的譜分解結果，提高儲層預測的精度。 具體的技術方案如下所述。
[0007] 本發(fā)明為了解決上述技術問題，采用以下技術方案：
[0008] 本發(fā)明中所述方法是為了克服上述現(xiàn)有技術的缺點，主要針對去除Wigner高階 矩譜的交叉項，并獲得高時頻聚集性的譜分解結果，提高儲層預測的精度的問題，提出了一 種基于匹配追蹤的Wigner高階譜地震信號譜分解方法，包括以下步驟：
[0009] 步驟1 :讀入地震剖面的q道數(shù)據(jù)，選定原子類型；
[0010] 步驟2 :設置初始取變量i = 1 ;
[0011] 步驟3 :讀取第i道地震數(shù)據(jù)Xl (t)，設置匹配追蹤分解次數(shù)N，將第η次分解的信 號記為Rn⑴，用Xi⑴對札⑴賦值，即札(t) = Xi (t)，n e [1，2,…，Ν];
[0012] 步驟4 :，設置初始取變量η = 1 ;
[0013] 步驟5 :對分解信號Rn(t)進行復數(shù)道分析，確定Rn(t)對應原子的初始時間延遲 u。、初始頻率ω。和初始相位1%三個參數(shù)；
[0014] 步驟6 :對尺度因子進行全局搜索，確定分解信號Rn(t)對應原子的初始尺度因子 σ。，由此得到初始參數(shù)集％: = 叫，％丨；
[0015] 步驟7 :對初始參數(shù)集Ad%，丨進行局部搜索，找到與信號Rn(t)最匹配 的原子f);
[0016] 步驟8 :計算信號Rn(t)的最匹配原子S';.,.⑴的k階Wigner高階譜的對角切片譜 SWkgJt,f)；
[0017] 步驟9 :計算信號比(〇在最匹配原子^⑴方向上投影后的殘差，并將殘差視為新 的信號Rn+1 (t)，令n = n+1，重復步驟5-8,直到達到最大分解次數(shù)N ;
[0018] 步驟10 :對η所有取值下的信號Rn⑴的最匹配原子:仏,⑴的k階Wigner高階譜 的對角切片譜求和，作為該道地震數(shù)據(jù)的Wigner高階譜時頻譜截取 單頻切片SFliV;
[0019] 步驟11 :令變量i = i+1，重復步驟3-10,直到取完整個地震剖面的數(shù)據(jù)，即直到i =q，由所有道地震數(shù)據(jù)的不同單頻切片可以構成整個二維剖面的一系列單頻屬性，即譜分 解的結果SFV= {SFliV，1彡i彡q}。
[0020] 上述技術方案中，所述步驟5包括以下幾個步驟：
[0021] 步驟 5. 1 :對 Rn(t)做 Hilbert 變換：
[0023] 其中，t是時間變量，P表示柯西主值，τ表示時間積分變量。
[0024] 步驟5. 2 :令X (t) = Rn⑴，確定Rn⑴的解析信號z⑴：
[0026] 其中，j是虛數(shù)單位，e是自然常數(shù)，a表示瞬時振幅；
[0027] 步驟5. 3 :確定瞬時振幅a (t):
[0029] 步驟5. 4 :選取瞬時振幅a (t)達到最大的時刻u。作為原子的時間延遲初始值：
[0030] iir, =argmaxi7(/)
[0031] 步驟5. 5 :確定原子的初始相位即該時刻的瞬時相位爐:
[0034] 步驟5. 6 :確定原子的初始頻率ω。，即該時刻的瞬時頻率ω (u。）：
[0037] 上述技術方案中，所述步驟6包括以下幾個步驟：
[0038] 步驟6. 1 :設置尺度因子〇的全局搜索范圍[σ σ 2]
[0039] 步驟6. 2 :根據(jù)尺度因子〇的數(shù)值計算原子：
L0041」 具中，gY⑴是通迓對基本汲形g⑴做Ζ 樹雙化而得到的原子波形表 達式，t是時間變量，j是虛數(shù)單位，e是自然常數(shù)；u為時間延遲，控制基本波形的時間中心； ω為頻率調制因子，控制基本波形的頻率中心；識為相位調制，控制波形的形狀；〇為尺度 因子，控制其在時間域的跨度；
[0042] 步驟6. 3 :計算信號Rn(t)與原子⑴的內積<Rn(t), (t)> :
[0044] 步驟6.4:找到使〈&(〇4"〇>達到最大值的尺度因子〇作為原子的初始尺度 因子〇。：
[0046] 步驟6. 5 :得到初始參數(shù)集γ。：
[0048] 其中，γ表示四個參數(shù)的集合用以描述一個原子。
[0049] 上述技術方案中，所述步驟7包括以下幾個步驟：
[0050] 步驟7. 1 :設置四參數(shù)集γ ξ的局部搜索范圍為[γ (J-Δ γ, yQ+A γ];
[0051] 步驟7. 2 :找到與信號Rn(t)最匹配的原子4,.(0:
[0053] 其中，Α ω是對應γ ξ參數(shù)集的原子，γ n表示第η次分解信號下的最佳參數(shù)集。
[0054] 上述技術方案中，所述步驟8包括以下步驟：
[0055] 步驟8. 1 :計算信號Rn(t)的最匹配原子及...⑴的k階Wigner高階譜
[0057] 其中，t是時間變量，f是頻率變量，j是虛數(shù)單位，e是自然常數(shù)，τ表示時間積 分變量，r、m和s表示整數(shù)變量。
[0058] 步驟8. 2 :計算信號Rn(t)的最匹配原子⑴的k階Wigner高階譜的對角切片譜
[0060] 上述技術方案中，所述步驟9中計算信號Rn(t)在最匹配原子方向上投影后的殘 差，并將殘差視為新的信號Rn+1 (t)，包括以下步驟：
[0061] 步驟9. 1 :計算信號Rn(t)在最匹配原子方向上投影后的殘差，并將殘差視為新的 十目號Rn+i⑴：
[0063] 上述技術方案中，所述步驟10包括以下步驟：
[0064] 步驟10. 1 :對第i道地震數(shù)據(jù)的所有原子的k階Wigner高階譜的對角切片譜求 和^ (〖,/):
[0066] 步驟10. 2 :將I (?,/)作為第i道地震數(shù)據(jù)的Wigner高階譜時頻譜；
[0067] 步驟10. 3 :截取單頻切片SFliV:
[0069] 其中，v表示單頻切片的頻率值。
[0070] 本發(fā)明是一種基于匹配追蹤的Wigner高階譜地震信號譜分解方法，并將該模型 應用于非平穩(wěn)信號時頻分析及地震信號處理領域。
[0071] 綜上所述，由于采用了上述技術方案，本發(fā)明的有益效果是：
[0072] 本發(fā)明將基于匹配追蹤的Wigner-Ville時頻分布擴展到基于匹配追蹤的Wigner 高階譜，利用匹配追蹤方法去除Wigner高階譜的交叉項，可以獲得時頻聚集性更高的地震 譜分解結果，為后續(xù)地震儲層預測和流體識別提供更精確的信息，具有較高的實用性。
【附圖說明】
[0073] 圖1為方法流程圖。
[0074] 圖2為一個原始地震剖面；
[0075] 圖3為一個Morlet原子；
[0076] 圖4為第15道地震數(shù)據(jù)；
[0077] 圖5為第一次分解的最佳匹配原子的波形；
[0078] 圖6為第一次分解的最佳匹配原子的時頻譜，6a為Wigner-Ville時頻譜，6b為 Wigner雙譜（k = 2)對角切片譜，6c為Wigner三譜（k = 3)對角切片譜；
[0079] 圖7為第15道地震數(shù)據(jù)的Wigner高階譜時頻譜，7a為Wigner-Ville時頻譜，7b 為Wigner雙譜（k = 2)對角切片譜，7c為Wigner三譜（k = 3)對角切片譜；
[0080] 圖8為地震數(shù)據(jù)的譜分解結果，8a為Wigner-Ville時頻譜得到的45Hz單頻剖面， 8b為Wigner雙譜（k = 2)對角切片譜得到的45Hz單頻剖面，8c為Wigner三譜（k = 3) 對角切片譜得到的45Hz單頻剖面。
【具體實施方式】
[0081] 本說明書中公開的所有特征，或公開的所有方法或過程中的步驟，除了互相排斥 的特征和/或步驟以外，均可以以任何方式組合。
[0082] 下