信號自適應衰減的多殼擴散磁共振成像方法

羅伶俐，王遠軍

(上海理工大學 醫學影像工程研究所，上海 200093)

1 引 言

擴散磁共振成像(Diffusion Magnetic Resonance Imaging，dMRI)技術被廣泛應用于諸如阿爾茨海默癥、精神分裂癥、腦損傷等諸多大腦疾病的研究中.dMRI的基礎原理是由于擴散加權MR信號對器官組織中的內生水分子的隨機運動十分敏感，這其中應用最為廣泛的擴散加權磁共振成像(Diffusion-weighted MRI，DW-MRI)是使用了兩種在180°附近的射頻脈沖自旋回波，其磁場梯度的方向和幅度在MR序列上均相同[1]，若在兩種梯度中氫核的位置不一致時，其磁矩會發生凈位移.若一群隨機運動的氫核在擴散過程中呈現的相位不一致，會導致整體信號的衰減，進而生成擴散MR信號.

此處用E(q)表示歸一化后的擴散MR信號，即是在波矢量q處所測的dMRI信號與不施加擴散靈敏梯度場時的原始MR信號之間的比值.從數學上講，E(q)與一個重要指標：總體平均擴散傳播算子(Ensemble Average diffusion Propagator，EAP)有關，其重建對于dMRI意義重大：不僅能由EAP推導出諸如用于描述纖維束方向特征的方向分布函數(Orientation Distribution Function，ODF)，其他相關細胞大小、方向，及跨膜交換的信息特征同樣可由此推導.

EAP的一種重建方法需要一個將粒子擴散過程與表示擴散信號的函數關聯起來的理想模型，但前提是信號的函數表示夠精確.這其中最經典的成像方法為擴散張量成像(Diffusion Tensor Imaging，DTI)[2]，DTI方法假設E(q)是一個以原點為中心的高斯函數，但這種過于簡化的假設對于含有復雜纖維結構(如交叉、相接的纖維束)的體素，建模局限性很大.相關研究表明，在當前dMRI成像分辨率的條件下，超過90%的體素存在復雜的多纖維群[3].為突破DTI的這一局限性，許多基于高角度分辨率擴散成像(High Angular Resolution Diffusion Imaging，HARDI)的方法被提出，這其中包括擴散光譜成像(Diffusion Spectrum Imaging，DSI)技術[4].DSI需在q空間中密集的笛卡爾網格點上進行數據采樣，在得出表示歸一化信號的連續函數E(q)后，對其進行傅里葉逆變換可得到相應的EAP[5,6]：

(1)

其中P(r)表示運動粒子凈位移為r的概率.q為q空間中的波矢量.由于擴散信號是掃描采樣所得數據，因此EAP的求解問題簡化成了如何根據所采樣的散點信號數據，估計出表示信號的連續函數E(q).

但過長的采樣時間使得DSI在臨床上并不實用.為解決這一問題，近年來的改進方法主要通過建立合適的信號模型以減少信號采樣，或用適當的插值基函數表示q空間中的擴散信號.例如，球極傅里葉基函數[7]，貝塞爾傅里葉基函數[8]，以量子力學簡單諧波振蕩器哈密頓量(亦稱為埃爾米特函數)，其本征函數的線性組合表示三維q空間中的擴散信號的MAP[9]方法.而MAPL[10]則是在MAP-MRI的基礎上，利用信號的拉普拉斯范數以對所求解的系數進行優化.但由于更多纖維組織的潛在結構信息在高b值條件下更易顯現，需將成像條件由單殼推廣到多殼條件(即多個b值條件下)時，以上方法均忽略了dMRI信號隨b值增大而衰減的相關性質，該性質為研究白質微結構提供了重要信息.Rathi等人[11]以球型脊波作為dMRI信號的插值基函數，用徑向函數項作為約束條件規范了整體信號的徑向衰減性，并結合信號函數的全變分(Total Variation，TV) 范數優化，該方法在稀疏的多殼采樣條件下實現了不錯的EAP重建效果.Ning等人[12]則提出用徑向基函數(Radial Basis Function，RBF)作為信號的插值基函數，該方法將dMRI信號表示為以q空間中不同位置為中心的各向異性高斯基函數的線性組合，并構建相應的約束條件以保證信號的衰減性，該方法也取得了理想的EAP重建效果并呈現了較強的魯棒性.然而上述方法中約束信號衰減的條件均為構建信號插值基函數的外部條件，并沒有將信號的衰減性考慮進基函數的構建中.為進一步提高EAP的重建精度，同時提升重建算法的計算效率，本文在以RBF作為頻域信號插值基函數方法的基礎上，添加了信號自適應衰減項以對信號衰減性進行建模，這不僅保留了原方法能夠明確估算EAP的二階、四階矩張量的優勢，且就體模數據的各項重建指標的結果來看，該方法能在提升計算效率的同時，進一步提高了EAP的重建精度.

