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

基于多特征融合的刀具磨損識別方法*

2014-02-19 04:18:42石志標
振動、測試與診斷 2014年3期
關鍵詞:特征提取特征融合

關 山, 石志標, 劉 炎

(東北電力大學機械工程學院 吉林,132012)

引 言

刀具狀態監測技術作為保障工件表面質量和尺寸精度,防止工件報廢和機床損壞,優化加工過程,降低成本,提高生產效率的有效手段,越來越引起人們的重視。許多學者提出了多種監測方法[1-4],其中,多傳感器融合法[5-8]在一定程度上克服了單一傳感器獲取信息量有限、抗干擾能力差的局限性,大大提高了監測結果的準確率,是最為有效的監測方法。但是,基于多傳感器融合的監測系統勢必帶來成本的增加,干涉機床操作的靈活性,甚至影響機床的性能,從而影響實際監測系統的推廣使用。基于采樣信號多特征融合[9]的監測技術克服了多傳感器融合的缺點,拓展了監測系統的可靠性和應用范圍。筆者通過采集刀具不同磨損階段的聲發射信號,分別從信號的時域、頻域、時-頻兩域進行特征提取,構造聯合多特征向量,然后采用核主元分析法對聯合多特征向量進行融合,剔除冗余特征,以最小二乘支持向量機作為分類器實現了刀具磨損狀態的識別,取得了較滿意的結果。

1 實驗系統及實驗方法

實驗系統如圖1所示,傳感器采用PXR30諧振式聲發射傳感器,諧振頻率為300kHz,帶寬為80kHz~400kHz。較高的頻率能有效地接收到感興趣信號的高頻成分,又可濾除低頻噪聲,較寬的頻帶有利于實驗過程中寬頻信號的采集。前置放大器為PXPAⅡ寬帶聲發射放大器,帶寬為15kHz~2MHz,增益為40dB。預處理實現對聲發射信號帶通濾波。采用PCI-1721數據采集卡,利用Lab-VIEW軟件編寫數據采集程序完成數據采集,采樣頻率為1MHz。實驗材料為高溫合金GH4169,刀片為肯納公司的KC9125硬質合金涂層刀片,采用CKA6136i數控車床進行車削實驗。

圖1 聲發射數據采集系統示意圖Fig.1 The diagram of acoustic emission data collection system

實驗的目的是為了研究變切削條件下刀具磨損狀態的分類問題。如果將實驗所選用的3種因素(切削速度、進給量、切削深度)的所有參數(水平)進行全面組合,形成多種切削條件,不僅會導致實驗量過大,而且由于切削速度對刀具壽命的影響要遠遠大于進給量和切削深度對刀具壽命的影響,雖然后續對采樣數據進行了去均值及歸一化處理,以削除切削條件變化對刀具磨損的影響,過大的切削速度變化范圍仍會給變切削條件下刀具狀態的分類及磨損量預測帶來較大的誤差。因此,將所選用的切削速度參數依據相近原則分成3組,如表1所示。對應于每一組各選取進給量和切削深度這兩個因素的3種切削用量,根據日本學者田口玄一提出的正交實驗設計法,設計了3組3因素3水平正交實驗。

表1 各組正交實驗所選用的切削速度Tab.1 Cutting speed of each orthogonal experiment r/min

根據正交實驗表,選定一切削條件,在這一確定的切削條件下,實驗方法如下:

1)取新刀片1進行切削實驗,切削10s后停車,僅采集切削過程中6~10s間的數據,取下刀片,測量磨損量VB值;

2)更換新刀片2,本次切削20s后停車,只記錄15~20s間的數據,取下刀片,測量VB值;

3)再次更換新刀片,切削時間比上一次再增加10s,僅記錄本次切削時間內最后5s的數據,停車,取下刀片,測量VB值;

4)步驟3反復進行,每次切削均更換新刀片,切削時間均比上次增加10s,直到切削時間累加到足夠長,新刀片在這個切削時間內進行切削能夠磨損為止,本切削條件下的切削實驗終止;

5)選取另一切削條件,重復1,2,3,4步實驗過程,直至完成全部選定切削條件下的切削實驗,切削實驗結束。

