黃敏，張璐堅，周到，陳軍波

(中南民族大學(xué)生物醫(yī)學(xué)工程學(xué)院，武漢 430074)

1 引 言

磁共振成像(MRI)技術(shù)已被廣泛應(yīng)用于臨床檢查，但MRI一次掃描只能得到T1、T2或質(zhì)子密度中的一種加權(quán)像，對醫(yī)生的經(jīng)驗要求較高。2013年,Griswold等人提出了“磁共振指紋”(magnetic resonance fingerprinting,MRF)成像方法[1]，可通過一次掃描，得到多參數(shù)的定量圖像。

經(jīng)典MRF成像框架需要建立一個包含被掃描部位所有組織的參數(shù)信息的字典,字典大小與參數(shù)匹配的準確度和效率密切相關(guān)。該框架存在匹配時間過長，字典不具有通用性等缺陷。基于分組匹配方法[2]和奇異值分解(SVD)[3-4]后再匹配的方法可以加快匹配速度。基于壓縮感知[5-6]的方法可加速匹配過程，提高準確度，但對不同的采樣序列，當改變TR和FA后，需重新生成字典，增加了額外的時間和空間成本。

目前提出的基于無字典的MRF參數(shù)量化框架的方法中，有采用深度學(xué)習得到參數(shù)圖像，代替基于字典的匹配過程[7-8]，但深度學(xué)習需要大量數(shù)據(jù)作為訓(xùn)練支撐，可獲得的MRF采集數(shù)據(jù)和標簽數(shù)據(jù)非常少。Zhang等人提出了使用卡爾曼濾波的方式對參數(shù)進行量化，無需訓(xùn)練數(shù)據(jù)，但量化時間較長[9-10]。為此，我們研究了基于擴展卡爾曼濾波的磁共振指紋參數(shù)量化方法，并對量化速度和精度進行了優(yōu)化。

2 參數(shù)量化方法

MRF技術(shù)采用偽隨機變化的TR和FA組合(如1 000組)進行掃描，獲得對量化參數(shù)敏感的獨特信號。為縮短掃描時間，在每組TR時間內(nèi)，可采用單支螺旋軌跡采樣，得到K空間欠采樣數(shù)據(jù)。對采樣數(shù)據(jù)進行逆NUFFT得到空間域圖像，將空間域圖像各點對應(yīng)位置的多幀信號連接起來，即得到每個體素的“磁共振指紋”信號。在經(jīng)典MRF參數(shù)量化中，需要根據(jù)所有組織參數(shù)組合與序列參數(shù)，通過Bloch方程建立字典庫。將指紋與字典信號進行逐條匹配，通過匹配結(jié)果得到參數(shù)像。

卡爾曼濾波(KF)可用于線性系統(tǒng)的信號跟蹤和觀測，是一種最優(yōu)估計方法[11]。而利用擴展卡爾曼濾波(EKF)，可實現(xiàn)對非線性系統(tǒng)(MRF屬于此類)的優(yōu)化估計問題。我們先建立磁共振指紋成像模型；對模型的MRF信號進行EKF跟蹤觀測，使估計值趨近于真實值，實現(xiàn)對該信號的輸入?yún)?shù)進行反推和估計。該方法簡稱為EKF-MRF方法。

2.1 模型建立

在MRF中，磁場不是絕對均勻的，存在偏振頻率df。若k時刻某組織的磁化矢量信號為Mk=[Mxk;Myk;Mzk], 待量化的參數(shù)組合為pk=[T1;T2;df]，定義此時的聯(lián)合狀態(tài)向量為：

Xk=[Mk;pk]=[Mxk;Myk;Mzk;T1;T2;df]

(1)

將MRF數(shù)據(jù)演變過程看作一個非線性系統(tǒng)，系統(tǒng)動態(tài)方程可以表示為：

Xk=f[k-1,Xk-1]+wk

(2)

Zk=h[k,Xk]+vk

(3)

其中Zk為觀測值，wk∈R6×1為過程噪聲(近似高斯白噪聲)，其協(xié)方差為Qk。vk∈R2×1為觀測噪聲，其協(xié)方差為Rk。

由MRF原理可知，f(*)為非線性變換，h(*)為線性變換。k時刻的測量值[Mxk;Myk]與系統(tǒng)變量之間是線性關(guān)系，觀測矩陣H表示為：

(4)

若采用bSSFP序列[1]進行MRF數(shù)據(jù)采集，根據(jù)Bloch方程，k時刻的系統(tǒng)變量與k-1時刻的系統(tǒng)變量之間是一個復(fù)雜的非線性關(guān)系。

已知k-1時刻的磁化矢量Mk-1，結(jié)合組織參數(shù)pk，要計算出下一時刻的磁化矢量Mk，需要經(jīng)過三步：

(1)首先，經(jīng)過dt1時間的自由進動后，磁化矢量可通過Bloch方程計算，即：

(5)

