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

基于EEMD-TEO熵的高速列車軸承故障診斷方法

2018-04-11 08:48:07林建輝伍川輝黃晨光
西南交通大學學報 2018年2期
關鍵詞:模態振動故障

靳 行,林建輝,伍川輝,鄧 韜,黃晨光

(西南交通大學牽引動力國家重點實驗室,四川 成都 610031)

確保車輛的穩定安全運行對軸承開展早期故障診斷的研究具有重要意義.走行部作為車輛最重要的部件之一對支撐車體和列車走形起著非常重要的作用,同時對車行駛安全性及穩定性也有重要影響[1].軸箱軸承作為最主要的重載軸承,在聯系構架與輪對、保證輪對的回轉運動中發揮關鍵作用.車輛在長期運行中,軸承受到蛇形、橫向振動、軸向動態沖擊力等影響,軸承的外圈、滾動體和保持架等部件極易發生故障,影響列車的運行品質,嚴重者危及列車運行安全[2].

滾動軸承故障振動信息多表現為非線性、非平穩的調制信號[3],尤其是滾動運行中萌生的早期故障,特征信息微弱,受設備運行產生的強噪聲干擾,給故障診斷帶來困難[4-5].

作為信號處理的一種新型自適應、近似正交的分解方法,經驗模態分解同時適用于非線性、非平穩信號的分析.熵作為信號不穩定程度的一種新的度量也已經在信號處理和特征提取方面取得了廣泛應用.Pincus[6]將近似熵的改進算法應用在生理信號處理上,并得到廣泛關注.趙志宏等[7]提出了一種新的軸承故障診斷方法,該方法是在原樣本熵的基礎上基于EEMD(ensemble empirical mode decomposition)與樣本熵相結合的非線性數據維數約簡法,能有效地發掘非線性高維數據的結構本質,對高維數約簡與非線性數據的分析有顯著成效.用信息熵描述定量信號其不確定性與復雜程度的統計特性參數[8]具有抗噪能力強、穩定性好的特點.何政友等[9]在電力系統故障信號的特征提取中應用了信息熵.秦娜等[10]針對故障時高速列車轉向架振動信號特點,提出基于聚合經驗模態分解和5種信息熵相結合的特征提取方法.孫暉[11]利用離散Teager能量算子(Teager energy operator,TEO),改善了求取幅度包絡和瞬時頻率的精度.

本文將EEMD[12-13]、 TEO、信息熵相組合,對比HHT (Hilbert-Huang transform)與信息熵組合模式,提出一種高速列車軸承故障自適應診斷方法.使用EEMD處理輪對軸箱位置的振動信號,利用EEMD的自適應性將原信號分解成不同尺度的細節信號,對拆分的細節信號應用改進的TEO算法進行特征熵的提取,最終得到清晰的由故障引起的特征信息熵構成高維特征矢量,最終用于故障類型的判斷.通過高速列車軸承振動試驗數據對該方法進行檢驗,識別了高速列車的軸承狀態與故障類型.

1 EMD分解及等效濾波特性分析

1.1 振動信號EMD分解

黃顎做出突破性假設:任何信號都是由數個固有模態(intrinsic mode function,IMF)組成;對于所有的IMF分量,可以是線性,也可以為非線性的;一個信號可以包含多個IMF分量;每個IMF之間相互重疊,組成復信號,這種“篩分”的處理方法根據信號本身的特點自適應地將復雜信號的非平穩信號分解為若干個IMF分量,這種自適應多分辨率持性改進了小波分析中的恒定性,從而能有效地識別故障信號.信號x(t)由n個IMF分量Ci(t)和一個殘余項Rt(t)的和組成,其中i表示n個IMF分量的第i層,

(1)

式中:t={t1,t2,…,tk,…,tN},k∈N.

但是由于EMD在分解過程中,存在模態混疊問題,因此在EMD過程中,添加高斯白噪聲得到IMF分量Ci(ne,t),ne表示第ne次添加噪聲,聚合M次添加噪聲的IMF分量求均值,得到聚合分解后的IMF分量

(2)

EEMD可以有效地抑制模態混疊問題,其分解流程如圖1.

1.2 EEMD分解的等效濾波特性分級

