基于奇異值分解的改進GRAPPA算法研究

魏小琴, 何汶靜, 李 楊， 杜 勇

(川北醫學院附屬醫院 a.醫學影像學院； b.影像科，四川 南充 637000)

0 引 言

磁共振成像(Magnetic Resonance Imaging,MRI)是利用磁共振原理，對物體內部結構進行成像，MRI由于其超高的軟組織對比度，在醫療診斷中有著舉足輕重的作用。磁共振成像中使用多通道相控陣線圈接收磁共振掃描信號，每個線圈含有整個成像空間中的所有信號，這種多通道采集的信號冗余使得加速成像成為可能。并行磁共振成像(parallel Magnetic Resonance Imaging, pMRI)技術正是利用這種信號冗余可壓縮的特性以及頻率空間(也被稱為K空間(k-Space))的相關性，通過減少信號的采集量(即減少相位編碼次數)來實現節約掃描時間提升掃描速度。影響pMRI成像的主要兩個因素是：K空間信號數據采樣掃描時間與由K空間信號數據進行重建的圖像質量[1]。pMRI在醫療圖像高速成像及圖像重建領域已有長足發展，該技術已應用于各大磁共振設備。

pMRI圖像重建技術目前主要分為兩類：一類為基于圖像域的重建，該類算法主要利用線圈靈敏度進行圖像重構，得到無偽影圖像，代表算法有敏感度編碼(Sensitivity Encoding，SENSE)及隨后發展出的擴展方法SC-SENSE(Self-calibrating)、PILS(Parallel Imaging Local Sensitivities)等，該類算法要求通道線圈有更為精確的靈敏度，線圈靈敏信息由采集低頻信號計算而得，當靈敏度已知，文獻[2-3]中將壓縮感知與SENSE 結合，增加正則項，取得了較好的重構效果，而準確估計往往非常困難。另一類是基于K空間數據重建技術，該類技術采集中間低頻信號行數據作為自校準信號( Auto Calibration Signals,ACS)，用于計算出多通道線圈的權重系數，用這些權重系數擬合欠采樣缺失數據，再將其重建為診斷圖像。這種多通道直接重構最終圖像，不需要估計線圈靈敏度，代表算法有AUTO-SMASH(Simulataneous Acquisition of Spatial Harmonics)、GRAPPA(Generalized Auto-calibrating Partially Parallel Acquisitions)和VD-AUTO-SMASH(Variable Density-AUTO-SMASH)算法[4]等。

目前GRAPPA算法對2倍降采樣能夠很好重建出沒有卷褶偽影的圖像，但對于3倍以及3倍以上的高倍數降采樣的重建，由于擬合精度不夠導致重建結果SNR降低或是仍然還有卷褶偽影殘留，重建效果欠佳。這對于磁共振三維采集、動態成像、實時成像以及快速成像而言,提高GRAPPA算法是一個亟需解決的問題，本文提出一種基于奇異值分解(Singular Value Decomposition,SVD)算法來改變GRAPPA算法中的權重計算方式，以此降低矩陣維度、降低計算權重方程的條件數，增加權重的穩定性，提升計算速度，提高多倍降采樣的情況下的圖像重建質量。

GRAPPA算法有著更好的魯棒性，對成像物體的運動不敏感，但GRAPPA算法中對多通道線圈的權重計算要求較高，其精確性關系到重建結果的質量。近年來并沒有提出新的針對平面內的GRAPPA算法的改進算法，大量的帶有GRAPPA、SENSE關鍵字的算法如slice-GRAPPA、split-slice-GRAPPA等算法都是進行同時多層成像的技術方法[5]，是基于GRAPPAS算法的核心層面思想在層面方向進行加速的算法，與平面內的并行成像算法GRAPPA并非相同用途。也有在平面內非笛卡爾采集方式(如radial、spiral等采集方式)下應用GRAPPA的研究[6-7]。本文提出的優化方法在普通的笛卡爾采集方式下即可完成，它在運用GRAPPA算法計算時對數據量進行優化和降維，雖然增加了算法復雜度，但它對參與計算的數據量進行了縮減，減少了計算時間，且不降低計算的精度，還提高了GRAPPA算法的重建效果，是更底層的算法改進。有報道提出關于GRAPPA的優化算法，例如采用偽隨機采樣模板進行欠采樣[8]，在擬合的過程加入正則項進行權重的計算[9]等，其中最有效的步驟仍然是基本的GRAPPA算法。

