帶有松弛因子和懲罰項(xiàng)的最小二乘PET圖像重建

劉 祎 桂志國

（中北大學(xué)電子測(cè)試技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，太原 030051）

引言

正電子放射斷層成像（positron emission tomography，PET）技術(shù)是一種放射性核素顯像技術(shù)，是核醫(yī)學(xué)成像技術(shù)和斷層成像技術(shù)的最新發(fā)展。PET的最大優(yōu)勢(shì)是能定量評(píng)價(jià)人體組織的生理、生化功能，被譽(yù)為活體的分子斷層圖像［1］，也使得PET成為核醫(yī)學(xué)診斷中的重要技術(shù)手段。

PET重建算法分為解析法和迭代法，經(jīng)典解析法是濾波反投影算法（filtered back projection，F(xiàn)BP）［2］。FBP法簡(jiǎn)單，重建速度快，但FBP算法過分簡(jiǎn)化了PET的物理模型，要求有完整的投影數(shù)據(jù)，無法有效地抑制噪聲［3］。而迭代重建算法適用于信息量相對(duì)不足的情況。經(jīng)典的迭代算法是基于泊松模型的最大似然算法（maximum likelihood，ML）和基于高斯模型的加權(quán)最小二乘算法（weighted least squares，WLS）。但這兩種算法都存在著由于噪聲而導(dǎo)致的圖像退化和收斂速度慢的問題，而且隨著迭代次數(shù)的增加，圖像質(zhì)量下降，噪聲非常嚴(yán)重［5］。為了加速迭代收斂，Hudson和Larkin提出了有序子集（ordered subset，OS）方法，應(yīng)用于MLEM算法，取得了較好的圖像效果。自此，OS算法不僅用于EM算法，更與其他的算法結(jié)合派生出新的“OS”型算法［6］，最典型的是OSLS算法。其他的加速收斂速度的算法還有Fessler提出的廣義的空間更新期望值最大化（space alternating generalized expectation maximization algorithm，SAGE）算法［7-8］。雖然這些加速方法可以加快算法的收斂速度，但依然具有較差的噪聲特性。為了克服噪聲的影響，許多學(xué)者針對(duì)算法引入了約束條件和懲罰函數(shù)，以改善圖像質(zhì)量。

為解決噪聲和收斂問題，本研究帶有松弛因子和懲罰項(xiàng)的最小二乘算法。首先從最小二乘目標(biāo)函數(shù)出發(fā)，為平滑圖像和抑制噪聲加入二次懲罰項(xiàng)，其次，為減少運(yùn)算量而采用修改的逐次超松弛（successive over-relaxation，SOR）迭代法求解最優(yōu)解［9］，并利用可變有序子集（MOS）對(duì)其進(jìn)行加速收斂，獲得基于可變有序子集的懲罰最小二乘重建算法，文中記為MOS-PLS+MSOR。仿真實(shí)驗(yàn)中，將MOS-PLS+MSOR算法與最小二乘（least squares，LS）、懲罰最小二乘（penalized least squares，PLS）等算法的重建結(jié)果作比較，實(shí)驗(yàn)證明，MOS-PLS+MSOR有更好的重建效果。

1 算法理論

在PET重建中，常采用的離散模型形式為

其中P=［p1，p2，…pN］T表示觀測(cè)的投影數(shù)據(jù)，G為N×M維的系統(tǒng)矩陣，λ=［λ1，λ2，…λM］T表示圖像向量。gij表示第j個(gè)像素被第i對(duì)探測(cè)器檢測(cè)到的概率。

基于最小二乘估計(jì)的目標(biāo)函數(shù)為

寫成矩陣的形式為

因此，最小二乘估計(jì)為：

顯然，對(duì)式（1）的求解相當(dāng)于解下面的方程組：

