2 2 2 2 2 2 2 数 学 杂 志 J . of Mat h. ( PRC) 3 Vol. 29 (2009) No . 4 图像恢复的一种快速迭代正则化方法 吕小红 , 吴传生 (武汉理工大学理学院 ,湖北武汉 , 430070) 　　摘要 : 本 文 研 究 了 将 图 像 恢 复 问 题 转 化 为 大 型 的 线 性 不 适 定 问 题 的 求 解. 利 用 由 Landweber 迭代正则化方法改进所得到的快速收敛的迭代正则化方法 ,处理具有可分离点扩散 函数的图像恢复问题. 图像恢复实验表明该方法可大大提高收敛速度 ,且在计算中只需要较少 的存储量. 　　关键词 : 不适定问题 ;图像恢复 ;迭代正则化 ; 点扩散函数 　　MR( 2000) 主题分类号 : 65F22 ; 68 U10 　　　　中图分类号 : TP391 　　文献标识码 : A 　　　　文章编号 : 0255 7797( 2009) 04 0563 04 1 　引言 图像恢复 ,是指去除或减轻在获取数字图像过程中发生的图像质量的退化. 这些退化 形式是多种多样的 ,如传感器噪声 、摄像机未聚焦 、物体与摄像设备间的相对移动 、随机大 气湍流 、光电系统的相差 、成像光源或射线的散射等 ,这些因素都会使成像的分辨率和对比 度退化 ,从而造成宇航 、遥感 、天文 、侦察照片 、电子显微图片 、医学射线图片等的退化. 因而 如何从模糊的 、不真实的退化图像得到逼真的图像具有重要的现实意义[ 1 ] . 图像恢复问题是一个二维反卷积问题 ,是典型的具有不适定性的反问题. 本文讨论的 范畴限于积分核可分离的二维卷积型问题上 ,即如下的第一类 Fredholm 积分方程 ∫1 0∫1 κ( x - x′)ω( y - y′) f ( x′, y′) dx′d y′= g ( x , y) 0 (1) 在实际情况中 ,大部分点扩散函数是可分离的. 对于不适定问题 (1) ,若不采取特殊的方法 求解 ,将得不到合理的结果 ,而正则化方法是解决此类问题的有效方法[ 2 6 ] ,另外图像复原问 题规模较大 ,迭代正则化方法是处理此类问题较好的方法. 但是迭代法通常收敛速度较慢 , 比如 Landweber 迭代法 ,因此在实际应用中 ,寻找一种快速的迭代方法至关重要. 本文将利用一种快速的迭代正则化方法求解点扩散函数可分离的图像复原问题 ,该方 法同传统的 Landweber 迭代正则化方法相比 ,可大大提高收敛速度 ,而且具体的计算中 ,只 需要较少的存储量. 2 　迭代正则化方法 利用矩形积分公式将 (1) 离散化处理 :首先 ,将对 y′的积分近似为 n 项求和 收稿日期 : 2008 作者简介 : 吕小红 (1969 08 12 　　　接收日期 : 2009 03 16 ) ,男 ,湖北黄陂 ,硕士生 ,讲师 ,从事反问题研究. E mail : lvxiaohong2008 @sina. com 

       465 数 　　学 　　杂 　　志 　　　　　　　　　　　　　Vol. 29 　 ∫1 ω( y - y′) f ( x′, y′) dy′≈ n- 1 ∑ l = 1 0 ω( y - y′l ) f ( x′, y′l ) = φ( x′, y) 这里 y′l是积分节点 , f 代表计算的近似解 ;其次 ,用另一组和逼近对 x′的积分 ∫1 κ( x - x′)φ( x′, y) dx′≈ m - 1 ∑ k = 1 0 κ( x - x′k)φ( x′k , y) = ψ( x , y) n m 这里 x′k也是积分节点. 并且在 m n 个点 ( x i , y j ) 上取 ψ( x i , y j ) = g ( x i , y j ) 　i = 1 , …, m 　j = 1 , …, n 这样 ,最终就可以构造成一个 m n 阶的线性方程组 ,包含有 m n 个未知量 四个矩阵 A , A , F , G: A ik = m - 1κ( x i - x′k) , A jl = n - 1ω( y j - y′l ) , Fkl = Gij = g ( x i , y j ) , ( i , k = 1 , …, m , j , l = 1 , …, n) . 则 (1) 离散化为 f ( x′k , y′l ) , f ( x′k , y′l ) . 引入 于是问题 (1) 就可转化为如下最小二乘问题 : A F A T = G A F A T - G = min F = min 其中 · 为 Euclid 范数 , (3) 的解为 A G ( A T ) 这里可以利用矩阵间的 Kronecker 乘积 . 将 (2) 转化为标准的线性方程组的形式[7 ] 其中 ,“vec”是这样定义的 ,如果 X = ( x1 , x2 , …, x n) 是 m ×n 矩阵 ,那么 ( A A) vec ( F) = vec ( G) vec ( X) = ( x1 T , x2 T , …, x n T ) T 是 m n ×1 矩阵. (2) (3) (4) 下面以标准的线性方程组的形式介绍我们的迭代方法[8 ] . 对于条件数非常大的病态线 性方程组 假设其右端项实际已知的是带有观测误差的 bδ ,且已知 b - bδ ≤δ. Cx = b 定义算子 R k R k = k- 1 ω[ ∑ j = 0 ( I - ωC T C) ( 1 a ) j ] a C T (5) (6) 其中 a 是正整数 ,0 < ω≤ 1 CT C 是松弛因子. 当 k →∞, R kb →x , x 是方程 (5) 的最小模最小 二乘解 , k 是正则化参数. 以 R k作为 C 当 a = 1 时 ,此方法即为 Landweber 迭代法 ,当 a 取大于 1 的整数时 , 将加快求正则解的收 敛速度. k = R kbδ作为 x b 的逼近. 的逼近 ,以 xδ = C 定理 1 　假设τ> 1 , x = ( CT C)νz ,选取正则化参数 k = k (δ) 使得 A R kbδ - bδ ≤τδ 第一次出现. 则有 2 =Ο(δ 2ν R kbδ - x 1) k =Ο(δ- (2ν+ 1) a) 　2) 定理证明见文献[8 ]. 迭代参数的选取是迭代正则化方法中至关重要的问题. 迭代参数的取法有先验和后验 两种方式 , 先验取法基于精确解的一些条件 ,这通常是未知的 ,故基于数据误差水平信息和 数据本身的后验取法更为实用 ,定理 1 解决了迭代参数的后验选取问题. 2ν+ 1 ) . 考虑迭代式 : D k + 1 = D k ( I - ( ωCT C) 1 a ) + I , D1 = I. 显然 D k = ∑ ( I - ( ωC T C) 1 a ) j . 而正 k- 1 j = 0 

