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

基于隨機森林的中小跨徑公路梁橋雙柱墩地震易損性模型

2022-03-08 12:20:06魯冠亞王克海邱文華張秉哲
地震工程與工程振動 2022年1期
關鍵詞:橋梁模型

魯冠亞,王克海,邱文華,張秉哲

(1.東南大學交通學院,江蘇 南京 210096;2.交通運輸部公路科學研究院,北京 100088)

引言

中小跨徑梁橋是最常見的橋梁類型,是我國高速公路網絡的重要組成部分。汶川地震震害顯示中小跨徑公路梁橋多采用雙柱墩,橋墩易出現開裂、壓潰和傾覆等破壞,導致交通系統的中斷,對經濟產生影響[1]。在對區域或線路中的每一座中小跨徑梁橋進行建模和地震易損性評估是不可行的,因此,有必要參數化表征橋梁群組并建立數值分析模型進行易損性評估。近年來,報道的文獻已經形成了不同區域的橋梁數據庫,例如美國加州[2],歐洲土耳其[3]、希臘[4]、意大利[5]等,雖然這些研究成果建立了區域橋梁的橋墩和支座等構件的地震易損性曲線,但在建立過程中采用的主觀分類和傳統的分析方法存在諸多不足,例如:(1)采用方差或協方差分析不同橋梁類型的地震響應時,必須滿足處理效應的獨立性和方差齊次性假設[6];(2)采用概率地震需求模型時,地震需求的對數正態性假設和恒定離差假設是地震易損性分析的主要誤差來源[7],并且在高強度地震動時,橋梁構件將經歷非線性行為,導致其概率需求模型的回歸效率降低。上述假定對橋梁地震易損性的準確性造成了一定的偏差,而沒有假定則無法做任何經典統計推斷;(3)傳統方法以單參數(地震動強度(IM))為條件回歸建立易損性,無法確定其他橋梁屬性在地震易損性中的相對重要性。

隨機森林(RF)是一種流行的機器學習技術,與其他機器學習技術相比,RF是依賴于一系列不等式規則的學習算法,不對映射函數進行任何強有力的假設,并且是非參數的[8]。近年來,在工程結構的地震響應分析和性能評估中已經有所應用。例如,Tesfamariam 等[9]研究了包括RF 在內的機器學習技術在過去地震對建筑物損傷狀態分類中的應用。Jia 等[10]使用了汶川地震的橋梁數據,根據隨機森林模型進行了結構、場地和支座類型等9 個特征的重要性排序,并判斷了橋梁的地震損傷。Mangalathu 等[11]在為加州三跨和四跨橋梁建立了構件易損性曲線的過程中,采用隨機森林建立了構件多維參數下的需求模型,但未考慮構件實際的能力模型,依然采用了傳統的對數正態分布。吳曉陽[12]采用回歸決策樹方法建立了場地地震動預警參數指標,指出該快速評估方法的總體成功率達到90%以上。

文中提出了一種結合增量分析法和RF技術為特定線路中小跨徑梁橋的雙柱式橋墩建立地震易損性曲線的方法。基于上述文獻,以下2 點值得關注:(1)已有研究采用RF 技術僅建立了構件的地震需求模型,未關注構件的抗震能力模型[11];(2)中國橋梁的結構體系、材料性質、構造設計與美國橋梁存在差異,直接引用推薦值是否會導致偏差[13]。文中以一條山區高速公路為研究對象,首先采用統計工具得到其橋墩的結構、幾何和材料特性的概率分布,由拉丁超立方體抽樣(LHS)為該線路中小跨徑梁橋生成橋墩屬性數據集,建立參數化的有限元模型。然后通過Pushover分析和基于非彈性需求譜的能力譜方法獲取雙柱墩的地震需求和抗震性能數據點,檢驗需求值和能力值的對數正態分布假設,進一步利用RF算法建立雙柱墩的地震易損性模型,并與傳統方法進行比較,說明RF建立易損性曲線的優點。

1 隨機森林

隨機森林是在決策樹為基學習器構建套袋集成的基礎上,進一步在決策樹的訓練過程中引入了隨機屬性選擇[8]。在套袋中使用訓練數據重復采樣獨立地構建每棵樹,最終輸出所有樹結果的平均值進行預測,一般算法如下:

(1)隨機且有放回地從訓練集中的抽取nt個訓練樣本。

(2)隨機選取包含k個屬性子集,得到最佳分割,從每個樣本生成決策樹。

