摘 要：低劑量輻射是CT成像中的研究熱點和難點，而不完全角度投影又是低劑量輻射中的一種有效手段。因此對不完全角度投影的重建問題也是一個研究熱點。建立在CS理論基礎上的TV-EM算法雖然能夠克服傳統算法的重建精度不高，去噪能力不明顯等特點，但也對邊緣進行了過度平滑。本文在懲罰正則項中加入中值先驗分布，在抑制噪聲的同時更能保留邊緣信息。加入中值先驗，形成新的TV-EM算法，并將其應用于不完全投影下的重建。通過仿真實驗，選擇合適的懲罰因子，驗證了其超強的抗噪能力和穩健的收斂性。

　　關鍵詞：不完全投影重建;稀疏懲罰;貝葉斯方法

　　1 引言

　　隨著科技的發展，CT在臨床診斷中的作用日顯重要。而人們對健康的重視程度也越來越高。因此CT在醫學診斷中的掃描所產生的射線輻射也受到人們的廣泛關注。畢竟，射線輻射對人體健康是有損害的。據報道，接受超過28次CT掃描的患者致癌幾率比平均水平高12%，而兒童在接受CT檢查時所受的輻射影響會更大。因此，人們都希望接受足夠小的輻射，而能得到精確的CT圖像輔助診斷。因此，低劑量輻射也成為醫學影像中亟待解決的熱點和難點。而有限角度和稀疏角度重建就是降低輻射劑量的一種有效手段。

　　在不完全角度投影數據條件下，無論是有限角度還是稀疏角度，經典FBP算法的重建結果總是不夠理想，雖然重建速度快，但是存在偽影，甚至有一部分信息嚴重缺失，導致重建圖像質量不盡如人意[1]。這就推動了迭代重建算法的研究。2006年，Candes等[2]提出了壓縮感知(CS)理論。隨后做了大量重建算法的研究，如Sidky等人提出全變分最小化算法[3，4]，Chen等[5]于2008年提出PICCS算法，而Wang等[6]于2010年提出RRD算法，Zhang等[7]于2013年提出ADTVM算法，都是基于CS理論，選取一種或幾種稀疏變換作為目標函數，再以成像模型AX=P作為約束條件，再求解出重建圖像[8]。在實際成像過程中，由于數據的隨機噪聲和成像模型的非隨機誤差，造成成像模型AX=P總是不相容的。所以用成像模型AX=P作為約束條件并不總是最合適的。本文研究的TV-EM(全變差-期望最大)算法[9]，是用貝葉斯方法建立目標函數，以ML-EM算法作為基礎，以圖像X的稀疏變換提取圖像的先驗知識，加入圖像的中值先驗到每一步迭代中，從而改善重建圖像的質量[10，11]。

　　2 獲取投影

　　重建物體的離散模型以網格像素(二維圖像)的形式進行刻畫，獲取投影數據模型如圖1所示[12]。探測器檢測結果為射線方向的累積值，并假設射線沒有寬度。探測器檢測的投影值與射線穿過體素的長度相關。實際上投影值是用像素內的線段長度加權的像素值的和。

　　這里，像素值xj為第j個像素的灰度值。探測器的第i個探測元發出的射線在每個像素內的線段長度記為aij。如果第i條射線與第j個像素不相交，則aij=0。而第i條射線的投影值pi是探測到的光子數。而投影方程可表示為：

　　pi=aijxj

　　3 貝葉斯圖像重建原理

　　對于一個成像問題來說，假定每一個圖像像素發出的光子數是一個獨立的泊松隨機變量。而實際上第i個探測元理論上探測到像素j發出的光子數的均值是aijxj。而所有的探測元探測到的pi組成一個樣本。那么可建立似然函數[17]：

　　prob=∏i(∑jaijxj) (2-1)

　　取對數似然

　　ln(prob)=∑i[-∑jaijxj+piln(∑jaijxj)-lnpi!] (2-2)

　　在圖像重建中，ML-EM方法就是尋找一個解，使得探測到的數據最有可能發生，即最大化概率Prob，其中pi為探測到的投影數據，而xj是待求圖像。可用如下優化目標來表示：

　　=argmaxx≥0ln(prob) (2-3)

　　為求解目標(2-3)，Shepp和Vardi提出ML-EM算法，迭代公式[16]為：

　　j=∑i (2-4)

　　其中pi中是投影數據，xj是圖像當前估計值。

　　4 TV-EM算法

　　TV約束是基于Donho，Candes和Tao等提出的壓縮感知理論，假定圖像具有分段光滑性質。選擇圖像具有稀疏性的TV范數作為正則項，加入(2-3)作為正則約束。

　　根據Candes和Donoho提出的理論，圖像重建目標變為：

　　=argmaxx≥0ln(prob)+?琢TV(x) (3-1)

　　?琢——正則化參數

　　TV(x)——圖像的TV范數

　　TV(全變差)范數是圖像的梯度的l1范數，但由于使用起來不方便，習慣用l2范數近似代替l1范數，這樣就可以對TV范數求導了。但是l2范數在抑制噪聲的同時，把邊緣信息也過度平滑了。引入TV范數為稀疏懲罰項的TV-EM算法迭代公式[9]為：

　　j=∑i (3-2)

　　其中，TV范數求偏導是離散化后的公式[12]：

　　=

　　+

　　- (3-3)

　　參數?著是一個很小的數，為了避免分母為0。

　　5 中值先驗

　　中值先驗在降低噪聲的同時，也保持了圖像的邊緣，它的分布如下[18]：

　　R(x)=aexp-∑j

　　a——常數;

　　?茁——正則化參數;

　　Mj——領域中像素的中間值。

　　推薦閱讀：低溫建筑技術發表論文要求