基于禁忌搜索的晶体位置谱分割与校准方法

邓贞宙，陈冠东，王 平，封子纪，韩春雷,2

(1.南昌大学 信息工程学院，江西 南昌 330000； 2.图尔库大学中心医院 芬兰国家PET中心，芬兰 图尔库 999018)

在核医学领域，正电子发射断层成像(PET)的关键是通过晶体阵列上的响应事件位置信息来确定湮灭事件响应线所在位置，响应事件位置信息的精确度影响着PET系统的能量、时间校正准确度。近年来，PET系统研发朝着高空间分辨率的方向发展，耦合的晶体阵列规模逐渐增大而晶体条尺寸不断减小，当晶体尺寸低至某个临界值时，固有空间分辨率便不会降低，于是位置数据获取及校正方法设计便成为了重要变量[1]。在现实情况下，由于硅光电倍增管(SiPM)耦合方式、探测器电路设计带来的噪声、位置谱生成程序的不足等因素使得位置谱呈现了大量非线性特征，甚至呈现枕形、桶形失真等问题，给后续的图像重建带来了误差甚至误判，因此需选择并设计一种位置谱分割方法来划分、归正晶体阵列中每个晶体条覆盖区的响应[2]。现有最普遍的分割方法为分水岭分割方法，其缺陷在于极易产生过分割，会误选原始位置谱中噪声信号带来的一些伪峰，且算法存在时间复杂度抛物线增长的问题[3-4]。

针对以上问题，本方法通过二维高斯模型将来自探测器的UDP数据包中的位置信息进行位置谱叠加，利用禁忌搜索(TS)算法的非局部搜索机制对位置谱进行全局搜索，将湮灭事件对应晶体上各响应区域的局部最大值记录到禁忌表并视为响应峰值点，同时将其附近禁忌邻域(以禁忌长度为边长的矩形，禁忌长度受折射率影响)内所有点视为禁忌，以此避免迭代搜索和过分割。随后将特赦准则作为程序的负反馈机制，将处于孤立状态的峰值点视为伪峰并解除其邻域的禁忌状态，继续迭代，以此过滤位置谱中的噪声，直至所有晶体条对应响应区域划分完毕，最后根据禁忌表和晶体条分布情况生成晶体位置查找表，实现湮灭事件位置校准。

1 数据获取系统架构

本方法的湮灭事件数据来源和实验载体为一套高性能PET数据获取系统，其硬件结构主要由4个PET基本探测器模组(BDM)、时钟系统、服务器等构成。作为PET成像系统中捕获γ射线事件信息的眼睛，单个BDM模块结构由4个硅酸钇镥(LYSO)闪烁晶体阵列和SiPM耦合构成的前端探头、模拟数字转换器(ADC)模块、稀疏量化电平(SQL)模块和时间数字转换器(TDC)电路组成。

PET数据获取系统模块结构如图1所示，18F-FDG药剂衰变后产生正电子，与负电子湮灭产生的γ射线在晶体阵列上某位置响应生成闪烁光子，SiPM将其转化为闪烁脉冲信号，随后ADC电路输出能量、位置信息，与SQL和TDC模块生成的时间信息送至FPGA电路，FPGA通过重心法和Anger Logic算法提取位置数据，并封装为单事件数据帧，以太网电路将其打包为UDP数据包，通过网线传输给本方法的位置谱分割校准模块进行处理[5-6]。

图1 PET数据获取系统模块结构Fig.1 PET data acquisition system module structure

2 位置谱生成程序

2.1 UDP数据帧提取模块

本文的数据帧封装方法为：通过FPGA软核Verilog代码将每16字节长度的16进制数据[7]封装为UDP数据包，包含了湮灭事件所有相关的位置、能量、时间数据和探测器相关信息。而本程序只提取第1、2、11和12字节的晶体编号、探测器编号和横纵坐标数据各响应事件的通道编号、探测器编号和横纵坐标数据。

