信號(hào)自適應(yīng)衰減的多殼擴(kuò)散磁共振成像方法

羅伶俐，王遠(yuǎn)軍

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

1 引 言

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

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

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

(1)

其中P(r)表示運(yùn)動(dòng)粒子凈位移為r的概率.q為q空間中的波矢量.由于擴(kuò)散信號(hào)是掃描采樣所得數(shù)據(jù)，因此EAP的求解問題簡(jiǎn)化成了如何根據(jù)所采樣的散點(diǎn)信號(hào)數(shù)據(jù)，估計(jì)出表示信號(hào)的連續(xù)函數(shù)E(q).

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

2 徑向基函數(shù)在dMRI中的應(yīng)用

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

(2)

其中s(x)作為待估計(jì)函數(shù)，可被表示成N個(gè)中心點(diǎn)在cn處，對(duì)應(yīng)權(quán)重為ωn的RBF的線性組合.相關(guān)RBF的基礎(chǔ)理論[14]表明：當(dāng)RBF的數(shù)量N足夠大，且中心點(diǎn)cn的選取適當(dāng)時(shí)，則RBF幾乎可以逼近任意連續(xù)函數(shù).

(3)

(4)

由于dMRI信號(hào)具有偶對(duì)稱性質(zhì)，即E(q)=E(-q)，因此對(duì)偶項(xiàng)的權(quán)重保持一致：

(5)

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

3 EAP及基于EAP計(jì)算的相關(guān)指標(biāo)

3.1 EAP的計(jì)算

由于將RBF作為q空間信號(hào)的插值基函數(shù)，則所估計(jì)的信號(hào)函數(shù)表達(dá)式E(q)是一系列高斯核函數(shù)的線性組合，因此對(duì)其進(jìn)行傅里葉逆變換后所得的EAP也可對(duì)應(yīng)地被寫作相關(guān)基函數(shù)的線性組合：

(6)

(7)

3.2 方向分布函數(shù)(ODF)

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

(8)

(9)

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

由于潛在纖維組織的重要特征可通過對(duì)EAP的高階矩張量計(jì)算所得，用高斯核RBF作為信號(hào)插值基函數(shù)的一大優(yōu)勢(shì)在于：可解析地計(jì)算EAP的高階矩統(tǒng)計(jì)量，這對(duì)于發(fā)掘白質(zhì)結(jié)構(gòu)微小異常的研究非常重要.因此衍生出了基于二階矩張量及四階矩張量計(jì)算的成像統(tǒng)計(jì)標(biāo)量：均方位移(MSD)，及平均四階位移(MFD).

其中MSD正比于測(cè)量時(shí)間內(nèi)的水分子平均擴(kuò)散量，MFD則對(duì)EAP的各向異性部分更敏感，因此可捕捉到更多潛在纖維組織的相關(guān)信息.二者計(jì)算方法為：

(10)

(11)

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

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

(12)

(13)

4 基于RBF的dMRI信號(hào)自適應(yīng)衰減模型

4.1 信號(hào)自適應(yīng)衰減項(xiàng)

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

(14)

(15)

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

4.2 參數(shù)估算

4.2.1 擴(kuò)散張量計(jì)算

在計(jì)算插值基函數(shù)之前，須先用DTI方法求解二階擴(kuò)散張量：

(16)

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

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

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

(17)

4.2.2 構(gòu)建基函數(shù)矩陣及約束條件

但當(dāng)所采集的信號(hào)數(shù)量較少時(shí)(例如小于基函數(shù)數(shù)量N的情況)，那么方程Av=E就成了一個(gè)欠定問題，因此需對(duì)所求解的系數(shù)向量進(jìn)行相關(guān)正則化以減小誤差，結(jié)合常用的正則化方法：l2范數(shù)正則(又稱蒂霍諾夫正則化)，目標(biāo)函數(shù)構(gòu)建為：

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

(18)

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

(19)

4.2.3 參數(shù)求解

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

算法1.求解系數(shù)向量v與衰減系數(shù)α

輸入：基函數(shù)矩陣A，信號(hào)向量E

輸出：系數(shù)向量v，衰減系數(shù)α

a) 初始化α

b) whileα>0

c) 根據(jù)目標(biāo)函數(shù)式(18)或式(19)求解系數(shù)向量v

d) 根據(jù)‖Av-E‖2對(duì)α求導(dǎo)，使用梯度下降法更新α

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

f) break;

g) end

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

(20)

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

(21)

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

(22)

(23)

乘子項(xiàng)的更新迭代公式為：

(24)

(25)

(26)

5 實(shí)驗(yàn)方法與結(jié)果

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

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

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

5.1 實(shí)驗(yàn)參數(shù)設(shè)置

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

5.2 相關(guān)評(píng)估指標(biāo)

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

(27)

(28)

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

5.2.2 估算角度(Estimated angle，EA)

每個(gè)體素的纖維主方向主要是由該處的ODF的峰值所對(duì)應(yīng)的方向所決定的.纖維束追蹤技術(shù)往往根據(jù)ODF的峰值來還原大腦中神經(jīng)纖維束的連接結(jié)構(gòu)，因此ODF的精確重建對(duì)于大腦白質(zhì)中神經(jīng)連接性的相關(guān)研究至關(guān)重要.本實(shí)驗(yàn)中所用體模的纖維交叉角度為45°，根據(jù)雙纖維交叉處的每個(gè)體素對(duì)應(yīng)的ODF，計(jì)算平均EA的公式為：

(29)

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

5.3 實(shí)驗(yàn)結(jié)果分析

5.3.1 ODF與估算角度

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

可見三種方法均能根據(jù)金標(biāo)準(zhǔn)數(shù)據(jù)成功地檢測(cè)出體素中的交叉纖維.然而三種方法對(duì)ODF的重建效果均隨著采樣數(shù)據(jù)的減少而逐漸變差：不難看出在交叉區(qū)域，三種方法均有未檢測(cè)到的交叉纖維.圖2(d)，圖2(f)中未顯示左上角和右下角區(qū)域是因?yàn)樵谶@些各向同性的區(qū)域處所估算的ODF出現(xiàn)了負(fù)值.

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

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

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

5.3.2 各標(biāo)量的NMSE

重建信號(hào)對(duì)應(yīng)的NMSE見圖4，在b={1000，3000}時(shí)，l1，l2正則化方法對(duì)應(yīng)的NMSE相較原方法分別下降了76.4%和64.9%，而在b={3000，5000}時(shí)，二者相較原方法NMSE分別下降了65.8%和68.1%，可見改進(jìn)方法大大提升了信號(hào)的重建精度.

圖4 重建信號(hào)的NMSEFig.4 NMSE of Signal on b-value shells

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

圖5 b={1000, 3000}時(shí)(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而言，改進(jìn)方法在K=15處NMSE取得了峰值，甚至RTAP對(duì)應(yīng)的NMSE最大值超過了原方法，這是由于梯度方向過少，采樣數(shù)據(jù)過稀疏所導(dǎo)致的高誤差.當(dāng)K>20，改進(jìn)方法所重建RTAP對(duì)應(yīng)的NMSE明顯降低，l1，l2正則化方法下的NMSE相較原方法分別下降了69.6%及60%.不過對(duì)于RTOP的重建，即使在K=15處，基于l1，l2正則化方法的NMSE相較原方法也分別降低了38.5%與44.3%，整體NMSE遠(yuǎn)遠(yuǎn)低于原方法(二者分別降低了88.2%與86.4%).

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

圖6 b={3000, 5000}時(shí)(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}

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

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

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

6 討論與總結(jié)

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

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