磁共振擴散張量成像中擴散敏感梯度磁場方向分布方案的研究進展*

劉良友 高嵩 ? 李莎 李兆同 夏一帆

1) (北京大學醫學部醫學技術研究院, 北京 100191)

2) (北京大學醫學技術研究發展中心, 北京 100191)

磁共振擴散張量成像可以定量無創研究人體內水分子在三維空間中的各向異性擴散規律, 進而獲取重要的病理及生理信息.為了得到水分子各向異性擴散信息, 需要按照一定的方案依次施加不同方向的擴散敏感梯度磁場, 測量水分子在這些方向上的擴散系數用以估算擴散張量.擴散張量成像測量結果的準確程度受梯度磁場方向分布方案的影響, 本文對擴散敏感梯度磁場方向分布方案進行綜述, 包括完全隨機方案、啟發式方案、規則多面體式方案和數值優化方案等, 分析這些方案的優勢與局限性, 并提出需進一步研究的問題.

綜述

1 引 言

在均勻介質中水分子的隨機擴散現象表現為各向同性, 其擴散程度可以用擴散系數D表示.在復雜的生物組織中, 水分子的擴散現象受組織結構的影響, 表現為各向異性擴散, 擴散程度用擴散張量D表示[1].擴散張量成像(diffusion tensor imaging,DTI)[2,3]技術是在標準磁共振成像(magnetic reso?nance imaging, MRI)脈沖序列中加入不同方向的擴散敏感梯度磁場(diffusion sensitive gradient,DSG), 測量不同方向上水分子的擴散情況, 最終獲得人體內大量的生理和病理信息.DTI技術提供了一種在微觀結構尺度上以非侵入性的方法來表征人體組織的結構和功能特性, 已經成為了解正常大腦組織結構[4,5]以及神經和精神性疾病病變過程的重要方式[6,7].

為了使DTI結果更加準確, 需要施加數量眾多、線性無關且空間均勻分布的DSG.近年來為了進一步獲取人體內復雜的微觀結構, 提出了高角分辨率擴散成像(high angular resolution diffusion imaging, HARDI)[8]技術, 對DSG的數量和方向分布的均勻性提出了更高的要求.本文從DTI基本原理出發, 介紹了多種不同的DSG方向分布方案及特點, 同時提出DSG方向分布方案應需進一步優化的問題.

2 DTI基本原理

2.1 Bloch-Torrey方程

為了描述擴散現象與磁化強度M隨時間的變化, Torrey在Bloch方程中加入擴散項得到關于擴散磁共振的方程[9]：

圖1 Stejskal?Tanner序列, 其中兩個同等的DSG脈沖置于180°RF脈沖兩側, 強度為G, 寬度為δ, 時間間隔為ΔFig.1.Stejskal?Tanner Scheme：Two diffusion sensitive gradients inserted before and after 180° RF refocusing pulse.G, amplitude;δ, duration of the DSG; Δ, time between the two sensitive gradient lobes.

其中B為磁場強度; γ為旋磁比; Mx, My, Mz是磁化矢量的3個分量; M0是平衡態下的磁化強度;T1, T2是縱向和橫向弛豫時間.(1)式右端的前三項是經典的Bloch方程, 當施加 90°射頻(radio frequency, RF)脈沖時, (1)式在橫向平面上的磁化強度Mxy的解為

其中b為描述施加的DSG強度(G)、持續時間等因素的擴散梯度因子, 單位是s/mm2, 值越大表示擴散權重越大, 對水分子擴散更加敏感, 可以表示為

Stejskal?Tanner序列是利用水分子信號衰減的特性來測量水分子的擴散[10,11], 它是將一對大小和持續時間完全相同的DSG沿著某一方向置于180°脈沖的兩側(圖1), 其信號衰減的大小可以由(1)—(3)式求解[12]：

其中S, S0表示施加和未施加DSG的回波信號值;Δ為兩個脈沖之間的時間間隔; δ為脈沖的持續時間; D單位是mm2/s, 值越大表示水分子擴散運動的能力越強.

2.2 擴散張量成像

引入擴散張量D描述人體內水分子沿各向異性擴散規律[1], 即：

擴散張量矩陣為對稱陣, 即Dij= Dji, i, j = x, y,z, 其中對角線元素Dxx, Dyy, Dzz和非對角線元素Dxy, Dxz, Dyz是擴散張量的6個獨立分量, 表示水分子在不同方向上的擴散系數[1].

