諶湘倩 馬紹惠 須文波

摘 要： 針對(duì)計(jì)算機(jī)斷層掃描（CT）重建過(guò)程中統(tǒng)計(jì)方法計(jì)算時(shí)間較長(zhǎng)的問(wèn)題，提出一種利用有序子集加速拆分算法的三維CT圖像重建方法。該方法充分利用線(xiàn)性約束凸優(yōu)化問(wèn)題的增廣拉格朗日（AL）方法在較弱條件下的收斂速度快的優(yōu)勢(shì)；同時(shí)針對(duì)內(nèi)部最小二乘問(wèn)題，使用AL方法的線(xiàn)性變形求解加權(quán)正則化最小二乘問(wèn)題，該方法使用可分離二次型代理函數(shù)代替縮放增廣拉格朗日中的二次型AL懲罰項(xiàng)，得到一種簡(jiǎn)單有序子集（OS）加速型拆分算法（OS?ASA），避免了繁瑣的參數(shù)調(diào)整，可快速收斂。實(shí)驗(yàn)結(jié)果表明，該文算法顯著加快了CT圖像重建的收斂速度，當(dāng)使用子集較多時(shí)，CT圖像重建可以減少OS偽影。

關(guān)鍵詞： 三維圖像重建； 計(jì)算機(jī)斷層掃描； 拉格朗日； 有序子集； 最小二乘

中圖分類(lèi)號(hào)： TN919?34； TP391 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 1004?373X（2016）06?0104?04

Three?dimensional CT image reconstruction based on accelerated splitting algorithm

of ordered subsets

CHEN Xiangqian1， MA Shaohui1， XU Wenbo2

（1. Department of Computer Science and Technology， Henan Institute of Technology， Xinxiang 453002， China；

2. College of Information Engineering， Jiangnan University， Wuxi 210034， China）

Abstract： Since the computing time of statistical method in CT （computed tomography） reconstruction process is long， a three?dimensional CT image reconstruction method based on accelerated splitting algorithm of ordered subsets （OS） is proposed. The method takes full advantage of the fast convergence speed when the augmented Lagrangian （AL） method of the linear constraint convex optimization is in weak condition. The weighted and regularized least square problems are solved with linear variant of AL method. The separable quadratic surrogate function is used in this method to replace the quadratic AL penalty term scaled in the AL to obtain a simple accelerated splitting algorithm of ordered?subsets （OS?ASA）， which can avoid the tedious parameter tuning and speed up the convergence rate. The experimental results show that the algorithm significantly accelerates the convergence speed of CT image reconstruction. The OS artifacts can be reduced by CT image reconstruction when the more subsets are used.

Keywords： three?dimensional image reconstruction； computed tomography； Lagrangian； ordered subset； least square

0 引 言