由于每次實驗都采用新刀片連續切削一定時間使刀片達到一定的磨損量,這與實際的切削過程相符。刀具磨損是漸變的過程,所以每次切削過程最后5s的實驗數據才能最真實地反應測量所得VB值所對應的刀具磨損狀態,并大大減少了采集的數據量,為后續數據處理提供了方便。采用每次切削實驗切削時間逐步遞增的方式,便于標定切削時間、刀具磨損量及信號特征之間的一一對應關系,為采用最小二乘支持向量機回歸算法實現刀具磨損量預測準備了數據條件。

圖2所示為在切削速度為560r/min、進給量為0.3mm/r、背吃刀量為0.4mm 時,刀具不同磨損階段采樣信號的時序圖(從上至下分別對應的VB值為0.11,0.13,0.17,0.24,0.26和0.31mm)。

圖2 刀具不同磨損階段采樣信號時序圖Fig.2 The time series diagram of different wear stages sample signal

根據VB值,筆者將刀具磨損分成如表2所示的4個等級。

表2 刀具磨損等級分類表Tab.2 The tool wear grade classification

2 聯合多特征向量的構造

2.1 基于EMD與AR模型的時域特征提取

經驗模態分解(empirical mode decomposition,簡稱EMD)是Huang[10]提出的一種自適應的、不需預先確定分解基的信號處理方法,通過EMD分解,可以將復雜的非平穩信號以有限個固有模態函數(intrinsic mode functions,簡稱IMF)之和的形式表示出來,各IMF分量包含了原信號中不同時間尺度的局部特征信息,實現了非平穩信號的平穩化處理。

時間序列的 AR(atuo-regressive)模型參數[11]凝聚了系統狀態的重要信息,對系統狀態變化的反映最為敏感。但AR模型僅對平穩過程具有較好的分析效果,刀具磨損過程中產生的聲發射信號表現為較強的非平穩特征,直接對其進行AR建模效果不理想,EMD分解正好解決了這一問題。因此,基于EMD與AR模型的信號時域特征提取步驟如下。

1)對采樣信號進行EMD分解。

圖3為在560r/min,0.3mm/r,0.4mm/r的切削條件下,刀具磨損量VB值分別為0.11(左列)、0.26(中列)、0.33(右列)時所采集的聲發射信號經EMD分解的結果。

圖3 不同磨損狀態信號經EMD分解結果Fig.3 The EMD decomposition results of different wear state

分析結果表明,不同切削條件、刀具不同磨損狀態下采集的數據經EMD分解后得到IMF分量的個數是不同的,如圖3所示。為使提取特征向量的維數一致,利用式(1)計算各IMF分量與原始信號的相關系數

其中:cov(X,Y)為兩序列的協方差;D(X),D(Y)為序列X和序列Y的方差。

計算結果表明,前7階IMF分量與原信號的相關性較好,攜帶原信號的主要信息,其余的IMF分量與原信號相關性較差,體現原信號中的噪聲成分,如圖4所示。筆者通過舍棄后面的IMF達到信號去噪、降維、歸范化的目的。

2)建立前7階IMF分量的AR模型。

采用最小信息準則函數,通過實例計算確定4階為最佳模型階次,分別建立各IMF分量的4階AR模型,提取模型的自回歸系數構成式(2)所示的28維特征向量。建模前先對各IMF分量進行歸一化處理,以消除切削條件變化的影響。

圖4 各階IMF分量與原始信號的相關性Fig.4 The correlation of each IMF component with the original signal

其中:φm,n(m=1~7;n=1~4)表示m階IMF分量4階AR模型的第n個系數。

2.2 基于雙譜分析的頻域特征提取

雙譜是高階譜分析的一個特例,它涵蓋了傳統功率譜分析不能表征的信息,這些信息完全可以顯示非高斯信號的特征,且均值為零的高斯過程,其雙譜為零,因此能有效地消除高斯噪聲。雙譜分析作為非平穩、非高斯信號特征提取的有效手段,已應用于刀具磨損狀態檢測中[12]。

設{x(t)}為均值為零的k階平穩隨機過程,其k階累積量ckx(τ1,τ2,…,τk-1)是絕對可和的,即

則{x(t)}的k階譜定義為k階累積量k-1維傅里葉變換

其中:ω=[ω1,ω2,…,ωk-1]T;τ=[τ1,τ2,…,τk-1]T。