由實(shí)際的系統(tǒng)模型可知G是大型的稀疏矩陣，直接求解相當(dāng)困難，可以通過共軛梯度法等方法求解，這里使用SOR迭代算法［10］，因?yàn)橥瑫r(shí)SOR算法高頻分量收斂速度比低頻分量快，本質(zhì)上也加快了整體的收斂速度［11］。而且OS方法雖然能加快重建速度，但OS重建的圖像往往不能收斂，隨著迭代次數(shù)的增加，噪聲大大增加。解決OS發(fā)散的一個(gè)方法就是使用松弛參數(shù)［12］。

記A=GTG，B=GTP，則式（5）可以寫成

很顯然，A=（aij）∈RM×M是對(duì)稱的非奇異矩陣，且滿足aii≠0，i=1，2，…，N，由SOR迭代公式得：

式中，ω為松弛因子，取值范圍為0＜ω＜2。將A=GTG、B=GTP帶入式（7）得：

由于系數(shù)矩陣的特殊性，對(duì)于較小的圖像，上式中的計(jì)算量已經(jīng)非常大，不僅耗費(fèi)的時(shí)間長(zhǎng)，而且還要占用很大的內(nèi)存空間。為了計(jì)算簡(jiǎn)便，本文考慮雅可比迭代的思想，將式（8）中的換成，即把式（8）中的后兩項(xiàng)合并起來，這樣就不用計(jì)算完后再計(jì)算，從而大大減少了迭代時(shí)間。從而得到簡(jiǎn)便的并行迭代公式，寫成矩陣的形式：

很顯然，式（9）屬于同時(shí)迭代重建法，有利于克服圖像對(duì)噪聲的敏感性。PET圖像重建實(shí)質(zhì)是一個(gè)病態(tài)解的求解問題，式（9）的方法雖然有利于克服圖像對(duì)噪聲的敏感性，但仍然不能克服求解病態(tài)解得困難。正則化是求解病態(tài)問題的常用方法。因此，本文中引入正則項(xiàng)（懲罰項(xiàng)），許多學(xué)者已經(jīng)提出很多懲罰函數(shù)［13］，一些懲罰函數(shù)目的在于保持邊緣銳化的同時(shí)平滑灰度一致的區(qū)域。本文引入簡(jiǎn)單的二次平滑懲罰項(xiàng)［14］，如下所示：

其中Sj表示第j個(gè)像素的8鄰近像素構(gòu)成的像素集合，ωjk表示權(quán)重量，水平方向和垂直方向上ωjk=1，對(duì)角線方向上ωjk=1/。即ωjk的值同像素b和j之間的距離成反比。不考慮邊緣效果，優(yōu)化矩陣R可以表示為：

此時(shí)目標(biāo)函數(shù)為

其中β是平滑參數(shù)，考慮二次懲罰項(xiàng)，式（9）的迭代公式變?yōu)?/p>

文中將這種方法記為PLS+MSOR方法。

2 基于可變有序子集的PLS算法

MSOR算法雖然減少了迭代過程中的計(jì)算量，使算法較之于傳統(tǒng)的SOR使用了較少的時(shí)間，但并不能改變收斂速度。OS方法作為加速收斂速度的重要方法，已經(jīng)得到廣泛的應(yīng)用。在理論上，隨著子集劃分個(gè)數(shù)增加，收斂速度增加，高頻分量得到恢復(fù)。但并不是子集個(gè)數(shù)越多越好，因?yàn)橐话阋裱白蛹胶狻钡臈l件，即每一個(gè)子集都含有相等的圖像放射性參數(shù)信息，子集個(gè)數(shù)增加使噪聲隨迭代次數(shù)的增加逐步放大［15-17］。然而，若采用較小的子集個(gè)數(shù)，圖像的高頻信息就難以恢復(fù)，并且收斂速度也會(huì)變慢。

