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

基于盲源信號分離的加工誤差分離方法研究

2016-11-10 08:01:18張發平吳迪張體廣張凌雲楊吉彬
兵工學報 2016年9期
關鍵詞:測量區域方法

張發平,吳迪,張體廣,張凌雲,楊吉彬

(北京理工大學機械與車輛學院,北京100081)

基于盲源信號分離的加工誤差分離方法研究

張發平,吳迪,張體廣,張凌雲,楊吉彬

(北京理工大學機械與車輛學院,北京100081)

針對當前加工誤差分離方法無法分離相近尺度系統誤差的缺點,提出了一種基于獨立成分分析的加工誤差分離方法。建立了各誤差源所導致的系統誤差和總的系統誤差之間的傳遞模型;根據盲源信號分離相關理論,建立了基于負熵固定點算法的加工誤差分離模型,實現了相近尺度的系統加工誤差分離,并利用主成分分析法給出了誤差源數量確定的方法。以某陀螺儀的某加工面加工誤差分離為例對所提出的方法進行了驗證,結果表明所提出的誤差分離方法可有效實現相近尺度的加工誤差分離。

機械制造工藝與設備;誤差源;誤差傳遞模型;獨立成分分析;主成分分析;誤差分離

0 引言

在武器裝備零件加工過程中,由于工藝系統諸多誤差源(如機床熱變形、主軸回轉誤差、刀具磨損、振動等)的存在,使得被加工工件表面存在加工誤差,而且不同誤差源所造成的誤差表現形式不同。如何根據加工誤差的不同表現形式,對加工誤差進行分離,確定加工誤差的組成和作用規律,進行消除誤差和進行誤差補償,是提高裝備精密零件加工精度的基礎。

誤差辨識模型精度直接決定了所得幾何誤差源的準確性,進而影響誤差補償效果[1]。目前加工誤差的分離辨識模型主要采用基于濾波的方法和統計的方法。根據加工誤差的尺度(頻率)大小,將工件表面誤差分為3類[2]:表面粗糙度、表面波紋度、形狀誤差。因此可以將加工誤差按照尺度的大小進行濾波,分析影響加工質量的因素并進行誤差溯源。基于濾波思想的誤差分離方法已有許多研究。文獻[3-5]都是基于小波變換的方法,將加工表面形貌誤差分解到不同的尺度空間進行誤差成分分析。Chi等[6]提出表面粗糙度的頻率成分構成,并利用多譜分析方法對測量的表面粗糙度進行分析。李世平等[7]通過仿真實驗,驗證了經驗模態分解(EMD)方法可以將線性誤差、指數誤差、冪函數誤差以及周期誤差分離開。李成貴等[8]利用EMD方法對超精密研磨表面的綜合誤差分解為不同頻率范圍內的分量,通過分析每個分量的Hilbert譜,進一步分析了加工過程中異常頻率產生的原因,實現對加工過程誤差的溯源。

基于統計的加工誤差分離方法根據加工誤差的表現形式和規律不同,將加工誤差分離為系統誤差和隨機誤差。陳岳坪等[9]利用空間統計學的方法,對加工表面回歸模型殘差進行空間獨立性分析,分解出系統誤差和隨機誤差,并對系統誤差進行在線補償。張磊等[10]針對傳統的誤差流模型僅適用于線性模型的缺陷,建立加工誤差流的半參數回歸模型,將工序系統誤差和隨機誤差分離出來。Suriano等[11]基于加工表面的高清晰度測量數據,提出一種改進的序貫策略檢測方法,對表面上局部的加工缺陷進行識別。

上述兩類加工誤差分離方法都存在一定的缺陷,當加工誤差是由相近頻率的誤差構成時,濾波方法難以將其分離;而基于空間統計學的誤差分離僅限于將加工誤差分離為系統誤差和隨機誤差,不能判斷系統誤差的構成成分。本文根據盲源信號分離的原理,提出一種基于獨立成分分析的誤差分離方法,建立了誤差源到加工誤差之間的傳遞方程,給出了利用基于負熵的固定點算法估計誤差源的方法,從而有效地分離了系統誤差。

1 加工誤差傳遞模型

1.1加工誤差的表示

工藝系統(機床、夾具、刀具、工件)在加工狀態下的誤差行為(幾何誤差、運動誤差、安裝誤差、受力變形、受熱變形)會通過各環節轉換為刀具與工件相對位置的偏差,使得實際加工表面和理論加工表面出現偏差,形成了加工誤差。加工誤差一般可定義為加工表面的設計參數與加工后測量值的偏離,滿足

