本發(fā)明屬于地震勘探識(shí)別技術(shù)領(lǐng)域，尤其與一種基于ceemd的高精度頻散avo屬性計(jì)算方法有關(guān)。
背景技術(shù)：
目前，大多數(shù)研究者都假設(shè)地震波速度會(huì)發(fā)生頻散，并利用頻譜分解技術(shù)把地震數(shù)據(jù)分解為窄帶地震信號(hào)，再利用地震波頻散理論及優(yōu)化算法提取地震波的頻散屬性。為了更有效地識(shí)別儲(chǔ)層流體，需要研究提高地震波頻散屬性計(jì)算精度與可靠性的方法，而最有效的提升計(jì)算精度方法，就是提高頻率域數(shù)據(jù)獲得準(zhǔn)確性。
頻譜分解技術(shù)基于各類時(shí)頻分析方法，時(shí)頻分析是針對(duì)非平穩(wěn)信號(hào)而逐漸發(fā)展起來(lái)的時(shí)間頻率域聯(lián)合分析方法。傅立葉變換只能給出全局信號(hào)各頻率的振幅譜，時(shí)間分辨率為0。時(shí)頻分析技術(shù)能夠直觀給出局部信號(hào)在各個(gè)頻率的能量分布，更適合分析地震信號(hào)等非平穩(wěn)信號(hào)。時(shí)頻分析方法主要分為線性和非線性方法兩類；線性時(shí)頻分析方法主要基于傅里葉變換的思想。法國(guó)地球物理學(xué)家于世紀(jì)年代提出了著名的小波變換，釆取低頻大尺度、高頻小尺度的分析方法，被譽(yù)為數(shù)學(xué)的“顯微鏡”。解決了時(shí)間頻率分辨率折中的難題，對(duì)信號(hào)具有更好的自適應(yīng)性和局部化特征。但小波變換存在幾個(gè)不足：由于小波變換時(shí)簡(jiǎn)諧波與高斯函數(shù)進(jìn)行同樣的伸縮和平移，小波變換是沿著時(shí)間方向絕對(duì)平移，小波變換相位局部化，各頻率相位基準(zhǔn)不一致；計(jì)算量較大，不便進(jìn)行數(shù)值計(jì)算。
經(jīng)驗(yàn)?zāi)B(tài)分解法(emd)是由huang等提出的一種信號(hào)分析方法，該方法通過(guò)提取復(fù)雜信號(hào)在每一個(gè)時(shí)刻局部的振蕩模式,按由高到低的自適應(yīng)頻率分解模式尋找信號(hào)內(nèi)蘊(yùn)的高頻信息，進(jìn)而分解得到若干個(gè)平穩(wěn)信號(hào)分量，即模態(tài)函數(shù)分量(imf)。emd分解法其應(yīng)用領(lǐng)域已經(jīng)遍及地震、雷達(dá)和語(yǔ)音信號(hào)處理及圖像分析等各個(gè)面。正是因?yàn)閑md的自適應(yīng)性，缺乏約束條件，使其不可避免的存在缺陷,這種缺陷稱為"模態(tài)混疊"。任何信號(hào)我們都可以看成其由若干個(gè)固有模態(tài)函數(shù)(imf)組成的,一個(gè)模態(tài)描述一個(gè)單一的震動(dòng)狀態(tài),而如果imf之間相互重疊，則形成復(fù)合信號(hào)。在經(jīng)驗(yàn)?zāi)B(tài)分析過(guò)程中,期望將這些單一的模態(tài)干凈的分離出來(lái),傳統(tǒng)的經(jīng)驗(yàn)?zāi)B(tài)分解方法(emd)由于算法本身的局限性，在分離出來(lái)的平穩(wěn)信號(hào)中會(huì)包含多個(gè)模態(tài)，從而造成模態(tài)混疊，其結(jié)果會(huì)造成頻譜分析的錯(cuò)誤。2014年，周竹生提出了利用小波變換進(jìn)行頻譜分解的分頻avo反演方法。該方法主要通過(guò)連續(xù)小波變換對(duì)疊前角度道集數(shù)據(jù)進(jìn)行頻譜分解,再在頻率域數(shù)據(jù)基礎(chǔ)上進(jìn)行頻散avo屬性計(jì)算。該方法不足之處在于，受限于小波變換自身算法原因，小波變換所獲得頻率域受限于小波窗的分布，往往包含其他頻率域數(shù)據(jù)，故而所得到的時(shí)頻譜在頻率域聚焦方面無(wú)法得到準(zhǔn)確的單一頻率數(shù)據(jù)，其時(shí)頻譜相對(duì)來(lái)說(shuō)過(guò)于粗化，直接影響后期頻散avo屬性計(jì)算的精度。針對(duì)上述缺陷，本專利申請(qǐng)人提供了一種基于ceemd的高精度頻散avo屬性計(jì)算方法，本發(fā)明主要為了解決在實(shí)際研究中，提出運(yùn)用ceemd方法進(jìn)行時(shí)頻分析，能夠較好克服模態(tài)混疊，消除端點(diǎn)效應(yīng)，并且頻率域上的時(shí)頻聚焦性更好，分辨率更高，其得到的頻率域信息的精確性明顯優(yōu)于小波變換，能夠在一定程度優(yōu)化人為后期的頻散avo計(jì)算結(jié)果，提高后期流體識(shí)別的準(zhǔn)確度。
技術(shù)實(shí)現(xiàn)要素：
針對(duì)上述背景技術(shù)存在的上述缺陷，本發(fā)明旨在提供一種基于ceemd的高精度頻散avo屬性計(jì)算方法，本發(fā)明能能夠較好克服模態(tài)混疊，消除端點(diǎn)效應(yīng)，并且頻率域上的時(shí)頻聚焦性更好，分辨率更高，其得到的頻率域信息的精確性明顯優(yōu)于小波變換，能夠在一定程度優(yōu)化人為后期的頻散avo計(jì)算結(jié)果，提高后期流體識(shí)別的準(zhǔn)確度。
為此，本發(fā)明采用以下技術(shù)方案：一種基于ceemd的高精度頻散avo屬性計(jì)算方法，包括以下步驟：
步驟一：首先定義算子ej(·)，當(dāng)給定一個(gè)信號(hào)，通過(guò)emd求得第j個(gè)模態(tài)；wi(n)表示單位方差的零均值高斯白噪聲n(0,1)；i＝1,...,i；εk系數(shù)允許在每個(gè)階段選擇信噪比，設(shè)疊前角度道集數(shù)據(jù)為輸入的目標(biāo)信號(hào)x(t)，使用不同的噪聲通過(guò)emd重復(fù)分解i次，計(jì)算總體平均值，并將其定義為目標(biāo)信號(hào)x(t)的imf1(t)，公式為
步驟二：對(duì)k＝1，計(jì)算一階殘差r1(t)，公式為
r1(t)＝x(t)-imf1(t)
步驟三：emd實(shí)現(xiàn)r1(t)+ε1e1(wi(t))，直到滿足第一個(gè)imf(t)條件，并定義總體平均值為imf2(t)，公式為
步驟四：對(duì)k＝2,...,k，計(jì)算k階殘差rk(t)，公式為
rk(t)＝rk-1(t)-imfk(t)
步驟五：提取rk(t)+εkek(wi(t))的imf1(t)分量，計(jì)算它們的總體平均值得到目標(biāo)信號(hào)的imf(k+1)(t)，公式為
步驟六：重復(fù)步驟(四)(五)，直到殘差不能再被分解為止，則得到最終殘差r(t)為
步驟七：實(shí)際的疊前角度數(shù)據(jù)x(t)包含n個(gè)角度的數(shù)據(jù)，所以x(t)可以按照角度數(shù)量寫(xiě)成x(t,n)，而經(jīng)過(guò)ceemd分解后會(huì)按照實(shí)際數(shù)據(jù)的頻率分量，分解成若干個(gè)imf個(gè)分量，所以同樣按照角度數(shù)量，可以將分解后得到的若干個(gè)imf分量寫(xiě)成s(t,n,f)，即，可以將對(duì)原始疊前角度數(shù)據(jù)進(jìn)行ceemd分解后得到不同頻率分量的imf下振幅譜s：
步驟八：由于地震記錄的振幅信息是地震子波與反射系數(shù)的褶積,振幅譜會(huì)受到“子波疊印”的影響,即能量在各個(gè)頻率分布不均衡,主要集中在主頻帶,因此要對(duì)不同頻率的振幅譜通過(guò)加權(quán)函數(shù)w進(jìn)行譜均衡：
sb(t,n,f)＝s(t,n,f)w(f)
步驟九：對(duì)于某一頻率分量f0，根據(jù)以下頻散公式：
可以得到以下關(guān)系式：
對(duì)上式，采用最小二乘法反演，可以計(jì)算頻譜振幅意義下f0頻率的縱橫波速度變化率；
步驟十：對(duì)步驟九的頻散公式在某一參考頻率f0處對(duì)縱橫波速度泰勒級(jí)數(shù)展開(kāi)，并舍去高階項(xiàng)，只保留一階導(dǎo)數(shù)得到：
其中ia和ib為縱橫波速度變化率關(guān)于頻率f的導(dǎo)數(shù)，即縱橫波速度變化率隨頻率變化的快慢，將其定義為頻散程度：
步驟十一：對(duì)步驟十的泰勒級(jí)數(shù)展開(kāi)后的頻散公式，求ia和ib，可以對(duì)其寫(xiě)為：
考慮m+1個(gè)頻率f1，f2，…，fm的情況，并定義列向量a為：
定義m×n行，2列的矩陣e如下：
將列向量a和矩陣e代入到上式中，可以寫(xiě)為:
于是，每一個(gè)采樣點(diǎn)t處的ia和ib可以通過(guò)最小二乘反演方法求得,最終得到頻散avo屬性：
使用本發(fā)明可以達(dá)到以下有益效果：本發(fā)明充分考慮地震數(shù)據(jù)的非平穩(wěn)性，利用ceemd經(jīng)驗(yàn)?zāi)B(tài)分解方法的算法優(yōu)勢(shì)，ceemd方法主要是通過(guò)向待分析信號(hào)中添加兩個(gè)相反的白噪聲信號(hào)，并分別進(jìn)行emd分解,分解的每一階段添加一個(gè)特定白噪聲，并計(jì)算一個(gè)唯一殘差以得到每個(gè)imf，該方法能夠在保證分解效果與eemd相當(dāng)?shù)那闆r下，減小了由白噪聲引起的重構(gòu)誤差，提供原始信號(hào)的精確重構(gòu)，同時(shí)更好的實(shí)現(xiàn)模態(tài)頻譜分離以及具有更低的計(jì)算成本,能夠有效將非平穩(wěn)信號(hào)分解為若干個(gè)平穩(wěn)信號(hào)，即imf分量。然后在得到頻率域imf分量的基礎(chǔ)上，進(jìn)行頻散avo屬性計(jì)算，這種方法相較于小波變換頻散avo計(jì)算來(lái)說(shuō)，得到的頻率域信息更為準(zhǔn)確，去掉了冗余的其他頻率域信息，以此計(jì)算的頻散avo更符合頻率域?qū)嶋H情況，精度相對(duì)更高。
具體實(shí)施方式
本發(fā)明是基于ceemd的頻散avo屬性計(jì)算方法，與emd、eemd一樣，ceemd可以將一個(gè)非平穩(wěn)信號(hào)分解成有限個(gè)平穩(wěn)信號(hào)，其中的分解得到的每個(gè)平穩(wěn)信號(hào)稱作為固有模態(tài)函數(shù)(intrinsicmodefunction,imf)。同時(shí)，與eemd一樣，ceemd也是一種噪聲輔助分析。首先定義算子ej(·)，當(dāng)給定一個(gè)信號(hào)，通過(guò)emd求得第j個(gè)模態(tài)；wi(n)表示單位方差的零均值高斯白噪聲n(0,1)；i＝1,...,i；εk系數(shù)允許在每個(gè)階段選擇信噪比。設(shè)疊前角度道集數(shù)據(jù)為輸入的目標(biāo)信號(hào)x(t)，則具體的實(shí)現(xiàn)步驟如下：
步驟一：imf1(t)的求取方法與eemd方法相同，即使用不同的噪聲實(shí)現(xiàn)通過(guò)emd重復(fù)分解i次，計(jì)算總體平均值，并將其定義為目標(biāo)信號(hào)x(t)的imf1(t)，公式為
步驟二：對(duì)k＝1，計(jì)算一階殘差r1(t)，公式為
r1(t)＝x(t)-imf1(t)
步驟三：emd實(shí)現(xiàn)r1(t)+ε1e1(wi(t))，直到滿足第一個(gè)imf(t)條件，并定義總體平均值為imf2(t)，公式為
步驟四：對(duì)k＝2,...,k，計(jì)算k階殘差rk(t)，公式為
rk(t)＝rk-1(t)-imfk(t)
步驟五：提取rk(t)+εkek(wi(t))的imf1(t)分量，計(jì)算它們的總體平均值得到目標(biāo)信號(hào)的imf(k+1)(t)，公式為
步驟六：重復(fù)步驟(四)(五)，直到殘差不能再被分解為止，則得到最終殘差r(t)為
步驟七：實(shí)際的疊前角度數(shù)據(jù)x(t)包含n個(gè)角度的數(shù)據(jù)，所以x(t)可以按照角度數(shù)量寫(xiě)成x(t,n)，而經(jīng)過(guò)ceemd分解后會(huì)按照實(shí)際數(shù)據(jù)的頻率分量，分解成若干個(gè)imf個(gè)分量，所以同樣按照角度數(shù)量，可以將分解后得到的若干個(gè)imf分量寫(xiě)成s(t,n,f)，即，可以將對(duì)原始疊前角度數(shù)據(jù)進(jìn)行ceemd分解后得到不同頻率分量的imf下振幅譜s：
步驟八：由于地震記錄的振幅信息是地震子波與反射系數(shù)的褶積,振幅譜會(huì)受到“子波疊印”的影響,即能量在各個(gè)頻率分布不均衡,主要集中在主頻帶附近。因此要對(duì)不同頻率的振幅譜通過(guò)加權(quán)函數(shù)w進(jìn)行譜均衡：
sb(t,n,f)＝s(t,n,f)w(f)
步驟九：對(duì)于某一頻率分量f0，根據(jù)以下頻散公式：
可以得到以下關(guān)系式：
對(duì)上式，采用最小二乘法反演，可以計(jì)算頻譜振幅意義下f0頻率的縱橫波速度變化率。
步驟十：對(duì)步驟九的頻散公式在某一參考頻率f0處對(duì)縱橫波速度泰勒級(jí)數(shù)展開(kāi)，并舍去高階項(xiàng)，只保留一階導(dǎo)數(shù)得到：
其中ia和ib為縱橫波速度變化率關(guān)于頻率f的導(dǎo)數(shù)，即縱橫波速度變化率隨頻率變化的快慢，將其定義為頻散程度：
步驟十一：對(duì)步驟十的泰勒級(jí)數(shù)展開(kāi)后的頻散公式，求ia和ib，可以對(duì)其寫(xiě)為：
考慮m+1個(gè)頻率f1，f2，…，fm的情況，并定義列向量a為：
定義m×n行，2列的矩陣e如下：
將列向量a和矩陣e代入到上式中，可以寫(xiě)為:
于是，每一個(gè)采樣點(diǎn)t處的ia和ib可以通過(guò)最小二乘反演方法求得,最終得到頻散avo屬性：
作為另一種實(shí)施方式，可以利用現(xiàn)有的基于小波變換頻譜分解方法對(duì)疊前角度數(shù)據(jù)進(jìn)行頻率域數(shù)據(jù)分解后，再進(jìn)行頻散avo屬性計(jì)算。
以上顯示和描述了本發(fā)明的基本原理和主要特征和本發(fā)明的優(yōu)點(diǎn)。本行業(yè)的技術(shù)人員應(yīng)該了解，本發(fā)明不受上述實(shí)施例的限制，上述實(shí)施例和說(shuō)明書(shū)中描述的只是說(shuō)明本發(fā)明的原理，在不脫離本發(fā)明精神和范圍的前提下，本發(fā)明還會(huì)有各種變化和改進(jìn)，這些變化和改進(jìn)都落入要求保護(hù)的本發(fā)明范圍內(nèi)。本發(fā)明要求保護(hù)范圍由所附的權(quán)利要求書(shū)及其等效物界定。