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

改進EEMD方法及混沌降噪應用研究

2017-09-25 06:02:07位秀雷林瑞霖劉樹勇楊慶超
振動與沖擊 2017年17期
關鍵詞:模態細節振動

位秀雷, 林瑞霖, 劉樹勇, 楊慶超

(1. 中國人民解放軍91404部隊,河北 秦皇島 066000; 2. 海軍工程大學 動力工程學院,武漢 430033)

改進EEMD方法及混沌降噪應用研究

位秀雷1, 林瑞霖2, 劉樹勇2, 楊慶超2

(1. 中國人民解放軍91404部隊,河北 秦皇島 066000; 2. 海軍工程大學 動力工程學院,武漢 430033)

在總體平均經驗模態分解(Ensemble Empirical Mode Decomposition,EEMD)降噪過程中,對本征模態分量(Intrinsic Mode Function,IMF)的有效處理一直是影響降噪效果的關鍵。為此,提出一種基于改進EEMD的去噪方法。基于“3σ”法則和奇異值分解(Singular Value Decomposition,SVD)提取第一個IMF分量中有用信號細節。利用連續均方誤差準則對剩余IMF分量進行高低頻區分,分別使用SVD和S-G算法提取高低頻分量的有用信號,可以有效避免了高頻部分有用信號的流失,同時剔除低頻分量中的部分噪聲,克服了EEMD去噪時IMFs難以有效處理的不足。為了驗證該方法的有效性,進行了數字仿真與雙勢阱混沌振動試驗,結果表明,該方法的降噪效果優于小波加權和EEMD去噪方法。

總體平均經驗模態分解;混沌信號;奇異值分解;降噪;S-G濾波

混沌振動信號在采集過程中不可避免受到復雜噪聲干擾,極大干擾了信號真實信息的解讀,從而影響了混沌振動信號控制、識別、預測等進一步分析,并且混沌復雜的動力學行為使其具有功率譜寬帶性和似噪聲性,其頻帶與疊加噪聲的頻帶往往全部或部分重疊,傳統的線性去噪方法很難實現有效濾波[1]。

Huang等[2]提出了一種新的信號處理方法-EMD(Empirical Mode Decompomposition)方法。與小波變換方法相比,EMD無需信號的先驗知識,其分解完全依賴信號本身,數據分解真實可靠,因此被廣泛應用于機械振動信號分析[3-6]。但是,在脈沖強干擾的影響下,EMD分解出來的本征模態分量(Intrinsic Mode Function,IMF)會發生畸變,導致信號失真[7],且EMD本身存在一些不足,如模式混疊、端點效應、停止條件等[8]。為了抑制模式混疊,Wu等[9]提出了總體平均經驗模態分解方法,有效地克服了這一缺陷。在IMF處理過程中,文獻[10]將EEMD 方法用于對疲勞應變信號降噪,其方法是對信號做EEMD 分解后,選取IMF 分量來重構信號以對疲勞信號降噪,但文中也未說明IMF分量的選取方法;陳仁祥等基于能量密度與其平均周期的乘積為一常量這一特點,設計了自動選擇IMF 分量重構信號的算法。但是,EEMD的時空域算法因簡單地去掉1個或多個IMF分量實現去噪,會導致相應分量上的有效信號一起被剔除,進而導致信號失真,為了在分離噪聲的過程中盡可能地保留有用信號,本文提出了改進EEMD降噪方法。首先對信號進行EEMD分解,分析各IMF分量特性,然后分別使用“3σ”法則、奇異值分解[11]和S-G濾波方法對IMF分量進行處理,克服了EEMD去噪時IMFS選擇的難題。為了驗證本文方法的有效性,進行了數字仿真與雙勢阱混沌振動試驗,結果表明,本文方法優于小波加權和傳統的EEMD去噪方法。

1 EEMD原理

