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

黃土高填方場地工后沉降預測新模型

2022-12-16 08:37:14于永堂鄭建國
西南交通大學學報 2022年6期
關鍵詞:工程模型

于永堂,鄭建國

(1.機械工業勘察設計研究院有限公司,陜西 西安 710043;2.西安建筑科技大學土木工程學院,陜西 西安 710055;3.中聯西北工程設計研究院有限公司,陜西 西安 710077)

我國西部黃土丘陵溝壑區地形起伏、溝壑縱橫、平地較少,為了增加建設用地,近年來實施了大量挖填交替、方量巨大的黃土高填方工程,一些工程的填方厚度達幾十米甚至上百米.黃土高填方場地的沉降變形除在施工期發生外,工后期還會繼續發生,且往往會持續很長時間.實際上,后續工程往往無法等到場地的工后沉降完全穩定時才開始修建,若建(構)筑物修建時間過早,剩余沉降或差異沉降超過其所能承受的限值時,將會導致其損毀破壞,嚴重時甚至無法使用.因此,準確地預測工后沉降,對確定后續工程的建設時機、保障工程安全和指導建(構)筑物的規劃布局等具有重要意義.

黃土高填方場地涉及地形地質條件復雜的原地基以及填料性質特殊、填土厚度大幅度變化的填筑體,受地下水位變動和地表水入滲等環境條件變化影響,準確預測工后沉降難度較大.沉降預測方法主要有兩類[1-3]:第一類是基于固結理論、本構模型的理論方法;第二類是基于前期實測數據外推預測的經驗方法.由于理論方法常存在大量簡化假定、模型參數不易準確測定等原因,導致理論預測值與實測值往往存在較大差異,實際工程中尚難以完全依靠理論方法進行沉降預測,而常采用基于前期實測數據外推預測的經驗方法,如雙曲線模型法[4]、指數函數模型法[5]、冪函數模型法[6]、星野法[7]、Weibull模型法[8]、Logistic模型法[9]、Gompertz模型法[10]、鄧英爾模型法[11]、灰色理論法[12]、神經網絡法[13]、數值反演分析法[14]等,其中回歸參數模型法因計算簡便、實用性強等特點,被工程技術人員廣泛采用.上述預測方法均有其適用范圍及適用條件,筆者在工程中實際應用時發現,一些預測模型存在收斂過早或發散嚴重等問題,并不適用于黃土高填方場地的工后沉降預測.

本文根據典型黃土高填方場地的工后沉降數據特點、曲線特征和發展演化規律,提出了收斂型和發散型兩種用于工后沉降預測的新回歸參數模型(簡稱新模型),介紹了新模型的基本性質與參數求解方法,并結合實測沉降數據檢驗了模型的預測精度和適應性.新模型可為黃土高填方場地的工后沉降預測提供參考.

1 典型黃土高填方工程概況

陜北某市地處黃土丘陵溝壑區,城市建在溝谷之中,因地形條件限制,建設用地緊張,近年來為增加建設用地,開展了多處挖填造地工程,其中工程Ⅰ、Ⅱ分別是該市在原中心主城區的北部和南部實施的兩處典型黃土高填方工程.現將兩處工程的基本情況簡要介紹如下:

工程Ⅰ的建設面積約10.5 km2,長度約5.5 km,寬度約2.0 km,最大填方厚度約112 m.本工程原始地形為溝谷形,通過開挖梁峁區土體填入溝谷中,形成大厚度填方場地.場地內出露的地層主要為第四系全新統洪積層粉土、上更新統馬蘭黃土、中更新統離石黃土、新近系棕紅色、暗紫色紅黏土和侏羅系砂巖、泥巖互層.工程區域分布有較大范圍的人工筑壩攔淤形成的淤積土,具有結構松散、含水率高,高填方荷載作用下沉降量大、變形穩定時間長等特點.地下水類型有第四系孔隙潛水和侏羅系基巖裂隙水兩大類,自周邊分水嶺地帶順地勢向溝谷徑流匯集,轉化為地表徑流排泄于區外.本工程建設過程在溝底設置了地下盲溝排水系統用于疏排地下水,對原地基采取強夯法處理,填筑體采用分層碾壓法(沖擊碾壓、振動碾壓)處理,填土的壓實系數要求不小于0.93 (重型擊實試驗控制).本工程自2012年4月17日起開工建設,至2013年10月底土方施工陸續完工,隨后布設地表沉降觀測點,觀測工作自2013年11月10日起,至2017年8月7日,累計觀測時間歷時1366 d.

