基于擴展卡爾曼濾波的磁共振指紋參數量化優化算法*

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

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

1 引 言

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

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

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

2 參數量化方法

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

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

2.1 模型建立

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

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

(1)

將MRF數據演變過程看作一個非線性系統，系統動態方程可以表示為：

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

(2)

Zk=h[k,Xk]+vk

(3)

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

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

(4)

若采用bSSFP序列[1]進行MRF數據采集，根據Bloch方程，k時刻的系統變量與k-1時刻的系統變量之間是一個復雜的非線性關系。

已知k-1時刻的磁化矢量Mk-1，結合組織參數pk，要計算出下一時刻的磁化矢量Mk，需要經過三步：

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

(5)

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

Mt+dt=g(Mt,dt)

(6)

可另寫為：

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

(7)

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

(8)

其中，

(9)

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

Mk=g(M2,k,dt2)

(10)

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

Xk=f(Xk-1)=

(11)

2.2 參數量化

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

(12)

(2)預測協方差矩陣Pk|k-1

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

(13)

(3)求解卡爾曼增益Kk

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

(14)

(4)狀態更新X(k|k)

(15)

(5)協方差更新Pk|k

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

(16)

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

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

(17)

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

2.3 EKF-MRF優化

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

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

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

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

3 實現與結果

3.1 單組織參數量化結果

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

圖1 翻轉角與重復時間Fig.1 Flip angle and repetition time

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

圖2 噪聲對結果的影響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。

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

圖3兩種組織的磁化矢量

Fig.3Magnetic vector of two tissues

圖4參數估計曲線

Fig.4Curves of parameter estimation

3.2 二維模型結果

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

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

另外，為了與字典量化法比較，我們建立了字典。參數選擇：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為間隔)；得到字典入口數為91 266個。由組織特性可知，無T1

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

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

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

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

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

圖5頭部模型的量化結果

Fig.5Quatification results for brain model

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

4 結論

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