基于稀疏編碼的彌散微循環(huán)模型參數(shù)估計神經(jīng)網(wǎng)絡

鄭天舒，顏國輝，葉初陽，吳 丹

（1.浙江大學生物醫(yī)學工程與儀器科學學院，杭州 310027；2.浙江大學醫(yī)學院附屬婦產(chǎn)科醫(yī)院，杭州 310003；3.北京理工大學集成電路與電子學院，北京 100081）

引 言

彌散磁共振成像（Diffusion magnetic resonance imaging，dMRI）是一種重要的醫(yī)學成像工具，它基于水分子在生物組織中受限彌散作用，可建立彌散信號的生物物理模型，無創(chuàng)地重建生物組織內(nèi)的微觀結(jié)構(gòu)［1］。常用的表觀彌散系數(shù)（Apparent diffusion coefficient，ADC）用一個單指數(shù)模型來計算，該單指數(shù)模型對腦卒中［2］、腫瘤等病理變化敏感［3-4］，但不能特異性地反映微觀層面病理變化。近30年來，領(lǐng)域內(nèi)提出了多種復雜的生物物理模型來表征特定的微結(jié)構(gòu)特性，如體素內(nèi)不相干運動模型（Intravoxel incoherent motion，IVIM）［5］、神經(jīng)元定向彌散和密度成像（Neurite orientation dispersionand density imaging，NODDI）［6］、用于腫瘤細胞術(shù)的血管、細胞外和限制性擴散模型（Vascular，extracellular and restricted diffusion for cytometry in tumors，VERDICT）［7］，胞體和神經(jīng)突密度成像（Soma and neurite density imaging，SANDI）［8］等。

準確估計微結(jié)構(gòu)模型參數(shù)對診斷具有重要意義。然而，大多數(shù)前沿的dMRI 模型由多個房室組成，并且這些房室在數(shù)學上復雜且高度非線性。用傳統(tǒng)的優(yōu)化方法例如最小二乘法，對這些模型進行擬合，容易產(chǎn)生估計誤差［9］。此外，從數(shù)據(jù)獲取的角度來看，先進的dMRI 模型往往需要在q空間中獲取多個b值和彌散梯度方向，這既耗時又容易受到運動偽影的影響。這對于運動的被試，例如腹部器官、胎兒和胎盤等，是一個不小的挑戰(zhàn)。

為了減少這些估計誤差和加速q空間采集，研究者們提出了許多方法。如Nedjati-Gilani等［10］首先提出了基于隨機森林的方法來對K?rger模型進行擬合；同一時間，Alexander等［11］通過隨機森林方法實現(xiàn)了在僅用單個b值的情況下，對NODDI 模型和球面平均技術(shù)（Spherical mean technique，SMT）模型的參數(shù)估計。

隨著深度學習技術(shù)的發(fā)展，深度學習的技術(shù)為dMRI 模型擬合開辟了新的途徑。Golkov 等［12］首先提出了基于q空間的深度學習方法，這是一種基于多層感知器（Multilayer perceptron，MLP）的方法，它利用q空間數(shù)據(jù)的一個子集來估計擴散峰度成像（Diffusion kurtosis imaging，DKI）參數(shù)。Barbieri 等［13］首次使用基于q空間的深度學習方法來對IVIM 模型進行參數(shù)估計。然而，傳統(tǒng)的基于q空間的深度學習方法與生物物理模型無關(guān)，因此很難被解釋。在深層神經(jīng)網(wǎng)絡中引入領(lǐng)域知識作為先驗信息被認為是提高網(wǎng)絡性能和可解釋性的有效途徑。Ye［14］首次將通過凸優(yōu)化加速微結(jié)構(gòu)成像（Accelerated microstructure imaging via convex optimization，AMICO）算法中字典構(gòu)建的過程與深度學習結(jié)合，摒棄之前預先定義字典的方法，通過深度學習方法來得到字典，實現(xiàn)了對NODDI 模型參數(shù)的準確估計。

