鄧貞宙，陳冠東，王 平，封子紀(jì)，韓春雷,2

(1.南昌大學(xué) 信息工程學(xué)院，江西 南昌 330000； 2.圖爾庫大學(xué)中心醫(yī)院 芬蘭國家PET中心，芬蘭 圖爾庫 999018)

在核醫(yī)學(xué)領(lǐng)域，正電子發(fā)射斷層成像(PET)的關(guān)鍵是通過晶體陣列上的響應(yīng)事件位置信息來確定湮滅事件響應(yīng)線所在位置，響應(yīng)事件位置信息的精確度影響著PET系統(tǒng)的能量、時間校正準(zhǔn)確度。近年來，PET系統(tǒng)研發(fā)朝著高空間分辨率的方向發(fā)展，耦合的晶體陣列規(guī)模逐漸增大而晶體條尺寸不斷減小，當(dāng)晶體尺寸低至某個臨界值時，固有空間分辨率便不會降低，于是位置數(shù)據(jù)獲取及校正方法設(shè)計便成為了重要變量[1]。在現(xiàn)實情況下，由于硅光電倍增管(SiPM)耦合方式、探測器電路設(shè)計帶來的噪聲、位置譜生成程序的不足等因素使得位置譜呈現(xiàn)了大量非線性特征，甚至呈現(xiàn)枕形、桶形失真等問題，給后續(xù)的圖像重建帶來了誤差甚至誤判，因此需選擇并設(shè)計一種位置譜分割方法來劃分、歸正晶體陣列中每個晶體條覆蓋區(qū)的響應(yīng)[2]?，F(xiàn)有最普遍的分割方法為分水嶺分割方法，其缺陷在于極易產(chǎn)生過分割，會誤選原始位置譜中噪聲信號帶來的一些偽峰，且算法存在時間復(fù)雜度拋物線增長的問題[3-4]。

針對以上問題，本方法通過二維高斯模型將來自探測器的UDP數(shù)據(jù)包中的位置信息進(jìn)行位置譜疊加，利用禁忌搜索(TS)算法的非局部搜索機(jī)制對位置譜進(jìn)行全局搜索，將湮滅事件對應(yīng)晶體上各響應(yīng)區(qū)域的局部最大值記錄到禁忌表并視為響應(yīng)峰值點，同時將其附近禁忌鄰域(以禁忌長度為邊長的矩形，禁忌長度受折射率影響)內(nèi)所有點視為禁忌，以此避免迭代搜索和過分割。隨后將特赦準(zhǔn)則作為程序的負(fù)反饋機(jī)制，將處于孤立狀態(tài)的峰值點視為偽峰并解除其鄰域的禁忌狀態(tài)，繼續(xù)迭代，以此過濾位置譜中的噪聲，直至所有晶體條對應(yīng)響應(yīng)區(qū)域劃分完畢，最后根據(jù)禁忌表和晶體條分布情況生成晶體位置查找表，實現(xiàn)湮滅事件位置校準(zhǔn)。

1 數(shù)據(jù)獲取系統(tǒng)架構(gòu)

本方法的湮滅事件數(shù)據(jù)來源和實驗載體為一套高性能PET數(shù)據(jù)獲取系統(tǒng)，其硬件結(jié)構(gòu)主要由4個PET基本探測器模組(BDM)、時鐘系統(tǒng)、服務(wù)器等構(gòu)成。作為PET成像系統(tǒng)中捕獲γ射線事件信息的眼睛，單個BDM模塊結(jié)構(gòu)由4個硅酸釔镥(LYSO)閃爍晶體陣列和SiPM耦合構(gòu)成的前端探頭、模擬數(shù)字轉(zhuǎn)換器(ADC)模塊、稀疏量化電平(SQL)模塊和時間數(shù)字轉(zhuǎn)換器(TDC)電路組成。

