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

基于狀態觀測器的空間機械臂關節故障診斷

2021-03-28 02:32:46賈慶軒符穎卓陳鋼徐文倩
航空學報 2021年1期
關鍵詞:機械數據庫故障

賈慶軒,符穎卓,陳鋼,徐文倩

北京郵電大學 自動化學院,北京 100876

空間機械臂通常被用以協助或替代宇航員執行大量繁重復雜的在軌操作任務,如懸停飛行器捕獲、載荷大范圍轉移、艙體組裝等[1-2]。然而,受太空環境惡劣、操作任務繁重和自身結構復雜等因素的影響,空間機械臂極易發生關節故障[3]。當關節故障發生時,若不及時診斷并獲得相關故障信息,不僅導致空間機械臂無法執行原定任務,且可能使得故障影響擴大,威脅機械臂的運行安全。因此,亟需開展空間機械臂關節故障診斷的研究,及時檢測并診斷出故障信息,以提出相應的故障處理應對策略,保障任務的順利進行和機械臂的安全運行。

在各類故障診斷技術中,基于解析模型的故障診斷技術是較常用的方法,其包含殘差生成和診斷決策兩個過程[4]。殘差用來反映被診斷系統的狀態,當系統無故障時,殘差趨于零值,而當有故障發生時,殘差偏離零值。現有研究常利用狀態估計法實現殘差的生成,如狀態觀測器方法[5]、粒子濾波器方法[6]、卡爾曼濾波器方法[7]等,其基本思想是利用實際系統的輸出與觀測器或濾波器的估計輸出進行比較形成殘差。考慮到狀態觀測器可以有效實時估計系統狀態或輸出值,本文將通過設計狀態觀測器跟蹤空間機械臂運行軌跡,生成空間機械臂各運行狀態的殘差。

在狀態觀測器穩定跟蹤空間機械臂運行軌跡后,若觀測殘差偏離零值,則表示發生了故障。然而,受太空中高能粒子、高層大氣、航天器振動等因素的干擾[8-9],空間機械臂在軌運行時不可避免地受到外界擾動的影響,造成殘差會在零值附近波動,但并不代表發生了關節故障。因此需事先為殘差序列設定閾值,只有當殘差超出閾值時,才可判定故障發生。當故障發生時,空間機械臂關節同時受到故障信號與外界擾動的作用,基于狀態觀測器的殘差信息實則為兩者共同作用的結果,若此時殘差信息仍保持在所設定的閾值之內,則診斷系統無法給予故障的警報。針對此情況,有較多學者利用神經網絡[10]、集合論[11]、特征結構配置[12]等方法對狀態觀測器進行干擾解耦,以消除外界擾動對殘差的影響,使得殘差只對故障信號敏感。然而,實際空間機械臂系統的復雜結構很難滿足完全解耦的條件,且外界擾動形式多樣,無法得知其變化特性,導致難以實現對外界擾動的完全解耦。借鑒文獻[13]將歷史數據作為衡量標準以避免誤判或漏判的思想,本文利用系統大量的歷史運行數據進行殘差閾值的選取,若殘差不超過設定的閾值,則將故障信號視為外界擾動的一種而不給出故障發生的警報,此時也不會對系統安全及任務執行效果產生影響。若故障信號與外界擾動共同作用下造成殘差超出了閾值,才給予故障發生的警報。

在實現故障檢測后,需對故障殘差信息進行特征分析,利用合適的決策規則識別出故障模式(包括故障發生位置及故障程度),用以實現故障系統的恢復或容錯控制。故障殘差分析方法有神經網絡法、模糊推理法、統計決策法等[14]。王建國[15]針對機器人殘差信息與故障間的不確定性關系,基于模糊推理法獲得系統每個運行狀態的殘差與不同故障模式間的隸屬度,進而利用有限個運行狀態的殘差實現故障模式的識別。然而,空間機械臂是一個強耦合系統,任何一種故障模式的發生均會造成機械臂所有運行狀態的殘差發生較大的變化,此時各運行狀態的殘差與不同故障模式間均有可能存在較高的隸屬度,若只選取有限個狀態信息進行故障模式識別,容易造成其它重要信息缺失而導致故障模式的識別準確率降低。由于空間機械臂系統內部強耦合關系的存在,不同故障模式所對應的故障后系統整體運行狀態存在較大的差異,考慮空間機械臂這一特點,本文通過模擬不同故障模式下機械臂運行狀況,獲取運動過程所有運行狀態的故障殘差信息,建立包含多種故障模式的故障數據庫,進而采用相似度算法[16-17]判斷實際故障導致的機械臂故障殘差信息與故障數據庫各類別間的近似程度,實現故障模式的識別。

