尹愛軍, 張 泉,2, 戴宗賢
(1.重慶大學機械傳動國家重點實驗室 重慶,400044) (2.天津送變電工程公司 天津,300161)(3.重慶市計量質(zhì)量檢測研究院第一分院 重慶,402260)
?
任意階PDE降噪特性分析*
尹愛軍1, 張 泉1,2, 戴宗賢3
(1.重慶大學機械傳動國家重點實驗室 重慶,400044) (2.天津送變電工程公司 天津,300161)(3.重慶市計量質(zhì)量檢測研究院第一分院 重慶,402260)
通過對分數(shù)階微積分原理的研究,提出了任意階偏微分方程(partial differential equations, 簡稱PDE)降噪的統(tǒng)一模型,實現(xiàn)了基于任意階PDE降噪的數(shù)值化方法,并分析了任意階PDE降噪特性。該數(shù)值化方法能夠快速實現(xiàn)信號降噪,耗時少。通過仿真實驗,分析了PDE降噪性能的影響因素,與其他去噪方法進行了對比分析,并對現(xiàn)場實測信號進行了降噪分析。結(jié)果表明,PDE數(shù)值求解降噪方法性能優(yōu)良,算法簡單。
偏微分方程; 分數(shù)階; 降噪; 振動信號
信號降噪在信號特征提取和參數(shù)化中具有重要的意義。濾波降噪理論與方法獲得了廣泛深入的研究,其理論從一般的平穩(wěn)信號降噪方法發(fā)展到各種非平穩(wěn)信號降噪方法,如小波去噪、經(jīng)驗模態(tài)分解(empirical mode decomposition, 簡稱EMD)、奇異值分解降噪及滑動平均濾波等降噪方法[1]。小波變換是一種很好的非平穩(wěn)信號降噪方法,得到了廣泛的研究和應用[2]。不同的小波基及閾值對去噪效果有很大影響[3]。EMD具有良好的自適應時頻分析能力,但EMD需要先驗知識,降低了該算法的魯棒性[4-5]。奇異值分解降噪是一種非線性濾波方法,奇異值數(shù)目選取是其降噪效果好壞的關(guān)鍵,但有效階次的選擇并沒有明確方法[6]?;瑒悠骄鶠V波能夠有效抑制隨機誤差,提高信噪比,但高頻變化的確定性成分會因平均而被削弱[7]。偏微分方程(PDE)在圖像去噪、圖像增強等方面取得了很好的應用效果[8-9]。PDE在一維振動信號也得到一些應用,Baravdish等[10]通過矩陣奇異值分解與PDE實現(xiàn)對音頻信號進行降噪。Yin等[11]研究了基于PDE的信號修補方法,在EMD端點效應處理取得了很好的效果。文獻[12]將PDE用于振動信號去噪取得了很好的實際應用效果。和整數(shù)階PDE一樣,分數(shù)階PDE受到廣泛的關(guān)注。Pu等[13]通過分數(shù)階的PDE降噪模型對二維圖像信號進行降噪,取得了更好的效果。Bai等[14]利用非線性的各向異性分數(shù)階擴散方程獲得更自然的影像。然而,分數(shù)階 PDE的階次被限制在一個很小的范圍內(nèi),這影響了PDE降噪效果。
筆者根據(jù)分數(shù)階微積分原理,從傳統(tǒng)濾波器的幅頻特性角度出發(fā),提出了任意階PDE降噪的統(tǒng)一模型,并研究了PDE濾波器的設(shè)計方法以及基于PDE降噪的數(shù)值化方法,建立了基于數(shù)值解降噪的一般過程。通過仿真實驗,分析了PDE降噪性能的影響因素,并與小波去噪等方法進行了對比分析。結(jié)果表明,PDE數(shù)值求解降噪方法性能優(yōu)良,算法簡單。
1.1 分數(shù)階微積分定義
從不同的應用角度去分析問題可以得到不同的分數(shù)階微積分定義。至今為止,分數(shù)階微積分仍沒有統(tǒng)一的定義表達式[15]。目前,主要有4種經(jīng)典的分數(shù)階微積分定義:Grünwald-Letnikov,Capotu,Riemann-Liouville和Cauchy定義。其中,后3種定義使用了Cauchy積分公式,計算復雜度較高,不利于分數(shù)階微分數(shù)值計算。
Grünwald-Letnikov定義是對整數(shù)階微分的差分近似遞推而來。式(1)為Grünwald-Letnikov分數(shù)階微積分定義表達式
(1)
Riemann-Liouville分數(shù)階的微積分定義是對Grünwald-Letnikov定義的改進,有
(2)
在一定條件下,可由Grünwald-Letnikov型分數(shù)階導數(shù)定義獲得Riemann-Liouville型分數(shù)階導數(shù)的數(shù)值化解法[16]。
Riemann-Liouville型分數(shù)階導數(shù)具有如下復合運算法則
(3)
1.2 任意階PDE濾波原理
筆者根據(jù)熱傳導方程,從濾波器角度推導了整數(shù)階偏微分方程的濾波原理和方法[11-12]。在這些研究基礎(chǔ)上,分析任意階偏微分方程的濾波特性。定義分數(shù)階偏微分方程式為
(4)
當f(x,t)=0時,利用傅里葉變換法,式(4)為
(5)
解得
(6)
其中:U(ω,t)=F(u(x,t)),Φ(ω)=F(φ(x))分別為u(x,t),φ(x)的關(guān)于x的傅里葉變換。
K(ω,t)為
(7)
當φ(x)表示原始信號時,式(4)的解即為對信號φ(x)的濾波過程,K(ω,t)即為濾波器頻率響應。
因在復數(shù)域中,(iω)v可表示為
則式(7)可以表示為
(8)
若偏微分方程的階次v為奇數(shù),幅頻特性為
(9a)
(9b)
進一步可得到任意階PDE降噪模型
(10)
定義信號φ(x)具有0邊界條件,即有u(k)(0,t)=0,(k=0,1,…,n-1),且a=0,則式(3)可變?yōu)?/p>
式(10)變?yōu)?/p>
(11)
其中:x∈R;t>0;n∈Z;v=r+n;0 此時 (12) 若v>0時,K(ω,t)為低通濾波器;若v<0時,為高通濾波過程。階次越高,濾波器衰減的越快。不同階次對應的濾波器幅頻特性如圖1所示。 圖1 濾波器幅頻特性Fig.1 Filter magnitude-frequency curves 1.3 PDE濾波器參數(shù)確定 設(shè)Y為濾波器截止頻率對應的衰減,由式(12)得到 即 (13) 其中:τ為演化時間。 由式(13)可知,當階次v及a確定后,通過確定濾波范圍即可確定演化時間。當v>0時,若截止頻率無窮大,τ趨近于0,此時輸入信號相當于通過了一個全通濾波器;而當截止頻率趨近于0,τ趨近于無窮大,PDE濾波后輸出直流;當v<0時,情形與v>0相反。 對于式(7),任意階PDE濾波器系數(shù)k(x,t)可以通過經(jīng)典的濾波器設(shè)計方法得到, 則濾波后的信號為φ′(x)=φ(x)*k(x,t)。然而在濾波器數(shù)值化及截斷過程中,實際PDE濾波器與理想PDE濾波器之間存在差異。因此,通過對分數(shù)階PDE的數(shù)值化求解可更準確地實現(xiàn)PDE濾波。 當n=1時,階次1≤v<2, 式(11)變?yōu)?/p> (14) 將u(x,t)變量分別替換成t和τ,將t和τ離散,并根據(jù)中心差分法、向后Euler格式和式(1),式(14)變?yōu)?/p> (15) (16) (17) 式(16)變?yōu)?/p> 即 (18a) 記C={c1,c2,…,cm,cm+1,cm+2},其中c1=-qg0,c2=-qg1+1。當3≤j≤m+1,cj=-q(gj-1-gj-3),cm+2=qgm-1,則 (18b) (19a) (19b) 式(18)可以記為 (20) 式(15)中,分數(shù)階微分方程向后Euler格式迭代過程對網(wǎng)格比q沒有限制[17]。對于式(19),當演化步長l=τ,總的迭代次數(shù)s=τ/l=1,即形成了PDE降噪的快速算法。 (21) 輸入噪聲信號u0,截止頻率為fn,采樣頻率為fs,PDE階次為v。輸出降噪后信號u。任意分數(shù)階PDE降噪算法如下: 2) 根據(jù)PDE的階次v、采樣頻率fs和截止頻率fn計算網(wǎng)格比q; 3) 根據(jù)gb和網(wǎng)格比q計算集合C; 4) 計算矩陣A并得到濾波矩陣A-1; 5) 得到降噪后信號u。 3.1 不同階次降噪性能比較 圖2(a)為頻率分別為5,7,9,10和20 Hz的正弦信號疊加而成的仿真信號,各信號分量的幅值為1 V。對該信號分別疊加不同信噪比的高斯白噪聲(信噪比分別為-10,-5,0,5和10 dB),采樣頻率為1 kHz,取1 000個采樣點進行降噪對比實驗。其中,圖2(b)為信噪比為-10 dB的高斯白噪聲。 圖2 仿真信號Fig.2 The simulation signal 圖3為降噪效果隨階次變化曲線。從圖3中可以看出,隨著階次的增大,信噪比并未呈單調(diào)增加,在4階附近獲得較好的降噪效果。這是因為在濾波器設(shè)計過程中,應根據(jù)噪聲信號的類型不同,合理設(shè)置濾波器通帶截止頻率和阻帶截止頻率。對于本研究的仿真信號,由于噪聲為寬帶噪聲,當階次較高時,雖PDE濾波器符合優(yōu)良濾波特性要求,但低頻噪聲將會更多地保留,從而降低了濾波效果。所以信噪比并未呈單調(diào)增加,而是在4階附近獲得較好的降噪效果。當為奇次階時,PDE濾波器幅頻特性為1,濾波前后信號幅頻特性不變,因此信噪比不變,且奇次階附近階次降噪效果較差。 圖3 降噪效果(SNR)隨階次變化圖Fig.3 Noise reductions (SNR) with different order 為了對比實際濾波器跟理想濾波器之間的差異,定義相對誤差百分比A為 (22) 其中:B為實際值;C為理論值。 圖4為IIR數(shù)字濾波器和PDE濾波器的相對誤差百分比對比圖。其中,IIR數(shù)字濾波器和PDE濾波器的參數(shù)相同,通帶截止頻率為35Hz,阻帶截止頻率為67Hz,PDE濾波器的階次為4。從圖4可以看出,IIR數(shù)字濾波器的誤差比PDE濾波器要高且波動大,PDE實際濾波器與理論濾波器之間的差異很小。 圖4 IIR數(shù)字濾波器與PDE濾波器的相對誤差對比圖Fig.4 The percentage absolute relative error of IIR filter and the proposed PDE filter 3.2 與其他降噪方法的比較 采用小波去噪、IIR數(shù)字濾波以及4階PDE差分求解去噪等方法分別對不同信噪比的仿真信號進行去噪分析。表1為不同濾波方法下降噪效果的對比。 表1 不同濾波方法下的降噪效果 表1中,小波去噪的小波基為db8,分解層數(shù)為5,閾值函數(shù)為軟閾值,表中IIR濾波采用的是3階巴特沃斯數(shù)字濾波器。根據(jù)圖3,選擇4階PDE進行比較。其中,IIR數(shù)字濾波器法、PDE差分求法的截止頻率為35 Hz。從表1看出,在信噪比較大時,小波方法和PDE差分方法具有相近的降噪效果(>5 dB);而信噪比較低時,PDE差分方法有較明顯的降噪效果(-5 dB)。在實際應用中,小波去噪過程復雜,需要根據(jù)信號的特點選擇不同的小波基、閾值函數(shù)等。圖5為對信號比SNR=5 dB的信號降噪比較。另外,從圖5(c)中可以看出,IIR數(shù)字濾波會產(chǎn)生時移(紅色標記)。 圖5 不同方法的降噪效果圖Fig.5 Noise reduction of different filtering methods 表2為幾種去噪方法的時間性能比較。 CPU為Intel(R) Core(TM) i3,主頻為2.27GHz,內(nèi)存為4GB,每個方法運行8次,剔除最高和最低值,求剩余6次平均運行時間。 表2 時間性能比較 從表2可以看出,幾種去噪方法中IIR 濾波、PDE數(shù)值求解耗時相當,小波降噪方法耗時最長。 3.3 現(xiàn)場應用 軸心軌跡作為轉(zhuǎn)子振動信號的一類重要圖形征兆,其形狀和動態(tài)特性包含了豐富的故障征兆信息。然而在現(xiàn)場實際測量的振動信號中往往都會受到噪聲污染,軸心軌跡非常復雜,不易獲得清晰的特征。因此,對軸心軌跡降噪成為轉(zhuǎn)子特征提取的關(guān)鍵之一。圖6是在重慶大學測試中心的轉(zhuǎn)子實驗臺現(xiàn)場實驗圖。 圖6 現(xiàn)場實驗圖Fig.6 Field experiment 圖7(a)和7(b)分別為轉(zhuǎn)軸原始軸心軌跡圖和使用PDE降噪之后的軸心軌跡圖。其轉(zhuǎn)子轉(zhuǎn)速為6 970 r/min,采樣頻率為1 652 Hz。從圖7(a)中可以看出,原始軸心軌跡較混亂。如圖3所示,PDE濾波器在階次4附近區(qū)域的降噪效果相近,為了方便PDE濾波的數(shù)值計算,圖7(b)PDE降噪階次選取整數(shù)階次4,PDE降噪的截止頻率為85 Hz。由圖7可以看出,該軸的軸心軌跡為橢圓形,表明該軸存在嚴重的不平衡,而其他故障很小。 圖7 濾波前后的軸心軌跡圖Fig.7 Original and PDE filtered shaft centerline orbits PDE在圖像濾波、修補和增強等方面都取得了很好的效果。將PDE用于一維信號的降噪處理的研究才剛剛起步,主要集中在整數(shù)階的降噪方法。筆者從濾波器設(shè)計的角度,提出了任意階PDE濾波器統(tǒng)一模型,研究了PDE濾波的兩種數(shù)值化過程,建立了基于數(shù)字差分求解的降噪一般過程,分析了階次對降噪性能和時間性能的影響。通過仿真信號對小波閾值去噪,IIR數(shù)字濾波降噪及PDE數(shù)值求解降噪等多種去噪方法進行了對比分析。實驗表明,分數(shù)階PDE數(shù)值求解去噪方法算法簡單,效果良好,是一種有效的信號去噪方法。 [1] 錢征文,程禮,李應紅.利用奇異值分解的信號降噪方法[J]. 振動、測試與診斷,2011(4):459-463. Qian Zhengwen, Cheng Li, Li Yinghong. Noise reduction method based on singular value decomposition[J]. Journal of Vibration, Measurement & Diagnosis, 2011(4):459-463. (in Chinese) [2] Donoho D L. De-noising by soft-thresholding information theory[J]. IEEE Transactions on , 1995,41(3):613-627. [3] 秦毅,王家序,毛永芳.基于軟閾值和小波模極大值重構(gòu)的信號降噪[J]. 振動、測試與診斷,2011(5):543-547. Qin Yi, Wang Jiaxu, Mao Yongfang. Signal denoising based on soft thresholding and reconstruction from dyadic wavelet transform modulus maxima[J]. Journal of Vibration, Measurement & Diagnosis, 2011(5):543-547. (in Chinese) [4] Dragomiretskiy K, Zosso D. Variational mode decomposition signal processing[J]. IEEE Transactions on, 2014,62(3):531-544. [5] 徐仁林,安偉.小波降噪在信號基于EMD的Hilbert變換中的應用[J]. 噪聲與振動控制,2008(3):74-77. Xu Renlin, An Wei. Wavelet denoise application in the signal hilbert transform based on EMD[J]. Noise and Vibration Control, 2008(3):74-77. (in Chinese) [6] 張磊,彭偉才,原春暉,等.奇異值分解降噪的改進方法[J]. 中國艦船研究,2012(5):83-88. Zhang Lei, Peng Weicai, Yuan Chunhui, et al. An improved method for noise reduction based on singular value decomposition[J]. Chinese Journal of Ship Research, 2012(5):83-88. (in Chinese) [7] 劉藝柱,楊瑞蘭.采用滑動平均濾波法提高硬幣識別準確率的研究[J].制造業(yè)自動化, 2010(1):42-44. Liu Yizhu, Yang Ruilan. Moving average filtering method used to improve the accuracy of coins to identify research[J]. Manufacturing Automation, 2010(1):42-44. (in Chinese) [8] 楊柱中,周激流,晏祥玉,等. 基于分數(shù)階微分的圖像增強[J]. 計算機輔助設(shè)計與圖形學學報,2008(3):343-348. Yang Zhuzhong, Zhou Jiliu, Yan Xiangyu, et al. Image enhancement based on fractional differentials[J]. Journal of Computer-Aided Design & Computer Graphics, 2008(3):343-348. (in Chinese) [9] Machado B B, Goncalves W N, Bruno O M. Enhancing the texture attribute with partial differential equations: a case of study with Gabor filters[J]. Advanced Concepts for Intelligent Vision Systems, 2011,6915:337-348. [10]Baravdish G, Evangelista G, Svensson O, et al. PDE-SVD based audio denoising[C]∥Communications Control and Signal Processing (ISCCSP). Washington DC, USA: IEEE, 2012:1-6. [11]Yin Aijun, Zhao Lei, Yang Zhengyi, et al. Noise reduction method for vibration signals 2D time‐frequency distribution using anisotropic diffusion equation [J]. Mathematical Methods in the Applied Sciences, 2015,38(4):609-616. [12]尹愛軍,王璇. PDE信號修補方法及在EMD端點效應處理中的應用[J]. 振動與沖擊,2012(2):6-9,61. Yin Aijun, Wang Xuan. Signal restoring based on PDE and its application in end effect processing of EMD[J]. Journal of Vibration and Shock, 2012(2):6-9,61. (in Chinese) [13]Pu Yifei, Siarry P, Zhou Jiliu, et al. Fractional partial differential equation denoising models for texture image[J]. Science China-Information Sciences, 2014,57(7):1-19. [14]Bai J, Feng X C. Fractional-order anisotropic diffusion for image denoising[J]. Ieee Transactions on Image Processing, 2007,16(10):2492-2502. [15]蒲亦非. 分數(shù)階微積分在現(xiàn)代信號分析與處理中應用的研究[D].成都:四川大學,2006. [16]Chen Changming, Liu F, Turner I, et al. A Fourier method for the fractional diffusion equation describing sub-diffusio[J]. Journal of Computational Physics, 2007,227(2):886-897. [17]El-Mikkawy M E A. On the inverse of a gene tridiagonal matrix[J]. Applied Mathematics and Computation, 2004,150(3):669-679. 10.16450/j.cnki.issn.1004-6801.2016.06.006 *國家自然科學基金資助項目(51105396); 中央高?;究蒲袠I(yè)務費資助項目(CDJZR13 11 55 01) 2014-12-09; 2015-01-07 TN911 尹愛軍,男,1978年5月出生。教授、博士生導師。主要研究方向為智能測試與虛擬儀器、現(xiàn)代信號分析處理及故障檢測與診斷等。曾發(fā)表《Thermography spatial-transient-stage mathematical tensor construction and material property variation track》(《International Journal of Thermal Science》2014,Vol.85)等論文。 E-mail:aijun.yin@cqu.edu.cn2 任意階PDE降噪的數(shù)值化
3 實驗分析與應用
4 結(jié)束語