任意階PDE降噪特性分析*

尹愛軍, 張 泉,2, 戴宗賢

(1.重慶大學機械傳動國家重點實驗室 重慶,400044) (2.天津送變電工程公司 天津，300161)(3.重慶市計量質量檢測研究院第一分院 重慶，402260)

?

任意階PDE降噪特性分析*

尹愛軍1, 張 泉1,2, 戴宗賢3

(1.重慶大學機械傳動國家重點實驗室 重慶,400044) (2.天津送變電工程公司 天津，300161)(3.重慶市計量質量檢測研究院第一分院 重慶，402260)

通過對分數階微積分原理的研究，提出了任意階偏微分方程(partial differential equations, 簡稱PDE)降噪的統一模型，實現了基于任意階PDE降噪的數值化方法，并分析了任意階PDE降噪特性。該數值化方法能夠快速實現信號降噪，耗時少。通過仿真實驗，分析了PDE降噪性能的影響因素，與其他去噪方法進行了對比分析，并對現場實測信號進行了降噪分析。結果表明，PDE數值求解降噪方法性能優良，算法簡單。

偏微分方程； 分數階； 降噪； 振動信號

引 言

信號降噪在信號特征提取和參數化中具有重要的意義。濾波降噪理論與方法獲得了廣泛深入的研究,其理論從一般的平穩信號降噪方法發展到各種非平穩信號降噪方法，如小波去噪、經驗模態分解(empirical mode decomposition, 簡稱EMD)、奇異值分解降噪及滑動平均濾波等降噪方法[1]。小波變換是一種很好的非平穩信號降噪方法，得到了廣泛的研究和應用[2]。不同的小波基及閾值對去噪效果有很大影響[3]。EMD具有良好的自適應時頻分析能力，但EMD需要先驗知識，降低了該算法的魯棒性[4-5]。奇異值分解降噪是一種非線性濾波方法，奇異值數目選取是其降噪效果好壞的關鍵，但有效階次的選擇并沒有明確方法[6]。滑動平均濾波能夠有效抑制隨機誤差，提高信噪比，但高頻變化的確定性成分會因平均而被削弱[7]。偏微分方程(PDE)在圖像去噪、圖像增強等方面取得了很好的應用效果[8-9]。PDE在一維振動信號也得到一些應用，Baravdish等[10]通過矩陣奇異值分解與PDE實現對音頻信號進行降噪。Yin等[11]研究了基于PDE的信號修補方法，在EMD端點效應處理取得了很好的效果。文獻[12]將PDE用于振動信號去噪取得了很好的實際應用效果。和整數階PDE一樣，分數階PDE受到廣泛的關注。Pu等[13]通過分數階的PDE降噪模型對二維圖像信號進行降噪，取得了更好的效果。Bai等[14]利用非線性的各向異性分數階擴散方程獲得更自然的影像。然而，分數階 PDE的階次被限制在一個很小的范圍內，這影響了PDE降噪效果。

筆者根據分數階微積分原理，從傳統濾波器的幅頻特性角度出發，提出了任意階PDE降噪的統一模型，并研究了PDE濾波器的設計方法以及基于PDE降噪的數值化方法，建立了基于數值解降噪的一般過程。通過仿真實驗，分析了PDE降噪性能的影響因素，并與小波去噪等方法進行了對比分析。結果表明，PDE數值求解降噪方法性能優良，算法簡單。

1 任意階PDE濾波模型

1.1 分數階微積分定義

從不同的應用角度去分析問題可以得到不同的分數階微積分定義。至今為止，分數階微積分仍沒有統一的定義表達式[15]。目前，主要有4種經典的分數階微積分定義：Grünwald-Letnikov,Capotu,Riemann-Liouville和Cauchy定義。其中，后3種定義使用了Cauchy積分公式，計算復雜度較高，不利于分數階微分數值計算。

Grünwald-Letnikov定義是對整數階微分的差分近似遞推而來。式(1)為Grünwald-Letnikov分數階微積分定義表達式

(1)

Riemann-Liouville分數階的微積分定義是對Grünwald-Letnikov定義的改進，有

(2)

在一定條件下，可由Grünwald-Letnikov型分數階導數定義獲得Riemann-Liouville型分數階導數的數值化解法[16]。

Riemann-Liouville型分數階導數具有如下復合運算法則