由(4)式可知, 對于各向異性擴散, 可以得到DTI圖像在不同DSG下S與S0的關系：

其中gi= (gxigyigzi)是歸一化矢量, 表示3個正交方向上的DSG強度, 等式左端為表觀擴散系數(apparent diffusion coefficient, ADC), (6)式展開為

其中i代表DSG方向的數量, ADCi表示第i個方向的表觀擴散系數, 張量D有6個獨立分量, 因此需要在至少6個不共線的方向上施加DSG, 并代入方程(7)得到張量D.為了得到更加準確的結果, 并用于計算擴散張量D, 同時減少噪聲對結果的影響, 通常會在眾多的空間方向上施加DSG[13],而HARDI技術對DSG方向數量要求則更多.

通過對擴散張量矩陣對角化[14]可以得到擴散張量的特征值及特征矢量.特征值和特征矢量是旋轉不變量, 且每個特征值對應一個特征矢量, 特征矢量之間相互垂直, 因此對于擴散的三維特性可以使用橢球體來描述, 橢球體的三個軸由擴散張量的三個特征值 λ1, λ2, λ3和特征矢量 ε1, ε2, ε3確定[15].同時可以利用特征值定義多種各向異性參數指標[16,17], 包括各向異性分數(fractional anisotropy,FA)、相對各向異性(relative anisotropy, RA)和容積比(volume ratio, VR)等, 還包括評價整體擴散效應指數的張量跡(trace)、平均擴散系數(medium diffusion, MD)等.

3 DSG方向分布方案

在DTI中, 通常在半個或整個球面上沿著至少6個線性無關的方向施加DSG, 并且要求DSG方向分布盡可能均勻, 才能獲得較為準確的擴散張量D[18,19].一個空間均勻的DSG方向分布方案對于提高張量估算的精確性有著重要作用[20?22], 在一定程度上增加DSG的數量也有利于提高圖像的信噪比(signal noise ratio, SNR)[23].HARDI技術對DSG的分布方案有更高要求, 但隨著DSG數量的增多, 圖像的采集時間也在不斷增加, 因此大量研究提出了多種DSG方向分布方案[21,22].

3.1 隨機分布方案

一般來說, 對于感興趣區內組織結構的方向分布并不是已知的, 通常假設最佳的DSG方向是均勻分布在單位球面上[21], 且所有的DSG方向之間不共線.隨機分布的DSG方向分布方案[24]即是在單位球表面隨機分布一定數量的DSG方向, 按照隨機方案產生的次序依次施加DSG, 最終獲得隨機分布方案的擴散信號.該方案雖然可以獲得任意數量的DSG分布方向, 但是分布的均勻性較差(圖2), 不能滿足DSG方向分布盡可能均勻的要求, 并且隨機分布的方案存在無法復現的問題, 因此隨機分布是一種較差的DSG方向分布方案.

3.2 啟發式分布方案

啟發式DSG方向分布方案是基于MRI系統的x, y, z梯度方向來定義的一組基本方向, 對應著邊長等于2, 中心點在(0, 0, 0)的立方體的三個面中心, 可以生成最多13個非共線的DSG方向(gx, gy, gz), 分別對應于面中心

邊中心線

立方體對角線

圖2 隨機分布60個方向DSG方向分布方案Fig.2.DSG encoding scheme with random distribution in 60 directions.

啟發式DSG方向分布方案包括：1)正交DSG方向分布方案, 由三個正交軸和三個邊中心線組合成6個DSG分布方向, 即(G1+ G2), 稱為金字塔錐體編碼, 常用于化學屏蔽張量譜學中[25]; 2)傾斜雙DSG方向分布方案, 由Basser和Pierpaoli[26]提出將所有的邊中心線產生 6個DSG方向進行組合, 即(G2+ G3), 同樣可以應用于化學屏蔽張量譜學中; 3)正交軸和四面體組合DSG方向分布方案, 由三個正交軸梯度和立方體對角線[27]梯度組合成7個方向, 即(G1+ G4), 不僅可以用于DTI[28]還可以用于熱膨脹的張量測量[29];4)十面體DSG方向分布方案, Skare 和Nordell[30]將邊中心線和立方體對角線組合成10個DSG方向, 即(G2+ G3+ G4); 5)完全啟發式DSG方向分布方案, 將所有的面中心、邊中心和對角線組合所產生的13個方向的DSG方向分布方案, 即(G1+ G2+ G3+ G4), 如圖3所示.

圖3 完全啟發式13個方向DSG方向分布方案Fig.3.DSG encoding scheme with heuristic distribution in 13 directions.

