偏微分方程在图像去噪中的应用_王正明

应 用 数 学
M A T H EM AT ICA A P PL ICAT A
2005,18(2):219~224
*偏微分方程在图像去噪中的应用
王正明,谢美华
(国防科技大学理学院系统科学与数学系,湖南长沙410073)
庆贺陈庆益先生八十寿辰
摘要:本文介绍用于图像去噪的偏微分模型、方法的发展历程.从理论上分析了线性
模型、简单非线性模型、复杂非线性模型、多步处理模型出现的背景和优缺点,并从空
域和频域上对偏微分方程模型的去噪原理进行了分析.最后,指出了偏微分方程去噪
与小波去噪结合的途径,据此对偏微分方程未来的发展方向进行了展望.
关键词:偏微分方程;模型;扩散;正则化;去噪
中图分类号:T P391 AMS(2000)主题分类:35R
文献标识码:A 文章编号:1001 9847(2005)02 0219 06
光学图像成像过程中的噪声污染通常来自于CCD成像器件、图像传输设备及处理设备的背景光子散粒噪声、暗电流散射噪声、读出噪声、热噪声、放大器噪声等,其中最终起限制作用的是光子散粒噪声,统计上服从高斯或泊松分布.同时图像数字化过程中的量化以及其它人为的因素也会导致噪声的产生,这种CCD传感器噪声和量化噪声可以仿真成 加性的或乘性的 , 信号相关的或信号无关的 以及 有色的或白色的 ,它们的存在极大地影响了图像的质量.这里只考虑加性的与信号无关的白噪声.此时,图像的成像模型可描述为
g(x,y)=f(x,y)+ (x,y),(1)其中f(x,y)为不含噪的真实图像,g(x,y)为实际观测到的图像, ~N(0, 2).那么,图像去噪问题就相当于,寻找合适的算子F R R,使得
F(g(x,y))=f(x,y).(2)由于噪声是随机性的,事实上我们只能得到f(x,y)的近似估计,而不可能使(2)式完全成立.因此,在评价算法的优劣时,通常以下述峰值信噪比作为评价指标.
PSNR=10*lg
2552 m n
i,j
(f(i,j)-f*(i,j))2
,(3)
其中f*(x,y)为采用某算法去噪后的图像,m,n为图像的尺寸.传统的图像去噪方法,如中
*收稿日期:2004 05 08
基金项目:全国优秀博士论文作者专项基金(200140),国家自然科学基金资助项目(60272013)
作者简介:王正明,男,汉,湖南长沙人,教授,博导,国防科技大学理学院院长,主要研究方向为图像处理中的数学方法、装备试验分析与评估.
值滤波、均值滤波,主要将图像的高频成分滤除.由于图像的细节如边缘纹理等也分布在高频区域,所以总是在对噪声进行滤除的同时将图像的边缘部分模糊了.事实上,数字图像在本质上可看成是R 2 R 的以图像的边缘为边界的分片连续的映射.基于这一性质,可以以图像的边缘为边界采用分片连续的函数来逼近图像中的真实信号,抑制其中的随机性噪声,如二元多项式就是一种可选的基函数,由于逼近是在图像的内区域进行的,所以不会造成对边缘的模糊,图像去噪的偏微分方程(PDE)方法就是基于这一性质构建的.
令u 0 R 2 R 表示一个含噪的灰值图像.对u 0(x ,y )去噪相当于极小化能量函数
[1]E (g)= 2 (u -u 0)2d x d y + |u 2x +u 2y |d x d y ,(4)
其中 为正则化参数, 为图像的支撑域, 称为正则化参函数,是梯度的增函数,满足 0.
(1)式中的前一项用以约束去噪图像与原图像的逼近程度,而后一项用以约束图像的光滑程度.它的含义是将图像近似为分片连续的零阶多项式.对极值问题(4),应用Eluer 方程,并采样最速下降法引进时间参数t,可转化为求解如下偏微分方程
u t
=F[u(x ,y ,t)],(5)其中u(x ,y ,t) R 2 [0, ] R 是随时间演化的图像,F R R 表示一个特定算法对应的算子.根据F 定义的不同可分为线性扩散、非线性扩散、各向异性扩散等.其中,最早出现的是如下线性扩散模型 t u = u,
u(x ,y ,0)=u 0(x ,y ),
(6)它对应于方程(4)中 =0, (s)=s 2/2的情形.
(6)式是一个解热传导方程的Cauchy 问题,通过对其进行Fourier 变换,可求得其解为
u(x ,y ,t)=u 0(x ,y )(K *u 0)(x ,y) t =0t >0,K (x ,y ,t)=1(4 t )3/2ex p -x 2+y 24t ,(7)显然这就是对应图像的高斯光滑,即利用高斯函数对邻域内的点进行加权平均来实现去噪.
(7)式给出的解可从空域上理解(6)式去噪的含义,下面从分离变量法的角度给出方程(6)的解[7],从频域上理解(6)式去噪的含义.设初始图像u 0(x ,y )按正弦函数可以展开成如下级数
u 0(x ,y )=
M k =1 N l =1A
k,l sin k x M sin l y N ,(8)其中N ,,M 分别为u 0(x ,y )沿x ,y 方向的离散采样点的个数.则(6)式的解可表示为u(x 1,x 2,t)= M k=1 N l =1A k,l e -k 2 2t M 2+l 2 2t N 2sin k x 1M sin l x 2N .(9)
从(9)式可以看出,经过(6)式去噪后的图像的正弦级数系数等于u 0(x ,y )的正弦级数的展开系数乘以一个与扩散时间相关的压缩因子w (k,l)=e -k 2 2t M 2+l 2 2t N 2,因为随着k,l 的增大w (k ,l)不断减小,因此(9)式对u 0(x ,y )的高频成分保留很少,因而能够实现对噪声的抑制.
由于方程(6)采用的正则化参函数 (s)=s 2/2,这就要求图像的梯度在整幅图像上实现最小,当扩散时间T 趋向于无穷时,得到的是用一个常数对图像进行逼近的结果,这无疑会导致对图像如边缘等微细结构的模糊.为了避免这种模糊的产生,人们采用了一种非线性的正则化参函数,利用分片连续的平面对图像进行逼近.此时参函数 的选取可以有很多形式,如 (s)=(k 2/2)log (1+(s/k)2),或220应 用 数 学2005
当 (s)=(k 2/2)log (1+(s/k )2)时,(5)式对应典型的非线性扩散方程-P M 扩散
[2],此时 t u =div (g(| u |2) u),
g(| u |2)=11+| u |2/k 2
,u(x ,y ,0)=u 0(x ,y ), 或 t u =div (g(| u |2) u),g(| u |2)=11+| u |2/k 2,u(x ,y ,0)=u 0(x ,y ).
(10)本着在边缘处扩散系数g(s)小的准则,一般选取g(s)使之满足如下条件:g(0)=1,g(s)为减函数,且lim s
g(s)=0.而且,在参数的选取上一般选取 为严格凸函数,这是因为定理1 当势能函数 (| u |)为严格凸函数时,能量函数E(u)=
(| u |)d x 正好有一个最小值.
证 令 (| u |)= ( ( u))=g(| u |2) u,由 (| u |)为严格凸函数可知 严格单调递增,又 (0)=0,从而当| u |>0时,有 (| u |)>0,从而 (| u |)是严格单调递增的.对一个常值图像 u 而言,其梯度恒等于零,而 (| u |)是关于| u |的严格单调增函数,因此任给u,| u | 0,必有 (| u |)> (| u |),又因为 (| u |) 0,从而有E(|u |)>E (| u |),因而 u 是E (u)的唯一最小值.定理证毕.
当定理1成立时,不必进行非凸优化所涉及的复杂的计算,而可使用标准的有限元近似来获得稳定的解.(10)式所示的非线性方程能根据图像的梯度| u |来判断边缘位置,使边缘处的扩散系数小,降低对边缘的模糊程度.但是,在| u |> 的情况下,其势能函数是非凸的,从而使得边缘等处的处理表现出不稳定性.而且,由于边缘处扩散系数很小,边缘处的噪声得不到有效抑制.
由于正则化参函数 选取的不恰当,线性模型与非线性模型都有各自的缺点.为克服这一缺点,在后续的研究中人们将更多的精力投入到了 函数的选取上.其发展主要沿两个方向进行.第一,由简单方程到复杂方程的转变,这种复杂方程的复杂性主要体现在参函数的形式上,如考虑采用高阶方程[3]、逆扩散方程[4],以及增加约束条件[5]等;第二,由偏微分方程的一步实现到多重实现的转变,这种多重实现指的是多次运用简单方程来实现复杂的操作[6,7].目前,偏微分方程去噪研究的重点仍放在这两方面.但是,由于采样复杂的参函数会增加约束参数的个数和模型的复杂性,给处理带来较大的麻烦,所以其前景较偏微分方程的多重实现差.
高阶偏微分方程 在能量函数(4)中,考虑对图像的二阶导数的约束有如下能量泛函[3]
E(u)=
f (| u |)d (11)时它所对应的偏微分方程为
u t
= [g(| u |) u].(12) 这是一个四阶偏微分方程,它的含义是将图像近似为分片连续的一阶多项式.使用四阶方程较前面的二阶偏微分方程的好处是能克服二阶方程常出现的 块状 (blocky )的效果.
前向 后向扩散方程 考虑逆扩散方程
t u (x ,t)=-c 2u(x ,t),c >0,(13)
这等价于一个高斯去卷积过程,当其作用于边缘附近时具有锐化边缘的作用.这一方程由于其数值不稳定性而被认为是病态的,但是Gilboa,Zeevi,Sochen 认为在局部范围内使用逆扩散过程不会破坏其稳定性,因此,提出了一种前向 后向扩散方程[4]
.即取扩散系数g(s)满足221第2期
王正明等:偏微分方程在图像去噪中的应用
g(s)=1-(s/k f )n ,
[((s -k b )/w )2m -1],0, 0 s k f ,k b -w s k b +w ,else.
(14)前向 后向扩散方程的优点是能在光滑区域内部的同时锐化边缘,缺点是参数选取困难.除了这些形式的方程所对应的参函数以外,还有复扩散方程[5]、基于边缘检测的扩散方
程、约束曲面面积扩散方程、拐点增强扩散方程,以及基于特殊边缘增强的藕合模型、稳健扩散、自适应扩散、多重网格等等.偏微分方程多步实现的代表是各向异性扩散方程.其中具有代表性的有边缘增强扩散和相干增强扩散.边缘增强扩散是为了对边缘处的噪声进行处理,而相干增强扩散是为了凸现图像的线状结构.在引入各向异性扩散模型之前,先介绍如下定理.
定理2 设u 0为原始图像,D 为连续的2 2扩散矩阵,p 1,p 2为其规范正交的特征向量分别对应图像的梯度方向和边缘方向, 1, 2为其相应的特征根,考虑如下扩散平滑问题
t u =div [D u],u(0,x ,y )=u 0(x ,y),
(15)
则利用(7)式对u 0进行光滑近似于分别以 1, 2的速度沿p 1,p 2方向光滑.
证 设p 1=(p 11,p 12)T ,p 2=(p 21,p 22)T ,由于图像的边缘方向在小范围内基本不变,所以p 1,p 2在一定范围内取常值,此时
u t
= (p 211u xx +2p 11p 12u xy +p 212u yy )+ 2(p 221u xx +2p 21p 22u xy +p 222u yy )= 1p 11 (p 11u x +p 12u y ) x +p 12 (p 11u x +p 12u y ) y
+ 2p 21 (p 21u x +p 22u y ) x +p 22 (p 21u x +p 22u y ) y = 1 (p 11u x +p 12u y ) p 1+ 2 (p 21u x +p 22u y ) p 2= 1 2u p 21+ 2 2u p 22
.因此,u 沿p 1,p 2方向扩散的速度分别为 1, 2,定理证毕.
各向异性扩散方程采用(7)式所对应的模型,其中扩散张量D 的特征向量p 1平行于梯度矢量,p 2垂直于梯度矢量.边缘增强扩散以 u u = u u T 为边缘定向算子,D 与
u u 有相同的特征向量,其中u =u( ,t)*
K ,K ( )=12 2 ex p -| |22 2.而D 的特征根为[6] 1=1, 2=l 2/(l 2+| u |2
),其中 1对应垂直于梯度方向的特征向量, 2对应平行于梯度方向的特征向量.显然该方程是先利用了一次线性方程对边缘定向,然后再利用非线性扩散方程实现了沿边缘方向的扩散,是分两次利用偏微分方程.边缘增强扩散的优点是考虑到了扩散系数沿边缘方向和垂直边缘方向的不同,缺点是所采样的边缘定向算子 u u 不能正确的对边缘定向.
与边缘增强扩散不同,相干增强扩散采用了一个新的边缘定向算子J 来对图像边缘定向,其扩散张量D 的特征向量与J 相同,其中[7]
J ( u )=K *( u u ), 0.
(16)
扩散张量D 的特征值取为
1 2,(17)
222
应 用 数 学2005
1, 2为 u u 的特征值.C 为大于零的常数, (0,1),常取一个很小的值.对应每一个J ,其中 1对应垂直于梯度方向的特征向量, 2对应平行于梯度方向的特征向量.相干增强扩散的优点是所采样的边缘定向算子能准确的描述边缘方向,缺点是其特征根的选取不适合平坦区域的去噪,会导致虚假边缘的产生.改变结构描述算子(16)的特征值的取值(17)还可以实现许多其他结构的定向增强.
为了进一步完善算子,[8]对扩散模型的最优停止时间进行了讨论,采用了一种最小化噪声协方差的去相关准则来实现最优停止.
为解决去噪问题(1),人们从不同的角度提出了很多不同的方法,其中主要的有正则化方法、小波系数伸缩法、总变分方法以及偏微分方程方法,文献[10]的研究指出这些不同方法之间存在某种程度的统一,并指出这种统一有利于综合利用这些方法的优点获得更好的去噪效果.前面也说明了偏微分扩散方程可以通过Euler 方程与正则化方法联系在一起.下面对偏微分方程与小波伸缩之间的相互关系进行讨论.
在应用最速下降法引入时间参数t 后,偏微分方程去噪每一步迭代的值可以看成空间{u t }的一个元素,可以证明{u t }满足尺度空间性质,于是偏微分方程一步扩散的结果就可以对应小波去噪的一步伸缩.文献[9]证明了H aar 小波的伸缩函数与非线性扩散方程的扩散系数之间存在一对一的关系,这种相互关系可通过下述过程获得
(i)给出偏微分方程进一步扩散的离散表达式.为方便记,我们考虑一维处理的情况,设所采用的偏微分方程的扩散系数为g,则离散化后的偏微分方程进一步扩散的表达式为
u n +1(i)=u n (i)+ (g |u n (i +1)-u n (i)|)(u n (i+1)-u n (i))
-g(|u n (i)-u n (i -1)|)(u n (i)-u n (i -1)).
(18) (ii)给出H aar 小波去噪的一步分解与重构的离散表达式.设所采用的H aar 小波的伸缩
函数为S ,则在一维情况下,基于H aar 小波一步分解与重构的离散化后地的表达式为
u n +1(i)=u n (i -1)+2u n (i)+u n (i +1)4+122S u n (i)-u n (i +1)2-122
S u n (i -1)-u n (i)2.(19)
(iii)对比(i),(ii),可得到两者之间的相互转换关系.将(18)式重写为
u n +1(i)=u n (i -1)+2u n (i)+u n (i +1)4+u n (i)-u n (i +1)4-u n (i -1)-u n (i)4+ (g(|u n (i +1)-u n (i)|)(u n (i +1)-u n
(i))
-g(|u n (i)-u n (i -1)|)(u n (i)-u n (i -1)))
=u n (i -1)+2u n (i)+u n (i +1)4
+(u n (i)-u n (i +1))(1/4- (g(|u n (i)-u n (i +1)|)))
-(u n (i -1)-u n (i))(1/4- (g(|u n (i -1)-u n (i)|))).
(20)对比(19)与(20)有
1
22S (x /2)=x (1/4- g (|x |)).从而S 与g 之间的相互关系可表述为S (x )=x (1-4 g (21)
223第2期王正明等:偏微分方程在图像去噪中的应用
g(|x|)=1
4
-
2
4 x
S
x
2
,(22)
由(21),(22)两式可在小波伸缩函数与扩散系数之间进行转换,从而可综合利用二者.
参考文献:
[1] W eickert J.A R eview o f N onlinear Diffusion F ilter ing in Scale space T heor y in Computer V ision[M].
Berlin:Spring er,1997.
[2] Per ona P,M alik J.Scale space and edg e detect ion using a niso tro pic diffusion[J].IEEE T r ansactio n on
P attern A nalysis and M achine Intellig ence,1990,12(7):629~639.
[3] L ysaker M,L under vo ld A,T ai X C.N oise remov al using fo ur th or der partial differential equatio n w ith
applicat ions to medical magnetic resonance imag es in space and time[J].IEEE T ransactions o n Imag e P ro cessing,2003,12(12):1579~1590.
[4] G ilbo a G,Zeevi Y Y,So chen N A.F orw ard and backwar d diffusion pro cesses fo r adaptiv e image enhance
ment denosing[J].IEEE T ransaction on Image P rocessing,2002,11(7):689~703.
[5] G ilboa G,Zeev i Y Y,Sochen N plex diffusion pr ocess for imag e filter ing[C].Scale Space,2001,
L N CS2106,Ber lin:Spr inger Ver lag,2001.
[6] W eicker t J.T heoret ical foundatio ns of naisotr opic diffusio n in imag e pro cessing[J].Co mputing,1996,11:
221~236.
[7] Weicker t J.Coherence Enhancing diffusio n filter ing[J].Inter national Jo ur nal o f Comput er V ision,1999,
31(2/3):111~127.
[8] M r azek P,N avara M.Selectio n o f optimal st opping time fo r nonlinear diffusio n filter ing[J].International
Jour nal of Computer V isio n,2003,52(2/3):189~203.
[9] Steidl G,W eicker t J.Br ox T,M razek P,W elk M.O n the equivalence o f so ft wavelet shr inkag e,to tal var i
ation diffusio n,total var iatio n regularizatio n,and SIDEs[J].SIA M Journal on N umer ical Analysis,2004, 42(2):686~713.
A Review of the Application of
Partial Differential Equation in Image Denoising
WAN G Zheng ming,X I E M ei hua
(Dep ar tment of Math.,Science College,N ational Univ er sity of Def ense T echnology,H u nan Changsha410073,China)
Abstract:This paper g ives an o verview o f the development of partial differential equation models for image deno ising.We sketch basic ideas behind the different filter m odels,such as linear model,simple nonlinear model,complex no nlinear m odel and its m ulti steps realiza tion,and analysis the denoising pr inciple of partial differential equation fro m spatial dom ain and fr equency dom st,w e point out a strateg y of com bining w avelet and partial differ ential equation,and give a pr eview of its development.
Key words:Partial differential equation;M odel;Diffusion;Regularizatio n;Denoising
224应 用 数 学2005。