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

基于遞歸特征分析的NURBS曲線插補算法

2015-10-28 09:46:30江本赤田曉青
中國機械工程 2015年1期

韓 江 江本赤 夏 鏈 田曉青

合肥工業(yè)大學(xué),合肥,230009

基于遞歸特征分析的NURBS曲線插補算法

韓江江本赤夏鏈田曉青

合肥工業(yè)大學(xué),合肥,230009

在分析NURBS曲線公式中基函數(shù)遞歸特征的基礎(chǔ)上,提出了一種新的NURBS曲線插補算法。該算法通過前一個參數(shù)增量值和前一段插補弦長來確定下一個參數(shù)增量值,進而實現(xiàn)對參數(shù)的密化。仿真實驗結(jié)果表明,該算法在改善弦長誤差和減小計算量方面優(yōu)于Taylor展開插補算法。

NURBS曲線;遞歸特征;迭代算法;插補

0 引言

在CAD/CAM系統(tǒng)中,一般采用參數(shù)形式來表示復(fù)雜的曲線和曲面。NURBS具備很強的形狀控制能力,于1991年被ISO確定為定義產(chǎn)品形狀的唯一數(shù)學(xué)方法[1]。然而,傳統(tǒng)的CNC系統(tǒng)只能進行直線和圓弧插補,需要用大量微直線段和圓弧逼近自由曲線,不僅需要大量的存儲空間,而且需要很長的計算時間,這顯然不符合高速加工的實時性要求。更重要的一點是用多段直線和圓弧表示曲線,段間的不連續(xù)破壞了曲線或曲面的光滑度,同時會導(dǎo)致速度、加速度的不連續(xù)[2]。因此,學(xué)者們對多種NURBS曲線直接插補技術(shù)進行了研究[3-6]。通過將NURBS曲線參數(shù)傳遞給數(shù)控系統(tǒng),然后由CNC系統(tǒng)生成軌跡點,實現(xiàn)對NURBS曲線的直接插補。NURBS曲線插補就是在已知當前參數(shù)的基礎(chǔ)上,通過合理確定參數(shù)增量值來求取下一個參數(shù)值,然后從參數(shù)空間映射到軌跡空間,下一個插補點就是將新的參數(shù)值代入數(shù)學(xué)模型而得到的。因而,NURBS曲線插補的關(guān)鍵在于確定下一個插補點的參數(shù),實現(xiàn)參數(shù)的密化。

Bedi等[7]首先提出將曲線參數(shù)在定義域內(nèi)等分的方法,該方法會造成插補速度的波動。Shpitalni等[8]針對參數(shù)曲線提出了Taylor展開近似方法的插補算法。現(xiàn)有文獻中報道的NURBS曲線插補算法大多是基于Taylor展開公式逼近而實現(xiàn)插補點的參數(shù)密化的[9-14]。由于Taylor一階逼近會產(chǎn)生較大的截斷誤差,羅福源等[10]對基于Taylor一階展開插補法的平穩(wěn)性進行了研究。Shiuh等[11]在Taylor展開近似計算的基礎(chǔ)上考慮了曲線的幾何特性,提出了閉環(huán)控制NURBS插補算法。Farouki等[12]和Yeh等[13]分別采用Taylor二階逼近和改變速度的方法來減小插補誤差。Liu等[14]提出了考慮動力學(xué)的NURBS插補算法。Feng等[15]提出了考慮軸加速度限制的NURBS實時自適應(yīng)插補算法。上述考慮加減速和機床動力學(xué)的NURBS曲線插補算法也采用了Taylor二階逼近。無論采用哪一種算法,Taylor一階逼近都會產(chǎn)生較高的插補弦長誤差,而Taylor二階逼近將大大增加運算量。