2 徑向基函數在dMRI中的應用

對于徑向基函數φ(x)，x∈d，可被寫作φ(‖x‖)，‖x‖表示由原點到x點的歐式距離，即該函數在點x處的值只取決于‖x‖.有時也可用馬氏距離替代歐氏距離[13].RBF通常被寫作式(2)形式以進行函數估計：

(2)

其中s(x)作為待估計函數，可被表示成N個中心點在cn處，對應權重為ωn的RBF的線性組合.相關RBF的基礎理論[14]表明：當RBF的數量N足夠大，且中心點cn的選取適當時，則RBF幾乎可以逼近任意連續函數.

(3)

(4)

由于dMRI信號具有偶對稱性質，即E(q)=E(-q)，因此對偶項的權重保持一致：

(5)

式(5)中的第一項v0φ0(q)，可使該模型更精確地估計各方向信號衰減基本呈各向同性的信號(如大腦灰質和腦脊液區域)，而剩余項則可被視作對信號衰減呈各向異性部分的信號估計.從而更全面地對整體dMRI信號特征進行建模.此外，用于插值dMRI 頻域信號所用的RBF是一種常見的高斯核函數φ(x)=exp(-σ‖x‖2)，而使用高斯核函數作為插值基函數的優勢在于：由于高斯核函數經傅里葉變換之后依舊是高斯核函數，因此可更簡便解析地計算后續EAP、ODF以及其他相關成像統計標量.

3 EAP及基于EAP計算的相關指標

3.1 EAP的計算

由于將RBF作為q空間信號的插值基函數，則所估計的信號函數表達式E(q)是一系列高斯核函數的線性組合，因此對其進行傅里葉逆變換后所得的EAP也可對應地被寫作相關基函數的線性組合：

(6)

(7)

3.2 方向分布函數(ODF)

ODF常被用于可視化纖維組織的方向特征，可由EAP通過式(8)中的積分公式進行計算：

(8)

(9)

3.3 均方位移(Mean-Squared-Displacement，MSD)及平均四階位移(Mean Fourth order Displacement，MFD)

由于潛在纖維組織的重要特征可通過對EAP的高階矩張量計算所得，用高斯核RBF作為信號插值基函數的一大優勢在于：可解析地計算EAP的高階矩統計量，這對于發掘白質結構微小異常的研究非常重要.因此衍生出了基于二階矩張量及四階矩張量計算的成像統計標量：均方位移(MSD)，及平均四階位移(MFD).

其中MSD正比于測量時間內的水分子平均擴散量，MFD則對EAP的各向異性部分更敏感，因此可捕捉到更多潛在纖維組織的相關信息.二者計算方法為：

(10)

(11)

3.4 回歸原點概率(Return-to-the-Origin Probability，RTOP)、回歸軸向概率(Return-to-the-Axis Probability，RTAP)

用于衡量在兩個擴散靈敏梯度場間水分子凈位移的概率是dMRI中的重要成像指標，其中回歸原點概率(RTOP)的取值由P(0)決定，由式(6)、式(7)可知，RTOP的計算公式為：

(12)

(13)

4 基于RBF的dMRI信號自適應衰減模型

4.1 信號自適應衰減項

在Ning等人提出用RBF作為dMRI信號插值基函數的方法中，約束dMRI信號衰減性質的方法是：構建約束矩陣B={φKq-φK(q+1)|q=1，2，…，Q}，對于所求系數v，通過保證Bv>0以確保信號隨b值增大而衰減的性質.然而該約束矩陣的引入占據了整個求解框架90%以上的計算時間.因此為進一步提高重建算法的計算效率，本文在原有徑向插值基函數上添加了信號自適應衰減項以對信號幅度的衰減率進行建模：

(14)

(15)