其中，ω為偏振角頻率，dt為當前信號采集時刻到下次射頻施加之前的時間段dt1(bSSFP序列為0.5TR(k-1))。得到自由進動狀態(tài)下磁化矢量：

Mt+dt=g(Mt,dt)

(6)

可另寫為：

M1,k=g(Mk-1,dt1)

(7)

(2)施加下一RF脈沖后(翻轉(zhuǎn)角為αk，相位為0和180°相互交替的φk)，得到新的磁化矢量：

(8)

其中，

(9)

(3)再經(jīng)過dt2時間的自由進動，按照式(2)計算新的磁化矢量Mk。其中，dt為下次射頻施加時刻到下次信號采集時刻的時間段dt2(bSSFP序列為0.5TR(k))，可表示為：

Mk=g(M2,k,dt2)

(10)

由于每個體素的組織參數(shù)T1，T2和磁場不均度df是不變的。將磁化矢量與參數(shù)相結(jié)合，可得到k時刻的聯(lián)合狀態(tài)變量Xk。最終公式為：

Xk=f(Xk-1)=

(11)

2.2 參數(shù)量化

基于擴展卡爾曼濾波的MRF參數(shù)量化過程，主要是對五個關(guān)鍵方程進行迭代計算：

(12)

(2)預(yù)測協(xié)方差矩陣Pk|k-1

Pk|k-1=FkPk-1|k-1FkT+Q

(13)

(3)求解卡爾曼增益Kk

Kk=Pk|k-1HT(HPk|k-1HT+R)

(14)

(4)狀態(tài)更新X(k|k)

(15)

(5)協(xié)方差更新Pk|k

Pk|k=(I-KkH)Pk|k-1

(16)

其中，f(Xk-1)為EKF-MRF的非線性狀態(tài)方程。將k-1時刻狀態(tài)變量代入式(11)，得到k時刻狀態(tài)的估計值。

通過狀態(tài)方程的一階線性化得到轉(zhuǎn)移矩陣，由于狀態(tài)向量Xk含有[Mk;pk]兩部分，所以最終狀態(tài)轉(zhuǎn)移矩陣為：

(17)

隨機輸入待量化的三個參數(shù)T1，T2和df的初始值，通過擴展卡爾曼濾波不斷跟蹤預(yù)測信號，對方差進行估計，更新參數(shù)組合的值,最終使參數(shù)組合通過狀態(tài)方程后，得到的信號趨近于信號真實值。說明此時參數(shù)值已趨近于組織真實值，即得到了量化結(jié)果。

2.3 EKF-MRF優(yōu)化

為減少MRF掃描時間，使用500組TR和FA的脈沖序列進行數(shù)據(jù)采集，但經(jīng)過500次卡爾曼迭代后，未能得到收斂的參數(shù)值，需要對觀測值進行數(shù)據(jù)復(fù)制。

通過實驗發(fā)現(xiàn)，過程噪聲Q和觀測噪聲R對實驗結(jié)果的準確性影響較大，Q值直接迭代次數(shù)。在第一輪(500次)EKF迭代時，由于初始參數(shù)值是隨機的，計算的預(yù)測值偏離真實值較遠，我們選取較大的Q值，使卡爾曼運算中更偏向測量值。在第一輪卡爾曼迭代過程中，各參數(shù)值向真實值迅速趨近，但不能精確收斂。當參數(shù)值已經(jīng)趨近于真實值后(誤差變小)，再將Q值逐步調(diào)小，使其收斂在一個固定值附近。

實驗表明：當R值過小，收斂后量化的T1,T2參數(shù)值會較真實值偏小；當R值過大，各參數(shù)將不能準確收斂。經(jīng)過大量試驗后，我們將卡爾曼濾波算法中的數(shù)據(jù)重復(fù)次數(shù)控制在四輪，即2 000個點處，可以得到逼近真實值的量化結(jié)果。

相較于Zhang等人的方法[10]，本研究還在狀態(tài)方程的一階線性化過程中做了優(yōu)化。前者在一階線性化過程中將狀態(tài)方程分成兩個非線性方程和一個線性方程，做了兩次一階線性化過程。我們基于EKF原理，將狀態(tài)方程作為一個整體僅做一次一階線性化，計算速度大大加快。

3 實現(xiàn)與結(jié)果

3.1 單組織參數(shù)量化結(jié)果

翻轉(zhuǎn)角FA采用添加了柏林噪聲的偽隨機變化值，并加入隨機噪聲；TR取范圍為6.5～14 ms的隨機值，見圖1。模擬腦部組織參數(shù)值，在df為4 Hz，SNR為15 dB條件下，對EKF-MRF方法進行實現(xiàn)和分析。

圖1 翻轉(zhuǎn)角與重復(fù)時間Fig.1 Flip angle and repetition time

