CT统计重建论文

基于最大似然和罚似然估计的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 is n’t 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

第一章绪论

1.1课题研究背景及意义 (1)

1.2国内外研究现状 (1)

1.3论文的主要工作及结构安排 (4)

第二章CT统计重建的基础

2.1CT成像理论 (5)

2.2统计重建的数学模型 (7)

2.3优化变换原理 (8)

2.4图像重建算法评价 (9)

第三章基于最大似然准则的CT统计重建算法

3.1最大似然准则的基本理论 (11)

3.2ML-EM算法 (12)

3.3ML-SPS算法 (14)

3.4有序子集的统计重建方法 (17)

3.4.1有序子集ML-EM (17)

3.4.2有序子集ML-SPS (19)

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.2OSL-EM算法 (31)

4.3OS-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 自适应选择正则项参数的实验结果分析 (39)

第五章锥束轨迹下有序子集罚似然统计(PL)重建

5.1锥束扫描下的投影矩阵的生成 (42)

5.2锥束下的OS-OSL-EM算法重建结果及分析 (44)

第六章总结与展望

6.1研究主要内容及成果 (47)

6.2存在的问题及以后的工作展望 (48)

参考文献

攻读硕士学位期间发表的论文

致谢

第一章绪论

1.1 课题研究背景及意义

计算机层析成像技术是利用有某种能量的射线源对物体进行断层扫描，根据所获得的某种物理量(指物质对射线的衰减系数)，运用特定的重建算法，重建出物体某个断层的无影像重叠的二维图像。自1967年CT(Computed Tomography)在英国问世以来，技术迅速发展，被公认为是20世纪影响人类发展的十大技术之一，其应用涉及到医学、工业、石油、物探、材料、生物、安全等众多领域。在国防工业中，CT已成为航天器、固体火箭发动机、高新技术等产品关键部件检测的必不可少的工具。由于CT技术是发展国防的基础技术之一，被发达国家列为敏感技术，作为国家机密进行严格的控制。因此，CT技术被视为国家科技实力的标志之一。

从数学上说，CT是从投影到重建图像的反问题[1]，具有普遍性，在数学界已经引起了广泛的重视。重建算法是工业CT技术中较为关键的一部分，CT重建算法主要可以分为滤波反投影重建算法和迭代重建算法。当投影数据完备、噪声不很严重时，解析法(滤波反投影(FBP)[2][3]等)可以得到很好的重建图像，但当CT系统在数据采集过程中受到严重的噪声干扰，或者投影数据采集不完全时，用解析重建算法得到的图像有伪迹。而迭代重建算法(统计迭代算法和代数迭代算法等)可以用来重建不完全数据、动态数据和噪声数据。其中统计迭代重建方法的优势是其对噪声的鲁棒性，通过加入先验模型和统计规律可以有效的抑制噪声，提高图像质量。

本课题主要研究统计算法在CT重建中的应用，通过研究有序子集的方法来提高收敛速度，加入惩罚项，约束项来降低噪声的影响，提高图像质量，为国内三维工业CT 的研制提供理论基础和技术支持。

1.2 国内外研究现状

近十几年来，CT广泛的应用在医学和工业上，相应的CT算法也迅速发展。其中解析算法重建速度较快，易于实现，是目前CT图像重建技术中最常用的算法。但是该算法通常要求采样数据是完全的，而且对噪声比较敏感，这在一定程度上影响了重建图像的质量，限制了它在实际中的应用。相对解析算法、代数迭代算法而言，统计迭代算法能准确描述投影数据的物理模型，对误差不敏感，易于加入约束条件等优点逐渐引起人们的重视，然而它与代数重建算法[4]一样重建速度慢，运算时间长，这些缺点极大地限制了它的应用，但是随着加速算法和计算机硬件技术的发展，统计重建将在实际中得到了广泛地应用。

一个完整的统计重建框架包括①图像的离散化②系统的物理几何模型③测量的统计模型(泊松分布)④优化准则(最大似然准则和罚似然准则)⑤求解方法。第①步是考虑如何对重建的图像进行建模，将待重建的图像离散成二维或者三维的像素矩阵，第②步考虑如何确定投影矩阵{}ij

=，A与系统的几何结构和建模方法相关[5-10]，第③步确

