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

PCA與SVD信號(hào)處理效果相似性與機(jī)理分析

2016-07-26 09:04:40聶振國(guó)趙學(xué)智
振動(dòng)與沖擊 2016年2期

聶振國(guó), 趙學(xué)智

(華南理工大學(xué) 機(jī)械與汽車(chē)工程學(xué)院,廣州 510640)

?

PCA與SVD信號(hào)處理效果相似性與機(jī)理分析

聶振國(guó), 趙學(xué)智

(華南理工大學(xué) 機(jī)械與汽車(chē)工程學(xué)院,廣州510640)

摘要:將主成分分析(Principal Component Analysis,PCA)用于信號(hào)處理,并與奇異值分解(Singular Value Decomposition,SVD)方法比較。分析總結(jié)PCA及SVD信號(hào)處理原理,提出基于PCA的特征值差分譜理論用于信號(hào)消噪。結(jié)果表明,PCA與SVD的處理效果較相似,相似性原因?yàn)樵季仃囉移娈愊蛄考礊閰f(xié)方差矩陣特征向量。SVD較PCA的重構(gòu)誤差小,因SVD無(wú)需計(jì)算協(xié)方差矩陣,可避免舍入誤差產(chǎn)生。

關(guān)鍵詞:主成分分析;奇異值分解;消噪;相似性;誤差

主成分分析由Pearson[1]在生理學(xué)研究中用于分析數(shù)據(jù)及建立數(shù)理模型,通過(guò)對(duì)協(xié)方差矩陣進(jìn)行特征分解[2],獲得數(shù)據(jù)的主成分與權(quán)值。現(xiàn)PCA已用于地理、生物、經(jīng)濟(jì)、數(shù)理統(tǒng)計(jì)等眾多領(lǐng)域。

奇異值分解理論由Beltrami提出,之后圍繞奇異值分解算法已有諸多研究[3],如文獻(xiàn)[4]提出的針對(duì)大型矩陣奇異值分解的多次分割雙向收縮QR算法,能克服單向收縮QR算法收斂較慢、長(zhǎng)時(shí)停滯甚至不收斂等缺陷,處理大型矩陣數(shù)值計(jì)算時(shí)具有迭代次數(shù)少、過(guò)程無(wú)停滯及收斂迅速等優(yōu)點(diǎn),SVD已廣泛用于信號(hào)處理、圖像壓縮、數(shù)字水印、語(yǔ)音編碼、控制系統(tǒng)等。

SVD與PCA的概念不同,應(yīng)用領(lǐng)域、最終實(shí)現(xiàn)目標(biāo)亦不同。PCA通過(guò)對(duì)多個(gè)樣本構(gòu)成的矩陣提取主成分進(jìn)行降維,減小數(shù)據(jù)量以便于數(shù)據(jù)分析,且不會(huì)丟失原矩陣的主要特征;而SVD則通過(guò)對(duì)矩陣進(jìn)行分解,選較奇異值進(jìn)行重構(gòu)獲得原矩陣的主要信息,其中細(xì)節(jié)信息如噪聲被濾掉。

關(guān)于SVD信號(hào)處理的文獻(xiàn)較多[5-9],而關(guān)于PCA應(yīng)用的文獻(xiàn)尚少見(jiàn)。用信號(hào)構(gòu)造Hankel矩陣作SVD處理,可據(jù)重構(gòu)矩陣恢復(fù)出信號(hào);而經(jīng)PCA處理所得各主成分只代表原變量指標(biāo)在新變量指標(biāo)上所占比重,較原矩陣維數(shù)已降低,雖能保留原始數(shù)據(jù)的主要特征,無(wú)法從中恢復(fù)出原始數(shù)據(jù)。本文通過(guò)研究發(fā)現(xiàn),PCA獲得某一主成分后再乘以相應(yīng)特征向量可得類(lèi)似SVD的重構(gòu)矩陣,且重構(gòu)矩陣與SVD重構(gòu)矩陣間有一一對(duì)應(yīng)關(guān)系,信號(hào)處理效果相似,并對(duì)相似性進(jìn)行分析、探討。

