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

基于聚類與Transformer算法的上軟下硬地層盾構滾刀正常磨損預測*

2024-02-26 07:36:44常佳奇黃宏偉張東明李章林
施工技術(中英文) 2024年1期
關鍵詞:模型施工

常佳奇,黃宏偉,張東明,李章林

(1.同濟大學地下建筑與工程系,上海 200092; 2.同濟大學巖土及地下工程教育部重點實驗室,上海 200092;3.上海隧道工程股份有限公司,上海 200032)

0 引言

盾構法作為隧道施工的主要方法,因其安全性高、對周邊環境擾動小、施工速度快、地層適應性強,得到了廣泛應用。隨著中國隧道建設的發展,盾構隧道正朝著大直徑、大埋深、長距離的方向發展。隨之而來的則是富水砂卵石地層、上軟下硬地層、極硬巖地層等復雜地層條件給盾構施工帶來的巨大挑戰。刀具作為盾構機切削、破碎地層的主要工具,在不良地質條件下極易受到嚴重磨損,如果無法準確預測刀具磨損,往往會導致施工參數異常、刀具破巖率下降、刀盤損壞、盾構掘進困難等問題,甚至會影響盾構掘進安全,造成嚴重的生命財產損失。因此,復雜地層條件下盾構刀具的磨損預測是十分必要的。

盾構刀具根據破碎、切削的地層不同可分為破巖的滾刀和切削軟土的切削刀,其中滾刀往往磨損較為嚴重,是研究的重點。在磨損機理研究上,科羅拉多礦業大學的CSM模型[1]和挪威科技大學的NTNU模型[2]建立了滾刀受力模型,并提出磨蝕性系數和磨蝕值預測滾刀壽命;Xue等[3]建立了滾刀破巖的離散元模型,分析滾刀的破巖過程及滾刀磨損情況;李剛等[4]則通過盾構比能的變化,從能量角度對盾構滾刀的磨損進行了分析預測;趙海鳴等[5]則基于磨粒磨損機理提出了滾刀磨損量的預測模型;吳俊等[6]將滾刀的磨損分為磨粒磨損、黏著磨損和疲勞磨損,其中磨粒磨損和黏著磨損是刀具磨損的主要原因。在滾刀磨損的預測上,袁大軍等[7]通過對工程實際刀具磨損情況的現場測量結果進行統計分析,總結出了砂卵石地層盾構刀具磨損與切削距離間的關系;陳子義[8]借助神經網絡方法建立了施工參數與滾刀磨損量之間的關系;楊俊哲等[9]建立了兩種支持向量機模型,基于盾構施工參數對滾刀磨損進行預測。

然而上述研究大多針對全斷面硬巖地層,全斷面硬巖地層滾刀破巖過程相對穩定,而復合地層中掌子面上軟下硬地層彼此交互,滾刀受力機理更為復雜,因此針對上軟下硬復合地層中滾刀磨損的研究較少,且多為經驗總結。李強等[10]依據杭州地鐵2號線工程,對盾構穿越上軟下硬復合地層段時的滾刀磨損情況以及處理措施進行了介紹;黃和平[11]分析了深圳地鐵6號線盾構在上軟下硬花崗巖復合地層中的刀具磨損特性,發現切削軌跡與磨損量存在明顯關系,并據此建議主動檢查、換刀點。韓冰宇等[12]通過對深圳地鐵9號線現場監測的研究,對不同類型刀具的磨損系數進行了推算,并總結了影響刀具磨損的施工參數主要包括推進速度、推力、刀盤轉速、扭矩、貫入度等。

由于復合地層滾刀磨損受到施工參數、地層條件等多種因素的影響,難以建立準確的力學關系或經驗公式。而機器學習算法能夠通過挖掘數據間的關系建立起復雜的模型,其中的無監督學習方法在圍巖分級預測中得到了廣泛應用,有監督學習則被應用于盾構掘進參數、盾構姿態的預測[13]。也有學者嘗試將機器學習算法應用到滾刀的磨損預測中[8-9,12],但仍是以全斷面硬巖中的滾刀磨損預測為主,對于軟硬地層交互,掌子面上軟硬地層比例變化明顯的隧道掘進過程的滾刀磨損預測十分少見。