A a

定投影数据的模型，在统计重建中，一般认为投影数据服从泊松分布。第④步是引入一定的优化准则。第⑤步在④的基础之上进行问题的求解，该问题的求解实际上是一个参数估计的数学优化问题。

统计重建的优化方法中涉及到优化变换(optimization transfer，OT)，增量方法(incremental method)，有序子集加速(ordered subset，OS)，这三种优化方法是统计重建方法发展的主线[11]。

在统计重建中，最大似然估计(Maximum Likelihood，ML)是一种常用方法。1977年，A. P. Dempster 等人把期望值最大化(Expectation Maximization, EM)算法引入图像重建之中。1982年，Shepp[12]和Vardi[13]首次把期望最大迭代算法[14]应用到CT算法中，使似然函数最大化，之后ML-EM算法得到了广泛的应用。1997年，Fessler[15]等人提出了基于替代函数(surrogate function)的优化理论，其基本思想是在M步(maximization step)直接对E 步(expectation step)中所得到目标函数的替代函数作最大化处理，这一突破性的改进在很大程度上应当归功于De Pierro的凸性(convexity)算法[16][17]。当然替代函数的选取需要满足一定的条件，其它替代函数算法还包括Fessler[18-21]，Zheng[22]提出的算法以及全局收敛的增量优化传递算法。

1994年，Hudson和Larkin提出了有序子集(Ordered Subsets) [23]EM算法(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)[27][28]是一种贝叶斯重建算法,也可以称为罚似然重建(Penalized Likelihood，PL)[29-32]，该算法和ML重建算法可以说是一对关系亲密的姊妹，因为ML算法是寻求合适的图像的估计值使得图像得到的已知投影的概率最大；而贝叶斯重建算法是从已知的投影出发，要求在给定投影情况下所求图像的概率最大，通过选择合适的先验分布模型来提高重建结果的质量，其作用等同于优化理论中的正则项（惩罚项）。正则化技术最早由Tikhonov在著作中提出，他研究了求解不适定问题稳定解的基本理论，并提出了著名的Tikhonov正则化技术。先验模型的引入虽然能改善重建效果，但同时增加了估计的难度。就泊松观测模型而言，它本身属于非线性模型，因此它的ML估计和MAP估计均没有闭型解析式。针对这一问题，Herbert最先提出GEM(generalized EM)算法[33]，Green给出了经典的OSL(one-step-late)算法[34]，Fessler[35]提出了罚似然的SAGE(space-alternating generalized EM)算法，另外De Pierro的[36][37]凸性算法也涉及到了最大后验估计的问题。在以后的工作中，De Pierro 和Yamagishi[38]将RA(row-action)思想与最大后验估计的思想相结合，提出了BSREM(block sequential regularizedexpectation maximization)算法。Chung Chan, Roger Fulton[39]等人提出了自适应结构先验在CT中的应用，根据图像结构块中灰度的不同自适应的选择滤波。Xiaochuan Pan[40]阐述了最小化TV正则项(total variation)在图像重建中的应用。Jing Wang, Tianfang Li[41]等人把基于最小二乘的罚函数方法用于投影数据的去噪，提高了图像重建的质量。

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算进行比较分析。

第六章对本课题的研究内容进行了总结，指出本文的创新之处和解决的问题，并提出了尚未解决的问题，指明了今后研究的重点。

第二章CT 统计重建的基础

在过去的几十年里，统计迭代法一直是断层重建研究的热点，其重建图像的过程就是利用参数估计，优化变换原理等来逼近最优解一个过程。在重建过程中，首先确定能够正确反映成像系统的物理模型和统计模型以及输入输出关系，然后利用已知的条件确定求解方法。而且统计重建可以利用与物体相关的约束条件和其他的与特定系统相关的边界条件，利用精确的物理模型和适当的统计模型，可以在数据缺失和不规则采样的情况下重建出满意的图像。

2.1 CT 成像理论

CT 图像重建的原理是用穿透力强的射线扫描被测物体，当一束强度大致均匀的射线投照到物体上时，射线一部分被吸收和散射，另一部分透过物体沿原方向传播。由于物体的各种结构在密度、厚度等方面存在不同，它们对射线的吸收量也不相同，从而使透过物体的射线强度分布发生变化，通过CT 采集、转换得到投影数据，并运用相应的数学模型，进行图像重建出二维或三维图像，最后通过算法将射线强度分布转换成图像上的灰度分布，形成人眼可视的图像。