由于之前的基于q空間的深度學習方法［12］缺乏先驗知識引導，受此基于稀疏先驗的模型構(gòu)建方法啟發(fā)［14］，結(jié)合之前基于q空間的深度學習方法［12］可以學習到降采樣的q空間信號與模型參數(shù)映射關(guān)系的思想。本文將兩者融合，提出了一種基于稀疏編碼深度學習網(wǎng)絡（Sparse coding deep neural network，SCDNN）的IVIM 模型參數(shù)估計網(wǎng)絡。該網(wǎng)絡由稀疏表示編碼器和參數(shù)映射解碼器兩部分構(gòu)成，其中稀疏表示編碼器和參數(shù)映射解碼器，均通過數(shù)據(jù)驅(qū)動的方式學習得到。首先，假設(shè)IVIM 模型所需要的b值具有稀疏性，通過稀疏表示編碼器，低維空間中稠密采集的dMRI 信號可以被轉(zhuǎn)化為高維空間的稀疏離散信號；隨后，通過參數(shù)映射解碼器，稀疏離散信號可以與IVIM 模型參數(shù)進行體素對體素的映射。將這兩部分結(jié)合，可以實現(xiàn)對IVIM 模型在q空間的降采樣，在采集時間更少的情況下獲得近似于相同質(zhì)量的圖像信息，具有較高的準確度和精確性。本文還將所擬的網(wǎng)絡和其他已有的dMRI 參數(shù)估計網(wǎng)絡在胎盤數(shù)據(jù)上進行了測試，發(fā)現(xiàn)基于稀疏編碼深度學習網(wǎng)絡的方法估計效果優(yōu)于其他IVIM 模型估計方法，并且更具備模型的解釋性。

1 理 論

本節(jié)將會簡要介紹IVIM 模型、傳統(tǒng)IVIM 模型參數(shù)的估計方法以及如何通過稀疏表示的方法來估計IVIM 模型的參數(shù)。

1.1 IVIM 模型

IVIM 模型由Le Bihan 等［5］提出，基于組織中和血液中水分子擴散率的不同建立雙指數(shù)模型，可以同時獲取組織（細胞、軸突等）中的水分子彌散信息和微循環(huán)（毛細血管、血小管）中的血流灌注信息。它已廣泛應用于檢測組織中的灌注情況，如大腦、腎臟、肝臟和胎盤等器官的灌注異常的檢測［15-17］。然而IVIM 成像需要采集多b值（≥10 個b值）的擴散信號，采集較長，容易受到運動偽影的影響［18］。體部器官隨呼吸運動偽影較大，特別地，胎盤的成像不僅受到呼吸的影響還受到胎兒運動的影響。

常見的IVIM 模型可以表示為

式中：b為施加的彌散梯度場強度，Sb為該b值下的信號值，S0為沒有擴散加權(quán)時的信號值，f為微循環(huán)的體積分數(shù)，D為組織水分子的彌散系數(shù)，D*為微循環(huán)血液中水分子的偽彌散系數(shù)。

1.2 IVIM 模型的參數(shù)估計

常見的IVIM 模型參數(shù)估計方法有兩類，分別為基于非線性最小二乘法的參數(shù)估計和基于貝葉斯的參數(shù)估計。對于胎盤而言，基于非線性最小二乘法的參數(shù)估計先通過分段方法，通過相對高b值（＞200 s/mm2）的數(shù)據(jù)，采用單指數(shù)模型估計D值以及組織中水分子組分對應的非擴散加權(quán)信號；通過相對低b值（＜200 s/mm2）的數(shù)據(jù)通過外插方法估計組織加血液對應非擴散加權(quán)信號，f值可通過進行估計；D*值可估計為D×10 得到f、D、D*的估算值，再以此為初始條件和上下限條件，通過建立非線性最小二乘擬合方法來對其中的f、D、D*來進行求解。貝葉斯估計則通過結(jié)合先驗信息，通過貝葉斯公式，優(yōu)化最大化后驗概率的分布來求f、D、D*［19］。

1.3 IVIM 模型的稀疏表示

稀疏表示可以被認為是一種表示學習方法，它已被普遍用于磁共振（Magnetic resonance imaging，MRI）的加速問題中［20-21］。對于IVIM 模型，q空間的彌散信號可以用以下字典學習框架來表示。

式中：y=Sb/S0∈RK是一系列通過S0歸一化后的信號Sb，其中K為該樣本中的彌散梯度個數(shù)；φ∈RK×N是被定義的字典，其中N為稀疏系數(shù)向量的長度，在本文中，該字典通過神經(jīng)網(wǎng)絡學習得到；x為與字典φ相關(guān)的稀疏系數(shù)向量；η為與x對應的噪聲向量。結(jié)合IVIM 模型，φ和x可被拓展為

式中：φD、φD*分別與離散化的D和D*相對應；x1-f對應于組織成分（D）的組成比例；xf為灌注成分（D*）的組成比例。其中為了滿足各組分之和為1 的限制，稀疏系數(shù)向量x被歸一化到［0，1］之間。

如果φ已經(jīng)被預先定義好，則稀疏系數(shù)向量可以被通過求解如下的帶有l(wèi)1正則化的目標函數(shù)得到。

