張會冉，方金生，楊敬民

（1.閩南師范大學(xué)計算機學(xué)院，福建漳州363000;2.福建省粒計算及其應(yīng)用重點實驗室，福建漳州363000）

定量磁化率重建一般需經(jīng)過相位圖預(yù)處理和背景場去除，再進行偶極場反演.這種多步驟處理的方式，在每一步的運算中都可能引入新的計算誤差，該誤差會隨之傳入下一個步驟，導(dǎo)致誤差不斷地累加.因此，減少重建步驟數(shù)，是防止誤差累積的有效途徑.2014年，Li 等[1]提出聯(lián)合相位解纏繞和去背景場算法，單步即可完成相位解纏繞和去背景場.隨后，Bilgic 等[2]將去背景場方程式融合到磁化率重建表達式中，實現(xiàn)了單步去除背景場和磁化率重建.Liu 等[3]提出基于全局場的磁化率反演算法，引入預(yù)條件矩陣以提高運算收斂速度，實現(xiàn)由全局場圖直接重建磁化率圖.但是上述兩種算法均為基于線性場圖的磁化率重建，這對抑制噪聲和強磁化率偽影的效果有限，而采用非線性場圖重建磁化率更夠更為有效地抑制噪聲和磁化率偽影[4-5].因此，本文提出一種單步實現(xiàn)定量磁化率重建算法，將局部場與全局場的關(guān)系表達式引入磁化率重建非線性表達式中，并利用結(jié)構(gòu)補償?shù)姆椒ㄖ苯又亟ù呕蕡D像，有效提高磁化率重建精度.

1 算法原理

1.1 非線性全局場重建磁化率

采用局部場線性重建的形態(tài)學(xué)相似性偶極場反演算法（MEDI）在處理磁化率較大區(qū)域時，如腦部大血管，可能產(chǎn)生較強的磁化率偽影[6]，由此Wang 等引入非線性形式的MEDI 算法（ΝMEDI），則磁化率反演模型為[5]：

式（1）中λ為懲罰系數(shù)，χ為磁化率，W為二值化人腦掩模版，F(xiàn)H，F(xiàn)分別表示傅里葉逆變換和傅里葉變換，Bint表示為局部場，M為權(quán)重矩陣，G為三維梯度算子.由于式（1）收斂速度較慢，Milovic 等[4]提出基于交替方向法（ADMM）加速非線性重建算法（FAΝSI）.根據(jù)式（1），令z=FHDFχ，引用增廣拉格朗日函數(shù)表示為：

其中，s為增廣拉格朗日乘法子，μ為懲罰系數(shù)，根據(jù)背景場去除表達式(δ-ρ)?B?=(δ-ρ)?Bint，其中B?為全局場，對等式兩邊分別作傅里葉變換F，再令C=F(δ-ρ)則可得[7]：

將式（3）及偶極場與磁化率公式FBint=D?Fχ代入式（2）可得非線性全局場QSM 重建算法（ΝLTFI）的表達式：

式（4）采用ADMM方法求解得：

初始化為z0=FHDFχ+s，則上式可有效重建QSM.

1.2 結(jié)構(gòu)補償?shù)姆蔷€性全局場反演

QSM重建方法通常提取幅值圖像的梯度值作為先驗稀疏約束，但由于磁化率圖幅值和相位的邊界并非完全一致，這使得基于幅值形態(tài)先驗的磁化率重建結(jié)果的準(zhǔn)確性不足[8]；另一方面，偶極子在54.7°附近區(qū)域的雙錐面上存在0值或者接近于0的數(shù)值，這使得由局部場圖重建磁化率圖時易產(chǎn)生磁化率偽影，對這些區(qū)域進行k空間的數(shù)據(jù)補償可有效抑制偽影的產(chǎn)生[9].采用正則化的磁化率重建方法，則對組織結(jié)構(gòu)邊緣產(chǎn)生過平滑現(xiàn)象，而若用無約束條件的代價函數(shù)可提取組織結(jié)構(gòu)邊緣，可有效補償被平滑的組織結(jié)構(gòu)[10].基于此，本文對提出的ΝLTFI 算法進行結(jié)構(gòu)補償（Structural Compensation，SC），記為SC-ΝLTFI，設(shè)已經(jīng)提取組織結(jié)構(gòu)邊緣的磁化率為χs，則代價函數(shù)可表示為：

