基于深度學(xué)習(xí)的心血管血流動力學(xué)參數(shù)反演

周 陽 潘 怡 崔 暢 陳明龍 孫蓓蓓

(1東南大學(xué)機械工程學(xué)院，南京 211189)(2南京醫(yī)科大學(xué)第一附屬醫(yī)院心血管內(nèi)科，南京 210029)

血流動力學(xué)建模是研究心血管活動的重要技術(shù)手段，其通過計算機數(shù)值模擬，以一種快速、準(zhǔn)確的方式對血液流動現(xiàn)象進行表征和預(yù)測[1-2].不同形式的心血管模型已被開發(fā)并被用于各種臨床場景，如零維模型[3]、一維模型[4]以及三維模型[5]等.零維-一維耦合模型由于能以較低的計算成本準(zhǔn)確提供脈搏波形態(tài)以及局部位置的壓力和速度信息，得到學(xué)者們的廣泛關(guān)注，并已被成功應(yīng)用于周身循環(huán)和局部血流動力學(xué)的研究[6].為了提高物理模型的預(yù)測精度，根據(jù)患者臨床實際測量數(shù)據(jù)對其進行個性化建模尤為重要，相關(guān)的參數(shù)反演策略已成為心血管血流動力學(xué)領(lǐng)域的研究熱點[7-8].

目前心血管參數(shù)反演方法主要包括參數(shù)優(yōu)化法和濾波法.參數(shù)優(yōu)化法以模型輸出與實際測量的差值構(gòu)建損失函數(shù)，并利用牛頓法或LM等算法通過最小化損失函數(shù)來求出未知參數(shù).Itu等[9]對基于牛頓法的心血管模型參數(shù)反演進行了大量的研究；Zhang等[10]對真實病人的臨床數(shù)據(jù)采用LM算法進行了參數(shù)反演實驗，研究結(jié)果表明個性化模型的輸出與實際測量吻合度良好.濾波法是一種動態(tài)遞歸的參數(shù)估計方法，基于卡爾曼濾波(KF)原理，利用每個時間步的卡爾曼增益來進行參數(shù)迭代更新[11].根據(jù)心血管數(shù)值模型的特點，學(xué)者們相繼提出利用降價無跡卡爾曼濾波[12]和集合卡爾曼濾波[13]來求解血流動力學(xué)的反問題，并取得了積極的成效.然而，這些方法在實際使用中存在迭代計算量大、初始值設(shè)置不合理以及迭代發(fā)散等情況，往往需要人工干預(yù)來調(diào)整初始參數(shù)和收斂條件，難以取得明顯的效率提升.此外，當(dāng)存在多個心動周期的測量數(shù)據(jù)時，上述方法并不能對這些數(shù)據(jù)進行較好的綜合利用.

由于深度學(xué)習(xí)技術(shù)的迅速發(fā)展，利用人工神經(jīng)網(wǎng)絡(luò)進行非線性映射為參數(shù)反演提供了新的思路，并在不同領(lǐng)域的參數(shù)反演任務(wù)中得到了初步的應(yīng)用[14].臨床測量中常涉及包括離散數(shù)值以及連續(xù)波形在內(nèi)的不同形式的數(shù)據(jù)，傳統(tǒng)的神經(jīng)網(wǎng)絡(luò)對這些復(fù)雜的數(shù)據(jù)難以取得滿意的映射效果.因此，提出一種混合多源輸入的深度神經(jīng)網(wǎng)絡(luò)(DNN)模型，根據(jù)不同類型的測量數(shù)據(jù)，分別采用卷積神經(jīng)網(wǎng)絡(luò)與全連接神經(jīng)網(wǎng)絡(luò)來提取其對參數(shù)的依賴特征.同時，針對波形噪聲的干擾，采用一種集成網(wǎng)絡(luò)模型來進一步提高參數(shù)反演的精度.通過引入一個多尺度血流動力學(xué)模型對所提方法進行有效性驗證，并討論了不同水平的噪聲對反演精度的影響.

