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

分層優(yōu)化PF在六維力傳感器下E型膜中的應(yīng)用

2014-06-06 10:46:47朱文超許德章
計算機工程 2014年9期
關(guān)鍵詞:信號系統(tǒng)

朱文超,許德章,方 濤

(安徽工程大學(xué)機械與汽車工程學(xué)院,安徽蕪湖241000)

分層優(yōu)化PF在六維力傳感器下E型膜中的應(yīng)用

朱文超,許德章,方 濤

(安徽工程大學(xué)機械與汽車工程學(xué)院,安徽蕪湖241000)

針對動載環(huán)境下,噪聲污染導(dǎo)致六維力傳感器測量精度急劇下降的問題,提出一種具有分層優(yōu)化步驟的改進粒子濾波算法。以雙E型彈性體六維力傳感器下E型膜為研究對象,根據(jù)正弦激勵力響應(yīng)和應(yīng)變的關(guān)系,建立非線性系統(tǒng)模型。在粒子濾波的框架下,將樣本集按權(quán)值的蛻化程度分層,引入野草繁殖算法,將最新的觀測信息融入高權(quán)值子集。基于Thompson-Taylor算法,通過聚合重采樣將高、低權(quán)值粒子隨機組合,產(chǎn)生中權(quán)值粒子集。將優(yōu)化后的粒子濾波算法在六維力傳感器動態(tài)測試系統(tǒng)中進行仿真研究,結(jié)果表明,該算法能以更小的估計誤差貼近真實后驗概率密度,在保持實時性的同時,有效地提高六維力傳感器的測量精度。

雙E型彈性體;六維力傳感器;下E型膜;粒子濾波;野草繁殖算法;Thompson-Taylor算法

1 概述

電阻應(yīng)變片式六維力傳感器的測量原理是彈性體表面應(yīng)變轉(zhuǎn)換為應(yīng)變片阻值的變化,通過惠斯頓電橋輸出便于測量的電信號。雙E型彈性體六維力傳感器采用組合梁結(jié)構(gòu)來探測空間6個方向的力和力矩。但由于電阻應(yīng)變片熱噪聲、放大電路以及電路周圍的電磁干擾源等原因,導(dǎo)致原始信號在傳輸、轉(zhuǎn)換、采集過程中不可避免地混入噪聲信號,嚴重影響了傳感器的測量精度和分辨率[1]。再加上電橋輸出信號弱,放大電路放大倍數(shù)高,有用信號很容易被噪聲信號淹沒,故有效地濾除隨機干擾和測量噪聲對傳感器課題的研究至關(guān)重要。

借助有用信號的規(guī)律性、噪聲信號隨機性特征,運用統(tǒng)計學(xué)最優(yōu)估計理論,濾除噪聲信號是現(xiàn)代信號處理技術(shù)研究的熱點之一[2]。國內(nèi)外專家試圖運用粒子濾波適用性廣泛、多模態(tài)處理能力強的特點,解決強非線性動態(tài)系統(tǒng)的狀態(tài)估計問題。文獻[3-4]針對移動機器人的運動控制及無人機組合導(dǎo)航問題,橫向?qū)Ρ攘朔蔷€性卡爾曼濾波與粒子濾波的性能。結(jié)果顯示,對于強非線性系統(tǒng)模型來說,當粒子數(shù)目達到閾值時,粒子濾波的估計精度要遠高于非線性卡爾曼濾波。但隨著時間的增加,粒子蛻化與貧化問題可能導(dǎo)致濾波精度大幅度降低,甚至出現(xiàn)發(fā)散現(xiàn)象。文獻[5-6]分別利用擴展卡爾曼粒子濾波算法與自適應(yīng)Unscented粒子算法構(gòu)造重要性密度函數(shù),提高了算法的估計精度,減小了權(quán)值蛻化速度,并將其應(yīng)用到編隊衛(wèi)星相對軌道的確定以及紅外弱小目標跟蹤問題中,獲得了較好的估計效果。然而上述2種算法均遵循傳統(tǒng)重采樣原則,經(jīng)過多次迭代后樣本多樣性的喪失仍無法避免。文獻[7]對重采樣后的每一個粒子施以MCMC移動,歷經(jīng)足夠的轉(zhuǎn)移步數(shù)后構(gòu)建的Markov鏈確保了重采樣后的粒子均能異步收斂到貝葉斯后驗。但收斂過程需要較長的燒穿時間,不能保證力傳感器能夠?qū)崟r、高精度完成作業(yè)。為了有效地解決粒子高效性和多樣性之間存在的矛盾,提高算法的實時性,本文在粒子濾波的基礎(chǔ)上,提出具有分層優(yōu)化步驟的粒子重采樣方案。根據(jù)權(quán)值蛻化程度,將樣本集分層;借鑒野草繁殖思想,將新息量測值融入高權(quán)值子集;基于Thompson-Taylor算法,通過聚合重采樣將高、低權(quán)值粒子隨機組合獲得中權(quán)值集,將低權(quán)值信息融入新粒子集。

