一种基于改进多光谱的稀疏角CB-XLCT成像方法

张文元，海琳琦，刘英杰，张海波

(1.西安航空职业技术学院航空维修工程学院，陕西 西安 710089； 2.西北大学信息科学与技术学院，陕西 西安 710127)

0 引 言

X射线发光断层成像(X-ray Luminescence Computed Tomography, XLCT)作为一种新的三维分子双模成像技术，已成为生物光学与医用成像领域的前沿热点[1]。结合纳米发光材料的光学特性，利用CT技术中的X射线源激发，发光材料受照射会发出近红外光，进而借助高灵敏度电耦合器件(Charge Coupled Device, CCD)采集透射出生物体表的微弱光信号，基于逆问题建模与成像算法可在分子层面对早期肿瘤进行有效检测[2]。

根据CT中X射线源的激发方式，XLCT技术可分为笔束型XLCT(Pencil-Beam XLCT, PB-XLCT)、扇形束XLCT(Fan-Beam XLCT, FB-XLCT)以及锥束型XLCT(Cone-Beam XLCT, CB-XLCT)[1,3-5]。其中，PB-XLCT和FB-XLCT分别利用笔束形X射线源与扇形X射线源各自的特点，利用“激发先验”可获得较高的成像分辨率，但由于系统耗时过长，不利于肿瘤的快速成像。随后，CB-XLCT被提出，其利用锥束型射线源扫描时间短和X剂量利用率高的优点，为XLCT应用于早期肿瘤的实时成像提供了可行性[6-11]。

类似于传统的生物发光断层成像(BioLuminescence Tomography, BLT)[12-13]和荧光分子断层成像(Fluorescence Molecular Tomography, FMT)[14-15]，XLCT成像问题在数学上属于不适定逆问题[16]。BLT和FMT等经典的重建方法大都可直接应用于传统CB-XLCT成像中。Chen等人[17]提出一种结合靶标深度位置的自适应Bregman算法，在算法中融入X射线衰减先验信息，有效降低了成像逆问题病态性。Pu等人[11]将主成分分析应用于CB-XLCT成像，通过实验验证了所提方法的有效性。Zhang等人[5]提出了一种基于高斯马尔可夫随机场的贝叶斯成像方法，用于传统CB-XLCT的逆问题求解。曲璇[6]提出了一种基于多网格的可行区域策略，但多网格的多次重建间接增加了计算成本。随后，稀疏角CB-XLCT(Sparse-view CB-XLCT)技术的出现，为CB-XLCT技术的实时成像转化进程提供了进一步促进。然而，由于稀疏角度数据的约束，其成像逆问题病态性相比传统多角度CB-XLCT成像(一般24角度)的逆问题病态性更为严重，直接影响了传统重建方法的成像性能[18]。

为了解决此问题，本文提出一种改进多光谱策略的稀疏角CB-XLCT成像方法。首先，将经典的多光谱成像策略应用于稀疏角CB-XLCT成像中，构建高维系统矩阵并建模逆问题；其次，考虑高维系统矩阵的规模对成像质量的影响，基于谱回归方法对系统矩阵进行子空间学习，优化成新的系统矩阵；再次，对新的系统矩阵进行矩阵预处理并建模新的成像逆问题；最后，采用经典的成像算法进行逆问题的准确求解。本文分别设计稀疏角重建实验、多角度重建实验以及鲁棒性实验，测试所提方法的有效性和鲁棒性。

1 理论方法

1.1 CB-XLCT成像模型

CB-XLCT的前向模型分为3个部分：

1)锥束X射线源发出的X射线在生物组织中的传输，该过程可用下式描述：

(1)

其中，X(r)表示生物组织中位置r处的X射线强度，X0为X射线源的初始强度，μx为是X射线穿过的组织的衰减系数[1]。

2)X射线照射生物体内的纳米靶标，激发其产生近红外光，该过程可表述为：

S(r)=εX(r)ρ(r)

(2)

其中，S(r)表示光源，ρ为目标浓度，ε为发光产额。

3)近红外光子经过生物组织的散射、折射和吸收等过程，传输到生物体表的光子信号可被高灵敏CCD探测并采集。经典的扩散近似模型可被用来描述组织中的光传输过程：

(3)

其中，D(r)=(3(μa(r)+μ′s(r)))-1为散射系数，μa(r)和μ′s(r)分别为吸收系数和散射系数。Ω和∂Ω分别表示成像区域和区域边界。Φ(r)表示光子通量，ν表示单位法向量，γ的值为依赖光折射系数的边界不匹配因子。

基于有限元理论[7-11]，式(3)可被离散为如下的矩阵形式：

Φ=M-1·F·ε·X·ρ

(4)

其中，F代表光源权重矩阵，M为正定矩阵。