因此,本文基于某大直徑泥水平衡盾構隧道掘進過程中的滾刀磨損情況以及施工參數記錄、地質勘查資料,先借助無監督聚類算法kmeans對盾構施工狀態進行劃分,針對不同的施工狀態分配不同的磨損修正系數,然后將磨損修正系數連同施工參數、掌子面上地層條件以及每把滾刀的切削半徑與切削距離作為輸入參數,以滾刀的磨損量作為輸出參數,建立有監督的機器學習Transformer模型,并使用遺傳算法對磨損修正系數和Transformer模型的超參數進行優化,從而得到上軟下硬地層中的盾構滾刀磨損預測模型,借助模型實時預測滾刀磨損量,以便及時更換磨損量較大的滾刀。

1 施工參數、地層參數與滾刀參數

1.1 施工參數

根據相關研究,復合地層中影響滾刀磨損的主要施工參數包括刀盤扭矩、刀盤轉速、前艙壓力、總推力、掘進速度等。刀盤扭矩與轉速反映了掌子面前方地層對刀盤轉動的阻礙作用,一定程度上反映出滾刀與巖石的相互作用情況[14]。前艙壓力主要用于平衡掌子面上的地層壓力,其值的大小會影響刀具的破巖力,總推力與前艙壓力共同決定了刀具與巖石間的壓力[15]。推進速度則決定了滾刀的破巖效率,適當的掘進速度能夠確保滾刀高效率破巖,而推進速度不合理則會大大增加滾刀的磨損[16]。此外,在不同的滾刀磨損條件下,上述施工參數間的關系也會發生變化。例如正常情況下扭矩往往隨著總推力的增加而增加,而當刀具磨損嚴重時,即使增加總推力,扭矩也很難升高,推進速度也大大降低,盾構掘進變得十分困難。選取刀盤扭矩、刀盤轉速、前艙壓力、總推力、掘進速度作為輸入參數,其數據特征如表1所示。

1.2 地層參數

上軟下硬地層中隧道沿線掌子面上的軟硬地層比例變化往往比較大,而硬巖和軟土中滾刀的磨損速率相差很大[17],因此要把掌子面上各種地層的面積占比作為輸入參數以考慮其影響。而同一類地層的物理力學參數往往變化很小,可以看作常數,因此不作為機器學習模型的輸入參數。

1.3 滾刀參數

盾構機刀盤如圖1所示,滾刀安裝在盾構機刀盤上,隨刀盤旋轉擠壓、破碎巖體。不同滾刀的安裝半徑不同,在刀盤旋轉一圈時劃過的距離不同,安裝半徑越大,刀盤旋轉一圈時劃過的距離越長,相應的磨損也越大。因此將刀具的安裝半徑作為輸入參數。此外,刀盤的旋轉圈數決定了滾刀累積滑過的距離,因此將刀盤旋轉圈數乘以安裝半徑得到的刀具切削距離為輸入參數。

圖1 盾構刀盤與滾刀

刀具磨損量的測量如圖2所示,采用專用的測量卡尺可以得到滾刀的磨損量。滾刀的磨損分為正常磨損與偏磨,正常磨損指磨損量沿圓周均勻分布,如圖3所示,滾刀的磨損大多以正常磨損為主[6],因此輸出參數選用滾刀的正常磨損量。

圖2 滾刀磨損量測量

圖3 滾刀正常磨損

2 算法介紹

2.1 kmeans聚類算法

kmeans聚類算法是無監督機器學習算法,通過計算樣本之間的歐氏距離來對樣本集進行分類[18-19],其算法基本原理如圖4所示。

