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

不同震源參數下的圍巖PGV分布規律及預測研究

2020-09-10 07:22:44王璽李施慶李秋濤程力劉煥新
黃金 2020年10期

王璽 李施慶 李秋濤 程力 劉煥新

摘要:為全面了解不同震源機制下地震波傳播和分布規律及其對巷道圍巖的影響,基于現場監測數據研究了不同震源參數下PGV的分布規律。研究結果表明,震源距離、矩震級及斷層滑移角都會對PGV的分布產生不同程度的影響。通過比較經驗公式計算PGV值與實測值間的擬合相關系數,判斷各預測方法適用性。研究發現Sadov預測公式對爆破事件的PGV預測值與實測結果相關度更高。依據研究結果,計算了動態震源機制下PGV空間分布,并依此給出了不同震源位置和震級條件下的圍巖支護建議。

關鍵詞:金屬礦山;震源;震動監測;PGV;預測公式;擬合相關系數

中圖分類號:TD853.34文章編號:1001-1277(2020)10-0029-08

文獻標志碼:Adoi:10.11792/hj20201006

引 言

隨著地表淺部資源的開采殆盡,金屬礦開采開始轉向深部。而由于地下金屬礦復雜的開采技術條件及多變的礦巖受力狀態,其采場內圍巖會產生不同程度的開采擾動,同時也導致地下金屬礦采場中存在多種高危的安全隱患[1]。在地下金屬礦開采過程中,若對礦山井下強擾動的影響范圍、具體位置、能量大小等缺乏準確的了解和認知,礦山工人與生產設備將持續處于極其危險且頻率高發的災害威脅中,并有可能引發重大的財產損失與人員傷亡,最終帶來難以估量的損失。

在礦山震動的研究工作中,預測震動強度對震動安全的判據有著重要意義。震動強度(I)的預測方法依據表征震動強度物理量的不同而對應不同的形式,通常各類預測方法都符合以下函數形式[2]:

I=f(m,n,k)(1)

式中:m為震源;n為傳播途徑;k為儀器特性。

通過公式(1)發現,影響震動強度的因素多樣且復雜,在計算過程中需要對多因素的影響進行綜合考量[3-4]。由于通過多變量確認函數形式的難度很高,該過程一般會選擇某些獨立變量作為常數,在此基礎上對其他變量進行研究,找到它們之間存在的函數關系。運用公式(1)及現場監測所得的震動相關數據能夠對質點峰值震速、加速度、位移等參數進行預測,再由預測結果對震動強度進行評判。實際工作中運用最廣泛的是對PGV(Peak Ground Velocity,PGV)進行預測。但是,在運用PGV經驗公式對實測數據進行回歸分析時,計算得出的回歸系數離散性往往較大。所以經驗公式在實際運用時會出現較多問題和局限性。為此,研究者結合現有經驗將公式進行變形,得出了很多適用于各類情況的相關公式。例如:陳壽如等[5]在計算露天礦爆破產生的震動模型時加入了高差這一變量,得到了新的變形公式;韓子榮[6]基于對主頻與結構自震頻率影響的考慮得出了折合速度公式。在預測PGV時除運用傳統線性回歸方法以外,近些年研究者研發出了不少新的預測方法,比較典型的有:王民壽等[7]通過對雙隨機變量進行回歸分析來預測震速;徐全軍等[8-9]通過神經網絡對震動峰值進行預測;黃光球等[10]通過遺傳規劃法對震動峰值進行預測。

綜合上述研究,針對震動事件產生的應力波,采礦中經常使用PGV描述圍巖在應力波作用下的能量平衡狀態(能量需求與能量吸收能力)和動態變化,由此對圍巖的穩定性和完整性作出評估[11]。因此,本次研究基于現場震動監測數據,對礦震數據在震源參數影響下的分布特征進行研究;并選取某礦震事件的震源參數,運用國內外已有的PGV計算公式進行預測計算;通過研究計算結果與實測值間的擬合相關系數,判斷不同的計算公式在該事件PGV預測的適用性;最后運用優選計算公式的預測結果給出巖體支護建議。

