图1 MACT-MI原理
Fig.1 MACT-MI schematic
摘要 感应式磁声磁粒子浓度成像（MACT-MI）是一种基于磁声耦合效应的磁纳米粒子（MNPs）浓度成像新方法。针对MACT-MI逆问题成像速度较慢的问题，该文引入势函数构建声压与MNPs浓度的关系，提出一种基于无转置拟最小残差（TFQMR）算法的洛伦兹力势声源图像重建方法。该方法降低了逆问题理论公式的求解复杂度，在保证图像高分辨率的前提下，进一步提高了成像速度。首先，建立了多种尺寸、形状，以及噪声情况下的磁纳米粒子模型。其次，将获取的数据用于浓度计算公式中进行图像重建。最后，对重建结果进行质量分析，分别对比不同模型在不同方法下的重建分辨率和重建速度。仿真结果表明：在相同浓度条件下，该方法在无噪声干扰时，相关系数平均高于0.947 6、相对误差平均低于0.399 3、结构相似性平均高于0.95、平均图像重建时间缩短至39.84 s。同时，该方法在不同噪声模型下具有较强的抗噪性，为MACT-MI的临床应用提供了理论支撑。
关键词：感应式磁声磁粒子浓度成像 洛伦兹力 势声源 TFQMR算法 逆问题成像
纳米医学[1-4]是现代医学中的重要组成部分，它在提升疾病诊断和治疗效果中起着至关重要的作用。通过靶向药物递送[5]和智能诊断技术[6]，纳米医学能够精确地将治疗剂送到病灶部位，有效地减少了对健康组织的伤害。其中磁粒子成像（Magnetic Particle Imaging, MPI）[7]在纳米医学中使用广泛，但是在复杂病变区域存在成像分辨率较低的问题，因此有望进一步提升成像分辨率。感应式磁声磁粒子浓度成像（Magneto-Acoustic Concentration Tomo- graphy with Magnetic Induction, MACT-MI）[8]作为MPI的一个分支，通过磁场激发使得磁性纳米粒子（Magnetic Nanoparticles, MNPs）高频振动，并利用其产生的超声波信号进行成像。该过程融合了电磁技术、超声技术的优势，有效地减少了线圈之间的电磁干扰，极大地改善了成像的清晰度和准确度。
MACT-MI作为一种新型的成像方法，关于逆问题重建的研究并不多，但是磁感应磁声成像（Magneto- Acoustic Tomography with Magnetic Induction, MAT- MI）、MPI和热声成像等逆问题的研究为本文提供了许多指导性思路。He Bin团队在2005年首次提出了关于感应式磁声成像技术，其核心是结合静态磁场和脉冲磁场诱导涡流产生洛伦兹力，从而激发超声波[9]；2007年，该团队开发了柱面扫描系统，在凝胶模型中成功地重建了多层结构，并首次在生物组织中区分肌肉和脂肪的电导率差异[10]；2008年，该团队提出了多激励策略的创新算法，设计了不同线圈配置的多次激励，结合改进算法成功地重建内部电导率并进行对比[11]；2010年，该团队深入解析声信号产生机制，解释了早期MAT-MI图像中边界显著的现象[12]；2014年，该团队开发矢量波束成形算法重建涡流源，在仿真和实验中达到1.5 mm分辨率，验证了生物组织应用价值[13]。2018年，M. Ziolkowski等使用有限元法和时逆算法重建了电导率分布和洛伦兹力散度，为MAT-MI提供了理论支持[14]。2020年，A. R. Zywica等使用反投影算法和时间反演法，通过分析系统矩阵中的奇异值，将电导率的分布与超声波信号耦合起来，成功地重建了复杂形状物体的电导率分布[15]。2020年，Wang Lili等采用稀疏矩阵处理方法，结合代数重建技术从声压数据中重建声源图像[16]。2024年，淡新贤等使用鲸鱼优化算法对电流密度函数进行优化，并设计了一种改进梯度线圈，有效地提高了磁声信号的信噪比[17]。2024年，孔晓涵等基于低频表面超材料的特性，提出了一种电感-电容谐振阵列，提高了超低场磁共振成像的信噪比[18]。2025年，孟凡钦等采用间接边界元法并从逆问题角度求出线圈流函数分布，再利用线性规划寻找最优线圈结构，从而实现了高信噪比的超低场磁共振图像[19]。2025年，汤云东等根据X-space的图像重建算法理论框架和投影重建成像算法中MPI检测装置磁场均匀性的重要性提出了一种基于正方形亥姆霍兹线圈的改进开放式检测装置，使得磁声信号的检测结果更接近于理想情 况[20]。2024年，A. R. Zywica等提出了激发脉冲的形状和持续时间会影响MAT-MI的重建图像质量，其中1 ms的高斯脉冲（delta形状）提供最佳重建效果[21]。2024年，I. Woods提出了一个扩展的MAT- MI算法，能够有效处理完全未知的各向异性导电张量，并提供了收敛性分析确保算法的稳定性[22]。2024年，杨丹等提出了一种基于线性零磁场的动脉血管扫描成像方法，设计了零磁场线双向扫描的组合线圈结构，理论上实现了对血管断层图像的实时监测[23]。
关于MACT-MI逆问题的研究，2020年，本课题组许正阳等提出基于矩量法的感应式磁声磁粒子浓度成像方法，该方法使用时间反演法进行声源重建，但稳定性差，重建结果存在边界奇异性[24]。2021年，胡宇等提出了基于截断奇异值分解（Truncated Singular Value Decomposition, TSVD）的感应式磁声磁粒子浓度成像方法，该方法计算复杂度高，难以处理大规模问题，且重建结果内部信息重构性差[25]。2022年，许洪等提出基于最小二乘QR分解法（Least Squares QR, LSQR）-梯形公式法的磁声磁粒子浓度成像方法，但存在算法迭代次数多、重建图像抗噪性能差等问题[26]。2024年，王宇飞等提出了基于矢量声源和稳定双共轭梯度法（Biconjugate Gradient Stabilized method, BICGSTAB）的磁声磁粒子浓度重建方法，该方法在处理病态问题时，求解不够稳定且收敛速度较慢，需要有效的预处理来加速收敛[27]。
综上所述，本文针对目前MACT-MI研究存在的问题，提出了一种基于无转置拟最小残差法（Transpose- Free Quasi-Minimal Residual, TFQMR）[28-30]的洛伦兹力势声源图像重建方法。通过引入势函数获得声压矩阵，根据势函数与MNPs浓度的关系，进一步得到浓度与声压的关系表达式。本方法避免了传统声压-浓度公式的求解复杂性，使得成像速度显著提升。最后利用TFQMR算法进行求解，该算法对于大型稀疏系统具有高效的求解能力，解决了以往图像重建时出现的边界奇异性、内部信息重构性差、抗噪性差等问题。因此，本方法在保证成像质量的同时，进一步提高了成像速度。
MACT-MI的原理如图1所示，由亥姆霍兹线圈提供静磁场B1使MNPs达到饱和状态，由麦克斯韦线圈提供均匀梯度磁场B2（B2在x、y方向的分量远小于z方向，因此可以忽略不计）。MNPs磁化后与B2相互作用受磁力振动产生超声信号，组织中不同浓度的MNPs会产生不同的磁力和超声信号，通过超声换能器检测含有不同浓度信息的超声信号用于逆问题重建。
图1 MACT-MI原理
Fig.1 MACT-MI schematic
超顺磁性MNPs的电磁特性可以用磁偶极子模型[31]描述，磁场中MNPs受到的磁力[32]为
（1）
式中，f(r,t)为t时刻MNPs在r处受到的磁力；N(r)为r处MNPs的浓度；m为MNPs的磁矩；ez为z方向的单位向量。
饱和磁化状态下MACT-MI的声场正问题满足有源线性声压波动方程，即
（2）
式中，p1(r,t)为MNPs在t时刻r处受振动产生的声压；cs为生物组织内声速；磁力散度
为声源项。
在饱和磁化状态下MACT-MI的研究中，以麦克斯韦线圈中通入电流的瞬间作为初始时刻，此时研究区域中的梯度磁场还未建立，超顺磁性MNPs所受磁力的磁力密度为零，粒子不发生振动，即不存在声源，因此可认为初始条件为
（3）
在MACT-MI系统中，超声换能器在t=|r-r0|/cs时刻r0处接收到的磁声信号w(r0,t)可表示为
（4）
式中，“*”表示卷积运算；p(r0,t)为检测点r0处t时刻声压信号；h(t)为换能器的时域脉冲响应函数；g(t)为噪声信号。
在频域上可表示为
（5）
式中，W(r0,w)、P(r0,w)、H(w)、G(w)分别为w(r0,t)、p(r0,w)、h(t)、g(t)的傅里叶函数，由维纳滤波反卷积理论可得到超声换能器检测点的原始声场为
（6）
式中，FFT-1( · )为反傅里叶符号；D=1/(y|H(w)|)，D为一个与|H(w)|变化趋势相反的函数，用以减小噪声的影响，y为乘数因子，用以匹配|H(w)|周期内的变化。反卷积构造的声场数据矩阵包括声场幅值矩阵和声场位置矩阵。
假定洛伦兹力近似无旋，即Ñ×f(r,t)=0，引入势函数V，满足f(r,t)=-ÑV，将其用于声压波动方程得到
（7）
式中，
为拉普拉斯算子；k为波数；P(r,k)为频域中位置r处的声压。
根据MAT-MI的理论推导[33]，波数k=w/cs，由于V的声压波动方程为0，即(Ñ2+k2)V=0，由式（7）化简得到声压表达式为
（8）
式中，Gk(r,r0)为波数k时声源点r处和超声换能器r0处响应的格林函数值；W为整个研究区域；P(r0,k)表示换能器波数k时的磁声信号，得到球坐标系下V的重建表达式为
（9）
式中，
为格林函数的共轭函数；S0为积分区域。
将格林函数进行化简[34]，相应的重建公式为
（10）
式中，|r-r0|为声源点与超声换能器之间的距离；z0和z1分别为超声换能器和声源点在z轴截距的大小。
根据式（1）中磁力与浓度的关系，引入势函数可以将声压p(r,t)和浓度N(r)相互连接。通过超声换能器检测到的声压数据，求解得到浓度数据，从而进行图像重建。
本文采用有限元网格划分的方法将磁场数据和声场数据离散化，得到声压幅值矩阵和系统势函数矩阵，构建浓度矩阵表达式，采用TFQMR方法求解浓度方程，最终获得浓度分布。逆问题的研究路线如图2所示，具体的重建步骤如下：
图2 逆问题研究路线
Fig.2 Inverse problem research roadmap
1）有限元法离散原始声场信号
根据有限元离散处理方法，对包围在MNPs检测闭合曲面上接收到的声场信号进行有限元网格剖分。将连续积分区域S0分割成有限个小的离散元素，每个元素代表声源的一个小区域，并且每个小区域对观测点的影响可以单独计算累加，通过组装这些局部解获得全局解，忽略其他方向的影响，只考虑z方向声源向量，可以得到势函数的离散化方程为
（11）
由式（11）计算势函数的负梯度得到洛伦兹力，联立式（1）与式（11），得到计算公式为
（12）
式中，
为t时刻声源点位置的梯度大小；
为离散化系统位置矩阵。
2）获取重建区域数据
重建区域是以生物组织内部MNPs集群所在位置为中心，选取边长为L的正方形区域，并对其进行有限元划分为
大小的网格，得到每一个网格处的数据信息。具体数据信息包括梯度磁场信息
、相关系数
及系统位置矩阵。
3）构建矩阵关系式
根据式（12）中声压与浓度之间的关系，通过系统矩阵Q进行表述，将这种关系概括为一个矩阵方程，表达式为
（13）
式中，系统矩阵Q包括梯度磁场数据、系数数据、系统位置矩阵A；M为声压梯度参数矩阵；浓度矩阵N为待求矩阵x。
4）利用TFQMR方法求解浓度分布
（14）
式中，QTQ为大小为n2×n2的方阵，记为P；N为n2×1的矩阵，记为x；QTM为n2×1的矩阵，记为R，式（14）可进一步改写为
（15）
算法步骤具体如下：
（1）条件初始化
①选择初始解x0，计算初始残差r0=R-Px0，设
用作双共轭的向量。
②初始化变量：q0=
,
、a0=1、J0=1。（其中a0=1表示初始迭代步长为1，J0=1表示初始松弛因子为1）。
③初始化搜索方向：使用初始残差向量作为第一个搜索方向q0=r0和用系统矩阵P作用于q0得到初始向量v0=Pq0。
（2）迭代过程
①残差方向更新：计算新的方向向量va=Pqa-1；更新变量a0=qa-1/
,
；更新解量
= xa-1+aaqa-1；更新残差向量
=ra-1+aava-1。
②残差修正：TFQMR算法的每步迭代会进行两次更新，以近似最小化残差，进行如下修正Ja=
/
；修正后的解为xa=+
；修正后的残差为ra=+
。
③方向向量的更新：计算新的方向标量ba=aa-1qa/qa-1Ja-1；更新搜索方向向量qa= ra+ba(qa-1-Ja-1va-1)。
④收敛判断：检查当前残差的大小，设置||ra||≤1×10-3时满足收敛条件，算法终止，输出解xa；否则，继续进行下一次迭代。
（3）终止条件
TFQMR的终止条件是满足以下两个之一：
①残差收敛：残差小于给定的容差e=1×10-3，则终止迭代。
②最大迭代次数：设置最大迭代次数100，若达到最大迭代次数仍未达到预期精度，终止算法，返回当前解。
（4）利用所得数据进行浓度重建
通过得到的浓度数据，可以精确地反映MNPs在目标区域内的分布情况，从而得到重建图像。
为验证基于TFQMR的洛伦兹力势声源MACT- MI图像重建方法的正确性，借助COMSOL建立二维轴对称模型进行计算仿真研究。分别建立不同半径、形状以及不同噪声条件下的MNPs模型，并与TV-MoM法、正则化预优化LSQR法、BICGSTAB法等进行对比，最后对仿真结果进行分析。
二维模型在xOy平面上分布示意图如图3所示，外围蓝色圆形区域是生物组织，设置组织中声速均匀且大小为1 500 m/s，以原点为中心选取25 mm×25 mm的矩形区域为成像区域，在xOy平面上沿半径40 mm的圆周等距分布165个超声换能器，MNPs均匀分布在原点为中心的红色圆形区域内。
图3 仿真模型
Fig.3 Simulation model
本研究MNPs选用聚乙二醇（PEG）涂层的超顺磁性氧化铁（SPIO）纳米粒子[35]，该粒子可在外加磁场作用下快速响应，且在磁场移除后不会保留磁性，PEG涂层的SPIO纳米粒子规格见表1。
该模型存在于麦克斯韦线圈和亥姆霍兹线圈产生的饱和磁场下。设置亥姆霍兹线圈和麦克斯韦线圈半径大小相同为100 mm，亥姆霍兹线圈通入电流方向相同，大小为50 A的电流，且线圈间距与线圈半径相同为100 mm。麦克斯韦线圈通入相反电流，线圈间距为线圈半径的
倍，即
mm。
表1 PEG涂层的SPIO纳米粒子规格
Tab.1 Specifications of PEG-coated SPIO nanoparticles
参 数数 值 粒径/nm8～10 饱和磁化强度/(kA/m)370 不同分子量PEG涂层厚度333Da～3.5 nm20 000Da～6.5 nm 水波动尺寸333Da～20 nm20 000Da～80 nm 循环半衰期/min105
2.2.1 不同尺寸对重建结果的影响
为分析不同尺寸的MNPs区域对重建结果的影响，以图3的仿真模型为基础，设置MNPs浓度为1×1019/mL，分别研究了四种不同重建算法下半径分别为2、5、8、12.5 mm的MNPs集群重建结果如图4所示。由图4可知，本文提出的TFQMR法相较于全变分正则化矩量法（Total Variation regularization Method of Moments, TV-MoM）法克服了边界奇异性问题；相较于BICGSTAB法克服了重建内部存在条纹、浓度不均匀的问题，保证了较高的重建图像分辨率。
为了更加直观地观察目标浓度与重建浓度之间的区别，在模型点 (0, 0, -25 mm) 和点 (0, 0, 25 mm)之间取三维截线，分别对比四种重建方法浓度截线差异，如图5所示，其中红色实线表示目标浓度，蓝色虚线表示重建浓度。由图5可知，TV-MoM法重建结果由于边界奇异性，其重建浓度曲线在边界处存在较大误差，且误差随半径增大而增大。BICGSTAB法内部浓度重建存在误差且边界处存在不稳定波动；TFQMR法能够均匀地重建内部浓度，在半径尺寸较大时，仍保持较高的稳定性，仅在边缘处存在较小重建误差，整体上更符合目标浓度曲线。
（a）TV-MoM法重建图（b）正则化预优化LSQR法重建图（c）BICGSTAB法重建图（d）TFQMR法重建图（e）目标浓度图
图4 不同算法下不同半径圆模型的MNPs集群重建结果对比
Fig.4 Comparison of MNPs cluster reconstruction results for different radius circular models using different algorithms
（a）TV-MoM法重建图（b）正则化预优化LSQR法重建图 （c）BICGSTAB法重建图 （d）TFQMR法重建图
图5 不同算法下不同半径圆模型的MNPs集群截线浓度对比
Fig.5 Comparison of the concentration of MNPs cluster cutting lines for different radius circular models using different algorithms
为进一步评价并比较TFQMR法与其他三种算法得到的图像质量，引入相关系数CC、相对误差ER、结构相似性SSIM对图像质量进行评估。
相关系数的表达式为
（16）
相对误差表达式为
（17）
式中，Nn为目标粒子n浓度值；Nn,r为重建区域r中粒子n浓度值；
为目标粒子浓度平均值；
为重建粒子浓度平均值。
给定两个图像X和Y，结构相似度的表达式为
（18）
式中，
和
分别为图像像素X和Y的均值；
、
分别为图像X和Y的标准差；
、
分别为图像X和Y的方差；
为X和Y的协方差；c1=(x1 j)2、c2=(x2 j)2为用于维持稳定的常数，j为像素值的动态范围，其中x1=0.01、x2=0.03。
不同半径模型下重建结果CC、ER、SSIM对比见表2。由表2可知，四种不同半径模型下，TFQMR重建结果中CC的平均值为0.960 0，ER的平均值为0.416 0，SSIM均稳定在0.918 8以上。结果表明，TFQMR法与其他算法相比，重建图与目标图之间相关系数和结构相似度更高，重建相对误差更小，故TFQMR算法能够有效地进行成像。
为了进一步比较成像速度，计算每种算法的迭代次数、成像时间与相对残差，对比见表3。其中TV-MoM法平均需66.25次迭代得到重建图像、平均成像时间为639.98 s；正则化预优化LSQR法平均迭代次数为47.75、平均成像时间为308.3 s；BICGSTAB法平均迭代次数为9.625次、平均成像时间为43 s；TFQMR法相比其他三种方法平均迭代次数更少为8.5次、成像时间更短平均为38.5 s，即可得到重建图像，且相对残差更小平均为0.008 9。
表2 不同算法下不同半径圆模型下重建结果CC、ER、SSIM对比
Tab.2 Comparison of CC, ER and SSIM reconstruction results under different radius circular models for various algorithms
半径/mmCCERSSIM TV-MoM法正则化预优化LSQR法BICGSTAB法TFQMR法TV-MoM法正则化预优化LSQR法BICGSTAB法TFQMR法TV-MoM法正则化预优化LSQR法BICGSTAB法TFQMR法 20.886 20.917 80.901 90.928 60.489 10.416 40.459 50.431 00.789 50.981 40.953 20.986 5 50.894 50.956 30.926 10.956 20.491 50.384 20.412 70.385 00.663 20.979 00.937 50.978 5 80.896 20.949 60.914 20.941 50.493 20.408 60.447 40.415 00.593 20.958 60.922 90.956 2 12.50.875 20.954 60.912 10.947 20.498 30.423 80.461 90.432 80.576 30.914 80.893 80.918 8
2.2.2 不同形状下的重建结果
为进一步验证算法的几何适应性，本文构建了圆、椭圆、三角形及六边形四种典型MNPs分布组合模型，如图6所示。
表3 不同算法下不同半径圆模型下迭代次数、重建时间和相对残差对比
Tab.3 Different algorithms compare the number of iterations, reconstruction time and relative residuals under different radius circular models
半径/mm迭代次数重建时间/s相对残差 TV-MoM法正则化预优化LSQR法BICGSTAB法TFQMR法TV-MoM法正则化预优化LSQR法BICGSTAB法TFQMR法TV-MoM法正则化预优化LSQR法BICGSTAB法TFQMR法 2664614.513559.6262.364.958.90.0560.008 50.010.008 8 567489.58668.2316.442.335.60.0150.008 60.009 50.008 4 865497.57663.2327.633.431.20.0340.00880.009 30.008 6 12.5674876668.9326.931.426.90.0250.009 40.009 70.009 7
图6 不同形状仿真模型
Fig.6 Different shape simulation models
四种模型中MNPs集群的尺寸及浓度参数见表4，其中r1为圆半径；l1为椭圆长轴，s1为椭圆短轴；l2为等腰三角形底边长，h1为等腰三角形的高；l3为六边形的上下边长，l4为六边形高；N表示不同形状MNPs浓度大小。
表4 四种不同形状模型中磁性纳米粒子集群尺寸及浓度
Tab.4 Size and concentration of magnetic nanoparticles in four different shape models
模型尺寸/mm浓度N/(1019/mL) ar1=5; l1=7; s1=41 br1=5; l2=10; h1=71 cl1=7; s1=4; l3=81 dl1=7; s1=4; l2=10; h1=71
不同形状模型重建结果如图7所示。由重建结果可知，重建浓度图像边界清晰，对复杂的形状仍具有较好的重建效果，整体与目标浓度图相吻合。
在点 (0, 0, -25 mm) 和点 (0, 0, 25 mm) 之间取三维截线，不同形状模型在截线上的浓度如图8所示。由图8可知，该方法整体上重建误差较小，与目标浓度截线基本吻合，证明了TFQMR可以重建不同形状的MNPs，并能取得较佳的成像效果。
不同形状模型重建质量分析见表5。由表5可知，四种模型的相关系数平均值为0.951 8，相对误差平均值为0.382 6，结构相似性平均值为0.940 1，迭代次数在8.5次即可收敛，平均重建时间为41.525 s。因此在多种复杂形状模型重建上，该方法在保证图像分辨率的同时，仍能保持一定的快速性。
图7 不同形状模型重建结果
Fig.7 Reconstruction results of models with different shapes
图8 不同形状模型截线上浓度重建结果
Fig.8 Reconstruction results of concentration on the cutting lines of different shape models
2.2.3 噪声情况下重建结果分析
在实际生物医学成像过程中，噪声是不可忽略的因素。MACT-MI成像过程中，信号受到电磁干扰、组织吸收、环境噪声等因素的影响，会导致重建图像质量下降。因此，评估成像算法在噪声环境下的鲁棒性至关重要。
表5 不同形状模型重建质量分析
Tab.5 Analysis of reconstruction quality for different shape models
模型CCERSSIM迭代次数重建时间/s a0.952 80.364 20.954 5943.5 b0.948 60.391 50.935 3839.2 c0.949 70.394 00.928 2839.3 d0.956 10.380 60.942 3944.1
为了模拟实验环境中的信号干扰，本文在仿真数据中加入高斯白噪声，并设定不同的信噪比（Signal-to-Noise Ratio, SNR）进行评估，SNR的定义为
（19）
式中，P1为原始磁声信号功率；P2为噪声信号功率。
本文对比了TFQMR法与LSQR-梯形公式法、正则化预优化LSQR法、BICGSTAB方法在5、10、15、20 dB四种不同SNR水平的重建结果，如图9所示。不同算法下不同噪声水平重建结果CC、ER、SSIM对比见表6。LSQR-梯形法和BICGSTAB法在高噪声环境下的图像重建质量明显下降，出现明显的伪影、边界模糊、内部浓度分布不均匀的问题；TFQMR方法在不同噪声水平下均能保持稳定的重建质量、边界清晰、内部浓度均匀，图像质量优于其他方法。
图9 不同噪声条件下的重建结果图
Fig.9 Reconstructed results under different noise conditions
表6 不同算法下不同噪声水平重建结果CC、ER、SSIM对比
Tab.6 Comparison of CC, ER and SSIM of reconstruction results under different algorithms and different noise levels
噪声/dBCCERSSIM LSQR-梯形公式法BICGSTAB法正则化预优化LSQR法TFQMR法LSQR-梯形公式法BICGSTAB法正则化预优化LSQR法TFQMR法LSQR-梯形公式法BICGSTAB法正则化预优化LSQR法TFQMR法 50.368 00.374 30.508 20.501 20.836 50.862 80.603 20.656 30.406 60.403 60.569 30.562 9 100.594 00.898 40.905 30.906 20.594 20.489 10.417 50.415 90.437 50.526 70.650 80.654 1 150.801 30.917 20.926 70.925 30.667 10.450 00.434 90.436 10.498 10.684 10.706 80.702 5 200.902 30.920 80.937 50.931 20.480 80.443 90.412 20.413 70.629 40.837 20.897 90.907 5
由表6可知，不同算法在不同噪声水平下的表现具有明显差异，TFQMR方法的CC、SSIM平均值分别为0.815 9和0.706 8，均明显高于LSQR-梯形公式法和BICGSTAB法，且ER平均值相比较小为0.480 5。证明了本方法具有稳定的图像重建质量和较高的图像分辨率。
不同算法下不同噪声水平迭代次数、重建时间和相对残差对比见表7。由表7可知，TFQMR法平均迭代次数为12.25，平均重建时间为53.5 s，重建速度相比较快，并且相对残差相对较低为0.009 2，证明了TFQMR法具有较快的成像速度。
表7 不同算法下不同噪声水平迭代次数、重建时间和相对残差对比
Tab.7 Comparison of iteration times, reconstruction time and relative residuals under different noise levels for different algorithms
噪声/dB迭代次数重建时间/s相对残差 LSQR-梯形公式法BICGSTAB法正则化预优化LSQR法TFQMR法LSQR-梯形公式法BICGSTAB法正则化预优化LSQR法TFQMR法LSQR-梯形公式法BICGSTAB法正则化预优化LSQR法TFQMR法 52928.55414126.5123.9336.561.00.420 00.390 00.008 00.007 8 103113.55313135.258.9426.856.70.260 00.009 80.009 10.009 4 153110.55311135.545.8432.448.00.150 00.009 70.009 60.009 8 2032105211139.645.9411.748.30.092 00.009 60.009 90.009 9
为提高MACT-MI重建图像的分辨率与重建速度，本文提出了基于TFQMR算法的洛伦兹力势声源图像重建方法，得到浓度数据并成像。通过建模与仿真，验证了理论公式的正确性及算法的高效性，本文得出以下结论：
1）提出了基于TFQMR算法的洛伦兹力势声源MACT-MI图像重建方法。根据势函数建立了洛伦兹力势声源重建理论，并利用TFQMR算法进行求解。该方法通过理论公式优化避免以往求解声压导函数的复杂性而带来的数据误差，同时，TFQMR算法减少了计算成本。本方法在保证图像分辨率的前提下，提高了成像速度。
2）研究对比了不同尺寸和形状MNPs模型的重建结果。本方法与TV-MoM法相比，克服了边界奇异性问题；与正则化预优化LSQR法相比，提高了重建速度；与BICGSTAB法相比，解决了内部重建不均匀的问题。证明本方法具有较强的普适性。
3）研究了不同噪声模型下的重建结果。本方法可以清晰地重建浓度的变化情况，并且在5～20 dB噪声水平下，仍保持稳定的成像质量，相对残差较低，证明本方法具有较强的稳定性与抗噪性。
本文采用的二维轴对称模型在重建中取得了较好的结果，但是在三维模型中的图像重建，还需要进一步探索。
参考文献
[1] 何祥, 荆慧. 基于纳米医学的光热治疗在肿瘤治疗中的研究进展[J]. 中国肿瘤临床, 2022, 49(8): 422-427.
He Xiang, Jing Hui. Research progress in photo- thermal therapy based on nanomedicine in tumor treatment[J]. Chinese Journal of Clinical Oncology, 2022, 49(8): 422-427.
[2] Pei Zerong, Chen Shuting, Ding Liqin, et al. Current perspectives and trend of nanomedicine in cancer: a review and bibliometric analysis[J]. Journal of Controlled Release, 2022, 352: 211-241.
[3] Huang Xiangang, Kong Na, Zhang Xingcai, et al. The landscape of mRNA nanomedicine[J]. Nature Medicine, 2022, 28(11): 2273-2287.
[4] Yang Zhe, Gao Di, Zhao Jing, et al. Thermal immuno- nanomedicine in cancer[J]. Nature Reviews Clinical Oncology, 2023, 20(2): 116-134.
[5] Liang Yujie, Duan Li, Lu Jianping, et al. Engineering exosomes for targeted drug delivery[J]. Theranostics, 2021, 11(7): 3183-3195.
[6] Chaouch M. Loop-mediated isothermal amplification (LAMP): an effective molecular point-of-care technique for the rapid diagnosis of coronavirus SARS-CoV-2[J]. Reviews in Medical Virology, 2021, 31(6): e2215.
[7] Gleich B, Weizenecker J. Tomographic imaging using the nonlinear response of magnetic particles[J]. Nature, 2005, 435(7046): 1214-1217.
[8] Shi Xiaoyu, Liu Guoqiang, Li Yanhong, et al. Simu- lation research on magneto-acoustic concentration tomography of magnetic nanoparticles with magnetic induction[J]. Computers in Biology and Medicine, 2020, 119: 103653.
[9] Xu Yuan, He Bin. Magnetoacoustic tomography with magnetic induction (MAT-MI)[J]. Physics in Medicine and Biology, 2005, 50(21): 5175-5187.
[10] Xia Rongmin, Li Xu, He Bin. Magnetoacoustic tomographic imaging of electrical impedance with magnetic induction[J]. Applied Physics Letters, 2007, 91(8): 83903.
[11] Ma Qingyu, He Bin. Magnetoacoustic tomography with magnetic induction: a rigorous theory[J]. IEEE Transactions on Bio-Medical Engineering, 2008, 55(2 Pt 2): 813-816.
[12] Li Xu, He Bin. Multi-excitation magnetoacoustic tomography with magnetic induction for bioimpe- dance imaging[J]. IEEE Transactions on Medical Imaging, 2010, 29(10): 1759-1767.
[13] Mariappan L, Hu Gang, He Bin. Magnetoacoustic tomography with magnetic induction for high- resolution bioimepedance imaging through vector source reconstruction under the static field of MRI magnet[J]. Medical Physics, 2014, 41(2): 022902.
[14] Ziolkowski M, Gratkowski S, Zywica A R. Analytical and numerical models of the magnetoacoustic tomography with magnetic induction[J]. COMPEL- the International Journal for Computation and Mathematics in Electrical and Electronic Engineering, 2018, 37(2): 538-548.
[15] Zywica A R, Ziolkowski M, Gratkowski S. Detailed analytical approach to solve the magnetoacoustic tomography with magnetic induction (MAT-MI) problem for three-layer objects[J]. Energies, 2020, 13(24): 6515.
[16] Wang Lili, Liu Guoqiang, Chen Siping, et al. Novel reconstruction algorithm of magnetoacoustic tomo- graphy based on ring transducer array for acoustic speed inhomogeneous tissues[J]. Medical Physics, 2020, 47(8): 3533-3544.
[17] 闫孝姮, 淡新贤, 陈伟华, 等. 基于改进目标场法梯度线圈的感应式磁声磁粒子浓度成像研究[J]. 电工技术学报, 2024, 39(14): 4305-4316.
Yan Xiaoheng, Dan Xinxian, Chen Weihua, et al. Magneto-acoustic magnetic particle concentration imaging based on improved target field method gradient coil[J]. Transactions of China Electro- technical Society, 2024, 39(14): 4305-4316.
[18] 孔晓涵, 张雅娜, 吴嘉敏, 等. 基于低频表面超材料的超低场磁共振成像信噪比增强方法[J]. 电工技术学报, 2024, 39(13): 3917-3927.
Kong Xiaohan, Zhang Yana, Wu Jiamin, et al. Signal-to-noise ratio enhancement method for ultra- low field magnetic resonance imaging based on low-frequency surface metamaterials[J]. Transactions of China Electrotechnical Society, 2024, 39(13): 3917-3927.
[19] 孟凡钦, 郭轶, 何为, 等. 基于边界元逆问题法的超低场磁共振大鼠脑部成像双通道射频线圈优化设计[J]. 电工技术学报, 2025, 40(1): 25-35, 63.
Meng Fanqin, Guo Yi, He Wei, et al. Optimal design of rat brain dual-channel RF coils for ultra-low-field magnetic resonance imaging based on boundary element inverse problem method[J]. Transactions of China Electrotechnical Society, 2025, 40(1): 25-35, 63.
[20] 汤云东, 丁宇彬, 金涛. 考虑磁场均匀性优化的开放式磁粒子成像检测装置改进方法[J]. 电工技术学报, 2025, 40(6): 1718-1728.
Tang Yundong, Ding Yubin, Jin Tao. Research on the improved method of open detection device for magnetic particle imaging considering magnetic field uniformity optimization[J]. Transactions of China Electrotechnical Society, 2025, 40(6): 1718-1728.
[21] Zywica A R, Ziolkowski M. The influence of shape and duration of excitation pulse on the reconstruction in magnetoacoustic tomography with magnetic indu- ction[J]. International Journal of Applied Electro- magnetics and Mechanics, 2024, 76(3): 356-369.
[22] Woods I. Reconstruction and convergence of anisotropic conductivity tensors using magneto-acoustic tomo- graphy with magnetic induction[D]. Philadelphia: Drexel University, 2024.
[23] 杨丹, 王雨忱, 李天兆, 等. 一种基于线性零磁场的动脉血管扫描成像方法仿真[J]. 电工技术学报, 2024, 39(2): 343-355.
Yang Dan, Wang Yuchen, Li Tianzhao, et al. Simulation of an arterial scanning imaging method based on linear zero magnetic field[J]. Transactions of China Electrotechnical Society, 2024, 39(2): 343- 355.
[24] Yan Xiaoheng, Xu Zhengyang, Chen Weihua, et al. Implementation method for magneto-acoustic con- centration tomography with magnetic induction (MACT-MI) based on the method of moments[J]. Computers in Biology and Medicine, 2021, 128: 104105.
[25] Yan Xiaoheng, Hu Yu, Guang Sichen, et al. Simu- lation research on magneto-acoustic concentration tomography of magnetic nanoparticles based on truncated singular value decomposition (TSVD)[J]. Medical & Biological Engineering & Computing, 2021, 59(11): 2383-2396.
[26] Yan Xiaoheng, Xu Hong, Li Jun, et al. Inverse problem of magneto-acoustic concentration tomo- graphy for magnetic nanoparticles with magnetic induction in a saturation magnetization state based on the least squares QR factorization method–trapezoidal method[J]. Medical & Biological Engineering & Computing, 2022, 60(11): 3295-3309.
[27] Yan Xiaoheng, Wang Yufei, Chen Weihua, et al. Simulation study on magneto-acoustic concentration tomography of magnetic nanoparticles based on vectorial acoustic source and BICGSTAB method[J]. Physica Scripta, 2024, 99(12).
[28] Freund R W. A transpose-free quasi-minimal residual algorithm for non-Hermitian linear systems[J]. SIAM Journal on Scientific Computing, 1993, 14(2): 470- 482.
[29] Khalevitsky Y V, Konovalov A V, Burmasheva N V, et al. Linear solver performance in elastoplastic problem solution on GPU cluster[C]//AIP Conference Proceedings. AIP Publishing, 2017, 1915(1).
[30] Izadkhah M M. Transpose-free quasi-minimal residual method based on tensor format for genera- lized coupled Sylvester tensor equations[J]. Calcolo, 2024, 61(3): 41.
[31] Zhao Zhiyuan, TorresDiaz I, Velez C, et al. Brownian dynamics simulations of magnetic nanoparticles captured in strong magnetic field gradients[J]. Journal of Physical Chemistry C, 2017, 121(1): 801-810.
[32] Yan Xiaoheng, Gao Peng, Cai Mingchen, et al. Simulation research on forward problem of mag- netoacoustic concentration tomograghy of magnetic nanoparticles with magnetic induction based on multi-coils[J]. Progress in Electromagnetics Research M, 2021, 104: 223-233.
[33] 刘国强. 磁声成像技术. 上册, 超声检测式磁声成像[M]. 北京: 科学出版社, 2018.
[34] Xia Rongmin, Li Xu, He Bin. Comparison study of three different image reconstruction algorithms for MAT-MI[J]. IEEE Transactions on Bio-Medical Engineering, 2010, 57(3): 708-713.
[35] Arbab A S, Wilson L B, Ashari P, et al. A model of lysosomal metabolism of dextran coated superpara- magnetic iron oxide (SPIO) nanoparticles: impli- cations for cellular magnetic resonance imaging[J]. NMR in Biomedicine, 2005, 18(6): 383-389.
Abstract Magnetoacoustic concentration tomography with magnetic induction (MACT-MI) is a new method for imaging the concentration of magnetic nanoparticles (MNPs) based on magnetoacoustic coupling. Existing studies on the inverse problem of MACT-MI have insufficient image resolution and slow reconstruction speed for MNPs. Therefore, this study focuses on the theoretical simplification and algorithmic optimization to simultaneously improve the quality and efficiency of reconstructed images.
Firstly, this paper presents a Lorentz force potential acoustic source image reconstruction method based on the transpose-free quasi-minimal residual (TFQMR) algorithm. A possible function is introduced to construct the correlation between the acoustic pressure signal and MNP concentration, effectively circumventing the complexity of the traditional acoustic pressure-concentration solution. The acoustic field data are discretized using the finite element method, transforming the concentration reconstruction into a sparse system of equations. Subsequently, the TFQMR algorithm solves the concentration distribution due to its fast convergence, low memory usage, and the ability to avoid explicitly computing the matrix transpose when dealing with large sparse systems, thereby significantly improving solution efficiency.
In addition, a variety of 2D axisymmetric simulation models were constructed in COMSOL. These models cover the distribution of MNPs at different sizes (radii 2 mm to 12.5 mm), geometries (circular, elliptical, triangular, hexagonal), and noise levels (signal-to-noise ratios 5 dB to 20 dB). The TV-MoM method, the regularized pre-optimization LSQR method, and the BICGSTAB method are compared. In the absence of noise interference, the average correlation coefficient is higher than 0.947 6, the average relative error is lower than 0.399 3, the average structural similarity is higher than 0.95, and the average image reconstruction time is reduced to 39.84 s. The shape contours of MNPs can be clearly reconstructed under different noise models, with relative residuals below 0.009 9, thereby maintaining stable imaging quality.
The results show that the present method overcomes the boundary-singularity problem in conventional image reconstruction, reduces streak artifacts within the reconstruction, and improves uneven concentration distribution and imaging resolution. The solution complexity is reduced by the Lorentz force potential acoustic source theory, and the matrix equations can be solved efficiently using the TFQMR algorithm, which has strong stability and noise immunity. The computational cost is reduced, and the reconstruction speed is enhanced.
Keywords：Magnetoacoustic concentration tomography with magnetic induction (MACT-MI), Lorentz force, potential acoustic source, transpose-free quasi-minimal residual (TFQMR) algorithm, inverse problem
中图分类号：TM12
DOI: 10.19595/j.cnki.1000-6753.tces.250397
国家自然科学基金青年科学基金项目（52207008）、国家自然科学基金面上项目（52477007）、辽宁省教育厅科技创新团队项目（LJ222410147025）和中国-波兰测控技术“一带一路”联合实验室开放课题项目（MCT202305）资助。
收稿日期 2025-03-11
改稿日期 2025-06-19
闫孝姮 女，1984年生，博士，教授，博士生导师，研究方向为电磁探测与成像。
E-mail: xiaohengyan@163.com（通信作者）
付 鹏 男，2000年生，硕士研究生，研究方向为电磁探测与成像。
E-mail: fupeng021488@163.com
（编辑 郭丽军）