为了能方便、快速地提取UDP数据，在各探测器通过交换机连接服务器构建UDP网络后，本程序将获取并绑定各IP地址与端口，定时定量地批量采集并读取UDP数据包，提取出数据帧中的以上信息，随后按照二维高斯模型生成位置谱。程序的可视化界面及网络编程通过Qt Creator框架开发。

2.2 位置谱生成模块

本质上，位置谱是γ光子入射在探测器晶体条上的位置信息累加所绘成的图谱。图2为三维曲面形式的13×13晶体位置谱，将位置谱视为一个三维直方图，其坐标系x、y轴对应晶体平面上位置点的坐标，z轴对应闪烁事件计数量(即为此位置响应事件的数量)。

图2 三维曲面形式的13×13晶体位置谱Fig.2 13×13 crystal position spectrum in form of three-dimensional curved surface

单字节最多容纳256个数值，故将位置谱的位置信息按256×256来划分。对n×n晶体阵列而言，理想的位置谱上会清晰显示n2个光斑。理论上，对于n×n晶体探测器生成的位置谱中的每个晶体条而言，累计记录的γ射线响应位置P均服从二维正态分布(又名二维高斯分布)，则：

Pij(x,y)=(2πσijξij)-1·

(1)

其中：i、j分别为位置谱对应晶体条的行和列；μij、λij分别为x、y方向的期望；σij、ξij分别为x、y方向的标准差，(x，y)上相关系数为0。

在晶体条阵列中，响应强度从晶体中心到边缘呈现逐渐减弱趋势，晶体条上响应的局部最大值与二维正态分布的峰值相互映射。位置谱上所有位置点的分布均呈现混合高斯分布特点，并受到晶体条上产生闪烁光子的权重因子W影响[8-9]，则：

(2)

图3为计数值累加示意图，位置谱的关键步骤是将数据帧列表中每帧数据按照封帧格式读取出事件位置信息，并将位置信息叠加到256×256矩阵的相应直方图上(位置谱根据探测器、通道编号单独处理)。

图3 计数值累加示意图Fig.3 Schematic of counting value accumulation

位置谱获取的程序流程图如图4所示，主要是将γ光子打在对应探测器各通道上的位置信息绘成图谱，并按照二维高斯模型累加求和，得到最终的位置谱。

图4 位置谱获取的程序流程图Fig.4 Flow chart of position spectrum acquisition algorithm

3 TS分割算法及晶体位置查找表

3.1 TS分割算法的关键要素及参数计算

为了划分位置谱中响应区域的离散点并滤除伪峰，设计了一种TS分割算法程序，本程序的核心思想为：将整个256×256位置谱视为矩阵，将所有点按计数量排序，通过禁忌表机制记录峰值并将其邻域坐标赋值为1，从而快速而不重复地迭代搜索，将等同于晶体条数目的禁忌邻域分割出来[10]，可构成晶体覆盖区上γ射线作用的响应位置分割结果图，算法包含以下5大要素。

1) 禁忌表(Tabulist)，数值为1的坐标即为禁止状态，其数据内容(禁忌点坐标与顺序)存储在Tabulist(晶体位置查找表的基础)变量中。Tabulist是一含有n2个标记点的n2行2列的矩阵，按照对应计数值从上到下降序排列。

2) 终止准则为晶体条的数量，如13×13晶体，其终止准则为输出169个禁忌邻域后循环终止。

3) 禁忌长度TL是指每个方块禁忌邻域的边长(即所含位置点个数)，其长度取决于晶体与SiPM间的光导玻璃&光学胶水等介质折射率Nd和晶体阵列总行(列)数n，其计算公式为：

(TL+5)(n-1)Nd=256

(3)