取一理想的X 射线源，它发出极细的笔束X 射线，在被测物体对面置一探测器，如图2.1(a)所示。假设强度为0I 的X 射线穿过均匀分布衰减系数为u 的物体，行进了x 的距离，强度变为I ，按Beer 定理有

00ln(/)ux i i e ux I I -==或

（2.1）

如图2.1(b)所示，若物体是分段均匀的，衰减系数分别是1,23,...u u u ,相应的长度为1,23,...x x x 则下式成立：

1230...ln(/)ux ux ux I I +++= （2.2）

更一般，若物体在xy 平面内为不均匀介质，则在某一方向l ，射线沿某一路径L 的总衰减值为

{}

0(,)ln(/)()exp (,)L

x y dl I

I I L I x y dl μμ=?=-?? （2.3）

其中(,)x y μμ=表示物体在二维平面上的衰减系数。()I L 是检测到的强度；L 是射线经过的直线，dl 是直线的弧长。

若没有具体路径，只说了沿着某一个方向，那么dl ?μ是投影，由等式(2.3)可知，测得0I 与I ，即可知道dl ?μ。在CT 中，入射的射线强度与出射的射线强度之比经过取对数运算后，则成了沿射线路径上衰减系数的线积分。CT 重建问题实际上是给定一个待重建物体的被测量的投影即线积分dl ?μ，然后运用算法去计算物体的衰减分布即被积函数μ。

（a ）均匀物体（b ）非均匀物体

图2.1 单能X 射线束在物体中的衰减图示

1917年，丹麦数学家雷当(J ．Radon)研究已经为CT 技术建立了数学基础，从数学理论上证明了某种物理参量的二维分布函数，由该函数在其定义域内的所有线积分确定。该结论指出了如果没有无穷多个且积分路径互不完全重叠的线积分，只能确定实际分布的一个估计近似值。

有了上述数学、物理基础后，为了实现工程技术上的应用，还需要解决两个主要问题。首先是如何采集到检测断层衰减系数线积分的数据集，其次是如何充分利用该数据集确定出衰减系数的二维分布。解决第一个问题可采用不同的扫描方式，即用射线束穿过被测体所检测断层并相应进行射线强度测量的规律性，可采用各种不同的扫描检测模式围绕提高扫描检测效率。解决第二个问题则是选择合适的图像重建算法，一类是解析法，以滤波反投影算法为代表。另一类是迭代重建法，有代数重建和统计重建。统计重

建是本论文研究的主要方向，下面具体介绍该算法。

2.2 统计重建的数学模型

在CT 图像重建中，图像区域是重建区域，首先应当明确，图像区域是一个正方形，其中心在坐标原点；图像函数是一个二元函数，其值在图像区域之外为零。一个具有n 个元素的网格把这个图像区域分成2n 个相等的正方形，每个小正方形(像素)内图像函数的值是均匀相等的。一幅图像的n n ?数字化是这样一幅n n ?的数字化图像，使得原图像在2n 个元素网格的任何像素上的积分值等于在同一个像素上的数字化积分值。在某一点的图像密度值就是在该点物质组织在某种能量下的相对衰减系数值，将物体的重建区域离散化为一幅数字化图像后，直接从投影数据来估计数字化图像每个像素点的密度。

图2.2 CT 数据采集

如图2.2所示，将待重建物体划分成M 个互不重叠的子区域，像素序列依次从1到M ，射线源和探测器是点状的，它们之间的连线对应一条射线。

核物理研究表明，X 光子的辐射满足泊松随机过程，结合这一统计特性，在CT 中选用泊松统计模型，可以表示成如下：

{}[]i i i y poisson A r +x :(1,2...)i N = （2.4）

其中'12(,,...)M x x x =x 表示重建图像矢量，用12(,...,)'N y y y =y 表示随机变量的N 个测量值,i r 表示第i 个投影的测量误差，{}ij A a =为投影矩阵，ij a 表示第i 个射线穿过第j 个像

素的长度。1[]M

i ij j j A a x ==∑x 表示第i 条射线的投影值。

