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

Hilbert-Huang變換中的模態(tài)混疊問題*

2016-11-23 11:07:50段玉波劉繼承
振動、測試與診斷 2016年3期
關(guān)鍵詞:模態(tài)信號分析

曹 瑩, 段玉波, 劉繼承

(東北石油大學(xué)電氣信息工程學(xué)院 大慶,163318)

?

Hilbert-Huang變換中的模態(tài)混疊問題*

曹 瑩, 段玉波, 劉繼承

(東北石油大學(xué)電氣信息工程學(xué)院 大慶,163318)

希爾伯特-黃變換(Hilbert-Huang transform,簡稱HHT)存在的模態(tài)混疊現(xiàn)象嚴(yán)重影響了實際應(yīng)用效果。在分析研究HHT原理及模態(tài)混疊產(chǎn)生機理的基礎(chǔ)上,提出了基于形態(tài)濾波預(yù)處理與端點延拓相結(jié)合的方法抑制模態(tài)混疊現(xiàn)象。與集合經(jīng)驗?zāi)B(tài)分解(ensemble empirical mode decomposition,簡稱EEMD)方法比較,所提出的方法能夠更快速、準(zhǔn)確地分解出表征信號的本征模態(tài)函數(shù)(intrinsic mode function,簡稱IMF)分量。將該方法應(yīng)用于滾動軸承的實測信號分析,結(jié)果表明,該方法在實際應(yīng)用中同樣具有很好的模態(tài)混疊抑制效果。

經(jīng)驗?zāi)B(tài)分解 ; 模態(tài)混疊; 形態(tài)濾波; 端點延拓

引 言

旋轉(zhuǎn)機械振動故障診斷主要是通過對機械設(shè)備的振動信號進(jìn)行一系列處理,進(jìn)而提取出能夠表征機械故障的特征信息,最終實現(xiàn)機械故障的診斷。在工程實際中,旋轉(zhuǎn)機械的振動信號大多為非線性、非平穩(wěn)的隨機信號,而傳統(tǒng)的Fourier變換無法滿足對此類信號的分析需求。1998年,Huang等人提出了HHT這種新型的時頻分析方法,該方法具有分析非平穩(wěn)、非線性信號及自適應(yīng)性的特點,在機械故障診斷領(lǐng)域得到了廣泛應(yīng)用[1-3]。

隨著HHT的不斷推廣和應(yīng)用,也逐漸暴露了在實際應(yīng)用中的問題。筆者主要針對HHT中的模態(tài)混疊問題進(jìn)行研究,通過對HHT原理及模態(tài)混疊現(xiàn)象產(chǎn)生機理的分析,針對性地提出了基于形態(tài)濾波預(yù)處理與端點延拓相結(jié)合的模態(tài)混疊抑制方法,并仿真驗證其可行性。通過對比分析其與集合經(jīng)驗?zāi)B(tài)分解方法在抑制模態(tài)混疊現(xiàn)象的效果差異,表明所提方法能夠更快速、準(zhǔn)確地分解出代表信號特征信息的IMF,并在對滾動軸承的實測振動信號分析中同樣得到了良好的模態(tài)混疊抑制效果。

1 HHT及模態(tài)混疊產(chǎn)生機理

1.1 HHT基本原理

HHT主要包括經(jīng)驗?zāi)B(tài)分解(empirical mode

decomposition,簡稱EMD)和Hilbert變換(Hilbert transform,簡稱HT)兩部分內(nèi)容。其中:EMD分解是將原始信號分解成若干個IMF分量;HT是對EMD分解所得的各項IMF分量進(jìn)行Hilbert變換,得到相應(yīng)的Hilbert時頻譜和邊際譜以進(jìn)行相應(yīng)的分析[1-2]。

對于任一給定信號,HHT過程如下。

1) EMD分解。首先,確定信號所有的局部極值點;其次,用三次樣條插值函數(shù)構(gòu)造上、下包絡(luò)線,并計算其均值曲線;最后,求取信號與該均值曲線的差,判斷其是否滿足IMF定義,并根據(jù)篩選停止條件反復(fù)判斷和篩選,直到最后剩余部分為一個單調(diào)信號,則分解完畢。

2) HT。EMD分解結(jié)束后,對各個IMF分量進(jìn)行Hilbert變換即可得到解析信號(不含殘余趨勢項)的Hilbert時頻譜及邊際譜。

