董 建，畢丹陽，楊耿煌，秦轉(zhuǎn)萍

（1.天津職業(yè)技術(shù)師范大學(xué)天津市信息傳感與智能控制重點實驗室，天津 300222；2.唐山工業(yè)職業(yè)技術(shù)學(xué)院自動化工程系，唐山 063299）

計算機斷層成像（computed tomography，CT）技術(shù)具有適用面廣、分辨率高、對比度好、價格相對低廉等優(yōu)勢，已成為使用最廣泛的無創(chuàng)臨床成像技術(shù)之一[1-2]。自1979年我國引入第一臺醫(yī)用CT 以來，CT 設(shè)備已經(jīng)在我國各級醫(yī)院工作了40年，早期投入使用的設(shè)備正相繼進入設(shè)備老化和故障期[3-4]。X 線探測器是CT設(shè)備的核心部件，也是故障高發(fā)部件。探測器故障會導(dǎo)致CT 圖像出現(xiàn)環(huán)狀偽影，嚴重影響正常組織結(jié)構(gòu)和病變區(qū)域的刻畫[5-6]。而現(xiàn)階段我國CT 研發(fā)及設(shè)備維護技術(shù)團隊已經(jīng)無力應(yīng)對各醫(yī)療機構(gòu)集中出現(xiàn)的設(shè)備故障問題，研發(fā)具有容錯機制的CT 圖像重建算法已經(jīng)迫在眉睫。據(jù)本團隊了解，還未發(fā)現(xiàn)國內(nèi)做相關(guān)研究的報告，本文提出的探測器故障下的CT 圖像重建的容錯算法開發(fā)尚屬解決此類問題的首例研究。

本文提出一種在未知探測器受損數(shù)量和具體位置信息的故障條件下，依然能夠精準重建的CT 圖像重建算法。傳統(tǒng)的迭代式重建算法采用正向投影與實測投影數(shù)據(jù)的最小二乘差，即L2 范數(shù)差作為最小化評價函數(shù)[7-8]。眾所周知，L2 范數(shù)對實測投影數(shù)據(jù)的異常值尤其敏感，這決定了此傳統(tǒng)方案無法避免CT 圖像中環(huán)狀偽影的出現(xiàn)。而一種能有效預(yù)測異常值位置并將其精準濾除的方案，即采用L1 范數(shù)為最小化評價函數(shù)。基于這一想法，開發(fā)了一種類似于代數(shù)式重建法（algebraic reconstruction technique，ART）的行處理型（row-action-type）迭代式重建算法[9-10]，該算法由凸優(yōu)化領(lǐng)域的近端分裂理論（proximal splitting）推導(dǎo)構(gòu)建而成[11-12]。此外，通過向L1 范數(shù)追加全變分（total variation，TV）懲罰項[13-16]，新算法在去除環(huán)狀偽影方面更具魯棒性。在經(jīng)典Shepp-Logan 模型重建實驗以及實際臨床CT 圖像的重建實驗中，新算法較傳統(tǒng)L2 范數(shù)重建算法表現(xiàn)出顯著優(yōu)越性。

1 算法介紹

1.1 問題定義

CT 設(shè)備獲取人體斷面圖像的掃描過程如圖1所示。

圖1 CT 設(shè)備獲取人體斷面圖像的掃描過程

重建一張人體斷層圖像，CT 設(shè)備需要對該斷面進行1 000～1 500 次X 射線掃描。在X 線的某一通路上，令X 線源發(fā)出的X 線的強度為I0，探測器接收到的被人體部分吸收后的X 線的強度為I1。在不考慮噪聲影響的前提下，I1與I0存在如下等式關(guān)系

式中：b 為投影數(shù)據(jù)；γ 為探測器對X 線的靈敏度。探測器正常工作條件下γ = 1，若探測器出現(xiàn)故障，則γ→0。此時投影數(shù)據(jù)b 為

