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
賺錢方法
捕魚
主站蜘蛛池模板: 成人精品视频一区二区在线| 亚洲侵犯无码网址在线观看| 99热这里都是国产精品| 欧美a在线看| 欧美亚洲一区二区三区在线| 久久久噜噜噜| 伊人久久大香线蕉成人综合网| 天天操精品| 欧洲亚洲一区| 97国产精品视频自在拍| 在线观看国产精美视频| 97se亚洲综合在线天天| 欧美国产日韩在线| 91啪在线| 亚洲精品另类| 国产乱子伦精品视频| 亚洲午夜综合网| 国产9191精品免费观看| 亚洲国产精品久久久久秋霞影院| 久久女人网| www.亚洲一区| 中国一级特黄视频| 国产亚洲欧美在线人成aaaa| 久久国产精品夜色| 成人精品免费视频| 国产精品亚欧美一区二区| 性激烈欧美三级在线播放| 一本大道无码日韩精品影视| 狠狠五月天中文字幕| 欧美成a人片在线观看| 欧美在线综合视频| 污视频日本| 福利小视频在线播放| 97国产成人无码精品久久久| 色亚洲成人| 亚洲一级毛片在线观| 欧洲日本亚洲中文字幕| 亚洲高清日韩heyzo| 亚洲综合经典在线一区二区| 亚洲精品福利网站| 国产精品无码AV片在线观看播放| 亚洲成a人片77777在线播放| 亚洲欧美在线看片AI| 99久久人妻精品免费二区| 影音先锋丝袜制服| AV色爱天堂网| 毛片视频网址| 日本黄色不卡视频| 99精品视频播放| 国产欧美在线| 国产精品天干天干在线观看| 在线观看精品国产入口| 欧美综合区自拍亚洲综合天堂| 五月天综合网亚洲综合天堂网| 欧美激情二区三区| 亚洲天堂网在线观看视频| 久久精品无码国产一区二区三区| 中文纯内无码H| 四虎成人免费毛片| 亚洲精品欧美日韩在线| 欧美色丁香| 亚洲无码精品在线播放| 狠狠色丁香婷婷| 永久免费av网站可以直接看的| 日韩精品成人网页视频在线 | AV老司机AV天堂| av色爱 天堂网| 91麻豆国产在线| 成AV人片一区二区三区久久| 欧美日韩中文字幕在线| 日韩美女福利视频| 男人天堂亚洲天堂| 国产成人无码AV在线播放动漫 | 亚洲国产精品美女| 在线观看无码av五月花| 亚洲娇小与黑人巨大交| 18禁黄无遮挡免费动漫网站| 亚洲视频欧美不卡| 国产精品流白浆在线观看| 国产香蕉97碰碰视频VA碰碰看| 91久久青青草原精品国产| 国产无码网站在线观看|