基于自適應字典學習的動態磁共振并行重建

王 悅 蔣慧敏 汪 洋,2*

1(臺州學院 浙江 臺州 318000)2(河海大學 江蘇 南京 210098)

0 引 言

磁共振成像(MRI)技術具有無創傷、無輻射、分辨率高和可多維成像等優點，被廣泛應用于臨床醫學各個系統，是繼CT以后的又一重要臨床檢測方法。但MR成像速度慢是其一大缺點，尤其是動態磁共振成像(dynamic magnetic resonance imaging，dMRI)，需要在較短時間內獲得高時空分辨率的MRI圖像序列，目前是醫學界的一個難題[1]。過長的掃描時間加上病人的器官運動(如呼吸、吞咽等)，會導致成像模糊，同時也無法滿足動態實時成像和功能成像的需求。在k空間進行降采樣是提高成像速度的一種方法，但如果直接從k空間逆傅里葉變換重建圖像，根據奈奎斯特采樣定理，會導致重建圖像產生混疊效應。動態磁共振成像數據在時空域具有很強的稀疏特性，使得壓縮感知(compressive sensing，CS)技術被廣泛應用于MR圖像重建。CS理論指出，如果一個信號在變換域是稀疏的，且變換基和測量矩陣是不相關的(又稱有限等距性質，RIP)，則采用CS方法可以從降采樣(遠小于奈奎斯特采樣率條件下)的數據樣本當中，通過非線性重建算法完美重建該信號。

dMRI重建方法可以分為在線模式和離線模式。采用離線模式重建時，需要在重建之前獲得整個dMRI序列的采樣數據。常見的離線重建方法有運動校正[2-3]、字典學習[4-6]、低秩近似[7-8]等。這些方法充分利用整個dMRI序列稀疏特性進行高精度重建，其缺點是重建速度較慢，并且需要等待較長的掃描時間。采用在線方法重建時，每一幀的重建僅僅跟之前的幀有關，可以邊掃描邊重建，節省了等待時間，但由于缺乏整個序列的完整信息，重建精度無法保障，同時對重建速度也提出了更高的要求。

在線重建有兩種常用的方案：串行和并行。串行方案通常利用圖像或變換域中相鄰幀之間的稀疏性，這也是大多數現有在線方法中常用的策略。然而，這些方法往往會導致累積誤差。Chen等[9]采用了一種新的動態全變差(dTV)的并行重建方案來解決誤差累積問題。但該方法每次只能利用兩幀之間的稀疏作為重建的先驗知識，導致比離線方法精度低。

本文提出了一種基于并行自適應字典學習的實時動態磁共振重建方法，將原本在離線模式下運行的相對較慢的字典學習算法應用到在線模式當中來。如圖1所示，以高精度采樣的第一幀為參考，實現對任意n個相鄰幀MR圖像的實時在線重建。以三維圖像小塊作為重建對象，采用K-SVD和batch-OMP兩種算法優化重建速度。實驗結果表明，該算法在重建精度和速度上都具有一定優勢。

圖1 并行自適應字典學習dMRI重建原理(以3幀為例)

1 壓縮感知MRI理論

MRI采集k域(傅里葉變換域)信號，再通過逆傅里葉變換重建空域信號。傅里葉變換是一種線性變換，要求k域信號數必須等于圖像域的像素數。2007年Lustig等[10]提出了壓縮感知磁共振成像(CSMRI)的概念，以遠低于奈奎斯特頻率的采樣頻率對k空間信號進行采集，從而極大降低成像時間，并通過非線性重建算法完美重建原始圖像。

給定一個dMRI圖像序列表示為x，k-t空間的欠采數據為y，壓縮感知dMRI重建問題可以歸結為如下l0范數最小化問題。

(1)

CSMRI需要滿足3個條件：

(1) 稀疏矩陣Ψ，使得原始MR圖像在稀疏域中僅有少量非零系數；

(2) 矩陣FuΨ-1滿足RIP條件，一般采用隨機或偽隨機采樣矩陣，降低采樣數據的相關性；

(3) 非線性重建算法，能快速準確重建原始圖像。

2 自適應字典學習并行重建模型

2.1 算法總體架構

