图像退化-图像复原
4记录和整理实验报告。图像降质的数学模型
图像复原处理的关键问题在于建立退化模型。输入图像f(x, y)经过某个退化系统后输出的是一幅退化的图像。为了讨论方便， 把噪声引起的退化即噪声对图像的影响一般作为加性噪声 考虑， 这也与许多实际应用情况一致，如图像数字化时的量化 噪声、 随机噪声等就可以作为加性噪声，即使不是加性噪声而 是乘性噪声， 也可以用对数方式将其转化为相加形式。
原始图像f(x, y) 经过一个退化算子或退化系统H(x, y) 的作
用， 再和噪声n(x,y)进行叠加，形成退化后的图像g(x, y)。图2-1表示退化过程的输入和输出的关系，其中H(x, y)概括了退化系统的物理过程，就是所要寻找的退化数学模型。
图2-1 图像的退化模型
数字图像的图像恢复问题可看作是： 根据退化图像g(x , y)和退化算子H(x , y)的形式，沿着反向过程去求解原始图像f(x , y)， 或者说是逆向地寻找原始
图像的最佳近似估计。图像退化的过程可以用数学表达式写成如下的形式：
g(x, y)=H ［f(x, y)］+n(x, y) (2-1)
在这里，n(x, y)是一种统计性质的信息。在实际应用中， 往往假设噪声是白噪声，即它的频谱密度为常数，并且与图像不相关。
在图像复原处理中， 尽管非线性、 时变和空间变化的系统模型更具有普遍性和准确性，更与复杂的退化环境相接近，但它给实际处理工作带来了巨大的困难， 常常找不到解或者很难用计算机来处理。因此，在图像复原处理中， 往往用线性系统和空间不变系统模型来加以近似。这种近似的优点使得线性系统中的许多理论可直接用于解决图像复原问题，同时又不失可用性。
H (x , y )f (x , y )g (x , y )
n (x , y )
~
匀速直线运动模糊的退化模型
在所有的运动模糊中，由匀速直线运动造成图象模糊的复原问题更具有一般性和普遍意义。因为变速的、非直线运动在某些条件下可以被分解为分段匀速直线运动。本节只讨论由水平匀速直线运动而产生的运动模糊。
假设图象()y x f ,有一个平面运动，令()t x 0和()t y 0
分别为在x 和y 方向上运动的变化分量，T 表示运动的时间。记录介质的总曝光量是在快门打开后到关闭这段时间的积分。则模糊后的图象为：
[]dt t y y t x x f y x g T ?--=000)(,)(),( （2-2）
式中g(x,y)为模糊后的图象。以上就是由于目标与摄像机相对运动造成的图象模糊的连续函数模型。
如果模糊图象是由景物在x 方向上作匀速直线运动造成的，则模糊后图象任意点的值为：
()[]dt y t x x f g T y x ?-=00,,)( （2-3）
式中()t x 0是景物在x 方向上的运动分量，若图象总的位移量为a ，总的时
间为T ，则运动的速率为()t x 0=at/T 。则上式变为：
dt y T at x f y x g T ?-=???????
?0,),( （2-4） 以上讨论的是连续图象，对于离散图象来说，对上式进行离散化得：
t y T at x f y x g L i ?∑-=-=????? ?
?10,),( （2-5） -
其中L 为照片上景物移动的像素个数的整数近似值。是每个像素对模糊产生影响的时间因子。由此可知，运动模糊图象的像素值是原图象相应像素值与其时间的乘积的累加。
从物理现象上看，运动模糊图象实际上就是同一景物图象经过一系列的距
离延迟后再叠加，最终形成的图象。如果要由一幅清晰图象模拟出水平匀速运动模糊图象，可按下式进行：
∑=-=10
),(1),(L i y x f L y x g （2-6） 这样可以理解此运动模糊与时间无关，而只与运动模糊的距离有关，在这种条件下，使实验得到简化。因为对一幅实际的运动模糊图象，由于摄像机不同，很难知道其曝光时间和景物运动速度。
我们也可用卷积的方法模拟出水平方向匀速运动模糊。其过程可表示为：
),(),(),(y x h y x f y x g *= (2-7)
其中
?????-≤≤=其它0101),(L x L y x h (2-8)
h(x,y)称为模糊算子或点扩散函数，“*”表示卷积，),(y x f 表示原始(清晰)图象，),(y x g 表示观察到的退化图象。
如果考虑噪声的影响，运动模糊图象的退化模型可以描述为一个退化函数和一个加性噪声项),(y x n ，处理一幅输入图象),(y x f 产生一幅退化图象),(y x g 。
),(),(),(),(y x n y x h y x f y x g +*= (2-9)
$
由于空间域的卷积等同于频率域的乘积，所以式(2-9)的频率域描述为： ),(),(),(),(v u N v u F v u H v u G += (2-10)
式(2-9)中的大写字母项是式(2-10)中相应项的傅里叶变换。
点扩散函数的确定
不同的点扩散函数（PSF ）会产生不同的模糊图象。明确的知道退化函数是很有用的，有关它的知识越精确，则复原结果就越好。首先讨论几个典型的点扩散函数。
2.3.1典型的点扩散函数
运动模糊的点扩散函数：假设图象是通过一个具有机械快门的摄像机获得的。摄像机和拍摄物体在快门打开期间T 的相对运动引起物体在图象中的平滑。假设V 是沿x 轴方向的衡常速度，时间T 内PSF 的傅里叶变换H(u,v)由下式给出：
Vu
vtu ππ)sin(v)H(u,= (2-11) 离焦模糊的点扩散函数：由于焦距不当导致的图象模糊可以用如下函数表示：
ar J
v u ),(1v)H(u,= (2-12) 其中J 1是一阶Bessel 函数，v u 222r +=，a 是位移。该模型不具有空间不
变性。
)
大气扰动的点扩散函数：大气的扰动造成的图象模糊在遥感和天文中是需要复原的。它是由大气的不均匀性使穿过的光线偏离引起的，以下给出了数学模型，其表达式为：
()e v u c v u )(,H 2265+-= (2-13)
其中c 是一个依赖扰动类型的变量，通常通过实验来确定。幂
5/6有时用1代替。
当我们得到一幅退化图象的时候，首先要判断其退化类型然后通过已知的先验知识进行恢复。以下的讨论主要针对运动模糊PSF 进行。
2.3.2运动模糊点扩散函数的离散化
对于运动模糊而言，根据相机与目标的相对运动速度，相机的焦距以及相机相对目标的距离等就可以计算出PSF 。例如通过计算得到一幅模糊图片的模糊方向是x=6,y=4，连续的PSF 为见图(a)。
根据式(2-14)和式(2-15)，可以得到离散化以后的PSF 如图(b)。可以看到由于离散化的原因PSF 并非是直线。在本文中的PSF 均由此方法获得。
1,,1,0-==???
?????n x x n m y (2-14) 1,,1,0-==????????n y y m n y (2-15)
$
n ：y 方向的模糊点数；m ：x 方向的模糊点数。当n>=m 的时候用式(2-14)计算当n 3、 运动模糊图象的复原方法及原理 为了抑制退化而利用有关退化性质知识的预处理方法为图象复原。多数图象复原方法是基于整幅图象上的全局性卷积法。图象的退化可能有多种原因：光学透镜的残次、光电传感器的非线性、胶片材料的颗粒度、物体与摄像机间的相对运动、不当的焦距、遥感或天文中大气的扰动、照片的扫描等等。图象复原的目标是从退化图象中重构出原始图象。 运动模糊图象的恢复是图象复原的主要分支之一，它的恢复算法有很多种。有些算法虽然有很好的恢复效果，但算法复杂，恢复时间比较长（如最大熵法）。有些算法虽然计算速度较快，但恢复效果不尽人意（如空间域逆向恢复）。 下面介绍逆滤波、维纳滤波和有约束最小二乘滤波三种恢复方法的原理。 逆滤波复原原理 在六十年代中期，逆滤波（去卷积）开始被广泛地应用于数字图象复原。Nathan 用二维去卷积方法来处理由漫游者、探索者等外星探索发射得到的图象。由于和噪声相比，信号的频谱随着频率升高下降较快，因此高频部分主要是噪声。Nathan 采用的是限定逆滤波传递函数最大值的方法。 在同一时期，Harris 采用PSF 的解析模型对望远镜图象总由于大气扰动造成的模糊进行了逆滤波处理，Mcglamery 则采用由实验确定的PSF 来对大气扰动图象进行逆滤波。从此以后，逆滤波就成了模糊图象复原的一种标准技术。 恢复退化图象最简单的方法是直接逆滤波。在该方法中，用退化函数除退化 图象的傅里叶变换F(u,v)来计算原始图象的傅里叶变换估计),(?v u G ，由式3-1可以得到逆滤波退化公式： ),(),(),(),(?v u H v u N v u F v u G -= （3-1） 这个公式说明逆滤波对于没有被噪声污染的图象很有效，这里不考虑在u ，v 空间的某些位置上当H(u,v)接近0时可能遇到的计算问题，幸运的是忽略这些点在恢复结果中并不会产生可感觉到的影响。但是，如果出现噪声就会引起几个问题：第一，对于H(u,v)幅值比较小的频率处噪声的影响可能变得显著起来。这种状况通常对于高频u ，v 。在实际中，通常H(u,v)幅值衰减得比N(u,v)快得多，因此噪声的影响可能支配整个复原结果。将复原限定在H(u,v)足够大得u ，v 原点处得一个小邻域中，可以克服这个问题。第二个问题针对噪声本身的频谱，我们通常没有充分的有关噪声的信息来足够好地确定N(u,v)。为了克服H(u,v)接近0所引起的问题，在分母中加入一个小的常数k ，将式(3-1)修改为： ( k v u H v u N v u G v u F +-=),(),(),(),(? （3-2） 维纳滤波复原原理 在大部分图象中，邻近的像素是高度相关的，而距离较远的像素其相关性 较弱。由此，我们可以认为典型图象的自相关函数通常随着与原点的距离增加下降。由于图象的功率谱是图象本身自相关函数的傅里叶变换，我们可以认为的功率谱随着频域的升高而下降。 一般地，噪声源往往具有平坦的功率谱，即使不是如此，其随着频率的升而下降的趋势也要比典型图象的功率谱慢得多。因此，可以料想功率谱的低频分以信号为主，然而高频部分则主要被噪声所占据。由于逆滤波滤波器的幅值常随着频率的升高而升高，因此会增强高频部分的噪声。为了克服以上缺点，出了采用最小均方误差的方法（维纳滤波）进行模糊图象恢复。 维纳(wiener)滤波可以归于反卷积(或反转滤波)算法一类，它是由Wiener 首提出的，并应用于一维信号，并取得很好的效果。以后算法又被引入二维信号理，也取得相当满意的效果，尤其在图象复原领域，由于维纳滤波器的复原效良好，计算量较低，并且抗噪性能优良，因而在图象复原领域得到了广泛的应用并不断得到改进发展，许多高效的复原算法都是以此为基础形成的。 如果取 R R Q n g 2121 = (3-3) R g 和R n 分别是图象和噪声的自相关矩阵。即??????=g g E R T g ，? ?????=n n E R T n ，并且都是正定对称矩阵，则有 f H R R H H g T n g T 1)(?1-+=-γ (3-4) g R R g Q n g ??2121=的模方最小，实际上就意味着使噪声和信号的比对复原图象影响最小。因为图象和噪声的相关矩阵都是把图象当作随机过程来研究，从而描述其统计特性的量，在这里最小二乘方的最佳已经演变成均方误差最小准则下的最佳。 同样根据式(3-4)可求得频域维纳滤波公式如下 \ ),(),() ,(),(),(),(1),(?22v u F v u S v u S v u H v u H v u H v u G g n γ+= (3-5) γ=1时，为标准维纳滤波器；γ≠1时，为含参维纳滤波器。若没有噪声时 即),(n v u S =0,维纳滤波器则退化成理想反滤波器。 实际应用中必须调节γ以满足式(3-4)。因为),(n v u S ,),(g v u S 实际很难求得因此，可以用一个比值k 代替两者之比，从而得到简化的维纳滤波公式 ),(),(),(),(1),(?22 v u F k v u H v u H v u H v u G += (3-6) 有约束最小二乘复原原理 由于大多数图象恢复问题都不具有唯一解，或者说恢复具有病态特征。为了克服这一问题，通常需要在恢复过程中对运算施加某种约束。 设对图象施加某一线性运算Q ，求在约束条件 n g H f 22?=- （3-7） 下，使g Q ?2为最小的g ?作为原图g 的最佳估计。 利用拉格朗日乘数法，先构造一辅助函数： | )?(?),?(222n g H f g Q g j ---=λλ （3-8） 令0?),?(=??g g j λ可得： 0)?(2)?(2=--g H f H g Q Q T T λ （3-9） 解之得： f H Q Q H H g T T T 1)(?-+=γ （3-10） 式中λγ1=。把式(3-10)代入式(3-7)中可以证明，n 2是γ的单调递增函 数。因此可以用迭代法求出满足约束条件(3-7)式的待定系数γ，首先任取一个γ, 代入(3-10)，把求得的g ?再代入式(3-7)，若结果大于n 2时，便减少γ；反之增大γ，再重复上述过程，直到约束条件式(3-11)被满足为止（实际求解时，只要能n g H f 2?2--之差小于某一给定值就可以了）。把求得的γ代入，便最后求得 最佳估计g ?。 我们可以直接从空间域的有约束最小二乘方恢复式(3-10)得到它的频域解 ),(),(),(),(),(1),(?222 v u F v u C v u H v u H v u H v u G γ+= (3-11) 应用有约束最小二乘方恢复方法时，只需有关噪声均值和方差的知识就可对每幅给定的图象给出最佳恢复结果。 4、运动模糊图像复原的实现 数字图像在获取的过程中，由于光学系统的像差、光学成像衍射、成像系统的非线性畸变、摄影胶片的感光的非线性、成像过程的相对运动、大气的湍流效应、环境随机噪声等原因，图像会产生一定程度的退化.因此，必须采取一定的方法尽可能地减少或消除图像质量的下降，恢复图像的本来面目，这就是图像复原，也称为图像恢复。 图像复原是试图利用退化过程的先验知识使已退化的图像恢复本来面目，即根据退化的原因，分析引起退化的环境因素，建立相应的数学模型，并沿着使图像降质的逆过程恢复图像。前面讲述的三种运动模糊图像复原方法均可在matlab 环境下实现，下面将对每一种方法进行举例。