啟發式DSG方向分布方案由立方體對應的不同方向組合而成, 生成5種可供選擇的子方案, 能夠產生較為均勻的空間分布, 應用于張量測量的同時還在化學領域有較為廣泛的應用, 但是啟發式DSG方向分布方案所產生的方向數量最多僅為13個, 并且不能生成連續的方向數, 對于張量的精確估算來說, 該方案方向數量則略顯不足.

3.3 幾何多面體分布方案

利用高對稱性的幾何多面體生成DSG方向分布方案.幾何多面體是由等邊三角形、正方形或者正五邊形組成的立體幾何, 最常用的幾何多面體包括四面體、立方體、八面體、十二面體和二十面體等.幾何多面體的頂點數與面數之和比邊數多出2個, 立方體的頂點與八面體的面彼此對應, 稱為對偶, 十二面體與二十面體也具有類似的性質, 在這四種正幾何多面體中, 從一個頂點通過中心繪制的線經過另一個頂點, 從一個面通過中心繪制的線經過另一個面, 從一條邊通過中心繪制的線經過另一條邊.因此它們的頂點、面、邊可以定義一半的不共線方向, 如3, 4, 6, 10, 15個方向[31].而正四面體則是從一個頂點穿過中心經過另一個面的中心, 因此可以生成3或4個DSG方向, 稱為四面體梯度方案[27], 且3和4個DSG方向可以組合生成7個DSG方向[32,33], 這就對應了四面體所有的頂點和邊或者對應立方體、八面體的所有頂點和面, 并且可以調整其中一些方向的極性使其分布更加均勻[31].

立方體和八面體的十二條邊可以定義6個均勻分布的DSG方向用于計算擴散張量值[21,26,34?36].同樣立方體的六個面和十二條邊可以組成9個DSG方向, 并且立方體的十二條邊可以被二十面體的十二個頂點取代產生新的9個DSG方向, 后者的空間均勻性更好[37].13個DSG方向可以來自于立方體所有的面、頂點、邊的組合, 其中任何一組的方向也可以由其他多面體DSG方向所取代或改變方向極性產生新的更加均勻的DSG方向分布.

二十面體包含的十二個頂點可以得到6個DSG方向分布, 二十個面可以得到10個DSG方向[21,31,36], 三十條邊可以得到15個DSG方向[31,34].若6, 10, 15個梯度為單一方向, 則可以相互組合成為數量更多的分布均勻的DSG方案, 如6, 10,15, 16, 21, 25, 31等方向, 在這些方向組合中, 有的需要對一個或多個方向進行極性調整來獲取最佳的DSG分布.對于一些其他數量的方向如12個方向, 則可以通過沿著x, y或z軸旋轉90°來生成,或者通過刪除二十面體15個方向中的立方體面產生的3個方向來獲得.其中二十面體的31個方向如圖4所示.

圖4 二十面體31個方向DSG方向分布方案Fig.4.DSG encoding scheme with icosahedron distribution in 31 directions.

將二十面體每條邊的中點作為新的頂點, 并將每個新的頂點相連, 來生成Ne= 5n2+ 1個方向,其中 n = 1, 2, 3, ··[38,39], 該方案可以產生 6, 21,46, 81, 126等方向數量.幾何多面體能夠生成較為均勻的DSG方向分布方案, 但是方向數量不多且不連續.

3.4 DISCOBALL分布方案

與DSG方向均勻的分布在球面上相比, 將DSG方向均勻的分布在圓環上更為簡單.Stirnberg等[40]將三維的球面分布簡化為二維的圓環分布,首先利用恒定的天頂角增量將球面分成均勻的片層, 然后將片層以恒定的方位角增量進行切割.每個片層的DSG方向的分布數量為：2Nsin(θ), 其中N為預定義的片層數, θ為天頂角增量, 同時將數量的結果四舍五入, 并且可以由

根據預定義的DSG方向數量Ntotal生成片層數來計算各片層上的DSG方向數, 最終得到整個球面的DSG方向的分布, 如圖5所示.

該方案可以產生均勻分布的DSG方向, 但是子集分布的均勻性較差, 對于不完整掃描的情況容易對張量產生錯誤的估算結果.

圖5 DISCOBALL 30個方向DSG方向分布方案Fig.5.DSG encoding scheme with DISCOBALL distribu?tion in 30 directions.

3.5 球面螺旋分布方案