令A=M-1·F·ε·X，得到生物体表面的发射光光子通量流率Φ与纳米发光材料三维分布ρ之间的关系,其公式为：

Φ=Αρ

(5)

其中，ρ为待求解的纳米靶标，Φ为采集的生物体表光信号，A为系统矩阵。

1.2 改进多光谱策略逆问题求解

经CCD相机采集的生物体表微弱光学信号，根据经典的多光谱策略可对生物体表微弱光学信号进行滤波处理，将其划分为多个不同波长的子带信号[19]，进而明显增加了生物体表的边界测量数据。因此，引入多光谱策略将式(5)表示为：

ΦM=AMρ

(6)

其中：

(7)

通过式(7)可以看出，采用多光谱策略构建的ΦM直接增加了边界测量数据量，有效缓减了逆问题的病态性。然而，相应系统矩阵AM明显增加了原系统矩阵A的矩阵规模。为了保证成像质量同时减少高维矩阵规模的影响，本文借鉴流形学习的思想，基于谱回归(Spectral Regression, SR)方法对系统矩阵AM进行子空间学习。本质上，SR方法来源于一个通用的谱嵌入框架[20-21]，基于该框架可概括为如下4步来获得AM的低维投影矩阵R：

Step1构建权重矩阵Kij与近邻图G。Kij表示节点之间的相邻关系。近邻图G来表示AM的内在高维几何结构，将高维系统矩阵AM中的每一行看成是一个多元统计变量:

(8)

其中，节点i和j分别表示原始数据集中的行向量xi和xj。如果节点i和j属于同一个类标签，将其用边连接，同时赋以权重1/lk表示边的权值，lk为相应标签内样本数量[20]。

Step2根据文献[22]构建拉普拉斯矩阵L=D-K,其中D为一个基于权重矩阵Kij构造的对角矩阵且Dii=∑jKji。

Step3利用公式(9),计算特征向量y的值。

(9)

其中，T表示矩阵转置。

Step4计算投影矩阵R=[r1,r2,…,rc-1](c

rk=(XXT+αI)-1Xyk

(10)

其中X=[x1,x2,…,xm]T，y=[y1,y2,…,ym]T。

基于上述3步计算获得的投影矩阵R，由于是对高维系统矩阵AM的子空间投影，对初始逆问题式(6)进行优化可得：

(11)

随后，基于稀疏感知理论[22-23]，为了进一步保证逆问题成像的质量，对式(11)进行矩阵预处理(preconditioning)[10]，随后具体步骤如下：

(12)

2)设置预处理矩阵R：

R=(ΛΛT+λI)-1/2UT

(13)

其中，I为单位矩阵，λ为正则参数(λ>0),T表示对矩阵进行转置操作。

3)对公式(11)进行矩阵预处理：

(14)

4)可获得优化后的新的逆问题如下所示：

(15)

其中，σ为正则化参数。对于式(15)，基于贪婪方法可采用经典的正交匹配追踪(Orthogonal Matching Pursuit, OMP)算法进行求解[24]。

2 实验与结果

2.1 实验设置

为了验证本文所提方法的有效性和鲁棒性，借鉴文献[9]的工作，分别仿真模拟了稀疏角度4投影和全角度24投影的成像对比实验。仿真数据采用的是非匀质仿体模型，仿体半径为10 mm，高为30 mm。该仿体模性包含了6个器官，分别是心脏、左肺、右肺、肝脏、骨骼以及肌肉。在肺部置入1个半径为0.5 mm、高为1 mm的光源来模拟早期肿瘤目标，光源距离X轴方向最近表面深度为4 mm。

在具体对比中，分别引入经典的主成分分析(Principal Components Analysis, PCA)以及不采用本文方法直接重建进行成像对比，最终的重建算法均为OMP算法。在结果量化分析中，引入中心定位误差(Location Error, LE)[2]、恢复浓度(Relative Quantity, RQ)[7]以及相对定量误差(Relative Quantity Error, RQE)[7]来作为量化评价指标。在重建中，为了避免经验选取正则参数带来的主观影响，本文采用经典L曲线方法自适应选取正则参数[1]。所有仿真实验都是基于Matlab平台，在配置为3.60 GHz Intel® Xeon® Processor i7-7820X处理器和64 GB内存的台式工作站上进行计算。

2.2 稀疏角度与全角度重建实验对比