1PCA與SVD基本理論

1.1基于PCA的信號(hào)處理原理

(1)

X可寫(xiě)成X=[x1x2… xm]。為消除量綱影響,將X化為零均值矩陣。

定義[10]:若x1,x2,…,xm為原變量指標(biāo),y1,y2,…,yl(l≤m)為新變量指標(biāo),則主成分分析即將m個(gè)原變量指標(biāo)由l個(gè)新變量指標(biāo)表示,即

(2)

式中:αij為系數(shù),須滿(mǎn)足

(3)

系數(shù)αij確定的原則如下:

(1) yi與yj(i≠j;i,j=1,2,…,l)相互無(wú)關(guān);

(2) y1為x1,x2,…,xm的一切線性組合中方差最大者,y2為與y1不相關(guān)的x1,x2,…,xm所有線性組合中方差最大者。以此類(lèi)推,yl為與y1,y2,…,yl-1均不相關(guān)的x1,x2,…,xm所有線性組合中方差最大者。則新變量指標(biāo)y1,y2,…,yl分別稱(chēng)原變量指標(biāo)x1,x2,…,xm的第1,2,…,l主成分。可證明αi=(α1i,α2i,,…,αmi)T為X協(xié)方差矩陣第i個(gè)特征值(從大到小排列)對(duì)應(yīng)的特征向量[11]。可獲得各主成分y1,y2,…,yl,若只用于數(shù)據(jù)降維分析,目的已達(dá)到,但若用于信號(hào)處理,則需重構(gòu)矩陣,據(jù)前述有

yi=Xαi,(i=1,2,…,m)

(4)

(5)

由于X的協(xié)方差矩陣為半正定對(duì)稱(chēng)陣,其特征向量間相互正交,有

(6)

原矩陣X的重構(gòu)表達(dá)式為

(7)

若選前l(fā)個(gè)主成分進(jìn)行重構(gòu),則得近似矩陣為

(8)

1.2基于SVD的信號(hào)處理原理

奇異值分解指對(duì)任意矩陣A∈Rn×m總存在正交矩陣U∈Rn×n及V∈Rm×m使

A=USVT

(9)

式中:S=[diag(σ1,σ2, …,σr),0]或其轉(zhuǎn)置,取決于n>m或n

設(shè)U=[u1,u2, …,un],V=[v1,v2, …,vn],其中ui∈Rn×1,vi∈Rm×1,則A可寫(xiě)成

A=σ1u1v1T+σ2u2v2T+…+σrurvrT

(10)

令A(yù)i=σiuiviT(i=1,2,…,r),按奇異值差分譜理論[13]選前l(fā)個(gè)分量進(jìn)行重構(gòu),得近似矩陣為

(11)

獲得近似重構(gòu)矩陣后亦需據(jù)矩陣構(gòu)造方式恢復(fù)出信號(hào),如文獻(xiàn)[7]中構(gòu)造Hankel矩陣進(jìn)行消噪處理,并給出具體信號(hào)恢復(fù)方法。

2PCA與SVD信號(hào)處理效果對(duì)比

2.1模擬信號(hào)去噪效果對(duì)比

2.1.1不含直流分量情況

以不含直流分量的模擬信號(hào)比較PCA與SVD的消噪效果,設(shè)信號(hào)為

f(t)=sin(2πt)+sin(4πt)+sin(20πt)

(12)

