基于位移流U-Net 和變分自動編碼器的心臟電影磁共振圖像左心肌運動追蹤*

王甜甜 王慧 朱艷春 王麗嘉?

1) (上海理工大學醫療器械與食品學院,上海 200093)

2) (中國聯通醫療基地,官洲生命科學創新中心,廣州 510000)

心血管疾病(cardiovascular diseases,CVDs)的高發病率和高死亡率已經嚴重影響了人類的生存質量.如何評估心臟功能、輔助臨床CVDs 診療和預后評估,是一個迫切需要解決的問題.針對這個問題,本文在前期心臟電影磁共振(cardiac cine magnetic resonance,CCMR)圖像左心肌分割的基礎上,提出一種基于位移流U-Net(DispFlow_UNet)和生物力學變分自動編碼器(variational autoencoder,VAE)的左心肌運動追蹤方法:DispFlow_UNet_VAE.主要研究內容有:1) 搭建壓縮激勵殘差U-net 網絡精準分割左心肌,根據分割結果計算心室體積、心肌質量等,評估心臟整體功能;2) 根據DispFlow_UNet_VAE 估計CCMR 圖像連續幀之間的左心室運動,結合左心肌分割掩膜得到左心肌密集位移場;3)利用模擬數據真實位移場、臨床數據集對追蹤結果進行對比和評估.結果表明,本文追蹤算法具有較高的精度和泛化能力.

1 引言

據世界衛生組織(world health organization,WHO)估計,每年全球有1790 萬人死于心血管疾病(cardiovascular disease,CVDs)[1],如果按照目前的趨勢繼續下去,2030 年全球將有2360 萬人死于CVDs[2].因此,如何準確評估心臟功能對于疾病的早診斷、早治療具有重要臨床意義.為了評估心臟功能,使用了兩類指標:全局指標[3?5]和局部指標[6,7].全局指標包括射血分數、心肌質量、心室容積等,在CVDs 診斷中具有重要意義.然而,該類指標不能用于局部分析,無法反映細微的室壁運動異常,因此,在整體評價心臟功能的同時,需通過探索心肌運動規律,進一步分析心肌局部功能變化.

目前,借助醫學成像手段評估心臟結構及功能指標是臨床診斷CVDs 的主要依據.其中,臨床常用的技術為心動超聲圖和CCMR 成像.超聲心動圖已廣泛應用于應變分析,但其結果的準確性取決于操作技巧,并受到狹窄聲窗的限制,在評估胸骨后的心臟時效果欠佳[8].CCMR 成像由于無電離輻射、軟組織對比度高、可多參數和多方位成像等優點而被廣泛應用于臨床[9].其中,基于穩態自由進動(steady state free precession,SSFP)序列的CCMR 成像可動態、精確地描繪心肌、心室及其內部的復雜結構,是無創性量化心臟功能的金標準技術[10].一個典型的基于SSFP 序列的CCMR 圖像如圖1 所示,代表由舒張末期(end diastole,ED)至收縮末期(end systole,ES),再至ED 的一系列相位組成的心臟運動周期,其中心肌在很大程度上是同質的,不能直接提供運動模式的信息[11],因此,需要通過追蹤算法來獲取心肌的運動信息,從而得到運動位移場,以輔助醫生早診斷早治療CVDs.

圖1 基于SSFP 序列的CCMR 圖像Fig.1.CCMR images based on SSFP sequence.

由于心臟的收縮是非剛性,心臟追蹤通常在像素級水平上進行.在CCMR 圖像中,一個相位的運動通常是相對參考相位來估計的,設t 相位的圖像為I(x,y,t),參考相位的圖像為I(x,y,tref).運動追蹤的目的是找到映射 Fθ[12]:

其中,Fθ是參數為θ 的映射函數,Vx和 Vy分別是沿x 和y 方向的位移矢量.

運動跟蹤方法通常可以根據 Fθ的不同進行分類:基于強度的方法[13]、光流法[14]、基于心膜輪廓的算法[15]、基于非剛性配準算法[16?18]和基于深度學習(deep learning,DL)的方法[19?23]等.這些方法在運動估計中呈現出不同的計算負擔和精確度,其中傳統的追蹤算法缺乏完全自動化,需要加入符合心肌運動模式的約束,限制了復現性,而基于DL 的追蹤算法已被證明具有高效性、精確性和可重復性.

