具有噪声鲁棒性的三维磁性纳米粒子成像快速重构方法

史玥婷，任仕伟，王晓华

（1. 北京理工大学 集成电路与电子学院， 北京 100081；2. 北京理工大学 长三角研究院(嘉兴)， 浙江， 嘉兴314019；3. 北京理工大学 重庆微电子中心， 重庆 401332）

磁性纳米粒子成像（magnetic particle imaging, MPI）是一种全新的示踪剂成像技术[1]，在血管造影[2]、肿瘤成像[3]、体细胞追踪[4]和磁热疗[5]等领域具有重要的应用前景. 该技术利用超顺磁性纳米粒子对变化磁场会产生非线性响应的原理，完成对磁粒子浓度空间分布的成像. 如何快速且高质量地从接收到磁性粒子电信号重构出的MPI 图像，一直是国内外研究的热点与难点.

MPI 重构的主要方法可以分为系统矩阵法和XSpace 重构法. 系统矩阵重构算法[1]是将MPI 接收线圈检测到感应信号与磁性纳米粒子示踪剂（super paramagnetic iron Oxide，SPIO）空间分布之间的映射关系描述为一个系统矩阵，然后再利用校验好的系统矩阵对待检测的图像进行重构. 该方法虽然可以准确地描述系统特性，但校准数据采集时间较长、计算复杂度高和对内存要求高. X-Space 重构算法的主要思想是使编码直接在接收到信号的空间中进行[6]，直接快速的在信号接收空间完成空间编码并重建出MPI 图像，可适用于一维至多维度MPI 成像和重构.在进行三维MPI 成像时，通常需要获取符合奈奎斯特采样定律要求的投影，再通过滤波反投影等方法进行重构.

根据MPI 基本原理，可以通过提高硬件性能的方法来提升成像速度，比如提高梯度磁场强度等等，但由于硬件的成本过高和磁场强度大小的限制，导致亟需在成像和重构算法上进行研究改进.DONOHO[7]提出的压缩感知理论表明信号具有稀疏性，当信号在某个变换域下表示为稀疏时，通过求解一个非线性最优化问题，可以以远低于奈奎斯特采样定理要求的采样数据来重建信号. 随后，国内外研究人员相继提出可用于CT、MRI 医学图像的压缩感知重构算法[8−9]. 就稀疏变换和重构算法而言，现有算法在充分利用图像的几何正则性、选取最优和最稀疏表达、重构图像计算复杂等方面还存在进一步优化空间. 此外，这些方法没有考虑到MPI 的成像原理，无法直接用于MPI 图像的稀疏采样与重构. 2013年，KONKLE 等[10]提出了基于X-Space 理论的零磁场点磁性纳米粒子成像（field free point MPI，FFP-MPI）直流信号恢复凸优化的重构算法，引入了MPI 图像的连续性和非负性两个先验条件，为进一步研究零磁场线磁性纳米粒子成像（field free line MPI ，FFLMPI）凸优化重构算法提供了新的方向.

本文通过设计基于X-Space 理论的磁性纳米粒子成像稀疏采样与重构方法，提出具有噪声鲁棒性的MPI 成像模型正则优化函数，实现3D TV 稀疏算子及低运算复杂度重构的软件程序框架，充分利用MPI 先验信息解决传统方法对成像投影数量完备性的要求，大幅提升MPI 三维成像和重构速度，为今后MPI 三维成像在各种生命科学及应用研究中发挥更好作用提供了重要基础.

1 磁性纳米粒子成像模型

1.1 点扩散函数及基础成像模型

磁性纳米粒子成像系统与常见的医学成像方法同样具有线性移不变（linear shift invariant, LSI）特性，MPI 图像可表示为纳米粒子浓度与MPI 成像点扩散函数（point spread function, PSF）的卷积[6,11]. MPI 的点扩散函数为Langevin 响应的导数形式，也可以称为成像模糊方程，MPI 的分辨率由PSF 的宽度决定.

在MPI 系统中，零磁场区域（field free region, FFR）是由两个相向的磁场创建产生的，由一个偏置磁场驱动FFR 在视场内按设定的轨迹进行扫描，还有一个高频正弦波激励磁场用于激发SPIO 产生信号. 当FFR 经过SPIO 位置时，检测线圈即可以采集到SPIO 磁矩方向翻转的信号，由于SPIO 信号在高频谐波部分，因此需要滤除基频信号，然后使用滤波后的磁粒子信号进行图像重建.