對(duì)該信號(hào)以采樣頻率200 Hz、采樣512點(diǎn),并疊加服從N(0,1)的高斯白噪聲,信噪比為1.705 3 dB,時(shí)域圖見(jiàn)圖1。以此信號(hào)構(gòu)造257×256的Hankel矩陣X,矩陣行列數(shù)越相近,去噪效果越好[14]。PCA用于信號(hào)去噪時(shí)若據(jù)貢獻(xiàn)率選主成分個(gè)數(shù),由于矩陣維數(shù)較大,貢獻(xiàn)率較高則主成分個(gè)數(shù)較多,會(huì)使去噪效果差。因此對(duì)PCA主成分個(gè)數(shù)采用類(lèi)似SVD的奇異值差分譜理論選取,并稱(chēng)其為特征值差分譜理論。區(qū)別在于,PCA用協(xié)方差矩陣特征值,而SVD用矩陣本身奇異值。本例特征值序列及差分譜見(jiàn)圖2(a),而奇異值序列及差分譜見(jiàn)圖2(b)。由圖2看出,奇異值較特征值大,但變化趨勢(shì)相同。

圖1 原始含噪信號(hào)Fig.1Theoriginalnoisysignal圖2 PCA特征值曲線與SVD奇異值曲線及差分譜Fig.2EigenvaluescurveofPCAandsingularvaluescurveofSVDandthedifferencespectrum

據(jù)特征值差分譜理論選前6個(gè)主成分進(jìn)行重構(gòu),獲得去噪信號(hào)見(jiàn)圖3(a),據(jù)奇異值差分譜理論選前6個(gè)分量進(jìn)行重構(gòu),得去噪信號(hào)見(jiàn)圖3(b)。由圖3看出,二者消噪效果非常相似。

以模擬調(diào)幅信號(hào)為例對(duì)比兩者消噪效果,設(shè)調(diào)幅信號(hào)為

y(t)=6sin(120πt+0.56)sin(12πt+0.78)

(13)

疊加服從N(0,1)分布的高斯白噪聲,信噪比為-4.194 1 dB。采樣頻率1 024 Hz,采樣1 024點(diǎn),其時(shí)域圖見(jiàn)圖4。

圖3 模擬信號(hào)PCA與SVD去噪效果對(duì)比Fig.3ThePCAandSVDmethodde-noisingeffectofanalogsignal圖4 含噪調(diào)幅信號(hào)Fig.4Thenoisyamplitudemodulationsignal

按文獻(xiàn)[14],構(gòu)造513×512的Hankel矩陣分別進(jìn)行PCA、SVD處理,特征值序列、奇異值序列及差分譜見(jiàn)圖5,消噪結(jié)果見(jiàn)圖6。

圖5 PCA特征值、SVD奇異值曲線及差分譜Fig.5 Eigenvalues curve of PCA and singular values curveof SVD and the difference spectrum

由圖3、圖6看出,PCA與SVD去噪效果幾乎無(wú)差異,不易分辨,但實(shí)際上其處理結(jié)果誤差卻不同,對(duì)式(12)信號(hào),可求得PCA法誤差均值為0.146 7,SVD法誤差均值為 0.129 0;而對(duì)式(13)調(diào)幅信號(hào),PCA法誤差均值為0.418 8,SVD法誤差均值為0.389 1。可見(jiàn)SVD誤差稍小。

PCA與SVD重構(gòu)誤差比較結(jié)果見(jiàn)表1及圖7。可見(jiàn),對(duì)每個(gè)信號(hào),SVD重構(gòu)誤差均小于PCA。原因?yàn)镻CA需計(jì)算協(xié)方差矩陣,且會(huì)產(chǎn)生舍入誤差,而SVD直接對(duì)原矩陣進(jìn)行處理,無(wú)需計(jì)算協(xié)方差矩陣,因而避免舍入誤差。

表1 PCA與SVD重構(gòu)誤差

圖6 調(diào)幅信號(hào)PCA、SVD法去噪效果對(duì)比Fig.6ThePCAandSVDmethodde-noisingeffectoftheamplitudemodulationsignal圖7 PCA與SVD誤差比較Fig.7TheerrorofPCAandSVD

2.1.2含直流分量情況

以含直流分量的模擬信號(hào)比較PCA與SVD的信號(hào)消噪效果,在模擬信號(hào)式(12)基礎(chǔ)上疊加幅值為2的直流分量,得信號(hào)為

f=2+sin(2πt)+sin(4πt)+sin(20πt)

(14)

