999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

非均勻介質條件下礦震震波三維傳播模型構建及其應用

2014-09-11 08:45:20賈寶新李國臻
煤炭學報 2014年2期
關鍵詞:模型

賈寶新,趙 培,姜 明,李國臻,王 帥

(1.遼寧工程技術大學 土木與交通學院,遼寧 阜新 123000;2.遼寧工程技術大學 沖擊地壓研究院,遼寧 阜新 123000)

非均勻介質條件下礦震震波三維傳播模型構建及其應用

賈寶新1,2,趙 培1,姜 明1,李國臻2,王 帥1

(1.遼寧工程技術大學 土木與交通學院,遼寧 阜新 123000;2.遼寧工程技術大學 沖擊地壓研究院,遼寧 阜新 123000)

為了研究速度模型對礦震定位誤差的影響,基于惠更斯原理建立了礦震波在層狀非均勻介質中的三維傳播模型,通過現場爆破振動試驗,采用礦山微震定位監測系統進行監測。對礦震震波的正演建立了礦震震波傳播的三維波動方程,并給出了礦震波在介質分界面上傳播速度計算公式,由此對經典線性定位方法進行合理改進,使得基于均勻介質的線性定位方法適合于非均勻分層介質。將理論推導結果應用于礦山微震監測定位系統中,進行實際礦震定位計算。試驗表明,速度模型對定位的影響較大,運用所建立的三維傳播模型應用于礦震定位,有效提高了定位精度。

礦震震波;非均勻介質;三維傳播模型;定位

對于礦山地震定位的研究,國內外已經有了較多的研究成果。由于礦震震源淺,速度模型可以簡化,因此Geiger法及其各改進方法已被應用于微震定位的研究[1-2],在有足夠臺站數的情況下值得應用。非線性方法中的Powell直接搜索法,具有不需要求目標函數的偏導數或逆矩陣,且對初值要求不高的特點,也已經應用于微震定位[3-4]。

但是在實際情況中,由于臺網數量少、震源淺、信號的傳播距離短,而且信號能量小,衰減較快,臺站數量達不到產生超越方程組的個數,使得一些定位方法不能應用,S波讀取困難,這些因素給精確定位帶來了困難[5]。因此對于礦震定位精度的提高還有待進一步研究,本文旨在通過建立礦震震波的傳播模型,并將其應用于礦震的定位中,以提高礦震定位精度。

在震波傳播模型的研究中,吳永國等[6]基于惠更斯原理,利用下行波模擬炮點所激發出的地震波在地質介質中傳播過程,實現了單程波動方程模擬疊前正演。趙愛華等[7]基于惠更斯原理和費馬原理研究了求取地震波走時及其反射波射線路徑的新方法。趙國敏等[8]根據地震波射線追蹤時的需要,基于惠更斯原理,研究了確定射線出射角與衰減角的有效方法。阮穎錚[9]基于復惠更斯原理,對瞬變地震波反射場、地震波場的輻射線展開等綜合地震波反射場方面做了深入研究。郭繼亮等[10]根據巖石物理學中的臨界孔隙度模型構建了簡潔的均勻彈性流體飽和孔隙介質模型,借此進行地震波傳播的研究。張麗琴等[11]根據地震波在離散介質中的傳播理論,采用偽譜法對三維離散介質中地震波的傳播進行了正演數值模擬。范春方等[12]以小擾動為前提條件,得出了不同震源條件下地震波在 Kelvin-Voigt均勻黏彈性介質中的 Kelvin-Voigt 均勻黏彈性三階波動方程的平面波解。巴晶[13]基于Biot理論,考慮了多孔介質中的固體骨架的松弛特征,對復雜多孔介質中的地震波傳播機理進行了研究。鞏思園等[14]構建了礦井尺度的微震監測系統異向波速模型,通過現場實驗證明其震源定位精度得到大幅提高。隨著科技的發展,越來越多的理論模型被構建并予以研究,但是對于礦震震波的研究還有待進一步完善。

1 礦震震波傳播模型的建立

1.1 利用惠更斯原理對分層介質中波陣面方程的建立

圖1為惠更斯原理示意。圖中,Ω1,Ω2表示兩種不同介質;O為震源;x軸為介質分界面;z軸為過震源且與介質分界面垂直的方向;波陣面恰好到達介質分界面時為初始狀態,此時波陣面與介質分界面切點為P;考察時刻波陣面與介質分界面的交點為T,與z軸交點為Q;θ為OT與z軸的夾角;μ為參數方程的參數,其值對應PT之間任意一點M,M點發射的子波在介質Ω2實際到達W點。另外,OM與初始狀態波陣面交點為N,OT與初始狀態波陣面交點為S。r為震源到介質分界面距離,也等于初始時刻球面波半徑;R為考察時刻波陣面半徑。