EEMD可以有效地抑制模態混疊,但是由于添加了高斯白噪聲,第1層IMF分量為添加噪聲的高頻信號,該信號會導致原有信號的頻率成分變得異常復雜,同時也會混淆、淹沒信號本身的故障特征和信號的固有振動特性.這種混淆會對軸承的早期故障檢測、故障判定和故障溯源帶來不利影響.為了避免這種混淆發生,本文在基于EEMD分解理論的基礎上,將Teager能量算子與信息熵的判定結合到EEMD分解中,提出融合EEMD、TEO和熵的故障診斷方法.

圖1 EEMD分解流程Fig.1 EEMD decomposition flow chart

2 EEMD-TEO熵特征提取

根據軸承的振動信號特點,將EEMD分解、TEO算子及信息熵相結合,并對特征信號進行提取,可以得到關于信號的線性復雜度統計特性,即為EEMD-TEO算子熵.下面給出了TEO提取方法與EEMD分解后得到的4種EEMD-TEO熵的定義.

Teager能量算子是一個非線性算子,定義x(t)的非線性算子為

(3)

對該無衰減自由振蕩的線性振子,振動位移x(t)=Acos(ωt+θ),其中:A為振幅;θ為初相;角頻率ω={ω1,ω2,…,ωα},ωα為最大分析角頻率,則有

ψc(x(t))=(Aω)2.

(4)

該振子的瞬時總能量E是一個常數,即

(5)

式中:m為振子的質量.

比較式(4)、(5)的結果,二者只相差一個常數因子m/2.質量與速度平方的乘積為能量,因此,將ψc(·)算子稱為能量算子.將一階微分代入得到能量算子為

(6)

由式(4)、(6)得到瞬時包絡a(t)和瞬時頻率f(t)為

(7)

2.1 離散信號TEO算法

與連續Teager能量算子相對應,離散信號x(tk)的Teager能量算子定義為

ψd(x(tk))=x2(tk)-x(tk-1)x(tk+1).

(8)

根據差分方程的意義,連續到離散的映射關系表示如下:

(9)

式中:Δt=tk-tk-1.

由映射關系得到離散的時間能量分離算法(discrete-time energy separation algorithm,DESA).對式(1)中的Ci(t)逐一計算ψc(·),得到瞬時包絡幅值函數a(tk)及瞬時頻率函數f(tk),

(10)

式中:y(tk)=(x(tk)-x(tk-1))/Δt.

2.2 改進的離散TEO算法

為進一步減少誤差,對ψd(·)進行低通濾波,數字濾波器采用切比雪夫Ⅰ型,其單位脈沖響應為hi(n),則有

ψd1(·)=hi(n)*ψd(·).

(11)

用ψd1(·)代替ψd(·)代入式(10),便可得到改進DESA,即改進的離散TEO算法.

通常在求解信號的瞬時頻率時采用希爾伯特變化,如圖2所示,(a)為采用希爾伯特得到信號的瞬時頻率存在大量波紋及重疊的瞬時頻率,(b)為采用TEO算法求得的瞬時頻率波同樣存在大量波紋,(c)為采用改進后的TEO算法能準確地解調出原信號的瞬時頻率,并基本消除(a)、(b)中的波紋現象,使瞬時譜曲線更為平滑。

圖2 瞬時譜圖Fig.2 Improved TEO transient spectrum

2.3 EEMD-TEO譜分析

基于EEMD信號其局部特征時間尺度,自適應地將原始信號分解為若干個IMF分量之和的形式,瞬間頻率這一概念具有實際的物理意義,從而可以逐個計算IMF分量的瞬時頻率與包絡幅值.

(12)

這里省略了殘余項.展開公式稱為TEO希爾伯特譜,用角頻率表示瞬時頻率fi(t)=ωi(t)/2π,記作

(13)

2.4 EEMD能量熵

信號x(t)在進行M次經驗模態分解后,聚合得到i個IMF的均值Ci(t),計算第i層IMF的能量譜Ei和能量譜熵概率分布Pi為

(14)

(15)

聚合經驗模態能量譜熵EEEE為

(16)

2.5 EEMD奇異值熵

將EEMD分解結果Ci(t)組成n×N的矩陣Dn×N.應用奇異值對該矩陣進行分解,分解后的結果為d個非負奇異值λl(l=1,2,…,d),則有第l階增量EMD奇異熵為