1 多尺度血流動力學(xué)模型

建立了一個零維-一維耦合的多尺度心血管血流動力學(xué)模型.如圖1所示.圖中，R0和R分別為末梢血管近端和遠端阻力；C為末梢血管的順應(yīng)性；T為心動周期；T1為心臟射血時間；Qmax為入口流率最大值；QAV為主動脈入口流率.一維模型用于描述血流在55條動脈內(nèi)的速度和壓力情況，而零維模型則用于描述外周小血管及毛細血管的阻力和順應(yīng)性.假設(shè)動脈內(nèi)的血流是黏性不可壓縮的，一維模型的控制方程如下[15]：

(a)55條動脈示意圖

(1)

(2)

式中，t為時間；x為動脈長度方向的坐標(biāo)；A為血管的橫截面積；P為該橫截面的平均壓力；U為相應(yīng)的平均流速；Kr為血流的黏性阻力系數(shù)，這里將其設(shè)為 8πν，其中動力黏度ν為4.43 cm2/s；ρ為血液的密度，其值假定為常數(shù)1.06 g/cm3.

采用彈性本構(gòu)方程來定義血流與動脈壁的耦合關(guān)系，即

(3)

(4)

式中，P0為參考壓力；A0為壓力處于P0時的血管橫截面積；β為血管的剛度參數(shù)；E和h分別為血管壁的彈性模量和壁厚.模型的主動脈入口流率QAV采用如圖1(b)所示的近似流量曲線，其中，Qmax與心臟每搏輸出量SV有如下關(guān)系[10]：

(5)

出口邊界采用一個三元素Windkessel模型來描述血流在末端血管的行為，如圖1(c)所示，其控制方程如下：

(6)

式中，Q為血管內(nèi)截面的平均流量，且有Q=AU.采用Maccormack格式對上述一維模型的雙曲偏微分方程進行數(shù)值離散[16].根據(jù)Zhang等[10]提出的參數(shù)估計設(shè)置，將55條動脈分為3個集群，如圖1(a)所示，其中臂動脈集群包括標(biāo)號為3、4、7、8、9、15、17、18、19、43、44、45、46的動脈段，頸動脈集群包括標(biāo)號為5、6、11、16、39、40、47、48 的動脈段，其他動脈段則被歸為主動脈集群.分別采用3個縮放系數(shù)Cβ_arm、Cβ_car和Cβ_aor對上述3個集群的剛度參數(shù)在標(biāo)稱值的基礎(chǔ)上進行縮放.采用參數(shù)CR用于統(tǒng)一縮放邊界阻力R.

模型中涉及的動脈幾何參數(shù)、剛度以及外周血管參數(shù)來自文獻[10].執(zhí)行參數(shù)反演所需的臨床測量包括心率(RH)、右肱動脈(標(biāo)號7)、右股動脈(標(biāo)號38)的壓力波形以及右頸動脈(標(biāo)號5)的血流速度波形.心臟射血時間T1根據(jù)測得的心率和已有的回歸公式進行計算[17].

2 基于DNN的參數(shù)反演模型

2.1 卷積神經(jīng)網(wǎng)絡(luò)模型

采用卷積神經(jīng)網(wǎng)絡(luò)(CNN)來提取連續(xù)生理波形的特征信息，如圖2所示.如圖2(a)所示，將一維的壓力波形W1、W2以及血流波形W3堆疊為一個二維的波形特征矩陣W，即

(7)

對于獲得的特征矩陣W，采用連續(xù)的多個卷積層對其進行特征提取，每一層之間的卷積運算如圖 2(b)所示，r+1層上的第k個卷積圖的計算可以表示為[18]

(a)生理波形特征矩陣

(8)

對于每一層卷積層(CONL)網(wǎng)絡(luò)的輸出，采用Relu激活函數(shù)來增強其非線性，如圖2(c)所示，其函數(shù)表達式為