即探測器故障情況下，測得的投影數(shù)據(jù)值遠小于探測器正常工作條件下的測得值。重建圖像中環(huán)狀偽影的產(chǎn)生機理闡述如下：在圖像重建的反投影這一關(guān)鍵步驟中，受損害的投影數(shù)據(jù)使得該通路上的各像素由于得不到有效恢復(fù)而取值大幅減小。此外，根據(jù)重建規(guī)則，反投影運算會在各角度方向依次進行，進而使得無法恢復(fù)正常取值的像素在某一路徑進行重合。180°的反投影運算使得該路徑逐漸演變?yōu)榘雸A形結(jié)構(gòu)，而且重合次數(shù)足夠多的像素取值偏大，重合次數(shù)少的像素取值偏小。該半圓形結(jié)構(gòu)即對應(yīng)重建圖像中由相鄰的亮區(qū)和暗區(qū)所組成的環(huán)狀偽影。完整、1 個和2 個探測器故障的投影數(shù)據(jù)與相應(yīng)濾波反投影算法（filtered back-projection，F(xiàn)BP）的重建圖像如圖2 所示。

圖2 完整、1 個和2 個探測器故障的投影數(shù)據(jù)與相應(yīng)FBP 算法的重建圖像

探測器無故障和隨機性出現(xiàn)1 個和2 個探測器故障時，由一般商業(yè)CT 搭載的FBP 算法重建了經(jīng)典Shepp-Logan 模型。由圖2 可知，完好投影數(shù)據(jù)可重建出理想圖像，探測器故障可引起投影數(shù)據(jù)的某一列數(shù)據(jù)異常，進而導(dǎo)致重建圖像中的環(huán)狀偽影。偽影出現(xiàn)個數(shù)與出現(xiàn)故障探測器個數(shù)一致。

CT 圖像重建即利用重建算法將投影數(shù)據(jù)轉(zhuǎn)化為CT 圖像的過程。迭代式重建算法可將問題定義為求解多元線性方程組 Ax = b 的問題，此處 x =（x1，x2，…，xJ）T為J 維未知量即待重建的CT 圖像，b=（b1，b2，…，bI）T為I 維已知量即實測投影數(shù)據(jù)，A={aij}為I×J 維系統(tǒng)矩陣[17-18]。本文研究I≥J 的超定方程組求解問題。迭代式重建領(lǐng)域的標準評價函數(shù)為L2 范數(shù)差

本研究將傳統(tǒng)L2 范數(shù)置換為L1 范數(shù)，如式（4），并通過將f（x）最小化而推導(dǎo)出行處理型迭代式重建算法。

2 種范數(shù)差展開式表示為

1.2 近端分裂理論

在實施算法推導(dǎo)前，有必要先對其基礎(chǔ)內(nèi)容的近端分裂理論及其對算法推導(dǎo)的重大貢獻做簡要介紹。首先考慮以下凸優(yōu)化問題

假定式中f（x）為可能不可微分的下半連續(xù)凸函數(shù)。近端算子法（proximal operator）為求解此類最優(yōu)化問題提供了良好解決方案。對應(yīng)于f（x）的近端算子表達式為

式中：α 為步長參數(shù)。

近端算子具備2 種優(yōu)秀屬性。

（1）它可針對任何一個下半連續(xù)凸函數(shù)進行定義，即使該函數(shù)不可微（如L1 范數(shù)）；

（2）在步長參數(shù)α ＞0 的前提下，滿足x=proxαf（g）的定點x 與使凸函數(shù)f（x）取最小值時的變量x 相一致。

這2 個性質(zhì)決定了可以應(yīng)用如下迭代式求得凸函數(shù)f（x）取得最小值時的變量x

式中：k 為迭代次數(shù)。

該迭代算法叫做近端最小化算法，它為不可微凸函數(shù)的最優(yōu)化問題提供了強有力的解決方案。