(3)通過nt個決策樹的平均值預測新數據集的輸出。

RF預測的輸出可表示為:

RF 依據計算袋外錯誤率(OOB)選擇最優的特征子集,首先使用每次重復抽樣得到樣本的生長樹,對每棵樹的袋外樣本進行預測,然后計算預測產生的OOB。

圖1展示了采用2個變量x1和x2分配預測A,B,C,D或E的單顆樹。樹起初基于變量x1劃分區域,若x1≥a1,則樹的左支被激活,A,B或C的分配取決于x2的值a2和x1的值a4。若樹的右支被激活,根據變量x2預測變量D或E。RF由隨機選取特征創建了大量的決策樹,由每個決策樹的平均預測產生輸出。

圖1 單顆樹二變量示意Fig.1 Illustration of a tree with two variables

在特征自變量的重要性判斷中,RF 采用基尼指數(GI)的降低值作為分類問題中特征重要性排序的指標,式(2)計算數據集D的基尼指數,其意義為刪除某個自變量后所有RF樹中每個特征的節點分裂不純度變化的平均值,該指標變化越大說明變量越重要。

式中:pv為數據集D中第v類樣本所占的比例(v=1,2,…,V);Dv為屬于第v類的樣本子集;a和V表示離散特征a可能具有的V個取值。

2 雙柱墩的屬性統計和數值模擬

2.1 屬性統計

文中采用新建汶川到馬爾康高速公路中小跨徑梁橋的589 個橋墩信息,該線路橋梁主要為橋面連續的預應力混凝土簡支T梁橋和預應力混凝土連續箱梁橋,簡支梁橋跨徑一般為30 m或40 m,主梁寬度為12 m;連續梁橋主梁采用單箱多室,單跨跨徑分布于20~40 m,主梁寬度分布于8.5~16.5 m。2種橋型的橋墩采用了不同的配置,前者為有蓋梁的雙柱墩,后者中墩采用無蓋梁形式(排架墩),均為圓形截面。支座采用板式橡膠支座,基礎均采用柱式樁基礎,無承臺。

簡支梁橋在每片T梁下設置板式橡膠支座,連續梁橋在排架墩2個柱頂分別設置一個板式橡膠支座。由式(3)定義板式橡膠支座剛度幾何系數,統計該值與跨徑和主梁寬度的關系,簡支梁30 m跨徑時,為2.2~2.6;40 m跨徑時,為2.6~3.2。連續梁的支座剛度幾何系數與跨徑和主梁寬度有一定的相關性,可按主梁寬度保守地取為每沿米主梁寬度的0.3~0.5倍。

式中:Ar為支座剪切面積;Σt為橡膠層的總厚度。

統計所調查線路中橋梁群組的雙柱墩結構、幾何和材料等屬性,由于橋墩的設計指標會受其他因素的影響,需要對各指標的統計數據做相關性分析,計算各指標的Spearman 相關系數,該相關系數對變量的分布不作要求,可評估2個連續變量之間的單調關系,一般認為相關系數在0.5~0.8時,兩變量顯著相關[14]。統計橋墩的蓋梁長度、柱間距與主梁寬度,結果表明,蓋梁長度和主梁寬度具有較明顯的相關性,通常二者相等;當主梁寬度為標準12 m 時,柱間距為7 m,遇到較陡橫坡地形,柱間距采用4.75 m;當主梁寬度非12 m 時,柱間距為主梁寬度的0.4~0.55 倍。對跨徑、柱高、橫向柱高比、柱間距、柱徑以及縱筋率和配箍率進行相關性分析,得出的相關系數矩陣,如圖2所示,結果表明墩柱直徑與跨徑、柱高和柱間距存在較為明顯的相關關系,不能作為獨立變量,其他指標之間無明顯相關性,可作為獨立變量。

圖2 橋墩各指標相關系數矩陣Fig.2 Correlation matrix of the pier indices

根據樁基設計參數(樁徑與樁長)和土層信息,對統計的橋墩基礎采用“m”法[15]計算其平動剛度和轉動剛度,得到基礎剛度的分布。

對于統計的離散型數據,生成相應屬性的非參數概率質量函數;對具有連續型隨機變量特性的橋梁屬性,采用Kolmogorov-Smirnov(K-S)擬合優度檢驗經驗分布是否符合某種理論分布,即通過比較理論分布和經驗分布的累積密度函數之間的最大垂直距離Dn判斷兩分布是否一致,若Dn的p值大于顯著性水平(工程中一般采用0.05),則認為兩分布一致[14]。表1總結了可視為獨立變量的橋墩屬性的統計分布和K-S檢驗p值,材料等屬性的統計分布直接應用規范或已有文獻采納的統計分布值。

