单光子发射计算机断层扫描(single photon emission computed tomography，SPECT)^[1-2]是一种重要的分子影像技术，可以通过追踪放射性药物在人体内的分布反映生理代谢情况。目前临床主要使用双探头平行孔SPECT系统。准直器作为SPECT系统的重要部件，可以阻止非特定方向的光子入射到探测器表面，从而对人体内放射性药物进行定位。

对于甲状腺癌的诊断和治疗评估，碘-131 SPECT平片成像是一种常用且有效的技术手段^[3-4]。碘-131衰变纲图如图 1所示，临床SPECT显像主要利用364.5 keV光子。对于采用铅合金和六边形设计的平行孔准直器(见图 2a)，较高能量的Gamma光子容易穿透准直器隔壁到达探测器，从而形成放射状伪影^[5](见图 2b)。该伪影一方面会导致难以区分真实的碘-131富集区域，另一方面可能会掩盖周围较低摄取区域，例如淋巴结转移。因此，消除放射状伪影对于辨别真实的碘-131分布以及开展进一步的临床诊断、治疗评估都具有十分重要的意义。

针对放射状伪影问题，一种常用的解决方法是采用点扩展函数(point spread function，PSF)对准直器穿透效应进行建模，并结合散射校正和先验函数实现解卷积算法，从而消除放射状伪影^[5-6]。先验函数的选择和相应迭代算法的设计都会对最终效果产生影响。传统方法采用全变分先验改善重建结果，但是在迭代算法方面存在一定问题，迭代算法的收敛性无法保证。本文提出了一种新的解卷积算法用于解决碘-131成像中的放射状伪影问题，选取适当的先验函数并推导了单调收敛的迭代算法；对提出的算法，考察了不同大小PSF和不同正则化参数对结果的影响，并基于AMIDE软件^[7]对感兴趣区(region of interest, ROI)进行定量分析；对于优化后的参数，采用临床数据对算法效果进行了进一步评估。

1 SPECT系统简介

本研究主要针对北京永新医疗设备有限公司的NET632 SPECT系统^[8]进行探索。该系统的探测器采用碘化钠(NaI)晶体，固有空间分辨率和能量分辨率分别为3.4 mm和9.8%@140 keV。该系统配备的用于碘-131成像的高能通用准直器参数如表 1所示。

2 伪影消除算法

考虑碘-131成像时的准直器穿透效应和噪声特性，平片图像y和真实图像x之间的关系可表示为

$ \boldsymbol{y} \sim \operatorname{Poisson}(\boldsymbol{x} \otimes \boldsymbol{h}+\boldsymbol{s}). $

(1)

其中：$\otimes$表示二维卷积操作；h表示二维PSF图像，大小为N×N，且满足$\sum\limits_{a, b} \boldsymbol{h}(a, b)=1$，a、b分别表示二维PSF图像行和列的索引；s表示散射数据。因此，可采用极大似然估计方法求解原始图像。

目标函数L(x)可表示为

$ L(\boldsymbol{x})=\sum\limits_i y_i \ln \left[(\boldsymbol{x} \otimes \boldsymbol{h})_i+s_i\right]-\left[(\boldsymbol{x} \otimes \boldsymbol{h})_i+s_i\right] . $

(2)

其中: i表示像素编号, y_i、s_i分别为平片图像y、散射数据s中的元素。

此外，为了压低噪声并保持图像边界清晰，改善重建图像质量，引入先验函数

$ U(\boldsymbol{x})=\frac{1}{4} \sum\limits_j \sum\limits_{k \in N_j} \omega_{j k} \sqrt{\left(x_j-x_k\right)^2+\delta^2} . $

(3)

其中：N_j表示像素j的邻域，本研究选取邻域大小为5×5；ω_jk表示j、k 2个像素之间的相关性，通常为像素距离的倒数; x_j、x_k分别表示真实图像x分别在像素j、k处的值；δ表示超参数，本研究默认δ=0.001。

综上所述，最终基于目标函数(4)对原始图像进行重建，具体可表示为

$ \psi(\boldsymbol{x})=L(\boldsymbol{x})-\beta U(\boldsymbol{x}) . $

(4)

其中：β表示正则化参数，用于控制正则化强度。

2.1 迭代算法

对于目标函数(4)采用迭代算法求极大值。为保证迭代过程单调收敛，采用优化转移的方法^[9-10]，寻找目标函数的替代函数并推导迭代公式。L(x)的替代函数可表示为：

