基于ICEEMDAN和分布熵的SS-Y伸縮儀信號隨機噪聲壓制方法

吳林斌

1 中國地震局地震研究所,武漢市洪山側(cè)路40號, 430071

2 武漢地震科學儀器研究院有限公司,武漢市洪山側(cè)路40號,430071

SS-Y伸縮儀是由中國地震局地震研究所自主研制,并在大地測量領域廣泛應用的地殼表面兩點因變量觀測儀器,具有遙控標定、高精度、短基線等特點[1-2]。SS-Y伸縮儀映震能力較好,但易受到前置放大故障、供電故障、超量程等儀器內(nèi)部故障或氣壓、溫度、雷電等復雜外部環(huán)境因素的影響[3-8],出現(xiàn)數(shù)據(jù)曲線固體潮不清晰、局部不光滑等情況。自適應噪聲完備集合經(jīng)驗模態(tài)分解(CEEMDAN)在信號去噪領域已取得顯著成果。該方法通過自主分析信號質(zhì)量,在信號分解過程中引入相匹配的白噪聲,以改善傳統(tǒng)EMD過程中常見的模態(tài)混疊現(xiàn)像,信號重構質(zhì)量較高且去噪效果較好。基于CEEMDAN的改進的自適應噪聲完備集合經(jīng)驗模態(tài)分解(ICEEMDAN)[9-11]能有效減少模態(tài)混疊現(xiàn)象和分解后虛假分量的產(chǎn)生,且針對信號在CEEMDAN處理過程中的問題進行優(yōu)化改進,克服了某些本征模態(tài)函數(shù)中仍含有較多殘余分量的現(xiàn)狀,信號重構質(zhì)量得到明顯提升。

本文提出一種基于ICEEMDAN和分布熵(DistEn)[12]的SS-Y伸縮儀信號隨機噪聲壓制方法。根據(jù)ICEEMDAN得到的各分量整體混亂程度,判斷其是否為原始信號的有效成分或是否為含噪分量,進而采取處置手段;同時,設計仿真信號與SS-Y伸縮儀實測信號去噪對比及單輪與多輪去噪對比實驗,結合多個指標評價去噪模型的性能。

1 信號去噪模型

1.1 ICEEMDAN原理

ICEEMDAN的主要步驟為[9-11]:

1)對原始信號x進行分解,構造原始信號的加噪合成信號:

x(i)=x+β0E1(w(i))

(1)

式中,β0為首次分解時的信噪比,w(i)為第i個被添加的高斯白噪聲,E1(·)為EMD的首個IMF分量。

2)重復步驟1)N次,得到N個x(i),并求全部局部均值函數(shù)值的均值,即首個殘余分量:

(2)

式中,M(·)為局部均值函數(shù)。

3)在原始信號的基礎上減去首個殘余分量,得到首個模態(tài):

d1=x-r1

(3)

4)根據(jù)第k次分解后計算的局部均值函數(shù)值,得到第k個殘余分量,并用上一個殘余分量減去該殘差結果,得到第k個模態(tài)dm:

dm=rk-1-rk

(4)

(5)

5)持續(xù)執(zhí)行步驟4),直至滿足分解停止條件。

1.2 DistEn原理

DistEn[12]是對樣本熵和近似熵的改進,可解決部分參數(shù)的高度依賴和非魯棒性問題,一般計算步驟為:

(6)

3)利用直方圖估算經(jīng)驗概率密度函數(shù)(EPDF,記作fk)。定義直方圖柱數(shù)q,則直方圖頻數(shù)分布為N=[n1,n2,…,nk,…,nq],距離矩陣D的EPDF為:

(7)

4)計算分布熵:

(8)

1.3 ICEEMDAN-DistEn去噪模型

ICEEMDAN-DistEn去噪模型步驟如下:

1)對SS-Y伸縮儀信號進行ICEEMDAN處理,得到若干個IMF分量,并依據(jù)頻率高低進行排序和標記。

2)分別計算各IMF分量的DistEn,并設定閾值δ,對符合不同閾值區(qū)間的IMF分量作對應的差異化處理。

3)當某個IMF的DistEn∈[δ,1]時,表示該IMF分量內(nèi)部幾乎不包含隨機噪聲成分,屬于原始信號的有效成分,可以直接參與線性重構。

4)當某個IMF的DistEn∈[0,δ]時,表示該IMF分量基本為隨機噪聲成分,具有類似的隨機振蕩特性,可以直接舍棄。

5)將步驟3)和4)所保留的全部IMF分量進行線性重構,即可得到SS-Y伸縮儀去噪后的數(shù)據(jù)。

由于ICEEMDAN方法屬于CEEMDAN方法的改進,參考ICEEMDAN-DistEn去噪模型的核心原理得到CEEMDAN-DistEn去噪模型,作為SS-Y伸縮儀信號去噪實驗的對照模型之一。

2 實驗與分析

