劉 霞, 李 文
(東北石油大學 a.電氣信息工程學院;b.土木建筑工程學院, 黑龍江 大慶 163318)
在地震勘探中, 主要通過研究地震波在地層中的傳播規(guī)律, 從而獲取地下的巖性分類以及油氣藏等信息。而檢波器采集到的地震數(shù)據(jù)中包含的噪聲會影響地震資料的反演與解釋, 如何有效抑制噪聲的影響, 一直是地質(zhì)勘探領(lǐng)域的研究熱點和難點問題之一。地震信號是一種非平穩(wěn)信號, 由不同傳播時間的具有不同振幅和頻率的地震子波組成。正是基于地震信號的這一特性, 小波變換[1-3]、曲波變換[4-5]、經(jīng)驗模態(tài)分解(EMD: Empirical Mode Decomposition)[6-10]、時頻峰值濾波算法[11]、局部均值分解[12]和變分模態(tài)分解(VMD: Variational Mode Decomposition)[13-16]等現(xiàn)代信號處理方法可以將地震信號進行頻率分解, 通過去除高頻分量實現(xiàn)信號去噪的目的。在各類算法中如何區(qū)分噪聲和有效信號的頻率分量是去噪算法的關(guān)鍵。小波變換閾值去噪算法是采用將高頻小波系數(shù)與閾值進行比較, 由大于閾值的高頻小波系數(shù)和低頻小波系數(shù)重構(gòu)信號。EMD方法能自適應的將信號分解為一組從高頻到低頻的固有模態(tài)分量, 去除高頻分量, 重構(gòu)低頻分量可實現(xiàn)對信號去噪, 但可能損失高頻中的有效信息[17-20]。楊凱等[17]將EMD方法與小波變換模極大值去噪方法相結(jié)合, 較好地抑制了地震信號中的隨機噪聲。張杏莉等[18]通過計算每個變分模態(tài)分量的能量熵, 通過能量熵的最小值確定噪聲與信號分量的分界, 將信號分量進行重構(gòu)進行去噪。陳輝等[19]通過不同頻帶上分量的梯度值和閾值將各IMF分量分區(qū)進行去噪處理。筆者在對地震信號進行VMD分解過程中發(fā)現(xiàn), 各IMF分量中都會不同程度包含噪聲的影響, 而高頻分量中也會有有效信息。如何去除低頻分量中的噪聲和有效保留高頻分量中的有效信息, 是高分辨率地質(zhì)勘探中要解決的問題。筆者將VMD算法與相關(guān)性和能量熵相結(jié)合, 將地震信號通過VMD分解為若干IMF分量, 然后將各IMF分量進行相關(guān)處理, 區(qū)分出有效信號和噪聲的數(shù)據(jù)點位置, 將去除有效信息的IMF分量分成若干子區(qū)間, 分別計算各子區(qū)間噪聲的能量熵, 選取能量熵最大區(qū)間的IMF分量系數(shù)作為該分量的噪聲方差帶入閾值公式計算得到該分量的閾值, 再把經(jīng)閾值處理后的IMF分量重構(gòu)得到去噪信號。該方法能有效地保留信號中的高頻分量, 從而保留了巖性油氣藏的重要信息。
Dragomiretskiy等[20]在2014年提出了變分模態(tài)分解算法, 該方法將信號分解為具有一定帶寬的K個模態(tài)分量(IMF), 每個模態(tài)分量的中心頻率為ωk。將IMF分量中心頻率和帶寬的估計問題轉(zhuǎn)化為具有約束的變分問題, 通過引入懲罰因子和Lagrange函數(shù)將有約束問題轉(zhuǎn)化為無約束問題, 采用交替乘子算法求最優(yōu)解, 將信號分解為具有一定帶寬的IMF頻率分量, 最后將各IMF頻率分量進行傅里葉逆變換得到時域信號??蓪MD算法看作是維納濾波器組的推廣, 有效解決了EMD算法中的模態(tài)混疊現(xiàn)象。具體方法如下。
1)將信號x(t)分解為K個模態(tài)分量uk(t),k=1,2,…,K,uk(t)定義為一組調(diào)幅-調(diào)頻信號
uk(t)=Ak(t)cos(φk(t))
(1)
(2)
分解的核心問題為找到這K個IMF分量。
2)對每個IMF分量uk(t)進行Hilbert變換, 構(gòu)造解析信號
3)將每個IMF分量的頻譜調(diào)制到以ωk為中心頻率的頻帶上,有
(3)
4)在滿足式(2)的約束條件下, 按照式(3)信號的梯度平方和(L2范數(shù))最小估算各IMF分量的帶寬, 即
5)通過引入懲罰因子α和Lagrange算子λ(t), 將有約束的變分問題轉(zhuǎn)化為無約束的變分問題
6)采用傅里葉變換, 求取二次優(yōu)化問題的頻域解更新公式為
(4)
(5)
(6)
7)對式(4)IMF分量的頻域結(jié)果進行傅里葉逆變換, 得各IMF分量的時域信號。
VMD分解算法是將信號分解為具有一定帶寬、頻率由高到低的K個模態(tài)分量, 而信號中隨機噪聲的頻率較高, 主要分布在高頻部分, 因此將K個模態(tài)分量中的高頻部分去除, 將剩余的低頻IMF分量重構(gòu), 即可得到去噪后的信號
其中l(wèi)+1,…,K個分量為噪聲所引起的。
該去噪算法的關(guān)鍵是找到信號中的有效分量和噪聲分量的分界點l。但去除的高頻分量中可能包含了部分有用信號, 而保留的低頻分量中也可能殘留少量的噪聲。針對該問題, 筆者提出了基于VMD的相關(guān)能量熵閾值去噪算法。
地震信號中的高頻成分可能包含了巖性油氣藏的重要信息, 因此, 在對地震信號進行去噪處理過程中應盡可能的保留這些高頻有效信息。在VMD分解中, 地震信號的有效成分主要分布在低頻分量中, 噪聲主要分布在高頻分量中, 但低頻分量中也殘留一些噪聲, 高頻分量中也含有有效信息。筆者依據(jù)有效信息和噪聲對地震信號的相關(guān)性不同的特性, 對各IMF分量中的有效信息進行檢測與定位, 將有效信息引起的系數(shù)剔除, 將噪聲引起的系數(shù)保留, 得到新的分量系數(shù);這些系數(shù)大部分是由噪聲引起的, 找到這些系數(shù)中噪聲影響最大的區(qū)域, 只要計算出該區(qū)域系數(shù)的能量, 再利用噪聲的能量構(gòu)造閾值, 對各原IMF分量系數(shù)進行閾值處理, 即可達到去除每個IMF分量中噪聲的影響。根據(jù)信息熵理論, 系統(tǒng)越有序, 系統(tǒng)的信息熵就越低;系統(tǒng)越混亂, 系統(tǒng)的信息熵就越高。在一個信號中, 噪聲是無序的, 因此噪聲的信息熵最大, 筆者利用信息熵理論, 將保留的噪聲IMF分量系數(shù)分成若干子區(qū)間, 分別計算各子區(qū)間的噪聲能量熵, 能量熵最大的區(qū)域, 可以認為是全部由噪聲引起的, 利用噪聲能量熵構(gòu)造閾值。該方法可根據(jù)信號的噪聲能量自適應地確定各IMF分量的閾值, 剔除每個IMF分量中噪聲的影響。下面給出改進算法的關(guān)鍵計算步驟。
信號的IMF分量中的有效信息與信號具有較強的相關(guān)性, 而噪聲與信號的相關(guān)性較弱, 根據(jù)上述特點, 通過相關(guān)處理可突出各IMF分量中的有效信息, 弱化噪聲的影響, 實現(xiàn)有效信息的檢測及定位。
定義rk(n)為第k個IMF分量uk(n)與x(n)的相關(guān)系數(shù)
(7)
其中x(n)、uk(n)為信號x(t)和IMF分量uk(t)的離散化表示,n=1,…,N,N為采樣點數(shù)。
為使相關(guān)系數(shù)與IMF分量系數(shù)具有可比性, 定義規(guī)范化相關(guān)系數(shù)Rk(n)為
(8)
通過各IMF分量的規(guī)范化相關(guān)系數(shù)與各IMF分量系數(shù)的比較, 得到有效信號的IMF分量系數(shù)位置。
根據(jù)信號在時頻域具有能量守恒特性, 定義第k個IMF分量的能量
將uk(n)分成l等份, 每個小區(qū)間的能量記為Eki, 采樣點數(shù)為N/l, 則有
第k個IMF分量中第i個子區(qū)間的能量Eki在該分量的總能量Ek中的概率為
pki=Eki/Ek
則第i個子區(qū)間對應的能量熵為
Ski=-pkilnpki
(9)
得到第k個IMF分量的能量熵序列
Sk={Sk1,…,Skl}
搜索第k個IMF分量的熵值最大子區(qū)間, 記為
Skm=max{Sk1,…,Skl}
求第k個IMF分量第m個子區(qū)間的IMF分量系數(shù)的平均值
定義第k個IMF分量的閾值Tk為
(10)
其中σk為第k個IMF分量的熵值最大子區(qū)間的系數(shù)的平均值做為噪聲方差。該閾值的確定方法可根據(jù)信號中噪聲的能量特征自適應地確定各IMF分量的閾值。
Step1 對信號進行VMD分解, 得到各IMF分量。
2)n=n+1;
3)按照式(4)k從1~K迭代更新uk;
4)按照式(5)k從1~K迭代更新ωk;
5)按照式(6)計算λ;
6)重復步驟2)~5), 直到滿足停止條件
(11)
其中ε為分解精度,ε>0;
7)對K個IMF分量進行傅里葉逆變換得到IMF分量的時域分解信號。
Step2 相關(guān)檢測及定位:按照式(7)計算各IMF分量與信號的相關(guān)系數(shù), 按照式(8)計算規(guī)范化相關(guān)系數(shù), 將規(guī)范化相關(guān)系數(shù)大于各IMF分量系數(shù)的位置記錄下來, 這些位置對應了各IMF分量中的有效信號。
Step3 計算噪聲方差:把Step2中確定出的含有有效信號的位置處對應的IMF分量系數(shù)置0, 得到新的系數(shù), 這些新的系數(shù)的采樣點數(shù)保持不變;將得到的新的系數(shù)分成l等份, 分別計算各區(qū)間的能量熵, 并進行比較, 選取能量熵最大的區(qū)間系數(shù), 認為該區(qū)間系數(shù)是由噪聲引起的, 計算該區(qū)間系數(shù)的平均值σk, 作為第k個IMF分量的噪聲方差。
Step4 計算閾值, 按照式(10)計算第k個IMF分量的閾值Tk。
Step5 根據(jù)得到的閾值Tk, 利用軟閾值函數(shù)分別對各IMF分量系數(shù)進行閾值處理, 得到新的IMF分量系數(shù)。
Step6 重構(gòu):利用新的IMF分量系數(shù)重構(gòu)信號, 得到去噪后的信號。
圖1 Ricker子波及加噪信號
我國很多陸相油氣田都屬于薄儲層油氣藏, 為模擬薄儲層油氣藏, 構(gòu)造由2個Ricker子波合成的加噪信號, 分析筆者改進VMD去噪算法的性能, 并與直接去除高頻分量VMD去噪方法進行對比。Ricker子波及加噪信號(信噪比為2 dB)如圖1所示, 對加噪信號進行3層VMD分解的IMF分量如圖2所示, 2種方法的去噪效果對比如圖3所示。不同噪聲水平下的2種方法信噪比(SNR: Signal to Noise Ratio)和均方根誤差(RMSE: Root Mean Squared Error)對比如表1所示。
圖2 Ricker子波3層VMD分解 圖3 去噪效果對比
表1 不同噪聲水平下的去噪效果對比
從圖1原信號和加噪信號的波形可以看出, 第2個Ricker子波分量已大部分淹沒在噪聲中, 從圖2的VMD3層分解的IMF分量中可以看出, 原信號的2個Ricker子波分量和噪聲較好的分離出來, IMF3分量噪聲含量較高, 但前2個IMF分量中也都有噪聲的影響。采用筆者VMD改進算法, 對每層IMF分量都進行閾值處理, 從圖3的2種方法的去噪效果可以看出, 筆者算法重構(gòu)效果優(yōu)于直接去除高頻分量的VMD去噪效果, 去除了大部分噪聲, 較好地保留了2個Ricker子波的形狀;表1中給出了對不同信噪比的信號進行去噪處理的信噪比和均方誤差, 可以看出, 筆者方法在2個指標上都優(yōu)于直接去除高頻分量的VMD去噪效果, 在低信噪比較低時, 具有更好的去噪效果。
由一個低頻Ricker子波和一個高頻Ricker子波疊加構(gòu)造薄儲層楔形地震模型, 如圖4所示。反射層的上層速度為2 000 m/s, 下層速度為3 000 m/s, 深度為100 m, 最小偏移距離10 m, 記錄時間長度512 ms, 密度1, 采集道數(shù)30。對圖4模型加上信噪比為2 dB的高斯白噪聲, 構(gòu)造含噪地震模型如圖5所示。
圖4 地震信號剖面圖 圖5 含噪地震信號剖面圖
通過圖6采用直接去除高頻分量的VMD方法去噪效果中可以看出, 地震剖面中還殘留了部分噪聲, 去噪后信噪比為2.92 dB。從圖7中可以看出, 筆者方法很好的抑制了噪聲, 同相軸保留完整, 分辨率更好, 有效信號得到更好的保護, 去噪后信噪比為6.25 dB, 去噪效果更好。
圖6 傳統(tǒng)VMD去噪后的地震剖面 圖7 筆者方法去噪后的地震剖面
從某地區(qū)實際地震資料中選取100道信號, 每道采樣點為1 024, 采樣間隔為2 ms, 地震剖面如圖8所示。采用直接去除高頻分量的VMD方法和筆者方法進行去噪處理, 去噪剖面圖如圖9、圖10所示。
圖8 實際地震剖面圖 圖9 VMD方法去噪后地震剖面 圖10 筆者方法去噪后
從圖8的實際地震信號中可以看出, 剖面圖中含有大量的噪聲, 同相軸不清晰。通過圖9和圖10可以看到, 兩種方法都抑制了噪聲的影響。但直接去除高頻分量的VMD方法削弱了同相軸, 使有效信息損失較大, 而筆者方法去噪后的地震剖面同相軸連續(xù)、清晰, 很好的保留了有效信息。
筆者針對直接去除高頻分量的VMD去噪方法中, 不能有效區(qū)分各IMF分量中的有效信息和噪聲問題, 提出了基于變分模態(tài)分解的相關(guān)能量熵自適應閾值去噪方法, 該方法通過相關(guān)處理檢測有效信息并定位, 根據(jù)噪聲的能量熵自適應地確定各IMF分量的閾值。通過對不同噪聲水平下的地震模型測試和實際地震信號去噪處理表明, 筆者方法能更有效的從強噪聲背景下提取有效信息, 去除混疊噪聲的影響, 提高信噪比, 為地震資料的解釋等工作提供高分辨率地震數(shù)據(jù)。