定量磁化率成像（QSM）可以提供組織的磁化率信息，在臨床診斷和研究中有重要價值。目前9.4 T 下高分辨活體猴腦 QSM 工作尚未見報道。本研究獲得了 9.4 T 200 μm 各向同性的超高分辨率活體猴腦QSM成像，并發現在超高分辨磁共振圖像中，QSM 固有的奇點問題會在 L2 正則化處理中嚴重發散。本文引入了k空間閾值分割法來消除圖像偽影，并發現當閾值在 0.2 到 0.3 的區間內圖像重建效果最佳。高分辨率活體猴腦 QSM 可為相關腦科學研究提供一種新的手段。
引用本文： 溫清清, 楊鴻毅, 鐘凱. 9.4 T超高分辨活體猴腦磁化率成像研究. 生物醫學工程學雜志, 2019, 36(3): 349-355. doi: 10.7507/1001-5515.201809009 復制
版權信息： ?四川大學華西醫院華西期刊社《生物醫學工程學雜志》版權所有，未經授權不得轉載、改編
引言
定量磁化率成像（quantitative susceptibility mapping，QSM）是一種基于磁敏感加權成像（susceptibility weighted imaging，SWI）發展起來的新的成像方法[1]。QSM 從相位圖中提取信息，獲得組織固有的磁化率。磁化率信息與組織中的順磁物質，如鐵的沉積等密切相關，而很多神經系統疾病都表現出大腦中鐵的過度沉積現象[2]，因而腦 QSM 定量檢測已成為近些年的研究熱點。
QSM 分析主要包含三個步驟：① 相位解纏繞；② 背景場去除；③ 磁化率重建。這三個步驟又各自包含很多不同的算法。常用的相位解纏繞算法包括質量圖引導法[3]和基于拉普拉斯方法的解纏繞[4]。復雜諧波偽影去除法[5]和偶極場投影法[6]是目前效果最好的背景場去除方法。磁化率重建是一個不適定逆問題，傅里葉域的場圖 
可由磁化率 
的離散傅里葉變換與 k 空間卷積核 
的乘積計算所得，但由于傅里葉域中圓錐曲面
的存在，不能直接用場圖與卷積核相除來計算磁化率。如何將局部場圖反演計算出磁化率圖，是 QSM 算法的一個核心問題。近年來也有多種算法相繼提出，如：貝葉斯正則化方法[7]、多方向采樣磁化率計算方法[8]、基于形態學的偶極子反演算法[9]、k 空間加權微分法[10]和 k 空間閾值分割法（thresholded K-space division，TKD）[11]等。這些方法都有著各自的優缺點及其最佳適用條件，且場強、回波時間和空間分辨率等條件都會對磁化率的計算精度產生影響[12]，研究者們需要根據實驗需求選擇最佳的算法組合。
QSM 算法目前主要用于人腦的磁化率分析，對獼猴腦磁化率的研究尚未見報道。獼猴和人的遺傳物質極為相似，獼猴的神經系統也與人最為接近，包括大腦結構、細胞構架、神經投射和功能環路[13-14]。因此，獼猴是模擬人類腦功能和腦疾病的最佳模式動物。目前很多腦疾病都采用獼猴作為研究對象，對磁化率的分析可以定量猴腦中的鐵沉積情況，從而為猴腦模型相關神經系統疾病的觀察和治療提供定量的生物學標記。人腦的三維圖像分辨率通常都低于 500 μm，且很多實驗為了獲得較高的平面分辨率而選擇較大的層厚，而我們在 9.4 T 下三維猴腦成像的分辨率已達到 200 μm。超高分辨率圖像能夠提供更精細的解剖結構，減小部分容積效應并且提高估計磁化率值的精確度[15]。目前 QSM 算法的研究工作主要關注如何能提高磁化率圖像重建速度或者減小偽影的存在[11,16]，分辨率的倍數級提高是否會對 QSM 算法帶來影響還未見系統研究。分辨率的提高會使頻域中定義的圓錐表面
上值為 0 的奇點變得更加密集，在此區域，局部場的噪聲在磁化率中會被放大，因此，可能會加重反演后的磁化率圖像存在的偽影問題[8]。針對這種情況，需要尋找合適的磁化率重建算法，以便為后續的超高分辨率猴腦 QSM 研究提供可行的方法。
本研究采集了分辨率為 200 μm 的活體猴腦三維梯度回波序列圖像，采用基于拉普拉斯方法的解纏繞算法進行相位解纏繞，之后采用復雜諧波偽影去除法去除背景場，這兩個算法都是目前效果較好且使用率極高的方法。在磁化率反演重建這一步選用 L2 范數法和不同閾值的 TKD 方法來進行，以比較并獲得適用于這種高分辨率成像且效果良好的磁化率重建算法。
1 實驗方法
1.1 QSM 重建算法
假設樣品磁化率為
，放置于 z 方向磁場
中，則會引發場擾動 
，且[5]：
 '/> |
