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

寒區隧道凍結特性及凍結深度預測

2023-09-15 02:33:54黎忠灝邱軍領賴金星羅燕平馮志華
隧道建設(中英文) 2023年8期
關鍵詞:圍巖模型

馬 超, 黎忠灝, 邱軍領, *, 賴金星, 羅燕平, 曾 斌, 馮志華

(1. 長安大學公路學院, 陜西 西安 710064; 2. 四川川交路橋有限責任公司, 四川 廣漢 618300; 3. 河北省交通規劃設計研究院有限公司, 河北 石家莊 050011)

0 引言

隨著我國交通運輸行業的快速發展,隧道工程建設逐步向高海拔、高緯度地區發展,所面臨的環境、氣候等工程地質條件越來越復雜。在季節性凍土區,隧道施工和運營遇到了一系列嚴重的凍融病害。截至2018年底,我國共有15 117座鐵路隧道。其中,發生凍害的鐵路隧道有7 921座,占鐵路隧道總數的52.4%;發生嚴重凍害的鐵路隧道有3 855座,占鐵路隧道總數的25.5%。我國發生凍害的隧道主要集中于高緯度和高海拔地區,分布于青藏高原、東北和內蒙古等地區。寒區隧道凍害的發生與溫度場和圍巖的凍結深度有關。圍巖的最大凍結深度是影響凍脹力的重要因素,同時也是防排水設施埋設深度的重要依據[1]。

夏才初等[1]以準穩態假設為基礎,利用積分方法推導了在考慮襯砌、保溫層和凍結圍巖中未凍水含量情況下的最大凍結深度解析表達式。張晨曦等[2]對20座隧道的溫度場進行了數值模擬,結果表明,隨著地表溫度和洞外氣溫的上升,隧道的徑向凍結深度逐漸減小。丁云飛等[3]通過數值仿真得出隧道入口段凍結深度與凍結時間和外界溫度密切相關的結論,凍結時間越長凍結深度越大,在極端氣溫下洞口凍結深度可達到8.2 m。王建軍等[4]對豎井施工過程中的溫度變化進行了現場監測,得出豎井凍結深度的主要影響因素是通風時間,而風速對豎井凍結深度的影響不大。周元輔等[5]通過室內試驗對凍結深度的不同計算方法進行分析,發現采用魯基揚諾夫公式計算凍結深度比其他方法更加精確。趙鑫等[6-7]對興安嶺隧道進行了溫度實測,研究了溫度場的時空分布特征,并基于斯蒂芬法建立了凍結深度的計算方程。隨著計算機技術的飛速發展,機器學習方法已經應用到各行各業中[8-9]。一些學者將機器學習引入到隧道工程的施工與運營中。例如: 王冰泉等[10]利用集成模擬策略訓練了基準時期的支持向量機模型,用于模擬過去和未來季節凍土的最大凍結深度。

上述對寒區隧道凍結深度特性的研究,主要關注的是凍結深度與溫度、時間的關系,而對凍結深度沿隧道徑向及環向的變化規律及影響因素研究較少。針對凍結深度的預測問題,學者們多采用現場試驗與理論推導相結合的方法,而采用機器學習方法的研究較少。本文對凍結深度沿隧道徑向及環向的變化規律以及凍結深度各影響因素的顯著性及敏感性進行探究,并提出一種基于XGBoost(extreme gradient boosting)和LightGBM(light gradient boosting machine)的混合模型對凍結深度進行預測,以期為寒區隧道防凍害研究及設計提供參考。

1 工程背景

興安嶺隧道是內蒙古最寒冷地區呼倫貝爾市牙克石和阿榮旗交界處的一條重要公路隧道,其橫穿大興安嶺中南部的中、低洼山區。該地區最高海拔為1 210 m,山勢總體平緩。由于地處歐亞高緯度地區,該地區春季干燥風大,夏季溫涼短促,秋季氣溫聚降,冬季寒冷漫長。山區積雪深厚,部分地區霜凍深度超過3 m。歷年最低氣溫為-46.7 ℃,最大風速為29 m/s。隧道為雙向4車道,最高設計時速為100 km,左線隧道里程為ZK164+605~ZK168+565,全長3 960 m,交通方向由牙克石到阿榮旗。

2 溫度場及凍結深度特性分析

2.1 數值模型建立