為了解決這個(gè)問題，采用可變有序子集的方法（MOS），先選擇大的子集數(shù)，以恢復(fù)高頻信息，再適當(dāng)減小子集數(shù)，減輕噪聲對(duì)重建圖像的影響，然后再選擇稍大的子集數(shù)，以恢復(fù)沒有恢復(fù)的高頻信息，如此循環(huán)往復(fù)，即加快了收斂速度，同時(shí)也減少了噪聲的影響。MOS算法在迭代過程中改變子集的大小，克服了OS算法只采用固定子集造成的不能同時(shí)顧全噪聲和高頻信息的缺點(diǎn)。

3 實(shí)驗(yàn)方法

為驗(yàn)證MOS-PLS+MSOR算法的有效性，采用大小為128像素×128像素的Sheep-Logan頭部模型（λ*）進(jìn)行仿真，如圖1所示。圖像像素點(diǎn)取值范圍為0到8，模擬的投影參數(shù)為N=128×128，即在0°～180°內(nèi)均勻分布128個(gè)投影方向，每個(gè)方向有128個(gè)探測(cè)器對(duì)。系統(tǒng)矩陣G不考慮探測(cè)器的效率和衰減校正等因素，利用公式P*=G·λ*產(chǎn)生無噪聲數(shù)據(jù)，并以P*為泊松變量噪聲均值來生成所需要的觀測(cè)數(shù)據(jù)P。為了公平比較，設(shè)置所有的算法均使用相同的均勻正值起始圖像，初始值設(shè)為0.5，實(shí)驗(yàn)步驟如下：

步驟1：置k=0，給未知量λk賦初值；

步驟2：在t（t=1，...，T）個(gè)子集內(nèi)順序執(zhí)行下面的迭代更新公式

步驟3：T個(gè)子集都更新，即完成一次迭代；

步驟4：設(shè)置迭代停止條件，滿足條件則停止迭代，不滿足條件則將更新后的像素值λk再重復(fù)執(zhí)行步驟2。

圖1 Sheep-Logan模型Fig.1 Sheep-Logan model

4 實(shí)驗(yàn)仿真

4.1 圖像質(zhì)量評(píng)價(jià)參數(shù)

為評(píng)價(jià)重建圖像質(zhì)量，需對(duì)重建結(jié)果進(jìn)行定量描述。采用一些常見的質(zhì)量評(píng)價(jià)參數(shù)定量分析圖像的重建質(zhì)量，分別是相關(guān)系數(shù)、信噪比、歸一化均方誤差，定義如下：

1）相關(guān)系數(shù)（correlation coefficient，CORR）

2）信噪比（signal to noise ratio，SNR）

信噪比越高表示重建圖像質(zhì)量越好。

3）歸一化均方誤差（root mean squared error，RMSE）

表征重建圖像與原始圖像的近似程度，RMSE越小則表明重建結(jié)果越接近原始圖像。

4.2 幾種重建算法的比較

將本研究的PLS+MSOR算法與LS，PLS算法做比較，圖2分別給出了計(jì)量數(shù)為1×106時(shí)3種算法重建的結(jié)果。LS迭代21次，PLS迭代27次，PLS+MSOR迭代35次（ω=0.25β=0.002 12）。從圖2的結(jié)果可以看出，在視覺效果上，3種算法中PLS+MSOR算法重建的圖像較之于LS算法和PLS有較少的偽影，而且中間區(qū)域比較平滑，重建質(zhì)量最好。同時(shí)，圖3給出了3種算法的CORR曲線變化和SNR曲線變化。

表1中列出了各種算法的質(zhì)量參數(shù)和迭代時(shí)間。從表中可以看出，由于PLS中含有懲罰項(xiàng)，因此迭代時(shí)間比LS的迭代時(shí)間長(zhǎng)。但是PLS+SOR的迭代時(shí)間卻比PLS的迭代時(shí)間短，甚至比LS的迭代時(shí)間還要短，足以看出PLS+MSOR算法在時(shí)間上的絕對(duì)優(yōu)勢(shì)。