由于獲取較低劑量X射線(xiàn)的CT掃描具有保持圖像質(zhì)量的潛力，圖像重建的統(tǒng)計(jì)方法已廣泛應(yīng)用于計(jì)算機(jī)斷層掃描（Computed Tomography，CT）[1]，然而，統(tǒng)計(jì)方法計(jì)算時(shí)間長(zhǎng)是實(shí)踐應(yīng)用中的最大瓶頸。為了加速統(tǒng)計(jì)方法，已經(jīng)有很多學(xué)者研究其優(yōu)化技術(shù)，例如，文獻(xiàn)[2]的增廣拉格朗日（Augmented Lagrangian，AL）方法（包括交替方向變型）使用兩個(gè)拆分求解正則化逆問(wèn)題，AL方法通過(guò)引入輔助變量可以分離非光滑[l1]正則項(xiàng)，得到簡(jiǎn)單懲罰型最小二乘內(nèi)部問(wèn)題，再使用快速傅里葉變換（FFT）算法和軟閾值[l1]范數(shù)的近端映射能夠有效解決該問(wèn)題。然而，類(lèi)似于X射線(xiàn)CT圖像重建的應(yīng)用程序中，統(tǒng)計(jì)加權(quán)的巨大動(dòng)態(tài)范圍造成了Hessian矩陣高度平移，導(dǎo)致內(nèi)部最小二乘問(wèn)題非常具有挑戰(zhàn)性[1]。為了解決該問(wèn)題，文獻(xiàn)[3]引入了分離平移變量和加權(quán)二次數(shù)據(jù)擬合項(xiàng)的平移不變量，使得預(yù)處理共軛梯度（Preconditioned Conjugate Gradient，PCG）方法能有效解決寬松條件下的最小二乘問(wèn)題。實(shí)驗(yàn)結(jié)果表明，該方法顯著加速了2?D CT[3]。然而，在具有錐形束幾何形狀的3?D CT中，為內(nèi)部最小二乘問(wèn)題構(gòu)建預(yù)處理器更加困難。有序子集（Ordered?subsets，OS）算法[4]是一種具有對(duì)角預(yù)處理器的一階方法，使用較小的步長(zhǎng)，但非常簡(jiǎn)單，適用于3?D CT。通過(guò)集合投影為滿(mǎn)足“子集平衡條件”的M個(gè)有序子集，并使用M個(gè)子集梯度增量更新圖像，與標(biāo)準(zhǔn)梯度下降方法的每次外部迭代的圖像更新一樣多，從而導(dǎo)致標(biāo)準(zhǔn)梯度下降法早期迭代M次加速。當(dāng)在約束條件下隨機(jī)選擇子集使子集梯度無(wú)偏置且具有有限方差時(shí)，也可以稱(chēng)為隨機(jī)梯度方法[5]。當(dāng)M增加時(shí)，快速OS算法有“更大的”限制周期，并在重建圖像中出現(xiàn)偽影，但它可以通過(guò)使用松弛動(dòng)量減小，即增長(zhǎng)的對(duì)角優(yōu)化器（或等價(jià)地，遞減步長(zhǎng)），以較慢的收斂速度為代價(jià)[6]去除偽影。文獻(xiàn)[7]的研究表明，加速近端梯度方法在梯度誤差和近側(cè)映射計(jì)算中對(duì)誤差更敏感。

本文使用線(xiàn)性AL方法（Linear Augmented Lagrangian Method，LALM）求解正則化（加權(quán)）最小二乘問(wèn)題，優(yōu)化了二次AL懲罰項(xiàng)，代替了平滑數(shù)據(jù)擬合項(xiàng)，在縮放增廣拉格朗日中，使用固定對(duì)角線(xiàn)優(yōu)化器，產(chǎn)生更加簡(jiǎn)單的OS加速型拆分算法（OS?ASA）。為了進(jìn)一步加速，使用二階遞歸系統(tǒng)分析設(shè)計(jì)了確定性向下連續(xù)化方法，避免了繁瑣的參數(shù)調(diào)整，可快速收斂。實(shí)驗(yàn)結(jié)果表明，提出的算法在早期迭代中顯著加速了X射線(xiàn)CT圖像重建。

1 本文提出的OS?加速拆分算法

首先考慮一個(gè)普通復(fù)合凸優(yōu)化問(wèn)題：

[x∈argminx{g（Ax）+h（x）}] （1）

式（1）的等效約束最小化問(wèn)題為：

[（x，u）∈argminx，u{g（u）+h（x）} s.t. u=Ax] （2）

式中：[g]和[h]是封閉且恰當(dāng)?shù)耐购瘮?shù)；CT中，[A]表示系統(tǒng)（投影）矩陣；[x]表示待構(gòu)建圖像；[g]表示一個(gè)加權(quán)二次數(shù)據(jù)擬合項(xiàng)；[h]表示邊緣保持正則項(xiàng)。求解約束最小化問(wèn)題（2）的一種方式是使用AL方法（交替方向），交替最小化縮放后的增廣拉格朗日：

[LA（x，u，d；ρ）?g（u）+h（x）+ρ2Ax-u-d22] （3）

