一种饱和磁化状态下磁声磁粒子浓度图像重建方法

1.本发明属于浓度分布图像重建技术领域，具体涉及一种基于最小二乘qr分解法-梯形公式法的饱和磁化状态下磁声磁粒子浓度图像重建方法。


背景技术：

2.磁性纳米粒子(mnps)因其低毒性、良好的生物相容性、磁响应性以及在外加磁场作用下的可控性在生物医学领域得到广泛应用，包括：磁热疗、药物递送、靶向治疗、基因治疗等。磁粒子成像(magnetic particle imaging,mpi)是最早将磁纳米粒子应用于医学诊断的成像方法，2005年，gleich b等人首次在nature上对mpi成像方法进行报道。但其空间分辨率受理论和设备因素的影响，目前在1-5mm，为进一步提高空间分辨率，2020年，史晓玉等人首次提出了感应式磁声磁粒子浓度成像(magneto-acoustic concentration tomography with magnetic induction mact-mi)，该方法天然克服了驱动线圈和检测线圈之间的电磁干扰问题，融合了电磁技术、超声技术的优势，兼具无创、对比度好、灵敏度高以及空间分辨率高等优点。
3.对于mact-mi逆问题的研究，2019年12月25日，许正阳等人申请的《一种磁声耦合的磁性纳米粒子浓度图像重建方法》(专利号：201911020966.0)公开了一种磁声耦合的磁性纳米粒子浓度图像重建方法，使用时间反演法进行声源重建，涉及声压求导过程，放大了数据噪声对重建结果的影响，算法稳定性较差，且重建结果存在边界奇异性。2021年，胡宇等人提出了基于tsvd(截断奇异值分解)的磁声磁粒子浓度成像方法，该方法在对于系数矩阵为小型稀疏矩阵时成像效果较好，但应对大型稀疏系数矩阵成像效果不理想，且其理论公式仅适用于求解浓度均匀分布的情况，适用范围较窄，所以mact-mi的逆问题重建还需要进行进一步研究。


技术实现要素：

