基于SPECT探头的高灵敏编码孔径γ成像系统

刘鑫垚魏清阳许天鹏张朝晖

1（北京科技大学自动化学院北京市工业波谱成像工程技术研究中心 北京 100083）

2（国家核安保技术中心 北京 102445）

核辐射技术在工业、农业、医疗卫生等领域应用广泛，在推动了我国各行各业发展的同时，确保核应用安全也是十分重要的。此外，当前国际形势错综复杂，各国面临着核材料走私和核恐怖主义袭击的风险；放射源丢失和被窃事故也时有发生，因此远距离、大范围、快速实时检测放射源对国土安全具有重要意义［1］。γ相机可以远距离提取放射性物质位置、强度以及类别等信息，实现对放射性物质状态实时监控，极大地减少工作人员在工作过程中的辐射剂量，因而被广泛应用，其中编码孔径γ相机具有很大的潜力。但是，目前编码孔径γ相机趋于便携化［2−4］，有效探测面积小，且编码板准直器也限制了成像系统的视场，因此探测效率相对较低。如果对大范围内的放射源进行精确定位，需要借助多个设备的联合测量，不仅需增加设备安装与空间使用的成本，而且在设备后续维护过程中也会加大工作人员负荷量。

近些年，为了实现大范围内快速监测放射源，国际上主要是通过增大探测器面积实现高灵敏的编码孔径γ相机研制。2014年，美国为了打击海上非法核材料走私，研制了船载γ射线成像系统［5］，有效探测面积为0.74 m2，对于远距离放射源的成像时间约为1~120 s；2016年，美国研制一款车载核辐射检测系统［6］，有效探测面积为1.03 m2，角分辨率约为14°；2020年，国内高能物理研究所基于自主研发的大面积探测器研制了高灵敏度的γ相机［7］，有效探测面积为330 mm×330 mm，3.7×107Bq137Cs点源距离探测器为10 m处能够1 s有效探测，角度分辨率为5.7°。目前已经广泛应用的编码孔径γ相机RLS-1000［8］，有效视场为30°，探测器有效探测面积为46.4 mm×46.4 mm，角分辨率为1°~2°，其探测效率有限。为了有效提高γ成像系统的探测效率，采用国产的滨松BHP6601型临床单光子发射计算机断层成像（Single Photon Emission Computed Tomography，

SPECT）探 头［9］作 为 探 测 器，其 固 有 分 辨 率 为3.55 mm，有效探测面积为510 mm×390 mm。因此，本文将基于此探测器设计仿真一种具有高分辨率的高灵敏编码γ成像系统，其探测效率约相当于RLS-1000的100倍［8］。通过MATLAB解析模拟不同的放射源事件，对该γ成像系统的各项性能进行研究和评估，分析该系统的可行性，并挖掘该系统的潜能。

1 解析计算模型

本文将基于滨松BH6601型临床SPECT系统的探测器设计并仿真一种大面积高灵敏编码板γ相机，主要通过MATLAB对该系统的物理探测过程进行解析构造，模拟流程如图1所示。首先，基于探测器设计整个系统各项参数，模拟整个系统。其次，利用位于不同视场角度的放射源发射的γ射线经过编码板到达探测器的几何结构，从而得到系统传输矩阵A；然后解析构造不同的放射源事件，利用放射源与该系统的几何关系计算落入探测器上光子计数，利用系统传输矩阵的累积概率分布得到投影数据。最后采用最大似然估计期望最大化（Maximum Likelihood Expectation Maximization，MLEM）算法对投影进行重建。

图1 编码孔径γ成像系统物理探测过程的模拟流程图Fig.1 Flow chart of physical detection process simulation of coded aperture γ imaging system

1.1 系统设计