2 系統(tǒng)模型

2.1 雙E型彈性體六維力傳感器的結(jié)構(gòu)

本文以雙E型彈性體六維力傳感器為研究對象,其結(jié)構(gòu)如圖1所示。傳感器主要由外傳力環(huán)1、內(nèi)傳力環(huán)2、薄矩形板3、上E型膜4、中心支柱5、下E型膜6和底座7組成。其中,傳感器底座與下轉(zhuǎn)接板剛性連接;上下E型膜與中心支柱相連;上E膜周圍設(shè)置內(nèi)傳力環(huán),并通過矩形薄板連接外傳力環(huán)。上E型膜用于檢測MX,MY方向力矩;下E型膜用來檢測FX,FY方向的力(切向力)和FZ方向的力(法向力),下E型膜受力分析如圖2所示;矩形梁用來檢測MZ方向的力矩。傳感器材料采用LY12,其彈性模量為72 GPa,密度為2 700 kg/m3,泊松比為0.33。

圖1 雙E型彈性體六維力傳感器結(jié)構(gòu)

圖2 下E型膜受力分析

2.2 非線性系統(tǒng)模型

當傳感器受法向力FZ作用時,外力垂直加載在上傳力環(huán)螺孔上,通過中心傳力環(huán)傳遞到下E型膜。其受力情況以及電阻應(yīng)變片布置位置如圖3所示。由于薄板所受橫向載荷繞Z軸對稱,則其彈性曲面也繞Z軸對稱。因此分析時可認為下E型膜只發(fā)生垂直方向的位移。定義下E型膜的邊界條件為外邊界固定、內(nèi)邊界自由,且在內(nèi)邊界受到正弦激勵力P(r,φ,t)=Ksinωtδ(r-ri)作用。

基于模態(tài)疊加理論,結(jié)合零初始條件,下E型膜對正弦激勵力的響應(yīng)可按極坐標系下的主振型函數(shù)與反映其運動規(guī)律的時間函數(shù)展開為如下無窮級數(shù)[8]:

其中,ωj為j階固有頻率;Wj(r,φ)為圓環(huán)薄板的主振型函數(shù);Jn,Yn為第一、二類貝塞爾函數(shù);In,Kn為第一、二類虛宗量貝塞爾函數(shù);固有頻率系數(shù):

其中,ω為激振頻率;為單位面積質(zhì)量。文獻[9]給出待定系數(shù)Dj1~Dj4,D′j1~D′j4的矩形齊次線性方程組,結(jié)合固有頻率系數(shù)λ,可獲得最終解答。

仔細審視法向力FZ主應(yīng)變片的分布位置后,不難發(fā)現(xiàn)應(yīng)變片的中心軸與薄板的徑向軸重合;圓環(huán)薄板受激勵力作用產(chǎn)生軸對稱彎曲變形,其表面的徑向應(yīng)變將轉(zhuǎn)化為電阻應(yīng)變片的軸向應(yīng)變。圓環(huán)薄板的徑向應(yīng)變信息是構(gòu)建系統(tǒng)量測方程的重要因素之一。

運用Kirchhoff彈性薄板小撓度彎曲理論,結(jié)合圓環(huán)薄板的撓度表達式,進而獲得圓環(huán)薄板的徑向應(yīng)變:

根據(jù)零初始條件下,下E型膜對正弦激勵力的響應(yīng)wr,φ(t)與徑向應(yīng)變εr,φ(t)之間關(guān)系,六維力傳感器下E型膜非線性離散時間系統(tǒng)狀態(tài)方程的一般形式可以表示為:

根據(jù)惠斯頓電橋原理可知,輸出電壓信號與信號檢測位置處4枚應(yīng)變片的阻值變化趨勢有關(guān)。量測方程由電阻應(yīng)變片的應(yīng)變信息與輸出電壓信號的關(guān)系所決定。非線性離散時間量測方程的一般形式可以表示為:

其中,V(k+1)是均值為0,方差為Rσ的高斯白噪聲序列。

