一種新小波閾值算法的變形監測數據去噪方法

沙 磊

1 安徽省地質實驗研究所(國土資源部合肥礦產資源監督檢測中心),合肥市阜陽北路318號,230041

變形監測數據一般包含真值和噪聲2個部分,對變形監測數據進行去噪處理、提取變形特征、分析變形規律,是變形分析研究的重要內容[1]。小波分析具有良好的時頻特性,被廣泛應用于變形分析中,特別是小波閾值去噪在非平穩變形監測數據處理中取得了較好的效果[2-3]。眾多學者從影響小波去噪的關鍵因素出發,如小波分解層數[4-5]、去噪評價[6-7]、閾值函數[8-10]、閾值估計[8-11]等對傳統小波閾值去噪進行研究。在閾值估計方面,常用的4種閾值分別為通用閾值、啟發式閾值、自適應閾值和極大極小閾值。其中,自適應閾值和極大極小閾值選取規則比較保守,僅將部分系數置0;而通用閾值和啟發式閾值選取規則的去除噪聲效果更佳,但易將有用信號的高頻部分當作噪聲信號除去[9]?？紤]到噪聲在小波變換下隨尺度增大而減小的特性,學者們給出了改進閾值估計算法,取得了一定效果,但也存在一些不足。如李明峰等[8]給出的閾值算法首層閾值估計過小,去噪不徹底;甘若等[10]、趙瑞珍等[11]給出的閾值算法首層閾值估計過大,易產生“過扼殺”小波系數的現象,且上述文獻都沒有考慮帶噪小波系數含噪幅度隨尺度增加而減小的函數關系。基于以上問題,本文提出一種新的小波閾值估計算法,并通過理論分析、工程實例驗證其優越性和有效性。

1 小波閾值去噪法

小波閾值去噪是基于小波變換的“集中”或“聚焦”能力。通常認為,信號在經過小波變換后,由信號產生的小波系數包含有信號的重要信息,其幅值大、數目少,而噪聲對應的小波系數幅值小、數目多。通過在不同分解尺度上選取合適的閾值,并將小于該閾值的小波系數置0,保留大于閾值的小波系數,使信號噪聲得到有效抑制,最后經過小波逆變換,得到去噪后的重構信號。

1.1 ZR閾值法

考慮到噪聲在小波變換下隨尺度增大而減小的特性[10],給出小波閾值法的改進法(稱為ZR閾值法)[11]:

(1)

1.2 1/2閾值法

對于均值為0、方差為σ2的高斯白噪聲n(t),其二進小波變換為WTj,kn(t),若小波函數ψj,k(t)為實函數時,則有[4]:

|WTj,kn(t)|2=|?Rn(t)ψj,k(t)dt|2

(2)

由于E[n(t)n(t)]=σ2,小波函數ψj,k(t)為小波基ψ(t)伸縮平移所得,小波函數中每一個都有規范化的范數‖ψj,k(t)‖=‖ψ‖=1,故:

(3)

即隨著小波分解尺度的增加,含有噪聲的觀測序列在低頻信號中噪聲影響幅度按1/2的規律逐漸減小[12]。隨著分解尺度的增加,高一級尺度上的低頻信號中含有的噪聲影響幅度會越來越小,曲線也會也來越光滑。

對變形觀測數據序列進行如式(4)所示的小波分解與重構:

an+1=an+dn=an-1+dn-1+dn=

…=a0+d1+d2+…+dn

(4)

式中,an為低頻信號,表現為趨勢項,dn為高頻信號,含有噪聲、粗差等局部變形信息。

設an噪聲影響幅度為σ2,由以上分析規律可知,an-1噪聲影響幅度為σ2/2,由式(4)可知,dn-1=an-an-1,即an噪聲影響幅度σ2在高一級尺度上一分為二,也就是an-1分解到的噪聲影響幅度為σ2/2,dn-1分解到的噪聲影響幅度為σ2/2。同理,dn-2=an-1-an-2,dn-2分解到的噪聲影響幅度為σ2/4。因此,隨著分解尺度的增加,高頻信號分解到的噪聲影響幅度也是按1/2的規律逐漸減小。

綜上所述,噪聲在小波變換下隨著分解尺度的增大而減小,且高頻信號分解到的噪聲影響幅度也是按1/2的規律逐漸減小。據此,本文給出新閾值計算公式(稱1/2閾值法):

(5)

2 仿真實驗

