基于Moreau-包絡的近似平滑迭代磁共振圖像重建算法

劉曉暉，路利軍，馮前進，陳武凡

(南方醫科大學 生物醫學工程學院，廣州 510515)(*通信作者電子郵箱chenwf@fimmu.com)

0 引言

磁共振成像(Magnetic Resonance Imaging, MRI)因其成像中不使用電離輻射、可以更靈活地選擇成像參數、對軟組織結構及代謝功能的完美表達以及對許多特定參數更加敏感等優點成為臨床應用中性能最卓越的醫學成像方式之一[1-2]。然而，MRI特殊的成像機制導致的相對較慢的成像速度阻礙了MRI進一步的發展和應用：磁共振(MR)圖像對于運動非常敏感，因而成像時間長會產生更多的運動偽影，這種情況常常導致需要對病人進行鎮靜或麻醉；較低的空間分辨率限制了MR對某些和呼吸相關的部位的成像，如腹部成像和心臟成像；較長的掃描時間在增加成像成本的同時也會限制MRI的適用病人范圍[3]。出于這些原因，許多研究人員在尋求在減少采樣數據的同時又不降低圖像質量的成像方法。

從20世紀70年代MR還未被廣泛地應用于臨床以來，MRI的加速問題就被廣泛的探討和研究。MRI中成像數據為Fourier變換空間，通常稱為k空間。MRI的采樣時間主要由k空間的分段采樣模式所決定，在給定一個采樣序列以后，k空間中的分段數就決定了MR圖像的信號采集時間。基于這個原因，MRI的加速通常是采用減少k空間中的分段數，即對k空間進行降采樣。這種方法利用了k空間的冗余特性，即k空間中的各個點是由位于成像物體不同空間位置的信號疊加產生的[4]，主要代表技術是并行MRI(parallel MRI, pMRI)技術[5-8]。pMRI技術利用多組線圈同時接收MR信號，利用線圈的空間敏感度差異來代替部分梯度編碼，從而減少梯度編碼步數，節省掃描時間。pMRI無需修改成像序列，在保持掃描視野和空間分辨率不變的情況下減少成像時間，因此為現代商業MRI系統廣泛采用，成為目前的常規掃描方法。然而，由于欠采樣導致的混疊偽影，以及重建算法對噪聲的放大效果，pMRI方法在實際臨床應用中只能使用小加速倍數以保證理想的重建質量。

近年來在信號處理領域發展起來的壓縮感知(Compressed Sensing, CS)[9-12]理論為快速MRI技術提供了一種新的可能。CS理論認為，一個本身稀疏的或者可壓縮的未知信號，選擇一組合適的感知信號(即感知矩陣)，在感知信號個數遠小于原始未知信號維數的情況下，通過求解一個非線性最優化問題，可以精確恢復出原始信號。在MRI中：1) MR圖像被認為是稀疏(如血管MR圖像、三維MR圖像相鄰圖像層之間以及動態MR圖像在時域上的圖像序列)或是可壓縮的(一幅MR圖像在某個變換域中的系數矩陣具有可壓縮性)；2) MRI的感知信號空間是Fourier空間，其成像方式決定了可以有選擇的掃描k空間中的一部分數據。上述兩個條件決定了CS可以應用于快速MRI中，這種快速MRI技術一般稱為CS-MRI。

CS-MRI通常需要求解一個包含稀疏正則項的非線性凸優化問題。為了更好地表達MR圖像的稀疏性，CS-MRI的數學模型中往往包含兩個稀疏正則項，即優化問題中包含三個函數：數據保真項、1范數稀疏正則項和MR圖像的總變分(Total Variation, TV)正則項；而經典的近似迭代類算法求解的是包含兩個函數(通常為數據保真項和一個正則項)的最優化問題，所以在使用近似點迭代類算法的CS-MRI方法中，通常的做法是將初始優化問題根據正則項分解為兩個優化子問題，分別對這兩個子問題用近似迭代算法進行求解得到相應的解，再對求得的兩個解取平均值得到最終解。這類算法由于每個子問題的輸入變量相同可以稱為多正則項優化問題的并行求解方法。