相對(duì)于[x]和[u]，按照[d]的梯度增加，產(chǎn)生下列AL迭代[2]：

[x（k+1）∈argminxh（x）+ρ2Ax-u（k）-d（k）22u（k+1）∈argminug（u）+ρ2Ax（k+1）-u-d（k）22d（k+1）=d（k）-Ax（k+1）+u（k+1）] （4）

式中：[d]為拆分變量[u]縮放后的拉格朗日乘數(shù)；[ρ>0]為對(duì)應(yīng)的AL懲罰參數(shù)。

線(xiàn)性AL方法（LALM）[8]在式（4）的x更新中替代二次AL懲罰項(xiàng)：

[θk（x）?ρ2Ax-u（k）-d（k）22] （5）

可將二次代理（SQS）函數(shù)分離為：[θk（x；x（k））?θk（x（k））+?θk（x（k）），x-x（k）+ρL2x-x（k）22 =ρ2tx-（x（k）-tA′（Ax（k）-u（k）-d（k）））22+ （constant independent of x）] （6）

式中，[L>A22=λmax（A′A）]確保[LI-A′A>0]和[t?1L]，該函數(shù)滿(mǎn)足“優(yōu)化”條件：

[θk（x；x）≥θk（x）， ?x，x∈Dom θkθk（x；x）=θk（x）， ?x∈Dom θk] （7）

可很容易地將[L]泛化為對(duì)稱(chēng)半正定矩陣[S+]，例如，基于OS的算法中所用的對(duì)角矩陣[7][DL]，仍能確保式（7）。當(dāng)[S+=A′A]時(shí)，LALM則為標(biāo)準(zhǔn)AL方法，優(yōu)化對(duì)角矩陣產(chǎn)生更加簡(jiǎn)單的[x]?更新，對(duì)應(yīng)的LALM迭代如下[8]：

[x（k+1）∈argminx?k（x）?h（x）+θk（x；x（k））u（k+1）∈argminug（u）+ρ2Ax（k+1）-u-d（k）22d（k+1）=d（k）-Ax（k+1）+u（k+1）] （8）

將[g]限制為二次損失函數(shù)，即[g（u）?（12）y-u22]，則最小化問(wèn)題（1）變?yōu)檎齽t化最小二乘問(wèn)題：

[x∈argminxΨ（x）?12y-Ax22+h（x）] （9）

令[L（x）?g（Ax）]為式（9）的二次數(shù)據(jù)擬合項(xiàng)。假設(shè)[L]適合利用OS加速，即[L]可分解為[M]個(gè)較小的滿(mǎn)足“子集平衡條件”[7]的二次型函數(shù)[L1，L2，…，LM]：

[?L（x）≈ML1（x）≈ML2（x）≈…≈MLM（x）] （10）

有助于子集梯度近似[L]的梯度。由于[g]為二次型，其近端映射為線(xiàn)性，故簡(jiǎn)單閉合形式解可表示為：

[u（k+1）=ρρ+1（Ax（k+1）-d（k））+1ρ+1y] （11）

[d]的更新如下：

[u（k+1）+ρd（k+1）=y] （12）

若初始化[d]為[d（0）=ρ-1（y-u（0））]。令[u?u-y]為拆分剩余差，則簡(jiǎn)化的LALM迭代如下所示（求x的迭代）：

[s（k+1）=A′（ρ（Ax（k）-y）+（1-ρ）u（k））x（k+1）∈prox（ρ-1t）h（x（k）-（ρ-1t）s（k+1））u（k+1）=ρρ+1（Ax（k+1）-y）+1ρ+1u（k）] （13）

每次迭代式（13），其計(jì)算復(fù)雜度減少，為一個(gè)[A]的乘法、一個(gè)[A′]的乘法和一個(gè)[h]的近端映射，可非迭代求解或不使用[A]，[A′]迭代求解。由于[L]的梯度為[A′（Ax-y）]，令[g?A′u]（拆分剩余的逆投影）表示拆分梯度，式（13）可重寫(xiě)為：

