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

基于SSWD優化LMD的高速電梯滾動導靴振動信號特征提取

2021-06-03 03:24:18許有才魯云波喬王治楊春宇張俊喃
振動與沖擊 2021年10期
關鍵詞:特征提取電梯細節

陶 然, 許有才, 和 杰,3, 魯云波,3, 喬王治,3, 楊春宇,3, 張俊喃,3, 李 珺, 王 華

(1. 昆明理工大學 質量發展研究院,昆明 650500;2.云南省特種設備安全檢測研究院,昆明 650228;3.云南惠民勞務服務有限公司,昆明 650228)

導靴與導軌共同組成電梯的導向裝置,其主要作用是將轎廂固定在導軌上做上下運動,其狀態好壞直接影響電梯的運行狀況。當滾動導靴失效時,導靴滾輪不能跟隨轎廂移動同時滾動。主要原因有滾輪磨損、老化、腐蝕,彈簧失效,調節螺絲松動等。

2016年10月在云南省昆明市海倫國際小區發生“電梯門爆炸”事件,該事件被都市晚高峰曝光,引起了社會各界的廣泛關注。據云南省特種設備安全檢測研究院調查結果表明,導致該事故發生的原因是:電梯轎廂導靴故障導致轎廂門刀撞上層門地坎。故障診斷主要由3個步驟組成:故障特征提取、計算特征參數構建特征集合、模式識別。故障特征提取作為故障診斷的第一步,也是最重要的一步,完整地提取故障特征,是高效模式識別的前提和基礎。

當導靴發生故障時,其振動信號會表現出明顯的非線性、非平穩性。對于非線性、非平穩振動信號,需要采用合適的時頻分析方法提取故障特征分量。常用的時頻分析方法有小波變換、經驗模態分解(empirical model decomposition, EMD)、LMD(local mean decomposition)等[1-3]。但是小波變換的小波基一旦選定就無法改變,缺乏自適應性[4-5]。EMD方法雖然具有自適應性,但是存在過包絡、欠包絡、端點效應和無物理意義的負頻率等局限性。相比EMD而言,LMD雖然具有端點效應不明顯、虛假分量少,不會出現負頻率的優點,但其分解得到的乘積函數(product function, PF)分量中仍然存在相近的特征時間尺度分量被分解到兩個及兩個以上PF分量中,即模態混淆現象[6-7]。

文獻[8]通過對比研究指出,LMD在保留信號信息完整性、迭代次數少及抑制端點效應等方面要優于EMD方法。文獻[9]將LMD方法應用到齒輪故障診斷中,得到了故障特征信號;但是該故障特征信號分散在多個PF分量中,存在模態混淆現象,故障診斷精度較低。

小波分解(wavelet decomposition, WD)是Mallat 在Mallat 算法基礎上提出來的一種多分量信號分解方法,可以在強噪聲背景下準確捕捉到多組不同頻率、不同衰減系數的沖擊響應信號。該方法首先構造高通濾波器、低通濾波器;然后在頻率多分辨的基礎上通過帶通濾波器,按魚骨算法原則將原始信號分解成高頻細節信號與低頻近似信號的一種特征分離方法。廣泛應用于圖像處理、降噪、濾波等信號處理領域。

文獻[10]將小波分解的多分辨率濾波特性與Elman網絡相結合應用于地下水埋深預測耦合模型,將地下水埋深序列分解成高頻細節信號和低頻近似信號,再將低頻近似信號輸入到Elman網絡,有效提高了地下水埋深預測精度。文獻[11]利用小波分解、重構提取原始數據的形態特征,然后再進行樣本分類,將趨勢與階躍異常狀態轉化為近似線性可分的類別,大大降低了趨勢類與階躍類之間的誤判率。文獻[12]提出了一種EMD與Laplace小波相結合的方法應用于齒輪故障診斷,該方法采用EMD抑制了噪聲和其他干擾信號以突出故障特征信息,但是EMD不能一次提取到故障特征信息,仍存在模態混淆現象,即一個特征時間尺度分量被分解到兩個及兩個以上固有模態函數中。