(9)

式中，X為Relu函數(shù)的輸入值.在激活函數(shù)后需要對獲得的特征圖進行特征選擇與特征過濾，以降低其維度，常見的操作包括池化運算和增大卷積步長.考慮波形矩陣在長度方向的尺度明顯大于高度方向，故采用增大長度方向的卷積步長s1的方式來降低矩陣維度.

2.2 全連接神經(jīng)網(wǎng)絡(luò)

為綜合考慮臨床測量中的離散數(shù)值，采用圖3所示的全連接神經(jīng)網(wǎng)絡(luò)(FCNN)來提取心率、脈搏波傳輸時間(TPWT)與待估計參數(shù)之間的依賴性.采用一種切線法來自動計算壓力波形之間的TPWT值，如圖3(a)所示，取波形最低點切線與上升支最大一階導(dǎo)數(shù)點切線的交點為特征點，2個壓力波形特征點之間的時間差則為TPWT值[19].

單個神經(jīng)元對輸入的作用如圖3(b)所示，其中r+1層第g個神經(jīng)元的輸出Or+1,i可表示為

(10)

(a)特征點計算示意圖

2.3 血流動力學(xué)參數(shù)反演DNN模型

如圖4所示，建立了一個用于心血管血流動力學(xué)參數(shù)反演的混合多源輸入深度神經(jīng)網(wǎng)絡(luò)(DNN)模型.圖中，W為連續(xù)血壓、血流波形，V為離散特征向量；CONL1(16×3×10)表示第1層卷積層，且其卷積核個數(shù)為16，高度和長度尺寸分別為3和10；FCL1(6×1)表示第1層全連接層，其神經(jīng)元個數(shù)為6，其余卷積層和全連接層的參數(shù)可以此類推.卷積核在高度方向的步長為1，長度方向的步長為3.壓平層用于將卷積得到的特征圖重排列為一維向量，以便與全連接網(wǎng)絡(luò)連接.融合層用于將2部分輸入經(jīng)各自網(wǎng)絡(luò)處理后得到的向量串聯(lián)，用于提取它們對輸出的共同依賴性.采用預(yù)測值與實際值的均方誤差(MSE)作為損失函數(shù)，同時采用反向傳播算法來進行網(wǎng)絡(luò)權(quán)重的更新.

圖4 DNN模型結(jié)構(gòu)

采用正向運行心血管物理模型的方式來獲取可用樣本.用于采樣的參數(shù)如表1所示.其中Cβ_arm、Cβ_aor、Cβ_car和CR為待估計參數(shù)，均為無量綱值.RH為已知的測量值，在生理范圍內(nèi)改變RH和SV值以模擬不同的入口流量.

表1 參數(shù)抽樣范圍

在這些參數(shù)構(gòu)成的參數(shù)空間內(nèi)，采用Sobol法進行參數(shù)抽樣，并將得到的參數(shù)樣本輸入心血管物理模型中進行正向計算，一個參數(shù)樣本與其對應(yīng)的模型輸出波形便構(gòu)成一個有效樣本.共生成10 000個樣本用于DNN模型訓(xùn)練，其中1 000個用于驗證.額外生成200個樣本用于網(wǎng)絡(luò)的測試與評價.

2.4 集成網(wǎng)絡(luò)模型

考慮到實際測量的波形常會受到噪聲的干擾，采用一種集成網(wǎng)絡(luò)模型來減小噪聲對預(yù)測的影響.對于同一個DNN模型，重復(fù)對其進行基于隨機初始化和隨機批樣本的網(wǎng)絡(luò)訓(xùn)練，可得到一批權(quán)重參數(shù)收斂到不同值的DNN模型.在測試階段，同時采用這些DNN模型進行預(yù)測，并將各模型的輸出進行簡單平均，其所得結(jié)果即為集成網(wǎng)絡(luò)模型的輸出.該模型可表示為