(17)

于是得到聚合驗模態奇異熵EESE為

(18)

EESE可用來表示被分析信號的頻率組成成分與各頻率成分之間的分布特征情況.

2.6 Teager能量算子熵

將聚合經驗模態分解結果Ci(t)代入式(8)得到Teager能量算子ψd(x(tk)-x(tk-1)),則設第i層IMF的TEO概率為

(19)

Teager能量算子熵ETEOE為

(20)

2.7 EEMD-TEO時頻熵

EEMD得到的IMF分量經過TEO改進算法,得到的TEO希爾伯特譜H(ω,t)是二維矩陣,設能量在時域與頻域都是均勻分布,沿頻率變量ω和時間變量t分別得到概率矢量Pt(ω,t)和Pω(ω,t)為

(21)

(22)

Pt(ω,t)和Pω(ω,t)構成時間、頻率的近似熵,分別為EVTFEω(ω,t)、EVTFEk(ω,t)

(23)

(24)

將式(23)與式(24)合并求均值,得到EEMD-TEO時頻熵EETTFE.

(25)

3 軸承早期故障EMD-TEO熵檢測模型與試驗驗證

本文利用高速列車滾動試驗臺上的軸箱垂向振動信號作為樣本數據,對其信號進行EEMD分解,根據垂向振動信號被自適應地逐一分解成IMF,每個IMF均為一個近似穩態正交的簡單成分信號.對IMF矩陣計算4個EEMD-TEO熵值特征組成特征向量,將該向量作為支持向量機(support vector machine,SVM)的輸入.診斷結果由SVM給出.信號的處理過程如圖3所示.

圖3 信號的處理流程Fig.3 Signal processing flow

由圖3可知,軸承故障檢測模型有4個關鍵流程節點:第1,對原始振動信號EEMD分解得到的IMF分量;第2,改進的TEO算法;第3,進行4種EEMD-TEO熵的特征計算;第4,將4種EEMD-TEO熵作為SVM分類輸入,診斷軸承的狀態.

本文采用能量算子解調法,運用解調后的數據,根據故障特性,進行故障判斷.輪對滾動試驗臺試驗如圖4所示.

圖4 試驗現場Fig.4 Test site map

試驗臺的模擬列車運行速度為200 km/h.對輪對軸承做了人工傷處理,設置了滾動體故障與保持架故障,在兩個軸承的滾動體與保持架分別縱向燒蝕一條凹痕,深度1 mm,寬度1 mm.對正常軸承、保持架故障、滾動體故障3種工況分別進行試驗.圖5看不出因軸承故障產生的具有周期沖擊性特性的異常波形,無法判斷此軸承狀態.

圖5 軸承的振動試驗數據Fig.5 Bearing vibration test data

對圖5所示進行傅里葉變換,得到圖6所示頻域信號也很難區分正常軸承與故障軸承.

圖6 軸承的振動試驗數據Fig.6 Vibration data of bearing

參考文獻[14]所述雙列圓錐滾子軸承故障時,其故障頻率可以由轉速、軸承主要參數計算得到.根據表1參數,高速列車在200 km/h時單個滾動體故障特征頻率為68.08 Hz,保持架故障特征頻率8.86 Hz.

保持架故障頻率屬于低頻,極易被其他高頻信號淹沒,造成模態混疊,從信息熵特性來看,非線性信號越趨于線性,信息熵越小;線性信號中的低頻調幅波會出現的混疊失真導致故障信號熵變小.

表1 雙列圓錐滾子軸承的主要參數Tab.1 Main parameters of double row tapered roller bearings

將車速200 km/h時的3種工況的滾動試驗臺采集的軸承垂向振動信號進行截取,每1 s為1個樣本,采用簡單交叉驗證,通常預留1/4~1/3的樣本作為交叉驗證,剩下的作為訓練集使用[15],每種工況選用64個樣本作為訓練集,占總樣本的1/3,對所有300個樣本數據進行特征提取,在進行EEMD時,添加的白噪聲標準差為0.1,并將聚合次數M設置為1 000次.