1 計算方法與數據收集

1.1 GRAPPA算法求解

GRAPPA算法在全采樣數據ACS計算權重，在非ACS信號欠采樣數據獲取的條件下擬合欠采樣數據點的值，K空間數據欠采樣與計算權重擬合欠采樣數據點的過程如圖1所示，利用ACS信號對欠采樣數據進行權重求解[10]：

sn(kx,ky+rΔky)=

(1)

式中：(kx,ky+rΔky)為未采樣數據的坐標，r為采集點到前一已采集點的偏移量，r=1,2,…,R-1；R為欠采樣倍數，也是后面提到的加速倍數；Δkx、Δky分別為x、y方向的單位偏移量;l為線圈索引;L為線圈數;h為x方向索引；H1、H2為x方向插值窗界限；t為y方向索引；N1、N2為y方向插值窗界限；wj, r(l,h,t)為第l個線圈插值窗內的已采樣數據權重值。

圖1 K空間數據欠采樣與線性擬合欠采樣數據過程

設第ki行為n行線圈ACS采集數據，第n個欠采樣數據Dn(ki)可以寫成求和的形式:

(2)

式中:w(k)為權重數據;D(k)為采集數據;Dn(ki)為欠采樣數據插值目標點;kernel為插值窗內的已采樣數據點的集合。根據已采集數據D(k)及其權重w(k)來估計，上述計算權重w(k)的過程變成矩陣求解過程，Dn(ki)為列向量T，D(k)為矩陣S，w(k)為矩陣W，上式轉換為：

SW=T

(3)

對于不同的插值窗，預測單個未采集數據點所需要的源采樣數據點的個數

N=L(H2-H1)[(N2-N1)+1]

(4)

如果采集的ACS行的數據個數為Y，則S矩陣是Y×N矩陣，向量T的長度是Y。由GRAPPA算法采樣及權重計算過程可知，S矩陣為梯度線圈采集的K空間數據，根據權重矩陣大小和ACS的數量不同，大多數情況下S的每個單線圈權重計算不是一個方陣，所以式(3)是一個超定方程，GRAPPA算法的權重矩陣W使用普通最小二乘法求解得到：

W=(SHS)-1SHT

(5)

式中:SH是S的共軛轉置，(SHS)-1是對(SHS)求逆矩陣。普通最小二乘法有兩個基本假設是：①S非奇異或者列滿秩;② 向量T或者矩陣S存在加性噪聲[11]。但在工程應用中，采樣數據間相關性很強，矩陣S無法做到列滿秩，因此GRAPPA算法很難得到最優解，即使調整插值窗也難以解決秩虧缺。SHS矩陣的條件數:

cond(SHS)=[cond(S)]2

(6)

S的誤差δS和T的誤差δT對W的誤差δW的影響與S的條件數的平方成正比，影響到整個權重值求解[12]，權重值的計算精度降低會直接影響插值的準確性，降低重建后的圖像質量。

1.2 基于SVD的GRAPPA方法

為解決在其采樣數據后，多倍降采樣數據在采樣及計算過程中出現的噪聲及計算權重的超定方程病態問題，引入奇異值分解方法(Singular Value Decomposintion Technique，SVD)[13-15]，在降低權重維度，提高采樣倍數情況下，同時降低采樣時間，提高采樣效率，提升重構圖像質量。