式中：β為一個常數(shù)，用來控制稀疏系數(shù)向量x的稀疏性，β越大時，稀疏系數(shù)向量x越稀疏。l1正則化的優(yōu)化可以通過許多算法來解決。在MRI 領(lǐng)域，普遍使用迭代收縮閾值算法（Iterative shrinkage thresholding algorithm，ISTA）［22-23］。優(yōu)化該目標函數(shù)可以通過如下的迭代過程展開。

式中：I為單位矩陣；HM為非線性算子，可定義為

通過如上步驟，可以獲得稀疏表示的信號。隨后，可以通過如下的線性組合得到模型估計參數(shù)。

2 方 法

本節(jié)將從如何構(gòu)建基于IVIM模型的深度神經(jīng)網(wǎng)絡、訓練所需樣本以及訓練方法3個方面來進行闡述。

2.1 網(wǎng)絡構(gòu)建

如前所述，基于IVIM 模型的深度神經(jīng)網(wǎng)絡由稀疏表示編碼器和參數(shù)映射解碼器構(gòu)成，具體方法如圖1 所示。

圖1 基于IVIM 模型的SCDNN 結(jié)構(gòu)示意圖Fig.1 Structure of IVIM model based SCDNN

對于稀疏表示編碼器，它將該接受q空間降采樣的輸入信號，該輸入信號通過信號提取層后輸入字典層A后得到其對應的輸出信號，隨后將這一輸出信號輸入連續(xù)的迭代單元，迭代單元依照可以通過對式（6）字典構(gòu)建的方法進行設(shè)計，通過迭代過程實現(xiàn)字典的構(gòu)建，其中紅色虛線所示的即對應迭代過程。

每個迭代單元中輸入信號依次經(jīng)過閾值層和字典層B，閾值層為一非線性算子，對應于式（6）中的HM，字典層B 則對數(shù)據(jù)進行組合，之后將字典層B 輸出的信號與該迭代單元的輸入信號通過殘差連接相疊加，作為通過該迭代單元后的輸出信號并作為下一個迭代單元的輸入信號；第n個迭代單元中通過對信號的處理實現(xiàn)閾值迭代算法在IVIM 模型中的通項表達式：xn+1=HM[φHy+(I-φHφ)xn]，其中φH和φ對應于輸入的線性層和迭代單元中線性層所學習到的權(quán)重，φ為由離散化的組織水分子的彌散系數(shù)D和微循環(huán)血液中水分子的偽彌散系數(shù)D*所構(gòu)成的字典向量，I是單位矩陣，HM為非線性算子，y為輸入信號，xn和xn+1分別代表第n個迭代單元的輸入信號和輸出信號。如此進行反復迭代，得到通過所有迭代單元后的最終輸出信號，將最終輸出信號輸入歸一化層后，得到體積分數(shù)f，通過歸一化后，使信號限制在0～1，滿足實際生理狀況；再將通過歸一化層的輸出分為等長的兩部分，分別進行再次歸一化后輸入兩個不同線性組合層進行信號組合，一個線性組合層輸出由離散化的組織水分子的彌散系數(shù)D重新組合，另一個線性組合層輸出由離散化的微循環(huán)血液中水分子的偽彌散系數(shù)D*重新組合。

2.2 數(shù)據(jù)及預處理

采用29 名孕期在13～37 周的被試數(shù)據(jù)中進行了數(shù)據(jù)集的構(gòu)建。

2.2.1 序列參數(shù)

29 名磁共振掃描在通用電氣SIGNA HDXT 1.5 T 磁共振掃描儀上進行，采用擴散加權(quán)的平面回波成像序列從母體矢狀位采集胎盤圖像：回波時間（Echo time，TE）/重復時間（Repetition time，TR）=76/3 000 ms，視野（Field of view，F(xiàn)OV）=320 mm×320 mm，平面內(nèi)分辨率為1.25 mm×1.25 mm，層厚4 mm，共15 層，共采集10 個b值的數(shù)據(jù)，分別為0、10、20、50、80、100、150、200、300 和500 s/mm2。

2.2.2 配 準

對29 例被試中的每個胎盤數(shù)據(jù)進行自配準，以去除母體呼吸、胎兒運動導致的各個b值擴散加權(quán)圖像間的位移。在配準過程中采用循環(huán)的配準方式操作6～10 次之間不等。對每一位被試，將其各個b值的擴散加權(quán)圖像進行平均，得到平均模板。然后，通過剛體變換和仿射變換將各個b值的擴散加權(quán)圖像與平均模板進行比對配準。將配準之后得到的數(shù)據(jù)再次進行平均，得到第二次的平均模板，如此循環(huán)6～10 次，選取其中位移最小的圖像作為得到目標圖像的配準結(jié)果。

