基于贝叶斯框架的三维磁粒子成像重建方法、系统和设备与流程

1.本发明属于磁粒子成像领域，具体涉及了一种基于贝叶斯框架的三维磁粒子成像重建方法、系统和设备。


背景技术：

2.在临床上，一般需要借助于影像手段对病灶进行准确、客观的诊断和检测，因此，如何实时的、准确的获得病灶的分布信息，一直是国际上的研究热点。磁粒子成像技术(magnetic particle imaging,mpi)是一种新兴的断层成像技术。由于其只能追踪磁性纳米粒子信号，且不受成像深度限制，因此，mpi成像技术具有高特异性、高灵敏度、高时空分辨率、且成像深度不受限等优点，是一种非常具备临床应用潜力的新型成像技术。
3.目前，大多数的mpi图像重建方法要借助于系统矩阵来完成。因为基于系统矩阵的重建方法可以提供更高的成像分辨率。然而，更高的成像分辨率意味着更高的系统矩阵维度，即更高的计算复杂度。在实际的临床应用中，如血流成像、血管灌注成像等，实时的成像效果能更便捷的引导治疗策略的制定。因此，快速的mpi重建方法显得尤为重要。


技术实现要素：

4.为了解决现有技术中的上述问题，即现有的技术随着分辨率提高带来的系统矩阵复杂度过高难以实时进行成像的问题，本发明提供了一种基于贝叶斯框架的三维磁粒子成像重建方法，所述方法包括：
5.步骤s100，通过校准程序，获得系统矩阵；
6.步骤s200，基于测量信号、待重建的磁纳米粒子浓度和所述系统矩阵构建mpi重建逆问题；
7.步骤s300，在贝叶斯框架下将所述mpi重建逆问题转化为求解有关待重建的磁纳米粒子浓度的最大后验概率估计问题；
8.步骤s400，求解所述有关待重建的磁纳米粒子浓度的最大后验概率估计问题，获得磁纳米粒子的浓度分布；
9.步骤s500，基于所述磁纳米粒子的浓度分布，完成三维磁粒子成像重建，获得重建后的三维图像。
10.在一些优选的实施方式中，所述步骤s100，具体为：通过校准程序，获得系统矩阵其中，表示复数域，m表示系统矩阵的行数，n表示系统矩阵的列数。
11.在一些优选的实施方式中，所述mpi重建逆问题，具体为：s*c＝u v
12.其中，表示待重建的磁纳米粒子浓度，表示系统测量到的电压信号，表示噪声信号。
13.在一些优选的实施方式中，所述步骤s300，具体为：
14.在贝叶斯框架下，将所述mpi重建逆问题转化为求解有关待重建的磁纳米粒子浓度c的最大后验概率估计问题：
15.其中，表示待重建的磁纳米粒子浓度c的最大后验概率估计，p(u|c)表示有关待重建的磁纳米粒子浓度c的条件似然分布函数，p(c)表示待重建的磁纳米粒子浓度c的先验概率分布函数。
16.参数化的似然函数为：
17.其中，p(u|c；λ)表示参数化的有关待重建的磁纳米粒子浓度c的条件参数化似然分布函数，λ表示噪声方差。
18.在一些优选的实施方式中，所述步骤s400，具体为：
19.步骤s410，假设所述待重建的磁纳米粒子浓度c满足稀疏分布，进一步假设c满足零均值的高斯分布：
20.其中，表示变量服从高斯分布，α＝(α1，α2，α3，...，αn)
t
表示未知的超参数，αj表示向量α中的第j个元素，cj表示向量c中的第j个元素；
21.所述噪声的方差λ服从gamma分布，则有关方差λ的概率密度函数p(λ|a，b)：
22.其中，a表示预设的尺度参数，b表示预设的形状参数，a和b作为已知参数需要事先输入到算法中，其最优值需要根据实验结果来调节；γ()表示gamma函数；
23.步骤s420，通过ii型最大似然方法估计所述未知的超参数α，使边界似然函数l取最大值：
24.其中，中间变量b＝i sλ-1st
，i表示单位矩阵，λ-1
＝diag(α1，α2，α3，...，αn)，λ-1
表示由未知超参数组成的对角矩阵，当所述边界似然函数l取最大值时，得到最优的超参数
25.步骤s430，基于所述最优的超参数计算c的最大后验概率估计：
26.其中，均值μ为：μ＝∑s
tu27.方差∑为：∑＝(s
t
s ξ)-1
28.表示由最优超参数组成的对角矩阵；
29.上述计算过程需要不断的迭代优化求解，迭代结束后的均值μ
*
即为待重建的磁纳米粒子浓度c。
30.在一些优选的实施方式中，所述测量信号为通过mpi成像设备测得的电压信号。
31.本发明的另一方面，提出了一种基于贝叶斯框架的三维磁粒子成像重建系统，包括系统矩阵获取模块、重建逆问题构建模块、后验概率估计问题转化模块、磁纳米粒子的浓度分布计算模块和磁粒子成像模块；
32.所述系统矩阵获取模块，配置为通过校准程序，获得系统矩阵；
33.所述重建逆问题构建模块，配置为基于测量信号、待重建的磁纳米粒子浓度和所述系统矩阵构建mpi重建逆问题；
34.所述后验概率估计问题转化模块，配置为在贝叶斯框架下将所述mpi重建逆问题转化为有关待重建的磁纳米粒子浓度c的最大后验概率估计问题；
35.所述磁纳米粒子的浓度分布计算模块，配置为求解所述有关待重建的磁纳米粒子浓度的最大后验概率估计问题，获得磁纳米粒子的浓度分布；
36.所述磁粒子成像模块，配置为基于所述磁纳米粒子的浓度分布，完成三维磁粒子成像重建，获得重建后的三维图像。
37.在一些优选的实施方式中，所述mpi重建逆问题，具体为：s*c＝u v
38.其中，表示待重建的磁纳米粒子浓度，表示系统测量到的电压信号，表示噪声信号，表示复数域，m表示系统矩阵的行数，n表示系统矩阵的列数。
39.本发明的第三方面，提出了一种电子设备，包括：至少一个处理器；以及与至少一个所述处理器通信连接的存储器；其中，所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现上述的基于贝叶斯框架的三维磁粒子成像重建方法。
40.本发明的第四方面，提出了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于被所述计算机执行以实现上述的基于贝叶斯框架的三维磁粒子成像重建方法。
41.本发明的有益效果：
42.(1)本发明通过引入稀疏先验和贝叶斯框架，将mpi三维图像重建问题转化成求解磁纳米粒子分布的最大似然估计问题，极大地加快了超大系统矩阵维度下的重建速度和重建实时性，提高重建精度，进一步提升mpi成像系统的实用性和临床转化潜力。
附图说明
43.通过阅读参照以下附图所作的对非限制性实施例所作的详细描述，本技术的其它
特征、目的和优点将会变得更明显：
44.图1是本技术实施例中基于贝叶斯框架的三维磁粒子成像重建方法的流程示意图；
45.图2是本技术实施例中基于贝叶斯框架的三维磁粒子成像重建系统的系统框图；
具体实施方式
46.下面结合附图和实施例对本技术作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。
47.需要说明的是，在不冲突的情况下，本技术中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本技术。
48.本发明提供一种基于贝叶斯框架的三维磁粒子成像重建方法，本方法基于稀疏先验和贝叶斯框架展开求解，可以大大加快超大系统矩阵维度下的重建速度，提高重建精度，进一步提升mpi成像系统的实用性和临床转化潜力。
49.本发明的一种基于贝叶斯框架的三维磁粒子成像重建方法，所述方法包括：
50.步骤s100，通过校准程序，获得系统矩阵；
51.步骤s200，基于测量信号、待重建的磁纳米粒子浓度和所述系统矩阵构建mpi重建逆问题；
52.步骤s300，在贝叶斯框架下将所述mpi重建逆问题转化为有关待重建的磁纳米粒子浓度的最大后验概率估计问题；
53.步骤s400，求解所述有关待重建的磁纳米粒子浓度的最大后验概率估计问题，获得磁纳米粒子的浓度分布；
54.步骤s500，基于所述磁纳米粒子的浓度分布，完成三维磁粒子成像重建，获得重建后的三维图像。
55.本方法基于事先通过校准程序测得的系统矩阵和mpi成像系统测得的磁纳米粒子产生的感应电压信号构建重建逆问题。进而将该逆问题转到贝叶斯框架下进行求解，在求解过程中，引入磁纳米粒子分布的稀疏先验和噪声分布的gamma先验，通过求解未知量(在本技术中为待重建的磁纳米粒子浓度c)的最大后验概率分布，最终得到磁粒子的准确的三维空间分布状态。
56.为了更清晰地对本发明基于贝叶斯框架的三维磁粒子成像重建方法进行说明，下面结合图1对本发明实施例中各步骤展开详述。
57.本发明第一实施例的基于贝叶斯框架的三维磁粒子成像重建方法，包括步骤s100-步骤s500，各步骤详细描述如下：
58.步骤s100，通过校准程序，获得系统矩阵；其中系统矩阵是根据不同设备，不同测量物质测定的。
59.在本实施例中，所述步骤s100，具体为：通过校准程序，获得系统矩阵其中，表示复数域，m表示系统矩阵的行数，n表示系统矩阵的列数。
60.对人体或者说成像整体而言，病灶的体积都是很小的，也就是病灶的分布具有稀疏性。因此，可以考虑将稀疏性引入到mpi重建领域，加速mpi重建过程。
61.步骤s200，基于测量信号、待重建的磁纳米粒子浓度和所述系统矩阵构建mpi重建逆问题；
62.在本实施例中，所述测量信号为通过mpi成像设备测得的电压信号，是系统中的磁纳米粒子与感应设备产生的感应电压信号。
63.在本实施例中，所述mpi重建逆问题，具体为：s*c＝u v
64.其中，表示待重建的磁纳米粒子浓度，表示系统测量到的电压信号，表示噪声信号。
65.步骤s300，在贝叶斯框架下将所述mpi重建逆问题转化为有关待重建的磁纳米粒子浓度的最大后验概率估计问题；
66.所述步骤s300，具体为：
67.在贝叶斯框架下，将所述mpi重建逆问题转化为求解有关待重建的磁纳米粒子浓度c的最大后验概率估计问题：
68.其中，表示待重建的磁纳米粒子浓度c的最大后验概率估计，p(u|c)表示有关待重建的磁纳米粒子浓度c的条件似然概率分布函数，p(c)表示待重建的磁纳米粒子浓度c的先验概率分布函数；
69.参数化的似然函数为：
70.其中，p(u|c；λ)表示待重建的有关待重建的磁纳米粒子浓度c的条件参数化的似然分布函数，λ表示噪声方差。
71.步骤s400，求解所述有关待重建的磁纳米粒子浓度的最大后验概率估计问题，获得磁纳米粒子的浓度分布；
72.在本实施例中，所述步骤s400，具体为：
73.步骤s410，假设所述待重建的磁纳米粒子浓度c满足稀疏分布，进一步假设c满足零均值的高斯分布：
74.其中，表示变量服从高斯分布，α＝(α1，α2，α3，...，αn)
t
表示未知的超参数，αj表示向量α中的第j个元素，cj表示向量c中的第j个元素；
75.所述噪声的方差λ服从gamma分布，则有关方差λ的概率密度函数p(λ|a，b)：
76.其中，a表示预设的尺度参数，b表示预设的形状参数，a和b作为已知参数需要事先输入到算法中，其最优值需要根据实验结果来调节；γ()表示gamma函数；
77.步骤s420，通过ii型最大似然方法估计所述未知的超参数α，使边界似然函数l取
最大值：
78.其中，中间变量b＝i sλ-1st
，i表示单位矩阵，λ-1
＝diag(α1，α2，α3，...，αn)，λ-1
表示由未知超参数组成的对角矩阵，当所述边界似然函数l取最大值时，得到最优的超参数
79.步骤s430，基于所述最优的超参数计算c的最大后验概率估计：
80.其中，均值μ为：μ＝∑s
tu81.方差∑为：∑＝(s
t
s ξ)-1
82.表示由最优超参数组成的对角矩阵；
83.上述计算过程需要不断的迭代优化求解，迭代结束后的均值μ
*
即为待重建的磁纳米粒子浓度c。
84.步骤s500，基于所述磁纳米粒子的浓度分布，完成三维磁粒子成像重建，获得重建后的三维图像。
85.本发明第二实施例的基于贝叶斯框架的三维磁粒子成像重建系统，如图2所示，包括系统矩阵获取模块、重建逆问题构建模块、后验概率估计问题转化模块、磁纳米粒子的浓度分布计算模块、磁粒子成像模块；
86.所述系统矩阵获取模块，配置为通过校准程序，获得系统矩阵；
87.所述重建逆问题构建模块，配置为基于测量信号、待重建的磁纳米粒子浓度和所述系统矩阵构建mpi重建逆问题；
88.所述后验概率估计问题转化模块，配置为在贝叶斯框架下将所述mpi重建逆问题转化为有关待重建的磁纳米粒子浓度的最大后验概率估计问题；
89.所述磁纳米粒子的浓度分布计算模块，配置为求解所述有关待重建的磁纳米粒子浓度的最大后验概率估计问题，获得磁纳米粒子的浓度分布；
90.所述磁粒子成像模块，配置为基于所述磁纳米粒子的浓度分布，完成三维磁粒子成像重建，获得重建后的三维图像。
91.在本实施例中，所述mpi重建逆问题，具体为：s*c＝u v
92.其中，表示待重建的磁纳米粒子浓度，表示系统测量到的电
压信号，表示噪声信号，表示复数域，m表示系统矩阵的行数，n表示系统矩阵的列数。
93.所属技术领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统的具体工作过程及有关说明，可以参考前述方法实施例中的对应过程，在此不再赘述。
94.需要说明的是，上述实施例提供的基于贝叶斯框架的三维磁粒子成像重建系统，仅以上述各功能模块的划分进行举例说明，在实际应用中，可以根据需要而将上述功能分配由不同的功能模块来完成，即将本发明实施例中的模块或者步骤再分解或者组合，例如，上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块，以完成以上描述的全部或者部分功能。对于本发明实施例中涉及的模块、步骤的名称，仅仅是为了区分各个模块或者步骤，不视为对本发明的不当限定。
95.本发明第三实施例的一种电子设备，包括：至少一个处理器；以及与至少一个所述处理器通信连接的存储器；其中，所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现上述的基于贝叶斯框架的三维磁粒子成像重建方法。
96.本发明第四实施例的一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于被所述计算机执行以实现上述的基于贝叶斯框架的三维磁粒子成像重建方法。
97.所属技术领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的存储装置、处理装置的具体工作过程及有关说明，可以参考前述方法实施例中的对应过程，在此不再赘述。
98.术语“第一”、“第二”等是用于区别类似的对象，而不是用于描述或表示特定的顺序或先后次序。
99.术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备/装置不仅包括那些要素，而且还包括没有明确列出的其它要素，或者还包括这些过程、方法、物品或者设备/装置所固有的要素。
100.至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征做出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。