陳榮前，聶德明

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

?

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

陳榮前，聶德明

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

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

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

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

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

1 格子玻爾茲曼模型

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

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

宏觀流體密度和速度u的計(jì)算表達(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)相對(duì)應(yīng)的常量.表示無量綱松弛時(shí)間.

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

·u=0，

(4)

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

(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)做泰勒張開，對(duì)fi(x,t)按照Chapman-Enskog展開，然后對(duì)比方程兩邊相同階數(shù)的結(jié)果：

(8)

(9)

(10)

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

(11)

(12)

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

(13)

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

(14)

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

(15)

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

(16)

(17)

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

2 松弛因子設(shè)計(jì)

2.1 Sobel梯度算子

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

圖2 Sobel模板

(18)

(19)

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

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

(20)

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

圖3 Sobel算子邊緣檢測(cè)效果Figure 3 Results of Sobel operater

2.2 基于Sobel算子的松弛因子

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

(21)

上式中L為一個(gè)常數(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 實(shí)驗(yàn)分析

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

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

(23)

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

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

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

(24)

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

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

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

為了驗(yàn)證本文去噪算法的穩(wěn)定性，圖6為算法的迭代次數(shù)為16次時(shí)，高斯噪聲方差分別為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算子的格子玻爾茲曼模型去噪算法.實(shí)驗(yàn)結(jié)果顯示，本文提出的SLBM不僅能夠?qū)崿F(xiàn)圖像的快速去噪，又能有效保護(hù)圖像的邊緣信息，而且適用于不同強(qiáng)度高斯噪聲的快速去噪，因此具有一定的實(shí)際應(yīng)用價(jià)值.

[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 《中國計(jì)量大學(xué)學(xué)報(bào)》網(wǎng)址：zgjl.cbpt.cnki.net

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

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

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

TN911.7

A