工程Ⅱ的建設面積約0.093 km2,長度約360 m,寬度約260 m,最大填方厚度約42 m.本工程填方區原始地形為溝谷形,梁峁區主要地層為第四系上更新統馬蘭黃土及中更新統離石黃土、新近系紅黏土和侏羅系砂泥巖;沖溝區分布有第四系全新統沖積層.場地內的地下水主要為基巖裂隙水,以下降泉的形式在溝谷內出露但水量較少.高填方施工過程在填方區沿原溝谷溝底設置了地下排水盲溝,溝谷區原地基采用強夯法處理,填筑體采用分層碾壓法處理,填土壓實系數要求不小于0.95 (重型擊實試驗控制).本工程自2009年8月28日起開工建設,至2010年9月底土方施工陸續完工,隨后布設地表沉降觀測點,觀測工作自2010年10月6日起,至2012年9月17日,累計觀測歷時713 d,共獲得32期沉降數據[15].

2 黃土高填方場地的工后沉降特征

根據前述兩處典型黃土高填方工程的工后沉降數據,對黃土高填方場地的工后沉降特征進行分析.工程Ⅰ、Ⅱ的工后沉降數據均為非等時距,前期觀測周期短、后期觀測周期長.

圖1代表工程Ⅰ、Ⅱ中典型監測點的工后沉降變形.圖1(a)中監測點J1、J2處的填土厚度分別為103.8、64.2 m,圖1(b)中監測點S1、S2處的填土厚度分別為32.4、27.3 m,剔除了觀測歷時為469 d和510 d (當年11月至次年3月間冰凍期)受凍融影響的異常沉降數據.

圖1 典型黃土高填方場地的工后沉降曲線Fig.1 Duration curves of post-construction settlement of typical loess high fill site

由圖1可知:黃土高填方場地的工后沉降曲線主要呈“J”形和“S”形兩種形態,由于各監測點的填土厚度相差較大,使壓縮土層厚度及自重荷載相差大,不同監測點的沉降量也差異較大;工后期階段,填土施工剛完成,自重荷載不再增加,土方施工加載引起的瞬時沉降已經完成,絕大多數監測點的工后沉降曲線如圖1(a)所示,曲線形態呈現出由快速增長向平緩增長,沉降曲線總體呈緩變型,近似呈“J”形曲線形態,根據填方工程沉降曲線的發展特點,最終將逐漸趨于某一定值;部分工后沉降曲線如圖1(b)所示,表現出初始階段短時相對平緩地增長(初始沉降階段),然后進入快速增長階段(加速沉降階段),達到一定程度后又趨于平緩(趨穩沉降階段),最終也將趨于某一定值(極限沉降階段),近似呈“S”形曲線形態.

工后沉降的發展是事物發生、發展、成熟并到達一定極限的過程,因此可以借此對“S”形曲線的產生機理進行解釋.該類曲線的發展過程可分為以下4個階段:

1)第Ⅰ階段:初始沉降階段.在沉降變形初期,上部非飽和土層排氣條件好,下部非飽和土層的排氣條件差.已有研究表明,排氣條件是影響非飽和填筑體沉降發展的關鍵因素,排氣條件較好土層的工后初期沉降量甚至是排氣條件較差土層沉降量的2倍以上[16].此時,測點主要反映上部土層的彈性或近似彈性的沉降變形.

2)第Ⅱ階段:加速沉降階段.上部土層的荷載逐步傳遞到下部土層,引起下部土體進入彈塑性狀態.隨著塑性區的不斷開展,測點的沉降速率也在不斷地增加.

3)第Ⅲ階段:趨穩沉降階段.由于填土固結尚未完成以及土體的流變,測點的沉降將隨著時間的推移而繼續,但沉降速率遞減.

4)第Ⅳ階段:極限沉降階段.從極限理論上講,當時間為無窮大時,沉降最終將達極限狀態,此時沉降將不隨時間發生變化.

由于填土分層填筑荷載相對較小,加之施工結束后首先需埋設沉降觀測標點,無法立即進行沉降觀測,即存在一定的時間滯后性,導致代表“S”形沉降曲線“發生”的“初始沉降階段”在觀測到的全過程沉降曲線所占時間很短,沒有文獻[17]中的沉降曲線那么明顯.當未觀測得到該階段變形時,所沉降-時間曲線則呈現為“J”形.