采用有限元軟件ANSYS建立洞口段溫度場,得出襯砌及圍巖中溫度場的分布情況。隧道受環境影響的最大深度為20 m[11],因此所建立的模型上下邊界、左右邊界均距隧道襯砌20 m。模型的對稱面、恒溫邊界及對流換熱邊界如圖1所示。其中,為了清晰地顯示對流換熱邊界,將二維平面圖進行了旋轉。

圖1 計算模型(單位: m)

空氣與襯砌之間的對流換熱系數取10 W/(m2·k)[12-13]。洞內沿程的溫度荷載根據興安嶺隧道各斷面壁面溫度的監測數據確定,即將興安嶺隧道現場實測數據進行函數擬合,然后將所擬合的函數作為溫度荷載施加于隧道襯砌表面。考慮到數值模擬的結果主要用于定性研究,研究中并未具體到某一天,且為了方便建模并簡化計算,溫度荷載加載的時間均取360 d,加載步長均取1 d。模型的初始溫度是圍巖的初始地溫,綜合考慮當地氣溫條件并參考文獻[14-15],取興安嶺隧道全年實測溫度值的平均值并取整作為初始溫度,則得到模型的初始溫度為3.5 ℃。結合興安嶺隧道的相關資料及文獻資料[9-10],模型的熱力學參數如表1所示。

表1 模型的熱力學參數

2.2 數值模型合理性驗證

為驗證數值模型的合理性,提取各斷面徑深1 m處的溫度數值模擬結果與實測數據進行對比,如圖2所示。

圖2 各斷面徑深1 m處的溫度數值模擬結果與實測數據對比

通過圖2可知,溫度模擬結果與實測數據變化規律一致且變化曲線基本吻合。實測數據與數值模擬結果最大差值為1.51 ℃,平均差值為0.42 ℃,滿足工程中對數值模擬計算的誤差要求[3]。

2.3 數值模擬結果分析

2.3.1 溫度場分布規律

為進一步分析隧道圍巖徑向的溫度場分布規律,建立仰拱中心沿隧道徑向的求解路徑,提取數值模擬結果,并將各斷面徑向的溫度模擬結果繪制成如圖3所示的點線圖。

(a) 距洞口50 m

(b) 距洞口150 m

(c) 距洞口300 m

(d) 距洞口500 m

由圖3可以看出:

1)各個斷面沿徑向的溫度場分布規律與洞內氣溫緊密相關。對于氣溫高于初始地溫的8、9、10月,圍巖徑向溫度逐漸減小;對于氣溫低于初始地溫的1、2、3、4、5月,圍巖徑向溫度逐漸增加。但最終都趨于穩定,穩定后的圍巖溫度值與初始地溫相近。11月出現了反常的先增后減現象,6月、7月出現了反常的先減后增現象。這主要是由于6月、7月、11月處于正負溫交替時期,而圍巖中熱量的傳遞具有一定的滯后性,6月、7月、11月受到圍巖中原賦存熱量的影響出現了反常的先增后減或先減后增現象。

2)根據各斷面圍巖不同徑深溫度場的變化情況,并綜合考慮溫度場的變化對隧道的影響程度[2-3],可將圍巖溫度沿徑向的變化劃分為快速變化階段、緩慢發展階段、穩定階段3個階段。不同階段劃分的主要依據是圍巖溫度沿徑向的變化幅度。當每米范圍內圍巖溫度的變化幅度超過0.5 ℃時,將其界定為快速變化階段;當每米范圍內圍巖溫度的變化幅度在0.1~0.5 ℃時,將其界定為緩慢發展階段;當每米范圍內圍巖溫度的變化幅度小于0.1 ℃時,將其界定為穩定階段。則快速變化階段的范圍是沿徑深0~5 m,緩慢發展階段的范圍是沿徑深5~10 m,在沿徑深10 m之后,圍巖溫度穩定在初始地溫附近,這也表明了隧道內氣溫對圍巖溫度場的最大影響范圍為10 m。

2.3.2 凍結深度分布規律

在不考慮隧道進出口高程差、隧道內交通量以及車輛行駛速度的情況下,并以0 ℃作為圍巖是否發生凍結的判斷依據,各斷面不同位置凍結深度變化規律如圖4所示。

圖4 各斷面不同位置凍結深度變化規律(單位: m)

