MLEM算法全过程推导.pdf

1、写在前面 ： 本文 主要从一个初探着的角度，一步一步地推导出了 PET 图像重建 MLEM 算法的迭代公式。 认真读完这篇文章的朋友将 具有 自行推导 MLEM 算法的迭代公式 的能力 ，知道每一步推导的原因 。 推导过 程中 重要的 已 地方用红字标出。 希望本文对从事或正准备从事 PET 图像重建 工作 或对 PET 图像重建有兴趣 的 你有一定的帮助 _！ 一、 PET 图像重建 ML-EM 算法 从头到尾的推导 过程（ PET 统计迭代重建算法入门）： 先说说我对 PET 图像重建中 统计迭代 重建算法的理解。统计迭代重建算法实际上是对PET 的图像重建建立了一个统计学模型，然后用相应

2、的统计学原理或者统计学算法来解决或者优化最后得到结果，也就是图像。 而为什么强调 统计 二字，是因为正电子发射本身就是一种满足泊松分布的统计模型。因此采用统计模型对 PET 图像重建进行建模在 噪声、边界 上能更方便地加入一些信息来对其进行约束 如 MAP、 TV 等 。这就是 PET 图像重建往往选择统计迭代方法而不选择代数 解析 迭代 方法的原因。 而说到 PET 图像重建的统计学建模， ML-EM （ Maximum Likehood- Expectation Maximization）算法 中的统计学建模可以说是最经典，最被业界公认的一种建模方法。由于ML-EM 中对 PET 图像重建

3、统计学建模的优势，在这种建模方法的基础之上， 也 衍生 了许许多多的更适应实际的图像重建算法 ，如为了提高收敛速度采用有序子集 OSEM等方法重建 。 在 迭代算法中，通常将图像 表示为一维数组 X。 Xj 表示排列的第 j 个像素。 一） ML-EM 重建算法中的统计学建模 模型假设 1. 我们最后要得到的图像像素 值 f0j 表示对应像素位置的活度 （单位为 Beq 或 mCi） 。因为 （ 1）衰变的正电子发射本身是满足泊松分布的； （ 2）发射的正电子数的期望与活度 f0j 成正比； （ 3）默认正电子湮灭朝各个方向发射伽马光子的概率是相等的。 因此该像素位置上 在一定扫描时间内， 朝

4、各个方向发射的伽马光子 的总 数 ，记为 fj， fj被看成 满足 fjpoison(Xj)的随机变量， 其期望值 Xj=a*f0j， a 为 时间 。 2. 系统响应矩阵 Hij 则理解为，像素 f0j 内 发射的伽马光子被第 i 条 LOR 上探测器探测到的概率。 （ 1）每个像素之间正电子发射相互独立； （ 2）第 i 条 LOR（ line of response，响应线） 探测光子对 的期望值 E( Pi ) = (j)Hij*fj，的含义为 像素 f0j 发射光子总数 fj 中，被 第 i 条 LOR 上探测到的光子数； （ 3）由 1. fj 又是满足泊松分布的随机变量 ，因此

5、根据泊松分布的性质， 有 Pipoison( ( j)Hij*Xj )。 3. 系统响应矩阵式可以通过 Monte Carlo 仿真得到 ，或是通过一些 数学 模型计算。 它 是可以预先保存下来 ，或是在重建过程中 计算 的常量。 二） ML-EM 重建算法中的统计学建模 符号系统（仅限于本文） f0j： 通过一维数组排列的第 j 个像素的活度 （可以理解为正电子发射的速度） Xj： （ 把二维图像的像素排列成一维数组 ） ， 第 j 个像素朝各个方向发射的伽马光子 的总数 。是随机变量的 期望 与 f0j 成线性关系，也是最后要求的值。 fj： 第 j 个像素朝各个方向发射的 总 伽马光子

6、总 数 。 随机变量 Hij： 系统相应矩阵，表示第 j 个像素发射的伽马光子被第 i 条 LOR 探测到的概率 。常量 E( Pi )： 第 i 条 LOR 探测到的伽马光子数 。 期望 Pi： 第 i 条 LOR 探测到的伽马光子数 。 随机变量 ，由探测器探测得到 Zij： 第 i 条 LOR 上探测到的第 j 个像素发射出的光子数，随机变量 p( k | =0 )： 这里的 表示 期望为 0 的泊松分布在 变量取 k 的时候的概率 后文的 粗体 变量表示向量或数组。在我阅读的很多文献里，经常会把 随机变量 和 期望 的概念弄混淆。随机变量是不确定的，而期望 往往是取决于事物本身的性质，

7、 是一个确定的值。在此我特地把这两个概念区分开，用不同的符号表示，以便减少读者的疑惑。 三） ML（ Maximum Likehood）最大似然估计 最大似然估计 ML 本身是一个很简单的思想：对于 目标函数 为似然函数 g( f | X )，变量为 f，在参数 X 取何值的时候达到最大，通常通过 令偏导等于 0。在 PET 图像重建中： 1） f： 表示像素朝各个方向发射的伽马光子 总 数 随机变量 的向量组； 2） X： 表示像素朝各个方向发射的伽马光子数的 期望值 的向量组。 *PET 图像重建中的 ML* 在 PET 图像重建中。似然函数 g 是通过叫做“ 完全数据集 ”的概念来构造的