基于單位球面恒定角速度采樣的DSG方向分布方案會導致赤道位置DSG方向分布的密度較低, 而靠近兩極處則密度較高, 對于盡可能均勻分布的空間DSG方向來說, 恒定角速度的DSG分布方案顯然無法滿足要求.Wong和Roos[41]提出球面螺旋的分布方案, 該方法是利用求解球面的方位角及頂角的速率變化生成恒定速率的螺旋路徑,并掃過單位球面來產生均勻分布的DSG方向.DSG方向(gx, gy, gz)經歸一化后可以寫為

其中 N 為 DSG 方向數量, n = 1, 2, ··, N.此方法適用于DSG方向數量較多的分布方案, 可以提供覆蓋單位球面的方向螺旋, 由于相鄰點之間的距離相等, 因此DSG方向分布近似均勻[21,31], 如圖6所示.該方案的不足之處是容易在兩極附近產生空洞, 且DSG方向子集分布的均勻性較差.

圖6 球面螺旋分布60個方向DSG方向分布方案Fig.6.DSG encoding scheme with spherical spiral distribu?tion in 60 directions.

3.6 Jones方案

Jones方案[20,21,31,42?44], 又稱為靜電力排斥或者最小作用力算法的DSG方向分布方案, 其建立在Conturo等[27]的平衡化學sp3雜化軌道中靜電排斥模型的基礎之上, 是目前應用范圍最為廣泛的DSG方向分布方案之一.假設建立一種模型,在該模型中DSG方向都是經過球體中心的線, 在線與球體表面相交處的兩個點上放置單位點電荷,不斷改變線的方向, 使產生的每個DSG方向都用一對點來表示, 在相反方向上的DSG同樣可以用正方向上的DSG進行擴散衰減的測量.根據庫侖定律, 一對點電荷之間的排斥力與點電荷之間距離的平方成反比, 因此, 用于在三維空間中均勻排列DSG方向的算法被用來優化這些正交梯度的方向,直到所有可能的電荷對之間的作用力總和最小.該方案生成的方向見圖7(a).

該方法可以得到數量眾多方向連續且分布均勻的DSG方向分布方案, 并且對于各向異性程度較高的組織, 眾多不共線的DSG方向可以提供更為穩健的張量估算.但是由于初始方向的隨機性,該方案無法復現, 單位點電荷之間最小靜電力斥力算法也是一種迭代耗時的數值計算方法, 且很難得到一個簡單的解析式, 其最大的不足是該方案的子集分布的均勻性不高, 無法滿足一些臨床要求.

圖7 兩種Jones 60個方向DSG方向分布方案 (a) Jones方案; (b)排序的Jones方案Fig.7.DSG encoding scheme with Jones (a) and Ordered Jones (b) in 60 directions.

Dubois等[22,45]提出把圖像采集過程中位置接近的DSG方向進行分離, 以此來提高部分采集時張量估算的精確性, 稱為加權的Jones方案, 主要方法是在能量公式中引入權重的概念：

其中Eij為兩個不同方向i, j間相互作用的能量,ωij為相互作用權重, 在一個完全各向同性的方向分布方案中, ωij= 1.

Dubois等[22,45]提出三種不同的加權方案, 第一種產生DSG方向的子集, 使每個子集包含6個方向, 如果方向i, j在同一個子集中則ωij= 1, 否則, ωij= 0.6; 第二種方案是相鄰子集間的權重為ωij= 0.8; 第三種方案是不產生子集, 但是權重因子隨著序列的采集順序不斷下降, 如果|i — j| ≤10, 則 ωij= 1, 如果|i — j| ＞ 10, 則 ωij= |i — j|—α,α為常數值.

三種加權方案對應三種不同的采樣方式, 第一種采樣方式由不同方向的子集組成, 每個子集包含相同的方向數量(6個), 屬于同一個子集的方向之間是完全相互作用的, 即權重因子等于1, 而屬于不同子集的方向相互作用系數大小為0.6; 第二種采樣方式原理與第一種方式類似, 不同的是方向在不同子集間的相互作用都是不相同的, 并且后一個子集應包括前一個子集中的方向, 各個子集的時間間隔越長, 相互作用越小, 權重因子也越小, 這就要求后一個子集必須在之前的集合基礎上提供更加準確的空間方向信息; 第三種采樣方式是完成一次完整的DSG方向采樣, 不生成子集, 因此要求最少有6個DSG方向均勻分布, 在采樣模型中,某一個方向與它時間上相距最近的10個方向完全相互作用, 與其他方向的相互作用隨著時間的增加呈指數模型規律遞減.該方案采用的是模擬退火法來最小化全局能量.