式中:M表示實際加工后所測得的表面;N表示理論設計表面;E為表面加工誤差。

1.2加工誤差的傳遞

加工誤差是由工藝系統各種誤差行為(誤差源)產生的,這些誤差經過工藝系統多個部分傳遞,最終引起刀具-工件相對位置的偏差,產生了加工誤差。因此加工誤差可以以誤差源產生的誤差作為輸入,加工表面誤差為輸出的傳遞函數

式中:si表示第i個誤差源;n為誤差源的數目。不同誤差源所產生的誤差特性不同,例如,主軸回轉誤差的主要成分是周期性信號,低階諧波(1~4階)占有很大的比重[12]。本文把所研究的加工表面劃分為m個區域,這m個區域可以是根據走刀次數來劃分或者其他方式,并假設在每個區域內的誤差源相同,其作用模式也相同。那么在第i個區域時由第j個誤差源sj獨立產生的加工誤差Eij可以表示為

式中:aij為誤差傳遞系數,數值大小代表誤差源對表面加工誤差的貢獻;bij為常數,代表誤差傳遞的偏移量。則第i(i=1,…,m)個區域所有誤差源形成的加工誤差Ei可以表示為

工件所有表面加工區域的加工表面誤差E可表示為

誤差源的傳遞系數矩陣;S=[s1,s2,…,sm]T,為誤差源誤差特性矩陣為誤差傳遞偏移值矩陣。(5)式表示了誤差源產生的誤差成分S與總的加工誤差之間的關系,只有E能通過對走刀路徑的測量和計算得到,是各個誤差源綜合作用的結果,而其余都是未知量。S是由不同的誤差源產生的誤差成分,其誤差特性各不相同,并且S的分量之間是相互獨立的。本文的目的正是利用這些特性,設計算法從E的測量值來估計各誤差源的誤差sj,以及誤差源的數量,從而實現誤差的分離,為誤差溯源提供基礎。

2 加工誤差的獨立成分分析分離方法

在加工誤差E已知,誤差的傳遞系數矩陣A和誤差源S未知的條件下,加工誤差S分離可以用盲源信號分離(BBS)的方法來處理。BBS[13]方法是一種在信號源和混合過程未知的情況下,根據觀測數據來確定分離過程的參數,使分離后的各信號是各源信號的拷貝或估計的方法。獨立成分分析(ICA)是其中常用的一種方法,其原理是:如果信號源之間相互獨立,通過尋找合適的分離矩陣W,使分離后的變量之間相互獨立,那么他們就是源信號的拷貝(分離信號和源信號存在幅值和順序的差異)[14]。誤差源是由工藝系統不同的部位產生,因此誤差源產生的誤差是相互獨立的,可以將ICA方法應用于本文所述的加工誤差分離中。

2.1誤差源數量估計

為了保證分離結果的準確,一般要求測量樣本的數目大于誤差源的數目,因此在進行分離之前,需要根據樣本的信息來估計含有的誤差源數量。由于主成分分析能使輸出變量之間相互獨立,對測量樣本做主成分分析,一般而言,統計能夠表征原數據95%以上信息的主成分個數,就可認定樣本中誤差源的數目,而其余5%的信息為包含了噪聲等干擾信息。確定誤差源的數目對測量樣本的數量以及樣本測量方式的選擇、工件表面測量區域劃分都有一定的指導作用。在此詳細說明:

設E1,E2,…,Em為測量得到的每個區域的加工誤差,每個區域內的誤差幅值可認為是隨機過程,第i個區域內p個測量點的數值可以看作是長度為p的向量,記為xi=[xi1,xi2,…,xip],那么Ei=[x1,x2…xp].m個測量區域中全部測量值構成的樣本集合記為矩陣X=[E1,E2,…,Em]。首先,計算協方差矩陣,設每一列的均值為g,則樣本集合的協方差矩陣其次,計算協方差矩陣cov X的特征值λ1,λ2,…,λp,特征值的累計貢獻率貢獻率的大小表明了攜帶原始特征變異信息的大小,貢獻率越大,表明它對應的特征向量對原特征信息的解釋能力越強。當Ci大于0.95時,即前i個主成分能夠包含原始誤差數據95%的信息,即包含了i個不相關的主要誤差成分,因此可以推斷出原始測量數據中誤差源的數目有i個,即取n=i.