利用EMD進行信號處理時,由于異常事件干擾導致極值點分布不均勻產生了模式混疊現象。為此,Wu等將白噪聲加入待分解信號來抑制異常事件,利用白噪聲頻譜的均勻分布來使不同尺度的信號自動分布到合適的參考尺度上。同時,利用白噪聲的零均值特性,經過多次平均使噪聲相互抵消,從而抑制甚至完全消除噪聲的影響。EEMD的本質就是疊加高斯白噪聲的多次經驗模式分解,其步驟如下:

(1) 在原始信號x(t)中疊加均值為0,幅值和標準差為常數的高斯白噪聲ni(t),i=1~M,疊加次數為M(M>1),即:

xi(t)=x(t)+ni(t)

(1)

(2) 對xi(t)進行EMD分解,得到N個IMF記為aij(t),j=1~N,余項表示為ri(t)。其中aij(t)表示第i次疊加高斯白噪聲后,分解得到的第j個IMF分量。

(3) 由于不相關隨機序列的統計均值為0,所以將以上步驟所得的IMF進行平均運算,即可消除多次疊加高斯白噪聲對真實IMF的影響,平均后得到的IMF為

(2)

式中,aj(t)表示對原始信號進行EEMD分解后所得的第j個IMF分量。

2 閾值決策方法

2.1利用“3σ”法則和SVD提取IMF1中的有用信號

信號經EEMD去噪處理,通常認為第一層IMF分量全部由噪聲構成,但隨著學者深入研究發現,IMF1中仍含有一定量的信號成分[12]。提取IMF1中的信號細節成分,會提高去噪效果,減少信號失真。但由于缺乏先驗知識,因此,對IMF1處理仍是個難題。在IMF1中,噪聲占絕大部分,僅有少量的信號細節成分,且所含噪聲仍近似服從零均值正態分布,文獻[12]利用“3σ”法則進行細節信息提取。

(3)

由式(3)可以看出,“3σ”法則提取的細節分量是一種窄時間段的高幅值沖擊信號,圖1是利用“3σ”法則對信噪比為5 dB的Lorenz時間序列經EEMD分解后IMF1提取的細節信息,顯然,直接把其作為IMF1提取的有用信號是不合理的,并且經大量仿真實驗證明:信號含噪聲強度越大,“3σ”法則提取的細節分量就越多;反之,越少,甚至當噪聲強度小至系數幅值都小于3σ,就提取不到任何細節分量,圖2是利用“3σ”法則對信噪比為15 dB的Lorenz時間序列經EEMD分解后IMF1提取的細節信息。因此需分兩種情況對進行處理,第一是“3σ”法則提取的細節分量不為零,即對象是噪聲強度大的混沌信號;第二是是“3σ”法則提取的細節分量為零,即對象噪聲強度小的混沌信號。

圖1 SNR為5 dB的IMF1細節

圖2 SNR為15 dB的IMF1細節

對于第二種情況,噪聲強度較小,利用傳統的非線性去噪方法即可取得較好的降噪效果,本文不再詳述。第一種情況,即低信噪比、高強度的混沌含噪信號,先利用“3σ”法則提取IMF1的細節分量,再采用奇異值分解(SVD)方法對有用信號進行提取,SVD降噪的關鍵在于選擇合適的奇異值進行重構,由于“3σ”法則提取IMF1的細節分量可以作為IMF1所含有用信號估計的先驗知識,因此可以利用有用信號能量相近原則選取SVD重構的奇異值個數。具體的步驟如下:

1) 利用式(3)對進行IMF1細節提取;

2) 利用下式估計IMF1所含有用信號的能量;

(4)

式中:M為滿足“3σ”法則條件的本征模態分量個數。

3) 對IMF1進行奇異值分解,假設分解后奇異值為λ=[λ1,λ2,…,λn],若選擇前k個特征值重構得到降噪后的信號S=[S1,S2,…,Sk],則根據下式便可求出k的值:

(5)

2.2SVD提取高頻分量中的有用信號