　No . 4 　　　　　　　　吕小红等 　图像恢复的一种快速迭代正则化方法 565 则解 xδ k = R kbδ = D k C T bδ ,因此 ,当 C 的规模不是很大时 ,我们可使用下列迭代算法 : 算法 1 1) 将 CT C 正交对角化 , Q T C T CQ =Λ; 2) 取 3) 计算 D k + 1 = D k ( I - ( ωQ T C T bδ) ,若 CR kbδ - bδ ≤τδ终止 ,否则返回 3) . 4) 计算 R kbδ = QD a k ( ω为 C T C 的最大特征根的倒数 , D1 = I ; 1 a ) + I ; ωΛ) 3 　图像恢复的迭代算法 将图像恢复问题转化为标准的线性方程组的形式 (4) 后 , 我们可以利用算法 1 进行求 A 的规模型相当大 ,为了减少存储量和计算量 ,具体到方程 (4) ,还必须利 解 ,但是由于 A 用 Kronecker 乘积的性质对算法 1 作适当的改进. 首先注意到 Kronecker“乘积满足 ( A B) ( C D) = ( A C B D) 和 ( A B) T = ( A T B T) , 在算法 1 中的第 1 步 ,对于 CT C = ( A 别对 A T A 和 A T A 正交对角化来实现. 若 A T A 及 A T A 已正交对角化 ,有 Q T A T A Q =Λ1 , P T A T A P =Λ2 ,则 ( Q A) = ( A T A ) P) T ( ( A T A ) ( A T A ) ) ( Q A) T ( A P) =Λ1 Λ2 =Λ. ( A T A ) 的正交对角化可以分 其次 ,由于 vec( A X B) = ( B T A) vec( X) . 在算法 1 中的第 4 步 , R k vec ( G) = ω( Q ω( Q = ω( Q = P) D a P) D a P) D a P T ) ( A k ( Q T k ( Q T A T k vec ( P T A T G A Q) A) T vec ( G) P T A T ) vec ( G) 令 vec( S k) = D a 即第 k 次迭代得到的恢复结果为 F k = ωPS kQ T ,同时 k vec ( P T A T G A Q) , S k与 F 行列数相同 ,则 R k vec ( G) = ωvec( PS kQ T ) , ( A A) R k vec ( G) - vec ( G) = ωA PS kQ T A T - G F . 最终对于 (2) 的求解算法为 ω为 ( A T A ) 算法 2 1) 将 A T A 正交对角化 , Q T A T A Q =Λ1 ;将 A T A 正交对角化 , P T A T A P =Λ2 ;Λ1 2) 取 3) 计算 D k + 1 = D k ( I - ( 4) 计算 Fk = 4 　数值实验 ωPS kQ T ,若 Fk满足 A Fk A T - G F ≤τδ停止 ,否则返回 3) . ( A T A ) 的最大特征根 (Λ对角线上的最大数) 的倒数 , D1 = I ; a ) + I , vec( S k) = D a k vec ( P T A T G A Q) ; ωΛ) 1 Λ2 =Λ; 我们考虑 Matlab 中的“Cameraman”图像 (256 ×256) 被高斯点扩散函数 1 2πσ exp ( - u2 + v2 2σ2 ) 模糊并加了信噪比为 50 dB 的高斯白噪声 (均值为 0) 的复原问题. 实验中取σ= 2 ,τ= 1. 1. 实验中的相关结果及数据 (计算机用时 、迭代次数 、复原图像与原始图像的相对误差 Fk - F F F F ) 见表 1. 

