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

基于改進EMD的多繩摩擦提升機載荷信息特征提取

2014-06-07 05:55:04石瑞敏楊兆建
煤炭學報 2014年4期
關鍵詞:模態振動信號

石瑞敏,楊兆建

(太原理工大學機械工程學院,山西太原 030024)

基于改進EMD的多繩摩擦提升機載荷信息特征提取

石瑞敏,楊兆建

(太原理工大學機械工程學院,山西太原 030024)

為解決多繩摩擦提升機有效提升載荷的間接識別問題,提出了一種基于振動信號分析處理的特征提取方法。該方法以主軸裝置軸承處水平振動信號為對象,利用改進的經驗模式分解(EMD)方法將振動信號分解為若干有效固有模態函數(IMFs),改善了經典EMD方法模態混疊與端點效應現象,選取各階IMF的能量、方差貢獻率與能量矩作為特征值,探討了特征值與提升機有效載荷之間的內在聯系。結果表明,該方法中分析信號的獲取不改變提升機主軸系統結構,不影響提升機正常運行,易于實現長期在線監測,便于大量基礎數據的采集,為提升載荷定量識別積累了樣本數據。所選取的3個特征值從大小、權重及分布的角度較好地反映了提升載荷特征信息。

改進EMD;提升載荷;特征提取;振動分析

礦井提升設備作為煤礦生產中聯系井上井下的紐帶,是關系生產與人員安全的關鍵設備之一,一旦出現故障將造成巨大的經濟損失甚至人員的傷亡。在諸類安全事故中,由于過裝或未卸凈引起的過卷、滑動甚至墜罐,以及由于罐道劣變等引起的卡繩、松繩甚至斷繩等安全事故時有發生,而這些故障都可以直接或間接的由提升載荷反映出來。目前對提升載荷的動態檢測方法主要有壓輪法與直接測量法,前者通過將壓輪推向工作過程中的鋼絲繩,測得提升鋼絲繩對壓輪作用的壓力從而計算得到提升載荷,在調整鋼絲繩張力差方面得到了較好的應用效果[1]。后者將傳感器安裝于提升鋼絲繩與提升容器間[2]或安裝于提升容器器件間[3]直接測得提升載荷。兩類方法均需要對提升系統結構進行改造,對傳感器與監測系統安全可靠性要求較高。

振動信號是旋轉機械故障特征信息的載體[4],筆者在對提升機主軸兩側軸承進行監測的過程中發現,提升載荷作用于鋼絲繩,在提升過程中引發的振動將傳導至滾筒乃至主軸,使得軸承處振動信號包含著提升載荷的重要信息。與提升設備其他系統相比,提升載荷的特征信號比較微弱,信噪比低,具有很強的非線性、非平穩及噪聲特性,傳統的信號處理方法在處理這類信號方面存在缺陷。如短時傅里葉變換對不同時段變化了的信號只能加相同的窗,不適應信號頻率不斷變化的不同要求;小波分析克服了傅里葉變化窗口不可調的缺點,提高了分辨率,但對小波窗內的信號平穩性提出要求,因為沒有擺脫傅里葉變換的局限,并且存在信號能量泄露和非自適應性的缺點。經驗模式分解(empirical mode decomposition, EMD)是1998年由黃鍔等提出的可對非線性非平穩信號進行有效分析的信號處理方法[5]。由于其良好的局部特性與自適應性,EMD自提出以來被應用于地震信號[6-7]、語音識別[8]、故障診斷[9-11]等各個領域,并取得良好的效果。但是EMD方法仍舊存在一些不足,如模態混疊、端點效應、包絡線選擇及停止準則等[12-14]。本文旨在對EMD方法存在的問題進行改進,并在此基礎上提出一種多繩摩擦提升機載荷信息特征提取方法的新思路,為進一步實現載荷定量識別提供理論依據。

1 改進的EMD方法

1.1 EMD理論及存在問題

EMD方法是一種適用于非線性、非平穩性信號的時頻分析方法。此方法根據被分析信號本身特點構造基函數,自適應的選擇頻帶,從局部時間尺度入手,確定信號在不同頻段的分辨率,將復雜信號分解為具有不同尺度特征的固有模態函數(intrinsic mode functions,IMF),因此,該方法具有良好的局部適應性和時頻分辨率,較好地克服了小波變換存在的問題。EMD方法實質上可看作截止頻率與帶寬隨信號變化的帶通濾波器,從數據的特征時間尺度入手對信號進行分解,原始信號x(t)經EMD分解后,得到有限數目IMF和一個剩余分量,即