電梯導靴按照運動方式不同可以分為滑動導靴與滾動導靴。滑動導靴主要應用于低速貨梯、低速客梯轎廂側或對重側;滾動導靴主要應用于高速電梯的轎廂側。文獻[13]根據Kalker三維Hertz滾動接觸簡化理論建立了滾動導靴-導軌的三維滾動接觸模型及高速曳引電梯系統動力學方程,很好地預測了高速曳引電梯的振動響應;但是其沒有研究振動響應與設備運行狀況之間的關系。文獻[14]針對電梯的5種常見機械故障,提出了一種基于遺傳算法優化支持向量機的電梯故障診斷方法,采用最優小波包計算故障振動信號的能量分布,并將能量分布與時域指標相結合構造故障特征向量對故障類型進行識別,取得了較高的識別效率。文獻[15]利用大數據分析手段對電梯轎廂振動信號進行特征參數提取,綜合運用聚類分析與回歸分析方法,實現了對電梯導靴的故障診斷;但該方法直接從振動信號中提取特征參數構建特征向量,沒有對振動信號進行振動特征提取,只能識別正常和故障兩種類型,有待進一步研究細化故障類型。文獻[16]以速度小于或等于1 m/s的客梯、貨梯滑動導靴為研究對象,提出了一種基于SVD(singular value decomposition)優化LMD的電梯導靴振動信號故障特征提取方法,該方法采用SVD抑制原始信號中的噪聲、干擾信號及平穩信號以達到凸顯高頻細節特征分量的目的,再使用LMD進行故障特征分量提取,該方法取得了很好的效果;但是將該方法應用于高速電梯(6 m/s)滾動導靴振動信號故障特征分量提取,由于SVD在抑制原始信號中噪聲、干擾信號及平穩信號的同時也會對高頻細節特征分量(故障特征分量)產生一定程度的抑制,因此仍然存在模態混淆現象。綜上,電梯導靴運行狀態的診斷已取得很大進展,但是關于高速電梯滾動導靴振動信號的故障特征提取,需要采用合適的分析方法進行深入研究。

針對高速電梯滾動導靴振動信號故障特征分量提取,本文將WD多分辨率濾波特性、信號增強技術與LMD的特征提取優勢相結合,提出了一種基于SSWD優化LMD方法。該方法首先構造高、低通濾波器、小波函數、尺度函數,利用WD的多分辨率濾波特性將原始信號分解為高頻細節特征信號和低頻近似信號;然后對高頻細節特征信號進行增強,將增強后的高頻細節特征信號與低頻近似信號進行重構;最后采用LMD對重構信號進行特征提取。通過與陶然等研究中的方法進行對比,表明該方法完整地提取了滾動導靴振動信號的故障特征分量,避免了模態混淆現象出現,為高速電梯滾動導靴的故障診斷提供了一個有效的參考方法。

1 SSWD與LMD的基本理論

電梯在正常運行期間,需要諸多機械、電氣部件的協同配合(如,曳引機、曳引繩、排風機等),因此,采集的振動信號具有強噪聲背景的特點。小波分解由于其強大的多分辨率濾波特性,成為了故障特征提取的重要分析方法之一。小波分解通過濾波器將信號分解成不同尺度的一系列小波基函數之和,因此小波分析的效果依賴于小波基函數的選擇。

1.1 構造高、低通濾波器

dbN小波是世界著名的小波分析學者I.Daubechies構造的小波函數,簡寫為dbN。當N=1時,為Harr小波。除Harr小波之外,dbN不具有對稱性,其低通濾波器與高通濾波器表達式

(1)

(2)

(3)

(4)

式中:N=2L為濾波器長度;ω為頻率;e為自然底數;H為低通濾波器;G為高通濾波器。

H,G為一組頻率均分帶通濾波器,其頻率能夠覆蓋整個原始信號的頻譜范圍。即將原始信號在其頻率范圍內被均分為高頻細節特征信號與低頻近似信號兩部分。

1.2 小波分解

原始信號S(t)第j層的一維離散小波分解可以表示為

(5)

(6)

(7)

式中:cDj為高頻細節系數;cAj為低頻近似系數,(j=1,2,…,L;k=0,1,…,N/2j-1),N為S(t)的長度,L為分解層次;P0S=S;g(n)為高通濾波器;h(n)為低通濾波器;ψ(t)為小波基函數;φ(t)為尺度函數。

由于滾動導靴振動信號既含有從中間向兩邊衰減的“魚腹狀”波形信號又含有單邊衰減的沖擊響應信號,因此本文將dbN小波引入信號分析。除Harr小波之外,dbN的小波基函數和尺度函數沒有明確的表達式,且不具有對稱性。db10的小波函數和尺度函數,如圖1所示,由圖1可知小波函數表現出了從中間向兩邊衰減的“魚腹狀”波形。

