諶湘倩 馬紹惠 須文波

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

關鍵詞： 三維圖像重建； 計算機斷層掃描； 拉格朗日； 有序子集； 最小二乘

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

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

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

首先考慮一個普通復合凸優(yōu)化問題：

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

式（1）的等效約束最小化問題為：

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

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

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

相對于[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]為對應的AL懲罰參數(shù)。

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

[θ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ù)滿足“優(yōu)化”條件：

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

可很容易地將[L]泛化為對稱半正定矩陣[S+]，例如，基于OS的算法中所用的對角矩陣[7][DL]，仍能確保式（7）。當[S+=A′A]時，LALM則為標準AL方法，優(yōu)化對角矩陣產(chǎn)生更加簡單的[x]?更新，對應的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]，則最小化問題（1）變?yōu)檎齽t化最小二乘問題：

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

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

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

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

[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]為拆分剩余差，則簡化的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），其計算復雜度減少，為一個[A]的乘法、一個[A′]的乘法和一個[h]的近端映射，可非迭代求解或不使用[A]，[A′]迭代求解。由于[L]的梯度為[A′（Ax-y）]，令[g?A′u]（拆分剩余的逆投影）表示拆分梯度，式（13）可重寫為：

[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）稱為基于梯度的LALM，因為僅[L]的梯度用于執(zhí)行更新，每次迭代式（14）的計算復雜度變?yōu)閇L]的梯度評估和[h]的近端映射。

基于梯度的LALM為正則化最小二乘成本函數(shù)[ψ]的廣義近端梯度下降法，其中，步長為[ρ-1t]、搜索方向為[s（k+1）]，即[L]的梯度和拆分梯度的加權平均。[ρ]越小，產(chǎn)生的步長越大，當[ρ=1]時，式（14）變?yōu)榻颂荻确椒ɑ虻湛s/閾值算法（ISTA）。即利用LALM，通過減小[ρ]可任意增加近端梯度方法的步長。假設所有更新均準確，即對于所有[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））]，可重寫為[ρ]的函數(shù)：

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

當：

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

該常量獲得最小值；當[L?A22]，[ρopt?1]且滿足第一項控制[C]，使[ρopt<ρ≤1]。由于[ρρopt?1]，原始對偶差的上界變?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）

當[LM+1=L1]時，類似于典型OS算法，當[M=1]時該算法收斂，當[M>1]時，無法確保收斂。

2 OS?ASA重建CT圖像

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

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

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

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

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

3 實驗結果

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

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

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

圖1 參考圖像和重建圖像

圖2 差分圖像

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

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

4 結 語

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

參考文獻

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

[2] 石英，肖金生，劉春曉，等.配氣凸輪優(yōu)化設計的懲罰函數(shù)法和增廣拉格朗日乘子法[J].武漢理工大學學報（交通科學與工程版），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] 林少春.基于有序子集的懲罰加權最小二乘低劑量CT優(yōu)質(zhì)重建算法研究[J].中國組織工程研究與臨床康復，2008，12（4）：679?682.

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

[6] 何玲君.基于最大似然和罰似然估計的CT統(tǒng)計重建算法研究[D].太原：中北大學，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算法的編碼孔徑光譜圖像壓縮與復原系統(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.