研究表明,心臟運動估計可以重新表述為一個可學習的問題.Qin 等[19]提出一種結合分割和運動估計的多任務框架,學習的心臟運動場用于扭曲分割掩膜,并以半監督方式引導分割模塊,結果表明,分割和運動估計性能都有所提高.Zheng 等[20]提出一種表觀流網絡,它是一種改進的U 型網,在訓練中使用分割掩膜來改進運動估計.然而,DL追蹤算法使用常見的機器學習技術(如所有可學習參數的L2 正則化[20,22]、運動估計的直接正則化(如平滑懲罰[19]、解剖感知[23]))作為約束項,缺乏應用特定的先驗知識來指導優化.因此,本文提出一種DispFlow_UNet_VAE 左心肌追蹤算法,并使用基于生物力學的先驗知識來約束學習,其可以隱式地學習正則化行為以指導優化,提高算法的高效性及復現性.

2 方 法

2.1 基于壓縮激勵殘差U-net 的左心肌分割方法

精準的心肌分割是運動追蹤的基礎.Ronneberger 等[24]提出的U 形網絡(U-Net)憑借所需訓練集數量少、分割模糊邊界精度高等優點而廣泛應用于左心肌分割中,此后多種U-Net 改進算法應運而生,本文采用王慧[25]提出的壓縮激勵殘差U-net (squeeze and excitation residual U-net,SERU-net)用于全心肌的自動分割,如圖2 所示.該網絡融合了壓縮激勵(squeeze-and-excitation,SE)模塊和殘差模塊,其中SE 模塊使網絡在提取特征時能夠通過學習自動獲取每個特征通道的重要程度,更好提取左心肌有效特征;殘差模塊有效抑制了梯度消失和訓練過擬合問題,使得信息前后向傳播更加順暢,使網絡學習效果更好,分割精度更高.

圖2 壓縮激勵殘差U-net 網絡結構Fig.2.Squeeze-and-excitation residual U-shaped network.

圖3 給出了測試集中某一例數據從基底到頂端的分割結果,其中黃色為專家手動分割的金標準,紅色為SERU-net 分割結果.從圖3 可以看出,SERU-net 能夠更好適應左心肌形狀變化,精確勾畫出左心肌.

圖3 SERU-net 左心肌分割結果Fig.3.Results of left myocardium segmentation by SERU-net.

2.2 運動追蹤算法

深度學習方法通常可分為有監督訓練策略與無監督訓練策略.在醫學圖像中,存在大量未標記數據和少量已標記數據,半監督學習被認為是充分運用可獲得數據進行自動訓練的有效方法[26].因此本文采取半監督深度學習策略.圖4 給出了DispFlow_UNet_VAE 框架,其 由位移流U-Net (2.2.1 節)和VAE 正則化器(2.2.2 節)兩部分組成.首先,VAE 用于學習模擬變形的概率分布以捕捉潛在的生物力學特征;接著,將學習好的VAE 用作位移流U-Net 配準網絡的正則化函數,這樣做的優點是,使解空間正則化并且有助于預測生物力學上合理的運動特征,而不需要任何顯示的懲罰項.

圖4 DispFlow_UNet_VAE 運動追蹤框架Fig.4.The motion tracking architecture of DispFlow_UNet_VAE.

2.2.1 DispFlow_UNet 網絡

如圖5 所示,DispFlow_UNet 是U-net 的變體.在給定同一切片上的目標相位和參考相位作為輸入的情況下,DispFlow_UNet 產生兩者之間像素級的位移矢量場(displacement vector field,DVF),結合分割掩膜,從位移場中提取心肌運動特征.

圖5 DispFlow_UNet 網絡框架Fig.5.The network architecture of DispFlow_UNet.

由于無法獲取圖像之間的真實密集位移場,因此網絡通過時間跟蹤空間特征,依靠時間強度變化作為自我監督.將兩個輸入圖像在位置P=(x,y)處的像素強度記為Is(P)和It(P),DispFlow_UNet產生s(源)與t(目標)相位之間的位移圖記為Ft,Ft(P)=(Ftx(P),Fty(P)),通過空間變化器將Ft進行圖像重建,使得以下強度差異最小化:

2.2.2 VAE 正則化器

為了進一步確保DispFlow_UNet 產生的位移場符合心肌運動實際,本文提出學習模擬變形概率分布的VAE 正則化網絡,以隱含地學習生物力學上的形變規律,如圖6 所示.VAE[27,28]被用作一個基于學習的正則化器,以模擬生物力學似是而非的變形的概率分布,其編碼器用于將變形映射到潛在變量,通過高斯分布進行正則化,并通過解碼器解碼為變形場.VAE 正則化器被訓練于重建生物力學模擬變形的一階導數,以消除剛性平移的任何影響.用 Φ=[u,v]∈R2×H×W表示二維(2D)空間中的變形,其中,u,v 分別表示沿x 和y 方向的位移.變形場的一階梯度表示為