在前向仿真中，本节分别模拟稀疏角度4次X射线激发与全角度24次X射线激发，采用XLCT领域中5种常见的子带波长(585 nm, 627 nm, 667 nm, 720 nm 和 802 nm)进行逆向重建[1-2,7-10,25-29]，不同子带波长下的各个器官的光学参数根据Alexandrakis等人[30]的工作进行计算。仿真实验纳米发光材料模拟NaGdF4: Eu3，相应的X射线的衰减系数为0.054 mm-1[31]。目标的纳米发光产额设置为0.15 cm3/mg[7]。纳米靶标颗粒浓度设置0.3183μg/mm3[7]。经过前向仿真计算，对于单光源数字仿体实验的5种子带波长(585 nm, 627 nm, 667 nm, 720 nm 和 802 nm)对应各自的能量权重系数分别为0.06, 0.15, 0.22, 0.26和0.31[29]。图1为采用的数字仿体以及2种激发方案下各自前向仿真结果展示。基于有限元方法，前向仿真采用的仿体被离散为167663个四面体单元和29597个结点。

图1 实验采用的仿体以及2种激发方案下各自前向仿真结果展示

在逆向多光谱重建中，分别采用本文所提方法(SR+preconditioning)、结合PCA的方法(PCA+preconditioning)以及直接重建(with no processing)进行对比，采用的仿体包含47873个四面体单元和8674个结点。图2和图3分别是稀疏4角度投影和全24角度投影各自的重建结果展示。表1和表2分别为相应的量化重建结果。

表1 稀疏4角度投影下的量化重建结果

表2 全24角度投影下的量化重建结果

图2 稀疏4角度投影下，重建结果在Z=15 mm平面处的二维截面图与重建结果的三维展示图

通过上述结果可以看出，对于全24角度投影重建实验，所有方法的重建误差LE均在1 mm以内，恢复浓度均在0.265 μg/mm3～0.295 μg/mm3，相应的相对量化误差RQE均在20%以内。对于稀疏4角度投影重建实验，由于测量数据的明显不足，使得逆问题病态性明显加剧，直接重建的误差LE已大于1 mm，恢复浓度已降至0.235 μg/mm3，相应的相对量化误差RQE已增至26.26%。而对于本文提出的SR+preconditioning方法与PCA+preconditioning方法，重建的误差LE仍然在1 mm以内，恢复浓度相对于全24角度投影重建结果略有降低，但相应的相对量化误差RQE仍在20%以内。特别地，不论是重建误差LE还是恢复浓度RQ及相应的相对量化误差RQE，本文所提方法的成像质量最好。

2.3 鲁棒性测试

在上一节重建实验结果的基础上，本节针对稀疏4角度投影，进一步设计了2组鲁棒性实验。考虑影响成像质量常见的因素[6-7,19]，分别设计了靶标深度水平测试与噪声水平测试。其中，靶标深度水平测试中设置光源距X轴方向最近表面深度分别为8 mm、7 mm、6 mm、 5 mm和4 mm，噪声水平测试中分别设置噪声为高斯白噪声，该噪声服从均值为0的正态分布，方差以表面测量值均值的5%为间隔在10%到30%之间变换。实验结果如图4和表3、表4所示。

图4 稀疏4角度投影下深度水平测试与噪声水平测试中所提方法各自的重建结果展示

表3 稀疏4角度投影下深度水平测试中所提方法量化重建结果

表4 稀疏4角度投影下噪声水平测试中所提方法量化重建结果

通过重建结果可以看出，尽管目标深度的增加或者噪声水平的增加都会不同程度增加逆问题重建的难度，进而影响成像质量。但本文所提方法在上述因素的测试中，重建误差LE均在1 mm以内。对于深度水平测试，随着目标深度的减小，恢复浓度有所改善，相应的相对量化误差也逐渐减小。对于噪声水平测试，重建误差LE在噪声水平的变化下整体较稳定，而恢复浓度RQ及相应的相对量化误差RQE相对于重建误差LE更易受噪声的干扰。其中，当噪声水平设置为30%时，相对量化误差RQE达到了30.62%。总体上，本文所提方法对于深度因素和噪声因素的影响，成像质量是可接受的。

3 结束语

稀疏角CB-XLCT成像作为一种新型的光学断层成像技术，具有对早期肿瘤进行实时监测的应用潜力。但受限于测量数据的严重不足，相对于传统多角度成像，逆问题病态性明显加剧。本文提出了一种改进多光谱的稀疏角CB-XLCT成像方法，将谱回归与矩阵预处理先后对多光谱策略生成的高维系统矩阵进行子空间学习以及进一步优化重建逆问题，最终采用经典OMP算法对优化后的逆问题进行准确求解。本文设计了稀疏4角度投影仿真实验、全24角度投影仿真实验、深度水平测试以及噪声水平测试。实验结果表明，本文所提方法不仅可以准确定位目标和恢复浓度，同时对于常见干扰因素也具有较好的鲁棒性。

总之，作为医学分子影像的研究热点，CB-XLCT应用于早期肿瘤的实时监测还有很多问题需要解决，比如在成像中减少投影数量进行单视图重建以及在重建中对肿瘤形状进行准确勾勒；此外，对于核谱回归方法在CB-XLCT实时成像中的进一步探索也是笔者未来的研究重点。