對(duì)此信號(hào)以采樣頻率200 Hz、采樣512點(diǎn),并疊加服從N(0,1)的高斯白噪聲。其時(shí)域圖見(jiàn)圖8。構(gòu)造257×256的Hankel矩陣X,用PCA及SVD對(duì)此信號(hào)進(jìn)行消噪,所得特征值序列及差分譜圖見(jiàn)圖9(a),奇異值序列及差分譜見(jiàn)圖9(b)。對(duì)比圖9、圖2看出,加直流分量信號(hào)與未加直流分量信號(hào)的特征值差分譜圖相同,而奇異值差分譜圖變化較大。

圖8 原始含噪信號(hào)Fig.8Theoriginalnoisysignal圖9 PCA特征值曲線、SVD奇異值曲線及差分譜Fig.9EigenvaluescurveofPCAandsingularvaluescurveofSVDandthedifferencespectrum

添加直流分量后,雖Hankel矩陣每列均值增加2,但由于PCA進(jìn)行零均值化處理,在零均值化過(guò)程中每列每個(gè)元素減去均值,因此,無(wú)論添加的直流分量幅值大小,最終所得特征值差分譜均與無(wú)直流分量特征值差分譜圖相同。因此,利用特征值差分譜確定有用分量個(gè)數(shù)時(shí),無(wú)論信號(hào)是否含直流分量均可簡(jiǎn)單以特征值差分譜最大峰值點(diǎn)確定有用分量個(gè)數(shù),此為特征值差分譜理論。而按文獻(xiàn)[13]奇異值差分譜理論,當(dāng)用SVD進(jìn)行信號(hào)消噪時(shí),信號(hào)中若含直流分量,應(yīng)剔除奇異值差分譜圖中第1點(diǎn)影響,再取奇異值差分譜峰值最大點(diǎn)確定有用分量個(gè)數(shù),以此進(jìn)行矩陣重構(gòu)。

據(jù)以上理論,此例中按照特征值差分譜理論,應(yīng)選前6個(gè)分量進(jìn)行重構(gòu),所得PCA去噪效果見(jiàn)圖10(a);按奇異值差分譜理論用SVD去噪時(shí),應(yīng)選前7個(gè)分量進(jìn)行重構(gòu)方能獲得較好去噪效果,見(jiàn)圖10(b)。

2.2實(shí)際信號(hào)消噪效果對(duì)比

以某軸承振動(dòng)信號(hào)為例,用PCA法、SVD法進(jìn)行消噪。原始信號(hào)見(jiàn)圖11,可見(jiàn)信號(hào)中有較大噪聲污染。用PCA、SVD對(duì)其進(jìn)行消噪處理,PCA特征值序列及差分譜見(jiàn)圖12(a),SVD奇異值序列及差分譜見(jiàn)圖12(b)。據(jù)差分譜理論,分別取前2個(gè)特征值、奇異值進(jìn)行重構(gòu),所得重構(gòu)結(jié)果見(jiàn)圖13。可見(jiàn),經(jīng)PCA、SVD消噪后,信號(hào)已變干凈,噪聲完全消去。

圖10 模擬信號(hào)PCA法與SVD法去噪效果對(duì)比Fig.10ThePCAandSVDmethodde-noisingeffectofanalogsignal圖11 原始信號(hào)Fig.11Theoriginalsignal

圖12 PCA特征值、SVD奇異值曲線及差分譜Fig.12 Eigenvalues curve of PCA and singular values curve of SVD and the difference spectrum

圖13 實(shí)際信號(hào)PCA、SVD法去噪效果對(duì)比Fig.13 The PCA and SVD method de-noising effect of analog signal

3PCA、SVD信號(hào)處理效果相似原因分析

無(wú)論模擬信號(hào)或?qū)嶋H工程信號(hào),由PCA、SVD重構(gòu)結(jié)果看出兩者相差無(wú)幾。對(duì)該相似性的內(nèi)在機(jī)理分析如下:

PCA重構(gòu)的某分量可寫(xiě)為

Xi=yiαiT=XαiαiT

