基于奇異值分解的改進(jìn)GRAPPA算法研究

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

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

0 引 言

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

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

目前GRAPPA算法對(duì)2倍降采樣能夠很好重建出沒有卷褶偽影的圖像，但對(duì)于3倍以及3倍以上的高倍數(shù)降采樣的重建，由于擬合精度不夠?qū)е轮亟ńY(jié)果SNR降低或是仍然還有卷褶偽影殘留，重建效果欠佳。這對(duì)于磁共振三維采集、動(dòng)態(tài)成像、實(shí)時(shí)成像以及快速成像而言,提高GRAPPA算法是一個(gè)亟需解決的問題，本文提出一種基于奇異值分解(Singular Value Decomposition,SVD)算法來改變GRAPPA算法中的權(quán)重計(jì)算方式，以此降低矩陣維度、降低計(jì)算權(quán)重方程的條件數(shù)，增加權(quán)重的穩(wěn)定性，提升計(jì)算速度，提高多倍降采樣的情況下的圖像重建質(zhì)量。

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

1 計(jì)算方法與數(shù)據(jù)收集

1.1 GRAPPA算法求解

GRAPPA算法在全采樣數(shù)據(jù)ACS計(jì)算權(quán)重，在非ACS信號(hào)欠采樣數(shù)據(jù)獲取的條件下擬合欠采樣數(shù)據(jù)點(diǎn)的值，K空間數(shù)據(jù)欠采樣與計(jì)算權(quán)重?cái)M合欠采樣數(shù)據(jù)點(diǎn)的過程如圖1所示，利用ACS信號(hào)對(duì)欠采樣數(shù)據(jù)進(jìn)行權(quán)重求解[10]：

sn(kx,ky+rΔky)=

(1)

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

圖1 K空間數(shù)據(jù)欠采樣與線性擬合欠采樣數(shù)據(jù)過程

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

(2)

式中:w(k)為權(quán)重?cái)?shù)據(jù);D(k)為采集數(shù)據(jù);Dn(ki)為欠采樣數(shù)據(jù)插值目標(biāo)點(diǎn);kernel為插值窗內(nèi)的已采樣數(shù)據(jù)點(diǎn)的集合。根據(jù)已采集數(shù)據(jù)D(k)及其權(quán)重w(k)來估計(jì)，上述計(jì)算權(quán)重w(k)的過程變成矩陣求解過程，Dn(ki)為列向量T，D(k)為矩陣S，w(k)為矩陣W，上式轉(zhuǎn)換為：

SW=T

(3)

對(duì)于不同的插值窗，預(yù)測(cè)單個(gè)未采集數(shù)據(jù)點(diǎn)所需要的源采樣數(shù)據(jù)點(diǎn)的個(gè)數(shù)

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

(4)

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

W=(SHS)-1SHT

(5)

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

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

(6)

S的誤差δS和T的誤差δT對(duì)W的誤差δW的影響與S的條件數(shù)的平方成正比，影響到整個(gè)權(quán)重值求解[12]，權(quán)重值的計(jì)算精度降低會(huì)直接影響插值的準(zhǔn)確性，降低重建后的圖像質(zhì)量。

1.2 基于SVD的GRAPPA方法

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

基于SVD的GRAPPA算法步驟：

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

(2) 利用權(quán)重系數(shù)和已采樣K空間數(shù)據(jù)計(jì)算欠采樣數(shù)據(jù)；

(3) 對(duì)完整K空間數(shù)據(jù)進(jìn)行傅里葉逆變換，重構(gòu)各個(gè)線圈圖像；

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

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

S=UΣVT

(7)

(8)

故最小的奇異值σmin會(huì)增大條件數(shù)，出現(xiàn)比較大的擾動(dòng)。因此需要計(jì)算超定方程的有效秩r，可以采用歸一化奇異值的方法

σmin/σmax>ε

(9)

選擇一個(gè)閾值ε，ε是某個(gè)很小的正數(shù)，例如可以選擇ε=0.1，對(duì)較小奇異值進(jìn)行截?cái)?。舍去小于ε的奇異值及其?duì)應(yīng)特征向量，相當(dāng)于舍去了相關(guān)性較強(qiáng)的向量，降低了條件數(shù)。

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