綜上所述,本文通過兩個階段開展空間機械臂關節故障診斷的研究,第一是設計狀態觀測器跟蹤空間機械臂運行軌跡,獲得機械臂各運行狀態的殘差信息,進而將其與設定的閾值比較,實現故障的檢測并獲得故障發生時刻;第二是通過在故障時刻引入不同故障模式,獲得空間機械臂運行狀態故障殘差信息,用以構建故障數據庫,進而將實際關節故障所導致的機械臂故障殘差信息與故障數據庫對比,識別出故障發生位置及其故障程度。本文所提空間機械臂關節故障診斷方法,可實時檢測故障的發生并獲得有效的故障信息,為故障系統的恢復或容錯控制提供指導。

1 空間機械臂關節故障檢測

檢測關節故障是否發生的思路是將空間機械臂實時運行狀態與期望的狀態進行對比,若兩者間的差異超出了閾值,則可判定發生了故障。同時,采集系統各運行狀態的故障殘差信息可作為識別故障模式的依據。本節在描述關節故障數學模型的基礎上建立關節故障空間機械臂動力學模型,進而結合滑模變結構理論設計滑模狀態觀測器,通過跟蹤空間機械臂運行軌跡,實現對關節故障的檢測并生成故障殘差信息。

1.1 關節故障空間機械臂數學模型

關節故障空間機械臂數學模型是設計狀態觀測器的基礎。在給出關節故障數學描述的基礎上,結合常態下n自由度空間機械臂動力學模型,建立關節故障空間機械臂動力學模型。

1.1.1 關節故障數學描述

按照分類方式的不同,關節故障類型可分為加性故障和乘性故障[18-19]。對于加性關節故障,其故障信號可表示為

uf=us

(1)

式中:us為關節發生加性故障后的偏移輸出。

對于乘性關節故障,其故障信號可表示為

uf=(λ-I)uθ

(2)

式中:uθ為期望關節輸出;λ=diag(λ1,λ2, …,λn)為關節有效輸出因子矩陣,其取值范圍為0≤λk≤1(k=1, 2, …,n);I為單位矩陣。

當關節同時存在加性故障與乘性故障時,關節故障信號表示為

uf=(λ-I)uθ+us

(3)

則關節的實際輸出為

ua=uθ+uf=λuθ+us

(4)

1.1.2 關節故障空間機械臂動力學模型

對于n自由度自由漂浮空間機械臂,其動力學方程為[20]

(5)

當發生關節故障時,結合關節故障模型式(3)~式(4),可得關節故障空間機械臂動力學方程為

(6)

1.2 滑模觀測器結構設計

狀態觀測器的思想是通過構造與被控系統數學模型相似的觀測器模型,將系統實際輸出與觀測器輸出對比產生輸出誤差(即殘差),進而反饋至觀測器模型,使得正常狀態下觀測器殘差趨于零。本節基于空間機械臂動力學模型,利用滑模變結構控制理論[21]設計滑模狀態觀測器,用以生成空間機械臂各運行狀態的殘差。

(7)

利用空間機械臂狀態方程,設計滑模觀測器:

(8)

(9)

然而,符號函數極易導致觀測器輸出抖振[22],為消除抖振帶來的高頻干擾,本文采用以下雙曲正切函數[23]來代替式(9)符號函數:

(10)

定義觀測器殘差為

(11)

將式(7)與式(8)代入式(11),可得殘差方程:

(12)

選取狀態觀測器的滑模面為s=e1=0,根據滑模面存在條件[21]:

(13)

(14)

此時滑模觀測器滿足吸引性條件,即

(15)

則e2=-η2tanh(e1)-k2e1也將收斂至零,滿足漸進穩定的判定條件。

1.3 故障檢測及故障殘差生成

當空間機械臂正常工作時,隨著運行時間的增加,滑模觀測器將逐漸跟蹤機械臂運行軌跡,此時殘差只受外界擾動的影響,其值較小,會在零值附近波動。當關節故障發生后,機械臂關節輸入由uθ+ud變為λuθ+ud+us,此時式(7)的空間機械臂狀態方程轉化為