表1 橋墩屬性統計分布Table 1 Statistical distributions of the pier properties

由雙柱墩的系梁設置發現,對于有蓋梁橋墩,橋墩高度大于12 m 設置1 根系梁,大于20 m 設置2 根系梁。對于排架墩,第1根系梁頂部距離柱頂0.3~0.5 m,當墩高大于12 m時設置第2根系梁。以墩柱頂部(蓋梁底部)到系梁頂部的距離與柱高的比值為指標表示系梁在墩柱的位置,該指標可滿足正態分布,見表2。系梁采用矩形截面,高度一般小于柱徑0.2 m,高寬比約為1.15。

表2 系梁位置指標分布Table 2 Distributions of the tie beam location indices

2.2 建立總體和參數化有限元模型

根據統計的橋墩屬性,采用LHS 模擬概率分布函數中的可能值生成橋墩總體。LHS 是一種分層隨機過程,提供了一種從變量分布中抽樣的有效方法[16]。在k個獨立變量X1,X2,…,Xk中,LHS 對每個變量的累積分布分為N個等概率間隔,從每個間隔中隨機選擇一個值,每個變量獲得的N個值再與其他變量隨機配對,抽樣生成樣本。由LHS生成具有95%保證率的1 200個橋墩樣本,基于有限元平臺OpenSEES[17]建立參數化的雙柱墩空間有限元模型。由彈性梁單元模擬蓋梁和系梁,彈塑性纖維梁單元模擬墩柱。混凝土和鋼筋纖維分別采用Mander 模型和雙折線模型。由于僅研究雙柱墩的易損性,所以,以集中質量點的形式考慮上部結構質量,各質量點之間以及質量點與支座頂部均采用剛臂連接,支座和樁基礎采用線性彈簧模擬,有限元模型如圖3所示。

圖3 雙柱墩數值模型Fig.3 Numerical models of double column piers

3 基于RF的雙柱墩易損性模型

3.1 雙柱墩的地震需求和能力分析

由于中小跨徑梁橋主要以一階振型控制為主,高階振型對結構的影響可忽略不計,所以,在橋墩頂部位置施加固定荷載模式,采用位移控制的擬靜力加載方法[18-19],即在橋墩頂部開始施加使墩頂產生逐級位移的單調遞增側向荷載模式,如圖4所示。在整個Pushover的過程中,監測每級加載位移下混凝土和鋼筋的應力應變、沿墩身各截面的曲率和橋墩頂部的位移,從材料和構件層面進行對應,使橋墩經歷由彈性到屈服,進入塑性以及最終破壞的全過程,標記出控制橋墩各損傷極限狀態的墩頂位移,表征雙柱墩的地震需求和能力值。

圖4 Pushover加載模式Fig.4 Loading modes of Pushover

基于非彈性需求譜的能力譜方法[20]計算所需IM 下雙柱墩的位移需求,以位移延性系數作為工程需求參數,需求譜阻尼比為5%。能力譜法的步驟如下:

(1)通過Pushover分析確定基底剪力Vb和頂部位移Dt曲線,按式(4)將Pushover曲線轉化為能力譜。

式中:A和D分別為擬加速度譜和擬位移譜;M1*為基本振型的有效質量(模態質量);Γ1為基本振型參與系數;φ1為基本振型在頂部的振幅。

(2)將上步所得的力-位移關系轉化為理想的彈塑性形式,理想化雙線性系統的彈性周期T*為:

式中,Fy*和Dy*分別為屈服力和屈服位移。

(3)確定等效單自由度系統的地震需求。首先判斷橋墩是否進入非彈性階段,若交點的譜位移小于Dy*,說明橋墩處在彈性狀態,該交點由式(4)返回求得需求點;若交點的譜位移大于Dy*,說明橋墩進入非彈性階段,則按以下步驟對彈性需求譜進行折減:

1)計算等效彈性系統與彈性需求譜的交點得到Sae以及屈服點對應的譜加速度值Say,并求出延性折減系數Rμ=Sae/Say。

2)比較T*和Tc(地震動特征周期),按下式確定位移延性系數μ:

3)由式(7)得到非彈性需求譜。

