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

        ?

        MEMS陀螺隨機誤差趨勢項的支持向量回歸機預(yù)測補償算法

        2016-12-23 01:44:46孫兆偉
        中國慣性技術(shù)學(xué)報 2016年5期
        關(guān)鍵詞:陀螺趨勢補償

        成 雨,葉 東,孫兆偉

        (哈爾濱工業(yè)大學(xué) 衛(wèi)星技術(shù)研究所,哈爾濱 150001)

        MEMS陀螺隨機誤差趨勢項的支持向量回歸機預(yù)測補償算法

        成 雨,葉 東,孫兆偉

        (哈爾濱工業(yè)大學(xué) 衛(wèi)星技術(shù)研究所,哈爾濱 150001)

        為了適應(yīng)低成本快速響應(yīng)衛(wèi)星的發(fā)展趨勢,越來越多的衛(wèi)星采用了MEMS陀螺系統(tǒng)與信息融合方法相結(jié)合的設(shè)計方案。為提高傳統(tǒng)的基于自回歸滑動平均(Auto-Regressive and Moving Average,ARMA)模型和支持度的信息融合方法的精度,抑制隨機誤差趨勢項對衛(wèi)星穩(wěn)定控制的不利影響,提出了一種基于支持向量回歸機(Support Vector Regression,SVR)的預(yù)測補償算法。該算法對各個MEMS陀螺輸出數(shù)據(jù)進行濾波,并提取相應(yīng)的隨機誤差趨勢項,通過相空間重構(gòu)獲得訓(xùn)練樣本并進行SVR建模,用以實時補償。然后使用低成本商用器件搭建了MEMS陀螺系統(tǒng),并在單軸氣浮轉(zhuǎn)臺上進行了實驗。實驗結(jié)果表明,預(yù)測補償算法使得MEMS陀螺系統(tǒng)輸出數(shù)據(jù)的方差降為原先的31.39%,融合精度得到了顯著的提高。

        微機電陀螺;時間序列;支持向量機;相空間重構(gòu)

        傳統(tǒng)衛(wèi)星的設(shè)計制造過程具有成本高、耗時長的特點,不利于衛(wèi)星技術(shù)在快速響應(yīng)任務(wù)中的應(yīng)用。商用現(xiàn)貨MEMS(Micro Electro Mechanical System)元器件體積小、重量輕、成本低,能夠有效提高衛(wèi)星的功能密度,縮短研制周期,為衛(wèi)星快速設(shè)計制造提供了思路。在這一思路的指導(dǎo)下,各國均有采用MEMS陀螺作為姿態(tài)確定系統(tǒng)的衛(wèi)星問世,如斯坦福大學(xué)的QuakeSat,東京工業(yè)大學(xué)的 Cute-1.7+APD2,哈爾濱工業(yè)大學(xué)的紫丁香2號等。為了適應(yīng)快速響應(yīng)任務(wù),加快衛(wèi)星設(shè)計制造進程,以MEMS陀螺構(gòu)成的衛(wèi)星姿態(tài)確定系統(tǒng)越來越受到關(guān)注[1-2]。

        由于工藝水平的限制,MEMS陀螺在精度和可靠性方面尚存在劣勢,無法替代高精度姿控系統(tǒng)中的傳統(tǒng)陀螺。為了滿足衛(wèi)星姿態(tài)確定系統(tǒng)的設(shè)計要求,可以采取兩種策略:一是大幅提高 MEMS陀螺測量精度,使之能夠匹敵傳統(tǒng)陀螺;二是通過設(shè)計體系結(jié)構(gòu),使用多個MEMS陀螺,構(gòu)建聯(lián)合姿態(tài)測量系統(tǒng),并采用信息融合算法提高測量系統(tǒng)的精度。由于前者會大幅增加整星成本,降低任務(wù)預(yù)期收益,因此后一種策略得到了學(xué)者們更多的研究[3]。

        一般地,MEMS陀螺信息融合算法首先對隨機誤差進行建模,然后應(yīng)用所建模型從量測數(shù)據(jù)中估計出真實狀態(tài)。時間序列模型將隨機誤差看成是平穩(wěn)線性的時間序列,常利用自回歸滑動平均算法(Auto-Regressive and Moving Average,ARMA),假設(shè) MEMS陀螺量測噪聲為歷史數(shù)據(jù)和歷史白噪聲的線性組合。這種方法建立的線性模型可以方便地應(yīng)用于 Kalman濾波器中,能夠有效地抑制隨機誤差[4]。但這種建立在線性平穩(wěn)假設(shè)之上的方法,對MEMS陀螺隨機誤差的建模必然會存在誤差。為了克服這一缺陷,文獻[5]和文獻[6]提出了基于小波變換的數(shù)據(jù)處理方法。該方法由Fourier變換發(fā)展而來,通過小波基和mallat算法,能夠同時在頻域和時域中對指定信號進行細分,可以獲得更多的細節(jié),從而有利于數(shù)據(jù)降噪[7]。但該方法為事后處理,實時性次于 Kalman濾波。另一方面,文獻[8]和文獻[9]使用神經(jīng)網(wǎng)絡(luò)替代ARMA模型。神經(jīng)網(wǎng)絡(luò)理論上可以逼近任意的非線性函數(shù)[10],因此所建立的模型能夠較好地適應(yīng)隨機誤差的非線性。但是,神經(jīng)網(wǎng)絡(luò)遵循經(jīng)驗風(fēng)險最?。‥mpirical Risk Minimization,ERM)的原則,當(dāng)樣本較少時,由于未考慮模型的置信風(fēng)險,存在過學(xué)習(xí)問題。

        為了保證MEMS陀螺信息融合方法的實時性,本文使用了基于 ARMA模型和支持度的 Kalman濾波器,同時針對線性模型精度的不足,提出了利用泛化性能更優(yōu)的支持向量回歸機(Support Vector Regression,SVR)對MEMS陀螺隨機誤差趨勢項進行預(yù)測補償。

        1 基于ARMA模型的Kalman濾波

        補償確定性誤差后,MEMS陀螺輸出值可以看作是由真實角速率kΩ和隨機誤差kε所組成的。對隨機誤差kε進行ARMA建模:

        式中:kε為k時刻的隨機誤差,ak為k時刻的白噪聲,iφ為自回歸參數(shù),jθ為滑動平均參數(shù)。為了便于觀測真實角速率kΩ,本文將其擴充在濾波器狀態(tài)變量中。在衛(wèi)星姿態(tài)保持階段和勻速機動過程中,可以認為wk為白噪聲序列。因此Kalman濾波的系統(tǒng)方程可以寫為:

        2 基于支持度的權(quán)值分配方法

        本文研究的MEMS陀螺系統(tǒng)由6個陀螺構(gòu)成,每個陀螺的輸出可以看作是真實角速率和噪聲的疊加,即

        每個 MEMS陀螺對同一輸入角速率進行冗余測量,其輸出結(jié)果可以相互仲裁。當(dāng)某一陀螺量測值明顯異于其它陀螺,則說明該陀螺不可靠,支持度低,因此在融合過程中給予較小的權(quán)值,抑制其對整個量測系統(tǒng)的不利影響。本文使用支持度函數(shù)衡量各個陀螺量測值之間的差異。

        對于理想的支持度函數(shù),當(dāng)兩個樣本點距離接近時,即相互支持度較高時,支持度函數(shù)值應(yīng)當(dāng)快速上升,反之則快速下降。因此使用 KMOD核函數(shù)[11]作為支持度函數(shù)。k時刻第 j個敏感器對第i個敏感器量測值的支持度為

        式中:k、γ和σ為可調(diào)參數(shù),用于調(diào)節(jié)支持度對量測值差異的敏感度。由式(5)構(gòu)造k時刻批量 MEMS陀螺的支持度矩陣

        和綜合支持度函數(shù)

        為了簡化計算,對式(7)進行歸一化處理,定義一致性度量:

        陀螺數(shù)據(jù)的可靠性由當(dāng)前時刻和歷史時刻的一致性度量共同決定,當(dāng)MEMS陀螺歷史時刻的一致性度量均較高,說明陀螺可靠,數(shù)據(jù)穩(wěn)定,應(yīng)給予較大的權(quán)值。

        k時刻第i個陀螺量測值的一致性均值為

        一致性方差為

        由式(9)和式(10)可以得到k時刻第i個陀螺量測值的權(quán)值為

        圖1 融合算法流程圖Fig.1 Flow chart of fusion algorithm

        式中:λ為可變參數(shù),可以調(diào)節(jié)權(quán)值分配策略。根據(jù)式(11),狀態(tài)變量估計值為

        綜上所述,融合算法流程如圖1所示。

        3 基于SVR的趨勢項預(yù)測補償算法

        從仿真結(jié)果中可以看出,經(jīng)過基于ARMA模型的Kalman濾波器處理后的MEMS陀螺的角速率估計誤差由趨勢項和高頻噪聲構(gòu)成,其中,對角度積分的影響較大,需要將其消除,以降低MEMS陀螺輸出數(shù)據(jù)的方差。由于趨勢項真實角速率ω和高頻噪聲混雜在一起,即

        當(dāng)星上高精度角速率敏感器正常工作時,即存在參考角速率ωref,其值可以看作是真實角速率,即,能夠用來標(biāo)定MEMS陀螺。擬合曲線與真實角速率之差即為趨勢項。當(dāng)參考角速率ωref有效時,對進行SVR建模。

        當(dāng)高精度角速率敏感器短時失效時,批量陀螺系統(tǒng)進入預(yù)測校正階段,以保證短時量測精度。

        圖2 趨勢項的預(yù)測補償算法流程圖Fig.2 Predictive compensation algorithm of trend term

        在衛(wèi)星進行穩(wěn)定控制時,可以認為真實角速率ω保持不變。擬合值與真實角速率做差,得到趨勢項的當(dāng)前值,與歷史值一同代入 SVR模型中,對下一時刻的趨勢項進行預(yù)測,從而補償了趨勢項的影響,降低了MEMS陀螺輸出數(shù)據(jù)的方差。

        綜上所述,基于SVR的趨勢項預(yù)測補償算法流程如圖 2 所示。

        3.1 SVR建模

        傳統(tǒng)的神經(jīng)網(wǎng)絡(luò)方法以 ERM為基礎(chǔ),當(dāng)樣本較少時,即使經(jīng)驗風(fēng)險很小,訓(xùn)練結(jié)果很好,學(xué)習(xí)機的實際風(fēng)險也可能很大,這就是 ERM 方法存在的過學(xué)習(xí)問題。而SRM(Structural Risk Minimization,SRM)綜合考慮了訓(xùn)練精度和置信范圍,使得學(xué)習(xí)機的實際風(fēng)險最小化,避免了過學(xué)習(xí)問題。

        給定訓(xùn)練樣本集:

        對于給出的訓(xùn)練樣本集,回歸問題便是尋找一個實值映射 f(·)來描述yi與xi的關(guān)系,從而對于任意一個輸入x,通過映射 f(·)便可得到輸出y。在分類問題中,要求分類超平面使得正反兩類樣本點的間隔最大。而在回歸問題中,要求所有樣本點到分類超平面(即回歸函數(shù)f(x))的距離盡量小。

        對于線性SVR,考慮下面的回歸函數(shù):

        訓(xùn)練樣本為

        假設(shè)在精度ε下,訓(xùn)練樣本存在擬合函數(shù) f( x),那么尋找最小ω的問題可以寫為凸優(yōu)化問題

        對于某些樣本點,函數(shù)f( x)在精度ε下無法擬合該樣本點,即線性不可分。松弛變量的引入可以解決這一問題,式(17)寫為

        引入Lagrange函數(shù)和對偶變量:

        選擇iα中的一個分量由此計算

        非線性 SVR通過非線性映射Φ(·)將輸入樣本映射到高維空間,然后在此空間中進行線性回歸。由此解決了輸入空間不可線性回歸的問題。

        但是這種提高系統(tǒng)維度的方法增加了計算復(fù)雜度,由此引入KMOD核函數(shù),由此計算

        對于常用的基于距離的核函數(shù),輸入空間內(nèi)的兩點距離越近,其相關(guān)性越強,因此需要在計算空間里盡可能將兩個點去相關(guān)性,同時還要盡可能地保留其它信息。設(shè)計的核函數(shù)需要擁有如下兩種特性:一是在零點附近快速減小,二是在遠離零點的區(qū)域緩慢下降。為此,N. Ayat[11]等人提出了緩和下降核函數(shù):

        式中:k為常數(shù),可以改變核函數(shù)在零點的值;γ和σ為可調(diào)參數(shù),用于改變核函數(shù)在零點附近的減小速度和核函數(shù)寬度。Ayat[11]等人已經(jīng)證明,在二分類問題中,KMOD核函數(shù)能夠以最少的支持向量獲得最大的分類間隔,具有良好的泛化能力。因此本文使用KMOD核函數(shù)構(gòu)建本章的SVR。

        由于核函數(shù)的固有屬性,可以很好地完成映射,并顯著降低計算量,由此優(yōu)化問題轉(zhuǎn)化為在式(20)的約束下,由

        求取

        所以回歸函數(shù) f( x)就可以寫為

        3.2 相空間重構(gòu)

        式中:m為嵌入維數(shù),τ為延遲時間。

        4 仿真結(jié)果及分析

        為了驗證上述方法的有效性,本節(jié)使用低成本商用器件搭建的 MEMS陀螺硬件系統(tǒng)(圖3),在單軸氣浮轉(zhuǎn)臺(圖4)上進行了驗證。氣浮轉(zhuǎn)臺以xPC為處理單元,俄制μFOG-6型光纖陀螺為測量單元,飛輪為執(zhí)行機構(gòu)。使用設(shè)計的批量MEMS陀螺系統(tǒng)測量轉(zhuǎn)臺的角速率,并在計算機上進行數(shù)據(jù)處理。

        圖3 批量MEMS陀螺系統(tǒng)Fig.3 System of multiple MEMS gyros

        圖4 單軸氣浮轉(zhuǎn)臺Fig.4 Single-axis air bearing turntable

        4.1 單軸氣浮轉(zhuǎn)臺穩(wěn)定控制

        當(dāng)氣浮轉(zhuǎn)臺穩(wěn)定控制時,在1~308 sample時,對光纖陀螺輸出數(shù)據(jù)取平均,看作是轉(zhuǎn)臺的真實角速率,此時角速度可以將此段過程看作是高精度角速度敏感器正常的情況,即存在參考角速度,且,在這段時間里建立角速率估計誤差趨勢項的 SVR模型。第309~700 sample時,認為高精度角速度敏感器短時失效,但可以認為參考角速率不變,即MEMS陀螺系統(tǒng)進入預(yù)測校正模式,通過之前建立的SVR模型預(yù)測, 并對陀螺Kalman濾波后的結(jié)果進行補償,消除趨勢項的影響。最后,對補償后的各陀螺數(shù)據(jù)進行了融合。下面以陀螺2、陀螺3、陀螺6為例,展示了相關(guān)數(shù)據(jù)處理的結(jié)果,如圖5~圖10所示。

        圖5 陀螺2的補償結(jié)果Fig.5 Compensation result of gyro 2#

        圖6 陀螺2的補償誤差和加權(quán)系數(shù)Fig.6 Compensation errors and weighted coefficients of gyro 2#

        圖7 陀螺3的補償結(jié)果Fig.7 Compensation result of gyro 3#

        圖8 陀螺3補償誤差和加權(quán)系數(shù)Fig.8 Compensation errors and weighted coefficients of gyro 3

        圖9 陀螺6的補償結(jié)果Fig.9 Compensation result of gyro 6#

        圖10 陀螺6補償誤差和加權(quán)系數(shù)Fig.10 Compensation errors and weighted coefficients of gyro 6#

        4.2 結(jié)果分析

        圖5、圖7和圖9為轉(zhuǎn)臺穩(wěn)定控制過程中基于SVR的趨勢項預(yù)測補償?shù)慕Y(jié)果,各圖上半部分為趨勢項的預(yù)測結(jié)果,下半部分為各子濾波器趨勢項補償前后的結(jié)果對比。圖6、圖 8和圖10的上半部分為補償誤差,下半部分為各陀螺在數(shù)據(jù)融合過程中的加權(quán)系數(shù)。

        圖5、圖7和圖9分別代表了基于SVR的趨勢項預(yù)測補償算法的三類結(jié)果:一是預(yù)測與實際在某些樣本點處存在偏差;二是預(yù)測精度較高;三是預(yù)測后期發(fā)散。出現(xiàn)這三種結(jié)果的原因在于前308個樣本點建模結(jié)果的不同泛化能力,由于采用的是基于模型的預(yù)測補償算法,因此不能保證每一個子濾波器的都可以達到較好的補償結(jié)果。

        對于這些由模型誤差引起的補償誤差,本文利用各子濾波器補償后結(jié)果的差異性,通過基于支持度融合算法予以仲裁和修正。由圖6、圖8和圖10可以看出,融合算法對各補償結(jié)果差異性的敏感性。由圖 8中上半部分可以看出,陀螺3的補償效果較好,其誤差的方差較小,可以認為在大部分時間內(nèi),其加權(quán)系數(shù)在六個陀螺之中是較高的。但是從圖6中可以看出,在第550~600個樣本點時,補償結(jié)果出現(xiàn)了較大的誤差,因而使得其加權(quán)系數(shù)在第550~600個樣本點處有所下降。同樣,從圖10也可以看出,由于補償結(jié)果的發(fā)散,其第550個樣本點之后的加權(quán)系數(shù)迅速下降。由此便可以看出,本文設(shè)計的融合算法對補償誤差的敏感性和對預(yù)測誤差引起的個別差異較大的補償結(jié)果抑制的有效性。從方差分析來看,未使用預(yù)測補償算法為 3.6753×10-5,使用預(yù)測補償算法為 1.1537×10-5??梢钥闯觯诤纤惴@著降低了角速率估計值的方差。

        5 結(jié) 論

        本文針對多個同類 MEMS陀螺構(gòu)成的衛(wèi)星姿態(tài)確定系統(tǒng),設(shè)計了基于ARMA模型的Kalman濾波器,根據(jù)支持度函數(shù)實時分配各個陀螺的加權(quán)系數(shù),從而得到最優(yōu)的角速率估計值。為了抑制陀螺隨機誤差趨勢項對衛(wèi)星穩(wěn)定控制的不利影響,本文提出了基于SVR的趨勢項預(yù)測補償算法。半實物仿真結(jié)果表明,經(jīng)過預(yù)測補償算法融合后的角速率估計值的方差是未使用時的31.39%,顯著提高了MEMS陀螺系統(tǒng)的測量精度。

        本文提出的趨勢項預(yù)測補償算法是針對衛(wèi)星穩(wěn)定控制過程的,適用于衛(wèi)星姿態(tài)保持階段或勻速機動過程。因此,后續(xù)工作將考慮在衛(wèi)星變速機動過程的信息融合算法。

        (References):

        [1]李博文, 姚丹亞. 低成本車載 MEMS慣導(dǎo)導(dǎo)航定位方法[J]. 中國慣性技術(shù)學(xué)報, 2014, 22(6): 719-723.Li B, Yao D. Low-cost MEMS IMU navigation positioning method for land vehicle[J]. Journal of Chinese Inertial Technology, 2014, 22(6): 719-723.

        [2]Abhilash M, Kumar S, Sandya S, et al. Implementation of the MEMS-based dual-axis sun sensor for nano-satellites[C]//Metrology for Aerospace (MetroAeroSpace). IEEE,2014: 190-195.

        [3]王新龍, 李娜. MEMS 陀螺隨機誤差的建模與分析[J].北京航空航天大學(xué)學(xué)報, 2012, 38(2): 170- 174.Wang X, Li N. Error modeling and analysis for random drift of MEMS gyroscopes[J]. Journal of Beijing University of Aeronautics and Astronautics, 2012, 38(2):170-174.

        [4]Chen M M, Gao G W. Research on MEMS gyroscope random error compensation algorithm based on ARMA model[C]//Applied Mechanics and Materials. 2014, 602-605: 891-894.

        [5]Shi Y S, Gao Z F. Study on MEMS gyro signal de-noising based on improved wavelet threshold method[J]. Applied Mechanics and Materials, 2013, 433-435: 1558-1562.

        [6]Yuan J, Yuan Y, Liu F, et al. An improved noise reduction algorithm based on wavelet transformation for MEMS gyroscope[J]. Frontiers of Optoelectronics, 2015, 8(4):413-418.

        [7]Farahani M A, Wylie M T V, Castillo-Guerra E, et al.Reduction in the number of averages required in BOTDA sensors using wavelet denoising techniques [J]. Journal of Lightwave Technology, 2012, 30(8): 1134-1142.

        [8]徐曉蘇, 周峰, 張濤,等. 遺傳算法優(yōu)化的神經(jīng)網(wǎng)絡(luò)在SINS/GPS中的應(yīng)用[J]. 中國慣性技術(shù)學(xué)報, 2015,23(3): 322-327.Xu X, Zhou F, Zhang T, et al. Application of neural network by genetic algorithm optimization in SINS/GPS[J]. Journal of Chinese Inertial Technology, 2015, 23(3):322-327.

        [9]Xiao W X, Zhang H, Ma X Q. Using BP neural network and kalman filter to signal processing of MEMS inertial sensors[J]. Sensors & Transducers, 2013, 156(9), 81- 88.

        [10]姜宇, 金晶, 張迎春. 基于異方差分析的多 MEMS陀螺隨機誤差補償方法[J]. 宇航學(xué)報, 2012, 33(6): 776-780.Jiang Y, Jin J, Zhang Y. Compensation method for random error signals of multiple MEMS gyroscopes based on heteroscedasticity analysis[J]. Journal of Astronautics, 2012, 33(6): 776- 780.

        [11]Ayat N E, Cheriet M, Remaki L, et al. KMOD-A new support vector machine kernel with moderate decreasing for pattern recognition. Application to digit image recognition[C]//Proceedings of the Sixth International Conference on Document Analysis and Recognition. IEEE,2001: 1215-1219.

        Predictive compensation algorithm for trend terms of MEMS gyro random errors based on support vector regression

        CHENG Yu, YE Dong, SUN Zhao-wei
        (Research Center of Satellite Technology, Harbin Institute of Technology, Harbin 150001, China)

        In order to keep up with the development trend of low cost and fast response satellite, more and more satellites adopt such design schema as combining the MEMS (Micro Electro Mechanical System) gyro system with an information fusion method. In order to improve the precision of traditional information fusion method based on ARMA (Auto-Regression and Moving Average) model and support degree, and to inhibit the random error trend-term’s adverse influences on satellite stability control, a predictive compensation algorithm based on SVR (Support Vector Regression) is proposed. By filtering the output data of each MEMS gyro, the algorithm extracts the corresponding random error trend term. With the phase space reconstruction, the training samples are obtained and then SVR model is built for the real time compensation.The MEMS gyro system is built with low cost commercial devices, and the experimental results are carried out on a single-axis air bearing table. Experiment results show that, by applying the predictive compensation algorithm, the output data’s variance of MEMS gyro system is reduced to 31.39% of the original one,showing that the fusion precision is significantly improved.

        MEMS gyro; time series; support vector machine; phase space reconstruction

        V441

        A

        1005-6734(2016)05-0600-07

        10.13695/j.cnki.12-1222/o3.2016.05.008

        2016-06-28;

        2016-09-18

        國家自然科學(xué)基金資助項目(61603115);中國博士后科學(xué)基金資助項目(2015M81455);微小型航天器技術(shù)重點學(xué)科實驗室開放基金資助(HIT.KLOF.MST.201501)

        成雨(1990—),男,博士研究生,從事微慣性系統(tǒng)與信息融合算法研究。E-mail: 14B918050@hit.edu.cn

        聯(lián) 系 人:葉東(1985—),男,博士/講師。E-mail: yed@hit.edu.cn

        猜你喜歡
        陀螺趨勢補償
        趨勢
        無功補償電容器的應(yīng)用
        山東冶金(2019年5期)2019-11-16 09:09:38
        做個紙陀螺
        玩陀螺
        陀螺轉(zhuǎn)轉(zhuǎn)轉(zhuǎn)
        軍事文摘(2018年24期)2018-12-26 00:58:18
        我最喜歡的陀螺
        快樂語文(2018年36期)2018-03-12 00:56:02
        初秋唇妝趨勢
        Coco薇(2017年9期)2017-09-07 21:23:49
        解讀補償心理
        SPINEXPO?2017春夏流行趨勢
        植物補償和超補償作用
        一本色道久久88精品综合| 青青自拍视频成人免费观看| 亚洲精品熟女av影院| 亚洲久悠悠色悠在线播放| 四川丰满妇女毛片四川话| 国语自产偷拍精品视频偷| 最新精品国偷自产在线婷婷| av资源吧首页在线观看| 国产精品国产传播国产三级| 亚洲国产精品无码av| 免费国产裸体美女视频全黄| 国产一区曰韩二区欧美三区| av网址大全在线播放| 在线播放国产自拍av| 久久精品国产自在天天线| 欧美日韩不卡视频合集| 国产三级精品美女三级| 亚洲国产av一区二区三| 虎白m粉嫩小在线播放| 东北老女人高潮大喊舒服死了| 97久久久久人妻精品专区| 中文字幕亚洲精品人妻| 深夜一区二区三区视频在线观看 | 人妻丝袜中文字幕久久| 久久99热只有频精品8国语| 一本色道无码道在线观看| 天堂国精产品2023年| 99久久婷婷国产亚洲终合精品 | 69精品人人人人| 亚洲尺码电影av久久| 成人试看120秒体验区| 日韩人妻无码一区二区三区久久99 | 全免费a级毛片| 亚洲黄片久久| 熟女肥臀白浆一区二区| 最近在线更新8中文字幕免费| 一本无码人妻在中文字幕| 免费看黄片视频在线观看| 国产又猛又黄又爽| 五月天丁香久久| 青青草视频在线播放81|