根根《礦震震波傳播規律的三維模型及其應用》[15]中的結論,礦震波在某一種介質中的傳播速度為

(1)

式中,E為介質的彈性模量;μ為介質的泊松比;ρ為介質的線密度。

對于Ω1,有

(2)

對于Ω2,有

(3)

假設此時震源O出發的球面波已經擴散至T點,此時球面波半徑達到R,OT與OP交角為θ(圖1)。

根據惠更斯原理,介質中波陣面(波前)上的各點都可以看作是發射子波的波源,其后任一時刻這些子波的包跡就是新的波陣面。

由圖1可知,波陣面由S點到達T點時,水平方向上陣面由P點傳播到Q點,所經歷時間設為t,則

(4)

現考察PT上任一點M,設OM與OP交角為μθ(0<μ<1),當μ=0時,M點與P點重合;當μ=1時,M點與T點重合。

(5)

此后當波陣面水平方向傳播到Q點時,以M點為子波波源的子波水平方向傳播到W點,所經歷時間為tμ,則

(6)

以P點為坐標原點,水平方向為z軸,豎直方向為x軸,建立坐標系,則W點的坐標可表示為(xw,zw),其中

(7)

其中,r,θ,v1,v2為常量。F(x,z)是以μ為參數的參數方程。其表示的曲線就是波穿過介質分界面后形成的新波陣面的形狀,消去μ可得F(x,z)的隱式方程為

(8)

1.2 利用波陣面方程的參數形式對波陣面的繪制

為了直觀看到波陣面所表示的曲線形狀,取特殊的θ,r,v1,v2,利用以μ為參數的參數方程繪制波陣面圖形,并與波在均勻介質中的球狀波陣面對比。

(1)取r=10 m,v1=a,v2=2a(即v1=0.5v2)。

當θ=π/6時,參數方程為

(9)

當θ=π/4時,參數方程為

(10)

當θ=π/3時,參數方程為

(11)

將上述3種θ的波陣面繪制在同一幅圖中,并與均勻介質中的球形波陣面對比,如圖2所示。

圖2 Ω1中波速小于Ω2中波速的波陣面分布Fig.2 Distribution of wavefront when wave velocity in Ω1 less than that in Ω2

從圖2可以看出,當Ω1中的波速v1小于Ω2中的波速v2時,Ω2中的波陣面更尖銳,比均勻介質中的球形波陣面超前。因此該方程的描述符合實際情況。

(2)取r=10 m,v1=2a,v2=a(即v1=2v2)。

當θ=π/6時,參數方程為

(12)

當θ=π/4時,參數方程為

(13)

當θ=π/3時,參數方程為

(14)

將上述3種θ的波陣面畫在同一幅圖中,并與均勻介質中的球形波陣面對比,如圖3所示。

圖3 Ω1中波速大于Ω2中波速的波陣面分布Fig.3 Distribution of wavefront when wave velocity of Ω1 larger than that of Ω2

從圖3可以看出,當v1大于v2時,Ω2中的波陣面更平緩,比均勻介質中的球形波陣面落后,因此該方程的描述符合實際情況。

(3)當v1=v2時,參數方程為

(15)

式(15)表示的曲線為圓的一部分,并與均勻介質中相同θ值對應的球形波陣面重合。即當分界面兩側的介質為同一介質時,穿過分界面后礦震波的波陣面轉化為均勻介質中的球形波陣面。這與實際情況是吻合的。

1.3 三維波陣面方程

(16)

相應的參數方程為

(17)

式中,ψ為OM在x,y面上的投影與x軸的夾角。

1.4 經過介質變化面后任意點的波動方程

對于穿過介質分界面后的任意一點(x0,y0,z0),可由x0或y0求出參數u0,代入式(17)可得該點沿z軸方向傳播的子波的傳播距離為

(18)

根根《礦震震波傳播規律的三維模型及其應用》[15]中的結論,有

(19)

令r=0,得子波的初振動為

(20)

(21)

這就是震波經過介質后在某點形成的子波的波動方程。根據惠更斯原理,波陣面上某點發射的子波仍為球面波,即(x0,y0,z0)處子波的波陣面關于(x0,y0,z0)呈中心對稱,所以這個方程對于x,y,z三個方向的波動均適用。

