格子玻爾茲曼模型的圖像去噪方法

陳榮前，聶德明

(中國計量大學(xué) 計量測試工程學(xué)院，浙江 杭州 310018)

?

格子玻爾茲曼模型的圖像去噪方法

陳榮前，聶德明

(中國計量大學(xué) 計量測試工程學(xué)院，浙江 杭州 310018)

為了實現(xiàn)快速去噪以及更好地保護(hù)圖像的邊緣信息，從各向同性擴(kuò)散模型出發(fā)，提出了一種基于索貝爾(Sobel)梯度算子的格子玻爾茲曼去噪模型，通過該模型求解非線性擴(kuò)散方程實現(xiàn)圖像去噪.模型引入了Sobel算子作為檢測邊緣的工具，同時滿足前向擴(kuò)散的要求，保證了方程的穩(wěn)定性.實驗分析表明，應(yīng)用該方法進(jìn)行圖像去噪，能夠在保護(hù)圖像邊緣的情況下進(jìn)行快速去噪.

圖像去噪；擴(kuò)散方程；格子玻爾茲曼模型；索貝爾梯度算子

在圖像的獲取和傳輸過程中，由于圖像設(shè)備性能受到各種因素的影響，使得圖像不可避免會受到噪聲的污染.圖像的濾波去噪研究一直是圖像處理中的熱點問題[1].基于擴(kuò)散方程的圖像濾波方法是從20世紀(jì)90年代發(fā)展起來，并且表現(xiàn)出優(yōu)于傳統(tǒng)圖像濾波方法的特性[2].Koenderink[3]等較早地進(jìn)行了熱傳導(dǎo)方程的圖像去噪研究.1990年P(guān)erona和Malik[4]提出了各向異性的擴(kuò)散模型(PM模型)，其基本思路是對于邊緣和目標(biāo)區(qū)域采用不同的去噪策略，從而可以在去噪的同時保護(hù)圖像的邊緣信息.但是該模型存在逆向擴(kuò)散問題，從而導(dǎo)致該模型存在穩(wěn)定性差的缺點[5].Catté[6]等對模型進(jìn)行了改進(jìn)，解決了PM模型的逆向擴(kuò)散問題.然而，一般情況下非線性擴(kuò)散方程難以求得解析解，隨著計算機(jī)技術(shù)的發(fā)展，數(shù)值方法成為求解擴(kuò)散方程的有效手段.如有限元、有限差分等方法通過對原有的微分方程進(jìn)行離散，得到差分格式，進(jìn)而利用初始條件和邊界條件進(jìn)行迭代求解；雖然精度達(dá)到了求解問題的要求，但往往計算效率較低，特別是對于像素量在幾十萬甚至上百萬級別的圖片而言，傳統(tǒng)的數(shù)值求解方法無法滿足應(yīng)用要求.因此，選擇合適的數(shù)值方法來求解這些方程顯得至關(guān)重要.

格子玻爾茲曼方法(Lattice Boltzmann method, LBM)是一種用來模擬流體宏觀流動的方法.可以用于求解N-S方程，由于其物理意義清晰，易于編程實現(xiàn)，且具有良好的計算并行性等特性，因此該算法受到了廣泛的關(guān)注與應(yīng)用[7].近年來LBM已經(jīng)成功應(yīng)用于求解擴(kuò)散方程，Chang[8]等根據(jù)變分模型，提出了LBM去噪算法，但是其使用單一的松弛因子，實現(xiàn)的是各項同性的擴(kuò)散模型，該模型忽略了圖像的梯度信息，從而在多次迭代后容易造成圖像的模糊.本文在各向同性擴(kuò)散模型的基礎(chǔ)上，提出了一種基于Sobel梯度算子的格子波爾茲曼模型，模型將Sobel梯度算子作為邊緣截止函數(shù)，將其引入到格子波爾茲曼模型的松弛因子中，從而能夠在去噪的同時可以很好地保留圖像的邊緣信息，且算法簡單，去噪效果好，具有實際應(yīng)用價值.

1 格子玻爾茲曼模型