施加一个强磁梯度为 −µ0G的磁场，那么待成像样品中的SPIO 处于的磁场大小为

为了刺激粒子翻转，需要对样品再施加一个激励磁场Hs，该激励磁场一般为正弦激励场. 通过求解H(x)+Hs(t)=0，得到瞬时场FFP 信号为

如前所述，MPI 信号是通过接收线圈的磁通 ϕ变化的结果. 这种变化是由于纳米粒子的总磁化强度M的变化引起的. 在一维情况下，可以得出[6]

其中，B1只是一个常数，它模拟了接收线圈的灵敏度.朗之万理论用来描述SPIO 的磁化强度是其所受外加磁场的函数，磁化强度的大小可以看作SPIO 粒子中与外加磁场对齐的原子偶极子的数量. 如果所有的原子偶极子都对齐，则称之为磁饱和. 朗之万理论认为磁场H中单点x的磁化强度为

由式(2)和式(3)可进一步推导出

式中： γ ≜B1mG/Hsat；m[Am2]为磁性纳米粒子的磁矩；Hsat[A/m]为使纳米颗粒示踪物半饱和所需的磁场振幅；s(t) [m/s]为瞬时FFP 速度； L为表征SPIO 磁化特性的Langevin 函数.

利用X-Space 重构方法可以将时域信号方程直接解析地转换成本地MPI 图像. 这只需要将接收到的信号按瞬时速度归一化，然后网格化到FFP 的瞬时位置.

1.2 零磁场线投影成像及重构

在本文实验平台上，FFL 是沿y轴方向驱动场激励矢量沿z轴方向，移位矢量沿x轴方向，扫描系统[x y z]T是 非旋转框架，样品成像平 台 [x′y′z′]T是旋转平台. 为了获得投影重建所需的多个投影图像，将待成像样本沿z轴轴心旋转，这等效于FFL 磁场绕z轴旋转.

沿y轴线方向的FFL 磁场可以用梯度矩阵H(x)=Gx来 描述. FFL 由均匀移位磁场Hs=Hx控制在x轴方向移动. 粒子沿着z轴被激励场He=Hz激发. 这里同样假设SPIO 粒子响应是瞬时的，没有弛豫效应，叠加后的磁场可表示为[12]

式中：x为 扫描仪坐标系上的位置向量；Gxx=−Gzz为满足麦克斯韦方程组.

要确定FFL 在空间中的位置，需要计算磁场的幅值（为简化结果，公式中为该幅值的平方值）:

其中， θ为旋转角度. 对式（8）进行整理得到零磁场的条件是：z′=−Hz/Gzz和x′cosθ+y′sinθ=Hx/Gzz. 定义投影线:

由此，求解出FFL 处于

并由Hz激发场和梯度场决定其在z轴的位置. 注意，物理上l表示垂直于FFL 的移位位置. SPIO 信号沿着这条线进行积分，从而可以任何旋转角度获取样品的投影图像. 通过调节变量Hx,Hz和 θ以精确控制并收集一组完整的投影数据.

2 磁性纳米粒子成像稀疏采样及三维重构算法

2.1 模型假设

不同于系统矩阵法，X-Space 重构算法的主要思想是在接收到信号的空间中进行直接编码. MPI 需要强大的磁场梯度使FFP 以外所有位置的粒子处于饱和状态，进而进行空间编码. 换而言之，X-Space 重构依赖于3 个主要假设：1) 磁性纳米粒子绝热地与所施加的磁场对齐；2) FFP 位置是唯一的；3) 由于对基波进行滤波而导致的低频信息丢失是可恢复的[1,6].除此之外，MPI 图像还具有连续性和非负性两个先验条件[9].

2.2 3D TV 稀疏约束算子

由于医学图像一般具有平滑性，即图像多具有接近的灰度分布仅在某些结构的边缘产生灰度急剧变化，TV 稀疏的约束可以看作是一种自适应的正则化形式，3D TV 准则可表示为