盡管依據不同地熵定義,得到如表2所示不同軸承狀態下熵的均值,均值是存差異的.這種差異如圖7的盒圖表示尤為明顯,均值的大小與中值、上下四分位數分布一致,如本文所用的EEMD-TEO時頻熵與改進的TEO熵,均值位置都能反映不同工況下樣本數據的類內聚集位置.

表2不同工況下的各種熵均值
Tab.2Various entropy values under different operating conditionsbit

工況EEEEEESEEETOEEETTFEHHT無故障2.102.7514.953.666.28保持架故障2.012.7714.943.726.24滾動體故障2.352.8014.883.886.26

圖7 EEMD-TEO時頻熵與改進的TEO熵盒圖Fig.7 EEMD-TEO time-frequency entropy and Improved TEO entropy box diagram

由式(16)得到EEMD能量熵,在不同IMF分量尺度下,熵的大小如圖8所示.EEMD熵繼承了EEMD對信號的自適應性,當該尺度下IMF能量信息更明確時其能量熵也就越小;反之,存在故障時由于調制與模態混疊導致不同尺度下信號發生變化,熵值隨之變換.因此在不同工況中不同尺度的EEMD能量熵值均會發生改變,而不是簡單的增長.因此在不同工況下不同尺度的熵有明顯的區別,不是簡單的分離.由此IMF從信號范圍從高到低分解的頻率以及頻率較低的IMF往往表現出明顯的趨勢信號,因此其包含的信息也更為確定,熵更趨近于0.

圖8 EEMD能量譜熵的熵值分布Fig.8 EEMD energy spectrum entropy distribution

由式(18)得到EEMD奇異值熵特性是隨著信號的進一步分解而降低,但是不同軸承狀態的下降速度并不一致,如圖9所示,雖然第一層尺度上滾動體故障信息最低,但是其下降速度明顯不如其他兩種狀態快,對于正常軸承,其下降更為均勻.

圖9 EEMD奇異熵的熵值分布Fig.9 EEMD Singular entropy distribution

EEMD能量熵實質上是一種近似功率密度的熵,本文提出的TEO熵實質上是受時間尺度影響的功率熵.因此樣本時間越長熵值越大.如圖10所示,TEO算子熵每個尺度皆無明顯的變化趨勢,不同工況下,熵值大小有一定區別.

如圖11所示,EEMD-TEO熵是一種結合能量熵、TEO熵的近似熵,該熵包含能量熵的功率密度趨勢,又包含TEO熵中時間尺度的變化趨勢.因此在不同尺度下,不同工況的信息差異更大.

綜上所述,由于軸承狀態不同,振動信號熵的大小是不一致的.在噪聲干擾或強沖擊干擾情況下,僅僅依據輪對軸承振動信號的總熵值來辨別輪對軸承狀態是困難的[16-17].不同的熵有不同的意義,將不同意義的熵組成特征向量作為SVM的輸入,可以提高SVM的測試精度和泛化能力,對提高識別率有著較為積極的作用.

圖10 TEO算子熵分布Fig.10 TEO operator entropy distribution

圖11 EEMD-TEO熵分布Fig.11 EEMD-TEO time frequency entropy distribution

采用改進TEO算法對白噪音進行濾波后得到如圖12 不同種工況下樣本三維特征空間分布.

圖12 改進TEO時頻樣本熵向量Fig.12 Improved TEO time frequency sample entropy vector

可以看出在高維特征空間中不同工況的EEMD-TEO熵向量具已經有明顯的類間分離性與類內聚集性,沒有明顯的混疊現象.

為了驗證本文方法與以往方法的優越性,圖 13 所示給出由能量熵、奇異值熵、HHT時頻熵組成的傳統經驗模態熵向量的樣本分布.

圖13 HHT時頻熵向量Fig.13 HHT time frequency entropy vector

從圖12可以明顯看出,本方法具有更好的類間分離性.表3分別給出不同工況之間的各類熵向量方差,相對經驗模態熵向量,本文所提出的EEMD-TEO熵向量差異明顯大于經驗模態熵法.

表3 不同工況下熵值協方差Tab.3 Different conditions of entropy covariance