典型的D2Q9格子玻爾茲曼方法的格子結(jié)構(gòu)如圖1，每個格點具有9個方向的格子速度ci和速度分布函數(shù)fi(x,t)，其中x表示格子的位置，t表示時間，i表示速度方向.

圖1 D2Q9格子玻爾茲曼網(wǎng)格結(jié)構(gòu)圖Figure 1 Structure of D2Q9 LBM

宏觀流體密度和速度u的計算表達(dá)式為

(1)

格子玻爾茲曼碰撞流過程的演化方程為

fi(x+ciΔt,t+Δt)-fi(x,t)=

(2)

式(2)中的局部平衡態(tài)分布函數(shù)f(eq)由Bhatnager等人[9]給出

Ci(ci·u)2+Diu2).

(3)

其中Ai至Di是與網(wǎng)格結(jié)構(gòu)相對應(yīng)的常量.表示無量綱松弛時間.

假設(shè)流動的馬赫數(shù)充分小，根據(jù)Chapman-Enskog展開，可以導(dǎo)出不可壓N-S方程組

·u=0，

(4)

上式中，σ為應(yīng)力張量.當(dāng)網(wǎng)格步長的平方與時間步長比值為1時，即dx×dx/dt=1，松弛時間τ由流體的運動粘度ν決定

(5)

為了將LBM應(yīng)用于求解擴(kuò)散方程，將式(3)修改為

(6)

將式(6)代入式(1)，并進(jìn)行Chapman-Enskog分析，Chapman-Enskog是一種多尺度技術(shù)：

(7)

ε是努森數(shù)，表示分子平均自由程與宏觀長度的比值.將式(2)中fi(x+ciΔt,t+Δt)做泰勒張開，對fi(x,t)按照Chapman-Enskog展開，然后對比方程兩邊相同階數(shù)的結(jié)果：

(8)

(9)

(10)

fi滿足下面關(guān)系式：

(11)

(12)

對(10)式子所有i方向進(jìn)行累加得到

(13)

聯(lián)立式(6)、(9)、(12)、(13)得到拋物型擴(kuò)散方程如下

(14)

其中Dxy表示擴(kuò)散張量

(15)

若τ是常數(shù)，則可以得到拋物型擴(kuò)散方程

(16)

(17)

由式(17)可知，若格子玻爾茲曼模型中的松弛時間為常數(shù)，得到的擴(kuò)散方程的中的擴(kuò)散系數(shù)也為常數(shù)，這意味著模型的平滑效果在圖像的各個點都會相同，從而在去噪聲的同時也會損失邊緣的信息.所以需要重新設(shè)計松弛因子.

2 松弛因子設(shè)計

2.1 Sobel梯度算子

為了達(dá)到保護(hù)邊緣的目的，必須在邊緣處抑制平滑，在目標(biāo)區(qū)域加速平滑，即在邊緣附近需要用較小擴(kuò)散系數(shù)，而在目標(biāo)區(qū)域采用較大的擴(kuò)散系數(shù).圖像的梯度圖像經(jīng)常用于檢測圖像中的邊緣信息，其中，Sobel算子是一種離散型的差分算子，如圖2，該算子包含兩個坐標(biāo)軸方向的兩組3×3的矩陣wx,wy.

圖2 Sobel模板

(18)

(19)

其中,I(x+s,y+t)代表圖像的灰度在(x+s,y+t)處的灰度值，為了加快計算速度，通過絕對值近似得到圖像的梯度幅值

MSobel(x,y)≈|gx|+|gy|.

(20)

Sobel梯度算子檢測邊緣的效果如圖3，圖3(a)為原始圖像，圖3(b)為用Sobel梯度幅值圖像，可以看到Sobel算子能夠檢測到較為精確的圖像邊緣信息.而且Sobel算子對噪聲的敏感性較低[11]，因此適合用來設(shè)計松弛因子.

圖3 Sobel算子邊緣檢測效果Figure 3 Results of Sobel operater

2.2 基于Sobel算子的松弛因子