PET數(shù)據(jù)獲取系統(tǒng)模塊結(jié)構(gòu)如圖1所示，18F-FDG藥劑衰變后產(chǎn)生正電子，與負(fù)電子湮滅產(chǎn)生的γ射線在晶體陣列上某位置響應(yīng)生成閃爍光子，SiPM將其轉(zhuǎn)化為閃爍脈沖信號，隨后ADC電路輸出能量、位置信息，與SQL和TDC模塊生成的時間信息送至FPGA電路，F(xiàn)PGA通過重心法和Anger Logic算法提取位置數(shù)據(jù)，并封裝為單事件數(shù)據(jù)幀，以太網(wǎng)電路將其打包為UDP數(shù)據(jù)包，通過網(wǎng)線傳輸給本方法的位置譜分割校準(zhǔn)模塊進(jìn)行處理[5-6]。

圖1 PET數(shù)據(jù)獲取系統(tǒng)模塊結(jié)構(gòu)Fig.1 PET data acquisition system module structure

2 位置譜生成程序

2.1 UDP數(shù)據(jù)幀提取模塊

本文的數(shù)據(jù)幀封裝方法為：通過FPGA軟核Verilog代碼將每16字節(jié)長度的16進(jìn)制數(shù)據(jù)[7]封裝為UDP數(shù)據(jù)包，包含了湮滅事件所有相關(guān)的位置、能量、時間數(shù)據(jù)和探測器相關(guān)信息。而本程序只提取第1、2、11和12字節(jié)的晶體編號、探測器編號和橫縱坐標(biāo)數(shù)據(jù)各響應(yīng)事件的通道編號、探測器編號和橫縱坐標(biāo)數(shù)據(jù)。

為了能方便、快速地提取UDP數(shù)據(jù)，在各探測器通過交換機(jī)連接服務(wù)器構(gòu)建UDP網(wǎng)絡(luò)后，本程序?qū)@取并綁定各IP地址與端口，定時定量地批量采集并讀取UDP數(shù)據(jù)包，提取出數(shù)據(jù)幀中的以上信息，隨后按照二維高斯模型生成位置譜。程序的可視化界面及網(wǎng)絡(luò)編程通過Qt Creator框架開發(fā)。

2.2 位置譜生成模塊

本質(zhì)上，位置譜是γ光子入射在探測器晶體條上的位置信息累加所繪成的圖譜。圖2為三維曲面形式的13×13晶體位置譜，將位置譜視為一個三維直方圖，其坐標(biāo)系x、y軸對應(yīng)晶體平面上位置點的坐標(biāo)，z軸對應(yīng)閃爍事件計數(shù)量(即為此位置響應(yīng)事件的數(shù)量)。

圖2 三維曲面形式的13×13晶體位置譜Fig.2 13×13 crystal position spectrum in form of three-dimensional curved surface

單字節(jié)最多容納256個數(shù)值，故將位置譜的位置信息按256×256來劃分。對n×n晶體陣列而言，理想的位置譜上會清晰顯示n2個光斑。理論上，對于n×n晶體探測器生成的位置譜中的每個晶體條而言，累計記錄的γ射線響應(yīng)位置P均服從二維正態(tài)分布(又名二維高斯分布)，則：

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

(1)

其中：i、j分別為位置譜對應(yīng)晶體條的行和列；μij、λij分別為x、y方向的期望；σij、ξij分別為x、y方向的標(biāo)準(zhǔn)差，(x，y)上相關(guān)系數(shù)為0。

在晶體條陣列中，響應(yīng)強(qiáng)度從晶體中心到邊緣呈現(xiàn)逐漸減弱趨勢，晶體條上響應(yīng)的局部最大值與二維正態(tài)分布的峰值相互映射。位置譜上所有位置點的分布均呈現(xiàn)混合高斯分布特點，并受到晶體條上產(chǎn)生閃爍光子的權(quán)重因子W影響[8-9]，則：

(2)

圖3為計數(shù)值累加示意圖，位置譜的關(guān)鍵步驟是將數(shù)據(jù)幀列表中每幀數(shù)據(jù)按照封幀格式讀取出事件位置信息，并將位置信息疊加到256×256矩陣的相應(yīng)直方圖上(位置譜根據(jù)探測器、通道編號單獨處理)。

圖3 計數(shù)值累加示意圖Fig.3 Schematic of counting value accumulation

