基于小波稀疏的磁性纳米粒子成像算法研究

张玉录 柯丽 杜强 赵宇楠 祖婉妮

0 引言

磁性纳米粒子成像(magnetic particle imaging，MPI)是一种可实现快速、精准成像的新兴临床成像技术。该技术通过检测磁性纳米粒子在变化磁场中的非线性磁化特性来重构粒子在待测区域中的浓度分布状况。在该过程中，变化的磁场产生移动的零场区域，只有零场区域中粒子磁矩会发生变化，进而在接收装置中产生MPI信号[1]，这不仅有利于MPI的高分辨率成像，也可实现MPI的快速动态成像，所以MPI在临床上的快速诊断和实时监控中具有广泛的应用前景，已被证实可以应用于细胞成像[2]、血管成像[3]、肿瘤成像[4]、干细胞追踪[5]、热疗[6]和药物靶向[7]等多个领域。

MPI成像方法是实现快速、精准成像的关键。目前MPI常用的成像方法主要有系统矩阵(system matrix，SM)方法和X-space方法[1]。由于X-space方法对零场区域的磁场和粒子在磁场中的弛豫时间提出更高要求[8]，SM方法成为MPI成像方案的首要选择。在SM方法中，由SM和测量向量构成矩阵方程的求解算法决定了成像速度和质量。

在SM成像方案的早期研究过程中，由于噪声的存在，矩阵方程的求解问题被转化为最小二乘问题，采用正则化和矩阵求逆技术求解该最小二乘问题是求解粒子浓度信息的常用方法[1]，但是大型MPI矩阵的逆计算复杂。随着MPI技术的发展，Lübeck大学的Panagiotopoulos等[1]针对大型SM在矩阵求逆中面临的大计算量问题，提出一种采用收敛迅速的代数重建算法(algebraic reconstruction technique，ART)求解吉洪诺夫正则化最小二乘模型成像的方法，成像结果表明ART算法大大缩短了重建时间,但是所成图像边缘较为模糊。Heidelberg大学的Storath等[9]针对吉洪诺夫正则化最小二乘模型所成图像边缘模糊的问题，提出用交替方向乘子算法(alternating direction method of multipliers，ADMM)求解非负融合Lasso模型，成像结果表明该算法可以实现粒子分布的精准成像，但是成像耗时是ART算法的7倍。近年来，研究人员相继提出随机奇异值分解加速的ART算法[10]、逐级背景分割的图像重建算法[11]、编码校准场景框架重构SM方法[12]等多种成像方法，虽然在求解速度或者成像质量上有一定的优势，但难以满足MPI的快速精准成像要求。

在ART快速成像算法基础上通过简单的图像处理改善图像边缘模糊问题，是实现MPI快速精准成像的一种有效方法。小波变换具有良好的表征图像局部特征的能力，且运算耗时较少。就边缘模糊的图像而言，采用简单数学工具对小波变换提取的图像边缘特征进行处理，可以快速去除图像中的干扰信息，提升边缘区域的图像质量。

本文为了实现MPI中粒子浓度空间分布的快速精准成像，针对系统矩阵成像方法构建的矩阵方程求解问题，提出一种基于小波变换和阈值算子的迭代MPI算法。在ART算法每次迭代结束后，首先采用小波变换提取图像中粒子分布边界的非平稳特征，然后通过阈值算子稀疏运算去除图像中的干扰信号，增加少量计算时间的同时通过强化图像非平稳特征提高了粒子分布边缘图像的清晰度，结合采用的线性零磁场MPI电磁系统可在短时间内实现对粒子分布状况的高质量成像。

1 磁性纳米粒子成像基础研究

1.1 磁性纳米粒子成像原理

MPI通过检测磁性纳米粒子在变化磁场中的非线性磁化特性[图1(a)]来重构粒子在待测区域中的浓度分布状况。在MPI中，选择场产生零场区域，聚焦场控制零场区域的自由移动，当对成像区域施加一个时变驱动磁场[图1(c)]时，只有零场区域中的粒子磁化强度改变，从而产生磁化响应信号[图1(b)]。粒子磁化响应信号会在接收线圈中感应产生电压信号[图1(d)]，该电压信号即为MPI信号。

图1 磁性纳米粒子对外部磁场的响应Figure 1 The response of magnetic particles to an external magnetic field

随后，对成像区域进行空间编码[13]，通过成像方法建立接收线圈感应产生的MPI信号与各个编码点粒子浓度值间的关系，是实现磁性纳米粒子成像的关键。

1.2 系统矩阵成像方法

为了实现对粒子分布状况的精准成像，采用SM成像方法[14]。在SM成像方法中，接收线圈感应产生的MPI信号u(t)与各编码点rn(n=1,…,N)处的粒子浓度信息在时域的关系可用式(1)表示：