3 粒子濾波

3.1 基本粒子濾波

在實際應(yīng)用過程中,粒子濾波在經(jīng)過若干次迭代后會出現(xiàn)蛻化現(xiàn)象。 局部重采樣[11-12](Partial Resampling,PR)可以有效地防止蛻化問題,提高粒子濾波的性能。其基本思想是將粒子按權(quán)值大小分為高、中、低3種權(quán)值的粒子,在進行重采樣時,低權(quán)值粒子被直接舍棄,中權(quán)值粒子權(quán)值不變,高權(quán)值粒子被多次復(fù)制,雖然局部重采樣只對部分粒子進行運算,提高了計算速度,但其只保留了高、中權(quán)值粒子信息,舍棄了低權(quán)值粒子信息。所以仍存在粒子貧化問題,不能保證粒子的多樣性。

3.2 分層優(yōu)化粒子濾波算法

Thompson-Taylor算法[13]是一種新穎的生成隨機樣本的方法,它既不過分依賴樣本的狀態(tài)空間分布,也不需要進行高斯近似,而是通過隨機數(shù)的方法對樣本集中某個樣本最近鄰的m個樣本進行平滑操作進而生成新的樣本。

Step 2 對高權(quán)值樣本集進行野草繁殖。

Step 3 通過Thompson-Taylor算法產(chǎn)生新的粒子集。

(2)產(chǎn)生均勻分布的隨機數(shù)集合:

其中,m為平滑參數(shù),其值為暫存域中需要優(yōu)化的粒子數(shù),即m=N-M。

3.3 復(fù)雜度計算

假設(shè)粒子與閾值比較的計算量為c1,高斯采樣的計算量為c2,計算適應(yīng)度函數(shù)的計算量為c3,高權(quán)值樣本集粒子排序的計算量為c4,暫存域粒子加權(quán)平均的計算量為c5,產(chǎn)生均勻分布隨機數(shù)的計算量為c6,產(chǎn)生中權(quán)值粒子的計算量為c7。以單位粒子的所需計算量為度量單位,對于LOPF算法,Step1劃分高低權(quán)值子集需要計算量Nc1;Step2高斯采樣產(chǎn)生子代粒子需要計算量Mc3+Lc2,優(yōu)化高權(quán)值子集需要計算量(M+L)(c3+c4);Step3暫存域中所有粒子加權(quán)平均需要進行Nc5次乘法,產(chǎn)生N個隨機數(shù)所需計算量為Nc6;最終產(chǎn)生中權(quán)值子集的產(chǎn)生所需進行Nc7次乘法。綜上所述,k時刻LOPF的計算復(fù)雜度可以表示為:

其中,L為產(chǎn)生的所有子代粒子數(shù)。

以上運算過程中忽略了部分實數(shù)加法的計算量。由式(13)可知,對于k時刻,分層優(yōu)化粒子濾波算法與基本粒子濾波以及局部重采樣粒子濾波算法[15]相比,計算復(fù)雜度并未明顯增加,均為O(N)。

4 六維力傳感器下E型膜濾波實例

4.1 六維力傳感器動態(tài)測試系統(tǒng)

如圖3所示,六維力傳感器動態(tài)測試系統(tǒng)主要由數(shù)據(jù)分析PC機、動態(tài)標定實驗臺、虛擬儀器NIPXI-1042Q、控制器PXI-8196、板卡 PXI-4461、放大電路和壓電陶瓷驅(qū)動模塊等部分組成。測試系統(tǒng)的工作原理如下:將雙E型彈性體六維力傳感器(圖4所示)固定在動態(tài)標定試驗臺上,選擇壓電陶瓷驅(qū)動模式,產(chǎn)生幅頻可調(diào)的正弦激勵力信號;通過信號轉(zhuǎn)接板加載在傳感器上,輸出的六路信號經(jīng)過放大電路增幅和緩沖后被虛擬儀器的數(shù)據(jù)采集模塊(板卡PXI-4661)高速采集;利用信號調(diào)理模塊(控制器PXI-8196)對采集數(shù)據(jù)進行實時處理,并將其導(dǎo)入PC機中的LabView平臺進行分析。

圖3 六維力傳感器動態(tài)測試系統(tǒng)

圖4 雙E型彈性體六維力傳感器

4.2 信號濾波

4.2.1 模型簡化