滨松BH6601型临床SPECT系统的探测器采用一块连续NaI闪烁晶体与光电倍增管阵列构成的探测器，其有效探测面积为510 mm×390 mm，厚度为9.5 mm，固有空间分辨率为3.55 mm。该探测器对137Cs的检测效率为23%。本文使用的探测器具有较多优点：1）固有分辨率较高，可以有效提高系统的分辨率；2）有效探测面积很大，具有高灵敏度的特性；3）成本低，各个性能参数优越，性价比较高。因此，本文将采用滨松BH6601型临床SPECT系统的探测器开展大面积高灵敏编码孔径γ相机的设计工作。

由于放射源发射γ光子的过程是各向同性的，探测晶体每个像素位置对接收远距离放射源发射的γ光子的概率基本相同，因此，无法直接用探测晶体显示放射源图像。准直器通过一定的开孔规则改变光子传播轨迹，从而达到空间定位选择的作用，可以有效提高成像距离较远的效果。单孔准直器虽然有较好的分辨率，但是开孔率低，进光量少，探测效率相对较低，因此，多针孔准直器被研发使用；但是多针孔设计困难，且很难平衡探测效率和空间分辨率等之间的需求矛盾，因此编码孔径成像技术在多针孔准直器的基础上，提供了一种全新的设计理论和加工制作方法。编码板准直器采用修正冗余阵列（Modified Uniformly Redundant Arrays，MURA）的编码模式，因为该模式不仅具有50%的高开孔率，而且在有限的几何大小内具有理想的相关性，还具有反对称结构优势，因此，具有较大的光通量，可以提高射线收集效率和灵敏度，这更能满足低剂量污染的成像需求［10］。为了提高编码板成像系统的视野范围，采用编码板嵌套模式。主要是先对编码板的编码模式做中心化处理，得到编码板的中心基本单元。根据MURA编码图案设计规则［11］，采用MURA（23×23）作为编码板的中心基本单元，利用MURA阵列周期相关性，使4个MURA（23×23）循环周期排列，从而形成编码板嵌套模式MURA（45×45），如图2所示。

图2 编码孔径的编码模式：MURA(45×45)Fig.2 Pattern of the coded mask:MURA(45×45)

由于方孔具有透光率高、加工方便、重建图像受噪声影响较小等优点，因此，编码板的孔径选择方形孔。按照编码板与探测器尺寸为2:1计算，MURA（45×45）探测面积为1 020 mm×780 mm。编码板材料将采用钨合金（密度18.75 g∙cm−3），其吸收射线能力强且成本低。综合考虑编码板对γ射线的阻挡本领、设备重量，编码板的厚度设计为10 mm。针对不同能量的放射源，质量衰减系数不同。针对于放射能量为0.661 7 MeV的137Cs放射源，其质量衰减系数为0.094 6 cm2∙g−1。根据质量衰减系数和编码板材料密度，可知该编码板的线衰减系数为1.77 cm−1。

编码孔径γ相机系统的视场主要依赖于编码板与探测器之间的距离，其中，在无伪影视场的范围内的放射源可以将MURA（23×23）可以完全投影在探测器上，其几何关系如图3所示。

图3 成像系统几何设置示意图Fig.3 Illustration of the geometry setup of the imaging system

由于探测器为矩形，成像系统视场设计目标：在X方向上，非伪影视场不低于60°，在Y方向上，非伪影视场不低于48°。主要根据X方向的设计目标设计该系统的焦距，如式（1）：式中：dm是编码板的中心基本单元在X方向的长度；b是成像系统的焦距（即编码板到探测器的距离）；θ是成像系统视野角度。在X方向上，结合dm=510 mm和θ=30°，由式（1）可得焦距解析估计值为441 mm。

根据以上的设计要求和方法，利用MATLAB模拟系统时，可供参考的设计参数如表1所示。

表1 成像系统的参数Table 1 Parameters of imaging system

1.2 模拟放射源投影数据