(1)

式中：M(rn,t)为位于位置rn处的粒子产生的磁化响应信号；s(rn)为线圈灵敏度；Δv代表rn处粒子对应的体积，则：

(2)

用矩阵-向量形式表示线性系统方程组(1)如式(3)所示：

G×c=u

(3)

式中：u是接收线圈感应产生的MPI信号；c是待测的磁性纳米粒子的浓度；G是SM。对等式两侧的数据进行傅里叶变换，就可以得到方程组在频域的矩阵-向量形式：

(4)

展开表示为：

(5)

式中：m代表频率分量的个数；n代表成像区域中编码点的个数。通过求解式(5)所示的矩阵方程组即可获得各个编码点的粒子浓度信息，进而成像。

1.3 基于ART迭代的矩阵方程求解算法

在实际的MPI中，通常检测不同断层血管中粒子的浓度分布状况进行成像。对于单一断层而言，首先，总是存在粒子浓度为零的组织部分；其次，血管边界使得粒子在血管中的分布相对均匀，这样断层上的粒子分布图像就具有分段恒定区域。在对图像恒定区域精确成像前提下实现粒子分布边缘区域的良好成像，是实现高质量成像的关键，这也对矩阵方程(5)的求解算法提出了更高的要求。

在理想的MPI电磁系统中，基于定点迭代的经典ART算法是求解矩阵方程(5)的快速有效方法。而实际的MPI系统中，通常结合正则化运算对矩阵方程(5)进行扩展来减小外界噪声对系统的影响，从而实现粒子浓度分布的良好成像[15]。在含噪声的一维MPI电磁系统中，对图2(a)所示的粒子分布模型进行成像，采用经典ART算法求解扩展后的矩阵方程所成图像如图2(b)所示。

图2 一维MPI粒子浓度分布对比图Figure 2 Contrast diagram of one-dimensional MPI particle concentration distribution

尽管已采用正则化运算降低噪声的影响，但经典ART算法求得的粒子浓度值相较真实粒子浓度值仍然存在一定误差，特别是粒子分布边界附近的粒子浓度差异，弱化了粒子分布边缘处的非平稳特征，图像相应地变得模糊。

2 基于小波稀疏的MPI算法

2.1 小波稀疏处理

通过简单运算增强MPI中被弱化的非平稳特征是实现MPI快速、精准成像的关键。为确保MPI信号的完整性，采用平稳小波变换(stationary wavelet transform，SWT)对MPI信号进行处理，MPI信号中的非平稳特征在SWT后的小波域表现为数值较大的小波系数，而不含粒子区域表现为数值较小的小波系数，对变换后的小波系数做稀疏处理可以强化信号的非平稳特征，提升成像质量。

阈值算子可实现对信号的稀疏处理。软阈值(soft thresholding)、硬阈值(hard thresholding)、NNG(nonnegative garotte thresholding)阈值是3种最常见的阈值算子，3种阈值算子处理后的数据除了在阈值处及高于阈值部分产生一定差异性外，均可实现对信号的稀疏处理。SWT和阈值算子结合的小波稀疏处理可以降低MPI信号中干扰信号的强度，从而增强粒子分布图像中被弱化的非平稳特征。

以图2(b)中采用ART算法求得的一维MPI信号为例，采用不同阈值算子下的小波稀疏处理对该信号进行处理，得到不同状况下的重构信号如图3所示。

图 3 小波稀疏处理前后效果对比Figure 3 Comparison of effects before and after wavelet sparse processing

3种阈值算子下的小波稀疏处理均可实现对MPI信号的稀疏处理。不同的是，基于软阈值的小波稀疏处理重构信号中粒子浓度数值明显小于真实值，信号失真明显；基于NNG阈值的小波稀疏处理重构信号在粒子浓度数值稍微减小的基础上，实现对非平稳特征的良好重构；基于硬阈值的小波稀疏处理虽然可良好实现对粒子分布状况成像，但是在粒子分布边界处会出现振荡，从而影响重构信号的精度，据此本文采用NNG阈值算子进行稀疏处理。综上，基于NNG阈值的小波稀疏处理可以降低ART解中干扰信息的强度，进而实现对存在非平稳特征的MPI粒子信号的优化。

2.2 基于小波稀疏的MPI算法

将小波稀疏处理应用于经典ART算法的迭代过程中，可在每次迭代后强化所成图像中的粒子分布边界非平稳特征，从而有效解决粒子分布图像中的边界模糊问题。在这种情况下：令φ：Rn→RJ×n表示分解级数为J的平稳小波变换，MPI中粒子浓度的重构问题[式(5)]可以被转化为式(6)的最小化问题：