(11)

3 參數(shù)靈敏度分析及反演結(jié)果

3.1 參數(shù)靈敏度分析

圖5給出了3處待測量部位的生理輸出對不同參數(shù)的靈敏度.分別觀察了Cβ_arm擾動對W1的影響、Cβ_car擾動對W3的影響、Cβ_aor和CR的擾動對W1和W2的影響，參數(shù)擾動統(tǒng)一設(shè)置為從1減小到0.7和從1增加到2.

如圖5所示，Cβ_arm對肱動脈壓力波形的主峰影響不大，對次峰影響明顯，而Cβ_aor對2個壓力波形的主峰和次峰都有著明顯的影響.Cβ_car對頸動脈流速波形的震蕩范圍有較大影響，當(dāng)Cβ_car增大時，波峰和波谷都有明顯的減弱.CR與壓力呈現(xiàn)正相關(guān)，能夠使壓力波形整體沿著垂直軸發(fā)生移動，且變化較為顯著.由此可發(fā)現(xiàn)，Cβ_arm對于W1的影響相對較弱，而其余3個參數(shù)對于相應(yīng)波形的影響都較為明顯.

(a)不同Cβ_arm值下的W1

3.2 參數(shù)反演結(jié)果

在生理波形中添加白噪聲N(0,σ2)來模擬實際測量的誤差, 標(biāo)準(zhǔn)差σ依次從0增加到1.5.圖6給出了單個DNN模型和集成網(wǎng)絡(luò)模型在不同噪聲干擾下的各參數(shù)反演情況，縱坐標(biāo)為200個測試樣本的參數(shù)真實值與反演值的MSE值.由圖6(a)可見，當(dāng)噪聲標(biāo)準(zhǔn)差σ增加不超過0.5時，參數(shù)反演的精度變化相對平穩(wěn)，且各參數(shù)MSE值相差不大.當(dāng)σ進一步增加時，參數(shù)MSE值有明顯上升，且Cβ_arm的MSE值上升的幅值明顯高于其他3個參數(shù).當(dāng)σ增加到1.5時，各參數(shù)的MSE值達到最大，其中Cβ_arm的MSE值為0.028，遠高于其他3個參數(shù)，其次是Cβ_car的MSE值，為0.008，Cβ_aor和CR的MSE值相對較小，分別為0.004和0.006.

可發(fā)現(xiàn)，在測量噪聲增加時，參數(shù)Cβ_arm的反演精度會有顯著下降.結(jié)合3.1節(jié)靈敏度分析不難發(fā)現(xiàn)，Cβ_arm對壓力波形的影響較小，在噪聲的干擾下，波形次峰的特征可能會有所減弱，導(dǎo)致DNN網(wǎng)絡(luò)對其識別能力有所降低.

圖6(b)給出了集成網(wǎng)絡(luò)模型的反演結(jié)果.可以發(fā)現(xiàn)，各參數(shù)的變化趨勢與單個DNN網(wǎng)絡(luò)的結(jié)果基本相符，但在噪聲幅值較大時，集成模型預(yù)測結(jié)果的MSE值明顯低于單個DNN網(wǎng)絡(luò)結(jié)果.其中，當(dāng)σ=1.5時，Cβ_arm的MSE值減小為0.015，Cβ_aor、Cβ_car以及CR的MSE值也分別減小為0.003、0.005和0.003.

(a)單個DNN模型

進一步地，選取一組特定參數(shù)設(shè)置來分析所提DNN法的反演效果，其中，Cβ_arm、Cβ_aor、Cβ_car和CR的真實值分別為2.18、1.84、1.67和1.14，而RH和SV分別為70.87和72.65 mL.圖7給出了噪聲方差σ為1.5時集成網(wǎng)絡(luò)模型中各DNN模型反演值的箱線圖，圖中線框值從上到下分別代表最大值、上四分位數(shù)、中位數(shù)、下四分位數(shù)、最小值，星號標(biāo)記代表剔除的異常值.從圖7可看出，正常值的分布較為集中，中位數(shù)的位置相對合理，且異常值數(shù)量較少，說明各DNN網(wǎng)絡(luò)模型的預(yù)測性能基本一致，受波形噪聲的影響，估計值在真實值附近波動.