1 圍巖震動測定與分析

1.1 震源參數與PGV測定

山東金洲礦業集團有限公司金青頂礦區(下稱“金青頂礦區”)坐落在黃壘河南岸的金青頂上,位于山東省乳山市下初鎮。金青頂礦區整體屬于硬巖礦山,最大開采深度已超過1 000 m,大部分采場和巷道布置集中、開挖空間較小,且眾多采場與井巷工程的開挖空間結構參數和圍巖性質不盡相同,這些因素導致該礦區的震動事件錯綜復雜,嚴重影響礦區范圍內礦震監測和分析工作。

為解決該問題,需要對PGV在各項震源參數下的分布情況進行研究。利用局域微震監測系統(見圖1)對該礦區-1 145 m中段E8、E9、E10及E11穿脈內的活動事件進行了跟蹤記錄,記錄信息包括事件時間、矩震級、實測PGV、傳感器距震源距離及相對位置等參數。將記錄的監測數據與震源特征參數繪制成圖,即可得到PGV與相關震源參數的關系圖。本文分別討論了震源距離、矩震級及斷層滑移角對PGV分布的影響情況。

1.2 震源參數對PGV分布影響

1)震源距離對PGV的影響。隨機抽取8次震級不一的事件,繪制不同矩震級(Mw)下震源距離(d)與PGV(v)的關系,見圖2。

從圖2可以看出:隨著震源距離的增大,震動的峰值速度都呈指數型下降趨勢,說明傳播過程中,巖石介質將對PGV產生明顯的削弱作用。此外,PGV的衰減程度與震源距離也有著密切的聯系,當震源距離d<300 m時,事件8(Mw=0.542 397 9)PGV隨震源距離的衰減率為k1=4.776 784×10-6;當300 m≤d≤500 m時,衰減率k2=2.772 07×10-7;當d>500 m時,衰減率k3=5.883 29×10-8。由此可見,整體上PGV在傳播過程中的衰減作用受震源距離的影響,在距震源距離較小時,震源距離對PGV的衰減影響較大;隨著震源距離增大,PGV的衰減受震源距離的影響作用逐漸弱化。

2)矩震級對PGV的影響。矩震級是利用地震矩的大小確定震級,它由地震斷層的破裂面積、平均錯動量及巖石剪切模量的乘積來確定,是一個描述震動發生時力學強度的物理量。基于監測數據,取震源距離為100~200 m、200~300 m、300~400 m、400~500 m及500 m以上的5種情況下矩震級和PGV數值,分別繪制PGV(v)相對于矩震級(Mw)大小的散點分布圖,見圖3。

從圖3可以看出:矩震級增大時,PGV值高的散點數量增多,整體上PGV呈現增大趨勢;分析發現其原因可能是,隨著矩震級的增大,斷層的長度等尺寸參數變大,使得PGV隨矩震級增大也顯著變大[12]。但是,由PGV-Mw變化趨勢公式(見表1)可知:在不同震源距離下,PGV隨矩震級增大的程度有所差異,在近震源情況下,PGV受矩震級影響較大,隨矩震級呈近指數趨勢增長;但當監測點逐漸遠離震源時,PGV受矩震級的影響減小,增大趨勢逐漸減弱。由此說明,矩震級對PGV的分布有較大影響,且該作用效果隨震源距離的由近至遠逐漸衰減。

3)斷層滑移角對PGV的影響。在金屬礦山中,斷層滑移角通常指礦體上下盤相對運動的方向(斷層滑移方向)與斷層的走向之間的夾角。為探究不

同斷層滑移角對PGV分布的影響,在監測所得的數據中取滑移角(αr)為0°、45°、90°、-45°和-90°的5次震動事件,分別繪制PGV關于震源距離變化的曲線,見圖4。

從圖4可以看出:不同斷層滑移角下隨著震源距離的增大,所監測到的PGV值都在逐漸減小,整體上呈現出隨震源距離增大而衰減的趨勢,但在相同震源距離下,不同滑移角下的PGV值大小存在差異。