事實上,通過考察NURBS曲線的定義公式不難發(fā)現(xiàn),其中的B樣條基函數(shù)是由遞推公式定義的,存在遞歸特性。該特性的幾何意義體現(xiàn)在,當相鄰的兩個參數(shù)增量值接近時,通過插補所得到的兩個弦長值亦較接近。因此,一方面,可以通過預(yù)估的參數(shù)增量值Δu初步獲取所需弦長;另一方面,盡管參數(shù)與弦長之間不存在精確的解析關(guān)系,但仍可根據(jù)初步得到的弦長值反過來對參數(shù)增量值Δu進行適當?shù)奈⒄{(diào),經(jīng)過少量的迭代計算得到理想的參數(shù)增量值。

本文提出了一種基于NURBS自身遞歸特征的NURBS曲線簡易插補算法。該算法通過前一個參數(shù)增量值和前一個弦長值來獲取下一個插補點的參數(shù)增量值。為得到理想的弦長值,在算法中設(shè)置了過渡參數(shù)增量值,最終的有效參數(shù)增量值則是利用上述遞歸特征更新過渡增量值而得到的。對一條典型NURBS曲線進行的插補仿真實驗驗證了所提算法的可行性。

1 NURBS曲線的定義

一條p次NURBS曲線定義為

(1)

a≤u≤b

其中,Pi為控制頂點;wi為權(quán)因子;Ni,p(u)是定義在非周期節(jié)點矢量U上的p次B樣條基函數(shù),它是在區(qū)間[ui,ui+p+1)上的非零多項式,其數(shù)值為

(2)

(3)

節(jié)點矢量U中,a和b的重復(fù)度都為p+1,除非特別聲明,通常取a=0,b=1。

式(2)和式(3)為B樣條基函數(shù)定義公式,不難看出其具有遞歸特征。以定義在節(jié)點矢量U=(0,0,0,0,0.25,0.5,0.75,1,1,1,1)上的3次基函數(shù)為例,考察基函數(shù)的變化情況,如圖1所示。為便于觀察,圖1中只給出了其部分基函數(shù)。可見,各基函數(shù)是連續(xù)的,不存在突變。

圖1 某3次NURBS曲線的部分基函數(shù)

當Δu=0.005時,取相鄰兩基函數(shù)的差值ΔNi,p(u)=Ni,p(u)-Ni-1,p(u),可以更清楚地觀察基函數(shù)的變化情況,如圖2所示。可見,基函數(shù)的遞歸特征決定其變化的平滑性,因而利用遞歸特征求取均勻的插補弦長存在理論上的可行性。

圖2 相鄰基函數(shù)的變化量

2 基于遞歸特征的NURBS曲線插補算法

NURBS曲線插補的關(guān)鍵在于求取下一個插補點對應(yīng)的參數(shù),而軌跡插補算法就是實現(xiàn)參數(shù)u在節(jié)點矢量區(qū)間上的密化。針對Taylor展開算法的不足,為減小運算量和改善插補弦長誤差,本文提出了基于NURBS自身遞歸特征的插補算法。

2.1求取初始參數(shù)增量Δu1的算法

設(shè)節(jié)點區(qū)間為0≤u≤1,即umin=0,umax=1。按照遞歸特征求取參數(shù)增量值,需要前一個參數(shù)增量值和前一段插補弦長,但是用這種方法顯然無法求出第一個參數(shù)增量值Δu1。考慮到當參數(shù)u從0增大到1時,就完成了對整個NURBS曲線的插補,即當參數(shù)的變化量為Δut=umax-umin時,弧長的變化量為NURBS曲線的總弧長Lt。設(shè)期望弦長值為Ld,可嘗試用Δut∶Lt=Δu1∶Ld這個比例關(guān)系來初步求出Δu1。但是由于計算弧長Lt比較復(fù)雜,經(jīng)分析NURBS曲線的有關(guān)參數(shù)信息,本文擬利用控制點之間的距離總長Lk(即控制多邊形的總邊長)代替Lt來初步求出Δu1,所產(chǎn)生的誤差可以通過較簡單的迭代計算來消除。Lk的計算公式為

(4)

式中,xP(i)、yP(i)為控制頂點Pi的坐標值。

由于Lt與Lk之間存在較大誤差,將此處所求出的參數(shù)增量值稱為過渡參數(shù)增量值,記作Δu1,1,其計算公式為

Δu1,1=ΔutLd/Lk

(5)

同理,記過渡參數(shù)值u2,1=u1+Δu1,1,利用NURBS公式求出過渡插補點C(u2,1)之后,可得到首個過渡弦長L1,1=|C(u2,1)-C(u1)|。將L1,1與期望弦長Ld之間相對誤差率記為ε1,其表達式為

ε1=(Ld-L1,1)/Ld

(6)

記給定的許用相對誤差為[εL],當ε1>[εL]時,則需求出第2次過渡參數(shù)增量值,按照下式進行1次迭代計算(為了明確算法的計算量,將使用NURBS公式的次數(shù)稱作迭代次數(shù)):

Δu1,1=Δu1,1Ld/L1,1

(7)

同理,可得到過渡參數(shù)值u2,2=u1+Δu1,2和過渡弦長L1,2=|C(u2,2)-C(u1)|,并求出相對誤差ε2。則求取初始參數(shù)增量值的通式為

Δu1,k=Δu1,k-1Ld/L1,k-1k≥2

(8)

連續(xù)應(yīng)用式(8),直到弦長相對誤差滿足要求為止,而此時的過渡值即為最終的有效值。即當εk≤[εL]時,Δu1,k=Δu1,u2,k=u2,L1,k=L1。計算Δu1所需迭代次數(shù)的確定方法將在實驗仿真部分詳細說明。

2.2求取其余參數(shù)增量Δui(i>1)的算法

按照上述對遞歸特征的分析,當兩個參數(shù)增量值大小接近時,所得到的兩弦長值也比較接近,即當Δui-1≈Δui時,Li-1≈Li,反之亦然。因此,可以通過期望弦長值反求Δui。上文已經(jīng)求出了Δu1和u2,其余的插補點參數(shù)為ui+1=ui+Δui(i>1),其中的參數(shù)增量Δui,1的計算公式為

Δui,1=Δui-1Ld/Li-1

(9)

類似于式(8),通過m-1次迭代求取其余參數(shù)增量值的通式為

Δui,m=Δui,m-1Ld/Li,m-1m≥2

(10)

當εm≤[εL]時,Δui,m=Δui。特別地,當ui>umax時,取ui=umax。計算Δui(i>1)所需迭代次數(shù)的確定方法也將在實驗仿真部分說明。

3 仿真實驗與分析

為了驗證本文算法的可行性并評估其性能,現(xiàn)以節(jié)點矢量U=(0,0,0,0,0.1,0.2,0.3,0.4,0.5,0.6,0.8,0.9,1,1,1,1),控制頂點的坐標分別為(2.0,8)mm,(3.0,4.8)mm,(4.0,3.0)mm,(5.0,2.0)mm,(6.2,2.5)mm,(7.5,8.8)mm,(9.0,8.5)mm,(10.5,7.0)mm,(11.0,5.0)mm,(13.5,0)mm,(16.0,5.5)mm,(18.0,7.0)mm,wi=1,p=3的NURBS曲線為例進行仿真實驗。如果將期望弦長設(shè)為恒定值,則利用本文算法所獲得的插補點分布情況如圖3所示。

圖3 期望弦長設(shè)為恒定值時的插補點分布示意圖

3.1迭代次數(shù)的確定

3.1.1確定計算Δu1所需迭代次數(shù)

第一個參數(shù)增量Δu1的求解,由于沒有前一段弦長值可供參考,本文參照了NURBS曲線控制多邊形的邊長值,雖然會產(chǎn)生較大的弦長相對誤差,但可通過后續(xù)的迭代計算來改善。為了控制誤差范圍,設(shè)置許用弦長誤差率[εL]作為迭代計算的判斷依據(jù),即當前一次計算所得弦長的誤差率超過[εL]時,就需進行下一次迭代計算。下面的仿真實驗是按照2.1節(jié)所述的計算方法開展的。

實驗結(jié)果表明,當期望弦長Ld=0.1mm時,未經(jīng)迭代計算(僅根據(jù)控制多邊形的邊長值按比例求解)所得的過渡參數(shù)增量Δu1,1=2.9067×10-3,弦長誤差率ε1=186.46%;經(jīng)1次迭代計算后的Δu1,2=1.0147×10-3,弦長誤差率ε2=1.33%;迭代2次所得的Δu1,3=1.0013×10-3,弦長誤差率則降至ε3=0.01%。可見,對于該NURBS曲線,若設(shè)[εL]=2%作為控制條件,則只需進行1次迭代;若設(shè)[εL]=1%,則需進行2次迭代。

當期望弦長Ld=1μm時,未經(jīng)迭代計算所得的Δu1,1=2.91×10-5,ε1=192.29%;經(jīng)1次迭代計算后的Δu1,2=9.9445×10-6,而弦長誤差率降至ε2=0.01%;迭代2次所得的Δu1,3=9.9432×10-6,弦長誤差率則驟降至ε3=10-8。顯然,當期望弦長值為μm級時,若以[εL]=1%作為控制條件,則只需進行1次迭代,而即使[εL]=0.1%,迭代2次已經(jīng)足夠。

作為首個參數(shù)增量,Δu1的精確度將影響后續(xù)參數(shù)增量Δui(i>1)計算的迭代效率。因此,必要時可以考慮減小[εL]值以適當增加迭代計算的次數(shù),但是3次迭代已經(jīng)能獲得足夠高的精度。

3.1.2確定計算Δui(i>1)所需迭代次數(shù)

對照上文計算Δu1的方法,來探究求解Δui(i>1)所需的迭代次數(shù),仍分別取Ld=0.1mm和Ld=1μm進行實驗。為了表征弦長值的波動情況,特引入標準差σ,其計算公式如下:

(11)

圖4 Ld=0.1 mm時的插補弦長變化情況

圖5 Ld=1 μm時的插補弦長變化情況

綜合圖4和圖5可見,在同等弦長誤差率限制條件下,期望弦長Ld越小,求解參數(shù)增量Δui所需的迭代次數(shù)越少,當Ld為μm級時,未經(jīng)迭代計算而根據(jù)遞歸特征直接求解即可得到非常高的計算精度。

3.2誤差分析

曲線插補的誤差一般分為徑向誤差和弦高誤差。對于本文所提出的直接插補算法,由于插補點都是通過NURBS公式直接計算出來的,即所有插補點都落在曲線上,不存在徑向誤差。其弦高誤差Ei計算公式為

(12)

式中,ri為第i個插補點C(ui)處的曲率半徑;Li為第i段插補弦長。

由式(12)可以看出,當弦長值Li很小時,Ei也很小。以本文所用的NURBS曲線為例,其曲率分布情況如圖6所示,當u=0.224時,曲率半徑取得最小值rmin=0.559mm,此時的弦高誤差Ei將達到最大值Emax。當Li=Ld=0.1mm時,Emax=2μm;當Li=Ld=1μm時,Emax=0.2236nm。

圖6 NURBS曲線曲率分布圖

可見,當期望弦長為μm級時,本文算法產(chǎn)生的弦高誤差與目前數(shù)控機床的分辨率相比,完全可以忽略。

3.3與泰勒展開方法對比

為了進一步驗證本算法的性能,從運算量和插補弦長誤差兩個方面,將其與常用的Taylor展開系列算法進行對比。

3.3.1插補運算量的比較

Taylor二階展開插補算法常用的公式為

而本算法中當?shù)螖?shù)為2次以上(m≥2)時所用的公式為