由于含噪信號經EEMD分解,IMF1相當于一個高通濾波器,其頻譜占據了原始信號所含頻譜的一半,白噪聲的各階IMF的能量密度按照2倍關系逐漸遞減,即隨著IMF階數的增加白噪聲的比重大幅度減少,而有用信號的比重大幅度增加,因此對于高頻IMF,IMF1經奇異值分解,各個特征量比較平均,而其他IMF經奇異值分解后,信號的特征值會明顯大于噪聲的特征值,依據這一特性,便可利用差分譜的方法選擇合理的奇異值個數進行重構,達到去噪目的。

首先,需要對IMF分量進行高低頻區分,噪聲主要分布在高頻IMF分量上,而信號主要分布在低頻IMF分量上,因此可以利用連續均方誤差準則[13](Consecutive Mean Square Error,CMSE)對信號分量起主導作用模態與噪聲起主導作用模態進行區分,即找到一個索引值js,使得從該索引開始往后的IMF分量對信號進行重構的誤差最小。

(6)

式中:N為信號長度;n為IMF分量的個數;IMFk(ti)表示第k個IMF的第ti個分量的重構誤差,基于該準則,索引值js可由下式給出:

(7)

式中,arg min表示重構誤差取最小的函數。

2.3S-G濾波方法平滑低頻分量

由于低頻IMF分量中大部分都是有用信號,噪聲成分比較少,所以奇異值分解后特征值分布比較平均,不適合利用SVD提取有用信號。趙志宏等指出Savitaky-Golay濾波算法對低頻IMF分量有很好的降噪效果,因此,本文低頻IMF分量按式(8)做SG平滑處理。設ci是其中的一個IMF分量系數,在ci附近以nl+nr+1個點在最小二乘意義下擬合一個M次多項式pi(c),多項式pi(c)在ci的值,即光滑數值gi表示為

(8)

式中:nl為ci左邊點的個數;nr為ci右邊點的個數;bp為多項式的系數。設實測數據為yi,為了使pi(c)擬合測試數據,必須定義系數bp,使得式(9)達到最優。

(9)

3 閾值決策降噪仿真分析

實驗信號為Lorenz方程產生的混沌信號,如圖3所示。

(10)

其中參數分別為:a=10,r=34,c=8/3。

(a) 原始時間序列

(b) 信噪比SNR=5 dB的含噪序列

圖4是Lorenz時間序列的EEMD結果,包括7個IMF分量(IMF1~7)和一個余項(res),根據式(6),(7)計算索引值js=4,即IMF1~3為高頻噪聲主導,IMF4~7為低頻信號主導,因此,EEMD閾值決策降噪的步驟具體為:

1) 對IMF1依據“3σ”法則和SVD提取有用信號細節d1。圖5為IMF1經SVD后的特征值,可以看出,特征值分布比較均勻,利用差分譜選取特征值個數無法取得理想效果,而先利用“3σ”法則估計IMF1中信號分量的能量,為SVD特征值的選取提供先驗知識,則可達到較好的降噪效果。去噪后的信號d1如圖6所示,相比于圖1,提取的信號比較均勻,更為合理,由于提取的信號僅占第一個IMF分量的0.003,因此,信號幅度很小。

圖4 Lorenz信號EEMD結果

圖5 IMF1SVD的特征值

圖6 去噪后的IMF1

2) 利用SVD提取高頻分量的有用細節。由于分解得到的IMF分量頻率按2的指數降冪排列,IMF分量頻率驟降,利用連續均方誤差準則確定高低頻分界點js=4,即噪聲起主導作用的模態為IMF2、IMF3,圖7為IMF2、IMF3經SVD后的特征值,可以看出,前面幾個特征值比較大,后面的比較均勻,利用差分譜選取特征值個數即可取得理想效果,重構選定的奇異值,得到d2、d3,如圖8所示。

(a) IMF2 SVD的特征值

(b) IMF3SVD的特征值