$ Q_L\left(\boldsymbol{x} ; \boldsymbol{x}^n\right)=\sum\limits_j p_j \left( {\hat x_{j, {\rm{EM}}}^{n + 1}{\rm{ln}}{x_j} - {x_j}} \right), $

(5)

$ \hat{x}_{j, \mathrm{EM}}^{n+1}=\frac{x_j^n}{p_j}\left[\boldsymbol{h}^{\mathrm{T}} \otimes \bf{b} \bf{p}\right]_j, $

(6)

$ \bf{b p}_i=\frac{y_i}{[\boldsymbol{h} \otimes \boldsymbol{x}+\boldsymbol{s}]_i}, $

(7)

$ \boldsymbol{h}^{\mathrm{T}}(a, b)=\boldsymbol{h}(N-1-a, N-1-b), $

(8)

$ \boldsymbol{p}=\boldsymbol{h}^{\mathrm{T}} \otimes \boldsymbol{I} . $

(9)

其中：x_jⁿ表示第n轮迭代时得到的像素j的值；$\hat x_{j, {\rm{EM}}}^{n + 1}$表示第(n+1)轮迭代时根据式(6)期望最大化(expectation maximization, EM)算法更新得到的像素j的估计值；xⁿ表示第n轮迭代重建的图像；T表示翻转，h^T表示二维PSF图像h(沿横、纵方向)翻转的结果；bp表示用于反投影的图像矩阵；bp_i表示bp的第i个元素，是真实投影数据和当前模型计算结果的比值，用于对当前图像估计值进行修正；p_j表示卷积模型中像素j处的灵敏度因子，p表示灵敏度因子构成的图像；I表示所有元素均为1的图像。

U(x)的替代函数可表示为：

$ Q_U\left(\boldsymbol{x} ; \boldsymbol{x}^n\right)=\frac{1}{2} \sum\limits_j w_j^n\left(x_j-\hat{x}_{j, \text { reg }}^{n+1}\right)^2, $

(10)

$ w_j^n=\sum\limits_{k \in N_j} w_{j, k}^n, $

(11)

$ \hat{x}_{j, \text { reg }}^{n+1}=\frac{1}{2 w_j^n} \sum\limits_{k \in N_j} w_{j, k}^n\left(x_j^n+x_k^n\right), $

(12)

$ w_{j, k}^n=\frac{\omega_{j, k}}{\sqrt{\left(x_j^n-x_k^n\right)^2+\delta^2}} . $

(13)

其中：w_jⁿ表示根据第n轮迭代图像xⁿ计算得到的针对像素j的权重因子；$\hat x_{j, {\rm{reg}}}^{n + {\rm{1}}}$表示一个中间变量，对应当前图像估计值的平滑结果；w_{j, k}ⁿ表示平滑操作时的权重因子。

综上可得目标函数ψ(x)的替代函数，

$ Q_\psi\left(\boldsymbol{x} ; \boldsymbol{x}^n\right)=Q_L\left(\boldsymbol{x} ; \boldsymbol{x}^n\right)-\beta Q_U\left(\boldsymbol{x} ; \boldsymbol{x}^n\right) . $

(14)

可证Q_ψ(x; xⁿ)满足文[9]描述的替代函数的2个基本条件，即

$ Q_\psi\left(\boldsymbol{x} ; \boldsymbol{x}^n\right)-Q_\psi\left(\boldsymbol{x}^n ; \boldsymbol{x}^n\right) \leqslant \psi(\boldsymbol{x})-\psi\left(\boldsymbol{x}^n\right), $

(15)

$ \nabla Q_\psi\left(\boldsymbol{x}^n ; \boldsymbol{x}^n\right)=\nabla \psi\left(\boldsymbol{x}^n\right). $

(16)

此后，ψ(x)的每一轮优化都可以转变为求当前Q_ψ(x; xⁿ)的极值。最终的迭代公式可表示为：

$ x_j^{n+1}=\frac{\sqrt{\left(1-\beta_j^n \hat{x}_{j, \mathrm{reg}}^{n+1}\right)^2+4 \beta_j^n \hat{x}_{j, \mathrm{EM}}^{n+1}}-\left(1-\beta_j^n \hat{x}_{j, \mathrm{reg}}^{n+1}\right)}{2 \beta_j^n}, $

(17)

$ \beta_j^n=\frac{\beta w_j^n}{p_j} . $

(18)