(15)

SVD某分量形式為

(16)

(17)

對(duì)比式(15)、(17)知,PCA、SVD所得分量矩陣Xi具有相似形式,但存在向量αi與vi的區(qū)別。分析αi與vi的關(guān)系,X的協(xié)方差矩陣為

式中:cov(xi,xj)=E[(xi-E(xi))(xj-E(xj))]。

X的奇異值分解為

X=USVT

(19)

據(jù)式(19)可得

XTX=(USVT)TUSVT=

VSTUTUSVT=VSTSVT

(20)

(21)

據(jù)式(21)可知X的右奇異向量vi即為XTX的特征向量,而X奇異值的平方即為XTX的特征值。對(duì)X進(jìn)行零均值處理后,有E(xi)=0(i=1,2,…m),因此X的協(xié)方差矩陣可寫(xiě)為

(22)

αi=vi

(23)

結(jié)合式(15)、(17)、(23)可見(jiàn),PCA及SVD所得重構(gòu)矩陣Xi相等,此為PCA及SVD的消噪效果相似性的本質(zhì)原因;但PCA需計(jì)算協(xié)方差矩陣,存在數(shù)據(jù)舍入誤差,而SVD無(wú)需計(jì)算協(xié)方差矩陣,且無(wú)需進(jìn)行零均值化處理,SVD計(jì)算量少,能避免計(jì)算協(xié)方差矩陣的舍入誤差,相對(duì)誤差較小。而PCA在處理帶直流分量的信號(hào)時(shí),因其進(jìn)行零均值化處理時(shí)能消除直流分量影響,故無(wú)、有直流分量的特征值差分譜完全相同,即無(wú)需考慮有無(wú)直流分量。

由此可見(jiàn),雖PCA、SVD信號(hào)處理效果相似,但SVD處理結(jié)果誤差較PCA小,且SVD計(jì)算量少,故工程實(shí)際應(yīng)用中應(yīng)用SVD去噪。

4結(jié)論

通過(guò)分析PCA、SVD信號(hào)處理算法,并實(shí)例對(duì)比,結(jié)論如下:

(1) 兩者結(jié)果相似原因?yàn)樵季仃囉移娈愊蛄考礊閰f(xié)方差矩陣特征向量。信號(hào)中含直流分量、用PCA去噪時(shí),特征值差分譜與無(wú)直流分量相同,只需據(jù)特征值差分譜最大峰值確定有用分量個(gè)數(shù),無(wú)需考慮信號(hào)有無(wú)直流分量。

(2) SVD誤差較PCA小的原因?yàn)椋琍CA在計(jì)算協(xié)方差矩陣時(shí)會(huì)產(chǎn)生舍入誤差,而SVD無(wú)需計(jì)算協(xié)方差矩陣,避免舍入誤差產(chǎn)生。

參 考 文 獻(xiàn)

[1] Pearson K. On lines and planes of closest fit to systems of points in space[J]. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science,1901,2(11): 559-572.

[2] Abdi H, Williams L J. Principal component analysis[J]. Wiley Interdisciplinary Reviews: Computational Statistics, 2010, 2(4): 433-459.

[3] Stewart G W. On the early history of the singular value decomposition[J]. SIAMReview, 1993, 35(4): 551-566.

[4] 趙學(xué)智,葉邦彥,陳統(tǒng)堅(jiān). 大型矩陣奇異值分解的多次分割雙向收縮 QR算法[J]. 華南理工大學(xué)學(xué)報(bào):自然科學(xué)版, 2010, 38(1):1-8.

ZHAO Xue-zhi, YE Bang-yan, CHEN Tong-jian. Multi-partition and double-direction shrink QR algorithm for singular value decomposition of large-scale matrix[J]. Journal of South China University of Technology: Natural Science Edition, 2010, 38(1): 1-8.

[5] 趙學(xué)智,葉邦彥,陳統(tǒng)堅(jiān). 基于 SVD 的奇異性信號(hào)檢測(cè)原理及其應(yīng)用[J]. 振動(dòng)與沖擊, 2008, 27(6): 11-14.