高頻細節系數即信號的高頻細節特征,如果以此為中心,經過變換上式可表示為

S=Q1S+Q2S+…+QjS+PjS

(8)

上式中的相加,不是簡單的數值相加而是以截斷的形式相加。

圖1 db10小波Fig.1 db10 wavelet

1.3 信號增強

圖像增強是指采用某種方法將圖像中的某些感興趣的細節分量凸顯出來[17]。本文將圖像增強的方法應用到電梯導靴振動信號中以凸顯原始信號中的高頻細節特征分量。

信號增強的必要性:①對于非平穩、非線性信號是在復雜的噪聲背景下采集得到的,其中很多故障特征分量被噪聲及其它干擾信號淹沒,如果不對其進行信號增強,很難凸顯出來;②通過信號增強,高頻細節特征分量在原始信號中的幅值、相對質量將會增大,再采用LMD方法進行特征提取時,能夠使第一個PF分量中含有更多的高頻細節特征信號,改善模態混淆現象。

信號增強的基礎是確定信號各點鄰域強度的變化值,然后將信號點鄰域值有顯著變化的點凸顯出來。該算法的數學表達式為

X(n)=((x(n)-x(n-1))2+(x(n)-x(n+1))2)1/2

(9)

1.4 局部均值分解

LMD方法是科研工作者Jonathan S.Smith于2005年,在EMD方法的研究基礎上提出的。LMD方法是一種基于局部時間尺度的特征提取方法,即能從原始振動信號中提取包含局部振動特征的PF分量。

LMD的實質是從原始振動信號中分離出純調頻信號和包絡信號,將純調頻信號和包絡信號相乘便可以得到具有物理意義的PF分量,迭代處理至所有的PF分量分離出來。對于任意信號x(t),其處理結果如式(10),具體分解步驟參見程軍圣等的研究。

(10)

式中,u(t)為殘余分量。

1.5 瞬時Teager能量算子

傳統意義上的信號能量是幅值的平方,代表勢能或動能。其雖然能夠表征信號的瞬時幅值振動特征,但是對于沖擊幅值較小,則該振動模態很容易被其它模態淹沒[18]。

Teager-Kaiser能量算子(Teager-Kaiser energy operator, TKEO)是由Kaiser提出的一種非線性算子。該算子由信號瞬時幅值及其微分通過非線性組合成信號的總能量[19]。與傳統信號能量相比,該算子能更有效地提取信號的瞬時特征,并具有良好的時間分辨率,近年來被學者廣泛應用于求取信號的瞬時振動特征。

離散信號x(n)的TKEO定義為

ψd[x(n)]=x2(n)-x(n-1)x(n+1)

(11)

由式(11)可知,計算任意離散時間信號的TKEO,只要該信號具有三個樣本點即可計算任意點處的瞬時能量。因此,TKEO對信號的瞬時變化具有良好的時間分辨率,能夠高精度地檢測信號中的瞬態振動特征。

2 基于SSWD優化LMD的特征提取

基于以上論述,本文提出的基于SSWD優化LMD的高速電梯滾動導靴振動信號特征提取方法,首先構建低通濾波器、高通濾波器、小波基函數、尺度函數,利用WD的多分辨率濾波特性將原始信號分解為高頻細節特征信號和低頻近似信號;然后對高頻細節特征信號進行信號增強,將增強后的高頻細節特征信號與低頻近似信號進行重構;最后采用LMD從重構信號中提取能夠表征滾動導靴振動特征的PF分量,求取PF分量的瞬時Teager能量進行對比分析。

該方法具體步驟如下:

步驟1構建低通濾波器H、高通濾波器G、小波基函數ψ(t)、尺度函數φ(t);

步驟2對原始信號x(t)進行一層WD,將得到高頻細節特征分量D(t)和低頻近似分量A(t);由于H,G為一組頻率均分帶通濾波器,即D(t)的頻譜范圍為[1/2fmax,fmax],fmax為X(t)頻譜的最大值;則D(t)包含原始信號的高頻細節特征信號,因此,只需要進行一層WD即可滿足實驗分析要求;