2.2.3 金標準獲取

將2.2.2 節(jié)中配準好的圖像，通過貝葉斯方法進行參數(shù)估計［24］，得到對應的f、D、D*。由于通過貝葉斯方法得到的f、D、D*存在因噪聲造成的誤差，因此將對應的f、D、D*代入式（1）中，得到該組信號的無偏信號，并添加噪聲。噪聲的添加方式與文獻［25］相同。

式中：Snoisy為數(shù)據(jù)集的輸入信號；S為通過式（1）得到的信號；ξ1，ξ2～(0，σ2)，σ=S0/SNR，其中與文獻［25］相同，S0=1，訓練數(shù)據(jù)SNR=30。

2.3 訓 練

將預處理好的29 例孕周在13～37 周的正常孕婦樣本按照6∶4 的比例劃分為訓練集和測試集。其中訓練集中有15 例樣本（324 016 個體素，其中與文獻［12］相同，90%用于訓練，10%用于驗證），測試中有9 例樣本（248 059 個體素）。本文選擇使用Adam 作為優(yōu)化器，學習率設(shè)置為1×10-4，總計10 個歷時，批次大小為128。在網(wǎng)絡訓練的前3 個歷時中使用了余弦預熱方法，學習率在3 個歷時后衰減。

3 實 驗

對于2.2 節(jié)中獲取的數(shù)據(jù)集，采用基于IVIM 模型的深度神經(jīng)網(wǎng)絡進行了估計，并與其他用于dMRI 參數(shù)估計的方法進行了性能對比。同時，對基于IVIM 模型的深度神經(jīng)網(wǎng)絡參數(shù)設(shè)置進行了分析。

3.1 性能分析

對基于IVIM 模型的深度神經(jīng)網(wǎng)絡進行了性能分析。首先，對基于IVIM 模型的深度神經(jīng)網(wǎng)絡的在訓練集和驗證集上的損失函數(shù)衰減曲線分別進行了觀察，相關(guān)結(jié)果如圖2 所示。

圖2 基于SCDNN 的IVIM 參數(shù)估計Fig.2 Parameter estimation of IVIM model based on SCDNN

3.1.1 收斂性分析

結(jié)合圖2，可以觀察到，第1 次迭代訓練的損失為0.33 左右，驗證集誤差為0.27 左右，說明模型剛開始訓練時，損失較大，準確率較低。隨著迭代周期增加，訓練集誤差和驗證集誤差在5 個歷時后趨于穩(wěn)定，可以認為其收斂，此時得到的網(wǎng)絡參數(shù)就是最終優(yōu)化的網(wǎng)絡參數(shù)指標。

3.1.2 結(jié)果分析

首先對于一例樣本進行分析，結(jié)果如圖3 所示。對于該例樣本，計算各參數(shù)網(wǎng)絡輸出的結(jié)果與其對應的金標準之間的決定系數(shù)（Coefficient of determination）R2。R2可定義為

式中：yi為該參數(shù)的金標準；fi為該參數(shù)網(wǎng)絡估計值；yˉ為該參數(shù)金標準的平均值。對于3 個參數(shù)f、D、D*，其對應的R2如圖3（a）所示，分別為0.998 2、0.999 1、0.998 5。3 個參數(shù)的網(wǎng)絡估計結(jié)果示意圖及殘差圖（參數(shù)金標準-參數(shù)網(wǎng)絡輸出結(jié)果）如圖3（b）所示。由圖3（b）中的殘差圖可以看出，f、D、D*均可以達到較好的估計效果。

圖3 神經(jīng)網(wǎng)絡輸出與金標準的相關(guān)性圖及神經(jīng)網(wǎng)絡輸出的參數(shù)圖和誤差圖Fig.3 Correlation between neural network output and gold standard and parameter maps and error maps of neural network output

在全部測試集上比較了5 種不同算法，包括兩種傳統(tǒng)優(yōu)化方法（最小二乘法和貝葉斯估計）以及3種基于學習的方法（q空間MLP［12］、IVIM 網(wǎng)絡（IVIM-NET）［13］和本文提出的基于SCDNN 的IVIM參數(shù)估計），結(jié)果如表1 所示。與傳統(tǒng)優(yōu)化方法最小二乘法和貝葉斯方法相比，基于學習的方法具有更高的估計精度。在基于學習的方法中，模型驅(qū)動方法的性能優(yōu)于無先驗信息的q空間的深度學習。

