基于最大似然和罰似然估計的基于最大似然和罰似然估計的 CT 統計重建算法研究統計重建算法研究 摘 要 當 CT 系統在數據采集過程中受到嚴重的噪聲干擾，或者投影數據采集不完全時， 用解析重建算法重建出的圖像有偽跡。而統計重建算法具有物理模型準確、對噪聲不 敏感、易于加入約束條件等優點，其重建的圖像質量要優于傳統的 FBP 方法。針對 CT 統計重建，本文主要研究了以下內容： 研究了基于最大似然估計的統計重建算法：期望最大值法(Maximum Likelihood Expectation Maximization，ML-EM)，可分離拋物面型替代(Separable Paraboloidal Surrogate，ML-SPS)算法和這兩種算法的有序子集形式（OS-ML-EM 和 OS-ML-SPS） 。 仿真實驗結果表明 ML-SPS 的初始收斂速度比 ML-EM 快，但該兩種算法的收斂速度 遠慢于其對應的有序子集形式。用 OS-ML-EM 算法和 OS-ML-SPS 算法對 CT 采集的 實際投影數據進行重建，得到了較好的重建結果。 研究了基于罰似然(Penalized Likelihood，PL)的統計重建算法:OSL(One Step Late)- EM 算法和 PL-SPS 算法。重點討論了基于 Gibbs 分布的罰函數，從理論上分析了勢函 數需要滿足的條件以及勢函數的選取對圖像的影響，著重分析了二次勢函數和 Huber 勢函數的優缺點，通過仿真實驗對它們的重建結果進行誤差分析。 提出了罰似然重建中正則項參數的自適應選取的方法，該方法充分利用每一次迭 代的重建結果信息，不斷對正則項參數進行更新。并用于仿真實驗和重建 CT 實際投影 數據，重建結果表明該方法降低了噪聲，提高了圖像質量。 將 OS-OSL-EM 算法推廣到三維錐束軌跡下的圖像重建，取得了較好的仿真結果。 關鍵詞：最大似然估計，CT 重建，罰似然估計，有序子集 The Study of CT Statistical Reconstruction Algorithm Based on Maximum Likelihood and Penalized Likelihood Estimates Abstract When the CT data has serious noise in the process of acquisition system and projection data is incomplete , analytic reconstruction algorithm get images with artifact. The statistical reconstruction algorithm with a accurately physical model isnt sensitive to noise and easy to add constraints etc.Therefore, the reconstructed image quality is superior to conventional FBP methods. For the CT Statistics reconstruction, This paper mainly studies the following content: The maximum likelihood estimate for the statistical reconstruction algorithms is studied and mainly includes expectation maximum algorithm (ML-EM),separable paraboloidal surrogate(ML-SPS) algorithm and the both algorithms with ordered subset(OS-ML-EM and OS-ML-SPS). In simulation experiments, it shows that the initial convergence rate of ML- SPS is faster than the ML-EM ,however, the both algorithms is slower than OS-ML-EM and OS-ML-SPS algorithms. At last, OS-ML-SPS and OS-ML-EM algorithms are used to reconstruct actual CT projection data and get better image. Penalized maximum likelihood(PL) reconstruction for the statistical algorithms is based on the penalized likelihood estimates and adds a penalty term to suppress noise, which mainly includes OSL-EM algorithm and PL-SPS algorithm. It is focused on the penalty function which is based on Gibbs distribution and analyzes the need of potential functions to meet the conditions, the influence on image for the selection of potential functions and the advantages and disadvantages of the quadratic function and the Huber function. At last, simulation experiments give the error analysis of reconstructed image. This method which is an adaptively regularied CT image reconstruction can adaptively choose regularization parameters with making full use of the results of every iteration to update regularization parameters. Simulation results show that this method reduces the noise and improve image quality. The method is used to reconstruct the actual CT projection data to demonstrate the feasibility and practicality. OS-OSL-EM algorithm is applied to three-dimensional cone-beam image reconstruction and gets the better image. Key words: maximum likelihood estimate, CT Reconstruction, Penalized maximum likelihood, ordered subsets I 目目 錄錄 第一章 緒論 1.1 課題研究背景及意義 1 1.2 國內外研究現狀 1 1.3 論文的主要工作及結構安排 4 第二章 CT 統計重建的基礎 2.1 CT 成像理論.5 2.2 統計重建的數學模型 7 2.3 優化變換原理 8 2.4 圖像重建算法評價.9 第三章 基于最大似然準則的 CT 統計重建算法 3.1 最大似然準則的基本理論 11 3.2 ML-EM 算法.12 3.3 ML-SPS 算法14 3.4 有序子集的統計重建方法.17 3.4.1 有序子集 ML-EM17 3.4.2 有序子集 ML-SPS19 3.5 實驗結果及分析 20 3.5.1 ML-EM 算法與 ML-SPS 算法的重建結果21 3.5.2 OS-ML-EM 算法與 OS-ML-SPS 算法的重建結果.23 第四章 基于罰似然的 CT 統計重建算法 4.1 罰似然的基本理論.25 4.2 OSL-EM 算法 .31 II 4.3 OS-PL-SPS 算法.33 4.3.1 PL-SPS 算法.33 4.3.2 有序子集 PL-SPS 算法.34 4.4 自適應選擇正則項參數 35 4.5 實驗結果與分析.36 4.5.1 OSL-EM 算法實驗結果分析.37 4.5.2 OS-PL-SPS 算法實驗結果分析 38 4.5.3 自適應選擇正則項參數的實驗結果分析40 第五章 錐束軌跡下有序子集罰似然統計(PL)重建 5.1 錐束掃描下的投影矩陣的生成 42 5.2 錐束下的 OS-OSL-EM 算法重建結果及分析.44 第六章 總結與展望 6.1 研究主要內容及成果 47 6.2 存在的問題及以后的工作展望 48 參考文獻 攻讀碩士學位期間發表的論文 致 謝 1 第一章 緒論 1.1 課題研究背景及意義 計算機層析成像技術是利用有某種能量的射線源對物體進行斷層掃描，根據所獲 得的某種物理量(指物質對射線的衰減系數)，運用特定的重建算法，重建出物體某個斷 層的無影像重疊的二維圖像。自 1967 年 CT(Computed Tomography)在英國問世以來， 技術迅速發展，被公認為是 20 世紀影響人類發展的十大技術之一，其應用涉及到醫學、 工業、石油、物探、材料、生物、安全等眾多領域。在國防工業中，CT 已成為航天器、 固體火箭發動機、高新技術等產品關鍵部件檢測的必不可少的工具。由于 CT 技術是發 展國防的基礎技術之一，被發達國家列為敏感技術，作為國家機密進行嚴格的控制。 因此，CT 技術被視為國家科技實力的標志之一。 從數學上說，CT是從投影到重建圖像的反問題1，具有普遍性，在數學界已經引 起了廣泛的重視。重建算法是工業CT技術中較為關鍵的一部分，CT重建算法主要可以 分為濾波反投影重建算法和迭代重建算法。當投影數據完備、噪聲不很嚴重時，解析 法(濾波反投影(FBP)23等)可以得到很好的重建圖像，但當CT系統在數據采集過程中 受到嚴重的噪聲干擾，或者投影數據采集不完全時，用解析重建算法得到的圖像有偽 跡。而迭代重建算法(統計迭代算法和代數迭代算法等)可以用來重建不完全數據、動態 數據和噪聲數據。其中統計迭代重建方法的優勢是其對噪聲的魯棒性，通過加入先驗 模型和統計規律可以有效的抑制噪聲，提高圖像質量。 本課題主要研究統計算法在 CT 重建中的應用，通過研究有序子集的方法來提高收 斂速度，加入懲罰項，約束項來降低噪聲的影響，提高圖像質量，為國內三維工業 CT 的研制提供理論基礎和技術支持。 1.2 國內外研究現狀 2 近十幾年來， CT 廣泛的應用在醫學和工業上，相應的 CT 算法也迅速發展。其 中解析算法重建速度較快，易于實現，是目前 CT 圖像重建技術中最常用的算法。但是 該算法通常要求采樣數據是完全的，而且對噪聲比較敏感，這在一定程度上影響了重 建圖像的質量，限制了它在實際中的應用。相對解析算法、代數迭代算法而言，統計 迭代算法能準確描述投影數據的物理模型，對誤差不敏感，易于加入約束條件等優點 逐漸引起人們的重視，然而它與代數重建算法4一樣重建速度慢，運算時間長，這些缺 點極大地限制了它的應用，但是隨著加速算法和計算機硬件技術的發展，統計重建將 在實際中得到了廣泛地應用。 一個完整的統計重建框架包括圖像的離散化系統的物理幾何模型測量的統 計模型(泊松分布)優化準則(最大似然準則和罰似然準則)求解方法。第步是考慮 如何對重建的圖像進行建模，將待重建的圖像離散成二維或者三維的像素矩陣，第 步考慮如何確定投影矩陣，與系統的幾何結構和建模方法相關5-10，第步 ij AaA 確定投影數據的模型，在統計重建中，一般認為投影數據服從泊松分布。第步是引 入一定的優化準則。第步在的基礎之上進行問題的求解，該問題的求解實際上是 一個參數估計的數學優化問題。 統計重建的優化方法中涉及到優化變換(optimization transfer，OT)，增量方法 (incremental method)，有序子集加速(ordered subset，OS)，這三種優化方法是統計重建 方法發展的主線11。 在統計重建中，最大似然估計(Maximum Likelihood，ML)是一種常用方法。1977 年，A. P. Dempster 等人把期望值最大化(Expectation Maximization, EM)算法引入圖像 重建之中。1982年，Shepp12和Vardi13首次把期望最大迭代算法14應用到CT算法中， 使似然函數最大化，之后ML-EM算法得到了廣泛的應用。1997年，Fessler15等人提出 了基于替代函數(surrogate function)的優化理論，其基本思想是在M步(maximization step)直 接對E步(expectation step)中所得到目標函數的替代函數作最大化處理，這一突破性的改 進在很大程度上應當歸功于De Pierro的凸性(convexity)算法1617。當然替代函數的選取 需要滿足一定的條件，其它替代函數算法還包括Fessler18-21，Zheng22提出的算法以及 全局收斂的增量優化傳遞算法。 3 1994年，Hudson和Larkin提出了有序子集(Ordered Subsets) 23EM算法(OS-EM)，大 大提高了收斂速度。OS-EM算法將投影數據分解成有限個有序子集，每次迭代時只處 理其中一個子集的數據，由于加速效果明顯，OS-EM算法及其各種變形算法此后便被 廣泛地應用。其中RAMLA(row-action maximun likelihood algorithm)24重建算法的思想 來源于Herman的代數重建算法和Hudson等人的OS-EM算法，該算法也是將投影數據分 成一系列不正交的投影數據子集，并引入一個松弛因子，該算法可以從理論上證明其 收斂到ML的收斂點，并且其收斂速度和OS-EM算法差不多。雖然子集類算法收斂速度 快，但不能全局收斂，迭代到一定次數會出現極限環，為了解決該問題，S. Ahn 和J. A. Fessler提出了改進的塊迭代EM算法25，E .Quan和D .S Lalush等人提出了基于快速子 集的凸算法在CT中的應用26。 最大后驗概率重建算法(maximum a posteriori，MAP)2728是一種貝葉斯重建算法, 也可以稱為罰似然重建(Penalized Likelihood，PL)29-32，該算法和 ML 重建算法可以說 是一對關系親密的姊妹，因為 ML 算法是尋求合適的圖像的估計值使得圖像得到的已 知投影的概率最大；而貝葉斯重建算法是從已知的投影出發，要求在給定投影情況下 所求圖像的概率最大，通過選擇合適的先驗分布模型來提高重建結果的質量，其作用 等同于優化理論中的正則項（懲罰項） 。正則化技術最早由 Tikhonov 在著作中提出， 他研究了求解不適定問題穩定解的基本理論，并提出了著名的 Tikhonov 正則化技術。 先驗模型的引入雖然能改善重建效果，但同時增加了估計的難度。就泊松觀測模型而 言，它本身屬于非線性模型，因此它的 ML 估計和 MAP 估計均沒有閉型解析式。針對 這一問題，Herbert 最先提出 GEM(generalized EM)算法33，Green 給出了經典的 OSL(one-step-late)算法34，Fessler35提出了罰似然的 SAGE(space-alternating generalized EM)算法，另外 De Pierro 的3637凸性算法也涉及到了最大后驗估計的問題。在以后的 工作中，De Pierro 和 Yamagishi38將 RA(row-action)思想與最大后驗估計的思想相結合， 提出了 BSREM(block sequential regularizedexpectation maximization)算法。Chung Chan, Roger Fulton39等人提出了自適應結構先驗在 CT 中的應用，根據圖像結構塊中灰度的 不同自適應的選擇濾波。Xiaochuan Pan40闡述了最小化 TV 正則項(total variation)在圖 像重建中的應用。Jing Wang, Tianfang Li41等人把基于最小二乘的罰函數方法用于投影 數據的去噪，提高了圖像重建的質量。 4 1.3 論文的主要工作及結構安排 本文 CT 統計迭代重建算法為研究課題，主要做了以下工作： 第一章主要介紹了課題的研究背景、意義，綜述了 CT 及其發展歷史，CT 技術的 發展意義及統計重建算法的研究現狀。最后提出了本文的研究內容，確定了本文的研 究方向。 第二章首先給出了 CT 統計重建的物理和數學基礎，然后介紹了優化變換原理，最 后介紹了圖像重建質量的評價標準。 第三章首先詳細介紹了最大似然估計類的重建算法：ML-EM 算法和可分離拋物面 型替代(Separable Paraboloidal Surrogate，ML-SPS)算法，它們的推導過程和實現步驟， 然后對有序子集的方法進行闡述，并詳細介紹了 OS-ML-EM 算法和 OS-ML-SPS 算法 的實現步驟。實驗中對比了 ML-EM 算法和 ML-SPS 算法的實驗結果，并用 OS-ML- EM 算法和 OS-ML-SPS 算法對 CT 采集的實際投影數據進行重建，得到了較好的重建 結果。 第四章先介紹了罰似然估計的理論知識，重點介紹了基于 Gibbs 分布的罰函數， 從理論上分析了勢函數需要滿足的條件以及勢函數的選取對圖像的影響。然后介紹了 OSL (One Step Late)-EM 算法和 PL-SPS 算法的推導和實現，接著提出了自適應選擇正 則項參數的罰似然重建，通過對正則項參數的改進來提高圖像的質量。算法實現中選 用了二次勢函數和 Huber 勢函數來重建圖像，對它們的重建結果進行誤差分析，并用 提出的自適應選取正則項參數的方法來重建 CT 實際投影數據，取得了較好的結果，驗 證了方法的實用性。 第五章介紹了 OS-OSL-EM 算法在三維錐束軌跡下的應用，首先介紹了錐束掃描 結構下的投影矩陣的計算，然后用 OS-OSL-EM 算法對仿真投影數據進行重建，并與 OS-OSL-EM 算進行比較分析。 第六章對本課題的研究內容進行了總結，指出本文的創新之處和解決的問題，并 提出了尚未解決的問題，指明了今后研究的重點。 5 第二章 CT 統計重建的基礎 在過去的幾十年里，統計迭代法一直是斷層重建研究的熱點，其重建圖像的過程就 是利用參數估計，優化變換原理等來逼近最優解一個過程。在重建過程中，首先確定 能夠正確反映成像系統的物理模型和統計模型以及輸入輸出關系，然后利用已知的條 件確定求解方法。而且統計重建可以利用與物體相關的約束條件和其他的與特定系統 相關的邊界條件，利用精確的物理模型和適當的統計模型，可以在數據缺失和不規則 采樣的情況下重建出滿意的圖像。 2.1 CT 成像理論 CT 圖像重建的原理是用穿透力強的射線掃描被測物體，當一束強度大致均勻的射 線投照到物體上時，射線一部分被吸收和散射，另一部分透過物體沿原方向傳播。由 于物體的各種結構在密度、厚度等方面存在不同，它們對射線的吸收量也不相同，從 而使透過物體的射線強度分布發生變化，通過 CT 采集、轉換得到投影數據，并運用相 應的數學模型，進行圖像重建出二維或三維圖像，最后通過算法將射線強度分布轉換 成圖像上的灰度分布，形成人眼可視的圖像。 取一理想的 X 射線源，它發出極細的筆束 X 射線，在被測物體對面置一探測器， 如圖 2.1(a)所示。假設強度為的 X 射線穿過均勻分布衰減系數為的物體，行進了 0 Iu 的距離，強度變為，按 Beer 定理有xI 00 ln(/ ) ux ii euxII 或 （2.1） 如圖 2.1(b)所示，若物體是分段均勻的，衰減系數分別是,相應的長度為 1,23 ,.u u u 則下式成立： 1,23 ,.x x x 1230 .ln(/ )uxuxuxII （2.2） 6 更一般，若物體在平面內為不均勻介質，則在某一方向 ，射線沿某一路徑的總衰xylL 減值為 00 ( , )ln(/ )( )exp( , ) LL x y dlIII LIx y dl （2.3） 其中表示物體在二維平面上的衰減系數。是檢測到的強度；是射線經( , )x y( )I LL 過的直線，是直線的弧長。dl 若沒有具體路徑，只說了沿著某一個方向，那么是投影，由等式(2.3)可知，dl 測得與，即可知道。在 CT 中，入射的射線強度與出射的射線強度之比經過 0 IIdl 取對數運算后，則成了沿射線路徑上衰減系數的線積分。CT 重建問題實際上是給定一 個待重建物體的被測量的投影即線積分，然后運用算法去計算物體的衰減分布即dl 被積函數。 x o I 射線源 I 探測器 o I 射線源 I 探測器 1 x 2 x 3 x N x 1 2 3 N （a）均勻物體 （b）非均勻物體 圖 2.1 單能 X 射線束在物體中的衰減圖示 1917 年，丹麥數學家雷當(JRadon)研究已經為 CT 技術建立了數學基礎，從數學 理論上證明了某種物理參量的二維分布函數，由該函數在其定義域內的所有線積分確 定。該結論指出了如果沒有無窮多個且積分路徑互不完全重疊的線積分，只能確定實 際分布的一個估計近似值。 有了上述數學、物理基礎后，為了實現工程技術上的應用，還需要解決兩個主要 問題。首先是如何采集到檢測斷層衰減系數線積分的數據集，其次是如何充分利用該 數據集確定出衰減系數的二維分布。解決第一個問題可采用不同的掃描方式，即用射 7 線束穿過被測體所檢測斷層并相應進行射線強度測量的規律性，可采用各種不同的掃 描檢測模式圍繞提高掃描檢測效率。解決第二個問題則是選擇合適的圖像重建算法， 一類是解析法，以濾波反投影算法為代表。另一類是迭代重建法，有代數重建和統計 重建。統計重建是本論文研究的主要方向，下面具體介紹該算法。 2.2 統計重建的數學模型 在 CT 圖像重建中，圖像區域是重建區域，首先應當明確，圖像區域是一個正方形， 其中心在坐標原點；圖像函數是一個二元函數，其值在圖像區域之外為零。一個具有 個元素的網格把這個圖像區域分成 2 n個相等的正方形，每個小正方形(像素)內圖像函n 數的值是均勻相等的。一幅圖像的數字化是這樣一幅的數字化圖像，使得原nnnn 圖像在 2 n個元素網格的任何像素上的積分值等于在同一個像素上的數字化積分值。在 某一點的圖像密度值就是在該點物質組織在某種能量下的相對衰減系數值，將物體的 重建區域離散化為一幅數字化圖像后，直接從投影數據來估計數字化圖像每個像素點 的密度。 圖2.2 CT數據采集 如圖2.2所示，將待重建物體劃分成個互不重疊的子區域，像素序列依次從1到M ，射線源和探測器是點狀的，它們之間的連線對應一條射線。M 核物理研究表明，X光子的輻射滿足泊松隨機過程，結合這一統計特性，在CT中選 用泊松統計模型，可以表示成如下： 8 iii ypoissonArx:(1,2. )iN （2.4） 其中表示重建圖像矢量，用表示隨機變量的個測量 12 ( ,.) M x xxx 12 (,.,) N y yyyN 值, 表示第 個投影的測量誤差，為投影矩陣，表示第 個射線穿過第個 i ri ij Aa ij aij 像素的長度。表示第 條射線的投影值。 1 M iijj j Aa x xi 2.3 優化變換原理 優化變換42(optimization transfer，簡稱OT)的思想是將復雜的統計重建的優化問題 轉換成易于求解的優化問題，轉換后的問題一般具有易于優化、分離變量、降低維數 的特點。 OT的優化問題是尋求目標函數的最大值，表達式如下：( ) x argmax( ) x D xx （2.5） 由于該目標函數的最優解比較難求，可以選擇一個用來逼近目標函數極大值的替 代函數，通過不斷的迭代來求解使替代函數取得極大值。之所以要用迭代的 ( ) ( ,) n x x 方法來求解最小值問題是由替代函數的特性決定的，替代函數可能為非二次曲線，不 易求最大值，計算量非常大，這些都給求解帶來了困難，所以一種可行的方法就是選 擇一個可以用來逐步逼近替代函數極大值的替代函數，通過不斷的迭代計算來得到替 代函數的極值。但是替代函數一定要滿足一定的條件才可以得到穩定的結果，首先替 代函數的目的就是為了簡化求解極大點的計算量，如果替代函數是一個多維的非二次 函數，就難于計算。那么替代函數就要盡量將多維的求解極值問題變為一維的問題來 簡化求解過程。另外為了保證收斂的一致性，如圖 2.3 所示，要求替代函數在當前點與 替代函數相切并且位于替代函數的下部，則替代函數滿足的條件如下： ( )( )( ) (;)() nnn xxx 9 （2.6） ( ) ( ;)( ) n x xx （2.7） 由條件(2.6)，(2.7)可以得出 ( )( )( )( ) ( )()( ;)(;) nnnn xxx xxx （2.8） 換句話說，我們選擇 使得來使得x ( )( )( ) ( ;)(;) nnn x xxx ( ) ( )() n xx 圖 2.3 替代函數曲線圖 其中為替代函數，為第 次迭代后的結果。只有替代函數滿足了這上述條 ( ) ( ,) n x x ( )n xn 件才能在簡化計算的同時保證得到穩定一致的收斂結果。 選定了替代函數以后，求解最大值的過程就可以簡化成式(2.9)所示的迭代過程， 每次迭代都是以上次的結果為一點，得到在這一點上的替代函數并求解極大點作為本 次迭代的結果。則該優化問題的迭代過程可以簡化成： 1( ) 0 argmax ( ,) nn x xx x （2.9） 2.4 圖像重建算法評價 10 對于重建算法的可行性,一般按照下面三個方面來評價: 1) 算法是否易于實現。 2) 計算量的大小是衡量算法可行性的一個重要指標。計算時間是一個相對量, 在 同一計算機上，可以通過幾種算法計算時間來衡量。 3) 迭代算法的誤差分析。用歸一化均方差距離、歸一化平均絕對距離、相關系數、 信噪比等評價指標反映重建圖像的誤差和噪聲。 相關系數 CC(Correlation Coefficient)衡量： 1 1/2 22 11 ()() ()() P recrecrefref jjjj j PP recrecrefref jjjj jj xxxx CC xxxx （2.10） 歸一化平均絕對距離 NMAE (Normalized Mean Absolute Error)來衡量： 11 ()/() MM recrefref jjj jj MAExxx （2.11） 信噪比 SNR (Signal Noise Ratio)： 22 11 10log() /() ) MM refrefrecref jjjj jj SNRxxxx （2.12） 歸一化均方差距離 NMSE(Normalized Mean Square Error): 1 2 2 2 recref jj i refref jj i xx MSE xx （2.13） ，分別為被檢測圖像和重建圖像對應像素的值，為被檢測圖像的平均灰度 ref j x rec j x ref j x 值，為重建圖像的平均灰度值。相關系數反映了被測試圖像和重建圖像之間的相 rec j x 11 似程度，相關系數越大，NMAE 越小說明重建效果越好。用信噪比 SNR 來評價圖像的 噪聲的大小，SNR 越大，圖像所含的噪聲越小。NMSE 越小表示圖像與原始圖像的誤 差越小，圖像質量越好。 第三章 基于最大似然準則的 CT 統計重建算法 最大似然估計考慮了數據的統計特性，能夠充分利用系統的固有分辨率，從而重 建的圖像分辨率和噪聲特性均優于卷積反投影重建的重建結果43。但是在圖像重建優 化求解的過程中，似然函數的解析表達式很難找到，因而限制了 ML 的使用。但是期 望最大法通過構造“完全數據集”的似然函數，充分利用完全數據似然函數的簡單性 去求解實際投影數據的最大似然估計，從而解決了這個難題。 3.1 最大似然準則的基本理論 圖像重建的極大似然模型是建立在假定（X 射線發射的光子數是服從泊松 （Poisson）分布）的基礎上的。最大似然準則為：假設是相互獨立的，已知的條yy 件下尋求x使下面的對數似然函數最大 ( )lnLPxy x （3.1） 12 基于 Poisson 分布的似然函數為P y x 1 () ( ) ! i i y N Ax i i i A SPe y x xy x （3.2） 假定滿足以下兩個條件 ij a 01 ij a （3.3） 1 1, M ij i a (1,2,)jM （3.4） 則(3.2)可以簡化得到下式： 1 1 1 () ( ) ! iM ijj j M y ijj a x N j i i a x Se y x （3.5） 圖像重建的極大似然方法就是求出圖像矢量，使似然函數最大，也 12 , T M x xxxL 就是使對數似然函數最大。因此，ML 估計為：( )lnLPxy x argmax ( ) ML Lxx 1 argmaxln()ln(!) N ii ii i AyAy xx 11111 argmax(log()log(!) MNNMN jijiijji jiiji xaya xy （3.6 ） 另外，還可以根據判定極值的二階充分條件（即目標函數的 Hessian 矩陣是半正定 的）來說明似然函數的凸性，從而可知 ML 估計(3.6)存在且唯一。 13 3.2 ML-EM 算法 目標(3.6)的優化具有一定的難度，由于目標函數是非線性的，不存在解的閉型表 達式；而且在求解過程中隱含著非負性約束條件，為了解決該問題，Shepp44和Vardi45建 議使用EM算法對(3.6)進行優化。1994年，EM 算法被成功地應用在圖像重建中，成為 研究不完全的數據極大似然估計的一種方法。EM 算法的收斂解是非負的，迭代形式 便于計算機實現，在圖像重建算法研究中，EM算法的應用相當重要。 EM算法認為直接被觀測到的投影數據總是不完全(incomplete)的，依賴一個較為y 完全的數據，它們之間的關系由映射表征，即，一般 為已知。稱完全R( )yzxZ 數據，它包含了待求未知參數 的所有信息。基本思想就是通過重復處理完全數據的統x 計問題來解決不完全數據的統計問題。首先，假設初始圖像為；接著根據求一 (0) x (0) x 次近似圖像，再由求二次近似圖像，以此類推到第 次迭代的圖像。 (1) x (1) x (2) xn ( )n x 根據ML估計，優化的目標函數為，利用優化變換原理，需要尋找( )lnPxy x 一個替代函數，其必須滿足(2.9)中的單調性的條件。EM算法在理論上總能 ( ) ( ,) n EM x x 以概率收斂到某個的ML估計46。 ML-EM 算法的迭代過程分兩步： E步：在給定當前圖像估計及測量得到投影值的情況下，完全數據對應的對 ( )n xy 數似然函數的條件期望為， ( )( ) ( ,)log(), nn EM EZ x xxy x （3.7） M步：最大化替代函數對圖像向量 進行估計x 1( ) 0 argmax( ,) nn EM x xx x （3.8） 根據計算步驟，需要找出替代函數，推導過程中應用了算法中的凹性特征，假設，0 i r 當時，則有 ( ) 0 n x 1 ( ) M iijj i j lAa x xx 14 111 ( )log()() NMM ijjiijji ijj La xra xr x ( ) ( )( ) ( )( )( ) 111 log()()() ()() n NMM ijjj nni iiiijji nnn ijj iji a xx r ylla xr lxl xx xx gg （3.9 ） 利用凹性特征有： （3.10 ( ) ( )( ) ( )( )( ) 111 ( )log()log ()() ()() n NMM ijjj nni iiiijji nnn ijj iji a xx r Lyl xl xa xr l xxl x xgg ） 忽略(3.9)式中的常數項，即可得到替代函數 ( ) ( ,) n EM x x ( ) ( ) ( ) 11 ( ,)() log () n MN ijjn EMijjj n ji i a x yxa x l x x xg （3.11） 其中，根據 M 步對(3.11)最大化，求關于任意的偏導數并令其為零即 1 N jij i aa j x ，可以立即得到： ( ) ( ,) 0 n EM j x x x ( ) ( )11 1 0 n NN jiij ij M nii j ijj j xy a a x a x （3.12） 即推到出 ML-EM 算法的迭代公式： ( ) (1)(1) ( ) 1 11 ,1,2.,0 n N jijinn jjNM n i ijijj ij xa y xjM x aa x （3.13） 根據迭代公式看出 EM 算法有兩個明顯的特點： 1） 每一次迭代后均提高后驗密度函數值，即有。 (1)( ) ()() nn ll xx 15 2） 給定非負的初始估計，以后迭代所得的也都是非負的，滿足了 (0) x (1)(2) ,xx 圖像像素的非負性約束。 ML-EM 算法的具體步驟如下： 1) 設定圖像的初始估計值: (0)( 1,2.,); jj xxjM 2) 計算所有射線方程的估計值 (0) 1 M iijj j la x (1,2,)iN 3) 計算 i i i y l (1,2,)iN 4) 計算第個像素的修正值 j 1 1 1 N jiijN i ij i Ca a 5) 對的值進行修正 j x ; )0( jjj Cxx 通過該像素的所有射線對該像素進行修正，從而完成第一輪迭代。 6) 將以上一輪迭代的結果作為初值，重復 2)到 5)可得第 n 輪結果，就得到一序列 若符合收斂要求,即對事先給定的很小的正數存在正整數 NN，使 (0)(1)(2)( ) , n xxxx 當時，有 .nNN ( )(1) ,(1,2, ) nn jj xxjJ 3.3 ML-SPS 算法 SPS(separable paraboloidal surrogates)算法是全局收斂的，并且滿足像素的非負性 條件，迭代的初始收斂速度要比SAGE快。SPS 適用于ECT的ML具有Convex形式先驗 的MAP問題47-51，Erdogan52將其推廣到TCT的ML和MAP問題。該算法是從優化變換 的思想出發，通過構造易于最大化的二次型替代函數簡化EM算法的計算復雜度。 首先定義非負集合：:0, MM j RRxj x 定義集合 : M DRAiri xx 根據最大似然估計理論，最大似然SPS算法的目標函數表達式如下： 16 1 ( )( )( ) N pp ii i LhA xxx （3.14） 其中對數似然函數的表達式可以為： , 其中 ( )log()() iiii h lylrlr 1 ( ) M ijj j lAa x xx （3.15） 函數有如下性質：( ) i h l （1），其中 i lr ( )( ) ii h lh l ii lyr （2）在區間上，單調遞增，在區間上單調遞減。 , ii r l ( ) i h l , i l （3）是凹函數( ) i h l 在 ML-SPS 算法中，首先找到一維函數的拋物面型替代函數，替代( ) i h l ( ) ( ,) pn i ql l 函數需要滿足的條件：在當前迭代點處，它與函數相切且在函數 ( ) nn ii lAx( ) i h l 的下方。通過對做二階 Taylor 展開：( ) i h l( ) i h l ( )( )( )( )( )( )2 1 ( ,)()()()()() 2 pnnnnnn iiiiiiiii ql lh lh lllc lll )( ; ) N npn ii i QqAA x xxx （3.17） 利用的凹性，可以得到的可分離的二次替代函數 ( ) ( ;) n Q x x( )L x ( ) ( ,) pn qx x ( )( )( )( )( )( )( ) 1 ( ,)()() ()() ()() 2 pnnnnnnn j qdiag cLL x xxxxxxxxxx （3.18） 17 其中為的一階偏導， 為曲率，曲 ( ) () n Lx( )L x 1 ( )( ) N jijii i ca ac xx 1 M iij j aa ( ) i c x 率有以下幾種選擇： 最大曲率： 0 max( ) MC ii l ch l jj xxjM 2) 計算函數，第 條射線的累加和和曲率( ) ii h l jjj xxC 通過該像素的所有射線對該像素進行修正，從而完成第一輪迭代。 6）以上一輪迭代的結果作為初值，重復2)到5)進行新一輪迭代，直到取得符合收 斂要求的結果為止。 3.4 有序子集的統計重建方法 Hudson 和 Larkin 提出有序子集方法的核心思想是對投影數據進行子集劃分，在每 次迭代的過程中只使用一個子集的投影數據而不是全部的投影數據。該方法大大提高 了 ML-EM 重建速度，同樣 OS 也可應用到 MAP 類的算法中來提高重建速度。盡管 OS 算法能夠加快重建速度，但是 OS 重建的圖像常常不能收斂，隨著迭代次數的增加， 圖像的噪聲大大增加。解決 OS 發散的一個方法就是使用松弛參數即隨迭代次數增加而 逐漸減小步長。OS 方法還可以和其它算法進行結合派生出新的“OS”型算法5455， 如 OS-SPS。目前許多研究者作了大量的研究。這些研究有最優有序子集選取方法的討 論，OS 收斂性必須滿足的條件，有序子集平衡和一些推廣算法，例如 BI 迭代算法 RA（row action）算法等。這些推廣算法均具有良好的收斂性。 OS方法是把投影數據按照某種選擇規則分成個有序子集，子集水平(level)表示S 子集的數目。原算法中的一步計算被分為步，每一步是對其中的一個子集內的投影S 數據進行重建圖像。這樣所有的子集都對像素校正一次，即圖像更新次，完成一次S 19 迭代。以往的算法是所有的投影數據計算完一次之后，其像素值更新一次。在OS方法 中，所有的的投影數據計算完一次之后，其像素值已更新次。S 3.4.1有序子集ML-EM OS-EM 算法是將 OS 方法用于 EM 算法中，將投影數據分成個有順序 12 (,.) N y yyS 的子集（） 。接著在每個子集內使用標準的 EM 算法，其重建結果為下一個子 12 ,. S T TT 集的重建初始值，在第個子集對像素校正完后，則完成一次迭代。上次重建的結果S 為下次迭代的初值。在 EM 算法中，圖像的修正值是用所有的投影數據計算得到，而 在 OS-EM 中，圖像的修正值是用每個子集內的投影數據計算。則 OS-EM 的迭代公式 為： ( , ) (1,1)(1,1) ( , ) 1 ,1,2.,0 s s s n T jijisnsn jjTM s n i ijj ij j i xa y xjM x a xa （3.23） OS-EM 算法的實現步驟： 1) 設定圖像的初始估計值: (0)( 1,2.,); jj xxjM 2) 對第 個子集s a) 計算所有射線方程的估計值 (0) 1 M iijj j la x (1,2,)iN b) 計算 i i i y l (1,2,)iN c) 計算第個像素的修正值 j 1 1 1 N jiijN i ij i Ca a d) 對的值進行修正 j x (0) ;jjj xxC 3) 將重建出的代入下一個子集中，重復步驟(a)到(d)，這樣循環直到做完第個 j xS 20 子集，這就完成了第一輪迭代。 4）以上一輪迭代的結果作為初值，重復(2)到(3)進行新一輪迭代，直到有符合收 斂要求的結果為止。 子集的選擇一般應該遵循以下準則： 1）每個子集內的圖像信息和統計信息應該最大； 2）子集內的像素對每個子集的活動貢獻應該相同； 3