基于SVD的GRAPPA算法步驟：

(1) 將圖像進行多倍降采樣，如圖1所示，將已采集到K空間數據與ACS進行線性加權擬合，用SVD計算得到各線圈權重系數；

(2) 利用權重系數和已采樣K空間數據計算欠采樣數據；

(3) 對完整K空間數據進行傅里葉逆變換，重構各個線圈圖像；

(4) 將各個線圈圖像整合為一張輸出重構圖。

采樣矩陣S為奇異陣，由SVD方法，將矩陣S進行奇異值分解[16]：

S=UΣVT

(7)

(8)

故最小的奇異值σmin會增大條件數，出現比較大的擾動。因此需要計算超定方程的有效秩r，可以采用歸一化奇異值的方法

σmin/σmax>ε

(9)

選擇一個閾值ε，ε是某個很小的正數，例如可以選擇ε=0.1，對較小奇異值進行截斷。舍去小于ε的奇異值及其對應特征向量，相當于舍去了相關性較強的向量，降低了條件數。

利用壓縮后的矩陣再對降采樣缺失行進行填充，其間仍采用最小二乘法式(5)求解，獲得完整的K空間數據。

1.3 數據收集

實驗驗證過程中的驗證數據來自ALLTECH 1.5T磁共振掃描系統獲取的8通道頭線圈的全部原始數據(raw data)，采用梯度回波脈沖序列(Gradient Echo Pulse Sequence，GRE)，脈沖序列的參數信息為：回波時間，TE=9 ms，脈沖重復時間，TR=143 ms，水脂偏移系數，Water Fat Shift=0.8，矩陣大小為256×278。數據采集方式如圖1所示,并行磁共振成像沿著相位編碼方向上，以加速因子R進行周期性軌跡采樣。將相應相位編碼行進行置零以模擬降采樣過程，保留K空間中間低頻部分行作為計算卷積核權重的ACS部分。完全采樣的K空間數據重建出的圖像作為對比參考圖像。

本文將對加速因子R為2～5這4種情況進行數據降采樣，K空間低頻部分ACS部分取16和32兩種情況，為提高采樣時間，減少采樣引起的運動偽影，大部分ACS取16。根據采樣數據，將對文中所描述的改進算法進行圖像重建，與傳統的GRAPPA算法圖像重建進行對比分析,圖像重建仿真計算在Matlab軟件上進行。

1.4 算法評價

為評價改進算法在并行成像多倍降采樣下重構圖像質量，從圖像重構質量與算法運行時間評價算法優劣。

圖像質量主要從改進算法重構圖像與傳統算法重構圖像進行分析，將重建后的圖像與全采樣圖像進行比較，分別生成改進算法與傳統算法差異十倍顯示比較圖，改進算法與全采樣重建圖像誤差圖，傳統算法與全采樣重建圖像誤差圖。

時間評價算法采用卷積核權重算法時間及算法重建時間作為評價標準。

2 結果分析

2.1 圖像結果分析

本研究用8通道全采樣頭部圖像作為算法重建過程對比標準圖如圖2所示，實驗分別用加速因子R為2～5，ACS線為16、32進行圖像傳統算法與改進算法SVD+GRAPPA重建。圖3所示為在不同采樣條件下的改進算法SVD+GRAPPA、傳統GRAPPA算法重建圖與誤差圖，該誤差圖為兩種重建算法差值放大10倍后的效果圖。

圖2 8通道頭線圈全采樣圖像

圖3 在不同采樣條件下改進算法SVD+GRAPPA、傳統GRAPPA算法重建圖與誤差圖比較