本文設計一套基于并行自適應字典學習的實時動態磁共振重建方法，算法流程如圖2所示。

圖2 流程框圖

2.2 字典學習和稀疏編碼

過完備字典在信號的稀疏表示當中有著廣泛的應用。各種基于圖像庫訓練的非自適應圖像重建方法因其適應性較差、重建效果一般等缺點逐漸被基于自適應字典學習的圖像重建算法所取代。這里的自適應是指由壓縮感知系統采集到的測量值作為圖像觀測樣本數據來訓練字典。由于字典學習的計算量會隨著字典尺寸的增加而快速增加，我們通常將待處理的圖像分為重疊的小塊，來代替整幅圖像進行稀疏表示。字典學習問題可描述為：

(2)

式中：D為過完備字典；Ri為重疊取塊算子；αi為x的稀疏表示系數；Γ={α1,α2，…,αI}為稀疏稀疏αi的集合；T0為稀疏度閾值常數。式(2)中第一項和稀疏約束條件確保了過完備字典對每個圖像塊的最優稀疏逼近。第二項為數據保真項，v是一個與k空間采樣時疊加的高斯白噪聲的標準差σ有關(v=λ/σ)的參數。

由于dMRI動態掃描時，每一幀的主要器官相對位置近似不變，我們選取序列的第一幀圖像為高精度采樣(一般取大于50%的采樣率)，為后續幀的重建提供高精度參考，這可由磁共振機的一次預掃描完成。第一幀與其余任意相鄰的n-1幀圖像組成一系列n幀的子序列圖像Xs(j)，以3幀為例，如圖1所示，Xs(j)=[x1,x2j,x2j+1]，j=1,2,…，x2j表示原dMRI序列x中的第2j幀圖像。只要相鄰兩幀采樣完畢，即可進行實時并行重建，無需等待后續幀采樣完成。每個子序列的重建與其他子序列不相關。

式(2)可通過以下過程進行優化求解：

字典學習過程，固定子序列圖像Xs，將問題轉化為求解式(3)中D和ɑ最優解的子問題。

(3)

(4)

(5)

算法1Batch-OMP算法

1: 輸入：α0=DTx,0=xTx,G=DTD,重建誤差

2: 輸出：稀疏系數γ滿足x≈Dγ

3: Init: SetI:=(),L1:=[1],δ0:=0,n:=1

4:whilen-1>do

7:ifn>1then

8:w:=Solve forw{Ln-1w=gIn}

10:endif

13:βn=GIncn

14:αn:=α0-βn

17:n:=n+1

18:endwhile

表1 算法1參數說明

算法1中，上標n為迭代次數，D為初始化DCT字典，x為字典學習圖像小塊對應的列向量，α為稀疏系數。這里G在每次迭代當中保持不變，可以在迭代之前預先計算好，以節省重建時間。在每次迭代過程當中只更新α當中的一列(|αk|取值最大的那一列)。算法1采用丘拉斯基分解方法(8-9)，將正定矩陣GIn,In分解為下三角矩陣L和其轉置的乘積的分解，減少了矩陣求逆過程的迭代運算量。最后采用誤差約束模型(15-16)循環更新D和α。

2.3 圖像重構

圖像重構過程，固定過完備字典和稀疏系數，將問題轉化為求解式(5)中待重建子序列Xs最優解的子問題。

(6)

式(6)是一個普通的最小二乘法問題，對唯一的變量xs求導并令其等于0得到：

(7)

直接對式(7)進行求解，計算量相當大，因為要對一個P×P的矩陣進行求逆才可得到重建結果，這使得求解此問題的時間復雜度達到O(p3)。我們可將運算由圖像空間變換到傅里葉空間進行，對式(7)兩邊進行傅里葉變換得到：

(8)

(9)

3 實 驗

所有的重建實驗都在一臺帶有3.0 GHz Intel i7 CPU的筆記本電腦上使用MATLAB 2013a進行仿真。選取一組心肌灌注dMRI序列(分辨率為192×192×30幀)用于評估算法的性能。偽射線采樣矩陣用來對測試圖像在k-t空間進行欠采樣，模擬磁共振機加速成像。如圖3所示，第一幀采樣率為50%，其他幀采樣率為15%。