[s（k+1）=ρ?L（x（k））+（1-ρ）g（k））x（k+1）∈prox（ρ-1t）h（x（k）-（ρ-1t）s（k+1））g（k+1）=ρρ+1?L（x（k+1））+1ρ+1g（k）] （14）

式（14）稱(chēng)為基于梯度的LALM，因?yàn)閮H[L]的梯度用于執(zhí)行更新，每次迭代式（14）的計(jì)算復(fù)雜度變?yōu)閇L]的梯度評(píng)估和[h]的近端映射。

基于梯度的LALM為正則化最小二乘成本函數(shù)[ψ]的廣義近端梯度下降法，其中，步長(zhǎng)為[ρ-1t]、搜索方向?yàn)閇s（k+1）]，即[L]的梯度和拆分梯度的加權(quán)平均。[ρ]越小，產(chǎn)生的步長(zhǎng)越大，當(dāng)[ρ=1]時(shí)，式（14）變?yōu)榻颂荻确椒ɑ虻湛s/閾值算法（ISTA）。即利用LALM，通過(guò)減小[ρ]可任意增加近端梯度方法的步長(zhǎng)。假設(shè)所有更新均準(zhǔn)確，即對(duì)于所有[k]，[εk=0]。由式（12）得：

[-ρd（k）=u（k）-y→Ax-y=z，k→∞（-ρd（0））-z=u（0）-Ax]

合理的初始化，即[u（0）=Ax（0）]，因此，[g（0）=?L（x（0））]，可重寫(xiě)為[ρ]的函數(shù)：

[C（ρ）=x（0）-x22ρ-1t+A（x（0）-x）22ρ] （15）

當(dāng)：

[ρopt=A（x（0）-x）2x（0）-x2≤1] （16）

該常量獲得最小值；當(dāng)[L?A22]，[ρopt?1]且滿(mǎn)足第一項(xiàng)控制[C]，使[ρopt<ρ≤1]。由于[ρρopt?1]，原始對(duì)偶差的上界變?yōu)椋?/p>

[Ω（zk，x）-Ω（z，xk）≤C2k≈L2x（0）-x22ρ-1k] （17）

即相比近端梯度方法（[ρ=1]），利用[ρ-1]的因子提升了式（14）的收斂速度（邊界），[ρopt<ρ≤1]。

最后，由于式（14）僅要求[L]的梯度執(zhí)行更新，從而產(chǎn)生本文提出的OS?加速LALM（OS?ASA）：

[s（k，m+1）=ρM?Lm（x（k，m））+（1-ρ）g（k，m））x（k，m+1）∈prox（ρ-1t）h（x（k，m）-（ρ-1t）s（k，m+1））g（k，m+1）=ρρ+1M?Lm+1（x（k，m+1））+1ρ+1g（k，m）c（k，m+1）=c（k+1）=c（k+1，1）， c∈s，x，g] （18）

當(dāng)[LM+1=L1]時(shí)，類(lèi)似于典型OS算法，當(dāng)[M=1]時(shí)該算法收斂，當(dāng)[M>1]時(shí)，無(wú)法確保收斂。

2 OS?ASA重建CT圖像

X?ray CT圖像重建問(wèn)題是一種約束正則化加權(quán)最小二乘問(wèn)題，使用下式置換求解式（14）：

[A←W12Ay←W12yh←R+lΩ] （19）

式中：[lΩ]指凸集[C]的特征函數(shù)，因此，式（14）中內(nèi)部極小化問(wèn)題轉(zhuǎn)變?yōu)橐环N約束去噪問(wèn)題。使用快速迭代收縮/閾值算法（FISTA）[9]求解該內(nèi)部約束去噪問(wèn)題，以上一次更新作為熱啟動(dòng)的開(kāi)始。