式中:Sae和Sde分別是彈性譜中(偽)加速度值和位移值;Sa和Sd分別是非彈性譜的加速度值和位移值;T為周期。

4)計算非彈性需求譜與能力譜的交點,再由式(4)返回求出需求位移值。

為了檢查橋墩地震需求的數據分布,進行了K-S檢驗,其原假設為各數據服從對數正態分布。圖5給出了0.5 g時橋墩的位移延性需求值分布,可見,K-S檢驗的p值接近于零,因此數據不服從對數正態分布。

圖5 PGA=0.5 g時橋墩的位移延性需求值分布Fig.5 Distributions of seismic demand for column displacement ductility at PGA= 0.5 g

傳統方法通常采用Hwang 的建議以鋼筋屈服、保護層混凝土剝落、鋼筋拉斷或核心混凝土壓潰定義橋墩的損傷狀態(DS),分別記為DS1、DS2和DS3,假定各損傷狀態能力模型服從對數正態分布,對應的位移延性系數分別為1.0、1.76和4.76[21],其對數標準差取0.35[2]。表3顯示了1 200個橋墩樣本的屈服位移、各損傷極限狀態位移延性系數的數據分布,采用K-S 檢驗橋墩的屈服位移和各位移延性系數是否滿足對數正態假設,K-S檢驗結果表明,各損傷極限狀態指標不能完全服從對數正態分布,且與建議值存在明顯差異,倒塌狀態位移延性系數明顯低于建議值,且各位移延性系數的對數標準差隨著橋墩損傷的變化而變化,未產生一致規律,更不是常值。所以,既有建議值是不適用于特定線路橋梁群組的橋墩性能。

表3 Pushover分析的橋墩損傷狀態性能點Table 3 Performance points of pier damage states by Pushover analyses

3.2 雙柱墩易損性分析

在傳統的易損性方法中,地震動被縮放到相同的IM進行地震需求分析,得到構件需求值的概率分布,并與構件的能力模型進行卷積得到IM處的損傷概率,其易損性函數可寫為:

式中:P為在某級IM 下,指定損傷狀態Sc下需求響應Sd的累積概率;λ為需求的對數均值。不確定性由需求響應的對數標準差(βd|IM)和損傷狀態的不確定性(βc)決定。該方法能夠計算每級IM 上的損傷概率進而生成易損性曲線。但是,該方法應在每級IM上的地震需求以及能力模型必須滿足對數正態分布。考慮到傳統方法中地震需求模型和能力模型的缺點,文中提出基于RF建立雙柱墩易損性模型,步驟如下:

步驟1:使用LHS 生成雙柱墩樣本并建立有限元模型,對1 200 個樣本在縱橋向和橫橋向上分別進行Pushover分析。

步驟2:獲取雙柱墩的位移需求,以位移延性系數作為工程需求參數。

步驟3:訓練RF 在對數空間中建立預測模型。將樣本屬性作為RF 的輸入,比較步驟1 得到的各性能點對應的位移延性系數和步驟2 獲得的位移延性系數,以找到橋墩所屬的類,輸出1 和0 的二進制向量分別表示橋墩響應超過和不超過各性能點閾值。

步驟4:基于輸入變量的概率分布隨機生成N個值(文中N取10萬),使用各自的RF模型輸出結果,計算失效概率P,即P=Nf/N,其中Nf為RF預測結果為1的個數。

步驟5:對各級IM重復步驟2到4建立易損性曲線。

文中采用PGA 作為地震動強度指標,輸入RF 模型的橋梁屬性包括主梁的跨徑(L)和寬度(Dw),橋墩類型(Type),橋墩高度與截面直徑之比(Ha·d)、橫向柱高比(rh),系梁設置(Tieb),縱筋率(ρl)和配箍率(ρw),樁基的平動(KPT)和轉動剛度(KPR),混凝土和鋼筋的類型(Con、Ste)和強度(fcon、fste)以及上部結構質量系數(MF)。有蓋梁橋墩和無蓋梁橋墩分別記為A和B;根據系梁數量,將系梁設置記為a、b和c三類。由于橋墩的其他指標與上述指標具有相關性,所以,不作為RF模型的輸入變量。

