（应用数学专业论文）正则化模型下图像处理的算法设计与实现.pdf

国防科学技术大学研究生院硕士学位论文 摘要 正则化是求解反问题最具普适性且行之有效的方法。本文针对图像处理中正 则化方法的研究及应用进行了广泛深入的探讨与分析，主要工作如下： 第一，运用处理反问题的思想，来研究病态线性方程组的求解，提出了一种 基于解空间分解的加速g m r e s 算法，将解空间分解为k r y l o v 子空间和一个辅助 子空间，其中一部分解用一种加速g m r e s 法迭代得到，另一部分解用直接求解 的方法得到，数值实验和分析表明这种算法是行之有效的，在达到相同的估计精 度的条件下，迭代速度大大提高。大量仿真计算表明求解时间只有普通g m r e s 算法的五分之一，甚至更少；而且在迭代次数相同的情况下，解的精度更高，如 解的均方误差平均是普通g m r e s 算法的五分之三；将该算法应用到光学图像复 原问题中并得以实现，比较了图像经共轭梯度法求解复原和本文方法求解复原后 的效果，实验表明经本文方法复原后图像在信噪比和视觉效果方面都有明显的改 善。 第二，从s a r 图像提高分辨率方法出发，介绍了正则化模型下s a r 图像提高 分辨率的频域正则化方法、复图像域正则化方法及灰度图像正则化方法；针对点 散射模型，灵活应用数学理论合理有效地设计新的求解算法，将正则化理论、g a u s s 消去法、共轭梯度法、n e w t o n 二次迭代法等有效的结合起来并运用到s a r 图像提 高分辨率正则化模型的求解中，并分析了算法的收敛性、计算的复杂度。 主题词：反问题正则化g m r e s 图像复原点散射模型s a i l 第i 页 国防科学技术大学研究生院硕士学位论文 a bs t r a c t n 圮m o s te f f e c t i v em e t h o d st h a tc a nt a c k l et h ei 1 1 p o s e dp r o b l e m sa l et h es o c a l l e d r e g u l a r i z a t i o nm e t h o d sd u r i n gt h ei n v e r s ep r o c e s s t l l i st h e s i si n t r o d u c e sa n da n a l y s e s i nd e t a i lt h er e s e a r c ha n da p p l i c a t i o no fr e g u l a r i z a t i o ni ni m a g ep r o c e s s ，t h em a i n c o n t r i b u t i o n sa n dc r e a t i v i t i e sa r el i s t e da sf o l l o w s ： f i r s t l y ，t h i sp a p e rr e s e a r c h e st h es o l u t i o no fl i n e a ri l l - p o s e de q u a t i o n su s i n gt h e i d e ad u r i n gt h ei n v e r s ep r o c e s s ，t h e np r e s e n t sad e c o m p o s i t i o nm e t h o dt h a t s p l i tt h e s o l u t i o ns p a c ei n t oak r y l o vs u b s p a c eb ya na c c e l e r a t e dg m r e sm e t h o da n da n a u x i l i a r ys u b s p a c et h a tc a nb ec h o s e nt oh e l pr e p r e s e n tp e r t i n e n tf e a t u r e so ft h es o l u t i o n n u m e r i c a lr e s u l t ss h o wt h er e l i a b i l i t ya n de f f i c i e n c yo ft h en e wm e t h o d ：w i t ht h es a m e p r e c i s i o n , t h ec p u t i m ei 8o n l yo n ef i f t h so ft h a to ft h eg m r e sm e t h o d ；a n dw i t ht h e s a m ei t e r a t i v es t e p s ，t h ep r e c i s i o ni sh i g h e r , f o ri n s t a n c e ，t h em s ei st h r e ef i f t h so ft h a t o ft h eg m r e sm e t h o do na v e r a g e t h em e t h o di sa p p l i e dt os o l v et h ep r o b l e m so f r e g u l a r i z e di m a g er e s t o r a t i o n s i m u l a t i o nr e s u l t ss h o wt h a tt h es i 印a 1t on o i s er a t i o ( s n r ) i sh i g h e rt h a nt h ec gm e t h o da n dt h ev i s i o ne f f e c ta r ea l s oi m p r o v e d s e c o n d l y ，t h i sp a p e ri n t r o d u c e s t h e r e g u l a r i z a t i o nm e t h o df o rs a ri m a g e s u p e r r e s o l u t i o ni nf r e q u e n c yd o m a i n , c o m p l e xi m a g ed o m a i na n dp o w e ri m a g ed o m a i n b a s e do np o i n ts c a t t e r i n gm o d e l ，m a t h e m a t i c a lt h e o r i e s a r ea p p l i e dt od e s i g nn e w s o l v i n ga l g o r i t h m sa n dm a n ym e t h o d si n c l u d i n gg a u s se l i m i n a t i o nm e t h o d ，c o n j u g a t e g r a d i e n ta l g o r i t h m ，n e w t o nq u a d r a t i ci t e r a t i v em e t h o da r ec o m b i n e dt os o l v et h e r e g u l a r i z a t i o nm o d e lf o rs a ri m a g es u p e r e s o l u t i o n ，n l ec o n v e r g e n c ea n dc o m p l e x i t y o fa l g o r i t h m sa r ei n v e s t i g a t e d k e yw o r d s ：i n v e r s ep r o b l e m s ，r e g u l a r i z a t i o n , g m r e s ，i m a g er e s t o r a t i o n , p o i n t s c a t t e r i n gm o d e l ，s a r 第i i 页 国防科学技术大学研究生院硕士学位论文 表目录 表2 2 1 线性方程组中矩阵a 的性质1 4 表2 2 2 三个方程组收敛性质比较1 4 表2 3 1b a b a r a 退化图像复原的i s n r 比较1 7 表3 2 1 点目标位置与幅度参数2 4 表3 2 2 求逆法点目标位置与幅度参数2 5 表3 2 3g a u s s 消去法点目标位置与幅度参数2 5 表3 2 4c g 法点目标位置与幅度参数2 6 表3 2 5 三种方法重构时间比较2 6 表3 3 1 点目标位置与幅度参数3 1 表3 3 2 不动点迭代法点目标位置与幅度参数3 2 表3 3 3n e w t o n - c g 法点目标位置与幅度参数3 2 表3 3 4 不动点迭代法与n e w t o n - c g 法重构时间比较3 2 第1 i i 页 国防科学技术大学研究生院硕士学位论文 图目录 图1 1 本文的总体结构框图8 图2 3 1 基于解空间分解的加速g m r e s 算法的图像复原框图1 7 图2 3 2b a b a r a 模糊图像复原对比结果：1 8 图3 2 1 原始含噪图像( 仃= 5 ) 2 4 图3 2 2 求逆法含噪图像的重构结果( 仃= 5 ) 2 5 图3 2 3g a u s s 消去法含噪图像的重构结果( 仃= 5 ) 。2 5 图3 2 4c g 法含噪图像的重构结果( 盯= 5 ) 2 6 图3 3 1 原始含噪图像( 仃25 ) 3l 图3 3 2 不动点法含噪图像的重构结果( 盯= 5 ) 3 1 图3 3 3n e w t o n c g 法含噪图像的重构结果( 盯= 5 ) 3 2 独创性声明 本人声明所呈交的学位论文是我本人在导师指导下进行的研究工作及取得的 研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其 他人已经发表和撰写过的研究成果，也不包含为获得国防科学技术大学或其它教 育机构的学位或证书而使用过的材料。与我一同工作的同志对本研究所做的任何 贡献均已在论文中作了明确的说明并表示谢意。 学位论文题目：垂则焦搓型王圈倦缝堡煎篡洼遮i 士生塞塑 学位论文作者签名：堡盐 日期：2 0 0 7 年j j 月b1 3 学位论文版权使用授权书 本人完全了解国防科学技术大学有关保留、使用学位论文的规定。本人授权 国防科学技术大学可以保留并向国家有关部门或机构送交论文的复印件和电子文 档，允许论文被查阅和借阅；可以将学位论文的全部或部分内容编入有关数据库 进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。 ( 保密学位论文在解密后适用本授权书) 学位论文题目：垂则丝搓型王图像矬垄鲍篡洼遮盐盏塞理 学位论文作者签名：垫丑日期：2 口。7 年i1 月侈日 1n 作者指导教师签名：j 例7 日期：夕年月日 国防科学技术大学研究生院硕士学位论文 1 1 研究背景 第一章绪论 自2 0 世纪6 0 年代以来，在图像处理、模式识别、遥感技术、地球物理、生 命科学、材料科学、工业控制乃至经济决策等众多的科学技术领域中，都提出了 “由效果、输出反求原因、输入 的反问题，通称“数学物理反问题。由于此 类问题有着广泛而重要的应用背景，其理论又具有鲜明的新颖性与挑战性，反问 题已成为应用数学中发展和成长最快的领域之一【1 】【刁。 数字图像处理已日渐成为一门广受关注的学科。并且随着数字摄影技术和雷 达成像技术的不断发展，其应用范围和规模也发生了巨大的变化d 】- 铂。 光学图像复原就是通过对图像退化的过程模型化，对退化图像进行处理，使 其趋向复原为没有退化前的理想图像。由于在各个不同的领域中图像的来源可能 是千差万别的，导致图像降质的原因也可能不同。但是对于光学图像复原问题， 却有本质的共同点：图像复原与图像降质是一个互逆过程，图像降质可以用一个 卷积过程来描述，因而图像复原可以通过反卷积来实现1 。 在二维情况下，由于卷积机制而造成图像降质的直观表现就是图像的分辨率 降低。照相机光学系统不精良、拍摄时对焦不准、拍摄时相机有移动等，都会使 照片变得不清晰。这些降质现象可以用卷积过程来描述，在成像过程中，景物上 的一个点不是仅仅反映到图像上的一个点，而是被扩散成图像平面上的一个区域。 因此，图像上的每个点是景物的许多个点经过混合叠加而成，这在数学上可以用 一个叠加积分式来描述： y c u ，v ) = 双嚣嚣版m 墨f 为和，f ) 出奶+ 善 ，v ) ( 1 1 1 ) 式中，似，功和0 ，d 分别表示像平面和物平面上点的二维空间坐标，h ( u ，v , s ，) 描述 成像传输的点扩展函数性质，称为点扩展函数或降质函数，它一般表示了一个处 于物平面上( s ，f ) 位置的景物被成像一传输系统扩展成像平面( 1 ，) 上的一个二维函 数。s ) 是一个非线性算子，善( 甜，v ) 表示加性噪声。在常规的图像复原研究中，总 是将非线性的影响加以忽略，从而对( 1 1 ) 进行简化可得： y ( u ，) = 岱窿h ( u s ，一t ) x ( s ，t ) d s d t + 善( “，1 ，) ( 1 1 2 ) 图像复原的一个基本问题就是反降质或去模糊，即依据一个降质并受噪声污染的 观测值y ( u ，) 来估计原来的图像x ( s ，r ) 。如果成像系统的点扩展函数h ( u ，v ) 已知， 则反降质问题是一个常规的反卷积问题，否则它是一个盲反卷积问题。 但是，由于点扩展函数和输入信号的卷积等价于对输入信号的一种滤波，而 第l 页 国防科学技术大学研究生院硕士学位论文 实际的点扩展函数总是形成一个低通滤波器，使输入信号的高频成分受到抑制甚 至丢失。因此，反卷积就是一个逆过程，它要“找回 丢失了的高频成分。一般 的，在反卷积过程中不可避免的会将观测值中的高频噪声放大，使得反卷积的结 果有可能严重偏离真实的解。因此，为了获得尽可能精确的解，在反卷积的过程 中需要加进抑制噪声的考虑，在信号复原与噪声放大之间做出适当的折衷。这里 需要指出的是，即使系统函数h q ) 或点扩展函数h ( u ，v ) 己知，反卷积问题也不是一 个简单的问题。在这种情况下，( 1 1 2 ) 是第一类f r e d h o l m 积分方程，它的求解是 一个病态问题。一个方程的病态表现为方程的解不是连续依赖于观测数据。也就 是说，观测数据的微小变动可能导致解的很大变动。由于观测数据受噪声的污染， 而反卷积过程又伴随着对高频噪声的放大。因此，方程的解可能会偏离真解相当 远，这是关于反卷积病态的一种直观了解。另一方面，反卷积是属于数学物理问 题中的一类“反问题 。在数学物理问题中，信号和图像经过一个系统后得到输 出是正问题，在数学上就是卷积积分，它的反问题就是图像复原和系统辨识。 合成孔径雷达( s y n t h e t i ca p e r t u r er a d a r ，s a r ) 是一种全天候、全天时的高 分辨率微波遥感成像雷达，自发明以来，其理论和技术在世界范围内受到极大重 视并获得飞跃发展，已成为空间对地观测发展的“热点 口】阳】。同时，许多应用领 域对s a r 图像的分辨率提出了很高的要求。因此，高分辨率是s a r 技术发展的主 要目标之一，利用相关技术与方法进行提高s a r 图像分辨率处理必然成为s a r 的 研究重点和热点。 s a r 距离剖面频域成像模型为阳1 ： g = r f + w( 1 1 3 ) 其中，g 为频域复回波数据( 相位历史数据) ，厂是场景散射系数，w 表示噪 声。该模型描述了场景后向散射系数与回波数据的直接关系。式( 1 1 3 ) t i p 为所考虑 的s a r 图像的离散观测模型。希望根据观测量g ，利用可以获得的有关图像的先 验信息，得到尽可能不失真的原来量厂，从而提高其分辨率。s a r 图像提高分辨 率可以通过成像系统的点扩散函数的伪逆实现反卷积来得到，该逆问题是病态的， 从数学模型的特征分析来说s a r 图像提高分辨率问题就是一个反问题。 无论是光学图像复原问题，还是s a r 图像提高分辨率问题，在数学上都是一 个反问题。求解反问题面临的两个本质性的实际困难是n 3 ： ( 1 ) 原始数据可能不 属于所论问题精确解所对应的数据集合( 例如积分算子或微分算子的值域) ，因 而，在经典意义下的近似解可能不存在； ( 2 ) 近似解的不稳定性，即：原始资料 的小的观测误差( 这在实际中是不可避免的) 会导致近似解与真解的严重偏离。 概言之，反问题常常是不适定的，若不用特殊的方法求解，将得不到合理的解答 或答案。 第2 页 国防科学技术大学研究生院硕士学位论文 反问题的求解已发展了各种方法，诸如脉冲谱技术( p s t ) 、广义脉冲谱技术 ( g p s t ) 、最佳摄动量法n 们、蒙特卡罗方法( m o n t ec a r l om c 吐l o d ) n 、各种优化 方法和正则化方法n 2 1 等。其中，最具普适性、在理论上最完备且行之有效的方法， 就是由著名学者t i k h o n o v 以第一类算子( 特别是积分算子) 方程为基本数学框架， 于2 0 世纪6 0 年代创造性地提出、后来得到深入发展的正则化方法。其基本思想 是：用一族与原问题相邻近的适定问题的解去逼近原问题的解。从而，如何构造 “邻近问题”而获得所谓的正则算子和正则解、如何控制与原问题的“邻近程度 而决定与原始资料的误差水平相匹配的正则参数以及上述工作的快速数值实现， 就成为正则化理论和方法的三大核心问题。 正则化方法的数值求解往往采用迭代的方法，而且迭代方法还有这样的优点： 当将无穷维的问题转化为有限维问题后，迭代过程不会破坏系数矩阵的结构，而 且存储量较之直接求解方法要节省。这对于求解大规模的问题是相当有利的。作 为提高图像分辨率与清晰度的正则化方法，是一种常用的解决病态反问题的有效 方法，具有深刻的数学理论与应用前景，这就对基于图像处理中的正则化方法的 深入研究提供了良好的理论与工程背景。 1 2 国内外研究现状 用迭代方法处理各种正问题已有悠久的历史，同样，迭代法还可用来求解反 问题。研究发现，用迭代法求解反问题时会出现所谓的“半收敛”现象，即：在迭代 的早期阶段，近似解可稳定地得到改进，展现出“自正则化”的效应，而迭代次数超 过某个阈值后便会趋向于发散。因而关键问题是要寻找一个合适的终止原则，使 得迭代次数与原始数据的误差水平相匹配。研究表明，迭代指数，即迭代步数正 好起到正则化参数的作用，而这个终止准则对应着正则化参数的某种选择方法。 另外，使用迭代方法显然还有以下几个优点：( 1 ) 可以避免对算子求逆；( 2 ) 关于解的一些先验知识可以合并到迭代过程中；( 3 ) 在求解的过程中可以对解进 行监控；( 4 ) 可以加入约束来控制噪声的影响。因此在求解正则化问题时首选迭 代方法。 以下是几种正则化问题的迭代方法：l a n d w e b e r 迭代法、v a nc i r e r t 迭代方 法、最速下降方法、共轭梯度法和迭代t i k h o n o v 正则化的求解方法，以及正则化 方法的快速数值实现n 1 ： l a n d w e b e r 迭代法 考虑第一类算子方程 a x = y( 1 2 1 ) 的正则化求解问题，设真解为矿。 第3 页 国防科学技术大学研究生院硕士学位论文 l a n d w e b e r 建议使用如下的迭代格式 略= - l + 国么( y 一么_ 1 ) ( 1 2 2 ) 来求解第一类算子方程的近似解。其中0 缈去为松弛因子，表示算子彳 z 的伴随算子，x o = ，r 为初始猜测值，不失一般性，可取，= 0 。 通常在( 1 2 1 ) 式中，右端项y 由观测得到为，带有误差水平万：8 y 一忙艿。 这时l a n d w e b e r 迭代为 = 硪i + 刎一磁1 ) ( 1 2 3 ) 一般来说l a n d w e b e r 迭代序列的收敛速度较慢，常采用加速l a n d w e b e r 迭代 法，即半迭代法。半迭代法又名多项式加速方法，其基本思想是：每步迭代不仅 仅利用上一步的信息，而且要充分利用前面各步所得的信息。 记 r k l = 彳杪一以一1 ) ( 1 2 4 ) 半迭代点列以如下方式产生： 噍= 以七一l + + 心+ q ，j 一l ， k = l ，2 ， 一七+ + 船= 1 ，0 ( 1 2 5 ) 可见，l a n d w e b e r 迭代法是式( 1 2 5 ) 当c o k = 缈，1 2 i = = 心= 0 时的特例。 v a nc i t t e r t 迭代方法 考虑方程( 1 2 1 ) 的正则化求解问题。如果4 的尺寸太大，使得标准的代数方程 求解方法无法实施，则v a nc i t t e r t 迭代就是一个不错的选择。设初始值= 缈， 迭代格式为： 毫+ l = 筑f l ( y - a 主q )(124-a x k2 6 )黾+ l2 &)【i o j 可以看出，迭代的每一步工作量在于计算以。因为a 是一个循环矩阵，则以 可由卷积和包含卷积的运算来形成。假定真解为工，它使( 1 2 1 ) 成立，即x = a y 。 于是 x 一量= x 一圣- , 8 ( a x 一出)( 1 2 7 ) k + l七 令误差向量气= x 一氟，则有 e k + l = ( ，一f l a ) e 女= ( x - 肛) “1 e o( 1 2 8 ) 式中，e o = x - 磊是起始误差。设彳有个不同的特征值，则有相似变换 a = p 。1 人p ，其中人是彳的特征值丑o = 1 , 2 忉构成的对角阵，p 是一个满秩矩阵。 于是 e k + l = p q ( i - 肚) “1 p e o ( 1 2 9 ) 第4 页 国防科学技术大学研究生院硕士学位论文 对任何e o ，为了使e - - - h0 ，充分必要条件为：1 1 一跳l 1 ，对于任意的i ； 从而0 0 ， 只要 o u ( u 1 ，z 乞) 万( 占)( u l ，, 2 u ) ( 2 3 2 ) 就有 所( 五，7 , 2 ) 占( 彳毛= u i 彳乞= ”2 ) ( 2 3 3 ) 反之，若上述三个条件中至少有一个不能满足，则称其为不适定的。 适定性的三个条件无疑具有深刻的实际背景。首先，对实际问题而言，我们 自然期望其解是存在且唯一的。更重要的是，实际获取的数据资料都是近似的， 第9 页 国防科学技术大学研究生院硕士学位论文 即我们实际处理的是近似数据( 4 ，) ，而不是“精确”的数据( 爿，“) ，因为“精确 的数据往往是未知的。若原始数据的很小的误差将导致近似解对于真解的严重 偏离，则计算所得的数值结果将毫无意义。 不适定问题的求解过程从大的方面可归结为：重新定义原问题的近似解；近 似解确定后，还需要稳定的数值方法，来确定原问题的稳定的数值解。对于不适 定问题的求解，已经有文献作了大量的研究，其中比较普遍采用的框架就是正则 化方法，即用一簇与原问题相“邻近 的适定问题的解去逼近原问题的真解n 羽。 应当认识到，反问题的不适定性是问题本身所固有的一种特征：如果没有关 于欲求解的问题的附加信息( 先验信息) ，这一本质性的困难是无法克服的。反 过来说，就是要依据所能提供的关于解的先验信息，尽可能多、尽可能稳定的恢 复有关原问题的部分信息。 对于图像处理问题，由于涉及到大规模的方程组求解，法方程的维数太大， 因而代数方程求解的一些基本技术就变得不现实。而迭代方法处理各种正问题( 无 穷维或有限维的适定问题) 已有悠久的历史，并且迭代方法还有这样的优点：当 将无穷维的问题转化为有限维问题( 线性代数方程组) 后，迭代过程不会破坏系 数矩阵的结构，而且存储量较之直接求解方法要节省。这对于求解此类大规模的 问题是相当有利的。 2 1 2 光学图像复原问题的病态性 光学图像复原问题的有效性关键之一取决于描述图像退化过程模型的精确 性。但在实际的图像处理过程中，图像均需以数字离散函数表示，所以必须将退 化模型离散化。 对于下式 g ，力= ff ( a ，p ) h ( x 一口，y - f 1 ) d a d f l + n ( x ，力 ( 2 1 4 ) 如果将式中f ，h ，n ，g 按相同间隔采样，产生相应的阵列【f ( i ，纠叙口、 【 ( f ，朋c d 、【刀( f ，力】爿。口、 g ( f ，j ) 】彳。b ，然后将这些阵列补零增广得到大小为mxn 的周期延拓阵列，为避免重叠误差，这里m 么+ c 一1 ，n b + d 一1 。由此，当 七= 0 , 1 ，m 一1 ；l = o ，l ，一l 时，即可得到二维离散退化模型形式： ，一l 一i ( 七，z ) = 五( 1 ，歹) 魂( 七- i ，z 一歹) + ( 是，f )( 2 1 5 ) 1 = 0 = 0 如果用矩阵表示上式，则可写为： g = h f + n( 2 1 6 ) 其中，厂，g ，行为一个行堆叠形成的m n x l 列向量，日为m n x m n 阶的块循环矩阵。 第1 0 页 国防科学技术大学研究生院硕士学位论文 假设函数厂属于某一空间h ，函数g 属于是某一个空间皿，寻找一个变换算 子丁，对厂进行丁变换得n g ，即r = g ，图像复原问题就是找到一个逆变换， 使得丁- 1 ( g ) = 厂。 反问题就是讨论逆变换的存在性和唯一性。理论表明，成像的方程是第一类 f r e d h o l m 积分方程，它的求解是一个病态问题。因此在对图像复原问题求解时会 出现解不存在，或解不唯一，或解不连续依赖于观测数据。解不连续依赖于观测 数据表现为观察数据的微小变动可能导致解的很大的变动，数学上表示为 z 一1 ( g + s ) = f + 8 ，其中占为任意小的值，万为造成厂的相应扰动，当病态时万 占。 为了解决病态问题，即转化为良态问题，t i k h o n o v 提出正则化方法，即使下 面的泛函最小化： m ( a ，力= g ( f ，g ) + a n ( f )( 2 1 7 ) 式中第一项为数据逼近项，表示解与原数据之间的距离，即保证求得的解逼 近真解；第二项使为了克服病态问题，当数据带有误差时，保证解在真解附近连 续地依赖于数据，不会把误差或噪声过分地放大；仃是控制第一项和第二项权重 的正则化参数。在理想的情况下，如果上式达到最小，且第一项趋于o 时，参数 也趋于0 。 根据以上t i k h o n o v 正则化的观点，一般常规的正则化方法是把图像复原问题 规划为以下优化问题。即通过建立目标函数j ( 厂) ，使 f = a r g m i n j ( f ) ( 2 1 8 ) l 其中，j ( f ) = i i g 一矿肟+ 甲( 门。右式第一项为数据逼近项( 其最小化实际观测量 和由系统模型对重建域产生的观测值之间的平方误差) ，第二项反映我们将采用的 先验信息，重建图像依赖于我们所使用的限制条件，因而甲( 厂) 的选择非常重要。 在图像复原中，正则化的基本思路即是利用解的先验知识，构造附加的解的 某种先验信息作为约束条件，改变求解策略，限制解的数据，修改空间和拓扑， 使用统计估计方法等，使图像复原的不适定问题( 病态问题) 变为适定问题( 良 态问题) ，从而使解变得稳定和确定，这正是正则化方法发挥作用的地方。 2 2 基于解空间分解的g m r e s 算法 病态线性方程组的求解一直是一个活跃而重要的研究方向。从数值代数的观 点来看，有其独立的研究价值；但从反问题的角度来看，可将它看成一类特殊的 反问题。因而，现有的许多求解病态方程组的方法或技巧可用于反问题的数值求 解n 1 。 第1 l 页 国防科学技术大学研究生院硕士学位论文 从处理反问题的角度来讨论病态线性方程组的求解。考虑求解线性方程组 a x = b ，么r n x n , x e 础，b r 一( 2 2 1 ) 由于彳r 删，故它是有限维空间中的有界线性算子，从而是紧算子，故方程组 ( 2 2 1 ) 是第一类算子方程。方程组( 2 2 1 ) 可能会出现以下几种情况： ( 1 ) 当d o t ( a ) = 0 时，因为彳奇异，故其逆算子不存在；但对于任意的b 舭， 它有无穷多个解；这时显然是一个不适定问题。 ( 2 ) 当诎似) 0 时，彳有逆算子，从而方程组( 2 2 1 ) 的解存在且唯一。但是 根据其条件数c o n d ( a ) = r a i i i a 一10 大小的不同，又会有两种情形： 若它是良态的( w e l l - c o n d i t i o n e d ) ，即条件数较小时，方程组( 2 2 1 ) 的解对于右 端项b 的扰动是不敏感的( 稳定的) ，故求解式( 2 2 1 ) 是一个适定问题； 若它是病态的( i l l c o n d i t i o n e d ) ，即条件数甚大时，右端项b 的小扰动将会引起解 的大变化，这时求解式( 2 2 1 ) 是一个不适定问题n 1 。 2 2 1g m r e s 算法的基本理论 考虑求解线性方程组( 2 2 1 ) ，为了尽可能减少存储空间和计算开销，k r y l o v 子空间迭代法n 4 1 是求解方程组( 2 2 1 ) 的有效方法。当系数矩阵彳对称正定，共轭梯 度法( c g ) 或预共轭梯度法( p c g ) 可快速准确求解方程组( 2 2 1 ) 的近似解；当彳 对称但不正定时，极小残量法或预极小残量法则能有效求解方程组( 2 2 1 ) 。对于一 般的非对称矩阵，常采用广义极小残量法( g m r e s ) h s 、共轭梯度法( c g ) 来求解。 g m r e s 算法每次迭代产生的近似解z ，x o + k ，( t o ，彳) 满足最小化条件 忙一鸽i f ：2 难而m 蚂i n ( 叫) l l b 一出。其中= 6 一如，而为迭代初值。算法利用a m o l d i 过程产生k r y l o v 子空间k i 的正交基向量，a m o l d i 过程的每一次迭代需要利用所 有前面迭代生成的正交基计算产生下一个正交基n 引。 文献【1 7 】提到了一种加速的g m r e s 算法： 构造一个空间r 联p ，在此空间中找出a c o = d 叫的近似解万，从而获得方 程组( 2 2 1 ) 的新的近似解可= + y + 枥。其中，d 件1 是a m o l d i 过程产生的正交 基向量，巧= ( q ，呸9 o - 9 d ， 是利用a m o l d i 过程构造的子空间k ，的正交基矩阵。 根据砌= d ，+ 1 解的性质，文献【1 7 】给出了的一种构造方法。 2 2 2 基于解空间分解的加速g m r e s 算法 文献 1 8 1 提出了求解大型离散不适定问题的分解方法，该方法较之g m r e s 算 第1 2 页 国防科学技术大学研究生院硕士学位论文 法能够在相同估计精度条件下减少求解次数。本节将这种分解方法进行了改进， 在将k r y l o v 解子空间进行分解后，其中一部分解改用加速的g m r e s 算法n 迭代 得到，而文献 1 8 中贝l j 采用普通g m r e s 算法迭代。由于求解其中部分解的速度变 快，并且结合了文献【1 7 】和文献【1 8 】中两种算法的优势，数值实验表明本节算法在 计算速度和解的精度上优于这两种方法。 假设形r 蒯中标准正交化列向量张成的空间代表人为选定的线性空间。并且 引入直交投影易= w w r 和磷= i - 易。按照 = 弓+ 弓，弓= ，弓= 唠勺，j = l ，2 ， ( 2 2 2 ) 将计算得到的近似解x ，分解。在文献 1 8 1 中，是通过普通的g m r e s 算法迭 代来得到的。而本节将加速的g m r e s 算法应用到求解# 中，这样，在求解t 的 过程中迭代速度变快，从而整个算法迭代速度要快于文献【18 】中方法。而，通常选 得比较小，如1 ，3 ，则可以通过直接求解的方法来得到t 。 下面给出基于解空间分解的加速g m r e s 算法的详细实现过程： 首先，对矩阵进行舛分解n 卅 a w = 鲫( 2 2 3 ) 其中q r 州有标准正交化的列向量，尺r i 。i 是上三角矩阵。因为在这里选 取的，很小，所以分解( 2 2 3 ) 可快速实现。这里假定矩阵r 是非奇异的，可以通过 把，取得足够小来实现。并假设吃= q q r ，访= i - 吃。 其次，与( 2 2 2 ) 类似，把x 分解为易x 和唠x ，可得到如下分解表达式： p q a p 矽x + p q a e w x x = 印 ( 2 2 4 ) 呓彳磁x = 呓6 ( 2 2 5 ) 路4 = o ( 2 2 6 ) 再次，利用文献 1 7 中提到的加速g m r e s 迭代方法来解大型线性系统( 2 2 5 ) 。 从( 2 2 6 ) 中可以得到呓彳彤= 呓彳，因此可将这种加速的g m r e s 迭代方法应用于 呓= 眵 ( 2 2 7 ) 假设t 是( 2 2 7 ) 的最佳近似解。则： 。= 彤z(22x 8 ) = 馆z ， u 通过解方程( 2 2 4 ) 来得到方程( 2 2 1 ) 中近似解的一部分弓，这样，弓= 名勺 满足 钙= 矽一尼q ( 2 - 2 9 ) 而这个方程组和下面的小型线性方程组等价： 第1 3 页 国防科学技术大学研究生院硕+ 学位论文 r z = q 7 p q ) ( 2 2 1 0 ) 因此，可以通过直接解方程组的方法来得到方程( 2 2 9 ) 的解，于是： 弓= 嗡 ( 2 2 1 1 ) 最后，i 1 荩过( 2 2 2 ) 得到 当 嘻万 时迭代终止，艿为预先取定的收敛判断参数。 2 2 3 数值实验与结果分析 ( 2 2 1 2 ) 使用m a t l a b 对上节算法进行数值实验，比较g m r e s 算法、加速g m r e s 算 法和2 2 2 节提出的新算法的收敛性质。在实验中，初始向量为= 0 。算法l 代 表g m r e s 算法n 副；算法2 代表加速g m r e s 算法n7 l ；算法3 代表基于分解的 g m r e s 算法n 射；算法4 代表2 2 2 节提出的新算法。实验在w i n d o w sx p 环境下 运行，电脑c p u 为p 4 ，2 8 6 g h z ：内存为2 5 6 m b 。 以这种方式构造矩阵彳：4 = s d s 一1 ，其中丘s ，d r 1 0 0 。l o o ，d 为对角阵，s 为 双对角阵，其对角线元素都为l ，上次对角线元素都为实数，。 表2 2 1 线性方程组中矩阵彳的性质 线性方程组 ，d 矩阵条件数 l0 9 d i a g ( 1 ，2 ，10 0 )1 5 1 0 2 21 1 d i a g ( 1 ，2 ，1 0 0 )1 2 i 0 9 3o 9 d i a g ( - l o ，- 1 ，l ，l o o )4 2 1 0 2 表2 2 1 显示出实验三中所需的三个矩阵彳的特点，其中d i a g ( 1 ，2 ，1 0 0 ) 表示 对角线上元素依次为l ，2 ，1 0 0 。取占= 彳( 1 ，1 ，1 ) r 1 0 0 x l 的列向量；并且 添加的误差项p 是均值为零，标准偏差为0 0 0 1 的g a u s s 噪声；那么6 = 占+ e ；在 本实验中选定形2 i - 赤o o ( i ，l ，o r r 1 0 0 c l ；收敛判断参数艿= 1 胪。 表2 2 2 三个方程组收敛性质比较 方程组1算法1算法2算法3算法4 运行时间 0 1 8 so 1 0 s0 0 8 s0 0 4 s 迭代次数 1 3 753 迭代残量 9 4 1 0 92 3 1 0 98 9 1 0 - 95 1 1 0 - 9 第1 4 页 国防科学技术大学研究生院硕士学位论文 方程组2算法1 算法2算法3 算法4 运行时间 0 1 8 s0 1 3 so 1 0 s 迭代次数不收敛 1 4 l l 7 迭代残量 8 1 1 0 - 9 7 9 1 0 - 93 2 1 0 - 9 方程组3算法1算法2算法3算法4 运行时间 2 2 4 s0 1 7 so 1 s0 0 6 s 迭代次数 9 01 31 05 迭代残量 8 8 1 0 - 95 6 l i f o4 7 1 0 - 99 4 1 0 。1 0 从表2 2 2 可以看到，g m r e s 算法只对方程组1 有比较快的收敛速度，而对 方程组2 却不收敛，因为矩阵彳比较病态；对方程组3 收敛但速度很慢，因为矩 阵彳含有负的特征值。相比之下，加速的g m r e s 算法和基于分解的g m r e s 算法 都比g m r e s 算法收敛快，而本文算法将这两种方法结合，从计算速度和求解精 度来看都优于前三种方法。 2 2 2 节提出的基于解空间分解的加速g m r e s 算法为什么具有在迭代速度上 具有优势呢? 首先，在q r 分解中，由于所选的w r 利中，取得比较小，所以分解 ( 2 2 3 ) 可快速实现。其次，在求解大型线性系统( 2 2 5 ) 时，采用了一种加速的g m r e s 迭代方法，这本身就提高了算法的计算速度。最后，在求解线性方程组( 2 2 1 0 ) 时， 由于其规模较小而采用直接求解的方法也提高了整个算法的计算速度。 2 3 光学图像复原结果 图象退化模型，其矩阵形式为： g = h f + n ( 2 3 1 ) g ， 刀r ( n 是图像像素个数) 分别是c a - - 维离散观测退化图像、原始图像和加性 噪声按列拉长所得到的列向量。日是由空间退化的点扩展函数( p s f ) 生成的矩阵 2 0 1 【2 l 】 常见的t i k h o n o v 正则化的目标泛函形式为： m ( a ，力= i 矽一吖+ 口l c f ( 2 3 2 ) 这里i 可一g r 为数据的逼近项，它反映了观测图像g 对原始图像f 的逼近程度；口 为正则化参数，它用来平衡近似解的逼近程度和平滑性；l 明z 为正则项，它根据 图像的平滑性对解进行约束，此时c 称为正则算子，它是由一高通滤波算子g 生 成的大小与日相同的分块循环矩阵，通常巴取为二维l a p l a c i a n 算子，比如： 第1 5 页 国防科学技术大学研究生院硕士学位论文 0 1 0 1 阽扑- 4 1l 【o 1 o j 对于正则化目标泛函( 2 3 2 ) ，正则化图像复原可以看作约束最优化问题，即原 图像厂的最优估计值厂为： 夕= a r g n l i l l 0 彤一g i l 2 + 口l 盯1 2 ( 2 3 3 ) 对于给定的正则算子c ，并适当地选择正则化参数口，可得到正则解厂： 夕= ( 日r h + t z c 丁c ) - 1 日r g ( 2 3 4 ) 最小化目标泛函( 2 3 2 ) 则得到如下法方程： ( h t h + a c r c ) f = h r g ( 2 3 5 ) 令a = h t h + a c t c ，b = h r g ，x = f ，则上式可写为一般的代数方程形式 a x = b ( 2 3 6 ) 2 3 1 算法的设计 对于方程( 2 。3 6 ) ，利用基于解空间分解的加速g m r e s 算法迭代求解。实验证 明，用该方法求解可以有效改善收敛性，并且在较少的迭代步数内达到收敛。 计算步骤如下： s t e p l 置初始值= 0 ，并令8 = 1 0 - 8 s t e p 2 用基于解空间分解的加速g m r e s 算法迭代求解方程( 2 3 6 ) ，在第歹步 的值为； 蛳若皆蛾则灿玳渖则聊肼髅断嘶2 。 在迭代运算中，正则化参数吩司5 k 老芋随着每一步的迭代而自动更新。 本节算法设计的基本框图如下： 第1 6 页 国防科学技术大学研究生院硕士学位论文 图2 3 1 基于解空间分解的加速g m r e s 算法的图像复原框图 2 3 2 仿真实验结果 仿真实验中所处理的是2 5 6 2 5 6 尺寸的0 2 5 5 灰度级的b a b a r a 图像。 这里用改进信噪比( i s n r ) 来衡量算法的复原性能。i s n r 表征了复原前后图 像相对于原始图像地信噪比变