1.2 模態(tài)混疊產(chǎn)生機理

模態(tài)混疊現(xiàn)象最早是由Huang等通過對含有間斷信號的EMD分解時發(fā)現(xiàn)的,該現(xiàn)象的產(chǎn)生與EMD分解中求包絡(luò)均值的篩分過程有關(guān)。由于間斷信號等一系列不連貫信號的存在導(dǎo)致局部極值點分布異常,為了保證信號包絡(luò)線的柔性和光滑性,包絡(luò)將會產(chǎn)生失真而出現(xiàn)模態(tài)混疊現(xiàn)象。模態(tài)混疊的出現(xiàn)不僅會導(dǎo)致嚴(yán)重錯假的時頻分布,也使IMF失去了物理意義,嚴(yán)重影響了EMD分解的準(zhǔn)確性。后續(xù)研究表明,除間斷信號外,脈沖干擾和噪聲等信號也可引起模態(tài)混疊現(xiàn)象[3]。

模態(tài)混疊具體有以下兩種表現(xiàn)形式:a.單個IMF中包含不同尺度或頻率的多個信號;b.同一尺度或頻率的信號被分解到多個不同的IMF中。

假設(shè)原信號是由頻率為50Hz的正弦信號和10dB的白噪聲信號疊加而成,則仿真得到圖1所示的原始信號(混有噪聲的正弦信號時域波形),圖2為對該信號進(jìn)行EMD分解后得到的IMF分量。

圖1 原始信號Fig.1 Original signal

圖2 EMD分解后的IMF分量Fig.2 IMF components after EMD

通過對比圖2中IMF2,IMF3分量的波形可知,在IMF2中出現(xiàn)了頻率較小、幅值較大的IMF3分量中的信號。由此說明,IMF2已經(jīng)發(fā)生了明顯的模態(tài)混疊現(xiàn)象。

2 模態(tài)混疊抑制方法

目前,對于抑制由間斷信號和噪聲干擾引起的模態(tài)混疊現(xiàn)象的方法很多,如中斷檢測、基于獨立分量分析、EEMD以及通過小波變換對信號進(jìn)行預(yù)處理等[4-6]。其中,EEMD主要利用白噪聲頻譜的均勻分布特性,通過對原始信號多次加入不同的白噪聲后進(jìn)行EMD分解,將多次分解結(jié)果進(jìn)行平均而得到最終的IMF分量,是一種簡便易行的分解方法[7-9]。對于旋轉(zhuǎn)機械的振動信號而言,產(chǎn)生模態(tài)混疊的主要原因即為環(huán)境噪聲對極值點的干擾,若采取合理的方法對信號進(jìn)行降噪處理,將會有效地抑制模態(tài)混疊現(xiàn)象。筆者研究了一種組合形態(tài)濾波預(yù)處理與特征尺度匹配延拓相結(jié)合的模態(tài)混疊抑制方法,從降噪和波形延拓兩方面來改善信號的極值點分布情況,以實現(xiàn)抑制模態(tài)混疊的效果。

2.1 形態(tài)濾波降噪

與傳統(tǒng)的降噪處理方法相比,數(shù)學(xué)形態(tài)學(xué)具有計算簡單、實用性好及時延較小等優(yōu)點。它主要包含腐蝕、膨脹、開運算和閉運算共4種基本運算[10-12]。

腐蝕和膨脹是最基本的運算。假設(shè)輸入序列f(n)為Df=(0,1,…,N-1)上的離散函數(shù),序列結(jié)構(gòu)元素g(n)為Dg=(0,1,…,M-1)上的離散函數(shù),且N≥M,則具體的計算方法如式(1)、式(2)所示

(1)

(2)

其中:m∈0,1,…,M-1。

開運算、閉運算是由膨脹和腐蝕兩種運算組合后得到的具有濾波性質(zhì)的運算方法,具體計算方法如式(3)、式(4)所示

(3)

(4)

基于上述4種運算,筆者通過開、閉運算級聯(lián)和組合平均構(gòu)造出一種三角形組合形態(tài)濾波器,這是一種平均組合形式的濾波器,能夠有效抑制信號中的各種噪聲成分。此濾波器輸出信號y(n)的表達(dá)式為

(5)