步驟3對高頻細節特征分量D(t)進行信號增強,然后采用WD逆運算方法對增強后的高頻細節特征信號與低頻近似信號進行重構,得到凸顯高頻細節特征分量的重構信號X′(t);

步驟4對X′(t)進行LMD,得到能夠表征原始信號振動特征的分量PFi;

步驟5根據TKEO理論求取PFi的瞬時Teager能量,根據Teager能量波形分析故障特征提取效果。

3 實驗研究

為了驗證SSWD優化LMD方法的有效性,將采用上述原理和方法分別對高速電梯轎廂側的上導靴故障信號、下導靴故障信號進行分析。本文采用的數據來源于云南省特種設備安全檢測研究院。

實驗采用SD011型滾動導靴,如圖2所示,其詳細參數如表1所示。在不影響電梯正常運行的情況下,通過調節螺絲分別使上、下導靴滾輪與導軌之間的間隙為1.2 mm,轎廂為空載運行,運行區間為底層站至頂層站,共70層69站69門,約310 m。該電梯的額定載質量為1 800 kg,額定速度為6.0 m/s。

1.滾輪; 2.滾軸; 3.輪臂; 4.軸承; 5.彈簧; 6.靴座。圖2 SD011型滾動導靴Fig.2 SD011 type rolling guide shoe

表1 滾動導靴參數

滾動導靴的振動信號由安裝在轎廂地面的加速度傳感器MPU6050進行采集,采樣頻率為100 Hz,波特率為115 200即每隔10 ms輸出1幀數據,電梯結構簡圖如圖3所示。采用USB轉TTL電平的串口模塊與計算機連接,具體連接方法如圖4所示 。

4 數據分析

實驗數據采集均在轎廂以額定速度運行時,3種狀態振動信號時域圖如圖5所示,瞬時Teager能量波形如圖6所示。從圖5可以看出,導靴處正常狀態運行時,其振動信號在時域圖的振幅遠小于故障振動信號。從圖6也可以看出,導靴正常狀態運行時,其不僅沒有明顯的沖擊特征,而且與故障振動信號相比,其能量譜峰值相差約100倍。從理論角度分析,導靴處于正常狀態運行時,其振動信號主要由平穩振動信號、其它干擾信號、噪聲組成,不具備明顯的沖擊特性。

圖3 電梯結構簡圖Fig.3 Elevator structure diagram

圖4 振動信號采集模塊Fig.4 Vibration signal acquisition module

綜上所述,將理論分析與現場工況數據表現相結合可以得到,導靴處正常狀態運行時,其振動信號不具有沖擊響應特性;因此,如果將本文方法應用于導靴正常狀態振動信號特征提取無法體現本方法的有效性,也沒有實際意義。所以本文僅對兩種故障振動信號進行分析;同時由于篇幅限制,本文僅選取一組轎廂下導靴故障振動信號進行詳細分析。將分別采用SSWD優化LMD方法與陶然等提出的SVD優化LMD方法分別對故障信號進行特征提取,并進行對比分析。

采用db10小波對下導靴故障振動信號進行1層WD。得到的高頻細節特征分量及其瞬時Teager能量波形如圖7所示。從圖7可以發現,D1(t)在1.5~2.0 s和3.5~4.0 s處,時域圖和瞬時Teager能量波形圖都表現出了明顯的沖擊響應特性,即含有大量故障特征信息。

圖5 3種狀態振動信號時域波形圖Fig.5 Three kinds of state vibration signal time domain waveform

圖6 3種狀態振動信號瞬時Teager能量波形Fig.6 Three kinds of state vibration signal instant Teager energy waveform

圖7 DW分量Fig.7 DW component

將D1(t)進行信號增強,得到其時域圖與瞬時Teager能量波形如圖8所示。將圖8與圖7進行對比可以發現,D1(t)經信號增強后其時域圖中的沖擊成分的最大峰值由約1.0增大到約2.0,瞬時Teager能量波形圖中的沖擊成分的最大峰值由約0.4增大到約1.0;因此可以得到結論,D1(t)經過增強,高頻細節特征信號相對質量得到很大提升,沖擊特性得到凸顯;后續進行信號重構,高頻細節特征信號在重構信號中的相對質量也將得到提升,沖擊特性也將得到凸顯。

圖8 DW分量增強Fig.8 DW component enhancing