2.2預處理

中心化處理。ICA算法要求測量誤差的均值為0,因此需要在分離之前對測量誤差進行零均值處理,即中心化處理。對加工誤差進行中心化處理可按照(6)式進行:

式中:E[E(x)]為測量誤差E(x)的數學期望。~E(x)是測量數據中心化后的結果。在實際中,由于測量到的加工誤差數據樣本長度N是有限的,因此可以用樣本數據的平均值代替期望,即

白化處理。為了保證測量誤差之間不相關,減少待估計參數的個數,降低問題的分析難度,需要對中心化后的數據進行白化處理。白化處理是一種數據預處理方法,其通過將信號的協方差矩陣的對角化處理,使向量元素之間不相關。白化處理將去均值后的測量誤差~E(x)進行線性變換P,得到

式中:V(x)即白化后的結果,且滿足E[~E(x)~ET(x)]= I;P為變換矩陣。求解方法如下,設Cx為去除均值的測量誤差的協方差矩陣,Cx=E[~E(x)~ET(x)],將其進行特征值分解,可得Cx=FDFT,

那么

式中:F表示特征矢量矩陣;D表示對角線上為特征值的對角矩陣。F和D可通過對Cx對角化處理得到,在此不詳述。

2.3誤差分離

對測量誤差數據進行中心化和白化之后,根據誤差源的數量構造初始n×n階隨機分離矩陣W. ICA的目的是逐步尋優,迭代計算出分離矩陣W,使分離后的誤差源S=[s1,s2,…,sn]對應的估計向量Y(x)=[y1(x),y2(x),…,yn(x)],滿足

En是從Em中隨機取的n個誤差;Y(x)、E(x)為測量位置的函數。對(10)式進行變換可得

通過設定滿足Y(x)各變量之間獨立性的目標函數,優化分離矩陣W,實現Y(x)各分量之間的統計獨立性最大化,從而獲得對誤差源S的估計。而Y(x)各分量之間的統計獨立性最大化可通過基于Y(x)的負熵最大化來判定。

2.4基于負熵最大化的獨立性判定

根據(10)式中的誤差分離模型,要求分離后的誤差變量y1(x),y2(x),…,ym(x)盡可能獨立,因此需要一個準則來度量分離后誤差變量之間的獨立性。根據中心極限定理:一定條件下,多個獨立分布變量的和趨向于高斯分布。加工誤差是由各個誤差源獨立作用產生,測量得到的誤差是誤差源分量的線性組合,所以其高斯型要比誤差源的高斯型要強。在具有相同方差的隨機變量中,高斯分布的隨機變量具有最大的信息熵,而非高斯性越強,信息熵越小,從而負熵越大。因此,負熵最大化意味著分離后的誤差源估計值y1(x),y2(x),…,ym(x)相互獨立程度最大,根據BBS理論,可以認為此時的y1(x),y2(x),…,ym(x)是對誤差源(s1,s2,…,sm)的一個估計[15]。

設誤差源的估計矩陣Y(x)的概率密度為p(Y),那么Y(x)的信息熵H定義為

信息熵是誤差變量高斯性的信息度量,高斯性越大,熵值越大。負熵N(Y)定義為式中:YGauss是具有與Y相同方差的高斯變量,可以由Y的分布計算獲得。(11)式需要給出概率密度的估計函數p(Y),計算較為困難,為了避免該問題,可采用(13)式的近似形式進行計算。

式中:JC(Y)為N(Y)的近似值,表示分離后誤差源估計值的負熵,作為優化的目標函數;G(·)是任意的非二次函數;C是一任意常數;yGauss是均值為0,方差為1的高斯變量。

如前所述,誤差分離就是尋找合適的分離矩陣W,使Y=WV具有最大的非高斯性,即負熵最大化。根據Kuhn-Tucker條件,在‖W‖2=1的約束下,E{G(WV)}的最優解滿足

式中:β為常數;g(·)為G(·)的導數;β可用E{WV g(WV)}來近似得到。利用牛頓迭代法求解(15)式從而對W進行更新:

迭代執行(15)式,直到W收斂。此時的W矩陣即為加工誤差的分離矩陣,代入(10)式中,即可獲得誤差源S的估計值Y,從而實現誤差的分離。

3 算例驗證