信號自適應衰減項H(q)=exp(-αq2)，α>0，H(q)∈[0，1]對多殼dMRI信號進行單調遞減的指數衰減建模，這一項的引入可確保在高b值的低信噪比條件下整體信號依舊呈衰減趨勢.當α=0時，式(14)則是針對單殼條件成像的表達式.此外整個求解框架除了計算系數v，還有對H(q)中衰減系數α的估算.而少量的待估參數確保了整個算法的穩定性，以及更好的魯棒性.

4.2 參數估算

4.2.1 擴散張量計算

在計算插值基函數之前，須先用DTI方法求解二階擴散張量：

(16)

a)對式(16)兩邊取對數后，用加權線性回歸方法估算等式右邊的權重D，由此減少線性求解過程中噪聲帶來的影響[19].

b)為保證張量的正定性，對于上一步所求的權重進行平方根分解(又稱Cholesky分解)：D=U0TU0，其中U0為上三角陣.

c)將U0作為估算初值，基于Levenberg-Marquadt非線性擬合算法擬合式(17)以計算上三角陣U，最終擴散張量由D0=UTU計算所得.

(17)

4.2.2 構建基函數矩陣及約束條件

但當所采集的信號數量較少時(例如小于基函數數量N的情況)，那么方程Av=E就成了一個欠定問題，因此需對所求解的系數向量進行相關正則化以減小誤差，結合常用的正則化方法：l2范數正則(又稱蒂霍諾夫正則化)，目標函數構建為：

min‖Av-E‖2+λ‖v‖2

(18)

不過基于壓縮感知理論[19]：對于欠定問題Av=E，假設向量v的元素分布足夠稀疏，則由Av還原出的信號精度越高.因此l1范數正則化也常被用于求解該類線性逆問題.此處設η為待估計信號與實際信號之間的噪聲，則目標函數為：

(19)

4.2.3 參數求解

由于在本實驗中除了求解系數向量v，還需估算信號的衰減系數α，因此本文采用共軛梯度算法來交替地估算v和α，算法流程為：

算法1.求解系數向量v與衰減系數α

輸入：基函數矩陣A，信號向量E

輸出：系數向量v，衰減系數α

a) 初始化α

b) whileα>0

c) 根據目標函數式(18)或式(19)求解系數向量v

d) 根據‖Av-E‖2對α求導，使用梯度下降法更新α

e) ifα<0 或‖v‖1不再下降

f) break;

g) end

而對于式(19)中結合系數向量l1正則化的目標函數，本實驗中則運用交替方向乘子法(Alternating Direction Method of Multipliers，ADMM)[20]，可有效地將該目標函數分解成一系列更簡單的優化問題，此處將輔助變量Z引入式(19)，將其寫成另一種形式：

(20)

該算法求解過程中的每一次迭代均通過交替地更新參數v與Z來最小化以下拉格朗日增廣函數：

(21)

其中ρ1，ρ2，ρ3分別為對應約束項的懲罰因子，以保證所求解在可行域中.λ1，λ2，λ3則分別是約束條件AZ-E=0，cTZ-1=0，v-Z=0的乘子項.設在第i+1次迭代中，各參數的更新公式為：

(22)

(23)

乘子項的更新迭代公式為：

(24)

(25)

(26)

5 實驗方法與結果

本實驗中所用的掃描數據來源于圖1中的球型物理體模[21]，該體模與人腦內部組織具有相似的擴散特性.該體模上有兩道大小為1×0.7cm2，交叉角度為45°的溝槽.

(a) Spherical phantom (b) Baseline image of the dataset

體模數據由西門子3T掃描儀以2mm×2mm×7mm的空間分辨率對體模進行掃描所得.對于五個不同的b值，即b={1000，2000，3000，4000，5000}，每個b值條件下在81個梯度方向下掃描十次，并取十次掃描數據的均值，并用該數據計算所得的成像標量(如MSD等)，作為后續實驗數據所對比的金標準.而用于實驗的測試數據則是在同掃描條件下，對于單個b值分別取K個梯度方向，K={15，20，25，30，40，45，50，55，60}，按照掃描參數b值的高低分類，實驗數據被分為b={1000，3000}、b={3000，5000}這兩組.

5.1 實驗參數設置