凸函數(shù)f（x）可拆分成若干分量函數(shù)fi（x），即其中fi（x）（i=1，2，…，I）均為可能不可微但滿足下半連續(xù)的凸函數(shù)。此時，假定凸函數(shù)f（x）的近端算子求解難度巨大，而其各個分量函數(shù)fi（x）的近端算子能夠由低計算成本求得。近端分裂理論能夠為該條件下的凸優(yōu)化問題提供解決方案，它的處理方式為在求解由x(k)到x(k+1)的過程中，進一步將問題分解為按照i=1，2，…，I 的順序進行各分量函數(shù)所對應(yīng)的近端算子的計算，即

當 i=I 時，x(k+1)=x(k，I+1)。

式中：k 為主循環(huán)變量；i 為子循環(huán)變量。

步長參數(shù)α(k)值會隨k 值改變。為了使所求解序列x(k)收斂于函數(shù)f（x）取最小值時的變量x，α(k)需滿足如下條件[12]

近端分裂理論為本文L1 范數(shù)重建算法的構(gòu)建提供了理論依據(jù)。

1.3 L1范數(shù)重建算法描述

式中：ai為系統(tǒng)矩陣A 所對應(yīng)的第i 行分量。進而求解各分量函數(shù)所對應(yīng)的近端算子，即求解以下最優(yōu)化問題

該近端算子的計算可由拉格朗日乘子技術(shù)簡單求得。首先，引入替代變量則以上問題轉(zhuǎn)變?yōu)榧s束極小化問題

對應(yīng)于式（13）的拉格朗日函數(shù)可定義為

式中：λ 為拉格朗日乘子，也稱作對偶變量。

對應(yīng)于式（14）的對偶問題可定義為

通過對該對偶問題的求解，可得到未知解x 的更新公式為

同時可得到對偶變量的表達式及取值范圍

可將λ 分情況表示為

將式（18）代入式（16），即可得到 x(k,i)的更新解x(k,i+1)=proxα(k)fi（x(k,i)）。由于步長參數(shù)α(k)需滿足隨循環(huán)變量k 而遞減的條件，故將其設(shè)定為

式中：α0＞ 0，? ＞ 0 為提前給定的正數(shù)常量。

至此，L1 范數(shù)重建算法已經(jīng)構(gòu)建完成。

1.4 L1范數(shù)重建算法

步驟0設(shè)定步長參量α0和? 的初始值，為變量x 取初值 x(0，1)，以 k=0，1，2，…的順序執(zhí)行步驟 1 到步驟3；

步驟1計算步長參數(shù)

步驟2以i=1，2，…，I 的順序執(zhí)行以下操作：

①計算λ

②更新x

步驟3x(k+1,1)=x(k,I+1)

為了得到更好的圖像重建效果，進一步向L1 范數(shù)的評價函數(shù)中追加了全變分懲罰項‖x‖TV，定義如下

式中：xm，n為在一幅二維圖像中坐標為（m，n）的像素的像素值。全變分模型是一個依靠梯度下降流對圖像進行平滑的各向異性的模型，能對圖像起到良好的去噪和平滑作用[13-16]。因此，L1-TV 重建的最小化評價函數(shù)為

式中：β 為可調(diào)節(jié)全變分作用強弱的超參數(shù)。

2 實驗結(jié)果及分析

2.1 少數(shù)探測器故障下的圖像重建

本文對經(jīng)典Shepp-Logan 數(shù)字模型進行重建實驗。該數(shù)字模型所含像素數(shù)為256×256，像素值的分布范圍為[0，4.5]。圖像重建領(lǐng)域的常用做法是獲取與重建圖像的像素數(shù)等量的投影數(shù)據(jù)量，因此本文以[0°，180°]和 180°/256 的角度間隔進行投影數(shù)據(jù)采集，投影數(shù)據(jù)的像素數(shù)同樣為256×256。