位置譜獲取的程序流程圖如圖4所示，主要是將γ光子打在對應(yīng)探測器各通道上的位置信息繪成圖譜，并按照二維高斯模型累加求和，得到最終的位置譜。

圖4 位置譜獲取的程序流程圖Fig.4 Flow chart of position spectrum acquisition algorithm

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

3.1 TS分割算法的關(guān)鍵要素及參數(shù)計算

為了劃分位置譜中響應(yīng)區(qū)域的離散點并濾除偽峰，設(shè)計了一種TS分割算法程序，本程序的核心思想為：將整個256×256位置譜視為矩陣，將所有點按計數(shù)量排序，通過禁忌表機(jī)制記錄峰值并將其鄰域坐標(biāo)賦值為1，從而快速而不重復(fù)地迭代搜索，將等同于晶體條數(shù)目的禁忌鄰域分割出來[10]，可構(gòu)成晶體覆蓋區(qū)上γ射線作用的響應(yīng)位置分割結(jié)果圖，算法包含以下5大要素。

1) 禁忌表(Tabulist)，數(shù)值為1的坐標(biāo)即為禁止?fàn)顟B(tài)，其數(shù)據(jù)內(nèi)容(禁忌點坐標(biāo)與順序)存儲在Tabulist(晶體位置查找表的基礎(chǔ))變量中。Tabulist是一含有n2個標(biāo)記點的n2行2列的矩陣，按照對應(yīng)計數(shù)值從上到下降序排列。

2) 終止準(zhǔn)則為晶體條的數(shù)量，如13×13晶體，其終止準(zhǔn)則為輸出169個禁忌鄰域后循環(huán)終止。

3) 禁忌長度TL是指每個方塊禁忌鄰域的邊長(即所含位置點個數(shù))，其長度取決于晶體與SiPM間的光導(dǎo)玻璃&光學(xué)膠水等介質(zhì)折射率Nd和晶體陣列總行(列)數(shù)n，其計算公式為：

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

(3)

晶體探頭所用光導(dǎo)介質(zhì)為聚甲基丙烯酸甲脂(一般室溫25 ℃條件下折射率為1.49)，但由于粘合晶體的光學(xué)膠水(有機(jī)硅凝膠，其折射率為1.41)降低了折射率，且BDM在采集數(shù)據(jù)時會發(fā)出大量熱導(dǎo)致光導(dǎo)玻璃折射率下降[11]，故實際折射率約為1.41～1.49。對多種晶體陣列進(jìn)行評估后，得出晶體到SiPM的實際折射率Nd=256/175≈1.462 9，符合溫度及膠水帶來的影響。圖5為7種晶體陣列的禁忌長度變化曲線。在實際的PET探測器應(yīng)用中，常用的4×4、6×6、10×10、12×12、13×13、16×16、20×20等7種晶體陣列的TL分別為53、30、14、11、10、7、4個坐標(biāo)長度(256×256位置譜已是數(shù)據(jù)容量極限，位置點個數(shù)不存在小數(shù)，故取近似為整數(shù))。

圖5 7種晶體陣列的禁忌長度變化Fig.5 Tabu length change of seven kinds of crystal arrays

4) 禁忌鄰域，對應(yīng)湮滅事件響應(yīng)位置，以局部響應(yīng)峰值坐標(biāo)(xi±0.5TL，yi±0.5TL)區(qū)域的矩形。圖6為13×13晶體陣列的禁忌鄰域，13×13晶體的禁忌鄰域為11×11區(qū)域，含121個位置點坐標(biāo)。

圖6 13×13晶體陣列的禁忌鄰域Fig.6 Tabu neighborhood of 13×13 crystal array

5) 特赦準(zhǔn)則，檢測響應(yīng)標(biāo)記譜Ff中所有響應(yīng)標(biāo)記點Ki(xi,yi)是否滿足孤立狀態(tài)，具體實現(xiàn)為(xi+0.5TL,yi)=(xi-0.5TL,yi)=0或(xi,yi+0.5TL)=(xi,yi-0.5TL)=0，若是則實行特赦，解除此點及其禁忌鄰域的禁忌狀態(tài)，繼續(xù)進(jìn)行禁忌搜索更新分割。圖7為偽峰的特赦處理，峰值點延伸的水平/豎直線不與其他任何響應(yīng)標(biāo)記點相交。

