約束球面反卷積可以從大腦擴散磁共振成像數據中量化白質纖維取向分布信息，該方法僅適用于單殼的擴散磁共振成像數據；在包含各向同性擴散信號的白質組織中，該方法會提供錯誤的纖維方向信息。針對這一不足，本文在約束球面反卷積的基礎上，結合多殼數據和多種擴散模型下估計的響應函數，提出一種基于多模型響應函數的約束球面反卷積方法。多殼數據可以提高纖維方向估計的穩定性，多模型響應函數可以衰減腦白質中各向同性擴散信號，提供更加準確的纖維方向信息。為了驗證算法的有效性，利用模擬數據和來自公開數據庫的真實人腦數據進行對比實驗。結果表明，本文算法可以衰減白質組織中的各向同性擴散信號，克服部分容積效應的影響，能更準確地進行纖維方向估計；重建的纖維方向分布穩定，偽峰少，而且對交叉纖維的識別能力也更強，為纖維束追蹤技術的進一步研究奠定了基礎。
引用本文： 潘映鈺, 王遠軍. 基于多模型響應函數約束球面反卷積的纖維方向估計. 生物醫學工程學雜志, 2022, 39(6): 1117-1126. doi: 10.7507/1001-5515.202202034 復制
引言
擴散磁共振成像（difussion magnetic resonance imaging，dMRI）是研究腦白質（white matter，WM）微結構的成像技術。其中，擴散張量成像（difussion tensor imaging，DTI）在臨床上有較高的使用率，該方法使用擴散張量的主特征向量來描述WM纖維的方向，但不適用于多根纖維共存特別是存在部分容積效應的腦區[1]。在磁共振成像（magnetic resonance imaging，MRI）掃描過程中，因受到層厚內不同組織的影響，掃描結果出現偏差，MRI設備測算出的數值不能代表腦組織本身的MRI值，這種現象被稱為部分容積效應[2]。大腦中的部分容積效應主要由灰質（grey matter，GM）、腦脊液（cerebrospinal fluid，CSF）等非WM組織引起的，如果忽略這部分差異，會影響神經纖維方向估計的準確性[3]。為了解決上述的局限性，研究者開發了許多新的方法，如擴散峰度成像（difussion kurtosis imaging，DKI）[4]、q空間成像（q space imaging，QBI）[5]、約束球面反卷積（constrained spherical deconvolution，CSD）[6]等。
CSD方法通過單位球面上多個擴散梯度方向獲取的數據來估計WM纖維走向，該方法是一種基于單纖維響應函數和纖維方向分布函數的神經纖維重建技術[7-8]。CSD方法忽略了dMRI信號中的各向同性擴散信號，導致無法精確地重建交叉纖維束。Roine等[9]基于CSD方法研究了部分容積效應對WM纖維方向估計的影響。實驗結果表明，隨著體素中GM體積分數的增加，即體素中各向同性擴散信號占比增大，部分容積效應變得尤為嚴重，導致重建的交叉纖維出現嚴重偽峰，無法準確描述真實纖維的走向。Jeurissen等[10-11]在CSD的基礎上，提出了多殼（即多個b值）CSD（multi-shell CSD，MSCSD）方法，分別對WM、GM和CSF三種組織進行建模，旨在提高存在部分容積效應的WM組織纖維方向估計的準確性。
本文針對部分容積效應現象，提出一種基于多模型響應函數的CSD方法，旨在通過衰減WM組織中的各向同性擴散信號來提高WM神經纖維的重建精度。分別選取不同的擴散模型，如基于固定張量（fixed tensor，FT）的DTI模型、DKI模型、神經突方向離散度和密度成像（neurite orientation dispersion and density imaging，NODDI）模型來構建WM響應函數。通過模擬數據實驗驗證本文算法具有良好的抗噪性能，以期可有效區分WM、GM以及CSF三種組織，減少各向同性擴散信號對WM纖維重構的影響。進一步，本文通過來自公開數據庫的真實人腦數據進行實驗，驗證本文算法對神經纖維方向估計準確性是否提升并減少偽峰，以達到優化神經纖維的成像效果，給大腦纖維束追蹤提供準確的方向信息。
1 方法原理
1.1 CSD模型
球面反卷積（spherical deconvolution，SD）在研究大腦復雜結構時通常是以其中的一個體素為單位，目標體素的整體擴散信號由一個軸對稱的響應函數矩陣R（g，r）與單位球面上的方向概率密度函數F（r）來建模得出[12-13]，表達式如式(1)所示：
 |
其中，S（g）為擴散加權信號，S0為非擴散加權信號，R（g，r）為響應函數矩陣，F（r）為單位球面上的方向概率密度函數，R2為積分球域。表示纖維方向的權重系數求解[14]如式(2)、式(3)所示：
 |
 |
其中，y為S（gi）/S0，f為纖維方向的權重系數，Aij為響應函數矩陣Ri(g，r)與方向概率密度函數Fj(r)的卷積結果，A為多個擴散梯度方向下計算的響應函數矩陣。
磁共振設備采集數據時，使用的梯度方向數量難以滿足式（3）的最優求解運算，導致結果出現較大的偏差。為了克服這一缺點，Tournier等[15]提出了CSD方法，實現系數求解過程中只涉及一個非零元素的解集合，即在式（3）的基礎上添加由懲罰參數λ和平滑解集合的約束矩陣T的約束項，如式(4)所示：
 |