其中,三 階 譜S3,x(ω1,ω2)叫 做 雙 譜,記 為Bx(ω1,ω2),定義為

筆者采用參數化雙譜估計法[13]計算雙譜,基于雙譜分析的頻域特征提取步驟如下。

1)對采樣信號進行去均值及歸一化處理,以消除切削條件變化的影響

2)對歸一化后信號進行雙譜分析。

圖5所示為進給量(0.3mm/r)和切削深度(0.4mm)相同、切削速度不同時(a,b,c三列對應的切削速度分別為560,630和800r/min),刀具不同磨損階段的雙譜分析圖。

從雙譜分析等高線圖中看出:在相同磨損階段,磨損信號的雙譜差異較小;而在不同的磨損階段,采樣信號的雙譜存在明顯的差別,且在三種切削速度下均可以區分出來。這說明信號經去均值及歸一化處理后的雙譜特征可以有效地實現變切削條件下刀具磨損狀態的分類。

3)基于奇異值分解的雙譜特征提取。

奇異值是矩陣的固有特征,當信號無噪聲或具有較高的信噪比時,奇異值隨著秩的增大而迅速減小;而當信號完全由噪聲成分構成時,矩陣是列滿秩的,各奇異值幾乎相等。因此,矩陣的奇異值中,前面若干個比較大,反映信號中的特征成分,其余的值較小,對應信號中的噪聲成分。定義奇異譜為)

其中:λ1>λ2>…>λr≥0為矩陣的奇異值。

奇異譜分量pi表示各狀態變量在整個系統中所占能量的相對關系,通過保留前s個奇異譜的方法提取特征向量并降維。奇異譜數量s的選取一般根據式(8)確定

采樣數據經雙譜分析后,得到復數矩陣。對矩陣中的元素取模,構造一個以模為元素的實矩陣MR,對MR進行奇異分解。計算結果表明:在不同切削條件、不同磨損狀態下,能量累積貢獻率大于85%的奇異譜數量并不相同,一般在10~15之間。筆者統一選取前15個奇異值譜構造特征向量如下

圖6為p1-p2,p1-p3的二維分布散度圖,由圖6可以看出,所提取特征的聚類效果是很明顯的。

2.3 基于小波包變換的能量特征提取

對不同切削條件、不同磨損階段所采集的信號進行傅里葉分析,發現信號的頻率主要集中在200kHz以下,并且隨著刀具磨損的加劇,信號在該頻段內的幅值有較大的變化。圖7所示為在切削條件為140r/min,0.3mm/r,0.4mm,VB值分別為0.146和0.347mm時,磨損信號傅里葉變換的頻譜圖。

圖5 雙譜分析等高線圖Fig.5 The bispectral analysis contour map

圖6 p1-p2,p1-p3 的散度圖Fig.6 The divergence chart of p1 -p2and p1 -p3

同時計算信號從低頻到高頻能量的累積占總能量的百分比,發現絕大多數情況下信號在0~200kHz內的能量占總能量的百分比超過85%,個別情況下在0~230kHz內超過85%。本研究中信號的采樣頻率為1MHz。根據采樣定理,采樣頻率ωs和連續信號最高頻率ωmax之間的關系必須滿足ωs≥2ωmax,從而得出采集有效信號的最高頻率為500kHz。根據小波包分解理論,對采樣信號進行4層小波包分解,可將信號分解為16個頻段,那么前8個頻段將包括信號中頻率為0~256kHz的范圍。所以只取前8個頻段的信號作為研究對象,就可以包含原信號中的絕大多數有用信息。基于小波包變換的能量特征的提取步驟如下。

圖7 刀具2種磨損狀態信號的頻譜圖Fig.7 Two wear states signal spectrum

1)對采樣信號進行4層小波包分解,提取前8個頻段的小波包分解系數并重構,得到各頻帶內的時域信號S4,0~S4,7。

2)求各頻帶信號的能量e4,j(j=0,1,…,7)

其中:xjk(j=0,1,…,7;k=0,1,…,n)表示重構信號s4,j的離散點的幅值。

3)以8個頻段內信號的能量為元素構造特征向量T,則T表示為考慮到變切削條件下刀具磨損狀態的監測問題,筆者對能量特征進行了歸一化處理