圖7 偽峰的特赦處理Fig.7 Amnesty treatment of false peak

3.2 TS分割算法的實現(xiàn)步驟

圖8為TS分割算法流程圖，具體內(nèi)容如下。

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

S1賦值：將位置譜F1的三維位置信息進(jìn)行分割預(yù)處理并整合為256×256矩陣F2。

S2排序：取i=1、j=1，按數(shù)值從大到小排序?qū)2中的坐標(biāo)轉(zhuǎn)化為169行2列的索引矩陣Ix，同時把禁忌鄰域譜Ff賦值為256×256的0矩陣，取矩陣Tabulist賦值為169行2列的0矩陣。

S3禁忌檢查：為判斷峰值是否已被禁忌，記錄Ix索引中數(shù)值排第j行的坐標(biāo)(xj,yj)，判斷Ff矩陣內(nèi)的(xj,yj)坐標(biāo)數(shù)值是否為1，是則繼續(xù)，否則j=j+1并重復(fù)S3步驟。

S4添加禁忌表：將F2中的(xj,yj)坐標(biāo)視為峰值點Ki(即相應(yīng)晶體條的γ射線響應(yīng)區(qū)域局部最大值)，將坐標(biāo)賦值到Tabulist中第j行，將F2中的Ki點用紅色圓圈標(biāo)記。

S5禁忌鄰域：為劃分出晶體條的γ射線響應(yīng)區(qū)域，將矩陣Ff內(nèi)(xj±0.5TL,yj±0.5TL)方形鄰域內(nèi)所有坐標(biāo)賦值為1，視為1個禁忌鄰域Ti。

S6更新結(jié)果圖：輸出此時含所有Ki的響應(yīng)峰值圖，如圖9a所示；輸出此時含所有Ti的禁忌鄰域圖(最終結(jié)果為響應(yīng)位置分割圖)，如圖9b所示。

圖9 正在標(biāo)記響應(yīng)點的響應(yīng)峰值圖(a)和禁忌鄰域圖(b)Fig.9 Respons peak map (a) and tabu area map (b) of response point being marked

S7終止準(zhǔn)則：為確保所有晶體條映射區(qū)域的響應(yīng)皆被分割，判斷i

S8特赦：作為算法負(fù)反饋機(jī)制，排除噪聲導(dǎo)致誤選偽峰為響應(yīng)峰值的情況，判斷是否符合特赦準(zhǔn)則，若符合返回步驟S5，否則結(jié)束。

隨著尋峰迭代和禁忌操作，單次搜索的對象數(shù)會逐步減少，避免了搜索過程陷入循環(huán)或局部最優(yōu)，尋峰操作也會越來越快，時間復(fù)雜度大幅減??；此外，由步驟S8可知，特赦準(zhǔn)則濾除了出現(xiàn)沖突的偽峰情況，避免了過分割。

3.3 構(gòu)造晶體位置查找表

為確定各晶體單元覆蓋區(qū)域內(nèi)γ射線作用位置與實際位置譜光斑位置點的映射關(guān)系，構(gòu)建了晶體位置查找表以便更準(zhǔn)確地確定湮滅事件所在響應(yīng)線的位置。

本方法以晶體條分布情況為約束，將響應(yīng)位置分割圖的各禁忌鄰域編號為1～n2，將Tabulist中所有點按坐標(biāo)逐行、逐列排序填入其中，將各離散位置點歸屬到對應(yīng)禁忌鄰域，輸出所有晶體的位置查找表CPLT，輸出的結(jié)果如圖10所示，從而根據(jù)Tabulist確定每個晶體條對應(yīng)的區(qū)域編號并進(jìn)行湮滅事件響應(yīng)線位置校準(zhǔn)。

圖10 晶體響應(yīng)映射圖(a)及晶體查找表結(jié)果(b)Fig.10 Crystal response map (a) and crystal lookup table result (b)

4 實驗結(jié)果

