亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        聯(lián)合深度學(xué)習(xí)的通用血流向量成像方法

        2021-12-07 10:10:02羅婭茹謝盛華尹立雪
        計(jì)算機(jī)應(yīng)用 2021年11期
        關(guān)鍵詞:徑向速度心動(dòng)圖左心室

        彭 博,羅婭茹,謝盛華,尹立雪

        (1.西南石油大學(xué)計(jì)算機(jī)科學(xué)學(xué)院,成都 610500;2.四川省醫(yī)學(xué)科學(xué)院·四川省人民醫(yī)院心血管超聲及心功能科,成都 610072;3.超聲心臟電生理學(xué)與生物力學(xué)四川省重點(diǎn)實(shí)驗(yàn)室,成都 610072)

        0 引言

        心血管疾病是日常生活中嚴(yán)重威脅人體健康的疾病,其發(fā)病率、致死率一直高居首位。根據(jù)世界衛(wèi)生組織(World Health Organization,WHO)的統(tǒng)計(jì),全世界每年大約有1700萬(wàn)人死于心血管疾病,占全世界死亡人數(shù)的31%。心血管疾病之所以有如此高的致死率,是因?yàn)樵谛难芑疾≡缙谕ǔ2粫?huì)出現(xiàn)明顯的臨床表現(xiàn)。所以心臟疾病的防治關(guān)鍵在于早發(fā)現(xiàn)、早治療。其實(shí),病人在患病早期心臟形態(tài)結(jié)構(gòu)和血流模式已與正常人有所不同,可以通過(guò)觀察心臟流體力學(xué)狀態(tài)的變化及時(shí)有效反映心臟功能的異常[1]。因此,心臟流體力學(xué)運(yùn)動(dòng)狀態(tài)的有效可視化觀察和量化評(píng)價(jià)有助于心臟功能異常的早期篩查,為心血管疾病的早發(fā)現(xiàn)提供了一種新的技術(shù)手段。

        超聲具有無(wú)傷、無(wú)輻射、實(shí)時(shí)、價(jià)格低廉、操作簡(jiǎn)單等優(yōu)點(diǎn),因此在臨床上已廣泛使用超聲成像技術(shù)進(jìn)行血流可視化觀察及測(cè)量,例如連續(xù)多普勒、脈沖多普勒、彩色多普勒等。但是,這些方法只能在一定程度上反映流場(chǎng)的部分信息,不能全面反映流體的運(yùn)動(dòng)狀態(tài)。血流向量成像(Vector Flow Mapping,VFM)技術(shù)的出現(xiàn),突破了常規(guī)多普勒超聲的局限,實(shí)現(xiàn)了心臟全流場(chǎng)可視化描述,特別是在左心室運(yùn)動(dòng)的可視化觀察和量化評(píng)價(jià)研究方面被廣泛關(guān)注。VFM 技術(shù)由Ohtsuki 等[2]首次提出,是一種新穎的心功能檢測(cè)技術(shù)。Uejima 等[3]使用三維數(shù)值模型進(jìn)行驗(yàn)證,但是他們的流場(chǎng)速度估計(jì)方法存在理論缺陷。Garcia 等[4]提出將室壁運(yùn)動(dòng)信息結(jié)合多普勒速度的改進(jìn)方法,得到了有效的左心室流場(chǎng)速度矢量。Itatani 等[5]在Garcia 等[4]研究的基礎(chǔ)上進(jìn)行權(quán)重劃分的改進(jìn),并通過(guò)計(jì)算流體動(dòng)力學(xué)(Computational Fluid Dynamics,CFD)得到驗(yàn)證。謝盛華等[6]通過(guò)擴(kuò)展流函數(shù)來(lái)計(jì)算多普勒流量函數(shù),最終實(shí)現(xiàn)了二維平面流體運(yùn)動(dòng)的可視化描述。Asami 等[7]通過(guò)將VFM 與粒子圖像測(cè)速的血流測(cè)量值進(jìn)行比較,定性且定量地評(píng)估了VFM 在三維左心室血流場(chǎng)中的準(zhǔn)確性。Tanaka等[8]提出了一種稱(chēng)為后驗(yàn)VFM準(zhǔn)確度估計(jì)的新方法,并提高了VFM 的臨床有效性。Zhuang 等[9]將YOLO 深度學(xué)習(xí)模型與改進(jìn)的塊匹配算法相結(jié)合,對(duì)左心室壁進(jìn)行定位和跟蹤,并改進(jìn)權(quán)重函數(shù),提高了VFM 的準(zhǔn)確性。但由于現(xiàn)有VFM 技術(shù)首先需要依賴(lài)于特定超聲設(shè)備獲得原始射頻數(shù)據(jù),然后再采用商用散斑跟蹤軟件進(jìn)行流場(chǎng)可視化分析,并且在計(jì)算過(guò)程中需要對(duì)左心室壁進(jìn)行手工勾勒,所以其通用性不足,有一定的局限性和不可移植性。

        基于Garcia等[4]的研究,本文提出了一種聯(lián)合深度學(xué)習(xí)的通用VFM方法。該方法充分利用彩色多普勒超聲心動(dòng)圖像彩色血流信息,從數(shù)字圖像處理的角度出發(fā),聯(lián)合深度學(xué)習(xí)網(wǎng)絡(luò)模型,提取心臟左心室流場(chǎng)信息。相較于現(xiàn)有VFM方法,本文方法主要改進(jìn)包括:1)不依賴(lài)特定超聲設(shè)備導(dǎo)出原始射頻信號(hào)提取多普勒速度信息,直接利用任意能夠?qū)С霾噬嗥绽粘曅膭?dòng)圖的超聲設(shè)備,基于圖像數(shù)據(jù)獲取速度信息,打破了超聲設(shè)備的局限性;2)采用U-Net 模型自動(dòng)識(shí)別左心室壁輪廓,且采用經(jīng)遷移學(xué)習(xí)重新訓(xùn)練的PWC-Net(convolutional neural Networks for optical flow using Pyramid,Warping,and Cost volume)模型[10]對(duì)所識(shí)別左心室壁進(jìn)行位移估計(jì),減少耗時(shí),提高計(jì)算效率。本文方法作為一種完全基于圖像信息的快速通用型VFM 方法,不需要任何供應(yīng)商的技術(shù)支持和專(zhuān)有軟件,可以進(jìn)一步推進(jìn)VFM在臨床工作流程中的應(yīng)用。

        1 本文方法

        本文提出完全基于二維彩色多普勒超聲心動(dòng)圖像信息的VFM 方法。該方法使用速度標(biāo)尺提取沿聲束方向的徑向速度代替從特定設(shè)備的原始射頻數(shù)據(jù)獲取多普勒速度作為徑向速度;利用兩個(gè)深度學(xué)習(xí)網(wǎng)絡(luò)識(shí)別左心室壁位置并獲取室壁速度,根據(jù)左心室前后壁的運(yùn)動(dòng)模式與心臟流體運(yùn)動(dòng)的關(guān)聯(lián)性,結(jié)合流體力學(xué)中的連續(xù)性方程,以左心室前后壁的切向速度分量作為邊界條件,計(jì)算心臟流場(chǎng)各血液質(zhì)點(diǎn)的切向速度分量,最終合成速度流場(chǎng)。本文主要使用速度矢量圖和平面流線(xiàn)圖對(duì)左心室流體運(yùn)動(dòng)進(jìn)行可視化描述,反映流動(dòng)趨勢(shì)、速度分布特點(diǎn),以及整體血液流動(dòng)結(jié)構(gòu)。

        1.1 徑向速度的提取

        彩色多普勒血流信息包含血流運(yùn)動(dòng)的方向信息,即:朝向探頭方向的血流以紅色顯示,背離探頭方向的血流以藍(lán)色顯示。同時(shí)彩色多普勒血流信息包含了血流速度的變化信息,即:顏色亮度越大,沿超聲聲束方向的速度越大;顏色亮度越小,沿超聲聲束方向的速度越小。

        本文方法中提取彩色血流信息中的多普勒速度是進(jìn)一步研究左心室流體運(yùn)動(dòng)的至關(guān)重要的一個(gè)步驟。本文將彩色血流信息提取的速度作為徑向速度,成為最終形成的流場(chǎng)中的第一個(gè)速度分量。為獲得該速度信息,將彩色多普勒超聲心動(dòng)圖上生成的速度標(biāo)尺作為度量速度值的標(biāo)準(zhǔn)。通過(guò)彩色編碼后每個(gè)血液質(zhì)點(diǎn)對(duì)應(yīng)的彩色血流信息,搜索彩色多普勒?qǐng)D像中的血液質(zhì)點(diǎn)在速度標(biāo)尺中的最佳匹配點(diǎn),本文采用最小二乘法計(jì)算兩者的相似程度:

        其中:c表示血液質(zhì)點(diǎn)與速度標(biāo)尺中點(diǎn)的相似程度;Rb、Gb、Bb分別代表當(dāng)前計(jì)算的血液質(zhì)點(diǎn)的RGB 顏色三分量;Rs、Gs、Bs分別代表速度標(biāo)尺上點(diǎn)的RGB 顏色三分量。這一步主要為了找到當(dāng)前血液質(zhì)點(diǎn)在速度標(biāo)尺中的最優(yōu)匹配點(diǎn)位置,獲取該點(diǎn)在速度標(biāo)尺上的距離參數(shù)。然后,得到徑向速度與距離的分段線(xiàn)性函數(shù)[11](如圖1):

        圖1 徑向速度與距離的分段函數(shù)示意圖Fig.1 Schematic diagram of piecewise function of radial velocity and distance

        式中:d代表最優(yōu)匹配點(diǎn)與速度標(biāo)尺最底端的像素距離;d1代表速度標(biāo)尺藍(lán)色區(qū)域的像素高度;d2代表d1與標(biāo)尺中黑色區(qū)域像素高度和;dmax代表整個(gè)速度標(biāo)尺的像素高度;Vmax代表當(dāng)前速度標(biāo)尺測(cè)量范圍內(nèi)的最大速度。

        通過(guò)速度標(biāo)尺獲得的徑向速度值是不光滑的,無(wú)法滿(mǎn)足后續(xù)的連續(xù)性方程的求解,為確保有比較穩(wěn)定的積分環(huán)境,使用速度標(biāo)尺求出原始徑向速度后,對(duì)徑向速度數(shù)據(jù)進(jìn)行平滑處理。平滑處理的方法包括對(duì)徑向速度數(shù)據(jù)進(jìn)行中值濾波和高斯濾波。中值濾波的處理方法如圖2(a)所示:對(duì)于當(dāng)前所求血液質(zhì)點(diǎn)的徑向速度,以圖2(a)中實(shí)線(xiàn)圈內(nèi)值為例,沿著超聲聲束方向定位該點(diǎn)的前一血液質(zhì)點(diǎn)和后一血液質(zhì)點(diǎn)的徑向速度值;以圖2(a)中兩個(gè)虛線(xiàn)圈內(nèi)值為例,獲取三個(gè)血液質(zhì)點(diǎn)的中值作為當(dāng)前血液質(zhì)點(diǎn)的徑向速度值。高斯濾波的處理方法如圖2(b)所示:通過(guò)二維高斯函數(shù)(見(jiàn)式(3))生成窗口大小為5× 5高斯模板(如圖2(b)中箭頭所示),與當(dāng)前血液質(zhì)點(diǎn)(如圖2(b)中實(shí)心方框)為中心的相同窗口大小范圍的徑向速度進(jìn)行卷積計(jì)算,計(jì)算結(jié)果作為當(dāng)前血液質(zhì)點(diǎn)的徑向速度值,重復(fù)該操作直到完成所有血液質(zhì)點(diǎn)的計(jì)算。中值濾波方法可以消除異常的速度值,高斯濾波方法可以消除小規(guī)模的速度波動(dòng)。

        圖2 平滑方法示意圖Fig.2 Schematic diagram of smoothing method

        其中:G表示高斯模板矩陣;σ表示標(biāo)準(zhǔn)差,本文設(shè)置為8;(x,y)表示模板的二維位置坐標(biāo),由于本文設(shè)置模板大小為5× 5,所以x,y∈[-2,-1,0,1,2]。

        1.2 左心室流體切向速度的計(jì)算

        1.2.1 流體力學(xué)中連續(xù)性方程的應(yīng)用

        流體力學(xué)已成為心血管生物力學(xué)基礎(chǔ)研究領(lǐng)域的重要研究方向,其中連續(xù)性方程是質(zhì)量守恒定理在流體力學(xué)中的應(yīng)用,即流進(jìn)絕對(duì)坐標(biāo)系中任何閉合曲面內(nèi)的質(zhì)量等于從這個(gè)曲面流出的質(zhì)量。本文的研究中,將左心室血液看作一種不可壓縮流體,基于二維平面流動(dòng)假設(shè)下,建立起便于流場(chǎng)分析的極坐標(biāo)系(r,θ),以超聲虛擬探頭O的位置作為極坐標(biāo)的極點(diǎn),以水平于超聲探頭的直線(xiàn)OX作為極軸,角度從極軸的左邊線(xiàn)段出發(fā),取逆時(shí)針?lè)较驗(yàn)檎较?,?guī)定左心室前壁的角度為θ-,左心室后壁的角度為θ+,如圖3所示。

        圖3 極坐標(biāo)系建立示意圖Fig.3 Schematic diagram of polar coordinate system construction

        將彩色多普勒超聲心動(dòng)圖像提供的沿著超聲聲束方向的單一速度作為流體的徑向速度分量,結(jié)合連續(xù)性方程(式(4)),用于計(jì)算流體的切向速度分量。

        其中:r表示距離探頭的半徑大??;θ表示極坐標(biāo)下的角度;z表示垂直貫穿于采集的二維彩色多普勒超聲圖像平面的方向;Vr表示極坐標(biāo)平面的徑向速度;Vθ表示極坐標(biāo)平面的切向速度;Vz表示垂直貫穿平面的速度。本文基于二維平面流動(dòng)假設(shè),使用長(zhǎng)軸切面的超聲心動(dòng)圖作流場(chǎng)分析,忽略垂直貫穿于極坐標(biāo)平面速度帶來(lái)的細(xì)微影響,因此可以將Vz設(shè)置為0,那么式(4)則可以化簡(jiǎn)為如下形式:

        式(5)中將沿聲束方向的徑向速度(Vr)與切向速度(Vθ)聯(lián)系起來(lái),直觀地表示了通過(guò)徑向速度獲得血液質(zhì)點(diǎn)切向速度的計(jì)算方法。該公式在左心室區(qū)域內(nèi)進(jìn)行求解,本文中以左心室的前后壁的切向速度作為式(5)的邊界條件參與切向速度的計(jì)算,得出方程的兩個(gè)解:

        其中,Vθ-和Vθ+分別是左心室前壁與左心室后壁的切向速度,積分部分的計(jì)算方法已經(jīng)在式(5)中給出。第一個(gè)解(式(6))是從左心室前壁起沿第一條積分路徑(圖4 中的左心室前壁與所求血液質(zhì)點(diǎn)虛線(xiàn)相交位置所示)到所求血液質(zhì)點(diǎn)的切向速度積分結(jié)果,第二個(gè)解(式(7))是從左心室后壁起沿第二條積分路徑(圖4 中的左心室后壁與所求血液質(zhì)點(diǎn)虛線(xiàn)相交位置所示)到所求血液質(zhì)點(diǎn)的切向速度積分結(jié)果。

        圖4 連續(xù)性方程的兩個(gè)解Fig.4 Two solutions of continuity equation

        1.2.2 權(quán)重函數(shù)的計(jì)算

        利用權(quán)重函數(shù)w和線(xiàn)性組合連續(xù)性方程的兩個(gè)解,可以得到最終的切向速度,計(jì)算式如下:

        考慮第一個(gè)解(式(6))和第二個(gè)解(式(7))對(duì)于最終血液質(zhì)點(diǎn)的切向速度的貢獻(xiàn)值大小,盡可能地使兩個(gè)解在最終切向速度中所占計(jì)算比例均勻合理,減小結(jié)果的誤差。采用角度距離比值來(lái)定義權(quán)重函數(shù):

        即在極坐標(biāo)系中將血液質(zhì)點(diǎn)與左心室前壁角度差和左心室后壁與左心室前壁角度差的比值作為權(quán)重值。

        1.3 左心室壁切向速度的計(jì)算

        文獻(xiàn)[12-14]研究表明,心室壁心肌運(yùn)動(dòng)狀態(tài)與心臟疾病密切相關(guān),室壁的異常運(yùn)動(dòng)將直接導(dǎo)致心腔內(nèi)血液動(dòng)力學(xué)特征的改變,所以,在研究血液動(dòng)力學(xué)特征時(shí)結(jié)合心室壁的運(yùn)動(dòng)是非常有價(jià)值的。本文使用時(shí)間分辨率重建后的待分析彩色多普勒超聲心動(dòng)圖像和其相鄰的后一幀彩色多普勒超聲心動(dòng)圖像,通過(guò)U-Net 深度學(xué)習(xí)框架自動(dòng)識(shí)別左心室壁輪廓。再將識(shí)別后的左心室壁輪廓,采用重新訓(xùn)練的PWC-Net 模型估計(jì)像素位移。最終,使用超聲設(shè)備設(shè)定的時(shí)間分辨率和心室壁的位置信息獲得左心室壁的物理切向速度,用來(lái)作為連續(xù)性方程的邊界條件。

        1.3.1 左心室壁輪廓自動(dòng)識(shí)別方法

        現(xiàn)有的左心室血流向量成像技術(shù)研究中,絕大多數(shù)使用臨床醫(yī)生手動(dòng)勾勒或者手動(dòng)勾勒結(jié)合自動(dòng)識(shí)別的半自動(dòng)方法識(shí)別左心室壁輪廓[3-4,7-8],所以往往要依靠醫(yī)生的臨床經(jīng)驗(yàn)進(jìn)行判斷,具有較強(qiáng)的主觀性;而且這樣的手動(dòng)方法雖然精確,但是影響了計(jì)算效率,且不同醫(yī)生勾勒的左心室壁位置也會(huì)存在差異,會(huì)對(duì)計(jì)算的流場(chǎng)結(jié)果造成影響。

        本文利用深度學(xué)習(xí)方法實(shí)現(xiàn)左心室輪廓的自動(dòng)分割識(shí)別。對(duì)于醫(yī)學(xué)圖像,U-Net 模型表現(xiàn)出較好的分割性能。UNet模型[15]是一種編解碼器結(jié)構(gòu),編解碼器結(jié)構(gòu)的網(wǎng)絡(luò)非常適用于分割。U-Net遵循特定的方法,其中每個(gè)下采樣步驟和上采樣步驟都要經(jīng)過(guò)兩個(gè)3× 3的卷積層,而每個(gè)下采樣過(guò)程的特征數(shù)量增加一倍,上采樣過(guò)程的特征數(shù)量減少一半。具體來(lái)說(shuō),該模型左半邊為下采樣過(guò)程,右半邊為上采樣過(guò)程。下采樣過(guò)程主要連續(xù)地經(jīng)過(guò)3× 3 的卷積層并使用ReLU 激活層,2 × 2 的最大池化層,不斷提取圖像特征直至最高維。上采樣過(guò)程通過(guò)連續(xù)的2 × 2 的上卷積方式將高維特征向低維映射,并且在映射過(guò)程中連接來(lái)自下采樣過(guò)程中相應(yīng)裁剪的圖像特征,最終得到對(duì)圖像中每一個(gè)像素點(diǎn)的分類(lèi),輸出分割圖像。U-Net模型的訓(xùn)練過(guò)程如圖5所示。

        圖5 U-Net模型的訓(xùn)練與應(yīng)用過(guò)程Fig.5 Training and application process of U-Net model

        本文使用用于二維超聲心動(dòng)圖評(píng)估的最大并且公開(kāi)的全注釋CAMUS 數(shù)據(jù)集[16]。該數(shù)據(jù)集是帶有完整標(biāo)簽注釋的大型數(shù)據(jù)集,其中包含了10 個(gè)文件夾,每個(gè)文件夾中有50 個(gè)患者二維B 模式下的至少一個(gè)心動(dòng)周期的超聲采集序列,舒張末期和收縮末期的兩腔心圖、四腔心圖,以及其對(duì)應(yīng)的由三個(gè)心臟疾病專(zhuān)家在協(xié)定協(xié)議下標(biāo)注的左心室內(nèi)膜(即左心室壁)、左心室外膜、左心房和背景四種標(biāo)簽。

        本文關(guān)注的區(qū)域是左心室壁的區(qū)域,需要識(shí)別左心室壁的位置作為連續(xù)性方程的邊界條件。為了提高訓(xùn)練的速度,本文對(duì)原始數(shù)據(jù)集的標(biāo)簽進(jìn)行預(yù)處理,只對(duì)左心室內(nèi)膜和背景標(biāo)簽作保留,使用U-Net 模型對(duì)左心室進(jìn)行二分類(lèi)的訓(xùn)練。在實(shí)際的模型訓(xùn)練過(guò)程中,以Keras 為深度學(xué)習(xí)框架,使用8個(gè)文件夾下的總共400 個(gè)患者的數(shù)據(jù)作為訓(xùn)練集(共1600 張圖像),1 個(gè)文件夾下的總共50 個(gè)患者的數(shù)據(jù)(共200 張圖像)作為測(cè)試集。根據(jù)本文所提出的方法,針對(duì)左心室壁輪廓的識(shí)別需求,經(jīng)過(guò)大量的實(shí)驗(yàn)訓(xùn)練,最終確定訓(xùn)練參數(shù):為了獲得更高維度的特征圖,在下采樣路徑中使用卷積層和最大池化層,直至特征維數(shù)至16 × 16;使用批量歸一化方法;設(shè)置批尺寸為2;學(xué)習(xí)率設(shè)置為1E-04;使用隨機(jī)梯度下降方法作為最優(yōu)化策略;損失函數(shù)為交叉熵和正則化。最終得到左心室壁的分割模型。

        1.3.2 PWC-Net光流網(wǎng)絡(luò)計(jì)算室壁運(yùn)動(dòng)

        PWC-Net模型的思想與多尺度光流法非常接近。該模型的結(jié)構(gòu)如圖6 所示。在第i個(gè)尺度中,PWC-Net 模型首先通過(guò)輸入的兩張圖像I1和I2提取出圖像特征,可以記作fI1i和fI2i,使用先前尺度(即第i-1 尺度)的上采樣位移估算值對(duì)fI2i進(jìn)行卷繞,作為運(yùn)動(dòng)補(bǔ)償特征,記為。然后,將和fI1i匹配成對(duì)后,使用多層卷積神經(jīng)網(wǎng)絡(luò)(Convolutional Neural Network,CNN)估算位移。在CNN 之后,再增加一個(gè)可選擇使用的上下文網(wǎng)絡(luò)層(額外的卷積層),用于進(jìn)一步對(duì)位移估計(jì)值進(jìn)行優(yōu)化。重復(fù)第i尺度的整個(gè)訓(xùn)練過(guò)程,直到所有尺度均執(zhí)行完該過(guò)程為止,最終輸出位移數(shù)據(jù)。

        圖6 PWC-Net模型結(jié)構(gòu)Fig.6 Structure of PWC-Net model

        由于心臟運(yùn)動(dòng)是一個(gè)自由形變的非剛性運(yùn)動(dòng)過(guò)程,原來(lái)的PWC-Net 光流網(wǎng)絡(luò)模型的訓(xùn)練是通過(guò)自然圖像數(shù)據(jù)完成的,因此本文使用由有限元軟件ANSYS和超聲模擬器Field Ⅱ產(chǎn)生的超聲非剛性運(yùn)動(dòng)數(shù)據(jù)集。首先,通過(guò)有限元軟件構(gòu)建組織模型并施加載荷模擬超聲探頭擠壓組織的過(guò)程,保留產(chǎn)生的位移,并加到壓縮前的組織模型內(nèi)用于獲得壓縮后的圖像;然后,通過(guò)Field Ⅱ平臺(tái)對(duì)壓縮前后的組織模型進(jìn)行超聲模擬成像,獲得壓縮前后的B 型超聲圖像對(duì)。模擬產(chǎn)生的超聲數(shù)據(jù)集包括壓縮前后的4 個(gè)小球仿體數(shù)據(jù)對(duì)和壓縮前后的復(fù)雜乳腺結(jié)構(gòu)數(shù)據(jù)對(duì)兩種類(lèi)型,數(shù)據(jù)集信息如表1所示。

        表1 超聲模擬數(shù)據(jù)集概要信息Tab.1 Summary information of ultrasound simulation datasets

        使用模擬數(shù)據(jù)保留的位移作為標(biāo)簽,模擬的B 型超聲圖像對(duì)作為訓(xùn)練集,對(duì)原有的PWC-Net 模型進(jìn)行遷移學(xué)習(xí)得到新的PWC-Net 模型[17]。模型的訓(xùn)練過(guò)程中,使用原始PWCNet 模型參數(shù)進(jìn)行初始化,4 個(gè)小球仿體數(shù)據(jù)參與模型的第一階段的遷移訓(xùn)練。模型收斂后,使用第一階段獲得的模型初始化PWC-Net,將復(fù)雜乳腺仿體數(shù)據(jù)用于模型的第二階段的遷移訓(xùn)練。模型重要訓(xùn)練參數(shù)信息如表2 所示,其中Slong、Sfine與文獻(xiàn)[18]描述一致:第一階段為粗訓(xùn)練階段,選擇L2 loss作為損失函數(shù),收斂速度快且具有穩(wěn)定解,控制模型的復(fù)雜程度,一定程度避免過(guò)擬合現(xiàn)象;第二階段為精訓(xùn)練階段,選擇L1 loss作為損失函數(shù),具有更好魯棒性,處理模擬數(shù)據(jù)中的異常值[19]。

        表2 PWC-Net模型重要訓(xùn)練參數(shù)信息Tab.2 Important training parameter information of PWC-Net model

        本文將裁剪后的當(dāng)前圖像幀的左心室區(qū)域和相鄰的下一圖像幀的對(duì)應(yīng)區(qū)域作為輸入圖像,利用重新訓(xùn)練后的PWC-Net模型,最終獲得整個(gè)左心室區(qū)域的位移,包括橫向的像素位移值和縱向的像素位移值。然后,使用本文中訓(xùn)練的U-Net模型得到左心室壁位置信息,只保留左心室壁的位移數(shù)據(jù),并且通過(guò)超聲設(shè)備采集時(shí)設(shè)置的圖像幀間的時(shí)間間隔和所求室壁點(diǎn)在極坐標(biāo)系統(tǒng)上的角度位置,可以獲得左心室壁的物理切向速度。

        1.4 本文方法的完整步驟

        本文可視化描述方法的流程如圖7所示,具體步驟如下:

        圖7 左心室血流場(chǎng)可視化方法流程Fig.7 Flow chart of visualization method of left ventricular blood flow field

        輸入 彩色多普勒超聲心動(dòng)圖像;

        輸出 速度矢量圖與流線(xiàn)圖可視化結(jié)果。

        步驟1 合理選擇心動(dòng)周期中的典型階段的圖像,本文選擇兩組實(shí)驗(yàn)圖像,這兩組圖像均處于收縮期階段。裁剪出左心室范圍作為感興趣區(qū)域。

        步驟2 使用一維線(xiàn)性插值的方法補(bǔ)充感興趣區(qū)域的空隙處彩色血流信息(見(jiàn)2.2.2節(jié))。

        步驟3 結(jié)合式(2),利用速度標(biāo)尺提取出感興趣區(qū)域的全部血流質(zhì)點(diǎn)的徑向速度分量并作平滑處理。平滑處理包括中值濾波和高斯濾波。

        步驟4 使用訓(xùn)練的U-Net 模型自動(dòng)識(shí)別左心室壁位置,并重建圖像序列提高時(shí)間分辨率(見(jiàn)2.2.1節(jié))。

        步驟5 利用重新訓(xùn)練的PWC-Net 深度學(xué)習(xí)模型求解左心室壁的切向速度。

        步驟6 通過(guò)式(5)~(9)的計(jì)算,得到左心室區(qū)域內(nèi)所有血液質(zhì)點(diǎn)的切向速度分量。

        步驟7 得到徑向速度分量和切向速度分量,即可合成血流速度場(chǎng)。繪制速度矢量圖用來(lái)觀察血液速度分布的特點(diǎn)和大小,并在速度矢量圖的繪制基礎(chǔ)上繪制流線(xiàn)圖用以觀察左心室腔內(nèi)的流體的血液結(jié)構(gòu)特征。

        2 實(shí)驗(yàn)設(shè)置與評(píng)價(jià)方法

        2.1 數(shù)據(jù)收集

        為了驗(yàn)證本文所提出的聯(lián)合深度學(xué)習(xí)的通用血流向量成像方法,收集不同正常人的多個(gè)心動(dòng)周期的心臟長(zhǎng)軸切面的彩色多普勒超聲心動(dòng)圖。本文的數(shù)據(jù)來(lái)源于四川省人民醫(yī)院心血管超聲及心功能科。為觀察正常狀態(tài)下左心室血液動(dòng)力學(xué)的特征,本文研究一共選取了5 名正常人的超聲心動(dòng)圖數(shù)據(jù)。超聲數(shù)據(jù)采集設(shè)備為Hitachi Aloka Prosound F75 型彩色多普勒超聲診斷儀,在圖像采集的過(guò)程中均采用超聲心動(dòng)圖標(biāo)準(zhǔn)切面進(jìn)行圖像數(shù)據(jù)的采集。最終導(dǎo)出彩色多普勒超聲圖像進(jìn)行離線(xiàn)分析。

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

        2.2.1 重建彩色多普勒超聲心動(dòng)圖圖像序列

        通常情況下,彩色多普勒超聲心動(dòng)圖時(shí)間分辨率不高,圖像幀率大約為30 frame/s,在這種情況下,采集到的圖像相鄰幀之間時(shí)間間隔較長(zhǎng),直接導(dǎo)致相鄰圖像幀之間的相關(guān)性不強(qiáng),最終形成的超聲圖像序列顯然無(wú)法滿(mǎn)足后續(xù)的心室壁運(yùn)動(dòng)跟蹤計(jì)算分析要求。因此,為了獲取心室壁運(yùn)動(dòng)信息,必須提高彩色多普勒超聲心動(dòng)圖的時(shí)間分辨率。

        本文基于時(shí)間超分辨率重建的思想,根據(jù)彩色多普勒超聲心動(dòng)圖中的心電圖時(shí)相位置信息,對(duì)多個(gè)心動(dòng)周期圖像按心電圖時(shí)相先后順序重新進(jìn)行排序,搜索后續(xù)的心動(dòng)周期中時(shí)相在第一心動(dòng)周期中對(duì)應(yīng)相鄰幀之間的圖像,將搜索到的圖像插入對(duì)應(yīng)相鄰幀。直到第一心動(dòng)周期所有相鄰幀均被后續(xù)心動(dòng)周期圖像插入,即可將多個(gè)心動(dòng)周期圖像序列合成為時(shí)間分辨率提升后的一個(gè)心動(dòng)周期圖像序列,從而實(shí)現(xiàn)了提高時(shí)間分辨率。

        本文利用3 個(gè)心動(dòng)周期的彩色多普勒超聲心動(dòng)圖圖像序列信息,通過(guò)重建得到一個(gè)3 倍于原圖像時(shí)間分辨率的彩色多普勒心動(dòng)圖新序列,保證心室壁運(yùn)動(dòng)位移估計(jì)的計(jì)算條件。

        2.2.2 彩色血流信息補(bǔ)償

        臨床診斷中常規(guī)彩色多普勒超聲心動(dòng)圖成像技術(shù)基于脈沖多普勒原理,首先利用多普勒中頻移的信號(hào),獲得取樣面積內(nèi)的血流平均速度信息和速度變化信息,然后進(jìn)行彩色多普勒編碼處理,形成的彩色血流信息疊加到二維B 模式的灰階圖像上,構(gòu)成一幅完整的超聲彩色多普勒?qǐng)D像[20]。在心臟多普勒?qǐng)D像的采集過(guò)程中,心臟左心室某些區(qū)域的血流速度過(guò)低會(huì)導(dǎo)致編碼后的彩色血流信息暗淡甚至無(wú)彩色血流信息。使用速度標(biāo)尺提取徑向速度時(shí),對(duì)這些低速區(qū)域不敏感,導(dǎo)致丟失部分的血流的徑向速度細(xì)節(jié),可能使得最終的流場(chǎng)結(jié)果不連續(xù)。所以,需要在原有的彩色血流信息的基礎(chǔ)上補(bǔ)償原始彩色多普勒心動(dòng)圖像中部分重要的彩色血流信息,從而獲得更加豐富準(zhǔn)確的流場(chǎng)信息。本文使用一維線(xiàn)性插值方法,利用原始圖像數(shù)據(jù)中已有的彩色血流信息進(jìn)行色彩補(bǔ)償,主要補(bǔ)償區(qū)域?yàn)閮煞N不同顏色的彩色血流信息之間的空隙處。一維線(xiàn)性插值方法的計(jì)算式如下:

        其中:y代表所求的部分無(wú)彩色血流信息或彩色血流信息暗淡的血液質(zhì)點(diǎn)的顏色值;y0、y1分別代表當(dāng)前所求血液質(zhì)點(diǎn)的左右最鄰近的含有彩色血流信息的血液質(zhì)點(diǎn)的顏色值;k代表當(dāng)前所求血液質(zhì)點(diǎn)到左邊鄰近含有血流色彩的血液質(zhì)點(diǎn)的距離與當(dāng)前所求血液質(zhì)點(diǎn)到右邊鄰近含有血流色彩的血液質(zhì)點(diǎn)的距離的比值。

        2.3 效果評(píng)價(jià)

        本文實(shí)驗(yàn)選擇心功能分析的舒張期對(duì)本文所提方法進(jìn)行驗(yàn)證。為評(píng)估本文方法的有效性,主要使用兩種方法對(duì)本文實(shí)驗(yàn)結(jié)果進(jìn)行評(píng)價(jià)。首先,Hitachi Aloka 提供的VFM 系統(tǒng)工作站作為業(yè)界首個(gè)血流向量成像系統(tǒng),在心功能臨床分析中被廣泛關(guān)注并使用。對(duì)比Hitachi Aloka VFM 系統(tǒng)工作站離線(xiàn)分析的左心室血流速度矢量圖和平面流線(xiàn)圖結(jié)果與本文方法得到的結(jié)果的一致性,并將其作為直接評(píng)估參考。通過(guò)觀察速度矢量圖的速度分布特征、血液流動(dòng)走向,以及平面流線(xiàn)圖的渦流大小、位置等判斷兩者的一致性程度。其次,將由四川省人民醫(yī)院的兩名具有多年工作經(jīng)驗(yàn)的超聲醫(yī)生的主觀評(píng)價(jià)作為評(píng)估參考,根據(jù)左心室血液動(dòng)力學(xué)的特征,判斷本文實(shí)驗(yàn)結(jié)果的合理性和有效性。

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

        3.1 原始圖像彩色信息補(bǔ)充結(jié)果

        圖8 為一個(gè)正常心臟的超聲心動(dòng)圖,處于心臟臨床分析的典型時(shí)期——舒張期,每組測(cè)試圖像的右上角給出了對(duì)應(yīng)的彩色血流信息的補(bǔ)充結(jié)果。從兩組測(cè)試圖像結(jié)果可以看出,根據(jù)一維線(xiàn)性插值的方法,無(wú)顏色區(qū)域補(bǔ)償?shù)牟噬餍畔⒅?,越靠近有顏色區(qū)域的點(diǎn)顏色越亮,速度相較于遠(yuǎn)離有顏色區(qū)域的點(diǎn)大,能夠合理有效地補(bǔ)償出不同流向的血流之間空隙處的部分血流信息。這個(gè)補(bǔ)償結(jié)果為后續(xù)的左心室腔內(nèi)的血液流場(chǎng)分析提供了更加完整的信息。

        圖8 測(cè)試彩色多普勒?qǐng)D像與色彩補(bǔ)償結(jié)果Fig.8 Test color-Doppler images and color compensation results

        3.2 左心室輪廓自動(dòng)識(shí)別結(jié)果

        本文使用公開(kāi)數(shù)據(jù)集中50 個(gè)患者組成的測(cè)試集對(duì)訓(xùn)練后的U-Net 模型進(jìn)行性能評(píng)估。使用每一個(gè)患者兩腔心圖和四腔心圖的左心室壁的標(biāo)簽(總共200 張)與訓(xùn)練的U-Net 模型預(yù)測(cè)的左心室壁位置基于度量指標(biāo)計(jì)算平均值。通常針對(duì)醫(yī)學(xué)圖像分割模型,常用的度量指標(biāo)包括戴斯相似性系數(shù)DICE、體積重疊誤差(Volumetric Overlap Error,VOE)、相對(duì)體積差異(Relative Volume Difference,RVD)、精確性(Precision)和召回率(Recall)等。其中:DICE 表示兩標(biāo)簽結(jié)果相交的面積占總面積的比值大??;VOE 與DICE 的原理類(lèi)似,采用兩標(biāo)簽結(jié)果相減的面積占總面積的比值來(lái)表示錯(cuò)誤率;RVD 表示兩標(biāo)簽之間的體積差異;Precision 表示判斷為真的正例占所有判斷為真的樣例的比重;Recall 表示判斷為真的正例占總正例的比率。以上評(píng)估值均在0~1。評(píng)估模型性能時(shí),DICE、Precision、Recall的值越接近1,則模型性能越好;而VOE、RVD的值越接近0,模型的性能越好。

        對(duì)200張圖像的5個(gè)評(píng)估指標(biāo)值進(jìn)行平均值的計(jì)算,模型的度量指標(biāo)結(jié)果如表3 所示。從表3 可以看出,DICE、Precision、Recall 值均處于0.9 以上,在較大值范圍內(nèi);VOE、RVD 值均處于0.04 以下,在較小值范圍內(nèi),表明本文訓(xùn)練的U-Net模型在左心室壁輪廓分割中性能較好。

        表3 U-Net模型評(píng)估結(jié)果Tab.3 U-Net model evaluation results

        使用訓(xùn)練好的U-Net 模型對(duì)圖8 中的彩色多普勒超聲心動(dòng)圖測(cè)試組進(jìn)行左心室壁的分割預(yù)測(cè),圖9(a)對(duì)應(yīng)圖8(a)中的左心室壁輪廓,圖9(b)對(duì)應(yīng)圖8(b)中的左心室壁輪廓。兩組圖像數(shù)據(jù)顯示,左心室壁輪廓位置合理,分割的輪廓邊緣連續(xù)光滑,并且通過(guò)了四川省醫(yī)學(xué)科學(xué)院·四川省人民醫(yī)院心血管超聲及心功能科的醫(yī)生對(duì)分割結(jié)果的主觀認(rèn)定??傮w來(lái)說(shuō),本文訓(xùn)練的自動(dòng)分割模型能夠較好地對(duì)左心室壁輪廓進(jìn)行識(shí)別。

        圖9 左心室輪廓自動(dòng)識(shí)別結(jié)果Fig.9 Automatic recognition results of left ventricular contour

        3.3 室壁運(yùn)動(dòng)追蹤結(jié)果

        圖10 為原始的PWC-Net 和重新訓(xùn)練的PWC-Net 模型對(duì)彩色多普勒超聲心動(dòng)圖的左心室區(qū)域的心室壁運(yùn)動(dòng)追蹤結(jié)果對(duì)比。通過(guò)圖10(a)、(b)與圖10(c)、(d)形成橫向位移與縱向位移的對(duì)比可以看出,直接使用原始PWC-Net 模型的測(cè)試數(shù)據(jù)出現(xiàn)嚴(yán)重分層曲線(xiàn)與噪點(diǎn),如圖10(a)、(c)中小方框所示,重新訓(xùn)練的PWC-Net 模型的測(cè)試結(jié)果則表現(xiàn)得更加穩(wěn)定平滑,預(yù)測(cè)結(jié)果具有連續(xù)性,表明重新訓(xùn)練的PWC-Net模型能夠?qū)ψ笮氖冶诘倪\(yùn)動(dòng)進(jìn)行準(zhǔn)確估計(jì)。

        圖10(e)、(f)展示了兩種模型分別由圖10(a)、(c)和圖10(b)、(d)所示位移合成的位移矢量圖。由可視化的運(yùn)動(dòng)矢量圖結(jié)果看出,圖10(f)中左心室前壁與后壁均呈現(xiàn)向外擴(kuò)張的運(yùn)動(dòng)趨勢(shì),符合舒張期心室壁的動(dòng)力學(xué)特征。

        圖10 PWC-Net模型的左心室室壁位移估計(jì)結(jié)果Fig.10 Left ventricular wall displacement estimation results obtained by PWC-Net model

        綜上,位移估計(jì)結(jié)果驗(yàn)證了重新訓(xùn)練的PWC-Net 模型對(duì)于心臟的非剛性運(yùn)動(dòng)檢測(cè)的有效性,并且相較于傳統(tǒng)依賴(lài)迭代求解的非剛性配準(zhǔn)運(yùn)動(dòng)估計(jì)[21],重新訓(xùn)練的PWC-Net 模型計(jì)算幀率可達(dá)到45 frame/s,滿(mǎn)足實(shí)時(shí)成像的需求。

        3.4 血流速度矢量圖

        根據(jù)兩組速度矢量圖結(jié)果,可以有效地對(duì)左心室腔內(nèi)的血液運(yùn)動(dòng)趨勢(shì)以及血液流動(dòng)速度大小的分布進(jìn)行觀察。圖11(a)、(b)對(duì)應(yīng)圖8(a)中的速度矢量圖結(jié)果,圖11(c)、(d)對(duì)應(yīng)圖8(b)中的速度矢量圖結(jié)果。從圖11(a)的速度矢量圖可以看出,左心室腔內(nèi)的血流速度在流入道一側(cè)較高,血液流向心尖位置碰及心尖段室壁,血液運(yùn)動(dòng)方向轉(zhuǎn)向流出道位置。從圖11(c)的速度矢量圖可以看出,該組彩色多普勒心動(dòng)圖對(duì)應(yīng)心動(dòng)周期的階段,稍晚于圖11(a)的舒張期階段,且相較于圖11(a)的速度矢量,圖11(c)的速度矢量經(jīng)心尖位置后完成大量的血流轉(zhuǎn)向,整體的血液流動(dòng)速度分布更加遠(yuǎn)離心尖位置,朝著左心室基底段方向靠近,整體速度均勻,并在流入道一側(cè)速度相較于圖11(a)稍降低。

        圖11 速度矢量圖結(jié)果對(duì)比圖11 Velocity vector map result comparison

        將本文提出的血流向量成像方法與Aloka VFM 工作站的后處理分析方法進(jìn)行對(duì)比,結(jié)果顯示兩種方法的結(jié)果具有高度的一致性,進(jìn)行實(shí)驗(yàn)的兩組測(cè)試圖像中,血流區(qū)域速度流場(chǎng)呈現(xiàn)出合理的血流運(yùn)動(dòng)趨勢(shì)和血流速度大小分布,能夠有效地觀察到左心室速度分布情況與血流速度的運(yùn)動(dòng)走向。

        3.5 平面流線(xiàn)圖

        在繪制速度矢量圖的基礎(chǔ)上繪制平面流線(xiàn)圖,平面流線(xiàn)圖的繪制結(jié)果有助于更加直觀地觀察左心室血液流動(dòng)的運(yùn)動(dòng)結(jié)構(gòu)模式。圖12(a)、(b)和圖12(c)、(d)分別為對(duì)應(yīng)圖8(a)、(b)的平面流線(xiàn)圖結(jié)果。通過(guò)流線(xiàn)圖的可視化可以看出,本文的方法能夠清晰地觀察到左心室中的渦流結(jié)構(gòu)。渦流結(jié)構(gòu)是血液動(dòng)力學(xué)的重要指標(biāo)之一,可以為心臟疾病的輔助診斷提供更加直接的證據(jù)。

        從圖12(a)中可以觀察到,在這個(gè)心動(dòng)周期階段成像的血液的渦流狀態(tài),在靠近流入道一側(cè)的室壁位置出現(xiàn)一個(gè)較小的渦流,并在靠近流出道的一側(cè)形成一個(gè)較大的渦流結(jié)構(gòu)(箭頭所指處)。這種動(dòng)力學(xué)特征符合人體心臟運(yùn)動(dòng)的特征。這個(gè)階段,血流碰撞左心室心尖段的室壁轉(zhuǎn)向,所以形成一個(gè)較大的渦流狀態(tài)。

        從圖12(c)中可以觀察到,肉眼可見(jiàn)一個(gè)大型的渦流充盈整個(gè)左心室腔(箭頭所指處),同時(shí)在靠近流入道一側(cè)的室壁位置出現(xiàn)一個(gè)較小的渦流。從圖12(a)到圖12(c),渦流未消失,并且保持旋轉(zhuǎn),左心室血液借助渦流的動(dòng)能進(jìn)行傳遞與轉(zhuǎn)向,為收縮期的到來(lái)提供足夠的動(dòng)能支持。

        圖12 平面流線(xiàn)圖結(jié)果對(duì)比Fig.12 Plane streamline map result comparison

        將圖12(a)與圖12(b)進(jìn)行比較可以看出,血液渦流運(yùn)動(dòng)狀態(tài)相似,即本文方法和Aloka VFM 工作站提供的方法均能夠檢測(cè)到靠近流出道一側(cè)的較大渦流與靠近流入道一側(cè)的較小渦流。將圖12(c)與圖12(d)進(jìn)行比較可以看出,在彩色區(qū)域都可以觀察到充盈左心室的大型渦流與流入道一側(cè)的小渦流,有相似的血流運(yùn)動(dòng)結(jié)構(gòu)和相似的渦流結(jié)構(gòu)。本文方法得到的流線(xiàn)結(jié)果與Aloka VFM 工作站得到的結(jié)果高度一致,并且選擇的實(shí)驗(yàn)圖像結(jié)果表明,本文方法得到的速度矢量圖和平面流線(xiàn)圖結(jié)果能夠有效可視化主要血流區(qū)域的血流動(dòng)力學(xué)模式的特征。

        4 結(jié)語(yǔ)

        本文提出了一種心臟功能輔助評(píng)價(jià)的VFM 方法。本文方法完全基于彩色多普勒超聲心動(dòng)圖像信息,結(jié)合心臟室壁運(yùn)動(dòng)信息和流體力學(xué)連續(xù)性方程,實(shí)現(xiàn)心臟流體運(yùn)動(dòng)信息的提取和分析計(jì)算,從而實(shí)現(xiàn)了左心室流體運(yùn)動(dòng)的可視化描述。實(shí)驗(yàn)結(jié)果表明,本文所提出的基于彩色多普勒超聲心動(dòng)圖的分析通用方法,與Aloka VFM 工作站的結(jié)果呈現(xiàn)出一致性,擺脫了超聲設(shè)備的限制,有效實(shí)現(xiàn)了左心室流體運(yùn)動(dòng)狀態(tài)的可視化描述。該方法能夠更加廣泛地應(yīng)用于臨床診斷中,有較好的應(yīng)用前景。但由于左心室血液流體是一個(gè)三維流動(dòng)狀態(tài),而本文方法是基于二維平面假設(shè)的,忽略了垂直平面的血流運(yùn)動(dòng)影響,會(huì)存在一定的差異;同時(shí),本文方法利用速度標(biāo)尺提取速度信息,與原始射頻數(shù)據(jù)獲取速度信息的方法不同,導(dǎo)致了結(jié)果具有細(xì)微差異,未來(lái)可以在這些方面做進(jìn)一步的研究。

        猜你喜歡
        徑向速度心動(dòng)圖左心室
        超聲心動(dòng)圖診斷Fabry病1例
        王新房:中國(guó)超聲心動(dòng)圖之父
        心電向量圖診斷高血壓病左心室異常的臨床應(yīng)用
        早孕期超聲心動(dòng)圖在胎兒嚴(yán)重先心病中的應(yīng)用
        超聲心動(dòng)圖診斷Loffler心內(nèi)膜炎1例
        非圓形光纖研究進(jìn)展
        航空兵器(2017年6期)2018-01-24 15:00:10
        臺(tái)風(fēng)威馬遜造成云南文山州強(qiáng)降水天氣雷達(dá)回波分析
        初診狼瘡腎炎患者左心室肥厚的相關(guān)因素
        距離頻率ML方法無(wú)模糊估計(jì)動(dòng)目標(biāo)徑向速度
        基于改進(jìn)多尺度ASM和非剛性配準(zhǔn)的4D-CT左心室分割
        国产老熟女狂叫对白| 日韩在线精品视频一区| 国产精品福利一区二区| 中国熟妇人妻xxxxx| 亚洲人妻无缓冲av不卡| 玩弄丝袜美腿超短裙校花| 亚洲av不卡免费在线| 国产一极内射視颍一| 国产亚洲欧美在线| 久久少妇呻吟视频久久久| 97人妻精品一区二区三区男同| 中文字幕人妻被公上司喝醉| 欧美精品在线一区| 国内精品极品久久免费看| 伊人久久这里只有精品| 黄瓜视频在线观看| 337p日本欧洲亚洲大胆色噜噜| 国产精品自拍视频免费看| 欧美亚洲精品suv| 亚洲av无码一区二区二三区| 日韩亚洲欧美精品| 国产亚洲中文字幕一区| 无码国内精品久久人妻| 国产乱子伦露脸在线| 日本骚色老妇视频网站| 女同同性av观看免费| 色八区人妻在线视频免费| 国产一级淫片免费播放电影| 一区二区三区在线观看人妖| 含紧一点h边做边走动免费视频 | 久久亚洲av成人无码国产最大| 中文字幕人妻熟女人妻洋洋| 国产视频最新| 午夜少妇高潮在线观看视频 | 另类内射国产在线| 日韩美女高潮流白浆视频在线观看| 日本中文字幕乱码中文乱码| 成人无码网www在线观看| 国产在线一区观看| 视频在线播放观看免费| 极品尤物一区二区三区|