Δui,m=Δui,m-1Ld/Li,m-1

顯然本算法極大地減少了運算量,更易于滿足高速插補的實時性要求。

3.3.2插補弦長誤差的比較

仍以圖3所示的NURBS曲線為例,設(shè)置期望弦長Ld=0.1 mm,比較兩種算法插補運算的弦長誤差情況。

圖7 一階Taylor展開與一次迭代計算效果對比

圖8 二階Taylor展開與二次迭代計算效果對比

結(jié)果表明,與Taylor展開系列算法相比,本算法不僅運算量小,而且弦長誤差也更小,總體上優(yōu)于基于Taylor展開的插補算法。

4 結(jié)語

在分析NURBS自身遞歸特征的基礎(chǔ)上,本文提出了一種新的NURBS曲線插補算法。該算法通過前一段插補弦長和前一個參數(shù)增量的關(guān)系來獲取下一個插補點對應(yīng)的參數(shù)。實驗結(jié)果表明:①通過利用NURBS公式進行少量的迭代計算可顯著改善插補的弦長誤差,并可根據(jù)所設(shè)置的精度自動調(diào)整迭代次數(shù);②在運算量相當?shù)臈l件下,期望弦長越短,插補效率越高,當期望弦長值為μm級時,未經(jīng)迭代計算而根據(jù)遞歸特征直接求解亦可得到極高的插補精度,而且弦高誤差也足夠小。