圖6比較了使用傳統方法和基于RF 的方法建立的橋墩易損性曲線,可以看出,對于DS1和DS2,與傳統方法相比,RF 方法曲線在各級IM 上的斜率較大,說明易損性模型的標準差值較小,導致在地震動強度較小時,RF 模型的損傷概率小于傳統方法的,但隨著地震動強度的增大,傳統方法將會低估橋墩的損傷概率;對于DS3,傳統方法低估了橋墩的易損性。這說明直接采用對數正態分布假設和既有建議值建立特定線路橋梁群組的橋墩易損性模型會導致較大偏差,即傳統方法對易損性模型的中值和標準差值均會產生偏差。

圖6 傳統方法和RF方法的易損性曲線比較Fig.6 Comparisons of fragility curves established by RF and traditional method

根據提出方法的步驟3,采用RF的分類精度表明RF算法準確分類的能力,精度定義為準確預測的數據與總數據的比[8]。由表4 得到,RF 易損性模型在各級IM 上的精度均在83%以上,具有較高的預測能力。此外,RF 并未對數據進行任何假設,相比于傳統方法的對數正態性假設和能力模型的推薦值,RF 方法建立的易損性模型更具可靠性。

表4 RF模型的精度Table 4 Accuracy of RF models g

由于易損性函數在計算中比散點連成的折線更方便,可采用易損性均值和標準差值描述易損性模型。因此,使用正態累積分布函數擬合損傷概率點,獲得特定損傷極限狀態的平滑易損性曲線,擬合見式(9)。

式中,λ和β分別為橋墩達到指定損傷狀態時,地震動強度(PGA)的均值和標準差。表5列出了3個損傷狀態下橋墩易損性的均值和標準差值。可以看到,雙柱墩橫橋向的抗震性能優于縱橋向的,易損性均值和標準差值隨著損傷的累積均增大。該易損性模型可直接用于研究線路橋梁中橋墩的損失估算。

表5 基于RF的雙柱墩易損性模型Table 5 Fragility models of the double column piers using RF algorithm g

3.3 橋墩屬性的相對重要性分析

分類模型除了計算特定IM下橋墩的損傷概率外,還可以用來快速識別損傷狀態。若給出特定橋墩屬性的參數值,則可以通過訓練的RF模型快速判斷橋墩在某地震動強度下的損傷狀態,無須進行重復的有限元計算模擬。選擇該線路中典型橋墩的參數配置,如表6 所示,RF 模型識別該橋墩在PGA 等于0.2、0.4、0.8 g時,分別處于橋墩屈服、保護層混凝土剝落和核心混凝土壓潰狀態。

表6 案例橋墩參數Table 6 Case pier parameters

根據式(2)計算RF 模型的GI降低值,由該值的最大最小歸一化值[8]得出每個輸入屬性對橋墩易損性的相對重要性,圖7 給出了各損傷狀態3 種不同IM(PGA 分別為0.2、0.5、0.8 g)下輸入屬性的相對重要性,得到以下結論:

圖7 橋墩屬性對易損性的相對重要性Fig.7 Relative importance of the pier attributes on fragility

(1)橋墩的幾何特性(Ha·d、rh)和類型對不同損傷狀態的易損性有重大影響。

(2)屬性的相對重要性隨PGA 的變化而變化。例如,在DS1損傷狀態,PGA 為0.2 g時,rh的影響重大,隨著PGA的增大,Ha·d的影響更加明顯。

(3)不同的損傷狀態影響橋墩屬性的相對重要性。由于橋墩非線性程度的發展,墩柱的類型、上部結構尺寸、配箍率等對易損性的影響會變大。

(4)縱、橫橋向的重要屬性具有差異。例如,DS3損傷狀態,縱橋向的易損性對rh更敏感,橫橋向則對Ha·d敏感。這主要是因為雙柱墩在縱橫向上產生損傷時的內力分布有關,在縱橋向,柱底是薄弱區;而在橫橋向,不同類型的橋墩和系梁設置導致柱子的薄弱區的變化,雙柱墩的蓋梁根部以及系梁與柱子的節點區域是薄弱區,排架墩柱底的損傷則要遠超于其他區域的損傷。

由于橋墩的類型、幾何特性、上部結構尺寸、配箍率等是墩柱的內力分布以及損傷出現的位置和發展的顯著影響因素,進而決定橋墩的易損性。

上述可見,基于RF的雙柱墩易損性模型具有以下優點:(1)RF分類模型是非參數的,沒有對橋墩的需求值和能力值進行任何假設;(2)在訓練后的RF 模型輸入特定橋墩的屬性,可以快速識別橋墩的損傷狀態;(3)在訓練RF模型的過程中,可以確定影響易損性的重要屬性。

4 結論