表1 比較5 種方法在使用減少的b 值（20、50、150、300、500 s/mm2）估計彌散微循環(huán)模型參數(shù)的均方誤差Table 1 Mean squared error comparison of five algorithms using a subset of b-values（20，50，150，300，500 s/mm2）for IVIM parameter estimation

3.2 泛化性分析

本節(jié)對基于IVIM 模型的深度神經(jīng)網(wǎng)絡進行了泛化性分析。為了評估基于IVIM 模型的深度神經(jīng)網(wǎng)絡的泛化性，在另一臺通用電氣GE 750 W 的3.0 T 磁共振掃描儀上通過相同的掃描序列參數(shù)進行了數(shù)據(jù)采集。采用擴散加權(quán)的平面回波成像序列從母體矢狀位采集胎盤圖像：TE/TR=76/3 000 ms，F(xiàn)OV=320 mm×320 mm，平面內(nèi)分辨率為1.25 mm×1.25 mm，層厚4 mm，共15 層，共采集10 個b值數(shù)據(jù)，分別為0、10、20、50、80、100、150、200、300 和500 s/mm2。與3.1.2 節(jié)相同，比較了所有的5 種算法，結(jié)果如表2 所示。從表2 中可以觀察到，在同一品牌的3.0 T 數(shù)據(jù)上進行測試，發(fā)現(xiàn)其他基于學習的方法有不同程度上的誤差升高，其中，無監(jiān)督IVIM-NET 誤差升高最多，本文方法具有最高的準確度，證明本文方法具有更好的泛化性。

表2 泛化性測試比較5 種方法在使用減少的b 值（20、50、150、300、500 s/mm2）在另一個中心采集到的數(shù)據(jù)估計彌散微循環(huán)模型參數(shù)的均方誤差Table 2 Generalizability test on data acquired from another center，based on the mean squared error of five algorithms using a subset of b-values（20，50，150，300，500 s/mm2）for IVIM parameter estimation

3.3 網(wǎng)絡參數(shù)分析

本節(jié)對SCDNN 的參數(shù)設(shè)置進行了分析。首先，對于不同b值的選取進行了實驗，之后對于所建立字典大小進行了實驗。

3.3.1b值組合分析

本節(jié)對于如何選取IVIM 模型訓練時所需要的b值數(shù)量和b值組合進行了研究。在之前的實驗中，選取了減少到5 的b值數(shù)量（20、50、150、300、500 s/mm2）來進行估計。在這里對5 個b值數(shù)量的不同組合進行了估計，同時也對不同數(shù)量的b值選取進行了估計。其中對數(shù)量為5 個，b值大小為20～500 s/mm2的5 種組合進行了估計，結(jié)果如表3 所示。

結(jié)合表3，可以發(fā)現(xiàn)，在b值數(shù)量為3 時，參數(shù)f、D、D*的估計結(jié)果與b值數(shù)量為5 或7 時相比較差。當b值數(shù)量為7 時，參數(shù)f、D、D*的估計結(jié)果與b值數(shù)量為5 時最優(yōu)結(jié)果總體誤差相接近，且并不優(yōu)于b值數(shù)量為5 時的結(jié)果。因此，在最終的b值選取上選擇5 個b值組合為20、50、150、300、500 s/mm2的降采樣方式。

表3 不同數(shù)量及不同b 值組合得到的IVIM 參數(shù)估計值與金標準的均方誤差Table 3 Mean squared errors of IVIM parameter estimation using different numbers of b-values and different combinations of b-values

3.3.2 字典大小分析

本節(jié)對字典大小為200～600 的結(jié)果進行了分析，結(jié)果如圖4 所示。從圖4 可以觀察到，隨著字典大小從200 到600 的增加，驗證集均方誤差逐漸降低，字典大小從600 增加800，驗證集誤差升高；并且隨字典大小的增加，訓練時間顯著增長。因此選擇字典大小為600。

圖4 SCDNN 隨字典大小變化驗證集誤差及訓練時間示意圖Fig.4 Changes of the validation error and training time of SCDNN with dictionary size

4 結(jié)束語

本文通過提出了一種基于稀疏編碼深度學習網(wǎng)絡的彌散微循環(huán)模型參數(shù)估計方法來估計IVIM 模型參數(shù)f、D、D*，該方法借助深度學習方法，結(jié)合彌散微循環(huán)模型的物理原理，可以通過少數(shù)b值的q空間信息擬合，有效縮短了采集數(shù)據(jù)所消耗的時間，并且能夠比現(xiàn)有的參數(shù)估計方法有更小的誤差，同時從多中心比較實驗來看，該算法具有良好的泛化性能。