加權的Jones方案可以解決原Jones方案中子集方向分布不均勻的問題, 它的劣勢在于雖然采用最小化加權能量公式產生DSG方向分布方案,但該方案卻無法滿足未加權的能量最小化, 因此其整體的DSG方向分布的均勻性不如Jones方案,另外由于能量公式中權重的變大會導致局部最小值數量的增加, 提高了求解全局最小值的難度.

Cook等[46]提出將所有的靜電點集分為相同大小的子集, 并且使每個子集中的能量最小化以達到分布的各向同性, 具體來說就是最小化能量函數：

當 i, j在同一個子集中時, δij= 1, 否則 δij= 0.該方案的目的就是建立各向同性的子集, 利用每個子集來更加精確的估算擴散張量, 得到用于運動校正的張量信息.與未經排序的方向集相比, 提高了部分掃描的質量, 但是由于該方案是在獨立處理每個子集方向, 導致對于前P個方向的優化無法達到最佳, 比如對于第一個子集和第二個子集的前一半組成的掃描序列可能會呈現各向異性較大的結果.

對于上述方法的不足, Cook等[43]再次提出一種新的方案來優化全局的采集順序, 同時在不影響整體掃描質量的情況下, 提高了部分掃描的結果.該方案稱為排序的Jones方案, 與上述方案不同的是, 子集的劃分為嵌套式, 而后者將方向劃分為非重疊的子集.目標是同時最小化所有子集的靜電能量來求解最優排序：

EP為前P個方向子集的靜電能量, 當P個方向各向同性分布時, 其靜電能近似正比于P2, 因此歸一化因子P—2使得每個子集的對目標函數具有相近的作用, 求解最優排序采用模擬退火算法, 方向如圖7(b)所示[47].

4 總結與展望

由DTI原理可知, DSG分布越均勻, 對擴散張量的估算越精確, 同時對組織結構和纖維束走行分布的描述越真實.隨機分布的DSG方向分布方案雖然可以產生數量較多且連續的DSG方向, 但DSG方向分布不均勻, 臨床上很少采用這種方案;啟發式方案彌補了隨機方案方向分布不均勻的缺點, 但其產生的方向數量較少且方向數不連續, 應用范圍較窄; 球面螺旋和DISCOBALL分布方案可以產生分布均勻且數量連續的DSG, 但是子集分布的均勻性較差, 應用較少; 幾何多面體方案,能夠產生均勻的DSG方向空間分布, 但是采樣的數量和連續性受限; 目前應用最為廣泛的是Jones DSG方向分布方案, 它能夠產生任意數量且分布均勻的DSG方向, 并且加權和排序的Jones方案也解決了子集分布不均勻的問題, 可以應對臨床采集過程中由于患者不配合或者不自主運動導致的數據集損壞的狀況, 但是由于Jones方案的隨機性, 造成梯度方案無法復現, 同時該方案采用迭代的數值計算方法, 導致計算過程復雜且耗時.

單纖維取向擴散張量模型假設每個體素中只能有一個主擴散方向, 但當纖維束出現交叉、分支、匯聚等情況時, 該模型則顯得過于簡單, 同時對于張量的估算及復雜的神經元微結構精度的判斷變得不確定.HARDI[48?50]技術則沿著更多方向施加DSG, 因此它能夠更準確解決纖維束交叉等問題, 其采樣是分布于q空間的單球殼(single?shell)或者多球殼(multi?shell)上, 相較于 q空間的大量采樣, HARDI減少了采樣的數量, 降低了采集時間, 同時對b值的要求較低, 也得到了相對較好的SNR.但由于單球殼的HARDI (sHARDI)采樣重建出的擴散概率分布函數不存在q空間的徑向信息, 因此為了能夠從徑向信息中獲取更多的擴散方向信息, 有研究提出多球殼HARDI (mHARDI)采樣技術[51], 在每個球殼上DSG方向同樣是均勻分布, 其在解決纖維束交叉等問題上的結果要優于sHARDI采樣技術.

基于以上所有DSG方向分布方案所存在的不足和需要進一步探索的問題, 可考慮將黃金分割法應用于球面DSG方向分布方案, 并采用統計學方法驗證該方案方向分布的均勻性, 進而得到DSG的均勻性對DTI結果的影響, 并且該方案既有望能夠滿足HARDI的要求, 同時在面對臨床中可能碰到的數據集損壞問題時也可以得到較準確的張量估算結果.