• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于時間反轉(zhuǎn)法的相控?fù)Q能器聲場的仿真研究?

    2015-10-26 08:07:14韓珍珍丁鑫羅明艷菅喜岐
    應(yīng)用聲學(xué) 2015年4期
    關(guān)鍵詞:聲強旁瓣同心

    韓珍珍 丁鑫 羅明艷 菅喜岐

    (天津醫(yī)科大學(xué)生物醫(yī)學(xué)工程與技術(shù)學(xué)院天津300070)

    基于時間反轉(zhuǎn)法的相控?fù)Q能器聲場的仿真研究?

    韓珍珍丁鑫羅明艷菅喜岐?

    (天津醫(yī)科大學(xué)生物醫(yī)學(xué)工程與技術(shù)學(xué)院天津300070)

    相控?fù)Q能器具有焦距可調(diào)的優(yōu)勢。本文以82陣元相控?fù)Q能器建立的3D數(shù)值仿真模型為例,基于時間反轉(zhuǎn)法提取陣元的激勵信號,利用時域有限差分(FDTD)法對Westervelt聲波非線性傳播方程進(jìn)行聲場數(shù)值仿真,研究不同陣元分布、焦點偏離聲軸的距離、設(shè)定焦距大小對形成聲場的影響,可調(diào)控范圍及其消除旁瓣方法。研究結(jié)果表明,隨機分布相控陣可明顯降低聲場中的旁瓣;隨著偏離聲軸距離的增加,主瓣聲壓幅值逐漸減小,旁瓣與主瓣的最大聲強比值r逐漸增大,且沿聲軸的可調(diào)控范圍逐漸減??;隨聲軸方向上設(shè)定焦距的增加,主瓣聲壓幅值先增大后減小,r值先減小后增大;基于時間反轉(zhuǎn)法的高聲壓旁瓣消除法可在一定程度上擴大相控陣聲場的可調(diào)控范圍。

    高強度聚焦超聲,相控?fù)Q能器,陣元分布,時間反轉(zhuǎn)法,高聲壓旁瓣消除法

    1 引言

    高強度聚焦超聲腫瘤治療具有無創(chuàng)、可重復(fù)治療等優(yōu)勢,超聲換能器是HIFU腫瘤治療系統(tǒng)中最關(guān)鍵的部件之一。目前單陣元換能器已應(yīng)用于乳腺癌、子宮肌瘤、前列腺癌等實體軟組織腫瘤的臨床治療[1],但其焦距固定,腫瘤治療時需通過機械調(diào)整換能器與病患部位的相對位置進(jìn)行治療定位。而相控陣換能器可通過相位調(diào)控其焦距大小而實現(xiàn)動態(tài)聚焦,該研究受到眾多研究者的關(guān)注。

    1996年Goss等[2]提出球面陣元隨機分布可改善相控陣聲場分布的方法。2005年Lu等[3]采用遺傳算法對256矩形陣元相控陣中每個陣元的幅值和相位進(jìn)行調(diào)控以降低聲場中旁瓣。2009年Hand等[4]通過仿真和離體實驗對254陣元隨機相控陣單焦點和多焦點的動態(tài)調(diào)控范圍進(jìn)行了測定。2012年Ji等[5]設(shè)計了用于治療深部腫瘤的90陣元大開口半隨機分布的球面相控陣。

    本文根據(jù)周期、同心環(huán)和隨機分布設(shè)計的82陣元相控陣換能器建立3D數(shù)值仿真模型,并基于時間反轉(zhuǎn)法提取相位控制信號,利用Westervelt聲波非線性傳播方程進(jìn)行時域有限差分?jǐn)?shù)值仿真研究,以旁瓣與主瓣的最大聲強比值r≤0.25為判斷標(biāo)準(zhǔn)[6]討論不同聚焦位置時的聲場變化規(guī)律、焦點可調(diào)控范圍。并提出擴大焦域可調(diào)控范圍的高聲壓旁瓣消除法。

    2 方法與模型

    2.1聲波傳播方程式

    Westervelt聲波傳播方程式[7-8]為

    (1)式中:?為拉普拉斯算子;p為聲壓;c0和ρ0分別為介質(zhì)的聲速和密度;β=1+B/(2A)為聲波非線性系數(shù),B/A為流體介質(zhì)的非線性系數(shù);為聲波擴散系數(shù);α為吸收系數(shù);ω=2πf為角頻率,f為頻率。

    采用時域有限差分法[9-10]對(1)式進(jìn)行中心差分,差分方程式為

    i、j、k分別為直角坐標(biāo)系下x、y、z三個坐標(biāo)軸方向的坐標(biāo),dx、dy、dz分別表示x、y、z三個坐標(biāo)軸方向的空間步長,dt為時間步長,n為計算時刻。

    2.2陣元相位獲取法

    基于時間反轉(zhuǎn)法[11-12]在目標(biāo)焦點S處設(shè)置正弦函數(shù)點聲源S0(t),數(shù)值仿真得到傳到相控陣編號為m陣元中心點的聲壓信號S0m(t),將該信號按時間序列進(jìn)行反轉(zhuǎn)后,得到對應(yīng)陣元m的信號S0m(T-t)。利用最小二乘函數(shù)擬合法計算一段時間內(nèi)S0m(T-t)的相對初始相位延遲Δtm,然后以同一輸入聲強對正弦信號幅值進(jìn)行調(diào)制,陣元m的激勵信號為

    如果在(7)式激勵信號觸發(fā)相控陣換能器陣元形成的聲場中出現(xiàn)聲壓幅值較高的旁瓣,確定其位置R,采用同樣的方法計算聚焦在R處相控?fù)Q能器陣元m的激勵信號。對分別聚焦在S和R處的激勵信號線性疊加,調(diào)控聚焦在R處的輸入聲強使在R處聚焦時穩(wěn)態(tài)聲壓的最大幅值與聚焦在S處時R位置的穩(wěn)態(tài)聲壓幅值相等,且使聚焦在R處的相位相對聚焦在S處的相位180°延遲,從而達(dá)到降低R處高聲壓旁瓣幅值[13]的目的。

    2.3數(shù)值仿真模型和參數(shù)

    圖1為開口直徑100 mm、曲率半徑80 mm、曲率半徑與開口直徑之比F值為0.8的82陣元相控陣換能器和水體組成的數(shù)值仿真模型。其中數(shù)值仿真區(qū)域為邊長為100 mm的立方體,聲軸為x軸。圖2(a)、2(b)、2(c)分別為陣元按照周期、同心環(huán)和隨機分布構(gòu)成的82陣元相控?fù)Q能器,其中陣元直徑均為8 mm、陣元填充率為52.48%,周期分布相鄰陣元中心間距均為9.5 mm;同心環(huán)分布中陣元中心點所在環(huán)半徑由內(nèi)到外依次為9 mm、18 mm、27 mm、36 mm、45 mm,對應(yīng)的陣元數(shù)目分別為6、11、16、21、27;隨機分布通過Matlab隨機數(shù)生成器獲得,四個象限內(nèi)陣元數(shù)目依次為21、21、20、20,在相鄰陣元中心間距不小于8.7 mm的條件下,依次選定各象限中其余陣元的中心位置。

    水體的仿真參數(shù)為ρ0=998 kg/m3,c0= 1486 m/s,β=3.5,α=0.02 Np/(m·MHz)[14]。換能器激勵函數(shù)為頻率為0.7 MHz的連續(xù)正弦波,邊界采用Mur一階邊界吸收條件??臻g步長dx=dy=dz=0.25 mm,時間步長dt=10 ns。

    圖1 數(shù)值仿真模型圖(單位:mm)Fig.1 Numerical simulation model(Unit:mm)

    圖2 陣列結(jié)構(gòu)圖Fig.2 Array structure diagram

    3 結(jié)果

    3.1不同陣元分布形成的聲場

    當(dāng)輸入聲強為1 W/cm2,設(shè)定焦點為幾何焦點(80,0,0)mm時,周期分布、同心環(huán)分布和隨機分布相控陣在x軸上形成的聲壓變化曲線如圖3所示,其中黑色虛線為周期分布、灰色實線為同心環(huán)分布、黑色實線為隨機分布。由圖3可知,三種分布陣列形成的最大聲壓(主瓣)幅值分別為4.800 MPa、4.716 MPa、4.777 MPa,實際焦距分別為78.25 mm、77.75 mm、77.75 mm,均小于設(shè)定的80 mm幾何焦距,周期分布和同心環(huán)分布陣列在主瓣左側(cè)旁瓣多且幅值較大,隨機分布相控陣僅在靠近主瓣左側(cè)有一個較小旁瓣。圖4為三種不同分布陣列在幾何焦點聚焦時實際形成的焦平面聲場分布圖,由圖4可知周期分布陣列焦域周圍有六個明顯的旁瓣,同心環(huán)分布陣列焦域周圍有兩層旁瓣,隨機分布陣列僅在焦域附近有分布無規(guī)則的較小旁瓣。對上述三種不同分布相控陣在幾何焦點聚焦時實際形成的三維聲場進(jìn)行分析,周期、同心環(huán)、隨機分布陣列中旁瓣與主瓣的最大聲強比值分別為0.1400、0.2358、0.1247,隨機分布陣列的比值最低。

    圖3 不同陣元分布相控陣在x軸上形成的聲壓變化曲線圖Fig.3The acoustic pressure curve along the xaxis of phased arrays with different element distribution

    圖4 在幾何焦點聚焦時的焦平面聲場分布圖Fig.4 The acoustic field distribution of the focal plane when focusing at geometric focus

    3.2偏離聲軸聚焦的聲場分布

    當(dāng)輸入聲強為1 W/cm2,設(shè)定焦點位置為x軸80 mm、z軸0 mm、沿y軸正方向以1 mm為步長偏離x軸聚焦時,形成焦點位置處最大聲壓的變化曲線如圖5(a)所示,其中灰色實線為周期分布、黑色虛線為同心環(huán)分布、黑色實線為隨機分布,由圖5(a)可知隨著偏離x軸距離的增加三種分布相控陣最大聲壓值均逐漸減小,且變化趨勢基本一致。圖5(b)為與圖5(a)相應(yīng)的旁瓣與主瓣最大聲強比值r的變化曲線,由圖5(b)可知隨著偏離x軸距離的增加該比值均逐漸增大,同心環(huán)分布和周期分布陣列的r值均明顯高于隨機分布陣列。由圖3~5可知,與周期和同心環(huán)分布陣列相比,隨機分布陣列可以明顯降低聲場中的旁瓣,為此選用隨機分布相控?fù)Q能器進(jìn)行聲場討論分析。

    圖5 三種不同分布陣列沿y軸聚焦時的聲場變化Fig.5The acoustic variation as a function of focal spot location along the y-axis for three phased arrays with different element distribution

    以圖2(c)隨機分布相控陣為例,當(dāng)設(shè)定焦點為(80,5,0)、(80,10,0)、(80,15,0)時實際形成的焦平面聲場分布如圖6所示,由圖6可知隨著偏離x軸距離的增加,分布在焦域周圍的旁瓣逐漸移至焦域下方,并且旁瓣數(shù)量逐漸增多、幅值逐漸增大。由圖5(b)可知,隨機分布相控陣在偏離x軸8 mm聚焦形成的聲場中r=0.2664,超過了0.25,其xy平面聲場分布如圖7(a)所示,對圖7(a)所示的靠近焦域左下方的旁瓣進(jìn)行高聲壓旁瓣消除后,xy平面的聲場分布如圖7(b)所示,使其旁瓣幅值明顯降低,且r=0.1534。同時高聲壓旁瓣消除后,二次聲壓旁瓣幅值由1.963 MPa下降為1.192 MPa,下降率為39.28%,而主瓣聲壓幅值由3.803 MPa變?yōu)?.403 MPa,下降率僅為10.52%。當(dāng)對隨機分布相控陣進(jìn)行高聲壓旁瓣消除后,可使其實現(xiàn)最大偏離x軸12 mm的調(diào)控聚焦。

    圖6 隨機分布相控陣在不同焦點位置聚焦時歸一化后焦平面聲場分布圖Fig.6 The acoustic field distribution of the focal plane for the quasi-random phased array with different focal location

    圖7 隨機分布相控陣在(80,8,0)聚焦時歸一化后xy平面聲場分布圖Fig.7The acoustic field distribution in the xy plane for the quasi-random phased array when focusing at(80,8,0)

    3.3隨機分布相控陣沿聲軸聚焦

    當(dāng)輸入聲強為1 W/cm2,設(shè)定目標(biāo)焦點在z=0平面,y方向上分別沿x軸、偏離x軸5 mm、10 mm、12 mm,沿x軸方向的焦距范圍均在55~90 mm調(diào)控的條件下,形成聲場中最大聲壓值隨設(shè)定焦點位置的變化曲線如圖8(a)所示,其中灰色虛線為沿x軸、灰色實線為偏離x軸5 mm、黑色虛線為偏離x軸10 mm、黑色實線為偏離x軸12 mm,由圖8(a)可知在y軸方向偏離x軸距離相同的條件下,隨x軸方向上設(shè)定焦距的增加,形成的最大聲壓值先增大后減小,且偏離x軸的距離不同形成最大聲壓的位置也不同,分別為x=70 mm、75 mm、75 mm、80 mm。在x軸方向上調(diào)控位置相同時,隨著y軸方向上偏離x軸距離的增加,形成的最大聲壓值逐漸降低。圖8(b)為與圖8(a)相同條件下隨設(shè)定焦點位置的r值變化曲線,由圖8(b)可知隨著y方向上偏離x軸距離的增加,r值逐漸增大,隨著x軸方向上設(shè)定焦距值的增加r值先減小后增大,形成最小r值的位置有所不同,分別為x=65 mm、65 mm、75 mm、80 mm。在r≤0.25的判斷條件下,x軸上的可調(diào)控范圍為55~90 mm,y方向上偏離x軸5 mm、10 mm沿x軸方向聚焦的可調(diào)控范圍分別為55~90 mm、60~80 mm,y方向上偏離x軸12 mm僅有x=80 mm和x=90 mm兩點。隨著y方向上偏離x軸距離的增加,沿x軸方向的聲場可調(diào)控范圍逐漸減小。同時由圖8(b)可知,在偏離x軸10 mm、12 mm聚焦時,分別有x=55、60、85、90 mm和x=55、60、65、70、75、85 mm位置處r值均大于0.25,且在x=85 mm調(diào)控聚焦時均有一波峰值,分別進(jìn)行高聲壓旁瓣消除后的r值變化曲線如圖8(c)所示,由8(c)可知偏離x軸10 mm、12 mm分別對應(yīng)x軸方向可調(diào)控范圍分別不小于60~90 mm、80~90 mm。焦點在(80,12,0)位置調(diào)控聚焦時形成的xy和yz平面聲場分布如圖9所示。由圖9可知,在焦域附近沒有較大的聲壓旁瓣出現(xiàn)。

    圖8 隨機分布相控陣y方向上偏離x軸不同距離時沿x軸以5 mm步長調(diào)控聚焦時的聲場變化規(guī)律Fig.8 The acoustic variation as a function of focal spot location along the x-axis and off y-axis different distance from x-axis for the quasi-random phased array

    圖9 隨機分布相控陣在(80,12,0)位置聚焦時的歸一化聲場分布圖Fig.9 The acoustic field distribution when focusing at(80,12,0)for the quasi-random phased array

    4 討論與結(jié)論

    本文設(shè)計了周期、同心環(huán)、隨機三種不同分布的相控陣,基于時間反轉(zhuǎn)法進(jìn)行相位調(diào)控,分析比較其聲場分布特點和變化規(guī)律,并結(jié)合基于時間反轉(zhuǎn)法的高聲壓旁瓣消除法增大隨機分布相控陣聲場的可調(diào)控范圍。結(jié)論如下:

    (1)與周期和同心環(huán)分布陣列相比,隨機分布陣列可明顯降低聲場中的旁瓣,改善聲場分布,這與Gavrilov等[6]的研究結(jié)果一致。

    (2)隨聲軸方向上設(shè)定焦點和換能器間距離的增加,主瓣聲壓幅值先增大后減小,旁瓣與主瓣的最大聲強比值先減小后增大。同時隨著偏離聲軸距離的增加,主瓣聲壓幅值逐漸減小,旁瓣與主瓣的最大聲強比值逐漸增大,且沿聲軸的聲場可調(diào)控范圍逐漸減小,該結(jié)果與Raju等[15]的研究結(jié)果相近。

    (3)基于時間反轉(zhuǎn)法的高聲壓旁瓣消除法可在一定程度上增大相控陣聲場的可調(diào)控范圍。

    上述結(jié)論是在高強度聚焦聲場的仿真研究中得出的,對于聲壓幅值較小的超聲波而言,一般呈線性傳播,但對HIFU治療腫瘤而言,其焦點聲壓一般較高,因此需要考慮超聲傳播的非線性特性。

    在基于時間反轉(zhuǎn)法獲取陣元的激勵信號時,聲波的不同傳遞過程中具有一定非線性差異,2010年P(guān)inton等[16]的研究結(jié)果表明對相控陣的相位差計算時其影響可忽略不計。為此本文在獲取陣元激勵信號時未考慮其非線性的影響。

    本文僅以水體內(nèi)相控?fù)Q能器形成的聲壓場進(jìn)行分析討論。對于非均質(zhì)人體組織的仿真研究中,可利用2003年Aubry等[17]提出的基于CT數(shù)據(jù)獲取組織聲學(xué)參數(shù)的方法,建立結(jié)合臨床實際的數(shù)值仿真模型,研究其形成的聲場特性。基于人體CT數(shù)據(jù)的數(shù)值仿真正在進(jìn)行中。

    [1]KENNEDY J E.High-intensity focused ultrasound in the treatment of solid tumours[J].Nature Reviews Cancer,2005,5(4):321-327.

    [2]GOSS S A,F(xiàn)RIZZELL L A,KOUZMANOFF J T,et al.Sparse random ultrasound phased array for focal surgery[J].IEEE Trans.Ultrason.Ferroelectr.Freq. Control,1996,43(6):1111-1121.

    [3]LU Mingzhu,WAN Mingxi,XU Feng,et al.Focused beam control for ultrasound surgery with spherical-section phased array:sound field calculation and genetic optimization algorithm[J].IEEE Trans.Ultrason.Ferroelectr. Freq.Control,2005,52(8):1270-1290.

    [4]HAND J W,SHAW A,SADHOO N,et al.A random phased array device for delivery of high intensity focused ultrasound[J].Phys.Med.Biol.,2009,54(19):5675-5693.

    [5]JI Xiang,SHEN Guofeng,BAI Jingfeng,et al.The characterization of an ultrasound spherical phased array for the ablation of deep-seated tissue[J].Appl.Acoust.,2012,73(5):529-534.

    [6]GAVRILOV L R,HAND J W.A theoretical assessment of the relative performance of spherical phased arrays for ultrasound surgery[J].IEEE Trans.Ultrason.Ferroelectr. Freq.Control,2000,47(1):125-139.

    [7]HAMILTON M F,BLOCKSTOCK D T.Nonlinear Acoustics[M].Boston:Academic Press,1998.

    [8]BAILEY M,KHOKHLOVA V,SAPOZHNIKOV O,et al. Physical mechanism of the therapeutic effect of ultrasound(A review)[J].Acoustic Physics,2003,49(4):369-388.

    [9]YEE K S.Numerical solution of initial boundary value problems involving Maxwell's equations[J].IEEE Trans. Antennas Propag,1966,14(3):302-307.

    [10]RADZEVICIUS S J,CHEN C C,ETERS J L,et al.Nearfield dipole radiation dynamics through FDTD modeling[J].Journal of Applied Geophysics,2003,52(2):75-91.

    [11]THOMAS J L,F(xiàn)INK M A.Ultrasonic beam focusing through tissue inhomogeneities with a time reversal mirror:application to transskull therapy[J].IEEE Trans. Ultrason.Ferroelectr.Freq.Control,1996,43(6):1122-1129.

    [12]王意喆,張千,周文征,等.HIFU經(jīng)顱治療中降低顱骨聲壓的仿真研究[J].計算機仿真,2014,31(5):209-213. WANG Yizhe,ZHANG Qian,ZHOU Wenzheng,et al. Simulation study on reducing acoustic pressure on skull during HIFU transcranial therapy[J].Computer Simulation,2014,31(5):209-213.

    [13]LEDUC N,OKITA K,SUGIYAMA K,et al.Focus control in HIFU therapy assisted by time-reversal simulation with an iterative procedure for hot spot elimination[J]. Journal of Biomechanical Science and Engineering,2012,7(1):43-56.

    [14]GINTERS.Numericalsimulationofultrasoundthermotherapy combining nonlinear wave propagation with broadband soft-tissue absorption[J].Ultrasonics,2000,37(10):693-696.

    [15]RAJU B I,HALL C S,SEIP R.Ultrasound therapy transducers with space-filling non-periodic arrays[J].IEEE Trans.Ultrason.Ferroelectr.Freq.Control,2011,58(5):944-954.

    [16]PINTON G,AUBRY J F,F(xiàn)INK M,et al.Effects of nonlinear ultrasound propagation on high intensity brain therapy[J].Medical Physics,2011,38(3):1207-1216.

    [17]AUBRY J F,TANTER M,PERNOT M,et al.Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans[J]. The Journal of the Acoustical Society of America,2003,113(1):84-93.

    The simulation study of acoustic pressure field about HIFU phased array transducer based on time reversal method

    HAN ZhenzhenDING XinLUO MingyanJIAN Xiqi
    (Department of Biomedical Engineering and Technology,Tianjin Medical University,Tianjin 300070,China)

    Phased array transducer has the advantage of adjustable focal length.In this paper,the incentive signal of each element is abstracted by time reversal method and the acoustic pressure numerical simulation is undergone through finite difference time domain(FDTD)of Westervelt acoustic nonlinear transmission formula based on the 3D numerical simulation model of 82-element phased array transducer.The effect of different element distribution,the distance off the acoustic axis and the setting focal length on the forming acoustic pressure field is discussed.The adjustable range of acoustic pressure field and the method which eliminate sidelobe is also analyzed.The results show that the phased array transducer with quasi-random distribution can suppress the sidelobe significantly.As the distance off the acoustic axis increasing,the acoustic pressure amplitude of the main lobe decreases,the acoustic intensity ratio between the maximum sidelobe and the main lobe(r value)increases,and the adjustable range along the acoustic axis decreases.The acoustic pressure amplitude of the main lobe increases firstly and then decreases,while r value decreases firstly and then increases with the increase of setting focal length along the acoustic axis.To some extent,the acoustic adjustable range can be extended by the high acoustic pressure sidelobe elimination method based on time reversal method.

    High-intensity focused ultrasound,Phased array transducer,Element distribution,Time reversal method,High acoustic pressure sidelobe elimination method

    R318.5

    A

    1000-310X(2015)04-0344-07

    10.11684/j.issn.1000-310X.2015.04.009

    2014-10-29收稿;2015-02-16定稿

    ?國家自然科學(xué)基金項目(81272495)

    韓珍珍(1989-),女,山東濰坊人,碩士研究生,研究方向:超聲醫(yī)學(xué)。?

    E-mail:jianxiqi@tmu.edu.cn

    猜你喜歡
    聲強旁瓣同心
    三維聲強與振速波束形成結(jié)合的聲源識別方法研究
    超音速聲強近似測量方法
    同心戰(zhàn)"疫" 攜手前行
    福利中國(2022年3期)2022-07-20 08:05:00
    基于圓柱陣通信系統(tǒng)的廣義旁瓣對消算法
    同心逐夢
    一種基于線性規(guī)劃的頻率編碼旁瓣抑制方法
    勠力同心 堅韌前行
    河北畫報(2020年3期)2020-06-24 06:03:42
    勠力同心 堅韌前行
    河北畫報(2020年5期)2020-06-10 01:50:12
    基于加權(quán)積分旁瓣最小化的隨機多相碼設(shè)計
    基于四項最低旁瓣Nuttall窗的插值FFT諧波分析
    免费不卡的大黄色大毛片视频在线观看| 亚洲av欧美aⅴ国产| 热re99久久精品国产66热6| 成人黄色视频免费在线看| tube8黄色片| 成人影院久久| 日本免费在线观看一区| 国产成人欧美| 中文字幕最新亚洲高清| 国产成人精品福利久久| 亚洲精品中文字幕在线视频| 大片电影免费在线观看免费| 你懂的网址亚洲精品在线观看| 九九爱精品视频在线观看| 久久国产亚洲av麻豆专区| 免费人妻精品一区二区三区视频| 日韩在线高清观看一区二区三区| 久久99热6这里只有精品| 国产精品免费大片| 最黄视频免费看| 欧美最新免费一区二区三区| 成人午夜精彩视频在线观看| 性色av一级| 天天躁夜夜躁狠狠躁躁| 99热这里只有是精品在线观看| 中国国产av一级| 国产欧美日韩综合在线一区二区| 日本-黄色视频高清免费观看| 男人添女人高潮全过程视频| 青春草国产在线视频| 国产1区2区3区精品| av免费在线看不卡| 日韩一本色道免费dvd| 国产在线视频一区二区| 又黄又爽又刺激的免费视频.| 国产精品一二三区在线看| 日日摸夜夜添夜夜爱| 国产男人的电影天堂91| 如日韩欧美国产精品一区二区三区| 老司机影院毛片| 精品人妻在线不人妻| 精品99又大又爽又粗少妇毛片| 久久99热6这里只有精品| 亚洲国产精品一区三区| 精品一区二区免费观看| 日韩av不卡免费在线播放| 日本免费在线观看一区| 亚洲精品美女久久av网站| 亚洲精华国产精华液的使用体验| tube8黄色片| 飞空精品影院首页| 狂野欧美激情性bbbbbb| 国产亚洲精品久久久com| 成人漫画全彩无遮挡| 五月玫瑰六月丁香| 成人国产麻豆网| 久久久亚洲精品成人影院| 插逼视频在线观看| 一级片免费观看大全| 日本黄色日本黄色录像| 青春草视频在线免费观看| 韩国高清视频一区二区三区| 亚洲四区av| 成人18禁高潮啪啪吃奶动态图| 一本色道久久久久久精品综合| 国产精品偷伦视频观看了| 母亲3免费完整高清在线观看 | 中文字幕人妻丝袜制服| 老熟女久久久| 国产成人aa在线观看| 午夜福利,免费看| 亚洲人与动物交配视频| 高清不卡的av网站| √禁漫天堂资源中文www| 国产一区二区在线观看av| 欧美日韩成人在线一区二区| 亚洲第一区二区三区不卡| 国产在线视频一区二区| 午夜视频国产福利| 国产精品久久久av美女十八| 久久久久久人妻| 激情视频va一区二区三区| 国产精品一区www在线观看| 大片免费播放器 马上看| 丰满迷人的少妇在线观看| 国产成人免费观看mmmm| 少妇人妻 视频| 日本与韩国留学比较| 国产男女超爽视频在线观看| 捣出白浆h1v1| 亚洲性久久影院| 精品国产一区二区三区四区第35| 建设人人有责人人尽责人人享有的| 最近手机中文字幕大全| 精品少妇内射三级| 日本av手机在线免费观看| 欧美另类一区| 视频中文字幕在线观看| 亚洲精品国产色婷婷电影| 久久久久视频综合| 精品一区二区三区四区五区乱码 | 黑人巨大精品欧美一区二区蜜桃 | 亚洲av日韩在线播放| 欧美最新免费一区二区三区| 蜜臀久久99精品久久宅男| 美国免费a级毛片| 国产国语露脸激情在线看| 三级国产精品片| 国产黄色免费在线视频| 美女大奶头黄色视频| 免费女性裸体啪啪无遮挡网站| 午夜免费男女啪啪视频观看| 十八禁网站网址无遮挡| videosex国产| 一级毛片 在线播放| 啦啦啦中文免费视频观看日本| 成年人午夜在线观看视频| 十八禁高潮呻吟视频| 久久国产精品大桥未久av| 成人国产麻豆网| 亚洲成人手机| 国产精品久久久av美女十八| 蜜臀久久99精品久久宅男| av天堂久久9| 激情视频va一区二区三区| a级片在线免费高清观看视频| 日本与韩国留学比较| 精品少妇久久久久久888优播| 18禁动态无遮挡网站| 五月开心婷婷网| 精品一区二区免费观看| 亚洲欧美一区二区三区国产| 国产欧美亚洲国产| 国产男人的电影天堂91| 黄网站色视频无遮挡免费观看| 亚洲在久久综合| 中文欧美无线码| 亚洲欧美日韩卡通动漫| 日本欧美视频一区| 香蕉国产在线看| xxxhd国产人妻xxx| 午夜影院在线不卡| 国产一区二区三区av在线| 免费观看无遮挡的男女| 91精品国产国语对白视频| 在线观看人妻少妇| 1024视频免费在线观看| 男男h啪啪无遮挡| 熟妇人妻不卡中文字幕| 男的添女的下面高潮视频| 国产高清国产精品国产三级| 久久久久久人妻| 日韩,欧美,国产一区二区三区| 亚洲美女视频黄频| 九色亚洲精品在线播放| www.av在线官网国产| 日本vs欧美在线观看视频| 午夜91福利影院| 欧美亚洲日本最大视频资源| 国产高清不卡午夜福利| 一边亲一边摸免费视频| 中文字幕制服av| 国产av码专区亚洲av| 自拍欧美九色日韩亚洲蝌蚪91| 色婷婷久久久亚洲欧美| 国产精品人妻久久久久久| 亚洲色图 男人天堂 中文字幕 | av福利片在线| 免费日韩欧美在线观看| 不卡视频在线观看欧美| 中文精品一卡2卡3卡4更新| 亚洲精品第二区| 亚洲国产日韩一区二区| √禁漫天堂资源中文www| 欧美成人午夜精品| 国产熟女欧美一区二区| 免费大片18禁| 在线观看一区二区三区激情| 一区在线观看完整版| 国语对白做爰xxxⅹ性视频网站| 国产成人精品一,二区| 91午夜精品亚洲一区二区三区| 色视频在线一区二区三区| av线在线观看网站| 天天操日日干夜夜撸| 9色porny在线观看| 精品人妻一区二区三区麻豆| 大香蕉久久成人网| 欧美 亚洲 国产 日韩一| 午夜免费男女啪啪视频观看| www.色视频.com| 美女国产高潮福利片在线看| 尾随美女入室| xxx大片免费视频| 建设人人有责人人尽责人人享有的| 女性生殖器流出的白浆| 国产日韩一区二区三区精品不卡| 一个人免费看片子| 日本色播在线视频| 亚洲欧美中文字幕日韩二区| 人人妻人人澡人人爽人人夜夜| av片东京热男人的天堂| 香蕉精品网在线| 中文字幕另类日韩欧美亚洲嫩草| 免费观看无遮挡的男女| 国产午夜精品一二区理论片| 成人国产av品久久久| 色哟哟·www| 国产一级毛片在线| 久久久久久久久久久久大奶| 成人无遮挡网站| 久久婷婷青草| 美女国产视频在线观看| 蜜桃国产av成人99| 亚洲一码二码三码区别大吗| 91精品伊人久久大香线蕉| 青春草视频在线免费观看| 午夜av观看不卡| 一本大道久久a久久精品| 深夜精品福利| 丰满乱子伦码专区| 一个人免费看片子| 少妇人妻 视频| 国产 精品1| 久久女婷五月综合色啪小说| 久久久久人妻精品一区果冻| 午夜影院在线不卡| 精品一区二区免费观看| 欧美 日韩 精品 国产| 日本-黄色视频高清免费观看| 国产日韩一区二区三区精品不卡| 亚洲丝袜综合中文字幕| 日韩制服骚丝袜av| 婷婷色av中文字幕| 久久免费观看电影| 欧美人与性动交α欧美软件 | 黄色毛片三级朝国网站| 亚洲精华国产精华液的使用体验| 五月玫瑰六月丁香| 黑人巨大精品欧美一区二区蜜桃 | 18禁动态无遮挡网站| 精品福利永久在线观看| 亚洲国产最新在线播放| 成年人免费黄色播放视频| 久久久国产一区二区| 韩国高清视频一区二区三区| 交换朋友夫妻互换小说| 亚洲成人手机| 九草在线视频观看| 男女高潮啪啪啪动态图| 亚洲三级黄色毛片| 国产成人免费观看mmmm| 亚洲av综合色区一区| 国产免费福利视频在线观看| 性色avwww在线观看| 亚洲国产精品999| 爱豆传媒免费全集在线观看| 高清毛片免费看| 视频在线观看一区二区三区| 热99国产精品久久久久久7| 久久精品国产亚洲av涩爱| 自拍欧美九色日韩亚洲蝌蚪91| 国产爽快片一区二区三区| 天天躁夜夜躁狠狠久久av| 丝袜美足系列| 日韩,欧美,国产一区二区三区| 国产精品欧美亚洲77777| 久久久久久人妻| 只有这里有精品99| 久久久久网色| 熟女电影av网| 黄色怎么调成土黄色| 精品国产一区二区三区四区第35| 国产成人91sexporn| 久久精品久久久久久久性| 人成视频在线观看免费观看| 国产在视频线精品| 亚洲成色77777| 综合色丁香网| 18禁国产床啪视频网站| 国产精品三级大全| 丝袜喷水一区| www.色视频.com| 香蕉精品网在线| 男女边摸边吃奶| 欧美 亚洲 国产 日韩一| 性高湖久久久久久久久免费观看| 国产成人av激情在线播放| 国产色婷婷99| www.色视频.com| 亚洲国产欧美日韩在线播放| 亚洲精品,欧美精品| 五月天丁香电影| 欧美最新免费一区二区三区| av在线播放精品| 免费看不卡的av| 少妇人妻久久综合中文| 天天躁夜夜躁狠狠躁躁| 国产熟女午夜一区二区三区| 亚洲精品中文字幕在线视频| 免费黄频网站在线观看国产| 日韩伦理黄色片| 久久人人爽人人片av| 男女下面插进去视频免费观看 | 男女边摸边吃奶| 夜夜骑夜夜射夜夜干| 九草在线视频观看| 免费观看无遮挡的男女| 国产在线免费精品| 一二三四中文在线观看免费高清| 久久久久精品人妻al黑| 人人澡人人妻人| 午夜免费男女啪啪视频观看| 激情五月婷婷亚洲| 久热这里只有精品99| 久久精品久久久久久噜噜老黄| 人妻系列 视频| 国产女主播在线喷水免费视频网站| 国产免费一区二区三区四区乱码| 国产白丝娇喘喷水9色精品| 久久综合国产亚洲精品| 欧美激情 高清一区二区三区| 久久久国产一区二区| av在线老鸭窝| 成人国语在线视频| 91成人精品电影| 99热国产这里只有精品6| 久久99热6这里只有精品| 欧美精品人与动牲交sv欧美| av女优亚洲男人天堂| 在线天堂最新版资源| 精品卡一卡二卡四卡免费| 欧美日本中文国产一区发布| 日韩成人av中文字幕在线观看| 赤兔流量卡办理| 精品午夜福利在线看| 十八禁高潮呻吟视频| 少妇人妻 视频| 久久久久久人人人人人| 99香蕉大伊视频| 性色av一级| 一级,二级,三级黄色视频| 久久久久精品性色| 精品福利永久在线观看| 狂野欧美激情性bbbbbb| 男女啪啪激烈高潮av片| 欧美xxxx性猛交bbbb| 国产免费一级a男人的天堂| 日韩免费高清中文字幕av| 精品国产一区二区三区久久久樱花| 高清不卡的av网站| 五月开心婷婷网| 日日爽夜夜爽网站| 亚洲av电影在线进入| 日日爽夜夜爽网站| 亚洲av电影在线进入| av卡一久久| 岛国毛片在线播放| 久久国产亚洲av麻豆专区| 国产探花极品一区二区| 啦啦啦视频在线资源免费观看| 亚洲久久久国产精品| 国产日韩一区二区三区精品不卡| 亚洲久久久国产精品| 黄片播放在线免费| 久久精品国产亚洲av涩爱| 视频中文字幕在线观看| 尾随美女入室| 免费观看性生交大片5| 99国产精品免费福利视频| 精品第一国产精品| 亚洲国产日韩一区二区| 99久久综合免费| 看十八女毛片水多多多| 亚洲高清免费不卡视频| 国产av精品麻豆| 久久午夜福利片| 黑丝袜美女国产一区| 亚洲婷婷狠狠爱综合网| 精品熟女少妇av免费看| 国产 精品1| 亚洲激情五月婷婷啪啪| 欧美丝袜亚洲另类| 美女主播在线视频| 伊人久久国产一区二区| 日韩大片免费观看网站| 在线观看一区二区三区激情| 一级a做视频免费观看| 成人午夜精彩视频在线观看| 久久ye,这里只有精品| 国产午夜精品一二区理论片| 视频中文字幕在线观看| 国产免费现黄频在线看| 一区二区三区乱码不卡18| 99精国产麻豆久久婷婷| 在线 av 中文字幕| 男女啪啪激烈高潮av片| 欧美3d第一页| 久久精品国产自在天天线| 久热久热在线精品观看| 国产国语露脸激情在线看| 久久精品国产鲁丝片午夜精品| 菩萨蛮人人尽说江南好唐韦庄| 免费黄网站久久成人精品| a级毛片在线看网站| 国语对白做爰xxxⅹ性视频网站| 新久久久久国产一级毛片| 热re99久久精品国产66热6| 边亲边吃奶的免费视频| 青春草亚洲视频在线观看| 在线观看一区二区三区激情| 国产极品粉嫩免费观看在线| 韩国高清视频一区二区三区| 男女免费视频国产| 久久韩国三级中文字幕| 国产午夜精品一二区理论片| 成人国产av品久久久| 好男人视频免费观看在线| 热99久久久久精品小说推荐| 色吧在线观看| 久久久久久久大尺度免费视频| 国产成人精品无人区| 日日爽夜夜爽网站| 少妇的逼水好多| 成人国语在线视频| 18禁观看日本| 2022亚洲国产成人精品| 我要看黄色一级片免费的| 看免费成人av毛片| 国产日韩一区二区三区精品不卡| 免费少妇av软件| 久久久久久伊人网av| 欧美日韩亚洲高清精品| 蜜臀久久99精品久久宅男| 亚洲国产色片| 成人亚洲欧美一区二区av| 综合色丁香网| 热99久久久久精品小说推荐| 中文精品一卡2卡3卡4更新| 男女午夜视频在线观看 | 亚洲精品中文字幕在线视频| 久久免费观看电影| 成年人午夜在线观看视频| 22中文网久久字幕| 欧美3d第一页| 色婷婷av一区二区三区视频| 日韩视频在线欧美| 中国国产av一级| 国产成人精品福利久久| 一本久久精品| 国产精品 国内视频| 国产无遮挡羞羞视频在线观看| av福利片在线| 国产精品久久久av美女十八| 成人影院久久| 精品酒店卫生间| 一级毛片我不卡| 久久99一区二区三区| 韩国av在线不卡| av免费观看日本| 久久精品国产a三级三级三级| av片东京热男人的天堂| 午夜视频国产福利| 免费观看性生交大片5| 99国产综合亚洲精品| 黄色视频在线播放观看不卡| 日本午夜av视频| 国产乱人偷精品视频| xxx大片免费视频| 国产精品人妻久久久久久| 亚洲,一卡二卡三卡| 欧美日韩综合久久久久久| 久久久久网色| 男人操女人黄网站| videos熟女内射| 九色成人免费人妻av| av视频免费观看在线观看| 国产综合精华液| 999精品在线视频| 欧美日韩视频精品一区| 精品亚洲成a人片在线观看| 国产亚洲欧美精品永久| 午夜福利,免费看| 亚洲高清免费不卡视频| 一级片免费观看大全| 妹子高潮喷水视频| 最近最新中文字幕免费大全7| 日韩电影二区| 人人妻人人添人人爽欧美一区卜| 在线天堂最新版资源| 亚洲欧美日韩另类电影网站| 国产一区二区在线观看日韩| 一个人免费看片子| 美女xxoo啪啪120秒动态图| 日韩熟女老妇一区二区性免费视频| 成人无遮挡网站| 日韩一本色道免费dvd| 在线亚洲精品国产二区图片欧美| 啦啦啦在线观看免费高清www| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 国产熟女午夜一区二区三区| 国产精品麻豆人妻色哟哟久久| 久久午夜综合久久蜜桃| 99精国产麻豆久久婷婷| 99热全是精品| 国产精品秋霞免费鲁丝片| 日韩视频在线欧美| 久久青草综合色| 亚洲第一区二区三区不卡| 亚洲精品第二区| 国产精品秋霞免费鲁丝片| 久久这里有精品视频免费| 日日爽夜夜爽网站| 成人毛片a级毛片在线播放| 国产无遮挡羞羞视频在线观看| 菩萨蛮人人尽说江南好唐韦庄| 午夜免费观看性视频| 久久精品久久精品一区二区三区| 免费观看在线日韩| 婷婷成人精品国产| 免费不卡的大黄色大毛片视频在线观看| 另类亚洲欧美激情| 亚洲国产最新在线播放| 欧美亚洲 丝袜 人妻 在线| 精品亚洲乱码少妇综合久久| 狂野欧美激情性xxxx在线观看| 亚洲人成77777在线视频| 国产精品一区www在线观看| 免费观看a级毛片全部| 国产高清三级在线| 人人妻人人添人人爽欧美一区卜| 寂寞人妻少妇视频99o| 亚洲精品中文字幕在线视频| 精品午夜福利在线看| 十八禁网站网址无遮挡| 成年人午夜在线观看视频| 日韩中字成人| 秋霞伦理黄片| 国产成人免费无遮挡视频| 国产永久视频网站| 国产在线一区二区三区精| 国国产精品蜜臀av免费| 女性生殖器流出的白浆| 久久久久国产精品人妻一区二区| 搡女人真爽免费视频火全软件| 久久狼人影院| 男男h啪啪无遮挡| 亚洲精品自拍成人| 久久精品国产亚洲av天美| av一本久久久久| 色94色欧美一区二区| 日韩 亚洲 欧美在线| 最新中文字幕久久久久| 久久这里有精品视频免费| 国产精品久久久久久av不卡| 精品人妻偷拍中文字幕| 九九爱精品视频在线观看| 91久久精品国产一区二区三区| 国产成人精品在线电影| 亚洲五月色婷婷综合| 十分钟在线观看高清视频www| 伦精品一区二区三区| 欧美人与性动交α欧美软件 | 亚洲国产毛片av蜜桃av| 成人漫画全彩无遮挡| 在线天堂最新版资源| 国产精品久久久久久精品古装| 亚洲精品色激情综合| 久久精品国产亚洲av涩爱| 曰老女人黄片| 国产一区二区在线观看av| 亚洲国产精品999| 色婷婷久久久亚洲欧美| 国产国语露脸激情在线看| 免费看不卡的av| 国产精品久久久久成人av| 国产精品一国产av| 午夜福利乱码中文字幕| 午夜老司机福利剧场| av在线观看视频网站免费| 十八禁高潮呻吟视频| 观看美女的网站| 最近最新中文字幕大全免费视频 | 国产精品免费大片| 亚洲人成77777在线视频| 欧美丝袜亚洲另类| 少妇的逼好多水| videosex国产| 丰满饥渴人妻一区二区三| 少妇的丰满在线观看| 免费看光身美女| 亚洲欧美中文字幕日韩二区| 日本vs欧美在线观看视频| 亚洲久久久国产精品| 水蜜桃什么品种好| 99久久中文字幕三级久久日本| 捣出白浆h1v1| 黄色视频在线播放观看不卡| 日日撸夜夜添| 亚洲成av片中文字幕在线观看 | 国产一级毛片在线| 国国产精品蜜臀av免费| 一级毛片黄色毛片免费观看视频| 亚洲精品中文字幕在线视频| 多毛熟女@视频| 肉色欧美久久久久久久蜜桃| 18禁观看日本| 成人漫画全彩无遮挡| 天天操日日干夜夜撸| 97在线人人人人妻| 在线观看美女被高潮喷水网站| 母亲3免费完整高清在线观看 | 久久久欧美国产精品| 一本大道久久a久久精品| 99热国产这里只有精品6|