(a) 第一幀采樣矩陣(b) 其他幀采樣矩陣圖3 采樣矩陣

以3幀構成一個子序列為例，進行并行重建。圖像質量的客觀評價方法采用峰值信噪比(Peak Signal to Noise Rate,PSNR)，單位為dB。PSNR計算公式如下：

PSNR(x,y)=20lg(Imax/RMSE(x,y))

(10)

式中：x和y分別為原始圖像和重建圖像；Imax為圖像像素最大值，一般為了減少矩陣運算量，要對圖像進行歸一化操作，則Imax=1。RMSE為兩幅圖像的均方根誤差。

實驗參數設置：為了快速并行重建dMRI序列，本文算法將三維圖像塊大小設置為最小的3×3×3，經過實驗測試，重建圖像質量與圖像塊尺寸關系不大，但是重建速度會隨著圖像塊尺寸的增加而迅速下降。字典大小設為27×108，每次在線學習的原子個數為5 000。迭代次數默認設置為20次，迭代結束條件設置為兩次迭代的PSNR值誤差在0.02 dB以內。對比的算法有離線模式的k-t SLR[8]和在線模式的TV[9]、DTV[9]三種。三個對比實驗算法的參數均按其最優參數設置。

圖4給出了本文算法重建一個dMRI子序列(以3幀為例)的過程。經過14次迭代，耗時2.8秒，就將k空間欠采樣后零填充的低信噪比圖像(31.1 dB)，重建為高信噪比的圖像(36.3 dB)。圖5給出了重建前后的圖像視覺對比，可見，重建后的圖像已經完全去除了欠采樣造成的偽影。

圖4 迭代次數與PSNR關系

(a) 零填充圖像 (b) 重建后圖像圖5 重建前后圖像對比

圖6對比了無采樣噪聲條件下各種算法的重建效果，以列為單位，第一行為各自的重建圖像，第二行為重建圖像與原始圖像的誤差圖，為凸顯效果，我們將誤差放大5倍進行比較。顯然本文算法的誤差最小，尤其是關鍵的心臟部位，視覺效果最為逼真。

圖6 各算法重建MR圖像和誤差圖像比較

在有采樣噪聲情況下，假設k空間疊加的高斯白噪聲標準差為σ，式(9)中λ=q/σ的取值就與重建效果息息相關。圖7給出了參數q與重建圖像的RMSE之間的關系。采樣信噪比分別取23 dB(強噪聲環境，σ=0.05)、37 dB(中等噪聲環境，σ=0.01)和57 dB(低噪聲環境，σ=0.001)。可見，q=0.005是考慮了在不同的采樣噪聲環境下的一個最優選擇。尤其是在強噪聲環境下，q=0.005時重建效果最佳。

圖7 重建均方根誤差與q的關系

本文算法適用于任意相鄰的n幀子序列構成的并行重建，表2給出了子序列幀數n與重建PSNR值和算法重建時間之間的關系(20%的采樣率條件下)。可見，以3幀作為一個子序列重建單元重建的精度(PSNR值)是最高的，而重建的時間又是最少的。并且各種參數設置下，重建精度都在可控的誤差范圍內，體現了本文算法的穩定性。

表2 重建子序列幀數與重建效果之間關系

本文算法與三個對比算法的重建PSNR值對比結果如表3所示。在不同的采樣率下，本文算法均是最佳的。尤其在15%的低采樣率條件下，相比其他算法有0.4 dB以上的性能提升。以3幀為一個dMRI子序列為例，本文算法平均每1.4 s即可重建一幀圖像，重建速度遠遠快于離線模式(k-t SLR)，但慢于兩種采用TV方法的算法。

表3 不同采樣率下重建圖像的PSNR值

4 結 語

快速成像一直是MRI創新的重點，從最早期掃描一個序列要花費一個多小時，到如今最新的壓縮感知技術已然可以實現10倍以上加速的超快速成像。本文使用一種基于自適應字典學習的動態磁共振并行重建算法，實現了對dMRI圖像的實時高精度重建。與傳統的離線模式k-t SLR算法和在線模式DTV算法相比，本文算法在不同的采樣率下均有0.2 dB以上的性能提升，為實時在線重建dMRI提供了一種解決方案，在磁共振成像領域具有廣闊的應用前景。