3 模型的函數表達式與基本性質

本次針對黃土高填方場地的工后沉降數據特點、曲線特征和發展演化規律,提出了發散型和收斂型兩種用于預測工后沉降的新模型.

3.1 新模型Ⅰ的函數表達式與基本性質

新模型Ⅰ的函數表達式為

式中:st為預測沉降量;t為時間;a、b、c、d為待求模型參數,均大于0.

該新模型具有以下基本性質:

1)過原點:當t= 0時,s0= 0,模型不包含初始加載時的瞬時沉降.

2)無界性:當t→∞ 時,st→+∞,模型無法直接獲得最終沉降量,僅能以達到某一較小沉降速率時的沉降量作為最終沉降量.

3)單調性:對st求一階導數,可得沉降速率s′t的函數表達式(2),由s′t>0可知,模型單調遞增.

4)適應性:對式(1)求二階導數,可得

圖2 新模型Ⅰ的參數變化對曲線形態的影響Fig.2 Influences of parameter change of the new model I on the shape of curves

由圖2可知:通過調整模型參數,本模型能夠在較大范圍內描述幾何上為“J”形和“S”形的沉降曲線.當討論其他參數變化時的曲線形態,可采取類似方法進行計算和分析.

選取圖2(a)中“J”形曲線模型參數為a=2.000,b= 0.500,c= 0.100,d= 0.500以及圖2(b)中“S”形曲線模型參數為a= 0.200,b= 0.100×10?5,c= 0.500×10?3,d= 0.400的沉降曲線,根據式(2)、式(3)繪制沉降速率、沉降加速度的全程變化曲線,如圖3.由圖3可知:類型A的沉降速率和沉降加速度均持續降低,類型B的沉降速率存在一個峰值點,沉降加速度存在正負兩個峰值點;對于類型B,若令s′t′=0,則可確定本模型預測曲線拐點(沉降速率最大值點)所對應的時間,同時,若令s′t′′=0,可確定出沉降加速度函數的兩個峰值點所對應的時間.綜合上述分析可知,通過調整模型參數a、b、c、d值,可模擬相當大變化范圍內的沉降曲線,能夠描述幾何形態上為“J”形和“S”形的沉降曲線,表明新模型Ⅰ具有較強的適應性.

圖3 新模型Ⅰ的沉降速率和沉降加速度全程變化曲線Fig.3 Settlement velocity and acceleration curves of the new model I

3.2 新模型Ⅱ的函數表達式與基本性質

鑒于新模型Ⅰ屬發散型模型,無法直接獲得最終沉降量,故又提出一種收斂型新模型Ⅱ,其函數表達式為

該新模型具有以下基本性質:

1)過原點:當t=0時,s0=0,可見該模型通過原點,不包含初始加載時的瞬時沉降.

2)有界性:當t→∞ 時,s1t→a/b,表明該模型屬于收斂模型,a/d為模型的最終沉降量.

3)單調性:對s1t求一階導數如式(5)所示,由沉降速率s′1t>0,表明該模型單調遞增.

新模型Ⅱ中參數變化時的曲線形態如圖4所示.由圖4可知:固定其他參數不變,僅參數b在(0,1]范圍內變化,能夠在較大范圍內描述“J”形沉降曲線特征;固定其他參數不變,僅參數b在(1, + ∞)范圍內變化,能夠在較大范圍內描述“S”形曲線特征.

圖4 新模型Ⅱ的參數變化對曲線形態的影響Fig.4 Influences of parameter change of the new model Ⅱ on the shape of curves

當a= 1.200,b= 1.00,c= 0.600,d= 0.005(類型A':對應“J”形曲線)和a= 0.300,b= 2.200,c= 50.000,d= 0.001 (類型B':對應“S”形曲線)時沉降速率、沉降加速度的全程變化曲線如圖5.由圖5可知:類型A' 的模型參數b在(0, 1]范圍內,沉降速率和沉降加速度均持續降低;類型B' 的模型參數b在(1, + ∞)范圍內,沉降速率存在一個峰值點,沉降加速度存在正負兩個峰值點,峰值點處曲線斜率為0.與新模型Ⅰ類似,對于新模型Ⅱ中的類型B',若令=0及=0 則分別可確定本模型沉降加速度函數的兩個峰值點所對應的時間及模型預測曲線的拐點(沉降速率最大值點)所對應的時間.通過以上分析可以發現,同樣通過調整模型參數a、b、c、d值,可實現對不同類型沉降曲線的模擬.模型所反映出來的沉降量、沉降速率以及沉降加速度隨時間的變化規律均與實際黃土高填方場地工后沉降監測過程中遇到的“S”形或“J”形沉降曲線的特征相符合.

