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

基于傳感器集群BLSTM模型的結構損傷定位

2022-12-01 10:24:44韓慶華黨大智
振動與沖擊 2022年22期
關鍵詞:結構模型

韓慶華,馬 乾,黨大智,徐 杰

(1.天津大學 中國地震局地震工程綜合模擬與城鄉抗震韌性重點實驗室,天津 300350;2.天津大學 濱海土木工程結構與安全教育部重點實驗室,天津 300350;3.天津大學 建筑工程學院,天津 300350)

通過對結構表面布設的加速度傳感器的振動數據進行分析,識別結構早期損傷并及時做出維護策略,對保證大型橋梁、大跨度空間結構等公共建筑在服役期間的安全性具有重要意義。按照難易程度劃分,結構損傷識別大致分為三個層次[1]:判斷損傷與異常是否存在(L1)、識別損傷位置(L2)和損傷程度(L3)。其中對于L1層次,研究學者已開展了大量研究,并基于信號分析與統計學理論提出多種損傷判斷方法[2]。

相比于判斷損傷與異常是否存在,在無監督模式下識別損傷位置和程度(L2和L3)相對更為困難?;诠逃蓄l率和振型等模態特征的分析是目前損傷定位的主要方法。近年來研究學者已構造了模態曲率[3]、模態應變能[4]等多種損傷敏感特征用于結構損傷定位。盡管模態特征參數對損傷位置信息較為敏感,但由于模態分析過程復雜、對傳感器布置方式要求高,且固有頻率等模態參數極易受環境影響等因素[5],導致其在實際工程的復雜結構損傷定位中應用存在一定難度。

因此,直接從加速度時域數據中提取特征,從而發現損傷位置和程度得到了更多學者的關注。時間序列建模是目前影響較大、效果較好的時域特征提取方法之一。近年來,利用自回歸模型[6]、帶有外生變量的自回歸模型[7]、自回歸滑動平均模型[8]和帶有外生變量的自回歸滑動平均模型[9]等時間序列模型,研究學者從振動信號中提取不同的損傷敏感特征,用于結構損傷定位。但線性時間序列模型對數據的平穩性等要求較高,在長期監測中應用有場景有限。基于神經網絡模型對時序數據進行分析,可建立時序數據的內在非線性關聯,相較于線性時間序列模型具有更高的泛化性,Yan等[10]將子結構技術與有外生變量的非線性自回歸神經網絡相結合,并使用統計推斷標準F檢驗對結構損傷進行定位;Umar等[11]考慮鄰近傳感器信號之間的相關性,基于傳感器集群的有外生變量的非線性自回歸模型對振動數據進行預測,并基于預測殘差的標準差與概率密度峰值構造損傷指標,實現無監督的結構損傷定位。有外生變量的非線性自回歸模型的主要結構為前饋神經網絡,前饋神經網絡由于只考慮當前時刻輸入對輸出的影響,而不考慮其他時刻輸入的影響,無法有效發現較長時間序列數據的規律。

長短時記憶神經網絡是一種改進的循環神經網絡,其通過門控結構保持訓練時的梯度穩定,可以充分提取信號的時間特征。文獻[12-13]利用長短時記憶神經網絡對結構振動信號的時間依賴性進行學習,實現了不同結構損傷工況的分類。雙向長短時記憶(bidirectional long short-term memory,BLSTM)神經網絡是由前后兩個方向的長短時記憶神經網絡層對序列數據進行建模的網絡結構,在處理時序數據的分類和回歸問題中通常表現出比長短時記憶神經網絡更高的性能[14],但目前利用包括BLSTM在內的長短時記憶神經網絡結構進行損傷定位的研究均屬于有監督學習方法,需要大量有標簽樣本進行分類模型的訓練,在實際結構的損傷識別中具有很大局限性。