不同位置且不同活度的放射源通过编码板后落在探测器上的光子计数是不同的。放射源以自身为圆心向四周发射光子，形成一个球体空间，其中只有放射源与成像系统组成的立体角范围内的光子才有可能落在探测器上。因放射源距离探测器较远，在有效视场内，放射源与成像系统的垂直距离相同时，即便放射源相对成像系统中心偏离不同角度时，探测器检测到的光子计数几乎相同，因此我们将成像系统检测到的光子计数计算过程简单化，成像系统探测到的光子计数方法简化二维模型如图4所示。

图4 光子计数计算的二维模型Fig.4 Two-dimensional model of photon counting estimation for detection

由于编码板的开孔率为50%，因此，放射源与成像系统组成的立体角度范围内的光子只有50%通过编码板，落入成像系统。除此之外，该探测器对137Cs的有效检测效率为23%，因此探测器实际能够检测到光子数NS可由式（2）估算：

式中：Sd是编码板的中心基本单元的面积；R是以放射源为中心向四周发射光子，以编码板接收到探测光子为界限而形成的球体空间的半径，具体如图4所示；H是放射源的活度，Bq；η是NaI晶体对于137Cs的探测效率，为23%；T是成像系统采集光子的时间。

不同的编码板都对应不同的系统传输矩阵［12］。当放射源距离探测器较远时，将放射点源所在平行于编码板的平面，即视野平面，按照水平和垂直两个方向上，分割为一个P×P的矩阵，其中P为61。探测器平面分离为143×110的矩阵，则本文系统传输矩阵由式（3）表示为：

式中：μ是线衰减系数，主要由编码板的质量衰减系数和密度决定；lij是视野编号j处的放射源到达探测器编号i处所经过的编码板厚度，该厚度是由Siddon算法［13］解析计算。本文在编码板的厚度中选择10个点，从而更加精确判断放射源辐射的γ射线被吸收的程度，即更精确地计算了各个体素辐射的路径；ϕ和ψ是视野平面编号j处的射线与探测器平面X和Y轴的夹角。

第j个位置的放射源基于系统传输矩阵A的累积概率分布PC（i），如式（4）所示：

式中：Asj为第j个像素位置处的放射源对探测器第s个元素位置贡献粒子数的概率；N为探测器像素总个数。

本研究主要通过解析模拟特定事件的点源在探测器上形成的投影数据P。P的过程为：首先根据式（2），对不同位置和不同活度的放射源经过编码板后，计算探测器能够检测得到的光子数量为NS。其次基于累计概率分布（PC）执行均匀采样，该过程主要是通过选取随机函数中小于累计概率分布函数的索引值，将索引值中的最小值作为光子所在探测器的位置，这个过程循环NS次，即NS个光子根据概率分布在探测器上投影P；最后，将随机生成的均匀分布的伪随机整数作为探测器的索引值，在索引值对应的探测器位置处添加一个光子作为背景噪声，该过程循环1 000次，即每秒随机添加1 000个光子计数作为探测器的背景噪声，近似模拟真实的SPECT系统的探测器。

1.3 重建算法

本文主要通过MLEM迭代算法利用投影数据PC重建放射点源的位置。由于该模型建立在两个泊松统计过程：放射衰变点过程和探测器接收光子的过程，因此该算法能够有效抑制泊松噪声。在放射源平面的X方向上，无伪影视场角度约为60°，步长为0.979 2°；在Y方向上，无伪影视场角度约为47°，步长为0.759 8°。探测器平面分成Q个像素，即i=1，2，3，…，Q。MLEM算法迭代过程如下：

式中：i是探测器像素编号；j是视场位置编号；pi为探测器像素编号i上的光子计数；f kj是重建视场位置编号j处的第k次迭代结果；Aij是探测器平面的第i个像素对源平面的第j个像素的响应，相当于放射源所在的平面第j个像素对探测器第i个元素贡献粒子的概率，A可描述为系统传输矩阵。

由式（5），投影估计值∑jAij f k j和实际投影值pi相比较，将其比值通过传输矩阵进行相关运算对第k次图像重建的迭代结果进行修正，再利用系统传输矩阵A反投影到第k+1次的放射源强度分布图像，使它逼近放射源图像。