由表4可知,本文由于采用改進的TEO算子,避免了希爾伯特對稱性導致的偏差,層次更分明,在絕大多數情況下,改進TEO可提高希爾伯特變換得到的時頻熵的故障識別率.運行速度200 km/h時,EEMD-TEO熵提取法在對保持架故障識別率明顯優于傳統方法,對高速列車軸承振動信號的故障特征識別率方面也有很大提升.

表4兩種特征提取方法的識別率對比
Tab.4Comparison of recognition rates of two types of feature extraction methods%

特征提取方法正常軸承保持架故障滾動體故障傳統經驗模態熵95.3092.3098.46EEMD?TEO熵100.0098.4698.46

4 結 論

本文提出一種新型混合故障診斷模型,集成了聚合經驗模態、Teager能量算子、信息熵、SVM.

(1) 對EEMD分解后的IMF進行EEMD-TEO熵特征提取,將高維EEMD-TEO熵特征矢量作為SVM分類器的輸入,分類識別結果最終由SVM診斷得出,驗證了EEMD-TEO熵對高速列車早期軸承故障診斷的有效性.

(2) 對比EEMD-TEO熵與經驗模態熵,EEMD-TEO熵更有效的處理非線性信號,并有更好的類間分離性.

所用的樣本數據為輪對滾動平臺采集的人工傷的數據,而實際中軸承的故障常常是漸變的,因此研究當故障參數為漸變時,EEMD-TEO熵是否也可以自適應地有效識別工作狀態與故障類型將是下一步研究的重點和難點.

參考文獻:

[1]李熙.城市軌道交通車輛走行部安全評估方法研究[D].北京:北京交通大學,2011.

[2]唐德堯.廣義共振、共振解調故障診斷與安全工程:鐵路篇[M].北京:中國鐵道出版社,2006:11-13.

[3]胡愛軍,馬萬里,唐貴基.基于集成經驗模態分解和峭度準則的滾動軸承故障特征提取方法[J].中國電機工程學報,2012,32(11):106-111.

HU Aijun,MA Wanli,TANG Guiji.Rolling bearing fault feature extraction method based on ensemble empirical mode decomposition and kurtosis criterion[J].Proceedings of the CSEE,2012,32(11):106-111.

[4]陳志新,徐金梧,楊德斌.基于復小波塊閾值的降噪方法及其在機械故障診斷中的應用[J].機械工程學報,2007,43(6):200-204.

CHEN Zhixin,XU Jinwu,YANG Debin.Denoising method of block thresholding based on DT-CWT and its application in mechanical fault diagnosis[J].Chinese Journal of Mechanical Engineering,2007,43(6):200-204.

[5]李富才,何正嘉,陳進.小波域相關濾波法及其早期故障預示應用[J].振動工程學報,2005,18(2):145-148.

LI Fucai,HE Zhengjia,CHEN Jin.Wavelet transform domain correlation filter and its application in incipient fault prognosis[J].Journal of Vibration Engineering,2005,18(2):145-148.

[6]PUBCYS S M.Approximate entropy as a complexity measure[J].Chaos,1995,5(1):110-117.

[7]趙志宏,楊紹普.一種基于樣本熵的軸承故障診斷方法[J].振動與沖擊,2012,31(6):136-140

ZHAO Zhihong,YANG Shaopu.Sample entropy-based roller bearing fault diagnosis method[J].Journal of Vibration and Shock,2012,31(6):136-140.

[8]TANELLI M,PIRODDI L,SAVARESI S M.Real-time identification of tier-road friction conditions[J].IET Control Theory Applications,2009,3(7):891-906.

[9]HE Zhengyou,CHEN Xiaoqing,LUO Guoming.Wavelet entropy measure definition and its application for transmission line fault detection and identification;(part I:definition and methodology)[C]∥Power System Technology,2006.International Conference on Power System Technology.[S.l.]:IEEE,2006,82(10):1-6.

[10]秦娜,王開云,金煒東,等.高速列車轉向架故障的經驗模態熵特征分析[J].交通運輸工程學報,2014,14(1):57-64,74.

QIN Na,WANG Kaiyun,JIN Weidong,et al.Fault feature analysis of high-speed train bogie based on empirical mode decomposition entropy[J].Journal of Traffic and Transportation Engineering,2014,14(1):57-64,74.

[11]孫暉.經驗模態分解理論與應用研究[D].杭州:浙江大學,2005.