關于算法1：α的初始值設為10-6，而梯度下降法中用于更新α的單次下降步長d=2×10-10.對于結合系數l2正則化的實驗方法，設正則化參數λ=0.00025，單次迭代次數為1000次.對于結合ADMM算法的系數l1正則化參數求解方法，設各約束項的懲罰因子為ρ1=ρ2=ρ3=1×10-4，用于判斷是否終止迭代的極小值設為ε1=ε2=ε3=1×10-8，單次迭代次數為300次.

5.2 相關評估指標

5.2.1 歸一化均方差(Normalized Mean-Squared Error，NMSE)

(27)

(28)

同樣，對于RTAP，MSD，MFD等標量的NMSE計算方式與式(28)一致，評估方法也一致：即更低的NMSE意味著更高的重建精度.

5.2.2 估算角度(Estimated angle，EA)

每個體素的纖維主方向主要是由該處的ODF的峰值所對應的方向所決定的.纖維束追蹤技術往往根據ODF的峰值來還原大腦中神經纖維束的連接結構，因此ODF的精確重建對于大腦白質中神經連接性的相關研究至關重要.本實驗中所用體模的纖維交叉角度為45°，根據雙纖維交叉處的每個體素對應的ODF，計算平均EA的公式為：

(29)

其中ΩD表有雙峰值ODF的體素集合，且|ΩD|表示該體素集的體素總個數.dx，1，dx，2則為位于x處的體素，其ODF雙峰值所對應的兩個單位矢量.顯然當EA越接近實際值45°時，ODF的重建精度越高.

5.3 實驗結果分析

5.3.1 ODF與估算角度

三種方法所估算的ODF可視化結果見圖2，圖2左側一列是根據金標準數據計算所得ODF，圖2右側一列是根據在K=60，b={1000，3000}條件下(常用的臨床掃描方式)采樣數據所估計的ODF.

可見三種方法均能根據金標準數據成功地檢測出體素中的交叉纖維.然而三種方法對ODF的重建效果均隨著采樣數據的減少而逐漸變差：不難看出在交叉區域，三種方法均有未檢測到的交叉纖維.圖2(d)，圖2(f)中未顯示左上角和右下角區域是因為在這些各向同性的區域處所估算的ODF出現了負值.

圖2 ODF可視化Fig. 2 Visualization of ODF

估算角度結果見圖3，可見無論在低b值還是高b值條件下，三種方法均在K=15處取得了最大值，這是由于采樣數據過稀疏所導致的.當b={1000，3000}時，由原始方法，l1，l2正則化方法所求EA與實際值的平均角度誤差分別是：1.29°，0.46°與1.05°.在高b值處，三種方法所得EA與實際值的平均角度誤差分別是：0.43°，0.39°與1.14°.因此結合兩種結果來看，l1正則化方法在估算纖維交叉角度問題上更勝一籌.而在高b值處由于信號整體信噪比降低，l1正則化方法使得整體EA的平均值小于45°，但是整體誤差相較低b值時有所降低.

圖3 (a)b={1000，3000}(b)b={3000，5000}時的估算角度Fig.3 Estimated Angle on b-value shells with (a)b={1000，3000}(b)b={3000，5000}

5.3.2 各標量的NMSE

重建信號對應的NMSE見圖4，在b={1000，3000}時，l1，l2正則化方法對應的NMSE相較原方法分別下降了76.4%和64.9%，而在b={3000，5000}時，二者相較原方法NMSE分別下降了65.8%和68.1%，可見改進方法大大提升了信號的重建精度.

圖4 重建信號的NMSEFig.4 NMSE of Signal on b-value shells

b={1000，3000}條件下各標量評估結果見圖5，可見兩種改進方法重建的MFD對應的NMSE遠低于原方法，l2正則化方法對應的NMSE相較原方法降低了32%，效果略優于l1正則化方法(二者相差近10.2%).且l2正則化方法取得了最佳的MSD的重建效果，而l1正則化方法對MSD的重建效果卻不盡人意，隨著K的增大，NMSE顯著提高，在K=60處取得了最大值，相較原方法提升了44%.但MSD的整體NMSE保持在0.0012以下，依舊保留了較高的重建精度.

圖5 b={1000, 3000}時(a)MFD (b)MSD (c)RTAP (d)RTOP的NMSEFig.5 NMSE of (a)MFD (b)MSD (c)RTAP (d)RTOPon b-value shells with b={1000,3000}