1.3 數(shù)據(jù)收集

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

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

1.4 算法評(píng)價(jià)

為評(píng)價(jià)改進(jìn)算法在并行成像多倍降采樣下重構(gòu)圖像質(zhì)量，從圖像重構(gòu)質(zhì)量與算法運(yùn)行時(shí)間評(píng)價(jià)算法優(yōu)劣。

圖像質(zhì)量主要從改進(jìn)算法重構(gòu)圖像與傳統(tǒng)算法重構(gòu)圖像進(jìn)行分析，將重建后的圖像與全采樣圖像進(jìn)行比較，分別生成改進(jìn)算法與傳統(tǒng)算法差異十倍顯示比較圖，改進(jìn)算法與全采樣重建圖像誤差圖，傳統(tǒng)算法與全采樣重建圖像誤差圖。

時(shí)間評(píng)價(jià)算法采用卷積核權(quán)重算法時(shí)間及算法重建時(shí)間作為評(píng)價(jià)標(biāo)準(zhǔn)。

2 結(jié)果分析

2.1 圖像結(jié)果分析

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

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

圖3 在不同采樣條件下改進(jìn)算法SVD+GRAPPA、傳統(tǒng)GRAPPA算法重建圖與誤差圖比較

從圖像上分析，當(dāng)加速因子R=2，ACS=16時(shí)，改進(jìn)算法SVD+GRAPPA和傳統(tǒng)GRAPPA算法重建的圖像幾乎都能達(dá)到臨床診斷效果，兩算法間重建差異小。當(dāng)加速因子R=3時(shí)，傳統(tǒng)GRAPPA算法已經(jīng)出現(xiàn)異常增亮影，即卷褶偽影。當(dāng)R為5時(shí)，傳統(tǒng)算法已經(jīng)完全無法分辨大腦內(nèi)部結(jié)構(gòu)，白質(zhì)和灰質(zhì)為一整片區(qū)域，無密度比較。加速因子越大，傳統(tǒng)GRAPPA算法中偽影越多，可識(shí)別組織信息越少，欠采樣會(huì)導(dǎo)致數(shù)據(jù)采集完整性降低，傳統(tǒng)算法欠采樣程度越大，圖像質(zhì)量損失越嚴(yán)重。

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

當(dāng)R=4,ACS=16時(shí)，改進(jìn)算法SVD+GRAPPA重構(gòu)圖像質(zhì)量達(dá)到R最大化時(shí)，降采樣重建效果最優(yōu)。當(dāng)R更大時(shí)，需增加ASC采樣線，才能保證優(yōu)質(zhì)重構(gòu)圖像質(zhì)量，可應(yīng)用于臨床診斷。

2.2 重建程序運(yùn)算時(shí)間分析

重建時(shí)間為定量分析重要評(píng)價(jià)標(biāo)準(zhǔn)，主要從兩方面分析，一是從卷積核權(quán)重計(jì)算時(shí)間上比較，二是從整個(gè)算法采集、權(quán)重計(jì)算、圖像重建過程時(shí)間進(jìn)行比較。

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

表1 改進(jìn)算法SVD+GRAPPA與原算法GRAPPA卷積核權(quán)重計(jì)算時(shí)間比較

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

表2 改進(jìn)算法SVD+GRAPPA與原算法GRAPPA重構(gòu)時(shí)間比較

3 結(jié) 語(yǔ)

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

實(shí)驗(yàn)發(fā)現(xiàn)改進(jìn)算法相比原GRAPPA方法在加速因子R=3、4時(shí)表現(xiàn)出更細(xì)致的圖像信息，能更好地分辨出大腦組織解剖結(jié)構(gòu)。同時(shí)，MRI對(duì)人體部位掃描一次需要花費(fèi)很長(zhǎng)一段時(shí)間，肺部或心臟等部位的運(yùn)動(dòng)會(huì)降低圖像質(zhì)量，改進(jìn)算法增加采樣倍數(shù)同時(shí)優(yōu)化重建圖像算法，減少受檢者檢查時(shí)間。實(shí)驗(yàn)表明，基于奇異值分解的改進(jìn)算法重建效果優(yōu)于原GRAPPA方法。