(1)整體上看,走滑斷層(αr=0°)下的PGV值最大,且分布較為規律,其次是正斷層(αr=45°)、垂直斷層(αr=±90°),而逆斷層(αr=-45°)下產生的PGV值最小。

(2)從曲線的走向和趨勢上看,αr=-45°、αr=90°和αr=-90°的3條曲線較為吻合,呈現的規律性相似,特別是當d>600 m時,3條曲線幾乎重合,說明逆斷層和垂直斷層PGV隨震源距離分布具有一定的相似性,且這種相似性隨與震源距離的增加而愈發明顯。

2 PGV預測方法

對于爆破作業引起的震動事件,薩道夫斯基的經驗公式[13],又稱“Sadov公式”,是目前最為公認、應用面最廣的經驗公式,其具體表達式如下:

vs=K(Q1/3R)α(2)

式中:vs為測點最大震動速度(cm/s);K為與地形、爆破方式等因素相關的系數,堅硬巖石50~150,中等硬度巖石150~250,軟巖250~350;Q為裝藥量(kg);R為測點距爆心的距離(m);α為爆破地震波衰減指數。

董隴軍等[14]根據Sadov公式的原理,發現裝藥量對爆破威力的影響與地震動峰值速度反映出的爆破荷載的能量幅值具有相同趨勢,基于上述理論,用震源PGV替代裝藥量Q作為計算參數,對Sadov公式進行了變形,從而得到新的圍巖爆破震動峰值速度計算公式:

vb=m1v1/30/lm2(3)

式中:vb為質點震動峰值速度(m/s);v0為震源震動峰值速度(m/s);l為測點距爆心的距離(m);m1為受地形、爆破方式影響的相關系數;m2為震動波衰減系數。

對于非爆破事件引起的礦山震動事件,目前已有的PGV計算通用方程式大多是基于McGarr等人先前的工作:

lg Rvmax=AM+1(4)

式(4)變形可得:

vmax=C10AM/R(5)

式中:vmax為質點峰值速度(m/s);M為地震時間的震級;A和C為礦山特征參數。

1996年Kaiser對來自Brunswick礦山、El Teniente礦山和Creighton礦山的地震數據總體進行了95 %置信區間下的回歸分析,綜合McGarr 1984年的研究結果,確定參數A=0.5和C=0.25。因此,有:

vmax=C10a(Mw+1.5)R≈1.4×10(Mw+1.5)R(6)

式中:a為與礦山相關的特征參數。

由上述內容可知,目前國內外對PGV計算方法的研究已經較為成熟。但是,對經驗公式在不同震動事件下適用性的討論研究尚少,工程和研究人員在選擇計算公式時缺少理論依據,所以需要結合現場數據對上述公式進行計算討論,并將理論計算結果與實際監測值進行擬合,根據擬合的相關系數來確定不同公式的適用范圍和精確度[15-16]。

為了獲得礦山特征參數,取Mw=0.053 067 5事件中的震源距離和實測PGV作為數據來源(見表2),對公式(5)進行迭代擬合運算,從而獲得適用于金青頂礦區的特征參數:C=0.260 2和A=0.355 1。

此外,由于PGV和矩震級由震源能量釋放決定,同時能量與速度的平方呈正相關,所以可將公式(3)

中的震源v0替換為矩震級Mw,得到公式(3)的變形公式:

vb=m1M1/6w/Rym2(7)

式中:Ry為測點距震源的距離(m)。