仔細分析下E型膜非線性離散時間系統(tǒng)狀態(tài)方程的一般形式,不難發(fā)現(xiàn),隨著狀態(tài)方程階數(shù)的增加,系統(tǒng)干擾項及控制項相應(yīng)的增加。過多的系統(tǒng)參數(shù)項會同時影響濾波器的運行時間和濾波精度,降低力傳感器的實時性和高效性。為了更好驗證改進粒子濾波算法的性能,將下E膜系統(tǒng)量測方程進行簡化,取其二階系統(tǒng)模型進行研究,即:

結(jié)合下E型膜前兩階主振型函數(shù)與電阻應(yīng)變片的貼片位置,可以獲得應(yīng)變片k+1時刻應(yīng)變之間的關(guān)系,即εr,x1(k+1)=εr,x4(k+1);εr,x2(k+1)=εr,x3(k+1),從而將下E型膜離散時間量測方程簡化為:

其中,R為應(yīng)變片的初始阻值;Ks為應(yīng)變片靈敏系數(shù);c為應(yīng)變函數(shù)εr,k(x3)與εr,k(x1)的幅值比。

考察下E型膜的結(jié)構(gòu)尺寸,其內(nèi)徑為50 mm,外徑為100 mm,板厚h為2 mm;4枚電阻應(yīng)變片質(zhì)心的位置參數(shù)分別為rx1,rx4=40 mm;rx2,rx3=20 mm;φx1,φx2=π/4;φx3,φx4=5π/4;取收斂系數(shù)A=0.05;a=0.05;B=0.5;應(yīng)變片靈敏系數(shù)Ks=1.4,初始阻值R=50 Ω;正弦激勵力的幅值K=20;激振頻率ω=200 Hz;根據(jù)下E型膜前兩階主振型相應(yīng)的參數(shù)(表1所示)求得幅值比c=3.900 9,固有頻率系數(shù)λ1=99.396 6,λ2=184.648 0,結(jié)合電阻應(yīng)變片X1的相關(guān)信息,獲得簡化的二階系統(tǒng)量測數(shù)學(xué)模型:

表1 固有頻率與振型振幅函數(shù)

舍棄式(16)的后4個小系數(shù)項,作為系統(tǒng)干擾進行處理;將剩下的余弦函數(shù)項視為系統(tǒng)控制項,則最終下E型膜的二階系統(tǒng)模型可以表達為:

其中,W(k)是均值為0,方差為0.001的高斯白噪聲序列。

4.2.2 信號濾波

設(shè)置采樣周期T為0.05 s,首先在空載狀態(tài)下,采集傳感器輸出的六路時域數(shù)據(jù),統(tǒng)計其規(guī)律,獲得測量噪聲V(k+1)方差值Rσ。接著將正弦激勵力P(t)=20sin(200·t),作為Fz方向的動載荷,加載在傳感器上(如圖3所示),獲得下E型膜動態(tài)輸出信號的時域數(shù)據(jù)集。取粒子數(shù)目N=200,高權(quán)值子集高斯采樣所產(chǎn)生的子代數(shù)為0~3。隨機從時域數(shù)據(jù)集中取出50個量測值,依次利用的局部重采樣粒子濾波算法(PRPF)、擴展卡爾曼粒子濾波(EKPF)以及基于分層優(yōu)化粒子濾波算法(LOPF)進行50個時刻的Monte Carlo仿真。狀態(tài)估計曲線與絕對誤差曲線如圖5、圖6所示,采樣時間間隔為0.05 s。EKPF的濾波精度與PRPF相比,略有提升。然而LOPF算法的系統(tǒng)狀態(tài)估計更接近于真實值(曲線的相似程度高),其濾波絕對誤差更小,估計精度更高。所以本文的LOPF算法在濾波效果上明顯優(yōu)于PRPF與EKPF算法。

圖5 3種濾波算法的系統(tǒng)狀態(tài)估計

圖6 3種濾波算法的估計誤差

通過11次實驗評價3種濾波器狀態(tài)估計均方誤差MSE的性能,并記錄3種算法的平均運行時間。如圖7與表2所示。LOPF算法的狀態(tài)估計均方誤差要明顯低于PRPF與EKPF算法,估計誤差的波動小,能更準確地估計上E型膜非線性系統(tǒng)的狀態(tài)。從表2可以看出,雖然EKPF與LOPF在迭代過程中同時融入了觀測值,但后者的計算耗時明顯少于前者。與PRPF相比,LOPF算法的平均運行時間并沒有明顯增加。故本文算法在提高六維力傳感器量測精度的同時,可以保持其實時性。

圖7 3種濾波器狀態(tài)估計均方誤差