本部分以某陀螺儀的某加工面為研究對象,利用ICA算法分析其表面形貌誤差的特性。本實例將ICA算法應用于三維的加工表面,對其進行誤差分析和診斷。具體的操作步驟如下。

3.1加工數據準備

本實例的利用CMM對陀螺儀的某加工面進行測量,該面上X、Y兩個方向的采樣頻率分別是:X方向為1.363 6個點每毫米,Y方向為0.6個點每毫米。得到點云數量為150×30的平面,如圖1所示。

圖1顯示的是零件表面的測量值,要得到該表面的形貌誤差,需要用測量值減去基準。本文利用最小二乘法獲取該表面的基準,如圖2所示。利用測量值減去該基準之后,得到的表面形貌誤差如圖3所示。

在得到表面形貌誤差之后,將三維表面劃分成多個區域,每個區域作為測量信號。這里有兩個問題需要解決:一是區域的個數;二是如何劃分區域。對于第一個問題,通過主成分分析來確定獨立變量的個數,從而確定區域的個數。對于第二個問題,該表面是沿Y方向銑削加工得到的,被加工表面沿著Y方向記錄下了完整的加工信息。而X方向則沒有完整的記錄結果。因此,將表面形貌沿著X方向進行區域劃分。為了確定區域的個數,將區域的個數從2到150(X方向最大數目)逐一遞增,并計算每次劃分區域后的獨立成分個數,結果如圖4所示。

圖1 某表面的測量結果Fig.1 Measurement of a certain surface

圖2 最小二乘基準面Fig.2 Least square datum

圖3 表面形貌誤差Fig.3 Profile error of surface

圖4 表面形貌誤差中獨立成分個數隨區域數目變化圖Fig.4 Independent component number of surface profile error vs.region number

從圖4可以看出:不同的區域劃分方法,得到的是不一樣的樣本輸入以及輸入長度,因此對應不一樣的獨立成分個數。在所有的計算得到的獨立成分個數中,個數5是最多的,說明樣本中包含5個獨立成分的概率最大;另外,隨著區域的增加,每個區域在X方向越來越窄,包含的信息也就越少;因此,這里選擇第一次出現5個獨立成分時對應的區域數10,此時每個區域包含的信息量也最大。

3.2獨立成分分析

將這10個區域作為ICA的輸入,獨立成分個數為5,得到對5個獨立誤差成分的估計,結果圖5所示。

這5個誤差的估計在三維空間表達不夠清楚,不能在直觀上展示各自的特性以及彼此之間的不同。為克服這些缺點,將這5個獨立成分的估計沿Y方向展成一維模式。利用Fourier變換并計算每個獨立成分的頻譜,得到如圖6所示的結果。

從圖6中上面的一維展開圖可以明顯地看出每種獨立成分在形態上的不同,每種獨立成分對應了一種誤差源特征。獨立成分頻譜圖則更為清晰地表示了不同獨立成分的異同之處:從頻譜峰值讀數可以看出,全部獨立成分的主要頻率成分為0.09 Hz及其倍頻(2~7倍頻);雖然構成的主要成分相同,但是每一種獨立成分又具有自己固有特點:例如第3個獨立成分的頻譜可以看出,峰值較陡,能量集中;第2個獨立成分的頻譜則峰值處相對平緩,能量發散;其余各個獨立成分也具有各自固有而與其他成分區別明顯的頻譜特征。這些差異反映了不同的誤差源獨特的表現形式,印證了誤差源相互獨立作用這一特性。

本實例很好地說明了ICA能夠將三維加工表面根據誤差源特性的不同,實現分離。進一步分析分離出的結果,結合實際的加工狀況,可以分析出表面誤差中誤差產生的來源,這也是下一步研究的重點。

圖5 5個獨立誤差成分的估計Fig.5 Estimation of five independent error components

4 結論

本文基于BBS理論對加工表面誤差進行分析。通過對加工過程的誤差形成進行分析,建立了從誤差源到加工表面誤差的線性傳遞模型;根據主成分分析法給出了誤差源數目的估計方法;采用負熵最大化的判定準則,提出了基于負熵的誤差分離算法來求得誤差源的估計,從而實現了加工誤差的分離。由于該方法不涉及誤差間的尺度問題,只要相互獨立就可以把相近尺度誤差進行分離,所以在誤差分離方面有一定優勢。最后以某陀螺儀的某加工面誤差分離實驗驗證了該方法的有效性。該方法為實現加工誤差溯源奠定了基礎。