晶体探头所用光导介质为聚甲基丙烯酸甲脂(一般室温25 ℃条件下折射率为1.49)，但由于粘合晶体的光学胶水(有机硅凝胶，其折射率为1.41)降低了折射率，且BDM在采集数据时会发出大量热导致光导玻璃折射率下降[11]，故实际折射率约为1.41～1.49。对多种晶体阵列进行评估后，得出晶体到SiPM的实际折射率Nd=256/175≈1.462 9，符合温度及胶水带来的影响。图5为7种晶体阵列的禁忌长度变化曲线。在实际的PET探测器应用中，常用的4×4、6×6、10×10、12×12、13×13、16×16、20×20等7种晶体阵列的TL分别为53、30、14、11、10、7、4个坐标长度(256×256位置谱已是数据容量极限，位置点个数不存在小数，故取近似为整数)。

图5 7种晶体阵列的禁忌长度变化Fig.5 Tabu length change of seven kinds of crystal arrays

4) 禁忌邻域，对应湮灭事件响应位置，以局部响应峰值坐标(xi±0.5TL，yi±0.5TL)区域的矩形。图6为13×13晶体阵列的禁忌邻域，13×13晶体的禁忌邻域为11×11区域，含121个位置点坐标。

图6 13×13晶体阵列的禁忌邻域Fig.6 Tabu neighborhood of 13×13 crystal array

5) 特赦准则，检测响应标记谱Ff中所有响应标记点Ki(xi,yi)是否满足孤立状态，具体实现为(xi+0.5TL,yi)=(xi-0.5TL,yi)=0或(xi,yi+0.5TL)=(xi,yi-0.5TL)=0，若是则实行特赦，解除此点及其禁忌邻域的禁忌状态，继续进行禁忌搜索更新分割。图7为伪峰的特赦处理，峰值点延伸的水平/竖直线不与其他任何响应标记点相交。

图7 伪峰的特赦处理Fig.7 Amnesty treatment of false peak

3.2 TS分割算法的实现步骤

图8为TS分割算法流程图，具体内容如下。

图8 TS分割算法流程图Fig.8 TS segmentation algorithm flow chart

S1赋值：将位置谱F1的三维位置信息进行分割预处理并整合为256×256矩阵F2。

S2排序：取i=1、j=1，按数值从大到小排序将F2中的坐标转化为169行2列的索引矩阵Ix，同时把禁忌邻域谱Ff赋值为256×256的0矩阵，取矩阵Tabulist赋值为169行2列的0矩阵。

S3禁忌检查：为判断峰值是否已被禁忌，记录Ix索引中数值排第j行的坐标(xj,yj)，判断Ff矩阵内的(xj,yj)坐标数值是否为1，是则继续，否则j=j+1并重复S3步骤。

S4添加禁忌表：将F2中的(xj,yj)坐标视为峰值点Ki(即相应晶体条的γ射线响应区域局部最大值)，将坐标赋值到Tabulist中第j行，将F2中的Ki点用红色圆圈标记。

S5禁忌邻域：为划分出晶体条的γ射线响应区域，将矩阵Ff内(xj±0.5TL,yj±0.5TL)方形邻域内所有坐标赋值为1，视为1个禁忌邻域Ti。

S6更新结果图：输出此时含所有Ki的响应峰值图，如图9a所示；输出此时含所有Ti的禁忌邻域图(最终结果为响应位置分割图)，如图9b所示。

图9 正在标记响应点的响应峰值图(a)和禁忌邻域图(b)Fig.9 Respons peak map (a) and tabu area map (b) of response point being marked

S7终止准则：为确保所有晶体条映射区域的响应皆被分割，判断i

S8特赦：作为算法负反馈机制，排除噪声导致误选伪峰为响应峰值的情况，判断是否符合特赦准则，若符合返回步骤S5，否则结束。

随着寻峰迭代和禁忌操作，单次搜索的对象数会逐步减少，避免了搜索过程陷入循环或局部最优，寻峰操作也会越来越快，时间复杂度大幅减小；此外，由步骤S8可知，特赦准则滤除了出现冲突的伪峰情况，避免了过分割。