由圖4可知: 1)距洞口50、150、300、500 m處的最大凍結深度分別為3.99、2.61、2.42、2.35 m,相比于前一斷面,各斷面最大凍結深度降幅分別為34.6%、7.3%及2.9%。這表明距洞口越遠圍巖凍結深度越小,距洞口越近圍巖凍結深度變化越大,在距洞口300 m之后凍結深度趨于穩定。2)各個斷面不同部位凍結深度的變化規律基本一致。除仰拱部位凍結深度明顯較大外,其余各部位未出現明顯差別。這主要是由于仰拱部位相較于其他部位地下水含量高,圍巖中賦存孔隙水多;此外,地下水的流動會帶走熱量,促進圍巖的凍結。

3 凍結深度影響因素分析

3.1 正交試驗及影響因素確定

凍結深度受初始地溫、最冷月平均氣溫、襯砌和圍巖熱力學參數(熱導率、比熱容)等諸多因素的影響。為了探討影響凍結深度各種因素的顯著性、敏感性及其與凍結深度的關系,采用ANSYS軟件并結合正交試驗方法進行分析。

正交試驗本質上是一種局部試驗,在試驗中找到一組有代表性的點作為試驗對象,通過正交試驗可以減少試驗次數,降低工作量[16]。本次試驗研究的影響因素包括最冷月平均氣溫﹑初始地溫、襯砌導熱系數、圍巖導熱系數、襯砌比熱容、圍巖比熱容,每個因素設計5個水平。襯砌導熱系數、圍巖導熱系數、襯砌比熱容及圍巖比熱容以表1中的數據為基準參數,最冷月平均氣溫取距洞口50 m斷面處的監測結果,并在此基礎上對各參數分別增15%、30%和減15%、30%作為試驗的5個試驗水平,具體取值如表2所示。本試驗是以6因素5水平為標準進行正交試驗的,因此選擇L50(511)正交表進行試驗,共需進行50次試驗。

3.2 試驗結果分析

3.2.1 各影響因素敏感性與顯著性分析

正交試驗分析一般采用方差分析法或直觀分析法。直觀分析法計算簡單,但無法給出各因素的顯著性;方差分析法比較精細,計算量稍大,但可以根據方差分析結果給出各因素的敏感性及顯著性。因此,本文采用方差分析法進行分析。

方差分析主要包括計算離差平方和、計算自由度、計算平均離差平方F和顯著性檢驗4步[17]。F值越大表明該因素對試驗結果的影響越敏感。給定顯著性水平0.005,若F>F0.005(fi,fE)=F0.005(f4,f19),則該因素對結果影響顯著。其中,fi為因素自由度,fE為誤差自由度。

表2 影響因素參數取值

經過方差分析之后的試驗結果如圖5所示。從圖中F值的大小可以看出,不同因素對凍結深度的影響從大到小排序為: 初始地溫>最冷月平均氣溫>圍巖比熱容>圍巖導熱系數>襯砌導熱系數>襯砌比熱容。在顯著性水平為0.005的條件下,F(A)=345.602>F0.005(4,19)=5.27;F(B)=393.229>F0.005(4,19)=5.27;F(C)=86.119>F0.005(4,19)=5.27;F(D)=207.918>F0.005(4,19)=5.27;F(E)=35.181>F0.005(4,19)=5.27;F(F)=4.644

3.2.2 各影響因素與凍結深度之間的變化關系分析

根據所確定的影響因素以及相應的試驗水平,選用L50(511)正交表進行試驗,試驗次數共計50次。在得到試驗結果之后,求出各影響因素估算邊際平均值,以各因素的實際水平作為橫坐標,以凍結深度作為縱坐標,繪制各影響因素與凍結深度間的變化關系曲線圖,結果如圖6所示。

圖5 各影響因素F值

(a) 最冷月平均氣溫

(b) 初始地溫

(c) 圍巖導熱系數

(d) 圍巖比熱容

(e) 襯砌導熱系數

(f) 襯砌比熱容

從圖6中可以看出: 1)凍結深度與襯砌導熱系數、圍巖導熱系數呈線性正相關,凍結深度隨著這2種因素的增大而逐漸增大; 2)凍結深度與最冷月平均氣溫、初始地溫和圍巖比熱容呈線性負相關,凍結深度隨著這3種因素的增大而逐漸減小; 3)凍結深度與襯砌比熱容并未呈現出明顯的線性關系,且各水平之間的凍結深度最大差值僅為0.056 m,說明襯砌比熱容的變化對凍結深度影響不顯著。

4 凍結深度預測模型構建

由于凍結深度受多種因素綜合影響,且現場監測費用高、難度大,難以直接測得,常采用理論分析、經驗回歸公式進行預測。隨著計算機技術的不斷發展,對凍結深度的預測不再局限于統計回歸方法,而是采用SVR、XGBoost、MLP、LGBM等大數據算法進行預測。本節提出了一種基于權重分配組合模型的凍結深度預測方法。