近似迭代算法主要是針對優化問題中存在的非平滑函數，如1范數正則項、核范數以及TV項等，這些非平滑函數不能直接求導，因此在優化算法中需要求取它們的次梯度，在近似迭代類算法中往往使用近似算子來求解這個次梯度，近似算子本身即為一個最優化問題的求解過程，通過對輸入信號x0迭代地施加近似算子即可求解原始優化問題的最優解。在前面提到的多正則項優化問題的并行求解算法中，兩個子問題的解在同一個迭代步驟中是互相獨立的，即互相不以對方為輸入，這類算法在兩個子問題的求解過程中沒有任何優化步驟。針對這個問題，本文提出一種基于Moreau包絡的近似平滑迭代算法(Proximal Smoothing Iterative Algorithm, PSIA)，該算法首先將原始優化問題中的一個非平滑正則項的梯度用它的平滑函數的梯度近似表示，然后把數據保真項和平滑后的正則項的線性組合看作一個平滑可導函數，最后利用經典的近似迭代算法對優化問題進行求解。

1 磁共振圖像稀疏重建相關算法

在文獻[4]中，對CS-MRI的理論基礎及數學模型都作了詳細的探討，并證明了在附加一個1范數小波變換約束項和TV約束項的前提下，可以從降采樣的k空間數據中準確地重建出MR圖像，并提出如式(1)中的非約束凸優化問題：

α‖Wx‖1+β‖x‖TV

(1)

其中：x∈RN×N為待重建的MR圖像，N×N為MR圖像尺寸；Fu:RN×N→CM×N(M≤N)，為降采樣Fourier變換，是傳感矩陣作用于Fourier矩陣的結果，即Fu=SF，其中S∈RM×N，F:RN×N→CN×N分別為傳感矩陣和Fourier變換；b∈RM×N為k空間降采樣數據；W:RN×N→RN×N為小波變換；‖·‖TV:RN×N→R為MR圖像的TV值；α，β>0為正則項權重系數。式(1)中的幾個范數定義分別為式(2)～(4)：

(2)

(3)

(4)

式中D1、D2為分離梯度算子，分別表示對圖像x作水平、垂直兩個方向滿足Neumann邊界條件的前向有限差分：

(5)

并滿足下面的Neumann邊界條件：

(6)

式中n1、n2滿足n1,n2∈[1,N]。對于問題式(1)的求解，最主要的困難來自于第二項小波變換1范數和第三項TV項，這兩項都不能直接求取微分或導數，因而不能對問題式(1)直接使用梯度法求最小值。針對這個問題，文獻[4]提出了共軛梯度(Conjugate Gradient, CG)下降算法來求解問題式(1)，但是這種算法的缺點是收斂速度慢，對于一幅512×512像素大小的圖像，CG重建出視覺效果比較好的圖像所需要的時間是2 min～3 min，即使重建一幅256×256像素大小的圖像，也需要十幾秒，當對三維物體成像時，所需時間就要以小時計，這顯然不能滿足實際應用的需要。

為了提高重建速度，考慮到問題式(1)的正則項是兩個截然不同的范數表達式，數學優化算法中的分離算法被用來解決式(1)這樣的多個正則項優化問題。分離算法主要包括算子分離算法和變量分離算法。算子分離算法的核心問題是尋找最小化問題所對應的最大單調算子為零的解，經典的算子分離算法有forward-backward算法[13]、Douglas-Rachford算法[14]和Peaceman-Rachford算法[15]，它們都是求解目標函數為兩個函數和的最小化問題的算法。與之相對應的，變量分離算法則是使用變方向法(Alternating Direction Method, ADM)求解最小化問題的增廣拉格朗日形式。變量分離算法是Douglas-Rachford算法和Peaceman-Rachford算法的優化變形，最早在20世紀70年代被提出用來解決偏微分方程(Partial Differential Equation, PDE)問題[16]。近來，在文獻[17]中提出了變量分離算法可以用來解1范數和TV正則優化問題。文獻[18]和文獻[19]分別提出算子分離類方法TV1范數壓縮MRI(TVl1Compressed MR Imaging, TVCMRI)算法與變量分離類方法部分k空間重建算法(Reconstruction from Partial k space algorithm, RecPF)來求解問題式(1)。這兩種算法都比算法CG要快，重建一幅像素為256×256的磁共振圖像所需時間為2 s～3 s，在保證重建準確度的前提下，大幅提高了重建速度。然而由于分離算法的特點，這兩種算法人為添加了一定數量的輔助變量，使得算法代碼的復雜性較大，從而算法的重現具有一定的難度；另外，這兩種算法中，對于兩個正則項求導過程中的近似處理都是簡單地加上一個非負二次項，對于算法的簡化沒有貢獻。

