基于變分模態(tài)分解的相關(guān)能量熵自適應(yīng)閾值去噪

劉 霞, 李 文

(東北石油大學(xué) a.電氣信息工程學(xué)院；b.土木建筑工程學(xué)院, 黑龍江 大慶 163318)

0 引 言

在地震勘探中, 主要通過(guò)研究地震波在地層中的傳播規(guī)律, 從而獲取地下的巖性分類以及油氣藏等信息。而檢波器采集到的地震數(shù)據(jù)中包含的噪聲會(huì)影響地震資料的反演與解釋, 如何有效抑制噪聲的影響, 一直是地質(zhì)勘探領(lǐng)域的研究熱點(diǎn)和難點(diǎn)問(wèn)題之一。地震信號(hào)是一種非平穩(wěn)信號(hào), 由不同傳播時(shí)間的具有不同振幅和頻率的地震子波組成。正是基于地震信號(hào)的這一特性, 小波變換[1-3]、曲波變換[4-5]、經(jīng)驗(yàn)?zāi)B(tài)分解(EMD: Empirical Mode Decomposition)[6-10]、時(shí)頻峰值濾波算法[11]、局部均值分解[12]和變分模態(tài)分解(VMD: Variational Mode Decomposition)[13-16]等現(xiàn)代信號(hào)處理方法可以將地震信號(hào)進(jìn)行頻率分解, 通過(guò)去除高頻分量實(shí)現(xiàn)信號(hào)去噪的目的。在各類算法中如何區(qū)分噪聲和有效信號(hào)的頻率分量是去噪算法的關(guān)鍵。小波變換閾值去噪算法是采用將高頻小波系數(shù)與閾值進(jìn)行比較, 由大于閾值的高頻小波系數(shù)和低頻小波系數(shù)重構(gòu)信號(hào)。EMD方法能自適應(yīng)的將信號(hào)分解為一組從高頻到低頻的固有模態(tài)分量, 去除高頻分量, 重構(gòu)低頻分量可實(shí)現(xiàn)對(duì)信號(hào)去噪, 但可能損失高頻中的有效信息[17-20]。楊凱等[17]將EMD方法與小波變換模極大值去噪方法相結(jié)合, 較好地抑制了地震信號(hào)中的隨機(jī)噪聲。張杏莉等[18]通過(guò)計(jì)算每個(gè)變分模態(tài)分量的能量熵, 通過(guò)能量熵的最小值確定噪聲與信號(hào)分量的分界, 將信號(hào)分量進(jìn)行重構(gòu)進(jìn)行去噪。陳輝等[19]通過(guò)不同頻帶上分量的梯度值和閾值將各IMF分量分區(qū)進(jìn)行去噪處理。筆者在對(duì)地震信號(hào)進(jìn)行VMD分解過(guò)程中發(fā)現(xiàn), 各IMF分量中都會(huì)不同程度包含噪聲的影響, 而高頻分量中也會(huì)有有效信息。如何去除低頻分量中的噪聲和有效保留高頻分量中的有效信息, 是高分辨率地質(zhì)勘探中要解決的問(wèn)題。筆者將VMD算法與相關(guān)性和能量熵相結(jié)合, 將地震信號(hào)通過(guò)VMD分解為若干IMF分量, 然后將各IMF分量進(jìn)行相關(guān)處理, 區(qū)分出有效信號(hào)和噪聲的數(shù)據(jù)點(diǎn)位置, 將去除有效信息的IMF分量分成若干子區(qū)間, 分別計(jì)算各子區(qū)間噪聲的能量熵, 選取能量熵最大區(qū)間的IMF分量系數(shù)作為該分量的噪聲方差帶入閾值公式計(jì)算得到該分量的閾值, 再把經(jīng)閾值處理后的IMF分量重構(gòu)得到去噪信號(hào)。該方法能有效地保留信號(hào)中的高頻分量, 從而保留了巖性油氣藏的重要信息。

1 VMD去噪原理

1.1 變分模態(tài)分解原理

Dragomiretskiy等[20]在2014年提出了變分模態(tài)分解算法, 該方法將信號(hào)分解為具有一定帶寬的K個(gè)模態(tài)分量(IMF), 每個(gè)模態(tài)分量的中心頻率為ωk。將IMF分量中心頻率和帶寬的估計(jì)問(wèn)題轉(zhuǎn)化為具有約束的變分問(wèn)題, 通過(guò)引入懲罰因子和Lagrange函數(shù)將有約束問(wèn)題轉(zhuǎn)化為無(wú)約束問(wèn)題, 采用交替乘子算法求最優(yōu)解, 將信號(hào)分解為具有一定帶寬的IMF頻率分量, 最后將各IMF頻率分量進(jìn)行傅里葉逆變換得到時(shí)域信號(hào)?？蓪MD算法看作是維納濾波器組的推廣, 有效解決了EMD算法中的模態(tài)混疊現(xiàn)象。具體方法如下。