由于該算法避免了復(fù)雜的NURBS曲線求導(dǎo)運算,代之以簡單的四則運算,故其計算量遠小于基于Taylor展開的系列算法,改善了插補的實時性。而且該算法插補弦長的標準差也小于Taylor展開算法,可以更有效地抑制插補速度的波動。

[1]王田苗,曹宇男,陳友東,等.基于de Boor算法的NURBS曲線插補和自適應(yīng)速度控制研究[J].中國機械工程,2007,18(21):2608-2613.

Wang Tianmiao, Cao Yunan, Chen Youdong, et al. NURBS Interpolation and Feedrate Adaptive Control Based on de Boor Algorithm[J].China Mechanical Engineering, 2007, 18(21):2608-2613.

[2]葉佩青,趙慎良.微小直線段的連續(xù)插補控制算法研究[J].中國機械工程,2004,15(15):1354-1356.

Ye Peiqing, Zhao Shenliang. Study on Control Algorithm for Micro-line Continuous Interpolation[J].China Mechanical Engineering, 2004, 15(15):1354-1356.

[3]趙國勇,徐志祥,趙福令.高速高精度數(shù)控加工中NURBS曲線插補的研究[J].中國機械工程,2006,17(3):291-294.

Zhao Guoyong, Xu Zhixiang, Zhao Fuling.Study on NURBS Curve Interpolator in the High Speed and High Accuracy CNC Machining[J].China Mechanical Engineering, 2006, 17(3):291-294.