1.5 波陣面在介質變化面上的傳播速度

研究波陣面在介質分界面上覆蓋范圍的傳播速度,以利于在分界面處布設監測站來獲得參數θ,r。設某時球面波在介質分界面上的覆蓋范圍為2h,則

(22)

兩邊對t求導得

(23)

因為

(24)

所以

(25)

2 傳播模型在礦震定位中的應用

2.1 定位方法

通過情況下,不小于4個臺站的定位,稱為多臺定位。設礦震震源點E的礦區坐標為(x0,y0,z0),發震時刻為t0,假定P波在煤巖體介質中以常速度vp傳播,則礦震震源與第i個傳感器之間的走時方程[16-17]為

(26)

其中,(xi,yi,zi)為第i個傳感器的測量坐標;ti為監測到時;m為接收到信號的傳感器的個數;(x0,y0,z0,t0)為所要求的微震源的時空參數。用一般方法處理這個非線性問題比較困難,這就需要尋找一線性系統來代替,基本函數就是描述實際用來替代的方程組,由式(26)經代數變換得到。

如果用第i個測點的走時方程減去第k個測點的走時方程,將會得到SW-GBM算式的所有函數形式,即

2(xi-xk)x0+2(yi-yk)y0+2(zi-zk)z0-

(27)

2(xi-xi-1)x0+2(yi-yi-1)y0+2(zi-zi-1)z0-

(28)

式(28)以矩陣形式表示為

(29)

考慮本定位問題中m=n=4,通常這種情況是不能用走時方程定位的。但其中的一種特殊情況可以,即4個臺站處于同一水平面上。這時,礦震震源點E到4個臺站的鉛垂距離相等,相當于消去一個未知數。反映在式(28)上就是3個未知數,m-1=3個方程,方程組可解,但要注意系數矩陣A33為病態的可能。編程中要使用可以求解病態方程組的算法[18]。

通過爆破振動試驗對上述理論進行驗證。試驗地點為北京昊華集團木城澗煤礦,監測系統為遼寧工程技術大學研發的礦山微震監測定位系統[19-21]。表1為定位結果比較。由表1可知,第1組誤差最大,第3組誤差最小;4組結果中鉛垂向都較其他兩方向誤差要大,原因可能在于z值是回代求解,這也驗證了定位中震源深度誤差較大的共識。

表1實際震源與計算結果對比
Table1Correlationoftheactualsourceandcalculationresults

序號發震時刻爆破震源/m計算震源/m誤差/mx0T=-18000x0C=-18005 0373Δx=5 037312010-04-06T10:06:02y0T=4421000y0C=4421003 0013Δy=3 0013z0T=810z0C=802 9532Δz=-7 0468x0T=-16000x0C=-16004 0079Δx=4 007922010-04-06T11:12:41y0T=4420670y0C=4420673 0017Δy=3 0017z0T=750z0C=756 0576Δz=6 5760x0T=-16205x0C=-16201 9876Δx=-3 012432010-04-06T14:32:41y0T=4420550y0C=4420547 9917Δy=2 0083z0T=800z0C=794 7528Δz=-5 2472x0T=-17050x0C=-17053 0035Δx=3 003542010-04-06T15:20:33y0T=4420350y0C=4420345 9560Δy=-4 0440z0T=700z0C=693 6969Δz=6 3031

2.2 應用傳播模型改進后的定位方法

走時方程為

SW-GBM算式的所有函數形式為

2(xi-xk)x0+2(yi-yk)y0+2(zi-zk)z0-

(31)

選取應以每個傳感器貢獻的信息均等為準則,式(32)給出了一種特殊集合,它的優點是能夠抵消測量得到的各個傳感器的坐標和到時的誤差。

2(xi-xi-1)x0+2(yi-yi-1)y0+2(zi-zi-1)z0-

(32)

基于式(32),用VB編程進行求解。通過爆破振動試驗對上述理論進行驗證。表2為定位結果的比較。由表2可見,定位誤差較表1所示誤差相比有所降低,說明定位精度有了明顯提高。

3 結 論

(1)基于惠更斯原理,建立了礦震震波傳播的三維波動方程:

(r2+x2+y2)(v2cosθ)2=(v2r-zv1cosθ)2

(2)基于礦震震波的正演結果,以及得到的介質分界面上波速,對經典線性定位方法進行合理改進,使得基于均勻介質的線性定位方法適合于非均勻分層介質。