將增強后的D1(t)替換原分量進行重構,再進行LMD,得到的PFi分量中,含有原始信號振動特征最多的分量是PF1,其次是PF2,PF3,…,PFn。因此,本文主要對PF1進行處理。PF1的時域圖與瞬時Teager能量波形,如圖9所示。

下導靴故障振動信號采用構造吸引子軌跡矩陣方法搭建二維矩陣,采用陶然等提出的SVD優化LMD方法對二維矩陣進行分析處理,得到的PF1分量的時域圖與瞬時Teager能量波形,如圖10所示。

將圖9與圖10進行對比,可以發現兩圖的時域圖和瞬時Teager能量波形圖在3~4 s和7~8 s處均表現出明顯的故障沖擊特性,且振幅較大;圖9的時域圖最大振幅約2.0、瞬時Teager能量波形圖最大振幅約4.0,圖10的時域圖最大振幅約1.0、瞬時Teager能量波形圖最大振幅約1.0。綜上對比分析,與SVD優化LMD方法相比,SSWD對原始信號中的高頻細節特征分量進行增強、重構,再采用LMD方法進行振動特征提取,使高頻細節特征分量在PF分量中的相對質量得到提升,沖擊特性得到凸顯。

與圖9的瞬時Teager能量波形圖相比,圖10在3~4 s處僅提取到了次級沖擊成分和主沖擊成分,一級衰減成分幾乎沒有提取到;在7~8 s處僅提取到了次級沖擊成分和一級衰減成分,主沖擊成分幾乎沒有提取到。綜上對比分析,SVD優化LMD方法僅提取到了部分故障沖擊成分,還有部分故障沖擊成分殘留在其它分量中,存在一個故障沖擊成分(振動模態)被分解到兩個或兩個以上分量中的情況,即模態混淆現象。SSWD優化LMD方法通過提取原始信號中的高頻細節特征信號進行增強、將增強后的高頻細節特征信號與低頻近似信號進行重構,再采用LMD方法進行振動特征提取,完整地提取了原始信號中的故障沖擊成分,避免了模態混淆現象的出現。

圖9 下導靴信號SSWD優化LMD結果Fig.9 The SSWD optimization LMD result of down guide boots signal

圖10 下導靴信號SVD優化LMD結果Fig.10 The SVD optimization LMD result of down guide boots signal

當電梯轎廂上導靴發生故障時,采用SSWD優化LMD方法得到的結果,如圖11所示,采用SVD優化LMD方法得到的結果,如圖12所示。

將圖11與圖12進行對比,可以發現兩圖的時域圖和瞬時Teager能量波形圖在2~3 s和6~7 s處均表現出明顯的故障沖擊特性,且圖11的振幅峰值均大于圖12;可以得到,SSWD優化LMD方法使高頻細節特征分量在PF分量中的相對質量得到提升,沖擊特性得到凸顯。

將圖11與圖12的瞬時Teager能量波形圖進行對比可以發現,在2~3 s處圖12僅提取到了主沖擊成分,次級沖擊成分而一級衰減成分幾乎沒有提取到;在6~7 s處圖12僅提取到了主沖擊成分和一級衰減成分,二級衰減成分幾乎沒有提取到。可以得到,SVD優化LMD方法僅提取到了部分故障沖擊成分,還有部分故障沖擊成分殘留在其它分量中,存在模態混淆現象。 SSWD優化LMD方法完整地提取了能夠表征原始信號振動特征的故障沖擊成分,避免了模態混淆現象的出現。

圖11 上導靴信號SSWD優化LMD結果Fig.11 The SSWD optimization LMD result of up guide boots signal

圖12 上導靴信號SVD優化LMD結果Fig.12 The SVD optimization LMD result of up guide boots signal

將圖9與圖11的瞬時Teager能量波形進行類比,可以發現圖9在3~4 s和7~8 s處、圖11在2~3 s處的故障沖擊信號波形(均含有次級沖擊成分、主沖擊成分、一級衰減成分)與db10小波函數從中間向兩邊衰減的“魚腹狀”波形非常相似。綜上分析,對振動特征類似“魚腹狀”波形的信號,采用db10小波對其高頻細節特征信號進行提取,是一種有效的方法。

5 結 論

