本发明专利技术的ＰＥＴ成像中卡尔曼滤波图像重建方法，通过ＰＥＴ正电子发射断层扫描仪得到原始投影线的正弦图，然后建立状态空间体系，通过基于状态空间的卡尔曼滤波法得出放射性活度分布，重建图像。由于本发明专利技术利用基于状态空间体系的卡尔曼滤波法重建ＰＥＴ图像，有效地提高了重建图像的质量；通过与现有重建方法的实验比较，重建结果的偏差和方差都有了很大地提高。
【技术实现步骤摘要】
本专利技术涉及一种医学断层影像重建方法，尤其是涉及一种在PET成像 中的卡尔曼滤波图像重建方法。
技术介绍
正电子发射断层成像(positron emission tomography, PET)作为一种 生物医学研究技术和临床诊断手段是核医学成像装置中最为重要的应用 之一。PET的基本思想向生物体内部注射正电子同位素标记的化合物， 在体外通过图像重建技术测量它们的空间及时间分布。如今许多临床医学 领域已经开始广泛使用PET图像进行肿瘤、心脏疾病(心肌存活的检测、 心脏移植、介入治疗前后监测等)、神经和精神系统疾病(如癫痫病灶的 确定、帕金森病、老年痴呆和精神病等)的诊断，部分发达国家已经把部 分疾病的PET检查列入到医疗保险范围。PET探测人体内发出的放射性信号，经过符合和采集系统处理，形成 投影线，并以正弦图(Sinogram图)方式存放于计算机硬盘中。计算机以 Sinogram图为输入，调用图像重建模块，计算得出人体横切断层图像，反 映体内某组织器官的新陈代谢情况，其实这类重建逆问题都可以用方程 y(0:i^(0 + v(f)来表示，r为PET的测量数据，X是要重构的图像,D是 系统模型。基于含有噪声的数据Y估计X的方法，可以划分为两类解 析法和统计迭代法。相对于解析法而言，统计迭代算法包括基于泊私、漠型 的最大似然法、最大后验概率法和基于高斯模型的最小二乘法，是以统计 规律为&出，对不完全数据的适应性好，能够得到更精确的结果，也因此 受到研究人员的广泛关注。迭代法是从一幅假设的初始图像出发，采用逐步逼近的方法，将理论 投影值同实际测量投影值进行比较，在某种最优化准则指导下寻找最优 解。在PET中常用的迭代法包括最大似然法(Maximum Likelihood Expectation-Maximization, MLEM)和OSEM (有序子集最大似然法，4Ordered Subset Expectation Maximization)算法。迭^f^、法^尤,泉之一是可以才艮 据具体成像条件引入与空间几何有关的或与测量值大小有关的约束条件， 如可进行空间分辨不均匀性的校正，物体几何形状约束，平滑性约束等控 制迭代的操作，在某些场合下，比如在相对欠采样、低计数的核医学成像 中可发挥其高分辨的优势。迭代法最大的缺点是计算量大，计算速度慢。迭代法求解过程是a. 假定一初始图Y象；b. 计算该图像投影；c. 同测量投影值对比；d. 计算校正系数并更新初始图像值；e. 满足停步规则时，迭代中止，否则以新的重建图像作为初始图像 从b步开始。OSEM是近年来发展完善的快速迭代重建算法，每次重建时只使用一 个子集对投影数据进行校正，重建图像更新一次，这样所有的子集都对投 影数据校正一次，称为一次迭代。和传统的迭代算法MLEM相比，大大 加快了图像重建速度，缩短了重建时间。
技术实现思路
本专利技术提供一种正电子断层扫描中的基于状态空间理论的卡尔曼滤 波图像重建方法，大幅度提高PET重建图像的质量。本专利技术的卡尔曼滤波图像重建方法，包括以下步骤 (1)输入原始投影线的正弦，=Z)X(/) + v(/))建立状态空间体系;^+1) = ^^)+<《其中，t表示时间；Y就是正弦图数据；D为系统矩阵，表示发射的光 子被探测器接收到的概率；A是状态转移矩阵，在瞬时稳定状态下是一个 单位矩阵；X为》史射性浓度状态变量，为需要重建的对象；w是过程噪声， 服从正态高斯分布w(0 W(0，2(O); v为测量噪声是各种噪声的集合， 服从正态高斯分布v(,) ~ W(O，A(O);(3)利用卡尔曼滤波，根据下列方程得出放射性活度分布，重建图像5<formula>formula see original document page 6</formula>其中，z(0为测量值，A为浓度初始值，？。为初始浓度误差协方差，迭代从初始^。，尸。，込，A出发，通过测量值z(O，不断修正估计出来的浓度值I,最终给出放射性活度分布，重建图像。上述的正弦图(即Sinogram图)通过PET正电子发射断层扫描仪得到。PET正电子发射断层扫描仪进行透射扫描和发射扫描，透射扫描得到图像的衰减校正系数，发射扫描在180度内选取48个采样角度和34个径向单元。在上述步骤(3)中，具体地迭代重建包括以下步骤1) 首先设定状态的初始值和初始协方差义。,尸。；2) 利用时间更新方程方程 和②及时向前推算出当前状态和误差协方差的先验估计值；3) 利用离散采集时间点测得的观测值Z(O ，根据状态更新方程 、 来更新先验值，得到这一时刻的最优估计；4) 重复交替步骤2 )和步骤3 )直至获得最优重建结果。本专利技术的优点是通过基于状态空间体系的卡尔曼滤波法重建PET图像，有效地提高了重建图像的质量；通过与现有重建方法的实验比较，重建结果的偏差和方差都有了很大地提高。附图说明附图l是扫描图像数据；附图2a是MLEM法对附图1的重建结果示意附图2b是本专利技术对附图1的重建结果示意附图3是模型中间一列像素对应真值的放射性浓度分布值示意图；附图4是体模示意附图5是MLEM法对附图4的重建结果示意图；附图6是本专利技术对附图4的重建结果示意图。具体实施例方式正电子发射断层扫描仪探测人体内发出的放射性信号，经过符合和采集系统处理，形成投影线，并以Sinogram图(正弦图)方式存放于计算机硬盘中。计算机以Sinogram图为输入，调用图像重建模块，计算得出人体横切断层图像。在应用PET正电子发射断层扫描仪时，进行透射扫描和发射扫描。透射扫描得到图像的衰减校正系数。发射扫描在180度内选取48个采样角度和34个径向单元，系统矩阵由Fessler教4受研究小组提供的软件包生成。对原始采集到的Sinogram图数据进行各类校正，从而建立一个离散的测量方程其中，Y是Sinogram图数据，D是系统矩阵，v表示观测噪声。建立发展方程人体内的同位素分布在成像过程中暂时稳定时有》(0 = 0,通过对发展方程的离散化，并且加入过程噪声，可得到更一般的表达邵+ 1) = ^￥(/) +《)其中，A是状态转移矩阵，在同位素分布暂时稳定时是一个单位矩阵。构建状态空间体系把离散的测量方程和发展方程联合起来，建立一个状态空间体系r(f)二/)邻)+v(0其中，v服从正态高斯分布v(0 W(0^W)， w服从正态高斯分布iV(0,卯))。利用卡尔曼滤波，根据下列方程重建图像邵+①〖(r +1) = P(r +1, f)Z/ [D尸(f+1， f )Z/ +及)r1 ③Z(1) = Z(, + 1)-④A^ + l)=1,0 ++1)2(^ + 1) ⑤尸(/+1)=[/—〖0+1)]尸0+1,0 ⑥其中，z(0为测量值，1。为浓度初始值，《为初始浓度误差协方差，迭代从初始A，A，，A出发，通过测量值z(O,不断修正估计出来的浓度值义，最终得出放射性活度分布，重建图像。采用Kalman滤波进行图像重建时，主要分成四个步骤1) 首先设定状态的初始值和初始协方差％。，/>。；2) 利用时间更新方程①和②及时向前推算出当前状态和误差协方差的先验估计值；3) 利用离散采集时间本文档来自技高网...
【技术保护点】
一种ＰＥＴ成像中卡尔曼滤波图像的重建方法，其特征在于：包括以下步骤： （１）输入原始投影线的正弦图； （２）建立状态空间体系： Ｙ（ｔ）＝ＤＸ（ｔ）＋ｖ（ｔ） Ｘ（ｔ＋１）＝ＡＸ（ｔ）＋ｗ（ｔ） 其中，ｔ表示时 间；Ｙ就是正弦图数据；Ｄ为系统矩阵，表示发射的光子被探测器接收到的概率；Ａ是状态转移矩阵，在瞬时稳定状态下是一个单位矩阵；Ｘ为放射性浓度状态变量，为需要重建的对象；ｗ是过程噪声，服从正态高斯分布ω（ｔ）～Ｎ（０，Ｑ（ｔ））；ｖ为测量噪声是各种噪声的集合，服从正态高斯分布ｖ（ｔ）～Ｎ（０，Ｒ（ｔ））； （３）利用卡尔曼滤波，根据下列方程得出放射性活度分布，重建图像： Ｘ（ｔ＋１，ｔ）＝ＡＸ（ｔ） ① Ｐ（ｔ＋１，ｔ）＝ＡＰ（ｔ）Ａ↑［Ｔ］＋Ｑ（ｔ） ② Ｋ（ｔ＋１）＝Ｐ（ｔ＋１，ｔ）Ｄ↑［Ｔ］［ＤＰ（ｔ＋１，ｔ）Ｄ↑［Ｔ］＋Ｒ（ｔ）］↑［－１］ ③ ＊（ｔ＋１）＝Ｚ（ｔ＋１）－ＤＸ（ｔ＋１，ｔ） ④ Ｘ（ｔ＋１）＝Ｘ（ｔ＋１，ｔ）＋Ｋ（ｔ＋１）＊（ｔ＋１） ⑤ Ｐ（ｔ＋１）＝［Ｉ－Ｋ（ｔ＋１）Ｄ］Ｐ（ｔ＋１，ｔ） ⑥ 其中，Ｚ（ｔ）为测量值，Ｘ↓［０］为浓度初始值，Ｐ↓［０］为初始浓度误差协方差，迭代从初始Ｘ↓［０］，Ｐ↓［０］，Ｑ↓［０］，Ｒ↓［０］出发，通过测量值Ｚ（ｔ），不断修正估 计出来的浓度值Ｘ，最终得出放射性活度分布，重建图像。...
【技术特征摘要】
1、一种PET成像中卡尔曼滤波图像的重建方法，其特征在于包括以下步骤(1)输入原始投影线的正弦图；(2)建立状态空间体系Y(t)＝DX(t)+v(t)X(t+1)＝AX(t)+w(t)其中，t表示时间；Y就是正弦图数据；D为系统矩阵，表示发射的光子被探测器接收到的概率；A是状态转移矩阵，在瞬时稳定状态下是一个单位矩阵；X为放射性浓度状态变量，为需要重建的对象；w是过程噪声，服从正态高斯分布ω(t)～N(0，Q(t))；v为测量噪声是各种噪声的集合，服从正态高斯分布v(t)～N(0，R(t))；(3)利用卡尔曼滤波，根据下列方程得出放射性活度分布，重建图像X(t+1，t)＝AX(t) ①P(t+1，t)＝AP(t)AT+Q(t) ②K(t+1)＝P(t+1，t)DT[DP(t+1，t)DT+R(t)]-1 ③ ④ ⑤P(t+1)＝[I-K(t+1)D]P(t+1，t) ⑥其中，Z(t)为测...
【专利技术属性】
技术研发人员：沈云霞，刘华锋，施鹏程，
申请(专利权)人：刘华锋，
类型：发明
国别省市：86[中国|杭州]
还没有人留言评论。发表了对其他浏览者有用的留言会获得科技券。