其中W為與幅值有關(guān)的權(quán)重矩陣，F(xiàn)為傅里葉變換，C=F(δ?ρ)，B?為全局場，D為偶極子在k 空間的表示形式，式（6）采用共軛梯度法求解，迭代次數(shù)為ns=100時，可求得組織結(jié)構(gòu)的細(xì)節(jié)；定義數(shù)值大于0.15的區(qū)域雙錐面為掩模版Mcone，則聯(lián)合式（5）和（6）對ΝLTFI得到的磁化率圖進行結(jié)構(gòu)補償，得到如下表達式：

2 實驗方法

截斷k空間閾值法（TKD）[11]、MEDI和FAΝSI是經(jīng)典的定量磁化率重建算法，本文以這三種方法作為對比算法，通過數(shù)值仿真人腦的QSM 驗證本文提出的ΝLTFI 及SC-ΝLTFI 的有效性，以均方根誤差（RMSE）與均值結(jié)構(gòu)相似度（MSSIM）作為評估度量.三種對比算法均采用復(fù)雜調(diào)和相位去除方法（SHARP）做背景場去除[7]，卷積核半徑為6mm.實驗中所有的代碼均運行于Matlab 2013a，電腦主機配置為Intel(R)CoreTM i7-4790 CPU，8GB RAM.

2.1 數(shù)值仿真

本文實驗?zāi)M的磁共振主磁場強度為3.0T.如圖1所示，仿真人腦的矩陣大小為224×224×150，分辨率為(1×1×1)mm3，根據(jù)實際人腦組織的解剖結(jié)構(gòu)設(shè)定磁化率值：蒼白球（GP），尾狀核（CΝ），殼核（PT），紅核（RΝ）和黑質(zhì)（SΝ）分別為0.19×10-6，0.05×10-6，0.10×10-6，0.14×10-6和0.16×10-6；血管，灰質(zhì)和白質(zhì)分別為0.3×10-6，0.02×10-6和-0.015×10-6；其他腦組織的磁化率設(shè)0，同時，在腦組織外部增加磁化率值為9.4×10-6結(jié)構(gòu)用于模擬鼻竇產(chǎn)生背景場.全局場利用前向法計算磁化率與偶極子在k空間中的乘積求得，然后用一個二值化的人腦掩模版提取感興趣區(qū)域（VOI）內(nèi)的全局場圖，如圖1b所示，最后通過模擬的回波時間TE=24 ms將場圖轉(zhuǎn)換為相位圖.幅值圖根據(jù)大腦成像的T1，T*2及質(zhì)子加權(quán)密度圖計算得到，接下來將幅值圖與相位圖組合成復(fù)數(shù)形式的信號，同時加入標(biāo)準(zhǔn)差為一的零均值高斯隨機噪聲以仿真實際成像中可能出現(xiàn)的干擾.

圖1 仿真磁化率圖及對應(yīng)的全局場圖的三方向視圖Fig.1 Three views of simulated QSM and the corresponding total field maps

3 實驗結(jié)果與討論

首先，將TKD，MEDI，F(xiàn)AΝSI和ΝLTFI、SC-ΝLTFI的相關(guān)正則化參數(shù)進行優(yōu)化.圖2a顯示TKD算法的截斷閾值（Th）與RMSE、MSSIM的關(guān)系曲線，隨著Th值的增加，RMSE值減少而MSSIM相應(yīng)增加并趨于收斂，由兩個曲線圖可知Th=0.18為最優(yōu)值.MEDI的正則化參數(shù)λMEDI優(yōu)化如圖2b所示，當(dāng)λMEDI=250時，結(jié)果最優(yōu).FAΝSI與ΝLTFI具有相同的優(yōu)化參數(shù)λ和μ，由于二者為倍數(shù)關(guān)系，令μ=b*λ，進而將問題轉(zhuǎn)化為優(yōu)化為λ 和b.實驗中，我們分別取λ值為：[0.1，0.04，0.02，0.01，0.004，1.6×10-3，8.0×10-4，4.0×10-4，2.0×10-4，1.6×10-4，1.0×10-4，2×10-5，2.0×10-6]，b 為[1，10，50，100，200].由圖3c 可知，F(xiàn)AΝSI 中λ 和b 的最優(yōu)值為λ=8.0×10-4，b=10；根據(jù)圖2d，ΝLTFI及SC-ΝLTFI在λ=1.6×10-3，b=10時，取得最佳結(jié)果.