[4]游有鵬,王珉,朱劍英.NURBS曲線高速高精度加工的插補控制[J].計算機輔助設(shè)計與圖形學(xué)學(xué)報, 2001, 13(10): 943-947.

You Youpeng, Wang Min, Zhu Jianying. An Interpolator for NURBS Curve Machining with High Speed and High Accuracy[J].Journal of Computer-Aided Design & Computer Graphics, 2001, 13(10): 943-947.

[5]Nam S H, Yang M Y.A Study on a Generalized Parametric Interpolator with Real-time Jerk-limited Acceleration[J]. Computer-Aided Design, 2004, 36(4): 27-36.

[6]彭芳瑜,何瑩,羅忠誠,等.NURBS曲線機床動力學(xué)特性自適應(yīng)直接插補[J].華中科技大學(xué)學(xué)報(自然科學(xué)版),2005,33(7):80-83.

Peng Fangyu, He Ying, Luo Zhongcheng, et al. NURBS Curve Interpolation Algorithm with the Adaptation to Machine Tool Kinetic Character[J]. J. Huazhong Univ. of Sci. & Tech. (Nature Science Edition), 2005,33(7):80-83.

[7]Bedi S, Ali I, Quan N. Advanced Techniques for CNC Machines[J].Journal of Engineering for Industry, 1993, 115(3): 329-336.

[8]Shpitalni M, Koren Y, Lo C C.Real-time Curve Interpolators[J]. Computer-Aided Design, 1994, 26(11): 832-838.

[9]劉宇,戴麗,劉杰,等. 泰勒展開NURBS曲線插補算法[J].東北大學(xué)學(xué)報(自然科學(xué)版),2009,30(1):117-120.

Liu Yu, Dai Li, Liu Jie, et al. Study on NURBS Interpolator with Taylor Expansion[J].Journal of Northeastern University (Natural Science),2009,30(1):117-120.

[10]羅福源,游有鵬,尹娟.NURBS曲線泰勒展開插補法的平穩(wěn)性與改進研究[J].中國機械工程, 2012, 23(4):383-388.

Luo Fuyuan, You Youpeng, Yin Juan. Research on Stability and Improvement of Taylor-Expansion-Based Approach for NURBS Curve Interpolation[J].China Mechanical Engineering, 2012, 23(4):383-388.[11]Shiuh S, Hsu P. Adaptive-feedrate Interpolation for Parametric Curves with a Confined Chord Error[J].Computer-Aided Design, 2002, 34(3):229-237.

[12]Farouki R T, Tsai Y F.Exact Taylor Series Coefficients for Variable-feedrate CNC Curve Interpolators[J].Computer-Aided Design, 2001, 33(2): 155-165.

[13]Yeh S, Hsu P.The Speed-controlled Interpolator for Machining Parametric Curves[J]. Computer-Aided Design, 1999, 31(1): 349-357.

[14]Liu X B, Ahmad F, Yamazaki K, et al. Adaptive Interpolation Scheme for NURBS Curves with the Integration of Machining Dynamics[J].International Journal of Machine Tools and Manufacture, 2005, 45(4/5):433-444.

[15]Feng J C, Li Y H,Wang Y H,et al. Design of a Real-time Adaptive NURBS Interpolator with Axis Acceleration Limit[J].The International Journal of Advanced Manufacturing Technology,2010,48(1/4): 227-224.