1.4 图像评价参数

对比度噪声比（Contrast-to-noise Ratio，CNR）用于评价图像的对比度，CNR更好意味着更容易找到源的位置，可以定量评价重建点源图像的质量［14−16］。CNR的定义如下：

式中：Rmax为重建图像的最大值；Bm为重建图像的背景噪声，即去除放射源附近3×3像素后每个像素的平均值；Bij为除源区域外的每个像素的值；Bsum为除源区域外的像素总数。CNR越大表示恢复图像相对噪声对比度越大，重建图像的质量越好。

角度分辨率用于评价成像系统能够分辨两个放射源相距的最小视场角度。该系统的视场角度分辨率主要是由成像系统几何角度分辨率和探测器固有角度分辨率共同决定［17］。理论来讲，在X和Y方向上的空间角度分辨率相差不多，因为系统几何角度分辨率主要是受编码孔的单位孔径和系统焦距影响。根据系统设计参数，编码孔径是矩形，因此，成像系统的角度分辨率将从X和Y方向分别评估。角度分辨率的确定主要采用了两种评估方式：第一种是采用单点源的重建图像，根据单点源的重建信号在X和Y方向上的归一化强度进行高斯拟合，得到半高宽（Full Width at Half Maximum，FWHM）［18］，根据式（7）求出该系统的角度分辨率Δθ；第二种是通过模拟两个相距一定视场角度的放射源，观察其能够完全将两个点源分辨开时最小的视场角度距离，这个相差的角度距离即角度分辨率Δθ。

式中：FWHM为单点源信号在X或Y方向上像素值的拟合高斯曲线的半高宽；n是单点源在X或Y方向上的像素总数；NAFOV'是实际重建图像显示的无伪影视场。

2 性能参数测试与讨论

2.1 成像系统的灵敏度

灵敏度是编码孔径γ相机的重要参数，灵敏度主要评价了成像系统对于较远或较弱的放射源的探测能力。在无伪影视场范围内，该系统对不同位置的单点源进行重建，观察重建效果。本文模拟了3.7×107Bq137Cs点源距离大面积编码孔径γ相机分别为10 m、15 m、18 m的情况下，及其位于视野角度分别为（−29.64°，23.62°）、（0°，23.65°）、（0°，0°）、（28.66°，0°）、（29.64°，0°）、（−14.87°，−14.87°）处，共计18种不同位置的模拟情况，探测时间为1 s。如图5所示，每行分别是放射源与成像系统距离相同的情况下，不同视场角度的放射源重建图像；每列分别是放射源视场角度位置相同的情况下，放射源位于不同距离处的重建效果。其中每个图中圆圈代表放射源应该重建的位置，亮斑代表系统实际重建的放射源位置。图5中第2列表示放射源位于Y轴方向上的边缘位置时的重建图像，可知图像重建质量良好，无伪影；第4列和第5列表示放射源位于视场X轴方向上的边缘位置时的重建图像。在（28.66°，0°）处，单点源的重建图像质量较好，而在（29.64°，0°），单点源的重建图像出现了伪影，即关于Y轴的对称位置有一个容易发生误判的亮斑，这是因为在设计编码板的过程中，实际的编码板的中心基本单元大小的长度比探测器的有效面积的长度小，并没有与目标尺寸达到完全一致。因此在X方向上，该成像系统的实际有效全视场范围应该在−28.66°~28.66°，在Y方向上，实际有效视场范围应该在−23.65°~23.65°。

图5 单点源3.7×107 Bq137Cs与大面积编码孔径γ相机的距离分别为10 m、15 m、18 m，及其位于视野角度分别为(−29.64°,23.65°)、(0°,23.65°)、(0°,0°)、(28.66°,0°)、(29.64°,0°)、(−14.87°,−14.87°)的重建图像Fig.5 Reconstructed images of a point source at different distances and different FOV angles:the 1st,2nd,3rd,4th,5th and 6th columns are positions at(−29.64°,23.65°),(0°,23.65°),(0°,0°),(28.66°,0°),(29.64°,0°)and(−14.87°,−14.87°);the 1st,2nd and 3rd rows are the distance between the radiation source and the detector of 10 m,15 m and 18 m,respectively