圖6 獨立成分估計的一維表示及其頻譜Fig.6 One-dimensional representation of independent components and their frequency spectra

(References)

[1] 田文杰,牛文鐵,常文芬,等.數控機床幾何精度溯源方法研究[J].機械工程學報,2014,50(7):128-135.TIANWen-jie,NIU Wen-tie,CHANG Wen-fen,et al.Research on geometric error tracing of NC machine tools[J].Journal of Mechanical Engineering,2014,50(7):128-135.(in Chinese)

[2] 蔣莊德,趙卓賢.形狀誤差、波度和表面粗糙度劃分的譜分析方法[J].計量學報,1989,10(3):170-174. JIANG Zhuang-de,ZHAO Zhuo-xian.A spectrum analysis method to distinguish form error,waviness and surface roughness[J].Acta Metrologic Sinica,1989,10(3):170-174.(in Chinese)

[3] Rosenboom L,Kreis T,JüptnerW.Surface description and defect detection by wavelet analysis[J].Measurement Science and Technology,2011,22(4):45102-45110.

[4] 李海燕,劉國棟,劉炳國,等.雙密度小波在表面形貌信號分離中的應用[J].光學精密工程,2008,16(6):1094-1097. LI Hai-yan,LIU Guo-dong,LIU Bing-guo,et al.Application of double density wavelet transform to surface topographic signal separation[J].Optics and Precision Engineering,2008,16(6): 1094-1097.(in Chinese)

[5] Liew K M,Wang Q.Application of wavelet theory for crack identification in structures[J].Journal of Engineering Mechanics,1998,124(2):152-157.

[6] Chi FC,Lee W B.A multi-spectrum analysis of surface roughness formation in ultra-precision machining[J].Precision Engineering,2000,24(1):77-87.

[7] 李世平,付宇,張進.一種基于EMD的系統誤差分離方法[J].中國測試,2011,37(3):9-15,36. LI Shi-ping,FU Yu,ZHANG Jin.Systematic error separation method based on EMD[J].China Measurement and Test,2011,37(3):9-15,36.(in Chinese)

[8] 李成貴,李寶貴,孫丹.基于HHT變換的超光滑加工表面譜特征分析[J].北京航空航天大學學報,2006,32(11):1341-1344. LI Cheng-gui,LI Bao-gui,SUN Dan.Application of HHT analysis in surface roughness[J].Journal of Beijing University of Aeronautics and Astronautics,2006,32(11):1341-1344.(in Chinese)

[9] 陳岳坪,高健,鄧海祥,等.復雜曲面零件在線檢測與誤差補償方法[J].機械工程學報,2012,48(23):143-151. CHEN Yue-ping,GAO Jian,DENG Hai-xiang,et al.On-line inspection and machining error compensation for complex surfaces[J].Journal of Mechanical Engineering,2012,48(23):143-151.(in Chinese)

[10] 張磊,張志勝.基于半參數回歸模型的制造過程加工誤差流建模與分析[J].機械工程學報,2013,49(15):181-185. ZHANG Lei,ZHANG Zhi-sheng.Stream of variation modeling and analysis for manufacturing processes based on a semi-parametric regression model[J].Journal of Mechanical Engineering,2013,49(15):181-185.(in Chinese)

[11] Suriano S,Wang H,Hu S J.Sequential monitoring of surface spatial variation in automotive machining processes based on high definition metrology[J].Journal of Manufacturing Systems,2012,31(1):8-14.

[12] 韓良,芮小健,鐘秉林.主軸回轉誤差信號的周期性試驗研究[J].制造技術與機床,1994(2):26-29. HAN Liang,RUI Xiao-jian,ZHONG Bing-lin.Experimental study on periodicity of spindle rotation error signal[J].Manufacturing Technology and Machine Tool,1994(2):26-29.(in Chinese)

[13] Tong L,Liu R,Soon V C,et al.Indeterminacy and indentifiability of blind indetification[J].IEEE Transactions on Circuits and System,1991,38(5):499-506.

[14] Hyv?rinen A,Oja E.Independent component analysis:algorithms and applications[J].Neural Networks,2000,13(4/5): 411-442.

[15] 趙忠華,楊曉梅.基于FastICA的語音盲源分離方法[J].四川大學學報:自然科學版,2015,52(4):830-834. ZHAO Zhong-hua,YANG Xiao-mei.Blind separation of sound signal by using FastICA[J].Journal of Sichuan University:Natural Science Edition,2015,52(4):830-834.(in English)