[12]WU Z H,HUANG N E.Ensemble empirical mode decomposition:a noise-assisted data analysis method[R].Calcerton:Center for Ocean-Land-Atmosphere Studies,2009.

[13]丁建明,王晗,林建輝,等.基于 EMD-Hankel-SVD 的高速列車萬向軸動不平衡檢測[J].振動與沖擊,2015,34(9):164-170.

DING Jianming,WANG Han,LIN Jianhui,et al.Detection of dynamic imbalance due to cardan shaft in high-speed train based on EMD-Hankel-SVD method[J].Journal of Vibration and Shock,2015,34(9):164-170.

[14]RANDALL R B,ANTONI J.Rolling element bearing diagnostics-a tutorial[J].Mechanical Systems and Signal Processing,2011,25(2):485-520.

[15]Linkin.斯坦福大學機器學習——交叉驗證(cross validation)[EB/OL].[2017-07-04].http://blog.csdn.net/linkin1005/article/details/42869331.

[16]YI C,LIN J,RUAN T,et al.Real time cardan shaft state estimation of high-speed train based on ensemble empirical mode decomposition[J].Shock and Vibration,2015(4):1-12.

[17]易彩.高速列車輪對軸承狀態表征與故障診斷方法研究[D].成都:西南交通大學,2015.

猜你喜歡
模態振動故障
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
故障一點通
中立型Emden-Fowler微分方程的振動性
奔馳R320車ABS、ESP故障燈異常點亮
國內多模態教學研究回顧與展望
故障一點通
基于HHT和Prony算法的電力系統低頻振蕩模態識別
江淮車故障3例
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 亚洲日韩AV无码一区二区三区人| 国产成人综合亚洲欧美在| 亚洲人成网18禁| 国产成人欧美| lhav亚洲精品| 日韩精品毛片人妻AV不卡| 免费一级毛片在线播放傲雪网 | 欧美国产精品不卡在线观看| 天堂中文在线资源| 亚洲AⅤ永久无码精品毛片| av无码一区二区三区在线| av无码久久精品| 国国产a国产片免费麻豆| 久久午夜夜伦鲁鲁片不卡| 一级黄色网站在线免费看| 免费国产高清视频| 高清视频一区| 四虎影视8848永久精品| 国产情侣一区| 无码高清专区| 久久99精品久久久久纯品| 99热这里只有精品在线播放| 国产91在线|中文| 欧美一级黄色影院| 欧美一区中文字幕| 国产成人精品一区二区三在线观看| 国产精品第一区在线观看| 久久公开视频| 亚洲妓女综合网995久久| 国产乱人乱偷精品视频a人人澡| 六月婷婷激情综合| 中文成人在线| 国产精品对白刺激| 国产亚洲成AⅤ人片在线观看| 狠狠ⅴ日韩v欧美v天堂| 国产一区二区精品高清在线观看| 日韩麻豆小视频| 亚洲中文字幕日产无码2021| 无码一区18禁| 99热线精品大全在线观看| 男人天堂亚洲天堂| 亚洲天堂精品在线| 国产精品永久在线| jijzzizz老师出水喷水喷出| 亚洲乱码在线播放| 久久视精品| 一级毛片免费播放视频| 日日噜噜夜夜狠狠视频| 国产精品第页| 成人中文在线| 欧美不卡视频在线| 高潮毛片免费观看| 日韩欧美国产三级| 国产99视频在线| 91精品综合| 真实国产乱子伦视频| 亚洲一区二区黄色| 男女精品视频| 91精品伊人久久大香线蕉| 欧美精品1区| 欧美a级在线| av大片在线无码免费| 国产白浆在线观看| 国产精品视频导航| 国产一区二区三区日韩精品| 91亚洲视频下载| 亚洲三级片在线看| 欧美啪啪一区| 青青青国产视频手机| 亚洲国产第一区二区香蕉| 久青草国产高清在线视频| 99这里只有精品在线| 亚洲无码视频一区二区三区| 久久精品这里只有国产中文精品 | 欧美专区在线观看| 一本大道香蕉高清久久| 国产成人高清精品免费| 91精品啪在线观看国产60岁 | 久草中文网| 国产在线观看第二页| 欧美日韩激情| 亚洲成a人片77777在线播放|