表2是对应图5中每个重建图像的CNR。根据图5每列可知，随着放射源与大面积编码孔径γ相机的距离的增加，CNR下降，图像质量随之下降，但是，圆圈和亮斑基本吻合，CNR基本都大于10，说明该系统对单点源远距离的位置能够准确重建。当3.7×107Bq137Cs的单点源与该系统距离为18 m时，1 s的采集时间，该系统能够对无伪影视场内的点源位置准确成像。根据图5与表2可知，圆圈和亮斑基本吻合，在实际全视场范围内，能够对单点源的位置准确成像，因此大面积编码孔径γ相机远距离探测辐射源具有一定的可行性。其中当放射源位于中心位置时，CNR最低，随着单点源与编码板中心偏离的角度增加，图像的CNR增大，图像质量提高，但是放射源与成像系统距离相同时，不同角度的放射源重建图像的CNR相差不多，说明放射源的距离对成像的质量影响较大。

表2 不同位置的放射源重建图像的CNRTable 2 CNR of reconstructed images from different positions of radioactive sources

当1.11×106Bq137Cs低活度点源置于成像系统前4 m时，分别位于（−29.64°，23.65°）、（0°，23.65°）、（0°，0°）、（28.66°，0°）、（−14.87°，−14.87°）这5个视场角度处，采集时间仅需要2 s即可重建，如图6所示。其中圆圈表示放射源应该被重建的位置，利用方形对图像中强度最大的亮斑进行框选，即辅助定位，圆圈与方形基本吻合，且CNR基本大于10，因此该成像系统对低活度放射源在短时间内进行探测时能够准确重建放射源的位置。

图6 单点源1.11×106 Bq137Cs与大面积编码孔径γ相机的距离为4 m，及其位于视野角度分别为(−29.64°,23.65°)、(0°,23.65°)、(0°,0°)、(28.66°,0°)、(−14.87°,−14.87°)的重建图像Fig.6 Reconstruction images of 1.11×106 Bq137Cs at distance of 4 m from detector with different FOV angles:the 1st,2nd,3rd,4th and 5th columns are positions at(−29.64°,23.65°),(0°,23.65°),(0°,0°),(28.66°,0°)and(−14.87°,−14.87°),respectively

2.2 成像系统的空间角度分辨率

空间角度分辨率是成像系统性能的重要指标，代表了成像系统能够分辨两个放射点源相距的最小视场角度。本文采用两种方法评估角度分辨率。

第一种方法：采用了3.7×107Bq137Cs在10 m处，成像系统通过1 s的采集数据进行重建，为了获取更多的数据点，该方法将放射源所在的平面即重建的视野平面，分为221×221的阵列，重建图像如图7（a）所示。然后选取放射源分别在X和Y方向上的对称的数据点归一化后进行高斯拟合，如图7（b）和图7（c）所示。根据式（7），可知空间角度分辨率在X和Y方向上分别为0.86°和0.85°。

图7 (a)单点源3.7×107 Bq137Cs与大面积编码孔径γ相机的距离为10 m，视野角度为(0°,0°)的重建图像，单点源沿X(b)和Y(c)方向上归一化后的数据点的高斯拟合像素强度分布Fig.7(a)Reconstructed image of 3.7×107 Bq137Cs point source at distance of 10 m from detector with FOV angle(0°,0°),Gaussian fitting of normalized data points of pixel intensity distribution along the X direction(b)and the Y direction(c)of a single point source