圖4 kmeans聚類算法流程

1)定義類別數量n,從樣本集中隨機選取n個點作為聚類中心。

2)計算每個樣本xi={xi1,xi2,…,xin}到每個聚類中心ci={ci1,ci2,…,cin}的距離di:

(1)

并將該樣本劃分到距離最近的聚類中心所屬的類別。

3)計算每個類別的幾何中心c′i,作為新的聚類中心:

(2)

4)重復2)和3),直到聚類中心的位置變化小于指定值或達到最大迭代次數。

作為無監督機器學習方法,kmeans不需要標簽值,能夠根據施工參數間的相關關系自動為其分類。既有的盾構施工研究表明,不同的施工參數關系反映了刀具的磨損情況,例如黃繼輝等[15]指出當掘進速度與推力都較高時,滾刀磨損速度正常;前艙壓力較高、總推力較低時滾刀磨損較快。龔永旺等[20]指出:當總推力較低、刀盤扭矩較低而推進速度正常時滾刀磨損緩慢;總推力正常、刀盤扭矩正常而速度較高時,滾刀磨損較快;總推力較高、刀盤扭矩較高而速度較低時,滾刀磨損最快。可見不同的施工參數關系下,滾刀的磨損速率不同。基于kmeans對施工參數的聚類結果,為不同類別的施工區間分配不同的滾刀磨損修正系數,并將此系數作為輸入參數。修正系數的初始值根據既有研究成果進行賦值,值越高代表磨損越快,具體值如表2所示。由于既有研究并沒有對不同施工關系下的刀具磨損量進行定量研究,因此表2定義的磨損修正系數初始值未必準確,僅僅反映了磨損速率的相對大小。所以,需要借助遺傳算法對磨損修正系數的值進行優化,以得到準確率更高的預測模型。

表2 施工參數相關關系及對應的初始修正系數

2.2 Transformer算法

由于不同滾刀安裝在不同的刀盤輻條上,且同一輻條上相鄰滾刀的切削軌跡通常并不相鄰,因此滾刀間的相互影響是一個復雜的空間問題。而Transformer獨特的自注意力機制與位置編碼方式十分適合這一問題。Transformer算法由Vaswani等于2017年提出[21],主要用于自然語言識別問題。為了將其應用于滾刀磨損問題,對其結構稍作改進,改進后Transformer模型的結構如圖5所示:輸入參數與位置編碼結合后輸入編碼器,編碼器數量為N,每個編碼器由一層自注意力層和一層前饋層組成,并且每層后面都緊接著殘差與正則化操作。經過N個編碼器后進入解碼器,解碼器數量為M,每個解碼器包含一個自注意力層、一個編碼-解碼注意力層和一個前饋層,每一層后同樣緊接著殘差與正則化操作。經過M個解碼器后數據再次經過前饋層與歸一化,最終得到預測值。計算預測值與輸出參數實際值的誤差,進行反向傳播來學習編碼器和解碼器中的網絡權重。

圖5 適合于滾刀磨損預測的改進Transformer模型結構

假設輸入參數的維度為m×1,則第i個樣本的輸入參數為ai={a1,a2,…,am},位置編碼操作時將維度同樣為m的位置編碼向量與ai相加,得到加入位置編碼后的向量bi={b1,b2,…,bm}。自注意力層包含三個矩陣:查詢矩陣WQ,鍵矩陣WK和值矩陣WV,將bi分別與3個矩陣相乘得到該樣本的查詢向量qi,鍵向量ki和值向量vi。其中查詢向量qi和鍵向量ki相乘得到的矩陣Wi就是自注意力的權重,Wi第i行第j列的元素wi,ij就代表第j個特征對第i個特征的影響權重。值向量vi則表示每個特征的初始值與位置編碼的影響,將矩陣Wi與值向量vi相乘并進行歸一化就得到了最終的自注意力結果:

(3)