圖6 VAE 網絡Fig.6.The network architecture of VAE.

其中H 和W 代表空間的維度.VAE 的重構損失公式為

其中,?Φ′表示輸入 ?Φ 的重建項,z 代表由VAE編碼的潛在矢量,p(z) N(0,I) 代表先驗的高斯分布,qθ代表由θ 參數化的編碼器,DKL代表KL 散度(Kullback-Leibler divergence),α 是控制重建質量和潛在空間規則程度之間折衷的超參數.

VAE 損失 RVAE提供一個定量的度量指標,用于確定變形在生物力學上的可信性.靠近已學習的VAE 潛在流形的解將產生較低的 RVAE,而遠離流形的解將給出較高的 RVAE.

2.2.3 損失函數

密集追蹤網絡DispFlow_UNet_VAE 的目標損失函數由圖像強度差異 LIMG(Ft) 和VAE 重構損失 RVAE(?Φ ⊙M) 兩部分組成:

其中,M 是分割掩膜的二值圖,僅用于調整感興趣區域,⊙ 代表哈達瑪積(Hadamard product).γ 是超參數,權衡了圖像的相似性和變形的物理合理性.

3 實 驗

3.1 實驗數據

該方法在Automated Cardiac Diagnosis Challenge (ACDC)數據集[29]上進行訓練、驗證和預測;在臨床數據集上進行模型泛化能力評估;在模擬MRI 數據集[30]上進行對比實驗.實驗數據的具體情況如表1 所列.

表1 實驗數據Table 1.The experimental data.

ACDC 數據集包括100 名受試者(健康和患病病例)的短軸CCMR 圖像,使用兩種不同主磁場的磁共振掃描儀(1.5T Siemens Area 和3.0T Siemens Trio Tim)在六年內通過SSFP 序列收集得到.所有數據采集都是屏氣進行的,以確保在視頻中只能觀察到心臟運動.每個受試者包含9—10層圖像,相位數在12—35 之間變化.總體來說,該數據集有951 個,每個序列提供ED 和ES 相位的金標準,如圖7 所示.此外,100 名受試者被平均分為5 類,每類20 名受試者,分別為正常病例(normal subjects,NOR)、伴有埂塞的收縮性心力衰竭(myocardial infarction,MINF)、擴張性心肌病(dilated cardiomyopathy,DCM)、肥厚型心肌病(hypertrophic cardiomyopathy,HCM)和右心室異常(right ventricle abnormality,ARV).

圖7 ED (左)與ES (右)的原始圖像及金標準Fig.7.Original image and its ground truth of ED (left) and ES (right).

75 例臨床數據集在 GE 1.5T 磁共振掃描儀通過 SSFP 序列獲得,所有數據均符合倫理要求.采集參數是:圖像大小為256×256、視野大小(field of vision,FOV)為360 mm×360 mm、層厚為6—8 mm、層間距為2—4 mm、重復時間(repetition time,TR)為 3.5 ms、回波時間(echo time,TE)為1.4 ms、層數為6—10 層、相位數為20—28 個.

3.2 實驗流程

實驗將從數據分配、數據預處理、訓練與驗證三步驟進行.

3.2.1 數據分配

在數據驅動的機器學習中,常假設訓練樣本和測試樣本來自同一個分布(內部分布),然而測試集中的外部分布通常會給出較差的模型泛化.本文按照留下一種疾病的方法在ACDC 數據集上進行了五重交叉驗證,根據已知的疾病將內部分布和外部分布分開.由于不同疾病之間存在顯著的心臟解剖和動力學差異,因此與其他4 種疾病相比,一種疾病類別可被視為外部分布.

對于內部分布中的受試者,將其分為訓練集(80*80%=72 名)、驗證集(80*10%=8 名)和預測集(80*10%=8 名),外部分布中的所有受試者(20 名)被用作測試.在訓練過程中,排除了覆蓋LV 腔的所有切片的頂部20%和底部20%,而選擇了中間剩余的60%,該設計旨在進一步減少平面外運動的影響.表2 給出了詳細的數據分配情況.

表2 ACDC 數據集分配情況Table 2.ACDC data set allocation.

3.2.2 數據預處理

為了避免周圍結構對左心肌的影響,需要在原始圖像上確定感興趣區域(region of interest,ROI),所有訓練圖像被裁剪為96×96 大小,同時將強度歸一化到[0,1]的范圍,并通過隨機旋轉、平移和縮放等操作擴充數據,提高網絡精確性及魯棒性.