同樣采用Mw=0.334 718 2事件中的數據對公式(7)進行迭代擬合,求得m1、m2的取值為18和1.9。 公式(7)、公式(5)和公式(6)的實際監測值、公式計算值與震源距離的關系曲線見圖5(為表述方便,后文統稱公式(?。?、公式(ⅱ)、公式(ⅲ))。

對公式(?。?、公式(ⅱ)、公式(ⅲ)的計算值與其對應的實測值進行相關系數擬合計算,獲得3個公式的相關系數分別為r1=0.859、r2=0.798和r3=0.756。 結合實測值、計算值同震源距離的關系圖及二者間的擬合相關系數,對3種PGV計算公式的擬合情況和適用性進行分析:

1)在3種公式的關系圖中,PGV隨震源距離的變化趨勢皆與擬合曲線相同,說明3種公式中都體現了PGV隨震源距離的衰減規律。

2)通過比較實測值的擬合曲線與3條計算曲線可發現,3條曲線在與擬合曲線的吻合度上存在差異,公式(?。┯嬎闱€在d>500 m的區間內與擬合曲線有較高的擬合度,公式(ⅱ)和公式(ⅲ)整體上與擬合曲線有較大差異,說明3種PGV計算公式在當前情況下適用性不一。

3)比較三者的相關系數:r1=0.859>r2=0.798>r3=0.756,發現公式(ⅰ)的計算曲線和擬合曲線的吻合程度較高[17],證明公式(?。┑暮瘮登€和PGV的實際分布規律高度相關,而公式(ⅱ)、公式(ⅲ)的擬合度較低,說明對于金青頂礦區的震動事件Mw=0.053 067 5而言,公式(?。┑倪m用性優于另2種公式。

4)3種計算公式對該事件PGV計算所得結果仍存在較大誤差。

針對3種公式在該事件中體現出的適用性差異,結合公式原理和震動事件本身的特征對誤差進行分析:

(1)基于對原始數據的分析,了解到Mw=0.053 067 5 事件的震源以爆破作用為主導,而公式(?。┦怯蒘adov 公式推導變形產生,Sadov公式是以爆破信號為震源的PGV計算公式,其計算原理和有關參數是基于爆破震動的形成機理和傳播規律,所以理論上其變形公式更適用于圍巖爆破震動峰值速度的計算[18];而公式(ⅱ)、公式(ⅲ)是基于McGarr經驗公式的改進,其本身的參數A、C僅與礦山特征相關,而與爆破震動特征相關性較低,所以在計算非剪切作用為主導的震動事件時體現出較低的適用性。

(2)公式(ⅲ)中的礦山特征參數A、C源自于Kaiser對Brunswick礦山、El Teniente礦山和Creighton礦山的地震數據進行95 %置信區間下回歸分析所得的結果,但理論上金屬礦山間一般都存在較大的差異性,因此公式(ⅲ)中的特征參數缺乏普適性。

(3)因為深部金屬礦在地質條件、采場結構和作業情況等方面都極為復雜,所以通過迭代擬合計算和總結出的爆破相關系數m1、衰減系數m2及礦山特征參數C、A可能與礦山實際情況不符,從而導致公式(ⅰ)及公式(ⅱ)的計算結果與實測值間存在一定差異。

3 基于監測數據與預測結果的支護建議

根據金青頂礦區的震動監測數據和分析結果,可以發現不同的震源參數會對采場及巷道圍巖內的PGV分布產生明顯的影響,且在同一震源條件下,PGV的分布在空間上表現出較強的不均勻性。因此,若在支護工作中不考慮震源參數對PGV分布規律的影響,則易導致出現支護手段與圍巖穩定狀態不匹配的現象[19-20]。根據Ju Ma等[12]提出的在考慮震源機制下的圍巖支護設計建議,應基于該地區震動數據監測記錄,對地震事件進行震源機制的反演分析,獲得斷層面解和矩震級的分布規律,根據所獲取的斷層面解和矩震級分布,合成該地區在可能的震源模型下及可能的最大矩震級風險下的PGV分布圖,結合基于動態荷載和巖體質量等級的支護建議(見表3),采取適宜的支護手段。

由于公式(?。┰谠摰V區條件下的適用性優于另外2種,所以選用Sadov公式的變形公式對矩震級Mw=1.65和Mw=0.99的2次震動事件中監測所得數據進行2種震源假設位置下的PGV分布預測分析。

3.1 震源位于采場M8-3~M8-5區域

震源位于采場M8-3~M8-5區域內時所獲得的E8、E9、E10和E11穿脈內圍巖PGV預測結果分布情況見圖6。

1)由圖6-a)中的預測結果可知:在Mw=1.65震源條件下,E8、E9、E10和E11穿脈內圍巖PGV預測最大值分別為0.004 46 m/s、0.004 22 m/s、0.002 89 m/s和0.002 47 m/s,根據表3對各穿脈內的支護作業給出如下建議:

(1)E8、E9穿脈內預測結果較大,建議對2條穿脈靠近采場位置的圍巖進行噴射混凝土+普通錨桿支護,對存在節理裂隙的區域可利用錨網進行補充。

(2)E10、E11穿脈內圍巖PGV預測值偏低,故整體上采用噴射混凝土作為支護手段。但是,E10、E11穿脈與采場的交叉區域內巖石暴露面積較大,應對該區域的圍巖進行噴射混凝土+吸能螺栓+錨網的預支護。

2)由圖6-b)中的預測結果可知:在Mw=0.99震源條件下,E8、E9、E10和E11穿脈內圍巖PGV預測最大值分別為0.001 48 m/s、0.000 96 m/s、0.000 44 m/s 和0.000 19 m/s。通過與圖6-a)所示結果比較可知:

(1)Mw=0.99震源條件下PGV隨震源距離的變化趨勢與Mw=1.65大致相同,但臨近穿脈或采場間的PGV分布差異性較小,整體上PGV分布會更加均勻。

(2)根據上述現象,可將噴射混凝土作為各穿脈的通用支護手段,同時針對作業點附近的圍巖或存在局部破碎的圍巖,可考慮在噴射混凝土的基礎上安裝普通錨桿作為補充支護。

3.2 震源位于E10-1鉆孔硐室內

震源位于E10-1鉆孔硐室內時所獲得的E8、E9、E10和E11穿脈內圍巖PGV預測結果分布情況見圖7。

1)由圖7-a)中的預測結果可知:在該震源條件下,E8、E9、E10和E11穿脈內圍巖PGV預測最大值分別為0.001 84 m/s、0.002 71 m/s、0.004 49 m/s和0.003 22 m/s,根據表3對各穿脈內的支護作業給出如下建議:

(1)E10穿脈內預測結果最大,且由于穿脈與采場及鉆孔硐室相接,圍巖的暴露面積大。故建議對E10穿脈靠近采場位置的圍巖進行噴射混凝土+普通錨桿+錨網支護,對靠近主巷道的圍巖進行噴射混凝土支護。

(2)E11穿脈內圍巖PGV預測值偏低,建議對E11穿脈采用噴射混凝土的方式作為支護手段。由圖7可知,E11、E10穿脈間采場內的預測值偏大,且該采場內頂板圍巖破碎并已設置木樁支護,所以建議在木樁支護的基礎上添加噴射混凝土+普通錨桿作為補充支護。

(3)E8、E9穿脈內PGV預測值普遍低于0.002 5 m/s,故對E8、E9穿脈圍巖采取噴射混凝土支護即可。

2)由圖7-b)中的預測結果可知:在Mw=0.99震源條件下,E8、E9、E10和E11穿脈內圍巖PGV預測最大值分別為0.000 12 m/s、0.000 67 m/s、0.001 50 m/s 和0.001 48 m/s,根據表3對各穿脈內的支護作業給出如下建議:

(1)圖7-b)中各穿脈內PGV值都較小,說明該情況下震動對圍巖穩定性的影響較弱,且各穿脈間PGV分布的差異性較低,所以建議對E8、E9、E10和E11穿脈內的圍巖采用噴射混凝土支護。

(2)由于E10穿脈與采場及鉆孔硐室相接,其圍巖的暴露面積大。故建議對E10穿脈靠近采場位置的圍巖進行噴射混凝土+普通錨桿支護。

4 結 論

本文在研究國內外震動機理和相關理論的基礎上,采用理論分析、現場監測數據分析和擬合相關系數計算等手段對不同震源參數下PGV分布規律和預測技術進行了深入研究和分析,得出以下結論:

1)震源距離對PGV分布產生明顯作用;且PGV隨震源距離的衰減程度與震源距離之間存在密切的聯系,當距震源較遠時,震源距離對PGV的衰減作用逐漸弱化。

2)相同震源距離下的PGV會隨矩震級增大顯著變大,且這種影響會根據距震源遠近逐漸衰減。

3)不同的斷層滑移角下的PGV大小和分布情況存在較大差異,一般情況下走滑斷層(αr=0°)>正斷層(αr=45°)>垂直斷層(αr=±90°)≈逆斷層(αr=-45°)。

4)Sadov公式的變形公式在爆破條件下PGV計算中體現出良好的適用性,同時也說明運用與震動事件的類型和機理相符的公式,所得結果更能反映PGV的實際分布規律。

5)通過迭代擬合獲得的McGarr經驗公式中的特征值A、C不能完全適合該礦區,且Kaiser回歸分析獲得的特征值在特定礦山下缺乏適用性,導致計算結果與實測值間的相關性較低。

6)在制定采場或巷道的支護方式時,應結合不同震動事件下,支護區域及周邊圍巖的PGV分布情況,根據不同震源機制下的圍巖支護設計建議選擇更加合理和時效化的支護方式。

[參 考 文 獻]

[1] MORRISION R G K.Theory and the practical problem of rock bursts[J].Engineering and Mining Journal,1948,149(3):66-72.

[2] 馮叔瑜,呂毅,顧毅成.城市控制爆破[M].北京:中國鐵道出版社,1987.

[3] 何姣云.礦山采動災害監測及控制技術研究[D].武漢:武漢理工大學,2007.

[4] 熊仁欽.頂板大面積來壓破壞機理的研究[J].煤炭學報,1995(增刊1):38-41.

[5] 陳壽如,宋光明,史秀志,等.近河堤采礦爆破的震動監測與控制[J].中國有色金屬學報,2000,10(1):136-139.

[6] 韓子榮.金川礦區露天地下聯合開采的爆破震動安全性評價[J].礦冶工程,1985,5(1):6-11.

[7] 王民壽,郭慶海.用雙隨機變量回歸改進爆破震速回歸分析[J].爆炸與沖擊,1998,18(3):283-288.

[8] 徐全軍,劉強,聶渝軍,等.爆破地震峰值預報神經網絡研究[J].爆炸與沖擊,1999,19(2):133-138.

[9] 徐全軍,張慶明,揮壽榕.爆破地震峰值的神經網絡預報模型[J].北京理工大學學報,1998,18(14):472-475.

[10] 黃光球,桂中岳.確定爆破工程中真實經驗公式的遺傳規劃方法[J].工程爆破,1997,3(3):15-22.

[11] 熊代余,顧毅成.巖石爆破理論與技術新進展[M].北京:冶金工業出版社,2002.

[12] MA J,DONG L J,ZHAO G Y,et al.Qualitative method and case study for ground vibration of tunnels induced by faultslip in underground mine[J].Rock Mechanics and Rock Engineering,2018,52(3):1-15.

[13] 田蜜.爆破擾動誘發煤礦沖擊地壓的數值模擬研究[D].阜新:遼寧工程技術大學,2014.

[14] 董隴軍,王鈞暉,馬舉.不同微震震源機制下地下硐室圍巖響應及支護建議[J].隧道與地下工程災害防治,2019,1(3):68-76.

[15] 肖和平.煤礦構造礦震機理[J].湖南地質,1999,18(2/3):141-146.

[16] 肖和平.煤礦礦震應力窗口效應[J].華南地震,1999,19(1):85-90.

[17] 董隴軍,李夕兵,唐禮忠,等.無需預先測速的微震震源定位的數學形式及震源參數確定[J].巖石力學與工程學報,2011,30(10):2 057-2 067.

[18] 鄒景波.爆破地震波作用下結構的動力響應及安全評價研究[D].青島:青島理工大學,2012.

[19] 毛暉.建筑物爆破震動的安全控制技術研究[D].長沙:中南大學,2004.