其中softmax是歸一化函數。模型的訓練過程就是不斷調整查詢矩陣WQ,鍵矩陣WK和值矩陣WV中值的過程。編碼-解碼注意力層與自注意力層的區別在于其計算查詢向量qi使用的查詢矩陣WQ并不需要重新定義,而直接使用編碼器的查詢矩陣。殘差操作就是將經過自注意力層或前饋層之前的數據直接加到自注意力層或前饋層處理過后的數據上,目的是保留一定的初始信息,歸一化就是使向量的各個特征按照相對大小縮放到0~1,從而避免數值的絕對大小對結果的影響。前饋層是一個全連接的神經網絡層,將來自前一層的信息進行匯總,其權重同樣在訓練過程中更新。

從Transformer的算法原理可以看出其位置編碼步驟能夠考慮輸入參數之間的位置關系,而刀具本身的安裝半徑作為輸入參數之一,反映了滾刀間的空間位置,因此能夠借助位置編碼將滾刀的位置信息輸入進模型中。而自注意力機制可以考慮輸入參數間的相互影響,同樣能夠為對某一把滾刀影響較大的其他滾刀或施工參數分配較高的注意力值,因此Transformer模型十分適合于刀盤上滾刀磨損這類具有空間位置信息的問題預測。Transformer模型具有兩個超參數,分別為編碼器的數量N以及解碼器的數量M,這兩個超參數同樣由遺傳算法進行優化,以得到最優參數組合,確保模型的準確率最高。

2.3 遺傳算法

遺傳算法是一種高效的優化算法,通過定義染色體種群,將需要優化的參數作為染色體上的基因值,通過交叉互換和變異的方法生成新的個體,計算每個染色體的適應度,將適應度高的個體保留,將適應度低的個體剔除,經過若干代更新后,最終得到適應度最高的個體,其基因值就是優化后的參數結果。具體遺傳算法原理可以參考相關文獻[12,22]。滾刀磨損預測模型中需要優化的參數包括經過聚類得到的不同施工狀態對應的磨損修正系數以及Transformer模型的超參數,染色體的適應度為模型預測的準確率。

2.4 基于kmeans聚類與Transformer算法的預測模型

提出的基于kmeans聚類與Transformer算法的預測模型流程如圖6所示。輸入參數包括施工參數、滾刀參數、地層參數和磨損修正系數四類。其中施工參數包括總推力、刀盤扭矩、刀盤轉速、前艙壓力和推進速度這五項對滾刀磨損量影響明顯的施工參數;滾刀參數包括滾刀的安裝半徑和滾刀的切削距離;地層參數包括了掌子面上軟硬地層的面積占比;磨損修正系數由kmeans聚類算法根據施工參數聚類得到,為不同的聚類結果分配不同的磨損修正系數初始值。輸出參數為滾刀的正常磨損量。輸入參數輸入到Transformer模型中進行訓練,并借助遺傳算法對Transformer模型的超參數以及磨損修正系數進行優化,在準確率達到要求或達到最大迭代次數后得到滾刀磨損預測模型,從而對滾刀的正常磨損量進行預測。

圖6 基于kmeans聚類與Transformer算法的預測模型

3 樣本庫建立

3.1 工程信息

選取中國某東南沿海城市隧道盾構段作為本方法的應用案例。該段長2 060m,埋深約25.9~33.5m,沿線典型地質剖面如圖7所示,其中藍線為隧道結構外界線,紅線為軟土地層與硬巖地層的分界線。盾構直徑15.0m,所穿越的地層包括⑧1砂質黏土、⑩1全風化混合花崗巖、⑩2-1強風化混合花崗巖上層、⑩2-2強風化混合花崗巖下層、⑩3中風化混合花崗巖、⑩4微風化混合花崗巖。

圖7 所選工程案例地質剖面

3.2 樣本組成

