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

一種新的MEMS陀螺儀信號去噪方法

2017-11-08 02:07:12孫田川劉潔瑜
關(guān)鍵詞:信號方法

孫田川,劉潔瑜

(火箭軍工程大學(xué) 控制工程系,西安710025)

一種新的MEMS陀螺儀信號去噪方法

孫田川,劉潔瑜

(火箭軍工程大學(xué) 控制工程系,西安710025)

為提高M(jìn)EMS陀螺儀輸出信號的去噪效果,將稀疏分解(sparse decomposition)與提升小波變換(lifting wavelet transform)相結(jié)合,提出了一種新的信號去噪方法.首先,建立MEMS陀螺帶噪信號的誤差模型,并利用小波提升正變換計算帶噪信號的非稀疏的小波系數(shù);然后,利用稀疏分解理論恢復(fù)小波系數(shù)的稀疏性;最后,再通過小波提升反變換重構(gòu)信號,從而達(dá)到去噪的目的.考慮到梯度投影(gradient projection)算法具有全局最優(yōu)解,運(yùn)算效率更高,將梯度投影思想引入恢復(fù)信號稀疏性的過程中,提出了基于梯度投影的稀疏分解算法,給出了利用梯度投影算法進(jìn)行信號系數(shù)分解的具體步驟,大大簡化了計算復(fù)雜度,同時提升了算法的穩(wěn)定性.為驗證所提方法的性能,進(jìn)行了MEMS陀螺信號去噪的靜態(tài)實驗和跑車實驗.實驗結(jié)果表明,此種方法在動靜態(tài)條件下都可以有效地去除MEMS陀螺儀輸出信號中的噪聲,尤其是在靜態(tài)條件下的去噪效果要優(yōu)于小波閾值濾波方法. 同時采用的梯度投影算法相比于正交匹配追蹤算法和基追蹤算法具有更高的運(yùn)算效率.

MEMS陀螺儀;信號去噪;稀疏分解;提升小波變換;梯度投影;凸優(yōu)化

MEMS陀螺儀具有成本低、體積小、功耗低、抗沖擊能力強(qiáng)等優(yōu)點,近年來得到了快速的發(fā)展[1],然而由于存在較大的漂移,限制了它在很多領(lǐng)域的應(yīng)用[2-3].所以,使用合適的去噪方法對MEMS陀螺儀的輸出信號進(jìn)行去噪處理,可以有效降低漂移的影響,從而提高精度[4].

MEMS陀螺儀通常采用的去噪方法包括卡爾曼濾波、時間序列分析、小波去噪等[5].小波在時頻和細(xì)節(jié)描述的優(yōu)越性拓展了它在MEMS陀螺去噪領(lǐng)域的應(yīng)用[6].小波去噪方法可以分為:小波閾值法、模極大值法和平移不變量法等,而閾值法應(yīng)用則比較多.這種方法的主要依據(jù)是:在小波域中,有用信號可以用較大的小波系數(shù)表示,噪聲信號則相反,用較小的小波系數(shù)表示,那么通過設(shè)置一個合理的閾值,就可以將較小的小波系數(shù)濾除,保留較大的小波系數(shù),從而實現(xiàn)對噪聲信號的抑制[7].這種方法的難點就在閾值的選取上,由于閾值選擇的固定性限制了小波方法的自適應(yīng)性,所以小波去噪對多變信號的處理效果就較差[8].

為了克服小波閾值濾波的缺陷,稀疏分解去噪方法被提出并受到廣泛關(guān)注.稀疏分解去噪主要原理是干凈信號通過提升小波正變換后,得到的各個尺度下的小波系數(shù)都較小,也就是說小波系數(shù)是稀疏的;而當(dāng)信號中含有噪聲時,得到的小波系數(shù)中含有的較大幅值的小波系數(shù)就會變多,換言之,小波稀疏的稀疏性隨之降低[9].所以,利用稀疏分解理論對非稀疏小波系數(shù)進(jìn)行稀疏性恢復(fù),可以達(dá)到去噪的目的.