4.1 預測模型介紹

4.1.1 XGBoost預測模型

XGBoost算法是一種設計良好的梯度增強決策樹(gradient boosting decision tree,GBDT)算法。該方法在對數據進行逐個構造的基礎上,逐步積累和綜合多個弱學習器,從而獲得較好的回歸和分類性能,其計算流程如圖7所示。XGBoost算法運算速度快、訓練效果較好,且可以有效避免過度擬合[18]。

圖7 XGBoost算法流程圖

XGBoost算法預測函數如下:

(1)

式中:yi為整個模型在這個樣本上的預測結果;k為弱評估器的總數量;fk為第k棵決策樹;xi為樣本i對應的特征向量。

XGBoost的目標函數由2部分組成: 一是模型誤差,即樣本真實值和預測值之間的差值;二是模型的結構誤差,即正則項,用于限制模型的復雜度[19-20]。目標函數

(2)

(3)

為了使得目標函數O(θ)最小,對損失函數部分進行二階泰勒展開,并求導使其為0,得到每個節點的最優預測值為

(4)

式中:ωj為第t棵樹的葉節點j對應的最優預測值;Gj和Hj分別為所有樣本在ft的j節點上gi與hi的和。

將式(4)代入目標函數(最優評分函數)得

(5)

XGBoost模型可以選擇信息增益最大的節點分裂,假設分裂前的節點為j,該節點分裂為左子節點L和右子節點R,對目標函數(最優評分函數)的貢獻為:

(6)

(7)

節點j分裂前的目標函數(最優評分函數)的貢獻和為:

(8)

因此,節點j分裂后的信息增益為:

Gain=O(R+L)-O(R)-O(L)。

(9)

式中Gain為分裂前損失與分裂后損失的差值,差值越大,代表分裂后的損失越小。

XGBoost計算Gain時遍歷所有特征所有可能的分割點,選取Gain值最大的節點進行分割。該方法的優點是精度高,但計算量太大。

4.1.2 LightGBM預測模型

LightGBM是GBDT算法的一種實現方法,該算法將Leaf-wise算法和直方圖算法相結合,使其訓練速度快、內存消耗小、準確率高,且支持分布式處理大量的數據。LightGBM計算流程如圖8所示[21]。

4.1.2.1 直方圖算法

首先,對連續的浮點特征值進行離散,得到K個分段函數,并構造寬度為K的直方圖;然后,通過對各分塊進行迭代,得到各分塊在直方圖上的累計統計量。在進行特征選取時,僅需從直方圖上求取最佳分割點即可。因此,直方圖算法具有節約儲存空間、提高計算效率等優點[22]。直方圖算法如圖9所示。

圖8 LightGBM計算流程[21]

圖9 直方圖算法

4.1.2.2 帶深度限制的Leaf-wise算法

LightGBM在分裂過程中采取了Leaf-wise生長策略,每一次選擇1片具有最大分裂能力的葉片進行分裂,具體分裂過程如圖10所示。○代表最佳分裂節點,第1次分裂以節點B作為最佳分裂節點進行分裂,第2次分裂遍歷節點D、E、C后以節點E作為最佳分裂節點進行分裂。在分裂次數相同的情況下,Leaf-wise可以降低誤差,精度較高,且LightGBM會在Leaf-wise的基礎上增加一個最大深度的限制,在保證高效率的同時防止過度擬合[23]。

圖10 帶深度限制的Leaf-wise分裂過程

4.2 建模過程

4.2.1 模型混合方式

XGBoost預測模型具有兼容中小型數據集、處理缺失值的內置函數、可以在每次迭代后運行交叉驗證等優點;LightGBM預測模型具有更快的訓練速度和更高的效率,可降低內存使用率,比其他增強算法具有更高的準確性。為了將2種模型的優勢相結合,提出了一種混合模型算法。該混合模型算法的具體步驟是:

1)利用XGBoost模型對凍結深度進行預測得到ya;

2)利用LightGBM模型對凍結深度進行預測得到yb;

3)將XGBoost模型和LightGBM模型進行融合,利用融合模型Hyb對凍結深度進行預測得到yc;

4)自適應權重部分,通過批量梯度下降法迭代訓練模型權重α、β、γ;

5)將4)中得到的權重α、β、γ分別分配給混合模型中的單模型預測值ya、yb、yc,并通過線性公式y=αya+βyb+γyc得到最終預測值y。