3.2.3 訓練與驗證

網絡搭建采用Pytorch 框架,操作系統為Windows 10.通過(3)式中的損失函數來訓練網絡參數,選取Adam 作為優化器,學習率=0.0001,批量大小為16,訓練次數epoch=300,參數 γ=0.001,通過驗證集選擇得出.利用訓練好的模型對預測集進行預測,根據預測結果進行評估.

3.3 評價方法

從三個方面評估本文追蹤方法的精確性和泛化能力.

1)利用幾何指標值量化本文算法的追蹤精度.一般來說,手動制作心臟運動的參考標準是不可能的,為了進行定量分析,通常利用分割掩膜作為獨立的參考標準.應用訓練好的密集追蹤網絡去產生FES,使用FES來扭曲ES 幀的分割金標準MES,得到 MES°WFES,再利用Dice 系數、Hausdroff 距離(HD)和心肌輪廓距離(myocardial contour distance,MCD)測量其與ED 幀心膜掩膜圖之間的重疊程度,Dice 和HD 的公式見(4)式和(5)式.

其中,A 是變形得到的ED 幀心肌輪廓,G 是ED幀掩膜,D(.,.)為歐幾里得距離.Dice 系數取值越高、MCD 和Hausdroff 距離取值越低代表追蹤精度越高.

2)利用模擬數據位移場驗證本文方法的可行性.計算預測位移場與模擬數據位移場之間的平均角度誤差(average angle error,AAE)和平均終止點誤差(average endpoint error,AEPE),測量兩位移場之間的重合程度.

設圖像某像素點P 的角度誤差AE(P)是預測位移矢量 u0(P)=(x0,y0) 和位移真值矢量u1(P)=(x1,y1)在點P 二維空間的角度,則:

AAE 為計算區域中AE(P)的均值,其取值越小表明追蹤結果越準確.

圖像中某像素點P 的終止點誤差EPE(P)表示兩個位移終止點 u0(P)=(x0,y0)、u1(P)=(x1,y1)之間的歐氏距離:

AEPE 為計算區域中EPE(P)的均值,其取值越小表明追蹤結果越準確.

3)利用臨床數據集(即外部數據)驗證本文模型的泛化能力.使用訓練好的模型直接預測臨床數據集,通過對比模型在內部數據和外部數據集上的性能評估模型的泛化能力.

4 結果與討論

4.1 全局功能評估

根據SERU-net 方法分割左心肌,利用面積長度法(即:每一層面心室面積與層厚(包含層間距)的乘積)計算心室體積(EDV,ESV)、每搏輸出量(SV)、射血分數(EF)、左心室質量(ED_LVM,ES_LVM)等指標,將指標值與正常范圍值進行比較,可以輔助臨床醫師評估心臟整體功能.18 名受試者各功能指標數據的均值和標準差如表3 所列.

表3 左心室功能指參數 (均值±標準差)Table 3.Left ventricular function parameters (Mean±standard deviation).

4.2 追蹤結果評估

基于目標損失函數((3)式)對密集追蹤網絡進行訓練,網絡訓練 300 個epoch 后得到的訓練集與驗證集的損失曲線如圖8 所示,可以看出,所提出的網絡可以較好地擬合,驗證集的損失率可以降到 0.02 左右.

圖8 訓練集(藍)與驗證集(橙)的損失曲線Fig.8.Loss curves of training set (blue) and verification set(orange).

利用預測位移場(圖9 中DVF)將ES 幀掩膜重采樣并與ED 幀金標準掩膜進行比較,結果如圖9 所示,可以看出,扭曲后的圖像(圖9 中Warped_label_ES)與金標準(圖9 中label_ED)的LVC 和LVM 高度重合,表明了左心肌運動追蹤的精度.

圖9 利用預測位移場將ES 扭曲至ED 的示例圖Fig.9.Example diagram of warping ES to ED using the predicted displacement field.

此外,預測位移場可以輔助臨床醫生診斷不同類型的CVDs,以達到早診斷早治療的目的.圖10給出了幾個例子,不難發現,網絡追蹤得到的位移場確實足夠好,足以表征病理類別中典型病例的心臟運動.其中正常人(NOR)的位移場在整個LVM上有大致相同的振幅,表明MOR 病例LVM 的收縮和增厚是同步的;對于MINF 病例,可以看到LVM 上位移流不均勻,部分心肌節段收縮和增厚明顯小于其他節段,這是MINF 病例的典型癥狀;對于HCM 病例,LVM 上的位移流過大,這意味著LVM 收縮和增厚過度;相反,DCM 病例LVM 上的位移流很小,這是因為DCM 的心臟通常沒有足夠的收縮和增厚.