假定出故障的探測器個數(shù)為1 個和2 個，并且假定出現(xiàn)故障的探測器的位置具有隨機性。分別用L2范數(shù)重建算法、L1 范數(shù)重建算法和L1-TV 重建算法對含有異常值的投影數(shù)據(jù)進行重建實驗。經(jīng)實驗驗證，對于本實驗的Shepp-Logan 數(shù)據(jù)模型，在各參數(shù)取以下數(shù)值時可得到最佳重建效果。在各項實驗中α0=0.03，?=50，L1-TV 實驗中 β =260。各項實驗的迭代次數(shù)k 設(shè)定為100。在不同探測器故障個數(shù)條件下各算法的實驗結(jié)果如圖3所示。其中，圖3（a）、（c）、（e）為 1 個探測器故障時的重建結(jié)果，圖3（b）、（d）、（f）為2 個探測器故障時的重建結(jié)果。由結(jié)果圖像可知，L2范數(shù)重建的圖像仍含有嚴重的環(huán)狀偽影，并且在半圓環(huán)偽影的起點和終點處二次產(chǎn)生了直線型放射狀偽影。L2 范數(shù)的重建結(jié)果與FBP 算法的重建結(jié)果的圖像效果相似。L1 范數(shù)的重建效果較L2 范數(shù)的重建效果有了顯著提高，環(huán)狀偽影得以有效去除。然而，偽影處的黑色瘢痕仍有殘留，并且圖像全體可見微弱的直線型偽影。L1-TV 算法得到了最好的重建效果，圖像中的環(huán)狀偽影得以完全去除，并且各組織結(jié)構(gòu)的灰度得到了有效的平滑。

2.2 多數(shù)探測器故障下的圖像重建

本文進一步設(shè)計了條件嚴苛的驗證實驗，即假定探測器的故障個數(shù)為9 個。各項實驗中的參數(shù)取值不變。9 個探測器故障條件下各算法的無噪聲影響和有噪聲影響的重建結(jié)果如圖4所示。圖4（a）、（c）、（e）為在無噪聲影響下用各種范數(shù)進行圖像重建的實驗結(jié)果。由結(jié)果圖像可知，嚴重的環(huán)狀偽影分布于整幅L2范數(shù)的重建圖像，并且在環(huán)狀偽影的起點和終點處存在顯著的直線型偽影。圖像的某些細節(jié)信息被偽影覆蓋，并且高灰度值的環(huán)狀偽影降低了圖像原有信息的對比度。在L1 范數(shù)的重建圖像中，高灰度值的環(huán)狀偽影被有效去除，并且在環(huán)狀偽影的起點和終點處的直線型偽影也被完全去除。圖像的細節(jié)信息得到了較好的保存。然而，在L1 的結(jié)果圖像中，仍有暗黑色環(huán)狀偽影瘢痕殘留。L1-TV 算法得到了最佳實驗效果，偽影被完全去除，圖像細節(jié)保存完好，并且各組織結(jié)構(gòu)均勻一致。鑒于CT 圖像采集過程中會受到電子噪聲的影響，為了使該實驗更接近實際運行機制，在投影數(shù)據(jù)中加入了均值為0、方差值為6.0 的高斯噪聲和發(fā)射光子數(shù)為5.0×105條件下的泊松噪聲。各項實驗中其余參數(shù)取值不變，為控制統(tǒng)計噪聲將循環(huán)次數(shù)k設(shè)定為 200。圖4（b）、（d）、（f）為在 2 種噪聲影響下用各種范數(shù)進行圖像重建的實驗結(jié)果。在L2 范數(shù)的重建圖像中充滿了亮暗并存的環(huán)狀偽影以及統(tǒng)計噪聲，環(huán)狀偽影與物體邊緣交界處可見部分組織結(jié)構(gòu)變形。L1 范數(shù)重建充分去除了環(huán)狀偽影的明亮區(qū)域和部分統(tǒng)計噪聲。L2 范數(shù)重建中的部分組織結(jié)構(gòu)變形得以恢復(fù)。L1-TV 算法有效去除了環(huán)狀偽影和統(tǒng)計噪聲，然而由于為達到去除統(tǒng)計噪聲目的而進行的過度平滑使得圖像微小細節(jié)結(jié)構(gòu)有所消失。

圖3 不同數(shù)量探測器故障條件下各算法的實驗結(jié)果