圖7 各DNN模型預(yù)測值箱線圖

為進一步評估所提方法的參數(shù)估計效果，將DNN法與現(xiàn)有卡爾曼濾波(KF)法所得估計值進行對比，分別給出了各參數(shù)在不同噪聲下的誤差，如表2所示.可看出，噪聲水平對于2個方法的估計誤差有較大的影響，噪聲增大時，估計誤差總體呈現(xiàn)上升趨勢；DNN法的估計誤差明顯低于KF法，兩者的誤差對比在較低噪聲水平下并不顯著，但隨著噪聲的提高，DNN法的估計精度具有顯著的優(yōu)勢.

表2 DNN法與KF法反演誤差的對比 %

對于血流動力學(xué)參數(shù)反演，個性化后的模型輸出與實際測量波形的誤差是評估參數(shù)反演是否有效的重要標(biāo)準(zhǔn).因此，在噪聲方差σ為1.5的條件下，分析了模型輸出波形與實際測量波形的誤差，圖8給出了這組參數(shù)的模型仿真波形與測量波形的對比.可以發(fā)現(xiàn)，2種方法個性化后的仿真波形與測量波形基本一致.計算可得，DNN法模型的模擬輸出與測量波形均方根誤差(RMSE)分別為242 Pa、233 Pa和2.21 cm/s，而對應(yīng)的KF法波形RMSE為665 Pa、595 Pa和3.23 cm/s.可看出，所提DNN法可以更好地擬合實際測量結(jié)果.

(a)波形W1

4 討論

基于深度學(xué)習(xí)技術(shù)，提出了一種新的血流動力學(xué)參數(shù)反演方法，并通過數(shù)值模擬對其進行了初步驗證.結(jié)果表明，相比于現(xiàn)有技術(shù)，本方法同時提高了參數(shù)反演的效率與精度，能夠?qū)崿F(xiàn)持續(xù)的動態(tài)監(jiān)測，具有較高的實際意義.現(xiàn)階段的研究存在以下不足：① 參數(shù)的可識別性包括靈敏度分析和相關(guān)性分析，考慮到相關(guān)性分析在文獻[10]中已經(jīng)給出，故本文未對此部分展開敘述；② 利用臨床中較為常見的3處波形對包括剛度與阻力在內(nèi)的4個參數(shù)進行反演，缺少局部血管的參數(shù)反演測試.臨床應(yīng)用時可根據(jù)反演任務(wù)的不同建立多個網(wǎng)絡(luò)模型，以滿足周身及局部個體化建模的需要；③ 考慮到血管長度、直徑等參數(shù)也會對血流速度與壓力產(chǎn)生影響，在應(yīng)用所提方法時，需要對已經(jīng)訓(xùn)練好的模型進行少量的遷移學(xué)習(xí)，未來工作中將進一步研究遷移算法對反演精度的影響.

5 結(jié)論

1)基于深度學(xué)習(xí)技術(shù)提出了一種新的心血管血流動力學(xué)參數(shù)反演方法，并提出了一種集成網(wǎng)絡(luò)模型來降低測量噪聲對參數(shù)反演的影響.

2)對所提方法進行數(shù)值模擬驗證，結(jié)果表明，所提方法可以準(zhǔn)確地從帶有噪聲的測量中反演出目標(biāo)參數(shù)，且反演過程無需人工操作和模型正向計算，兼顧了反演效率與反演精度.

3)未來工作將對所提方法進行基于遷移學(xué)習(xí)算法的優(yōu)化與擴展，使該方法能夠有效應(yīng)用于不同個體.