2.3 优化变换原理

优化变换[42](optimization transfer ，简称OT)的思想是将复杂的统计重建的优化问题转换成易于求解的优化问题，转换后的问题一般具有易于优化、分离变量、降低维数的特点。

OT 的优化问题是寻求目标函数()Φx 的最大值，表达式如下：

?arg max ()x D

∈=Φx

x （2.5）由于该目标函数的最优解比较难求，可以选择一个用来逼近目标函数极大值的替代函数()(,)n φx x ，通过不断的迭代来求解使替代函数取得极大值。之所以要用迭代的方法来求解最小值问题是由替代函数的特性决定的，替代函数可能为非二次曲线，不易求最大值，计算量非常大，这些都给求解带来了困难，所以一种可行的方法就是选择一个可以用来逐步逼近替代函数极大值的替代函数，通过不断的迭代计算来得到替代函数的极值。但是替代函数一定要满足一定的条件才可以得到稳定的结果，首先替代函数的目的就是为了简化求解极大点的计算量，如果替代函数是一个多维的非二次函数，就难于计算。那么替代函数就要尽量将多维的求解极值问题变为一维的问题来简化求解过程。另外为了保证收敛的一致性，如图2.3所示，要求替代函数在当前点与替代函数相切并且位于替代函数的下部，则替代函数满足的条件如下：

()()()(;)()n n n φ=Φx x x （2.6） ()(;)()n φ≤Φx x x （2.7）

由条件(2.6)，(2.7)可以得出

()()()()()()(;)(;)n n n n φφΦ-Φ≥-x x x x x x （2.8）

换句话说，我们选择x 使得()()()(;)(;)n n n φφ≥x x x x 来使得()()()n Φ≥Φx x

图2.3 替代函数曲线图

其中()(,)n φx x 为替代函数，()n x 为第n 次迭代后的结果。只有替代函数满足了这上述条件才能在简化计算的同时保证得到稳定一致的收敛结果。

选定了替代函数以后，求解最大值的过程就可以简化成式(2.9)所示的迭代过程，每次迭代都是以上次的结果为一点，得到在这一点上的替代函数并求解极大点作为本次迭代的结果。则该优化问题的迭代过程可以简化成：

1()0

arg max (,)n n x φ+≥=x x x （2.9）

2.4图像重建算法评价

对于重建算法的可行性,一般按照下面三个方面来评价: 1) 算法是否易于实现。

2) 计算量的大小是衡量算法可行性的一个重要指标。计算时间是一个相对量, 在同一计算机上，可以通过几种算法计算时间来衡量。

3) 迭代算法的误差分析。用归一化均方差距离、归一化平均绝对距离、相关系数、信噪比等评价指标反映重建图像的误差和噪声。

相关系数CC(Correlation Coefficient)衡量：

11/2

2211()()

()()P

rec rec ref ref j j j j j P

rec rec ref ref j j j j

j j x

x x x CC x x x x ===--=

??--????

∑∑∑ （2.10）

归一化平均绝对距离NMAE (Normalized Mean Absolute Error)来衡量： 11

()/()M

rec

ref

ref j

j j j MAE x x x ===-∑∑ （2.11）

信噪比SNR (Signal Noise Ratio)： 2

211

10log(()/())M

ref

rec ref j

j j j j j SNR x x x x ===--∑∑ （2.12）

归一化均方差距离NMSE(Normalized Mean Square Error):

()()1

222

r e c

r e

f j j

ref ref j j i x x MSE x x ??

- ?= ?- ???

∑∑ （2.13）

ref j x ，rec j x 分别为被检测图像和重建图像对应像素的值，ref j x 为被检测图像的平均灰度值，rec j x 为重建图像的平均灰度值。相关系数反映了被测试图像和重建图像之间的相似程度，相关系数越大，NMAE 越小说明重建效果越好。用信噪比SNR 来评价图像的噪声的大小，SNR 越大，图像所含的噪声越小。NMSE 越小表示图像与原始图像的误差越小，图像质量越好。

第三章基于最大似然准则的CT 统计重建算法