4.本发明所要解决的技术问题在于针对上述现有技术的不足，提供了一种饱和磁化状态下磁声磁粒子浓度图像重建方法。该制备方法利用lsqr-梯形公式法，能够快速、稳定、高质量地重建出浓度分布，图像边界清晰，具有良好的应用前景。
5.为解决上述技术问题，本发明采用的技术方案是：一种饱和磁化状态下磁声磁粒子浓度图像重建方法，其特征在于，该方法包括以下步骤：
6.步骤一、设定饱和磁化状态下磁声磁粒子的仿真初始条件，基于预设的仿真模型，获取超声换能器处的声压数据和重建区域处的梯度磁场数据；所述重建区域是以磁纳米粒子群为中心选取wmm
×
wmm的区域作为重建区域，对重建区进行有限元划分，划分成m
×
m个网格；重建区域处的梯度磁场数据包括每个网格处的梯度磁场数据；
7.步骤二、利用每个超声换能器处获取的声压数据构建声压矩阵，利用所述重建区域处的梯度磁场数据构建系统矩阵，该系统矩阵用于描述声压与磁性纳米粒子的浓度偏导数之间存在直接的对应关系；
8.步骤三、利用步骤二中构建的声压矩阵和系统矩阵用lsqr方法获取磁性纳米粒子的浓度偏导数分布；
9.步骤四、利用磁性纳米粒子浓度偏导数分布获取磁性纳米粒子浓度分布，从而获得饱和磁化状态下磁声磁粒子浓度重建图像。
10.优选地，步骤一中设定饱和磁化状态下磁声磁粒子的仿真初始条件，基于预设的仿真模型，获取超声换能器处的声压数据和重建区域处的梯度磁场数据；以磁纳米粒子群为中心选取wmm
×
wmm的区域作为重建区域，对重建区域进行有限元划分，划分成m
×
m个网格；重建区域处的梯度磁场数据包括每个网格处的梯度磁场数据；包括：
11.初始化磁纳米粒子参数、helmholtz线圈和maxwell线圈电流的仿真条件，选用型号为emg308作为磁性纳米粒子，磁纳米粒子群置于线圈中心位置处，分别向helmholtz线圈和maxwell线圈中通入电流，由helmholtz线圈提供静磁场b
sat
使磁性纳米粒子达到饱和状态，由maxwell线圈提供均匀梯度磁场bg；以磁纳米粒子群为中心，以固定的扫描半径画圆圈，在该圆圈上离散设置多个超声换能器，获取每个超声换能器接收的多个时间点的声压数据，并获取每个超声换能器处的原始声场p(r,t)。
12.优选地，以磁纳米粒子群为中心，以固定的扫描半径画圆圈，在该圆圈上离散设置165个超声换能器，获取每个超声换能器接收的470个时间点的声压数据；以磁纳米粒子群为中心选取50mm
×
50mm的区域作为重建区域，对重建区域进行有限元划分，划分成250
×
250个网格。
13.优选地，利用每个超声换能器处获取的声压数据构建声压矩阵，利用所述重建区域处的梯度磁场数据构建系统矩阵，该系统矩阵用于描述声压与磁性纳米粒子的浓度偏导数之间存在直接的对应关系，包括：
14.磁性纳米粒子受到的磁力表示为
[0015][0016]
式(1)中，f(r
′
,t)为t时刻磁性纳米粒子受到的磁力；r
′
为所求声源点的位置，t为时间；n(r
′
)为磁性纳米粒子的浓度；m表示磁性纳米粒子的磁矩；为t时刻声源点位置的梯度磁场大小；ez表示z方向上的单位向量；
[0017]
声源与声压关系式如下：
[0018][0019]
式(2)中，p(r,t)为t时刻任意点的声压，单位为帕；r为任意点位置；cs为生物组织内的声速，单位是m/s；为声源项；t为时间；r
′
为所求声源点的位置，r＝|r-r
′
|；ω为整个研究区域；
[0020]
由于磁力的单向性，声源项可表示为
[0021]
[0022]
将式(3)带入式(2)，得到
[0023][0024]
由于在计算瞬时时刻空间分布问题时，可按照常量处理，故可得
[0025][0026]
由上式可得，声压与浓度偏导数之间存在直接的对应关系，将这种对应关系通过系统矩阵a描述，可以抽象得到矩阵关系式：
[0027]
ax＝b
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0028]
其中a为系统矩阵，x为浓度偏导数，b为声压矩阵，声压矩阵b中的声压数据对应p(r,t)；根据式(5)和式(6)并结合梯度磁场数据，从而得到系统矩阵a。
[0029]
优选地，步骤三中利用步骤二中构建的声压矩阵和系统矩阵用lsqr方法获取磁性纳米粒子的浓度偏导数分布，包括：
[0030]
引入lsqr方法，对作为大型稀疏系数矩阵的系统矩阵进行求解，最终求得浓度偏导数分布；
[0031]
求解过程包括：
[0032]
令标准正交矩阵uk＝[u1,u2,
…
,uk](uj∈r
mi
)和双对角矩阵为
[0033][0034]
其中(α1,α2,
…
,αk∈r；β2,β3,
…
,β
k 1
∈r)。
[0035]
迭代过程如下：迭代过程如下：
[0036]
步骤301，条件初始化
[0037]
β1＝||b||；α1＝||a
t
u1||2；β1u1＝b；
[0038][0039]
步骤302，对角化系数矩阵
[0040]
β
j 1
＝||av
j-αjuj||2；α
j 1
＝||a
tuj 1-β
j 1
vj||2；
[0041][0042]
(j＝1,2,...,k)。
[0043]
步骤303，计算qr分解中间变量
[0044][0045][0046]
步骤304，更新x和中间变量w的值
[0047][0048]
步骤305，判断迭代条件
[0049]
如果满足||ax
k-b||≤ε则终止迭代，其中ε是允许的误差，设定ε＝0.01；
[0050]
对于j＝1,2,...,k重复步骤302～305。
[0051]
优选地，步骤四中利用磁性纳米粒子浓度偏导数分布获取磁性纳米粒子浓度分布，从而获取磁声磁粒子浓度重建图像，包括：
[0052]
已知浓度偏导数与浓度分布关系式满足下列关系：
[0053][0054]
假设网格的边长为l，l＝0.2mm，则浓度偏导数x可表示为:
[0055][0056]
浓度分布n(r')可表示为
[0057][0058]
引入计算简单、精度较高的梯形公式法求解，可得
[0059][0060]
其中0≤c≤249,1≤d≤250，结合x(c,d)边界条件则可得到浓度分布。
[0061]
本发明与现有技术相比具有以下优点：
[0062]
1、本发明是为了解决现有的磁声磁粒子浓度成像方法中存在的重建速度慢、重建质量差以及存在边界奇异性等问题。本发明基于lsqr-梯形公式法求解磁声磁粒子浓度分布，可以高质量、快速、稳定地重建出浓度分布，图像边界清晰，不存在奇异性。
[0063]
2、本发明将计算区域按照有限元法进行规定的栅格的网格划分，将声场数据和磁场数据离散化，构造反映声压和浓度偏导数联系的系统矩阵，利用lsqr方法重建浓度偏导数分布，利用梯形公式法重建得到浓度分布图像。
[0064]
下面通过附图和实施例对本发明的技术方案作进一步的详细说明。
附图说明
[0065]
图1为本发明的一种基于最小二乘qr分解法-梯形公式法的饱和磁化状态下磁声磁粒子浓度图像重建方法流程图。
[0066]
图2为仿真模型。
[0067]
图3为helmholtz线圈电流。
[0068]
图4为声压数据及磁场数据获取原理图。
[0069]
图5为渐变浓度模型。
[0070]
图6a为仿真结果图。
[0071]
图6b为目标浓度图。
具体实施方式
[0072]
实施例1
[0073]
磁声磁粒子浓度成像方法
[0074]
如图1所示，将计算区域按照有限元法进行规定的栅格的网格划分，将声场数据和磁场数据离散化，构造反映声压和浓度偏导数联系的系统矩阵，利用lsqr方法重建浓度偏导数分布，利用梯形公式法重建得到浓度分布图像。如图1所示，本发明通过以下技术方案来实现：本发明提供一种基于lsqr-梯形公式法的饱和磁化状态下磁声磁粒子浓度图像重建方法，具体包括以下步骤：
[0075]
s1、设定饱和磁化状态下磁声磁粒子的仿真初始条件，基于预设的仿真模型，获取超声换能器处的声压数据和重建区域处的梯度磁场数据；
[0076]
进一步地，所述饱和磁化状态下磁声磁粒子的仿真模型搭建方式为：借助comsol multiphysics5.6建立二维轴对称模型进行仿真研究，如图2所示，线圈材料为铜，以z轴为中心轴放置，线圈半径取r＝150mm，以原点为线圈中心，中间圆柱区域为研究区域。
[0077]
1.仿真条件
[0078]
(1)磁性纳米粒子的参数
[0079]
取自emg 308(ferrotec(usa)corporation)，其规格如表1所示。
[0080]
表1 emg 308规格
[0081][0082]
(2)激励电流
[0083]
在helmholtz线圈的上下线圈同时通入逆时针方向大小为50a的恒定电流，产生的磁场强度约为1.93
×
105a/m,大于emg308的饱和磁场强度，磁纳米粒子达到饱和状态。在maxwell线圈的上线圈入逆时针方向的电流，下线圈中通入同样大小的顺时针方向的电流，线圈中均通入时间特性如图3所示的电流。
[0084]
具体地，初始化磁纳米粒子参数、helmholtz线圈和maxwell线圈电流的仿真条件，选用型号为emg308作为磁性纳米粒子，磁纳米粒子群置于线圈中心位置处，分别向helmholtz线圈和maxwell线圈中通入电流，由helmholtz线圈提供静磁场b
sat
使磁性纳米粒子达到饱和状态，由maxwell线圈提供均匀梯度磁场bg；如图4所示，以磁纳米粒子群为中心，以固定的扫描半径画圆圈，在该圆圈上离散设置165个超声换能器，获取每个超声换能器接收的470个时间点的声压数据，并获取每个超声换能器处的原始声场p(r,t)；以磁纳米粒子群为中心选取50mm
×
50mm的区域作为重建区域，对重建区域进行有限元划分，划分成250
×
250个网格，提取每个网格处的梯度磁场数据。
[0085]
s2、利用每个超声换能器处获取的声压数据构建声压矩阵，利用所述重建区域处的梯度磁场数据构建系统矩阵，该系统矩阵用于描述声压与磁性纳米粒子的浓度偏导数之间存在直接的对应关系；
[0086]
利用每个超声换能器处获取的声压数据构建声压矩阵，利用所述重建区域处的梯度磁场数据构建系统矩阵，该系统矩阵用于描述声压与磁性纳米粒子的浓度偏导数之间存在直接的对应关系，包括：
[0087]
磁性纳米粒子受到的磁力表示为
[0088][0089]
式(1)中，f(r
′
,t)为t时刻磁性纳米粒子受到的磁力；r
′
为所求声源点的位置，t为时间；n(r
′
)为磁性纳米粒子的浓度；m表示磁性纳米粒子的磁矩；为t时刻声源点位置的梯度磁场大小；ez表示z方向上的单位向量；
[0090]
声源与声压关系式如下：
[0091][0092]
式(2)中，p(r,t)为t时刻任意点的声压，单位为帕；r为任意点位置；cs为生物组织内的声速，单位是m/s；为声源项；t为时间；r
′
为所求声源点的位置，r＝|r-r
′
|；ω为整个研究区域；
[0093]
由于磁力的单向性，声源项可表示为
[0094][0095]
将式(3)带入式(2)，得到
[0096][0097]
由于在计算瞬时时刻空间分布问题时，可按照常量处理，故可得
[0098][0099]
由上式可得，声压与浓度偏导数之间存在直接的对应关系，将这种对应关系通过系统矩阵a描述，可以抽象得到矩阵关系式：
[0100]
ax＝b
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0101]
其中a为系统矩阵，x为浓度偏导数，b为声压矩阵，声压矩阵b中的声压数据来自p(r,t)。
[0102]
根据式(5)和式(6)并结合梯度磁场数据，从而得到系统矩阵a，该系统矩阵的大小为77550
×
62500；利用超声换能器处的声压数据构造声压矩阵b，该矩阵大小为77550
×
1。
[0103]
s3、利用步骤二中构建的声压矩阵和系统矩阵用lsqr方法获取磁性纳米粒子的浓度偏导数分布；
[0104]
进一步地，获取磁性纳米粒子的浓度偏导数分布包括：
[0105]
引入lsqr方法，对作为大型稀疏系数矩阵的系统矩阵进行求解，最终求得浓度偏导数分布；
[0106]
求解过程包括：
[0107]
令标准正交矩阵uk＝[u1,u2,
…
,uk](uj∈r
mi
)和双对角矩阵为
[0108][0109]
其中(α1,α2,
…
,αk∈r；β2,β3,
…
,β
k 1
∈r)。
[0110]
迭代过程如下：迭代过程如下：
[0111]
步骤301，条件初始化
[0112]
β1＝||b||；α1＝||a
t
u1||2；β1u1＝b；
[0113][0114]
步骤302，对角化系数矩阵
[0115]
β
j 1
＝||av
j-αjuj||2；α
j 1
＝||a
tuj 1-β
j 1
vj||2；
[0116][0117]
(j＝1,2,...,k)。
[0118]
步骤303，计算qr分解中间变量
[0119][0120][0121]
步骤304，更新x和中间变量w的值
[0122][0123]
步骤305，判断迭代条件
[0124]
如果满足||ax
k-b||≤ε则终止迭代，其中ε是允许的误差，设定ε＝0.01；
[0125]
对于j＝1,2,...,k重复步骤302～305。
[0126]
s4、利用磁性纳米粒子浓度偏导数分布获取磁性纳米粒子浓度分布，从而获得饱和磁化状态下磁声磁粒子浓度重建图像。
[0127]
进一步地，利用磁性纳米粒子浓度偏导数分布获取磁性纳米粒子浓度分布，从而获得饱和磁化状态下磁声磁粒子浓度重建图像，包括：
[0128]
已知浓度偏导数与浓度分布关系式满足下列关系：
[0129][0130]
假设网格的边长为l，l＝0.2mm，则浓度偏导数x可表示为:
[0131][0132]
浓度分布n(r')可表示为
[0133][0134]
引入计算简单、精度较高的梯形公式法求解，可得
[0135][0136]
其中0≤c≤249,1≤d≤250，结合x(c,d)边界条件则可得到浓度分布。
[0137]
通过以上步骤便可重建出浓度分布，考虑到mnps在生物组织中是弥散渐变的，建立浓度渐变模型，如图5所示，外围圆形区域为仿生物区域，内部区域为仿mnps区域，mnps选用emg 308，重建结果及浓度重建目标如图6a及图6b中显示，本发明所提出的方法能够高质量地重建出浓度分布，不存在奇异性，为进一步评价本发明适用性，引入相关系数cc及相对
误差re，对目标浓度分布与重建浓度分布进行比较，具体采用方法采用如下公式：
[0138][0139][0140]
式中nn为目标纳米粒子浓度值，n
n,r
为重建纳米粒子浓度值，为目标纳米粒子浓度平均值，为目标纳米粒子浓度平均值。
[0141]
通过对本发明的评价，可以发现本发明的相关系数和相对误差分别为0.9816、0.2761，重建精度较高。
[0142]
以上所述，仅是本发明的较佳实施例，并非对本发明作任何限制。凡是根据发明技术实质对以上实施例所作的任何简单修改、变更以及等效变化，均仍属于本发明技术方案的保护范围内。