(a) 去噪后的IMF2

(b) 去噪后的IMF3

3) 對剩余IMF分量利用Savitaky-Golay濾波算法進行平滑處理。SG算法的擬合階數取3,數據窗口為15,得到d4~d8。

4) 重構信號x_denoised(t)=d1+d2+…+d8。

為了對比降噪效果,對Lorenz含噪序列分別采用小波加權閾值去噪方法[14]、傳統EEMD和所提方法進行降噪處理,小波加權閾值去噪方法基小波選取“sym8”小波,分解層數為5;傳統EEMD降噪中,采用硬閾值方法,3種方法對5 dB的Lorenz含噪序列的去噪效果如圖9所示。對于模型已知的Lorenz時間序列,采用信噪比SNR和均方誤差MSE來評價降噪效果,SNR越大,MSE越小,降噪效果越好。為了全面檢驗3中方法對不同噪聲強度的含噪信號的去噪效果,分別計算對SNR=0 dB、5 dB、10 dB、15 dB的Lorenz含噪序列降噪后的SNR和MSE,詳見表1。

由圖9可以看出,小波加權閾值法和EEMD方法去噪后的信號含有很多毛刺,整體不夠平滑,后者尤為突出,而本文方法去噪后的信號較為光滑平整更接近原始信號。由表1的計算結果看出,本文所提方法整體上要優于小波加權閾值和EEMD方法,并且對較低信噪比的Lorenz含噪序列尤為突出,隨信噪比逐漸增大,3種方法的差距減小,這和2.1節中的分析一致,也是由于低強度噪聲對混沌時間序列的軌跡造成干擾較小,傳統的非線性降噪方法都比較有效。

(a) 小波加權閾值法

(b) EEMD方法

(c) 本文所提方法

降噪方法SNR/MSE0dB5dB10dB15dB小波閾值加權法14.273/0.225716.536/0.167820.305/0.07624.578/0.032EEMD14.095/0.237516.475/0.176220.257/0.08324.557/0.037本文方法15.237/0.209717.335/0.109420.828/0.05524.732/0.025

4 機械式混沌振動信號處理

混沌振動是混沌科學研究的重要課題,在混沌隔振、混沌振動壓路機等方面有著廣泛的應用,國內外學者圍繞混沌振動的應用、控制和識別等問題開展了深入的研究。但是,在混沌振動信號采集過程中,難免受到周圍環境噪聲的干擾,嚴重影響了關聯維數、Lyapunov指數等混沌特征量的計算,為混沌參數區域確定和混沌控制的研究帶來了困難。

為了進一步驗證本文所提方法對于模型未知的混沌信號的去噪效果,本文基于雙勢阱理論的單端磁吸式混沌振動試驗裝置產生的振動信號為對象,雙勢阱單端磁吸式混沌振動裝置如圖10所示。

實驗本質為正弦信號的慢速頻率掃描實驗,掃描頻率范圍為5~25 Hz,采樣頻率為2 kHz,數據采集時長為5 s。調節功率放大器增益為1、激勵頻率為13 Hz時,重構相空間參數中嵌入維數為3,延遲時間為10,得到的重構相圖如圖11(a)所示。分別使用小波閾值加權法、EEMD和本文方法對采集的振動信號進行降噪處理。

圖10 雙勢阱單端磁吸式混沌振動裝置

以及本文方法對原始信號進行降噪處理,去噪后的二維相圖如圖11(b)、(c)、(d)所示,通過比較可以看出,本文方法降噪后的二維相圖曲線更加光滑,更能清晰地展現原信號吸引子的幾何結構。由于原始信號為機械式振動信號,其信號和噪聲未知,不能利用SNR和MSE進行定量比較,但是對于混沌信號,其自相關函數值比較大,且遠遠大于噪聲的自相關函數值,因此,可以利用3種方法去噪后的自相關函數值以評價去噪效果。

(a)原始信號(b)小波加權閾值方法(c)EEMD方法(d)本文方法