設計3個實驗分析ICEEMDAN-DistEn去噪模型的有效性,包括仿真信號的高、低頻噪聲去噪實驗和SS-Y伸縮儀實測信號去噪實驗,并利用相關系數(shù)、均方誤差(MSE)及信噪比(SNR)等評價指標衡量模型的去噪效果。

2.1 仿真信號的高頻噪聲去噪實驗與分析

設計仿真信號去噪實驗,初始化一段純凈信號x1(t)和兩段位置不同且幅值差異較大的噪聲成分x2(t)和x3(t)(圖1),線性合成后構成需要處理的仿真含噪信號X(t),采樣率為10 Hz。ICEEMDAN及ICEEMDAN-DistEn去噪波形如圖2和3所示:

圖1 利用MATLAB生成仿真純凈信號與加噪信號

圖2 仿真含噪信號經(jīng)過ICEEMDAN處理后得到的IMF

圖3 基于ICEEMDAN-DistEn模型的仿真信號去噪結果

X(t)=x1(t)+x2(t)+x3(t)

(9)

x1(t)=sin(0.08πt),1≤t≤1 000

(10)

x2(t)=sin(0.5π(t-201)),201≤t≤300

(11)

x3(t)=sin(0.8π(t-301)),501≤t≤700

(12)

圖2為仿真含噪信號的ICEEMDAN處理結果,包含8個IMF分量。可以看出,IMF1分量即可從原始含噪信號中分離出2個發(fā)生位置及持續(xù)時長完全不同的噪聲成分,而IMF2分量和 IMF3分量是原始純凈信號的主要幅值成分,其他IMF分量幅值較小,對去噪結果影響不大。

表1為原始含噪信號經(jīng)過ICEEMDAN處理后各IMF分量的分布熵值。可以看出,各IMF分量在最大幅值、混亂程度、局部變異、對稱性等多個維度均存在極大差異,其中IMF1的分布熵值最小,IMF7的分布熵值最大,說明IMF1的局部噪聲成分具有顯著規(guī)律性和周期性,很容易被ICEEMDAN算法識別出。結合圖2和表1可以發(fā)現(xiàn),有規(guī)律性的高頻噪聲成分分布熵值較小,而作為原始含噪信號的有效成分由于已被ICEEMDAN處理為不同頻率下的IMF分量,其分布熵值均較大,均可被用于線性重構。

表1 仿真含噪信號經(jīng)過ICEEMDAN處理后得到的各IMF分量的分布熵值

結合圖1和3可知,基于仿真含噪信號的ICEEMDAN-DistEn模型去噪效果較顯著。為證明該模型的優(yōu)越性,對比分析ICEEMDAN-DistEn、CEEMDAN-DistEn、卡爾曼濾波、小波去噪、70階低通FIR濾波和Savitzky-Golay濾波等6個模型的去噪效果。表2為6種去噪模型的含噪信號處理結果,表中相關系數(shù)、MSE、SNR等評價指標基于原始純凈信號和去噪后信號進行比較。

表2 不同去噪模型下的仿真含噪信號去噪結果

由表可知,原始含噪信號相對于純凈信號的信噪比為5.228 8,6個模型中ICEEMDAN-DistEn、小波去噪、70階低通FIR濾波和Savitzky-Golay濾波的去噪效果明顯,SNR值大幅提高;綜合相關系數(shù)和MSE來看,70階低通FIR濾波模型的信號重構損失稍大,只有ICEEMDAN-DistEn模型的各指標均有相對均衡的去噪效果;卡爾曼濾波受制于狀態(tài)方程參數(shù)設置的不確定性和隨機性,去噪結果變化較大,本實驗中沒有得到預期效果。

2.2 仿真信號的低頻噪聲去噪實驗與分析

設計仿真含噪信號的相對低頻噪聲剔除實驗(圖4),如式(13)～(16)所示,低頻噪聲信號的頻率分別為0.005 Hz和0.002 5 Hz,遠低于純凈信號的0.04 Hz:

圖4 利用MATLAB生成仿真純凈信號與低頻噪聲

X′(t)=x′1(t)+x′2(t)+x′3(t)

(13)

x′1(t)=sin(0.08πt),

1≤t≤1 000

(14)

x′2(t)=sin(0.01π(t-201)),

201≤t≤300

(15)

x′3(t)=sin(0.005π(t-301)),

501≤t≤700

(16)

圖5為含低頻噪聲信號的ICEEMDAN處理結果。可以看出,IMF2分量與原始純凈信號的波形特征和幅值最為接近;IMF1分量的幅值整體較小且呈現(xiàn)無規(guī)律的尖刺現(xiàn)象,可判斷為噪聲成分;其他IMF分量的曲線較為光滑,需借助評價指標來判斷優(yōu)劣程度。由表3可知,各IMF按照分布熵數(shù)值降序排列依次為IMF7、IMF8、IMF6、IMF2、IMF5、IMF3、IMF4、IMF1,分布熵值大于0.9的IMF分量有4個,含有極少的噪聲,可以參與到后續(xù)的信號重構中。最終得到的去噪結果如圖6所示。

圖5 含低頻噪聲信號經(jīng)過ICEEMDAN處理后得到的IMF