通常，F(xiàn)ISTA迭代次數(shù)越多，本文算法收斂速度越快，但是，當(dāng)[n]較大時(shí)，迭代內(nèi)更新的開(kāi)銷(xiāo)不可忽略。在典型X?ray CT圖像重建問(wèn)題中，優(yōu)化一般較容易，因?yàn)榻y(tǒng)計(jì)學(xué)加權(quán)[W]的巨大動(dòng)態(tài)范圍。因此，大部分情況下[t?1]，極大地消除了約束去噪問(wèn)題中的正則化功率。實(shí)踐中，約束去噪問(wèn)題可解決，一至兩次迭代后，正則化問(wèn)題就可滿(mǎn)足容差條件。

本文使用具有對(duì)角Hessian的SQS函數(shù)[DL?diag{A′WA1}]優(yōu)化加權(quán)二次型數(shù)據(jù)擬合項(xiàng)[8][L]，還使用具有Hessian [DL]的SQS函數(shù)優(yōu)化縮放增廣拉格朗日中的二次型懲罰參數(shù)，并使用具有位反轉(zhuǎn)順序的子集梯度增量更新圖像。

3 實(shí)驗(yàn)結(jié)果

針對(duì)具有各種幾何形狀的真實(shí)CT掃描的3?D X?ray CT圖像，對(duì)幾種基于OS算法的重建結(jié)果進(jìn)行比較，算法包括：OS?SQS?M，具有[M]個(gè)子集的標(biāo)準(zhǔn)OS算法[7]；OS?Nes05?M：基于[M]個(gè)子集的Nesterov快速梯度方法，即OS+momentum算法[10]；OS?ASA?M?[ρ]?[n]，利用[M]個(gè)子集和[n]次FISTA迭代，提出固定的AL懲罰參數(shù)用于求解內(nèi)部約束去噪問(wèn)題。

OS?SQS是層析重建的標(biāo)準(zhǔn)迭代方法，OS?Nes05是使用Nesterov技術(shù)的X?ray CT圖像快速重建方法。不同于其他幾種基于OS算法，由于迭代內(nèi)更新，本文算法具有額外開(kāi)銷(xiāo)，然而，當(dāng)[n=1]時(shí)，即約束去噪問(wèn)題具有單梯度降落，每次迭代一個(gè)向前/向后投影對(duì)和[M]個(gè)正則化梯度評(píng)估，故幾種算法具有相同的計(jì)算復(fù)雜度。

本文從肩區(qū)域螺旋CT掃描重建一幅512[×]512[×]109圖像，其中，正弦圖大小為888[×]32[×]7 146，間隙為0.5，子集最大數(shù)約為40。圖1為初始FBP圖像中央橫斷平面的裁剪圖像、參考重建以及使用本文算法（OS?ASA?40?c?1）在第30次迭代時(shí)的重建圖像。圖（a）為肩掃描從初始FBP圖像的中央軸平面裁剪的圖像[x（0）]（顯示從800～1 200 HU）；圖（b）為參考重建x*；圖（c）為使用提出的算法（OS?ASA?40?c?1）在第30次迭代處的重建圖像[x（30）]。圖2表示差分圖像，即[x（30）-x*]，針對(duì)幾種基于OS的算法，標(biāo)準(zhǔn)OS算法（20和40子集）在差分圖像中展示了可見(jiàn)條紋偽影和結(jié)構(gòu)化高頻噪聲。圖2中肩掃描：使用基于OS的算法，從[x（30）]-x*的中央軸平面裁剪的差分圖像（顯示從-30～30 HU）。當(dāng)M=20時(shí)，OS+Momentum算法的差分圖像更加結(jié)構(gòu)化和不均勻，但從整體上看，OS+Momentum算法與本文算法的差分圖像相似。當(dāng)M=40時(shí)，本文算法的差分圖像仍然均勻，而OS+Momentum算法的差分圖像帶有類(lèi)似噪聲的OS偽影，當(dāng)M增加時(shí)，M=40時(shí)，OS+Momentum算法并沒(méi)有偽影惡化。顯然，OS?ASA比OS+Momentum方法具有更好的梯度容差。

圖1 參考圖像和重建圖像

圖2 差分圖像