開-閉Foc、閉-開Fco的組合運算如下

(6)

2.2 端點延拓

目前,針對數(shù)據(jù)延拓的方法有很多,如鏡像延拓、波形特征匹配延拓、支持向量回歸機和神經(jīng)網(wǎng)絡(luò)延拓等。在具體應(yīng)用時,對端點處的數(shù)據(jù)不能盲目延拓,需要根據(jù)信號的特點,選擇合適的延拓方法以保證延拓后的波形要符合原始信號在端點處的變化趨勢[13-14]。

筆者采用的是一種基于特征尺度匹配的延拓方法,在信號內(nèi)在規(guī)律性較強的情況下,通過采用信號內(nèi)部和邊緣處變化趨勢最為相似的子波來對端點處數(shù)據(jù)進(jìn)行延拓;在信號內(nèi)在規(guī)律性較弱的情況下,只需考慮邊緣處的局部信息,根據(jù)邊緣局部極值點的特征,在信號邊緣兩側(cè)各添加一對極大值點和極小值點,對延拓的極值點序列進(jìn)行包絡(luò)擬合,估計出均值曲線,即得到完整的延拓波形。下面以信號右邊界點數(shù)據(jù)Sr的延拓為例,說明其具體延拓步驟。

1) 設(shè)Mr,Nr為Sr前的第1個極大值點和極小值點,則以Sr-Mr-Nr為邊界特征波形,在全部數(shù)據(jù)中找到與其構(gòu)成的三角形最接近的波形Sri-Mri-Nri(i=1,2,…)。

2) 計算波形Sri-Mri-Nri與邊界特征波形的匹配誤差Eri。此時,Eri存在以下兩種情況:a.若Eri在誤差允許范圍內(nèi),則取Eri最小的波形為匹配波形,從Sri的后一點數(shù)據(jù)開始,向后延拓波形數(shù)據(jù),使延拓數(shù)據(jù)符合信號的自然走向;b.若Eri不滿足允許范圍,即表明信號波形的內(nèi)在規(guī)律不明顯,難以找到與端點處變化趨勢最為相似的子波,此時,需要根據(jù)端點局部極值點的特征進(jìn)行延拓,即在信號右端添加一對極大值點、極小值點。其中,極值點的幅值根據(jù)與端點Sr臨近的3(或4)個極大、極小值點的平均幅值確定,對應(yīng)的添加位置根據(jù)與端點Sr最靠近的極大或是極小值點與臨近3(或4)個極值點的平均時間間隔確定。添加好極值點后,對新的極值點序列進(jìn)行包絡(luò)擬合,進(jìn)而估計出相應(yīng)的均值曲線。

該方法在具體實現(xiàn)過程中,僅需一次延拓即可完成,同時能夠使延拓后的數(shù)據(jù)與原始信號特征保持良好的一致性,很好地反映出信號的實際特征。

在HHT過程中,完成對原信號的端點延拓后,需進(jìn)行正常的EMD分解,對得到的IMF分量做Hilbert變換,然后按原始信號的長度及位置截取有效數(shù)據(jù)進(jìn)行相應(yīng)的結(jié)果分析。

2.3 仿真驗證

為了驗證本方法的可行性,對1.2節(jié)中的原信號進(jìn)行了仿真驗證。同時,為了對比說明本方法與傳統(tǒng)EEMD的優(yōu)劣,將兩種方法仿真得到的IMF分量進(jìn)行了對比,具體如圖3、圖4所示。

圖3 EEMD分解后的IMF分量Fig.3 IMF components after EEMD

圖4 本方法處理后的IMF分量Fig.4 IMF components after processing

在對比圖2~圖4的同時,為了進(jìn)一步分析兩種方法在抑制模態(tài)混疊方面的優(yōu)劣,筆者從運算時間、IMF分量個數(shù)及抑制效果上進(jìn)行了差異對比,見表1。

表1 模態(tài)混疊抑制效果差異對比

Tab.1 The differences contrast in the effects of mode-mixing restrain

措施前后方法運算時間/sIMF數(shù)量抑制效果措施前—1.026—措施后EEMD4.387很好形態(tài)濾波+端點延拓2.164很好