其中，沿 z 方向單元偶極子
，
是
與
之間的角度。將式（1）變換到傅里葉域[12]：
 |
其中，
為傅里葉域坐標。
。擾動場由背景場校正之后的相位信息 
（亦可稱為局部場）決定，令
，相位信息與磁化率之間的關系[11]可以表示成如下形式：
 |
其中
，
為旋磁比，值為 42.58 MHz/T，
為梯度回波序列的回波時間，
是局部場的傅里葉變換。因此磁化率可由如下公式計算得到：
 |
由于存在
的情況，直接求解
并不可行，因此如何求解這種不適定逆問題，是磁化率重建算法的關鍵。
由于單方向采樣在實際應用中更為廣泛，本實驗中猴腦的數據為單方向采樣獲得，因此雖然多方向采樣磁化率計算方法可以有效地壓制偽影，但并不適用于本實驗。目前的單方向采樣 QSM 重建算法大致可以分為兩類：一類是基于傅里葉空間的閾值處理方法，一類是增加約束項引入先驗信息并利用最優化算法來進行求解[17]。通常認為引入平滑的先驗信息，可以使所獲得的磁化率圖像在偽影的抑制上優于基于傅里葉空間的閾值處理方法。L2 范數法作為正則化算法中的一種，屬于重建算法的第二類，這種算法所需的運算時間短且結果精確，已成功應用于人腦的磁化率圖像重建。因此，本文首先考慮選擇這種方法來處理猴腦磁化率圖像。另一方面，TKD 算法作為一種簡單直接的磁化率計算方法，屬于基于傅里葉空間的閾值處理方法，計算過程中不需要參數值擬合求解，但需要根據實驗條件選擇合適的閾值[11]。
1.1.1 L2 范數方法
將位于
處的體素看做一個具有強度
的磁偶極子，一個特定體素產生的頻率響應是所有體素的場偏移之和，可以寫成如下的矩陣形式：
 |
其中，
指標準化場移動的列向量[12]，
是磁化率，
即為單元偶極子 
在傅里葉域中的表達。采用標準 Tikhonov 正則化方法[7]求解磁化率的核心公式為：
 |
為線性算子，
是相應的正則化參數。正則化方法包含貝葉斯正則化重建方法、L1-正則化重建方法（L1 范數法）和 L2-正則化重建方法（L2 范數法），它們的核心思想基本一致，區別主要在于約束項（式（6）的第二項）的不同。L2 范數法[16]的核心公式為：
 |
其中，
為單位矩陣 
或者是從幅值圖中推導出的對角矩陣，
為三維梯度算子，
是相應的正則化參數。
1.1.2 TKD 方法
TKD 作為一種簡單直接的磁化率重建算法，采用一個常數值
來替代傅里葉域中偶極子
中 0 值附近的小的絕對值情況，即：
 |
此時，
求解公式如下：
 |