1.2 多模型響應函數估計
單殼CSD方法的dMRI數據通常是在b值約為3 000，60個梯度方向的條件下采集的[16]。雖然理論上采用高b值的dMRI數據重建神經纖維角度分辨率更高，但信噪比（signal to noise ration，SNR）降低也會導致結果出現偏差，因此在實際應用中要權衡SNR與角度分辨率之間的關系以使結果達到最優。隨著MRI設備的改進，數據采集時間明顯縮短，同時使得多b值dMRI數據的采集成為可能[17-18]，多b值的dMRI數據可以進一步提高纖維方向估計的穩定性。基于傳統的CSD方法，為了適應多b值的dMRI數據，本文在原來響應函數的基礎上增加行數來擴展矩陣的大小，如式(5)所示：
 |
其中，n代表不同b值的個數，擴展的響應函數矩陣RN實際上是由每個在不同b值下估計的響應函數矩陣疊加而成。為了衰減WM組織中各向同性擴散信號，響應函數矩陣需要同時結合各向異性擴散和各向同性擴散兩部分，即在式（5）的基礎上增加列數來描述不同的組織類型。根據大腦實際的組織結構，本文考慮了GM、CSF兩種各向同性擴散組織[19]，根據WM、GM以CSF的擴散特性不同，響應函數矩陣的改進RM如式（6）所示：
 |
利用表觀擴散系數（apparent difussion coefficient，ADC）方程（以符號ADC表示），如式(7)所示，可以計算各向同性擴散組織（GM、CSF）的響應函數。
 |
其中，S（gi）為第i個擴散梯度方向gi的擴散加權信號，S0為非擴散加權信號，b為擴散敏感因子。根據兩種組織的擴散特性，本文將GM擴散張量的3個特征值設為［0.000 7，0.000 7，0.000 7］ mm2/s，CSF設為［0.003，0.003，0.003］ mm2/s。對于主要呈各向異性擴散的WM組織，下文討論了構造WM響應函數矩陣的一些擴散模型。
1.2.1 基于FT的DTI模型
傳統的CSD方法計算響應函數矩陣需要先選出各向異性分數（anisotropy fraction，FA）值大于預先設定閾值的所有體素，然后通過選出的體素進行計算[20]。本文使用基于模型的方法來構建WM響應函數矩陣，根據WM組織的擴散特性將張量的3個特征值設為［0.0017，0.000 2，0.000 2］ mm2/s，如式(8)所示。這種通過設定固定特征值來構建響應函數矩陣的方法是將人腦中水分子的擴散近似為呈高斯分布，即非高斯擴散效應對于每個纖維群都可以忽略不計。
 |
其中，S（gi）為第i個擴散梯度方向gi的擴散加權信號，S0為非擴散加權信號，D代表各向異性擴散張量，DTI模型是描述大腦神經纖維束最簡單的擴散模型。
1.2.2 DKI模型
大腦中水分子的擴散由于受到細胞膜和生物大分子的影響而呈非高斯分布，且高b值條件下采集的dMRI數據呈非高斯擴散程度越明顯。DKI是基于DTI的改進方法，通過研究大腦水分子的非高斯擴散程度來分析人腦微結構的變化。Jensen等[21]的研究引入DKI模型來描述高b值dMRI數據的水分子運動受限情況，如式(9)所示：
 |
其中，SDWI、S0分別為擴散加權信號和非擴散加權信號，D為特定擴散梯度方向的擴散系數，K為峰度系數。本文先從dMRI數據中生成一個FA值大于0.7的體素模板，從這些體素中估計出張量值和峰度系數值，然后取平均值用于WM響應函數矩陣的計算。
1.2.3 NODDI模型
WM組織響應函數矩陣還可以通過微觀結構模型來構建。本文選取了NODDI模型，該模型綜合分析了神經突的密度以及纖維方向的離散度[22-23]。考慮到水分子在細胞內外的擴散程度不同，該模型將細胞內、細胞外以及CSF三種類型的微觀結構環境區分開來，每一種環境都以獨特的方式影響著組織中水分子的擴散[24]。本文根據大腦組織的擴散特性，將細胞內的體積分數設為1，自由擴散系數設為0.001 7 mm2/s，濃度參數設為3.5，各向同性體積分數設為0來構造響應函數矩陣，這些參數都是描述WM擴散特性的典型參數。
1.3 改進的CSD模型
在WM組織和非WM組織的交界處，WM信號多少都會受到如GM、CSF等各向同性擴散信號的影響[25]，這就是上文提到的部分容積效應。由于WM組織的各向異性擴散程度遠大于各向同性擴散，通常認為其總體擴散是呈各向異性的，因此研究人員通常忽略WM組織中存在著各向同性擴散信號這部分差異，沒有將總體呈各向同性擴散的GM、CSF組織考慮在SD模型當中，這給存在明顯部分容積效應區域的纖維方向估計帶來許多不確定性。本文在傳統CSD方法的基礎上，嘗試將各向異性和各向同性兩種不同的擴散加權信號分離以衰減WM組織中的各向同性信號，減少部分容積效應對WM纖維重構的影響，以期產生更好的纖維束成像效果。
基于傳統的CSD方法，改進的CSD模型表達式如式(10)所示：
 |
其中，S（g）為擴散加權信號，S0為非擴散加權信號，
為各向同性擴散信號，RM，WM（g，r）為響應函數矩陣RM的WM響應函數列，F（r）為單位球面上的方向概率密度函數，R2為積分球域。
改進的CSD模型求解表達式如式(11)～(13)所示：
 |
 |
 |
