亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        一種基于橢球約束的加速度計(jì)粒子群優(yōu)化靜態(tài)修正方法

        2014-06-07 10:02:35張康寧于國(guó)輝李文超
        計(jì)量學(xué)報(bào) 2014年3期
        關(guān)鍵詞:橢球參數(shù)估計(jì)加速度計(jì)

        金 梅, 張康寧, 金 菊, 于國(guó)輝, 李文超

        (1.燕山大學(xué)河北省測(cè)試計(jì)量技術(shù)與儀器重點(diǎn)實(shí)驗(yàn)室,河北秦皇島 066004;

        2.河北工業(yè)大學(xué)土木工程學(xué)院,天津 300401; 3.秦皇島視聽(tīng)機(jī)械研究所,河北秦皇島 066000)

        一種基于橢球約束的加速度計(jì)粒子群優(yōu)化靜態(tài)修正方法

        金 梅1, 張康寧1, 金 菊2, 于國(guó)輝3, 李文超1

        (1.燕山大學(xué)河北省測(cè)試計(jì)量技術(shù)與儀器重點(diǎn)實(shí)驗(yàn)室,河北秦皇島 066004;

        2.河北工業(yè)大學(xué)土木工程學(xué)院,天津 300401; 3.秦皇島視聽(tīng)機(jī)械研究所,河北秦皇島 066000)

        研究了微機(jī)電系統(tǒng)加速度計(jì)基于橢球約束的靜態(tài)修正法。根據(jù)微機(jī)電系統(tǒng)加速度計(jì)輸出數(shù)學(xué)模型,采用標(biāo)準(zhǔn)粒子群優(yōu)化算法估計(jì)微機(jī)電系統(tǒng)加速度計(jì)的輸出數(shù)學(xué)模型中的刻度因數(shù)、零偏及正交誤差。實(shí)驗(yàn)結(jié)果表明,該方法有效剔除了采樣過(guò)程中的隨機(jī)誤差,粒子群優(yōu)化方法的引入,使得微機(jī)電系統(tǒng)加速度計(jì)輸出能夠準(zhǔn)確地反映其真實(shí)值。

        計(jì)量學(xué);微機(jī)電系統(tǒng);加速度計(jì);慣性運(yùn)動(dòng)跟蹤;橢球約束;靜態(tài)修正;粒子群優(yōu)化算法

        1 引 言

        在慣性運(yùn)動(dòng)跟蹤系統(tǒng)中,主要使用慣性微機(jī)電系統(tǒng)(MEMS:micro-electro-mechanical system)傳感器和磁力計(jì)實(shí)現(xiàn)角度估計(jì),因?yàn)镸EMS陀螺儀存在明顯的固有偏差和噪聲誤差,因此姿態(tài)角的累積誤差會(huì)隨時(shí)間而不斷增長(zhǎng)。為了修正陀螺儀的累積誤差,人們引入了加速度計(jì)和磁力計(jì)作為參考。由于傳感器的制造原因和使用環(huán)境,加速度計(jì)的輸出存在一定的誤差,限制了其使用效果,需要對(duì)其進(jìn)行修正。根據(jù)重力加速度為固定矢量,人們目前先后提出了多種修正方法,主要有3類:第一類是根據(jù)加速度的幅值為1g,典型的有工程中普遍使用的最大值最小值法[1],但該方法僅能估計(jì)標(biāo)度因子誤差和零位偏差,而且容易受到隨機(jī)誤差的干擾;第二類是根據(jù)已知姿態(tài)下加速度矢量值,進(jìn)行靜態(tài)標(biāo)定,主要有Seong-hoon[2]、秦永元[3]等提出了六位置測(cè)試法,該方法計(jì)算簡(jiǎn)單,修正精度高,但需要外部旋轉(zhuǎn)平臺(tái)提供精確的空間姿態(tài),限制了使用范圍。第三類就是和陀螺儀組成慣性測(cè)量單元進(jìn)行測(cè)試,北京航空航天大學(xué)的李建利[4]、孫宏偉[5]提出了組合標(biāo)定補(bǔ)償方法,可以同時(shí)校準(zhǔn)靜態(tài)和動(dòng)態(tài)誤差,數(shù)據(jù)利用率高,但其輸出數(shù)學(xué)模型復(fù)雜而且需要高精度旋轉(zhuǎn)平臺(tái)。

        在實(shí)際修正過(guò)程中,六位置測(cè)試法和組合標(biāo)定法,均需要依靠旋轉(zhuǎn)平臺(tái)提供參考數(shù)據(jù),這限制了其使用范圍;最大最小值法,通過(guò)加速度計(jì)勻速旋轉(zhuǎn)過(guò)程中的輸出重力場(chǎng)矢量,不需要外部參考設(shè)備,更適合普通用戶的使用,但是,該方法只能估計(jì)部分誤差參數(shù),包括零偏和量化誤差兩種,所以如何根據(jù)輸出的重力場(chǎng)矢量提高估計(jì)精度,擴(kuò)展估計(jì)的誤差參數(shù)種類,成為提高其估計(jì)精度的主要改進(jìn)手段。通過(guò)分析MEMS加速度計(jì)的誤差來(lái)源及其數(shù)學(xué)模型,在最大最小值法定基礎(chǔ)上,本文引入了基于橢球約束的最小二乘靜態(tài)修正方法,擴(kuò)展了修正的誤差參數(shù)種類。為了得到更為完整的加速度計(jì)的誤差修正模型,本文引入了(PSO:particle swarm optimization)算法,提出了基于橢球約束的加速度計(jì)PSO靜態(tài)修正方法,通過(guò)實(shí)驗(yàn)驗(yàn)證所提方法的有效性。

        2 MEMS加速度計(jì)的誤差模型

        實(shí)驗(yàn)中使用ADI公司生產(chǎn)的三軸加速度計(jì)ADXL345,該傳感器主要用于消費(fèi)類電子的微慣性傳感測(cè)量,最大可以測(cè)量±16g,測(cè)量精度為3.9g(每最低有效位),具體的現(xiàn)場(chǎng)采集模塊見(jiàn)圖1,三軸加速度計(jì)安裝誤差見(jiàn)圖2。本文設(shè)置其為13位A/D轉(zhuǎn)換,靈敏度為3.9g,量程為-2g~2g,故1g為255 Counts,通過(guò)MSP430F169的I2C總線進(jìn)行數(shù)據(jù)采集,采樣速度40 Hz。

        根據(jù)加速計(jì)的測(cè)量原理,MEMS加速度計(jì)的主要誤差包括零偏、量化、對(duì)準(zhǔn)、正交誤差和軸間干擾等。按照其誤差特點(diǎn),可得某MEMS陀螺捷聯(lián)慣導(dǎo)系統(tǒng)中MEMS加速度計(jì)的誤差模型方程為:

        圖1 加速度計(jì)現(xiàn)場(chǎng)采集模塊

        圖2 三軸加速度計(jì)安裝誤差

        式中:HS為傳感器載體坐標(biāo)下的加速度矢量,Gp為MEMS平臺(tái)載體坐標(biāo)下的加速度矢量,B為零偏,S為各軸的量化誤差,M為傳感器坐標(biāo)和MEMS平臺(tái)的正交誤差,?i,j為兩個(gè)坐標(biāo)系的對(duì)準(zhǔn)誤差角度。其中,

        2 基于橢球約束的最小二乘靜態(tài)修正

        因?yàn)橹亓铀俣葹榇笮『头较蚬潭ǖ目臻g矢量,當(dāng)重力加速度計(jì)圍繞中心點(diǎn)勻速旋轉(zhuǎn)時(shí),加速度矢量應(yīng)該是分布在中心點(diǎn)為(0,0,0)、半徑為1g的球面上。但因?yàn)闇y(cè)量誤差的存在,實(shí)際的加速度矢量分布在中心為(bx,by,bz)的橢球面上,因此加速度計(jì)的靜態(tài)修正轉(zhuǎn)換為橢球面約束下的參數(shù)估計(jì)[6]。最小二乘法就是根據(jù)重力加速度矢量的大小固定,將重力加速度真值作為參考值,進(jìn)行參數(shù)估計(jì),并建立目標(biāo)函數(shù),見(jiàn)式(2)。

        基于橢球約束的最小二乘修正方法。將式(1)代入式(2)可得目標(biāo)函數(shù),因此需要確定矩陣K和向量C中,共9個(gè)參數(shù),加速度計(jì)的誤差模型為:

        將式(3)代入式(2),可得式(4),因此誤差參數(shù)的估計(jì)就轉(zhuǎn)化為對(duì)式(4)中參數(shù)mi的最小二乘估計(jì),根據(jù)計(jì)算所得mi,即可得到相應(yīng)的誤差參數(shù)。

        3 基于橢球約束的PSO靜態(tài)修正

        在文獻(xiàn)[6]中,提出了更為全面的誤差模型。因?yàn)榧铀俣扔?jì)的各軸間存在相互干擾,有時(shí)可達(dá)到測(cè)量值的5%,因此將矩陣S修改為

        式(7)中,待估計(jì)的參數(shù)共有12個(gè),如果仍舊采用第2節(jié)提出的最小二乘法推導(dǎo)出和式(5)類似的方程組,很難直接計(jì)算出對(duì)應(yīng)的誤差參數(shù)。為了解決式(7)中參數(shù)估計(jì)的問(wèn)題,本文引入了PSO算法(粒子群優(yōu)化算法)。

        PSO算法[7]具有收斂速度快,計(jì)算模型簡(jiǎn)單等優(yōu)點(diǎn),隨著其誕生就在工程應(yīng)用中得到廣泛的推廣,因此也衍生了很多粒子群優(yōu)化的改進(jìn)算法。在分析了眾多改進(jìn)的粒子群算法后,標(biāo)準(zhǔn)粒子群算法分別在拓?fù)浣Y(jié)構(gòu)、慣性權(quán)重、粒子數(shù)目、粒子初始化和邊界條件等方面提出了相應(yīng)的參考標(biāo)準(zhǔn),并在大量的函數(shù)上驗(yàn)證了該算法的準(zhǔn)確性和穩(wěn)定性。

        將式(7)代入式(2),可以確定PSO估計(jì)算法的適應(yīng)度函數(shù),由此可知,待估計(jì)的參數(shù)共有12個(gè),則定義粒子維數(shù)為12。根據(jù)重力場(chǎng)的橢球分布特點(diǎn)分布和第2節(jié)所述,確定PSO算法的適應(yīng)度函數(shù)為:

        根據(jù)標(biāo)準(zhǔn)粒子群算法[8],定義D=12維的搜索空間,粒子數(shù)目為50,采用環(huán)形拓?fù)浣Y(jié)構(gòu)。第k次循環(huán)迭代中,粒子的位置和速度為Xi(k)=(Xi1(k),Xi2(k),Xi3(k),…,XiD(k))和Vi(k)=(Vi1(k),Vi2(k),Vi3(k),…,ViD(k)),相應(yīng)的各粒子個(gè)體最優(yōu)位置為Pi(k)=(Pi1(k),Pi2(k),Pi3(k),…,PiD(k)),個(gè)體最優(yōu)適應(yīng)度y=(y1,y2,y3,…,y50)和局部最優(yōu)位置gi(k)=(gi1(k),gi2(k),gi3(k),…,giD(k)),局部最優(yōu)適應(yīng)度依次為yg=(yg1,yg2,yg3,…,ygn)。各粒子的速度和位置更新如下式

        粒子速度更新中,右邊第1項(xiàng)是粒子飛行的慣性速度,來(lái)保證粒子處于飛行狀態(tài);第2項(xiàng)是自我學(xué)習(xí)部分,表示粒子從自身的飛行經(jīng)驗(yàn)學(xué)習(xí),向自己曾經(jīng)找到的最優(yōu)位置靠近,c1為其學(xué)習(xí)因子,r1為0~1的隨機(jī)數(shù);第3項(xiàng)是社會(huì)經(jīng)驗(yàn)部分,從鄰域中的其他粒子的飛行經(jīng)驗(yàn)中學(xué)習(xí),向鄰域粒子的最優(yōu)位置靠近,c2為其學(xué)習(xí)因子,r2為0~1的隨機(jī)數(shù)。在粒子的飛行過(guò)程中,限制了各粒子的飛行速度,同樣對(duì)粒子的搜索空間也進(jìn)行了相應(yīng)的限制,使各粒子分布在最優(yōu)位置的附近區(qū)域內(nèi),不僅可以提高搜索速度,而且可以提高準(zhǔn)確性。

        按照標(biāo)準(zhǔn)粒子群算法,其計(jì)算步驟:(1)設(shè)置相應(yīng)的參數(shù),粒子數(shù)目為N=50,各粒子為D=12;在可能的搜索區(qū)域內(nèi)對(duì)粒子群進(jìn)行隨機(jī)初始化,包括位置和速度。(2)計(jì)算各粒子的適應(yīng)度。按照索引號(hào)的環(huán)形拓?fù)浣Y(jié)構(gòu)內(nèi),取鄰域半徑為1,計(jì)算各粒子的個(gè)體最優(yōu)位置和局部最優(yōu)位置的適應(yīng)度。(3)在[-Vmax,Vmax]和[-Pmax,Pmax]內(nèi),對(duì)群體中各個(gè)粒子的速度和位置進(jìn)行更新。(4)更新粒子的個(gè)體最優(yōu)位置和其適應(yīng)度,并和鄰域范圍內(nèi)的局部最優(yōu)位置的適應(yīng)度比較后,更新其局部最優(yōu)位置和其適應(yīng)度。(5)如果迭代結(jié)束,則輸出最優(yōu)位置和其適應(yīng)度,否則返回步驟(3)。按照前面的方法,得到誤差修正模型結(jié)果:

        4 仿真實(shí)驗(yàn)及分析

        通過(guò)I2C總線,按照40 Hz進(jìn)行重力加速度值原始數(shù)據(jù)采樣,共采樣1936次,采樣點(diǎn)基本覆蓋了球面的所有點(diǎn),以保證參數(shù)估計(jì)的精度和準(zhǔn)確性。按照前文介紹的最小二乘法和PSO法對(duì)采樣所得數(shù)據(jù)進(jìn)行靜態(tài)修正處理,分別代入式(6)和式(11),2種方法修正后重力加速度幅值見(jiàn)圖3。

        圖3 重力加速度修正前后的幅值分布

        實(shí)驗(yàn)結(jié)果如表1所示,證明了橢球約束的有效性,但PSO法則能實(shí)現(xiàn)更好的誤差參數(shù)估計(jì),給出更為完善和全面的誤差模型,所以其參數(shù)估計(jì)效果更為理想。

        表1 重力加速度修正前后的幅值比較

        因此PSO靜態(tài)修正法明顯減小了重力加速度的測(cè)量誤差。造成這種區(qū)別的主要原因是最小二乘法對(duì)數(shù)學(xué)模型進(jìn)行了簡(jiǎn)化處理,只能對(duì)部分誤差進(jìn)行估計(jì)。而PSO靜態(tài)修正法,估計(jì)的誤差參數(shù),包括了加速度的主要誤差參數(shù),給出了整個(gè)空間內(nèi)的最優(yōu)估計(jì),更接近理想的誤差模型。

        在N-E-D坐標(biāo)系下,當(dāng)加速度計(jì)處于靜止時(shí),X、Y、Z三軸的測(cè)量值依次為X=-sin(θ)、Y=cos(θ)sin(φ)和Y=cos(θ)cos(φ),因此可知θ=asin(X)、φ=atan(Y/Z)(θ和φ分別為俯仰角和滾轉(zhuǎn)角)。為了驗(yàn)證文中靜態(tài)修正的效果,設(shè)計(jì)2組實(shí)驗(yàn)來(lái)測(cè)試角度測(cè)量的精度,分布為:實(shí)驗(yàn)1是θ=0,φ=[-180°,165°];實(shí)驗(yàn)2是φ=90°,θ=[-75°,75°];其角度間隔均為15°。

        采樣過(guò)程中,每個(gè)靜態(tài)位置間隔30 s采樣一次,共采樣15次,將各軸向平均值作為該姿態(tài)的重力加速度矢量,以消除隨機(jī)誤差的影響。為了避免歐拉角坐標(biāo)系中萬(wàn)向節(jié)死鎖的問(wèn)題,本文選擇θ的測(cè)量范圍為[-75°,75°]。

        根據(jù)式(6)和式(11)分別修正靜態(tài)位置下的加速度計(jì)各軸分量,見(jiàn)圖4,X、Y、Z各軸的標(biāo)準(zhǔn)差見(jiàn)表2。

        表2 標(biāo)準(zhǔn)姿態(tài)下各軸的標(biāo)準(zhǔn)差

        從修正結(jié)果可以看到,PSO的修正效果比最小二乘法有一定程度的改善,2組實(shí)驗(yàn)證明,在大部分情況下,PSO法的修正精度要好于最小二乘法。因?yàn)镻SO法和最小二乘法都是根據(jù)采樣的加速度計(jì)輸出矢量進(jìn)行統(tǒng)計(jì)分析,是針對(duì)重力場(chǎng)矢量在橢球面的分布進(jìn)行參數(shù)估計(jì),會(huì)出現(xiàn)某些姿態(tài)下PSO法的效果不如最小二乘法,如實(shí)驗(yàn)1中Z軸誤差分布圖。

        將各個(gè)靜態(tài)位置下測(cè)得的加速度計(jì)各軸分量,分別代入式(6)和式(11),完成靜態(tài)靜態(tài)修正。根據(jù)修正值,計(jì)算加速度計(jì)的俯仰角和滾轉(zhuǎn)角,并將其與真值進(jìn)行比較,得到各角誤差,如圖5所示。從圖5結(jié)果可知,經(jīng)過(guò)修正后,誤差角度小于5°,滿足了運(yùn)動(dòng)捕捉的使用要求;PSO的修正效果明顯提高了角度估計(jì)的精度,其角度誤差平均值為2°以內(nèi)。

        圖4 標(biāo)準(zhǔn)姿態(tài)下各軸的誤差

        圖5 標(biāo)準(zhǔn)姿態(tài)下各角的誤差分布

        5 結(jié) 論

        本文提出的基于橢球約束的加速度計(jì)PSO靜態(tài)修正方法針對(duì)原有的加速度計(jì)誤差模型,引入了新的智能算法,能夠有效地估計(jì)加速度計(jì)的主要誤差參數(shù),包括軸間誤差,保證了誤差模型的完整,給出了整個(gè)空間內(nèi)的最優(yōu)估計(jì),同時(shí)將該方法與最小二乘靜態(tài)修正法進(jìn)行了對(duì)比,實(shí)驗(yàn)結(jié)果證明,PSO靜態(tài)修正方法能實(shí)現(xiàn)更好的誤差參數(shù)估計(jì),給出更為完善和全面的誤差模型,所以其參數(shù)估計(jì)效果更為理想。在最大最小值的基礎(chǔ)上,本文提出了基于橢球約束的種修正方法,分別是最小二乘靜態(tài)修正法和PSO靜態(tài)修正法,2種方法的誤差模型和目標(biāo)函數(shù)基本一致,都是根據(jù)重力場(chǎng)矢量的橢球分布,是根據(jù)待估計(jì)矢量的強(qiáng)度完成參數(shù)估計(jì),所以2種方法都不能估計(jì)加速度計(jì)坐標(biāo)和載體坐標(biāo)間的對(duì)準(zhǔn)誤差,所以有必要對(duì)對(duì)準(zhǔn)誤差專門的誤差修正。另一方面,參數(shù)估計(jì)的準(zhǔn)確度和采樣點(diǎn)的分布情況有直接關(guān)系,在橢球面上,應(yīng)該實(shí)現(xiàn)采樣點(diǎn)的均勻分布,所以下一步的研究將集中在如何確定采樣過(guò)程中,加速度計(jì)的旋轉(zhuǎn)操作流程,減少采樣點(diǎn)的數(shù)目,并提高估計(jì)精度。

        [1] 黃徽,劉建業(yè),李榮冰.六余度配置MEMS-IMU誤差標(biāo)定[J].應(yīng)用科學(xué)學(xué)報(bào),2009,(27):106-110.

        [2] Won SP,Golnaraghi F.A Triaxial Accelerometer Calibration Method Using a Mathematical Model[J].Instrumentation and Measurement,2010,59(8):2144-2153

        [3] 宋麗君,秦永元.MEMS加速度計(jì)的六位置測(cè)試法[J].測(cè)控技術(shù),2009,28(7):11-17.

        [4] 李建利,房建成,盛蔚.微小型捷聯(lián)慣性測(cè)量單元標(biāo)定及補(bǔ)償方法[J].宇航學(xué)報(bào),2008,29(3):947-951.

        [5] 孫宏偉,房建成,盛蔚.一種基于MEMS的微慣性測(cè)量單元標(biāo)定補(bǔ)償方法[J].北京航空航天大學(xué)學(xué)報(bào),2008,34(4):439-442.

        [6] Skog I,Handel P.Calibration of A MEMS Inertial Measurement Unit[C]//XVII IMEKO WORLD CONGRESS,Rio de Janeiro,Brazil,2006:17-22.

        [7] 黃席樾,向長(zhǎng)城,殷禮勝.現(xiàn)代智能算法理論及應(yīng)用[M].北京:科學(xué)出版社,2009:192-203.

        [8] Bratton D,Kennedy J.Defing a Standard for Particle Swarm Optimization[C]//Swarm Intelligence Symposium,Harbin,Heilongjiang,China,2007:120-127.

        A Particle Swarm Optim ization Static Calibration Method of Micro-electro-mechanical System Accelerometer Based on Ellipsoidal Restriction

        JIN Mei1, ZHANG Kang-ning1, JIN Ju2, YU Guo-hui3, LIWen-chao1
        (1.Measurement Technology and Instrumentation Key Lab of Hebei Province,Yanshan University,Qinhuangdao,Hebei066004,China;2.School of Civil Engineering,Hebei University of Technology,Tianjin 300401,China;3.Audio-Visual Machinery Research Institute,Qinhuangdao,Hebei066000,China

        The particle swarm optimization static calibrationmethod ofmicro-electro-mechanical system accelerometer based on ellipsoidal restriction is studied.According to the output mathematical model of the micro-electro-mechanical system accelerometer,the standard particle swarm optimization algorithm is taken to estimate scale factors,bias,and installation errors in the outputmathematical model.The experimental results show that the method can eliminate the random errors in the sampling process effectively.With the introduction of particle swarm optimization algorithm,the outputs of themicro-electro-mechanical system accelerometer can reflect their true values.

        Metrology;Micro-electro-mechanical system;Accelerometer;Inertial motion tracking;Ellipsoidal restriction;Static calibration;Particle swarm optimization algorithm

        TB934

        A

        1000-1158(2014)03-0230-06

        10.3969/j.issn.1000-1158.2014.03.08

        2013-06-08;

        2013-08-09

        國(guó)家自然科學(xué)基金(61077071)

        金梅(1977-),女,黑龍江明水人,燕山大學(xué),博士,主要研究方向?yàn)橐曈X(jué)伺服、智能控制及智能信號(hào)處理。meijin297@126.com

        猜你喜歡
        橢球參數(shù)估計(jì)加速度計(jì)
        基于新型DFrFT的LFM信號(hào)參數(shù)估計(jì)算法
        獨(dú)立坐標(biāo)系橢球變換與坐標(biāo)換算
        基于加速度計(jì)的起重機(jī)制動(dòng)下滑量測(cè)量
        橢球槽宏程序編制及其Vericut仿真
        智能制造(2021年4期)2021-11-04 08:54:44
        橢球精加工軌跡及程序設(shè)計(jì)
        基于外定界橢球集員估計(jì)的純方位目標(biāo)跟蹤
        基于遺傳算法的加速度計(jì)免轉(zhuǎn)臺(tái)標(biāo)定方法
        Logistic回歸模型的幾乎無(wú)偏兩參數(shù)估計(jì)
        常溫下硅微諧振加速度計(jì)零偏穩(wěn)定性的提高
        基于向前方程的平穩(wěn)分布參數(shù)估計(jì)
        国产伦精品免编号公布| 亚洲一区有码在线观看| 久久精品熟女亚洲av香蕉| 成人中文乱幕日产无线码| 四房播播在线电影| 亚洲人成网站久久久综合| 白白色发布视频在线播放| 国产精品美女久久久网av| 国产大陆亚洲精品国产| 色播中文字幕在线视频| 亚洲av一区二区在线| 免费国产黄网站在线观看视频| 精品福利视频一区二区三区| 丰满熟妇人妻无码区| 亚洲视频在线观看第一页| 天天摸夜夜摸摸到高潮| 岛国AV一区二区三区在线观看| 精品中文字幕手机在线 | 美女下蹲露大唇无遮挡| 欧美内射深喉中文字幕| 国产成人亚洲综合无码精品| 久久精品国产亚洲av专区| 日本一区二区三区爆乳| 99久久国产露脸精品竹菊传媒| 成人免费无码视频在线网站| 人妖在线一区二区三区| 狠狠噜天天噜日日噜无码| 国产午夜在线观看视频播放| 日本高清一区二区三区色| 四虎成人精品在永久免费| 成人做爰高潮尖叫声免费观看| 无码一区二区丝袜| 日本超级老熟女影音播放| 久久精品国产久精国产| 国产精品青草久久久久婷婷| 亚洲av男人的天堂在线| 中文字幕网伦射乱中文| 欧美第五页| 日本超骚少妇熟妇视频| 在线中文字幕乱码英文字幕正常 | 精品视频一区二区杨幂|