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

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

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

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

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

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

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

(a)55條動(dòng)脈示意圖

(1)

(2)

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

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

(3)

(4)

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

(5)

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

(6)

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

模型中涉及的動(dòng)脈幾何參數(shù)、剛度以及外周血管參數(shù)來自文獻(xiàn)[10].執(zhí)行參數(shù)反演所需的臨床測(cè)量包括心率(RH)、右肱動(dòng)脈(標(biāo)號(hào)7)、右股動(dòng)脈(標(biāo)號(hào)38)的壓力波形以及右頸動(dòng)脈(標(biāo)號(hào)5)的血流速度波形.心臟射血時(shí)間T1根據(jù)測(cè)得的心率和已有的回歸公式進(jìn)行計(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堆疊為一個(gè)二維的波形特征矩陣W，即

(7)

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

(a)生理波形特征矩陣

(8)

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

(9)

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

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

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

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

(10)

(a)特征點(diǎn)計(jì)算示意圖

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

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

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

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

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

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

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

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

(11)

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

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

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

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

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

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

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

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

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

(a)單個(gè)DNN模型

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

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

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

表2 DNN法與KF法反演誤差的對(duì)比 %

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

(a)波形W1

4 討論

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

5 結(jié)論

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

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

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