針對以上問題，近似迭代算法被引進非平滑正則約束的優化問題中。近似迭代算法的基本原理是給出非平滑函數的近似算子[20]，這個過程就是一個求取子問題的最小化的過程。對給定一個初始點的原問題式(1)，迭代地算出它的近似點，即可求解得出理論上的最優解。這類算法和前述幾類算法不同的地方在于：直接對非平滑項的近似算子進行操作，而不是對非平滑項求取梯度或是微分[21]。兩種最著名的近似迭代算法是迭代收縮閾值算法(Iterative Shrinkage Thresholding Algorithm, ISTA)[22]和它的加速形式——快速迭代收縮閾值算法(Fast Iterative Shrinkage Thresholding Algorithm, FISTA)[23]，這兩種算法被提出用來求解兩個函數和的最小化問題。在快速復合分離算法(Fast Composite Splitting Algorithm, FCSA)[24-25]中將問題式(1)分解為兩個子問題，即數據保真項與小波變換1范數的和的最小化問題以及數據保真項與TV項的和的最小化問題，分別對這兩個問題使用快速迭代收縮閾值算法，將求得的解取平均，即認為所得的值為問題式(1)的解。在FCSA中，兩個最小化子問題的權重被簡單地取值為0.5，然而在實際應用中，這兩個子問題的權重則根據實際所要重建的磁共振圖像的內容而有所不同，這就向算法的魯棒性提出了挑戰。本文提出的重建算法則避免了這個困擾，即不對問題式(1)進行子問題分離，利用近似平滑方法的數學特性求解。

2 快速近似平滑迭代算法

2.1 相關快速重建算法

在給出本文算法之前，先討論目標函數為兩個凸函數線性組合的優化問題的求解。近似迭代類算法FISTA被用來求解式(15)所示兩個函數和的最小化問題：

Γ(x)=p(x)+q(x):x∈RN×N

(7)

其中:p(x)為一個凸的Lp平滑函數;q(x)為一個不可微的連續凸函數。FISTA提出通過迭代地求解函數q(x)的近似點便可得到問題式(7)的最優解，即：

xk+1=proxq(xk-η▽p(xk));

(8)

其中：xk+1、xk分別為第k+1與第k次迭代求得的解；▽p(x)為函數p(x)的導數；proxq(x)為函數q(x)在x處的近似值函數。

FISTA算法偽代碼為：

步驟1 參數及變量初始化：

參數：正則化參數α,β，迭代計數k=1，近似算子參數

η=1/Lp。

變量：圖像變量x0=0，y1=0，加速步長變量θ1=1。

步驟2 變量迭代：

Fork=1,2,… do

xk=proxq(yk-η▽p(yk))

yk+1=xk+((θk-1)/θk+1)(xk-xk-1)

End for

算法FCSA將問題式(1)分解為兩個子問題，分別是小波1范數正則約束問題和TV正則約束問題，利用FISTA對這兩個問題求解后所得解為x1、x2，再取平均即為問題式(1)的解。本文將提出另外一種利用算法FISTA求解問題式(1)的近似平滑迭代算法(PSIA)。

2.2 近似平滑迭代算法

本文需要求解的問題式(1)是由三個函數線性組合而成，將問題式(1)簡化表達為：

Γ(x)=f(x)+g(x)+h(x):x∈RN×N

(9)

其中f(x)、g(x)和h(x)的表達式分別為：

(10)

如果要對問題式(9)使用算法FISTA進行求解，則對照式(8)可得其變量近似迭代步驟為：

xk+1=proxh(xk-η▽(f+g)(xk));

(11)

式中的proxh(x)在已知x時很容易求得，而導數▽(f+g)(x)由于函數g(x)的非平滑性不能直接求取。為了得到▽(f+g)(x)的表達式，首先要對函數(f+g)(x)進行平滑近似，由于函數f(x)是平滑函數，因此平滑近似主要是對函數g(x)進行。在凸優化理論[26]中，可以通過求解非平滑函數的Moreau包絡來對其進行平滑近似，其定義為：

定義1 Moreau包絡[26]。對于一個非平滑不可微函數g(x):RN×N→R∪{+∞}，它的Moreau包絡為：

(12)

定理1[27]給定函數g(x)為一個下連續的非平滑凸函數且μ>0，則由定義1給出的Moreau包絡為一個1/μ平滑函數，且其導數可由g(x)的近似算子得到：