對滾刀磨損量的預測采用Transformer模型,屬于有監督學習,每個樣本包含輸入參數和輸出參數(標簽)。如前所述,影響滾刀磨損的因素有很多,包括滾刀安裝半徑、切削距離、地層參數、施工參數以及kmeans聚類得到的施工狀態對應的磨損修正系數。所使用的數據來源于某大直徑泥水平衡盾構隧道前288環的掘進過程。實際施工過程中滾刀磨損量的測量方式是每掘進一段距離停機抽刀檢查,因此輸出參數就選為兩次抽刀檢查之間滾刀的正常磨損增量,則輸入參數也要選取兩次抽刀檢查之間各個因素的代表值。切削距離由兩次抽刀檢查之間刀盤旋轉的圈數乘以滾刀安裝半徑算得。

對于地層參數,以兩次抽刀檢查之間隧道沿線掌子面上各類地層的面積平均值作為輸入參數。該隧道沿線掌子面上的地層類型包括⑤1黏土、⑥1淤泥質黏土、⑥2細砂、⑥3黏土、⑧1砂質黏土、⑩1全風化混合花崗巖、⑩2強風化混合花崗巖、⑩3中風化混合花崗巖。為了減少輸入參數維度,降低模型復雜度,避免過擬合,將⑩1全風化混合花崗巖及以上地層統一定義為軟土,則地層輸入參數包括掌子面上軟土面積、⑩2強風化混合花崗巖面積和⑩3中風化混合花崗巖面積。根據地勘報告,本項目隧道沿線掌子面上3種地層的面積分布如圖8所示,圖中虛線為掌子面面積,可見隧道沿線掌子面上的三類地層占比變化明顯,既存在全斷面軟土的情況,也存在全斷面硬巖的情況。

圖8 隧道沿線掌子面上各類地層面積

對于施工參數,選取了刀盤扭矩、刀盤轉速、前艙壓力、總推力、掘進速度作為輸入參數,以兩次抽刀檢查之間的平均值作為代表值。在施工參數的聚類分析中,kmeans算法以每一環各個施工參數的平均值作為聚類點,聚類結果如圖9所示。圖9a繪制了每一環施工過程施工參數間關系的聚類結果,根據隧道沿線的地層情況以及換刀記錄,將聚類類別定為4類。圖9b則繪制了不同類別中5種施工參數的相互關系,4種顏色的圓分別對應于圖9a中的顏色相同4類聚類結果。聚類結果3,即黑色點,存在于盾構剛剛開始掘進的階段,這個階段屬于駕駛員對盾構機性能的適應階段,因此各項參數設定的相對保守,總推力、前艙壓力和刀盤扭矩都較低,刀盤轉速和掘進速度正常,滾刀屬于正常磨損情況,取磨損修正系數為1.0;聚類結果1,即紅色點,此時駕駛員對盾構機性能基本熟悉,選取了正常的推力、前艙壓力,對應的扭矩和推進速度也比較正常,滾刀屬于正常磨損情況,取磨損修正系數為1.1;聚類結果2,即藍色點,此時為了加快施工,駕駛員選取了相對激進的掘進策略,設定了較高的總推力與前艙壓力,對應的扭矩也較高,推進速度有所提高,滾刀屬于較快磨損情況,取磨損修正系數為2.5;聚類結果4,即綠色點,此時掌子面前方地層中硬巖比例迅速上升,導致刀具磨損急劇加速,駕駛員選用了正常的總推力與前艙壓力,但是對應的刀盤扭矩和推進速度都很低,滾刀屬于快速磨損情況,取磨損修正系數為3.0。

圖9 聚類結果及每一類別施工參數間相關關系

3.3 數據預處理與樣本集構建

由圖9b所示,不同施工參數由于量綱的問題導致其絕對數值相差很大,為了避免絕對數值較大的參數遮蓋了絕對數值較小的參數對滾刀磨損的影響,需要對輸入參數進行歸一化處理,選用最大值-最小值歸一化:

(4)