最大似然估计考虑了数据的统计特性，能够充分利用系统的固有分辨率，从而重建的图像分辨率和噪声特性均优于卷积反投影重建的重建结果[43]。但是在图像重建优化求解的过程中，似然函数的解析表达式很难找到，因而限制了ML 的使用。但是期望最大法通过构造“完全数据集”的似然函数，充分利用完全数据似然函数的简单性去求解实际投影数据的最大似然估计，从而解决了这个难题。

3.1 最大似然准则的基本理论

图像重建的极大似然模型是建立在假定（X 射线发射的光子数是服从泊松（Poisson ）分布）的基础上的。最大似然准则为：假设y 是相互独立的，已知y 的条件下寻求x 使下面的对数似然函数最大

()()l n L P =x y x （3.1）基于Poisson 分布的似然函数为()P y x

()[][]1

()()!

y N

Ax i i i A S P e

y -===∏x x y x （3.2）

假定ij a 满足以下两个条件

01ij a ≤≤ （3.3）

1,M

i a

==∑(1,2,,)j M = （3.4）

则(3.2)可以简化得到下式：

()()!

ij j

j M

y ij j a x N

j i i a x S e

y =-

==∑

=∑∏x （3.5）

图像重建的极大似然方法就是求出图像矢量()12,,,T

M x x x =x L ，使似然函数最大，也就

是使对数似然函数()()ln L P =x y x 最大。因此，ML 估计为：

?a r g m a x ()

ML L =x x [][]{}1arg max ln()ln(!)N

i i i i i A y A y ==-+-∑x x

arg max(log()log(!)M N N M N

j ij i ij j i j i i j i x a y a x y ======-+-∑∑∑∑∑ （3.6）

另外，还可以根据判定极值的二阶充分条件（即目标函数的Hessian 矩阵是半正定的）来说明似然函数的凸性，从而可知ML 估计(3.6)存在且唯一。

3.2 ML-EM 算法

目标(3.6)的优化具有一定的难度，由于目标函数是非线性的，不存在解的闭型表达式；而且在求解过程中隐含着非负性约束条件，为了解决该问题，Shepp [44]和Vardi [45]建议使用EM 算法对(3.6)进行优化。1994年，EM 算法被成功地应用在图像重建中，成为研究不完全的数据极大似然估计的一种方法。EM 算法的收敛解是非负的，迭代形式便于计算机实现，在图像重建算法研究中，EM 算法的应用相当重要。

EM 算法认为直接被观测到的投影数据y 总是不完全(incomplete)的，依赖一个较为完全的数据R ，它们之间的关系由映射χ表征，即()y z χ=，一般x 为已知。Z 称完全数据，它包含了待求未知参数x 的所有信息。基本思想就是通过重复处理完全数据的统计问题来解决不完全数据的统计问题。首先，假设初始图像为(0)x ；接着根据(0)x 求一次近似图像(1)x ，再由(1)x 求二次近似图像(2)x ，以此类推到第n 次迭代的图像()n x 。

根据ML 估计，优化的目标函数为()()ln P Φ=x y x ，利用优化变换原理，需要寻找一个替代函数()(,)n EM φx x ，其必须满足(2.9)中的单调性的条件。EM 算法在理论上总能以概率收敛到某个的ML 估计[46]。

ML-EM 算法的迭代过程分两步：

E 步：在给定当前图像估计()n x 及测量得到投影值y 的情况下，完全数据对应的对数

似然函数的条件期望为，

()()

(,)log()),n n EM E Z φ??=??

x x x y x （3.7） M 步：最大化替代函数对图像向量x 进行估计

1()0

a r g m a x (,)

n n EM x φ+≥=x x x （3.8）根据计算步骤，需要找出替代函数，推导过程中应用了算法中的凹性特征，假设0i r >，当()

0n >x

时，[]1

()M

i ij j i j l A a x ===∑x x ，则有

111()log()()N

M M ij j i ij j i i j j L a x r a x r ===??

=+-+????

∑∑∑x

()()()