其中,n為分解后IMF數量;ci(t)為分解得到的n組IMF分量,各分量均滿足平穩性條件;r(t)為分解后的剩余分量。EMD方法因其在處理非線性非平穩性信號時表現出的良好特性,在實踐中有大量應用,但是理論本身仍有以下需要完善的不足之處。

(1)在EMD分解過程中存在模態混疊現象,即同一頻率信號被分解在不同的IMF中,或一個IMF中包含多個不同頻率的信號。當原始信號中存在間隙性成分或脈沖成分時,模態混疊現象更為明顯。

(2)EMD分解過程中采用極值點構造包絡曲線,而被分析信號端點處往往非極值點,這就導致每個分解過程中都會在端點處產生擬合誤差,隨著分解進行,誤差影響會逐步擴大。尤其是分解后期的低頻IMF分量,端點效應現象更為明顯,誤差的累積和傳播會導致分解結果的失真,產生虛假的IMF分量。

1.2 改進的EMD

針對以上提出的模態混疊、端點效應及虛假IMF分量的缺陷,提出采用集合經驗模式分解(ensemble empirical mode decomposition,EEMD)方法改善模態混疊現象,利用鏡像邊界延拓與灰色預測模型延拓方法相結合來抑制端點效應,最后采用相關系數法剔除假IMF分量。

(1)EEMD算法是Zhaohua Wu和Huang提出的一種噪聲輔助的數據分析方法[15],可以在不影響EMD自適應性的前提下減少模態混疊對EMD分解結果的影響。EEMD分解過程為:將高斯白噪聲與原始信號疊加后進行多次EMD分解,由于高斯白噪聲具有頻率均勻分布的統計特性,使得不同尺度的信號成分自動分布到合適的參考尺度上,從而降低各IMF分量的模態混疊程度;并且由于零均值噪聲的特性,經過足夠多次平均運算處理后所加入的高斯白噪聲將相互抵消,原始信號得到還原。

在應用EEMD方法時,有2個重要參數需要考慮:首先,對于所加入的高斯白噪聲的幅值選取十分重要,如果幅值過小,可能不足以引起信號局部極值點的變化,使得白噪聲的加入對EMD無作用;如果幅值過大,則會產生過多虛假IMF分量,使分解誤差增大,甚至淹沒原始信號特征使分解失真。其次,隨著分解次數的增加所加入白噪聲的幅值會隨之降低,理論上,分解次數越大則分解誤差越小,但是過大的分解次數代價是計算效率的降低,故應選取適當值。Huang[15]建議所加入白噪聲的幅值可以由原始信號的標準差乘以一個系數(一般取0.2)來確定,信號中高頻成分多時系數適當減小,反之則適當增大。分解次數可通過設置分解誤差來確定。

(2)抑制端點效應多采用在原始信號兩端通過一定方法增加信號長度,經HHT變換后再從結果中除去延長部分方法。目前所提出波形延拓、神經網絡延拓、支持向量回歸預測及邊界極值延拓等方法,在一定程度上消除了端點效應的影響,但仍存在局限性。本文提出利用灰色預測模型與鏡像延拓相結合的方法來抑制端點效應。根據被分析信號的特性,考慮到其發展過程與各影響因素是一種灰色不確定性關系,選取灰色動態預測模型GM(1,1)對信號進行延拓,此方法具有要求樣本數據少、原理簡單、計算量小、可檢驗等優點。具體延拓過程如下:

①以信號左端延拓為例,取靠近左端點處的5個極大(小)值組成原始序列,建立原始信號x(t)的灰色微分方程:

其中,a為發展灰數;b為內生控制灰數。二者由已知數據及最小二乘法確定,在此不贅述。求出a,b后可得預測模型:

②先用原始序列建立的GM(1,1)模型預測得到一個值,然后把這個值加入到原始序列中,同時剔除一個最遠端數據,以新的原始序列重復步驟①,再預測下一個數值,循環進行直至達到要求的預測精度。為滿足信號延拓要求,需要對數據極大值和極小值分別預測兩個值。