在 切 削 條 件 為 220r/min,0.3mm/r,0.51mm,VB 值分 別為 0.083,0.143,0.252 和0.324mm時,小波包分解前8個頻段的能量圖如圖8所示。

圖8 刀具不同磨損狀態AE信號的頻帶-能量圖Fig.8 The band-energy bar of different tool wear AE signal

2.4 構造聯合多特征向量

1)采集對應不同切削條件、不同磨損狀態下的聲發射信號,讀取VB值,建立信號與刀具磨損狀態的對應關系。

2)構造如式(2)所示的28維特征行向量Tar,如式(9)所示的15維特征行向量Tho和如式(12)所示的8維特征向量Twa。

3)將Tar,Tho和Twa依次首尾相連,構造如式(13)所示的51維聯合多特征向量t

3 基于核主元分析的融合特征構造

核主元分析[14]是一種常用的高維數據線性降維和特征提取方法,它通過構造一組新的潛隱變量來降低原始數據空間的維數,再從新的映射空間中抽取主要變化的信息,提取統計特征。新映射空間中的變量由原始數據變量的線性組合構成,彼此正交,消除了變量間的關聯性,最大限度地攜帶原變量的有用信息,從而降低了投影空間的維數。

本研究以表3所列的9種切削條件及表2所示的刀具磨損等級為例說明融合特征的構造方法。

表3 實驗所選切削參數Tab.3 Experimental cutting parameters

在同一切削條件下,采集刀具從初期磨損到劇烈磨損的樣本50組,9種切削條件下共選取450組樣本。分別構造聯合多特征向量,組成一個450×51的聯合特征向量矩陣M,矩陣的每行代表一個樣本的聯合特征。對矩陣M進行核主元分析,得到51×51的主成分系數矩陣Acoeff

Acoeff矩陣中的第i列是第i個主成分的系數向量。同時計算樣本協方差矩陣的特征值向量,它是由51個特征值構成的列向量,特征值按降序排列,對應每一主成元對整體貢獻的大小。由圖9可以看出,前10個主元的累積貢獻率超過85%,因此,選取前10個主元,這樣主成分系數矩陣Acoeff變為如下的51×10的矩陣

對于一個新的樣本,按照聯合特征向量的構造方法,構造出51維聯合特征向量tnew=[t1,t2,…,t51],然后向新的主元投影,就得到降維后的融合特征向量Tnew∈R1×10,這樣特征向量從51維降為10維,投影算法如下

圖9 主元的貢獻率及累積貢獻率Fig.9 Contribution rate of the principal component and cumulative contribution rate

圖10為融合特征向量Tnew中不同主元的二維散度圖。由圖10可以看出,經KPCA降維后的融合特征具有較好的聚類能力。

4 基于LS-SVM的刀具狀態分類

將融合特征向量分為二組,一組用于訓練,一組用于驗證。兩組樣本數都為200個,其中A類樣本20個,其余每種類別各60個。

采用 LS-SVM[6,15]進行分類,徑向基核函數參數σ和懲罰參數γ的選擇直接影響所建模型的識別精度。筆者采用交叉驗證法優化參數,表4列出了選用不同的σ和γ值對分類準確率的影響。在選定優化參數的情況下,各類別的識別率及總體識別率如表5有淺色底紋行所示。

表4 采用不同的參數σ和γ對分類準確率的影響Tab.4 The impact on the classification accuracy with different parametersσandγ

為了驗證融合特征的優劣,筆者分別采用EMD法、雙譜分析法、小波分析法等提取的單一特征以及聯合特征送入LS-SVM進行識別,結果如表5所示。識別結果表明:a.采用融合特征進行識別具有最高的識別準確率;b.比較融合特征與聯合特征的識別率可以得出,主元分析可以有效地消除各變量間的冗余,去除次要特征,有效減少了特征的維數,減少了后續訓練和識別的時間,避免了“維數災難”,提高了識別的準確率;c.A類樣本的識別率均最低,這是由于在切削過程中,刀具在初始階段磨損較快,采集到的A類樣本較少,且磨損量多分布在0.067~0.083mm之間,樣本覆蓋范圍較窄,導致訓練階不夠充分,從而影響整體識別率,但是這類樣本的誤識別并不影響刀具的使用性能。進一步比較LSSVM與SVM及BP神經網絡、RBF神經網絡在分類方面的性能,采用融合特征分別送入以上分類器,分類結果如表6所示。