圖2 LS，PLS，PLS+MSOR算法重建圖像。（a）LS算法；（b）PLS算法；（c）PLS+MSOR算法（計(jì)量數(shù)1×106）Fig.2 Reconstruction images of LS，PLS，PLS+MSOR.（a）LS；（b）PLS；（c）PLS+MSOR

圖3 LS，PLS，PLS+MSOR算法的質(zhì)量參數(shù)曲線。（a）CORR曲線；（b）SNR曲線Fig.3 Quality parameter curve of LS，PLS，PLS+MSOR.（a）CORR；（b）SNR

然而在實(shí)際的PET檢測(cè)中，光子數(shù)量與攝入體內(nèi)的放射核素的濃度有關(guān)，光子計(jì)數(shù)值不可能很大，如果要得到較大的計(jì)數(shù)值，就需要注射過量的放射性核素，這將給患者造成傷害。另外可獲得的最大的計(jì)數(shù)值受探測(cè)器晶體的計(jì)數(shù)能力和后項(xiàng)處理電路的限制，因此在圖4中給出了低計(jì)數(shù)情況下的重建結(jié)果。從圖2和圖3對(duì)比可以看出，低計(jì)數(shù)量情況下，LS和的PET圖像重建質(zhì)量遠(yuǎn)遠(yuǎn)不如高計(jì)數(shù)量情況下的質(zhì)量，但是PLS+MSOR算法在低劑量情況下仍然有較好的重建效果。在表2中列出了計(jì)數(shù)量在1.26×105情況下3種算法重建圖像的質(zhì)量參數(shù)，可以看出，PLS+MSOR算法在低級(jí)數(shù)量時(shí)的重建效果要遠(yuǎn)好于LS算法和PLS算法。

圖4 LS，PLS，PLS+MSOR算法重建圖像。（a）LS算法；（b）PLS算法；（c）PLS+MSOR算法（計(jì)量數(shù)為1.26×105）Fig.4 Reconstruction images of LS，PLS，PLS+MSOR.（a）LS；（b）PLS；（c）PLS+MSOR

MOS-PLS+MSOR的子集個(gè)數(shù)序列為16，4，4，8，4，4，OS-LS和OS-PLS的子集數(shù)為4，并使相鄰子集滿足正交關(guān)系，重建結(jié)果如圖5所示。表3中列出了3種算法的質(zhì)量參數(shù)。相對(duì)于OS-LS和OS-PLS算法，本算法在有限的迭代中表現(xiàn)出較好的噪聲特性，說明了OS-PLS不僅能夠重建出質(zhì)量較好的圖像，而且具有收斂速度上的優(yōu)勢(shì)。為了更能說明本算法的可靠性，圖6給出了3種算法的RMSE曲線和SNR曲線。

圖5 OS-LS，OS-PLS和MOS-PLS+MSOR算法的重建圖像。（a）OSLS算法；（b）OS-PLS算法；（c）MOS-PLS+MSOR算法；Fig.5 Reconstruction images of OS-LS，OS-PLS and MOS-PLS+MSOR.（a）OS-LS；（b）OS-PLS；（c）MOS-PLS+MSOR

表1 LS，PLS和PLS+SOR的質(zhì)量參數(shù)（計(jì)量數(shù)為1×106）Tab.1 The quality parameter of LS PLS and PLS+SOR algorithm（The number of measurement：1×106）

表2 LS，PLS和PLS+SOR的質(zhì)量參數(shù)（計(jì)量數(shù)為1.26×105）Tab.2 The quality parameter of LS PLS and PLS+SOR algorithm（The number of measurement：1.26×105）

圖6 OS-LS，OS-PLS，MOS-PLS+MSOR算法的質(zhì)量參數(shù)曲線。（a）RMSE曲線；（b）SNR曲線Fig.6 Quality parameter curve of OS-LS，OSPLS，MOS-PLS+MSOR.（a）curveof RMSE；（b）curve of SNR

表3 OS-LS，OS-PLS和MOS-PLS+MSOR的質(zhì)量參數(shù)（計(jì)量數(shù)為1×106）Tab.3 The quality parameter of OS-LS，OS-PLS and MOS-PLS+SOR algorithm（The number of measurement：1×106）