從圖像上分析，當加速因子R=2，ACS=16時，改進算法SVD+GRAPPA和傳統GRAPPA算法重建的圖像幾乎都能達到臨床診斷效果，兩算法間重建差異小。當加速因子R=3時，傳統GRAPPA算法已經出現異常增亮影，即卷褶偽影。當R為5時，傳統算法已經完全無法分辨大腦內部結構，白質和灰質為一整片區域，無密度比較。加速因子越大，傳統GRAPPA算法中偽影越多，可識別組織信息越少，欠采樣會導致數據采集完整性降低，傳統算法欠采樣程度越大，圖像質量損失越嚴重。

改進算法SVD+GRAPPA重構圖中，當加速因子為R=2，ACS=16時，改進算法重構圖與傳統算法差值圖幾乎為全黑圖像，即兩圖無明顯差異。當加速因子R=3、4時，改進算法重構圖與傳統算法差值圖為非全黑圖像，即兩圖具有明顯差異，SVD+GRAPPA重構圖有較為細致的解剖結構，未出現異常增亮影，有清晰的白質和灰質密度影。當加速因子R=5時，圖像有較為模糊的解剖結構，出現噪聲影，但相較于傳統算法噪聲量少。當加速因子R逐漸增大時，差值圖像包含的信息量也越來越多，相較于傳統GRAPPA算法，改進算法圖像質量損失較少。

當R=4,ACS=16時，改進算法SVD+GRAPPA重構圖像質量達到R最大化時，降采樣重建效果最優。當R更大時，需增加ASC采樣線，才能保證優質重構圖像質量，可應用于臨床診斷。

2.2 重建程序運算時間分析

重建時間為定量分析重要評價標準，主要從兩方面分析，一是從卷積核權重計算時間上比較，二是從整個算法采集、權重計算、圖像重建過程時間進行比較。

表1為改進算法SVD+GRAPPA與原算法GRAPPA卷積核權重計算時間比較，在 GRAPPA算法基礎上解病態方程，相較于原算法來說，改進算法理論復雜度更高，但欠采樣對矩陣進行壓縮使得計算時間下降，運行效率更高。由表1可見，加速因子R=2、ACS=16時，改進算法比原算法用時減少79 ms，在原基礎上效率提高39.6%，當加速因子為R=3時，改進算法用時減少67 ms，效率提高38.3%，當加速因子為R=4時，改進算法用時減少60 ms，效率提高38.6%。

表1 改進算法SVD+GRAPPA與原算法GRAPPA卷積核權重計算時間比較

表2為改進算法SVD+GRAPPA與原算法GRAPPA在圖像采集、權重計算及重建過程等MRI并行計算整個過程時間比較。由表2可見，加速因子R=2、ACS=16時，改進算法比原算法用時減少0.61 s，在原基礎上效率提高了51.3%，當加速因子R=3時，改進算法用時減少0.48 s，效率提高49.6%，當加速因子R=4時，改進算法比原算法用時減少0.49 s，效率提高53.4%。

表2 改進算法SVD+GRAPPA與原算法GRAPPA重構時間比較

3 結 語

GRAPPA方法基于K空間，利用多通道相位控陣線圈冗余數據還原未采樣數據行，其插值權重計算的準確性對重建算法的結果至關重要。原GRAPPA算法采用最小二乘法，但未考慮采樣矩陣中向量的相關性問題，導致算法對插值窗大小和自動校準線數量的選取非常敏感，在加速因子R為3、4時，重建效果不佳。本文提出基于奇異值分解的自校準并行采集成像方法，相比傳統算法更多的考慮了各已知采樣數據向量之間的相關性。算法計算采樣矩陣的奇異值后，去掉部分較小奇異值以達到降低條件數的效果。

實驗發現改進算法相比原GRAPPA方法在加速因子R=3、4時表現出更細致的圖像信息，能更好地分辨出大腦組織解剖結構。同時，MRI對人體部位掃描一次需要花費很長一段時間，肺部或心臟等部位的運動會降低圖像質量，改進算法增加采樣倍數同時優化重建圖像算法，減少受檢者檢查時間。實驗表明，基于奇異值分解的改進算法重建效果優于原GRAPPA方法。