通過圖2~圖4,并結(jié)合表1數(shù)據(jù)的分析可知,經(jīng)過兩種方法處理后,原IMF2分量中的模態(tài)混疊現(xiàn)象均得到了很好的抑制。但由于EEMD涉及多次迭代運算,實時性較差,相較而言,本研究方法大大縮減了運算時間,在時效性方面更具優(yōu)勢。此外,根據(jù)IMF分量個數(shù)的對比可知,EEMD在分解過程中會產(chǎn)生許多無意義的IMF虛假分量,而經(jīng)本方法處理后的IMF分量個數(shù)明顯減少,為后續(xù)對IMF虛假分量的識別和剔除提供了便利,是較EEMD而言的另一優(yōu)勢所在。

3 實測信號分析

為了進(jìn)一步驗證筆者提出的方法的實操性,在仿真基礎(chǔ)上,對型號為HRB-N205EM的滾動軸承進(jìn)行了實測信號分析。實驗原始信號為滾動軸承的振動加速度信號,電機轉(zhuǎn)速為1 450 r/min,采樣頻率為5 kHz,采樣點數(shù)為5 000點。

實驗測得原始信號的時域波形如圖5所示,圖6為未采取模態(tài)混疊抑制方法所得到的IMF分量,圖7為采用筆者提出的模態(tài)混疊抑制方法得到的各IMF分量。

圖5 原始實測信號Fig.5 Time domain waveform of the original signal

圖6 未經(jīng)處理得到的IMF分量Fig.6 IMF components before processing

圖7 經(jīng)本方法處理后得到的IMF分量Fig.7 IMF components after processing

根據(jù)圖6、圖7的對比可知,筆者提出的方法在抑制模態(tài)混疊現(xiàn)象及減少無意義IMF分量上均有很好的效果。同時,在仿真結(jié)果的基礎(chǔ)上,綜合本次實測信號的分析結(jié)果可知,對于滾動軸承故障診斷而言,應(yīng)優(yōu)先考慮采用組合形態(tài)濾波方法與特征尺度匹配端點延拓方法來抑制模態(tài)混疊現(xiàn)象,既可以實現(xiàn)模態(tài)混疊現(xiàn)象的抑制,同時為后續(xù)的Hilbert變換及相應(yīng)的故障特征提取打下良好的基礎(chǔ),以實現(xiàn)對故障類型準(zhǔn)確、高效的判斷。

4 結(jié)束語

目前,HHT方法被廣泛用于旋轉(zhuǎn)機械故障的分析和診斷。該方法雖然對非平穩(wěn)、非線性信號的處理具有很大的優(yōu)勢,但仍存在端點效應(yīng)、模態(tài)混疊等問題。筆者通過對HHT原理及模態(tài)混疊現(xiàn)象產(chǎn)生機理的分析研究,有針對性地提出了基于形態(tài)濾波預(yù)處理與端點延拓相結(jié)合的模態(tài)混疊抑制方法。在仿真驗證該方法可行性的同時,從時間、IMF數(shù)量等方面對比分析了該方法與傳統(tǒng)EEMD的優(yōu)劣。分析表明所提出的方法能夠在抑制模態(tài)混疊的同時,更快速、準(zhǔn)確地分解出代表信號特征的IMF分量。同時,通過對滾動軸承實測信號的應(yīng)用分析,進(jìn)一步驗證了該方法在抑制模態(tài)混疊實際應(yīng)用中的良好效果。

筆者給出的方法在進(jìn)行端點延拓時,需要根據(jù)實際情況對匹配誤差Eri的限值Lm進(jìn)行調(diào)整。若對信號的規(guī)律性要求較低,可將限值設(shè)置較大;反之,需要將限值調(diào)小。但Lm不可以設(shè)置過大,否則延拓后的信號可能會與信號實際趨勢相差太遠(yuǎn);同樣,Lm也不可設(shè)置太小,否則會使對信號的規(guī)律性要求過于苛刻,導(dǎo)致不能找到合適的匹配波形。

[1] 雷亞國.基于改進(jìn)Hilbert-Huang變換的機械故障診斷[J].機械工程學(xué)報,2011,47(5):71-77.

Lei Yaguo.Machinery fault diagnosis based on improved hilbert-huang transform[J].Journal of Mechanical Engineering,2011,47(5):71-77. (in Chinese)

[2] 徐曉剛,徐冠雷,王孝通,等.經(jīng)驗?zāi)J椒纸?EMD)及其應(yīng)用[J].電子學(xué)報,2009,37(3):581-585.