[20] 許紅濤,盧文波.幾種爆破震動安全判據[J].爆破,2002,19(1):8-10.

Study on the distribution and prediction of PGV in surrounding

rock under different seismic source parameters

Wang Xi1,Li Shiqing2,Li Qiutao2,Cheng Li1,Liu Huanxin1

(1.Deep Mining Laboratory Subsidiary of Shandong Gold Mining Technology Co.,Ltd.;

2.Shandong Jinzhou Mining Group Co.,Ltd.)

Abstract:In order to totally understand the law of seismic wave propagation and distribution in the surrounding rock under different seismic source mechanisms,the paper studies the distribution of PGV under different seismic source parameters based on field monitoring data.The results show that the distance between stations and the source,moment magnitude scale,and the rake of fault slip have different effects on the distribution of PGV.The applicability of each prediction methods is analyzed by comparing the fitting correlation coefficient between PGV calculated from empirical formulas and others from field monitoring.The study finds that the predicted value of PGV calculated from the Sadovs prediction formula has a higher correlation to the measured result with blasting vibration event.Based on the above research results,the spatial distribution of PGV under dynamic seismic mechanisms is calculated before suggestions on the surrounding rock support under different conditions of the seismic source location and magnitude are given.

Keywords:metal mine;seismic source;vibration monitoring;PGV;prediction formula;fitting correlation coefficient

主站蜘蛛池模板: 午夜国产精品视频黄 | 毛片网站观看| 亚洲欧洲自拍拍偷午夜色无码| 成人欧美日韩| jizz国产视频| 亚洲综合经典在线一区二区| 欧美日韩资源| 日本黄网在线观看| 久久精品中文字幕免费| 免费a级毛片18以上观看精品| 欧美午夜在线观看| 国产精品网址你懂的| 亚洲日本中文字幕乱码中文 | 理论片一区| 欧美一级专区免费大片| 亚洲香蕉在线| 国产又粗又爽视频| 亚洲欧美日本国产综合在线 | 国产激爽大片高清在线观看| 国产女人在线| 中文字幕无码电影| 久久男人视频| 女人18毛片久久| 亚洲国产天堂久久综合226114| 亚洲精品波多野结衣| 无码精品国产dvd在线观看9久| 在线另类稀缺国产呦| 国产精品私拍在线爆乳| 18禁黄无遮挡网站| 爱爱影院18禁免费| 嫩草在线视频| 国产黄色免费看| 国产综合色在线视频播放线视| 亚洲天堂免费观看| 国产成人福利在线| 久久午夜夜伦鲁鲁片无码免费| 免费一级无码在线网站| 国产一区免费在线观看| 国内精品视频| 午夜成人在线视频| 色网站在线视频| 国产精品理论片| 国内精品91| 永久免费无码日韩视频| 国产流白浆视频| 欧美三级不卡在线观看视频| 一级毛片免费的| 综合天天色| 99爱在线| 一本视频精品中文字幕| 5388国产亚洲欧美在线观看| 日韩av无码精品专区| 欧美一道本| 国产免费网址| a级毛片网| 亚洲日本精品一区二区| AV天堂资源福利在线观看| 国产精品人莉莉成在线播放| 找国产毛片看| 九九热在线视频| 在线精品视频成人网| 亚洲国产精品日韩专区AV| 熟女视频91| 中文字幕天无码久久精品视频免费| 四虎永久在线精品国产免费| 国产精品9| 青青青国产视频| 免费高清a毛片| 91探花在线观看国产最新| 中文字幕久久波多野结衣| 日韩精品毛片人妻AV不卡| 国产麻豆精品手机在线观看| 久久精品91麻豆| 波多野结衣中文字幕一区二区| 亚洲精品色AV无码看| 亚欧乱色视频网站大全| 五月婷婷丁香综合| 日韩A∨精品日韩精品无码| 国产午夜一级毛片| 精品无码国产一区二区三区AV| 国产亚洲精品97AA片在线播放| 午夜视频www|