針對(duì)X?ray CT 圖像重建，圖3給出了使用FISTA求解內(nèi)部約束去噪問(wèn)題。圖3中肩掃描為本文提出求解內(nèi)約束的去噪問(wèn)題。其中，n取1，2或5以重建圖像x（k）和參考圖像之間x*的均方差作為迭代的函數(shù)。從圖3可看出，當(dāng)使用多次FISTA迭代求解內(nèi)部約束去噪問(wèn)題時(shí)，收斂速度沒(méi)有明顯提升。因此，只需1次FISTA迭代（n=1），每個(gè)子集更新即可滿(mǎn)足快速精確的X?ray CT圖像重建。

圖3 內(nèi)部約束去噪

4 結(jié) 語(yǔ)

增廣拉格朗日（AL）方法和有序子集（OS）是兩種分別使用分解和近似的加速優(yōu)化算法，通過(guò)考慮AL方法的線(xiàn)性變形，結(jié)合這兩種技術(shù)，本文提出了一種快速OS?加速型拆分算法（OS?ASA），解決了正則化（加權(quán)）最小二乘問(wèn)題。利用OS?ASA對(duì)X射線(xiàn)計(jì)算機(jī)斷層掃描（CT）圖像進(jìn)行重建，并將其與幾種先進(jìn)的OS算法進(jìn)行比較，使用具有不同幾何形狀的真實(shí)CT掃描圖像作為數(shù)據(jù)集。實(shí)驗(yàn)結(jié)果表明，OS?ASA具有較快的收斂速度和優(yōu)異梯度容差。今后的研究工作是確定性向下連續(xù)化方法的收斂速度分析，以及M大于1時(shí)OS?ASA的更嚴(yán)格收斂分析。

參考文獻(xiàn)

[1] 趙杰，龔碩然，王龍.腦部CT圖像的三維重建及紋理特征研究[J].激光雜志，2014，35（6）：62?65.

[2] 石英，肖金生，劉春曉，等.配氣凸輪優(yōu)化設(shè)計(jì)的懲罰函數(shù)法和增廣拉格朗日乘子法[J].武漢理工大學(xué)學(xué)報(bào)（交通科學(xué)與工程版），2002，26（3）：365?368.

[3] RAMANI S， FESSLER J A. A splitting?based iterative algorithm for accelerated statistical X?ray CT reconstruction [J]. IEEE transactions on medical imaging， 2012， 31（3）： 677?688.

[4] 林少春.基于有序子集的懲罰加權(quán)最小二乘低劑量CT優(yōu)質(zhì)重建算法研究[J].中國(guó)組織工程研究與臨床康復(fù)，2008，12（4）：679?682.

[5] 張華軍，趙金，羅慧，等.基于梯度投影法與隨機(jī)優(yōu)化算法的約束優(yōu)化方法[J].控制與決策，2014，27（10）：1777?1782.

[6] 何玲君.基于最大似然和罰似然估計(jì)的CT統(tǒng)計(jì)重建算法研究[D].太原：中北大學(xué)，2011.

[7] SCHMIDT M， ROUX N L， BACH F. Convergence rates of inexact proximal?gradient methods for convex optimization [J]. Advances in neural information processing systems， 2011， 24： 1458?1466.

[8] LIN Z， LIU R， SU Z. Linearized alternating direction method with adaptive penalty for low?rank representation [J]. Advances in neural information processing systems， 2011， 24： 612?620.

[9] 孫念，胡炳樑，王爽，等.基于FISTA算法的編碼孔徑光譜圖像壓縮與復(fù)原系統(tǒng)[J].紅外與激光工程，2014，19（1）：238?242.

[10] KIM D， RAMANI S， FESSLER J A. Accelerating X?ray CT ordered subsets image reconstruction with Nesterovs first?order methods [C]// Proceedings of the 12th International Meeting on Fully 3?D Image Reconstruction in Radiology and Nuclear Medicine. [S.l.： s.n.]， 2013： 22?25.