表2 3種濾波方法的MSE與計算耗時比較

5 結(jié)束語

本文提出一種具有分層優(yōu)化步驟的改進重采樣粒子濾波算法。該算法在粒子濾波的基礎(chǔ)上,將野草繁殖思想和Thompson-Taylor算法有機地結(jié)合起來。根據(jù)正弦激勵力響應(yīng)和應(yīng)變的關(guān)系,得到下E型膜非線性系統(tǒng)模型。通過野草繁殖法將最新的觀測信息融入高權(quán)值子集,并保證其產(chǎn)生的子代仍分布在高似然區(qū)域附近。基于Thompson-Taylor算法,通過聚合重采樣將低權(quán)值粒子信息融入中權(quán)值粒子集,在保證粒子高效性的基礎(chǔ)上,實現(xiàn)了粒子集合的多樣性。理論分析和實驗結(jié)果表明,LOPF算法在解決粒子退化問題的同時,較好地避免了粒子匱乏,提高了樣本的多樣性。該算法的濾波精度及穩(wěn)定性明顯優(yōu)于PRPF與EKPF。由于其計算耗時較少,估計精度較高,因此可以保證六維力傳感器實時、高精度地完成作業(yè)。

[1] 梁康橋.特殊應(yīng)用的多維力/力矩傳感器研究與應(yīng)用[D].合肥:中國科學(xué)技術(shù)大學(xué),2009.

[2] Gao J B,Harris C J.Some Remarks on Kalman Filters for the Multi-sensor Fusion[J].Information Fusion, 2002,7(3):191-201.

[3] Rigatos G G.Nonlinear Kalman Filters and Particle Filters for Integrated Navigation of Unmanned Aerial Vehicles[J].Robotics and Autonomous Systems,2012, 60(2):978-995.

[4] Rigatos G G.Extended Kalman and Particle Filtering for Sensor Fusion in Motion Control of Mobile Robots[J]. Mathematics and Computers in Simulation,2010,81 (3):590-607.

[5] 張劍鋒,曾國強.擴展卡爾曼粒子濾波在編隊衛(wèi)星相對軌道確定中的應(yīng)用[J].航天控制,2010,28(4): 40-43.

[6] 康 莉,謝維信,黃敬雄.基于unscented粒子濾波的紅外弱小目標跟蹤[J].系統(tǒng)工程與電子技術(shù),2007,29 (1):1-4.

[7] 周 航,葉俊勇.運用聚類方法的分層采樣粒子濾波算法[J].計算機應(yīng)用,2013,33(1):69-71.

[8] 徐芝綸.彈性力學(xué)(下冊)[M].北京:知識產(chǎn)權(quán)出版社,2011.

[9] 汪志紅.電阻應(yīng)變片式六維力傳感器彈性體力學(xué)特性的研究[D].蕪湖:安徽工程大學(xué),2013.

[10] Arulampalam S,Maskell S R,Gordon N J.A Tutorial on ParticleFiltersforOn-line Nonlinear/Non-Gaussian Bayesian Tracking[J].IEEE Transactions on Signal Processing,2002,50(2):174-188.

[11] 胡士強,敬忠良.粒子濾波原理及原理[M].北京:科學(xué)出版社,2010.

[12] Ryan A,Hedrick J K.Particle Filter Based Informationtheoretic Active Sensing[J].Robotics and Autonomous Systems,2010,58(2):574-584.

[13] 常天慶,李 勇,劉忠仁,等.一種改進的重采樣的粒子濾波算法[J].計算機應(yīng)用研究,2013,30(3): 748-750.

[14] 陳 歡,周永權(quán),趙光偉,等.基于混沌序列的多種群入侵雜草算法[J].計算機應(yīng)用,2012,32(7): 1958-1961.

[15] Ghirmai T,Bugallo M F,Miguez J,et al.A Sequential MonteCarloMethodforAdaptiveBlindTimeing Estimation and Data Detection[J].IEEE Transactions on Signal Processing,2005,53(8):2855-2865.

編輯 顧逸斐

Application of Hierarchical Optimal Particle Filtering in Lower E-type Membrane of Six-axis Force Sensor

ZHU Wen-chao,XU De-zhang,FANG Tao
(College of Mechanical and Automotive Engineering,Anhui Polytechnic University,Wuhu 241000,China)

