一體化 TOF-PET/MR 是獲得高精度正電子斷層掃描（PET）圖像同時一并獲得磁共振（MR）圖像的多模態成像設備，具備飛行時間（TOF）功能。其中 PET 系統利用了正電子湮滅產生的兩個方向 180° 的光子被探測器接收的時間差，來縮小響應線上湮滅發生位置的范圍。由于每個探測晶體自身硬件條件不同和外部噪聲影響，不同晶體的時間偏移不同，所以需要對一體化 TOF-PET/MR 進行精確的時間校正來使得系統正常工作。我們采取三種不同原理的方法對系統進行時間校正，第一種為利用幾何方法構建扇形束模型，來擬合時間偏移的高斯分布，從而求晶體時間偏差的迭代方法，簡稱扇形束法；第二種為構建超定方程組的解，并利用 L1 范數最小化求解的方法，簡稱 L1-norm 法；第三種利用大量數據構建直方圖，擬合尋峰后的結果構成超定方程組，利用 L2 范數最小化進行求解，簡稱 L2-norm 法。本文對這三種方法所需數據量和計算時間進行了比較，對時間校正后系統重建圖像進行了分析。為了減小采集數據時放射源位置偏移對校正結果的干擾，我們設計了位置校正算法對采集數據進行預處理，該算法能直接計算出筒源擺放位置，并減少因為偏心產生的圖像偽影。實驗結果證明，L2-norm 法噪聲比較小，但是計算速度慢；L1-norm 法擁有最快的計算速度，但是圖像質量較差；扇形束法較其他兩種方法擁有更高的圖像質量，尤其是對于微小病灶的探查能力更強，因此在一體化 TOF-PET/MR 中使用扇形束法進行時間校正最佳。
引用本文： 曾天翼, 楊卉, 曹拓宇, 胡凌志, 褚旭, 呂新宇, 陳群. 一體化 TOF-PET/MR 時間校正方法的對比研究. 生物醫學工程學雜志, 2019, 36(6): 1003-1011. doi: 10.7507/1001-5515.201809044 復制
引言
正電子斷層掃描/磁共振成像一體化設備（positron emission tomography/magnetic resonance imaging，PET/MR）是一種新的融合兩種領先成像方法的多模態成像系統，可以綜合 PET 系統對示蹤劑的高敏感度和 MR 系統對軟組織的高對比度[1]。目前，PET/MR 在乳腺癌、腹部及盆腔腫瘤和淋巴系統腫瘤等腫瘤診斷，神經系統疾病研究，以及心血管系統疾病中發揮著重要作用[2]。傳統 PET 圖像重建中認為湮滅均勻發生在響應線上，通過濾波反投影或正向反向迭代投影重建算法得到圖像，對比度較低，圖像空間分辨率較差。PET 飛行時間技術（time of flight，TOF）通過獲得一個符合事件內兩個光子被探測器接收的時間差，確定發生湮滅事件的位置區域。TOF 技術提高了定位精度，縮小了定位范圍，從而提高圖像對比度，改善圖像質量，在三維采集模式下能夠降低隨機計數率和散射計數率[3-4]。TOF-PET 技術是否可行，很大程度上取決于飛行時間的探測是否準確。對于 TOF-PET 系統，若有一正電子點源在視野中心位置，理論上任一響應線的 TOF 值應該呈現中心為 0 的高斯分布。受定時誤差影響，實際 TOF 差分布中心存在偏差，因此必須通過時間校正技術修正其影響。偏差產生的原因由晶體閃爍特性、脈沖識別差異、信號放大過程、光路長度、電子學延遲等因素共同作用，因此每個晶體的偏差都有可能不同。隨著探測器使用時長和環境改變，偏差數值也發生改變[5]，因此要周期性地對 TOF-PET 系統進行時間校正，而一體化 TOF-PET/MR 系統由于 MR 的高磁場會對 PET 穩定性造成影響，因此需要更加頻繁地進行時間校正。為了獲得更高質量的 PET 圖像以及 PET-MR 融合圖像，有必要對系統進行速度快、準確性高的時間校正。
目前，TOF-PET 系統時間校正的方法可以分為兩類。第一類方法把正電子源放在探測器中間，并用一個能夠探測正電子湮滅前逸出能量的探測器，探測正電子損失的動能以確定電子湮滅時間，作為校正參考時間。單光子采集模式下，對于在探測器環中間的正電子源，所有晶體到它的距離相同，光子被 PET 探測器接收后，晶體反饋回飛行時間。統計每一個晶體收到單光子的飛行時間，結合參考時間，就可以計算出每一個晶體的時間偏差[6]。
第二類方法同樣把正電子源放在探測器中央位置，采集模式為符合模式，即收集響應線（line of response，LOR）數據，不需要使用額外的探測器。將得到的 LOR 數據進行處理計算得到每一個晶體的時間偏差。中心放置源的形式有：探測器中心軸位置的線源，點源，以及中心軸與探測器中心軸重合的均勻筒源[7-9]。處理數據計算偏差的算法有很多種，包括扇形束計算法、L1 范數最小化法（以下簡稱 L1-norm 法）和 L2 范數最小化法（以下簡稱 L2-norm 法）等。
因為第一類方法需要額外增加探測逸出能量的探測器，且額外探測器擺放位置的準確程度對時間校正影響很大，因此本文應用了第二類方法中的扇形束法、L1-norm 法和 L2-norm 法對一體化 TOF-PET/MR 進行時間校正，并提出了一種對筒源放置偏差的校正方法。評估三種時間校正方法的優劣使用了美國電器制造商協會（National Electrical Manufactures Association，NEMA）圖像質量（image quality，IQ）模體的重建圖像和復旦大學附屬中山醫院的十位腫瘤患者的全身圖像。
1 方法
1.1 扇形束法
在扇形束法中假設，PET 探測器收集到的中央筒狀正電子源的所有光子符合事件，都是真實發生的符合事件，所以扇形束即指代每一個晶體與有可能和它一起產生符合事件的所有晶體組成的集合。如圖 1a 所示，以晶體連線與筒源相切為范圍，可知虛線對應范圍內的晶體形成扇形束，晶體的時間偏移量應該是統計其扇形束中探測到的所有 TOF 值的高斯分布中心距離 0 時間的偏差。計算得出每個晶體的時間偏移后，因為每個晶體受多條 LOR 影響，因此該晶體的時間偏移變化會讓扇形束中多條 LOR 的 TOF 改變，將改變后的 TOF 值替換原來實際測量得到的 TOF 值，因為考慮了晶體時間偏移，會更逼近真實的 TOF 值。迭代進行該過程，可以逐次逼近所有晶體真實的相對時間偏移。
圖1
扇形束法以及位置校正原理示意圖
a. 構建扇形面示意圖；b. 筒源位置偏差對 TOF 造成的影響
Figure1. Schematic diagram of fan-beam method and position calibrationa. schematic diagram of fan-beam area；b. influence of cylindrical source position deviation on TOF
1.2 L1-norm 法原理
L1-norm 法認為所有探測器獲得的符合事件的信息可以用一個超定方程組來歸納，超定方程組指方程個數大于未知量個數的方程組。對于放置在探測器中心的點源，兩個晶體 l、k 之間 LOR 的理論 TOF 差為 0。這時兩個晶體的實際 TOF 差可以表示為：
 |