其中：β_jⁿ表示像素j在灵敏度因子p_j、正则化参数β和权重因子w_jⁿ下的综合作用。由于参数β、δ的选择与实际的像素取值范围有关，为避免受到影响，因此实际处理过程中会预先将待处理数据按照最大值进行归一化，同时将此结果作为迭代重建的初始值，待重建完成后再等比例恢复。

2.2 相关工作对比

一种传统的PSF解卷积算法是基于采集数据的Poisson分布模型并结合全变分先验构建目标函数，具体可表示为：

$ J_{\mathrm{reg}}(\boldsymbol{x})=\int|\nabla \boldsymbol{x}(\boldsymbol{r})| \mathrm{d} \boldsymbol{r}, $

(19)

$ \frac{\partial}{\partial \boldsymbol{x}} J_{\mathrm{reg}}^n=-\operatorname{div}\left(\frac{\nabla \boldsymbol{x}^n(\boldsymbol{r})}{\left|\nabla \boldsymbol{x}^n(\boldsymbol{r})\right|}\right). $

(20)

其中：J_reg(x)表示针对真实图像x的全变分先验；r对应空间位置；$\frac{\partial}{\partial \boldsymbol{x}} J_{\mathrm{reg}}^n$表示全变分项对当前估计图像xⁿ的偏导数，其在像素j对应位置处的值记为$\left[\frac{\partial}{\partial \boldsymbol{x}} J_{\mathrm{reg}}\right]_j^n$。

采用一步延迟法(one-step-late, OSL)解决目标函数难以优化的问题^[6]。在具体考虑散射数据的影响后^[5]，可将最终迭代公式表示为：

$ x_j^{n+1}=\frac{x_j^n}{1+\beta\left[\frac{\partial}{\partial \boldsymbol{x}} J_{\mathrm{reg}}\right]_j^n}\left[\boldsymbol{h}^{\mathrm{T}} \otimes \bf{b}^{\prime}\right]_j, $

(21)

$ {\bf{b}}_i^{\prime}=\frac{\left[y_i-s_i\right]_{+}}{[\boldsymbol{h} \otimes \boldsymbol{x}]_i} . $

(22)

其中：bp′表示反投影因子，bp′_i表示bp′的第i个元素，是散射校正后真实投影数据和当前模型计算结果的比值；[*]₊表示取括号内数值的非负部分。

对比式(6)和(21)可以发现，该算法与传统EM算法相比，主要有3种不同做法：1) 将不同像素位置处的灵敏度因子p_j全部视为1；2) 引入全变分先验的导数对p_j进行修正，以体现先验函数对重建的影响；3) 基于扣除散射数据后的投影进行重建。

由于实际考察图像的尺寸有限，导致不同像素位置的灵敏度因子p_j存在一定差异。第1种做法会导致图像边缘处出现一定差异；第2种做法导致选择参数β时要十分谨慎，并且无法保证迭代过程的非负性和单调收敛性。直接扣除散射数据的投影，理论上不再服从Poisson分布；而且由于噪声影响，直接扣除散射数据的做法可能导致某些像素位置出现负值，式(22)中的[*]₊操作会强制将这些像素值置为0，这些都会对最终结果造成一定影响。

2.3 散射校正

对于散射数据可基于多能窗法^[11-12]进行估计，

$ \boldsymbol{s}=\left(\frac{C_{\text {lower }}}{W_{\text {lower }}}+\frac{C_{\text {upper }}}{W_{\text {upper }}}\right) \times \frac{W_{\text {main }}}{2}. $

(23)

其中：C_lower、C_upper分别表示低能窗、高能窗计数，W_lower、W_upper、W_main分别表示低能窗、高能窗和主能窗的宽度。

对于本研究的SPECT系统，主能窗设置为328~400 keV, 相邻的低能窗和高能窗分别设置为304~ 328 keV和400~424 keV。

2.4 PSF获取

获取准确的PSF是伪影消除算法的关键。本文采用Monte Carlo模拟方法，基于GATE v8.0平台^[13]对SPECT探测器和高能通用准直器进行模拟，并采集了碘-131点源的投影数据。考虑到SPECT平片成像时甲状腺到准直器前表面的典型距离，模拟时点源与准直器前表面的典型距离设置为10 cm。

对于模拟采集的数据，提取主能窗中未经散射的数据用于生成PSF。为了降低PSF图像的噪声，采用未经散射的数据约976×10⁴计数，获得的碘-131点源PSF图像如图 3所示，PSF图像大小为255×255，像素为2.2 mm×2.2 mm。