圖10 不同主元的散度圖Fig.10 The divergence chart of different principal component

表5 選用不同特征的識別率比較Tab.5 The recognition rate comparison of different characteristics

表6 選用不同識別方法的識別率比較Tab.6 The recognition rate comparison of different methods

結果表明,采用不同的識別方法,B,C和D類的識別結果大至相近,但是LS-SVM,SVM對A類的識別率明顯高于RBF和BP網絡。這表明LSSVM,SVM在小樣本的識別上要優于神經網絡。通過對不能正確分類樣本的進一步分析發現,這類樣本多集中在磨損分類的過渡階段,不論采用哪種特征提取方法,這類樣本的特征都比較接近。

5 結 論

1)筆者提出的融合特征與單一特征或聯合特征相比,可以更有效地刻畫刀具磨損過程產生的聲發射信號,得到最高的識別率。在不增加監測成本及設備安裝復雜性的基礎上,采用融合特征實現刀具狀態監測,為實用刀具磨損監測系統的開發研制提供了一條可借鑒的途徑。

2)徑向基核函數參數σ和懲罰參數γ的選擇直接影響所建LS-SVM模型的分類識別精度。在小樣本的情況下,LS-SVM方法用于刀具磨損狀態識別在實例計算中要優于神經網絡識別方法,得到較滿意的結果。

3)錯誤識別的樣本大多集中在刀具磨損的過渡階段,為了提高識別的準確率,找到更為有效的信號特征提取方法將是下一步研究的主要目標。

[1] Zhu Kunpeng,Wong Y S,Hong G S.Wavelet analysis of sensor signals for tool condition monitoring:a re-view and some new results[J].International Journal of Machine Tools & Manufacture,2009,49(7):537-553.

[2] Dutta S,Datta A,Das Chakladar D,et al.Detection of tool condition from the turned surface images using an accurate grey level co-occurrence technique[J].Precision Engineering,2012,36(3):458-466.

[3] Zhu Kunpeng,Wong Y S,Hong G S.Multi-category micro-milling tool wear monitoring with continuous hidden Markov models[J].Mechanical Systems and Signal Processing,2009,23(2)547-560.

[4] 康晶,馮長建,胡紅英.刀具磨損監測及破損模式識別[J].振動、測試與診斷,2009,29(1):5-9.

Kang Jing,Feng Changjian,Hu Hongying.Tool wear monitoring and pattern recognition of tool failure[J].Journal of Vibration, Measurement & Diagnosis,2009,29(1):5-9.(in Chinese)

[5] Cuneyt A,Metin Ertunc H,Hasan O.Tool wear condition monitoring using a sensor fusion mode based on fuzzy inference system[J].Mechanical Systems and Signal Processing,2009,23(2):539-546.

[6] Sultan B,Shihab A,Sohyung C,et al.Machine ensemble approach for simultaneous detection of transient and gradual abnormalities in end milling using multisensor fusion[J].Journal of Materials Processing Technology,2009,209(10):4728-4738.

[7] Ghosh N,Ravi Y B,Patra P,et al.Estimation of tool wear during CNC milling using neural network-based sensor fusion[J].Mechanical Systems and Signal Processing,2007,21(1):466-479.

[8] Ibrahim D,Khaled A,Firas H.On modeling of tool wear using sensor fusion and polynomial classifiers[J].Mechanical Systems and Signal Processing,2009,23(5):1719-1729.

[9] 鄭建明.基于HMM的多特征融合鉆頭磨損監測技術的研究[D].西安:西安理工大學,2004.

[10]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 A:Mathematical,Physical and Engineering Sciences,1998,454:903-995.

[11]王海麗,馬春翔,邵華,等.車削過程中刀具磨損和破損狀態的自動識別[J].上海交通大學學報,2006,40(12):2057-2062.

Wang Haili,Ma Chunxiang,Shao Hua,et al.The tool wear and breakage monitoring in turning using neural network[J].Journal of Shanghai Jiaotong University,2006,40(12):2057-2062.(in Chinese)

[12]李錫文,楊明金,謝守勇,等.銑削加工刀具磨損過程雙譜分析[J].農業機械學報,2007,38(9):143-146.