其中， D1x、 D2x、 D3x分别为沿i方向、j方向以及k方向的差分算子，约束并在保障三维图像边缘的稀疏度，对于位于远离边缘的平滑区域的体素，强正则化可以滤除振荡和高频分量. 但是对于位于边缘的体素，TV 正则化滤除很弱. 因此，对于边缘上的体素，高频分量不会被过滤掉，在重建过程中保留了边缘.与l2正则化相比，TV 正则化的这种边缘保持特性实现了更高的空间分辨率[13]，且具有较好的噪声抑制能力.

为了实现快速重构，本文设计了稀疏矩阵形式的3D TV 变换函数，可以通过Matrix Free 方法和一阶快速邻近梯度重构算法软件实现. 稀疏3D TV 算子的示意图如图1 所示，分别为x、y、z三个维度和完整3D 算子的可视化图，黑色和白色位置为进行全变分计算应对的稀疏点位置.

图1 稀疏3D TV 算子示意图Fig. 1 Diagram of the sparse 3D TV operator

2.3 具有噪声鲁棒性的稀疏采样与重构算法

本文将MPI 投影重构转换为凸优化重构形式，3D NRSS-MPI 优化问题可以描述为

在MPI 中，函数g(x)可以是空操作，也可以选择利用一些先验信息，如图像连续性等. 在本文3D NRSS-MPI 算法中，为了实现对MPI 含噪数据进行优化重构，设置了正则化函数为

其中，TV 范数具有n维度，D 表示一阶有限差分，那么具体的目标函数可定义为

进一步通过任意邻近梯度算法[14]求解重构图像:

换而言之，NRSS-MPI 方法解决了一个凸优化问题，试图最小化两个数据一致性项惩罚（模拟扫描和采集数据之间）和全变分稀疏域惩罚，其中MPI 成像模型用Matrix Free 线性变换表示. 在进行三维MPI成像时，根据如图2 所示，本文方法提供了基于广义Radon 变换的MPI 正向模型和3D TV 稀疏算子，可通过邻近梯度优化算法快速实现三维图像重构.

3 模型实验与性能分析

3.1 欠采样下三维MPI 图像重构实验与分析

在磁性纳米粒子成像系统中的重构误差通常表现为条状伪影. 伪影会降低图像的有效分辨率，影响图像的清晰度. 为了避免条纹伪影和分辨率损失，传统方法通常选择符合Nyquist 标准的完备投影采样集，但MPI 扫描和重构所需时间却大幅提升，这限制了许多MPI 应用场景，如活体MPI 成像，在漫长扫描过程中任何移动误差都将严重影响成像质量.

为了分析本文3D NRSS-MPI 算法和传统FBP 算法的性能，以及三维重构结果对空间分辨率的影响，进行了点阵模型实验. 点阵模型一组竖条纹中间间隔为1.5 mm，每组条纹之间间隔为12 mm，两组之间最邻近距离为7.5 mm. 根据奈奎斯特采样定律可以推导出避免混叠所需的投影数量[12]

式中：FOVxy为FOV 最小值，根据待成像实验设定；Kmax[m−1]为系统中最大的空间频率，由系统PSF决定.

通过滤波反投影（filtered back projection，FBP）方法和NRSS-MPI 方法对点阵模型成像并欠采样重构，实验结果如图3 所示，在1/3 完全投影数情况下，FBP 算法重构图像质量较低，有星状伪影，但依旧可以分清两条细竖线. 在1/4 完全投影数情况下，FBP 方法重构图虽然能看清细竖线，但图像质量已经不能用于诊断等后续应用了，NRSS-MPI 算法重构图虽然没有明显星状伪影，但细竖线间存在模糊.

图3 不同采样数据量下本文NRSS-MPI 方法与传统FBP 方法重构结果Fig. 3 The reconstruction results of NRSS-MPI and FBP with different sampling projections

在成像和重构时间方面，由于NRSS-MPI 方法可以以1/4 欠采样投影数进行重构，这意味着获取图像时间缩小了4 倍. 另外在重构的计算上，NRSS-MPI方法花了10 s 重建了一个3D 版本的模型，在实际应用中证明了该算法的快速性.

3.2 噪声鲁棒性实验与分析

为了定量评价重构图像的相似性和分析采样数据含噪声情况下CSMPI 算法的鲁棒性质，用SSIM指标对RSN分别为5、10、20 下的含噪声采样数据通过CS-MPI 方法和FBP 方法重构图进行评价测试.