其中，用 
表示 
/S0項，為WM中常被忽略的各向同性擴散信號，S（g）/S0為歸一化的dMRI信號，用字母y表示；f為纖維方向的權重系數，
為第j次迭代計算的各向同性擴散信號，f j為第j次迭代計算的纖維方向的權重系數，T為一個單位矩陣，λ為懲罰參數；RM，WM為響應函數矩陣RM的WM響應函數列，AM，WM是經過一定映射變換的多b值WM響應函數，Y是一個大小為n × 2的響應函數矩陣，描述GM、CSF兩種各向同性擴散組織，
表示球形卷積運算，
是初始f的進一步約束，即將小于f中值的所有系數設為0，旨在進一步減少神經纖維重建時偽峰的數量；fWM、fGM 、fCSF分別代表WM、GM、CSF的體積分數，
、
、
分別代表第j次迭代計算的WM、GM、CSF體積分數、
為第j次迭代計算的WM纖維方向的權重系數，j為迭代次數，本文設置最大迭代次數為50，但當相鄰兩次迭代所得的體積分數解不再變化時，提前結束迭代過程，如式(14)所示：
 |
其中，ε為一個閾值參數，本文設為0.001，在第一次迭代時，fGM 、fCSF的值都設為0。
2 實驗結果與分析
為了驗證部分容積效應下，改進的CSD方法的有效性，本文以算法開發及數值計算軟件Matlab R2016a（MathWorks Inc.，美國）為平臺，首先利用模擬的dMRI數據進行實驗，測試算法對噪聲的魯棒性以及分離各向同性信號的有效性，通過角度誤差以及產生偽峰的數量進行精確度評價；然后利用真實的大腦數據進行實驗，選取不同的大腦區域，將傳統的CSD方法與本文算法重建的纖維束對比分析。
2.1 模擬數據實驗
模擬的dMRI信號可以由單體素隔室的信號和來表示[26]。通常將擴散方向與程度大小具有相同性質的區域歸為一個“隔室”。在真實的人腦數據中，體素中除了WM隔室外還包括了GM、CSF等。Dell’acqua等[1]不是單純地用二階擴散張量來描述一個隔室內的纖維束，而是疊加多個隔室的貢獻。假設在沒有噪聲的情況下，疊加信號的表達式如式(15)所示：
 |
其中，S（gi）為第i個擴散梯度方向gi的擴散加權信號，S0為非擴散加權信號，Di為第i個擴散梯度方向對應的各向異性擴散張量，DGM為GM組織的擴散張量，DCSF為CSF組織的擴散張量，fj代表第j根神經纖維所占的體積分數，N代表WM中纖維束的總數，fGM 、fCSF分別為GM、CSF的體積分數，體素內三種組織的體積分數之和為1，如式(16)所示：
 |
在模擬數據實驗中，為了方便計算通常將的S0設為1。本文使用的模擬數據是在b = {0，1 000，3 000}，對應的擴散梯度向量數為g = {1，70，70}的條件下生成的。基于WM、GM、CSF三隔室模型模擬夾角為60°的交叉纖維，將WM擴散張量的特征值設置為[0.001 7，0.000 2，0.000 2] mm2/s，而各向同性組織GM和CSF平均擴散系數值分別設置為0.001 7 mm2/s和0.003 mm2/s。三種組織的信號體積分數滿足fWM + fGM + fCSF = 1。實驗進行了以下3種假設：
模擬數據1：fCSF = 0，fWM以0.1的步長從0增大到1；
模擬數據2：fGM = 0，fWM以0.1的步長從0增大到1；
模擬數據3：fWM以0.1的步長從0增大到1，fGM = fCSF = 0.5 × （1 ? fWM）。
上述模擬數據均是在假設無噪聲的條件下生成的，為了驗證改進的CSD算法的抗噪性能，實驗過程中向模擬數據加入了不同強度的萊斯噪聲，最終得到SNR分別為{80，100，1 000}的模擬數據。
2.2 實驗結果
模擬數據實驗研究了本文算法在不同SNR水平下分離WM、GM、CSF三種組織的性能。三組模擬數據估計的信號分數圖如圖1所示，每一行依次為WM體積分數fWM、GM體積分數fGM、CSF體積分數fCSF，每一列依次為模擬數據1、模擬數據2、模擬數據3的實驗結果。由此可以看出本文提出的方法可以從混合信號中準確地分離出擴散程度不同的信號，且隨著SNR增加，估計的組織信號分數值越接近真實值，說明本文算法對噪聲具有較強的魯棒性。模擬數據1和模擬數據2中的混合信號比較容易分離，而模擬數據3（即同時包含GM、CSF兩種各向同性擴散信號）估計的fWM偏高，fGM偏低，可能是由于GM組織中少量水分子的各向異性擴散。但估計值與真實值的線性變化呈相同的趨勢，并且可以較準確地估計出模擬的交叉纖維角度，因此認為可以將本文算法應用到真實人腦數據實驗中。
圖1
改進的CSD方法估計的三種組織信號分數圖
Figure1.
Three tissue signal fraction maps estimated by improved CSD method
為了展示纖維方向估計的準確性，本文通過計算預測的纖維交叉角度和估計的偽峰值數量來進行精確度的評價。兩根纖維方向的交叉角度計算公式如式(17)所示[27]：
 |