(3)

1.2 任意階PDE濾波原理

筆者根據熱傳導方程，從濾波器角度推導了整數階偏微分方程的濾波原理和方法[11-12]。在這些研究基礎上，分析任意階偏微分方程的濾波特性。定義分數階偏微分方程式為

(4)

當f(x,t)=0時，利用傅里葉變換法，式(4)為

(5)

解得

(6)

其中:U(ω,t)=F(u(x,t)),Φ(ω)=F(φ(x))分別為u(x,t)，φ(x)的關于x的傅里葉變換。

K(ω,t)為

(7)

當φ(x)表示原始信號時，式(4)的解即為對信號φ(x)的濾波過程，K(ω,t)即為濾波器頻率響應。

因在復數域中，(iω)v可表示為

則式(7)可以表示為

(8)

若偏微分方程的階次v為奇數，幅頻特性為

(9a)

(9b)

進一步可得到任意階PDE降噪模型

(10)

定義信號φ(x)具有0邊界條件,即有u(k)(0,t)=0,(k=0,1,…,n-1)，且a=0,則式(3)可變為

式(10)變為

(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濾波器參數確定

設Y為濾波器截止頻率對應的衰減，由式(12)得到

即

(13)

其中:τ為演化時間。

由式(13)可知，當階次v及a確定后，通過確定濾波范圍即可確定演化時間。當v>0時，若截止頻率無窮大，τ趨近于0，此時輸入信號相當于通過了一個全通濾波器；而當截止頻率趨近于0，τ趨近于無窮大，PDE濾波后輸出直流;當v<0時，情形與v>0相反。

2 任意階PDE降噪的數值化

對于式(7),任意階PDE濾波器系數k(x,t)可以通過經典的濾波器設計方法得到, 則濾波后的信號為φ′(x)=φ(x)*k(x,t)。然而在濾波器數值化及截斷過程中，實際PDE濾波器與理想PDE濾波器之間存在差異。因此，通過對分數階PDE的數值化求解可更準確地實現PDE濾波。

當n=1時,階次1≤v<2, 式(11)變為

(14)

將u(x,t)變量分別替換成t和τ，將t和τ離散，并根據中心差分法、向后Euler格式和式(1)，式(14)變為

(15)

(16)

(17)

式(16)變為

即

(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)中，分數階微分方程向后Euler格式迭代過程對網格比q沒有限制[17]。對于式(19)，當演化步長l=τ，總的迭代次數s=τ/l=1，即形成了PDE降噪的快速算法。

(21)

輸入噪聲信號u0,截止頻率為fn,采樣頻率為fs,PDE階次為v。輸出降噪后信號u。任意分數階PDE降噪算法如下：

2) 根據PDE的階次v、采樣頻率fs和截止頻率fn計算網格比q;

3) 根據gb和網格比q計算集合C;

4) 計算矩陣A并得到濾波矩陣A-1;

5) 得到降噪后信號u。

3 實驗分析與應用

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中可以看出，隨著階次的增大，信噪比并未呈單調增加，在4階附近獲得較好的降噪效果。這是因為在濾波器設計過程中，應根據噪聲信號的類型不同，合理設置濾波器通帶截止頻率和阻帶截止頻率。對于本研究的仿真信號，由于噪聲為寬帶噪聲，當階次較高時，雖PDE濾波器符合優良濾波特性要求，但低頻噪聲將會更多地保留，從而降低了濾波效果。所以信噪比并未呈單調增加，而是在4階附近獲得較好的降噪效果。當為奇次階時，PDE濾波器幅頻特性為1，濾波前后信號幅頻特性不變，因此信噪比不變，且奇次階附近階次降噪效果較差。

圖3 降噪效果(SNR)隨階次變化圖Fig.3 Noise reductions (SNR) with different order

為了對比實際濾波器跟理想濾波器之間的差異，定義相對誤差百分比A為

(22)

其中：B為實際值；C為理論值。

圖4為IIR數字濾波器和PDE濾波器的相對誤差百分比對比圖。其中，IIR數字濾波器和PDE濾波器的參數相同，通帶截止頻率為35Hz，阻帶截止頻率為67Hz，PDE濾波器的階次為4。從圖4可以看出，IIR數字濾波器的誤差比PDE濾波器要高且波動大,PDE實際濾波器與理論濾波器之間的差異很小。