圖5 新模型Ⅱ的沉降速率和沉降加速度全程變化曲線Fig.5 Settlement velocity and acceleration curves of the new model Ⅱ

4 模型參數的求解與估計

由式(1)、(4)可知,新模型Ⅰ、Ⅱ均為4參數非線性方程,很難采用解析法直接求解,故采用數值法求解.設有m組實測沉降數據樣本 (,,···,)用于建模,由最小二乘法原理建立第i期的預測值sti與實測值之間的目標函數Q,如式(7),對目標函數求極小值即為對待定參數a、b、c、d的尋優過程.

式中:eti為第i期預測值與實測值之差,i= 1,2,…,m;ti為第i期數據對應的時間.

采用Levenberg-Marquardt優化算法[18]對模型中的4個參數進行尋優.該方法是用模型函數對待估參數向量在其領域內做線性近似,利用泰勒展開,忽略二階以上的導數項,將優化目標方程轉化為線性最小二乘法問題進一步求解,它具有收斂速度快、適應性較強的特點,易于編程實現.本次利用MATLAB軟件中成熟的非線性復雜模型參數估計求解功能,確定模型參數.當求解模型參數時,由軟件隨機自動給出初始值進行迭代運算,最終找出最優解,確定出模型參數a、b、c、d的值.

5 模型預測效果的評價方法

當采用回歸參數模型進行沉降預測時,常會出現內擬合誤差低而外推預測誤差高,或者內擬合誤差高而外推預測誤差低的情況.因此,需要從內擬合誤差和外推預測誤差兩方面綜合評價模型的預測效果.

5.1 內擬合誤差的評價方法

本次對模型擬合效果的評價指標為決定系數(R2),R2可采用式(8)計算.R2越接近于1,則表明實測數據通過模型的解釋性就越強.

5.2 外推預測誤差的評價方法

對模型預測效果的評價指標為平均絕對百分比誤差(MAPE,MAPE)、平均預測誤差(MFE,MFE),分別采用式(9)、式(10)計算.MAPE能夠較好地衡量模型的預測精度,其值越小精度越高,當MAPE<10%時,是一個比較好的預測模型[19];MFE能較好體現模型的無偏性(如預測值相對于實測值的正負偏差).

6 工程實例分析與效果檢驗

本次為了綜合評價模型的預測效果,將已知沉降數據分為前后兩部分,若實測數據共n期,前一部分m期數據來識別和估計模型參數,后一部分(n?m)期數據來評價模型的外推預測效果.本次從前文所述工程Ⅰ、Ⅱ中共分別選取2個典型監測點的沉降數據,采用本文提出的兩個新模型預測其工后沉降,檢驗模型的實際預測效果.

工程Ⅰ中典型監測點JC6、T13的工后沉降曲線呈“J”形,其中監測點JC6觀測歷時642 d,共有32期實測數據,將實測數據分為前22期和后10期,利用前22期實測數據求解模型參數,然后采用向后預測的10期數據(外推預測時長/總數據時長 =49.5%)與后10期實測數據比較,檢驗模型的預測效果.監測點T13觀測歷時1366 d,共有70期實測數據,將實測數據分為前60期和后10期,利用前60期實測數據求解模型參數,然后采用向后預測的10期數據(外推預測時長/總數據時長 = 23.6%)與后10期實測數據比較,檢驗模型的預測效果.

工程Ⅱ中典型監測點S1、S2的工后沉降曲線呈“S”形,兩個監測點觀測歷時712 d,均有32期實測數據,將實測數據分為前22期和后10期,利用前22期實測數據求解模型參數,然后采用向后預測的10期數據(外推預測時長/總數據時長 = 75.6%)與后10期實測數據比較檢驗模型的預測效果.兩種新模型的回歸模型參數及預測效果評價指標結果如表1所示.

表1 新模型的回歸模型參數及預測效果評價結果Tab.1 Regression parameters and evaluation results of the new models