Xu Xiaogang,Xu Guanlei,Wang Xiaotong,et al.Empirical mode decomposition and its application[J].Acta Electronica Sinica,2009,37(3):581-585.(in Chinese)

[3] 胡愛軍,孫敬敬,向玲.經(jīng)驗?zāi)B(tài)分解中的模態(tài)混疊問題[J].振動、測試與診斷,2011,31(4):429-434.

Hu Aijun,Sun Jingjing,Xiang Ling.Mode mixing in empirical mode decomposition [J].Journal of Vibration,Measurement & Diagnosis,2011,31(4):429-434.(in Chinese)

[4] Antonio H C,Stephan H.Adaptive time-frequency analysis based on autoregressive modeling[J].Signal Processing,2011,91(4):740-749.

[5] 王晶,陳果,郝騰飛.滾動軸承早期故障的多源多方法融合診斷技術(shù)[J].振動、測試與診斷,2013,33(5):868-874.

Wang Jing,Chen Guo,Hao Tengfei.Multiple sources and multiple methods about integration of diagnostic techniques based on ball bearing of vibration [J].Journal of Vibration,Measurement & Diagnosis,2013,33(5):868-874.(in Chinese)

[6] 湯寶平,董紹江,馬靖華.基于獨立分量分析的EMD模態(tài)混疊消除方法研究[J].儀器儀表學(xué)報,2012,33(7):1477-1482.

Tang Baoping,Dong Shaojiang,Ma Jinghua.Study on the method for eliminating mode mixing of empirical mode decomposition based on independent component analysis[J].Chinese Journal of Scientific Instrument,2012,33(7):1477-1482.(in Chinese)

[7] Zhang Jian,Yan Ruqiang,Gao R X,et al.Performance enhancement of ensemble empirical mode decomposition[J].Mechanical Systems and Signal Processing,2010,24(7):2104-2123.

[8] 胡愛軍,馬萬里,唐貴基.基于集成經(jīng)驗?zāi)B(tài)分解和峭度準(zhǔn)則的滾動軸承故障特征提取方法[J].中國電機工程學(xué)報,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].Proceeding of the CSEE,2012,32(11):106-111.(in Chinese)

[9] 陳仁祥,湯寶平,楊黎霞,等.自適應(yīng)參數(shù)優(yōu)化EEMD機械故障特征提取方法[J].振動、測試與診斷,2014,34(6):1065-1071.

Chen Renxiang,Tang Baoping,Yang Lixia,et al.An EEMD-feature extraction method of mechanical fault based on adaptive parameter optimum [J].Journal of Vibration,Measurement & Diagnosis,2014,34(6):1065-1071.(in Chinese)

[10]杜必強,唐貴基,石俊杰.旋轉(zhuǎn)機械振動信號形態(tài)濾波器的設(shè)計與分析[J].振動與沖擊,2009,28(9):79-81.

Du Biqiang,Tang Guiji,Shi Junjie.Design and analysis of morphological filter for vibration signals of a rotating machinery[J].Journal of Vibration and Shock,2009,28(9):79-81.(in Chinese)

[11]宋平崗,周軍,陳建亨.形態(tài)濾波優(yōu)化算法用于滾動軸承故障診斷[J].振動、測試與診斷,2013,33(5):756-762.

Song Pinggang,Zhou Jun,Chen Jianheng.Fault diagno-

sis method of rolling bearings based on optimized morphological filter algorithm [J].Journal of Vibration,Measurement & Diagnosis,2013,33(5):756-762.(in Chinese)

[12]胡振邦,張東升,章云,等.?dāng)?shù)學(xué)形態(tài)學(xué)濾波器在轉(zhuǎn)子失衡識別中的應(yīng)用[J].振動、測試與診斷,2014,34(6):1038-1044.

Hu Zhenbang,Zhang Dongsheng,Zhang Yun,et al.Research of rotor unbalance recognition based on mathematical morphology filter[J].Journal of Vibration,Measurement & Diagnosis,2014,34(6):1038-1044.(in Chinese)

[13]Wu Qin,Sherman D R.Boundary extension and stop criteria for empirical mode decomposition[J].Advances in Adaptive Data Analysis,2010,2(2):157-169.