圖4 IIR數字濾波器與PDE濾波器的相對誤差對比圖Fig.4 The percentage absolute relative error of IIR filter and the proposed PDE filter

3.2 與其他降噪方法的比較

采用小波去噪、IIR數字濾波以及4階PDE差分求解去噪等方法分別對不同信噪比的仿真信號進行去噪分析。表1為不同濾波方法下降噪效果的對比。

表1 不同濾波方法下的降噪效果

表1中，小波去噪的小波基為db8，分解層數為5，閾值函數為軟閾值，表中IIR濾波采用的是3階巴特沃斯數字濾波器。根據圖3，選擇4階PDE進行比較。其中，IIR數字濾波器法、PDE差分求法的截止頻率為35 Hz。從表1看出，在信噪比較大時，小波方法和PDE差分方法具有相近的降噪效果(>5 dB)；而信噪比較低時，PDE差分方法有較明顯的降噪效果(-5 dB)。在實際應用中，小波去噪過程復雜，需要根據信號的特點選擇不同的小波基、閾值函數等。圖5為對信號比SNR=5 dB的信號降噪比較。另外，從圖5(c)中可以看出，IIR數字濾波會產生時移(紅色標記)。

圖5 不同方法的降噪效果圖Fig.5 Noise reduction of different filtering methods

表2為幾種去噪方法的時間性能比較。 CPU為Intel(R) Core(TM) i3，主頻為2.27GHz，內存為4GB，每個方法運行8次，剔除最高和最低值，求剩余6次平均運行時間。

表2 時間性能比較

從表2可以看出，幾種去噪方法中IIR 濾波、PDE數值求解耗時相當，小波降噪方法耗時最長。

3.3 現場應用

軸心軌跡作為轉子振動信號的一類重要圖形征兆，其形狀和動態特性包含了豐富的故障征兆信息。然而在現場實際測量的振動信號中往往都會受到噪聲污染，軸心軌跡非常復雜，不易獲得清晰的特征。因此，對軸心軌跡降噪成為轉子特征提取的關鍵之一。圖6是在重慶大學測試中心的轉子實驗臺現場實驗圖。

圖6 現場實驗圖Fig.6 Field experiment

圖7(a)和7(b)分別為轉軸原始軸心軌跡圖和使用PDE降噪之后的軸心軌跡圖。其轉子轉速為6 970 r/min，采樣頻率為1 652 Hz。從圖7(a)中可以看出，原始軸心軌跡較混亂。如圖3所示，PDE濾波器在階次4附近區域的降噪效果相近，為了方便PDE濾波的數值計算，圖7(b)PDE降噪階次選取整數階次4，PDE降噪的截止頻率為85 Hz。由圖7可以看出，該軸的軸心軌跡為橢圓形，表明該軸存在嚴重的不平衡，而其他故障很小。

圖7 濾波前后的軸心軌跡圖Fig.7 Original and PDE filtered shaft centerline orbits

4 結束語

PDE在圖像濾波、修補和增強等方面都取得了很好的效果。將PDE用于一維信號的降噪處理的研究才剛剛起步，主要集中在整數階的降噪方法。筆者從濾波器設計的角度，提出了任意階PDE濾波器統一模型，研究了PDE濾波的兩種數值化過程，建立了基于數字差分求解的降噪一般過程，分析了階次對降噪性能和時間性能的影響。通過仿真信號對小波閾值去噪，IIR數字濾波降噪及PDE數值求解降噪等多種去噪方法進行了對比分析。實驗表明，分數階PDE數值求解去噪方法算法簡單，效果良好，是一種有效的信號去噪方法。

[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] 秦毅,王家序,毛永芳.基于軟閾值和小波模極大值重構的信號降噪[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].制造業自動化, 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] 楊柱中,周激流,晏祥玉,等. 基于分數階微分的圖像增強[J]. 計算機輔助設計與圖形學學報,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]蒲亦非. 分數階微積分在現代信號分析與處理中應用的研究[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); 中央高校基本科研業務費資助項目(CDJZR13 11 55 01)

2014-12-09；

2015-01-07

TN911

尹愛軍，男，1978年5月出生。教授、博士生導師。主要研究方向為智能測試與虛擬儀器、現代信號分析處理及故障檢測與診斷等。曾發表《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.cn