表2實際震源與計算結果對比
Table2Correlationoftheactualsourceandcalculationresults

序號發震時刻爆破震源/m計算震源/m誤差/mx0T=-18000x0C=-18002 1245Δx=2 124512010-06-06T9:01:05y0T=4421000y0C=4421002 6618Δy=2 6618z0T=810z0C=805 9578Δz=-4 0422x0T=-16000x0C=-16001 2517Δx=1 251722010-06-06T10:21:32y0T=4420670y0C=4420672 8142Δy=2 8142z0T=750z0C=753 1625Δz=3 1625x0T=-16205x0C=-16202 8645Δx=-2 135532010-06-06T14:20:11y0T=4420550y0C=4420548 8411Δy=-1 1589z0T=800z0C=796 5489Δz=-3 4511x0T=-17050x0C=-17052 4123Δx=2 412342010-06-06T15:15:26y0T=4420350y0C=4420348 1625Δy=-1 8375z0T=700z0C=695 7748Δz=-4 2252

(3)將理論推導所建立的礦震震波傳播的三維模型及其他推導結果應用于遼寧工程技術大學研發的礦山微震監測定位系統中,進行實際礦震定位計算,試驗證明理論推導結果有效地提高了定位精度。

對于礦震定位的研究要同時將速度模型的研究、到時拾取、監測臺站的空間布局進行綜合分析,才能更有效提高礦震定位精度。這也是今后重點的研究工作。

[1] 胡新亮,馬勝利,高景春,等.相對定位方法在非完整巖體聲發射定位中的應用[J].巖石力學與工程學報,2004,23(2):277-283. Hu Xinliang,Ma Shengli,Gao Jingchun,et al.Location of acoustic emissions in non-integral rock using relative locating method[J].Chinese Journal of Rock Mechanics and Engineering,2004,23(2):277-283.

[2] 郭 飚,劉啟元,陳九輝,等.首都圈數字地震臺網的微震定位實驗[J].地震地質,2002,24(3):453-460. Guo Biao,Liu Qiyuan,Chen Jiuhui,et al.Test of epicenter determination of microearthquakes recorded by the digital seismicnetwork in capital circle[J].Seismology and Geology,2002,24(3):453-460.

[3] 汪素云,許忠淮,俞言祥,等.北京及鄰區現代微震重新定位及其構造含義[J].中國地震,1995,11(3):222-230. Wang Suyun,Xu Zhonghuai,Yu Yanxiang,et al.Relocation of microearthquakes in Beijing and its neighbouring areas and its tectonic implications[J].Earthquake Research in China,1995,11(3):222-230.

[4] 汪素云,許忠淮,俞言祥,等.北京西北地區現代微震重新定位[J].地震學報,1994,16(1):24-31. Wang Suyun,Xu Zhonghuai,Yu Yanxiang,et al.Relocation of microearthquakes in the northwest region of Beijing[J].Acta Seismologica Sinica,1994,16(1):24-31.

[5] 梁寶霞.礦震信號識別和定位及其在木城澗煤礦的應用[D].阜新:遼寧工程技術大學,2007. Liang Baoxia.Thesignal identification and localization of rockburst and their application at Muchengjian Colliery[D].Fuxin:Liaoning Technical University,2007.

[6] 吳永國,賀振華,黃德濟.基于惠更斯原理的波動方程共炮點道集地震正演[J].成都理工大學學報(自然科學版),2008,35(3):275-278. Wu Yongguo,He Zhenhua,Huang Deji.Wave-field modeling of common shot gatherrecords based on Huygens principle[J].Journal of Chengdu University of Technology(Science & Technology Edition),2008,35(3):275-278.

[7] 趙愛華,丁志峰,孫為國,等.復雜介質地震定位中震源軌跡的計算[J].地球物理學報,2008,51(4):1188-1195. Zhao Aihua,Ding Zhifeng,Sun Weiguo,et al.Calculation of focal loci for earthquake location in complex media[J].Chinese Journal of Geophysics,2008,51(4):1188-1195.

[8] 趙國敏,張 中,滕吉文,等.線性連續變化非彈性介質中地震波射線軌跡與走時研究[J].華北地震科學,1994,12(2):49-54. Zhao Guomin,Zhang Zhong,Teng Jiwen,et al.A study of raypath and traveltime of seismic waves in linear continuous varyingvis coelastic media[J].North China Earthquake Sciences,1994,12(2):49-54.