恢復(fù)小波系數(shù)稀疏性過程中所采用的算法對于恢復(fù)后的系數(shù)的稀疏度有很大的影響,文獻(xiàn)[10-14]采用貪婪算法中的正交匹配追蹤算法(orthogonal matching pursuit, OMP),取得較好的效果.文獻(xiàn)[15]采用凸松弛方法中的基追蹤算法(basis pursuit, BP),也驗證了稀疏分解去噪方法的優(yōu)勢.本文采用的梯度投影算法(gradient projection, GP)是凸松弛算法中的一種,相對于貪婪算法,此種算法具有全局最優(yōu)解,而相對于基追蹤算法運(yùn)算效率更高.

1 MEMS陀螺儀數(shù)學(xué)模型

低精度的MEMS陀螺儀輸出模型可描述如下:

W(t)=ω(t)+ω0.

式中:W(t)為輸出信號,ω0為漂移,°/s.

MEMS陀螺儀的漂移通常包括:常值漂移、周期漂移和白噪聲.由于MEMS陀螺儀性能較傳統(tǒng)陀螺儀性能普遍不高,通常在連續(xù)工作時間較短的場合得到廣泛應(yīng)用,那么對其短時漂移特性的研究就顯得更為重要,根據(jù)其短時漂移強(qiáng)周期性和規(guī)律性的特點,可以將其漂移的模型描述如下:

ω0=εb+Ωdsin(2πfd+θ0)+n(t).

式中:εb為陀螺儀的零偏,短時間內(nèi)可以近似為一個常數(shù);Ωd為周期分量的幅值;fd為周期分量的頻率;θ0則為初始相位;n(t)為零均值高斯白噪聲.

2 稀疏分解去噪與提升小波變換基本理論

2.1 稀疏分解去噪的數(shù)學(xué)模型

本文先對一維信號進(jìn)行研究,將信號和噪聲的關(guān)系表示如下:

y=s+n.

式中:s為有用信號;n為高斯白噪聲;y為帶噪信號.y、s、n分別為一維實向量,且長度是M,假設(shè)n~N(0,σ2).

本文可以用下式來表示對y的稀疏分解

y=s+n=Ax+n,

(1)

假設(shè)n=Az,則式(1)可以表示為

y=s+n=A(x+z).

值得注意的是,有用信號s在A上的分解系數(shù)x是稀疏的,或者說系數(shù)中的非零系數(shù)只有K個,而噪聲信號n在A上的分解系數(shù)z是非稀疏的,而且非零的稀疏分布在整個域.本文要做的就是首先得到y(tǒng)在A上的分解系數(shù),然后將其中幅值較大的K個系數(shù)選取出來,最后利用這K個系數(shù)所對應(yīng)的原子的線性組合來近似本文需要的有用信號s.

(2)

式中,spa(·)為問題的稀疏解.

從式(2)可以看出想要(P0,δ)滿足條件很困難,所以本文要解決這個問題,就需要尋求其他的表達(dá)方式,也就是基于l0范數(shù)的(P0,δ)問題的變體.

2.2 提升小波變換

提升小波變換是基于提升方案的小波變換.由于提升小波只在時域上進(jìn)行分析,而且它的反變換具有簡單直接和高效等優(yōu)點,所以本文采用這種變換方式.正向變換的過程包括分解、預(yù)測和更新.反向變換則相反,正向變換和反向變換的步驟如圖1所示.

圖1 提升小波變換分解與重建

本文采取的是5/3提升小波變換,其中正反變換的公式如下:

3 基于梯度投影的稀疏分解算法(gradient projection for sparse reconstruction, GPSR)

對于稀疏分解去噪與提升小波變換基本理論所述的(P0,δ)問題,考慮凸優(yōu)化方法,用l1范數(shù)代替l0范數(shù),即有

(3)

求解(P1,δ)問題,可以將其與基追蹤去噪(basis pursuit denoising, BPDN)算法關(guān)聯(lián)起來,并修改式(3)如下

首先,(P1,δ)問題可得變式如下:

(4)

式中,A為一個k×n矩陣,而τ為非負(fù)參數(shù),y=Ax+n.

為求解此問題,通過引入向量u,v對x做如下替換:

x=u-v,u≥0,v≥0.

那么,問題(4)可被轉(zhuǎn)換成一個帶約束的二次線性規(guī)劃問題如下:

(5)

問題(5)又可以寫成如下的標(biāo)準(zhǔn)形式:

(6)

其中:

(7)

表面看來,轉(zhuǎn)化后的問題維數(shù)是問題(4)的2倍,而實際上這種維數(shù)上的增大對于問題求解影響很小.由式(7)可以得到

