[理學(xué)]CT統(tǒng)計重建論文.doc

基于最大似然和罰似然估計的基于最大似然和罰似然估計的 CT 統(tǒng)計重建算法研究統(tǒng)計重建算法研究 摘 要 當(dāng) CT 系統(tǒng)在數(shù)據(jù)采集過程中受到嚴(yán)重的噪聲干擾，或者投影數(shù)據(jù)采集不完全時， 用解析重建算法重建出的圖像有偽跡。而統(tǒng)計重建算法具有物理模型準(zhǔn)確、對噪聲不 敏感、易于加入約束條件等優(yōu)點，其重建的圖像質(zhì)量要優(yōu)于傳統(tǒng)的 FBP 方法。針對 CT 統(tǒng)計重建，本文主要研究了以下內(nèi)容： 研究了基于最大似然估計的統(tǒng)計重建算法：期望最大值法(Maximum Likelihood Expectation Maximization，ML-EM)，可分離拋物面型替代(Separable Paraboloidal Surrogate，ML-SPS)算法和這兩種算法的有序子集形式（OS-ML-EM 和 OS-ML-SPS） 。 仿真實驗結(jié)果表明 ML-SPS 的初始收斂速度比 ML-EM 快，但該兩種算法的收斂速度 遠慢于其對應(yīng)的有序子集形式。用 OS-ML-EM 算法和 OS-ML-SPS 算法對 CT 采集的 實際投影數(shù)據(jù)進行重建，得到了較好的重建結(jié)果。 研究了基于罰似然(Penalized Likelihood，PL)的統(tǒng)計重建算法:OSL(One Step Late)- EM 算法和 PL-SPS 算法。重點討論了基于 Gibbs 分布的罰函數(shù)，從理論上分析了勢函 數(shù)需要滿足的條件以及勢函數(shù)的選取對圖像的影響，著重分析了二次勢函數(shù)和 Huber 勢函數(shù)的優(yōu)缺點，通過仿真實驗對它們的重建結(jié)果進行誤差分析。 提出了罰似然重建中正則項參數(shù)的自適應(yīng)選取的方法，該方法充分利用每一次迭 代的重建結(jié)果信息，不斷對正則項參數(shù)進行更新。并用于仿真實驗和重建 CT 實際投影 數(shù)據(jù)，重建結(jié)果表明該方法降低了噪聲，提高了圖像質(zhì)量。 將 OS-OSL-EM 算法推廣到三維錐束軌跡下的圖像重建，取得了較好的仿真結(jié)果。 關(guān)鍵詞：最大似然估計，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 國內(nèi)外研究現(xiàn)狀 1 1.3 論文的主要工作及結(jié)構(gòu)安排 4 第二章 CT 統(tǒng)計重建的基礎(chǔ) 2.1 CT 成像理論.5 2.2 統(tǒng)計重建的數(shù)學(xué)模型 7 2.3 優(yōu)化變換原理 8 2.4 圖像重建算法評價.9 第三章 基于最大似然準(zhǔn)則的 CT 統(tǒng)計重建算法 3.1 最大似然準(zhǔn)則的基本理論 11 3.2 ML-EM 算法.12 3.3 ML-SPS 算法14 3.4 有序子集的統(tǒng)計重建方法.17 3.4.1 有序子集 ML-EM17 3.4.2 有序子集 ML-SPS19 3.5 實驗結(jié)果及分析 20 3.5.1 ML-EM 算法與 ML-SPS 算法的重建結(jié)果21 3.5.2 OS-ML-EM 算法與 OS-ML-SPS 算法的重建結(jié)果.23 第四章 基于罰似然的 CT 統(tǒng)計重建算法 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 自適應(yīng)選擇正則項參數(shù) 35 4.5 實驗結(jié)果與分析.36 4.5.1 OSL-EM 算法實驗結(jié)果分析.37 4.5.2 OS-PL-SPS 算法實驗結(jié)果分析 38 4.5.3 自適應(yīng)選擇正則項參數(shù)的實驗結(jié)果分析40 第五章 錐束軌跡下有序子集罰似然統(tǒng)計(PL)重建 5.1 錐束掃描下的投影矩陣的生成 42 5.2 錐束下的 OS-OSL-EM 算法重建結(jié)果及分析.44 第六章 總結(jié)與展望 6.1 研究主要內(nèi)容及成果 47 6.2 存在的問題及以后的工作展望 48 參考文獻 攻讀碩士學(xué)位期間發(fā)表的論文 致 謝 1 第一章 緒論 1.1 課題研究背景及意義 計算機層析成像技術(shù)是利用有某種能量的射線源對物體進行斷層掃描，根據(jù)所獲 得的某種物理量(指物質(zhì)對射線的衰減系數(shù))，運用特定的重建算法，重建出物體某個斷 層的無影像重疊的二維圖像。自 1967 年 CT(Computed Tomography)在英國問世以來， 技術(shù)迅速發(fā)展，被公認(rèn)為是 20 世紀(jì)影響人類發(fā)展的十大技術(shù)之一，其應(yīng)用涉及到醫(yī)學(xué)、 工業(yè)、石油、物探、材料、生物、安全等眾多領(lǐng)域。在國防工業(yè)中，CT 已成為航天器、 固體火箭發(fā)動機、高新技術(shù)等產(chǎn)品關(guān)鍵部件檢測的必不可少的工具。由于 CT 技術(shù)是發(fā) 展國防的基礎(chǔ)技術(shù)之一，被發(fā)達國家列為敏感技術(shù)，作為國家機密進行嚴(yán)格的控制。 因此，CT 技術(shù)被視為國家科技實力的標(biāo)志之一。 從數(shù)學(xué)上說，CT是從投影到重建圖像的反問題1，具有普遍性，在數(shù)學(xué)界已經(jīng)引 起了廣泛的重視。重建算法是工業(yè)CT技術(shù)中較為關(guān)鍵的一部分，CT重建算法主要可以 分為濾波反投影重建算法和迭代重建算法。當(dāng)投影數(shù)據(jù)完備、噪聲不很嚴(yán)重時，解析 法(濾波反投影(FBP)23等)可以得到很好的重建圖像，但當(dāng)CT系統(tǒng)在數(shù)據(jù)采集過程中 受到嚴(yán)重的噪聲干擾，或者投影數(shù)據(jù)采集不完全時，用解析重建算法得到的圖像有偽 跡。而迭代重建算法(統(tǒng)計迭代算法和代數(shù)迭代算法等)可以用來重建不完全數(shù)據(jù)、動態(tài) 數(shù)據(jù)和噪聲數(shù)據(jù)。其中統(tǒng)計迭代重建方法的優(yōu)勢是其對噪聲的魯棒性，通過加入先驗 模型和統(tǒng)計規(guī)律可以有效的抑制噪聲，提高圖像質(zhì)量。 本課題主要研究統(tǒng)計算法在 CT 重建中的應(yīng)用，通過研究有序子集的方法來提高收 斂速度，加入懲罰項，約束項來降低噪聲的影響，提高圖像質(zhì)量，為國內(nèi)三維工業(yè) CT 的研制提供理論基礎(chǔ)和技術(shù)支持。 1.2 國內(nèi)外研究現(xiàn)狀 2 近十幾年來， CT 廣泛的應(yīng)用在醫(yī)學(xué)和工業(yè)上，相應(yīng)的 CT 算法也迅速發(fā)展。其 中解析算法重建速度較快，易于實現(xiàn)，是目前 CT 圖像重建技術(shù)中最常用的算法。但是 該算法通常要求采樣數(shù)據(jù)是完全的，而且對噪聲比較敏感，這在一定程度上影響了重 建圖像的質(zhì)量，限制了它在實際中的應(yīng)用。相對解析算法、代數(shù)迭代算法而言，統(tǒng)計 迭代算法能準(zhǔn)確描述投影數(shù)據(jù)的物理模型，對誤差不敏感，易于加入約束條件等優(yōu)點 逐漸引起人們的重視，然而它與代數(shù)重建算法4一樣重建速度慢，運算時間長，這些缺 點極大地限制了它的應(yīng)用，但是隨著加速算法和計算機硬件技術(shù)的發(fā)展，統(tǒng)計重建將 在實際中得到了廣泛地應(yīng)用。 一個完整的統(tǒng)計重建框架包括圖像的離散化系統(tǒng)的物理幾何模型測量的統(tǒng) 計模型(泊松分布)優(yōu)化準(zhǔn)則(最大似然準(zhǔn)則和罰似然準(zhǔn)則)求解方法。第步是考慮 如何對重建的圖像進行建模，將待重建的圖像離散成二維或者三維的像素矩陣，第 步考慮如何確定投影矩陣，與系統(tǒng)的幾何結(jié)構(gòu)和建模方法相關(guān)5-10，第步 ij AaA 確定投影數(shù)據(jù)的模型，在統(tǒng)計重建中，一般認(rèn)為投影數(shù)據(jù)服從泊松分布。第步是引 入一定的優(yōu)化準(zhǔn)則。第步在的基礎(chǔ)之上進行問題的求解，該問題的求解實際上是 一個參數(shù)估計的數(shù)學(xué)優(yōu)化問題。 統(tǒng)計重建的優(yōu)化方法中涉及到優(yōu)化變換(optimization transfer，OT)，增量方法 (incremental method)，有序子集加速(ordered subset，OS)，這三種優(yōu)化方法是統(tǒng)計重建 方法發(fā)展的主線11。 在統(tǒng)計重建中，最大似然估計(Maximum Likelihood，ML)是一種常用方法。1977 年，A. P. Dempster 等人把期望值最大化(Expectation Maximization, EM)算法引入圖像 重建之中。1982年，Shepp12和Vardi13首次把期望最大迭代算法14應(yīng)用到CT算法中， 使似然函數(shù)最大化，之后ML-EM算法得到了廣泛的應(yīng)用。1997年，F(xiàn)essler15等人提出 了基于替代函數(shù)(surrogate function)的優(yōu)化理論，其基本思想是在M步(maximization step)直 接對E步(expectation step)中所得到目標(biāo)函數(shù)的替代函數(shù)作最大化處理，這一突破性的改 進在很大程度上應(yīng)當(dāng)歸功于De Pierro的凸性(convexity)算法1617。當(dāng)然替代函數(shù)的選取 需要滿足一定的條件，其它替代函數(shù)算法還包括Fessler18-21，Zheng22提出的算法以及 全局收斂的增量優(yōu)化傳遞算法。 3 1994年，Hudson和Larkin提出了有序子集(Ordered Subsets) 23EM算法(OS-EM)，大 大提高了收斂速度。OS-EM算法將投影數(shù)據(jù)分解成有限個有序子集，每次迭代時只處 理其中一個子集的數(shù)據(jù)，由于加速效果明顯，OS-EM算法及其各種變形算法此后便被 廣泛地應(yīng)用。其中RAMLA(row-action maximun likelihood algorithm)24重建算法的思想 來源于Herman的代數(shù)重建算法和Hudson等人的OS-EM算法，該算法也是將投影數(shù)據(jù)分 成一系列不正交的投影數(shù)據(jù)子集，并引入一個松弛因子，該算法可以從理論上證明其 收斂到ML的收斂點，并且其收斂速度和OS-EM算法差不多。雖然子集類算法收斂速度 快，但不能全局收斂，迭代到一定次數(shù)會出現(xiàn)極限環(huán)，為了解決該問題，S. Ahn 和J. A. Fessler提出了改進的塊迭代EM算法25，E .Quan和D .S Lalush等人提出了基于快速子 集的凸算法在CT中的應(yīng)用26。 最大后驗概率重建算法(maximum a posteriori，MAP)2728是一種貝葉斯重建算法, 也可以稱為罰似然重建(Penalized Likelihood，PL)29-32，該算法和 ML 重建算法可以說 是一對關(guān)系親密的姊妹，因為 ML 算法是尋求合適的圖像的估計值使得圖像得到的已 知投影的概率最大；而貝葉斯重建算法是從已知的投影出發(fā)，要求在給定投影情況下 所求圖像的概率最大，通過選擇合適的先驗分布模型來提高重建結(jié)果的質(zhì)量，其作用 等同于優(yōu)化理論中的正則項（懲罰項） 。正則化技術(shù)最早由 Tikhonov 在著作中提出， 他研究了求解不適定問題穩(wěn)定解的基本理論，并提出了著名的 Tikhonov 正則化技術(shù)。 先驗?zāi)Ｐ偷囊腚m然能改善重建效果，但同時增加了估計的難度。就泊松觀測模型而 言，它本身屬于非線性模型，因此它的 ML 估計和 MAP 估計均沒有閉型解析式。針對 這一問題，Herbert 最先提出 GEM(generalized EM)算法33，Green 給出了經(jīng)典的 OSL(one-step-late)算法34，F(xiàn)essler35提出了罰似然的 SAGE(space-alternating generalized EM)算法，另外 De Pierro 的3637凸性算法也涉及到了最大后驗估計的問題。在以后的 工作中，De Pierro 和 Yamagishi38將 RA(row-action)思想與最大后驗估計的思想相結(jié)合， 提出了 BSREM(block sequential regularizedexpectation maximization)算法。Chung Chan, Roger Fulton39等人提出了自適應(yīng)結(jié)構(gòu)先驗在 CT 中的應(yīng)用，根據(jù)圖像結(jié)構(gòu)塊中灰度的 不同自適應(yīng)的選擇濾波。Xiaochuan Pan40闡述了最小化 TV 正則項(total variation)在圖 像重建中的應(yīng)用。Jing Wang, Tianfang Li41等人把基于最小二乘的罰函數(shù)方法用于投影 數(shù)據(jù)的去噪，提高了圖像重建的質(zhì)量。 4 1.3 論文的主要工作及結(jié)構(gòu)安排 本文 CT 統(tǒng)計迭代重建算法為研究課題，主要做了以下工作： 第一章主要介紹了課題的研究背景、意義，綜述了 CT 及其發(fā)展歷史，CT 技術(shù)的 發(fā)展意義及統(tǒng)計重建算法的研究現(xiàn)狀。最后提出了本文的研究內(nèi)容，確定了本文的研 究方向。 第二章首先給出了 CT 統(tǒng)計重建的物理和數(shù)學(xué)基礎(chǔ)，然后介紹了優(yōu)化變換原理，最 后介紹了圖像重建質(zhì)量的評價標(biāo)準(zhǔn)。 第三章首先詳細介紹了最大似然估計類的重建算法：ML-EM 算法和可分離拋物面 型替代(Separable Paraboloidal Surrogate，ML-SPS)算法，它們的推導(dǎo)過程和實現(xiàn)步驟， 然后對有序子集的方法進行闡述，并詳細介紹了 OS-ML-EM 算法和 OS-ML-SPS 算法 的實現(xiàn)步驟。實驗中對比了 ML-EM 算法和 ML-SPS 算法的實驗結(jié)果，并用 OS-ML- EM 算法和 OS-ML-SPS 算法對 CT 采集的實際投影數(shù)據(jù)進行重建，得到了較好的重建 結(jié)果。 第四章先介紹了罰似然估計的理論知識，重點介紹了基于 Gibbs 分布的罰函數(shù)， 從理論上分析了勢函數(shù)需要滿足的條件以及勢函數(shù)的選取對圖像的影響。然后介紹了 OSL (One Step Late)-EM 算法和 PL-SPS 算法的推導(dǎo)和實現(xiàn)，接著提出了自適應(yīng)選擇正 則項參數(shù)的罰似然重建，通過對正則項參數(shù)的改進來提高圖像的質(zhì)量。算法實現(xiàn)中選 用了二次勢函數(shù)和 Huber 勢函數(shù)來重建圖像，對它們的重建結(jié)果進行誤差分析，并用 提出的自適應(yīng)選取正則項參數(shù)的方法來重建 CT 實際投影數(shù)據(jù)，取得了較好的結(jié)果，驗 證了方法的實用性。 第五章介紹了 OS-OSL-EM 算法在三維錐束軌跡下的應(yīng)用，首先介紹了錐束掃描 結(jié)構(gòu)下的投影矩陣的計算，然后用 OS-OSL-EM 算法對仿真投影數(shù)據(jù)進行重建，并與 OS-OSL-EM 算進行比較分析。 第六章對本課題的研究內(nèi)容進行了總結(jié)，指出本文的創(chuàng)新之處和解決的問題，并 提出了尚未解決的問題，指明了今后研究的重點。 5 第二章 CT 統(tǒng)計重建的基礎(chǔ) 在過去的幾十年里，統(tǒng)計迭代法一直是斷層重建研究的熱點，其重建圖像的過程就 是利用參數(shù)估計，優(yōu)化變換原理等來逼近最優(yōu)解一個過程。在重建過程中，首先確定 能夠正確反映成像系統(tǒng)的物理模型和統(tǒng)計模型以及輸入輸出關(guān)系，然后利用已知的條 件確定求解方法。而且統(tǒng)計重建可以利用與物體相關(guān)的約束條件和其他的與特定系統(tǒng) 相關(guān)的邊界條件，利用精確的物理模型和適當(dāng)?shù)慕y(tǒng)計模型，可以在數(shù)據(jù)缺失和不規(guī)則 采樣的情況下重建出滿意的圖像。 2.1 CT 成像理論 CT 圖像重建的原理是用穿透力強的射線掃描被測物體，當(dāng)一束強度大致均勻的射 線投照到物體上時，射線一部分被吸收和散射，另一部分透過物體沿原方向傳播。由 于物體的各種結(jié)構(gòu)在密度、厚度等方面存在不同，它們對射線的吸收量也不相同，從 而使透過物體的射線強度分布發(fā)生變化，通過 CT 采集、轉(zhuǎn)換得到投影數(shù)據(jù)，并運用相 應(yīng)的數(shù)學(xué)模型，進行圖像重建出二維或三維圖像，最后通過算法將射線強度分布轉(zhuǎn)換 成圖像上的灰度分布，形成人眼可視的圖像。 取一理想的 X 射線源，它發(fā)出極細的筆束 X 射線，在被測物體對面置一探測器， 如圖 2.1(a)所示。假設(shè)強度為的 X 射線穿過均勻分布衰減系數(shù)為的物體，行進了 0 Iu 的距離，強度變?yōu)椋?Beer 定理有xI 00 ln(/ ) ux ii euxII 或 （2.1） 如圖 2.1(b)所示，若物體是分段均勻的，衰減系數(shù)分別是,相應(yīng)的長度為 1,23 ,.u u u 則下式成立： 1,23 ,.x x x 1230 .ln(/ )uxuxuxII （2.2） 6 更一般，若物體在平面內(nèi)為不均勻介質(zhì)，則在某一方向 ，射線沿某一路徑的總衰xylL 減值為 00 ( , )ln(/ )( )exp( , ) LL x y dlIII LIx y dl （2.3） 其中表示物體在二維平面上的衰減系數(shù)。是檢測到的強度；是射線經(jīng)( , )x y( )I LL 過的直線，是直線的弧長。dl 若沒有具體路徑，只說了沿著某一個方向，那么是投影，由等式(2.3)可知，dl 測得與，即可知道。在 CT 中，入射的射線強度與出射的射線強度之比經(jīng)過 0 IIdl 取對數(shù)運算后，則成了沿射線路徑上衰減系數(shù)的線積分。CT 重建問題實際上是給定一 個待重建物體的被測量的投影即線積分，然后運用算法去計算物體的衰減分布即dl 被積函數(shù)。 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 年，丹麥數(shù)學(xué)家雷當(dāng)(JRadon)研究已經(jīng)為 CT 技術(shù)建立了數(shù)學(xué)基礎(chǔ)，從數(shù)學(xué) 理論上證明了某種物理參量的二維分布函數(shù)，由該函數(shù)在其定義域內(nèi)的所有線積分確 定。該結(jié)論指出了如果沒有無窮多個且積分路徑互不完全重疊的線積分，只能確定實 際分布的一個估計近似值。 有了上述數(shù)學(xué)、物理基礎(chǔ)后，為了實現(xiàn)工程技術(shù)上的應(yīng)用，還需要解決兩個主要 問題。首先是如何采集到檢測斷層衰減系數(shù)線積分的數(shù)據(jù)集，其次是如何充分利用該 數(shù)據(jù)集確定出衰減系數(shù)的二維分布。解決第一個問題可采用不同的掃描方式，即用射 7 線束穿過被測體所檢測斷層并相應(yīng)進行射線強度測量的規(guī)律性，可采用各種不同的掃 描檢測模式圍繞提高掃描檢測效率。解決第二個問題則是選擇合適的圖像重建算法， 一類是解析法，以濾波反投影算法為代表。另一類是迭代重建法，有代數(shù)重建和統(tǒng)計 重建。統(tǒng)計重建是本論文研究的主要方向，下面具體介紹該算法。 2.2 統(tǒng)計重建的數(shù)學(xué)模型 在 CT 圖像重建中，圖像區(qū)域是重建區(qū)域，首先應(yīng)當(dāng)明確，圖像區(qū)域是一個正方形， 其中心在坐標(biāo)原點；圖像函數(shù)是一個二元函數(shù)，其值在圖像區(qū)域之外為零。一個具有 個元素的網(wǎng)格把這個圖像區(qū)域分成 2 n個相等的正方形，每個小正方形(像素)內(nèi)圖像函n 數(shù)的值是均勻相等的。一幅圖像的數(shù)字化是這樣一幅的數(shù)字化圖像，使得原nnnn 圖像在 2 n個元素網(wǎng)格的任何像素上的積分值等于在同一個像素上的數(shù)字化積分值。在 某一點的圖像密度值就是在該點物質(zhì)組織在某種能量下的相對衰減系數(shù)值，將物體的 重建區(qū)域離散化為一幅數(shù)字化圖像后，直接從投影數(shù)據(jù)來估計數(shù)字化圖像每個像素點 的密度。 圖2.2 CT數(shù)據(jù)采集 如圖2.2所示，將待重建物體劃分成個互不重疊的子區(qū)域，像素序列依次從1到M ，射線源和探測器是點狀的，它們之間的連線對應(yīng)一條射線。M 核物理研究表明，X光子的輻射滿足泊松隨機過程，結(jié)合這一統(tǒng)計特性，在CT中選 用泊松統(tǒng)計模型，可以表示成如下： 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 優(yōu)化變換原理 優(yōu)化變換42(optimization transfer，簡稱OT)的思想是將復(fù)雜的統(tǒng)計重建的優(yōu)化問題 轉(zhuǎn)換成易于求解的優(yōu)化問題，轉(zhuǎn)換后的問題一般具有易于優(yōu)化、分離變量、降低維數(shù) 的特點。 OT的優(yōu)化問題是尋求目標(biāo)函數(shù)的最大值，表達式如下：( ) x argmax( ) x D xx （2.5） 由于該目標(biāo)函數(shù)的最優(yōu)解比較難求，可以選擇一個用來逼近目標(biāo)函數(shù)極大值的替 代函數(shù)，通過不斷的迭代來求解使替代函數(shù)取得極大值。之所以要用迭代的 ( ) ( ,) n x x 方法來求解最小值問題是由替代函數(shù)的特性決定的，替代函數(shù)可能為非二次曲線，不 易求最大值，計算量非常大，這些都給求解帶來了困難，所以一種可行的方法就是選 擇一個可以用來逐步逼近替代函數(shù)極大值的替代函數(shù)，通過不斷的迭代計算來得到替 代函數(shù)的極值。但是替代函數(shù)一定要滿足一定的條件才可以得到穩(wěn)定的結(jié)果，首先替 代函數(shù)的目的就是為了簡化求解極大點的計算量，如果替代函數(shù)是一個多維的非二次 函數(shù)，就難于計算。那么替代函數(shù)就要盡量將多維的求解極值問題變?yōu)橐痪S的問題來 簡化求解過程。另外為了保證收斂的一致性，如圖 2.3 所示，要求替代函數(shù)在當(dāng)前點與 替代函數(shù)相切并且位于替代函數(shù)的下部，則替代函數(shù)滿足的條件如下： ( )( )( ) (;)() 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 替代函數(shù)曲線圖 其中為替代函數(shù)，為第 次迭代后的結(jié)果。只有替代函數(shù)滿足了這上述條 ( ) ( ,) n x x ( )n xn 件才能在簡化計算的同時保證得到穩(wěn)定一致的收斂結(jié)果。 選定了替代函數(shù)以后，求解最大值的過程就可以簡化成式(2.9)所示的迭代過程， 每次迭代都是以上次的結(jié)果為一點，得到在這一點上的替代函數(shù)并求解極大點作為本 次迭代的結(jié)果。則該優(yōu)化問題的迭代過程可以簡化成： 1( ) 0 argmax ( ,) nn x xx x （2.9） 2.4 圖像重建算法評價 10 對于重建算法的可行性,一般按照下面三個方面來評價: 1) 算法是否易于實現(xiàn)。 2) 計算量的大小是衡量算法可行性的一個重要指標(biāo)。計算時間是一個相對量, 在 同一計算機上，可以通過幾種算法計算時間來衡量。 3) 迭代算法的誤差分析。用歸一化均方差距離、歸一化平均絕對距離、相關(guān)系數(shù)、 信噪比等評價指標(biāo)反映重建圖像的誤差和噪聲。 相關(guān)系數(shù) 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） ，分別為被檢測圖像和重建圖像對應(yīng)像素的值，為被檢測圖像的平均灰度 ref j x rec j x ref j x 值，為重建圖像的平均灰度值。相關(guān)系數(shù)反映了被測試圖像和重建圖像之間的相 rec j x 11 似程度，相關(guān)系數(shù)越大，NMAE 越小說明重建效果越好。用信噪比 SNR 來評價圖像的 噪聲的大小，SNR 越大，圖像所含的噪聲越小。NMSE 越小表示圖像與原始圖像的誤 差越小，圖像質(zhì)量越好。 第三章 基于最大似然準(zhǔn)則的 CT 統(tǒng)計重建算法 最大似然估計考慮了數(shù)據(jù)的統(tǒng)計特性，能夠充分利用系統(tǒng)的固有分辨率，從而重 建的圖像分辨率和噪聲特性均優(yōu)于卷積反投影重建的重建結(jié)果43。但是在圖像重建優(yōu) 化求解的過程中，似然函數(shù)的解析表達式很難找到，因而限制了 ML 的使用。但是期 望最大法通過構(gòu)造“完全數(shù)據(jù)集”的似然函數(shù)，充分利用完全數(shù)據(jù)似然函數(shù)的簡單性 去求解實際投影數(shù)據(jù)的最大似然估計，從而解決了這個難題。 3.1 最大似然準(zhǔn)則的基本理論 圖像重建的極大似然模型是建立在假定（X 射線發(fā)射的光子數(shù)是服從泊松 （Poisson）分布）的基礎(chǔ)上的。最大似然準(zhǔn)則為：假設(shè)是相互獨立的，已知的條yy 件下尋求x使下面的對數(shù)似然函數(shù)最大 ( )lnLPxy x （3.1） 12 基于 Poisson 分布的似然函數(shù)為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） 圖像重建的極大似然方法就是求出圖像矢量，使似然函數(shù)最大，也 12 , T M x xxxL 就是使對數(shù)似然函數(shù)最大。因此，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 ） 另外，還可以根據(jù)判定極值的二階充分條件（即目標(biāo)函數(shù)的 Hessian 矩陣是半正定 的）來說明似然函數(shù)的凸性，從而可知 ML 估計(3.6)存在且唯一。 13 3.2 ML-EM 算法 目標(biāo)(3.6)的優(yōu)化具有一定的難度，由于目標(biāo)函數(shù)是非線性的，不存在解的閉型表 達式；而且在求解過程中隱含著非負(fù)性約束條件，為了解決該問題，Shepp44和Vardi45建 議使用EM算法對(3.6)進行優(yōu)化。1994年，EM 算法被成功地應(yīng)用在圖像重建中，成為 研究不完全的數(shù)據(jù)極大似然估計的一種方法。EM 算法的收斂解是非負(fù)的，迭代形式 便于計算機實現(xiàn)，在圖像重建算法研究中，EM算法的應(yīng)用相當(dāng)重要。 EM算法認(rèn)為直接被觀測到的投影數(shù)據(jù)總是不完全(incomplete)的，依賴一個較為y 完全的數(shù)據(jù)，它們之間的關(guān)系由映射表征，即，一般 為已知。稱完全R( )yzxZ 數(shù)據(jù)，它包含了待求未知參數(shù) 的所有信息?；舅枷刖褪峭ㄟ^重復(fù)處理完全數(shù)據(jù)的統(tǒng)x 計問題來解決不完全數(shù)據(jù)的統(tǒng)計問題。首先，假設(shè)初始圖像為；接著根據(jù)求一 (0) x (0) x 次近似圖像，再由求二次近似圖像，以此類推到第 次迭代的圖像。 (1) x (1) x (2) xn ( )n x 根據(jù)ML估計，優(yōu)化的目標(biāo)函數(shù)為，利用優(yōu)化變換原理，需要尋找( )lnPxy x 一個替代函數(shù)，其必須滿足(2.9)中的單調(diào)性的條件。EM算法在理論上總能 ( ) ( ,) n EM x x 以概率收斂到某個的ML估計46。 ML-EM 算法的迭代過程分兩步： E步：在給定當(dāng)前圖像估計及測量得到投影值的情況下，完全數(shù)據(jù)對應(yīng)的對 ( )n xy 數(shù)似然函數(shù)的條件期望為， ( )( ) ( ,)log(), nn EM EZ x xxy x （3.7） M步：最大化替代函數(shù)對圖像向量 進行估計x 1( ) 0 argmax( ,) nn EM x xx x （3.8） 根據(jù)計算步驟，需要找出替代函數(shù)，推導(dǎo)過程中應(yīng)用了算法中的凹性特征，假設(shè)，0 i r 當(dāng)時，則有 ( ) 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)式中的常數(shù)項，即可得到替代函數(shù) ( ) ( ,) n EM x x ( ) ( ) ( ) 11 ( ,)() log () n MN ijjn EMijjj n ji i a x yxa x l x x xg （3.11） 其中，根據(jù) M 步對(3.11)最大化，求關(guān)于任意的偏導(dǎo)數(shù)并令其為零即 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） 根據(jù)迭代公式看出 EM 算法有兩個明顯的特點： 1） 每一次迭代后均提高后驗密度函數(shù)值，即有。 (1)( ) ()() nn ll xx 15 2） 給定非負(fù)的初始估計，以后迭代所得的也都是非負(fù)的，滿足了 (0) x (1)(2) ,xx 圖像像素的非負(fù)性約束。 ML-EM 算法的具體步驟如下： 1) 設(shè)定圖像的初始估計值: (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) 將以上一輪迭代的結(jié)果作為初值，重復(fù) 2)到 5)可得第 n 輪結(jié)果，就得到一序列 若符合收斂要求,即對事先給定的很小的正數(shù)存在正整數(shù) NN，使 (0)(1)(2)( ) , n xxxx 當(dāng)時，有 .nNN ( )(1) ,(1,2, ) nn jj xxjJ 3.3 ML-SPS 算法 SPS(separable paraboloidal surrogates)算法是全局收斂的，并且滿足像素的非負(fù)性 條件，迭代的初始收斂速度要比SAGE快。SPS 適用于ECT的ML具有Convex形式先驗 的MAP問題47-51，Erdogan52將其推廣到TCT的ML和MAP問題。該算法是從優(yōu)化變換 的思想出發(fā)，通過構(gòu)造易于最大化的二次型替代函數(shù)簡化EM算法的計算復(fù)雜度。 首先定義非負(fù)集合：:0, MM j RRxj x 定義集合 : M DRAiri xx 根據(jù)最大似然估計理論，最大似然SPS算法的目標(biāo)函數(shù)表達式如下： 16 1 ( )( )( ) N pp ii i LhA xxx （3.14） 其中對數(shù)似然函數(shù)的表達式可以為： , 其中 ( )log()() iiii h lylrlr 1 ( ) M ijj j lAa x xx （3.15） 函數(shù)有如下性質(zhì)：( ) i h l （1），其中 i lr ( )( ) ii h lh l ii lyr （2）在區(qū)間上，單調(diào)遞增，在區(qū)間上單調(diào)遞減。 , ii r l ( ) i h l , i l （3）是凹函數(shù)( ) i h l 在 ML-SPS 算法中，首先找到一維函數(shù)的拋物面型替代函數(shù)，替代( ) i h l ( ) ( ,) pn i ql l 函數(shù)需要滿足的條件：在當(dāng)前迭代點處，它與函數(shù)相切且在函數(shù) ( ) 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） 利用的凹性，可以得到的可分離的二次替代函數(shù) ( ) ( ;) n Q x x( )L x ( ) ( ,) pn qx x ( )( )( )( )( )( )( ) 1 ( ,)()() ()() ()() 2 pnnnnnnn j qdiag cLL x xxxxxxxxxx （3.18） 17 其中為的一階偏導(dǎo)， 為曲率，曲 ( ) () 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) 計算函數(shù)，第 條射線的累加和和曲率( ) ii h l jjj xxC 通過該像素的所有射線對該像素進行修正，從而完成第一輪迭代。 6）以上一輪迭代的結(jié)果作為初值，重復(fù)2)到5)進行新一輪迭代，直到取得符合收 斂要求的結(jié)果為止。 3.4 有序子集的統(tǒng)計重建方法 Hudson 和 Larkin 提出有序子集方法的核心思想是對投影數(shù)據(jù)進行子集劃分，在每 次迭代的過程中只使用一個子集的投影數(shù)據(jù)而不是全部的投影數(shù)據(jù)。該方法大大提高 了 ML-EM 重建速度，同樣 OS 也可應(yīng)用到 MAP 類的算法中來提高重建速度。盡管 OS 算法能夠加快重建速度，但是 OS 重建的圖像常常不能收斂，隨著迭代次數(shù)的增加， 圖像的噪聲大大增加。解決 OS 發(fā)散的一個方法就是使用松弛參數(shù)即隨迭代次數(shù)增加而 逐漸減小步長。OS 方法還可以和其它算法進行結(jié)合派生出新的“OS”型算法5455， 如 OS-SPS。目前許多研究者作了大量的研究。這些研究有最優(yōu)有序子集選取方法的討 論，OS 收斂性必須滿足的條件，有序子集平衡和一些推廣算法，例如 BI 迭代算法 RA（row action）算法等。這些推廣算法均具有良好的收斂性。 OS方法是把投影數(shù)據(jù)按照某種選擇規(guī)則分成個有序子集，子集水平(level)表示S 子集的數(shù)目。原算法中的一步計算被分為步，每一步是對其中的一個子集內(nèi)的投影S 數(shù)據(jù)進行重建圖像。這樣所有的子集都對像素校正一次，即圖像更新次，完成一次S 19 迭代。以往的算法是所有的投影數(shù)據(jù)計算完一次之后，其像素值更新一次。在OS方法 中，所有的的投影數(shù)據(jù)計算完一次之后，其像素值已更新次。S 3.4.1有序子集ML-EM OS-EM 算法是將 OS 方法用于 EM 算法中，將投影數(shù)據(jù)分成個有順序 12 (,.) N y yyS 的子集（） 。接著在每個子集內(nèi)使用標(biāo)準(zhǔn)的 EM 算法，其重建結(jié)果為下一個子 12 ,. S T TT 集的重建初始值，在第個子集對像素校正完后，則完成一次迭代。上次重建的結(jié)果S 為下次迭代的初值。在 EM 算法中，圖像的修正值是用所有的投影數(shù)據(jù)計算得到，而 在 OS-EM 中，圖像的修正值是用每個子集內(nèi)的投影數(shù)據(jù)計算。則 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 算法的實現(xiàn)步驟： 1) 設(shè)定圖像的初始估計值: (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) 將重建出的代入下一個子集中，重復(fù)步驟(a)到(d)，這樣循環(huán)直到做完第個 j xS 20 子集，這就完成了第一輪迭代。 4）以上一輪迭代的結(jié)果作為初值，重復(fù)(2)到(3)進行新一輪迭代，直到有符合收 斂要求的結(jié)果為止。 子集的選擇一般應(yīng)該遵循以下準(zhǔn)則： 1）每個子集內(nèi)的圖像信息和統(tǒng)計信息應(yīng)該最大； 2）子集內(nèi)的像素對每個子集的活動貢獻應(yīng)該相同； 3