信息熵是表征系統信號狀態的紊亂程度和無序性的度量[15]。各種故障診斷和損傷識別的研究表明,隨著故障和損傷的發生,結構振動響應復雜性和離散程度增加,信息熵值隨之增加。陳仁祥等[16]利用信息熵與融合無量綱指標,對旋轉機械故障進行了診斷;項長生等[17]基于模態曲率效用信息熵,對梁結構損傷是否發生進行識別。

本文在上述研究背景的基礎上,提出了一種基于傳感器集群BLSTM神經網絡模型的結構損傷定位方法,以進一步提高無監督模式下的結構損傷定位精度。首先通過定義傳感器集群,利用基準工況虛擬脈沖響應函數(virtual impulse response function,VIRF)數據對不同傳感器集群的BLSTM模型進行訓練,進而通過訓練好的模型對未知工況數據進行預測,通過基準工況和未知工況的模型預測殘差的信息熵,提出了一種新的損傷識別指標,實現了結構的損傷定位與損傷程度判別。

1 本文提出的損傷定位方法

本文提出方法的核心是基于傳感器集群BLSTM模型和信息熵的結構進行損傷定位。過程中首先對不同傳感器集群的時間序列數據建立不同的BLSTM模型,然后利用基準工況的數據對BLSTM模型進行訓練,進而利用訓練好的BLSTM對未知工況的各個測點的數據進行預測,最后結合信息熵理論從模型的預測殘差中提取損傷敏感特征來進行損傷定位。方法具體實施流程如圖1所示,包括以下步驟:

步驟1對各工況下傳感器采集的加速度數據進行降噪、去離群值等預處理,并根據多自由度結構的傳感器布置特點,將傳感器分為若干集群,其中每一個集群包括一個待測傳感器以及一組與該傳感器相鄰的傳感器;

步驟2選取某參考傳感器的加速度數據作為虛擬激勵,求解各傳感器對應的虛擬脈沖響應函數,對其進行歸一化處理后作為BLSTM模型的輸入與輸出;對每個傳感器集群進行BLSTM建模,其中模型輸入為該傳感器集群中待測傳感器及與該傳感器相鄰一組傳感器的虛擬脈沖響應函數;輸出為該傳感器集群中待測傳感器的虛擬脈沖響應函數;設置合理的BLSTM參數,利用基準工況數據進行模型訓練,基于訓練好的BLSTM模型,分別對參考工況(健康工況)和未知工況(待測工況)下各個傳感器的虛擬脈沖響應函數進行預測;

步驟3計算參考工況(健康工況)和未知工況(待測工況)下各個傳感器BLSTM模型預測殘差,并基于預測殘差的信息熵構造損傷指標,基于損傷指標確定損傷位置與損傷程度。

圖1 本文方法的主要流程圖Fig.1 The main flow chart of the proposed method

下面就本文損傷定位方法中所涉及的主要理論進行介紹。

1.1 傳感器集群

在多自由度結構的損傷識別中,直接將所有測點的傳感器數據作為輸入BLSTM模型進行訓練,往往導致結構復雜、計算效率低且最終的損傷識別效果較差。由于相鄰測點的傳感器數據往往表現出相似的特征,因此可以將相鄰測點具有相似數據特征的傳感器定義為一個集群[18],通過對傳感器集群建立BLSTM模型,提高模型訓練效率和損傷定位精度。以如圖2所示的3自由度質量彈簧系統為例,簡單說明傳感器集群概念。該系統共設置了A1~A3三個測點,因此傳感器被分為三個集群,其中每個集群包含一個待測傳感器(BLSTM模型的輸出)以及一組和該待測傳感器相鄰的其他傳感器(BLSTM模型的輸入)。

圖2 3自由度質量彈簧系統的傳感器集群Fig.2 Sensor clustering of 3DOF mass-spring system

1.2 雙向長短時記憶神經網絡