圖11 振動信號去噪前后二維相圖

Fig.11 The phase diagrams of vibration signals

3種方法去噪前后的部分自相關函數值如表3所示,可以看出,原始信號受到噪聲的干擾,其自相關函數值相比去噪后要小很多,而本文所提方法降噪后序列的自相關函數值最大,進一步展現了其降噪的優越性。

通過掃頻試驗,不斷調節激振器激勵頻率(以步長0.5 Hz)和功率放大器的增益(以步長0.1 V),利用本文方法對采集到的信號進行降噪處理,然后分別計算降噪前后的信號的Lyapunov指數(LE)和關聯維數。當激勵增益為1時,部分去噪前后信號特征量如表3所示,混沌運動的LE為正值且關聯維數為分數維,去噪后信號的兩個特征量都證實系統運動為混沌運動,而f=14 Hz、15 Hz時,原始信號的關聯維數接近整數維,說明了噪聲干擾了系統特征量的計算,也證實了本文方法為實現混沌振動信號進一步分析及其應用提供了重要的前提條件。通過對本文方法降噪后信號特征量計算判斷振動的運動形式,并以激勵頻率為橫坐標,激勵增益為縱坐標作出系統不同響應的參數區域圖,如圖12所示,從圖中可以明顯獲得混沌運動的參數區域。

表2 3種方法去噪后部分自相關函數值

表3 部分實測信號的特征指數

5 結 論

基于總體平均經驗模態分解,提出一種閾值決策的去噪方法。信號經總體模態分解得到固有模態函數分量,首先根據“3σ”法則和SVD提取第一個IMF分量中有用信號細節。然后,利用連續均方誤差準則對剩余IMF分量進行高低頻區分,分別使用SVD和S-G算法提取高低頻分量的有用信號,可以有效避免了高頻部分有用信號的流失,同時剔除低頻分量中的部分噪聲,克服了EEMD去噪時IMFs選擇的難題。最后通過對多普勒仿真信號和實測振動信號進行去噪分析,結果表明,本文方法是非常有效的。

圖12 系統不同響應的參數區示意圖

致謝:感謝楊慶超講師和Southampton University Solent Institution of Acoustic: Jian Jiang James, Chris, Lee, Lawrance的討論。

[1] 趙志宏,楊紹普,申永軍. 一種改進的EMD降噪方法[J]. 振動與沖擊,2009,28(12):35-37.

ZHAO Zhihong,YANG Shaopu,SHEN Yongjun. Improved EMD based de-noising method J]. Journal of Vibration and Shock,2009,28(12):35-37.

[2] HUANG N E,SHEN Z,LONG S R,et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings of the Royal Society of London Series A-Mathematical Physical and Engineering Science,1998,454:903-995.

[3] MANUEL B V,WENG B W,BARNER K E. ECG signal denoising and baseline wander correction based on the empirical mode decomposition[J]. Computers in Biology and Medicine,2008,38(1):1-13.

[4] 王紅軍,萬鵬.基于EEMD和小波包變換的早期故障敏感特征獲取[J]. 北京理工大學學報,2013,33(9):945-950.

WANG Hongjun,WAN Peng. Sensitive features extraction of early fault based on EEMD and WPT[J].Transactions of Beijing Institute of Technology,2013,33(9):945-950.

[5] 甘雨,隋立芬,王冰. 經驗模態分解閾值消噪方法及其在慣性導航系統數據處理中的應用[J].測繪學報,2012,41(4):504-509.

GAN Yu,SUI Lifen,WANG Bing. EMD threshold de-noising and its applications in INS data processing[J].Acta Geodaetica et Cartographica Sinica,2012,41(4):504-509.

[6] 陳忠,符和超. 基于EEMD方法的永磁同步電機電磁噪聲診斷[J]. 振動與沖擊,2013,32(20):124-128.