(8)

那么本文要計算Bz只需要計算向量u-v和ATA.

由式(6)得其投影為

F(z)=c+Bz,

亦可由式(8)得zTBz如下

本文將描述使用GP算法求解該問題具體步驟:

在GP算法中通過迭代來確定k.

首先,選擇標(biāo)量參數(shù)α(k)>0,并令

w(k)=(z(k)-α(k)F(z(k)))+.

然后,選擇參數(shù)λ(k)∈[0,1],令

z(k+1)=z(k)+λ(k)(w(k)-z(k)).

再定義向量g(k)如下:

并選擇α0的初始值為

).

具體計算過程如下:

(9)

那么,GP算法解決此問題的具體過程如下.

Step2通過式(9)計算α0,并用mid(αmin,α0,αmax)替換α0.其中[αmin,αmax]是本文為了防止α0的值過大或過小而定義的一個區(qū)間,0<αmin<α0<αmax,而mid(αmin,α0,αmax)則是α0,αmin,αmax3個元素的中間值.

Step3(回溯行搜索法) 選擇αk作為序列α0,βα0,β2α0,…的第1個數(shù),那么就有

并令z(k+1)=(z(k)-α(k)F(z(k)))+.

Step4進(jìn)行收斂性測試,如果滿足那么就在z(k+1)處終止算法,否則令k←k+1,并回到Step 2.

基于上述分析,可以得出利用稀疏分解和提升小波變換求解此問題的算法流程如圖2所示.

圖2 算法流程

4 實驗與分析

將MEMS陀螺儀ADXRS300放置于安裝在隔離地基的溫控轉(zhuǎn)臺上,使用PC104工控機(jī)對輸出數(shù)據(jù)進(jìn)行采集處理.

待陀螺儀工作一段時間,輸出穩(wěn)定后,以200 Hz的采樣頻率采集一組MEMS陀螺儀ADXRS300的靜態(tài)數(shù)據(jù),作為陀螺儀隨機(jī)漂移的時間序列.取該時

間序列中長度為3 600的一段序列作為實驗數(shù)據(jù).然后分別采用本文方法和小波軟閾值濾波方法對陀螺儀輸出信號進(jìn)行處理,然后比較兩種方法的去噪效果.本文采用了db4作為小波基,對輸出信號進(jìn)行3尺度分解,小波軟閾值濾波方法選用了“Birge Massart”閾值.

對采集的數(shù)據(jù)分別用兩種方法進(jìn)行處理的結(jié)果如圖3、4所示.

從圖3、4可以看出,相對于小波閾值去噪方法,使用本文方法處理后的信號幅值下降更明顯.另外,利用本文方法進(jìn)行處理后的信號毛刺更少,也就是說處理后信號的零偏穩(wěn)定性更好.

為更好的驗證本文所提出去噪方法在MEMS陀螺信號處理中的可行性,設(shè)計了相應(yīng)的跑車試驗,具體試驗條件描述如下:將MEMS陀螺儀ADXRS300固定放置在車載平臺上,并使用PC104工控機(jī)對輸出信號進(jìn)行實時采集.實驗的流程為先150 s預(yù)熱,然后200 s的跑車,跑車路段包括直線和彎道,且不保持勻速行駛.

選擇一組實地跑車實驗中的MEMS陀螺儀數(shù)據(jù),分別使用本文方法和小波閾值去噪方法對其進(jìn)行分析處理.為了定量分析此兩種方法的性能,選擇均方差和信噪比兩項指標(biāo)來衡量其去噪效果.由均方差和信噪比的定義可知,同一信號去噪處理后,均方差越小,信噪比越大則去噪效果越好.

處理結(jié)果見表1.

圖3 小波濾波去噪效果

圖4 本文方法去噪效果

去噪方法SNR/dBMSE/10-4小波閾值法64.85643.3589本文方法70.36953.0125

從表1中可以看出,本文方法處理后的MEMS陀螺輸出信號,無論是信噪比還是均方差都優(yōu)于小波閾值去噪法.

圖5是本文所使用的梯度投影算法與基追蹤算法及l(fā)1_ls算法在解決此類問題時性能的對比圖.由圖5可以看出梯度投影算法更穩(wěn)定,且達(dá)到收斂所需迭代次數(shù)更少.