(16)

由式(11)可知,觀測器滑模面被破壞,狀態觀測器將不再跟蹤實際運行軌跡,即殘差將不收斂于零,進而將殘差與閾值相比較即可實現故障的檢測。

選取系統運行狀態的角度殘差作為待比較的對象信息。首先選取角度殘差閾值rth=(rth_1,rth_2,…,rth_n+6),其中rth_j(j=1,2,…,n+6)表示空間機械臂第j個運行狀態的角度殘差閾值。當無故障時,殘差閾值的選擇取決于系統外界擾動τd的邊界,設計為

rth=f(τd)

(17)

式中:f(·)為外界擾動與殘差間的函數關系,其解析式由關節故障空間機械臂數學模型確定。

殘差閾值的選取依賴于設計人員的經驗,選取過大時會使得系統對故障的反應遲鈍,導致故障檢測漏報或者延遲,選取過小又會使得系統對故障過于敏感而導致誤報。為此,本文通過采集大量機械臂正常工作時的歷史運行數據,確定外界隨機擾動(范數有界)作用下機械臂各運行狀態殘差的極值,從而將該值設定為閾值。閾值選取的具體流程為

(18)

(19)

由式(7)與式(8),此時有:

f(ud+uf)>rth

(20)

式(20)物理意義即為當故障信號與外界擾動共同作用所造成的系統運行狀態殘差大于設定的閾值時,關節故障可被檢測。

即使系統存在故障,若故障信號與外界擾動的共同作用使得殘差仍保持在閾值以內,此時診斷系統依然視機械臂工作正常,而將故障信號視為擾動的一種,只表現為可能引發故障的潛在因素。當故障信號與外界擾動共同造成的殘差超過閾值時,則故障可被診斷系統檢測。

2 基于殘差分析的關節故障模式識別

在完成故障檢測并獲得空間機械臂運行狀態的故障殘差信息后,通過分析殘差特征可實現故障模式的識別。通過模擬不同故障模式下空間機械臂運行狀況,采集各故障模式對應的系統運行狀態故障殘差信息,建立包含多類故障模式的故障數據庫;利用相似度算法,判斷實際故障導致的機械臂故障殘差信息與故障數據庫中各類別的近似程度,篩選出與實際故障殘差相似度最大的一類,完成故障發生位置及故障程度的診斷。

2.1 故障數據庫建立

通過在仿真平臺上于故障時刻引入不同故障模式,學習各類故障模式下空間機械臂運行狀態,進而采集這些運行狀態數據作為故障樣本,用以構建故障數據庫,數學描述為

定義故障位置集為P={Jk|k=1,2,…,n},Jk表示關節編號;故障程度集為D={φl|φ=(λ,us),l=1,2,…,m},m表示故障程度總個數。取該兩個集合的元素,倆倆組合構成不同的故障模式,可表示為

S={(P1,D1),(P1,D2),…,(Pi,Dj)}

(21)

式中:Pi、Dj分別表示故障位置集P的第i個元素和故障程度集D的第j個元素,則故障模式有n×m個類別。

對于故障模式(Jk,φl),由滑模觀測器獲得該模式下空間機械臂關節與基座的角度殘差:

k_le1=[Δk_lqb,Δk_lqθ]= [Δk_lxb,Δk_lyb,

Δk_lzb,Δk_lαb,Δk_lβb,Δk_lγb, Δk_lqθ1,

Δk_lqθ2,…,Δk_lqθn]

(22)

而每個狀態分量是運行時間t的函數,即在關節故障發生后在每個運行時刻ti的采樣數據:

k_le1j=[k_le1j_t1,k_le1j_t2, …,k_le1j_tN]

(23)

其中i滿足:

(24)

式中:tf、te分別表示故障發生時刻和仿真平臺上停止數據采集的時刻;Δt表示機械臂運動規劃步長;N代表故障數據庫數據的長度。

通過遍歷故障模式集S中的每一組元素,模擬不同模式下空間機械臂運行狀況,由滑模觀測器獲得對應的空間機械臂狀態運行狀態角度殘差,即可建立包含故障發生位置和故障程度信息的空間機械臂運行狀態故障數據庫。

2.2 基于殘差相似度的故障模式識別