[9] 阮穎錚.用復惠更斯原理綜合地震波反射場[J].電子科技大學學報,1990,19(4):336-341. Ruan Yingzheng.Synthesis of reflected seismic wavefield by complex Huyghens’ principle[J].Journal of University of Electronic Science and Technology of China,1990,19(4):336-341.

[10] 郭繼亮,牛濱華,孫春巖,等.基于臨界孔隙度模型的地震波傳播[J].地球物理學報,2012,55(11):3813-3820. Guo Jiliang,Niu Binhua,Sun Chunyan,et al.Seismic wave propagation based on critical porosity model[J].Chinese J.Geophys.,2010,55(11):3813-3820.

[11] 張麗琴,於文輝,王家映.離散介質中地震波傳播的數值模擬[J].西北地震學報,2004,26(3):89-92. Zhang Liqin,Yu Wenhui,Wang Jiaying.Numerical simulation of seismic waves propagation in discrete media[J].Northwestern Seismological Journal,2004,26(3):89-92.

[12] 范春方.Kelvin-Voigt均勻黏彈性介質中傳播的地震波[J].中國科學D輯:地球科學,2005,35(10):957-962. Fan Chunfang.Seismic wave propagation in homogeneous viscoelastic media[J].Science in China Ser.D:Earth Sciences,2005,35(10):957-962.

[13] 巴 晶.復雜多孔介質中的地震波傳播機理研究[D].北京:清華大學,2008. Ba Jing.Research on the seismic wave propagation in complex porous media[D].Beijing:Tsinghua University,2008.

[14] 鞏思園,竇林名,馬小平,等.煤礦礦震定位中異向波速模型的構建與求解[J].地球物理學報,2012,55(5):1757-1763. Gong Siyuan,Dou Linming,Ma Xiaoping,et al.Study on the construction and solution technique of anisotropic velocity model in the location of coal mine tremor[J].Chinese J.Geophys.,2012,55(5):1757-1763.

[15] 潘一山,賈寶新,王 帥,等.礦震震波傳播規律的三維模型及其應用[J].煤炭學報,2012,37(11):1810-1814. Pan Yishan,Jia Baoxin,Wang Shuai,et al.Three-dimensional of model and its application mine seismic wave propagation[J].Journal of China Coal Society,2012,37(11):1810-1814.

[16] 王煥義.巖體微震事件的精確定位研究[J].工程爆破,2001,7(3):5-8. Wang Huanyi.Study on precise localization of micro-seismic events in a rock-mass[J].Engineering Blasting,2001,7(3):5-8.

[17] 周長發.科學與工程數值算法(Visual Basic 版)[M].北京:清華大學出版社,2002. Zhou Changfa.Science and engineering numerical algorithm(Visual Basic Version)[M].Beijing:Tsinghua University Press,2002.

[18] 丁紅旗,李國臻,賈寶新.微震振源激振模型及振動周期與震級的關系[J].遼寧工程技術大學學報(自然科學版),2009,28(3):352-354. Ding Hongqi,Li Guozhen,Jia Baoxin.The microseism source excitation model and the relation between the magnitude and vibration cycle[J].Journal of Liaoning Technical University(Natural Science),2009,28(3):352-354.

[19] 丁紅旗,李國臻,賈寶新,等.微震振源激振模型及震級計算公式的建立[J].科學技術與工程,2009,9(14):4130-4133. Ding Hongqi,Li Guozhen,Jia Baoxin,et al.An excitation model of the microseismicity source and a formula to calculate its vibration magnitude[J].Science Technology and Engineering,2009,9(14):4130-4133.

[20] 潘一山,趙揚鋒,官福海,等.礦震監測定位系統的研究及應用[J].巖石力學與工程學報,2007,26(5):1002-1071. Pan Yishan,Zhao Yangfeng,Guan Fuhai,et al.Study on rockburst monitoring and orientation system and its application[J].Chinese Journal of Rock Mechanics and Engineering,2007,26(5):1002-1071.

[21] 賈寶新,李國臻.礦山地震監測臺站的空間分布研究與應用[J].煤炭學報,2010,35(12):2045-2048. Jia Baoxin,Li Guozhen.The research and application for spatial distribution of mines seismic monitoring stations[J].Journal of China Coal Society,2010,35(12):2045-2048.

Three-dimensionalpropagationmodelbuildinginheterogeneousmediumofmineearthquakeshockwaveanditsapplication

JIA Bao-xin1,2,ZHAO Pei1,JIANG Ming1,LI Guo-zhen2,WANG Shuai1