▽(gμ)(x)=(x-proxg(x))/μ

(13)

其中proxg(x,μ)的表達式為：

(14)

此時式(11)中的導數▽(f+g)(x)中g(x)部分可由式(13)給出。

本文提出求解優化問題式(9)的近似迭代算法為：首先對非平滑函數g(x)使用Moreau包絡進行平滑近似，則得到一個新的優化問題：

Γμ(x)=f(x)+gμ(x)+h(x):x∈RN×N

(15)

在優化問題式(15)中，由于函數f(x)、gμ(x)都是定義域為x∈RN×N的平滑凸函數，因此可以將這兩個函數的和看成一個新的函數K(x)，則式(15)可以表達為：

K(x)=f(x)+gμ(x):x∈RN×N

(16)

對問題式(16)，使用算法FISTA進行求解，其變量迭代步驟為式(11)。

本文算法偽代碼為：

步驟1 參數及變量初始化：

參數：正則化參數α,β，迭代計數k=1，近似算子參數μ=1/Lf,η=1/(Lf+1/μ)。

變量：所求取圖像變量x0=0，中間變量y1=0，加速步長變量：θ1=1。

步驟2 變量迭代求解：

Fork=1，2，… do

xk+1←proxh(yk-η▽(f+gμ)(yk));

yk+1=xk+((θk-1)/θk+1)(xk-xk-1)

End for

迭代步驟1中▽(f+gμ)的計算如下：

(17)

▽(gμ)(x)=(x-proxg(x))/μ=

(sign(Wx).×max{0,Wx-μ/Lf})/μ

(18)

(19)

由此可知上述算法PSIA中的所有計算單元都已知，給定一個初始點，即可迭代地計算出重建圖像。

3 實驗結果與分析

3.1 實驗設置

為了驗證算法的有效性，將本文提出的近似平滑迭代算法(PSIA)與第1章中提到的4種經典稀疏MR重建算法的重建結果進行對比，這4種重建算法分別為：CG、TVCMRI、RecPF和FCSA，選取這些算法進行對比是因為它們求解的目標函數和本文一致，且它們都是經典的磁共振稀疏重建算法。本文算法與對比算法都是在Inter Core i5- 2400、3.10 GHz CPU的PC上進行，編程環境使用Matlab 2015b，每個算法的迭代次數都是50。使用4種圖像數據進行重建結果對比，分別是：經典的Shepp-Logan體模圖像、人腦MR圖像、手部MR圖像，以及腦部血管MR圖像；為了方便對比，所有圖像尺寸都設為256×256像素大小，灰度取值范圍為0～255；實驗降采樣率約為20%，采樣矩陣與實驗圖像數據如圖1所示；在實驗中，對采樣數據添加標準差為σ=0.01的高斯白噪聲以模擬采樣中產生的噪聲。

圖1 采樣矩陣與四種參考實驗圖像

算法中的參數設置為：正則化參數α、β分別設為0.035和0.001，這個取值與4種對比算法中相同，以便于算法效果的比較；近似算子參數的設置根據Lipschitz常數分別設置為：1，0.5。

3.2 實驗結果評估標準

本文主要給出兩大類實驗結果評估標準，分別為定性評估即視覺效果評估，以及定量評估包括信噪比(Signal-to-Noise Ratio, SNR)、相對誤差(Relative Error, RE)和結構相似度指數(Structural SIMilarity, SSIM)[28]，另外還給出了算法復雜度對比(CPU時間)。其中SNR、RE和SSIM的表達式為：

(20)

3.3 實驗結果及討論

四種實驗圖像在五種不同算法(CG、TVCMRI、RecPF、FCSA和PSIA)中的重建結果與重建結果對應的差值圖像如圖2～5所示。在差值圖像的表示中，為了比較的公平性，在每個圖像的結果表示中，將重建結果最差的算法對應的差值圖像灰度值作為灰度直方圖范圍。

圖2 Shepp-Logan體模的重建結果

圖3 腦部MR圖像的重建結果

圖4 手部MR圖像的重建結果

