面結(jié)合圖1對本發(fā)明作詳細說明。
[0083] 為詳細說明本發(fā)明的技術(shù)內(nèi)容、構(gòu)造特征、所實現(xiàn)目的及效果，以下結(jié)合實施方式 并配合附圖詳予說明。
[0084] 本發(fā)明提出了一種基于匹配追蹤的Wigner高階譜地震信號譜分解方法，該模型 應(yīng)用于非平穩(wěn)信號時頻分析及地震信號處理取得良好的效果。整個算法的實現(xiàn)示意圖如圖 1，包括步驟：
[0085] 步驟1 :如圖2所示，讀入地震剖面的q道數(shù)據(jù)；選定原子類型，如圖3所示的 Morlet 原子；
[0086] 步驟2 :設(shè)置初始取變量i = 1 ;
[0087] 步驟3 :如圖4所示，讀取第i道地震數(shù)據(jù)Xl (t)，設(shè)置匹配追蹤分解次數(shù)N，將第η 次分解的信號記為Rn⑴，用Xl (t)對札(t)賦值，即札(t) = Xl (t)，n e [1，2,…，Ν];
[0088] 步驟4 :，設(shè)置初始取變量η = 1 ;
[0089] 步驟5 :對分解信號Rn(t)進行復(fù)數(shù)道分析，確定Rn(t)對應(yīng)原子的初始時間延遲 u。、初始頻率ω。和初始相位蚪三個參數(shù)；
[0090] 步驟 5. 1 :對 Rn(t)做 Hilbert 變換：
[0092] 其中，t是時間變量，P表示柯西主值，τ表示時間積分變量。
[0093] 步驟5. 2 :令X (t) = Rn⑴，確定Rn⑴的解析信號z⑴：
[0095] 其中，j是虛數(shù)單位，e是自然常數(shù)，a表示瞬時振幅；
[0096] 步驟5. 3 :確定瞬時振幅a (t):
[0098] 步驟5. 4 :選取瞬時振幅a (t)達到最大的時刻u。作為原子的時間延遲初始值：
[0100] 步驟5. 5 :確定原子的初始相位％，即該時刻的瞬時相位爐(〃。）：
[0103] 步驟5. 6 :確定原子的初始頻率ω。，即該時刻的瞬時頻率ω (u。）：
[0106] 步驟6 :對尺度因子進行全局搜索，確定分解信號Rn(t)對應(yīng)原子的初始尺度因子 〇 〇,由此得到初始參數(shù)集魏
[0107] 步驟6. 1 :設(shè)置尺度因子。的全局搜索范圍[0 u σ 2]
[0108] 步驟6. 2 :根據(jù)尺度因子〇的數(shù)值計算原子：
[0109]
[0110] 其中，gY(t)是通過對基本波形g(t)做r =hf7,嘆㈧變化而得到的原子波形表 達式，t是時間變量，j是虛數(shù)單位，e是自然常數(shù)；u為時間延遲，控制基本波形的時間中心； ω為頻率調(diào)制因子，控制基本波形的頻率中心；P為相位調(diào)制，控制波形的形狀；〇為尺度 因子，控制其在時間域的跨度；
[0111] 步驟6. 3 :計算信號Rn(t)與原子gY (t)的內(nèi)積<Rn(t), gY (t)> :
[0113] 步驟6.4:找到使〈比(〇^"〇>達到最大值的尺度因子〇作為原子的初始尺
[0114] 度因子
[0115] 步驟6. 5 :得到初始參數(shù)集γ。：
[0117] 其中，γ表示四個參數(shù)的集合用以描述一個原子。
[0118] 步驟7 :對初始參數(shù)集進行局部搜索，找到與信號Rn(t)最匹配 的原子心,夕)；
[0119] 步驟7. 1 :設(shè)置四參數(shù)集γ ξ的局部搜索范圍為[γ (J-Δ γ, yQ+A γ];
[0120] 步驟7. 2 :如圖5所示，找到與信號Rn(t)最匹配的原子：
[0122] 其中，心_/〇是對應(yīng)γ ξ參數(shù)集的原子，γ n表示第η次分解信號下的最佳參數(shù)集。
[0123] 步驟8 :計算信號Rn(t)的最匹配原子仏（/)的k階Wigner高階譜的對角切片譜 如圖6所示；
[0124] 步驟8. 1 :計算信號Rn(t)的最匹配原子SV, (0的k階Wigner高階譜 為：
[0126] 其中，t是時間變量，f是頻率變量，j是虛數(shù)單位，e是自然常數(shù)，τ表示時間積 分變量，r、m和s表示整數(shù)變量。
[0127] 步驟8. 2 :計算信號Rn(t)的最匹配原子gJO的k階Wigner高階譜的對角切片譜
[0129] Wigner高階譜的對角切片譜的計算主要為Wigner雙譜的對角切片譜（k = 2):
[0131] 和Wigner三譜的對角切片譜（k = 3):
[0133] 步驟9 :計算信號比⑴在最匹配原子^(0方向上投影后的殘差，并將殘差視為新 的信號Rn+1 (t)，令n = n+1，重復(fù)步驟5-8,直到達到最大分解次數(shù)N ;
[0134] 步驟9. 1 :計算信號Rn(t)在最匹配原子方向上投影后的殘差，并將殘差視為新的 十目號Rn+i⑴：
[0136] 步驟10 :對η所有取值下的信號Rn(t)的最匹配原子心,(〇的k階Wigner高階譜 的對角切片譜仏/)求和，作為該道地震數(shù)據(jù)的Wigner高階譜時頻譜G,/)，截取 單頻切片SFliV;
[0137] 步驟10. 1 :對第i道地震數(shù)據(jù)的所有原子的k階Wigner高階譜的對角切片譜求 和:?仏/):
[0139] 步驟10. 2 :將^；, (^/)作為第i道地震數(shù)據(jù)的Wigner高階譜時頻譜，如圖7所示；
[0140] 步驟10. 3 :截取單頻切片SFliV:
[0142] 其中，v表示單頻切片的頻率值。
[0143] 步驟11 :令變量i = i+Ι，重復(fù)步驟3-10,直到取完整個地震剖面的數(shù)據(jù)，即直到i =q，由所有道地震數(shù)據(jù)的不同單頻切片可以構(gòu)成整個二維剖面的一系列單頻屬性，即譜分 解的結(jié)果SFV= {SF liV，1彡i彡q}，如圖8所示。
【主權(quán)項】
1. 一種基于匹配追蹤的Wigner高階譜地震信號譜分解方法，其特征在于包括以下步 驟： 步驟1 :讀入地震剖面的q道數(shù)據(jù)，選定原子類型； 步驟2 :設(shè)置初始取變量i = 1 ; 步驟3 :讀取第i道地震數(shù)據(jù)X1 (t)，設(shè)置匹配追蹤分解次數(shù)N，將第η次分解的信號記 為 Rn ⑴，用 Xi ⑴對 R1 ⑴賦值，即 R1 (t) = Xi (t)，n e [1，2,…，Ν]; 步驟4 :，設(shè)置初始取變量η = 1 ; 步驟5 :對分解信號Rn(t)進行復(fù)數(shù)道分析，確定Rn(t)對應(yīng)原子的初始時間延遲u。、初 始頻率ω。和初始相位羚三個參數(shù)； 步驟6 :對尺度因子進行全局搜索，確定分解信號Rn (t)對應(yīng)原子的初始尺度因子σ。， 由此得到初始參數(shù)集：. 步驟7 :對初始參數(shù)集：進行局部搜索，找到與信號Rn (t)最匹配的原 子⑴； 步驟8 :計算信號Rn⑴的最匹配原子馬"_(:0:的k階Wigner高階譜的對角切片譜步驟9 :計算信號Rn(t)在最匹配原子心,.(0方向上投影后的殘差，并將殘差視為新的 信號Rn+1 (t)，令n = n+1，重復(fù)步驟5-8,直到達到最大分解次數(shù)N ; 步驟10 :對η所有取值下的信號Rn(t)的最匹配原子義,⑴的k階Wigner高階譜的對 角切片譜1求和，作為該道地震數(shù)據(jù)的Wigner高階譜時頻譜截取單 頻切片SFliv; 步驟11 :令變量i = i+Ι，重復(fù)步驟3-10,直到取完整個地震剖面的數(shù)據(jù)，即直到i = q，由所有道地震數(shù)據(jù)的不同單頻切片可以構(gòu)成整個二維剖面的一系列單頻屬性，即譜分解 的結(jié)果 SFv= {SFliV，1 彡 i 彡 q}。2. 根據(jù)權(quán)利要求1所述基于匹配追蹤的Wigner高階譜地震信號譜分解方法，其特征在 于，所述步驟5包括以下幾個步驟： 步驟5. 1 :對Rn⑴做Hilbert變換：其中，t是時間變量，P表示柯西主值，τ表示時間積分變量。 步驟5. 2 :令X (t) = Rn (t)，確定Rn⑴的解析信號ζ (t):其中，j是虛數(shù)單位，e是自然常數(shù)，a表示瞬時振幅； 步驟5. 3 :確定瞬時振幅a (t):步驟5. 4 :選取瞬時振幅a (t)達到最大的時刻u。作為原子的時間延遲初始值：步驟5. 5 :確定原子的初始相位熟即該時刻的瞬時相位0(?):步驟5. 6 :確定原子的初始頻率ω。，即該時刻的瞬時頻率ω (U。）：ω0= ω (u 〇)3. 根據(jù)權(quán)利要求1所述基于匹配追蹤的Wigner高階譜地震信號譜分解方法，其特征在 于，所述步驟6包括以下幾個步驟： 步驟6. 1 :設(shè)置尺度因子。的全局搜索范圍[〇 i，σ 2] 步驟6. 2 :根據(jù)尺度因子〇的數(shù)值計算原子：其中，gY (t)是通過對基本波形g(t)做變化而得到的原子波形表達 式，t是時間變量，j是虛數(shù)單位，e是自然常數(shù)；u為時間延遲，控制基本波形的時間中心； ω為頻率調(diào)制因子，控制基本波形的頻率中心；#為相位調(diào)制，控制波形的形狀；〇為尺度 因子，控制其在時間域的跨度； 步驟6. 3 :計算信號Rn(t)與原子gY (t)的內(nèi)積<Rn(t), gY (t)> :步驟6.4:找到使〈比(〇^"〇>達到最大值的尺度因子〇作為原子的初始尺 度因子σ。：步驟6. 5 :得到初始參數(shù)集γ。：其中，γ表示四個參數(shù)的集合用以描述一個原子。4. 根據(jù)權(quán)利要求1所述基于匹配追蹤的Wigner高階譜地震信號譜分解方 法，其特征在于，所述步驟7包括以下幾個步驟： 步驟7. 1 :設(shè)置四參數(shù)集γ ξ的局部搜索范圍為步驟7. 2 :找到與信號Rn (t)最匹配的原子其中，.I 是對應(yīng)γξ參數(shù)集的原子，γ n表示第η次分解信號下的最佳參數(shù)集。5. 根據(jù)權(quán)利要求1所述基于匹配追蹤的Wigner高階譜地震信號譜分解方法，其特征在 于，所述步驟8包括以下步驟： 步驟8. 1 :計算信號Rn(t)的最匹配原子(0的k階Wigner高階譜為：其中，t是時間變量，f是頻率變量，j是虛數(shù)單位，e是自然常數(shù)，τ表示時間積分變 量，r、m和s表示整數(shù)變量。 步驟8. 2 :計算信號Rn(t)的最匹配原子S；..⑴的k階Wigner高階譜的對角切片譜6. 根據(jù)權(quán)利要求1所述基于匹配追蹤的Wigner高階譜地震信號譜分解方法，其特征在 于，所述步驟9中計算信號Rn (t)在最匹配原子方向上投影后的殘差，并將殘差視為新的信 號Rn+1(t)，包括以下步驟： 步驟9. 1 :計算信號Rn(t)在最匹配原子方向上投影后的殘差，并將殘差視為新的信號 Rn+i(t)：7. 根據(jù)權(quán)利要求1所述基于匹配追蹤的Wigner高階譜地震信號譜分解方法，其特征在 于，所述步驟10包括以下步驟： 步驟10. 1 :對第i道地震數(shù)據(jù)的所有原子的k階Wigner高階譜的對角切片譜求和步驟10. 2 :將作為第i道地震數(shù)據(jù)的Wigner高階譜時頻譜； 步驟10. 3 :截取單頻切片SFliv:其中，V表示單頻切片的頻率值。
【專利摘要】一種基于匹配追蹤的Wigner高階譜地震信號譜分解方法，首先讀入地震剖面，選定原子類型；然后選取一道地震數(shù)據(jù)；接著對信號進行復(fù)地震道分析以及對尺度因子進行全局搜索，確定原子的初始參數(shù)集；對參數(shù)集進行局部搜索，找到與信號最匹配的原子；計算最匹配原子的Wigner高階譜的對角切片譜；計算該信號在最匹配原子方向上投影后的殘差，并將其視為新的分解信號；對所有分解得到的原子的Wigner高階譜的對角切片譜求和，作為該道地震數(shù)據(jù)的Wigner高階譜時頻譜，截取單頻切片；對所有地震數(shù)據(jù)采用相同方法以得到譜分解的結(jié)果。本發(fā)明利用匹配追蹤方法去除Wigner高階譜的交叉項，可以獲得時頻聚集性更高的地震譜分解結(jié)果，為后續(xù)地震儲層預(yù)測和流體識別提供更精確的信息。
【IPC分類】G01V1/30
【公開號】CN105353408
【申請?zhí)枴緾N201510815579
【發(fā)明人】彭真明, 王雨青, 李新彥, 王曉陽, 孔德輝, 何艷敏, 田琳
【申請人】電子科技大學(xué)
【公開日】2016年2月24日
【申請日】2015年11月20日