公式（1）可以表示成超定方程組形式：
 |
其中，符合事件總數為 m，晶體總數為 n，C為系數矩陣，每一行對應一個符合事件，每一行是公式 1 中 tl 和 tk 前面的系數，因此 C 內每行元素值幾乎都為 0，除了兩個元素分別為 1 和 ? 1。利用 L1 范數將方程轉化為：
 |
L1 范數最小化方法已經被廣泛應用在圖像去噪、信號壓縮和監督學習領域[10]，原因是 L1 范數的特征值往往較為稀疏，這使得該方法求解超定方程時，對噪聲敏感度較低。交替方向乘子法（alternating direction method of multipliers，ADMM）適用于求解約束問題，已被證明可以解決 L1 范數最小化問題[11]。ADMM 法迭代收斂后解出的答案即為每一個晶體的時間偏移。
1.3 L2-norm 法原理
L2-norm 法中，對采集筒源獲得的每一種 LOR 的 TOF 差進行統計，每一種 LOR 構建一個 TOF 差直方圖，找到其峰值與 0 時間的偏移距離，構建方程個數等于 LOR 種數的方程組。利用極大似然估計法可將方程組轉換為公式（4）。
 |
其中 c 為所求每個晶體的時間偏移，A 是 LOR 方程組的系數矩陣，與 L1 norm 方法中的 C 矩陣類似，取值為 ? 1、0、1，
是實際測量 TOF 差，
是實際測量噪聲。解該方程組實際上是求二范數 
的最小值，一般解法是
 |