文中介紹了利用隨機森林算法建立特定線路中小跨徑梁橋群組的雙柱墩的易損性模型的方法。考慮了橋梁的結構、幾何和材料屬性的不確定性,根據統計分布和拉丁超立方抽樣對線路橋梁的雙柱墩進行參數化表征,批量建立了有限元模型,并與傳統條帶方法對比,結論如下:

(1)提出基于RF建立雙柱墩易損性模型的步驟。采用Pushover分析和能力譜法計算橋墩的地震需求和能力值,由隨機森林建立輸入的橋梁屬性和輸出1和0的二進制向量之間的易損性模型,該模型具有較高的預測能力,克服了傳統方法中對數正態性假設和失實的能力模型建議值的應用。

(2)經過訓練的RF 模型可快速地估計線路橋墩的損傷概率,并且在各級IM 上的精度均在83%以上。此外,若給定橋墩的屬性參數值,則可以快速識別橋墩的損傷狀態,避免重復的有限元計算,提高效率。

(3)相比于以單參數IM 建立的易損性曲線,在訓練RF 模型的過程中根據GI 的降低可直接得出橋梁輸入屬性的重要性,結果表明幾何特性和類型、上部結構尺寸、配箍率對本線路橋墩的易損性具有顯著影響。

由于該模型是完全基于指定線路建立的,因此,文中之外的橋梁樣本使用該模型時需要注意,應進一步根據提出的研究步驟訓練和評估RF模型在其他線路橋墩中的性能。

猜你喜歡
橋梁模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
手拉手 共搭愛的橋梁
句子也需要橋梁
加固技術創新,為橋梁健康保駕護航
中國公路(2017年11期)2017-07-31 17:56:30
無人機在橋梁檢測中的應用
中國公路(2017年10期)2017-07-21 14:02:37
高性能砼在橋梁中的應用
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 亚洲男人天堂2020| 综合色在线| 国产人前露出系列视频| 欧美在线伊人| 久久精品人人做人人爽电影蜜月| 国产乱人视频免费观看| 99伊人精品| 不卡无码网| 国产精品网址你懂的| 欧美精品成人一区二区视频一| 免费无遮挡AV| 国产麻豆va精品视频| 99久久99这里只有免费的精品| 亚洲精品动漫在线观看| 欧洲高清无码在线| 国产91在线|中文| 夜夜爽免费视频| 国产在线第二页| 亚洲Av综合日韩精品久久久| 97综合久久| 在线va视频| 视频国产精品丝袜第一页| 欧美精品亚洲精品日韩专| 欧美性爱精品一区二区三区 | 福利片91| 成人在线不卡视频| 91久久精品日日躁夜夜躁欧美| 国产精品成人久久| 亚洲AⅤ综合在线欧美一区| 欧美一区二区三区欧美日韩亚洲| 国产一区二区在线视频观看| 波多野结衣久久精品| www.狠狠| 国内精品视频| 国产va在线观看免费| 日韩精品无码免费一区二区三区| 精品人妻系列无码专区久久| 日本人妻丰满熟妇区| 91在线一9|永久视频在线| 国产又色又爽又黄| 国产大片喷水在线在线视频| 久久国产精品电影| 国产精品网址你懂的| 精品1区2区3区| 欧美日韩福利| 欧美成人影院亚洲综合图| 日韩欧美国产三级| 精品一区二区三区中文字幕| 国产理论精品| 99久久婷婷国产综合精| 一级毛片视频免费| 在线观看国产小视频| 日韩高清欧美| 天堂在线www网亚洲| 国产97视频在线观看| 在线视频精品一区| 麻豆精品国产自产在线| 国产成人久久综合一区| 亚洲 欧美 偷自乱 图片 | 97国产成人无码精品久久久| 91色国产在线| 国产第一页亚洲| 波多野结衣一区二区三区四区视频| 欧美午夜网站| 麻豆精品在线| 九九九精品成人免费视频7| 国产在线观看第二页| 亚洲成a人片在线观看88| 波多野结衣久久精品| 国产SUV精品一区二区6| 日日噜噜夜夜狠狠视频| 国产成人1024精品下载| 久久婷婷色综合老司机| 国产全黄a一级毛片| 国产原创第一页在线观看| 日本久久久久久免费网络| 亚洲国产亚综合在线区| 青青草原国产av福利网站| 91福利国产成人精品导航| 超薄丝袜足j国产在线视频| 久久成人国产精品免费软件| 国产美女91呻吟求|