4.2.2 數據歸一化與反歸一化

由于數據特征因素之間的基本單位不同,需要進行歸一化處理,通過離差標準化將不同特征統一取值為0~1,計算公式為:

(10)

得到預測結果后再通過反歸一化還原為原本的量綱,其計算公式為:

y=xa(xmax-xmin)+xmin。

(11)

式(10)—(11)中:xb為原始數據;xa為歸一化后的量綱為1的數據;xmax為x的最大值;xmin為x的最小值;y為反歸一化后的數據。

4.2.3 模型評估方法

為了評估各模型的預測精度,衡量預測值與真實值之間的偏差,采用平均絕對百分比誤差(MAPE)、平方相關系數(R2)以及均方根誤差(RMSE)對模型的可靠性進行評估,其計算公式分別為:

(12)

(13)

(14)

由式(12)—(14)可看出,MAPE值越小模型的預測準確度越高;R2越接近于1,模型的擬合效果越好;RMSE越小模型離散度越低。

4.2.4 預測模型結構

凍結深度預測模型預測流程如圖11所示。計算特征選取對凍結深度影響顯著的初始地溫、最冷月平均氣溫、圍巖比熱容、圍巖導熱系數及襯砌導熱系數。相關特征融合是指將可能存在交互作用的特征進行融合,并作為新的特征參與訓練,例如: 最冷月平均氣溫與初始地溫均為環境特征,將最冷月平均氣溫與初始地溫進行融合形成最冷月平均氣溫×初始地溫特征參與訓練。

圖11 凍結深度預測模型預測流程

機器學習將數據集劃分為訓練集與測試集,訓練集用于學習,測試集用于檢驗預測效果。訓練集與測試集的數據來源于第3章正交試驗所獲取的50組試驗數據。考慮到機器學習需要的數據量較大,在此基礎上又進行了50次試驗,總計100組數據。正交試驗與所補充的50次試驗地溫、圍巖熱力學條件等參數是不相同的,以保證基礎數據的多樣性。訓練集與測試集按照7∶3的比例進行劃分,即訓練集為70組數據,測試集為30組數據。評估模型采用的指標為4.2.3節所述的MAPE、R2以及RMSE。

5 凍結深度預測結果與工程驗證

5.1 預測結果分析

XGB(XGBoost)、LGB(LightGBM)、Hyb及其混合模型預測結果如圖12所示。由圖可以看出: 與真實凍結深度值相比,XGB、LGB以及Hyb均有較好的預測結果,變化趨勢一致性較高;而LGB+XGB+Hyb模擬結果相較于其他模型來說更加穩定,擬合程度更高,幾乎與真實的凍結深度值重合。

各模型預測結果對比如圖13所示。圖中,預測結果越靠近對角線表明預測模型精度越高。則從圖中可以看出: LGB+XGB+Hyb的預測結果明顯優于其他模型,這表明提出的混合模型相較于XGB、LGB及Hyb模型適用性更強。

圖12 各模型預測結果

圖13 各模型預測結果對比

各預測方法評估指標值如表3所示。由表3可以看出: XGB+LGB+Hyb的R2值高于其他預測方法,MAPE值及RMSE值明顯小于其他預測方法,這進一步證明了所構建的預測方法具有更強的適用性。

表3 各預測方法評估指標值

5.2 工程驗證

將10余座寒區隧道的凍結深度實際值與預測值進行對比,如圖14所示。考慮到工程實踐里對凍結深度的預測可接受1個數量級范圍內的誤差[1],故在圖14中加入距實際值1個數量級的2條誤差線,認為在該誤差線范圍內的預測值達到了精度要求。從圖14中可以看出: 1)當凍結深度在1.5~2.0 m時,預測效果較好,均接近于實際值。2)當凍結深度值小于1.5 m或大于2.0 m時,預測值出現了偏差,這主要是由于在訓練所構建的XGB-LGB預測模型時,輸入的數據有73%的凍結深度值都處于1.5~2.0 m,因此使得在預測凍結深度值小于1.5 m或大于2.0 m時出現了偏差。

圖14 模型預測值與實測值對比

6 結論與討論

1)各個斷面沿徑向的溫度場分布規律與洞內氣溫密切相關。根據各斷面圍巖不同徑深溫度場的變化情況,可將圍巖溫度沿徑向的變化劃分為3個階段,即快速變化階段、緩慢發展階段以及穩定階段。