公式（5）中 t 為所求晶體偏移時間，A 和 
與公式（4）相同。利用迭代單邊坐標下降法去估計逆矩陣，可以快速解 L2 范數最小化問題[12-13]。得到的初步時間校正結果進行迭代重新尋峰后，達到收斂即為最終晶體時間偏移結果。
1.4 位置校正方法原理
在實際測量中，筒源不可避免與探測器軸心有肉眼難以發現的偏差，圖 1b 表示了筒源偏心或者歪斜產生的結果。筒源偏心時，LOR 與其截交線中心點發生移動，而中心點處即為 TOF 高斯分布中心。所以，筒源偏心影響了 LOR 對應兩方晶體的時間偏移。為了消除偏差對時間校正的影響，需要計算響應線中點和截交線中點的距離 d。
我們提出了一種利用最小二乘法擬合得到實際筒源擺放軸心位置的方法，以期減小位置偏差對校正結果的影響。在該方法中，首先將筒源沿軸向平分為 14 個薄層，在每個薄層內，讀取探測器接收的兩個晶體都落在該薄層的筒源符合事件。使用 Fong 等[14]提出的最小平方殘差（least squares minimal residual，LSMR）算法，其核心是利用最小二乘法迭代求解。該算法收斂結果為該薄面的中心點（x，y）坐標。利用空間直線擬合的方法，得到軸心直線在空間的表達式。
2 實驗結果及分析
我們使用上海聯影醫療科技有限公司生產的 uPMR790 一體化 TOF-PET/MR 設備來獲得原始數據。將灌有放射性活度 3 mCi 氟代脫氧葡萄糖溶液、半徑為 106.5 mm、長度為 319.2 mm 的筒狀體模放置在視野中心，筒源中軸與探測器中心軸重合，探測器接收光子模式為符合模式，符合時間窗寬為 3 ns，時間戳長度為 10 ns。L2-norm 法需要采集 40 min，得到大小 40 GB 的原始數據；扇形束法需要采集 90 s，得到大小 1.5 GB 的原始數據；L1-norm 法需要采集 60 s，得到大小 1 GB 左右的原始數據。
2.1 晶體時間偏移隨迭代次數的變化情況
L2-norm 法由于需要每種 LOR 大量計數獲得統計直方圖，所需數據量非常大，一般需要 3.03 × 109 個真符合事件。扇形束法需要的符合事件數量在 108 量級，根據經驗選取了 1.15 × 108 個真符合事件進行計算，結果與預期符合。L1-norm 法因為對噪聲不敏感，且能自動選擇特征值，因此僅需要 7 × 107 個符合事件，比扇形束法所需數據量更小。
任意選取晶體，對其每一次迭代的時間偏移進行追蹤，根據時間偏移變化值隨迭代次數的變化情況判斷收斂性。結果表明，在迭代初期，時間偏移變化劇烈，之后在某一范圍內波動。扇形束法的收斂結果最終在[? 6，+ 6] ps 范圍內波動，L1-norm 法和 L2-norm 法的收斂結果最終在[? 3，+ 3] ps 范圍內波動。與扇形束法對比，可知 L1-norm 法和 L2-norm 法收斂性更強，波動范圍更小。
我們用 CPU 主頻 2.53 GHz、內存容量 14 GB、硬盤容量 1 TB 的計算機運行程序。扇形束法取 10 次迭代后結果輸出，采集時間加上計算時間共需 214.68 s；L1-norm 法迭代次數達到 50 次后，可以認為達到收斂，采集時間加上計算時間共需 295.75 s；L2-norm 法迭代次數 5 次達到收斂，采集時間加上計算時間共需 45 min。
2.2 三種方法時間校正結果分析
如圖 2a、2b、2c 分別是扇形束法、L1-norm 法、L2-norm 法計算得到的每個晶體的時間偏移，三種方法計算所得結果在晶體分布圖上無直觀區別。圖 2d 是扇形束法與 L2-norm 法計算差異的直方圖，對于絕大多數晶體，兩種方法計算結果接近，但整體存在一個偏移量，在個別晶體上有較大差異。分析直方圖分布可知，結果差異呈正態分布，擬合后期望為 ? 259.95 ps，正態分布半高寬（full width at half maxima，FWHM）為 112.88 ps。圖 2e 是扇形束法與 L1-norm 法結果差異直方圖，二者計算結果最接近，是以 64.59 ps 為中心、FWHM = 181.65 ps 的正態分布。圖 2f 是 L1-norm 法與 L2-norm 法結果差異直方圖，分析兩者整體偏差，是以 ? 222.49 ps 為中心、FWHM = 238.01 ps 的正態分布。不同方法計算結果存在整體偏移的原因，是由于三種方法所求的并不是真的時間偏移值，而是使兩兩晶體時間校準的值，即兩晶體分別加上各自時間偏移后相減時間差等于實際 TOF 差。由于一體化 TOF-PET/MR 使用 TOF 差推測湮滅發生位置進行重建，所以晶體時間整體偏移不會對重建圖像產生影響。
圖2
三種方法計算的時間偏移與互相比較結果
a. 扇形束法按每個晶體顯示的計算結果；b. L1-norm 法按每個晶體顯示的計算結果；c. L2-norm 法按每個晶體顯示的計算結果；d. 扇形束法結果減去 L2-norm 法結果的分布直方圖，為類高斯分布；e. 扇形束法結果減去 L1-norm 法結果的分布直方圖，為類高斯分布；f. L1-norm 法結果減去 L2-norm 法結果的分布直方圖，為類高斯分布
Figure2. Timing offsets and comparison results calculated by three methodsa. calculation results of fan-beam method；b. calculation results of L1-norm method；c. calculation results of L2-norm method；d. distribution histogram of the calculation results of fan-beam subtracting L2-norm；e. distribution histogram of the calculation results of fan-beam subtracting L1-norm；f. distribution histogram of the calculation results of L1-norm subtracting L2-norm
2.3 位置校正準確性分析
為了確定位置校正算法準確性，我們使用基于有序子集最大期望值法（ordered subset expectation maximum，OSEM）的 TOF 重建獲得高清圖像[15]。將得到的圖像進行閾值分割得到筒源輪廓，在輪廓內部用兩條線逼近的方法，找到對應切面的圓心坐標，作為參考值。利用一個肉眼幾乎不可見的偏心筒源進行位置校正準確性驗證，結果見表 1。
表1
位置校正準確性驗證結果
Table1.
Accuracy of position correction
由于重建圖像在筒源首尾存在邊界模糊，導致首尾偏差比較大，我們取中間層進行分析，可知軸心位置計算偏差百分比在 15% 以內。軸心在 x 坐標偏心 1 cm 情況下，誤差在 1.5 mm 內。對于 TOF 校正結果的偏差為 10 ps，這個誤差是可以接受的。由于在實際采集中無法定量描述可能發生的肉眼擺放誤差，我們在實驗中采集了一組偏心較大筒源數據，偏移范圍在 4 cm 左右，在未進行位置校正情況下，對系統進行時間校正。得到校正結果后對該偏心筒源進行重建，圖像產生很大偽影。再使用本文提出的位置校正算法后進行時間校正，同樣重建條件下，圖像偽影有很大改善。結果表明我們所使用的位置校正算法，能減小因為位置偏差對時間校正的誤差。
2.4 NEMA IQ 實驗圖像質量分析
NEMA 定義了圖像質量的測試標準，該標準給出了臨床成像條件下圖像質量的分析方法。對比度復原系數（contrast recovery coefficient，CRC）表示 NEMA IQ 中球體在背景下的對比度還原情況，通過計算 NEMA IQ 模體內 4 個熱球和 2 個冷球的 CRC，可以反映 PET 圖像質量。
我們對同一個 NEMA IQ 模體采集的符合數據，對未進行時間校正以及分別采用 L2-norm 法、扇形束法和 L1-norm 法時間校正后的數據進行重建。重建使用前述 OSEM 算法，參數設置相同的情況下，取第三次迭代圖像結果的同一切面進行分析。同一機器同樣條件重復進行三次實驗，取平均值作統計圖，誤差條為正負標準差，結果如圖 3 所示。圖中黑色圓形表示熱球，半徑從 10 mm 到 22 mm 逆時針增加；白色圓形表示冷球，半徑為 28 mm 和 37 mm，逆時針增加；中間圓形為空氣，在圖表中記為 Lung。三種時間校正方法所得圖像均清晰無偽影，可以清楚看到熱球和冷球位置，與未校正結果對比明顯。定量分析圖 3 中各方法圖像的 CRC 值，未時間校正圖像的 CRC 與三種不同方法時間校正圖像 CRC 有明顯差異，同時噪聲也更大。扇形束法噪聲平均值為 6.85%，略大于 L1-norm 法的 5.92% 和 L2-norm 法的 5.49%。如圖 3 中右邊圖表可知，扇形束法的熱球 CRC 平均值 68.62% 略高于 L1-norm 法的熱球 CRC 平均值 67.58% 和 L2-norm 法的熱球 CRC 平均值 67.28%。這表明在半徑較小的熱球區域，扇形束法圖像質量較好，顯示了對小病灶的探查潛力。
圖3
NEMA IQ 模體重建圖像以及 CRC 值與噪聲值統計結果
Figure3.
Reconstructed images of NEMA IQ phantom and its statistical results of CRC and noise
2.5 腫瘤患者全身 PET 圖像的比較
隨機選擇了 10 例復旦大學附屬中山醫院 2018 年接受一體化 TOF-PET/MR 掃描的腫瘤患者，根據復旦大學附屬中山醫院和上海聯影醫療科技有限公司臨床合作協議，所有數據使用均已獲授權。用不同方法時間校正后的系統分別重建患者原始數據得到全身 PET 圖像，采用 PET 圖像中腫瘤區域基于去脂體重的標準攝取值（standard uptake value based on lean body mass，SUV lbm）比較結果。在同一正切面圖像中圈出感興趣區域（region of interest，ROI），比較 ROI 中 SUV 的平均值。SUV 在圖像中表現為體素的信號強度值，與病灶濃度 C病灶（kBq/mL）和注射劑量 D注射（MBq/kg）有關[16]，具體關系為：
 |