2 665 数 　　学 　　杂 　　志 　　　　　　　　　　　　　Vol. 29 　 算法 2 中 取值 Landweber 迭代法 ( a = 1) 算法 2 ( a = 2) 算法 2 ( a = 3) 算法 2 ( a = 4) 表 1 　实验中的相关数据 CPU 时间 (秒) 441. 579352 25. 579249 11. 294643 6. 770146 迭代次数 748 43 17 10 相对误差 0. 0967 0. 0952 0. 0937 0. 0932 该方法采用迭代法只需要较少的存储量 ,同时实验表明我们的迭代复原算法通过对 a 取值的调整可显著减少迭代次数 ,有效节约算法时间. 参考文献 : [ 1 ] 　邹谋炎. 反卷积和信号复原[ M ]. 北京 : 国防工业出版社 , 2001. [ 2 ] 　Engl H. W , Hanke M. and Neubauer. A Regularization of Dordrecht , 1996. Inverse Problems [ M ]. Kluwer , [ 3 ] 　Ramlau R. Modified Landweber Methods for Inverse Problems [J ]. Number. Funct . Anal. Optimiz. 1999 20 : 79 98. [ 4 ] 　Hanke M. Accelerated Landweber Iterations for the Solution of Ill posed Equations [ J ]. Numer. Math. 1991 60 :341 373. [ 5 ] 　曾三友 ,康立山 ,丁立新等. 一种基于正则化方法的准最佳图像复原技术 [J ]. 软件学报 , 2003 , 14 (3) : 689 696. [ 6 ] 　张海燕 ,闵 　涛 ,艾克锋. 二维截断奇异值分解方法在图像恢复中的应用 [J ]. 计算机工程与应用 , 2008 , 44 (1) : 60 62. [ 7 ] 　Loan. C. V. Comp utational Frameworks for the Fast Fourier Transform[ M ]. SIAM , Philadelp hia , PA , 1992. [ 8 ] 　吕小红 ,吴传生 ,周 　俊. 一种快速收敛的迭代正则化方法[J ]. 吉林大学学报 (理学版) ,2009 ,47 (2) : 269 272. A FAST ITERATIVE METHOD OF REGULARIZATION FOR IMAGE RESTORATION ( S chool of S ciences , W uhan U ni versit y of Technolog y , W uhan 430070 , Chi na) L Xiao hong , WU Chuan sheng Abstract : Image restoration is usually converted to large Iterative method of regularization is an effective method for large In this paper , we apply an imp roved fast iterative method of regularization based on Landweber iterative method to image restoration problems with separable PSF. The image restoration experiment shows that our method is efficient and less storage is needed. po sed problem. scale linear ill posed p roblem. scale linear ill Keywords : ill 2000 MR Subject Classif ication : 65F22 ; 68 U10 posed problems ;image restoration ; iterated regularization ;point sp read f unction