由2.1節(jié)的分析可知，Sobel梯度幅值越大的位置越有可能是邊緣位置，因此在這個位置的擴(kuò)散系數(shù)需要適當(dāng)減小，而在幅值越小的位置越有可能是平坦區(qū)域需要適當(dāng)增加擴(kuò)散系數(shù).另外，從(17)式可以看出，松弛因子越大，擴(kuò)散系數(shù)越大，因此設(shè)計的松弛因子函數(shù)應(yīng)該隨著Sobel梯度幅值越大而減小，且為了保證方程的穩(wěn)定性，擴(kuò)散系數(shù)必須大于0,所以松弛因子τ≥0.5.將文獻(xiàn)[12]中用邊緣檢測函數(shù)替換為Sobel梯度幅值函數(shù)，得到的松弛因子τs為

(21)

上式中L為一個常數(shù)，將式(21)代入式(18)得到擴(kuò)散系數(shù)與梯度幅值的關(guān)系

(22)

從式(22)可以看出圖像中Sobel梯度算子越大的位置，擴(kuò)散系數(shù)越小，因此從理論上可以達(dá)到保護(hù)邊緣的效果;而在目標(biāo)區(qū)域，即圖像灰度幅值較小的位置，擴(kuò)散系數(shù)越大，從而由加速去噪的效果.而且上式滿足擴(kuò)散系數(shù)d≥0的條件，從而保證了方程的穩(wěn)定性.

3 實驗分析

為了驗證基于Sobel梯度算子的格子玻爾茲曼方法(SLBM)進(jìn)行圖像去噪的有效性，對圖像進(jìn)行去噪實驗.首先需要對圖像添加噪聲，圖像的噪聲一般有沖擊噪聲、加性噪聲等.本文選取高斯噪聲作為噪聲源進(jìn)行去噪仿真.

在圖像的獲取及傳輸階段，常常會由于低照明度或傳感器的溫度變化以及電子電路干擾產(chǎn)生高斯噪聲.均值為μ,方差為σ高斯噪聲的概率密度函數(shù)G(x,y)為

(23)

在測試圖像中添加μ=0，σ2=0.01.的高斯噪聲，圖4(c)和圖4(d)分別為采用各向同性擴(kuò)散模型和SLBM進(jìn)行圖像去噪的實驗結(jié)果，從圖中可以看到，各向同性擴(kuò)散模型，在去噪的同時會噪聲圖像邊緣信息的損失.而在本文的算法中，由于擴(kuò)散系數(shù)灰度變化相關(guān)，因此屬于各向異性擴(kuò)散模型，與在去噪的的同時能夠很好的保留圖像中的邊緣信息，更加逼近于原圖像.

圖4 SLBM和各向同性擴(kuò)散模型去噪效果Figure 4 Denoise results of SLBM and isotropic diffusion model

為了定量評估本文算法的去噪效果，計算圖像的峰值性噪比隨迭代次數(shù)的關(guān)系如圖5，其中峰值性噪比(PNSR)的計算公式為

(24)

式(24)中，n為圖像的深度，MSE為圖像間的均方差.

從圖5中可以看出，隨著迭代次數(shù)的增加，PSNR先迅速增加后緩慢減小，在30次附近PSNR達(dá)到最大值，說明本文算法能夠迅速提升圖像的圖像質(zhì)量，而且多次迭代后，圖像質(zhì)量的損失速度緩慢，不易造成圖像的模糊現(xiàn)象，從而有效地保護(hù)了圖像中的邊緣信息.

圖5 峰值性噪比與迭代次數(shù)的變化曲線Figure 5 Relationship between PSNR and Iteration

為了驗證本文去噪算法的穩(wěn)定性，圖6為算法的迭代次數(shù)為16次時，高斯噪聲方差分別為2=0.02,0.04,0.08的去噪結(jié)果.從圖中看出，隨著噪聲強(qiáng)度的增加，在迭代次數(shù)一定的情況下，雖然SLBM的去噪效果有所減弱，但是總體的去噪結(jié)果還是較為令人滿意的，特別是圖6(c)中，圖像受到嚴(yán)重的噪聲污染，憑借人眼已然很難分辨出圖像細(xì)節(jié)，通過SLBM去噪，使得圖像的質(zhì)量有了很大的提升.說明本文算法適用于不同噪聲強(qiáng)度下的快速去噪.