圖10 不同病例類型預測的位移圖Fig.10.Predicted displacement fields of different case types.

4.3 對比結果

本文通過兩個方面展開對比實驗,首先將本文方法與其他追蹤算法所得的結果進行對比,其次將預測得到的位移場與模擬數據位移場進行對比.

1)為了驗證本文追蹤方法的準確性,將本文方法與二維B 樣條自由變形(free-form deformation,FFD)配準方法[31]、Qin 等[19]提出的基于深度學習的L2 范數正則化配準網絡(DL+L2)、未引入VAE 正則化器的追蹤方法(DispFlow_UNet)得到的結果進行對比.在5 重交叉驗證的結果中,LVC 和LVM 獲得的Dice 系數、MCD 和HD 的平均值(標準差)如表4,加粗字體表示結果更好,結果顯示DispFlow_UNet 優于DL+L2,FFD 法,而DispFlow_UNet_VAE 相比于DispFlow_UNet進一步取得更好的指標值.同時,各方法計算得到LVM 的Dice 和MCD 箱形圖如圖11 所示,更加直觀地看出本文方法取得最高的Dice 值和最低的MCD 值,表明追蹤精度最優.

圖11 本文方法與其他方法的DM (a)和MCD (b)指標箱形圖Fig.11.Box chart of DM (a) and MCD (b) indicators of the method presented in this paper and other methods.

表4 不同追蹤方法Dice 系數、MCD 和HD 的對比Table 4.Comparison of Dice coefficients,MCD and HD of different tracking methods.

2)將預測得到的左心肌位移矢量u0=(x0,y0)和位移真值矢量 u1=(x1,y1) 進行對比,結果如圖12 所示,可以直觀看出,兩者在方向和幅值上都很相近.為了更加清晰地比較兩者的相似性,圖13 將預測位移矢量場(紅色)與位移真值矢量場(綠色)疊加顯示,放大區域很好地展示了兩者的相似性.除了視覺上直觀的對比,本文還計算了定量評估運動追蹤精度的評價指標AEE 和AEPE.其中,AEE=9.79±7.43,AEPE=2.28±1.19 (均值±偏差)單位分別為度和像素,說明預測位移矢量和位移值真值矢量重合程度較高.

圖12 預測位移矢量與位移真值矢量對比Fig.12.Comparison between the predicted displacement field (left) and the true displacement field (right).

圖13 預測位移矢量(紅色)與位移真值矢量(綠色)對比Fig.13.Predicted displacement field (left) and the true displacement field (right).

4.4 模型泛化能力評估

在對模型進行泛化性研究中,利用ACDC 數據集訓練的模型對15 例臨床數據集直接進行預測,測試結果如表5 所列,其中比較了Dice 和MCD,結果表明本文所提出的方法在ACDC 數據集(內部預測集)和臨床數據集(外部預測集)上都取得了較滿意的結果,表明該方法在未知區域的數據上具有很好的通用性,具備較高的配準精度,這可能是由于生物力學規范的好處,它強制生成的變形場在生物力學上是合理的,并且對域轉移問題不敏感.

表5 臨床數據集上的模型泛化性能Table 5.Model generalization performance on clinical datasets.

5 結論

本文提出一種基于位移流U-Net 和生物力學引導的VAE 左心肌運動追蹤方法DisFlow_UNet_VAE.首先針對CCMR 圖像進行了數據預處理,包括數據歸一化、數據增強和感興趣區域的提取;其次,搭建壓縮激勵殘差U-Net 網絡精準分割ED 和ES 相位的左心肌,為后續的定量計算服務.接著運用DisFlow_UNet_VAE 追蹤心肌運動得到密集位移場,利用位移場將ES 相位的心肌掩膜變形至ED 相位,并計算其與ED 相位掩膜的重疊程度以評估追蹤精度;最終對預測集及臨床數據集進行預測,結果表明,DisFlow_UNet_VAE 能夠精確地追蹤左心肌運動,具有較強的模型泛化能力.密集位移場能夠輔助醫師診斷不同類型的CVDs,并為后續應變分析提供了很好的支持,以實現對心肌組織的活性以及局部心臟功能的評估.該研究將進一步將 2D 網絡推廣到 三維(3D) 網絡,對 3D 數據進行追蹤,更好地滿足臨床需求.