3.3 构造晶体位置查找表

为确定各晶体单元覆盖区域内γ射线作用位置与实际位置谱光斑位置点的映射关系，构建了晶体位置查找表以便更准确地确定湮灭事件所在响应线的位置。

本方法以晶体条分布情况为约束，将响应位置分割图的各禁忌邻域编号为1～n2，将Tabulist中所有点按坐标逐行、逐列排序填入其中，将各离散位置点归属到对应禁忌邻域，输出所有晶体的位置查找表CPLT，输出的结果如图10所示，从而根据Tabulist确定每个晶体条对应的区域编号并进行湮灭事件响应线位置校准。

图10 晶体响应映射图(a)及晶体查找表结果(b)Fig.10 Crystal response map (a) and crystal lookup table result (b)

4 实验结果

在具体实验中，放置装有活度为1.11×105Bq的18F-FDG注射器在距探测器模块1.5 cm处位置，使用探测器扫描10 s，每个探测器模块平均捕捉到约350 000个光子事件数。程序输出位置谱耗时0.31 s，输出位置谱结果如图11所示。

a——13×13位置谱；b——6×6位置谱图11 两种晶体位置谱结果Fig.11 Two crystal position spectrum results

为评估位置数据准确度，在溢出率为5%的条件下，将本方法获取的13×13晶体的位置谱按照水平、竖直方向进行了计数值的投影，所得结果如图12所示，通过计算获得位置谱的13个峰值与12个谷值的峰值比并求均值后可得：水平峰谷比为12.18，竖直峰谷比为9.45，高于现有的大多数方法生成的位置谱效果[12-13]，通过拟合方法计算投影的半高宽(FWHM)获得1.7 mm的最佳空间分辨率。

平均每执行70.43 ms完成1个BDM的晶体位置谱分割，分割结果如图13所示，使用的13×13和6×6晶体阵列分割出169个、36个晶体标记区域，其中每个晶体区域的响应光斑均能清晰显示，保证了所有晶体条覆盖区的γ射线响应位置得到映射，无响应区域遗漏，且无噪声引起伪峰，从结果上进一步证实了算法避免过分割的功能。

a——水平方向投影；b——竖直方向投影图12 13×13晶体位置谱的峰值及投影曲线Fig.12 Peak value and projection curve of 13×13 crystal position spectrum

a——13×13响应峰值图；b——13×13响应位置分割图；c——6×6响应峰值图；d——6×6响应位置分割图图13 位置谱分割结果Fig.13 Segmentation result of position spectrum

为从数据上进一步证实本方法的效率优势，将TS分割算法与分水岭分割算法[14-15]的时间复杂度进行了对比(N为阵列中晶体条数目)，结果列于表1。分水岭分割算法理论上运算次数为N(1+N)+p次(p对应后续待处理的伪位置点个数)，时间复杂度为O(N2)，运行时间约为1～3 s。本算法的运算次数为N(1+k)次(k为排序N以内被特赦准则滤除的伪峰数量)，即时间复杂度为O(N)，对应的13×13晶体阵列运算时间为0.070 43 s。

表1 两种方法的运行时间对比Table 1 Comparison of running time of two methods

5 结论

本文开发了一种以禁忌搜索算法为核心，结合UDP数据帧、二维高斯模型及光导折射率的位置信息分割校准方法，该方法获得的位置谱峰谷比优于大多数其他方法，分割效率高于分水岭算法，且避免了过分割现象，滤除了噪声引发的伪峰，最终实现了系统空间分辨率的提高。因此，本方法在自动分割强度分布不均或形变失真程度高的位置谱具有较大优势，具有校准响应事件位置准确性的负反馈机制，分割与校准可指导后续符合处理，从而准确地确定湮灭事件所在响应线的位置，同时可作为测试大量晶体阵列乃至PET系统的评估标准。本方法已在南昌大学PET实验室的智能PET影像系统中表现出良好的适用性。