在具體實驗中，放置裝有活度為1.11×105Bq的18F-FDG注射器在距探測器模塊1.5 cm處位置，使用探測器掃描10 s，每個探測器模塊平均捕捉到約350 000個光子事件數(shù)。程序輸出位置譜耗時0.31 s，輸出位置譜結(jié)果如圖11所示。

a——13×13位置譜；b——6×6位置譜圖11 兩種晶體位置譜結(jié)果Fig.11 Two crystal position spectrum results

為評估位置數(shù)據(jù)準(zhǔn)確度，在溢出率為5%的條件下，將本方法獲取的13×13晶體的位置譜按照水平、豎直方向進(jìn)行了計數(shù)值的投影，所得結(jié)果如圖12所示，通過計算獲得位置譜的13個峰值與12個谷值的峰值比并求均值后可得：水平峰谷比為12.18，豎直峰谷比為9.45，高于現(xiàn)有的大多數(shù)方法生成的位置譜效果[12-13]，通過擬合方法計算投影的半高寬(FWHM)獲得1.7 mm的最佳空間分辨率。

平均每執(zhí)行70.43 ms完成1個BDM的晶體位置譜分割，分割結(jié)果如圖13所示，使用的13×13和6×6晶體陣列分割出169個、36個晶體標(biāo)記區(qū)域，其中每個晶體區(qū)域的響應(yīng)光斑均能清晰顯示，保證了所有晶體條覆蓋區(qū)的γ射線響應(yīng)位置得到映射，無響應(yīng)區(qū)域遺漏，且無噪聲引起偽峰，從結(jié)果上進(jìn)一步證實了算法避免過分割的功能。

a——水平方向投影；b——豎直方向投影圖12 13×13晶體位置譜的峰值及投影曲線Fig.12 Peak value and projection curve of 13×13 crystal position spectrum

a——13×13響應(yīng)峰值圖；b——13×13響應(yīng)位置分割圖；c——6×6響應(yīng)峰值圖；d——6×6響應(yīng)位置分割圖圖13 位置譜分割結(jié)果Fig.13 Segmentation result of position spectrum

為從數(shù)據(jù)上進(jìn)一步證實本方法的效率優(yōu)勢，將TS分割算法與分水嶺分割算法[14-15]的時間復(fù)雜度進(jìn)行了對比(N為陣列中晶體條數(shù)目)，結(jié)果列于表1。分水嶺分割算法理論上運算次數(shù)為N(1+N)+p次(p對應(yīng)后續(xù)待處理的偽位置點個數(shù))，時間復(fù)雜度為O(N2)，運行時間約為1～3 s。本算法的運算次數(shù)為N(1+k)次(k為排序N以內(nèi)被特赦準(zhǔn)則濾除的偽峰數(shù)量)，即時間復(fù)雜度為O(N)，對應(yīng)的13×13晶體陣列運算時間為0.070 43 s。

表1 兩種方法的運行時間對比Table 1 Comparison of running time of two methods

5 結(jié)論

本文開發(fā)了一種以禁忌搜索算法為核心，結(jié)合UDP數(shù)據(jù)幀、二維高斯模型及光導(dǎo)折射率的位置信息分割校準(zhǔn)方法，該方法獲得的位置譜峰谷比優(yōu)于大多數(shù)其他方法，分割效率高于分水嶺算法，且避免了過分割現(xiàn)象，濾除了噪聲引發(fā)的偽峰，最終實現(xiàn)了系統(tǒng)空間分辨率的提高。因此，本方法在自動分割強(qiáng)度分布不均或形變失真程度高的位置譜具有較大優(yōu)勢，具有校準(zhǔn)響應(yīng)事件位置準(zhǔn)確性的負(fù)反饋機(jī)制，分割與校準(zhǔn)可指導(dǎo)后續(xù)符合處理，從而準(zhǔn)確地確定湮滅事件所在響應(yīng)線的位置，同時可作為測試大量晶體陣列乃至PET系統(tǒng)的評估標(biāo)準(zhǔn)。本方法已在南昌大學(xué)PET實驗室的智能PET影像系統(tǒng)中表現(xiàn)出良好的適用性。