1)將信號(hào)x(t)分解為K個(gè)模態(tài)分量uk(t),k=1,2,…,K,uk(t)定義為一組調(diào)幅-調(diào)頻信號(hào)

uk(t)=Ak(t)cos(φk(t))

(1)

(2)

分解的核心問(wèn)題為找到這K個(gè)IMF分量。

2)對(duì)每個(gè)IMF分量uk(t)進(jìn)行Hilbert變換, 構(gòu)造解析信號(hào)

3)將每個(gè)IMF分量的頻譜調(diào)制到以ωk為中心頻率的頻帶上,有

(3)

4)在滿足式(2)的約束條件下, 按照式(3)信號(hào)的梯度平方和(L2范數(shù))最小估算各IMF分量的帶寬, 即

5)通過(guò)引入懲罰因子α和Lagrange算子λ(t), 將有約束的變分問(wèn)題轉(zhuǎn)化為無(wú)約束的變分問(wèn)題

6)采用傅里葉變換, 求取二次優(yōu)化問(wèn)題的頻域解更新公式為

(4)

(5)

(6)

7)對(duì)式(4)IMF分量的頻域結(jié)果進(jìn)行傅里葉逆變換, 得各IMF分量的時(shí)域信號(hào)。

1.2 VMD去噪算法

VMD分解算法是將信號(hào)分解為具有一定帶寬、頻率由高到低的K個(gè)模態(tài)分量, 而信號(hào)中隨機(jī)噪聲的頻率較高, 主要分布在高頻部分, 因此將K個(gè)模態(tài)分量中的高頻部分去除, 將剩余的低頻IMF分量重構(gòu), 即可得到去噪后的信號(hào)

其中l(wèi)+1,…,K個(gè)分量為噪聲所引起的。

該去噪算法的關(guān)鍵是找到信號(hào)中的有效分量和噪聲分量的分界點(diǎn)l。但去除的高頻分量中可能包含了部分有用信號(hào), 而保留的低頻分量中也可能殘留少量的噪聲。針對(duì)該問(wèn)題, 筆者提出了基于VMD的相關(guān)能量熵閾值去噪算法。

2 基于VMD的相關(guān)能量熵閾值去噪算法

2.1 算法思想

地震信號(hào)中的高頻成分可能包含了巖性油氣藏的重要信息, 因此, 在對(duì)地震信號(hào)進(jìn)行去噪處理過(guò)程中應(yīng)盡可能的保留這些高頻有效信息。在VMD分解中, 地震信號(hào)的有效成分主要分布在低頻分量中, 噪聲主要分布在高頻分量中, 但低頻分量中也殘留一些噪聲, 高頻分量中也含有有效信息。筆者依據(jù)有效信息和噪聲對(duì)地震信號(hào)的相關(guān)性不同的特性, 對(duì)各IMF分量中的有效信息進(jìn)行檢測(cè)與定位, 將有效信息引起的系數(shù)剔除, 將噪聲引起的系數(shù)保留, 得到新的分量系數(shù)；這些系數(shù)大部分是由噪聲引起的, 找到這些系數(shù)中噪聲影響最大的區(qū)域, 只要計(jì)算出該區(qū)域系數(shù)的能量, 再利用噪聲的能量構(gòu)造閾值, 對(duì)各原IMF分量系數(shù)進(jìn)行閾值處理, 即可達(dá)到去除每個(gè)IMF分量中噪聲的影響。根據(jù)信息熵理論, 系統(tǒng)越有序, 系統(tǒng)的信息熵就越低；系統(tǒng)越混亂, 系統(tǒng)的信息熵就越高。在一個(gè)信號(hào)中, 噪聲是無(wú)序的, 因此噪聲的信息熵最大, 筆者利用信息熵理論, 將保留的噪聲IMF分量系數(shù)分成若干子區(qū)間, 分別計(jì)算各子區(qū)間的噪聲能量熵, 能量熵最大的區(qū)域, 可以認(rèn)為是全部由噪聲引起的, 利用噪聲能量熵構(gòu)造閾值。該方法可根據(jù)信號(hào)的噪聲能量自適應(yīng)地確定各IMF分量的閾值, 剔除每個(gè)IMF分量中噪聲的影響。下面給出改進(jìn)算法的關(guān)鍵計(jì)算步驟。

2.2 相關(guān)性檢測(cè)及定位

信號(hào)的IMF分量中的有效信息與信號(hào)具有較強(qiáng)的相關(guān)性, 而噪聲與信號(hào)的相關(guān)性較弱, 根據(jù)上述特點(diǎn), 通過(guò)相關(guān)處理可突出各IMF分量中的有效信息, 弱化噪聲的影響, 實(shí)現(xiàn)有效信息的檢測(cè)及定位。

