高斯模糊图像的正则逆扩散方程复原方法摘要：利用高斯卷积和线性扩散的等价性，从偏微分方程逆问题的角度，提出了一种针对高斯模糊图像的复原方法：RBD-PDE（Regularized Backward Diffusion） ；从频率域角度分析了逆扩散方程的正则化表达式和正则滤波之间的关系；得出正则滤波器最佳截止频率和反向扩散时间之间的关系，为以实验的方式进行盲反卷积提供便利。较传统的基于能量范涵的复原方法，如维纳滤波或 TV 模型，RBD-PDE 方法具有最佳复原效果（在高斯核标准方差已知或未知的情况下，RBD 的结果均优于传统能量泛函方法的最佳结果） 。关键词：偏微分方程，逆问题，正则化，图像复原1 引言图像复原是图像处理中的经典问题，对于线性系统，图像的模糊过程可以看作原始的清晰图像 与核函数（本文假设高斯核 ）的卷积，而图像复原或反卷积（去卷积）是0uG从模糊图像 复原清晰图像的过程，数学形式为：t(1)00tuudx=y图像反卷积包括核函数已知与核函数未知（盲反卷积）的两种情况，已有大量的研究文献提出了多种方法，如文献13,7为核函数已知情况，文献46,10,12为核函数未知的情况等。大多数方法都基于能量泛函理论，通过加入约束条件建立优化模型，如维纳滤波方法、有约束的最小二乘法、整体变分（TV）模型 7 等，或使用自然图像的统计特性取代梯度 4,6,10,11,12，以实现稳定和准确的进行求解。对于基于能量泛函的方法，准确知道核函数对于复原效果起着至关重要的作用 13。当高斯核函数的标准方差未知时，有无数组 满足式(1)，因此，需要加入对0,Gu的假设（先验知识） 。稀疏性是最常用的假设，对于主要应用于运动模糊的盲卷积能0,Gu取得较好的效果 4,5,6,10。但是当稀疏性不满足时，例如高斯核函数，传统的基于稀疏先验的方法难以取得较好的效果。不同于传统的基于能量泛函的方法，本文从偏微分方程和逆问题的角度出发，提出一种全新的针对高斯模糊图像的复原方法：RBD-PDE（Regularized Backward Heat Diffusion） 。较之于传统的基于能量泛函的方法，RBD-PDE 在高斯核标准方差未知的情况下，仍然能够有效地实现图像复原，性能优于传统方法。RBD-PDE 容易和现有的线性或非线性偏微分方程方法相结合，构成新的复原模型，因此具有更大的灵活性和方法的可拓展性。对于运动模糊图像已有许多有效的复原方法 4,5,10，而复杂的图像模糊可以分解为运动模糊和高斯模糊 ？ ，并分别进行复原。因此，高斯模糊的复原方法具有很重要的研究和实用价值。2 正向和逆扩散方程线性热扩散方程的解为高斯核（热核）函数与初始条件函数的卷积，热扩散过程等价于高斯模糊过程。自然的，从偏微分方程角度看，图像复原可视为正向热扩散的逆过程。2.1 正向扩散方程对于线性热扩散偏微分方程：（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,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）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表示无噪声精确解， 为噪声放大项：,exactvn（12 ）2,ptt服从期望为 0 方差为 的高斯分布，由式 (12)可以看出，噪声中的高频分量将被迅速放2大（以指数速度） ，覆盖真实（希望得到的复原）结果，使得实际应用中，无法直接利用式(6)进行图像复原。3 正则逆扩散方程利用逆扩散方程进行图像复原最早由 Garbor 提出。由于逆扩散方程的病态性，噪声被迅速放大，从而使得方程(6)的应用受到极大的制约。从逆问题理论出发，通过正则化方法控制噪声的传播和放大是解决式(6)的病态性的有效途径。3.1 正则化方法如果要求复原结果中的噪声放大项 满足：,nt（13 ）2expt对于任意 都成立，那么逆扩散的极限时间 需要满足：(无推导)( 随即变量不能取具体mat值, ?)2（14 ）max21lnt由式(14)可以看出， 由噪声的方差 决定，与正向扩散时间 无关。当 远小maxt2 1tmaxt于 时，无法用方程式(6) 进行有效复原。因此，增加逆扩散方程实用性的关键在于降低1t对 的限制，使 接近 ，解决的途径之一是强制使得式(12) 中较大 所对应的值为maxat1t 0，即放弃对原始图像中高频分量的恢复来换取计算过程的稳定，而这正是正则化的基本思想。采用低通滤波器去除图像中高频分量，如使用截断窗口函数， （15 ）if 0,式(15)可滤除图像中的高频分量，称为正则滤波器， 为对应的截止频率。对 进行滤波，逆扩散结果（式 11）变为：,vt（16 ）, , ,exactvvnt其中 ， 分别为：,exactnt（17 ）2, 01pexact ut（18 ）2,ex对于扩散方程，正则化等价于低通滤波过程，正则化后，逆扩散的极限时间 变为：maxt（19 ）max21lnt提高了逆扩散时间，扩大了逆扩散方程的适用范围。如同所有的正则化一样，正则化要付出一定的代价，正则滤波器 是以丢失高频信息为代价，换取逆扩散极限时间的增加，因此，其复原效果自然要受到一定的制约。下面给出 5050 的二值图像“矩形”的实验，在这里的计算采用有限差分法的显式格式，高斯白噪声的方差 ，实验结果如同图 1 所示：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) 的 DCT 系数（以对数形式显示） 。(f) 为正则化前图像的 DCT 系数（对数形式显示） ，可以看出，能量主要集中在高频（右下角） ，说明高频噪声被迅速放大。实验中采用归一化均方差：（20 ）0kEvu作为(c) 选取的标准。由于截断窗口函数式(15)不连续，滤波会产生吉布斯现象（边缘附近出现振荡，产生重影） 。选取其它函数来优化窗口函数可以解决该问题，例如，选取亚高斯窗口函数：（21 ）4exp作为低通滤波器，并对复原后的图像进行投影：（22 ）1,if,1, 0,i,vtvtttx可提高图像处理效果。图 2 是在相同实验条件下，优化后得到的实验结果：2040204000.5120 40204000.51(a) (b) (c)图 2 亚高斯窗口滤波及投影后的结果图 2 中(a)为亚高斯窗口（式(21)）正则滤波后经过式(22)进行投影后的结果；(b)为(a)的灰度值所对应的曲面；(c)为图 1(d)的灰度值所对应的曲面。对比(b)和(c) 可以看出，通过调整窗函数和投影，吉布斯现象得到了有效的抑制，图像效果有所提高。3.2 正则逆扩散框架从频率域分析的角度考虑，正则化方法实际上是对逆扩散方程的解进行低通滤波，以抑制噪声中高频分量的传播和放大，降低对逆扩散时间的限制，达到更好的复原效果。即： （23 ）21, 0expvFtFvx如果令：（24 ）100可以反解出正则复原方程，写为：（25 ）200, ,vtttRLvxx其中 表示低通滤波器 所对应的滤波算子。说明对于 LTI（线性时不变）的情况，理L论上，以 为初始条件的正则逆扩散方程等价于 【 为初始条件】的逆扩0vx0x0v散方程，即逆扩散方程的正则化等价于对扩散方程初始条件进行低通滤波。在数值计算中，由于存在机器误差，当迭代次数较大时，仍然需要对结果再进行低通滤波，以消除误差的积累。这种情况下，正则复原过程等价于：（26 ）21 0, 1expvFtFvx作为改进，可以考虑使用时变正则低通滤波器：（27 ）4,exptt抑制迭代过程中计算误差的传播和放大，由式(8)和式(17)可以求得对应的时变 RBD PDE:为什么 （28）2200, 0,vttvttRLxx其中 表示四阶偏导数算子。为了减小式(27) 中时间 变化的影响，可以采24211xx t用：为什么（29）22001, 0,vttvttRLxxRBD-PDE(式 29)所对应的时变滤波器为：取何值 （30 ）4,eplntt在 较小的条件下式(30)近似于式(21)。 什么时候停止t可以看出，RBD-PDE 是一个理论框架，可以根据本文提出的