循環神經網絡是以時間序列數據為輸入,在序列方向按鏈式連接的神經網絡。標準的循環神經網絡具有一種簡單的鏈式重復網絡模塊形式,例如圖3(a)所示的tanh層。圖3中:x(t)為當前輸入信息;h(t-1)與h(t)分別為上一時刻和當前時刻的輸出。循環神經網絡中簡單的鏈式結構在處理較長時序數據時,容易產生梯度消失及爆炸問題,長短時記憶神經網絡通過引入具有多個門控結構的記憶單元,使網絡具有控制內部信息積累的能力,有效解決了這一問題,其結構如圖3(b)所示。

圖3 循環神經網絡結構與長短時記憶神經網絡結構對比Fig.3 Comparison of recurrent neural network and long short-term memory neural network structures

長短時記憶神經網絡的記憶單元中核心的部分為3種門控結構,分別為:決定當前記憶單元接受或遺忘上一記憶單元信息大小的遺忘門f(t)、決定當前記憶單元儲存輸入信息大小的輸入門i(t)、決定當前記憶單元輸出信息大小的輸出門o(t),表達式分別為

f(t)=σ{Wf·[h(t-1),x(t)]+bf}

(1)

i(t)=σ{Wi·[h(t-1),x(t)]+bi}

(2)

o(t)=σ{Wo·[h(t-1),x(t)]+bo}

(3)

式中:W為權重;b為偏置項;σ為Sigmiod激活函數。記憶單元狀態通過門控值進行更新

(4)

(5)

h(t)=o(t)·tanh[c(t)]

(6)

BLSTM的網絡結構如圖4所示,前向層從上一時間步迭代到當前時間步,后向層則與之相反。圖4中,Cf(t)與Cb(t)分別為時間步t的前向層與后向層記憶單元狀態。

圖4 BLSTM網絡結構示意圖Fig.4 Schematic diagram of BLSTM network structure

向前和向后兩層長短時記憶神經網絡單元的狀態信息共同決定輸出

(7)

式中,W與b分別為權重和偏置。

通過這種雙向結構,BLSTM每個時刻的記憶單元狀態都包含前后輸入的信息,從而使網絡可以更充分的學習序列的時間依賴性。在本文中,采用BLSTM對各個傳感器集群的時域振動數據進行建模,模型的輸入為各個傳感器集群中所有測點的虛擬脈沖響應函數序列,模型的輸出為傳感器集群中待測傳感器的虛擬脈沖響應函數序列。選取基準工況的數據對各傳感器集群的BLSTM模型進行訓練,進而根據訓練好的模型,對未知工況數據進行預測。

1.3 基于信息熵的損傷指標

(8)

式中:e∈n×1為為參考工況(健康工況)或未知工況(待測工況)的BLSTM模型預測殘差;n為采樣點數。

信息熵是依賴于振動信號概率密度的統計特征,該特征根據PDF的直方圖計算,用來表征信號的不確定性和隨機性,對任意測點的模型預測殘差,其信息熵表示為

(9)

式中,p(ej)為e中第j個信息出現的概率。當結構發生損傷時,模型預測殘差的無序性和離散性增大,其信息熵增大,因此可進一步利用參考工況和未知工況的模型殘差的信息熵構造損傷指標

(10)

式中,EH和EU分別為參考工況和未知工況BLSTM模型預測殘差的信息熵。

2 試驗分析與驗證

2.1 8自由度質量彈簧系統試驗驗證

2.1.1 模型介紹