111log(()()()()()n N

M M ij j j n n i i i i ij j i n n n i j j i j i a x x r y l l a x r l x l ===??=+-+??????

∑∑∑x x x x g g （3.9）

利用凹性特征有：

()()()()()()

111()log(()log ()()()()n N

M M ij j j n n i i i i ij j i n n n i j j i j i a x x r L y l x l x a x r l x x l x ===??????≥+-+????????????

∑∑∑x g g （3.10）忽略(3.9)式中的常数项，即可得到替代函数()(,)n EM φx x

()()

()

11(,)()log ()n M

N ij j n EM i j j j n j i i a x y x a x l x φ==????=-??????

∑∑x x g （3.11）其中1

j i j i a a ==∑，根据M 步对(3.11)最大化，求关于任意j x 的偏导数并令其为零即

()(,)

0n EM j

x φ?=?x x ，可以立即得到：

()()

0n N

j i ij

ij M

n i i j

ij j

j x y a a x a x

===-+

=∑∑

∑

（3.12）即推到出ML-EM 算法的迭代公式：

()

(1)

(1)()

,1,2...,0n N

j ij i

n n j

j N

n i ij

ij j

i j x a y x j M x a

a x

++====

=>∑

∑∑ （3.13）

根据迭代公式看出EM 算法有两个明显的特点：

1）每一次迭代后均提高后验密度函数值，即有(1)()()()n n l l +≥x x 。

2）给定非负的初始估计(0)x ，以后迭代所得的(1)(2),,x x 也都是非负的，满足了图像像素的非负性约束。

ML-EM 算法的具体步骤如下：

1) 设定图像的初始估计值: (0)(1,2...,);

j j x x j M == 2) 计算所有射线方程的估计值 (0)1M

i i j j j l a x ==∑(1,2

,,i N = 3) 计算 i

i i y l ?=

(1,2,,)i N = 4) 计算第j 个像素的修正值 1

j i

i j

i ij

i C a

===

?∑∑

5) 对j x 的值进行修正 ;)0(j j j C x x ?=

通过该像素的所有射线对该像素进行修正，从而完成第一轮迭代。

6) 将以上一轮迭代的结果作为初值，重复2)到5)可得第n 轮结果，就得到一序列

(0)(1)(2)(),,,,,n x x x x ?????若符合收敛要求,即对事先给定的很小的正数ε存在正整数NN ，使当n NN >时，有 ()(1),(1,2,,)n n j j x x j J ε--<=???.

3.3 ML-SPS 算法

SPS(separable paraboloidal surrogates)算法是全局收敛的，并且满足像素的非负性条件，迭代的初始收敛速度要比SAGE 快。SPS 适用于ECT 的ML 具有Convex 形式先验的

MAP 问题[47-51]，Erdogan [52]将其推广到TCT 的ML 和MAP 问题。该算法是从优化变换的思想出发，通过构造易于最大化的二次型替代函数简化EM 算法的计算复杂度。

首先定义非负集合：{}:0,M

M j R R x j +

=∈≥?x 定义集合{}

:[]M D R A i r i =

∈>-+x x

根据最大似然估计理论，最大似然SPS 算法的目标函数表达式如下：

1()()([])N

i i i L h A =Φ==∑x x x

（3.14）其中对数似然函数的表达式可以为：

()log()()i i i i h l y l r l r =+-+, 其中[]1()M

ij j j l A a x ===∑x x （3.15）

函数()i h l 有如下性质：

（1）i l r ?≥-，?()()i i h l h l ≤，其中?i i l y r =-

（2）在区间?,i i r l ??-??上，()i h l 单调递增，在区间?,i l ??∞??上单调递减。

（3）()i h l 是凹函数

在ML-SPS 算法中，首先找到一维函数()i h l 的抛物面型替代函数()(,)p n i q l l ，替代函数需要满足的条件：在当前迭代点()[]n n i i l A =x 处，它与函数()i h l 相切且在函数()i h l 的下方。通过对()i h l 做二阶Taylor 展开：

()()()()()()21(,)()()()()()2

p n n n n n n i i i i i

i i i i q l l h l h l l l c l l l =+---& （3.16）其中()h l &为函数()h l 一阶微分。把公式（3.16）带入公式（3.14）中可以得到函数

()

()1(;)([];[])N

n p n i i i Q q A A ==∑x x x x （3.17）

利用()(;)n Q x x 的凹性，可以得到()L x 的可分离的二次替代函数()(,)p n q x x

{}()()'()()()'()()1

(,)()()()()()()2

p n n n n n n n j q diag c L L =---+?-+x x x x x x x x x x x （3.18）