③對信號右端進行相同處理,最后將延拓得到的端點值與原信號連接起來形成新的信號序列。在此基礎上應用鏡像延拓法,以預測得到的各端點作為鏡像點進行延拓,從而得到較理想的延拓數據。

(3)EMD分解后得到的各IMF分量中,真實IMF分量與原始信號在理論上滿足正交性,具有較大的相關系數,而虛假IMF分量與原始信號相關系數較小。因此,可通過計算IMF分量與原始信號的相關系數來識別并剔除虛假IMF分量。

2 改進EMD方法的信號消噪分頻處理

運用改進的HHT方法對原始信號進行處理流程如圖1所示,其步驟包括:

圖1 改進的EMD處理流程Fig.1 Flowchart of improved EMD

(1)將原始信號x(t)加入一個正態分布的白噪聲序列nj(t),使得x~(t)=x(t)+nj(t)。

(2)將得到的信號x~(t)進行EMD分解,過程中在求取信號極值點后,采用鏡像與灰色預測模型相結合的方式對信號進行延拓,最終得到各IMF分量ci(t),i=1,2,3,…,n。

(3)重復步驟(1)2N次,每次加入新的隨機正態分布的白噪聲序列nj(t)。

(4)將N次分解得到的固有模態分量IMF做集成平均處理,當N足夠大時,可使添加白噪聲的影響趨近于0,進而得到IMF的最終結果。

(5)計算各IMF分量與原始信號的相關系數,以此為依據篩選有效IMF集。

(6)將有效IMF集重構信號,進行Hilbert變換。

為驗證此方法的有效性,將3個頻率分別為3, 11,16 Hz,幅值分別為1,0.4,0.9的正弦信號疊加,如下式:x(t)=sin(2π×3t)+0.4sin(2π×11t)+ 0.9sin(2π×16t),采樣頻率為500 Hz,采樣長度1 000點,時域信號如圖2所示。

圖2 仿真信號時域Fig.2 Time-domain chart of simulation signal

對數據分別采用原始與改進EMD兩種方法進行分析。圖3為原始EMD分解結果與信號的Hilbert時頻譜圖,從圖3(a),(b)可看出,IMF1由兩種頻率的正弦信號疊加,出現了明顯的模態混疊現象,反映到圖3(c)上,11與16 Hz譜線對應位置出現較大變形,未正確區分。另外從Hilbert時頻譜中可看出,信號兩端出現了明顯的頻率發散,并向信號內部傳播,靠近邊界處尤為嚴重,只有信號中間部分可以反映信號的特征。采用本文改進EMD方法對信號重新分解結果如圖4所示,從圖4(a),(b)EMD分解結果及頻譜可以看出,相近頻率11與16 Hz成功分解到兩個IMF分量中,模態混疊現象得以解決,從圖4(c)亦可看出,頻率發散現象趨勢變緩,端點效應得到抑制。

圖3 原始EMD分解結果Fig.3 Decompose results of original EMD

圖4 改進EMD分解結果Fig.4 Decompose results of improved EMD

3 基于IMF能量信息的提升載荷信息特征提取

反映提升機載荷信息的特征參數有電流、振動或直接測量鋼絲繩張力等。當提升機上提和下放載荷時,提升容器的振動會通過鋼絲繩作用于提升機從而反應在包括軸承與減速器的主軸系統上,由于載荷變化會引起軸承水平方向振動信號的波動,故本文拾取主軸裝置兩側軸承水平方向振動信號,采用改進的EMD方法對振動信號中包含的載荷特征信息進行提取。經改進EMD處理后得到的有效IMF集中,各階IMF分量所含的頻率成分不同,階數越低所含信號的高頻成分則越多。各階IMF分量的能量以及在整個頻帶上的能量分布能夠刻畫出信號的整體特征,因此,本文提出以IMF分量的能量變化作為參考依據提取載荷信息特征。