為研究損傷識別問題,洛斯·阿拉莫斯國家實驗室設計了一個如圖5所示的8自由度質量彈簧系統。該系統由7個彈簧(k1~k7)和8個鋁合金質量塊(m1~m8)連接組成,其中:k1~k7=56.7 kN/m;m1=559.3 g;m2~m8=419.4 g。所有質量塊中間開孔,可沿高度拋光的圓柱形鋼棒軸向水平滑動。使用峰值輸出為215 N的激振器對質量塊1施加水平方向隨機激勵,設置采樣頻率為500 Hz,通過加速度傳感器記錄每個質量塊在各個設計工況下的加速度響應,每組試驗采樣時間8 s,共采集4 096個數據點。選取如表1所示的6種工況的試驗數據驗證本文方法的有效性。其中BC為基準(健康)工況,激勵水平為5 V;HC1和HC2分別為激勵水平為4 V和6 V的健康工況;DC1~DC3分別是不同激勵水平(4~6 V)下的線性損傷工況,通過將k5彈簧剛度的降低14%模擬。

圖5 8自由度線性質量彈簧系統Fig.5 The 8DOF mass-spring system

表1 6種工況試驗數據Tab.1 Test data of the 6 conditions

2.1.2 傳感器集群BLSTM模型建立

如第1章所述的損傷識別流程,首先定義該系統的傳感器集群,根據傳感器數量共創建了8個不同的傳感器集群:集群1包括前兩個傳感器,其中第1個傳感器為待測傳感器,第2個傳感器為其相鄰傳感器;集群2~7由第2~7個傳感器中相鄰的3個傳感器依次組成,每個集群中間的傳感器為待測傳感器;群集8由最后兩個傳感器組成,其中最后一個傳感器為待測傳感器,第7個傳感器為其相鄰傳感器。

基于振動時序提取損傷特征時,若直接將加速度原始數據作為BLSTM模型輸入,可能會因外部激勵的不確定性給識別結果帶來干擾。為了排除這種干擾,可以利用激勵與響應數據構造脈沖響應函數,進而基于脈沖響應函數進行時間序列建模。但是在實際結構監測過程中,外界激勵通常無法獲得,從而無法求解脈沖響應函數。為了在環境激勵未知條件下獲得更好的損傷特征,選取多自由度結構某參考點的加速度數據為虛擬激勵,計算其他測點響應與虛擬激勵的之間的頻率響應函數[19],并對此頻響函數做逆傅里葉變換從而得到虛擬脈沖響應函數。理論上,計算虛擬脈沖響應函數時可以選擇任意傳感器作為參考傳感器,本文選取加速度數值相對較小的第8個傳感器為參考傳感器,將其加速度響應作為虛擬激勵,設置傅里葉變換長度為4 096,計算各傳感器的虛擬脈沖響應函數,并截取前1 000個數據,將其歸一化到[-1,1]內,其中BC工況的虛擬脈沖響應函數如圖6所示。

圖6 BC工況歸一化后的虛擬脈沖響應函數Fig.6 Virtual impulse response function of BC after normalization

分別對8個傳感器集群建立相應的BLSTM模型,所有模型均由4層網絡組成,第一層為序列輸入層,其中輸入為該待測傳感器及其相鄰傳感器的虛擬脈沖響應函數序列,即集群中所有傳感器的虛擬脈沖響應函數序列。第二層為BLSTM層,采用Adam優化算法求解,具有100個隱含單元,BLSTM層后連接全連接層與回歸層,最終將序列數據進行輸出,輸出設置為每個集群中待測傳感器的虛擬脈沖響應函數序列。設置最大迭代輪數為500,每次最小批處理數30,初始學習率0.01,利用BC工況數據對模型進行訓練。

2.1.3 損傷識別結果分析

通過BC工況數據訓練得到的各傳感器集群的BLSTM模型后,對各個工況不同傳感器的虛擬脈沖響應函數進行預測。以傳感器1和傳感器5的結果為例,說明BLSTM模型在損傷前后的數據預測效果。如圖7所示,其中每幅圖的上部分是損傷前(BC工況)的預測結果,下部分是損傷后(DC3工況)的結果。由圖7可知,對于損傷前的數據,本文的BLSTM模型預測值與實際值吻合較好,說明訓練得到的BLSTM模型能夠較高精度地預測和表征未損傷結構的振動特征;對于損傷后的數據,本文的BLSTM模型預測值與實際值雖然存在一定的偏差,但這種偏差并不明顯,通過損傷前后的預測值與實際值的直觀對比并不能識別損傷是否發生。