ZHAO Xue-zhi, YE Bang-yan, CHEN Tong-jian. Principle of singularity detection based on SVD and its application[J]. Journal of Vibration and Shock, 2008, 27(6): 11-14.

[6] 段向陽(yáng),王永生,蘇永生.基于奇異值分解的信號(hào)特征提取方法研究[J]. 振動(dòng)與沖擊, 2009,28(11): 30-33.

DUAN Xiang-yang, WANG Yong-sheng, SU Yong-sheng. Feature extraction methods based on singular value decomposition[J]. Journal of Vibration and Shock, 2009,28(11): 30-33.

[7] 張波,李健君. 基于Hankel 矩陣與奇異值分解 (SVD) 的濾波方法以及在飛機(jī)顫振試驗(yàn)數(shù)據(jù)預(yù)處理中的應(yīng)用[J]. 振動(dòng)與沖擊, 2009, 28(2): 162-166.

ZHANG Bo, LI Jian-jun. De-noising method based on hankel matrix and SVD and its application in flight flutter testing data preprocessing[J]. Journal of Vibration and Shock, 2009,28(2): 162-166.

[8] 陳恩利,張璽,申永軍,等. 基于SVD降噪和盲信號(hào)分離的滾動(dòng)軸承故障診斷[J]. 振動(dòng)與沖擊, 2012, 31(23): 185-190.

CHEN En-li, ZHANG Xi, SHEN Yong-jun, et al. Fault diagnosis of rolling bearings based on SVD denoising and blind signals separation[J]. Journal of Vibration and Shock, 2012,31(23): 185-190.

[9] 王建國(guó),李健,劉穎源. 一種確定奇異值分解降噪有效秩階次的改進(jìn)方法[J]. 振動(dòng)與沖擊, 2014(12): 176-180.

WANG Jian-guo, LI Jian, LIU Ying-yuan. An improved method for determining effective order rank of SVD denoising[J]. Journal of Vibration and Shock, 2014(12): 176-180.

[10] 王峻峰. 基于主分量獨(dú)立分量分析的盲信號(hào)處理及應(yīng)用研究[D]. 武漢:華中科技大學(xué), 2005.

[11] Liang Y C, Lee H P, Lim S P, et al. Proper orthogonal decomposition and its applications-part I: theory[J]. Journal of Sound and Vibration, 2002, 252(3): 527-544.

[12] Golub G H, Van Loan C F,袁亞湘. 矩陣計(jì)算[M]. 北京:人民郵電出版社, 2011.

[13] 趙學(xué)智,葉邦彥,陳統(tǒng)堅(jiān). 奇異值差分譜理論及其在車(chē)床主軸箱故障診斷中的應(yīng)用[J].機(jī)械工程學(xué)報(bào),2010(1):100-108.

ZHAO Xue-zhi, YE Bang-yan, CHEN Tong-jian. Difference spectrum theory of singular value and its application to the fault diagnosis of headstock of lathe[J]. Journal of Mechanical Engineering,2010 (1): 100-108.

[14] 趙學(xué)智,葉邦彥,陳統(tǒng)堅(jiān). 矩陣構(gòu)造對(duì)奇異值分解信號(hào)處理效果的影響[J]. 華南理工大學(xué)學(xué)報(bào):自然科學(xué)版, 2008, 36(9):86-93.

ZHAO Xue-zhi, YE Bang-yan, CHEN Tong-jian. Influence of matrix creation way on signal processing effect of singular value decomposition [J]. Journal of South China University of Technology: Natural Science Edition, 2008, 36(9): 86-93.

基金項(xiàng)目:國(guó)家自然科學(xué)基金項(xiàng)目(51375178);廣東省自然科學(xué)基金項(xiàng)目(S2012010008789)

收稿日期:2014-11-19修改稿收到日期:2015-01-20

通信作者趙學(xué)智 男,教授,博士生導(dǎo)師,1970年11月生