定義rk(n)為第k個(gè)IMF分量uk(n)與x(n)的相關(guān)系數(shù)

(7)

其中x(n)、uk(n)為信號(hào)x(t)和IMF分量uk(t)的離散化表示,n=1,…,N,N為采樣點(diǎn)數(shù)。

為使相關(guān)系數(shù)與IMF分量系數(shù)具有可比性, 定義規(guī)范化相關(guān)系數(shù)Rk(n)為

(8)

通過(guò)各IMF分量的規(guī)范化相關(guān)系數(shù)與各IMF分量系數(shù)的比較, 得到有效信號(hào)的IMF分量系數(shù)位置。

2.3 噪聲能量熵計(jì)算

根據(jù)信號(hào)在時(shí)頻域具有能量守恒特性, 定義第k個(gè)IMF分量的能量

將uk(n)分成l等份, 每個(gè)小區(qū)間的能量記為Eki, 采樣點(diǎn)數(shù)為N/l, 則有

第k個(gè)IMF分量中第i個(gè)子區(qū)間的能量Eki在該分量的總能量Ek中的概率為

pki=Eki/Ek

則第i個(gè)子區(qū)間對(duì)應(yīng)的能量熵為

Ski=-pkilnpki

(9)

得到第k個(gè)IMF分量的能量熵序列

Sk={Sk1,…,Skl}

2.4 閾值計(jì)算

搜索第k個(gè)IMF分量的熵值最大子區(qū)間, 記為

Skm=max{Sk1,…,Skl}

求第k個(gè)IMF分量第m個(gè)子區(qū)間的IMF分量系數(shù)的平均值

定義第k個(gè)IMF分量的閾值Tk為

(10)

其中σk為第k個(gè)IMF分量的熵值最大子區(qū)間的系數(shù)的平均值做為噪聲方差。該閾值的確定方法可根據(jù)信號(hào)中噪聲的能量特征自適應(yīng)地確定各IMF分量的閾值。

2.5 筆者方法實(shí)現(xiàn)步驟

Step1 對(duì)信號(hào)進(jìn)行VMD分解, 得到各IMF分量。

2)n=n+1;

3)按照式(4)k從1～K迭代更新uk；

4)按照式(5)k從1～K迭代更新ωk；

5)按照式(6)計(jì)算λ；

6)重復(fù)步驟2)～5), 直到滿足停止條件

(11)

其中ε為分解精度,ε>0；

7)對(duì)K個(gè)IMF分量進(jìn)行傅里葉逆變換得到IMF分量的時(shí)域分解信號(hào)。

Step2 相關(guān)檢測(cè)及定位：按照式(7)計(jì)算各IMF分量與信號(hào)的相關(guān)系數(shù), 按照式(8)計(jì)算規(guī)范化相關(guān)系數(shù), 將規(guī)范化相關(guān)系數(shù)大于各IMF分量系數(shù)的位置記錄下來(lái), 這些位置對(duì)應(yīng)了各IMF分量中的有效信號(hào)。

Step3 計(jì)算噪聲方差：把Step2中確定出的含有有效信號(hào)的位置處對(duì)應(yīng)的IMF分量系數(shù)置0, 得到新的系數(shù), 這些新的系數(shù)的采樣點(diǎn)數(shù)保持不變；將得到的新的系數(shù)分成l等份, 分別計(jì)算各區(qū)間的能量熵, 并進(jìn)行比較, 選取能量熵最大的區(qū)間系數(shù), 認(rèn)為該區(qū)間系數(shù)是由噪聲引起的, 計(jì)算該區(qū)間系數(shù)的平均值σk, 作為第k個(gè)IMF分量的噪聲方差。

Step4 計(jì)算閾值, 按照式(10)計(jì)算第k個(gè)IMF分量的閾值Tk。

Step5 根據(jù)得到的閾值Tk, 利用軟閾值函數(shù)分別對(duì)各IMF分量系數(shù)進(jìn)行閾值處理, 得到新的IMF分量系數(shù)。

Step6 重構(gòu)：利用新的IMF分量系數(shù)重構(gòu)信號(hào), 得到去噪后的信號(hào)。

3 數(shù)值仿真

3.1 Ricker子波模擬地震信號(hào)仿真

圖1 Ricker子波及加噪信號(hào)