式中:y′為歸一化后的參數值;ymin為所有樣本中該參數的最小值;ymax為所有樣本中該參數的最大值。最終得到的樣本集尺寸為217×12,共有217個樣本,每個樣本的維度為12,包含2個滾刀參數(安裝半徑、切削距離)+3個地層參數(軟土面積、⑩2面積、⑩3面積)+5個施工參數(刀盤扭矩、刀盤轉速、前艙壓力、總推力、掘進速度)+1個磨損修正系數共11個輸入參數,和1個輸出參數(滾刀磨損量)。

4 模型建立與結果分析

4.1 模型建立

樣本總數為217,相對較少,為了更好地評價模型的性能,采用10次10折交叉驗證[22],即隨機將樣本集劃分為10份,輪流用其中的1份作為測試集,其余9份作為訓練集,計算10個測試集上的準確率的平均值,重復以上步驟10次,得到最終的準確率平均值,這個過程就稱為10次10折交叉驗證。以最終得到的準確率平均值為優化目標,采用遺傳算法對磨損修正系數和Transformer模型的超參數進行優化,得到準確率最高的結果。其中一次Transformer模型的訓練時間約為0.1s,遺傳算法優化過程約為20 000s。

4.2 模型評價指標

由于刀具的磨損屬于回歸問題,而且其影響因素復雜,當前工程中不要求嚴格的準確預測,一般在15%以內的預測誤差是可以的[13],因此選用最大容許磨損量40.0mm的15%作為容許誤差,對于測試集上的每一個樣本的真實值v和與預測值vp,如果二者滿足:

v(1-15%)≤vp≤v(1+15%)

(5)

則認為模型在這個樣本上的預測結果準確。此外,還選用平均預測絕對誤差MSE作為模型的評價指標:

(6)

式中:k為測試集中樣本數量;vi為測試集中第i個樣本輸出參數的真實值;vpi為測試集中第i個樣本輸出參數的模型預測值。準確率越高,MSE越小,模型對滾刀的磨損預測越準確。

4.3 模型預測結果

模型準確率隨遺傳算法的迭代逐漸提高,最終準確率為75%,如圖10所示,黑色線表示遺傳算法種群中表現型最好的樣本對應的準確率,紅色線表示全部種群樣本的平均準確率,可見在第4 966次迭代后最優樣本的準確率不再提高。最終最優樣本對應的4類施工狀態的磨損修正系數如表3所示,Transformer超參數取值為N=2,M=2。由表3可知,優化前后4種聚類結果的修正系數的大小關系并沒有變化,只是具體數值發生了變化,由此可見聚類結果確實反映出不同施工參數關系下刀具的磨損速度存在差異。

圖10 準確率隨遺傳算法迭代不斷提高

表3 優化前后磨損修正因子數值

將該大直徑泥水平衡盾構隧道施工過程中的刀具磨損檢測結果用于模型的驗證,最終模型的預測結果如圖11所示,模型在15%的容許誤差下的準確率為75%,平均預測絕對誤差為3.6mm。在本項目中,滾刀的容許磨損量為40mm,因此3.6mm 的誤差對于復合地層滾刀磨損預測這一復雜問題而言是可以接受的,模型的預測結果能夠為當前滾刀的磨損量提供參考,以便及時安排滾刀的更換,從而保證掘進效率與施工安全。

圖11 模型預測結果

5 結語

提出了一種結合機器學習中無監督kmeans聚類算法、有監督Transformer算法和遺傳算法的上軟下硬地層中盾構滾刀磨損預測方法,并基于收集到的某復合地層中的盾構換刀數據與相關參數,建立了滾刀磨損預測模型,模型能夠對復合地層中的盾構滾刀磨損量進行準確預測,得到如下結論。

1)無監督聚類方法能夠基于施工參數間的相關關系對施工參數進行分類,聚類結果反映了盾構的施工狀態,不同施工狀態下刀具的磨損速度有所不同,通過引入磨損修正系數的方式,為不同施工狀態分配不同的修正系數,從而將施工狀態的影響引入到預測模型中,幫助模型準確預測滾刀的磨損。