4 結(jié)論與討論

MOS-PLS+MSOR算法具有加速收斂和抑制噪聲的作用，較之于傳統(tǒng)的最小二乘算法有較好的重建效果。然而本研究中使用的二次先驗(yàn)懲罰項(xiàng)雖然比較簡(jiǎn)單且在理論上更加易于分析，但是對(duì)圖像邊緣區(qū)域造成過度平滑。同時(shí)因此使用更為合理的先驗(yàn)懲罰項(xiàng)，這也是進(jìn)行下一步研究的內(nèi)容。另外需要研究關(guān)于迭代過程中如何選取松弛因子的問題，但至今仍未能有效地解決，這也是下一步要研究的內(nèi)容。

［1］張澤寶，醫(yī)學(xué)影像物理學(xué)［M］.北京：人民衛(wèi)生出版社，2000：189-191.

［2］Nam-Yong Lee，Yong Choi.A modified OSEM algorithm for PET reconstruction using wavelet processing［J］.Computer Methods and Programs in Biomedicine，2005，80：236-245.

［3］趙凌云.三維正電子發(fā)射斷層成像理論及實(shí)現(xiàn)［D］.浙江：浙江大學(xué)，2001.

［4］Shepp SA，Vardi Y.Maximum likelihood restoration for emission tomography［J］.IEEE Trans Med Imageing，1982，1：113-122.

［5］Fessler J.A.，Hakan E.A paraboloidal surrogates algorithm for convergent penalized-likelihood emission image reconstruction［C］//Fessler J.A..IEEE Nuclear Science Symposium and Medical Imaging Conference Record.Toronto：IEEE，1998：1132-1135.

［6］周健.正電子發(fā)射斷層的圖像重建方法研究［D］.南京：東南大學(xué)，2006.

［7］Fessler JA，Hero AQ.Space-alternating generalized expectation maximization algorithm［J］.IEEE Trans on Signal Processing，1994，42（10）：2664-2676.

［8］Fessler JA，Hero AQ.Penalized maximum-likelihoodimage reconstruction using space-alternating generalized EM algorithms［J］.IEEE Trans Imag Proc，1995，4（10）：1417-1429.

［9］孔慧華，潘晉孝.序列子集聯(lián)合代數(shù)重建技術(shù)［J］.CT理論與應(yīng)用研究，2008，17（2）：40-45.

［10］封建胡，車剛明.數(shù)值分析原理［M］.北京：科學(xué)出版社，2001.

［11］Fessler JA.Penalized weighted least-squares image reconstruction for positron emission tomography［J］.IEEE Medical Imaging，1994，13（2）：290-300.

［12］顏建華.正電子發(fā)射斷層圖像重建算法研究［D］.武漢：華中科技大學(xué)，2007.

［13］Geman S，McClure D.E.，Statistical methods for tomographic image reconstruction［C］//Geman S.Processings 46th Session.ISI，Bulletin of the ISI.Tokyo：Legal Information Institute，1987：5-21.

［14］Sangtae A，F(xiàn)essler JA.Globally convergent ordered subsets algorithms：Application to tomography［J］.IEEE Transactions on Medical Imaging，2003，22（5）：613-626.

［15］張書文，陳英茂，田嘉禾.OSEM子集及迭代次數(shù)對(duì)PET圖像質(zhì)量的影響［J］.中國醫(yī)學(xué)影像學(xué)雜志，2003，11（3）：199-201.

［16］Liang Z，Jaszczak R，Creer K.On Bayesian image reconstruction from projections：uniform and nonuniform a priorisource information［J］.IEEE Transactions on Medical Imaging，1989，8（3）：227-235.

［17］Hudson HM，Larkin RS.Accelerated image reconstruction using ordered subsets of projection data［J］.IEEE Trans Medical Imaging，1994，13（4）：601-609.