(1)為凸顯信號中的高頻細節特征信號,常用的方法是抑制信號中的噪聲、干擾信號和平穩信號,但同時也會對信號中的高頻細節特征分量(故障特征分量)產生一定程度的抑制。本文提出的SSWD方法借鑒圖像增強方法,將信號分解為高頻細節特征信號和低頻近似信號;然后對高頻細節特征信號進行增強,將增強后的高頻細節特征信號與低頻近似信號進行重構,從而使故障特征分量得到凸顯,該方法可推廣應用到強噪聲背景下的故障信息提取。

(2)在獲取滾動導靴故障信息的過程中,SSWD只對目標信號進行了一層WD、LMD也只對目標信號進行一次分解就得到了滾動導靴故障信息,很大程度優化了算法的復雜度,減少了計算量。

(3)本文提出的基于SSWD優化LMD的高速電梯滾動導靴振動信號特征提取方法,完整地從高速電梯滾動導靴振動信號中提取了能夠表征原始信號振動特征的故障特征分量,避免了模態混淆現象出現,為高速電梯滾動導靴的故障診斷提供了參考。

猜你喜歡
特征提取電梯細節
以細節取勝 Cambridge Audio AXR100/ FOCAL ARIA 906
基于Gazebo仿真環境的ORB特征提取與比對的研究
電子制作(2019年15期)2019-08-27 01:12:00
留心細節處處美——《收集東·收集西》
被困電梯以后
一種基于LBP 特征提取和稀疏表示的肝病識別算法
細節取勝
Coco薇(2016年10期)2016-11-29 19:59:58
電梯不吃人
乘電梯
小說月刊(2015年4期)2015-04-18 13:55:18
基于MED和循環域解調的多故障特征提取
Walsh變換在滾動軸承早期故障特征提取中的應用
軸承(2010年2期)2010-07-28 02:26:12
主站蜘蛛池模板: 亚洲天堂日韩在线| 在线欧美日韩| 欧美69视频在线| 手机在线免费毛片| 91精品久久久无码中文字幕vr| 欧美综合一区二区三区| 久久综合丝袜长腿丝袜| 666精品国产精品亚洲| 国产乱视频网站| 欧美日韩中文国产va另类| 欧美精品黑人粗大| 亚洲精品自在线拍| 亚洲欧洲天堂色AV| 欧美中出一区二区| AV天堂资源福利在线观看| 日韩欧美网址| 国产成人亚洲精品无码电影| 国产浮力第一页永久地址| 国产麻豆aⅴ精品无码| 国产一在线观看| 国产人成在线观看| 国产又黄又硬又粗| 久久亚洲欧美综合| 国产91九色在线播放| 午夜激情婷婷| 亚洲精品无码久久毛片波多野吉| 2021精品国产自在现线看| 在线色综合| 91午夜福利在线观看精品| 激情在线网| 国产精品美女免费视频大全| 亚洲国产精品成人久久综合影院| 成人欧美日韩| 亚洲国产中文综合专区在| 亚洲第一综合天堂另类专| 久久黄色视频影| 国产精品流白浆在线观看| 国产一区二区三区夜色| 成年免费在线观看| 麻豆国产在线不卡一区二区| 精品少妇人妻av无码久久 | 四虎综合网| 国产精品福利在线观看无码卡| 2020最新国产精品视频| 欧美a级完整在线观看| 亚洲91精品视频| a级毛片网| 国产一区在线观看无码| 亚洲欧美在线精品一区二区| 亚洲国产天堂在线观看| 亚洲精品黄| 久久精品人妻中文系列| 麻豆精品国产自产在线| 国产传媒一区二区三区四区五区| 九九视频在线免费观看| 午夜爽爽视频| 国产高清无码麻豆精品| 无码一区18禁| 在线看片免费人成视久网下载 | 亚洲av综合网| 蜜桃臀无码内射一区二区三区| 中文字幕第4页| 玖玖精品视频在线观看| 在线欧美日韩国产| 蜜臀av性久久久久蜜臀aⅴ麻豆| 中文国产成人精品久久一| 91外围女在线观看| 欧美成人午夜视频免看| 91福利免费视频| 久久国产精品麻豆系列| 国产视频a| 久久婷婷六月| 亚洲色图综合在线| 成人字幕网视频在线观看| AV不卡国产在线观看| 国内精品自在自线视频香蕉| 亚洲成a人片| 亚洲永久精品ww47国产| 欧日韩在线不卡视频| 色婷婷电影网| 国产AV无码专区亚洲精品网站| 日韩欧美国产另类|