圖7 不同測點損傷前后實際值與模型預測值對比Fig.7 Identification results of different damage conditions

進一步求解BLSTM模型的預測殘差,并對預測殘差進行統計分析,繪制模型預測殘差的分布直方和概率密度。其中傳感器1和傳感器5的在損傷前后的模型預測殘差概率密度,如圖8所示。由圖8可知,相對于健康工況,損傷發生后模型預測殘差的無序性增大,概率分布更為離散,且在靠近損傷位置處的傳感器5的變化相對更為顯著。

按照式(9)求解各個工況的不同傳感器的模型預測殘差的信息熵,并通過參考工況與未知工況的數據求解損傷指標ID,最終結果如圖9所示。由圖9可知,在HC1和HC2健康工況下,各傳感器位置處的ID值均維持在較小水平(均值分別為0.13和0.06),說明通過BC工況數據訓練得到各傳感器集群的BLSTM模型,可對環境激勵變化下的健康工況進行正確識別;在DC1~DC3工況下,各傳感器位置處的ID相對健康工況明顯增大,且損傷附近的傳感器5和傳感器6位置處的ID明顯高于其他傳感器,本文方法可對損傷位置進行較好定位。此外對比DC1~DC3工況可以看出,隨著激勵水平的增大ID隨之增大,說明在確定的損傷程度下,ID對外荷載的較大變化也具有較好的敏感性。

圖8 不同測點損傷前后模型預測殘差的概率分布對比Fig.8 Comparison of probability distribution of residual prediction before and after damage at different measuring points

圖9 不同工況的損傷識別結果Fig.9 Damage identification results of different conditions

2.2 工字鋼梁模型試驗驗證

2.2.1 試驗過程與數據采集

為進一步研究本文損傷識別方法的適用性,設計了如圖10所示跨度為1 500 mm的工字型截面鋼梁試驗模型。鋼梁橫截面高140 mm,翼緣寬80 mm,腹板厚4 mm,翼緣厚5 mm。鋼梁兩端下翼緣底部通過螺栓固定,將鋼梁沿跨度方向劃分為(1)~(9) 9個單元,在除支座外的2~9節點上翼緣頂部布置了8個加速度器(A1~A8),傳感器靈敏度為0.5 V/g;采用激振器對鋼梁施加0.5~2 500 Hz的限帶高斯白噪聲激勵,激振器出力峰值50 N,通過M5頂桿與4節點下翼緣底部相連。設置采樣頻率為4 000 Hz,分別設置健康工況與若干損傷工況,對各工況重復試驗,采集簡支梁上A1~A8傳感器的加速度時程數據。通過在工字鋼梁下翼緣底部分別切割長30 mm,寬為5 mm,10 mm和15 mm縫隙,模擬結構輕度、中度和重度3種不同程度的損傷,設置了4種損傷工況如圖11所示。圖11中:DC1~DC3為單元(3)位置處的單一損傷工況;DC4為單元(3)和單元(5)位置處的多損傷工況。

圖10 工字鋼梁尺寸及傳感器布置圖(mm)Fig.10 Dimensions and sensor layout of I-beam test(mm)

圖11 模擬的損傷工況Fig.11 Simulated damage conditions

2.2.2 損傷識別結果分析