(編輯陳勇)

NURBS Curve Interpolation Algorithm Based on Recursive Feature Analysis

Han JiangJiang BenchiXia LianTian Xiaoqing

Hefei University of Technology,Hefei,230009

Based on the recursive feature analysis of the base function in NURBS formula, a novel interpolation algorithm was put forwards. The method obtained the next parameter increment by a previous parameter increment and a previous chord length, and then the parameter interpolation was achieved.Simulation results show that, the proposed interpolation algorithm has obvious advantages in improving chord length errors and lessening computation load compared with those using Taylor’s equation.

NURBS curve;recursive feature;iterative algorithm;interpolation

2013-09-12

國家自然科學(xué)基金資助項目(51275147)

TP273< class="emphasis_italic">DOI

:10.3969/j.issn.1004-132X.2015.01.019

韓江,男,1963年生。合肥工業(yè)大學(xué)機械與汽車工程學(xué)院教授、博士研究生導(dǎo)師。主要研究方向為先進制造技術(shù)、數(shù)控技術(shù)與數(shù)控裝備。發(fā)表論文100余篇。江本赤,男,1979年生。合肥工業(yè)大學(xué)機械與汽車工程學(xué)院博士研究生。夏鏈,女,1964年生。合肥工業(yè)大學(xué)機械與汽車工程學(xué)院教授。田曉青,女,1987年生。合肥工業(yè)大學(xué)機械與汽車工程學(xué)院講師。

主站蜘蛛池模板: 看你懂的巨臀中文字幕一区二区| 2020久久国产综合精品swag| 免费观看成人久久网免费观看| 少妇露出福利视频| 国产成人超碰无码| 粉嫩国产白浆在线观看| 国产麻豆福利av在线播放| 午夜不卡视频| 91无码视频在线观看| 日韩精品一区二区深田咏美| 国产电话自拍伊人| 91精品国产福利| 国产精品爽爽va在线无码观看| 99久久精品国产自免费| 91欧美亚洲国产五月天| 国产成人高清精品免费软件 | 亚洲高清资源| 欧美久久网| 欧美国产综合色视频| 精品人妻一区无码视频| 久久久久无码精品| 国产成人艳妇AA视频在线| 国产办公室秘书无码精品| av在线人妻熟妇| 国产日本欧美在线观看| 久久久久亚洲Av片无码观看| 亚洲欧洲日本在线| 九九九国产| 福利在线一区| 国产精品视频导航| 成人午夜视频网站| 亚洲国语自产一区第二页| 免费Aⅴ片在线观看蜜芽Tⅴ| 四虎精品国产永久在线观看| 国产成人综合久久| 免费人成网站在线观看欧美| 国产成人亚洲综合a∨婷婷| 亚洲精品777| 在线毛片免费| 精品一区二区三区视频免费观看| 国产精品嫩草影院av| 亚洲狠狠婷婷综合久久久久| 国产精品太粉嫩高中在线观看| 天天色天天综合网| 一级毛片免费观看不卡视频| 制服丝袜无码每日更新| 人妻少妇久久久久久97人妻| 国产乱子伦视频在线播放| 亚洲欧美在线看片AI| 精品91视频| 国产无码网站在线观看| 亚洲综合中文字幕国产精品欧美| 毛片最新网址| 免费国产在线精品一区| 国产99欧美精品久久精品久久| 欧美一道本| 亚洲国产欧美目韩成人综合| 日韩乱码免费一区二区三区| 激情网址在线观看| 97超爽成人免费视频在线播放| 国产成熟女人性满足视频| 亚洲视频二| 免费一看一级毛片| 丰满人妻久久中文字幕| 亚洲精品少妇熟女| 亚洲国产第一区二区香蕉| 毛片免费观看视频| 69综合网| 日韩免费无码人妻系列| 国产sm重味一区二区三区| 亚洲欧洲日本在线| 在线亚洲精品福利网址导航| 久久国产免费观看| 国产精品无码在线看| 波多野结衣久久高清免费| 亚洲AⅤ综合在线欧美一区| 一级毛片在线播放免费观看 | 日韩国产一区二区三区无码| 亚洲欧美日韩成人在线| 国产精品理论片| 国产一级在线播放| 一区二区三区四区在线|