(1.InstituteofCivilandTransportation,LiaoningTechnicalUniversity,Fuxin123000,China;2.ResearchInstituteofRockBurst,LiaoningTechnicalUniversity,Fuxin123000,China)

In order to study the influence of velocity model on the deviation of shock bump earthquake location,three dimensional propagation model of mine quake waves in the layered inhomogeneous media was built based on Huygens principle.Through blasting vibration tests,monitoring was made by microseismic system of mine.The three-dimensional wave equation was built for seismic waves forward,and propagation velocity’s calculation formula on the dielectric interface was also given,on which the classical linear positioning method was reasonably improved so that the linear positioning method based on uniform medium suits for inhomogeneous medium.The actual calculation for shock bump can be made through the application of deduction results to mine microseismic monitoring and positioning system.The results show that velocity model greatly affects positioning.By using 3D propagation model in the shock bump positioning,location accuracy is effectively improve.

mine earthquake shock wave;heterogeneous medium;three dimensional propagation model;positioning

10.13225/j.cnki.jccs.2013.2014

國家重點基礎研究發展計劃(973)資助項目(2010CB226803);遼寧省自然科學基金資助項目(201204414);遼寧省教育廳自然科學資助項目(L2012107)

賈寶新(1978—),男,遼寧撫順人,副教授。E-mail:jbx_811010@126.com

TD324

A

0253-9993(2014)02-0364-07

賈寶新,趙 培,姜 明,等.非均勻介質條件下礦震震波三維傳播模型構建及其應用[J].煤炭學報,2014,39(2):364-370.

Jia Baoxin,Zhao Pei,Jiang Ming,et al.Three-dimensional propagation model building in heterogeneous medium of mine earthquake shock wave and its application[J].Journal of China Coal Society,2014,39(2):364-370.doi:10.13225/j.cnki.jccs.2013.2014

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 欧美国产日产一区二区| 精品欧美一区二区三区久久久| 日韩视频免费| www.av男人.com| 久996视频精品免费观看| 国产精品丝袜视频| 91区国产福利在线观看午夜| 九九这里只有精品视频| 久久黄色免费电影| 日韩专区欧美| 国产亚洲日韩av在线| 欧美午夜在线观看| 亚洲美女久久| 欧美激情视频一区二区三区免费| 白浆视频在线观看| AV网站中文| 一本大道香蕉中文日本不卡高清二区| 亚洲人成在线免费观看| 日韩精品成人在线| 漂亮人妻被中出中文字幕久久| 任我操在线视频| 成人午夜视频在线| 久久久久久国产精品mv| 国产精品开放后亚洲| 91青草视频| 露脸国产精品自产在线播| 亚洲区第一页| 国产乱子精品一区二区在线观看| 色偷偷男人的天堂亚洲av| 亚洲国产欧美中日韩成人综合视频| 免费一级毛片在线播放傲雪网| 直接黄91麻豆网站| 在线看片免费人成视久网下载| 日韩一区二区在线电影| 99re精彩视频| 国产无码性爱一区二区三区| 综合久久久久久久综合网| 亚洲va视频| 一区二区三区四区精品视频 | 99热这里只有精品久久免费| 综合色亚洲| 亚洲激情99| 69精品在线观看| 亚洲另类国产欧美一区二区| 亚洲AV色香蕉一区二区| 亚洲男人在线| 国产乱子伦视频在线播放| 免费A∨中文乱码专区| 试看120秒男女啪啪免费| 国产白浆视频| 青青草原国产精品啪啪视频| 亚洲高清在线天堂精品| 91色爱欧美精品www| 98精品全国免费观看视频| 日本免费福利视频| 精品少妇人妻一区二区| 国产丝袜无码一区二区视频| 毛片免费试看| 久久综合丝袜日本网| 黄色三级毛片网站| 欧美在线视频不卡第一页| 成人第一页| 青青国产在线| 欧美区国产区| 国内精品九九久久久精品| 一级毛片在线播放| 一级毛片中文字幕| 亚洲国产精品日韩av专区| 欧美不卡二区| 人人澡人人爽欧美一区| a级毛片免费网站| 亚洲黄网视频| 四虎在线高清无码| 91欧美在线| 九九久久99精品| 久久99精品国产麻豆宅宅| 尤物视频一区| 国产成在线观看免费视频| 国产一级小视频| 天天色综网| 国产综合亚洲欧洲区精品无码| 亚洲欧美人成电影在线观看|