6.1 新模型對“J”形沉降曲線的預測效果

采用新模型對工程場地Ⅰ中典型監測點JC6、T13進行擬合及預測,沉降曲線如圖6所示.表1中MFE計算值結合圖6可知:新模型Ⅰ的預測值較實測值總體呈偏高(正偏差)趨勢,新模型Ⅱ的預測值較實測值總體呈偏低(負偏差)趨勢;根據檢驗數據計算獲得監測點JC6、T13后10期工后沉降預測結果的MAPE值,其中新模型Ⅰ分別為4.2%和1.2%,新模型Ⅱ分別為1.9%和0.7%,表明兩種新模型的預測精度均較高,此時新模型Ⅰ對“J”形沉降曲線的預測精度低于新模型Ⅱ;在檢驗數據區段內,絕大多數實測值基本處于兩種新模型預測曲線包絡帶范圍內.因此,可通過新模型Ⅰ和新模型Ⅱ來預測未來沉降區間.

圖6 采用新模型進行擬合及預測的沉降曲線(工程Ⅰ)Fig.6 Settlement curves fitted and predicted by the new models ( project Ⅰ)

6.2 新模型對“S”形沉降曲線的預測效果

新模型對監測點S1、S2的擬合及預測曲線如圖7所示.表1中MFE結合圖7可知:新模型Ⅰ的預測值較實測值總體呈正偏差(偏高),新模型Ⅱ的預測值較實測值總體呈負偏差(偏低).由表1可知,監測點S1、S2后10期工后沉降預測值的MAPE值,當采用新模型Ⅰ時分別為3.5%和4.6%,采用新模型Ⅱ時分別為11.2%和6.2%,表明兩種新模型對該類沉降曲線的預測精度均較高,此時新模型Ⅰ對“S”形沉降曲線的預測精度總體高于新模型Ⅱ.

圖7 采用新模型進行擬合及預測的沉降曲線(工程Ⅱ)Fig.7 Settlement curves fitted and predicted by the new models ( project Ⅱ)

6.3 新模型與傳統回歸參數模型的預測效果比較

在傳統回歸參數模型中,雙曲線模型、指數函數模型、對數函數模型、冪函數模型、平方根函數模型和星野法等較適合預測“J”形曲線;Weibull模型、Morgan-Mercer-Flodin (MMF)模型、改進Knothe模型、鄧英爾模型、Logistic模型、Gompertz模型、Usher模型、Spillman模型及Janoschek模型等較適合預測“S”形曲線.采用本文提出的兩種新模型及傳統回歸參數模型分別對具有“J”形(監測點JC6)及“S”形(監測點S2)曲線特征的沉降數據進行建模和外推預測后,得到各模型的擬合與預測曲線如圖8所示,文中所列的傳統回歸參數模型出現了收斂過早或發散嚴重的問題.在圖8(a)中:新模型Ⅰ的R2= 1.000,MAPE = 4.2%;新模型Ⅱ的R2= 0.999,MAPE = 1.9%;傳統回歸參數模型R2的變化范圍為0.951 ~ 0.999,MAPE的變化范圍為6.0% ~ 12.2%.新模型Ⅰ、新模型Ⅱ相較于傳統預測模型的外推預測誤差分別降低了30.0% ~ 65.6%和68.3% ~ 84.4%.在圖8(b)中:新模型Ⅰ的R2= 0.998,MAPE = 4.6%;新模型Ⅱ的R2= 0.998,MAPE = 6.2%;傳統回歸參數模型R2的變化范圍為0.990 ~ 0.998,MAPE的變化范圍為7.5% ~ 38.0%.新模型Ⅰ、新模型Ⅱ相較于傳統預測模型的外推預測誤差分別降低了78.7% ~95.8%和17.3% ~ 83.7%.由以上分析可知,兩種新模型與傳統回歸參數模型對“J”形和“S”形沉降曲線的擬合精度指標R2相差不大,但兩種新模型的預測精度指標MAPE均明顯小于傳統回歸參數模型,表明兩種新模型的工后沉降預測精度明顯優于文中所列的傳統回歸參數模型,且具有較強的適用性和穩定性.

圖8 新模型與傳統模型的工后沉降擬合與預測曲線Fig.8 Fitting and prediction curves of post-construction settlement between the new models and some other conventional models

7 結 論