圖5為某礦副立井型號為JKM 2.8×6(I)A的多繩摩擦提升機主軸系統,傳感器布置在軸承座頂端,采集每個提升過程中水平方向(以圖中X方向為正向)振動信號,采樣頻率500 Hz,截取勻速提升時6 s數據作為待分析信號。副立井主要負責運送人員、設備以及提升矸石等,故載荷情況較為復雜,每次提升載荷情況差別較大,本文以1號罐提升工況選取5組數據作為示例,載荷數據見表1。

圖5 多繩摩擦提升機主軸裝置傳感器布置Fig.5 Sensors placement on spindle device of multi-rope hoist

表1 示例載荷數據Table 1 Load data of sampleskN

首先采用改進的EMD方法對被分析信號進行消噪分頻處理,得到5階有效IMF集,如圖6所示。

圖6 實測信號的改進EMD分解Fig.6 Improved EMD of actual signal

針對被測振動信號的非平穩性特點,當載荷發生變化時,其振動信號在各個頻段的能量必將發生變化,而隨著工況不同,經改進的EMD分解后的各有效IMF分量包含了不同頻段的信息,為了更好地提取隱藏在信號中的載荷信息,選取IMF分量能量值,方差貢獻率及能量矩從能量大小、權重和分布的角度對信號本質特征進行分析。其定義分別為

圖7為各階IMF分量特征值,各特征值隨載荷的變化。從圖7(a)可以看到,隨著提升載荷總重的增加,第1階IMF能量值呈上升趨勢,但當提升載荷總重相近時,載荷總重大但載荷差小的工況第1階IMF能量值反而小,需要從第2階及以后的IMF能量上反映真實載荷情況;圖7(b)為各階IMF分量的方差貢獻率,體現了各階能量在總能量中的權重,可以看到能量主要集中于前兩階IMF分量中;還選取了能量矩作為表征載荷信息的特征值,如圖7(c)所示,能量矩指標考慮到了IMF能量隨時間t的分布情況,可以體現能量的分布特征。從上述分析可以得到,相對于不能完整快捷反映載荷信息的單一的能量指標,將可以體現能量的大小、權重及分布的指標結合起來,可以更好地揭示IMF集能量特征以及與載荷的對應關系,有利于后續的載荷識別與故障診斷。

圖7 各階IMF分量特征值Fig.7 Characteristic parameters of IMFs

4 結 論

(1)多繩摩擦提升機主軸裝置軸承處水平振動信號中蘊含了提升載荷的一些特征信息。

(2)提出了一種改進的經驗模態分解方法,仿真結果表明此方法改善了原始EMD方法存在的模態混疊、端點效應以及虛假IMF分量問題,快速準確的得到了有效IMF集,可用于其他類型信號的消噪分頻處理。

(3)利用固有模態函數的能量值、方差貢獻率與能量矩從能量大小、權重及分布3方面相結合作為提升載荷的主要特征,可以更好地反映提升載荷與振動信號的對應關系。

(4)特征提取是提升載荷識別和診斷中的重要一環,本文提出的方法可以在不改變提升主軸系統結構,不影響提升機正常運行的前提下有效提取多繩摩擦提升機提升載荷特征,可作為進一步建立提升載荷識別模型的依據,為提升載荷定量識別與監測提供了新思路。

[1] 楊兆建,賈長喜,王勤賢,等.多繩提升機鋼絲繩張力動態測試調整研究[J].煤炭學報,1995,20(5):519-524.

Yang Zhaojian,Jia Changxi,Wang Qinxian,et al.Dynamic test and adjustment of wire-rope tension of multiple rope friction winder[J].Journal of China Coal Society,1995,20(5):519-524.

[2] 騰孝來,肖興明,胡 明,等.多繩摩擦提升鋼絲繩張力動態無線監測[J].煤礦機械,2008,29(5):207-208.

Teng Xiaolai,Xiao Xingming,Hu Ming,et al.Dynamic wireless tension monitoring of wire-rope of multiple friction winder[J].Coal Mine Machinery,2008,29(5):207-208.

[3] 邵海燕.多繩提升機載荷動態監測與提升安全研究[D].青島:山東科技大學,2003.

Shao Haiyan.Study of the real-time monitoring and security of multiwire hoisting load[D].Qingdao:Shandong University of Science and Technology,2003.

[4] 鞠萍華,秦樹人,秦 毅,等.多分辨EMD方法與頻域平均在齒輪早期故障診斷中的研究[J].振動與沖擊,2009,28(5):97-101.