故障數據庫包含的故障模式有多類,當實際的空間機械臂運行狀態故障殘差與其中某一類故障模式所對應的殘差最為相似,即可識別出實際的故障模式。

設空間機械臂運行狀態實際故障殘差為

kd_lde1j=[kd_lde1j_t1,kd_lde1j_t2, …,kd_lde1j_tM]

(25)

式中:kd、ld表示待識別的實際故障發生位置和故障程度。殘差序列維度M(即數據長度)與數據采集的時間有關。在實際工程應用中,需使機械臂在發生關節故障后繼續運行一段時間,以便采集一定大小的故障信息,此時實際故障殘差的數據長度M與機械臂實際運行終止時間有關,即

(26)

式中:ts為故障后空間機械臂停止運行的時刻。一般地,實際故障殘差的數據長度較故障數據庫的短,即M≤N。M越大,代表采集的數據所包含的故障信息越豐富,更有利于實現故障模式的識別。

式(25)所顯示的殘差序列是一種特殊的既有大小又有方向的高維向量,在評估其與故障數據庫各故障模式間的相似程度時,需同時考慮這兩方面的特點進行計算近似度。為此,利用向量夾角的余弦值和歐式距離來評估實際故障殘差與故障數據庫各故障模式間的相似度,實現kd、ld的識別。

余弦相似度表征了對象在方向上的相似度,對于向量a、b,定義其間的余弦相似度為[16]

(27)

歐式距離相似度表征了對象在絕對數值上的相似度,定義其間的歐式距離相似度為[17]

(28)

dist(a,b)的值越小,表示a、b的相似度越高。

由于機械臂各運行狀態的殘差具有不同數量級的大小,需對數據進行歸一化處理以便在同一數量級下分析殘差特征。對于序列k_le1j、kd_lde1j,按式(29)進行歸一化:

(29)

利用歸一化處理后的殘差信息,基于余弦相似度和歐式距離相似度,定義實際故障殘差與故障數據庫各類故障模式間的相似度為

(30)

式(30)是基于機械臂所有運行狀態的角度殘差所構建的相似度指標,由此即可獲得維度為n×m的相似度矩陣Φ。通過尋找相似度矩陣Φ中最大的元素,將對應的行列序號k與l作為實際故障發生的位置和故障程度信息,完成故障模式的識別。

3 仿真研究

以7自由度空間機械臂為研究對象,驗證基于狀態觀測器的空間機械臂關節故障診斷方法。7自由度空間機械臂結構如圖1所示,其D-H參數和動力學參數分別如表1和表2所示。實驗仿真中,機械臂運動規劃步長為Δt=0.05 s。

圖1 7自由度空間機械臂結構Fig.1 Structure of 7-DOF space manipulator

表1 空間機械臂D-H參數Table 1 D-H parameters of space manipulator

表2 空間機械臂動力學參數Table 2 Dynamics parameters of space manipulator

3.1 關節故障檢測及故障殘差生成

利用隨機數法在0~1之間產生隨機值作為隨機變量,設外界隨機擾動力矩為

τd=(-1.5+3rand(·))sin(rand(·)Δt)

(31)

式中:rand(·)表示在[0,1)范圍內取值的隨機函數。

根據觀測器的收斂條件,可令所設計的滑模觀測器增益矩陣為

(32)

設空間機械臂初始狀態為[0 m, 0 m, 0 m, 0°, 0°, 0°, -50°, -170°, 150°, -60°, 130°, 170°, 0°],而滑模觀測器初始觀測值為[0 m, 0 m, 0 m, 0°, 0°, 0°, -45°, -165°, 155°, -55°, 15°, 175°, 5°],利用滑模觀測器跟蹤空間機械臂運行軌跡。仿真結果如圖2所示。

圖2 基于滑模觀測器的空間機械臂軌跡跟蹤Fig.2 Trajectory tracking of space manipulator based on sliding-mode observer

從圖2可以看出,空間機械臂開始運行時,由于機械臂初始狀態與觀測器初始觀測值不一致,兩者的軌跡不重合,但隨著運行時間的增加,觀測器逐漸跟蹤上空間機械臂運行軌跡,殘差逐漸趨于零,結果如圖3所示。

圖3 空間機械臂正常運行時滑模觀測器觀測殘差Fig.3 Residuals of sliding-mode observer with manipulator being at normal state