選取T1=500 ms，T2=70 ms，df=6 Hz，分析觀測噪聲和測量噪聲對實驗結(jié)果的影響，見圖2。圖2(a)中，當R值偏小時，df雖然能收斂，但收斂值偏??；(b)中，當R值偏大時，df不能收斂；(c)中，當Q值固定時，df在較大范圍內(nèi)震蕩，不能準確的收斂；(d)為本研究的df結(jié)果，能較快并準確的收斂。

圖2 噪聲對結(jié)果的影響Fig.2 Influence of noise on results

選取T1值最高的腦脊液(T1∶2 569 ms、T2∶329 ms)，以及T1值與T2值差距最大的肌肉(T1∶350 ms、T2∶70 ms)進行實驗，得到仿真的磁化矢量Mx和My，見圖3。

腦脊液和肌肉的參數(shù)估計值收斂曲線見圖4。測量數(shù)據(jù)在重復(fù)1輪至1 000個點時，估計值已收斂到較小范圍；到1 500個點時，能得到穩(wěn)定收斂且誤差很小的參數(shù)值，證明優(yōu)化方法可得到逼近真實值的參數(shù)估計值。

圖3兩種組織的磁化矢量

Fig.3Magnetic vector of two tissues

圖4參數(shù)估計曲線

Fig.4Curves of parameter estimation

3.2 二維模型結(jié)果

在配置為Intel Core(TM) i5-8300H 2.30 GHz 的PC上，在Matlab R2017b環(huán)境下，對上述算法進行實現(xiàn)。

為了將EKF-MRF優(yōu)化法與Zhang等人的方法[10]進行比較，使用相同的BrainWeb數(shù)據(jù)庫中人腦模型進行仿真。BrainWeb的模型大小為217×181×181，將每層擴充為256×256，在1T條件下，設(shè)置每層組織的T1、T2以及質(zhì)子密度值[4]。模擬信號采集的SNR為15 dB，df在范圍為0～12 Hz內(nèi)呈指數(shù)趨勢中心對稱變化，并對非感興趣區(qū)域添加掩膜。

另外，為了與字典量化法比較，我們建立了字典。參數(shù)選擇：T1值范圍為100～5 000 ms(100～2 000 ms之間，以20 ms為間隔；2 000 ms以上，300 ms為間隔)；T2范圍為20～3 000 ms(20～100 ms之間，以5 ms為間隔；100～200 ms之間，以10 ms為間隔；200 ms以上，200 ms為間隔)；df范圍為-20～20 Hz(2 Hz為間隔)；得到字典入口數(shù)為91 266個。由組織特性可知，無T1

EKF-MRF及字典法的參數(shù)量化圖像見圖5，匹配誤差及量化時間見表1。準確度采用參數(shù)量化結(jié)果與金標準之間的均方根誤差作為衡量標準。

結(jié)果表明，EKF-MRF法在量化的準確度和效率上均優(yōu)于字典直接匹配法。直接匹配法的最大缺陷是T2量化圖失真嚴重，df量化出現(xiàn)環(huán)狀偽影。EKF-MRF法T2量化精度明顯增強，誤差明顯縮小，df量化無環(huán)狀偽影。

表1 匹配誤差及量化時間Table1 Matching error and quantification time

Zhang等人的方法[10]，對BrainWeb數(shù)據(jù)進行64×64點量化，用時42 min；可估算該方法量化256×256點的圖像需要640 min。本研究的優(yōu)化方法，參數(shù)量化256×256點僅需約38 min，時間縮短了1個數(shù)量級，效率明顯提高。提高的關(guān)鍵在于：對測試數(shù)據(jù)的重復(fù)次數(shù)僅需要4次即可達到迭代的精度要求；雅克比矩陣僅用直接一次性計算，大大節(jié)省了時間。

為進一步測試量化的準確度，三個參數(shù)組合中保持兩個不變，觀測對剩余一個參數(shù)的量化影響。在T2=60 ms和df=6 Hz時，分析T1量化的準確性；在T1=600 ms和df=6 Hz時，分析T2量化的準確性；在T1=600 ms和T2=60 ms時，分析df量化的準確性。結(jié)果見圖6，表明EKF-MRF在T1,T2,df具有很好的量化效果。字典直接匹配法T1參數(shù)估計較準確，df估計誤差較大，T2誤差尤其大。

圖5頭部模型的量化結(jié)果

Fig.5Quatification results for brain model

圖6 參數(shù)值變化對準確度的影響Fig.6 Influence of parameter value change on accuracy

4 結(jié)論

基于擴展卡爾曼濾波的磁共振指紋參數(shù)量化優(yōu)化法，相較于字典量化法和Zhang等人的卡爾曼濾波法，在速度上有了很大提升，在T2參數(shù)量化的質(zhì)量上也有顯著改善。這種無字典的磁共振指紋參數(shù)量化方法，解決了不同部位、不同序列需要建立不同字典的問題。也使采樣序列的設(shè)計更加靈活，并且節(jié)省了儲存字典所需的空間。