圖4 9 個探測器故障條件下各算法的無噪聲影響和有噪聲影響的重建結(jié)果

2.3 臨床CT圖像的重建

為進一步驗證算法有效性，選取了臨床中常用的上腹部CT 圖像進行了重建實驗。該圖像的像素數(shù)為512×512，像素間距為0.625 mm，因此圖像的采集范圍為32 cm×32 cm。各像素的取值為對應(yīng)區(qū)域的X 線衰減系數(shù)值，像素值分布范圍為[0，0.387]。采集的投影數(shù)據(jù)像素數(shù)為512×512。各實驗參數(shù)設(shè)置如下：α0=0.002，?=50，L1-TV 實驗中 β =8 600（β 值越大，圖像的平滑效果越?。?。含3 個探測器故障的上腹部CT 斷層圖像的重建結(jié)果如圖5 所示，實際圖像的重建效果與上述數(shù)字模型重建效果保持一致。L2 范數(shù)重建結(jié)果中有明顯的環(huán)狀偽影出現(xiàn)，而L1 范數(shù)重建有效去除了環(huán)狀偽影的明亮區(qū)域，L1-TV 重建得到了良好質(zhì)量的重建圖像。

圖5 含3 個探測器故障的上腹部CT 斷層圖像的重建結(jié)果

實驗結(jié)果表明，在單個及多數(shù)探測器故障條件下L2 范數(shù)重建不能消除圖像中的環(huán)狀偽影，而L1 范數(shù)重建能有效去除高灰度值的環(huán)狀偽影，L1-TV 重建算法能夠得到良好的圖像重建效果。在受高斯噪聲和泊松噪聲影響的條件下，L1 算法依然展現(xiàn)出對環(huán)狀偽影去除的有效性。深究3 種算法的不同發(fā)現(xiàn)，L1 算法較L2 算法具有顯著優(yōu)越性的原因在于圖像更新過程具備異常值濾除程序。式（16）中，λ 的取值被限定在的取值被限定在一定合理區(qū)間，該過程即決定了若投影數(shù)據(jù)bi存在異常值則會被濾除的運行機制。而L1-TV 算法較L1 算法具有優(yōu)越性，原因在于全變分懲罰項能夠?qū)D像起到有效的平滑和去噪作用。

3 結(jié) 語

本文研究了在醫(yī)用CT 設(shè)備的核心部件X 線探測器部分損壞的條件下，進行具有容錯機制的圖像重建算法的開發(fā)問題。以經(jīng)典Shepp-Logan 數(shù)字模型作為重建對象，設(shè)計了單個及多個X 線探測器出現(xiàn)故障情況下的重建實驗。此外，還討論了在電子噪聲影響下的重建情況以及實際的上腹部CT 圖像的重建情況。分別對比了傳統(tǒng)L2 范數(shù)重建算法和本文提出的L1范數(shù)重建和L1-TV 重建算法。實驗結(jié)果表明，L2 范數(shù)重建不具備去除探測器故障引起的環(huán)狀偽影的功能，其重建效果接近于當前醫(yī)用商業(yè)CT 所搭載的FBP 算法的重建效果。而L1 范數(shù)重建算法較傳統(tǒng)L2 范數(shù)重建在去除探測器故障引起的環(huán)狀偽影方面具有顯著優(yōu)越性，能有效去除CT 圖像中光亮的環(huán)狀偽影，然而該算法在去除暗黑色環(huán)狀偽影方面還存在不足。增強算法L1-TV 算法有效彌補了L1 算法的不足，具備完全去除環(huán)狀偽影和有效保護圖像細節(jié)和對比度的能力。本文提出的探測器故障下醫(yī)用CT 圖像重建的容錯算法具有較強的容錯能力，能夠?qū)崿F(xiàn)多數(shù)探測器同時出現(xiàn)故障的情況下的精準圖像重建。因此，該技術(shù)為我國各臨床機構(gòu)中已經(jīng)進入老化和故障期的CT 器械能夠繼續(xù)服役提供了有效的技術(shù)支撐。