式中： µx和 µy分 别代表x、y的平均值； σ和 σ分别代表x、y的方差； σxy为协方差；c1和c2为常数，对比结果如表1 所示.

表1 不同含噪声数据在本文NRSS-MPI 方法和传统FBP 方法重构图结构相似度比较表Tab. 1 The comparison of SSIM between the proposed NRSSMPI and FBP with different noise level

从SSIM 水平而言，NRSS-MPI 优于FBP 重建，尤其是在RSN较低的情况下. 当噪声水平增加时，可以观察到NRSS-MPI 和FBP 的SSIM 收敛水平之间的差距增大，提出的CS-MPI 方法在重构图像的结构相似性分析中表现优越.

此外，对不同信噪比下FBP 算法和3D NRSSMPI 算法1D 重构轮廓进行了对比，分析本文NRSSMPI 算法对含噪声采样数据的重构能力，仿真实验结果如图4 所示. 对冠状血管模型的水平位置（y=40）的最强一帧进行1D 轮廓勾画，图中（x=19,x=40）位置是冠状动脉的两条分支血管. 在RSN=20 情况下，3D NRSS-MPI 和FBP 算法总体上都能重构出两处极大值位置，但是在1/4 欠采样时，FBP 方法出现了比较明显的尾部波动，即产生了背景伪影，而3D NRSSMPI 方法依然可以准确重构. 在RSN=5 的采样噪声较大的情况下，FBP 在不同采样投影数的实验中都存在比较明显的波动，3D NRSS-MPI 算法明显比FBP稳定、准确. 1D 轮廓的重构比较更直观的看出FBP算法在背景处噪声波动大于3D SS-MPI 算法，尤其在含噪声数据重构实验时FBP 算法重构数据与真值存在明显差距. 再一次验证了3D NRSS-MPI 算法能够快速、稳健地重建欠采样、噪声较大的FFL-MPI数据.

图4 不同含噪数据下本文NRSS-MPI 方法与传统FBP 方法重构结果1D 轮廓图Fig. 4 The comparison of 1D profile of the reconstruction results between the proposed NRSS-MPI and FBP with different noise level

3.3 三维成像与重构实验

用未经稀释的Resovist 纳米粒子示踪剂制作了5 点源模型，然后将它们固定在设计定制的3D 支架中，使得在3D 点源模型在x、y、z视角方向具有不同间隔值. 此外，设计并用激光切割机（universal laser systems, Inc., Scottsdale, AZ, USA）制作雕刻了3 mm厚的点阵模型和冠状动脉血管样待测模型，用胶带进行密封用针头注入稀释后的Resovist 纳米粒子示踪剂，针孔处用透明指甲油密封防止SPIO 蒸发或流出污染扫描平台. 用3D 打印的模型支架将血管模型呈60°角放置形成3D 模型，如图5（a）（c）所示. 最终待测模型通过3D 打印模型底托准确固定在成像孔，并加载到高分辨率FFL-MPI 扫描仪中，成像结果如图5（b）（d）所示，实验证明了本文算法在实际应用中是快速有效的.

图5 三维成像模型及重构图Fig. 5 3D Phantoms and the reconstruction results

4 结 论

本文构建了FFL-MPI 前向模型并提出了一种新的具有噪声鲁棒性的稀疏采样MPI 三维重构算法，该算法通过最小化模拟扫描和采集数据之间的差异性和3D 全变分域上的差异，然后通过邻近优化算法求解重构出MPI 三维图像. 所提出的MPI 成像及重构模型和快速稀疏的3D TV 算子充分利用了MPI 成像的先验信息，具有良好的噪声鲁棒性. 通过成像实验和含噪数据鲁棒性关键指标测试分析了滤波反投影FBP 方法和本文NRSS-MPI 方法在不同投影数下的重构图像质量，结果表明，当FBP 对约1/4 欠采样投影进行图像重构时会产生明显的伪影，但NRSSMPI 方法能正确地重建冠状动脉体模. 在不同含噪水平的数据重建中，本文NRSS-MPI 方法重构结果的结构相似性优于传统方法. 此外，还通过三维成像实验证明了本文算法在实际应用中是快速有效的.