2)不同因素對凍結深度的影響程度由大至小依次是: 初始地溫>最冷月平均氣溫>圍巖比熱容>圍巖導熱系數>襯砌導熱系數>襯砌比熱容。凍結深度與圍巖導熱系數、襯砌導熱系數呈線性正相關關系,而凍結深度與最冷月平均氣溫、初始地溫、圍巖比熱容呈線性負相關關系。

3)提出的XGBoost-LightGBM模型與傳統單一模型相比預測精度較高,具有更強的適用性。

值得注意的是,不同的寒區隧道會因為氣候等原因導致凍結深度趨于穩定時距洞口的距離不同,在距洞口300 m之后凍結深度趨于穩定這一結論在類似的地質、氣候和工程條件下才具有一定的適用性。另外,在探究凍結深度各影響因素敏感性及顯著性時并未涉及到各個影響因素間的相互作用(即交互性)。此外,XGBoost-LightGBM混合模型預測凍結深度的最佳適用范圍為1.5~2.0 m,超過該范圍的凍結深度預測值會出現一定的偏差。

猜你喜歡
圍巖模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
隧道開挖圍巖穩定性分析
中華建設(2019年12期)2019-12-31 06:47:58
軟弱破碎圍巖隧道初期支護大變形治理技術
江西建材(2018年4期)2018-04-10 12:37:22
3D打印中的模型分割與打包
復雜巖層大斷面硐室群圍巖破壞機理及控制
煤炭學報(2015年10期)2015-12-21 01:55:09
滑動構造帶大斷面弱膠結圍巖控制技術
山西煤炭(2015年4期)2015-12-20 11:36:18
采空側巷道圍巖加固與巷道底臌的防治
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 无码一区二区波多野结衣播放搜索| 高清欧美性猛交XXXX黑人猛交| 日本不卡视频在线| 亚洲欧美日本国产专区一区| 无码中文字幕加勒比高清| 久久久久国产一级毛片高清板| 免费va国产在线观看| 国产成人久久综合777777麻豆| 欧美精品二区| 久久亚洲国产最新网站| 成人精品午夜福利在线播放| 谁有在线观看日韩亚洲最新视频| 激情爆乳一区二区| av尤物免费在线观看| 亚洲二区视频| 狠狠做深爱婷婷综合一区| 亚洲精品无码久久毛片波多野吉| a毛片基地免费大全| 热久久这里是精品6免费观看| 亚洲天堂日本| 亚洲欧美不卡中文字幕| 国产高清在线精品一区二区三区| 嫩草国产在线| 日韩视频免费| 1769国产精品视频免费观看| 久久天天躁狠狠躁夜夜躁| 丁香亚洲综合五月天婷婷| 三上悠亚精品二区在线观看| www.99精品视频在线播放| 中文毛片无遮挡播放免费| 日韩精品专区免费无码aⅴ| 综合五月天网| 日韩精品无码免费一区二区三区 | 午夜性爽视频男人的天堂| 亚洲欧美综合精品久久成人网| 999国产精品永久免费视频精品久久 | 亚洲精品成人福利在线电影| 在线国产91| 99精品福利视频| 狠狠色噜噜狠狠狠狠色综合久| 国产亚洲欧美在线人成aaaa| 在线无码九区| 亚洲成a∧人片在线观看无码| 先锋资源久久| 四虎亚洲国产成人久久精品| 久久国产免费观看| 成人福利在线看| 国产极品嫩模在线观看91| 国产精品久久久久婷婷五月| 蜜桃臀无码内射一区二区三区| 亚洲欧美综合在线观看| 黄色网在线| 在线日韩一区二区| 亚洲国产成人久久77| 国产成人a在线观看视频| 成人av专区精品无码国产| 永久成人无码激情视频免费| 少妇精品在线| 免费三A级毛片视频| 国产精品亚洲va在线观看| 国产精品无码影视久久久久久久| 日韩AV手机在线观看蜜芽| 亚洲精品你懂的| 欧美不卡视频在线| 国产精品视频a| 亚洲国产精品久久久久秋霞影院| 欧洲高清无码在线| 国产在线观看一区精品| 成人免费网站久久久| 制服无码网站| 日韩精品一区二区三区swag| 四虎永久在线精品国产免费| 无码'专区第一页| 亚洲美女视频一区| 九色国产在线| 囯产av无码片毛片一级| 午夜精品久久久久久久无码软件| 亚洲最新在线| 囯产av无码片毛片一级| 2019国产在线| 国产精品自在线拍国产电影| 欧美日韩在线第一页|