從圖3可以看出,約從6 s開始,空間機械臂各狀態故障殘差二范數的數量級保持在10-2(單位:m或(°)),表明此時狀態觀測器已穩定跟蹤空間機械臂運行軌跡。

通過多次模擬仿真空間機械臂正常運行狀況,可以獲得觀測器穩定跟蹤空間機械臂運行軌跡后(6 s后)各狀態觀測殘差信息。設置仿真實驗次數為K=50,依據式(18)選取故障殘差閾值為:

rth=[0.014 0 m, 0.011 4 m, 0.013 9 m, 0.000 2°,

0.002 2°, 0.002 1°, 0.003 9°, 0.001 1°,

0.006 0°, 0.006 4°, 0.027 6°, 0.046 0°,

0.040 2°]

(33)

假設在8 s時空間機械臂關節2發生故障,不考慮發生加性故障的情況,設其故障程度為φ=[0.003,0]。則空間機械臂各運行狀態故障殘差如圖4所示。

圖4 關節2故障時滑模觀測器殘差Fig.4 Residuals of sliding-mode observer with failure occurring at joint 2

從圖4可以看出,在8 s之前,基于滑模觀測器的空間機械臂軌跡跟蹤曲線與圖2基本一致。而從8 s開始,隨著時間的積累,空間機械臂實際運行軌跡逐漸偏離觀測器所觀測的軌跡,各運行狀態的殘差不再收斂。在8.15 s時,有

(34)

關節1與關節2的角度殘差超出了設定的閾值,表明空間機械臂在8.15 s發生了關節故障。所檢測的故障發生時刻與實驗所設置的故障時刻大體一致,說明故障檢測方法是有效的。然而,此時仍無法識別故障位置及故障程度。3.2節可通過在8.15 s引入不同故障程度的不同關節故障,構建對應的空間機械臂運行狀態故障數據庫,基于此實現故障模式識別。

3.2 關節故障模式識別

當關節3在8.15 s發生故障時,空間機械臂運行狀態的故障殘差如圖5所示。與圖4對比,關節故障后的機械臂各狀態有很大的差異,說明不同故障模式所造成的機械臂運行狀態是不同的。從這個角度出發,若能夠建立包含不同故障模式對應的空間機械臂運行狀態故障數據庫,則當實際發生某種故障時,將實際故障數據與故障數據庫對比,即可識別故障發生位置與故障程度。

圖5 關節3故障時滑模觀測器殘差Fig.5 Residuals of sliding-mode observer with failure occurring at joint 3

不考慮發生加性故障的情況,即故障程度的偏移分量us為零,選取發生乘性故障的有效因子λ范圍為[0.001,0.005],在建立故障數據庫時,設置故障程度有[0.001,0],[0.002,0],[0.003,0], [0.004,0],[0.005,0]5種可能。而故障關節發生位置均有可能發生在各關節上,即有7種可能,因此故障數據庫的故障模式有35(=5×7)類,且故障數據庫每一類模式下由13(=6+7)個運行狀態的殘差組成。取故障殘差庫的數據長度為

考慮外界擾動的影響,在構建故障數據庫每類故障模式對應的運行狀態殘差時,設置模擬實驗次數為100次,然后在每個時刻取100次模擬實驗的各運行狀態故障殘差的平均值,作為對應運行狀態最終的故障殘差,進而通過遍歷各類故障模式,完成故障數據庫的建立。

設3.1節中空間機械臂在關節2發生故障后1 s停止運動,則采集的實際故障殘差的數據長度為

令權值系數ρ=0.8,將長度為20的實際故障殘差與故障數據庫對比,根據式(30)求解其間的相似度矩陣為

Φ=

(35)

矩陣Φ中,第2行第3列的元素值最大,說明與故障數據庫中關節2發生故障程度為[0.003,0]的這一故障模式最為相似,同本文仿真實驗的初始設置一致,驗證了本文所提診斷算法的正確性。

根據式(30),不同的權值系數ρ將會影響實際故障殘差與故障數據庫各類故障模式間的相似度,從而影響診斷的結果。為探究故障診斷正確率與所選的ρ間的關系,計算不同的權值系數ρ∈[0,1]下相似度矩陣。設置每個不同的ρ對應的實驗次數為50次,實驗結果如圖6所示。

圖6 故障診斷正確率與權值系數ρ間的關系Fig.6 Relationship between accuracy of fault diagnosis and weight coefficient ρ