(6)

在该最小化问题中，函数f代表通过阈值算子促进平稳小波变换后小波系数φc的稀疏处理过程。

综上，对于式(6)的最小化问题，可应用基于小波稀疏的MPI算法求解，算法实现过程为：首先用正则化运算对线性系统[式(5)]进行扩展，然后采用经典ART算法对扩展后的线性系统进行求解，之后对求解得到的粒子浓度c执行SWT，从而得到粒子浓度c的近似分量和细节分量，然后用阈值算子对小波变换后的各分量进行稀疏处理，再将阈值处理后的各分量通过小波反变换得到新的粒子浓度c，将所得到的粒子浓度重新带入到算法中进行下次迭代过程，重复迭代直到所得到的粒子浓度满足某种迭代终止条件算法终止。单次迭代过程可用式(7)表示：

(7)

式中：gi表示系统矩阵G的行向量；α为正则化参数。基于小波稀疏的MPI算法伪代码如下所示：

通过以上的迭代过程，可以实现对待测区域中各编码点的粒子浓度求解。需要注意的是，粒子浓度值是一个非负的值，这要求对迭代后的粒子浓度数值进行非负运算。在提出的基于小波稀疏MPI算法中，算法的收敛性与系统噪声水平有关，噪声水平越大，算法收敛性越差。

3 仿真实验与结果分析

3.1 仿真实验

采用图4所示开放式线性零磁场MPI电磁系统模型进行仿真，获取MPI信号。

图4 线性零磁场MPI电磁系统Figure 4 MPI electromagnetic system of linear zero magnetic field

在线性零磁场MPI电磁系统中：采用上下各2个矩形线圈组成的双平面梯度线圈结构组成一组选择场线圈，两组选择场线圈共同作用就可实现零场线在待成像平面的自由旋转；采用同样由矩形线圈组成的双平面梯度线圈组成聚焦场线圈，主要作用是控制选择场产生的零场线在待测平面平行移动。对选择场和聚焦场同时施加不同的电流激励，就可使零场线按照不同的扫描轨迹对待测区域进行扫描，本文采用成像效果良好的径向轨迹(图5)对待测区域扫描。亥姆霍兹结构的驱动场线圈在25 kHz激励频率的电流作用下可在成像区域z方向上产生磁通密度振幅为2.6 T的均匀变化磁场，该变化磁场可使所成图像的分辨率达到1 mm以下[16]。

图5 零场线径向扫描轨迹Figure 5 Radial scan trace of free field line

系统搭建完毕后，对两个接收线圈中间面积为3.1 cm×3.1 cm大小的成像区域按照1 mm的上下间距进行空间编码，得到31×31共961个编码点。通过模型与计算相结合的方法确定各编码点单独存在粒子时在接收线圈中感应产生的MPI信号，提取各编码点处粒子产生MPI信号频谱的1 000个频率分量得到大小为961×1 000的系统矩阵。

实际MPI成像过程中，受外界噪声的影响，接收线圈感应到的MPI信号与真实MPI信号间存在随机性的偏差，该偏差可用高斯白噪声表示。添加高斯白噪声的MPI信号与SM组成矩阵方程，求解矩阵方程即可得到各编码点的粒子浓度信息，进而成像。

所成图像中的粒子浓度值与真实粒子浓度值间的差异可用峰值信噪比参数(peak signal-to-noise ratio，PSNR)、图形相似性参数(structural similarity，SSIM)来衡量，所成图像中的粒子浓度值与真实粒子浓度值间的差异性越小，峰值信噪比参数值越大、图形相似性参数值也越大。两个参数的计算方法由式(8)和式(9)决定：

式中：x代表理想的粒子分布图像；y代表通过算法所成图像；xmax是x中最大像素值；μx、μy分别表示x、y中各像素值的均值；δx、δy分别表示x、y中各像素值的方差；δxy表示x和y像素值间的协方差；c1和c2是两个常数，主要的作用是避免式子中的分母为零。

本文仿真实验分别采用图6所示的几何形粒子分布模型和血管模型进行成像实验，其中点状编码点对应区域的粒子浓度为1 mol/L，“x”状编码点对应区域的粒子浓度为0.5 mol/L。

图6 粒子分布模型示意图Figure 6 Schematic diagram of particle distribution model

实验过程中设置算法的最大迭代次数为5 000，当两个连续迭代之间的粒子浓度间的差异参数er[式(10)]小于10-5时，迭代终止。

(10)

在相同的终止条件下，绘制出PSNR随阈值λ的变化曲线，当PSNR达到最大值时，选择对应的阈值λ为最终的阈值参数值。算法实现过程中采用基于Haar小波的2级SWT对迭代变量进行处理。

