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

        ?

        基于聲輻射模態(tài)的聲場重建中的測點優(yōu)化方法

        2017-03-09 08:29:28蘇俊博朱海潮毛榮富蘇常偉
        振動與沖擊 2017年3期
        關(guān)鍵詞:階數(shù)聲壓聲場

        蘇俊博, 朱海潮, 毛榮富, 郭 亮, 蘇常偉

        (海軍工程大學 船舶振動噪聲重點實驗室,武漢 430033)

        基于聲輻射模態(tài)的聲場重建中的測點優(yōu)化方法

        蘇俊博, 朱海潮, 毛榮富, 郭 亮, 蘇常偉

        (海軍工程大學 船舶振動噪聲重點實驗室,武漢 430033)

        少測點條件下,利用聲輻射模態(tài)理論重建聲場時,測點布置方式是決定聲場重建精度的一個關(guān)鍵因素。為求得最佳測點布置方式,提出了一種基于聲輻射模態(tài)的最優(yōu)測點選擇方法,即基于聲輻射模態(tài)矩陣的奇異值分解,采用循環(huán)迭代的方式,逐次去除對其最小奇異值最敏感的測點,從而得到了一組使重建方程條件數(shù)最小的測點。實驗結(jié)果表明,利用文中提出的最優(yōu)測點選擇方法布置測點,能夠?qū)β晥鲞M行有效的重建,重建效果優(yōu)于均勻布置方式,顯著提高了聲場重建精度。

        聲輻射模態(tài);聲場重建;測點優(yōu)化

        近場聲全息技術(shù)自20世紀80年代由MAYNARD等[1]提出后,許多學者進行了研究,目前已經(jīng)在噪聲源識別與定位、聲場可視化等方面取得了良好的效果。但近場聲全息技術(shù)要求在全息面上布置大量測點,特別是應用于大型結(jié)構(gòu)時,這限制了近場聲全息技術(shù)在工程實際中的應用。因此如何減少測點數(shù)目是一個需要解決的重要問題。HELS(Helmholtz Equation-Least Squares)方法[2-4]是一種所需測點數(shù)目遠少于近場聲全息技術(shù)的聲場重建方法,其原理是是將輻射聲場表示為Helmholtz方程一組特解的線性組合,利用聲場測量數(shù)據(jù)基于最小二乘法求解組合系數(shù),然后再求解聲場中的場點聲壓。HELS方法對類似球形聲源的聲場能夠?qū)崿F(xiàn)較為精確的重建,但對長寬比較大的結(jié)構(gòu)聲源聲場難以進行有效的重建。聲輻射模態(tài)理論的提出為解決少測點條件下任意結(jié)構(gòu)聲源的聲場重建問題提供了新的途徑[5-6]。

        聲輻射模態(tài)理論是BORGIOTTI等[7-8]學者在九十年代初提出,目前已在結(jié)構(gòu)聲輻射主動控制以及輻射聲功率的計算方面取得了許多成果[9-13]。聲輻射模態(tài)就是輻射體表面一種可能的輻射形式,任何表面法向速度或聲壓都可通過這些聲輻射模態(tài)的線性組合來表示,利用聲輻射模態(tài)進行聲場重建的整個過程可分為求解模態(tài)展開系數(shù),聲場聲壓求解兩個步驟。由此可見,聲輻射模態(tài)理論與HELS方法具有相似性,但是由于聲輻射模態(tài)中包含了結(jié)構(gòu)表面的幾何形狀信息,因此聲輻射模態(tài)理論適用于任意結(jié)構(gòu)輻射聲場的重建。而在少測點條件下利用聲輻射模態(tài)進行聲場重建時,測點布置方式是一個決定聲場重建精度的關(guān)鍵問題,如果測點的位置分布不當會使各測點之間有較強的相關(guān)性,在測點較少的情況下就無法采集到足夠的有效數(shù)據(jù),從而使求解精度降低。此外,在已知測點布置方式的情況下,聲輻射模態(tài)階數(shù)的選擇也對求解精度具有重要影響,這是由于利用已知測點求解模態(tài)展開系數(shù)是一個數(shù)學上的逆問題,模態(tài)階數(shù)過少不足以表示,模態(tài)階數(shù)過多又會使受到噪聲污染的高階模態(tài)參與到求解過程中,從而使求解精度降低。

        宋艷華等[14]研究了利用聲輻射模態(tài)進行聲源識別時測點位置對識別效果的影響,并得出結(jié)論:測點位置分布不均勻時,識別效果較好。但沒有對測點分布方式作進一步研究。KAMMER[15]提出了一種有效獨立方法(Effective Independence Method, EFI)。KIM[16]將該方法用于基于邊界元方法的聲場重建中,在聲場中選取了一組線性獨立的測點,但該方法的局限性在于選取的測點數(shù)量不能少于結(jié)構(gòu)表面離散的單元數(shù)。針對EFI方法的這一局限性,ZHANG[17]提出了一種基于奇異值分解的循環(huán)迭代方法。為解決利用聲輻射模態(tài)重建聲場時所面臨的測點布置問題,本文提出了一種基于聲輻射模態(tài)的最優(yōu)測點選擇方法。針對聲輻射模態(tài)階數(shù)選取問題,將一種基于測量的迭代方法[18]引入到重建過程中,確定了重建過程中所需選取的聲輻射模態(tài)階數(shù)。文章最后開展了雙層圓柱殼體實驗研究,通過與測點均勻布置方式下聲場重建效果的比較,驗證了本文提出的基于聲輻射模態(tài)的最優(yōu)測點選擇方法的有效性。

        1 基于聲輻射模態(tài)的聲場重建理論

        將輻射體表面離散成一系列等面積的輻射單元,由文獻[7]可知,輻射聲功率可表示為

        (1)

        R=ΦΛΦH

        (2)

        式中:Φ為N×N維矩陣,Φ的列向量φi(i=1,2,…,N)即為輻射體的聲輻射模態(tài)向量。

        在得到了輻射體的聲輻射模態(tài)以后,輻射體表面的速度向量可表示成如下形式

        v=Φc

        (3)

        式中:c為輻射模態(tài)展開系數(shù)。假設將輻射體表面離散成N個輻射單元,并計算得到了N階聲輻射模態(tài),由于式(3)具有良好的收斂性[14],表面速度可用截斷的聲輻射模態(tài)表示為

        v≈Φ(N×N1)c(N1×1)

        (4)

        式中:矩陣Φ(N×N1)表示只取了前N1階聲輻射模態(tài),當輻射體表面速度向量中有N2個元素是已知的,可以得到由這N2個元素組成的方程組

        v′=Φ′(N2×N1)c(N1×1)

        (5)

        由式(4)和式(5)可以求得表面速度向量

        (6)

        利用聲輻射模態(tài)由少量測點聲壓值重建輻射體表面聲壓具有相同的過程。在得到了輻射體表面的聲壓或振速值后,就可利用邊界元等方法方便地計算出整個聲場分布。

        2 最優(yōu)測點選擇方法

        針對測點選擇問題,文獻[17]用一種基于奇異值分解的循環(huán)迭代方法,基于使輻射體表面離散點與聲場中測點之間的傳遞函數(shù)矩陣條件數(shù)最小的準則,求得了聲場中的一組測點。本文將其思想應用于基于聲輻射模態(tài)的聲場重建中,提出一種基于聲輻射模態(tài)的最優(yōu)測點選擇方法。假定輻射體表面候選測點數(shù)量為N,從這N個候選測點中選取N2個作為實際布置傳感器的位置。

        假設輻射體表面聲壓向量為p0(N×1),可用聲輻射模態(tài)表示為

        p0(N×1)=Φ(N×N)c(N×1)

        (7)

        對Φ(N×N)進行奇異值分解

        Φ(N×N)=USVH

        (8)

        式中:S為奇異值矩陣,S對角線元素為聲輻射模態(tài)矩陣Φ(N×N)的奇異值。利用式(7)、式(8)求得的模態(tài)展開系數(shù)重建輻射體表面聲壓pr(N×1)

        pr(N×1)=Φ(N×N)c(N×1)=ΦΦ+p0=

        USVHVS+UHp0(N×1)=USS+UHp0(N×1)

        (9)

        將找到的對最小奇異值最敏感的測點從p0(N×1)中去除,并將模態(tài)矩陣Φ(N×N)中與該測點對應的行去除,并此時式(7)變?yōu)?/p>

        p0[(N-1)×1)]=Φ[(N-1)×N)c(N×1)

        (10)

        利用式(10)求解模態(tài)展開系數(shù)時,式(10)為一欠定方程,為使方程有解,去除模態(tài)矩陣Φ[(N-1)×N]中的最高階模態(tài)向量,也即模態(tài)矩陣的最后一列,得到新的模態(tài)矩陣Φ1[(N-1)×(N-1)],則式(10)變?yōu)?/p>

        p01[(N-1)×1]=Φ1[(N-1)×(N-1)]c1[(N-1)×1]

        (11)

        然后以式(11)代替式(7)。至此,由式(7)得到式(11)并以式(11)代替式(7)的過程完成了一次循環(huán)迭代,以相同過程迭代N-N2+1次,得到如下方程

        p0(N-N2)(N2×1)=Φ(N-N2)(N2×N2)c(N-N2)(N2×1)

        (12)

        則p0(N-N2)(N2×1)就是我們需要的N2個測點,式(12)即為由此N2個測點組成的重建方程。

        由以上推導過程可知,該方法是基于模態(tài)矩陣的奇異值分解,采用了一種循環(huán)迭代的方式,逐次去除對最小奇異值最敏感的測點,直至達到要求的測點數(shù)目。經(jīng)過此循環(huán)迭代過程直接得到了求解模態(tài)展開系數(shù)的方程,確保了重建方程系數(shù)矩陣的條件數(shù)最小。整個測點選擇流程如圖1所示。

        上述介紹的最優(yōu)測點選擇方法同樣適用于對某一頻段范圍內(nèi)的噪聲進行分析,僅需要在每一次的循環(huán)迭代過程中增加一個對角線元素求和的步驟,具體方法如下:根據(jù)式(9)可知,在每一次的循環(huán)迭代過程中,均有下式出現(xiàn):

        pr=Φc=ΦΦ+p0=

        USVHVS+UHp0=USS+UHp0

        (13)

        式(13)是由表面聲壓重建表面聲壓,因此矩陣T=USS+UH為單位矩陣,T的對角線元素與輻射體表面離散測點是一一對應的關(guān)系。為得到對最小奇異值最為敏感的測點,只需將最小奇異值置零。假設分析的對象為某個頻段范圍內(nèi)的噪聲,則在此頻段內(nèi)均勻選取一組離散的頻率點f1,f2, …fn,根據(jù)本文的最優(yōu)測點選擇方法得到一組矩陣T1,T2, …Tn,將每個矩陣的對應的S和S+中的最小奇異值置零,并將對應位置的對角線元素相加,得到一組新的對角線元素,則將這組新的對角線元素的最小值所對應的測點確定為此次迭代過程所要去除的測點。

        圖1 最優(yōu)測點選擇流程圖Fig.1 The process ofdetermining the optimum measurement points

        3 實驗研究

        為了驗證本文提出的最優(yōu)測點選擇方法的有效性,在消聲水池中進行了雙層圓柱殼體實驗。雙層圓柱殼體高度為2 m,半徑0.5 m,內(nèi)殼厚度0.008 m,外殼厚度0.002 m,近場水聽器陣列布置固定于旋轉(zhuǎn)裝置上,通過旋轉(zhuǎn)近場水聽器陣列達到測量四周近場水聲的目的,近場水聽器線陣由11個水聽器組成,水聽器間距0.22 m,線陣貼近外殼表面,距離外殼表面0.14 m,圓柱殼體內(nèi)部裝有激振器,采用頻率為168 Hz的單頻信號對圓柱殼體進行激勵。圖2為實驗裝置簡圖,圖3為實驗現(xiàn)場圖。

        圖2 實驗裝置示意圖Fig.2 The sketch of experimental device

        圖3 實驗現(xiàn)場Fig.3 The scene of experiment

        當近場水聽器陣列以20°為步進角度旋轉(zhuǎn)一周時,在圓柱全息面測得198個測點聲壓值,取其中66個點對這198測點進行重建,用本文提出的最優(yōu)測點選擇方法求得這66個測點的最優(yōu)分布方式,為便于觀察測點的分布情況,將圓柱面沿圓周方向展開成平面形式,如圖4所示。

        選取168 Hz為分析頻率,當激振器工作時測得圓柱全息面聲壓幅值如圖5所示。根據(jù)本文采用的最佳聲輻射模態(tài)階數(shù)選取方法,另外增加一組測點,根據(jù)圓柱的形狀特點,另外增加的測點可選中間的一行和一列,其示意圖如圖6所示。正則化方法對于逆問題的求解具有非常重要的作用,本文在對圓柱全息面聲壓值進行重建時均采用了Tikhonov正則化方法。

        圖4 最優(yōu)測點分布圖Fig.4 The optimum distribution of measurement points

        圖5 圓柱全息面聲壓測量值Fig.5 Measurement value of the pressure on cylindrical holographic surface

        圖6 增加測點示意圖Fig.6 The sketch of the increased measurement points

        重建中選取的聲輻射模態(tài)階數(shù)與增加測點的重建誤差之間的關(guān)系如圖7(a)所示,選取階數(shù)與所有測點的重建誤差之間的關(guān)系如圖7(b)所示。從圖7(a)與圖7(b)可觀察到二者具有相同的變化趨勢,重建誤差隨著選取的輻射模態(tài)階數(shù)的增加先減小后增大,圖7(a)在28階處重建誤差取得最小值14.57%,圖7(b)在29階處重建誤差取得最小值13.01%,這說明將這種基于測量的迭代方法引入到基于聲輻射模態(tài)的聲場重建中來確定聲輻射模態(tài)階數(shù)是有效的。此處確定的重建階數(shù)為28,其重建誤差為13.2%,重建效果如圖8所示。

        (a) 增加測點 (b)所有測點圖7 階數(shù)選取與重建誤差關(guān)系Fig.7 Relationship between number of modes and reconstruction error

        圖8 圓柱全息面聲壓重建值Fig.8 Reconstruction value of the pressure on cylindrical holographic surface

        為驗證本文中最優(yōu)測點選擇方法的有效性,在其他條件不變的情況下,采取下面兩種均勻布置的方式來重建圓柱全息面聲壓值,布置測點數(shù)目仍為66,其布置方式如圖9、圖10所示。

        圖9 均勻布置方式1Fig.9 The firsteven arrangement

        圖10 均勻布置方式2Fig.10 The second even arrangement

        兩種布置方式條件下選取的階數(shù)與所有測點的重建誤差之間的關(guān)系如圖11所示,從圖11中可知,均勻布置方式1條件下,重建誤差在46階處取得最小值64.32%,均勻布置方式2條件下,重建誤差在34階處取得最小值23.22%,重建效果如圖12所示。三種布置方式下圓柱全息面聲壓值重建誤差如表1所示。由圖11、12及表1可以看出,測點分布方式對最終的重建結(jié)果具有非常重要的影響,當測點按均勻布置方式1布置時,重建誤差達到64.32%,重建已經(jīng)失效;當測點按本文提出的方法進行布置時,重建誤差僅為13.20%,優(yōu)于另外兩種均勻布置方式。對比結(jié)果表明,利用本文提出的最優(yōu)測點選擇方法布置測點,能夠?qū)碗s結(jié)構(gòu)聲場進行有效的重建,提高了聲場重建精度。

        (a) 均勻布置方式1 (b) 均勻布置方式2圖11 階數(shù)選取與重建誤差關(guān)系Fig.11 Relationship between number of modes and reconstruction error

        (a) 均勻布置方式1 (b) 均勻布置方式2圖12 圓柱全息面聲壓重建值Fig.12 Reconstruction value of the pressure on cylindrical holographic surface

        最優(yōu)布置方式均勻布置方式1均勻布置方式213.20%64.32%23.22%

        上述過程利用全息測量面部分測點聲壓值重建了圓柱全息面所有測點的聲壓值,重建誤差在工程實際可接受范圍內(nèi),取得了較為理想的效果。求得圓柱全息面聲壓值后進而可以求解出聲場中任意一點的聲壓值,此計算過程本文不再敘述。

        4 結(jié) 論

        研究結(jié)果表明,利用聲輻射模態(tài)理論由輻射體表面少量測點重建聲場時,測點布置方式是決定聲場重建精度的一個關(guān)鍵因素。基于重建方程條件數(shù)最小準則,本文提出了一種基于聲輻射模態(tài)的最優(yōu)測點選擇方法,利用該方法布置測點,能夠?qū)β晥鲞M行有效的重建,重建效果優(yōu)于均勻布置方式,聲場重建精度明顯提高,并為實驗驗證。

        [ 1 ] MAYNARD J D, WILLIAMS E G, LEE Y. Near field acoustic holography Ⅰ. theory of generalized holography and the development of NAH[J]. J. Acoust. Soc. Am. , 1985, 78(4): 1395-1413.

        [ 2 ] WANG Zhaoxi, WU S F. Helmholtz equation-least-squares method for reconstructing the acoustic pressure field[J]. J. Acoust. Soc. Am. , 1997, 102(4): 2020-2032.

        [ 3 ] WU S F. On reconstruction of acoustic pressure fields using the Helmholtz equation least squares method[J]. J. Acoust. Soc. Am. , 2000, 107(5): 2511-2522.

        [ 4 ] 張海濱, 萬泉, 蔣偉康. HELS法在循環(huán)平穩(wěn)聲場全息重建中的理論與試驗研究[J]. 物理學報, 2009, 58(1):333-340. ZHANG Haibin, WAN Quan, JIANG Weikang. Theorical and experimental research on the HELS method in the cyclo stationary acoustic field[J]. Acta Physica Sinica, 2009, 58(1):333-340.

        [ 5 ] 姜哲. 聲輻射問題中的模態(tài)分析:Ⅲ.聲場重構(gòu)[J]. 聲學學報, 2005, 30(3): 242-248. JIANG Zhe. A modal analysis for the acoustic radiation problem: III. Reconstruction of acoustic fields [J].Acta Acustica, 2005, 30(3):242-248.

        [ 6 ] 聶永發(fā), 朱海潮. 利用源強密度聲輻射模態(tài)重建聲場[J]. 物理學報, 2014, 63(10):104303. NIE Yongfa, ZHU Haichao.Acoustic field reconstruction using source strength density acoustic radiation modes[J]. Acta Physica Sinica, 2014,63(10):104303.

        [ 7 ] BORGIOTTI G V. The power radiated by a vibrating body in an acoustic fluid and its determination from boundary measurements[J]. J. Acoust. Soc. Am. , 1990, 88(4): 1884-1893.

        [ 8 ] CUNEFARE K A. The design sensitivity and control of acoustic power radiated by three-dimensional structures[D]. Ph.D., The Pennsylvania State University, 1990.

        [ 9 ] ELLIOTT S J, JOHNSON M E. Radiation modes and the activecontrol of sound power[J]. J. Acoust. Soc. Am. , 1993, 94(4): 2194-2204.

        [10] PETERS H, KESSISSOGLOU N. Enforcing reciprocity in numerical analysis of acoustic radiation modes and sound power evaluation[J]. Journal of Computational Acoustics, 2012, 20(3): 5-33.

        [11] WU Haijun,JIANG Weikang, ZHANG Yilin. A method to compute the radiated sound power based on mapped acoustic radiation modes[J]. Journal of Computational Acoustics, 2014, 135(2): 679-692.

        [12] 杜向華,朱海潮,毛榮富. 利用聲輻射模態(tài)進行聲功率的靈敏度分析[J]. 振動與沖擊,2011,30(11):183-185. DU Xianghua, ZHU Haichao, MAO Rongfu. Acoustic sensitivity analysis of sound power using sound radiation modes[J]. Journal of Vibration and Shock, 2011, 30(11):183-185.

        [13] 袁國清, 姜哲. 基于聲輻射模態(tài)模型求解聲功率靈敏度[J]. 振動與沖擊,2009,28(8):109-112. YUAN Guoqing, JIANG Zhe. Method for solving sound power sensitivity based on the theory of acoustic radiationmodes[J]. Journal of Vibration and Shock, 2009, 28(8):109-112.

        [14] 宋艷華, 姜哲. 測量點位置對聲源識別的影響[J]. 振動與沖擊,2008,27(3):35-37. SONG Yanhua, JIANG Zhe. Effect of measurement points distribution on sound souce identification [J]. Journal of Vibration and Shock, 2008, 27(3):35-37.

        [15] KAMMER D C. Sensor placement for on-orbit model identification and correlation of large space structures[J]. J. Guid. Control Dyn. 1991,14: 251-259.

        [16] KIM B K, IH J G. On the reconstruction of the vibro-acoustic field over the surface enclosing an interior space using the boundary element method[J]. J. Acoust. Soc. Am. ,1996, 100(5): 3003-3016.

        [17] ZHANG Zhidong, VLAHOPOULOS N. A computational acoustic field reconstruction process based on an indirect boundary element formulation[J]. J. Acoust. Soc. Am. 2000,108 (5):2167-2178.

        [18] 雷宣揚,陳進,張桂才,等. 基于Helmholtz方程最小二乘法的聲場重構(gòu)[J]. 上海交通大學學報,2006,40(1):129-132. LEI Xuanyang, CHEN Jin, ZHANG Guicai, et al.The reconstruction of sound field using helmholtz equation-least squares method[J]. Journal of Shang Hai Jiao Tong University, 2006, 40(1):129-132.

        Optimization of measurement points in reconstruction of acoustic fieldbased on acoustic radiation modes

        SU Junbo, ZHU Haihao, MAO Rongfu, GUO Liang, SU Changwei

        (National Key Laboratory on Ship Vibration & Noise, Naval University of Engineering, Wuhan 430033, China)

        For an acoustic field reconstruction using the acoustic radiation mode theory, the distribution of measurement points is a key factor, especially, when the number of measurement points is small. In order to determine the optimal arrangement of measurement points, a method based on acoustic radiation modes was proposed. A loop iteration process was adopted in this proposed method. Singular value decomposition (SVD) of the matrix of acoustic radiation modes was employed in each loop iteration, and the measurement point that was the most sensitive to the minimum singular value was removed. Finally, those points that minimize the condition number of reconstruction equations were obtained. Furthermore, compared with the situation when measurement points distribute evenly, the test results showed that the acoustic field of sound sources can be reconstructed effectively and the precision of reconstruction is much better when measurement points distribute according to their optimal arrangement obtained with the proposed method.

        acoustic radiation modes; acoustic field reconstruction; optimal arrangement of measurement points

        國家自然科學基金(51305452)

        2015-09-28 修改稿收到日期:2016-01-14

        蘇俊博 男,博士生,1985年生

        朱海潮 男,教授,博士生導師,1963年生

        O423;O429

        A

        10.13465/j.cnki.jvs.2017.03.023

        猜你喜歡
        階數(shù)聲壓聲場
        基于嘴唇處的聲壓數(shù)據(jù)確定人體聲道半徑
        關(guān)于無窮小階數(shù)的幾點注記
        確定有限級數(shù)解的階數(shù)上界的一種n階展開方法
        基于BIM的鐵路車站聲場仿真分析研究
        車輛結(jié)構(gòu)噪聲傳遞特性及其峰值噪聲成因的分析
        汽車工程(2018年12期)2019-01-29 06:46:36
        探尋360°全聲場發(fā)聲門道
        基于GIS內(nèi)部放電聲壓特性進行閃絡定位的研究
        電測與儀表(2016年9期)2016-04-12 00:30:02
        一種新的多址信道有效階數(shù)估計算法*
        關(guān)于動態(tài)電路階數(shù)的討論
        板結(jié)構(gòu)-聲場耦合分析的FE-LSPIM/FE法
        91中文在线九色视频| 亚洲中文字幕午夜精品| 亚洲人成在久久综合网站| 亚洲国产精品综合久久网络| 免费看又色又爽又黄的国产软件 | 无码精品人妻一区二区三区av | 88久久精品无码一区二区毛片| 日韩国产一区| 中文字幕日产人妻久久| 日韩精品成人一区二区在线观看| 中国少妇和黑人做爰视频| 日日麻批视频免费播放器| 国产一区二区视频免费在线观看| 欧美日韩午夜群交多人轮换| 影音先锋男人站| 99热久久这里只精品国产www| 日韩精品网| 亚洲青青草视频在线播放| 蜜桃视频羞羞在线观看| 日本va欧美va精品发布| 女人被狂躁c到高潮视频| 久无码久无码av无码| 人妻少妇精品无码专区app| 青青久久精品一本一区人人 | 一区五码在线| 国产自拍精品在线视频| 国产av无毛无遮挡网站| 国产亚洲精品av一区| 亚洲av成人噜噜无码网站| 女人被做到高潮免费视频| 久久久久国产精品四虎| 日韩美女人妻一区二区三区| 亚洲人不卡另类日韩精品| 亚洲av无码一区二区三区观看| 亚洲最新版无码AV| 色妞一区二区三区免费视频| 日本不卡高字幕在线2019| 国产麻豆精品久久一二三| 国产精品丝袜美女在线观看| 日本免费播放一区二区| 日韩日韩日韩日韩日韩日韩|