其中，v1為第一根纖維的方向向量，v2為另一根纖維的方向向量。計算結果如圖2所示，子圖依次為預測角度圖、估計的偽峰數量圖。由此可以看出，雖然在重建交叉纖維時出現了少量的偽峰，但三組模擬數據均能較準確地估計出主纖維的方向。隨著SNR增大，預測的纖維交叉角度更接近真實值60°，出現偽峰的數量減少，說明本文算法對噪聲具有較強的魯棒性。同一SNR水平下，模擬數據3較前兩組數據估計的交叉纖維不僅角度誤差更小，而且檢測到偽峰的數量都不高于前兩組數據，說明本文算法適用于同時包含WM、GM和CSF三種擴散程度不同的實際人腦數據。
圖2
模擬數據的預測角度與偽峰數量估計
Figure2.
Prediction angle and false peaks number estimation for synthetic data
2.3 真實大腦數據實驗
本實驗使用的真實大腦數據來自爬蟲擴散圖像（diffusion image in python，DIPY）數據庫公開的功能整合神經科學中心的醫學圖像（center of functionally integrative neuroscience，CFIN）數據集（網址為：https://github.com/dipy），該數據是在多個擴散敏感梯度強度下采集的，參數如下：體素量為96 × 96 × 19，體素大小為2.5 × 2.5 × 2.5，選取b值為{0，1 000，3 000}的圖像數據，對應的擴散梯度向量數為{1，33，33}。
實驗將三種不同擴散模型下構造的WM響應函數矩陣分別運用到改進的CSD方法中，分別命名為FT-CSD、DKI-CSD、NODDI-CSD，進行組織信號分數及纖維方向分布估計，并對比分析了本文算法與MSCSD方法之間的纖維重建精度。四種方法估計的大腦組織信號分數圖如圖3所示，每一行依次為估計的WM信號分數圖fWM、GM信號分數圖fGM、CSF信號分數圖fCSF；每一列依次為MSCSD、FT-CSD、DKI-CSD、NODDI-CSD四種方法的估計結果；由此可以看出本文算法的估計結果基本符合真實的人腦解剖結構；且構造響應函數時擴散模型的選擇也會影響組織信號分數的估計，FT-CSD方法估計的GM的信號分數在WM區域的取值為零，與傳統的MSCSD方法相似；而DKI-CSD、NODDI-CSD兩種方法在該區域產生非零值，說明WM組織中除了各向異性擴散信號外還包含了各向同性擴散信號，即腦WM組織存在部分容積效應。如果忽略這部分差異，會影響神經纖維方向估計的準確性，導致較差的纖維成像效果。
圖3
大腦組織信號分數圖
Figure3.
Brain tissue signal fraction map
腦部神經纖維束的重建結果如圖4所示，每幅子圖依次為MSCSD、FT-CSD、DKI-CSD、NODDI-CSD四種方法的重建結果；綠色矩形框代表腦部典型的單纖維區域，由此可以看出，FT-CSD、DKI-CSD、NODDI-CSD三種方法重建的纖維與MSCSD方法走向一致，方向估計準確，符合真實的大腦纖維走向，纖維成像效果好。黃色矩形框代表腦部交叉纖維區域，由此可以看出，MSCSD方法對交叉纖維的識別能力比較差，容易產生偽峰；而FT-CSD、DKI-CSD、NODDI-CSD三種方法在減少偽峰出現的同時也能保持較好的交叉纖維角度分辨率，證明本文算法可以提高WM纖維的重建精度。此外，DKI-CSD和NODDI-CSD方法較FT-CSD方法，重建的纖維方向更準確，偽峰數量較少。這可能是由于FT模型是基于水分子擴散呈高斯分布的假設，而DKI研究的是大腦水分子的非高斯擴散程度，后者符合大腦真實的解剖結構特性，表明構造WM響應函數時擴散模型的選擇會影響纖維方向估計的準確性。
圖4
單纖維和交叉纖維區域的纖維重建結果
Figure4.
Fiber reconstruction results in single fiber and cross fiber region
為了進一步驗證本文算法可以減少部分容積效應對WM纖維重建的影響，實驗選取了存在明顯部分容積效應的區域（WM和GM的交界處）進行纖維束重建，重建結果如圖5所示，每幅子圖依次為MSCSD、FT-CSD、DKI-CSD、NODDI-CSD四種方法的重建結果；黃色矩形框代表腦部存在明顯部分容積效應的兩個區域，由此可以看出，在模型求解過程中衰減WM組織中的各向同性擴散信號，重建的神經纖維穩定，方向估計更加準確，在一定程度上減少了偽峰的數量。DKI-CSD、NODDI-CSD兩種方法較FT-CSD方法在重建的纖維視覺上更清晰平滑，偽峰數量少。說明本文算法可以解決由于部分容積效應造成纖維估計誤差大的問題，使得估計的纖維走向更符合真實的大腦纖維情況，具有一定的可行性。
圖5
WM和GM交界處的纖維重建結果
Figure5.
Fiber reconstruction results at the interface of WM and GM
為了比較本文算法與MSCSD在重建全腦WM纖維數量的差異，分別計算了四種方法在整個大腦WM中檢測到的峰值數量，如圖6所示，每幅子圖依次為MSCSD、FT-CSD、DKI-CSD、NODDI-CSD四種方法的計算結果，縱軸表示其中某類纖維數量占總纖維數量的比率，用歸一化頻率（normalized frequency，NF）表示。考慮到大腦中80%的體素交叉纖維的數量小于4[28]，因此本文將數量閾值設為3。由此可以看出，FT-CSD與MSCSD估計的WM纖維方向數量相近，而NODDI-CSD估計的具有3個或更多交叉纖維的體素數量減少；這可能是由于NODDI模型考慮了軸突擴散的影響[23]，從而降低了交叉纖維的角度分辨率。DKI-CSD方法在減少偽峰數量的同時，保留了原始CSD方法有效分辨多交叉纖維的能力。
圖6
檢測到的WM神經纖維數量與總數的占比
Figure6.
The proportion of the number and total number of nerve fibers detected
3 結論
本文提出一種基于多模型響應函數的CSD方法，并分別對算法抗噪性能和區分混合信號的有效性進行了分析，模擬數據實驗結果表明，本文算法對噪聲具有較強的魯棒性，且可以區分擴散程度不同的混合信號，實現各向同性擴散信號與各向異性擴散信號的有效分離。進一步，通過來自公開數據庫的真實人腦數據進行實驗，結果表明，本文算法可以減小部分容積效應對WM纖維重建效果的影響，實現纖維方向的準確估計，使得重建的纖維視覺上更平滑穩定，偽峰較少，且可以以較高的角度分辨率識別交叉纖維。因此，本文所提算法能為纖維束追蹤的研究提供更準確可靠的方向信息。但NODDI模型由于考慮了軸突的離散度影響，導致交叉纖維的角度分辨率降低，之后可以進一步優化該模型的初始化參數，在提高纖維方向估計準確性的同時提高交叉纖維的識別能力，獲得最優的纖維成像效果，這將對纖維束跟蹤技術在腦認知研究和神經退行性疾病診斷中的實踐應用非常有利。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：潘映鈺是本研究的實驗設計者和實驗研究的執行人，完成實驗數據的收集，論文初稿的寫作；王遠軍是該研究的負責人，指導實驗過程、論文寫作與修改。二位作者都閱讀并同意最終的文本。
引言
擴散磁共振成像（difussion magnetic resonance imaging，dMRI）是研究腦白質（white matter，WM）微結構的成像技術。其中，擴散張量成像（difussion tensor imaging，DTI）在臨床上有較高的使用率，該方法使用擴散張量的主特征向量來描述WM纖維的方向，但不適用于多根纖維共存特別是存在部分容積效應的腦區[1]。在磁共振成像（magnetic resonance imaging，MRI）掃描過程中，因受到層厚內不同組織的影響，掃描結果出現偏差，MRI設備測算出的數值不能代表腦組織本身的MRI值，這種現象被稱為部分容積效應[2]。大腦中的部分容積效應主要由灰質（grey matter，GM）、腦脊液（cerebrospinal fluid，CSF）等非WM組織引起的，如果忽略這部分差異，會影響神經纖維方向估計的準確性[3]。為了解決上述的局限性，研究者開發了許多新的方法，如擴散峰度成像（difussion kurtosis imaging，DKI）[4]、q空間成像（q space imaging，QBI）[5]、約束球面反卷積（constrained spherical deconvolution，CSD）[6]等。
CSD方法通過單位球面上多個擴散梯度方向獲取的數據來估計WM纖維走向，該方法是一種基于單纖維響應函數和纖維方向分布函數的神經纖維重建技術[7-8]。CSD方法忽略了dMRI信號中的各向同性擴散信號，導致無法精確地重建交叉纖維束。Roine等[9]基于CSD方法研究了部分容積效應對WM纖維方向估計的影響。實驗結果表明，隨著體素中GM體積分數的增加，即體素中各向同性擴散信號占比增大，部分容積效應變得尤為嚴重，導致重建的交叉纖維出現嚴重偽峰，無法準確描述真實纖維的走向。Jeurissen等[10-11]在CSD的基礎上，提出了多殼（即多個b值）CSD（multi-shell CSD，MSCSD）方法，分別對WM、GM和CSF三種組織進行建模，旨在提高存在部分容積效應的WM組織纖維方向估計的準確性。
本文針對部分容積效應現象，提出一種基于多模型響應函數的CSD方法，旨在通過衰減WM組織中的各向同性擴散信號來提高WM神經纖維的重建精度。分別選取不同的擴散模型，如基于固定張量（fixed tensor，FT）的DTI模型、DKI模型、神經突方向離散度和密度成像（neurite orientation dispersion and density imaging，NODDI）模型來構建WM響應函數。通過模擬數據實驗驗證本文算法具有良好的抗噪性能，以期可有效區分WM、GM以及CSF三種組織，減少各向同性擴散信號對WM纖維重構的影響。進一步，本文通過來自公開數據庫的真實人腦數據進行實驗，驗證本文算法對神經纖維方向估計準確性是否提升并減少偽峰，以達到優化神經纖維的成像效果，給大腦纖維束追蹤提供準確的方向信息。
1 方法原理
1.1 CSD模型
球面反卷積（spherical deconvolution，SD）在研究大腦復雜結構時通常是以其中的一個體素為單位，目標體素的整體擴散信號由一個軸對稱的響應函數矩陣R（g，r）與單位球面上的方向概率密度函數F（r）來建模得出[12-13]，表達式如式(1)所示：
|
其中，S（g）為擴散加權信號，S0為非擴散加權信號，R（g，r）為響應函數矩陣，F（r）為單位球面上的方向概率密度函數，R2為積分球域。表示纖維方向的權重系數求解[14]如式(2)、式(3)所示：
|
|
其中，y為S（gi）/S0，f為纖維方向的權重系數，Aij為響應函數矩陣Ri(g，r)與方向概率密度函數Fj(r)的卷積結果，A為多個擴散梯度方向下計算的響應函數矩陣。
磁共振設備采集數據時，使用的梯度方向數量難以滿足式（3）的最優求解運算，導致結果出現較大的偏差。為了克服這一缺點，Tournier等[15]提出了CSD方法，實現系數求解過程中只涉及一個非零元素的解集合，即在式（3）的基礎上添加由懲罰參數λ和平滑解集合的約束矩陣T的約束項，如式(4)所示：
|
1.2 多模型響應函數估計
單殼CSD方法的dMRI數據通常是在b值約為3 000，60個梯度方向的條件下采集的[16]。雖然理論上采用高b值的dMRI數據重建神經纖維角度分辨率更高，但信噪比（signal to noise ration，SNR）降低也會導致結果出現偏差，因此在實際應用中要權衡SNR與角度分辨率之間的關系以使結果達到最優。隨著MRI設備的改進，數據采集時間明顯縮短，同時使得多b值dMRI數據的采集成為可能[17-18]，多b值的dMRI數據可以進一步提高纖維方向估計的穩定性。基于傳統的CSD方法，為了適應多b值的dMRI數據，本文在原來響應函數的基礎上增加行數來擴展矩陣的大小，如式(5)所示：
|
其中，n代表不同b值的個數，擴展的響應函數矩陣RN實際上是由每個在不同b值下估計的響應函數矩陣疊加而成。為了衰減WM組織中各向同性擴散信號，響應函數矩陣需要同時結合各向異性擴散和各向同性擴散兩部分，即在式（5）的基礎上增加列數來描述不同的組織類型。根據大腦實際的組織結構，本文考慮了GM、CSF兩種各向同性擴散組織[19]，根據WM、GM以CSF的擴散特性不同，響應函數矩陣的改進RM如式（6）所示：
|
利用表觀擴散系數（apparent difussion coefficient，ADC）方程（以符號ADC表示），如式(7)所示，可以計算各向同性擴散組織（GM、CSF）的響應函數。
|
其中，S（gi）為第i個擴散梯度方向gi的擴散加權信號，S0為非擴散加權信號，b為擴散敏感因子。根據兩種組織的擴散特性，本文將GM擴散張量的3個特征值設為［0.000 7，0.000 7，0.000 7］ mm2/s，CSF設為［0.003，0.003，0.003］ mm2/s。對于主要呈各向異性擴散的WM組織，下文討論了構造WM響應函數矩陣的一些擴散模型。
1.2.1 基于FT的DTI模型
傳統的CSD方法計算響應函數矩陣需要先選出各向異性分數（anisotropy fraction，FA）值大于預先設定閾值的所有體素，然后通過選出的體素進行計算[20]。本文使用基于模型的方法來構建WM響應函數矩陣，根據WM組織的擴散特性將張量的3個特征值設為［0.0017，0.000 2，0.000 2］ mm2/s，如式(8)所示。這種通過設定固定特征值來構建響應函數矩陣的方法是將人腦中水分子的擴散近似為呈高斯分布，即非高斯擴散效應對于每個纖維群都可以忽略不計。
|
其中，S（gi）為第i個擴散梯度方向gi的擴散加權信號，S0為非擴散加權信號，D代表各向異性擴散張量，DTI模型是描述大腦神經纖維束最簡單的擴散模型。
1.2.2 DKI模型
大腦中水分子的擴散由于受到細胞膜和生物大分子的影響而呈非高斯分布，且高b值條件下采集的dMRI數據呈非高斯擴散程度越明顯。DKI是基于DTI的改進方法，通過研究大腦水分子的非高斯擴散程度來分析人腦微結構的變化。Jensen等[21]的研究引入DKI模型來描述高b值dMRI數據的水分子運動受限情況，如式(9)所示：
|
其中，SDWI、S0分別為擴散加權信號和非擴散加權信號，D為特定擴散梯度方向的擴散系數，K為峰度系數。本文先從dMRI數據中生成一個FA值大于0.7的體素模板，從這些體素中估計出張量值和峰度系數值，然后取平均值用于WM響應函數矩陣的計算。
1.2.3 NODDI模型
WM組織響應函數矩陣還可以通過微觀結構模型來構建。本文選取了NODDI模型，該模型綜合分析了神經突的密度以及纖維方向的離散度[22-23]。考慮到水分子在細胞內外的擴散程度不同，該模型將細胞內、細胞外以及CSF三種類型的微觀結構環境區分開來，每一種環境都以獨特的方式影響著組織中水分子的擴散[24]。本文根據大腦組織的擴散特性，將細胞內的體積分數設為1，自由擴散系數設為0.001 7 mm2/s，濃度參數設為3.5，各向同性體積分數設為0來構造響應函數矩陣，這些參數都是描述WM擴散特性的典型參數。
1.3 改進的CSD模型
在WM組織和非WM組織的交界處，WM信號多少都會受到如GM、CSF等各向同性擴散信號的影響[25]，這就是上文提到的部分容積效應。由于WM組織的各向異性擴散程度遠大于各向同性擴散，通常認為其總體擴散是呈各向異性的，因此研究人員通常忽略WM組織中存在著各向同性擴散信號這部分差異，沒有將總體呈各向同性擴散的GM、CSF組織考慮在SD模型當中，這給存在明顯部分容積效應區域的纖維方向估計帶來許多不確定性。本文在傳統CSD方法的基礎上，嘗試將各向異性和各向同性兩種不同的擴散加權信號分離以衰減WM組織中的各向同性信號，減少部分容積效應對WM纖維重構的影響，以期產生更好的纖維束成像效果。
基于傳統的CSD方法，改進的CSD模型表達式如式(10)所示：
|
其中，S（g）為擴散加權信號，S0為非擴散加權信號， 為各向同性擴散信號，RM，WM（g，r）為響應函數矩陣RM的WM響應函數列，F（r）為單位球面上的方向概率密度函數，R2為積分球域。
改進的CSD模型求解表達式如式(11)～(13)所示：
|
|
|
其中，用 表示 /S0項，為WM中常被忽略的各向同性擴散信號，S（g）/S0為歸一化的dMRI信號，用字母y表示；f為纖維方向的權重系數， 為第j次迭代計算的各向同性擴散信號，f j為第j次迭代計算的纖維方向的權重系數，T為一個單位矩陣，λ為懲罰參數；RM，WM為響應函數矩陣RM的WM響應函數列，AM，WM是經過一定映射變換的多b值WM響應函數，Y是一個大小為n × 2的響應函數矩陣，描述GM、CSF兩種各向同性擴散組織， 表示球形卷積運算， 是初始f的進一步約束，即將小于f中值的所有系數設為0，旨在進一步減少神經纖維重建時偽峰的數量；fWM、fGM 、fCSF分別代表WM、GM、CSF的體積分數，、、 分別代表第j次迭代計算的WM、GM、CSF體積分數、 為第j次迭代計算的WM纖維方向的權重系數，j為迭代次數，本文設置最大迭代次數為50，但當相鄰兩次迭代所得的體積分數解不再變化時，提前結束迭代過程，如式(14)所示：
|
其中，ε為一個閾值參數，本文設為0.001，在第一次迭代時，fGM 、fCSF的值都設為0。
2 實驗結果與分析
為了驗證部分容積效應下，改進的CSD方法的有效性，本文以算法開發及數值計算軟件Matlab R2016a（MathWorks Inc.，美國）為平臺，首先利用模擬的dMRI數據進行實驗，測試算法對噪聲的魯棒性以及分離各向同性信號的有效性，通過角度誤差以及產生偽峰的數量進行精確度評價；然后利用真實的大腦數據進行實驗，選取不同的大腦區域，將傳統的CSD方法與本文算法重建的纖維束對比分析。
2.1 模擬數據實驗
模擬的dMRI信號可以由單體素隔室的信號和來表示[26]。通常將擴散方向與程度大小具有相同性質的區域歸為一個“隔室”。在真實的人腦數據中，體素中除了WM隔室外還包括了GM、CSF等。Dell’acqua等[1]不是單純地用二階擴散張量來描述一個隔室內的纖維束，而是疊加多個隔室的貢獻。假設在沒有噪聲的情況下，疊加信號的表達式如式(15)所示：
|
其中，S（gi）為第i個擴散梯度方向gi的擴散加權信號，S0為非擴散加權信號，Di為第i個擴散梯度方向對應的各向異性擴散張量，DGM為GM組織的擴散張量，DCSF為CSF組織的擴散張量，fj代表第j根神經纖維所占的體積分數，N代表WM中纖維束的總數，fGM 、fCSF分別為GM、CSF的體積分數，體素內三種組織的體積分數之和為1，如式(16)所示：
|
在模擬數據實驗中，為了方便計算通常將的S0設為1。本文使用的模擬數據是在b = {0，1 000，3 000}，對應的擴散梯度向量數為g = {1，70，70}的條件下生成的。基于WM、GM、CSF三隔室模型模擬夾角為60°的交叉纖維，將WM擴散張量的特征值設置為[0.001 7，0.000 2，0.000 2] mm2/s，而各向同性組織GM和CSF平均擴散系數值分別設置為0.001 7 mm2/s和0.003 mm2/s。三種組織的信號體積分數滿足fWM + fGM + fCSF = 1。實驗進行了以下3種假設：
模擬數據1：fCSF = 0，fWM以0.1的步長從0增大到1；
模擬數據2：fGM = 0，fWM以0.1的步長從0增大到1；
模擬數據3：fWM以0.1的步長從0增大到1，fGM = fCSF = 0.5 × （1 ? fWM）。
上述模擬數據均是在假設無噪聲的條件下生成的，為了驗證改進的CSD算法的抗噪性能，實驗過程中向模擬數據加入了不同強度的萊斯噪聲，最終得到SNR分別為{80，100，1 000}的模擬數據。
2.2 實驗結果
模擬數據實驗研究了本文算法在不同SNR水平下分離WM、GM、CSF三種組織的性能。三組模擬數據估計的信號分數圖如圖1所示，每一行依次為WM體積分數fWM、GM體積分數fGM、CSF體積分數fCSF，每一列依次為模擬數據1、模擬數據2、模擬數據3的實驗結果。由此可以看出本文提出的方法可以從混合信號中準確地分離出擴散程度不同的信號，且隨著SNR增加，估計的組織信號分數值越接近真實值，說明本文算法對噪聲具有較強的魯棒性。模擬數據1和模擬數據2中的混合信號比較容易分離，而模擬數據3（即同時包含GM、CSF兩種各向同性擴散信號）估計的fWM偏高，fGM偏低，可能是由于GM組織中少量水分子的各向異性擴散。但估計值與真實值的線性變化呈相同的趨勢，并且可以較準確地估計出模擬的交叉纖維角度，因此認為可以將本文算法應用到真實人腦數據實驗中。
圖1
改進的CSD方法估計的三種組織信號分數圖
Figure1.
Three tissue signal fraction maps estimated by improved CSD method
為了展示纖維方向估計的準確性，本文通過計算預測的纖維交叉角度和估計的偽峰值數量來進行精確度的評價。兩根纖維方向的交叉角度計算公式如式(17)所示[27]：
|
其中，v1為第一根纖維的方向向量，v2為另一根纖維的方向向量。計算結果如圖2所示，子圖依次為預測角度圖、估計的偽峰數量圖。由此可以看出，雖然在重建交叉纖維時出現了少量的偽峰，但三組模擬數據均能較準確地估計出主纖維的方向。隨著SNR增大，預測的纖維交叉角度更接近真實值60°，出現偽峰的數量減少，說明本文算法對噪聲具有較強的魯棒性。同一SNR水平下，模擬數據3較前兩組數據估計的交叉纖維不僅角度誤差更小，而且檢測到偽峰的數量都不高于前兩組數據，說明本文算法適用于同時包含WM、GM和CSF三種擴散程度不同的實際人腦數據。
圖2
模擬數據的預測角度與偽峰數量估計
Figure2.
Prediction angle and false peaks number estimation for synthetic data
2.3 真實大腦數據實驗
本實驗使用的真實大腦數據來自爬蟲擴散圖像（diffusion image in python，DIPY）數據庫公開的功能整合神經科學中心的醫學圖像（center of functionally integrative neuroscience，CFIN）數據集（網址為：https://github.com/dipy），該數據是在多個擴散敏感梯度強度下采集的，參數如下：體素量為96 × 96 × 19，體素大小為2.5 × 2.5 × 2.5，選取b值為{0，1 000，3 000}的圖像數據，對應的擴散梯度向量數為{1，33，33}。
實驗將三種不同擴散模型下構造的WM響應函數矩陣分別運用到改進的CSD方法中，分別命名為FT-CSD、DKI-CSD、NODDI-CSD，進行組織信號分數及纖維方向分布估計，并對比分析了本文算法與MSCSD方法之間的纖維重建精度。四種方法估計的大腦組織信號分數圖如圖3所示，每一行依次為估計的WM信號分數圖fWM、GM信號分數圖fGM、CSF信號分數圖fCSF；每一列依次為MSCSD、FT-CSD、DKI-CSD、NODDI-CSD四種方法的估計結果；由此可以看出本文算法的估計結果基本符合真實的人腦解剖結構；且構造響應函數時擴散模型的選擇也會影響組織信號分數的估計，FT-CSD方法估計的GM的信號分數在WM區域的取值為零，與傳統的MSCSD方法相似；而DKI-CSD、NODDI-CSD兩種方法在該區域產生非零值，說明WM組織中除了各向異性擴散信號外還包含了各向同性擴散信號，即腦WM組織存在部分容積效應。如果忽略這部分差異，會影響神經纖維方向估計的準確性，導致較差的纖維成像效果。
圖3
大腦組織信號分數圖
Figure3.
Brain tissue signal fraction map
腦部神經纖維束的重建結果如圖4所示，每幅子圖依次為MSCSD、FT-CSD、DKI-CSD、NODDI-CSD四種方法的重建結果；綠色矩形框代表腦部典型的單纖維區域，由此可以看出，FT-CSD、DKI-CSD、NODDI-CSD三種方法重建的纖維與MSCSD方法走向一致，方向估計準確，符合真實的大腦纖維走向，纖維成像效果好。黃色矩形框代表腦部交叉纖維區域，由此可以看出，MSCSD方法對交叉纖維的識別能力比較差，容易產生偽峰；而FT-CSD、DKI-CSD、NODDI-CSD三種方法在減少偽峰出現的同時也能保持較好的交叉纖維角度分辨率，證明本文算法可以提高WM纖維的重建精度。此外，DKI-CSD和NODDI-CSD方法較FT-CSD方法，重建的纖維方向更準確，偽峰數量較少。這可能是由于FT模型是基于水分子擴散呈高斯分布的假設，而DKI研究的是大腦水分子的非高斯擴散程度，后者符合大腦真實的解剖結構特性，表明構造WM響應函數時擴散模型的選擇會影響纖維方向估計的準確性。
圖4
單纖維和交叉纖維區域的纖維重建結果
Figure4.
Fiber reconstruction results in single fiber and cross fiber region
為了進一步驗證本文算法可以減少部分容積效應對WM纖維重建的影響，實驗選取了存在明顯部分容積效應的區域（WM和GM的交界處）進行纖維束重建，重建結果如圖5所示，每幅子圖依次為MSCSD、FT-CSD、DKI-CSD、NODDI-CSD四種方法的重建結果；黃色矩形框代表腦部存在明顯部分容積效應的兩個區域，由此可以看出，在模型求解過程中衰減WM組織中的各向同性擴散信號，重建的神經纖維穩定，方向估計更加準確，在一定程度上減少了偽峰的數量。DKI-CSD、NODDI-CSD兩種方法較FT-CSD方法在重建的纖維視覺上更清晰平滑，偽峰數量少。說明本文算法可以解決由于部分容積效應造成纖維估計誤差大的問題，使得估計的纖維走向更符合真實的大腦纖維情況，具有一定的可行性。
圖5
WM和GM交界處的纖維重建結果
Figure5.
Fiber reconstruction results at the interface of WM and GM
為了比較本文算法與MSCSD在重建全腦WM纖維數量的差異，分別計算了四種方法在整個大腦WM中檢測到的峰值數量，如圖6所示，每幅子圖依次為MSCSD、FT-CSD、DKI-CSD、NODDI-CSD四種方法的計算結果，縱軸表示其中某類纖維數量占總纖維數量的比率，用歸一化頻率（normalized frequency，NF）表示。考慮到大腦中80%的體素交叉纖維的數量小于4[28]，因此本文將數量閾值設為3。由此可以看出，FT-CSD與MSCSD估計的WM纖維方向數量相近，而NODDI-CSD估計的具有3個或更多交叉纖維的體素數量減少；這可能是由于NODDI模型考慮了軸突擴散的影響[23]，從而降低了交叉纖維的角度分辨率。DKI-CSD方法在減少偽峰數量的同時，保留了原始CSD方法有效分辨多交叉纖維的能力。
圖6
檢測到的WM神經纖維數量與總數的占比
Figure6.
The proportion of the number and total number of nerve fibers detected
3 結論
本文提出一種基于多模型響應函數的CSD方法，并分別對算法抗噪性能和區分混合信號的有效性進行了分析，模擬數據實驗結果表明，本文算法對噪聲具有較強的魯棒性，且可以區分擴散程度不同的混合信號，實現各向同性擴散信號與各向異性擴散信號的有效分離。進一步，通過來自公開數據庫的真實人腦數據進行實驗，結果表明，本文算法可以減小部分容積效應對WM纖維重建效果的影響，實現纖維方向的準確估計，使得重建的纖維視覺上更平滑穩定，偽峰較少，且可以以較高的角度分辨率識別交叉纖維。因此，本文所提算法能為纖維束追蹤的研究提供更準確可靠的方向信息。但NODDI模型由于考慮了軸突的離散度影響，導致交叉纖維的角度分辨率降低，之后可以進一步優化該模型的初始化參數，在提高纖維方向估計準確性的同時提高交叉纖維的識別能力，獲得最優的纖維成像效果，這將對纖維束跟蹤技術在腦認知研究和神經退行性疾病診斷中的實踐應用非常有利。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：潘映鈺是本研究的實驗設計者和實驗研究的執行人，完成實驗數據的收集，論文初稿的寫作；王遠軍是該研究的負責人，指導實驗過程、論文寫作與修改。二位作者都閱讀并同意最終的文本。