就重建RTAP、RTOP而言，改進方法在K=15處NMSE取得了峰值，甚至RTAP對應的NMSE最大值超過了原方法，這是由于梯度方向過少，采樣數據過稀疏所導致的高誤差.當K>20，改進方法所重建RTAP對應的NMSE明顯降低，l1，l2正則化方法下的NMSE相較原方法分別下降了69.6%及60%.不過對于RTOP的重建，即使在K=15處，基于l1，l2正則化方法的NMSE相較原方法也分別降低了38.5%與44.3%，整體NMSE遠遠低于原方法(二者分別降低了88.2%與86.4%).

b={3000，5000}條件下各標量評估結果見圖6，由圖6(c)可知改進方法在K≤30時RTAP的重建效果不如原始方法，可見改進方法在稀疏采樣的條件下對RTAP的重建并不理想.但當K>30時，基于l1，l2正則化的改進方法對應的NMSE相較原方法分別下降了29.7%與42.7%.此外改進方法下其他指標對應NMSE的數據趨勢與低b值條件下類似，均達到了相對理想的重建效果，且l1正則化方法對MSD的重建效果相較低b值時顯著提升，l1，l2約束方法對應的NMSE相較原方法平均分別下降了84.5%和79%，重建效果顯著提升.但在高b值條件下標量的整體NMSE高于低b值同條件下的NMSE，這是由于隨著b值增大，信噪比隨之降低，噪聲對信號的影響程度變大所導致的.

圖6 b={3000, 5000}時(a)MFD (b)MSD (c)RTAP (d)RTOP的NMSEFig.6 NMSE of (a)MFD (b)MSD (c)RTAP (d)RTOP on b-value shells with b={3000,5000}

本文所有實驗數據的運行條件如下：中央處理器為Intel(R)Core(TM) i7-7700 CPU，運行內存與主頻分別為8GB和3.60GHz，測試平臺是Windows 8環境下的MATLAB 2018b.三種方法的計算時間對比見圖7.在b={1000，3000}，同采樣數據條件下，l1，l2正則化方法所需平均計算時間分別是原方法的6%與33.8%，大大提升了計算效率.而在b={3000，5000}的同數據計算條件下，l1，l2正則化方法所需計算時間更短，平均計算時間分別為2.72秒及8.32秒.

圖7 (a) b={1000,3000} (b) b={3000,5000}時的計算時間Fig.7 Caculating time on b-value shells with (a)b={1000,3000}(b)b={3000, 5000}

結合兩種方法對應的NMSE表現來看，在低b值條件下l2正則化方法的重建效果相較l1正則化方法更穩定，能夠有效提升各標量重建精度，因此更適用于低b值條件下的成像.而在高b值條件下，雖然兩種改進方法對各標量重建精度的提升相差無幾，但l1正則化方法對稀疏采樣數據的重建相較l2正則化方法略勝一籌，而在采樣數據充分(K≥30)的條件下，選取l1正則化方法可兼顧計算效率及各標量重建精度.

6 討論與總結

本文在用高斯核RBF作為頻域dMRI信號插值基函數的方法基礎上，引入了信號自適應衰減項，以實現在多殼條件下對每個體素進行信號衰減建模，同時改進了擴散張量的求解方式，運用基于l1，l2正則化的方法分別求解參數以進行對比.該改進算法不僅保留了原方法中以RBF做基函數的優勢，即以穩定解析的方式計算對受限擴散信號更敏感的高階統計量(如MSD、MFD等)，并保證了ODF的計算方式在任何纖維束追蹤算法下的普遍適用性，還能在提升整體重建精度的同時大大縮短所需計算時間，減少了原方法中相關約束矩陣帶來的計算負擔.

但是本文中的實驗數據只局限于體模仿真數據，沒有將實際掃描過程中的噪聲因素引入到重建過程中，且因缺少實際腦成像數據，無法分析該算法對腦成像的實際影響.同時在稀疏采樣條件下該改進方法就個別成像指標的重建效果仍有待提升.此外，由于dMRI技術的普遍缺陷：采集信號所需的長時間掃描仍是dMRI領域的一大難題，因此為減少掃描時間，提升dMRI的臨床實用性，如何在多殼條件下根據少量梯度方向下采集的稀疏數據恢復頻域信號，并進一步重建出相關成像指標是接下來的工作.