截取每個工況試驗數據的中間階段5 s時長加速度數據用于分析。與2.1.2節分析過程類似,按照本文的損傷識別流程,首先定義該結構的8個傳感器集群,進而選取靠近支座的A1傳感器為參考傳感器,計算不同工況各個傳感器的虛擬脈沖響應函數序列,并對其數據進行歸一化,利用健康工況虛擬脈沖響應函數對8個傳感器集群的BLSTM模型進行訓練;通過訓練好的模型對未知工況的虛擬脈沖響應函數數據進行預測并計算預測殘差,最后利用預測殘差信息熵構造的ID對損傷進行定位,結果如圖12所示。由圖12可以看出,對于不同的損傷工況,損傷發生處傳感器的ID值明顯高于其他傳感器,說明本文提出的傳感器集群BLSTM模型對單一位置損傷(DC1~DC3)和多位置損傷(DC4)均能較好的進行定位,且隨著損傷程度的增加ID隨之增大,說明本文提出的損傷識別指標同時也可較好的對損傷程度進行表征。

圖12 不同損傷工況的識別結果Fig.12 Identification results of different damage conditions

3 桁架數值模型驗證

3.1 模型介紹和數值模擬

圖13 鋼桁架試驗模型Fig.13 Model of the steel truss

使用ANSYS建立該結構的有限元模型,對本文方法的有效性進行進一步驗證。桿件采用BEAM188單元模擬彈性模量E=2.05×105MPa,材料密度ρ=7 850 kg/m3,泊松比μ=0.3。結構損傷通常表現為單元剛度的降低,本文中通過將原桿件替換成14.0×3.0 mm鋼管,模擬50%的單元剛度損失。共模擬了8種損傷工況,如表2和圖14所示。其中:DC1~DC4為單一損傷工況;DC5~DC8為多損傷工況。設置采樣頻率為200 Hz,每個工況持時20 s,共4 000個數據點。對結構整體施加垂直方向白噪聲加速度,模擬自然環境下的地脈沖激勵。采集如圖14所示的結構下弦1~13號節點處的加速度響應用于后續損傷定位。

圖14 模擬的損傷位置及傳感器位置示意圖Fig.14 Simulated damage location and sensor location schematic

3.2 損傷定位結果分析

按照第1章損傷定位流程,首先定義傳感器集群。共創建了13個不同的傳感器集群,集群1包括前兩個傳感器,其中第1個傳感器為待測傳感器;群集13由最后兩個傳感器組成,其中最后一個傳感器為待測傳感器;其余傳感器集群為中間相鄰3個傳感器。選取第1個傳感器為參考傳感器,將其加速度響應作為虛擬激勵,計算各傳感器的虛擬脈沖響應函數。分別對13個傳感器集群建立相應的BLSTM模型,建模過程與第2.1.2節所述類似,不再贅述。利用預測殘差信息熵構造的ID對損傷進行定位,損傷定位結果如圖15和圖16所示。由圖15可以看出,對于單一損傷工況,本文提出的損傷定位方法可以對損傷發生位置進行精確定位;對于多損傷工況,本文設計了兩種損傷類型,即相鄰桿件和不相鄰桿件同時發生損傷的工況。由圖16也可以看出,當損傷發生在相鄰桿件或不相鄰桿件時,本文方法均可進行較好的定位。

圖16 多損傷工況的識別結果Fig.16 Identification results of multiple damage conditions

3.3 抗噪性分析

以損傷工況DC2為例,進一步說明本文方法的抗噪性??紤]3種噪聲工況,信噪比β分別為60 dB(低噪聲)、40 dB(高噪聲)和20 dB(超高噪聲),在監測得到的原始加速度響應數據中按照如下公式添加不同程度的高斯白噪聲

(11)

為了便于分析不同噪聲水平對本文方法的影響,對每個噪聲工況不同傳感器的ID進行歸一化處理

(12)

(13)

式中,id為編號,對于DC2工況,id=4~5。

圖17 不同噪聲水平下DC2工況的識別結果Fig.17 Identification results of DC2 under different noise levels

4 結 論

提出了一種基于傳感器集群BLSTM模型的結構損傷定位方法,通過不同傳感器集群的BLSTM模型的建立,對不同位置處的虛擬脈沖響應函數進行了預測,并通過基準工況與未知工況模型預測殘差的信息熵構造損傷識別指標,依此進行損傷定位。通過8自由度質量彈簧系統試驗、工字鋼梁模型試驗,以及桁架結構數值模擬對方法進行了驗證,得到了以下結論:

(1) 通過基準(健康)工況的虛擬脈沖響應函數訓練得到的BLSTM模型,能夠較高精度地預測和表征未損傷結構的振動特征,可對環境激勵變化下的健康工況進行正確判斷,具有較好的環境適應性。

(2) 相比于健康工況,損傷發生后模型預測殘差的無序性增大,概率分布更為離散,且在靠近損傷位置處的傳感器的變化相對更為顯著,通過本文提出的基于模型預測殘差信息熵的損傷指標,可對單一位置和多位置的損傷進行準確定位,而且可通過同一位置處損傷指標的變化對損傷程度進行判別。

(3) 隨著噪聲水平的增大,損傷發生處傳感器位置的損傷指標減小,但即使對于超高噪聲水平,損傷發生處傳感器位置的損傷指標仍明顯高于其他位置,本文方法仍能較好的識別損傷位置,具有良好的抗噪聲干擾能力。

猜你喜歡
結構模型
一半模型
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
論《日出》的結構
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 91麻豆精品视频| 精品三级在线| 国产99视频精品免费视频7| aa级毛片毛片免费观看久| 国产女人水多毛片18| 91亚洲视频下载| 日韩高清中文字幕| 一级毛片在线播放| 精品久久国产综合精麻豆| 不卡的在线视频免费观看| 国产凹凸一区在线观看视频| AV老司机AV天堂| 欧美一级在线| 亚洲男人天堂2020| 免费不卡在线观看av| 婷婷六月激情综合一区| 国产爽歪歪免费视频在线观看| 成年看免费观看视频拍拍| 波多野结衣爽到高潮漏水大喷| 黄色网站在线观看无码| 特级做a爰片毛片免费69| 超清人妻系列无码专区| 亚洲综合香蕉| 国产精品黑色丝袜的老师| 精品亚洲国产成人AV| 精品国产美女福到在线直播| 亚洲一区免费看| 国产精品亚洲一区二区三区在线观看| 亚洲永久视频| 人妻无码一区二区视频| 精品福利网| 久久久久中文字幕精品视频| 欧美日韩高清在线| 中国国产高清免费AV片| 亚洲国产精品日韩专区AV| 国产精品性| 国产一二三区视频| 夜精品a一区二区三区| 国产精品极品美女自在线| 色网站免费在线观看| 亚洲精品图区| 免费国产一级 片内射老| 免费人成黄页在线观看国产| 国产手机在线小视频免费观看 | 内射人妻无套中出无码| 国产欧美综合在线观看第七页| 亚洲欧洲日韩综合色天使| 亚洲欧美一区二区三区蜜芽| 亚洲综合中文字幕国产精品欧美| 国产乱人激情H在线观看| 一本色道久久88| 72种姿势欧美久久久久大黄蕉| 黄色网站不卡无码| 午夜国产小视频| 理论片一区| 99久久精品免费看国产电影| 亚洲视频无码| 国产亚洲精品91| 99视频国产精品| 亚洲欧美成aⅴ人在线观看| 久久综合五月| 91精品国产情侣高潮露脸| 日韩毛片在线播放| 热九九精品| 日韩人妻精品一区| 亚洲日韩AV无码一区二区三区人| 国产精品亚欧美一区二区| 人禽伦免费交视频网页播放| 国产成人免费手机在线观看视频| 无码专区在线观看| 精品丝袜美腿国产一区| 国产手机在线小视频免费观看| 国产一级做美女做受视频| 婷婷激情五月网| 国产超碰在线观看| 香蕉在线视频网站| 91网在线| 丰满人妻一区二区三区视频| 中文字幕不卡免费高清视频| 性欧美久久| 国产成人综合亚洲网址| 最新国产高清在线|