從圖2～5中可以看出：本文提出的重建算法PSIA在重建的視覺效果上比之前的幾種經典算法都要好，具體體現在差值圖像中的邊緣細節的重建以及對噪聲的消除上。算法CG的重建效果最差，不僅沒有很好地消除引進的高斯白噪聲，在圖像的邊緣和解剖學結構上都丟失了很多細節；算法TVCMRI和RecPF的重建結果在視覺效果上相近，且同時在細節的恢復和噪聲的消除上優于CG；算法FCSA比TVCMRI和RecPF得到更好的重建結果，具體體現在消除了更多的噪聲，使重建圖像看起來更“干凈”；本文提出的算法PSIA比起FCSA在去噪方面表現更好，具體表現在差值圖中圖像背景部分的噪聲更小，在圖像細節的恢復上也更好，在視覺上與參考圖像以及真實MR圖像更接近。

四幅圖像在五種重建算法中的SNR與迭代次數的關系如圖6所示，用來比較算法的收斂性。

從圖6可以看出，對于不同的重建圖像，算法CG的效果是最差的，體現為SNR值最低；TVCMRI和RecPF的重建效果相近，這可以從算子分離和變量分離這兩類算法的關系推測出原因：變量分離類算法是算子分離類算法的優化變形；FCSA和本文中提出的PSIA算法的重建效果比前三種算法要好，而本文提出的算法PSIA在四個重建圖像上的重建結果都好過FCSA，這驗證了第1章中的討論，即：簡單地將兩個子問題的權重設為0.5帶來魯棒性問題，導致了雖然都是使用快速迭代收縮閾值算法(FISTA)，PSIA比FCSA的重建效果要更好。

表1～4給出了四種圖像不同算法50次迭代后的重建性能，包括SNR、RE以及SSIM；同時為了比較算法復雜度，給出了每個算法運行的CPU時間。SNR給出了重建圖像與參考圖像的信噪比，SNR的值越大，認為重建算法的重建效果越好；RE給出了重建圖像和參考圖像間的差值，RE的值越小，認為重建算法的重建效果越好；SSIM給出了重建圖像和參考圖像間的相似度測量，SSIM的值越大，認為重建效果越好。

表1 Sheep-Logan體模的重建性能

從表1～4可以看出，在這三個評估參數方面，PSIA都是最好的，同時也可以看出算法CG的結果是最差的，算法TVCMRI和RecPF表現相近，FCSA優于前三種算法，比PSIA表現差，這與圖2～5中的視覺評估以及圖6中的SNR隨迭代次數變化的結果一致。在CPU時間上，算法CG運行時間最長，從TVCMRI開始算法運行時間即加速到十倍數量級，RecPF的運行時間和TVCMRI相當，FCSA比以上兩個算法都要快，這得益于快速迭代收縮閾值算法的使用，PSIA在時間上和FCSA相當，這也與兩者使用的核心算法都是同一種算法所得出的推斷一致。

圖5 腦部血管MR圖像的重建結果

圖6 四種MR圖像重建結果的SNR

算法SNR/dBRE/%SSIMCPU時間/sCG10.1826.910.594311.97TVCMRI14.4213.890.62551.83RecPF15.3713.230.64381.81FCSA18.099.100.78221.36PSIA19.617.640.85281.37

表3 手部MR圖像的重建性能

表4 腦部血管MR圖像的重建性能

4 結語

針對現有的MR圖像稀疏重建算法中不能在一個迭代計算步驟中同時對兩個正則約束項同時施加近似算子的問題，提出了一種基于Moreau-包絡的近似平滑迭代算法，該算法具有以下特點：

1)使用基于近似平滑迭代的算法，使算法結構簡單、易于重現；

2)與同樣基于近似平滑迭代的算法FCSA相比，解決了將原始最小化問題分解為兩個子問題時面臨的權重設置問題；

3)此算法可以應用在不同正則約束稀疏化重建問題，在算法PSIA中，只要給出相對應的正則項的近似算子，即可直接使用該算法。

對仿真體模以及真實MR圖像的實驗結果表明，本文提出的近似迭代算法在重建精度以及算法復雜性方面都有很好的表現，重建效果比當前經典的CG、TVCMRI、RecPF和FCSA都要好，在算法復雜度方面和最快的算法FCSA相當。

下一步的工作主要有兩方面：一是對算法中參數的優化，在本文算法PSIA中，近似算子參數μ、η取的是固定的常數，在優化算法文獻中，提出可以對這類參數進行優化從而得到更好的重建結果；二是提出新的正則約束項并利用算法PSIA，尤其是針對TV約束項，目前已有多種對于總變分的變形以及新的定義，將這些新的約束用在重建中預期將得到更好的重建結果。