由于OMP并不是和GP同類型的凸松弛方法,不是帶約束的線性規(guī)劃問題,也就不存在目標(biāo)函數(shù),所以不能以目標(biāo)函數(shù)收斂所需時間和迭代次數(shù)作為指標(biāo)和GP對比.所以為充分對比這兩種算法的性能如何,分別使用這兩種算法稀疏重建一組信號x,x是在長度為n的全零信號中隨機(jī)插入m個非零值得到的仿真信號.y=Ax+N為觀測信號,N為白噪聲.

圖6分別是用這兩種方法對信號進(jìn)行處理并達(dá)到相同稀疏度后重構(gòu)信號的MSE和CPU時間的對比圖.

圖5 不同算法計算效率對比

圖6 GP與OMP性能對比

從圖6中可以看出:當(dāng)信號x的稀疏度比較小(x中的非零元素個數(shù)m>220)時,采用GP算法處理后的信號的MSE更小,而且很明顯GP的速度比OMP要快(除了m<50這種極端的情況).

而在分別使用這兩種方法對本文實測靜態(tài)數(shù)據(jù)進(jìn)行處理時,所需要的CPU時間分別為:tGPSR=1.248 0 s,tOMP=12.214 9 s,可見GP方法還是明顯要快于OMP算法.

5 結(jié) 論

1)本文將稀疏分解與提升小波變換相結(jié)合,并采用梯度投影算法恢復(fù)信號的稀疏性,進(jìn)而提出了基于梯度投影的改進(jìn)稀疏分解去噪算法,并用于MEMS陀螺信號的處理.

2)本文算法與小波軟閾值濾波方法對比結(jié)果表明,本文算法在動靜態(tài)條件下都可以有效地去除MEMS陀螺儀輸出信號中的噪聲,尤其是在靜態(tài)條件下的去噪效果要明顯優(yōu)于小波閾值濾波方法.

3) 本文采用的梯度投影算法與正交匹配追蹤算法和基追蹤算法的對比結(jié)果表明,梯度投影算法具有更高的運(yùn)算效率,適合在線濾波,同時具有良好的穩(wěn)定性和全局優(yōu)化能力.

[1] BARBOURN, CONNELLY J, GILMORE J, et al. Micromechanical silicon instrument and system development at draper laboratory [C]//AIAA Guidance Navigation and Control Conference. San Diego: AIAA, 1996: 1-7.DOI: 10.2514/6.1996-3709.

[2] 王昊,王俊璞, 田蔚風(fēng),等.梯度RBF神經(jīng)網(wǎng)絡(luò)在MEMS陀螺儀隨機(jī)漂移建模中的應(yīng)用[J].中國慣性技術(shù)學(xué)報,2006,14(4):44-48. DOI:10.13695/j.cnki.12-1222/o3.2006.04.010.

WANG Hao, WANG Junpu, TIAN Weifeng, et al. Application of gradient radial basis function network in the modeling of MEMS gyro’s random drift[J].Journal of Chinese Intertial Technology, 2006,14(4):44-48.DOI:10.13695/j.cnki.12-1222/o3.2006.04.010.

[3] 王新龍,李娜. MEMS陀螺儀隨機(jī)誤差的建模與分析[J].北京航空航天大學(xué)學(xué)報,2012, 38(2):170-174. DOI:10.13700/j.bh.1001-5965.2012.02.020.

WANG Xinlong, LI Na. Error modeling and analysis for random drift of MEMS gyroscope[J].Journal of Beijing University of Aeronautics and Astronautics, 2012, 38(2):170-174. DOI:10.13700/j.bh.1001-5965.2012.02.020.

[4] 李澤民,段鳳陽,馬佳智. 基于支持向量機(jī)的MEMS陀螺儀隨機(jī)漂移補(bǔ)償[J].傳感技術(shù)學(xué)報,2012,25(8):1084-1087. DOI: 10.3969/j.issn.1004-1699.2012.08.013.

LI Zemin, DUAN Fengyang, MA Jiazhi. Random drift compensation of MEMS gyroscope based on support vector machine[J].Chinese Journal of Sensors and Actuators, 2012, 25(8):1084-1087. DOI: 10.3969/j.issn.1004-1699.2012.08.013.