多期的變形監測信號一般包括趨勢項、周期項及多種噪聲,為保證與實際觀測數據相符合,構造復合信號進行仿真模擬實驗。如圖1所示,仿真信號時長為120 s,其中包括2個趨勢項、2個周期項及隨機噪聲項,實驗采用MATLAB中自帶的awgn函數,添加信噪比為db10的高斯白噪聲,其表達式為:

圖1 仿真數據Fig.1 Simulation data

Y=0.000 6t2-0.185 2t-

0.649 9sin(2πt/12+0.153 7)+

0.147 7sin(2πt/20+0.424 0)-

2.038 0+noise

為驗證本文1/2閾值法的降噪效果,采用分解能力較好的sym8小波基,使用文獻[4-5]的方法確定該信號最佳分解尺度為3層。使用硬閾值方法、軟閾值法、ZR閾值法進行對比,使用信噪比(SNR)、均方根誤差(RMSE)作為小波去噪的評價指標。

消噪的信號如圖2所示?？梢钥闯?軟閾值法去噪后曲線過于光滑,未能有效提取周期性信號,去噪效果最差;硬閾值法由于閾值函數的不連續性,在采樣序號75處有較大折線現象,即偽吉布斯現象,對周期性信號的提取有所提高,但在采樣序號5～10處未能提取到信息;ZR閾值法去噪信號曲線光滑,無多余折線現象,在采樣序號5～10處也未能提取到信息;1/2閾值法去噪信號曲線光滑程度最好,無多余振蕩,對周期性信號提取最好。從表1可知,1/2閾值法去噪的SNR最大,RMSE最小,去噪效果最佳,說明1/2閾值法去噪是可行的。

表1 4種閾值法的去噪指標Tab.1 Denoising indicators of four threshold methods

圖2 4種方法去噪結果Fig.2 Denoising results of four methods

3 工程實例

為進一步驗證本文小波閾值法去噪的有效性和優越性,對文獻[10]中某江岸危巖上H監測點10 a的沉降監測數據進行實例分析,如圖3所示?？梢钥闯?該危巖H點的高程隨時間的增加呈現出下降趨勢,同時含有周期性變形信息。由于各種不確定的因素,該信號含有很多毛刺,嚴重影響變形監測的準確分析和后續預測,因此需要對其進行去噪處理。

圖3 原始變形監測數據Fig.3 Original measuring deformation data

為了驗證1/2閾值法去噪的優越性,采用分解能力較好的sym8小波基,使用文獻[4-5]的方法確定該信號最佳分解尺度為3層,同時使用硬閾值法、軟閾值法、ZR閾值法進行對比,去噪結果如圖4所示。可以看出,軟閾值法由于估計小波系數和帶噪小波系數間有恒定的偏差,導致重構的信號失真比較嚴重,去噪后曲線過于光滑,未能有效保留周期性細節變形信息,去噪效果最差。硬閾值法由于閾值函數的不連續性,在采樣序號43處有較大折線現象,即偽吉布斯現象,去噪后曲線光滑性最差,對周期性細節變形信息保留有提高,如在采樣序號20～40處保留了周期性細節變形信息。ZR閾值法克服了軟、硬閾值法的缺點,也考慮到噪聲在小波變換下隨尺度增大而減小的特性,去噪后曲線的光滑性好,對周期性細節變形信息的保留程度有所提高,如在采樣序號50～60處保留了周期性細節變形信息。1/2閾值法克服了軟、硬閾值法的缺點,不僅與噪聲在小波變換各尺度上的傳播特性相一致,也更符合帶噪小波系數隨分解尺度增加而減小的函數關系,去噪后曲線的光滑性好,沒有多余振蕩,特別是在采樣序號110～120處,對周期性細節變形信息保留最好。從表2可知,1/2閾值法去噪的SNR最大,RMSE最小,與其他閾值法去噪對比,效果更好,對周期性細節變形信息提取更準確。

表2 4種閾值法的去噪指標Tab.2 Denoising indicators of four threshold methods

圖4 4種方法去噪結果Fig.4 Denoising results of four methods

4 結 語

通過理論分析與工程實例可知,本文1/2閾值法去噪效果優于軟、硬閾值法及ZR閾值法,能更有效地恢復被噪聲湮沒的周期性細節變形信息,進一步提高去噪的精度與可靠性,更好地反映危巖變形規律。1/2小波閾值法不僅與噪聲在小波變換各尺度上的傳播特性相一致,也更符合帶噪小波系數隨分解尺度的增加而減小的函數關系。需要指出的是,變形監測數據具有受噪聲干擾相對較小的特點,對噪聲干擾較大的復雜信號,本文方法的去噪效果有待進一步研究。