我國(guó)很多陸相油氣田都屬于薄儲(chǔ)層油氣藏, 為模擬薄儲(chǔ)層油氣藏, 構(gòu)造由2個(gè)Ricker子波合成的加噪信號(hào), 分析筆者改進(jìn)VMD去噪算法的性能, 并與直接去除高頻分量VMD去噪方法進(jìn)行對(duì)比。Ricker子波及加噪信號(hào)(信噪比為2 dB)如圖1所示, 對(duì)加噪信號(hào)進(jìn)行3層VMD分解的IMF分量如圖2所示, 2種方法的去噪效果對(duì)比如圖3所示。不同噪聲水平下的2種方法信噪比(SNR: Signal to Noise Ratio)和均方根誤差(RMSE: Root Mean Squared Error)對(duì)比如表1所示。

圖2 Ricker子波3層VMD分解 圖3 去噪效果對(duì)比

表1 不同噪聲水平下的去噪效果對(duì)比

從圖1原信號(hào)和加噪信號(hào)的波形可以看出, 第2個(gè)Ricker子波分量已大部分淹沒(méi)在噪聲中, 從圖2的VMD3層分解的IMF分量中可以看出, 原信號(hào)的2個(gè)Ricker子波分量和噪聲較好的分離出來(lái), IMF3分量噪聲含量較高, 但前2個(gè)IMF分量中也都有噪聲的影響。采用筆者VMD改進(jìn)算法, 對(duì)每層IMF分量都進(jìn)行閾值處理, 從圖3的2種方法的去噪效果可以看出, 筆者算法重構(gòu)效果優(yōu)于直接去除高頻分量的VMD去噪效果, 去除了大部分噪聲, 較好地保留了2個(gè)Ricker子波的形狀；表1中給出了對(duì)不同信噪比的信號(hào)進(jìn)行去噪處理的信噪比和均方誤差, 可以看出, 筆者方法在2個(gè)指標(biāo)上都優(yōu)于直接去除高頻分量的VMD去噪效果, 在低信噪比較低時(shí), 具有更好的去噪效果。

3.2 合成地震信號(hào)去噪效果對(duì)比

由一個(gè)低頻Ricker子波和一個(gè)高頻Ricker子波疊加構(gòu)造薄儲(chǔ)層楔形地震模型, 如圖4所示。反射層的上層速度為2 000 m/s, 下層速度為3 000 m/s, 深度為100 m, 最小偏移距離10 m, 記錄時(shí)間長(zhǎng)度512 ms, 密度1, 采集道數(shù)30。對(duì)圖4模型加上信噪比為2 dB的高斯白噪聲, 構(gòu)造含噪地震模型如圖5所示。

圖4 地震信號(hào)剖面圖 圖5 含噪地震信號(hào)剖面圖

通過(guò)圖6采用直接去除高頻分量的VMD方法去噪效果中可以看出, 地震剖面中還殘留了部分噪聲, 去噪后信噪比為2.92 dB。從圖7中可以看出, 筆者方法很好的抑制了噪聲, 同相軸保留完整, 分辨率更好, 有效信號(hào)得到更好的保護(hù), 去噪后信噪比為6.25 dB, 去噪效果更好。

圖6 傳統(tǒng)VMD去噪后的地震剖面 圖7 筆者方法去噪后的地震剖面

4 實(shí)際地震信號(hào)去噪

從某地區(qū)實(shí)際地震資料中選取100道信號(hào), 每道采樣點(diǎn)為1 024, 采樣間隔為2 ms, 地震剖面如圖8所示。采用直接去除高頻分量的VMD方法和筆者方法進(jìn)行去噪處理, 去噪剖面圖如圖9、圖10所示。

圖8 實(shí)際地震剖面圖 圖9 VMD方法去噪后地震剖面 圖10 筆者方法去噪后

從圖8的實(shí)際地震信號(hào)中可以看出, 剖面圖中含有大量的噪聲, 同相軸不清晰。通過(guò)圖9和圖10可以看到, 兩種方法都抑制了噪聲的影響。但直接去除高頻分量的VMD方法削弱了同相軸, 使有效信息損失較大, 而筆者方法去噪后的地震剖面同相軸連續(xù)、清晰, 很好的保留了有效信息。

5 結(jié) 語(yǔ)

筆者針對(duì)直接去除高頻分量的VMD去噪方法中, 不能有效區(qū)分各IMF分量中的有效信息和噪聲問(wèn)題, 提出了基于變分模態(tài)分解的相關(guān)能量熵自適應(yīng)閾值去噪方法, 該方法通過(guò)相關(guān)處理檢測(cè)有效信息并定位, 根據(jù)噪聲的能量熵自適應(yīng)地確定各IMF分量的閾值。通過(guò)對(duì)不同噪聲水平下的地震模型測(cè)試和實(shí)際地震信號(hào)去噪處理表明, 筆者方法能更有效的從強(qiáng)噪聲背景下提取有效信息, 去除混疊噪聲的影響, 提高信噪比, 為地震資料的解釋等工作提供高分辨率地震數(shù)據(jù)。