此外，应用时还分别截取了大小为127×127、63×63、31×31的PSF，用于评估不同大小的PSF对于伪影消除效果的影响。这3种大小的PSF用到的事件数(指用于生成PSF的数据)分别占总计数的96.9%、84.0%和69.8%，占比越小表明舍弃的有效信息越多。

3 模拟数据评估 3.1 Monte Carlo模拟

为了评估算法性能并优化参数，本文基于GATE v8.0平台模拟了美国国家电气制造商协会(National Electrical Manufacturers Association，NEMA)躯干模型，如图 4所示, 在直径分别为13、17和22 mm的3个球体内充满碘-131放射源。图 5展示了模拟采集的数据能谱和各能窗投影图。主能窗、高能窗和低能窗的光子计数分别为22.8×10⁴、2.1×10⁴和4.3×10⁴。

3.2 收敛性验证

为了验证算法收敛性，分别选用4组不同参数进行重建，并计算相对函数值增量和相对图像差异随迭代次数的变化：

$ \text { 相对函数值增量 }=\frac{\left|\psi\left(\boldsymbol{x}^n\right)-\psi\left(\boldsymbol{x}^0\right)\right|}{\left|\psi\left(\boldsymbol{x}^*\right)-\psi\left(\boldsymbol{x}^0\right)\right|} \text {, } $

(24)

$ \text { 相对图像差异 }=\frac{\left|\boldsymbol{x}^{n+1}-\boldsymbol{x}^n\right|}{\left|\boldsymbol{x}^{n+1}\right|} \text {. } $

(25)

其中x⁰和x^*分别表示初始图像和最终收敛图像。本次工作中用500次迭代结果作为x^*进行后续评估。图 6展示了基于不同正则化参数重建的算法收敛性验证(前100次迭代结果)。由图可知，对于4组不同的重建参数，随着迭代次数增加，相对函数值增量快速增加并趋于1，而相对图像差异快速减小并趋于0，从而证实算法的单调收敛性。50次迭代后，4组参数对应的相对函数值增量均大于0.980，而相对图像差异均小于0.003。为了节省计算时间，本文后续工作中均采用50次迭代。

3.3 结果分析和参数优化

图 7展示了不同大小PSF的重建结果。由图可知，直接扣除散射数据可一定程度降低放射状伪影，由于准直器穿透效应，热球附近区域仍有明显计数，而使用PSF解卷积可以明显改善放射状伪影，降低热球附近区域计数。但是PSF过小时，仍会存在一定的放射状伪影，且显著放大背景区域噪声。PSF大小为127×127和255×255时的结果相对更好，背景噪声较低，无明显放射状伪影。

为了对结果进行定量评估，在重建图像的放射状条纹、条纹间隙和背景区域分别勾画了ROI，如图 8所示。通过计算大小为160 mm×50 mm的矩形ROI(蓝色)的标准差，对背景噪声进行评估。此外，为了定量描述放射状伪影的改善情况，在放射状条纹区域放置4个直径为11 mm的圆形ROI(黄色)，在条纹间隙放置3个直径为11 mm的圆形ROI(绿色)，并根据式(26)计算归一化伪影计数，具体可表示为

$ \text { 归一化伪影计数 }=\frac{(\text { 条纹区均值 }- \text { 条纹间隙均值 })_{\text {校正后 }}}{(\text { 条纹区均值 } - \text { 条纹间隙均值 })_{\text {校正前 }}}. $

(26)

校正图像定量分析结果如表 2所示。由表 2可知，随着PSF增大，归一化伪影计数逐渐降低。但是相对于仅做散射校正的结果，大小为31×31的PSF处理结果的归一化伪影计数更大。对比图 7中编号2和3可知，基于大小为31×31的PSF解卷积虽然可以压低条纹间隙计数，但是对于条纹区域计数的压低效果并不理想，导致归一化伪影计数增大。另外，虽然更大的PSF可以降低背景噪声，但是相比于原始数据，过小的PSF反而会导致背景噪声更大。综合来看，PSF大小为127×127和255×255时效果十分接近, 明显优于另外2组选择。因此下一步将针对大小为127×127的PSF引入先验函数进行进一步优化。