Independent Component Analysis Based on Machining Error Separation

ZHANG Fa-ping,WU Di,ZHANG Ti-guang,ZHANG Ling-yun,YANG Ji-bin
(School of Mechanical Engineering,Beijing Institute of Technology,Beijing 100081,China)

For the unsuccessful separation of the multiple systematic errors with similar scales by current machining error separation method,a new method to separate the machining errors is proposed based on independent component analysis.An error transfer model is built to describe the relationship between the systematic error caused by individual error source and the final systematic error measured from the machining surface.Then according to the theory of blind signal separation,an optimization model is used for the machining error separation where the negative entropy of the estimated error is used as the optimization objective function,and the fixed point algorithm is used as the optimization method.A method to determine the number of machining error sources is also given by means of the principal component analysis.A study case of a certain gyroscope surface is tested to verify the efficiency of the error separation method.The proposed method provides a new way for separation of machining errors and tracing of error sources.

manufacturing technology and equipment;error source;error propagation model;independent component analysis;principal component analysis;error separation

TH165+.3

A

1000-1093(2016)09-1692-08

10.3969/j.issn.1000-1093.2016.09.020

2015-11-18

國家自然科學基金項目(51275049);國防“973”計劃項目(613243)

張發平(1970—),男,副教授,碩士生導師。E-mail:zfpnew@163.com

猜你喜歡
測量區域方法
把握四個“三” 測量變簡單
滑動摩擦力的測量和計算
滑動摩擦力的測量與計算
關于四色猜想
分區域
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
測量
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
基于嚴重區域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
主站蜘蛛池模板: 久久精品国产免费观看频道| 久久黄色一级片| 久久鸭综合久久国产| 亚洲狼网站狼狼鲁亚洲下载| 色婷婷在线播放| 中文字幕波多野不卡一区| 国内精品久久久久久久久久影视| 国产95在线 | 原味小视频在线www国产| 最新加勒比隔壁人妻| 日本不卡在线| 制服丝袜一区二区三区在线| 午夜丁香婷婷| 国产日韩久久久久无码精品| 激情乱人伦| 日本在线免费网站| 亚洲精品波多野结衣| 一级一级一片免费| 超碰精品无码一区二区| 一区二区三区毛片无码| 女同国产精品一区二区| 亚洲天堂网站在线| 呦女精品网站| 亚洲精品久综合蜜| 久久美女精品| 人人澡人人爽欧美一区| 久久99热66这里只有精品一| 亚洲无码视频喷水| 色丁丁毛片在线观看| 蜜桃视频一区二区三区| 人妻精品全国免费视频| 日韩专区欧美| 99久久精品美女高潮喷水| 真实国产精品vr专区| AV色爱天堂网| 亚洲欧美精品在线| 亚洲天堂在线免费| 丰满人妻被猛烈进入无码| 无码一区二区波多野结衣播放搜索 | 国产综合日韩另类一区二区| 国产一区二区三区日韩精品| 亚洲中文字幕av无码区| 亚洲午夜片| 亚洲视频在线观看免费视频| 亚洲第一国产综合| 伊人色在线视频| 一级一级一片免费| 最新亚洲人成网站在线观看| 日韩成人在线视频| 69综合网| 51国产偷自视频区视频手机观看| 五月激情综合网| 亚洲色中色| 国产亚洲欧美在线中文bt天堂| 国产精品视频白浆免费视频| 丁香婷婷激情综合激情| 国产综合色在线视频播放线视| 狠狠色香婷婷久久亚洲精品| 久久久久久久久亚洲精品| 精品国产三级在线观看| 2021精品国产自在现线看| 日本高清免费不卡视频| 极品性荡少妇一区二区色欲| 免费一级无码在线网站| 国产在线八区| 综合色亚洲| 97精品国产高清久久久久蜜芽| 国产精品页| 欧美激情伊人| 欧美啪啪视频免码| 中文字幕一区二区人妻电影| 2020国产免费久久精品99| 2021亚洲精品不卡a| 免费国产高清视频| 亚洲日韩第九十九页| 国产91蝌蚪窝| 国内a级毛片| 在线精品视频成人网| 素人激情视频福利| 中国成人在线视频| 91美女视频在线观看| 亚洲av无码牛牛影视在线二区|