本文在分析典型黃土高填方場地工后沉降數據特點、曲線特征和發展演化規律后,提出了發散型和收斂型兩種回歸參數預測新模型,并對模型的預測效果進行了工程實例分析和檢驗,主要結論如下:

1)黃土高填方場地工后實測沉降曲線的幾何形態主要為“J”形和“S”形兩類,以“J”形曲線為主,代表“S”形沉降曲線“發生階段”的“初始沉降階段”在全過程沉降曲線所占時間很短.

2)對模型預測效果評價時,建議將已知實測數據分為前后兩部分,分別檢驗擬合精度和外推預測精度.擬合精度可采用決定系數評價,外推預測精度可采用平均絕對百分比誤差和平均預測誤差綜合評價.

3)兩種新模型均能較好地反映“J”形和“S”形兩類工后沉降曲線的變化規律,表明新模型具有較好的適應性和通用性,其中發散型模型對“S”形沉降曲線的預測效果較好,收斂型模型對“J”形沉降曲線的預測效果較好.

4)兩種新模型對實測數據均具有較高的擬合精度,其中發散型新模型的預測值較實測值總體呈偏高趨勢,收斂型新模型的預測值較實測值總體呈偏低趨勢,用于模型檢驗的實測值多分布于由兩新模型預測曲線構成的包絡帶范圍內.

5)兩種新模型未考慮黃土濕陷變形的影響,主要反映壓縮變形和蠕變變形的長期發展趨勢,在一定程度上克服了傳統回歸參數模型收斂過早或發散嚴重的問題,在黃土高填方場地工后沉降預測的適用性、穩定性和準確性等方面,明顯優于文中所列的傳統回歸參數模型,具有較好的工程實用價值.

猜你喜歡
工程模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
子午工程
太空探索(2016年6期)2016-07-10 12:09:06
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
工程
工程
工程
工程
主站蜘蛛池模板: 国产精品视频导航| 亚洲Av激情网五月天| 日韩黄色在线| 狠狠色噜噜狠狠狠狠奇米777| 99视频在线免费观看| 黄网站欧美内射| 视频一区视频二区日韩专区| 精品無碼一區在線觀看 | 91最新精品视频发布页| 久久综合丝袜长腿丝袜| 第一页亚洲| 毛片卡一卡二| 日本久久久久久免费网络| 国产一在线| 天天视频在线91频| 中文无码精品A∨在线观看不卡 | 亚洲人成网18禁| 伊人色天堂| 无码人中文字幕| 一本久道久久综合多人| 亚洲男人天堂2020| 国产精品视频白浆免费视频| 久久精品一卡日本电影| 亚洲一区二区日韩欧美gif| 国产精品女主播| 国产一级在线播放| 狠狠色综合网| 日韩精品少妇无码受不了| 亚洲色无码专线精品观看| 国产精品 欧美激情 在线播放 | 亚洲最新地址| 亚洲人成高清| 欧美成人一区午夜福利在线| 亚洲人成日本在线观看| 精品福利网| 综合人妻久久一区二区精品| 六月婷婷综合| 欧美日韩精品一区二区在线线 | 999福利激情视频| 性网站在线观看| 91成人精品视频| 国产伦片中文免费观看| 毛片免费在线| 偷拍久久网| 国产真实乱了在线播放| 国产呦视频免费视频在线观看| 5555国产在线观看| 国产精品偷伦视频免费观看国产| 六月婷婷精品视频在线观看| 香蕉网久久| 久久国产乱子| 欧美日韩另类国产| 久久精品66| 国产97视频在线| 幺女国产一级毛片| 四虎综合网| 91亚洲精选| 欧美在线观看不卡| 国产毛片不卡| 啦啦啦网站在线观看a毛片| 999国产精品| 毛片免费视频| 一级毛片免费不卡在线| 99国产精品免费观看视频| 日韩国产精品无码一区二区三区| 九色在线观看视频| 欧美色99| 午夜啪啪网| 国产一级视频久久| 波多野结衣中文字幕一区二区| 亚洲天堂精品在线观看| 成人在线欧美| 日韩欧美在线观看| 免费不卡在线观看av| 青青热久免费精品视频6| 亚洲国产欧美国产综合久久 | 97国产精品视频自在拍| 欧美日韩一区二区在线免费观看| av一区二区无码在线| 亚洲精品动漫在线观看| 国产精品欧美在线观看| 999国产精品|