[理学]CT统计重建论文.doc

基于最大似然和罚似然估计的基于最大似然和罚似然估计的 CT 统计重建算法研究统计重建算法研究 摘 要 当 CT 系统在数据采集过程中受到严重的噪声干扰，或者投影数据采集不完全时， 用解析重建算法重建出的图像有伪迹。而统计重建算法具有物理模型准确、对噪声不 敏感、易于加入约束条件等优点，其重建的图像质量要优于传统的 FBP 方法。针对 CT 统计重建，本文主要研究了以下内容： 研究了基于最大似然估计的统计重建算法：期望最大值法(Maximum Likelihood Expectation Maximization，ML-EM)，可分离抛物面型替代(Separable Paraboloidal Surrogate，ML-SPS)算法和这两种算法的有序子集形式（OS-ML-EM 和 OS-ML-SPS） 。 仿真实验结果表明 ML-SPS 的初始收敛速度比 ML-EM 快，但该两种算法的收敛速度 远慢于其对应的有序子集形式。用 OS-ML-EM 算法和 OS-ML-SPS 算法对 CT 采集的 实际投影数据进行重建，得到了较好的重建结果。 研究了基于罚似然(Penalized Likelihood，PL)的统计重建算法:OSL(One Step Late)- EM 算法和 PL-SPS 算法。重点讨论了基于 Gibbs 分布的罚函数，从理论上分析了势函 数需要满足的条件以及势函数的选取对图像的影响，着重分析了二次势函数和 Huber 势函数的优缺点，通过仿真实验对它们的重建结果进行误差分析。 提出了罚似然重建中正则项参数的自适应选取的方法，该方法充分利用每一次迭 代的重建结果信息，不断对正则项参数进行更新。并用于仿真实验和重建 CT 实际投影 数据，重建结果表明该方法降低了噪声，提高了图像质量。 将 OS-OSL-EM 算法推广到三维锥束轨迹下的图像重建，取得了较好的仿真结果。 关键词：最大似然估计，CT 重建，罚似然估计，有序子集 The Study of CT Statistical Reconstruction Algorithm Based on Maximum Likelihood and Penalized Likelihood Estimates Abstract When the CT data has serious noise in the process of acquisition system and projection data is incomplete , analytic reconstruction algorithm get images with artifact. The statistical reconstruction algorithm with a accurately physical model isnt sensitive to noise and easy to add constraints etc.Therefore, the reconstructed image quality is superior to conventional FBP methods. For the CT Statistics reconstruction, This paper mainly studies the following content: The maximum likelihood estimate for the statistical reconstruction algorithms is studied and mainly includes expectation maximum algorithm (ML-EM),separable paraboloidal surrogate(ML-SPS) algorithm and the both algorithms with ordered subset(OS-ML-EM and OS-ML-SPS). In simulation experiments, it shows that the initial convergence rate of ML- SPS is faster than the ML-EM ,however, the both algorithms is slower than OS-ML-EM and OS-ML-SPS algorithms. At last, OS-ML-SPS and OS-ML-EM algorithms are used to reconstruct actual CT projection data and get better image. Penalized maximum likelihood(PL) reconstruction for the statistical algorithms is based on the penalized likelihood estimates and adds a penalty term to suppress noise, which mainly includes OSL-EM algorithm and PL-SPS algorithm. It is focused on the penalty function which is based on Gibbs distribution and analyzes the need of potential functions to meet the conditions, the influence on image for the selection of potential functions and the advantages and disadvantages of the quadratic function and the Huber function. At last, simulation experiments give the error analysis of reconstructed image. This method which is an adaptively regularied CT image reconstruction can adaptively choose regularization parameters with making full use of the results of every iteration to update regularization parameters. Simulation results show that this method reduces the noise and improve image quality. The method is used to reconstruct the actual CT projection data to demonstrate the feasibility and practicality. OS-OSL-EM algorithm is applied to three-dimensional cone-beam image reconstruction and gets the better image. Key words: maximum likelihood estimate, CT Reconstruction, Penalized maximum likelihood, ordered subsets I 目目 录录 第一章 绪论 1.1 课题研究背景及意义 1 1.2 国内外研究现状 1 1.3 论文的主要工作及结构安排 4 第二章 CT 统计重建的基础 2.1 CT 成像理论.5 2.2 统计重建的数学模型 7 2.3 优化变换原理 8 2.4 图像重建算法评价.9 第三章 基于最大似然准则的 CT 统计重建算法 3.1 最大似然准则的基本理论 11 3.2 ML-EM 算法.12 3.3 ML-SPS 算法14 3.4 有序子集的统计重建方法.17 3.4.1 有序子集 ML-EM17 3.4.2 有序子集 ML-SPS19 3.5 实验结果及分析 20 3.5.1 ML-EM 算法与 ML-SPS 算法的重建结果21 3.5.2 OS-ML-EM 算法与 OS-ML-SPS 算法的重建结果.23 第四章 基于罚似然的 CT 统计重建算法 4.1 罚似然的基本理论.25 4.2 OSL-EM 算法 .31 II 4.3 OS-PL-SPS 算法.33 4.3.1 PL-SPS 算法.33 4.3.2 有序子集 PL-SPS 算法.34 4.4 自适应选择正则项参数 35 4.5 实验结果与分析.36 4.5.1 OSL-EM 算法实验结果分析.37 4.5.2 OS-PL-SPS 算法实验结果分析 38 4.5.3 自适应选择正则项参数的实验结果分析40 第五章 锥束轨迹下有序子集罚似然统计(PL)重建 5.1 锥束扫描下的投影矩阵的生成 42 5.2 锥束下的 OS-OSL-EM 算法重建结果及分析.44 第六章 总结与展望 6.1 研究主要内容及成果 47 6.2 存在的问题及以后的工作展望 48 参考文献 攻读硕士学位期间发表的论文 致 谢 1 第一章 绪论 1.1 课题研究背景及意义 计算机层析成像技术是利用有某种能量的射线源对物体进行断层扫描，根据所获 得的某种物理量(指物质对射线的衰减系数)，运用特定的重建算法，重建出物体某个断 层的无影像重叠的二维图像。自 1967 年 CT(Computed Tomography)在英国问世以来， 技术迅速发展，被公认为是 20 世纪影响人类发展的十大技术之一，其应用涉及到医学、 工业、石油、物探、材料、生物、安全等众多领域。在国防工业中，CT 已成为航天器、 固体火箭发动机、高新技术等产品关键部件检测的必不可少的工具。由于 CT 技术是发 展国防的基础技术之一，被发达国家列为敏感技术，作为国家机密进行严格的控制。 因此，CT 技术被视为国家科技实力的标志之一。 从数学上说，CT是从投影到重建图像的反问题1，具有普遍性，在数学界已经引 起了广泛的重视。重建算法是工业CT技术中较为关键的一部分，CT重建算法主要可以 分为滤波反投影重建算法和迭代重建算法。当投影数据完备、噪声不很严重时，解析 法(滤波反投影(FBP)23等)可以得到很好的重建图像，但当CT系统在数据采集过程中 受到严重的噪声干扰，或者投影数据采集不完全时，用解析重建算法得到的图像有伪 迹。而迭代重建算法(统计迭代算法和代数迭代算法等)可以用来重建不完全数据、动态 数据和噪声数据。其中统计迭代重建方法的优势是其对噪声的鲁棒性，通过加入先验 模型和统计规律可以有效的抑制噪声，提高图像质量。 本课题主要研究统计算法在 CT 重建中的应用，通过研究有序子集的方法来提高收 敛速度，加入惩罚项，约束项来降低噪声的影响，提高图像质量，为国内三维工业 CT 的研制提供理论基础和技术支持。 1.2 国内外研究现状 2 近十几年来， CT 广泛的应用在医学和工业上，相应的 CT 算法也迅速发展。其 中解析算法重建速度较快，易于实现，是目前 CT 图像重建技术中最常用的算法。但是 该算法通常要求采样数据是完全的，而且对噪声比较敏感，这在一定程度上影响了重 建图像的质量，限制了它在实际中的应用。相对解析算法、代数迭代算法而言，统计 迭代算法能准确描述投影数据的物理模型，对误差不敏感，易于加入约束条件等优点 逐渐引起人们的重视，然而它与代数重建算法4一样重建速度慢，运算时间长，这些缺 点极大地限制了它的应用，但是随着加速算法和计算机硬件技术的发展，统计重建将 在实际中得到了广泛地应用。 一个完整的统计重建框架包括图像的离散化系统的物理几何模型测量的统 计模型(泊松分布)优化准则(最大似然准则和罚似然准则)求解方法。第步是考虑 如何对重建的图像进行建模，将待重建的图像离散成二维或者三维的像素矩阵，第 步考虑如何确定投影矩阵，与系统的几何结构和建模方法相关5-10，第步 ij AaA 确定投影数据的模型，在统计重建中，一般认为投影数据服从泊松分布。第步是引 入一定的优化准则。第步在的基础之上进行问题的求解，该问题的求解实际上是 一个参数估计的数学优化问题。 统计重建的优化方法中涉及到优化变换(optimization transfer，OT)，增量方法 (incremental method)，有序子集加速(ordered subset，OS)，这三种优化方法是统计重建 方法发展的主线11。 在统计重建中，最大似然估计(Maximum Likelihood，ML)是一种常用方法。1977 年，A. P. Dempster 等人把期望值最大化(Expectation Maximization, EM)算法引入图像 重建之中。1982年，Shepp12和Vardi13首次把期望最大迭代算法14应用到CT算法中， 使似然函数最大化，之后ML-EM算法得到了广泛的应用。1997年，Fessler15等人提出 了基于替代函数(surrogate function)的优化理论，其基本思想是在M步(maximization step)直 接对E步(expectation step)中所得到目标函数的替代函数作最大化处理，这一突破性的改 进在很大程度上应当归功于De Pierro的凸性(convexity)算法1617。当然替代函数的选取 需要满足一定的条件，其它替代函数算法还包括Fessler18-21，Zheng22提出的算法以及 全局收敛的增量优化传递算法。 3 1994年，Hudson和Larkin提出了有序子集(Ordered Subsets) 23EM算法(OS-EM)，大 大提高了收敛速度。OS-EM算法将投影数据分解成有限个有序子集，每次迭代时只处 理其中一个子集的数据，由于加速效果明显，OS-EM算法及其各种变形算法此后便被 广泛地应用。其中RAMLA(row-action maximun likelihood algorithm)24重建算法的思想 来源于Herman的代数重建算法和Hudson等人的OS-EM算法，该算法也是将投影数据分 成一系列不正交的投影数据子集，并引入一个松弛因子，该算法可以从理论上证明其 收敛到ML的收敛点，并且其收敛速度和OS-EM算法差不多。虽然子集类算法收敛速度 快，但不能全局收敛，迭代到一定次数会出现极限环，为了解决该问题，S. Ahn 和J. A. Fessler提出了改进的块迭代EM算法25，E .Quan和D .S Lalush等人提出了基于快速子 集的凸算法在CT中的应用26。 最大后验概率重建算法(maximum a posteriori，MAP)2728是一种贝叶斯重建算法, 也可以称为罚似然重建(Penalized Likelihood，PL)29-32，该算法和 ML 重建算法可以说 是一对关系亲密的姊妹，因为 ML 算法是寻求合适的图像的估计值使得图像得到的已 知投影的概率最大；而贝叶斯重建算法是从已知的投影出发，要求在给定投影情况下 所求图像的概率最大，通过选择合适的先验分布模型来提高重建结果的质量，其作用 等同于优化理论中的正则项（惩罚项） 。正则化技术最早由 Tikhonov 在著作中提出， 他研究了求解不适定问题稳定解的基本理论，并提出了著名的 Tikhonov 正则化技术。 先验模型的引入虽然能改善重建效果，但同时增加了估计的难度。就泊松观测模型而 言，它本身属于非线性模型，因此它的 ML 估计和 MAP 估计均没有闭型解析式。针对 这一问题，Herbert 最先提出 GEM(generalized EM)算法33，Green 给出了经典的 OSL(one-step-late)算法34，Fessler35提出了罚似然的 SAGE(space-alternating generalized EM)算法，另外 De Pierro 的3637凸性算法也涉及到了最大后验估计的问题。在以后的 工作中，De Pierro 和 Yamagishi38将 RA(row-action)思想与最大后验估计的思想相结合， 提出了 BSREM(block sequential regularizedexpectation maximization)算法。Chung Chan, Roger Fulton39等人提出了自适应结构先验在 CT 中的应用，根据图像结构块中灰度的 不同自适应的选择滤波。Xiaochuan Pan40阐述了最小化 TV 正则项(total variation)在图 像重建中的应用。Jing Wang, Tianfang Li41等人把基于最小二乘的罚函数方法用于投影 数据的去噪，提高了图像重建的质量。 4 1.3 论文的主要工作及结构安排 本文 CT 统计迭代重建算法为研究课题，主要做了以下工作： 第一章主要介绍了课题的研究背景、意义，综述了 CT 及其发展历史，CT 技术的 发展意义及统计重建算法的研究现状。最后提出了本文的研究内容，确定了本文的研 究方向。 第二章首先给出了 CT 统计重建的物理和数学基础，然后介绍了优化变换原理，最 后介绍了图像重建质量的评价标准。 第三章首先详细介绍了最大似然估计类的重建算法：ML-EM 算法和可分离抛物面 型替代(Separable Paraboloidal Surrogate，ML-SPS)算法，它们的推导过程和实现步骤， 然后对有序子集的方法进行阐述，并详细介绍了 OS-ML-EM 算法和 OS-ML-SPS 算法 的实现步骤。实验中对比了 ML-EM 算法和 ML-SPS 算法的实验结果，并用 OS-ML- EM 算法和 OS-ML-SPS 算法对 CT 采集的实际投影数据进行重建，得到了较好的重建 结果。 第四章先介绍了罚似然估计的理论知识，重点介绍了基于 Gibbs 分布的罚函数， 从理论上分析了势函数需要满足的条件以及势函数的选取对图像的影响。然后介绍了 OSL (One Step Late)-EM 算法和 PL-SPS 算法的推导和实现，接着提出了自适应选择正 则项参数的罚似然重建，通过对正则项参数的改进来提高图像的质量。算法实现中选 用了二次势函数和 Huber 势函数来重建图像，对它们的重建结果进行误差分析，并用 提出的自适应选取正则项参数的方法来重建 CT 实际投影数据，取得了较好的结果，验 证了方法的实用性。 第五章介绍了 OS-OSL-EM 算法在三维锥束轨迹下的应用，首先介绍了锥束扫描 结构下的投影矩阵的计算，然后用 OS-OSL-EM 算法对仿真投影数据进行重建，并与 OS-OSL-EM 算进行比较分析。 第六章对本课题的研究内容进行了总结，指出本文的创新之处和解决的问题，并 提出了尚未解决的问题，指明了今后研究的重点。 5 第二章 CT 统计重建的基础 在过去的几十年里，统计迭代法一直是断层重建研究的热点，其重建图像的过程就 是利用参数估计，优化变换原理等来逼近最优解一个过程。在重建过程中，首先确定 能够正确反映成像系统的物理模型和统计模型以及输入输出关系，然后利用已知的条 件确定求解方法。而且统计重建可以利用与物体相关的约束条件和其他的与特定系统 相关的边界条件，利用精确的物理模型和适当的统计模型，可以在数据缺失和不规则 采样的情况下重建出满意的图像。 2.1 CT 成像理论 CT 图像重建的原理是用穿透力强的射线扫描被测物体，当一束强度大致均匀的射 线投照到物体上时，射线一部分被吸收和散射，另一部分透过物体沿原方向传播。由 于物体的各种结构在密度、厚度等方面存在不同，它们对射线的吸收量也不相同，从 而使透过物体的射线强度分布发生变化，通过 CT 采集、转换得到投影数据，并运用相 应的数学模型，进行图像重建出二维或三维图像，最后通过算法将射线强度分布转换 成图像上的灰度分布，形成人眼可视的图像。 取一理想的 X 射线源，它发出极细的笔束 X 射线，在被测物体对面置一探测器， 如图 2.1(a)所示。假设强度为的 X 射线穿过均匀分布衰减系数为的物体，行进了 0 Iu 的距离，强度变为，按 Beer 定理有xI 00 ln(/ ) ux ii euxII 或 （2.1） 如图 2.1(b)所示，若物体是分段均匀的，衰减系数分别是,相应的长度为 1,23 ,.u u u 则下式成立： 1,23 ,.x x x 1230 .ln(/ )uxuxuxII （2.2） 6 更一般，若物体在平面内为不均匀介质，则在某一方向 ，射线沿某一路径的总衰xylL 减值为 00 ( , )ln(/ )( )exp( , ) LL x y dlIII LIx y dl （2.3） 其中表示物体在二维平面上的衰减系数。是检测到的强度；是射线经( , )x y( )I LL 过的直线，是直线的弧长。dl 若没有具体路径，只说了沿着某一个方向，那么是投影，由等式(2.3)可知，dl 测得与，即可知道。在 CT 中，入射的射线强度与出射的射线强度之比经过 0 IIdl 取对数运算后，则成了沿射线路径上衰减系数的线积分。CT 重建问题实际上是给定一 个待重建物体的被测量的投影即线积分，然后运用算法去计算物体的衰减分布即dl 被积函数。 x o I 射线源 I 探测器 o I 射线源 I 探测器 1 x 2 x 3 x N x 1 2 3 N （a）均匀物体 （b）非均匀物体 图 2.1 单能 X 射线束在物体中的衰减图示 1917 年，丹麦数学家雷当(JRadon)研究已经为 CT 技术建立了数学基础，从数学 理论上证明了某种物理参量的二维分布函数，由该函数在其定义域内的所有线积分确 定。该结论指出了如果没有无穷多个且积分路径互不完全重叠的线积分，只能确定实 际分布的一个估计近似值。 有了上述数学、物理基础后，为了实现工程技术上的应用，还需要解决两个主要 问题。首先是如何采集到检测断层衰减系数线积分的数据集，其次是如何充分利用该 数据集确定出衰减系数的二维分布。解决第一个问题可采用不同的扫描方式，即用射 7 线束穿过被测体所检测断层并相应进行射线强度测量的规律性，可采用各种不同的扫 描检测模式围绕提高扫描检测效率。解决第二个问题则是选择合适的图像重建算法， 一类是解析法，以滤波反投影算法为代表。另一类是迭代重建法，有代数重建和统计 重建。统计重建是本论文研究的主要方向，下面具体介绍该算法。 2.2 统计重建的数学模型 在 CT 图像重建中，图像区域是重建区域，首先应当明确，图像区域是一个正方形， 其中心在坐标原点；图像函数是一个二元函数，其值在图像区域之外为零。一个具有 个元素的网格把这个图像区域分成 2 n个相等的正方形，每个小正方形(像素)内图像函n 数的值是均匀相等的。一幅图像的数字化是这样一幅的数字化图像，使得原nnnn 图像在 2 n个元素网格的任何像素上的积分值等于在同一个像素上的数字化积分值。在 某一点的图像密度值就是在该点物质组织在某种能量下的相对衰减系数值，将物体的 重建区域离散化为一幅数字化图像后，直接从投影数据来估计数字化图像每个像素点 的密度。 图2.2 CT数据采集 如图2.2所示，将待重建物体划分成个互不重叠的子区域，像素序列依次从1到M ，射线源和探测器是点状的，它们之间的连线对应一条射线。M 核物理研究表明，X光子的辐射满足泊松随机过程，结合这一统计特性，在CT中选 用泊松统计模型，可以表示成如下： 8 iii ypoissonArx:(1,2. )iN （2.4） 其中表示重建图像矢量，用表示随机变量的个测量 12 ( ,.) M x xxx 12 (,.,) N y yyyN 值, 表示第 个投影的测量误差，为投影矩阵，表示第 个射线穿过第个 i ri ij Aa ij aij 像素的长度。表示第 条射线的投影值。 1 M iijj j Aa x xi 2.3 优化变换原理 优化变换42(optimization transfer，简称OT)的思想是将复杂的统计重建的优化问题 转换成易于求解的优化问题，转换后的问题一般具有易于优化、分离变量、降低维数 的特点。 OT的优化问题是寻求目标函数的最大值，表达式如下：( ) x argmax( ) x D xx （2.5） 由于该目标函数的最优解比较难求，可以选择一个用来逼近目标函数极大值的替 代函数，通过不断的迭代来求解使替代函数取得极大值。之所以要用迭代的 ( ) ( ,) n x x 方法来求解最小值问题是由替代函数的特性决定的，替代函数可能为非二次曲线，不 易求最大值，计算量非常大，这些都给求解带来了困难，所以一种可行的方法就是选 择一个可以用来逐步逼近替代函数极大值的替代函数，通过不断的迭代计算来得到替 代函数的极值。但是替代函数一定要满足一定的条件才可以得到稳定的结果，首先替 代函数的目的就是为了简化求解极大点的计算量，如果替代函数是一个多维的非二次 函数，就难于计算。那么替代函数就要尽量将多维的求解极值问题变为一维的问题来 简化求解过程。另外为了保证收敛的一致性，如图 2.3 所示，要求替代函数在当前点与 替代函数相切并且位于替代函数的下部，则替代函数满足的条件如下： ( )( )( ) (;)() nnn xxx 9 （2.6） ( ) ( ;)( ) n x xx （2.7） 由条件(2.6)，(2.7)可以得出 ( )( )( )( ) ( )()( ;)(;) nnnn xxx xxx （2.8） 换句话说，我们选择 使得来使得x ( )( )( ) ( ;)(;) nnn x xxx ( ) ( )() n xx 图 2.3 替代函数曲线图 其中为替代函数，为第 次迭代后的结果。只有替代函数满足了这上述条 ( ) ( ,) n x x ( )n xn 件才能在简化计算的同时保证得到稳定一致的收敛结果。 选定了替代函数以后，求解最大值的过程就可以简化成式(2.9)所示的迭代过程， 每次迭代都是以上次的结果为一点，得到在这一点上的替代函数并求解极大点作为本 次迭代的结果。则该优化问题的迭代过程可以简化成： 1( ) 0 argmax ( ,) nn x xx x （2.9） 2.4 图像重建算法评价 10 对于重建算法的可行性,一般按照下面三个方面来评价: 1) 算法是否易于实现。 2) 计算量的大小是衡量算法可行性的一个重要指标。计算时间是一个相对量, 在 同一计算机上，可以通过几种算法计算时间来衡量。 3) 迭代算法的误差分析。用归一化均方差距离、归一化平均绝对距离、相关系数、 信噪比等评价指标反映重建图像的误差和噪声。 相关系数 CC(Correlation Coefficient)衡量： 1 1/2 22 11 ()() ()() P recrecrefref jjjj j PP recrecrefref jjjj jj xxxx CC xxxx （2.10） 归一化平均绝对距离 NMAE (Normalized Mean Absolute Error)来衡量： 11 ()/() MM recrefref jjj jj MAExxx （2.11） 信噪比 SNR (Signal Noise Ratio)： 22 11 10log() /() ) MM refrefrecref jjjj jj SNRxxxx （2.12） 归一化均方差距离 NMSE(Normalized Mean Square Error): 1 2 2 2 recref jj i refref jj i xx MSE xx （2.13） ，分别为被检测图像和重建图像对应像素的值，为被检测图像的平均灰度 ref j x rec j x ref j x 值，为重建图像的平均灰度值。相关系数反映了被测试图像和重建图像之间的相 rec j x 11 似程度，相关系数越大，NMAE 越小说明重建效果越好。用信噪比 SNR 来评价图像的 噪声的大小，SNR 越大，图像所含的噪声越小。NMSE 越小表示图像与原始图像的误 差越小，图像质量越好。 第三章 基于最大似然准则的 CT 统计重建算法 最大似然估计考虑了数据的统计特性，能够充分利用系统的固有分辨率，从而重 建的图像分辨率和噪声特性均优于卷积反投影重建的重建结果43。但是在图像重建优 化求解的过程中，似然函数的解析表达式很难找到，因而限制了 ML 的使用。但是期 望最大法通过构造“完全数据集”的似然函数，充分利用完全数据似然函数的简单性 去求解实际投影数据的最大似然估计，从而解决了这个难题。 3.1 最大似然准则的基本理论 图像重建的极大似然模型是建立在假定（X 射线发射的光子数是服从泊松 （Poisson）分布）的基础上的。最大似然准则为：假设是相互独立的，已知的条yy 件下寻求x使下面的对数似然函数最大 ( )lnLPxy x （3.1） 12 基于 Poisson 分布的似然函数为P y x 1 () ( ) ! i i y N Ax i i i A SPe y x xy x （3.2） 假定满足以下两个条件 ij a 01 ij a （3.3） 1 1, M ij i a (1,2,)jM （3.4） 则(3.2)可以简化得到下式： 1 1 1 () ( ) ! iM ijj j M y ijj a x N j i i a x Se y x （3.5） 图像重建的极大似然方法就是求出图像矢量，使似然函数最大，也 12 , T M x xxxL 就是使对数似然函数最大。因此，ML 估计为：( )lnLPxy x argmax ( ) ML Lxx 1 argmaxln()ln(!) N ii ii i AyAy xx 11111 argmax(log()log(!) MNNMN jijiijji jiiji xaya xy （3.6 ） 另外，还可以根据判定极值的二阶充分条件（即目标函数的 Hessian 矩阵是半正定 的）来说明似然函数的凸性，从而可知 ML 估计(3.6)存在且唯一。 13 3.2 ML-EM 算法 目标(3.6)的优化具有一定的难度，由于目标函数是非线性的，不存在解的闭型表 达式；而且在求解过程中隐含着非负性约束条件，为了解决该问题，Shepp44和Vardi45建 议使用EM算法对(3.6)进行优化。1994年，EM 算法被成功地应用在图像重建中，成为 研究不完全的数据极大似然估计的一种方法。EM 算法的收敛解是非负的，迭代形式 便于计算机实现，在图像重建算法研究中，EM算法的应用相当重要。 EM算法认为直接被观测到的投影数据总是不完全(incomplete)的，依赖一个较为y 完全的数据，它们之间的关系由映射表征，即，一般 为已知。称完全R( )yzxZ 数据，它包含了待求未知参数 的所有信息。基本思想就是通过重复处理完全数据的统x 计问题来解决不完全数据的统计问题。首先，假设初始图像为；接着根据求一 (0) x (0) x 次近似图像，再由求二次近似图像，以此类推到第 次迭代的图像。 (1) x (1) x (2) xn ( )n x 根据ML估计，优化的目标函数为，利用优化变换原理，需要寻找( )lnPxy x 一个替代函数，其必须满足(2.9)中的单调性的条件。EM算法在理论上总能 ( ) ( ,) n EM x x 以概率收敛到某个的ML估计46。 ML-EM 算法的迭代过程分两步： E步：在给定当前图像估计及测量得到投影值的情况下，完全数据对应的对 ( )n xy 数似然函数的条件期望为， ( )( ) ( ,)log(), nn EM EZ x xxy x （3.7） M步：最大化替代函数对图像向量 进行估计x 1( ) 0 argmax( ,) nn EM x xx x （3.8） 根据计算步骤，需要找出替代函数，推导过程中应用了算法中的凹性特征，假设，0 i r 当时，则有 ( ) 0 n x 1 ( ) M iijj i j lAa x xx 14 111 ( )log()() NMM ijjiijji ijj La xra xr x ( ) ( )( ) ( )( )( ) 111 log()()() ()() n NMM ijjj nni iiiijji nnn ijj iji a xx r ylla xr lxl xx xx gg （3.9 ） 利用凹性特征有： （3.10 ( ) ( )( ) ( )( )( ) 111 ( )log()log ()() ()() n NMM ijjj nni iiiijji nnn ijj iji a xx r Lyl xl xa xr l xxl x xgg ） 忽略(3.9)式中的常数项，即可得到替代函数 ( ) ( ,) n EM x x ( ) ( ) ( ) 11 ( ,)() log () n MN ijjn EMijjj n ji i a x yxa x l x x xg （3.11） 其中，根据 M 步对(3.11)最大化，求关于任意的偏导数并令其为零即 1 N jij i aa j x ，可以立即得到： ( ) ( ,) 0 n EM j x x x ( ) ( )11 1 0 n NN jiij ij M nii j ijj j xy a a x a x （3.12） 即推到出 ML-EM 算法的迭代公式： ( ) (1)(1) ( ) 1 11 ,1,2.,0 n N jijinn jjNM n i ijijj ij xa y xjM x aa x （3.13） 根据迭代公式看出 EM 算法有两个明显的特点： 1） 每一次迭代后均提高后验密度函数值，即有。 (1)( ) ()() nn ll xx 15 2） 给定非负的初始估计，以后迭代所得的也都是非负的，满足了 (0) x (1)(2) ,xx 图像像素的非负性约束。 ML-EM 算法的具体步骤如下： 1) 设定图像的初始估计值: (0)( 1,2.,); jj xxjM 2) 计算所有射线方程的估计值 (0) 1 M iijj j la x (1,2,)iN 3) 计算 i i i y l (1,2,)iN 4) 计算第个像素的修正值 j 1 1 1 N jiijN i ij i Ca a 5) 对的值进行修正 j x ; )0( jjj Cxx 通过该像素的所有射线对该像素进行修正，从而完成第一轮迭代。 6) 将以上一轮迭代的结果作为初值，重复 2)到 5)可得第 n 轮结果，就得到一序列 若符合收敛要求,即对事先给定的很小的正数存在正整数 NN，使 (0)(1)(2)( ) , n xxxx 当时，有 .nNN ( )(1) ,(1,2, ) nn jj xxjJ 3.3 ML-SPS 算法 SPS(separable paraboloidal surrogates)算法是全局收敛的，并且满足像素的非负性 条件，迭代的初始收敛速度要比SAGE快。SPS 适用于ECT的ML具有Convex形式先验 的MAP问题47-51，Erdogan52将其推广到TCT的ML和MAP问题。该算法是从优化变换 的思想出发，通过构造易于最大化的二次型替代函数简化EM算法的计算复杂度。 首先定义非负集合：:0, MM j RRxj x 定义集合 : M DRAiri xx 根据最大似然估计理论，最大似然SPS算法的目标函数表达式如下： 16 1 ( )( )( ) N pp ii i LhA xxx （3.14） 其中对数似然函数的表达式可以为： , 其中 ( )log()() iiii h lylrlr 1 ( ) M ijj j lAa x xx （3.15） 函数有如下性质：( ) i h l （1），其中 i lr ( )( ) ii h lh l ii lyr （2）在区间上，单调递增，在区间上单调递减。 , ii r l ( ) i h l , i l （3）是凹函数( ) i h l 在 ML-SPS 算法中，首先找到一维函数的抛物面型替代函数，替代( ) i h l ( ) ( ,) pn i ql l 函数需要满足的条件：在当前迭代点处，它与函数相切且在函数 ( ) nn ii lAx( ) i h l 的下方。通过对做二阶 Taylor 展开：( ) i h l( ) i h l ( )( )( )( )( )( )2 1 ( ,)()()()()() 2 pnnnnnn iiiiiiiii ql lh lh lllc lll )( ; ) N npn ii i QqAA x xxx （3.17） 利用的凹性，可以得到的可分离的二次替代函数 ( ) ( ;) n Q x x( )L x ( ) ( ,) pn qx x ( )( )( )( )( )( )( ) 1 ( ,)()() ()() ()() 2 pnnnnnnn j qdiag cLL x xxxxxxxxxx （3.18） 17 其中为的一阶偏导， 为曲率，曲 ( ) () n Lx( )L x 1 ( )( ) N jijii i ca ac xx 1 M iij j aa ( ) i c x 率有以下几种选择： 最大曲率： 0 max( ) MC ii l ch l jj xxjM 2) 计算函数，第 条射线的累加和和曲率( ) ii h l jjj xxC 通过该像素的所有射线对该像素进行修正，从而完成第一轮迭代。 6）以上一轮迭代的结果作为初值，重复2)到5)进行新一轮迭代，直到取得符合收 敛要求的结果为止。 3.4 有序子集的统计重建方法 Hudson 和 Larkin 提出有序子集方法的核心思想是对投影数据进行子集划分，在每 次迭代的过程中只使用一个子集的投影数据而不是全部的投影数据。该方法大大提高 了 ML-EM 重建速度，同样 OS 也可应用到 MAP 类的算法中来提高重建速度。尽管 OS 算法能够加快重建速度，但是 OS 重建的图像常常不能收敛，随着迭代次数的增加， 图像的噪声大大增加。解决 OS 发散的一个方法就是使用松弛参数即随迭代次数增加而 逐渐减小步长。OS 方法还可以和其它算法进行结合派生出新的“OS”型算法5455， 如 OS-SPS。目前许多研究者作了大量的研究。这些研究有最优有序子集选取方法的讨 论，OS 收敛性必须满足的条件，有序子集平衡和一些推广算法，例如 BI 迭代算法 RA（row action）算法等。这些推广算法均具有良好的收敛性。 OS方法是把投影数据按照某种选择规则分成个有序子集，子集水平(level)表示S 子集的数目。原算法中的一步计算被分为步，每一步是对其中的一个子集内的投影S 数据进行重建图像。这样所有的子集都对像素校正一次，即图像更新次，完成一次S 19 迭代。以往的算法是所有的投影数据计算完一次之后，其像素值更新一次。在OS方法 中，所有的的投影数据计算完一次之后，其像素值已更新次。S 3.4.1有序子集ML-EM OS-EM 算法是将 OS 方法用于 EM 算法中，将投影数据分成个有顺序 12 (,.) N y yyS 的子集（） 。接着在每个子集内使用标准的 EM 算法，其重建结果为下一个子 12 ,. S T TT 集的重建初始值，在第个子集对像素校正完后，则完成一次迭代。上次重建的结果S 为下次迭代的初值。在 EM 算法中，图像的修正值是用所有的投影数据计算得到，而 在 OS-EM 中，图像的修正值是用每个子集内的投影数据计算。则 OS-EM 的迭代公式 为： ( , ) (1,1)(1,1) ( , ) 1 ,1,2.,0 s s s n T jijisnsn jjTM s n i ijj ij j i xa y xjM x a xa （3.23） OS-EM 算法的实现步骤： 1) 设定图像的初始估计值: (0)( 1,2.,); jj xxjM 2) 对第 个子集s a) 计算所有射线方程的估计值 (0) 1 M iijj j la x (1,2,)iN b) 计算 i i i y l (1,2,)iN c) 计算第个像素的修正值 j 1 1 1 N jiijN i ij i Ca a d) 对的值进行修正 j x (0) ;jjj xxC 3) 将重建出的代入下一个子集中，重复步骤(a)到(d)，这样循环直到做完第个 j xS 20 子集，这就完成了第一轮迭代。 4）以上一轮迭代的结果作为初值，重复(2)到(3)进行新一轮迭代，直到有符合收 敛要求的结果为止。 子集的选择一般应该遵循以下准则： 1）每个子集内的图像信息和统计信息应该最大； 2）子集内的像素对每个子集的活动贡献应该相同； 3