CHEN Zhong,FU Hechao. Electromagnetic noise diagnosis for a permanent magnet synchronous motor based on EEMD[J]. Journal of Vibration and Shock,2013,32(20):124-128.

[7] 柏林,劉小峰,秦樹人. 小波-形態-EMD綜合分析法及其應用[J]. 振動與沖擊,2008,27(5):1-4.

BO Lin,LIU Xiaofeng,QIN Shuren. Hybrid wavelet-morphology-EMD analysis and its application[J]. Journal of Vibration and Shock,2008,27(5):1-4.

[8] 陳仁祥,湯寶平,馬靖華. 基于EEMD的振動信號自適應降噪方法[J]. 振動與沖擊,2012,31(15):82-86.

CHEN Renxiang,TANG Baoping,MA Jinghua. Adaptive de-noising method based on ensemble empirical mode decomposition for vibration signal[J]. Journal of Vibration and Shock,2012,31(15):82-86.

[9] WU Z,HUANG N E. A study of the characteristics of white noise using the empirical mode decomposition method[J]. Proc. R. London A,2004,460:1597-1611.

[10] 陳雋,李想. 運用總體經驗模態分解的疲勞信號降噪方法[J].振動、測試與診斷,2011,31(1):15-19.

CHEN Juan, LI Xiang. Application of Ensemble Empirical Mode decomposition to noise reduction of fatigue signal[J].Journal of Vibration, Measurement & Diabnosis,2011,31(1):15-19.

[11] 雷達,鐘詩圣.基于奇異值分解和經驗模態分解的航空發動機健康信號降噪[J].吉林大學學報(工學版),2013,43(3):764-770.

LEI Da,ZHONG Shisheng. Aircraft engine health signal denoising based on singular value decomposition and empirical mode decomposition methods[J]. Journal of Jilin University (Engineering and Technology Edition),2013,43(3):764-770.

[12] 王文波,張曉東,汪祥莉. 基于主成分分析的經驗模態分解消噪方法[J]. 電子學報,2013,41(7):1425-1430.

WANG Wenbo,ZHANG Xiaodong,WANG Xiangli.Empirical mode decomposition de-noising method based on principal component analysis[J]. Acta Electronica Sinica,2013,41(7):1425-1430.

[13] 張文忠,周蓉,武旭紅,等. 利用白噪聲分解特征的EEMD閾值降噪方法[J].測繪科學技術學報,2013,30(3):255-259.

ZHANG Wenzhong,ZHOU Rong,WU Xuhong,et al. A threshold de-noising method based on the characteristics of white noise decomposed by EEMD[J]. Journal of Geomatics Science and Technology,2013,30(3):255-259.

[14] 位秀雷,林瑞霖,劉樹勇,等.基于改進小波變換的混沌信號去噪研究[J].武漢理工大學學報(交通科學與工程版),2013,37(5):1062-1065.

WEI Xiulei,LIN Ruilin,LIU Shuyong,et al. Denoising for chaotic signal based on an improved wavelet transform method[J].Journal of Wuhan University of Technology(Transportation Science & Engineering),2013,37(5):1062-1065.

Ade-noisingmethodforchaoticsignalsbasedonimprovedEEMD

WEI Xiulei1, LIN Ruilin2, LIU Shuyong2, YANG Qingchao2

(1. NO.91404 troops of PLA , Qinhuangdao, 066000, China; 2. College of Power Engineering,Naval University of Engineering,Wuhan 430033, China)

In de-noising process using the ensemble empirical mode decomposition (EEMD), the effective treatment of intrinsic mode functions (IMFs) is a key affecting noise reduction effect. Here, an improved EEMD de-noising method was proposed. Firstly, the useful signal details of the first IMF were extracted based on the “3σ” criterion and the singular value decomposition (SVD). Then the remaining IMFs were divided into higher frequency components and lower ones based on the consecutive mean square error (CMSE). Secondly, useful signals in higher frequency components and lower ones were extracted based on SVD and Savitzky-Golay(S-G) filtering method, respectively. Thus, the loss of useful signals in higher frequency components was avoided while parts of noise in lower frequency components were removed effectively to overcome the shortcoming of IMFs being difficult to treat with EEMD to do de-noising. In order to evaluate the effectiveness of the proposed method, the test rig of leaf spring based on the double-potential well theory was made and the test results showed that the proposed method is better than the wavelet weighted parameters method and the EEMD de-noising one.