基于大小为127×127的PSF，图 9展示了不同β的重建图像，编号1的2条线表示剖线位置，用于后续分析。可以发现，随着β增大，背景噪声得到一定程度的压制，但是β过大时会导致过平滑。图 10展示了β为2.5×10^-3、1.0×10^-2、4.0×10^-2共3组重建图像剖线图，剖线位置如图 9中编号1所示，图 10a对应水平剖线，图 10b对应竖直剖线。由图 10a可知，随着β增大，相邻的充满放射源的玻璃管中间背景区域的像素值会明显增大，导致对比度降低。由图 10b可知，随着β增大，条纹伪影区域的值会逐步增大。这是由于式(3)中的先验函数在降低背景噪声的同时还要保持边界，因此放射状伪影区域的条纹特征也会被当作边界进行一定程度的强化。图 11给出了不同β(0, 1.25×10^-3，2.5×10^-3，5.0×10^-3，1.0×10^-2，2.0×10^-2，4.0×10^-2)重建的定量分析结果。前期随着β增大，归一化伪影计数并无明显变化，背景噪声明显减小。但是随着β进一步增大，归一化伪影计数也会随之增大，这与图 10展示的条纹伪影被强化的现象一致。

综上所述，最终选择PSF为127×127、β为1.0×10^-2用于临床评估。此时基于Monte Carlo模拟数据的归一化伪影计数为0.04，背景噪声为0.96。

3.4 传统算法与本文算法对比

对于2.2节所述传统算法和本文算法的差异，基于Monte Carlo模拟数据对2种算法的重建结果进行对比评估。选择同样的初始值、PSF(大小为127×127)和迭代次数(50次)，考察不同正则化强度、不同算法的重建结果，不同解卷积算法效果对比如图 12所示，β相同时传统算法背景噪声相对更大，且传统算法对β值较为敏感，β显著增大时重建结果会明显出现问题。相比之下，本文算法综合表现更好，β增大后图像变得更加平滑，符合预期，且没有出现传统算法结果中的点阵伪影。

4 临床数据评估

对于临床采集的点源数据和人体数据，分别采用本文算法进行处理。其中点源数据基于连续床位运动采集。表 3展示了临床数据各能窗计数的情况。图 13展示了临床数据的处理结果。采集的真实投影数据大小为1 024×256，图 13仅展示关注部分，大小为256×256。不论是点源数据还是人体数据，引入解卷积算法后放射状伪影均得到了有效压制。对于点源数据，图像分辨率得到了明显改善，并且由于计数较多，点源数据解卷积后并未明显导致背景噪声放大。相对而言，由于人体数据计数低、背景噪声大，单纯的解卷积算法会导致背景噪声显著增大，而引入先验函数可使背景噪声变小，在图像上表现为背景区域更加“干净”。

图 14和表 4分别展示了临床数据ROI勾画示意图和临床数据定量分析结果。条纹伪影和条纹间隙区域ROI直径为11 mm, 背景区域ROI大小为160 mm×50 mm。由表 4可知，相对于仅做散射校正，解卷积校正结果的归一化伪影计数明显更小。对于低计数的人体数据，先验正则化可以显著改善解卷积导致的高背景噪声，与图 13观测结果一致，代价是略微提升归一化伪影计数。此外，点源数据本身放射状伪影区域较小，ROI区域条纹特征相对已不是特别明显，因此引入先验函数后反而会被进一步平滑，导致归一化伪影计数进一步降低。

5 结论

针对碘-131 SPECT平片成像中的放射状伪影问题，本文提出了一种PSF解卷积结合先验正则化的方法，构建目标函数并推导了迭代算法。在Monte Carlo模拟实验中，对算法的单调收敛性进行了验证，研究了不同大小PSF、不同正则化参数β对重建结果的影响，并验证了本文算法相对于传统算法的优势。结果表明：虽然PSF可以有效降低放射状伪影，但也会引入背景噪声，PSF越小背景噪声越明显。而先验约束在降低背景噪声的同时也会一定程度保持放射状伪影特征。综合考虑选择PSF大小为127×127、正则化参数β为1.0×10^-2时，Monte Carlo模拟数据处理结果的归一化伪影计数仅为0.04。临床数据处理结果也证明了算法的有效性，点源数据和人体数据处理后的归一化伪影计数分别为0.35、0.28。

本研究基于Monte Carlo模拟方法获取临床典型距离(10 cm)的PSF，而临床采集时实际位置会有一定偏差，这也可能是临床数据结果相对模拟数据较差的原因。此外，先验函数中的超参数δ用于区分边界和噪声，像素差异远大于该参数时，会被认为是图像边界进行维持，而像素差异远小于该参数时，会被当作背景噪声进行平滑。2组实验数据测试中，由于伪影特征强度不同，导致引入先验函数对伪影区域计数的影响也存在一定差异。下一步工作一方面要考察不同位置的PSF对算法结果的影响，另一方面也要基于更多临床数据对算法进行进一步的参数优化和有效性验证。