利用公式 
可以直接對式（4）進行求解。研究認為，當閾值 
取在 0.2～0.5 之間，可以有效地減小偽影[11]。
1.2 活體猴腦磁共振實驗
本實驗使用了 3 只 6 歲的健康雄性獼猴（安徽旌德縣皖南獼猴馴養繁殖專業合作社）。相關研究通過了本單位的動物實驗倫理委員會審核（實驗動物許可證號 SYXK（皖）2018-005）。所有磁共振成像實驗在 Bruker 40 cm 大口徑動物磁共振成像系統（Karlsruhe，德國）上開展，異氟烷麻醉后的獼猴以俯臥的姿勢放入磁體中心，身體與主磁場
平行，使用自制的磁共振射頻線圈進行三維梯度回波序列成像。成像參數如下：重復時間 TR = 40 ms，回波時間 TE = 15 ms，翻轉角 FA = 10°，矩陣大小為 384 × 384 × 256，視野大小 FOV 為 90 mm × 90 mm × 60 mm。
1.3 圖像處理
預處理：使用 FSL（Oxford，英國）軟件處理幅值圖獲得二進制掩膜圖。二進制掩膜圖與解纏繞之后的相位圖相乘，可去除頭骨以及大腦外部噪聲，使圖片只保留大腦部分。
QSM 的主體程序部分均在 Matlab 中實現。首先將三維相位圖與 FSL 所得的掩膜圖均轉換為 MAT 格式，以便 Matlab 可以直接調用，然后依次進行拉普拉斯相位解纏繞，復雜諧波偽影去除背景場和磁化率重建。磁化率重建采用了兩種算法：L2 范數法和 TKD，其中 TKD 算法的閾值選取 0.2、0.3、0.4、0.5 四種，以便比較獲得算法的最佳閾值。
2 實驗結果
基于拉普拉斯方法的解纏繞算法和復雜諧波偽影去除法可以得到高分辨率猴腦的相位解纏繞及局部場圖，且效果良好（見圖 1）。局部場圖中猴腦的解剖結構清晰，灰白質對比顯著。
圖1
原始圖像（幅值和相位圖）及其對應的相位解纏繞和背景場去除結果圖
Figure1.
The raw data (magnitude and phase) and corresponding calculated unwrapping phase image and background removal image
圖 2a 展示了 TKD 方法四種不同閾值情況下，傅里葉域中閾值修正后的偶極子 
的 y 方向中心層面（式（8）），修正后的偶極子 
在包含 384 × 384 × 256 個體素的網格上生成，此時的 
可由空間位置坐標計算所得，即 
，空間的中心點坐標為（192,192,128）。在選定閾值的情況下，將空間中每個點的絕對值求和并除以點的總個數 384 × 384 × 256 即可得到該閾值修正后的 
均值，實驗中共計算了閾值 
為 0.05、0.1、0.2、0.3、0.4、0.5 和 0.6 時 
的均值，均值與 
的關系如圖 2b 所示。從圖 2a 中可以看出隨著閾值 
的增大，
中的數值逐漸向 
和 
兩端聚集，0 附近的極小值數目減小，導致 
的均值增大，與圖 2b 結果相一致，而這種情況的出現會導致計算所得圖像整體的磁化率值范圍變小。因此，在閾值選擇時，要盡量控制所選閾值不要過大以避免給磁化率的精確度帶來較大誤差。
圖2
閾值
對公式（8）的影響
a. 修正偶極子的
on the equation (8)
a. the central slices of the
采用 L2 范數法重建分辨率較低（940 μm × 940 μm × 2 500 μm）的人腦磁化率圖像[16]時，磁化率圖像內部結構清晰，磁化率強的位置對比度明顯。然而，在處理分辨率達 200 μm 的高分辨率猴腦圖像時，L2 范數法計算出的圖像有著明顯的條紋偽影存在（見圖 3a），條紋偽影會影響組織結構的清晰度且會給磁化率值的計算帶來誤差。但由于 L2 范數法一直被認為結果較為準確，在本實驗中雖然有偽影的影響我們也以該方法計算的磁化率圖像作為一個對比度參照。圖 3b 展示了采用 TKD 方法進行磁化率重建，選取四種不同閾值（0.2、0.3、0.4、0.5）情況下的猴腦磁化率圖像。4 種不同閾值的 TKD 重建，圖像偽影都得到了很好的抑制，但當閾值取 0.2 或 0.3 時，圖像的對比度優于 0.4 和 0.5 的情況，且磁化率的數值與 L2 范數法結果比較接近（見表 1）。從表 1 中可以看出 TKD 算法閾值取較大時，磁化率計算值偏低，與圖 2 結果相一致。總體來看，閾值取到 0.4、0.5 時，磁化率值偏低；閾值取 0.2 時，TKD 算法對偽影的抑制表現略差于閾值為 0.3，但算法所計算的磁化率圖像對比度更接近 L2 范數法。綜合偽影抑制效果和磁化率準確度兩種因素考慮，TKD 方法計算本實驗高分辨率猴腦磁化率圖像的最佳閾值范圍為 0.2～0.3。
圖3
9.4 T 高分辨率猴腦磁化率圖像[磁化率單位 ppm（parts per million）]
a. L2 范數法所得磁化率圖；b. TKD 方法計算所得磁化率圖像
Figure3. The susceptibility maps of high resolution monkey brains at 9.4 T [susceptibility is presented in ppm (parts per million)]a. susceptibility map by L2-regularization; b. susceptibility maps by TKD
表1
不同重建算法計算所得 3 只猴腦的蒼白球與黑質的磁化率均值（ppm）
Table1.
The mean susceptibility of globus pallidus and substantia nigra of three monkey brains calculated by different reconstruction algorithms (ppm)
圖 4 展示了不同閾值下的 TKD 算法重建出的三個不同層面的猴腦磁化率圖像。不同閾值計算的磁化率圖像采用相同數值范圍來展示，隨著閾值的增大，圖像的對比度逐漸減弱，可以看出圖像的整體磁化率范圍在減小，與圖 2 的結果相符。四種不同閾值情況下偽影均得到了有效的抑制，不同結構之間的磁化率對比度明顯，超高分辨率磁化率圖像可以觀察到更精細的核團結構，例如在圖 3 所選的層面中能夠清晰地看到內蒼白球和外蒼白球之間的分隔。圖 4 中可以清晰地觀察到黑質和基底節結構。蒼白球和黑質的磁化率明顯高于周圍組織，灰白質對比明顯，這與以往人腦磁化率研究結果相一致[18]。
圖4
TKD 方法選擇不同閾值時，活體猴腦 3 個不同層面的磁化率圖
Figure4.
The representative slices of susceptibility maps of monkey brain in vivo calculated by TKD with different thresholding values
3 討論
目前 QSM 多用于處理 1.5 T 或 3 T 下的人腦圖像，這類圖像的分辨率基本上大于 500 μm，有的甚至大于 1 mm，例如 Persson 等[18]研究中的分辨率為 469 μm × 469 μm × 3 000 μm，Liu 等[19]研究中的分辨率為 500 μm × 500 μm × 500 μm，而本實驗中的猴腦分辨率接近 200 μm × 200 μm × 200 μm，高分辨率會增強磁化率成像中逆問題的不適定性帶來的偽影問題。通常認為引入平滑的先驗信息，可以更好地抑制磁化率圖像的偽影。然而本實驗發現，現有的 L2 范數法雖然引入了先驗知識，但計算出的磁化率圖像出現了嚴重的條紋偽影現象。我們推測該偽影的產生可能與高分辨率情況下奇點數目增加導致正則化方法最小二乘法擬合惡化有關，因而在后續算法選擇時考慮非擬合優化的算法。TKD 算法作為一種簡單直接的磁化率計算方法，計算過程中不需要參數值擬合求解，只需要選擇合適的閾值[11]。閾值選擇中要注意，閾值選擇過大會使得整體的磁化率值變小，而過小的閾值可能達不到足夠好的偽影抑制效果。本實驗中分別選擇 0.2、0.3、0.4 和 0.5 作為閾值進行了磁化率重建，圖像條紋偽影基本消除。同時發現當閾值取 0.2 或 0.3 時，圖像的磁化率信息更接近 L2 范數法且對比度更好。L2 范數重建算法在 9.4 T 超高分辨率圖像中的表現，說明圖像分辨率是定量磁化率成像中一個不可忽略的因素，分辨率的提高使得不適定逆問題中的奇點變得密集，可能導致很多算法需要進一步改進。
組織邊緣的一個體素可能同時包含兩種不同的組織類型[20]，這種部分容積效應會影響大腦中各種核團結構的分割和磁化率分析。高分辨率的圖像可以看到更多的解剖結構，例如低分辨率圖中無法觀察到的內蒼白球、外蒼白球和它們之間的分隔，有助于提高手工精細核團分割的準確度，從而得到更精確的核團磁化率值，因此高分辨率的磁化率成像是未來的一個重要研究方向。
很多中樞神經性系統疾病都會表現出磁化率值的異常，建立一個正常猴腦鐵沉積（磁化率）與年齡的關系圖，可以為猴腦的神經系統疾病研究提供一個基線參考，同時對磁化率的監控也對疾病研究過程中的診斷和治療有著重要的意義。這也是本課題組未來的一個研究目標。
4 結論
本研究對 9.4 T 下的高分辨率活體猴腦進行了定量磁化率成像研究，發現高分辨率會導致定量磁化率成像 L2 范數法的磁化率圖像條紋偽影增大，而 TKD 方法可以有效地消除偽影，其中閾值選擇 0.2～0.3 時重建效果最好。本文為高分辨率猴腦磁化率改變相關的疾病評估提供了一種可行的分析方法。
引言
定量磁化率成像（quantitative susceptibility mapping，QSM）是一種基于磁敏感加權成像（susceptibility weighted imaging，SWI）發展起來的新的成像方法[1]。QSM 從相位圖中提取信息，獲得組織固有的磁化率。磁化率信息與組織中的順磁物質，如鐵的沉積等密切相關，而很多神經系統疾病都表現出大腦中鐵的過度沉積現象[2]，因而腦 QSM 定量檢測已成為近些年的研究熱點。
QSM 分析主要包含三個步驟：① 相位解纏繞；② 背景場去除；③ 磁化率重建。這三個步驟又各自包含很多不同的算法。常用的相位解纏繞算法包括質量圖引導法[3]和基于拉普拉斯方法的解纏繞[4]。復雜諧波偽影去除法[5]和偶極場投影法[6]是目前效果最好的背景場去除方法。磁化率重建是一個不適定逆問題，傅里葉域的場圖 可由磁化率 的離散傅里葉變換與 k 空間卷積核 的乘積計算所得，但由于傅里葉域中圓錐曲面的存在，不能直接用場圖與卷積核相除來計算磁化率。如何將局部場圖反演計算出磁化率圖，是 QSM 算法的一個核心問題。近年來也有多種算法相繼提出，如：貝葉斯正則化方法[7]、多方向采樣磁化率計算方法[8]、基于形態學的偶極子反演算法[9]、k 空間加權微分法[10]和 k 空間閾值分割法（thresholded K-space division，TKD）[11]等。這些方法都有著各自的優缺點及其最佳適用條件，且場強、回波時間和空間分辨率等條件都會對磁化率的計算精度產生影響[12]，研究者們需要根據實驗需求選擇最佳的算法組合。
QSM 算法目前主要用于人腦的磁化率分析，對獼猴腦磁化率的研究尚未見報道。獼猴和人的遺傳物質極為相似，獼猴的神經系統也與人最為接近，包括大腦結構、細胞構架、神經投射和功能環路[13-14]。因此，獼猴是模擬人類腦功能和腦疾病的最佳模式動物。目前很多腦疾病都采用獼猴作為研究對象，對磁化率的分析可以定量猴腦中的鐵沉積情況，從而為猴腦模型相關神經系統疾病的觀察和治療提供定量的生物學標記。人腦的三維圖像分辨率通常都低于 500 μm，且很多實驗為了獲得較高的平面分辨率而選擇較大的層厚，而我們在 9.4 T 下三維猴腦成像的分辨率已達到 200 μm。超高分辨率圖像能夠提供更精細的解剖結構，減小部分容積效應并且提高估計磁化率值的精確度[15]。目前 QSM 算法的研究工作主要關注如何能提高磁化率圖像重建速度或者減小偽影的存在[11,16]，分辨率的倍數級提高是否會對 QSM 算法帶來影響還未見系統研究。分辨率的提高會使頻域中定義的圓錐表面上值為 0 的奇點變得更加密集，在此區域，局部場的噪聲在磁化率中會被放大，因此，可能會加重反演后的磁化率圖像存在的偽影問題[8]。針對這種情況，需要尋找合適的磁化率重建算法，以便為后續的超高分辨率猴腦 QSM 研究提供可行的方法。
本研究采集了分辨率為 200 μm 的活體猴腦三維梯度回波序列圖像，采用基于拉普拉斯方法的解纏繞算法進行相位解纏繞，之后采用復雜諧波偽影去除法去除背景場，這兩個算法都是目前效果較好且使用率極高的方法。在磁化率反演重建這一步選用 L2 范數法和不同閾值的 TKD 方法來進行，以比較并獲得適用于這種高分辨率成像且效果良好的磁化率重建算法。
1 實驗方法
1.1 QSM 重建算法
假設樣品磁化率為，放置于 z 方向磁場中，則會引發場擾動 ，且[5]：
| '/> |
其中，沿 z 方向單元偶極子，是與之間的角度。將式（1）變換到傅里葉域[12]：
|
其中，為傅里葉域坐標。。擾動場由背景場校正之后的相位信息 （亦可稱為局部場）決定，令，相位信息與磁化率之間的關系[11]可以表示成如下形式：
|
其中，為旋磁比，值為 42.58 MHz/T，為梯度回波序列的回波時間，是局部場的傅里葉變換。因此磁化率可由如下公式計算得到：
|
由于存在的情況，直接求解并不可行，因此如何求解這種不適定逆問題，是磁化率重建算法的關鍵。
由于單方向采樣在實際應用中更為廣泛，本實驗中猴腦的數據為單方向采樣獲得，因此雖然多方向采樣磁化率計算方法可以有效地壓制偽影，但并不適用于本實驗。目前的單方向采樣 QSM 重建算法大致可以分為兩類：一類是基于傅里葉空間的閾值處理方法，一類是增加約束項引入先驗信息并利用最優化算法來進行求解[17]。通常認為引入平滑的先驗信息，可以使所獲得的磁化率圖像在偽影的抑制上優于基于傅里葉空間的閾值處理方法。L2 范數法作為正則化算法中的一種，屬于重建算法的第二類，這種算法所需的運算時間短且結果精確，已成功應用于人腦的磁化率圖像重建。因此，本文首先考慮選擇這種方法來處理猴腦磁化率圖像。另一方面，TKD 算法作為一種簡單直接的磁化率計算方法，屬于基于傅里葉空間的閾值處理方法，計算過程中不需要參數值擬合求解，但需要根據實驗條件選擇合適的閾值[11]。
1.1.1 L2 范數方法
將位于處的體素看做一個具有強度的磁偶極子，一個特定體素產生的頻率響應是所有體素的場偏移之和，可以寫成如下的矩陣形式：
|
其中， 指標準化場移動的列向量[12]， 是磁化率， 即為單元偶極子 在傅里葉域中的表達。采用標準 Tikhonov 正則化方法[7]求解磁化率的核心公式為：
|
為線性算子， 是相應的正則化參數。正則化方法包含貝葉斯正則化重建方法、L1-正則化重建方法（L1 范數法）和 L2-正則化重建方法（L2 范數法），它們的核心思想基本一致，區別主要在于約束項（式（6）的第二項）的不同。L2 范數法[16]的核心公式為：
|
其中， 為單位矩陣 或者是從幅值圖中推導出的對角矩陣， 為三維梯度算子， 是相應的正則化參數。
1.1.2 TKD 方法
TKD 作為一種簡單直接的磁化率重建算法，采用一個常數值來替代傅里葉域中偶極子中 0 值附近的小的絕對值情況，即：
|
此時， 求解公式如下：
|
利用公式 可以直接對式（4）進行求解。研究認為，當閾值 取在 0.2～0.5 之間，可以有效地減小偽影[11]。
1.2 活體猴腦磁共振實驗
本實驗使用了 3 只 6 歲的健康雄性獼猴（安徽旌德縣皖南獼猴馴養繁殖專業合作社）。相關研究通過了本單位的動物實驗倫理委員會審核（實驗動物許可證號 SYXK（皖）2018-005）。所有磁共振成像實驗在 Bruker 40 cm 大口徑動物磁共振成像系統（Karlsruhe，德國）上開展，異氟烷麻醉后的獼猴以俯臥的姿勢放入磁體中心，身體與主磁場平行，使用自制的磁共振射頻線圈進行三維梯度回波序列成像。成像參數如下：重復時間 TR = 40 ms，回波時間 TE = 15 ms，翻轉角 FA = 10°，矩陣大小為 384 × 384 × 256，視野大小 FOV 為 90 mm × 90 mm × 60 mm。
1.3 圖像處理
預處理：使用 FSL（Oxford，英國）軟件處理幅值圖獲得二進制掩膜圖。二進制掩膜圖與解纏繞之后的相位圖相乘，可去除頭骨以及大腦外部噪聲，使圖片只保留大腦部分。
QSM 的主體程序部分均在 Matlab 中實現。首先將三維相位圖與 FSL 所得的掩膜圖均轉換為 MAT 格式，以便 Matlab 可以直接調用，然后依次進行拉普拉斯相位解纏繞，復雜諧波偽影去除背景場和磁化率重建。磁化率重建采用了兩種算法：L2 范數法和 TKD，其中 TKD 算法的閾值選取 0.2、0.3、0.4、0.5 四種，以便比較獲得算法的最佳閾值。
2 實驗結果
基于拉普拉斯方法的解纏繞算法和復雜諧波偽影去除法可以得到高分辨率猴腦的相位解纏繞及局部場圖，且效果良好（見圖 1）。局部場圖中猴腦的解剖結構清晰，灰白質對比顯著。
圖1
原始圖像（幅值和相位圖）及其對應的相位解纏繞和背景場去除結果圖
Figure1.
The raw data (magnitude and phase) and corresponding calculated unwrapping phase image and background removal image
圖 2a 展示了 TKD 方法四種不同閾值情況下，傅里葉域中閾值修正后的偶極子 的 y 方向中心層面（式（8）），修正后的偶極子 在包含 384 × 384 × 256 個體素的網格上生成，此時的 可由空間位置坐標計算所得，即 ，空間的中心點坐標為（192,192,128）。在選定閾值的情況下，將空間中每個點的絕對值求和并除以點的總個數 384 × 384 × 256 即可得到該閾值修正后的 均值，實驗中共計算了閾值 為 0.05、0.1、0.2、0.3、0.4、0.5 和 0.6 時 的均值，均值與 的關系如圖 2b 所示。從圖 2a 中可以看出隨著閾值 的增大， 中的數值逐漸向 和 兩端聚集，0 附近的極小值數目減小，導致 的均值增大，與圖 2b 結果相一致，而這種情況的出現會導致計算所得圖像整體的磁化率值范圍變小。因此，在閾值選擇時，要盡量控制所選閾值不要過大以避免給磁化率的精確度帶來較大誤差。
圖2
閾值對公式（8）的影響
a. 修正偶極子的
on the equation (8)
a. the central slices of the
采用 L2 范數法重建分辨率較低（940 μm × 940 μm × 2 500 μm）的人腦磁化率圖像[16]時，磁化率圖像內部結構清晰，磁化率強的位置對比度明顯。然而，在處理分辨率達 200 μm 的高分辨率猴腦圖像時，L2 范數法計算出的圖像有著明顯的條紋偽影存在（見圖 3a），條紋偽影會影響組織結構的清晰度且會給磁化率值的計算帶來誤差。但由于 L2 范數法一直被認為結果較為準確，在本實驗中雖然有偽影的影響我們也以該方法計算的磁化率圖像作為一個對比度參照。圖 3b 展示了采用 TKD 方法進行磁化率重建，選取四種不同閾值（0.2、0.3、0.4、0.5）情況下的猴腦磁化率圖像。4 種不同閾值的 TKD 重建，圖像偽影都得到了很好的抑制，但當閾值取 0.2 或 0.3 時，圖像的對比度優于 0.4 和 0.5 的情況，且磁化率的數值與 L2 范數法結果比較接近（見表 1）。從表 1 中可以看出 TKD 算法閾值取較大時，磁化率計算值偏低，與圖 2 結果相一致。總體來看，閾值取到 0.4、0.5 時，磁化率值偏低；閾值取 0.2 時，TKD 算法對偽影的抑制表現略差于閾值為 0.3，但算法所計算的磁化率圖像對比度更接近 L2 范數法。綜合偽影抑制效果和磁化率準確度兩種因素考慮，TKD 方法計算本實驗高分辨率猴腦磁化率圖像的最佳閾值范圍為 0.2～0.3。
圖3
9.4 T 高分辨率猴腦磁化率圖像[磁化率單位 ppm（parts per million）]
a. L2 范數法所得磁化率圖；b. TKD 方法計算所得磁化率圖像
Figure3. The susceptibility maps of high resolution monkey brains at 9.4 T [susceptibility is presented in ppm (parts per million)]a. susceptibility map by L2-regularization; b. susceptibility maps by TKD
表1
不同重建算法計算所得 3 只猴腦的蒼白球與黑質的磁化率均值（ppm）
Table1.
The mean susceptibility of globus pallidus and substantia nigra of three monkey brains calculated by different reconstruction algorithms (ppm)
圖 4 展示了不同閾值下的 TKD 算法重建出的三個不同層面的猴腦磁化率圖像。不同閾值計算的磁化率圖像采用相同數值范圍來展示，隨著閾值的增大，圖像的對比度逐漸減弱，可以看出圖像的整體磁化率范圍在減小，與圖 2 的結果相符。四種不同閾值情況下偽影均得到了有效的抑制，不同結構之間的磁化率對比度明顯，超高分辨率磁化率圖像可以觀察到更精細的核團結構，例如在圖 3 所選的層面中能夠清晰地看到內蒼白球和外蒼白球之間的分隔。圖 4 中可以清晰地觀察到黑質和基底節結構。蒼白球和黑質的磁化率明顯高于周圍組織，灰白質對比明顯，這與以往人腦磁化率研究結果相一致[18]。
圖4
TKD 方法選擇不同閾值時，活體猴腦 3 個不同層面的磁化率圖
Figure4.
The representative slices of susceptibility maps of monkey brain in vivo calculated by TKD with different thresholding values
3 討論
目前 QSM 多用于處理 1.5 T 或 3 T 下的人腦圖像，這類圖像的分辨率基本上大于 500 μm，有的甚至大于 1 mm，例如 Persson 等[18]研究中的分辨率為 469 μm × 469 μm × 3 000 μm，Liu 等[19]研究中的分辨率為 500 μm × 500 μm × 500 μm，而本實驗中的猴腦分辨率接近 200 μm × 200 μm × 200 μm，高分辨率會增強磁化率成像中逆問題的不適定性帶來的偽影問題。通常認為引入平滑的先驗信息，可以更好地抑制磁化率圖像的偽影。然而本實驗發現，現有的 L2 范數法雖然引入了先驗知識，但計算出的磁化率圖像出現了嚴重的條紋偽影現象。我們推測該偽影的產生可能與高分辨率情況下奇點數目增加導致正則化方法最小二乘法擬合惡化有關，因而在后續算法選擇時考慮非擬合優化的算法。TKD 算法作為一種簡單直接的磁化率計算方法，計算過程中不需要參數值擬合求解，只需要選擇合適的閾值[11]。閾值選擇中要注意，閾值選擇過大會使得整體的磁化率值變小，而過小的閾值可能達不到足夠好的偽影抑制效果。本實驗中分別選擇 0.2、0.3、0.4 和 0.5 作為閾值進行了磁化率重建，圖像條紋偽影基本消除。同時發現當閾值取 0.2 或 0.3 時，圖像的磁化率信息更接近 L2 范數法且對比度更好。L2 范數重建算法在 9.4 T 超高分辨率圖像中的表現，說明圖像分辨率是定量磁化率成像中一個不可忽略的因素，分辨率的提高使得不適定逆問題中的奇點變得密集，可能導致很多算法需要進一步改進。
組織邊緣的一個體素可能同時包含兩種不同的組織類型[20]，這種部分容積效應會影響大腦中各種核團結構的分割和磁化率分析。高分辨率的圖像可以看到更多的解剖結構，例如低分辨率圖中無法觀察到的內蒼白球、外蒼白球和它們之間的分隔，有助于提高手工精細核團分割的準確度，從而得到更精確的核團磁化率值，因此高分辨率的磁化率成像是未來的一個重要研究方向。
很多中樞神經性系統疾病都會表現出磁化率值的異常，建立一個正常猴腦鐵沉積（磁化率）與年齡的關系圖，可以為猴腦的神經系統疾病研究提供一個基線參考，同時對磁化率的監控也對疾病研究過程中的診斷和治療有著重要的意義。這也是本課題組未來的一個研究目標。
4 結論
本研究對 9.4 T 下的高分辨率活體猴腦進行了定量磁化率成像研究，發現高分辨率會導致定量磁化率成像 L2 范數法的磁化率圖像條紋偽影增大，而 TKD 方法可以有效地消除偽影，其中閾值選擇 0.2～0.3 時重建效果最好。本文為高分辨率猴腦磁化率改變相關的疾病評估提供了一種可行的分析方法。
