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