Ju Pinghua,Qin Shuren,Qin Yi,et al.Research on earlier fault diagnosis of gear by method of multi-resolution empirical mode decomposition and frequency domain averaging[J].Journal of Vibration and Shock,2009,28(5):97-101.

[5] Huang N E,Shen Z,Long S R,et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis[J].Proceedings of the Royal Society of London, 1998,454(A):903-995.

[6] Huang N E,Chern C C,Huang K,et al.A new spectral representation of earthquake data:Hilbert spectral analysis of Station TCU129, Chi-Chi,Taiwan,21 September 1999[J].Bulletin of the Seismological Society of America,2001,91(5):1310-1338.

[7] 石春香,羅奇峰.時程信號的Hilbert-Huang變換與小波分析[J].地震學報,2003,25(4):398-405.

Shi Chunxiang,Luo Qifeng.Hilbert-Huang transform and wavelet analysis of time history signal[J].Acta Seismologica Sinica,2003, 25(4):398-405.

[8] Aicha Bouzid,Noureddine Ellouze.Empirical mode decomposition of voiced speech signal[A].First International Symposium on Control, Communications and Signal Processing[C].2004:603-606.

[9] 程軍圣,于德介,楊 宇,等.EMD方法在轉子局部碰摩故障診斷中的應用[J].振動、測試與診斷,2006,26(1):24-27.

Cheng Junsheng,Yu Dejie,Yang Yu,et al.Application of EMD to local rub-impact fault diagnosis in rotor systems[J].Journal of Vibration Measurement&Diagnosis,2006,26(1):24-27.

[10] Bassiunya A M,Lib Xiaoli.Flute breakage detection during end milling using Hilbert-Huang transform and smoothed nonlinear energy operator[J].International Journal of Machine Tools and Manufac-ture,2007,47(6):1011-1020.

[11] 時培明,丁雪娟,李 庚,等.一種EMD改進方法及其在旋轉機械故障診斷中的應用[J].振動與沖擊,2013,32(4):185-190.

Shi Peiming,Ding Xuejuan,Li Geng,et al.An improved method of EMD and its application in rotating machinery fault diagnosis[J].Journal of Vibration and Shock,2013,32(4):185-190.

[12] 李方熙,陳桂明,劉希亮,等.希爾伯特-黃變換端點效應的自適應端點相位正弦延拓方法[J].上海交通大學學報,2013, 47(4):594-601.

Li Fangxi,Chen Guiming,Liu Xiliang,et al.Processing method for Hilbert-Huang transform end effects self-adaptive endpoint-phase sinusoidal extension[J].Journal of Shanghai Jiaotong University, 2013,47(4):594-601.

[13] 鐘佑明,金 濤,秦樹人.希爾伯特一黃變換中的一種新包絡線算法[J].數據采集與處理,2005,20(1):13-18.

Zhong Youming,Jin Tao,Qin Shuren.New envelope algorithm for Hilbert-Huang transform[J].Journal of Data Acquisition&Processing,2005,20(1):13-18.

[14] 王增才,王樹梁,任鍇勝,等.基于EEMD的提升機天輪軸承故障診斷方法[J].煤炭學報,2012,37(4):659-694.

Wang Zengcai,Wang Shuliang,Ren Kaisheng,et al.Research on the method of hoist head sheave bearing fault diagnosis based on EEMD[J].Journal of China Coal Society,2012,37(4):659-694.

[15] Wu Zhaohua,Norden Eh Huang.Ensemble empirical mode decomposition:A noise assisted data analysis method[J].Advances in Adaptive Data Analysis,2008,1(1):1-41.

Feature extraction for hoisting load of multiple rope friction hoist based on improved EMD

SHI Rui-min,YANG Zhao-jian

(School of Mechanical Engineering,Taiyuan University of Technology,Taiyuan 030024,China)