8、。即通常我们先得到一组随机变量 f 的值 ，然后根据 f 的独立性，将每个 fj 的似然函数相乘来构造整体似然函数： g( f | X ) = (j)g( fj | Xj )= (j)p( fj | =Xj ) 通过取对数求偏导 计算 X 取何值的时候， g( f | X )能取最大（或极大），这个 X 的值就是最大似然估计求得的结果。但是这个式子构造的“完全数据集“ f 我们并不能事先得知，因此不能利用 已知的 PET 探测器得到的 信息 P。 也就没法对 Xj 求极大。 因此我们要利用可探测到的 P 作为完全数据集的变量 ： E( Pi )= (j)Hij*fj ； Pipoison( (

9、j)Hij*Xj ) 所以我们调整 ”完全数据集“ 的对象为 P 似然函数为： g( P | X )= (i) p( Pi | =(j)Hij*Xj ) 因为泊松概 率乘性因子不好求导，所以取对数转换为加性因子，且 取对数不改变其 极值点 。 取对数化简后： ln |ln!lnj i j ji i j ij jPi HXj ij jiii j ij j i j ij jp P H XH X ePH X P H X co ns t const 是与 X 无关的量 看到在对数 ln 内还有“ “符号，对接下来的 ML 求偏道很不利。 *引入变量 Zij* 因此在 ML-EM 算法中， 引入完全观测

10、数据 Zij，它表示第 i 条 LOR 上探测到的第 j 个像素发射出的光子数。它也是服从泊松独立分布 Zijpoison( Hij*Xj )。 它满足 Pi=(j) Zij。 所以由 Zij 构成的完全数据集的似然函数为： g( Z | X )= (i) (j) p( Zij | =Hij*Xj ) 取对数化简后： g( Z | X )= ln |ln!lnij i j ji j ij ij jz HXij jijiji j ij j ij ij jp Z H XH X ezH X Z H X co ns t （ 3-1） 此式为一凸函数，最大值一定在极值点，即偏导为 0。特别指出 下式 ：

11、 lni j i j j i j i j jH X Z H X （ 此式极为 重要 ，很多文献跳过之前步骤直接写出此式 ） 再对 Xj 求偏导： i i ji i jjjZg HXX （ 3-2） 算到这一步似乎没法继续进行下去了 ， 偏导等于 0 后 求的 Xj 没什么 实际意义 ，因为 Zij这个变量我们不知道。我们只知道 Pi=(j) Zij 如何解决 Zij 的问题，在 EM 步中给出了解释。 四） EM（ Expectation Maximization）最大 化 期望值 1. EM 算法分为 两步 ： （假设式中有两个量 a、 b 未知） 通常先将一个量赋予一个初值，如令 b=1。

12、 步 1） Expectation： 求得 a 的期望用 E( a )。 步 2） Maximization：用 E( a )替代 a 来求得我们需要的最大化的值 b。 步 3）求得的值 b 又会继续影响 E( a )， 故 再重复步骤 1） 直至收敛。 2. 在 PET 图像重建 ML-EM 重建算法中，对应的 EM 步骤为： “E”步 1）：因为我们最后要通过 ML 得到所有像素朝各个方向发射的光子数的期望 X。 为了构造交替迭代，用 Zij 在当前所有像素朝各个方向发射的光子 总 数的期望 X 已知的情况下的 条件期望 代替随机变量 Zij 本身。即令 Zij=E( Zij | X 当前

13、 )。 即且已知 常量 Pi=(j)Zij=Zij+( j0 != j )Zij0。 可以求得 ： | kk ij jij ij i kj ij jHXZ E Z X PHX （ X 当前用 Xk 表示） 上式推倒过程如 下 （若读者不感兴趣可以直接看 “M“步 ）： 00|ii j i j iPjjkijki j ijP Z k P Z P kE Z X kP Z P 由泊松分布的性质知：00 ijjjZ也是满足泊松分布的， 为后文描述方便起见记： ij ij jHX i ij jj HX ，则上式等于： ! !|!ijjiik P kjjPikij Pkieek PkE Z X keP !

14、 !iiik P kjjPiPkik PkkP ， 再化简得 !ii k P kPjjik iPkk P k 1iiik P kPk jjPk kC 这是一个 二项分布 的求期望公式，可看出 n=Pi， p= j 。所以最后值等于： | k jijE Z X n p P i | k ij jij iij jjHXE Z X PHX “M”步 2）：将求得的 Zij 带入式（ 3-2）， 令偏导等于零 ： 0 0 00ki j ji ki j i j ji i jjHXPHXH X 0 0 0k ijji ki j ij ji ijjHXPHXHX0 0 01ii i j kj i j jkkj

15、 j ji i jPHHXX X XH 步 3）得到的更新图像 1kjX 又会在重复步 1）时作为 kjX 影响 E( Zij | X )。 综（一）（二）（三）（四）就得到了 ML-EM 重建算法的迭代公式： 0 0 01ii i j kj i j jkkjji i jPHHXXXH 小结 ： 本文 主要从一个初探着的角度，一步一步地推导出了 PET 图像重建 MLEM 算法的迭代公式。 认真读完这篇文章的朋友应该具有 自行推导 MLEM 算法的迭代公式 的能力 ，知道每一步推导的原因 。 至于 对 MLEM 迭代公式的理解，在本文就不详细叙述了，和其他迭代算法一样，大致分为 四 个步骤：投影、比较、反投影 、更新图像 。 关于 MLEM 算法的收敛性，本文 也 不做 详细 讨论 了，因为网上 已经 有很多 相关的博客和文章。 希望本文对从事或正准备从事 PET 图像重建 工作 或对 PET 图像重建有兴趣 的 你有一定的帮助！ 其实此文早在几个月前就写完了，但是过程比较粗糙，到现在才 想起把中间的过程做的更饱满、易懂， 但是我不是很拘谨与格式，把格式做好要浪费太多的时间了， 所以我在重要的地方用红字标出， 以便 让读者 看的清重点 。 全文由本人纯手打，若有笔误敬请谅解 若读者对本文的观点另有高见，欢迎指正！ 致谢 2017 年 12 月 5 日星期二