数字图像处理第5章图像复原北邮出版社200810

1、12v5.1 5.1 图像降质的数学模型图像降质的数学模型v5.2 5.2 无约束图像复原无约束图像复原v5.3 5.3 有约束图像复原有约束图像复原3高斯噪声引起模糊高斯噪声引起模糊4运动引起的模糊运动引起的模糊光学聚焦引起的模糊光学聚焦引起的模糊5 图像复原的关键：建立退化模型。图像复原的关键：建立退化模型。 常用图像退化和复原的模型：仅限于讨论线性、空不变系统常用图像退化和复原的模型：仅限于讨论线性、空不变系统 f(x,y)是原始图像，是原始图像，h(x,y) 是退化函数，是退化函数， n(x,y) 是加性噪声，是加性噪声，(x,y) 是复原是复原函数，函数， 是复原后的图像是复原后的图

2、像。 退化模型的关键：对退化系统退化模型的关键：对退化系统 h(x,y) 的了解。的了解。图图5.1 5.1 图像的退化及复原模型图像的退化及复原模型),(yxf6v5.1.1 5.1.1 连续图像退化的数学模型连续图像退化的数学模型 一幅连续图像一幅连续图像 f(x,y) 可以通过点源函数的卷积来表示可以通过点源函数的卷积来表示 对于线性空间不变系统，输入图像经退化后为对于线性空间不变系统，输入图像经退化后为 其中其中h(x-,y-)称为该退化系统的点扩展函数，或叫系统的冲激响应。称为该退化系统的点扩展函数，或叫系统的冲激响应。( , )( ,) (,)f x yfxyd d ( , ) (

3、 , )( ,) (,)( ,)(,)( ,) (,)g x yH f x yHfxyd dfHxyd dfh xyd d （5.15.1）（5.65.6）7 加性白噪声加性白噪声n(x,y)引起引起图像退化：图像退化： 频域表达：频域表达： 图像复原的实质：已知图像复原的实质：已知g(x,y)、h(x,y)，求求f(x,y)， 或或 已知已知G(u,v)、H(u,v)，求求F(u,v)。 关键是降质模型：时域的冲激响应函数关键是降质模型：时域的冲激响应函数h(x,y) ， 或频域的传递函数或频域的传递函数H(u,v) 。( , )( ,) (,)( , )( , )( , )( , )g x

4、 yfh xyd dn x yf x yh x yn x y ( , )( , )( , )( , )G u vF u v H u vN u v（5.85.8）（5.95.9）8v5.1.2 5.1.2 几个典型的退化模型几个典型的退化模型 1. 1. 受到孔径衍射造成的图像退化受到孔径衍射造成的图像退化l 光学孔径衍射效应是造成图像模糊的主要原因之一。光学孔径衍射效应是造成图像模糊的主要原因之一。l 相干光的光学成像系统的递函数由下式给出：相干光的光学成像系统的递函数由下式给出：l 圆形孔径的点扩展函数为圆形孔径的点扩展函数为图图5.2 5.2 小孔衍射造成的模糊小孔衍射造成的模糊22( ,

5、 )(,)H u vpd ud v( , ) ( , )h x yF p 92. 2. 目标相对运动造成的图像模糊目标相对运动造成的图像模糊l 在获取图像时，由于景物和摄像机之间的相对运动，造成图像的模糊。在获取图像时，由于景物和摄像机之间的相对运动，造成图像的模糊。l 设物体设物体f(x,y)在一平面作在一平面作匀速直线匀速直线运动，令运动，令x(t)和和y(t)分别是物体在分别是物体在x、y方向方向上的分量，上的分量，t 表示运动的时间，表示运动的时间，l 记录介质的总曝光量在快门打开到关闭这曝光瞬间内的积分，记录介质的总曝光量在快门打开到关闭这曝光瞬间内的积分，l 设设T为曝光时间，则曝

6、光成像后的降质图像为：为曝光时间，则曝光成像后的降质图像为： 图图5.5 5.5 目标相对运动造成的图像模糊目标相对运动造成的图像模糊 dttyytxxfyxgT000,（5.175.17）10l 对式（对式（5.175.17）两边进行傅立叶变换，得）两边进行傅立叶变换，得 dxdyvyuxjdttyytxxfdxdyvyuxjyxgvuGT)(2exp,)(2exp,000dtdxdyvyuxjtyytxxfvuGT 000)(2exp)(),(),(dttvytuxjvuFdttvytuxjvuFvuGTT000000)()(2exp,)()(2exp,dttvytuxjvuHT000)(

7、)(2exp,交换积分次序FT位移性质（5.185.18）（5.195.19）（5.215.21）),(),(,vuHvuFvuH（5.205.20）11l 如果只有如果只有x方向的匀速运动，即方向的匀速运动，即y0(t)=0，则模糊后图像，则模糊后图像g(x,y)可以简化为：可以简化为：l 如果在如果在T时间里物体运动的总位移量为时间里物体运动的总位移量为a，则在任意，则在任意t时间里物体在时间里物体在x方向上方向上的分量：的分量： 于是，于是，模糊系统的传递函数：模糊系统的传递函数：0002( , )exp2( )exp2(1)2TTjuaH u vjuxtdtatjudtTjTeuadt

8、ytxxfyxgT00),(),(（5.235.23）Tattx)(0（5.245.24）12 3. 3. 大气湍流造成的图像降质大气湍流造成的图像降质l 航空、卫星、天文图像中，由于受大气湍流的影响，使图像产生退化。航空、卫星、天文图像中，由于受大气湍流的影响，使图像产生退化。l 长时间作用的情况下，大气湍流降质图像的系统传递函数为：长时间作用的情况下，大气湍流降质图像的系统传递函数为： 式中式中C为与湍流性质有关的常数。为与湍流性质有关的常数。22( , )exp()H u vC uv（5.265.26）13v5.1.3 5.1.3 离散图像退化的数学模型离散图像退化的数学模型 1. 1.

9、 一维离散情况退化模型一维离散情况退化模型 设：离散输入函数设：离散输入函数f(x)，有有A个采样，退化系统的冲激响应个采样，退化系统的冲激响应h(x) ，有有B个个采样，则：经退化系统后的离散输出函数采样，则：经退化系统后的离散输出函数g(x)： 分别对分别对f(x)和和h(x)用添零延伸的方法扩展成周期为用添零延伸的方法扩展成周期为M M= =A A+ +B B-1-1的周期函数：的周期函数： 此时输出为：此时输出为： ( )( )( )g xf xh x（5.275.27） 1010MxAAxxfxfe 1010MxBBxxhxhe10( )( )( )( )()Meeeeemgxfxh

10、 xf m h xm（5.285.28）（5.295.29）14矩阵表示矩阵表示因为因为he(x)的周期为的周期为M，所以，所以he(x)=he(x+M) ，MM 阶矩阵阶矩阵H为：为：H是循环矩阵：每一行元素相同，每行以循环方式右移一位，是循环矩阵：每一行元素相同，每行以循环方式右移一位，循环矩阵相加还是循环矩阵，循环矩阵相乘还是循环矩阵。循环矩阵相加还是循环矩阵，循环矩阵相乘还是循环矩阵。 （5.315.31）（5.305.30）) 1()2() 1 ()0()0()2() 1()3() 1 ()2()2()0() 1 () 1() 1()0() 1()2() 1 ()0(MffffhMh

11、MhMhhhMhhhMhhhMggggeeeeeeeeeeeeeeeefhg)0()2()1()3()1 ()2()2()0()1 ()1()1()0(eeeeeeeeeeeehMhMhMhhhMhhhMhhhH15 2. 2. 二维离散模型二维离散模型推广到二维，设输入的数字图像推广到二维，设输入的数字图像f(x,y)大小为大小为AB，点扩展函数点扩展函数h(x,y)被均匀采样为被均匀采样为CD。用添零的方法将它们扩展成用添零的方法将它们扩展成M=A+C-1和和N=B+D-1的周期函数的周期函数。 输出降质数字图像输出降质数字图像其它且01010,ByAxyxfyxfe其它且01010,Dy

12、Cxyxhyxhe（5.335.33）（5.345.34） 1100,MNeeeeemngx yfm n hxm ynfx yhx y 16矩阵形式为：矩阵形式为： 式中式中g、f是是MN1维列向量，维列向量，H是是MNMN维矩阵。维矩阵。将将f(x,y)中的元素排成的列向量：中的元素排成的列向量： 将将g(x,y)中的元素排成的列向量：中的元素排成的列向量：fhg（5.355.35）（5.365.36）TM 1 (0,0), (0,1),(0,1), (1,0), (1,1),(1,1),(1,0)(1,1)fffNfffNf Mf MNf 第一行元素第二行元素第 行元素， ，T1 (0,0

13、), (0,1), (0,1), (1,0), (1,1), (1,1), (1,0), (1,1)MgggNgggNg Mg MNg 第一行元素第二行元素第行元素，17Hj (j=0,1,M-1)为子矩阵，大小为为子矩阵，大小为NN，H矩阵是由矩阵是由MM个大小为个大小为NN的子矩阵组成分块循环矩阵，的子矩阵组成分块循环矩阵，每个子矩阵每个子矩阵Hj是由延拓函数的第是由延拓函数的第 j 行构成的，构成方法如下：行构成的，构成方法如下：032121011210HHHHHHHHHHHHHMMMMMM( ,0)( ,1)( ,2)( ,1)( ,1)( ,0)( ,1)( ,2)( ,1)( ,2

14、)( ,3)( ,0)eeeeeeeejeeeehjhj Nhj Nhjhjhjhj Nhjhj Nhj Nhj NhjH （5.385.38）（5.375.37）18加上噪声影响，一个更加完整的离散图像退化模型加上噪声影响，一个更加完整的离散图像退化模型写成矩阵形式写成矩阵形式 1100( ,)(,)(,) +( ,)MNeeeemngx yfm n hxm ynnx y g = H f + n （5.395.39）（5.405.40）19v5.1.4 5.1.4 循环矩阵的对角化循环矩阵的对角化 1. 1. 循环矩阵的对角化循环矩阵的对角化对于循环矩阵对于循环矩阵H，标量函数，标量函数(k

15、)和和w(k)分别是它的特征值和特征向量分别是它的特征值和特征向量，由矩阵乘法可得由矩阵乘法可得( )( ) ( )0,1,1kkkkMH ww222( )(0)(1)exp(2)exp2(1)exp1eeeekhh Mjkh MjkhjMkMMM222( )1expexp2exp1TkjkjkjMkMMMw（5.425.42）（5.415.41）（5.435.43）20H的的M个特征向量组成一个个特征向量组成一个MM的矩阵的矩阵W ：把把H写成：写成： D是一个对角矩阵，其对角元素正是是一个对角矩阵，其对角元素正是H的特征值，的特征值，D(k,k)= (k)， 1H = WDW（5.455.

16、45）(0)000(1)000000(1)MD（5.465.46）T (0),(1),(1)MWwww（5.445.44）21 2. 2. 分块循环矩阵的对角化分块循环矩阵的对角化 推广到二维离散退化模型，可以使得块循环矩阵对角化。推广到二维离散退化模型，可以使得块循环矩阵对角化。 定义一个定义一个MNMN（包含（包含MM个个NN的块）的矩阵的块）的矩阵W ，W的第的第i行第行第m列列个子块为：个子块为： 其其中子块中子块WN为为NN矩阵，第矩阵，第k行第行第n列位置列位置的元素为：的元素为：2( ,)expNi mjimMWW 2( , )expNwk njknN （5.485.48）（5.

17、495.49）22 由上可知，由上可知，矩阵矩阵W是由分块循环矩阵是由分块循环矩阵H的特征向量组成，必定有逆矩阵的特征向量组成，必定有逆矩阵W-1，逆矩阵，逆矩阵W-1的的MN个特征向量是线性无关的，个特征向量是线性无关的，逆矩阵逆矩阵W-1的与的与W相似，是分块循环矩阵，第相似，是分块循环矩阵，第i行第行第m列个子块为：列个子块为：子块子块WN-1仍为仍为NN的矩阵，的矩阵，第第k行行n列元素为列元素为：-1-1H = WDWD = W HW或1112,expNi mjimMMWW112,expNwk njknNN（5.505.50） （5.535.53）(5.52)(5.52)（5.515

18、.51）分块循环矩阵分块循环矩阵H可以写成可以写成：23上式中矩阵上式中矩阵D是一个是一个MNMN维的矩阵，它是由维的矩阵，它是由MM个对角子矩阵组成的个对角子矩阵组成的分块循环矩阵，它的分块循环矩阵，它的MN个对角元素是个对角元素是H的特征值。即：的特征值。即：子矩阵子矩阵Dk是是NN的对角阵，对角元素是对应的特征值的对角阵，对角元素是对应的特征值：实际上，实际上，H的转置矩阵的转置矩阵HT 可用可用D的复共扼的复共扼D* 来表示，即来表示，即0M 1M 2110M 12M 1M 2M 30DDDDDDDDDDDDD T1*H = WD W(k)0000(k)00000(k)kD(5.54)

19、(5.54)24v5.1.5 5.1.5 对角化在降质模型中的应用对角化在降质模型中的应用先考虑一维情况。把式先考虑一维情况。把式(5.42)(5.42)中的中的 代入到式代入到式(5.29)(5.29)可得：可得：用用 左乘上式两边，得左乘上式两边，得1H =WDW-1g = WDW f1W-1-1-1-1W g = W WDW f = DW f111112221expexp2exp(1)12221exp2exp4exp(1)22221 exp(1) exp(1)2exp(1)(1)jjjMMMMjjjMMMMMjMjMjMMMMMW(5.55)(5.55)(5.56)(5.56)(5.57

20、)(5.57)25g是是M维的列向量，乘积维的列向量，乘积W-1g也是一个也是一个M维列向量，其第维列向量，其第k项为项为G(k)， 同理，同理，W-1 f的第的第k项记为项记为F（k）：它们分别是扩展序列它们分别是扩展序列ge(x)和和fe(x)的傅立叶变换，的傅立叶变换，D矩阵的主对角线元素是矩阵的主对角线元素是H矩阵的特征值矩阵的特征值(k)，将将(k)记为记为MH(k)，根据式，根据式（5.415.41）有：有： 1012( )( )expMeiG kg ijikMM1012( )( )expMeiF kf ijikMM 102( )expMeikh ijikMH kM(5.58)(5

21、.58)(5.59)(5.59)(5.60)(5.60)26其中其中H(k)便是扩展序列便是扩展序列he(x)的傅立叶变换。的傅立叶变换。综合上述分析，可将式综合上述分析，可将式（5.515.51）简化成一维傅立叶变换序列的对应项之积，简化成一维傅立叶变换序列的对应项之积，即：将上式简化成一维傅立叶变换序列的对应项之积即：将上式简化成一维傅立叶变换序列的对应项之积推广到二维退化情况。考虑到噪声项，有推广到二维退化情况。考虑到噪声项，有g=H f+n，用，用W-1 左乘得左乘得 W-1 g可用下式表示：可用下式表示：( )( ) ( )0,1,2,1G kMH k F kkM-1-1-1-1-1

22、-1-1-1-1W g = W (Hf +n) = W Hf +W n= W WDW f +W n= DW f +W n(5.61)(5.61)(5.62)(5.62)(5.63)(5.63)TT0,00,11,01,11,01,1GGNGGNG MG MN -1W gW-1 27对于任一对变量对于任一对变量（u,v），G(u,v)可表示成：可表示成： 此式即为此式即为 的傅立叶变换。的傅立叶变换。F(u,v)、N(u,v) 和和H(u,v)对应于对应于 的傅立叶变换。的傅立叶变换。MNMN维的对角矩阵维的对角矩阵D，其，其MN个对角元素可用下式表示：个对角元素可用下式表示：综合利用上述关系，

23、可得综合利用上述关系，可得可见：退化模型的庞大方程组可简化为计算大小为可见：退化模型的庞大方程组可简化为计算大小为MN的离散傅立叶变换，的离散傅立叶变换，如用如用FFT算法，可方便地实现。算法，可方便地实现。11001( , )( , )exp2MNexyuxvyG u vgx yjMNMN( , )egx y( , ),( , ),( , )eeefx y n x y h x yD = MNH( , )( , ) ( , )( , )G u vH u v F u vN u v(5.64)(5.64)28v5.2.1 5.2.1 逆滤波逆滤波由退化模型可知：由退化模型可知： 最优化：寻找一个最

24、优化：寻找一个 ，使得，使得g与与 之偏差之偏差 在最小二乘方意义上最在最小二乘方意义上最小，即使得小，即使得 的范数的平方最小，这等效于噪声项的范数最小。的范数的平方最小，这等效于噪声项的范数最小。 准则函数用准则函数用 表示：表示： 根据范数的定义可知：根据范数的定义可知： 为使得准则函数为使得准则函数 最小，实际上就是求它的极小值。最小，实际上就是求它的极小值。 在求极小值时，不受任何其它条件的约束，因此也称为无约束复原。在求极小值时，不受任何其它条件的约束，因此也称为无约束复原。n=g-HffHfg - Hf(fJ)22(fng-HfJ)=22TTn= n n,g-Hf= (g-Hf)

25、 (g-Hf)(fJ)(5.65)(5.65)g-Hf(5.66)(5.66)(5.67)(5.67)29根据向量微分的性质：根据向量微分的性质： 内积内积TX对对X的偏导的偏导 内积内积XT对对X的偏导的偏导 二次型二次型XTAX对对X的偏导的偏导可得：可得：TXXTXX2TX AXAXX(5.70)(5.70)(5.69)(5.69)(5.68)(5.68) 20TTTTTTTJ fg gg Hff H gf H HfHgHfff(5.71)(5.71)30 即：即： 可以解出：可以解出： 两边乘上两边乘上 ，可得：可得： 上二式表明在最小二乘方准则下寻找出的最优估计图像上二式表明在最小二

26、乘方准则下寻找出的最优估计图像 ，可由降质图像，可由降质图像g和和降质系统的冲激响应的逆矩阵降质系统的冲激响应的逆矩阵 或者其对角化矩阵或者其对角化矩阵D得出。得出。T-1T-1-1-1-1-1f = (H H) H g = H g = (WDW ) g = WD W g(5.72)(5.72)-1DW-1-1DW f = W gf-1HTTH Hf = H g 31 由前可知，由前可知，W-1g的各元素可表示成的各元素可表示成G(u,v)， 的各元素可表示成的各元素可表示成 ，而，而 D(u,v) = MNH(u,v)， 于是可以得到对应的元素之间的关系为：于是可以得到对应的元素之间的关系为

27、： MNH(u,v)是滤波函数，则它与是滤波函数，则它与 的乘积是退化图像的乘积是退化图像g(x,y)的傅立叶变换。的傅立叶变换。 复原后的图像复原后的图像 在考虑噪声的情况下，上式可写成：在考虑噪声的情况下，上式可写成：11( , )( , ) ( , )( , )N u vf x ysF u vsMNH u vfW1(,)F u v( , )( , ) ( , )G u vMNH u v F u v11( , )( , )( , )( , )G u vf x yF u vMNH u v ( , )F u v(5.74)(5.74)(5.76)(5.76)(5.77)(5.77)-1-1DW

28、 f = W g32v5.2.2 5.2.2 无约束图像复原的病态性质无约束图像复原的病态性质 病态性质：逆滤波函数中分母病态性质：逆滤波函数中分母H(u,v)有可能有可能出现零点，就会导致不定解。出现零点，就会导致不定解。l 即使没有噪声，一般也不可能精确地复原即使没有噪声，一般也不可能精确地复原f(x,y)；l 如果考虑噪声项如果考虑噪声项N(u,v) ，则出现零点时，噪声项将被放大，零点的影响将会，则出现零点时，噪声项将被放大，零点的影响将会更大，对复原的结果起主导地位。更大，对复原的结果起主导地位。 解决办法：解决办法：l （1 1）采用的有约束图像复原（下一节）；）采用的有约束图像复

29、原（下一节）；l （2 2）利用利用噪声一般在高频范围且衰减速度较慢，而信号一般在低频且衰减较噪声一般在高频范围且衰减速度较慢，而信号一般在低频且衰减较快的性质快的性质，在复原时，只限制在频谱坐标离原点不太远的有限区域内运行；，在复原时，只限制在频谱坐标离原点不太远的有限区域内运行；33l （3 3）为了避免为了避免H(u,v)的值太小，一种改进方法是在的值太小，一种改进方法是在H(u,v)=0的那些频谱点及其附近，的那些频谱点及其附近，人为地设置人为地设置 的值；（见后图）的值；（见后图）l （4 4）考虑）考虑到退化系统的传递函数到退化系统的传递函数H(u,v)的带宽比噪声的带宽要窄得多，

30、其频率特性的带宽比噪声的带宽要窄得多，其频率特性具有低通性质，取恢复转移函数具有低通性质，取恢复转移函数M(u,v)为：为：l （5 5）其中）其中 的选取原则是将的选取原则是将H(u,v)为零的点除去。这种方法的缺点是复原后图像为零的点除去。这种方法的缺点是复原后图像的振铃效果较明显。的振铃效果较明显。1( , )Hu v222022201( , )( , )1uvH u vM u vuv0(5.78)(5.78)34 例：人为地设置例：人为地设置 的零点（极点）值的零点（极点）值(a) 退化系统的传递函数H(u,v)频率(b) 逆滤波器传递函数(c) 改进的逆滤波器传递函数图图 5.7 5

31、.7 逆滤波器零点的影响及其改进逆滤波器零点的影响及其改进H-1(u,v)频率M(u,v)频率1( , )Hu v35v5.2.3 5.2.3 匀速直线运动引起的图像模糊的复原匀速直线运动引起的图像模糊的复原 假设：对平面匀速运动的景物拍摄一幅图像照片假设：对平面匀速运动的景物拍摄一幅图像照片f(x,y)， x0(t) 和和y0(t) 分别是景物在分别是景物在x和和y方向的运动分量，方向的运动分量，T为曝光时间，为曝光时间， 忽略其它因素，则实际所采集的模糊图像忽略其它因素，则实际所采集的模糊图像g(x,y) ： 如只在如只在x方向上运动，则该退化系统的传递函数如下：方向上运动，则该退化系统的

32、传递函数如下： 问题：当问题：当n为整数时，为整数时，H(u,v)=0出现零点，出现零点，u=n/a 解决：避免零点，在解决：避免零点，在 1/u 范围内进行复原运算，对复原结果影响不大。范围内进行复原运算，对复原结果影响不大。 000( ,)( ),( )Tg x yf xxtyytdt2( , )(1)sin ()2j uaj uajTH u veTecuaua(5.80)(5.80)36用递推的方法复原图像用递推的方法复原图像f(x,y) ：在只有在只有x方向上（方向上（0 xL ）的运动：）的运动：令令 代入上式，（代入上式，（0 xL），得：），得：为简单起见，忽略前面的系数，并对为

33、简单起见，忽略前面的系数，并对x求导数，得到：求导数，得到：设设L=Ka，K为正整数。设为正整数。设m为为L/a的整数部分，其值为的整数部分，其值为0，1，2，K-1，余，余数用数用z表示，用表示，用a表示图像内景物移动的总距离。则：表示图像内景物移动的总距离。则： x=z+ma 当当x=L时，时，m=K-1，z = a。0( )Tatg xfxdtTatxT( )( )xxaTg xfda( )( )()f xgxf xa(5.83)(5.83)(5.82)(5.82)(5.81)(5.81)(5.84)(5.84)37 把式（把式（5.84）代入式（）代入式（5.83）得：）得： 设设(z

34、)=f(z-a) 0z(z)代表了在曝光期间景物移入代表了在曝光期间景物移入0za的部分。此时，的部分。此时，式（式（5.85）可通过）可通过(z) 用递推的方式表示为另一种形式。用递推的方式表示为另一种形式。 当当m=0，有，有 当当m=1，有，有 当当m=2，有，有 1fzmagzmafzma(5.85)(5.85) f zgzf zagzzfzagzafz22fzagzafza38 以此类推，可以得到：以此类推，可以得到： 即：即： 其中其中 为正整数，为正整数，g(x)是已知的，要求是已知的，要求f(x)，就只需要估计，就只需要估计(x)。0( )()()mkf xg xkaxma 0

35、mjfzmagzjaz(5.86)(5.86)modxma 39 从模糊图像直接估计从模糊图像直接估计(x)的方法的方法: 当当x从从0变化到变化到L时，时，m从从0变化到变化到K-1， 的自变量为的自变量为x-ma，该变量总是在，该变量总是在0到到a之间变化，即之间变化，即0 x-maa， 因此，在因此，在0 xL 区间内计算区间内计算f(x)。 由于由于重复了重复了k次，定义次，定义 代入到式（代入到式（5.86），可写成），可写成 如果对每个如果对每个kax(k+1)a进行计算，并把进行计算，并把k=0,1,2，,k-1的结果加起来，得到的结果加起来，得到0()()mjfxgxja xm

36、afxfx(5.88)(5.88) 110011KKkkxf xkaf xkaKK(5.87)(5.87)40 上式右边第一项为未知项，但当上式右边第一项为未知项，但当K很大时接近很大时接近f(x)的平均值，设它为常数的平均值，设它为常数A， 将式（将式（5.87）代入到上式，得到在）代入到上式，得到在0 xL 区间内区间内 将上式将上式重新代入（重新代入（5.86）,得出只有得出只有x方向匀速直线运动模糊图像复原后图像表方向匀速直线运动模糊图像复原后图像表达式达式 0 xL 把变量把变量y加入，即可得到最终复原图像加入，即可得到最终复原图像 0 xL1010011KkKkkjxmaAfxka

37、maKAgxkamajaK 10001KkmkjjfxAgxmakj agxjaK10001,Kkmkjjf x yAgxmakj a ygxja yK(5.91)(5.91)(5.90)(5.90)(5.89)(5.89)41 无约束图像复原：无约束图像复原：（如逆滤波）满足最优准则（如逆滤波）满足最优准则 ， 只要了解退化系统的传递函数，无其它约束；只要了解退化系统的传递函数，无其它约束； 有约束图像复原有约束图像复原： 还需知道某些噪声的统计特性，或噪声与图像的某些相关情况。还需知道某些噪声的统计特性，或噪声与图像的某些相关情况。 根据不同的先验知识，采用不同的约束条件，得到不同的复原技

38、术。根据不同的先验知识，采用不同的约束条件，得到不同的复原技术。 最常见的是有约束的最常见的是有约束的最小二乘方最小二乘方图像复原技术，图像复原技术， 此外还有：功率谱均衡复原、几何均值滤波器复原等方法。此外还有：功率谱均衡复原、几何均值滤波器复原等方法。22(fng-HfJ)=42 有约束最小二乘方复原问题：有约束最小二乘方复原问题：l 令令Q为为 的线性算子，要设法寻找一个最优估计的线性算子，要设法寻找一个最优估计 ，使形式为，使形式为 、 服从约束条件服从约束条件 的准则函数最小化。的准则函数最小化。 l 采用拉格朗日乘子算法，使得准则函数最小：采用拉格朗日乘子算法，使得准则函数最小：l

39、 对上式求导并使其为零，求解得到：对上式求导并使其为零，求解得到：l 选择不同的选择不同的Q，就可得到不同类型的有约束的最小二乘方图像复原方法。，就可得到不同类型的有约束的最小二乘方图像复原方法。 222JfQfg - HfngHQQHHfT1TTf2fQ22nfHgf(5.92)(5.92)(5.93)(5.93)43当当QI（MNMN单位阵），即单位阵），即Q不起作用，这时：不起作用，这时：如如=0，上式退化为无约束情况的图像复原。，上式退化为无约束情况的图像复原。注意：是一个在空间域上进行图像复原的理论方法。注意：是一个在空间域上进行图像复原的理论方法。由于矩阵或向量的维数巨大，实际上很

40、难直接用上述方法求解。由于矩阵或向量的维数巨大，实际上很难直接用上述方法求解。通常都需要转换到频率域上进行处理，类似于逆滤波方法，通常都需要转换到频率域上进行处理，类似于逆滤波方法，在有约束条件下称为维纳滤波复原方法。在有约束条件下称为维纳滤波复原方法。 r-1TTf = H H + IH g(5.94)(5.94)44v5.3.1 5.3.1 维纳滤波维纳滤波 条件：图像为平稳随机过程，独立于噪声；噪声为加性噪声，是平稳过程。条件：图像为平稳随机过程，独立于噪声；噪声为加性噪声，是平稳过程。 准则：准则：使使f(x,y) 和和 之间的均方误差达到最小：之间的均方误差达到最小： 维纳滤波选用的

41、变换矩阵：维纳滤波选用的变换矩阵： 其中其中Rf图像图像f的相关矩阵，的相关矩阵，Rn是噪声是噪声n的相关矩阵，且是实对称矩阵：的相关矩阵，且是实对称矩阵：T-1fnQ Q = R REETTfnR =ffR=nn( , )f x y22min( , )( , )eEf x yf x y(5.97)(5.97)(5.96)(5.96)(5.95)(5.95)(5.98)(5.98)45l 像素点越近相关度越高，自相关函数随着与远点距离增加而下降；像素点越近相关度越高，自相关函数随着与远点距离增加而下降；l 功率谱是其自相关函数的傅立叶变换，图像功率谱随着频率的升高而下降；功率谱是其自相关函数的

42、傅立叶变换，图像功率谱随着频率的升高而下降；l 图像的平稳性：像素的相关性只是相互距离的函数，和位置无关；图像的平稳性：像素的相关性只是相互距离的函数，和位置无关；l 可将可将Rf 和和Rn 都用块循环矩阵表示，利用循环矩阵的对角化，可以写成：都用块循环矩阵表示，利用循环矩阵的对角化，可以写成： Rf = W A W-1Rn = W B W-1 l 其中其中A和和B为对角矩阵，为对角矩阵，根据循环矩阵对角化的性质可知，根据循环矩阵对角化的性质可知，A和和B中的诸元素中的诸元素分别为分别为Rf 和和Rn 中诸元素的傅立叶变换，并用中诸元素的傅立叶变换，并用Sf(u,v)和和Sn(u,v)表示。表

43、示。l 由由（5.895.89）式可得：式可得： 上式两边乘以上式两边乘以W-1，得到：，得到： gWWDBWWADWWDf*11111gWDBADDfW*1111(5.99)(5.99)(5.100)(5.100)(5.101)(5.101)(5.102)(5.102)46对应的频域表达式：对应的频域表达式：(1)(1) 如果如果 ，方括号内的项被称为维纳滤波器；，方括号内的项被称为维纳滤波器；(2)(2)无噪声时，无噪声时， ，上式退化成逆滤波器；，上式退化成逆滤波器；(3)(3)如果不知道噪声的统计性质，即当如果不知道噪声的统计性质，即当 和和 未知时，上式可未知时，上式可近似（近似（K

44、为噪声对信号频谱密度之比）：为噪声对信号频谱密度之比）：222( , )( , )( , )( , )( , )/( , )( , )1( , )( , )( , )( , )/( , )nfnfHu vF u vG u vH u vS u vSu vH u vG u vH u vH u vS u vSu v1( , )0nS u v ( , )nS u v( , )fSu v2( , )( , )( , )( , )Hu vF u vG u vH u vK (5.103)(5.103)(5.104)(5.104)47v5.3.2 5.3.2 功率谱均衡复原功率谱均衡复原 另设准则函数，如使估

45、计图像另设准则函数，如使估计图像 的功率谱的功率谱 与原图像与原图像f(x,y)功功率谱率谱 相等来导出功率谱均衡复原滤波器。相等来导出功率谱均衡复原滤波器。 复原滤波器复原滤波器M(u,v)、估计图像的功率谱、估计图像的功率谱 、退化图像的功率谱、退化图像的功率谱Sg(u,v)之间关系之间关系： 根据图像退化的降质模型和功率谱的定义可知：根据图像退化的降质模型和功率谱的定义可知： ( , )( , )ffSu vSu v( , )f x y( , )fSu v( , )fSu v( , )fSu v2( , )( , )( , )gfSu vM u vSu v*2,gfnSu vE G u

46、v Gu vEH u v F u vN u vH u v F u vN u vH u vSu vSu v(5.107)(5.107)(5.106)(5.106)(5.105)(5.105)48把上式代入到式把上式代入到式（5.1065.106），），可得可得根据功率均衡滤波器的要求根据功率均衡滤波器的要求功率谱均衡复原图像的频域表达式可以写成：功率谱均衡复原图像的频域表达式可以写成：1221( , )( , )( , )( , )/( , )nfF u vG u vH u vS u vSu v ( , )( , )ffSu vSu v(5.110)(5.110)(5.109)(5.109)(5

47、.108)(5.108)1221( , )( , )( , )( , )nfM u vH u vSu vSu v11222,ffgfnSu vSu vM u vSu vH u vSu vSu v49 维纳滤波器和功率均衡滤波器比较：维纳滤波器和功率均衡滤波器比较：两者基本相似，分子相差一项两者基本相似，分子相差一项H*(u,v)；都要预先知道功率谱都要预先知道功率谱Sf(u,v)和和Sn(u,v)；无噪声时，两者都简化为逆滤波器；无信号时，两者都完全截止；无噪声时，两者都简化为逆滤波器；无信号时，两者都完全截止； 不同在不同在H(u,v)0处，维纳滤波器的恢复转移函数强迫响应为零，处，维纳滤波

48、器的恢复转移函数强迫响应为零， 功率谱均衡滤波器的恢复转移函数则不等于零；功率谱均衡滤波器的恢复转移函数则不等于零；功率谱均衡滤波器复原能力强，在某些情况下其性能优于维纳滤波器。功率谱均衡滤波器复原能力强，在某些情况下其性能优于维纳滤波器。图图5.8 5.8 维纳滤波器和功率谱均衡滤波器比较维纳滤波器和功率谱均衡滤波器比较50v5.3.3 5.3.3 有约束最小平方复原有约束最小平方复原 维纳滤波维纳滤波和和功率谱均衡滤波功率谱均衡滤波基于统计特性的图像复原方法，以图像和噪基于统计特性的图像复原方法，以图像和噪声的相关矩阵为基础，是平均的意义上的最佳；声的相关矩阵为基础，是平均的意义上的最佳；

49、 维纳滤波要求图像和噪声为平稳随机场，相互独立，已知它们的频谱密度，维纳滤波要求图像和噪声为平稳随机场，相互独立，已知它们的频谱密度，但是在实际中，往往很难得到这方面的先验知识。但是在实际中，往往很难得到这方面的先验知识。 有约束最小平方复原有约束最小平方复原以平滑度为基础的图像复原方法，以平滑度为基础的图像复原方法，如使得某个函数如使得某个函数的二阶导数为最小；的二阶导数为最小； 采用该方法复原，对每个给定的图像都是最佳的；采用该方法复原，对每个给定的图像都是最佳的； 只需要知道有关噪声的均值和方差的先验知识，就可对每个给定的图像得到只需要知道有关噪声的均值和方差的先验知识，就可对每个给定的

50、图像得到最优结果。最优结果。51 最小二乘方滤波复原以式最小二乘方滤波复原以式(5.89)(5.89)为基础，关键是选择变换矩阵为基础，关键是选择变换矩阵Q Q； 一维情况：一维情况：给定一维离散函数给定一维离散函数f(x)，该函数在某一点，该函数在某一点x处的二阶导数可近似表示为：处的二阶导数可近似表示为： 有约束最小平方复原的最佳准则是使有约束最小平方复原的最佳准则是使 在所有的在所有的x处的和为最小，即处的和为最小，即 或用矩阵形式表示为：或用矩阵形式表示为：min fT CT C f 式中：式中：C为平滑矩阵，为平滑矩阵，f 为图像向量为图像向量22(1)2 ( )(1)ff xf x

51、f xx222fx 2min121xf xf xf x121121121121121121C(5.113)(5.113)(5.112)(5.112)(5.111)(5.111)(5.114)(5.114)52 二维情况，二维情况，f(x,y)在在(x,y)处的二阶导数可用下式近似处的二阶导数可用下式近似上式可用上式可用f(x,y)与下面的与下面的p(x,y)算子卷积得到：算子卷积得到：有约束最小平方复原的最佳准则：有约束最小平方复原的最佳准则：f(x,y)与与p(x,y)卷积的平方和的最小卷积的平方和的最小22224,1,1,1,1fff x yf xyf xyf x yf x yxy010,

52、141010p x y 22222,minx yffxy(5.117)(5.117)(5.116)(5.116)(5.115)(5.115)53 离散卷积计算过程中，用添零的方法扩展离散卷积计算过程中，用添零的方法扩展f(x,y)与与p(x,y)。 如果如果f(x,y)为为AB，而，而p(x,y)为为33，扩展后，扩展后 MA+3-1=A+2，NB+3-1=B+2，即，即 扩展后的卷积为扩展后的卷积为 由二维离散退化模型可知，式（由二维离散退化模型可知，式（5.119）可以写成矩阵形式）可以写成矩阵形式 g = C f 其中其中f是是MN维的列矢量，维的列矢量，C是是MNMN维的分块循环矩阵。

53、维的分块循环矩阵。,01,01,01,1ef x yxAyBfx yAxMByN ,02,02,031,31ep x yxyp x yxMyN1100,MNeeemngx yfm n pxm yn(5.118)(5.118)(5.119)(5.119)(5.120)(5.120)54 平滑矩阵平滑矩阵C重写：重写： 其中每个子矩阵其中每个子矩阵Cj是由是由pe(x,y)的第的第j行组成的行组成的NN维的循环矩阵，即维的循环矩阵，即0M-1M-2110M-122103M-1M-2M-30CCCCCCCCCCCCCCCCC,0,1,2,1,1,0,1,2C,2,1,0,3,1,2,3,0eeeee

54、eeejeeeeeeeepjpj Npj Npjpjpjpj Npjpjpjpjpjpj Npj Npj Npj(5.122)(5.122)(5.121)(5.121)55根据循环矩阵的对角化性质，利用式（根据循环矩阵的对角化性质，利用式（5.44）的矩阵）的矩阵W对对C进行对角化：进行对角化： E = W-1 C W 式中式中E为对角矩阵，其元素是为对角矩阵，其元素是C中元素的二维傅立叶变换。中元素的二维傅立叶变换。约束最小平方复原的最佳准则就是式（约束最小平方复原的最佳准则就是式（5.119）卷积平方和最小，表示成矩阵形式：）卷积平方和最小，表示成矩阵形式： min fT CT C f 变

55、换矩阵变换矩阵Q相当于相当于有约束最小平方复原有约束最小平方复原中的平滑矩阵中的平滑矩阵C，由此可得：，由此可得：利用（利用（5.45）（）（5.123）代入上式：）代入上式： 两边左乘以两边左乘以W-1：gHCCHHfT1TTgWWDEWWEDWWDf*1111gWDEEDDfW*111(5.127)(5.127)(5.124)(5.124)(5.125)(5.125)(5.126)(5.126)(5.123)(5.123)56 对应的有约束最小平方复原滤波器对应的有约束最小平方复原滤波器 在形式上与维纳滤波器相似，在形式上与维纳滤波器相似， 主要区别：只要知道噪声均值和方差，不需其它统计参

56、数，主要区别：只要知道噪声均值和方差，不需其它统计参数， 与维纳滤波器一样，与维纳滤波器一样，是一个调节参数，是一个调节参数， 当当满足满足 时，才能达到最优。时，才能达到最优。*22,eHu vF u vG u vH u vpu v22nfHg(5.128)(5.128)57 估计估计的方法：把估计图像代回到退化系统中，得到的输出会和已退化的图像有差异，的方法：把估计图像代回到退化系统中，得到的输出会和已退化的图像有差异，该差异定义为残差该差异定义为残差r: 由于噪声存在，残差不会等于零，因此残差的范数应该反映出噪声的特征，一般它是由于噪声存在，残差不会等于零，因此残差的范数应该反映出噪声的

57、特征，一般它是的函数，当调整的函数，当调整以达到以达到 时，即为所求的时，即为所求的。 如，对于如，对于NN图像，图像，Sawchuck导出了：导出了： 式中式中 和和 分别表示噪声的方差和均值。分别表示噪声的方差和均值。 可见，图像复原最后可归结为在满足约束的前提下，使得以平滑度为基础的最佳准则可见，图像复原最后可归结为在满足约束的前提下，使得以平滑度为基础的最佳准则达到最优。达到最优。gHCCHHHgfHgrTTT122nr2n2n22mNn2nnm(5.129)(5.129)(5.130)(5.130)58 有约束最小平方复原过程：有约束最小平方复原过程： (1)(1)选一个初始值赋给选一个初始值赋给 ，算得算得 的估计；的估计； (2)(2)计算计算 ； (3)(3)计算残差矢量，并计算残差矢量，并 ； (4)(4)根据根据 的差值决定进一步调整的差值决定进一步调整 ，直到直到 ； (5)(5)此时，将此时，将 代入，得到复原的图像代入，得到复原的图像 。 重申：上