[5] 于湘濤,張?zhí)m.基于小波最小二乘支持向量機(jī)的加速度計溫度建模和補(bǔ)償[J].中國慣性技術(shù)學(xué)報,2011,19(1):95-98. DOI:10.13695/j.cnki.12-1222/o3.2011.01.014.

YU Xiangtao, ZHANG Lan. Temperature modeling and compensation of accelerometer based on least squares wavelet support vector machine[J]. Journal of Chinese Inertial Technology,2011, 19(1):95-98.DOI:10.13695/j.cnki.12-1222/o3.2011.01.014.

[6] 秦偉偉,鄭志強(qiáng), 劉剛,等.基于小波分析與LSSVM的陀螺儀隨機(jī)漂移建模[J]. 中國慣性技術(shù)學(xué)報, 2008, 16(6): 721-724, 729. DOI:10.13695/j.cnki.12-1222/o3.2008.06.013.

QIN Weiwei, ZHENG Zhiqiang, LIU Gang, et al. Modeling method of gyroscope’s random drift based on wavelet analysis and LSSVM[J]. Journal of Chinese Inertial Technology, 2008, 16(6):721-724, 729. DOI:10.13695/j.cnki.12-1222/o3.2008.06.013.

[7] 石光明,劉丹華, 高大化,等. 壓縮感知理論及其研究進(jìn)展[J]. 電子學(xué)報,2009,37(5):1070-1081.DOI: 10.3321/j.issn:0372-2112.2009.05.028.

SHI Guangming, LIU Danhua, GAO Dahua, et al. Advances in theory and application of compressed sensing [J]. Acta Electronica Sinica, 2009, 37(5):1070-1081.DOI: 10.3321/j.issn:0372-2112.2009.05.028.

[8] TSAIG Y, DONOHO D L. Extensions of compressed sensing [J].Signal Processing, 2006, 86(3):549-571.DOI: 10.1016/j.sigpro.2005.05.029.

[9] FIGUEIREDO M A T, NOWAK R D, WRIGHT S J. Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems[J]. IEEE Journal of Selected Topics in Signal Processing,2007, 1(4): 586-594. DOI: 10.1109/JSTSP.2007.910281.

[10]程承,潘泉, 王申龍,等.基于壓縮感知理論的MEMS陀螺儀信號降噪研究[J]. 儀器儀表學(xué)報,2012, 33(4): 769-773. DOI:10.3969/j.issn.0254-3087.2012.04.008.

CHENG Cheng, PAN Quan, WANG Shenlong, et al. Research on MEMS gyroscope signal denoising based compressed sensing theory[J].Chinese Journal of Scientific Instrument, 2012,33(4):769-773. DOI:10.3969/j.issn.0254-3087.2012.04.008.

[11]張海鵬,房建成. MEMS陀螺儀短時漂移特性實驗研究[J].中國慣性技術(shù)學(xué)報, 2007,15(1):100-104. DOI:10.3969/j.issn.1005-6734.2007.01.026.

ZHANG Haipeng, FANG Jiancheng. Short-time drift characteristic of MEMS Gyroscope[J]. Journal of Chinese Inertial Technology, 2007, 15(1):100-104. DOI:10.3969/j.issn.1005-6734.2007.01.026.

[12]DAI Wei, MILENKOVIC O. Subspace pursuit for compressive sensing signal reconstruction[J]. IEEE Transactions on Information Theory, 2009, 55(5): 2230-2249. DOI: 10.1109/TIT.2009.2016006.

[13]LIVSHITZ E. On efficiency of orthogonal matching pursuit[M]. [S.L.]: Preprint, 2010.

[14] MALLAT S G, ZHANG Zhifeng. Matching pursuits with time-frequency dictionaries[J]. IEEE Transactions on Signal Processing, 1993, 41(12):3397-3415. DOI: 10.1109/78.258082.

[15]ZHU Wenjie, WANG Guanglong, QIAO Zhongtao, et al. A novel noise reduction algorithm of mems gyroscope based on compressive sensing and lifting wavelet transform[J]. Key Engineering Materials, 2014, 609-610: 1138-1143. DOI:10.4028/www.scientific.net/KEM.609-610.1138.

AnovelnoisereductionmethodforMEMSgyroscope

SUN Tianchuan, LIU Jieyu

(Dept. of Control Engineering, Rocket Force University of Engineering, Xi’an 710025, China)