In order to indirectly identity the payload of multi-rope friction hoist,a feature extraction method was proposed based on analysis of horizontal vibration signal which was collected from bearing caps of spindle device.Firstly, improved empirical mode decomposition(EMD)that can restrain the mode mixing phenomenon and reduce the boundary effect was used to decompose the signal into several intrinsic mode functions(IMFs),then the energy,contribution rate of variance and energy moment of IMFs were calculated as characteristic parameters.The intrinsic relationship between the characteristic parameters and the payload of hoist were discussed.The method of signal acquisition didn’t change the structure of spindle system and had no effect on normal operation of hoist.A long-term online monitoring was apt to accomplish,the basic data was easy to obtain,furthermore,a large amounts of sample datum were accumulated for quantitative identification of actual load.It is shown that the characteristic parameters proposed could sensitively reflect the load characters in terms of magnitude,weight and distribution.

improved EMD;hoisting load;feature extraction;vibration analysis

TD53

A

0253-9993(2014)04-0782-07

石瑞敏,楊兆建.基于改進EMD的多繩摩擦提升機載荷信息特征提取[J].煤炭學報,2014,39(4):782-788.

10.13225/j.cnki.jccs.2013.1502

Shi Ruimin,Yang Zhaojian.Feature extraction for hoisting load of multiple rope friction hoist based on improved EMD[J].Journal of China Coal Society,2014,39(4):782-788.doi:10.13225/j.cnki.jccs.2013.1502

2013-10-18 責任編輯:許書閣

山西省科技重大專項資助項目(20111101040)

石瑞敏(1983—),女,山西文水人,博士研究生。Tel:0351-6010414,E-mail:srm0018@link.tyut.edu.cn。通訊作者:楊兆建(1955—),男,河北高陽人,教授,博士生導師。Tel:0351-6010414,E-mail:yangzhaojian@tyut.edu.cn

猜你喜歡
模態振動信號
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
中立型Emden-Fowler微分方程的振動性
基于LabVIEW的力加載信號采集與PID控制
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: www.精品国产| 国产精品天干天干在线观看| 97视频免费看| 五月婷婷伊人网| 91www在线观看| 高清欧美性猛交XXXX黑人猛交 | 国产理论精品| 免费一极毛片| 日韩精品亚洲一区中文字幕| 久久夜色精品国产嚕嚕亚洲av| 色偷偷综合网| 久久久久人妻一区精品| 亚洲欧洲日产国产无码AV| 国产精品第5页| 国产成人超碰无码| 依依成人精品无v国产| 性喷潮久久久久久久久| 日韩福利在线视频| 在线观看免费国产| 26uuu国产精品视频| 国产欧美日韩综合在线第一| 国产成人免费视频精品一区二区| 丁香婷婷激情网| 久久公开视频| 国产剧情无码视频在线观看| 免费毛片a| 曰AV在线无码| 国产精品不卡片视频免费观看| 成人综合网址| 亚洲最新网址| 国产大全韩国亚洲一区二区三区| 亚洲二区视频| 色成人亚洲| 日本国产精品一区久久久| 亚洲美女一级毛片| 中文字幕1区2区| 国内精品一区二区在线观看| 色有码无码视频| www.狠狠| 色综合婷婷| 日韩精品无码免费一区二区三区| 国产网站免费| 国产成人综合网| 综合五月天网| 黄色福利在线| 欧美成人手机在线视频| 免费看一级毛片波多结衣| 91色综合综合热五月激情| 免费观看无遮挡www的小视频| 高潮毛片免费观看| 国产靠逼视频| 宅男噜噜噜66国产在线观看| 国产一二三区在线| 欧美一区精品| 国产精品lululu在线观看| 91亚洲影院| 国产91九色在线播放| 三上悠亚一区二区| 欧美另类视频一区二区三区| 亚洲AV无码一区二区三区牲色| hezyo加勒比一区二区三区| 精品国产香蕉在线播出| 久久9966精品国产免费| 色首页AV在线| 又大又硬又爽免费视频| 日韩AV无码免费一二三区| 国产精品理论片| 亚洲第一区精品日韩在线播放| 最新精品久久精品| 伊人成人在线| 成人午夜久久| 免费看a毛片| 亚洲一区二区三区香蕉| 秘书高跟黑色丝袜国产91在线| 中文国产成人精品久久一| 亚洲va在线∨a天堂va欧美va| 香蕉eeww99国产精选播放| 色老二精品视频在线观看| 国产特级毛片| 欧美无遮挡国产欧美另类| 国产精品福利社| 天天综合网站|