結合圖6與式(30)可以看出,當側重于使用余弦相似度進行故障診斷時,診斷正確率比較高;當側重于使用歐式距離相似度進行故障診斷時,則診斷正確率比較低。

仿真實驗在建立故障數據庫時只考慮了乘性故障,但并不影響故障診斷方法的有效性,即實際應用過程可以利用本文方法建立同時考慮乘性故障和加性故障的故障數據庫,以實現相關故障信息的診斷。

4 結 論

1) 在建立關節故障空間機械臂數學模型的基礎上,設計狀態滑模觀測器實現對空間機械臂運行軌跡的跟蹤,以獲得系統各運行狀態的殘差信息。

2) 通過多次模擬空間機械臂正常工作狀況選取殘差閾值,進而比較系統運行狀態殘差與閾值,實現了關節故障的檢測。

3) 通過在故障時刻引入不同故障模式,構建故障數據庫,并將實際關節故障殘差與故障數據庫進行對比,實現了故障發生位置及其故障程度的識別,完成空間機械臂關節故障診斷。

猜你喜歡
機械數據庫故障
調試機械臂
當代工人(2020年8期)2020-05-25 09:07:38
故障一點通
簡單機械
數據庫
財經(2017年2期)2017-03-10 14:35:35
機械班長
奔馳R320車ABS、ESP故障燈異常點亮
數據庫
財經(2016年15期)2016-06-03 07:38:02
數據庫
財經(2016年3期)2016-03-07 07:44:46
數據庫
財經(2016年6期)2016-02-24 07:41:51
按摩機械臂
主站蜘蛛池模板: 高清色本在线www| 在线免费无码视频| 九九九国产| 狠狠色噜噜狠狠狠狠色综合久 | 青青青国产免费线在| 久久久久无码国产精品不卡| 色成人亚洲| 五月婷婷综合色| 91福利免费| 國產尤物AV尤物在線觀看| 国产精品偷伦视频免费观看国产 | 国产香蕉97碰碰视频VA碰碰看| 成人亚洲天堂| 日日碰狠狠添天天爽| 免费毛片全部不收费的| 国产中文一区二区苍井空| 色九九视频| 欧美在线一级片| 91精品国产91久久久久久三级| 亚洲a级在线观看| 日本一区中文字幕最新在线| 欧美日韩国产在线人| 色欲不卡无码一区二区| 精品成人一区二区三区电影| 亚洲另类色| 免费99精品国产自在现线| 91精品小视频| 中文无码精品a∨在线观看| 国产成人综合欧美精品久久| 囯产av无码片毛片一级| 国产激情第一页| 小蝌蚪亚洲精品国产| 国产成年女人特黄特色毛片免 | 久久精品午夜视频| 亚洲日韩国产精品综合在线观看| 日韩资源站| 亚洲天堂网在线播放| 色综合五月婷婷| 亚洲成人在线免费观看| 欧美综合中文字幕久久| 精品第一国产综合精品Aⅴ| 爽爽影院十八禁在线观看| 久久精品视频亚洲| 亚洲Aⅴ无码专区在线观看q| 伊人久综合| 亚洲精品第一页不卡| 亚洲AV一二三区无码AV蜜桃| 激情六月丁香婷婷四房播| 亚洲 欧美 日韩综合一区| 日本欧美精品| 国产精品亚洲一区二区在线观看| 97国产在线视频| 亚洲欧美日韩另类在线一| 欧美激情成人网| 91小视频在线观看免费版高清| 国产成人精品日本亚洲| 国产黑丝一区| 国产xxxxx免费视频| 国产日韩欧美精品区性色| 98精品全国免费观看视频| 极品私人尤物在线精品首页| 国产色爱av资源综合区| 国产 日韩 欧美 第二页| 国产91av在线| 亚洲视频色图| 国产第一页第二页| 免费国产小视频在线观看| 欧美性精品| 天天摸天天操免费播放小视频| 日本在线免费网站| 国产69囗曝护士吞精在线视频| 国产欧美日韩精品综合在线| 九色视频最新网址| 亚洲男人的天堂在线| 沈阳少妇高潮在线| 毛片大全免费观看| 黄色a一级视频| 国产特级毛片aaaaaa| 国产女同自拍视频| 国产网友愉拍精品| 免费又黄又爽又猛大片午夜| 人人91人人澡人人妻人人爽|