第二种方法：该方法的放射源所在的平面即重建视野平面是61×61的阵列，采用两个3.7×107Bq137Cs点源，在X方向上间隔2.94°，置于编码孔径γ相机前10 m处，采集时间为1 s，该成像系统能够对5个不同的视场角度处的两个单点源进行清晰成像，如图8（a）所示，采用两个3.7×107Bq137Cs点源，在Y方向上间隔2.28°，置于编码孔径γ相机前10 m处，采集时间为1 s，在5个不同的视场角度处对双点源能够清晰成像，该系统能够基本区分这两个点源，如图8（b）所示。

图8 两个相同的3.7×107 Bq137Cs与大面积编码孔径γ相机的距离为10 m处的重建图像(a)X方向上，两个点源方位角差为2.94°的重建图像，左点源的视场角度为(−29.64°,23.65°)、(0°,23.65°)、(0°,0°)、(26.70°,0°)、(−14.87°,−14.87°)，(b)Y方向上，两个点源方位角差为2.28°的重建图像，上点源的视场角度为(−29.64°,23.65°)、(0°,23.65°)、(0°,0°)、(26.70°,0°)、(−14.87°,−14.87°)Fig.8 Reconstructed images of two 3.7×107 Bq137Cs sources at the distance of 10 m from detector with acquisition time of 1 s(a)The 1st,2nd,3rd,4th and 5th columns are left and upper point source position at(−29.64°,23.65°),(0°,23.65°),(0°,0°),(26.70°,0°)and(−14.87°,−14.87°),respectively,in the X-direction with azimuth difference of 2.94°,(b)The 1st,2nd,3rd,4th and 5th columns are upper point source position at(−29.64°,23.65°),(0°,23.65°),(0°,0°),(26.70°,0°)and(−14.87°,−14.87°),respectively,in the Ydirection with azimuth difference of 2.38°

2.3 讨论

基于SPECT探测器的高灵敏编码孔径γ相机具有较好的视场和角度分辨率。其无伪影的视场在X和Y方向上分别为57.32°和47.3°，根据几何结构，具有广阔的探测范围。根据双点源成像，该系统角度分辨率在X和Y方向上分别为2.94°和2.28°。该系统能够对18 m远处的3.7×107Bq的137Cs仅需采集1 s，CNR均大于10，能够有效图像重建，说明了该系统对于远距离的放射源的快速成像具有可行性。当1.11×106Bq的137Cs距离成像系统4 m时，能够准确地定位放射源的位置，采集时间仅需要2 s，可见该系统对于较弱的放射源可以实现快速定位。综上所述，该系统在灵敏度和空间分辨率上都有了新的突破，对于大范围内快速搜索丢失放射源具有很大的优势。在后续工作中，将通过Gate进行系统仿真以及放射源的物理模拟过程，进一步分析该系统的可行性和潜能。在此基础上进行一些优化：在图像重建的数据处理中，如何在保证系统性能指标的前提下，可以进行分区域的计算，缩小系统传输矩阵，可以有效提高工作效率，因此还需要进一步优化算法，提高算法执行速度，这样可以有效提高系统探测效率；在编码板的结构设计过程中还可以进行优化，使编码板的设计需要更加匹配探测器，从而在保证具有较高分辨率的条件下，提高系统灵敏度。

3 结语

创新点主要是基于滨松BH6601型临床SPECT系统的探测器进行系统的研制，该探测器不仅有较高的固有分辨率，而且有较大的有效探测面积，有效提高系统的空间角度分辨率和灵敏度；该探测器成熟且成本低，性价比较高，因此也大大缩短了该成像系统的开发周期。提出一种基于SPECT系统探测器的大面积高灵敏编码孔径γ相机的设计方法，并对该系统性能进行了测试和评价，结果表明：该系统实现大视野范围内探测，具有较高灵敏度的同时空间角度分辨率较好，弥补市场此类产品的空缺。

作者贡献声明刘鑫垚：负责设计详细方案、系统仿真、图像重建和数据分析工作，并完成论文的撰写；魏清阳：提出方案、论文修改、研究基金的获取；许天鹏：负责提供技术支持和指导、研究资金的获取；张朝晖：负责论文修改，工作进度的监督。