圖6 基于ICEEMDAN-DistEn模型的仿真低頻噪聲信號去噪結果

表4為利用ICEEMDAN-DistEn、CEEMDAN-DistEn、卡爾曼濾波、小波去噪及Savitzky-Golay濾波等5個模型的去噪結果。可以看出,ICEEMDAN-DistEn模型在相關系數(shù)、MSE及SNR上皆顯著優(yōu)于其他去噪模型,能較好地抑制原始含噪信號中的低頻噪聲成分,一定程度上還原了初始純凈信號。

表4 不同去噪模型下的仿真含低頻噪聲信號去噪結果

2.3 SS-Y伸縮儀實測信號的去噪實驗與分析

應用實測含噪信號進行ICEEMDAN-DistEn去噪實驗,數(shù)據(jù)來自SS-Y伸縮儀,地點位于武漢地震中心站,時間為2022-05-31,文件格式為EPD,模型測試環(huán)境為MATLAB 2023a,實驗結果如圖7～9和表5～6所示。

表5 SS-Y伸縮儀信號經(jīng)過ICEEMDAN處理后得到的各IMF分布熵值

圖7 SS-Y伸縮儀信號經(jīng)過ICEEMDAN處理后得到的IMF

圖7為SS-Y伸縮儀信號經(jīng)ICEEMDAN處理后得到的IMF結果。可以看出,IMF6分量擁有初始信號的主要能量;IMF1～IMF4分量的能量幅值極小,且波形極為雜亂,在零值上下劇烈振蕩;IMF5分量的波形較光滑,最大幅值相對較小,對信號重構的影響不顯著。表5為各IMF分量的分布熵值,其中IMF4～IMF6分量的分布熵值皆超過0.8,而IMF1～IMF3分量分布熵值差異明顯。

圖8為SS-Y伸縮儀單條信號的ICEEMDAN-DisEn去噪結果。從局部放大結果可以看出,該模型能有效地光滑曲線,且具有極好的信號還原性,線性重構結果較理想。引入CEEMDAN-DistEn、卡爾曼濾波、70階低通FIR濾波等模型進行信號去噪對比,表6為不同模型的信號去噪效果。可以看出,ICEEMDAN-DistEn模型的去噪效果相對最好;雖然CEEMDAN-DistEn模型的MSE值僅有0.800 5,表明其信號損失極小,但其SNR和相關系數(shù)都不理想,說明去噪過程中未處理真正的噪聲成分;卡爾曼濾波和70階低通FIR濾波模型的信噪比較高,但都出現(xiàn)極大的重構波形變異,結果完全失真。

表6 不同去噪模型下SS-Y伸縮儀信號去噪結果

圖8 基于ICEEMDAN-DisEn模型的單條SS-Y伸縮儀信號去噪結果

2022-05的SS-Y伸縮儀信號去噪結果如圖9所示,其中柱狀圖為相關系數(shù)的分布統(tǒng)計,折線圖為信噪比曲線。可以看出,幾乎全部信號的信噪比都超過了70 dB,去噪前后的相關系數(shù)值也都接近于1。

圖9 基于ICEEMDAN-DistEn模型的30 d SS-Y伸縮儀信號去噪結果

3 結 語

本文提出一種結合ICEEMDAN和分布熵的SS-Y伸縮儀信號去噪方法,能夠更精確地識別出噪聲和有效信號,線性重構誤差較小,無明顯的失真現(xiàn)象,綜合性能優(yōu)于CEEMDAN和分布熵組建的去噪模型,并得到以下結論:

1)作為2種完全不同的信號處理技術,ICEEMDAN和分布熵能分別提取和篩選出含噪信號中的噪聲成分和有效成分,適當?shù)貙S-Y伸縮儀含噪信號進行分解及熵值計算,剔除或保留信號并進行線性重構,可以有效地對信號進行去噪處理。

2)ICEEMDAN-DisEn模型有較明顯的曲線光滑效果,不會造成嚴重的信號失真,同時可以將各種噪聲的能量和位置準確地分離出來,便于后續(xù)開展SS-Y伸縮儀信號的異常檢測,并與其他形變觀測儀器進行同震響應分析。此外,本文設計的仿真信號去噪實驗證明該模型既能剔除高頻噪聲,也能有效抑制低頻噪聲,是性能較為全面的去噪模型。

3)由于觀測環(huán)境和儀器布設地點的影響,ICEEMDAN-DisEn模型中的分布熵閾值并非固定值,應根據(jù)待去噪樣本的數(shù)據(jù)特性預先設定一個值并檢驗其合理性,但閾值的選取原則應保持不變:較大分布熵值指示的IMF分量是有效成分的概率較大,較小分布熵值指示的IMF分量是噪聲的概率較大。

針對SS-Y伸縮儀信號去噪結果的可靠性檢驗仍有一定的改進空間。后續(xù)將從固體潮調(diào)和分析入手,通過計算O1波和M2波的振幅因子和相位滯后等參數(shù)來佐證去噪結果的合理性。