2)基于Transformer算法對復合地層中的盾構滾刀磨損量進行預測時,除了需要考慮施工參數外,還要考慮掌子面前方的軟硬地層的厚度和刀具的安裝半徑、切削距離,再加上無監督聚類得到的施工狀態修正系數,得到的模型對滾刀磨損量的預測誤差只有3.6mm,在15%的容許誤差下準確率為75%,施工時可以根據預測結果制定相應的換刀方案。

本研究仍存在不足之處:在全斷面硬巖地層中的滾刀磨損機理已有較多研究,設法將這些研究成果融入機器學習算法中,建立融合物理力學信息的機器學習預測模型,是提高模型準確率與泛化性能的關鍵,也是未來研究的主要方向。

猜你喜歡
模型施工
一半模型
后澆帶施工技術在房建施工中的踐行探索
后澆帶施工技術在房建施工中的應用
土木工程施工技術創新探討
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
土木工程施工實習的探討與實踐
扶貧村里施工忙
河南電力(2016年5期)2016-02-06 02:11:34
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 日韩免费毛片| aⅴ免费在线观看| 特级aaaaaaaaa毛片免费视频 | 日韩免费无码人妻系列| 天天躁日日躁狠狠躁中文字幕| 无码专区在线观看| 精品伊人久久久香线蕉| 亚洲 欧美 偷自乱 图片| 色综合网址| 色亚洲成人| 啪啪永久免费av| 在线观看网站国产| 中文字幕佐山爱一区二区免费| 人妻少妇乱子伦精品无码专区毛片| 中文字幕不卡免费高清视频| 国产精品久久久久久久久kt| 日韩视频福利| 国产91色在线| 国产免费高清无需播放器| 久夜色精品国产噜噜| 国产日韩欧美在线视频免费观看| 狠狠色狠狠综合久久| 欧美三级日韩三级| 日本亚洲欧美在线| 日本妇乱子伦视频| 精品自拍视频在线观看| 老色鬼欧美精品| 重口调教一区二区视频| 啦啦啦网站在线观看a毛片| 国产日韩欧美黄色片免费观看| 国产网站免费观看| 男女性色大片免费网站| 麻豆精品在线播放| 国产视频自拍一区| 色综合久久久久8天国| 超碰精品无码一区二区| 国产91无毒不卡在线观看| 一边摸一边做爽的视频17国产 | a亚洲视频| 狠狠五月天中文字幕| 欧美日韩理论| jizz在线免费播放| 色有码无码视频| 精品国产香蕉在线播出| 国产午夜一级毛片| 欧美亚洲激情| 亚洲国产精品成人久久综合影院| 亚洲综合色吧| 亚洲高清中文字幕| 国产AV毛片| 午夜小视频在线| 国产爽歪歪免费视频在线观看| 玖玖免费视频在线观看| 亚洲人成网站日本片| 国产三级视频网站| 婷婷综合亚洲| 亚洲乱伦视频| 久久精品aⅴ无码中文字幕| 强奷白丝美女在线观看| 国产毛片不卡| 日本在线视频免费| 婷婷色狠狠干| 亚洲妓女综合网995久久| 狠狠色狠狠色综合久久第一次| 久久综合一个色综合网| 亚洲人成网站观看在线观看| 国产在线专区| 成人免费一区二区三区| 六月婷婷激情综合| 欧美在线天堂| 国产自在线播放| 伊人精品视频免费在线| 久久精品国产精品青草app| 国产亚洲欧美在线视频| 97影院午夜在线观看视频| 亚洲视频一区| 国产另类视频| 国产色婷婷视频在线观看| 亚洲精品欧美日本中文字幕| 97视频在线观看免费视频| 精品午夜国产福利观看| 国产成人综合网|