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

        ?

        基于EGSnrc程序的錐形束CT模擬仿真研究

        2019-10-18 06:23:32孟慧鵬艾克文王克強張艷龍丁紅軍段敬豪
        醫(yī)療衛(wèi)生裝備 2019年10期
        關(guān)鍵詞:離軸錐形子程序

        付 娟,孟慧鵬,*,艾克文,顧 欣*,王克強,張艷龍,丁紅軍,段敬豪,孫 倩

        (1.武警特色醫(yī)學(xué)中心,天津 300162;2.天津大學(xué)精密儀器與光電子工程學(xué)院,天津 300072)

        0 引言

        機載式錐形束CT已經(jīng)廣泛用于調(diào)強放射治療(intensity modulated radiotherapy,IMRT)的圖像引導(dǎo),成為位置驗證及跟蹤腫瘤變化等方面的有力工具[1-3]。隨著放療技術(shù)向自適應(yīng)放射治療(adaptive radiotherapy,ART)發(fā)展,越來越多的學(xué)者開始關(guān)注基于錐形束CT進(jìn)行形變配準(zhǔn)、劑量計算等新的需求[4-6]。然而為方便此類研究的順利開展,首先需要解決的問題是如何對現(xiàn)有錐形束CT進(jìn)行有效的模擬仿真。因此本研究擬通過Monte Carlo模擬程序EGSnrc對錐形束CT進(jìn)行模塊化建模,并對所建立的系統(tǒng)進(jìn)行驗證,進(jìn)而探討所建立的仿真系統(tǒng)是否能夠準(zhǔn)確模擬錐形束CT的主要技術(shù)指標(biāo),及其是否具有進(jìn)一步研究其他錐形束CT相關(guān)問題的潛能。

        1 材料和方法

        1.1 儀器設(shè)備

        實驗設(shè)備為Varian IX直線加速器(美國,瓦里安公司);三方能譜文件由開源程序Spektr(Version 3.0)生成;測量時采用型號為MT-VW的紅固體水,密度為1.03 g/cm3;百分深度劑量測量使用德國PTW公司的MP3-M大水箱,離軸劑量測量使用該公司的OCTAVIUS Detector 729二維電離室矩陣;模擬仿真程序的版本為EGSnrc v2016(2016年1月發(fā)布),在Ubuntu 16.04的64位系統(tǒng)下運行[7]。

        1.2 方法

        1.2.1 主要部件參數(shù)

        擬模擬仿真的錐形束CT系統(tǒng)集成于Varian IX直線加速器上,主要組成構(gòu)件為X射線源和非晶硅平板探測器。X射線源球管為Varian公司生產(chǎn),型號為G242,峰值電壓和陽極最大熱容量分別為125 kVp和445 kJ,陽極材質(zhì)采用鎢合金,靶光束角為14°,初級濾波為2 mm厚的鋁合金;平板探測器同樣為Varian公司生產(chǎn),型號為4030CB,有效探測面積為39.7 cm×29.8 cm,矩陣尺寸為 1 024×768,像素尺寸為0.388 mm×0.388 mm,探測器可沿探測平面橫向向外偏移 16 cm,最大掃描視野(field of view,F(xiàn)OV)可達(dá)到46.6 cm[8]。

        1.2.2 錐形束CT的模擬仿真

        本研究通過EGSnrc的BEAMnrc子程序?qū)㈠F形束CT系統(tǒng)分為X射線源、射線源組件、探測器、Monte Carlo模擬用模體共計4個模塊建模,所建立的系統(tǒng)示意圖如圖1所示。根據(jù)掃描時采用的領(lǐng)結(jié)濾波不同,分為全扇束模式和半扇束模式。實驗使用的錐形束CT系統(tǒng)掃描協(xié)議共有3種掃描電壓,而125 kVp對應(yīng)2種領(lǐng)結(jié)濾波模式,因此共建立4種源模型:100 kVp全領(lǐng)結(jié)濾波、110 kVp半領(lǐng)結(jié)濾波、125 kVp全領(lǐng)結(jié)濾波和半領(lǐng)結(jié)濾波。

        圖1 EGSnrc模擬的CBCT系統(tǒng)示意圖

        EGSnrc的BEAMnrc子程序可提供多種建模組件模塊,通過不同組件的組合可實現(xiàn)錐形束CT的模擬仿真。可用于模擬錐形束CT系統(tǒng)的組件有XTUBE、CONSTAK、SLABS、BLOCK、JAWS和PYRAMIDS等,模擬仿真時,由SLABS組件構(gòu)建初級濾波、射束硬化濾波、機載影像系統(tǒng)窗和探測器4種結(jié)構(gòu);領(lǐng)結(jié)濾波既可以用PYRAMIDS組件構(gòu)建,也可以用JAWS組件構(gòu)建,本研究中采用PYRAMIDS組件。每個組件的最終輸入?yún)?shù)如下:

        (1)X射線源。

        ①XTUBE組件構(gòu)建靶:由Varian公司提供的球管參數(shù)文檔可知,陽極靶物質(zhì)為鎢合金,其中主要成分為金屬鎢(占比95%)、金屬銠(占比5%),等效物理密度為18.68 g/cm3。此外,靶光束角為14°,入射電子的焦距為0.4~0.8 mm,本實驗?zāi)M定義的焦斑直徑為0.5mm。②CONSTAK組件構(gòu)建球管出窗:玻璃窗的密度設(shè)置為2.23 g/cm3,厚度設(shè)置為0.9 mm;窗口與靶之間的距離設(shè)置為2.8cm,出窗與靶之間設(shè)置為真空。

        (2)射線源組件。

        ①SLABS組件構(gòu)建初級濾波:初級濾波為2 mm厚的鋁合金,與靶的距離為5.72 cm。其中鋁合金的等效物理密度為3.23 g/cm3,包含96.5%的鋁、1.3%的硅、1.2%的鎂和1%的錳。②BLOCK組件構(gòu)建初級準(zhǔn)直器:初級準(zhǔn)直器的材質(zhì)為2 cm厚的鉛,其上端開口(靠近源方向)和下端開口(遠(yuǎn)離源方向)分別為3.2cm×3.2cm和4.2cm×4.2cm,與靶的距離為6.52cm。③JAWS組件構(gòu)建擋片:擋片的材質(zhì)為0.2 mm厚的鉛,同準(zhǔn)直器的設(shè)置類似,上下端開口分別設(shè)置為4.72 cm×4.72 cm和4.86 cm×4.86 cm,與靶的距離為9.59 cm。④SLABS組件構(gòu)建射束硬化濾波:射束硬化濾波為0.5 mm厚的鈦,物理密度為4.54 g/cm3,與靶的距離設(shè)置為11.19 cm。⑤SLABS組件構(gòu)建機載影像系統(tǒng)窗:機載影像系統(tǒng)窗為1 mm厚的透光有機物聚碳酸酯,與靶的距離設(shè)置為14.9 cm。⑥PYRAMIDS組件構(gòu)建領(lǐng)結(jié)濾波:構(gòu)建全扇束領(lǐng)結(jié)濾波和半扇束領(lǐng)結(jié)濾波,領(lǐng)結(jié)濾波的材質(zhì)為鋁合金,與初級濾波的鋁合金成分相同。建模時由6層鋁層組成(每層中心開孔不同),從上到下的鋁層厚度分別為 1.5、1.5、7、13、3 和 2 mm,與靶的距離為 15 cm。

        (3)探測器。

        實驗采用的是Varian 4030CB平板探測器,BEAMnrc子程序建模時主要考慮兩個部分:一個是裝配于探測器陣列上方的10∶1的抗散射柵,另一個是平板探測器銅層和由無定形硅構(gòu)成的平板薄膜晶體管陣列(閃爍體層)。擺位采集到的錐形束CT圖像原始投影信號實際上就是探測器陣列對到達(dá)探測器粒子的響應(yīng),建模時X射線的響應(yīng)可通過提取該層的能量沉積實現(xiàn)。對于千伏級別的X射線,抗散射柵自身和平板探測器銅層與粒子的相互作用不可忽略,因此將探測器建模成1 mm厚的銅層,采用SLABS組件構(gòu)建,物理密度設(shè)置為8.9 g/cm3,與靶的距離為150cm,設(shè)置平板探測器的大小為39.7cm×29.8cm,矩陣尺寸為1024×768,像素尺寸為0.388mm×0.388mm。

        (4)Monte Carlo模擬用模體。

        為了驗證BEAMnrc子程序建立系統(tǒng)的準(zhǔn)確性及后續(xù)實驗,建立的固體水模尺寸為30 cm×30 cm×30 cm,同時建立相同尺寸的標(biāo)準(zhǔn)水模體。

        1.2.3 仿真系統(tǒng)的驗證

        為了驗證建立系統(tǒng)的幾何和光束特性是否正確,通過X射線能譜、百分深度劑量、離軸劑量等進(jìn)行分析比較。Monte Carlo模擬上述參數(shù)的結(jié)果分兩步生成:第一,在BEAMnrc子程序中相應(yīng)平面生成相空間文件;第二,由相應(yīng)子程序分析相空間文件得到對應(yīng)結(jié)果,其中BEAMDP子程序可獲得X射線能譜,DOSXYZnrc子程序可獲得百分深度劑量和離軸劑量。比對數(shù)據(jù)的生成如下:由開源程序Spektr生成X射線能譜,由實測數(shù)據(jù)獲取百分深度劑量和離軸劑量。

        (1)X射線能譜。

        用于比對的能譜文件由開源程序Spektr(Version 3.0)生成,分別生成100、110和125 kVp的能譜文件。具體設(shè)置時,在Spektr程序的UI界面中選擇Tube_09,并在tubeSettings函數(shù)中將Tube_09的值自定義為[13 2;74 0],在“ADDED FILTERATION”模塊中,增加Ti(鈦,Z=22),“Thickness”設(shè)置為0.5 mm,以確保生成的能譜與本實驗設(shè)備的參數(shù)盡可能一致。

        Monte Carlo模擬X射線能譜時,對應(yīng)錐形束CT系統(tǒng)掃描協(xié)議的3種X射線管電壓,實驗驗證進(jìn)行了3次主要的Monte Carlo模擬。每次運行均得到對應(yīng)的相空間文件,該文件生成于模擬仿真中的一個虛擬平面,可存儲所有通過該平面的粒子的能量、位置、方向等相關(guān)數(shù)據(jù),由BEAMDP子程序分析后得到X射線能譜。由于Spektr程序只能生成純鎢靶X射線能譜,且只能生成距離源100 cm處的能譜文件,為了便于比較,在生成分析X射線能譜的相空間文件模擬中,本實驗將所建立系統(tǒng)中的靶材料替換為純鎢(其余模擬均為鎢合金),相空間文件記錄平面設(shè)置為距離靶100 cm處。X射線能譜的Monte Carlo模擬采用的粒子數(shù)為1×109,在距離靶100 cm處模擬的射野大小為50 cm×50 cm,模擬時未加載領(lǐng)結(jié)濾波組件。

        (2)百分深度劑量和離軸劑量。

        實測數(shù)據(jù)均需錐形束CT射線源在0°處,因此選擇設(shè)備的服務(wù)模式,該模式可自定義相關(guān)掃描參數(shù),如管電壓、管電流、掃描時間、有無領(lǐng)結(jié)濾波、擋片開口大小等,所有數(shù)據(jù)均在透視模式下源皮距(source skin distance,SSD)100 cm處測量得到。為驗證所建立的仿真系統(tǒng)的準(zhǔn)確性,分別實測了100kVp全領(lǐng)結(jié)濾波、110 kVp半領(lǐng)結(jié)濾波、125 kVp全領(lǐng)結(jié)濾波和半領(lǐng)結(jié)濾波的百分深度劑量和離軸劑量。百分深度劑量測量時設(shè)置擋片野大小為26 cm×26 cm,在射束穩(wěn)定狀態(tài)下,采集射束中心軸0~20 cm范圍內(nèi)的劑量;離軸劑量測量時將0.5 cm厚的固體水放置于電離室矩陣之上,采集固體水下0.75 cm平面內(nèi)等中心位置橫向和縱向的劑量,射野大小為20 cm×20 cm。

        Monte Carlo模擬百分深度劑量和離軸劑量時,對應(yīng)錐形束CT系統(tǒng)實測數(shù)據(jù)的4種電壓-濾波模式進(jìn)行4次主要的Monte Carlo模擬,共記錄4組相空間文件,記錄位置為源到等中心距離(source isocentre distance,SID)100 cm 處,Monte Carlo 模擬使用的粒子數(shù)為1×109。記錄的相空間文件使用DOSXYZnrc子程序分析計算,使用的模體大小為30 cm×30 cm×30 cm,計算體素大小為 0.5 cm×0.5 cm×0.5 cm,粒子數(shù)為5×108,光子分裂數(shù)設(shè)置為100,對所有感興趣區(qū)域的體素計算中統(tǒng)計不確定性控制在2%以內(nèi)。DOSXYZnrc子程序的輸出文件格式為*.3ddose,該文件可由STATDOSE子程序分析。提取出射束中心軸0~20 cm范圍內(nèi)的劑量和模體深度為1.25 cm平面內(nèi)等中心位置橫向和縱向的劑量,繪制生成百分深度劑量和橫向、縱向離軸劑量曲線與實測結(jié)果進(jìn)行比較。

        2 結(jié)果

        在SID 100 cm平面處,3種管電壓100、110和125 kVp的Monte Carlo模擬結(jié)果與Spektr程序生成的對應(yīng)比較結(jié)果如圖2所示,所有能譜信息均與峰值信息做了歸一化處理。結(jié)果表明,2種能譜具有良好的一致性,差異為0.6%~3.0%,平均差異小于2.0%。由圖2還可知,3種管電壓的最大差異區(qū)域均出現(xiàn)在30~55 keV范圍內(nèi),其中125 kVp管電壓的36 keV能量處差異達(dá)到了3.0%。

        圖2 仿真系統(tǒng)生成的能譜信息與Spektr程序生成的能譜信息的比較

        Monte Carlo模擬的百分深度劑量和離軸劑量結(jié)果與實測的結(jié)果一致性良好,如圖3~5所示,所有的數(shù)據(jù)均對最大值做了歸一化處理。圖3為百分深度劑量分布的比較,4種電壓-濾波模式的射束中心軸劑量差異均較小,劑量差異平均值在±2.0%以內(nèi)。其中100 kVp全領(lǐng)結(jié)濾波差異最大,且表現(xiàn)為全領(lǐng)結(jié)濾波模式[如 3 圖(a)和(c)所示)]的差異大于半領(lǐng)結(jié)濾波模式[如圖 3(b)和(d)所示],且水下 5 cm厚度以內(nèi)的Monte Carlo模擬值與實測值的差異大于5 cm以外的。差異最大處在100 kVp全領(lǐng)結(jié)濾波的水下5 cm范圍內(nèi),達(dá)到了2.7%;4種模式水下5 cm以外差異值均逐漸減小,均小于1.0%。圖4為橫向離軸劑量分布,用于驗證系統(tǒng)加載全領(lǐng)結(jié)濾波和半領(lǐng)結(jié)濾波過濾器的幾何設(shè)計,4種電壓-濾波模式橫向離軸劑量的差異平均值在±1.1%以內(nèi)。由圖4可知,110 kVp半領(lǐng)結(jié)濾波模式Monte Carlo模擬值與實測值在4種模式中差異最大,最大處達(dá)到了4.0%,該模式中的差異平均值為1.3%。圖5為縱向離軸劑量分布比較,4種電壓-濾波模式縱向離軸劑量的差異平均值在±0.8%以內(nèi),最大差異為2.7%。由圖5可知,4種模式的曲線有類似的形狀,其中125 kVp全領(lǐng)結(jié)濾波和半領(lǐng)結(jié)濾波模式的曲線肉眼無法區(qū)分其差異。

        圖3 仿真系統(tǒng)中心軸百分深度劑量曲線的Monte Carlo模擬數(shù)據(jù)與實測數(shù)據(jù)比較

        圖4 仿真系統(tǒng)等中心位置橫向離軸劑量曲線的Monte Carlo模擬數(shù)據(jù)與二維電離室矩陣實測數(shù)據(jù)比較

        3 討論

        ART技術(shù)最終向劑量引導(dǎo)放射治療(dose-guided radiotherapy,DGRT)[9]發(fā)展,實現(xiàn)實際患者體內(nèi)劑量沉積的實時追蹤[4]。通過治療時采集的患者的擺位錐形束CT圖像可間接實現(xiàn)患者體內(nèi)累積劑量的追蹤,是ART技術(shù)當(dāng)前較易實現(xiàn)的一種方法。然而,基于錐形束CT進(jìn)行形變配準(zhǔn)、劑量計算的研究[5-6]必須基于一套準(zhǔn)確高效的模擬仿真系統(tǒng)。本研究在充分了解Varian OBI系統(tǒng)(錐形束CT系統(tǒng))的構(gòu)造及特點的基礎(chǔ)上,基于Monte Carlo模擬軟件的EGSnrc程序?qū)崿F(xiàn)了錐形束CT系統(tǒng)的模擬仿真。所建立系統(tǒng)計算得到的百分深度劑量、離軸劑量等物理參數(shù)與實測值誤差平均值均在±2.0%以內(nèi),與文獻(xiàn)[8,10-11]中其他學(xué)者的建模方法相比,各種幾何結(jié)構(gòu)及物質(zhì)組成更接近真實情況,如采用鎢合金和鋁合金代替其他學(xué)者采用純鎢和純鋁的模擬、考慮平板探測器的1 mm銅層的模擬等。此外,通過改變各組件有限的參數(shù),所建立的系統(tǒng)可對其他主流錐形束CT系統(tǒng)進(jìn)行模擬,具有良好的通用性,可用于模擬研究其他錐形束CT系統(tǒng)及錐形束CT系統(tǒng)的其他問題。

        值的一提的是,實驗中模擬生成的能譜未能與實際測量的能譜數(shù)據(jù)進(jìn)行比對,主要原因是能譜的測量復(fù)雜,實際測量較困難,且測量精度受多種因素制約,因此本實驗采用很多研究使用的成熟的第三方能譜生成軟件來生成比對能譜信息。能譜生成程序采用的是Spektr最新版本,其基于三次樣條插值算法的鎢陽極光譜模型(tungsten anode spectral model using interpolating cubic splines,TASMICS)的算法,相對于之前采用多項式插值的鎢陽極光譜模型(tungsten anode spectral model using interpolating polynomials,TASMIP)算法的老版本進(jìn)行了改進(jìn),3.0版本允許用戶自定義球管參數(shù)、不同厚度材料的過濾器及對生成的能譜文件優(yōu)化,可生成20~150 keV能量范圍內(nèi)的X射線能譜,其模擬出的能譜信息與實際數(shù)據(jù)的差異在4.0%以內(nèi),平均差異在2.7%以內(nèi)[10]。本中心Ming等[11-12]的前期工作表明,10%的能譜差異可導(dǎo)致約4%的射束中心軸百分深度劑量差異和小于1%的器官劑量差異,因此仿真系統(tǒng)生成的能譜文件中射束中心軸百分深度劑量差異均較小,可滿足錐形束CT系統(tǒng)的相關(guān)研究。對于百分深度劑量曲線,在5 cm以內(nèi)深度處差異較大的原因與模擬時未考慮實際照射的多向散射及電離室的精度有關(guān),其均會帶來淺表區(qū)域的差異增加;對于能譜信息30~55 keV能量范圍內(nèi)的差異大于其他區(qū)域的可能原因是初級濾波和射束硬化濾波對低能射線的過濾放大了該區(qū)域的差異,此外本實驗建模的1 mm銅層探測器也對低能部分影響較大。

        圖5 仿真系統(tǒng)等中心位置縱向離軸劑量曲線的Monte Carlo模擬數(shù)據(jù)與二維電離室矩陣實測數(shù)據(jù)比較

        橫向離軸劑量曲線全領(lǐng)結(jié)濾波模式在離軸距離[-10 cm,10 cm]區(qū)間內(nèi)差異略大于[-13 cm,-10 cm)和(10 cm,13 cm]區(qū)間,可能的原因是模塊化建模時忽略支撐領(lǐng)結(jié)濾波的鋼支架(厚度小于1 mm)導(dǎo)致的。全領(lǐng)結(jié)濾波結(jié)構(gòu)為中間開孔,兩邊厚,鋼支架對中心位置的劑量影響被放大,從而導(dǎo)致領(lǐng)結(jié)濾波開孔區(qū)域內(nèi)的劑量差異變大。半領(lǐng)結(jié)濾波模式在離軸距離[-13 cm,6 cm]區(qū)間內(nèi)的差異略大于(6 cm,13 cm]區(qū)間,原因與全領(lǐng)結(jié)濾波類似。盡管存在差異,但是橫向離軸劑量的差異整體在較小的范圍內(nèi),對使用該仿真系統(tǒng)進(jìn)行錐形束CT的其他研究影響甚微。

        縱向離軸劑量曲線中,125 kVp全領(lǐng)結(jié)濾波和125 kVp半領(lǐng)結(jié)濾波模式的曲線肉眼無法區(qū)分其差異,究其原因與領(lǐng)結(jié)濾波的結(jié)構(gòu)有關(guān)。領(lǐng)結(jié)濾波的中心開孔方向與縱向平行,導(dǎo)致125 kVp全領(lǐng)結(jié)濾波采集數(shù)據(jù)時領(lǐng)結(jié)濾波對其造成的影響與半領(lǐng)結(jié)濾波相同(僅與鋼支架發(fā)生作用)。領(lǐng)結(jié)濾波在[-10 cm,10 cm]區(qū)間內(nèi)縱向離軸劑量模擬值與實測值吻合性更好,其余區(qū)域差異略大,究其原因是實際測量使用的OCTAVIUS Detector 729二維電離室矩陣相鄰2個電離室的中心間距為1 cm,而模擬時的體素計算間距為0.5 cm,分辨力的差異造成劑量變化較快的區(qū)域出現(xiàn)更大誤差,此外差異可能還與二維電離室矩陣的測量精度有關(guān)。分析該組數(shù)據(jù)還發(fā)現(xiàn),Monte Carlo模擬值與實測值差異最大的模式為125 kVp全領(lǐng)結(jié)濾波,出現(xiàn)在離軸距離約-12 cm處,差異值為2.7%。同樣,縱向離軸劑量的差異整體上也在較小的范圍內(nèi),4種電壓-濾波模式縱向離軸劑量的差異平均值在±0.8%以內(nèi),對使用該仿真系統(tǒng)進(jìn)行錐形束CT的其他研究時的影響比橫向的影響更小。

        綜上所述,本研究建立的模擬仿真系統(tǒng)可實現(xiàn)錐形束CT的準(zhǔn)確模擬,可為研究錐形束CT的相關(guān)問題奠定基礎(chǔ)。下一步將基于該系統(tǒng)對如何提高錐形束CT圖像質(zhì)量及基于該系統(tǒng)的形變配準(zhǔn)算法等展開研究。

        猜你喜歡
        離軸錐形子程序
        離軸超構(gòu)透鏡設(shè)計與特性分析
        基于自由曲面雙波段離軸三反光學(xué)系統(tǒng)的優(yōu)化設(shè)計
        下頜管在下頜骨內(nèi)解剖結(jié)構(gòu)的錐形束CT測量
        反射式紅外多波段準(zhǔn)直投影光學(xué)系統(tǒng)設(shè)計
        航空兵器(2019年2期)2019-05-30 00:00:00
        寬譜段、動態(tài)局部高分辨離軸主動反射變焦系統(tǒng)
        錐形束CT結(jié)合顯微超聲技術(shù)診治老年鈣化根管的應(yīng)用
        宮頸錐形切除術(shù)后再次妊娠分娩方式的探討
        錐形流量計尾流流場分析
        淺談子程序在數(shù)控車編程中的應(yīng)用
        子程序在數(shù)控車加工槽中的應(yīng)用探索
        亚洲av无码专区亚洲av网站| 91精品国产91| 国产在线不卡免费播放| 亚洲av中文aⅴ无码av不卡| 国产一区二区三区日韩精品 | 免费一级欧美大片久久网| 日韩中文字幕无码av| 日韩一二三四区免费观看| 日本不卡不二三区在线看 | 日日碰日日摸日日澡视频播放| 亚洲日韩精品一区二区三区无码| 久久久精品一区aaa片| 午夜三级a三级三点| 亚洲 欧美 唯美 国产 伦 综合| 国产亚洲女在线线精品| 亚洲女同系列高清在线观看| 日韩字幕无线乱码免费| 精品国产一区二区三区av免费| 免费a级毛片18禁网站免费| 国产午夜伦鲁鲁| 亚洲а∨天堂久久精品2021| 亚洲男人天堂2019| 国产精品一区二区三区精品| 亚洲综合一区二区三区蜜臀av| 日本在线一区二区免费| 日本人妻伦理在线播放| 久久精品国产亚洲av麻豆图片| 国产成人综合久久亚洲精品| 国产精品丝袜黑色高跟鞋| 久久尤物av天堂日日综合| 成人综合亚洲国产成人| 国产91在线播放九色快色| 极品嫩模大尺度av在线播放| 免费人成激情视频在线观看冫| 18成人片黄网站www| 国产亚洲精品看片在线观看| 久久久久久免费播放一级毛片| 日本黄色特级一区二区三区| 成人久久久精品乱码一区二区三区| 小说区激情另类春色| 在线欧美中文字幕农村电影|