Li Xiwen,Yang Mingjin,Xie Shouyong,et al.Bispectral analysis of cutter wear process in milling operation[J].Transactions of the Chinese Society for Agricultural Machinery,2007,38 (9):143-146.(in Chinese)

[13]李勇,徐震.Matlab輔助現代工程數字信號處理[M].西安:西安電子科技大學出版社,2002:219-238.

[14]Elangovan M,Devasenapati S B,Sakthivel N R,et al.Evaluation of expert system for condition monitoring of a single point cutting tool using principle component analysis and decision tree algorithm[J].Expert Systems with Applications,2011,38(4):4450-4459.

[15]Elangovan M,Sugumaran V,Ramachandran K I,et al.Effect of SVM kernel functions on classification of vibration signals of single point cutting tool[J].Expert Systems with Application,2011,38(12):15202-15207.

猜你喜歡
特征提取特征融合
村企黨建聯建融合共贏
今日農業(2021年19期)2022-01-12 06:16:36
融合菜
從創新出發,與高考數列相遇、融合
《融合》
現代出版(2020年3期)2020-06-20 07:10:34
如何表達“特征”
基于Gazebo仿真環境的ORB特征提取與比對的研究
電子制作(2019年15期)2019-08-27 01:12:00
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
抓住特征巧觀察
一種基于LBP 特征提取和稀疏表示的肝病識別算法
基于MED和循環域解調的多故障特征提取
主站蜘蛛池模板: 毛片基地美国正在播放亚洲| 女人18毛片水真多国产| 国产精品手机视频| 欧美成人精品一级在线观看| 亚洲无码电影| 色呦呦手机在线精品| 欧美激情,国产精品| 亚洲国产成人精品青青草原| 国产熟女一级毛片| 日本人妻丰满熟妇区| 天天躁夜夜躁狠狠躁图片| 被公侵犯人妻少妇一区二区三区| 久久这里只有精品66| 国产福利在线免费观看| 欧洲高清无码在线| 91在线一9|永久视频在线| 3D动漫精品啪啪一区二区下载| 日韩小视频在线观看| 特级aaaaaaaaa毛片免费视频| 日本不卡视频在线| 亚洲成人在线免费| 91原创视频在线| 亚洲成a人片| 在线永久免费观看的毛片| 99激情网| 久久男人资源站| 天天做天天爱夜夜爽毛片毛片| 欧美yw精品日本国产精品| 精品久久久久久中文字幕女| AV不卡在线永久免费观看| 欧美19综合中文字幕| 99热国产在线精品99| 亚洲人成人无码www| 又大又硬又爽免费视频| 丰满人妻一区二区三区视频| www亚洲天堂| 国产玖玖玖精品视频| 久久精品丝袜| 欧美日韩导航| 欧美午夜理伦三级在线观看| 日日拍夜夜嗷嗷叫国产| 亚洲色无码专线精品观看| 国产SUV精品一区二区6| 无码有码中文字幕| 日本草草视频在线观看| 国产嫩草在线观看| 毛片基地视频| 色噜噜在线观看| 天堂va亚洲va欧美va国产| 在线观看国产黄色| 国产三级韩国三级理| 日韩欧美中文在线| 国产麻豆永久视频| 国产午夜精品鲁丝片| 亚洲欧美自拍视频| 乱色熟女综合一区二区| 欧美不卡视频一区发布| 无码国产伊人| 亚洲精品视频在线观看视频| 在线高清亚洲精品二区| 久久一级电影| 国产精品欧美激情| 粗大猛烈进出高潮视频无码| 亚洲午夜久久久精品电影院| 久热精品免费| 福利在线不卡一区| 亚洲国产看片基地久久1024| 国产精品久久精品| 国产亚洲美日韩AV中文字幕无码成人| 亚洲狠狠婷婷综合久久久久| 91精品国产一区| 亚洲 欧美 日韩综合一区| 国产精品尤物铁牛tv| 97国产精品视频自在拍| 国产精品自在在线午夜区app| 四虎成人精品| 91精品综合| 99热这里只有精品5| 日韩a级毛片| 中文字幕在线永久在线视频2020| 日韩中文无码av超清| 亚洲九九视频|