圖 4a 為未校正圖像，其他三幅分別為使用不同時間校正方法的晶體時間偏移數據后重建獲得的患者圖像。從圖 4 分析可知，經過時間校正后的重建圖像，對湮滅發生的位置定位更加準確，可以清晰地看到腫瘤區域，腫瘤區域的 SUV 平均值增大。
圖4
患者腫瘤圖像 SUV 對比
a. 未校正時患者全身冠狀面圖像；b. L2-norm 法校正后患者全身冠狀面圖像以及 SUV 值；c. 扇形束法校正后患者全身冠狀面圖像以及 SUV 值；d. L1-norm 法校正后患者全身冠狀面圖像以及 SUV 值
Figure4. Comparison of SUV calculated from patient’s tumor imagesa. uncalibrated whole-body coronal image；b. whole-body coronal image with L2-norm calibration；c. whole-body coronal image with fan-beam calibration；d. whole-body coronal image with L1-norm calibration
將全部 10 例患者腫瘤區域 SUV 平均值以散點圖形式表示，并對每種方法的平均值進行直方圖表示，誤差條是十個患者數據的標準差，如圖 5 所示。L2-norm 法由于使用較多數據，更加真實地反映了時間偏差，使得腫瘤區域 SUV 平均值稍高。利用單因素方差分析四組結果發現總體差異有統計學意義，其中未校正圖像和三種不同方法時間校正后圖像均有顯著差異，而三種校正方法兩兩進行方差分析其差異無統計學意義。檢驗結果表明，未校正圖像 SUV 平均值與經過時間校正后圖像 SUV 平均值差異顯著，而三種校正方法之間 SUV 平均值差異不明顯。
圖5
十位患者腫瘤圖像 SUV 平均值統計圖
直方圖是不同方法校正時所有患者 SUV 的平均值，散點圖為每個患者的具體腫瘤區域 SUV 平均值。*：
histogram indicates average value of ten patients’ SUV with different calibration，and scatter plot indicates specific SUV mean value of every patient. *：
3 討論與結論
TOF-PET 自 20 世紀 80 年代被提出，隨著新型閃爍晶體的應用，時間分辨率越來越好，TOF-PET 重建技術也隨之發展。TOF 系統在探測過程中能定位湮滅發生的區域，提高反投影的空間精度[17]。TOF 技術最大程度地發揮了摻鈰硅酸镥（lutetium oxyorthosilicate，LSO）閃爍晶體在 PET 中的表現，在減少隨機事件和系統噪聲以及加快掃描時間上有重要的價值[18]。對于一體化 TOF-PET/MR 系統，由于需要更加頻繁的時間校正，算法的優劣將對系統 TOF 表現造成很大影響。
我們在一體化 TOF-PET/MR 系統中利用扇形束法、L2-norm 法和 L1-norm 法進行時間校正，并考慮了采集時筒源擺放偏心情況的位置校正。GE 醫療在 2014 年研發出帶有 TOF 功能的一體化 PET/MR，其中使用了磁場條件下電子性能更好的硅光電二極管，其時間分辨率達到 390 ps。Grant 等[19]對該系統測試表明，TOF 能有效增加重建圖像的信噪比，具有更好圖像質量、更短掃描時間和更少注射劑量的潛力。Minamimoto 等[20]對一體化 TOF-PET/MR 的臨床表現進行了評估，擁有良好時間分辨率的 PET/MR 與無 TOF 功能 PET/MR、TOF-PET/CT 以及無 TOF 功能 PET/CT 相比，減少了圖像偽影，提高了 PET 臨床圖像質量。隨著 TOF 技術的發展，高時間分辨率的一體化 TOF-PET/MR 將會成為主流發展方向，因此能夠發揮硬件極限時間分辨率的時間校正算法愈發重要。本文比較了三種不同方法的效果，證明了時間校正的必要性，對時間校正方法的選取具有重要意義。
自 Luo 等[9]提出時間校正問題，已經有多種時間校正的方法被提出。主流方法之一是 Lenox 等[8]提出的扇形束法，通過計算每個晶體的時間偏差直方圖，對于圓柱狀模體，形成類似扇形面的 LOR 集合。另一主流方法是構建超定方程組解每個晶體的時間偏差，由 Reynolds 等[13]最先提出，將問題轉化為最小二乘問題，并利用 L2-norm 最小化的方法求解。后來，Freese 等[10]基于 L2-norm 最小化法，提出用 L1-norm 替代，計算得到魯棒性好的結果。然而，目前并沒有針對 TOF-PET 的各種主流時間校正方法效果的對比。本文對比了三種時間校正方法的整體流程所需時間和圖像質量，分析了各種方法的優劣，對一體化 TOF-PET/MR 時間校正方法提出了建議，且本文方法也可以用于其他 TOF-PET 系統。但本文所提三種方法，仍存在散射和隨機事件對校正結果影響較大的問題，在未來的研究中，應該去探究如何減少數據采集中的散射和隨機事件，或者研究新的時間校正算法，完全消除散射和隨機產生的影響。本文采取的位置校正方法不僅可以減少位置擺放偏差帶來的圖像偽影，而且能根據計算出的筒源位置，判斷筒源擺放是否準確，如果偏差非常大，產生錯誤，可以重新擺位進行采集。但是位置校正算法在部分切面上誤差大于 10%，這樣的誤差也使得位置校正無法完全消除偏心帶來的圖像偽影，如何提高該方法的準確性，獲取更真實的放射源位置，需要進一步探索。
綜上所述，扇形束法和 L1-norm 法與 L2-norm 法相比，所用數據量大幅減少，計算速度提升，使商品化 TOF-PET 周期性時間校正的工作量大幅降低，尤其對于需要更頻繁進行時間校正的一體化 TOF-PET/MR 來說具有很大優勢。扇形束法與其他方法相比有更佳的圖像質量，尤其是在較小病灶的探查上表現出更大潛力，因此綜合考慮時間成本和圖像質量，一體化 TOF-PET/MR 最適合選用扇形束法進行時間校正。本文在時間校正過程中提出了位置校正算法，對數據進行預處理，能發現肉眼無法分辨的擺放偏差并有效改善筒源擺放誤差對時間校正結果的影響。位置校正算法結合扇形束法將會使系統的時間分辨率提升，可廣泛應用在商品化和研究用途的一體化 TOF-PET/MR 中。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
正電子斷層掃描/磁共振成像一體化設備（positron emission tomography/magnetic resonance imaging，PET/MR）是一種新的融合兩種領先成像方法的多模態成像系統，可以綜合 PET 系統對示蹤劑的高敏感度和 MR 系統對軟組織的高對比度[1]。目前，PET/MR 在乳腺癌、腹部及盆腔腫瘤和淋巴系統腫瘤等腫瘤診斷，神經系統疾病研究，以及心血管系統疾病中發揮著重要作用[2]。傳統 PET 圖像重建中認為湮滅均勻發生在響應線上，通過濾波反投影或正向反向迭代投影重建算法得到圖像，對比度較低，圖像空間分辨率較差。PET 飛行時間技術（time of flight，TOF）通過獲得一個符合事件內兩個光子被探測器接收的時間差，確定發生湮滅事件的位置區域。TOF 技術提高了定位精度，縮小了定位范圍，從而提高圖像對比度，改善圖像質量，在三維采集模式下能夠降低隨機計數率和散射計數率[3-4]。TOF-PET 技術是否可行，很大程度上取決于飛行時間的探測是否準確。對于 TOF-PET 系統，若有一正電子點源在視野中心位置，理論上任一響應線的 TOF 值應該呈現中心為 0 的高斯分布。受定時誤差影響，實際 TOF 差分布中心存在偏差，因此必須通過時間校正技術修正其影響。偏差產生的原因由晶體閃爍特性、脈沖識別差異、信號放大過程、光路長度、電子學延遲等因素共同作用，因此每個晶體的偏差都有可能不同。隨著探測器使用時長和環境改變，偏差數值也發生改變[5]，因此要周期性地對 TOF-PET 系統進行時間校正，而一體化 TOF-PET/MR 系統由于 MR 的高磁場會對 PET 穩定性造成影響，因此需要更加頻繁地進行時間校正。為了獲得更高質量的 PET 圖像以及 PET-MR 融合圖像，有必要對系統進行速度快、準確性高的時間校正。
目前，TOF-PET 系統時間校正的方法可以分為兩類。第一類方法把正電子源放在探測器中間，并用一個能夠探測正電子湮滅前逸出能量的探測器，探測正電子損失的動能以確定電子湮滅時間，作為校正參考時間。單光子采集模式下，對于在探測器環中間的正電子源，所有晶體到它的距離相同，光子被 PET 探測器接收后，晶體反饋回飛行時間。統計每一個晶體收到單光子的飛行時間，結合參考時間，就可以計算出每一個晶體的時間偏差[6]。
第二類方法同樣把正電子源放在探測器中央位置，采集模式為符合模式，即收集響應線（line of response，LOR）數據，不需要使用額外的探測器。將得到的 LOR 數據進行處理計算得到每一個晶體的時間偏差。中心放置源的形式有：探測器中心軸位置的線源，點源，以及中心軸與探測器中心軸重合的均勻筒源[7-9]。處理數據計算偏差的算法有很多種，包括扇形束計算法、L1 范數最小化法（以下簡稱 L1-norm 法）和 L2 范數最小化法（以下簡稱 L2-norm 法）等。
因為第一類方法需要額外增加探測逸出能量的探測器，且額外探測器擺放位置的準確程度對時間校正影響很大，因此本文應用了第二類方法中的扇形束法、L1-norm 法和 L2-norm 法對一體化 TOF-PET/MR 進行時間校正，并提出了一種對筒源放置偏差的校正方法。評估三種時間校正方法的優劣使用了美國電器制造商協會（National Electrical Manufactures Association，NEMA）圖像質量（image quality，IQ）模體的重建圖像和復旦大學附屬中山醫院的十位腫瘤患者的全身圖像。
1 方法
1.1 扇形束法
在扇形束法中假設，PET 探測器收集到的中央筒狀正電子源的所有光子符合事件，都是真實發生的符合事件，所以扇形束即指代每一個晶體與有可能和它一起產生符合事件的所有晶體組成的集合。如圖 1a 所示，以晶體連線與筒源相切為范圍，可知虛線對應范圍內的晶體形成扇形束，晶體的時間偏移量應該是統計其扇形束中探測到的所有 TOF 值的高斯分布中心距離 0 時間的偏差。計算得出每個晶體的時間偏移后，因為每個晶體受多條 LOR 影響，因此該晶體的時間偏移變化會讓扇形束中多條 LOR 的 TOF 改變，將改變后的 TOF 值替換原來實際測量得到的 TOF 值，因為考慮了晶體時間偏移，會更逼近真實的 TOF 值。迭代進行該過程，可以逐次逼近所有晶體真實的相對時間偏移。
圖1
扇形束法以及位置校正原理示意圖
a. 構建扇形面示意圖；b. 筒源位置偏差對 TOF 造成的影響
Figure1. Schematic diagram of fan-beam method and position calibrationa. schematic diagram of fan-beam area；b. influence of cylindrical source position deviation on TOF
1.2 L1-norm 法原理
L1-norm 法認為所有探測器獲得的符合事件的信息可以用一個超定方程組來歸納，超定方程組指方程個數大于未知量個數的方程組。對于放置在探測器中心的點源，兩個晶體 l、k 之間 LOR 的理論 TOF 差為 0。這時兩個晶體的實際 TOF 差可以表示為：
|
公式（1）可以表示成超定方程組形式：
|
其中，符合事件總數為 m，晶體總數為 n，C為系數矩陣，每一行對應一個符合事件，每一行是公式 1 中 tl 和 tk 前面的系數，因此 C 內每行元素值幾乎都為 0，除了兩個元素分別為 1 和 ? 1。利用 L1 范數將方程轉化為：
|
L1 范數最小化方法已經被廣泛應用在圖像去噪、信號壓縮和監督學習領域[10]，原因是 L1 范數的特征值往往較為稀疏，這使得該方法求解超定方程時，對噪聲敏感度較低。交替方向乘子法（alternating direction method of multipliers，ADMM）適用于求解約束問題，已被證明可以解決 L1 范數最小化問題[11]。ADMM 法迭代收斂后解出的答案即為每一個晶體的時間偏移。
1.3 L2-norm 法原理
L2-norm 法中，對采集筒源獲得的每一種 LOR 的 TOF 差進行統計，每一種 LOR 構建一個 TOF 差直方圖，找到其峰值與 0 時間的偏移距離，構建方程個數等于 LOR 種數的方程組。利用極大似然估計法可將方程組轉換為公式（4）。
|
其中 c 為所求每個晶體的時間偏移，A 是 LOR 方程組的系數矩陣，與 L1 norm 方法中的 C 矩陣類似，取值為 ? 1、0、1， 是實際測量 TOF 差， 是實際測量噪聲。解該方程組實際上是求二范數 的最小值，一般解法是
|
公式（5）中 t 為所求晶體偏移時間，A 和 與公式（4）相同。利用迭代單邊坐標下降法去估計逆矩陣，可以快速解 L2 范數最小化問題[12-13]。得到的初步時間校正結果進行迭代重新尋峰后，達到收斂即為最終晶體時間偏移結果。
1.4 位置校正方法原理
在實際測量中，筒源不可避免與探測器軸心有肉眼難以發現的偏差，圖 1b 表示了筒源偏心或者歪斜產生的結果。筒源偏心時，LOR 與其截交線中心點發生移動，而中心點處即為 TOF 高斯分布中心。所以，筒源偏心影響了 LOR 對應兩方晶體的時間偏移。為了消除偏差對時間校正的影響，需要計算響應線中點和截交線中點的距離 d。
我們提出了一種利用最小二乘法擬合得到實際筒源擺放軸心位置的方法，以期減小位置偏差對校正結果的影響。在該方法中，首先將筒源沿軸向平分為 14 個薄層，在每個薄層內，讀取探測器接收的兩個晶體都落在該薄層的筒源符合事件。使用 Fong 等[14]提出的最小平方殘差（least squares minimal residual，LSMR）算法，其核心是利用最小二乘法迭代求解。該算法收斂結果為該薄面的中心點（x，y）坐標。利用空間直線擬合的方法，得到軸心直線在空間的表達式。
2 實驗結果及分析
我們使用上海聯影醫療科技有限公司生產的 uPMR790 一體化 TOF-PET/MR 設備來獲得原始數據。將灌有放射性活度 3 mCi 氟代脫氧葡萄糖溶液、半徑為 106.5 mm、長度為 319.2 mm 的筒狀體模放置在視野中心，筒源中軸與探測器中心軸重合，探測器接收光子模式為符合模式，符合時間窗寬為 3 ns，時間戳長度為 10 ns。L2-norm 法需要采集 40 min，得到大小 40 GB 的原始數據；扇形束法需要采集 90 s，得到大小 1.5 GB 的原始數據；L1-norm 法需要采集 60 s，得到大小 1 GB 左右的原始數據。
2.1 晶體時間偏移隨迭代次數的變化情況
L2-norm 法由于需要每種 LOR 大量計數獲得統計直方圖，所需數據量非常大，一般需要 3.03 × 109 個真符合事件。扇形束法需要的符合事件數量在 108 量級，根據經驗選取了 1.15 × 108 個真符合事件進行計算，結果與預期符合。L1-norm 法因為對噪聲不敏感，且能自動選擇特征值，因此僅需要 7 × 107 個符合事件，比扇形束法所需數據量更小。
任意選取晶體，對其每一次迭代的時間偏移進行追蹤，根據時間偏移變化值隨迭代次數的變化情況判斷收斂性。結果表明，在迭代初期，時間偏移變化劇烈，之后在某一范圍內波動。扇形束法的收斂結果最終在[? 6，+ 6] ps 范圍內波動，L1-norm 法和 L2-norm 法的收斂結果最終在[? 3，+ 3] ps 范圍內波動。與扇形束法對比，可知 L1-norm 法和 L2-norm 法收斂性更強，波動范圍更小。
我們用 CPU 主頻 2.53 GHz、內存容量 14 GB、硬盤容量 1 TB 的計算機運行程序。扇形束法取 10 次迭代后結果輸出，采集時間加上計算時間共需 214.68 s；L1-norm 法迭代次數達到 50 次后，可以認為達到收斂，采集時間加上計算時間共需 295.75 s；L2-norm 法迭代次數 5 次達到收斂，采集時間加上計算時間共需 45 min。
2.2 三種方法時間校正結果分析
如圖 2a、2b、2c 分別是扇形束法、L1-norm 法、L2-norm 法計算得到的每個晶體的時間偏移，三種方法計算所得結果在晶體分布圖上無直觀區別。圖 2d 是扇形束法與 L2-norm 法計算差異的直方圖，對于絕大多數晶體，兩種方法計算結果接近，但整體存在一個偏移量，在個別晶體上有較大差異。分析直方圖分布可知，結果差異呈正態分布，擬合后期望為 ? 259.95 ps，正態分布半高寬（full width at half maxima，FWHM）為 112.88 ps。圖 2e 是扇形束法與 L1-norm 法結果差異直方圖，二者計算結果最接近，是以 64.59 ps 為中心、FWHM = 181.65 ps 的正態分布。圖 2f 是 L1-norm 法與 L2-norm 法結果差異直方圖，分析兩者整體偏差，是以 ? 222.49 ps 為中心、FWHM = 238.01 ps 的正態分布。不同方法計算結果存在整體偏移的原因，是由于三種方法所求的并不是真的時間偏移值，而是使兩兩晶體時間校準的值，即兩晶體分別加上各自時間偏移后相減時間差等于實際 TOF 差。由于一體化 TOF-PET/MR 使用 TOF 差推測湮滅發生位置進行重建，所以晶體時間整體偏移不會對重建圖像產生影響。
圖2
三種方法計算的時間偏移與互相比較結果
a. 扇形束法按每個晶體顯示的計算結果；b. L1-norm 法按每個晶體顯示的計算結果；c. L2-norm 法按每個晶體顯示的計算結果；d. 扇形束法結果減去 L2-norm 法結果的分布直方圖，為類高斯分布；e. 扇形束法結果減去 L1-norm 法結果的分布直方圖，為類高斯分布；f. L1-norm 法結果減去 L2-norm 法結果的分布直方圖，為類高斯分布
Figure2. Timing offsets and comparison results calculated by three methodsa. calculation results of fan-beam method；b. calculation results of L1-norm method；c. calculation results of L2-norm method；d. distribution histogram of the calculation results of fan-beam subtracting L2-norm；e. distribution histogram of the calculation results of fan-beam subtracting L1-norm；f. distribution histogram of the calculation results of L1-norm subtracting L2-norm
2.3 位置校正準確性分析
為了確定位置校正算法準確性，我們使用基于有序子集最大期望值法（ordered subset expectation maximum，OSEM）的 TOF 重建獲得高清圖像[15]。將得到的圖像進行閾值分割得到筒源輪廓，在輪廓內部用兩條線逼近的方法，找到對應切面的圓心坐標，作為參考值。利用一個肉眼幾乎不可見的偏心筒源進行位置校正準確性驗證，結果見表 1。
表1
位置校正準確性驗證結果
Table1.
Accuracy of position correction
由于重建圖像在筒源首尾存在邊界模糊，導致首尾偏差比較大，我們取中間層進行分析，可知軸心位置計算偏差百分比在 15% 以內。軸心在 x 坐標偏心 1 cm 情況下，誤差在 1.5 mm 內。對于 TOF 校正結果的偏差為 10 ps，這個誤差是可以接受的。由于在實際采集中無法定量描述可能發生的肉眼擺放誤差，我們在實驗中采集了一組偏心較大筒源數據，偏移范圍在 4 cm 左右，在未進行位置校正情況下，對系統進行時間校正。得到校正結果后對該偏心筒源進行重建，圖像產生很大偽影。再使用本文提出的位置校正算法后進行時間校正，同樣重建條件下，圖像偽影有很大改善。結果表明我們所使用的位置校正算法，能減小因為位置偏差對時間校正的誤差。
2.4 NEMA IQ 實驗圖像質量分析
NEMA 定義了圖像質量的測試標準，該標準給出了臨床成像條件下圖像質量的分析方法。對比度復原系數（contrast recovery coefficient，CRC）表示 NEMA IQ 中球體在背景下的對比度還原情況，通過計算 NEMA IQ 模體內 4 個熱球和 2 個冷球的 CRC，可以反映 PET 圖像質量。
我們對同一個 NEMA IQ 模體采集的符合數據，對未進行時間校正以及分別采用 L2-norm 法、扇形束法和 L1-norm 法時間校正后的數據進行重建。重建使用前述 OSEM 算法，參數設置相同的情況下，取第三次迭代圖像結果的同一切面進行分析。同一機器同樣條件重復進行三次實驗，取平均值作統計圖，誤差條為正負標準差，結果如圖 3 所示。圖中黑色圓形表示熱球，半徑從 10 mm 到 22 mm 逆時針增加；白色圓形表示冷球，半徑為 28 mm 和 37 mm，逆時針增加；中間圓形為空氣，在圖表中記為 Lung。三種時間校正方法所得圖像均清晰無偽影，可以清楚看到熱球和冷球位置，與未校正結果對比明顯。定量分析圖 3 中各方法圖像的 CRC 值，未時間校正圖像的 CRC 與三種不同方法時間校正圖像 CRC 有明顯差異，同時噪聲也更大。扇形束法噪聲平均值為 6.85%，略大于 L1-norm 法的 5.92% 和 L2-norm 法的 5.49%。如圖 3 中右邊圖表可知，扇形束法的熱球 CRC 平均值 68.62% 略高于 L1-norm 法的熱球 CRC 平均值 67.58% 和 L2-norm 法的熱球 CRC 平均值 67.28%。這表明在半徑較小的熱球區域，扇形束法圖像質量較好，顯示了對小病灶的探查潛力。
圖3
NEMA IQ 模體重建圖像以及 CRC 值與噪聲值統計結果
Figure3.
Reconstructed images of NEMA IQ phantom and its statistical results of CRC and noise
2.5 腫瘤患者全身 PET 圖像的比較
隨機選擇了 10 例復旦大學附屬中山醫院 2018 年接受一體化 TOF-PET/MR 掃描的腫瘤患者，根據復旦大學附屬中山醫院和上海聯影醫療科技有限公司臨床合作協議，所有數據使用均已獲授權。用不同方法時間校正后的系統分別重建患者原始數據得到全身 PET 圖像，采用 PET 圖像中腫瘤區域基于去脂體重的標準攝取值（standard uptake value based on lean body mass，SUV lbm）比較結果。在同一正切面圖像中圈出感興趣區域（region of interest，ROI），比較 ROI 中 SUV 的平均值。SUV 在圖像中表現為體素的信號強度值，與病灶濃度 C病灶（kBq/mL）和注射劑量 D注射（MBq/kg）有關[16]，具體關系為：
|
圖 4a 為未校正圖像，其他三幅分別為使用不同時間校正方法的晶體時間偏移數據后重建獲得的患者圖像。從圖 4 分析可知，經過時間校正后的重建圖像，對湮滅發生的位置定位更加準確，可以清晰地看到腫瘤區域，腫瘤區域的 SUV 平均值增大。
圖4
患者腫瘤圖像 SUV 對比
a. 未校正時患者全身冠狀面圖像；b. L2-norm 法校正后患者全身冠狀面圖像以及 SUV 值；c. 扇形束法校正后患者全身冠狀面圖像以及 SUV 值；d. L1-norm 法校正后患者全身冠狀面圖像以及 SUV 值
Figure4. Comparison of SUV calculated from patient’s tumor imagesa. uncalibrated whole-body coronal image；b. whole-body coronal image with L2-norm calibration；c. whole-body coronal image with fan-beam calibration；d. whole-body coronal image with L1-norm calibration
將全部 10 例患者腫瘤區域 SUV 平均值以散點圖形式表示，并對每種方法的平均值進行直方圖表示，誤差條是十個患者數據的標準差，如圖 5 所示。L2-norm 法由于使用較多數據，更加真實地反映了時間偏差，使得腫瘤區域 SUV 平均值稍高。利用單因素方差分析四組結果發現總體差異有統計學意義，其中未校正圖像和三種不同方法時間校正后圖像均有顯著差異，而三種校正方法兩兩進行方差分析其差異無統計學意義。檢驗結果表明，未校正圖像 SUV 平均值與經過時間校正后圖像 SUV 平均值差異顯著，而三種校正方法之間 SUV 平均值差異不明顯。
圖5
十位患者腫瘤圖像 SUV 平均值統計圖
直方圖是不同方法校正時所有患者 SUV 的平均值，散點圖為每個患者的具體腫瘤區域 SUV 平均值。*：
histogram indicates average value of ten patients’ SUV with different calibration，and scatter plot indicates specific SUV mean value of every patient. *：
3 討論與結論
TOF-PET 自 20 世紀 80 年代被提出，隨著新型閃爍晶體的應用，時間分辨率越來越好，TOF-PET 重建技術也隨之發展。TOF 系統在探測過程中能定位湮滅發生的區域，提高反投影的空間精度[17]。TOF 技術最大程度地發揮了摻鈰硅酸镥（lutetium oxyorthosilicate，LSO）閃爍晶體在 PET 中的表現，在減少隨機事件和系統噪聲以及加快掃描時間上有重要的價值[18]。對于一體化 TOF-PET/MR 系統，由于需要更加頻繁的時間校正，算法的優劣將對系統 TOF 表現造成很大影響。
我們在一體化 TOF-PET/MR 系統中利用扇形束法、L2-norm 法和 L1-norm 法進行時間校正，并考慮了采集時筒源擺放偏心情況的位置校正。GE 醫療在 2014 年研發出帶有 TOF 功能的一體化 PET/MR，其中使用了磁場條件下電子性能更好的硅光電二極管，其時間分辨率達到 390 ps。Grant 等[19]對該系統測試表明，TOF 能有效增加重建圖像的信噪比，具有更好圖像質量、更短掃描時間和更少注射劑量的潛力。Minamimoto 等[20]對一體化 TOF-PET/MR 的臨床表現進行了評估，擁有良好時間分辨率的 PET/MR 與無 TOF 功能 PET/MR、TOF-PET/CT 以及無 TOF 功能 PET/CT 相比，減少了圖像偽影，提高了 PET 臨床圖像質量。隨著 TOF 技術的發展，高時間分辨率的一體化 TOF-PET/MR 將會成為主流發展方向，因此能夠發揮硬件極限時間分辨率的時間校正算法愈發重要。本文比較了三種不同方法的效果，證明了時間校正的必要性，對時間校正方法的選取具有重要意義。
自 Luo 等[9]提出時間校正問題，已經有多種時間校正的方法被提出。主流方法之一是 Lenox 等[8]提出的扇形束法，通過計算每個晶體的時間偏差直方圖，對于圓柱狀模體，形成類似扇形面的 LOR 集合。另一主流方法是構建超定方程組解每個晶體的時間偏差，由 Reynolds 等[13]最先提出，將問題轉化為最小二乘問題，并利用 L2-norm 最小化的方法求解。后來，Freese 等[10]基于 L2-norm 最小化法，提出用 L1-norm 替代，計算得到魯棒性好的結果。然而，目前并沒有針對 TOF-PET 的各種主流時間校正方法效果的對比。本文對比了三種時間校正方法的整體流程所需時間和圖像質量，分析了各種方法的優劣，對一體化 TOF-PET/MR 時間校正方法提出了建議，且本文方法也可以用于其他 TOF-PET 系統。但本文所提三種方法，仍存在散射和隨機事件對校正結果影響較大的問題，在未來的研究中，應該去探究如何減少數據采集中的散射和隨機事件，或者研究新的時間校正算法，完全消除散射和隨機產生的影響。本文采取的位置校正方法不僅可以減少位置擺放偏差帶來的圖像偽影，而且能根據計算出的筒源位置，判斷筒源擺放是否準確，如果偏差非常大，產生錯誤，可以重新擺位進行采集。但是位置校正算法在部分切面上誤差大于 10%，這樣的誤差也使得位置校正無法完全消除偏心帶來的圖像偽影，如何提高該方法的準確性，獲取更真實的放射源位置，需要進一步探索。
綜上所述，扇形束法和 L1-norm 法與 L2-norm 法相比，所用數據量大幅減少，計算速度提升，使商品化 TOF-PET 周期性時間校正的工作量大幅降低，尤其對于需要更頻繁進行時間校正的一體化 TOF-PET/MR 來說具有很大優勢。扇形束法與其他方法相比有更佳的圖像質量，尤其是在較小病灶的探查上表現出更大潛力，因此綜合考慮時間成本和圖像質量，一體化 TOF-PET/MR 最適合選用扇形束法進行時間校正。本文在時間校正過程中提出了位置校正算法，對數據進行預處理，能發現肉眼無法分辨的擺放偏差并有效改善筒源擺放誤差對時間校正結果的影響。位置校正算法結合扇形束法將會使系統的時間分辨率提升，可廣泛應用在商品化和研究用途的一體化 TOF-PET/MR 中。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