圖6 不同高斯噪聲強(qiáng)度下SLBM的去噪效果Figure 6 Results of SLBM for different intensity of Guass noise

4 結(jié) 語

本文在各向同性擴(kuò)散模型的基礎(chǔ)上，提出了基于Sobel算子的格子玻爾茲曼模型去噪算法.實驗結(jié)果顯示，本文提出的SLBM不僅能夠?qū)崿F(xiàn)圖像的快速去噪，又能有效保護(hù)圖像的邊緣信息，而且適用于不同強(qiáng)度高斯噪聲的快速去噪，因此具有一定的實際應(yīng)用價值.

[1] WEICKERT J. A review of nonlinear diffusion filtering[J].Lecture Notes in Computer Science,1999(1):1-28.

[2] AUBERT G, KORNPROBST P. Mathematical problems in image processing: partial differential equations and the calculus of variations[J]. Applied Mathematical Sciences,2006,147:84-107.

[3] KOENDERINK J J. The structure of images[J]. Biological Cyberernetics,1984,50(5):363-370.

[4] PERONA P, MALIK J. Scale-space and edge detection using anisotropic diffusion[J]. Pattern Analysis & Machine Intelligence IEEE Transactions on,1990,12(7):629-639.

[5] KICHENASSAMY S. The perona-malik paradox[J]. Siam Journal on Applied Mathematics,1997,57(5):1328-1342.

[6] CATTé F, LIONS P L, MOREL J M, et al. Image selective smoothing and edge detection by nonlinear diffusion[J].Siam Journal on Numerical Analysis,1992,9(1):845-866.

[7] 郭照立,鄭楚光.格子Boltzmann方法的原理及應(yīng)用[M].北京:科學(xué)出版社,2009:1-15.

[8] CHANG Qianshun, YANG Tong. A Lattice Boltzmann method for image denoising[J]. IEEE Transactions on Image Processing a Publication of the IEEE Signal Processing Society,2009,18(12):2797-2802.

[9] BHATNAGAR P L, GROSS E P, KROOK M. A model for collision processes in gases. I. small amplitude processes in charged and neutral one-component systems[J]. Physical Review,1954,94(3):511-525.

[10] QIAN Y H, D’HUMIRES D, LALLEMAND P. Lattice BGK models for Navier-Stokes equation[J]. Europhysics Letters,1992,17(6):479-484.

[11] GONZALEZ R C, WINTZ P. Digital image processing[J]. Prentice Hall International,2001,28(4):484 -486.

[12] ZHAN Ye. Lattice Boltzmann based PDE solver on the GPU[J]. Visual Computer International Journal of Computer Gaphics,2008,24(5):323-333.

A lattice Boltzmann model for image denoising

CHEN Rongqian，NIE Deming

(College of Metrology and Measurement Engineering, China Jiliang University, Hangzhou 310018, China)

A lattice Boltzmann model for fast denoising and image edge information protection was proposed based on the Sobel gradient operator. By solving the nonlinear diffusion equation, the model realized image denoising. The adoption of the Sobel operator to detect edges met the requirement of forward diffusion and guaranteed the stability of the equation. Experimental results show that the method is good for fast image denoising while protecting the edge information.

image denoising; diffusion equation; lattice Boltzmann model; Sobel gradient operator

2096-2835(2016)03-0345-05

10.3969/j.issn.2096-2835.2016.03.018

2016-04-19 《中國計量大學(xué)學(xué)報》網(wǎng)址：zgjl.cbpt.cnki.net

國家自然科學(xué)基金資助項目(No. 11272302),浙江省自然科學(xué)基金資助項目(No.LY15A020004).

陳榮前(1991- )，男，浙江省溫州人，碩士研究生，主要研究方向為顆粒兩相流模擬等.E-mail:Chrischenwork@163.com

聶德明，男，副教授.E-mail: nieinhz@cjlu.edu.cn

TN911.7

A