中圖分類(lèi)號(hào):TN91117; TH16513

文獻(xiàn)標(biāo)志碼:A

DOI:10.13465/j.cnki.jvs.2016.02.003

Similarity of signal processing effect between PCA and SVD and its mechanism analysis

NIE Zhen-guo, ZHAO Xue-zhi

(School of Mechanical and Automotive Engineering, South China University of Technology, Guangzhou 510640, China)

Abstract:The principal component analysis (PCA) was applied to signal processing and its effect was compared with that of the singular value decomposition (SVD). The signal processing principles of PCA and SVD were analyzed and summarized, and the theory of eigenvalues difference spectrum based on PCA for signal denoising was introduced. It is pointed out that PCA has a very similar signal processing effect to that of SVD when applied to signal de-noising. The reason for this similarity was analyzed theoretically, and it is found that this is because the right singular vectors of the original matrix are just the eigenvectors of its covariance matrix, and leads to the similarity between PCA and SVD in signal processing. It is also pointed out that the reconstruction error of SVD is smaller than that of PCA, and the reason is that SVD does not need to compute the covariance matrix, so the rounding errors are avoided.

Key words:principal component analysis; singular value decomposition; de-noising; similarity; error

第一作者 聶振國(guó) 男,碩士生,1991年10月生

主站蜘蛛池模板: 久久无码av一区二区三区| 国产情精品嫩草影院88av| 四虎成人精品| 97成人在线视频| 国产va视频| 亚洲一区毛片| 国产91蝌蚪窝| 伊在人亞洲香蕉精品區| 精品成人一区二区| 亚洲精品国产综合99| 亚洲国产亚综合在线区| 国产特级毛片| 免费国产高清精品一区在线| 日韩欧美中文在线| 亚洲AV成人一区国产精品| 国产成年女人特黄特色毛片免| 亚洲无卡视频| 国产成人AV综合久久| 国产精品欧美激情| 亚洲精品无码在线播放网站| 国产交换配偶在线视频| 在线亚洲小视频| 久草视频一区| 国产91熟女高潮一区二区| 九九久久99精品| 99热这里只有成人精品国产| 3p叠罗汉国产精品久久| 性色在线视频精品| 日本午夜视频在线观看| 色综合五月婷婷| 成人伊人色一区二区三区| 亚洲综合亚洲国产尤物| 亚洲欧美一级一级a| 亚洲精品无码日韩国产不卡| 亚洲AV一二三区无码AV蜜桃| 亚洲av无码专区久久蜜芽| 国产福利在线免费观看| 伊人中文网| 国产福利一区视频| 久久国产精品麻豆系列| 国产95在线 | 亚洲欧美日本国产综合在线 | 伊人网址在线| 国产精品女熟高潮视频| 看av免费毛片手机播放| 亚洲成人网在线观看| 青青国产在线| 超碰精品无码一区二区| 亚洲日韩精品伊甸| 亚洲男人天堂久久| 亚洲国产精品日韩av专区| 亚洲精品免费网站| 天堂网亚洲综合在线| 欧美不卡二区| 喷潮白浆直流在线播放| 国产精品久久精品| 欧美另类精品一区二区三区| 婷婷亚洲综合五月天在线| 99热这里只有精品久久免费| 亚洲无线一二三四区男男| 草草线在成年免费视频2| 国产精品自在在线午夜区app| 亚洲大尺度在线| 久久这里只有精品66| 亚洲成a人片77777在线播放| 国产亚洲视频播放9000| 欧美中文字幕第一页线路一| 国产成人在线无码免费视频| 蜜芽国产尤物av尤物在线看| 亚洲精品自在线拍| 激情国产精品一区| 国产精品蜜臀| 99视频免费观看| 麻豆精品久久久久久久99蜜桃| 国产综合网站| 九九热视频在线免费观看| 亚洲第一色网站| 欧美成人看片一区二区三区| 亚洲成人精品在线| 韩国v欧美v亚洲v日本v| 欧美特黄一免在线观看| 99久久国产综合精品女同|