ensemble empirical mode decomposition (EEMD); chaotic signal; singular value decomposition (SVD); denoising; Savitzky-Golay (S-G) filtering

國家自然科學基金(51579242;51179197);國家自然科學基金青年基金(51509253);海軍工程大學科研基金(425517K143)

2016-04-18 修改稿收到日期:2016-07-16

位秀雷 男,博士生,1988年生

劉樹勇 男,博士,副教授,碩士生導師,1975年生 E-mail:wxlcln@163.com

TN911

: A

10.13465/j.cnki.jvs.2017.17.006

猜你喜歡
模態細節振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
以細節取勝 Cambridge Audio AXR100/ FOCAL ARIA 906
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
留心細節處處美——《收集東·收集西》
中立型Emden-Fowler微分方程的振動性
細節取勝
Coco薇(2016年10期)2016-11-29 19:59:58
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 黄色网址免费在线| 婷婷六月综合| 免费无码AV片在线观看国产| 国产95在线 | a国产精品| 自拍亚洲欧美精品| 久久国产V一级毛多内射| 99无码中文字幕视频| 少妇极品熟妇人妻专区视频| 一级毛片高清| 成人国产精品2021| 久久青草精品一区二区三区 | 国产免费精彩视频| 无码'专区第一页| 无码精品国产dvd在线观看9久 | 亚洲h视频在线| 激情爆乳一区二区| 国产在线欧美| www欧美在线观看| 日本欧美视频在线观看| 精品国产www| av在线5g无码天天| 狠狠亚洲五月天| 日本不卡在线视频| 中文字幕精品一区二区三区视频 | 亚洲一级毛片免费观看| 色久综合在线| 国产v欧美v日韩v综合精品| 亚洲专区一区二区在线观看| 中文毛片无遮挡播放免费| 国产精品尹人在线观看| 无码又爽又刺激的高潮视频| 老司机久久99久久精品播放| 一本大道视频精品人妻| 99re热精品视频国产免费| 99视频在线精品免费观看6| 久久精品欧美一区二区| 亚洲综合婷婷激情| 日韩精品成人在线| 欧美亚洲网| 黄色在线不卡| 欧美久久网| 2020最新国产精品视频| 国产波多野结衣中文在线播放| 国产H片无码不卡在线视频| 91亚瑟视频| 看国产毛片| AV片亚洲国产男人的天堂| 国产一在线观看| 国产精品永久免费嫩草研究院| 亚洲v日韩v欧美在线观看| 国产乱人伦精品一区二区| 国产精品亚欧美一区二区| 欧美成人午夜影院| 91极品美女高潮叫床在线观看| 美女亚洲一区| 99在线观看国产| 欧美日韩动态图| 国产在线精彩视频论坛| 国产一级在线观看www色 | 亚洲va欧美ⅴa国产va影院| 日韩高清成人| 亚洲无码91视频| 欧美福利在线观看| 国产永久免费视频m3u8| 欧美亚洲国产精品久久蜜芽| 亚洲无码精彩视频在线观看| 亚洲美女久久| 欧美一级高清片久久99| 欧美日本在线观看| 国产精品天干天干在线观看| 毛片最新网址| 国产日韩精品欧美一区灰| 国产一级视频在线观看网站| 欧美日韩精品一区二区在线线| 91免费国产高清观看| 伊人AV天堂| 亚洲妓女综合网995久久| 国产精品9| 亚洲aaa视频| 国产99视频精品免费视频7| 精品国产成人国产在线|