[14]時培明,蔣金水,劉彬,等.基于邊界特征尺度匹配延拓的EMD改進(jìn)方法及應(yīng)用[J].中國機械工程,2014,25(12):1616-1623.

Shi Peiming,Jiang Jinshui,Liu Bin,et al.Improved method of EMD and its applications based on boundary characteristic scale matching extension method [J].China Mechanical Engineering,2014,25(12):1616-1623.(in Chinese)

10.16450/j.cnki.issn.1004-6801.2016.03.018

*黑龍江省長江學(xué)者后備支持計劃資助項目(2012CJHB005);黑龍江省教育廳科學(xué)技術(shù)研究資助項目(12531063)

2015-02-01;

2015-03-20

TH165.3; TP206.3

曹瑩,女,1987年2月生,博士生。主要研究方向為油氣田信息控制、信號處理及故障診斷、電力電子與電力拖動技術(shù)。曾發(fā)表《電壓自平衡式不對稱多電平逆電器對電機調(diào)速的研究》(《電工電能新技術(shù)》2011年第30卷第4期)等論文。

E-mail:cy1987@sina.cn

猜你喜歡
模態(tài)信號分析
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
隱蔽失效適航要求符合性驗證分析
完形填空二則
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
基于FPGA的多功能信號發(fā)生器的設(shè)計
電子制作(2018年11期)2018-08-04 03:25:42
電力系統(tǒng)及其自動化發(fā)展趨勢分析
基于LabVIEW的力加載信號采集與PID控制
國內(nèi)多模態(tài)教學(xué)研究回顧與展望
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
由單個模態(tài)構(gòu)造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 免费国产在线精品一区| 美女高潮全身流白浆福利区| 久久国产精品麻豆系列| 国产一区二区精品福利| 亚洲国产亚综合在线区| 亚洲一区二区三区中文字幕5566| 欧美精品1区| 少妇高潮惨叫久久久久久| 国产1区2区在线观看| 在线观看免费人成视频色快速| 亚洲精品国产首次亮相| 亚洲国产精品日韩专区AV| 国产在线观看人成激情视频| 国产成人AV大片大片在线播放 | 国产精品福利导航| 日韩无码黄色网站| 国产成人精品一区二区秒拍1o| 伊人色天堂| 天天爽免费视频| 伊人精品视频免费在线| 亚洲性视频网站| 亚洲天堂免费观看| 思思99思思久久最新精品| 99在线视频精品| 成人一级黄色毛片| v天堂中文在线| 男女男精品视频| 极品国产在线| 在线观看无码av五月花| 久久窝窝国产精品午夜看片| 青青极品在线| 深爱婷婷激情网| 国产欧美成人不卡视频| 国产网友愉拍精品| 丰满人妻被猛烈进入无码| 欧美性猛交xxxx乱大交极品| 伊人狠狠丁香婷婷综合色| 三级视频中文字幕| 天天综合网色| 亚洲第一极品精品无码| 97综合久久| 久久96热在精品国产高清| 亚洲欧洲综合| 无码专区国产精品一区| 国产欧美另类| 91精品国产情侣高潮露脸| 国产一区在线观看无码| 日韩国产亚洲一区二区在线观看| 中文字幕天无码久久精品视频免费| 中国一级特黄视频| 国产成人精品日本亚洲77美色| 国产日韩AV高潮在线| 99在线观看视频免费| 波多野结衣无码视频在线观看| 欧美精品影院| 欧美日韩在线成人| 日韩一级毛一欧美一国产| 亚洲AⅤ综合在线欧美一区| 国产精品久久国产精麻豆99网站| 黄色成年视频| 欧美不卡二区| 91精品国产一区| 日韩福利在线视频| 无套av在线| 色久综合在线| 婷婷综合色| 九九九精品成人免费视频7| 韩国福利一区| 一本一本大道香蕉久在线播放| 国产精品流白浆在线观看| 99久久精品无码专区免费| 99热这里只有精品5| 国产欧美精品一区aⅴ影院| 日韩欧美中文字幕一本| 午夜福利网址| 亚洲精品国产自在现线最新| 国产亚洲视频播放9000| 老司机午夜精品网站在线观看 | 浮力影院国产第一页| 中文精品久久久久国产网址| 国产欧美日韩在线一区| 91精选国产大片|