The measurement accuracy of the sensor which works on the environment of the dynamic load can be seriously affected by the pollution of noise signal.A new improved particle filtering which owns hierarchical optimal steps is proposed.This algorithm takes the rectangular thin plate of dual-E elastic body six-axis force sensor as the research object.The nonlinear state-space model based on the relationship between the response of sinusoidal excitation force and the strain is established.According to degenerate level,the sample sets can be divided into two parts.Based on the weeds breeding algorithm,the new measurement can be transferred to high likelihood region.Based on the Thompson-Taylor algorithm,the new particles set produced by random combinations of particles are achieved through polymerization resample of transferred particles.Simulation results indicate that the new algorithm can adjoin the real posterior probability density with smaller estimated error.It can effectively enhance the measurement accuracy of six-axis force sensor and maintain the real-time performance.

dual-E-type elastic body;six-axis force sensor;lower E-type membrane;particle filtering;weeds breeding algorithm;Thompson-Taylor algorithm

1000-3428(2014)09-0257-06

A

TP18

10.3969/j.issn.1000-3428.2014.09.052

國家自然科學(xué)基金資助項目(51175001);安徽省自然科學(xué)基金資助項目(11040606m144)。

朱文超(1989-),男,碩士研究生,主研方向:機器人信息感知;許德章,教授、博士;方 濤,碩士研究生。

2013-08-20

2013-10-15E-mail:zhuwenchao5102951@126.com

猜你喜歡
信號系統(tǒng)
Smartflower POP 一體式光伏系統(tǒng)
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
WJ-700無人機系統(tǒng)
ZC系列無人機遙感系統(tǒng)
北京測繪(2020年12期)2020-12-29 01:33:58
完形填空二則
基于PowerPC+FPGA顯示系統(tǒng)
半沸制皂系統(tǒng)(下)
孩子停止長個的信號
連通與提升系統(tǒng)的最后一塊拼圖 Audiolab 傲立 M-DAC mini
基于LabVIEW的力加載信號采集與PID控制
主站蜘蛛池模板: 亚洲日本精品一区二区| 97国产在线播放| 欧美色伊人| 久久公开视频| 亚洲天堂啪啪| 国产福利拍拍拍| 色香蕉影院| 日韩av电影一区二区三区四区| 国产成人91精品免费网址在线| 蝌蚪国产精品视频第一页| 精品综合久久久久久97超人| 亚洲成aⅴ人片在线影院八| 国产交换配偶在线视频| 国产午夜无码片在线观看网站| 91精品国产综合久久不国产大片| 亚洲人成人无码www| 国产极品嫩模在线观看91| 国产在线视频欧美亚综合| 亚洲欧美日韩色图| 亚洲一级毛片免费观看| 99热这里只有精品久久免费| 在线毛片免费| 成人无码区免费视频网站蜜臀| 六月婷婷精品视频在线观看| 亚洲色图在线观看| 另类重口100页在线播放| 国产精品极品美女自在线看免费一区二区| 国产成人一区在线播放| 亚欧成人无码AV在线播放| 国产福利微拍精品一区二区| 色天堂无毒不卡| 欧美一区二区三区不卡免费| 欧美福利在线观看| 亚洲伊人天堂| 国产精品黑色丝袜的老师| 日本免费高清一区| 国产精品天干天干在线观看| 午夜视频在线观看区二区| 五月丁香在线视频| 免费无码AV片在线观看国产| 国产午夜精品一区二区三| 久久综合九色综合97婷婷| 国产男女XX00免费观看| 在线看国产精品| 热思思久久免费视频| 欧美激情首页| 国产一级二级在线观看| 久久精品人人做人人综合试看| 无码一区二区三区视频在线播放| 国产主播在线观看| 国产精品视频导航| 久久国产成人精品国产成人亚洲| 精品丝袜美腿国产一区| 久久综合九九亚洲一区| 久久一级电影| 园内精品自拍视频在线播放| 成人免费黄色小视频| 国产精选自拍| www.亚洲一区二区三区| 亚洲精品日产AⅤ| 一级全黄毛片| 亚洲欧州色色免费AV| 高清无码手机在线观看| 在线观看的黄网| 国产成人综合网| 国产污视频在线观看| 欧美日韩激情在线| 亚洲精品人成网线在线 | 亚洲一区二区黄色| 福利一区在线| 国产精品福利社| 亚洲精品成人片在线观看| 欧美伦理一区| 青青草国产免费国产| 一级香蕉视频在线观看| 青青草国产免费国产| 国产噜噜噜| 亚洲无码高清一区二区| 久青草免费在线视频| 亚洲一区精品视频在线| 国产精品欧美激情| 日本高清视频在线www色|