1、高斯模糊图像的正则逆扩散方程复原方法摘要：利用高斯卷积和线性扩散的等价性，从偏微分方程逆问题的角度，提出了一种针对高斯模糊图像的复原方法：RBD-PDE（Regularized Backward Diffusion） ；从频率域角度分析了逆扩散方程的正则化表达式和正则滤波之间的关系；得出正则滤波器最佳截止频率和反向扩散时间之间的关系，为以实验的方式进行盲反卷积提供便利。较传统的基于能量范涵的复原方法，如维纳滤波或 TV 模型，RBD-PDE 方法具有最佳复原效果（在高斯核标准方差已知或未知的情况下，RBD 的结果均优于传统能量泛函方法的最佳结果） 。关键词：偏微分方程，逆问题，正则化，图像复原
2、1 引言图像复原是图像处理中的经典问题，对于线性系统，图像的模糊过程可以看作原始的清晰图像 与核函数（本文假设高斯核 ）的卷积，而图像复原或反卷积（去卷积）是0uG从模糊图像 复原清晰图像的过程，数学形式为：t(1)00tuudx=y图像反卷积包括核函数已知与核函数未知（盲反卷积）的两种情况，已有大量的研究文献提出了多种方法，如文献13,7为核函数已知情况，文献46,10,12为核函数未知的情况等。大多数方法都基于能量泛函理论，通过加入约束条件建立优化模型，如维纳滤波方法、有约束的最小二乘法、整体变分（TV）模型 7 等，或使用自然图像的统计特性取代梯度 4,6,10,11,12，以实现稳定和
3、准确的进行求解。对于基于能量泛函的方法，准确知道核函数对于复原效果起着至关重要的作用 13。当高斯核函数的标准方差未知时，有无数组 满足式(1)，因此，需要加入对0,Gu的假设（先验知识） 。稀疏性是最常用的假设，对于主要应用于运动模糊的盲卷积能0,Gu取得较好的效果 4,5,6,10。但是当稀疏性不满足时，例如高斯核函数，传统的基于稀疏先验的方法难以取得较好的效果。不同于传统的基于能量泛函的方法，本文从偏微分方程和逆问题的角度出发，提出一种全新的针对高斯模糊图像的复原方法：RBD-PDE（Regularized Backward Heat Diffusion） 。较之于传统的基于能量泛函的方
4、法，RBD-PDE 在高斯核标准方差未知的情况下，仍然能够有效地实现图像复原，性能优于传统方法。RBD-PDE 容易和现有的线性或非线性偏微分方程方法相结合，构成新的复原模型，因此具有更大的灵活性和方法的可拓展性。对于运动模糊图像已有许多有效的复原方法 4,5,10，而复杂的图像模糊可以分解为运动模糊和高斯模糊 ？ ，并分别进行复原。因此，高斯模糊的复原方法具有很重要的研究和实用价值。2 正向和逆扩散方程线性热扩散方程的解为高斯核（热核）函数与初始条件函数的卷积，热扩散过程等价于高斯模糊过程。自然的，从偏微分方程角度看，图像复原可视为正向热扩散的逆过程。2.1 正向扩散方程对于线性热扩散偏微分
5、方程：（2）20, 0utttRxx其中 ，是一个二维变量。 为热扩散方程的初始条件，对于图像问21,xR=题，表示原始的清晰图像。假设 定义在区间 上，式(2) 的解为 9：0ux2,1(3)00,tGud=xy表示二维高斯核函数：2t(4)221221expep的标准差。在已知 的条件下，图像复原问题等价于求解式（3）的 （第一类 Fredholm0u积分方程） 。现有许多求解方法，例如 LTI（线性时间不变）维纳滤波：(5)*10 2, tFGux其中 分别表示傅里叶变换及逆变换算子， “ ”表示算子的作用， “*” 表示算子的1,F伴随（共轭转置） ， 为正则参数，一些改进复原方法 2
6、,3，可以看作对 的优化。若式(2)中加入边界条件的约束，式(5)中的 可为傅里叶变换及逆变换的特殊形式。例如对于1,F第二类边界条件（本文中使用的边界条件） ， 表示余弦变换及逆变换。1,2.2 逆扩散方程在核函数未知的情况下，无法直接通过式(5)求解 。但是可以从式(2) 出发，通0,ux过逆过程，得到 的估计值，即将盲反卷积问题转换为一个偏微分方程逆问题。0ux引理 1：假设模糊图像是经过线性扩散方程 (2)（式(3)高斯核卷积）得到的，在理想情况下（没有噪声和计算误差） ，总可以通过逆扩散方程：（6）201, 0,vtttRuxx得到式(2)中的初始条件 。0u证明：对式(2)和式(6
7、) 两端做傅里叶变换，可得：（7）20,0dtutu（8）201,0dvtvtu，是一个二维频率域变量。 表示 的傅里21,R=,tt,utvtx叶变换；求解常微分方程(7)和(8) ，最终可得：（9）201,expvt t因此，当 时， ，从而 。并且，在满足1t,tu0,vt的条件下， ，即在满足引理 1 的假设条件下，总可以通过实验的1t0,v方法稳定地求得原始的清晰图像。2.3 逆扩散方程的病态性在噪声存在的情况下（本文假设噪声 是方差为 加性高斯白噪声） ，即式(6)中：2（10 ）01,vutxx进而可以得到噪声情况下式(6)的结果为：（11 ）,exacttnt表示无噪声精确解，
8、 为噪声放大项：,exactvn（12 ）2,ptt服从期望为 0 方差为 的高斯分布，由式 (12)可以看出，噪声中的高频分量将被迅速放2大（以指数速度） ，覆盖真实（希望得到的复原）结果，使得实际应用中，无法直接利用式(6)进行图像复原。3 正则逆扩散方程利用逆扩散方程进行图像复原最早由 Garbor 提出。由于逆扩散方程的病态性，噪声被迅速放大，从而使得方程(6)的应用受到极大的制约。从逆问题理论出发，通过正则化方法控制噪声的传播和放大是解决式(6)的病态性的有效途径。3.1 正则化方法如果要求复原结果中的噪声放大项 满足：,nt（13 ）2expt对于任意 都成立，那么逆扩散的极限时间
9、 需要满足：(无推导)( 随即变量不能取具体mat值, ?)2（14 ）max21lnt由式(14)可以看出， 由噪声的方差 决定，与正向扩散时间 无关。当 远小maxt2 1tmaxt于 时，无法用方程式(6) 进行有效复原。因此，增加逆扩散方程实用性的关键在于降低1t对 的限制，使 接近 ，解决的途径之一是强制使得式(12) 中较大 所对应的值为maxat1t 0，即放弃对原始图像中高频分量的恢复来换取计算过程的稳定，而这正是正则化的基本思想。采用低通滤波器去除图像中高频分量，如使用截断窗口函数， （15 ）if 0,式(15)可滤除图像中的高频分量，称为正则滤波器， 为对应的截止频率。对
10、 进行滤波，逆扩散结果（式 11）变为：,vt（16 ）, , ,exactvvnt其中 ， 分别为：,exactnt（17 ）2, 01pexact ut（18 ）2,ex对于扩散方程，正则化等价于低通滤波过程，正则化后，逆扩散的极限时间 变为：maxt（19 ）max21lnt提高了逆扩散时间，扩大了逆扩散方程的适用范围。如同所有的正则化一样，正则化要付出一定的代价，正则滤波器 是以丢失高频信息为代价，换取逆扩散极限时间的增加，因此，其复原效果自然要受到一定的制约。下面给出 5050 的二值图像“矩形”的实验，在这里的计算采用有限差分法的显式格式，高斯白噪声的方差 ，实验结果如同图 1 所
11、示：220.1(a).原始图像“矩形” (b).模糊加噪图像 (c).式(6) 的最佳复原结果 (d).正则化后的结果10 20 30 40 501020304050 -6-4-20210 20 30 40 501020304050 -505(e).图像(a)的 DCT 系数( 对数显示) (f)逆扩散后的图像的 DCT 系数( 对数显示)图 1 正则化过程和结果在图 1 中，(c) 为未正则化时的最佳复原结果，可以看出，没有明显的复原效果，相反，噪声被迅速放大；(d)为正则滤波后的逆扩散结果，截至频率 ，迭代时间与正向扩散相同，和(b) 相比有明显的复原效果，噪声得到抑制；(e) 为(a)
12、的 DCT 系数（以对数形式显示） 。(f) 为正则化前图像的 DCT 系数（对数形式显示） ，可以看出，能量主要集中在高频（右下角） ，说明高频噪声被迅速放大。实验中采用归一化均方差：（20 ）0kEvu作为(c) 选取的标准。由于截断窗口函数式(15)不连续，滤波会产生吉布斯现象（边缘附近出现振荡，产生重影） 。选取其它函数来优化窗口函数可以解决该问题，例如，选取亚高斯窗口函数：（21 ）4exp作为低通滤波器，并对复原后的图像进行投影：（22 ）1,if,1, 0,i,vtvtttx可提高图像处理效果。图 2 是在相同实验条件下，优化后得到的实验结果：2040204000.5120 40
13、204000.51(a) (b) (c)图 2 亚高斯窗口滤波及投影后的结果图 2 中(a)为亚高斯窗口（式(21)）正则滤波后经过式(22)进行投影后的结果；(b)为(a)的灰度值所对应的曲面；(c)为图 1(d)的灰度值所对应的曲面。对比(b)和(c) 可以看出，通过调整窗函数和投影，吉布斯现象得到了有效的抑制，图像效果有所提高。3.2 正则逆扩散框架从频率域分析的角度考虑，正则化方法实际上是对逆扩散方程的解进行低通滤波，以抑制噪声中高频分量的传播和放大，降低对逆扩散时间的限制，达到更好的复原效果。即：（23 ）21, 0expvFtFvx如果令：（24 ）100可以反解出正则复原方程，写
14、为：（25 ）200, ,vtttRLvxx其中 表示低通滤波器 所对应的滤波算子。说明对于 LTI（线性时不变）的情况，理L论上，以 为初始条件的正则逆扩散方程等价于 【 为初始条件】的逆扩0vx0x0v散方程，即逆扩散方程的正则化等价于对扩散方程初始条件进行低通滤波。在数值计算中，由于存在机器误差，当迭代次数较大时，仍然需要对结果再进行低通滤波，以消除误差的积累。这种情况下，正则复原过程等价于：（26 ）21 0, 1expvFtFvx作为改进，可以考虑使用时变正则低通滤波器：（27 ）4,exptt抑制迭代过程中计算误差的传播和放大，由式(8)和式(17)可以求得对应的时变 RBD PD
15、E:为什么 （28）2200, 0,vttvttRLxx其中 表示四阶偏导数算子。为了减小式(27) 中时间 变化的影响，可以采24211xx t用：为什么（29）22001, 0,vttvttRLxxRBD-PDE(式 29)所对应的时变滤波器为：取何值 （30 ）4,eplntt在 较小的条件下式(30)近似于式(21)。 什么时候停止t可以看出，RBD-PDE 是一个理论框架，可以根据本文提出的原理和实际图像复原特点，建立新的偏微分方程。同时，可以和一些现有的非线性方法，例如文献14等相结合。因此，RBD-PDE 具有很强的灵活性和可扩展性。3.3 频率域分析逆扩散过程中，噪声项的高频分
16、量被迅速放大，正则化的过程是通过低通滤波抑制噪声。由于图像的边缘也是高频成分，因此，正则化的过程同时也抑制了图像边缘的复原，换句话说，正则化需要在问题解的稳定性和精确性之间寻求一种平衡。正则化方法的优化问题等价于滤波器的设计问题，我们通过频率域分析，得到了逆扩散时间和正则滤波器最优截止频率的关系，为卷积实验提供了便利。为了便于讨论，后面的分析中，正则滤波器采用式(21)，于是，式(17) 可以进一步整理为：（31 ）224, 0 011expexactvutu 从频率域的角度，由原始清晰图像，经过高斯模糊（正向扩散） ，再经过正则逆扩散，得到正则复原结果中的整个过程可以看作一个等效滤波器 ，将
17、逆扩散的极限时间(式19)代人， 可以整理为：（32 ）221 1242412expln, ifln , iflt tt 在式(32)中， 的形状由正则滤波器的截止频率 唯一确定，合理的 可以增大的宽度，达到尽可能多地恢复原始图像信息的目的。图 3 给出了滤波器(式 32)随截止频率的变化曲线。1020304050102030400.20.40.60.8|w|pE(|w|,p)图 3 滤波器( 式 32)随截止频率的变化从图 3 中可以看出，当 较小时（ 的宽度加宽） ，丢失了大量的高频信息，当 较11大时，由于式（19）的约束，反向扩散作用微小， 近似于正向扩散的窗函数（高斯函数） 。只有当
18、（33 ）12lnopt取得最佳效果（ ） 。因此，对于 的信息，无法在噪声约束条件(12)下1opt t进行有效恢复。式(33)刻画了正则复原方法的极限，只有当原始清晰图像的频率分量集中在 的范围内时，正则复原方法才是有效的。1opt由于正向扩散时间 未知，实际情况下，无法精确得到 。但是，数值实验中，反1t 1opt向扩散时间 ，因此，式(33)约束了实验中反向扩散时间 和截止频率的关系，为盲卷t积实验提供了一种快速的搜索途径。从图 3 中可以看出， 被低估的影响小于被高估，因此，截止频率可以适当选择稍小一些。4 实验考虑到本文提出的 RBD-PDE 属于线性方法，而维纳滤波是噪声特性已知
19、情况下的最优线性反卷积方法 1，也是常用的图像复原方法。因此，本文的实验采用 LTI RBD-PDE 和 LTI维纳滤波方法进行比较。另外，虽然非线性方法，如 T-V 模型 7及其改进方法 16比维纳滤波复原效果要理想，但对于模糊程度较大的图像，与 LTI RBD-PDE 相比复原效果较差。实验中，采用有限差分法的显式格式对 RBD-PDE 进行离散化数值求解，为了说明 RBD -PDE 方法的通用性，正向扩散时，时间步长与空间间隔平方之比取为 ；210.1rtx反向扩散时，取 。对于 1 的信号或 的图像块，220.rtx。由于高斯核的标准差 ，真实高斯核和估计的高斯核的标准差1xN2t满足
20、：2,(34)1122trn其中 为正向和反向扩散时间， 为对应的迭代次数。由于 ， 表12,t 1,nmxR示像素点数目，可以进一步得到 和迭代次数 之间的关系： 。m21,2iirn4.1 二值图像实验对于二维图像，我们仍以 5050 的图像“矩形” （图 2(a)）进行实验，正向扩散时，；对应于 个像素点的高斯核。逆扩散时， ，分别取噪声方110.,8rn420.r差 。实验结果如图 4 所示：222,0.5(a).模糊图像( ) (b).式(24)得到的结果 (c).LTI RBD 结果( 对(a) (d).维纳滤波结果(对(a)0(e).模糊图像( ) (f).LTI RBD 结果(
21、对(e) (g).维纳滤波结果(对 (e) (h).模糊图像( )0.1 0.5(i).LTI RBD 结果(对(h) (j).维纳滤波结果(对(h) (k).维纳滤波结果(m=5) (l).LTI RBD 结果(n 2=63)(m).T-V 模型结果( 对(a) (n). T-V 模型结果 (对(e) (o). T-V 模型结果(对(h) (p). T-V 模型结果(m=5)图 4 正则锐化方程实验结果图 4 中(a)为 时的模糊图像，(b)为直接由式(24)得到的复原结果，由于迭代次数较多，20机器误差的积累使得结果中出现震荡失真；(c)为对(b) 正则滤波后的结果，取 ，震120荡得到有
22、效抑制；(d)为对应的维纳滤波方法的最优结果（(2.8)为标准） ；(e)为时的模糊图像，(f)为(e)正则复原后的结果（ ） ；(g)为对应的维纳22.1 15滤波方法的最优结果；(h)为 时的模糊图像，(i)为(h)正则复原的结果（220.5）;(j)为对应的维纳滤波方法的最优结果；(k)为高斯估计过大时（m=5 ） ，(h)经10过维纳滤波后的最佳结果；(l) 为相同条件下(t=(5/4) 2t1)，LTI RBD-PDE 的复原结果。实验表明，对于高斯核已知情况，LTI RBD-PDE 结果优于维纳滤波结果，偏离已知的高斯核对维纳滤波结果影响较大，但 LTI RBD-PDE 在相同条件
23、下仍能取得较好结果。作为对比，(m),(n),(o)为 T-V 模型的复原结果（分别对(b),(e),(h)） 。可以看出，对于模糊程度较大的情况，T-V 模型虽然有较清晰的效果，但是从图像复原角度看，结果却不理想（方形变成了圆形） 。(p)为高斯估计过大时（m=5） ，T-V 模型对于(h)的复原结果。可以看出，对于噪声特别小，模糊程度大的情况，T-V 模型与线性方法相比，没有优势，但是，对于噪声较大的情况，T-V 模型的非线性特性使得其对噪声放大作用的抑制效果比较好。因此可以借用 T-V 模型的非线性原理和思想，对线性 RBD-PDE 进行改进，加入非线性正则约束项，以更好地抑制噪声，从而
24、提高图像的复原效果。4.2 自然图像实验T-V 模型的“分片常数”特点，使得其对二值图像处理效果理想，但是对于自然图像复原效果并不理想，尤其是对模糊程度较大的图像。本节中，分别采用 LTI RBD-PDE，维纳滤波和 T-V 模型对自然图像进行实验。原始清晰图像选取 512512 的自然图像 “Boat”，正向扩散取 ，对应 个110.,rn20m像素点的高斯核，逆扩散时， ，实验结果如图 5 所示：20.r(a).原始清晰图像 (b).模糊图像(n 1=100)(c).LTI RBD 复原结果 (n2=50; t1= t2) (d).维纳滤波结果(核宽度估计准确 )20m(e).LTI RB
25、D 锐化效果 (n2=63; t1=1.26t2) (f). 维纳滤波结果 (核宽度估计过大 m=5)(g).T-V 模型结果(核宽度估计准确 ) (h). T-V 模型结果(核宽度估计过大 m=5)20m图 5 自然图像的正则复原结果从图 5 中可以看出，对于自然图像，LTI RBD-PDE 方法取得了较为理想的效果，当高斯核估计正确时，复原效果优于维纳滤波方法的结果（对比图 5 (c)和(d)） ；当高斯核估计较大时（m=5,对应的反向扩散时间为正向扩散时间的 1.25 倍，实验中取逆扩散迭代次数 n2=63） ，图像得到很好复原。由于逆扩散迭代次数大于造成模糊的正向扩散迭代次数，使得图像
26、获得 A 了进一步锐化的效果，看起来甚至比原始图像(a)更清晰。而在相同情况下，维纳滤波的结果则差很多。从图 5 (g)和(h)可以看出，T-V 模型由于其“分片常数”性，图像中很多细节和纹理在锐化过程中被“抹平” 。虽然 T-V 模型是一个很有影响，且有较好图像复原效果，但对于模糊程度大的图像，T-V 模型并不适用。实验说明 RBD-PDE 是一种有效的针对模糊程度较大图像的复原和锐化方法，相对于维纳滤波，RBD-PDE 有更高的灵活性和更好的效果，相对于 T-V 模型，有更广泛的适用性。5 结论处理图像的线性扩散方程等价于对图像的高斯模糊，于是，高斯模糊图像的复原可以看作一个偏微分方程的逆
27、问题，而逆扩散方程是线性扩散的逆过程。由于逆扩散方程的病态性，噪声和计算误差被迅速放大，使得用逆扩散方程来实现模糊图像复原并无实用性。本文从逆问题的角度出发，通过对图像逆扩散频率域分析，得到了逆扩散方程的正则化方法，给出了正则逆扩散方程（RBD-PDE） 。利用正则逆扩散方程可有效地对高斯模糊图像进行复原。实验表明，本文提出的 RBD-PDE 不仅在恢复图像的效果上有显著的提升，而且相比于传统的基于能量泛函的方法，在核函数大小未知的条件下，本文的方法具有更好的稳定性和更佳的效果。参考文献1 M. Bertero， P. Boccacci, Introduction to inverse pro
28、blems in imaging, Philadelphia, Pa.: Institute of Physics Pub., 19982 Mesarovic V.Z., Galatsanos N.P., Katsaggelos A. K. “Regularized Constrained Total Least-Squares Image Restoration”, IEEE Trans. on Image Processing, 4(8), 1096-1108, 1995.3 Banhan M, Katsaggelos AK, “Spatially adaptive wavelet-bas
29、ed multiscale image restoration.” IEEE Trans. on Image Processing, 5(4), 619634, 1996.4 Qi Shan, Jiaya Jia, and Aseem Agarwala, “High-quality Motion Deblurring from a Single Image”, SIGRAPH, 2008.5 Amit A, Yi Xu, Ramesh R, “Invertible Motion Blur in Video”, SIGGRAPH, 20096 Levin A., “Blind motion de
30、blurring using image statistics”, In NIPS, 2006.7 L. Rudin, S. Osher, E. Fatemi, “Nonlinear Total Variation based noise removal algorithms”, Physica D, 60, 259-268, 19928 Isakov V., Inverse problems for partial differential equations, New York : Springer, 2006.9 G. Aubert, P. Kornprobst, Mathematica
31、l problems in image processing: partial differential equations and the calculus of variations, New York: Springer, 200610 Fergus, R., Singh, B., Hertzmann A., Roweis S. T., Freeman W. “Removing camera shake from a single photograph”. ACM Trans. on Graphics, 787-794. 2006. 11 Y. Weiss, W. T. Freeman,
32、 “What makes a good model of natural images?”, In CVPR, 2007.12 N. Joshi, R. Szeliski, D. Kriegman, “Psf estimation using sharp edge prediction”. In CVPR, 2008.13A. Levin, Y. Weiss, F. Durand, W. T. Freeman “Understanding and evaluating blind deconvolution algorithms”, In CVPR, 2009.14 Perona P, Malik J., “Scale-space and edge detection using anisotropic diffusion”, IEEE Transactions on Pattern Analysis and Machine Intelligence, 12 (7). 629-639,199015 R. C. Gonzalez, R. E. Woods. Digital Image Processing, Prentice Hall, 2002.