To get a better de-noising effect, a novel noise reduction method combining the sparse decomposition with lifting wavelet transform is proposed. Firstly, the error model is established for the MEMS gyroscope output signal with noise, and wavelet coefficient of signals with noise can be obtained by lifting wavelet transform. Then the sparsity of the coefficient is recovered according to sparse decomposition theory. Finally, signals are reconstructed by lifting wavelet inverse transform, i.e. the de-noised signal is thus obtained. In addition, since the gradient projection algorithm is global optimal algorithm with high computational efficiency, the theory of gradient projection is used in the restoration of sparse signal. Specifically, a sparse decomposition based on gradient projection is designed to simplify the algorithm complexity and improve the stability of the algorithm. To verify the performance of the proposed algorithm, the static experiment and dynamic car test on MEMS gyroscope are implemented. The results show that the denoising performance of the new method is better than that of wavelet filter either under the static or dynamic condition, especially under the latter condition. Meanwhile, the CPU time of gradient projection is less than orthogonal matching pursuit (OMP) and basis pursuit (BP).

MEMS gyroscope; signal denoising; sparse decomposition; lifting wavelet transform; gradient projection; convex optimization

10.11918/j.issn.0367-6234.201606079

V241.5

A

0367-6234(2017)10-0060-06

2016-06-22

國家自然科學(xué)基金(61304001)

孫田川(1993—),男,博士研究生;

劉潔瑜(1970—),女,教授,博士生導(dǎo)師

劉潔瑜,jieyu.liu@gmail.com

(編輯張 紅)

猜你喜歡
信號方法
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
學(xué)習(xí)方法
孩子停止長個的信號
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 精品国产电影久久九九| 免费观看欧美性一级| 中文字幕无码av专区久久| 亚洲色图另类| 国产成人精品一区二区三在线观看| 九九九精品视频| 日本手机在线视频| 天天躁夜夜躁狠狠躁躁88| 欧美精品亚洲日韩a| 成人在线观看不卡| 久久永久视频| 免费国产福利| 亚洲欧美另类日本| 干中文字幕| 日本国产一区在线观看| 四虎影院国产| 男人天堂亚洲天堂| 欧美成人精品一级在线观看| 香蕉久久国产超碰青草| 国产成人综合日韩精品无码不卡| 国产精品浪潮Av| 小说 亚洲 无码 精品| 国产 在线视频无码| 青青青国产免费线在| 久久无码高潮喷水| 女人18毛片水真多国产| 青草视频在线观看国产| 97se亚洲综合| 在线观看的黄网| 亚洲日韩精品无码专区| a天堂视频在线| 亚洲精品成人片在线播放| 第一页亚洲| 亚洲精品天堂自在久久77| 欧美亚洲国产精品久久蜜芽| 久久精品国产亚洲麻豆| 亚洲一道AV无码午夜福利| 精品国产电影久久九九| 欧美一级高清片欧美国产欧美| 女人18毛片一级毛片在线 | 国产精品亚欧美一区二区| 国产91丝袜在线播放动漫 | 在线免费a视频| 亚洲国产精品日韩av专区| 国产精品视频猛进猛出| 五月天久久综合| 亚洲综合色区在线播放2019| 国产精品19p| 日韩最新中文字幕| 欧美一级在线| 怡春院欧美一区二区三区免费| 成年片色大黄全免费网站久久| 欧美一级黄片一区2区| 粉嫩国产白浆在线观看| 亚洲第一视频区| 青青热久免费精品视频6| 欧美亚洲国产一区| 成人午夜久久| 91在线一9|永久视频在线| 国模粉嫩小泬视频在线观看| 青青草原国产一区二区| 视频在线观看一区二区| 天天综合网色| 国产产在线精品亚洲aavv| 特级aaaaaaaaa毛片免费视频| 亚洲天堂在线免费| 久99久热只有精品国产15| 久久99精品久久久久纯品| 宅男噜噜噜66国产在线观看| 久久久久88色偷偷| 伦精品一区二区三区视频| 欧美一级夜夜爽www| 亚洲国产清纯| 日韩欧美国产区| 国产哺乳奶水91在线播放| a欧美在线| 精品国产黑色丝袜高跟鞋| 91偷拍一区| AV不卡在线永久免费观看| 精品自窥自偷在线看| 青草免费在线观看| 色综合狠狠操|