3.2 结果分析

在信噪比分别为30 dB、20 dB、10 dB的噪声下，采用基于小波稀疏的MPI算法和经典ART算法分别对图6所示几何形粒子分布模型和血管模型进行成像，得到的成像结果如图7和图8所示。

图7 基于小波稀疏的MPI算法和经典ART算法成像结果(几何模型)对比Figure 7 Comparison of imaging results (geometric model) between wavelet sparse MPI algorithm and classic ART algorithm

图8 基于小波稀疏的MPI算法和经典ART算法成像结果(血管模型)对比Figure 8 Comparison of imaging results (blood vessel model) between wavelet sparse MPI algorithm and classic ART algorithm

对代表两种算法所成图像的质量差异参数值进行显示，如图9所示。

图9 基于小波稀疏的MPI算法和经典ART算法成像质量对比Figure 9 Comparison of imaging quality between MPI algorithm based on wavelet sparseness and classic ART algorithm

横向对比发现，两种成像模型下，基于小波稀疏的MPI算法所成图像均表现出更好的成像质量。由于阈值算子统一将阈值以下的SWT各分量置零，结果在每次迭代结束后都对得到的粒子浓度信息进行稀疏处理，不含粒子编码点对应的粒子浓度值被有效降低，无粒子区域与粒子均匀分布区域间的浓度梯度相应地增加，实现了对MPI所成图像中粒子分布边界非平稳特征的良好成像，所以无论采用哪种阈值算子，基于小波稀疏的MPI算法所成图像都表现出更好的成像质量。

纵向对比发现，两种成像模型下，当系统信噪比为30 dB、20 dB、10 dB时，基于小波稀疏的MPI算法与经典ART算法相比所成图像的PSNR参数平均提升了67.83%、18.66%、8.05%。在基于小波稀疏的MPI算法中，小波稀疏处理对图像质量的提升幅度取决于阈值算子对SWT各分量的稀疏程度，随着噪声水平的增加，部分不含粒子编码点处对应的分量值接近于存在粒子编码点处对应的分量，这时小波稀疏处理也无法将这些过大的干扰信号有效除去，所以图像质量提升的幅度会随之下降。

分别以几何形粒子分布模型和血管模型为研究对象，在不同噪声水平下对基于小波稀疏的MPI算法和经典ART算法收敛过程进行分析，得到的图像如图10所示。

图10 基于小波稀疏的MPI算法和经典ART算法收敛过程Figure 10 The convergence process diagram of the MPI algorithm based on wavelet sparseness and classic ART algorithm

对比发现，基于小波稀疏的MPI算法收敛速度相较经典ART算法较慢，成像时间会有所增加。在基于小波稀疏的MPI算法迭代过程中，单次迭代所得的粒子浓度信号是由上次迭代所得粒子浓度信号经过ART迭代与小波稀疏处理后产生，相邻的两次迭代过程所得粒子浓度信号之间的差异参数相较经典ART算法更大，所以当设置的终止条件一致时，基于小波稀疏的MPI算法达到稳定就需要更多次的迭代，但是小波稀疏处理对粒子信号的优化作用是有限的，所以两种算法迭代次数的差异不是特别明显。这样，基于小波稀疏的MPI算法就可在稍微增加成像时间的前提下实现更高质量的成像。

4 讨论与总结

本文为了实现MPI中粒子浓度空间分布的快速精准成像，提出在耗时较短的经典ART算法每次迭代结束后对当前解进行小波稀疏处理的MPI算法。在小波稀疏处理过程中，SWT在保证信号完整性的前提下可实现对粒子分布边缘非平稳特征的良好提取，阈值算子对SWT后各分量的稀疏处理在确保含粒子区域粒子浓度值变化不大的前提下，有效降低了不含粒子区域的干扰信息强度，无粒子区域与粒子均匀分布区域间的浓度梯度相应增加，实现了对MPI所成图像中粒子分布边界非平稳特征的良好成像。

在提出的基于小波稀疏的MPI算法中，将小波稀疏处理融入到经典ART迭代过程中，可在每次迭代后耗费极少时间强化所成图像中的粒子分布边缘非平稳特征，当算法达到收敛时，计算得到的各编码点粒子浓度值就更接近真实粒子浓度值。算法收敛速度与系统中的噪声水平有关，系统噪声水平越低，算法达到收敛所需的时间就越少，所成图像的质量也越高。综上，在较低的噪声水平下，基于小波稀疏的MPI算法可以在短时间内实现对粒子分布状况的高质量成像。

在今后的工作中，将着重研究基于小波稀疏的MPI算法与MPI成像系统的兼容性，从而推动本文方法在MPI实时精准成像方面的实际应用。