圖2 仿真實驗中TKD，MEDI，F(xiàn)AΝSI與ΝLTFI算法的相關(guān)正則化參數(shù)優(yōu)化Fig.2 Optimization of regularization parameters of TKD，MEDI，F(xiàn)AΝSI and ΝLTFI

圖3b-f顯示相關(guān)算法的重建結(jié)果和對比實驗；圖3g顯示相應(yīng)的差值圖，采用將重建結(jié)果減去仿真圖的方法.如圖3b中紅色箭頭所示，由于TKD采用簡單的閾值截斷，導(dǎo)致k空間中數(shù)據(jù)的不連續(xù)，重建的磁化率在基底神經(jīng)節(jié)和血管附近具有很強的偽影.MEDI，F(xiàn)AΝSI及ΝLTFI三種方法均采用正則化方法，有效地抑制了磁化率偽影（如圖3c-e），但對大腦組織結(jié)構(gòu)均有過平滑現(xiàn)象，如3h中的輪廓曲線所示，而MEDI和FAΝSI在經(jīng)過背景場去除后，重建的磁化率值則略低于ΝLTFI.本文對ΝLTFI方法進行結(jié)構(gòu)補償后，SC-ΝLTFI的重建結(jié)果最接近于仿真值，而且在大血管位置處的偽影抑制要優(yōu)于其他幾種方法，如圖3中黃色箭頭所示.

圖3 仿真實驗重建結(jié)果及對比實驗Fig.3 Comparsion of reconstructed results with the simulated QSM

4 討論與總結(jié)

本文提出一種基于結(jié)構(gòu)補償?shù)姆蔷€性全局場反演算法，有效實現(xiàn)了由局場到磁化率的重建.通過數(shù)值仿真實驗驗證其有效性，將重建結(jié)果分別與TKD，MEDI及FAΝSI三種經(jīng)典算法作比較.我們以模擬的磁化率圖像基準(zhǔn)圖，根據(jù)相應(yīng)算法的磁化率重建結(jié)果繪制RMSE 和MSSIM 曲線進行相關(guān)參數(shù)的優(yōu)化.TKD 法易在血管周圍及基底核的部分區(qū)域留下嚴(yán)重的磁化率重建偽影；MEDI與FAΝSI對磁化率偽影的抑制效果遠(yuǎn)優(yōu)于TKD，但由于幅值圖像的組織邊緣與磁化率圖有所差異以及背景場去除進行的卷積運算可能對局部信號具有放大作用，因而仍可能產(chǎn)生磁化率偽影.本文提出的ΝLTFI 直接將局部場與全局場關(guān)系式代入磁化率重建函數(shù)中，省略背景場去除步驟，對磁化率偽影的抑制作用要優(yōu)于TKD，MEDI 和FAΝSI，且重建精度優(yōu)于這三種算法.但ΝLTFI 的重建結(jié)果與真實值之間仍存誤差，采用結(jié)構(gòu)補償法的SC-ΝLTFI 通過提取組織邊界的高頻信息可有效補償ΝLTFI 的缺陷，重建結(jié)果最接近于真實值.因此，本文提出的SC-ΝLTFI 算法，單步實現(xiàn)磁化率重建，簡化QSM 的重建過程，通過定量的計算和評估，重建效果優(yōu)于其它算法，具有潛在的臨床應(yīng)用價值.

當(dāng)然，本文算法仍存在有待改進的問題.在磁化率重建時，需將受強磁化率干擾的腦組織區(qū)域去除，這就使得大腦組織中有用的信息不夠完整，這將是我們未來研究的重點.