本發(fā)明涉及地球物理勘探領(lǐng)域,具體地說,涉及一種頻率域的子波分解方法。
背景技術(shù):
頻譜成像技術(shù)是近年來發(fā)展起來的一項基于頻率譜分解的儲層特色解釋技術(shù),是地震屬性分析中重要組成部分。該方法已逐漸成為研究復(fù)雜油氣區(qū)域的一種有價值的后期處理技術(shù),對分析薄儲層、描述沉積特性具有顯著的效果。頻譜成像技術(shù)具有在空間橫向上分辨率高的特點,是一種利用三維地震資料的多尺度信息對儲層進行高分辨率成像、檢測儲層時間厚度變化的工具。
頻譜成像技術(shù)的核心為信號的時頻分析方法,時頻分析是分析非平穩(wěn)信號的傳統(tǒng)方法。時頻分析方法一般分為線性時頻分析方法、雙線性時頻分析方法和參數(shù)化時頻分析方法等。
典型的線性時頻分析方法包含短時傅里葉變換、小波變換、S變換等。由于線性時頻分析方法受海森堡(Heisenberg)測不準(zhǔn)原理的約束,信號的時間分辨率和頻率分辨率為一對矛盾體,不能同時達到較高的時間分辨率和較高的頻率分辨率。典型的雙線性時頻分析方法為Wigner-Ville分布(即WVD)。由于WVD中不含任何窗函數(shù),避免了線性時頻分析中時間分辨率和頻率分辨率矛盾的問題,在分解單分量平穩(wěn)信號時具有很高的時頻分辨率。但是WVD為雙線性的,在分解多分量非平穩(wěn)信號時會出現(xiàn)嚴重的交叉項。
最典型的參數(shù)化時頻分析方法為匹配追蹤法(即MP算法),亦稱為子波分解法。匹配追蹤法為S.Mallat與Z.Zhang于1993年提出的,它是通過對時頻原子進行伸縮、時移、頻移和相移得到一個原子庫,同時根據(jù)最大匹配投影原理尋找最佳時頻原子的線性組合,以對原始信號分解為多個子波,且具有較高的時頻分辨率。但是MP算法是一種貪婪算法,它是通過不斷地迭代尋找局部最佳時頻原子的方法來實現(xiàn)對原信號的不斷逼近過程,因此,運算量較大且計算效率較低。
因此,亟需一種能夠快速且準(zhǔn)確地搜索時間域子波分解中的時頻原子,使得 根據(jù)時域原子確定的多個子波的疊加結(jié)果更接近于原始信號,且獲得的各時頻原子更可靠的方法。
技術(shù)實現(xiàn)要素:
本發(fā)明的目的之一在于解決現(xiàn)有技術(shù)在進行頻率域的子波分解的過程中,對時頻原子的搜索過程運算量較大,且搜索得到的時頻原子精度較低的技術(shù)缺陷。
本發(fā)明提供一種頻率域的子波分解方法,包括:初始參數(shù)確定步驟,由地震信號x(t)確定初始時移un、初始相位φn和初始頻率fn,并獲得振幅譜的實部XRe(ω)和虛部XIm(ω),搜索地震信號的初始尺度σn;
局部優(yōu)化步驟,根據(jù)初始時移un、初始相位φn、初始頻率fn和初始尺度σn設(shè)定局部搜索范圍,在局部搜索范圍之內(nèi)尋找局部最優(yōu)的控制參數(shù)r′n={u′n,σ′n,f′n,φ′n},其中,σ′n為局部最優(yōu)尺度,u′n為局部最優(yōu)時移,φ′n為局部最優(yōu)相位,f′n為局部最優(yōu)頻率;
子波形態(tài)確定步驟,由局部最優(yōu)的控制參數(shù)r′n={u′n,σ′n,f′n,φ′n}獲得子波振幅an,進而確定子波形態(tài)信號m(t);
殘差確定步驟,將地震信號x(t)和子波形態(tài)信號m(t)的差值做為地震信號殘差x’(t);
循環(huán)迭代步驟,以地震信號殘差x’(t)代替地震信號x(t),循環(huán)執(zhí)行初始參數(shù)確定步驟、局部優(yōu)化步驟、子波形態(tài)確定步驟和殘差確定步驟獲得其他子波。
在一個實施例中,初始時移un為地震信號x(t)的瞬時振幅最大值對應(yīng)的時間,初始相位φn為該時間對應(yīng)的瞬時相位,初始頻率fn為該時間對應(yīng)的瞬時頻率。
在一個實施例中,在搜索地震信號的初始尺度σn的步驟中,
將初始時移un、初始相位φn和初始頻率fn設(shè)定為固定數(shù)值;
選擇實部的初始尺度σn-Re,使得根據(jù)實部的初始尺度σn-Re得到的子波振幅譜的實部MRe(ω)與地震信號振幅譜的實部XRe(ω)的匹配程度最高;
選擇虛部的初始尺度σn-Im,使得根據(jù)虛部的初始尺度σn-Im得到的子波振幅譜的虛部MIm(ω)與地震信號振幅譜的虛部XIm(ω)的匹配程度最高;
為實部的初始尺度σn-Re和虛部的初始尺度σn-Im分別賦予權(quán)重并進行求和,確定初始尺度σn。
在一個實施例中,在局部優(yōu)化步驟中,
先設(shè)定時間偏移量△u、尺度偏移量△σ、頻率偏移量△f、相位偏移量△φ,確定局部搜索范圍為[rn-△r,rn+△r],其中rn={un,σn,fn,φn},△r=(△u,△σ,△f,△φ);
再在局部搜索范圍內(nèi),選擇實部的控制參數(shù)u′n-Re,σ′n-Re,f′n-Re,φ′n-Re,使得根據(jù)實部的控制參數(shù)得到的子波振幅譜的實部MRe(ω)與地震信號的實部XRe(ω)的匹配程度最高;以及
在局部搜索范圍內(nèi),選擇虛部的控制參數(shù)u′n-Im,σ′n-Im,f′n-Im,φ′n-Im,使得根據(jù)虛部的控制參數(shù)得到的子波振幅譜的虛部MIm(ω)與地震信號的虛部XIm(ω)的匹配程度最高;
為實部和虛部的控制參數(shù)分別賦予權(quán)重并進行求和,確定局部最優(yōu)的控制參數(shù)r′n={u′n,σ′n,f′n,φ′n}。
在一個實施例中,在所述局部最優(yōu)的控制參數(shù)中,
局部最優(yōu)時移u′n=αu′n-Re+(1-α)u′n-Im;
局部最優(yōu)尺度σ′n=βσ′n-Re+(1-β)σ′n-Im;
局部最優(yōu)頻率f′n=γf′n-Re+(1-γ)f′n-Im;
局部最優(yōu)相位φ′n=δφ′n-Re+(1-δ)φ′n-Im;
其中,權(quán)重參數(shù)α,β,γ,δ的取值范圍均為(0,1)。
在一個實施例中,在獲得子波振幅an的步驟中,
根據(jù)局部最優(yōu)的控制參數(shù)r′n={u′n,σ′n,f′n,φ′n}分別計算實部的子波振幅an-Re和虛部的子波振幅an-Im;
為實部的子波振幅an-Re和虛部的子波振幅an-Im分別賦予權(quán)重并進行求和,獲得子波振幅an。
在一個實施例中,所述子波振幅表示為:
an=εan-Re+(1-ε)an-Im
其中,權(quán)重參數(shù)ε的取值范圍為(0,1)。
在一個實施例中,所述子波形態(tài)信號表示為:
其中,σ′n為局部最優(yōu)尺度,u′n為局部最優(yōu)時移,φ′n為局部最優(yōu)相位,ω′n=2πf′n,f′n為局部最優(yōu)頻率,an為子波振幅。
在一個實施例中,在循環(huán)步驟中,
根據(jù)預(yù)設(shè)的迭代次數(shù)閾值或者殘差能量閾值判斷是否滿足迭代條件,在滿足迭代條件的情況下,循環(huán)執(zhí)行初始參數(shù)確定步驟、局部優(yōu)化步驟、子波形態(tài)確定步驟和殘差確定步驟。
在一個實施例中,所述殘差能量閾值設(shè)定為原始地震信號能量的5%。
本發(fā)明的實施例在頻率域?qū)Φ卣鹦盘栒穹V的實部和虛部分別進行子波分解??紤]到控制參數(shù)大小對實部和虛部形態(tài)影響不同,對子波形態(tài)控制參數(shù)采取實部和虛部分別匹配策略,分別設(shè)定一個系數(shù)決定實部和虛部的權(quán)重,從而能夠更精確的獲得局部子波形態(tài)。
本發(fā)明的其它特征和優(yōu)點將在隨后的說明書中闡述,并且,部分地從說明書中變得顯而易見,或者通過實施本發(fā)明而了解。本發(fā)明的目的和其他優(yōu)點可通過在說明書、權(quán)利要求書以及附圖中所特別指出的結(jié)構(gòu)來實現(xiàn)和獲得。
附圖說明
附圖用來提供對本發(fā)明的進一步理解,并且構(gòu)成說明書的一部分,與本發(fā)明的實施例共同用于解釋本發(fā)明,并不構(gòu)成對本發(fā)明的限制。在附圖中:
圖1是根據(jù)本發(fā)明實施例的頻率域子波分解方法的步驟流程圖;
圖2是一個實際地震信號測試結(jié)果的對比圖;
圖3是各局部子波形態(tài)展示圖;
圖4是各局部子波的振幅譜與原始地震信號振幅譜的對比圖;
圖5a是所有局部子波頻譜的實部相加與原始地震信號頻譜的實部對比圖;
圖5b是所有局部子波頻譜的虛部相加與原始地震信號頻譜的虛部對比圖;
圖6是所有子波的相加后的振幅譜與原始地震信號振幅譜的對比圖。
具體實施方式
為使本發(fā)明的目的、技術(shù)方案和優(yōu)點更加清楚,以下結(jié)合附圖對本發(fā)明作進一步地詳細說明。
以下結(jié)合說明書附圖對本發(fā)明的實施例進行說明,應(yīng)當(dāng)理解,此處所描述的優(yōu)選實施例僅用于說明和解釋本發(fā)明,并不用于限定本發(fā)明。并且在不相沖突的情況下,本發(fā)明的實施例中的特征可以相互結(jié)合。
本發(fā)明的實施例在迭代過程中采用全局粗粒度預(yù)測和局部最優(yōu)相結(jié)合的方法,考慮了尺度的重要性,提出了尺度濾波的思想,提高了算法的計算效率。在時間域搜索局部最佳時頻原子,將時頻原子作為一個整體去尋找最優(yōu)的一組參數(shù)。由于時頻原子的各參數(shù)對子波頻譜的實部和虛部影響程度是不同的,因此,本發(fā)明的實施例在匹配過程中為時頻原子的實部、虛部設(shè)定不同的權(quán)重,從而保證搜索到的局部時頻原子更加準(zhǔn)確。
實施例
圖1為本實施例的頻率域子波分解方法的步驟流程圖。
首先,執(zhí)行初始參數(shù)確定步驟。在本步驟中,先根據(jù)時域地震信號x(t)確定初始時移un、初始相位φn和初始頻率fn,并獲得振幅譜的實部XRe(ω)和虛部XIm(ω),再根據(jù)初始時移un、初始相位φn和初始頻率fn搜索地震信號的初始尺度σn。
如圖1所示,先對時域地震信號x(t)進行希爾伯特變換,計算瞬時振幅、瞬時相位和瞬時頻率(步驟S110),獲得初始時移un、初始相位φn和初始頻率fn(步驟S120)。其中,初始時移un為地震信號x(t)的瞬時振幅最大值對應(yīng)的時間,初始相位φn為該時間對應(yīng)的瞬時相位,初始頻率fn為該時間對應(yīng)的瞬時頻率。
然后,對時域地震信號x(t)進行傅里葉變換,獲得地震信號振幅譜的實部XRe(ω)和虛部XIm(ω)(步驟S130)。
隨后,將初始時移un、初始相位φn和初始頻率fn設(shè)定為固定數(shù)值,使得子波振幅譜的實部和虛部分別去匹配地震信號振幅譜的實部和虛部,搜索地震信號的初始尺度σn(步驟S140)。
具體來說,分別根據(jù)表達式(1)和(2)計算子波振幅譜的實部MRe(ω)和虛部MIm(ω):
其中,σn為初始尺度,un為初始時移,φn為初始相位,ωn=2πfn,fn為初 始頻率。
初始尺度σn是在一組固定數(shù)值的un、φn和fn情況下,通過計算最優(yōu)化公式(3)得到的。
其中,D={Mr(t)}為時頻原子的頻譜字典,是函數(shù)R(n)X和的內(nèi)積,且
需要說明的是,表達式(3)僅僅為優(yōu)化處理的通式,本實施例對實部的初始尺度σn-Re和虛部的初始尺度σn-Im分別進行計算。將初始時移un、初始相位φn和初始頻率fn設(shè)定為固定數(shù)值,根據(jù)最大投影匹配原理確定實部的初始尺度σn-Re和虛部的初始尺度σn-Im,具體說明如下。
實部的初始尺度表示為:
其中,R(n)表示初始控制參數(shù)rn={un,σn,fn,φn},XRe表示地震信號振幅譜的實部XRe(ω),表示由初始控制參數(shù)rn得到的子波振幅譜的實部MRe(ω),為R(n)XRe與的內(nèi)積,表示R(n)XRe在方向的投影,表示向量的投影長度,表示時頻原子的頻譜詞典中的實部成分。
則表達式(4)表明,σn-Re為在DRe范圍內(nèi)使取得最大值的初始尺度的實部數(shù)值。也就是說,根據(jù)實部的初始尺度σn-Re得到的子波振幅譜的實部MRe(ω)與地震信號振幅譜的實部XRe(ω)的匹配程度最高。
虛部的初始尺度表示為:
其中,XIm表示地震信號振幅譜的虛部XIm(ω),表示由初始控制參數(shù)rn得到的子波振幅譜的虛部MIm(ω),為R(n)XIm與的內(nèi)積,表 示R(n)XIm在防線的投影,表示向量的投影長度。表示時頻原子的頻譜詞典中的虛部成分。
則表達式(5)表明,σn-Im為在DIm范圍內(nèi)使取得最大值的初始尺度的虛部數(shù)值。也就是說,根據(jù)虛部的初始尺度σn-Im得到的子波振幅譜的虛部MIm(ω)與地震信號振幅譜的虛部XIm(ω)的匹配程度最高。
隨后,為實部的初始尺度σn-Re和虛部的初始尺度σn-Im分別賦予權(quán)重并進行求和,確定初始尺度σn,表示為:
σn=βσn-Re+(1-β)σn-Im,
其中,β和(1-β)分別為實部σn-Re和虛部σn-Im的權(quán)重,β∈(0,1)。
至此為止,在步驟S110至步驟S140中確定子波形態(tài)的初始控制參數(shù)rn={un,σn,fn,φn}的數(shù)值,這里的控制參數(shù)也就是時頻原子的控制參數(shù)。
隨后,執(zhí)行局部優(yōu)化步驟,根據(jù)初始時移un、初始相位φn、初始頻率fn和初始尺度σn設(shè)定局部搜索范圍,在局部搜索范圍之內(nèi)尋找局部最優(yōu)的控制參數(shù)r′n={u′n,σ′n,f′n,φ′n}(步驟S150)。
在局部優(yōu)化的過程中,先設(shè)定時間偏移量△u、尺度偏移量△σ、頻率偏移量△f、相位偏移量△φ,確定局部搜索范圍為[rn-△r,rn+△r],其中rn={un,σn,fn,φn},△r=(△u,△σ,△f,△φ)。
再在局部搜索范圍內(nèi),利用表達式(1)計算一組子波振幅譜的實部選擇實部的控制參數(shù)u′n-Re,σ′n-Re,f′n-Re,φ′n-Re,使得根據(jù)實部的控制參數(shù)得到的子波振幅譜的實部M′Re(ω)與地震信號的實部XRe(ω)的匹配程度最高。其中,匹配的過程可采用類似于表達式(4)的計算過程。
進而,在局部搜索范圍內(nèi),利用表達式(2)計算一組子波振幅譜的虛部選擇虛部的控制參數(shù)u′n-Im,σ′n-Im,f′n-Im,φ′n-Im,使得根據(jù)虛部的控制參數(shù)得到的子波振幅譜的虛部M′Im(ω)與地震信號的虛部XIm(ω)的匹配程度最高。其中,匹配的過程可采用類似于表達式(5)的計算過程。
為實部和虛部的控制參數(shù)分別賦予權(quán)重并進行求和,確定局部最優(yōu)的控制參數(shù)r′n={u′n,σ′n,f′n,φ′n}。
各個局部最優(yōu)控制參數(shù)分別為:
局部最優(yōu)時移u′n=αu′n-Re+(1-α)u′n-Im;
局部最優(yōu)尺度σ′n=βσ′n-Re+(1-β)σ′n-Im;
局部最優(yōu)頻率f′n=γf′n-Re+(1-γ)f′n-Im;
局部最優(yōu)相位φ′n=δφ′n-Re+(1-δ)φ′n-Im;
其中,權(quán)重參數(shù)α,β,γ,δ的取值范圍均為(0,1)。
隨后,執(zhí)行子波形態(tài)確定步驟,由局部最優(yōu)的控制參數(shù)r′n={u′n,σ′n,f′n,φ′n}獲得子波振幅an,進而確定子波形態(tài)信號m(t)(步驟S160)。
在步驟S160中,先根據(jù)局部最優(yōu)的控制參數(shù)r′n={u′n,σ′n,f′n,φ′n}分別計算實部的子波振幅an-Re和虛部的子波振幅an-Im。其中,子波振幅an的計算通式如下:
具體來說,實部的子波振幅an-Re的表達式為:
在表達式(7)中,R′(n)表示局部最優(yōu)的控制參數(shù)r′n={u′n,σ′n,f′n,φ′n},XRe表示地震信號振幅譜的實部XRe(ω),表示由局部最優(yōu)的控制參數(shù)r′n得到的子波振幅譜的實部M′Re(ω),為R′(n)XRe與的內(nèi)積,表示R′(n)XRe在方向的投影,表示向量的投影長度。
虛部的子波振幅an-Im的表達式為:
其中,XIm表示地震信號振幅譜的虛部XIm(ω),表示由局部最優(yōu)控制參數(shù)r′n得到的子波振幅譜的虛部M′Im(ω),為R′(n)XIm與的內(nèi)積,表示R′(n)XIm在方向的投影,表示向量的投影長度。
再為實部的子波振幅an-Re和虛部的子波振幅an-Im分別賦予權(quán)重并進行求和,獲得子波振幅an,表示為:
an=εan-Re+(1-ε)an-Im
其中,權(quán)重參數(shù)ε的取值范圍均為(0,1)。
需要強調(diào)的是,本實施里充分考慮了各控制參數(shù)對子波頻譜實部和虛部的影響的不同,引入了五個系數(shù)分別調(diào)節(jié)各控制參數(shù)實部和虛部的權(quán)重,能夠更準(zhǔn)確地定位時頻原子的時移、頻率、相位、尺度及振幅等信息,獲得局部最優(yōu)的時頻原子。
具體來說,時頻原子的實部相對于虛部而言,對振幅譜形態(tài)影響更大,所以各個系數(shù)的數(shù)值一般大于0.5;且每個參數(shù)對信號振幅譜影響程度又不同,如頻率的大小對時頻原子形態(tài)影響相對其他四個參數(shù)更大,所以頻率的系數(shù)γ較大,其他系數(shù)次之。
這樣,由局部最優(yōu)的控制參數(shù)r′n和子波振幅an根據(jù)表達式(9)得到時間域子波形態(tài)信號m(t):
其中,σ′n為局部最優(yōu)尺度,u′n為局部最優(yōu)時移,φ′n為局部最優(yōu)相位,ω′n=2πf′n,f′n為局部最優(yōu)頻率,an為子波振幅。
接下來,執(zhí)行殘差確定步驟,將地震信號x(t)和子波形態(tài)信號m(t)的差值做為地震信號殘差x’(t)(步驟S170),即地震信號殘差x′(t)=x(t)-m(t)。
隨后,執(zhí)行循環(huán)迭代步驟,以地震信號殘差x’(t)代替地震信號x(t),循環(huán)執(zhí)行初始參數(shù)確定步驟、局部優(yōu)化步驟、子波形態(tài)確定步驟和殘差確定步驟獲得其他子波。
其中,根據(jù)預(yù)設(shè)的迭代次數(shù)閾值或者殘差能量閾值判斷是否滿足迭代條件(步驟S180),在滿足迭代條件的情況下,回到步驟S110循環(huán)執(zhí)行初始參數(shù)確定步驟、局部優(yōu)化步驟、子波形態(tài)確定步驟和殘差確定步驟。
在這里,滿足迭代的條件可以是當(dāng)前的迭代次數(shù)小于預(yù)設(shè)的迭代次數(shù)閾值,或者當(dāng)前的殘差能量小于預(yù)設(shè)的殘差能量比重。
這樣,在各個迭代過程中將地震信號x(t)精確的分解為多個子波,且所確定子波的疊加能夠無限接近最初的地震信號x(t)。
應(yīng)用示例
以下通過一個具體的應(yīng)用示例說明本發(fā)明實施例的有益效果。在本示例中,設(shè)定迭代條件為迭代次數(shù)小于10次,各權(quán)重參數(shù)的數(shù)值分別設(shè)定為α=0.5、 β=0.6、γ=0.9、δ=0.7、ε=0.6。
圖2為本示例中用到的一道地震信號,實線為原始地震信號,圓圈表示利用本發(fā)明的實施例對地震信號進行子波分解得到的各子波相加的結(jié)果。從圖2中可以看出,兩條曲線基本重合,誤差較小。
圖3為采用本發(fā)明實施例獲得的各子波形態(tài)。本示例中的退出條件設(shè)置為迭代次數(shù)達到10次即退出,所以分解出10個局部最優(yōu)的子波形態(tài),并且子波的能量是在不斷降低的,第10個子波的能量已經(jīng)非常小了。
圖4展示了分解出的10個子波的振幅譜曲線與原始地震信號振幅譜對比圖。細實線為原始地震信號的振幅譜,方框、菱形、三角形等等分別為各子波的振幅譜,可以看出振幅譜能量在逐漸減小,與子波能量變化一致。圖4也反映出原始地震信號中能量較強的頻率并非最強信號的頻率,而是多個頻率信號疊加的結(jié)果。
圖5a展示了原始地震信號實部與各子波實部相加后的對比圖,圖5b展示了原始地震信號虛部與各子波虛部相加后的對比圖。圖5a和圖5b中,方框為信號的實部或虛部,星號為各子波頻譜實部或虛部的和。從中可以看出,實部和虛部匹配得都較好,所以圖6中的所有子波的相加后的振幅譜與原始地震信號振幅譜想比,也匹配得較好。
總之,本發(fā)明實施例的方法充分考慮了各參數(shù)對子波頻譜實部和虛部的影響程度不同,通過調(diào)節(jié)各參數(shù)實部和虛部的權(quán)重來獲得子波形態(tài)決定參數(shù),得到的局部子波形態(tài)更精確,分解結(jié)果更合理。
上述技術(shù)方案只是本發(fā)明的一種實施方式,對于本領(lǐng)域內(nèi)的技術(shù)人員而言,在本發(fā)明公開了應(yīng)用方法和原理的基礎(chǔ)上,很容易做出各種類型的改進或變形,而不僅限于本發(fā)明上述具體實施方式所描述的方法,因此前面描述的方式只是優(yōu)選的,而并不具有限制性的意義。
任何本發(fā)明所屬技術(shù)領(lǐng)域內(nèi)的技術(shù)人員,在不脫離本發(fā)明所公開的精神和范圍的前提下,可以在實施的形式上及細節(jié)上作任何的修改與變化,但本發(fā)明的專利保護范圍,仍須以所附的權(quán)利要求書所界定的范圍為準(zhǔn)。