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

動態(tài)總方差及其在陀螺振動信號分析中的應(yīng)用*

2015-08-24 02:53:41朱戰(zhàn)輝汪立新
傳感技術(shù)學(xué)報 2015年12期
關(guān)鍵詞:振動信號分析

朱戰(zhàn)輝,汪立新,李 燦

(第二炮兵工程大學(xué)控制科學(xué)與工程系,西安710025)

動態(tài)總方差及其在陀螺振動信號分析中的應(yīng)用*

朱戰(zhàn)輝,汪立新*,李燦

(第二炮兵工程大學(xué)控制科學(xué)與工程系,西安710025)

動態(tài)Allan方差是分析動態(tài)環(huán)境下陀螺儀隨機誤差變化規(guī)律的一種新方法。針對其運用窗函數(shù)截取信號導(dǎo)致方差估計值置信度降低,尤其是長相關(guān)時間上估計誤差大的問題,提出了動態(tài)總方差法。用矩形窗分段截取陀螺儀量測信號,再對截斷窗內(nèi)數(shù)據(jù)進行倒像映射延拓以增加方差估計的實際自由度,最后計算延拓后樣本的Allan方差并提取其噪聲系數(shù),將其按時間順序分別以三維和二維的形式表征出來,以細化和辨識陀螺輸出的動態(tài)特性。從對半球諧振陀螺線振動試驗實測數(shù)據(jù)處理結(jié)果來看,新算法不僅能及時跟蹤信號的非平穩(wěn)變化,而且有效提高了振動信號方差估計值的置信度。

陀螺儀;動態(tài)特性;動態(tài)Allan方差;總方差;振動信號

EEACC:7230;7310Gdoi:10.3969/j.issn.1004-1699.2015.12.010

陀螺儀是慣性導(dǎo)航系統(tǒng)的核心部件,其靜動態(tài)性能的優(yōu)劣直接影響到導(dǎo)彈、火箭及其它飛行器的工作精度。在慣導(dǎo)系統(tǒng)工作過程中,機械振動是不可避免且長期存在的,陀螺能否在振動環(huán)境下保持高穩(wěn)定性與高精度成為影響整個慣導(dǎo)系統(tǒng)精度的重要因素[1]。Allan方差是一種從時域數(shù)據(jù)中分析振蕩器頻率穩(wěn)定性的方法,能夠?qū)ν勇葺敵鼋撬俾手写嬖诘母鱾€噪聲項和整個噪聲系統(tǒng)進行細致的表征和辨識[2]。但該方法只能對陀螺靜態(tài)信號隨機誤差進行分析,不能表征其在動態(tài)激勵下的時域非平穩(wěn)特征。因此近年來有學(xué)者將用于分析原子時鐘頻率非平穩(wěn)性的動態(tài)Allan方差DAVAR(Dynamic Allan Variance)引入到陀螺隨機噪聲的分析中來[3-4]。文獻[5-6]將該方法應(yīng)用到光纖陀螺振動和變溫條件下的性能評價上,有效地表征出了輸出信號在動態(tài)試驗中所表現(xiàn)出的時變特性。然而由于動態(tài)Allan方差法引入了窗函數(shù)截斷的理論,使截斷窗內(nèi)參與Allan方差計算的數(shù)據(jù)相對整個序列來說大幅減少,這對短相關(guān)時間下噪聲系數(shù)的擬合來說影響較小,但中、長相關(guān)時間下噪聲系數(shù)(低頻噪聲系數(shù))的擬合精度將大大降低。如果一味的通過增大窗寬的方法來增加方差估計的置信度,必然會降低動態(tài)跟蹤效果,與動態(tài)Allan方差方法的設(shè)計原則相背離。因此,要想提高低頻噪聲系數(shù)的擬合精度就必須立足于有限數(shù)據(jù)量去改進方差估計算法,以提高中、長相關(guān)時間下的估計置信度。

實際上,學(xué)者們已經(jīng)對有限數(shù)據(jù)量下提高方差估計置信度的方法開展了大量的研究。比較成熟的方法主要有重疊Allan方差,是通過重疊采樣形成所有可能的相關(guān)時間為τ的子序列,來最大限度的利用現(xiàn)有數(shù)據(jù)估計Allan方差[7]。Allan總方差是采用倒鏡像映射的方法,對數(shù)據(jù)兩段進行延伸,然后再對延伸后的數(shù)據(jù)計算其Allan方差,可以有效提高長相關(guān)時間下方差的估計精度[8]。#1理論方差是模仿Allan方差的較新的統(tǒng)計方法,其相關(guān)時間可以從10到0.75 N,具有更好的置信度,但#1理論方差是對Allan方差的有偏估計,因此需要對其做偏差修正[9]。ThêoH方差是去除了#1理論方差與Allan方差之間偏差,再與Allan方差合成的較新的統(tǒng)計方法,尤其適合于相關(guān)時間較長和多種噪聲混合存在的情況,可以從相同的數(shù)據(jù)序列中最大程度的獲得更多的統(tǒng)計信息[10]。韓軍良、石國祥等將總方差法引入到光纖陀螺靜態(tài)信號的隨機誤差分析中,有效的提高了長相關(guān)時間下噪聲系數(shù)的估計值置信度[11-12]。本文將動態(tài)總方差法提出的倒鏡像延伸數(shù)據(jù)的思想運用到DAVAR算法截斷窗窗內(nèi)數(shù)據(jù)估計置信度的提高上,有效改善了中、長相關(guān)時間下隨機噪聲系數(shù)的擬合精度。通過對線振動試驗中半球諧振陀螺量測信號的分析,驗證了本文算法的有效性。

1 動態(tài)Allan方差

DAVAR法是Allan方差法的擴展,其原理是用固定長度的矩形窗函數(shù)分別截取信號在不同時段內(nèi)的數(shù)據(jù),計算其Allan方差,并把方差的集合、時間和相關(guān)時間(平均因子)在同一幅三維圖上表示出來,以表征信號方差在時間域的變化特征。進一步還可以用最小二乘法提取出每個截斷窗內(nèi)樣本的噪聲系數(shù),再分別按時間順序上排列出來,繪制出各個噪聲系數(shù)的時間變化曲線。

2 總方差

2.1總方差法的原理

Allan方差的估計是基于有限長度的數(shù)據(jù),估計的可信度依賴于數(shù)據(jù)的獨立組數(shù)。對于給定的隨機序列,如果數(shù)據(jù)樣本越少,可劃分的組數(shù)就越少,Allan方差的估計置信度較差,估計誤差就越大[11]。針對此種情況,總方差法通過對原始數(shù)據(jù)進行倒鏡像映射延拓的辦法來增加自由度,其原理如下:

設(shè)有一個時間序列信號xi(i=1,2,…,N),采樣時間為τ0,將這個序列通過映射延伸成一個新的、更長的虛擬序列x?i,映射的規(guī)則為

其經(jīng)過橫軸平移后的映射示意圖如圖1。

圖1 R總方差定義下的數(shù)據(jù)映射示意圖

得到虛擬序列的長度接近于原始序列的三倍,式(1)為映射后序列的計算公式,其中m為平均因子,1≤m≤N-1

對比二者的定義可以看出,Allan方差估計的相關(guān)時間τ只能達到,而總方差估計的相關(guān)時間則達到。

圖2為對同一組陀螺輸出進行Allan方差分析和總方差法分析的雙對數(shù)曲線圖,可見在長相關(guān)時間,也就是大的平均因子下,總方差法是收斂的,且具有更大的τ值,而Allan方差是發(fā)散的。

圖2 RAllan方差和總方差雙對數(shù)曲線比較圖

2.2總方差法的估計誤差

用總方差法分析陀螺數(shù)據(jù)時,僅給出方差的估計值是不完整的,還應(yīng)該給出這個估計值的置信度,置信度取決于自由度。在樣本數(shù)據(jù)量大小不能改變前提下,文獻[10]運用交疊Allan方差、總方差和ThêoH方差對頻率穩(wěn)定度進行了分析,結(jié)論是在τ=T/2情況下,ThêoH方差的自由度大于總方差的自由度,總方差的自由度又大于交疊Allan方差的自由度。但是ThêoH方差的計算量特別大,不適合實時在線分析。因此我們選用了總方差來替換Allan方差進行陀螺隨機誤差系數(shù)的動態(tài)提取。文獻[7]通過對陀螺和加速度計的分析,也證明了在長相關(guān)時間下總方差的估計準(zhǔn)確度要大于交疊Allan方差和Allan方差。

3 動態(tài)總方差算法設(shè)計

動態(tài)總方差DTVAR(Dynamic Total Variance)借用DAVAR法用窗函數(shù)截斷數(shù)據(jù)的思想,分段截取數(shù)據(jù)后按照一定規(guī)則進行數(shù)據(jù)延拓,計算其Allan方差和標(biāo)準(zhǔn)差并將其按時間序列排列,DTVAR方法的具體實現(xiàn)步驟如下:

圖3 RDTVAR算法設(shè)計流程圖

①確定隨機信號x(t)分析時間起始點t1。

②用中心點為t1,寬度為L的窗函數(shù)截斷隨機信號x(t),獲得窗口截斷信號yT(t1,t′),支撐變量t′描述窗內(nèi)漸逝的時間

PW(t′)為長度L的窗函數(shù),截取得到的信號為

③將窗口截斷信號yT(t1,t′)按照式(1)進行倒像映射延伸,得到虛擬序列

④將y*T(t1,t′)同Allan窗口hτ(t′)做卷積建立增量過程Δ(t1,t′,τ)

這時變量t′的范圍變?yōu)?/p>

0≤τ≤τmax,通常取τmax=L/3,則動態(tài)總方差在t1時刻的估計值為

Allan方差可以被定義為上式的總體期望值

對比Allan方差和總方差的定義可以看出,Allan方差估計的相關(guān)時間τ只能達到Nτ0/2,而總方差估計得相關(guān)時間τ則達到Nτ0,估計范圍得到了擴展。

⑤陀螺的部分主要隨機噪聲系數(shù)可以通過最小二乘擬合分離出來,公式如下

其中Q、N、B、K、R分別代表了陀螺的噪聲系數(shù):量化噪聲(quantization noise)、隨機游走系數(shù)(angle random walk)、零偏不穩(wěn)定性(bias instability)、速率隨機游走(rate random walk)和速率斜坡(the rate slope)。

⑥將窗口滑動到t2,用窗函數(shù)PW(t′)截取x(t),獲得窗口截斷信號yT(t2),重復(fù)步驟(2)~步驟(6),得到σ(t2,τ),以此類推,可以得到時間域的總方差序列σ(tN,τ)。其中tn+1截斷的窗口最好與以tn為中心截斷的窗口相重疊。

⑦最終將動態(tài)Allan方差序列σ2(tN,τ)或動態(tài)Allan標(biāo)準(zhǔn)差DADEV(Dynamic Allan Deviation)σ(tN,τ)按時間順序繪制在σ(tN,τ)-τ-t的三維圖上,從宏觀上表征隨機誤差的動態(tài)特征。

4 動態(tài)試驗驗證

半球諧振陀螺(HRG)是一種基于哥氏效應(yīng)的陀螺儀,具有啟動時間短,長期穩(wěn)定性好等特點[13]。本文以71號陀螺在線振動試驗中采集到的角速率量測信號為原始數(shù)據(jù),振動試驗初始加速度為5 g,線振動固定頻率為1.5 Hz,采樣時間為0.1 s,采集數(shù)據(jù)6 173個。對該信號進行預(yù)處理,得到隨機誤差信號如圖5所示。

圖4 R陀螺線振動實驗臺

圖5 R振動條件下半球諧振陀螺量測信號

分別運用DTVAR和DAVAR方法對半球諧振陀螺量測信號進行三維分析的結(jié)果如圖6和圖7所示,可以看到,動態(tài)過程在三維圖中得到了表征,顏色亮度的變化同標(biāo)準(zhǔn)差數(shù)值大小成正比。此外,圖7比圖6在長相關(guān)時間下方差估計值的波動要大的多,這正是Allan方差在長相關(guān)時間下自由度低于總方差法造成的。比較看來,DTVAR方法能更好的對HRG進行動態(tài)特征分析。

圖6 R半球諧振陀螺振動試驗信號的DTVAR分析

圖7 R半球諧振陀螺振動試驗信號的DAVAR分析

將隨機噪聲系數(shù)中比較重要的角度隨機游走和零偏不穩(wěn)定性分別用DTVAR和DAVAR方法對其進行辨識,再用最小二乘法擬合出來,其中DTVAR用的是401窗寬,DAVAR分別用401窗寬、801窗寬和1 201窗寬,比較結(jié)果如圖8、圖9所示。

圖8 R零偏不穩(wěn)定性的DAVAR和DTVAR比較分析

圖9 R角度隨機游走的DAVAR和DTVAR比較分析

DAVAR方法提出者Galleani L指出,截斷窗越短,動態(tài)跟蹤能力越強,但方差估計值的置信度較低;截斷窗越長,方差估計值的置信度較高,動態(tài)跟蹤能力越差,兩者是一組矛盾,通常很難平衡。而DAVAR方法在一定程度上解決了這個矛盾,即用較短的窗寬截取信號依然可以保持較高的置信度。為驗證這一點,我們分別用DAVAR和DTVAR方法對振動段信號進行分析,截斷窗長同為401。從圖8和圖9中可以觀測到,無論是角度隨機游走還是零偏不穩(wěn)定性,DAVAR方法的方差估計值總是大于DTVAR。要想使DAVAR與DTVAR的估計值相同就需要增加DAVAR截斷窗長,零偏不穩(wěn)定性噪聲系數(shù)的分析窗長需要增加到1 201,而角度隨機游走系數(shù)的分析窗長需要增加到801。然而,隨之而來的是動態(tài)跟蹤能力的下降。

從圖8和圖9還可以看到:DAVAR方法1 201窗長和801窗長所定位的突變發(fā)生點都較大的偏離了真實的突變發(fā)生點。這是因為長窗寬截斷窗包含了更多的數(shù)據(jù),過早的把振動信號包含了進來;振動結(jié)束后,窗口又遲遲不能離開振動發(fā)生點。此外,還可以看到1 201窗長使整個振動過程都發(fā)生了嚴重的畸變,振動過渡段被大幅拉長,而實際上,無論振動開始,還是振動結(jié)束,其過渡過程都是很短暫的。DTVAR的估計值比DAVAR的波動性小,更適合分析動態(tài)特性及動態(tài)變化時間點的準(zhǔn)確定位。

HRG動態(tài)試驗是固定頻率的正弦線振動試驗,所以當(dāng)振動進入平穩(wěn)振動段,也就是2 000采樣點到3 500采樣點時,我們可以認為此時間段里陀螺的輸出也是平穩(wěn)的,DTVAR三維分析圖也驗證了這個推斷。為更準(zhǔn)確的比較這兩種算法的準(zhǔn)確度,我們參照文獻[14]的方法,分別對平穩(wěn)振動段(2 000~3 500采樣點)和靜態(tài)段(4 300~5 800采樣點)的信號展開分析,分別繪制出陀螺5個主要隨機噪聲系數(shù)在這兩個階段的動態(tài)變化曲線。同時用這兩個數(shù)據(jù)段的Allan方差值作為標(biāo)準(zhǔn)值,截斷窗內(nèi)樣本數(shù)據(jù)為401個屬于小樣本,而該數(shù)據(jù)段包含了1 500個數(shù)據(jù),樣本數(shù)據(jù)量大,顯然置信度要高于截斷窗內(nèi)數(shù)據(jù)的置信度,因此將其作為大樣本下的標(biāo)準(zhǔn)值。對兩種算法展開比較,結(jié)果如表1和表2所示。可以看出:動態(tài)總方差法的噪聲系數(shù)擬合準(zhǔn)確度要高于DAVAR方法。尤其在中、長相關(guān)時間下的噪聲系數(shù),如速率隨機游走和速率斜坡的擬合精度上,有了較大改善。

圖10是分別用DTVAR401窗長,DAVAR401窗長和DTVAR1201窗長對信號進行辨識獲得的零偏穩(wěn)定性時間變化曲線。48.67是振動段大數(shù)據(jù)樣本(2 000~3 500點)下用Allan方差辨識出來的高置信度基準(zhǔn)值。可以看到,用DTVAR401窗長提取的振動條件下零偏穩(wěn)定性與基準(zhǔn)值最為接近。

圖10 R振動條件下零偏不穩(wěn)定性比較分析圖

表1 R振動條件下噪聲系數(shù)的方差估計值

表2 R靜態(tài)條件下噪聲系數(shù)的方差估計值

5 結(jié)論

針對動態(tài)Allan方差法用窗函數(shù)截斷數(shù)據(jù)造成方差估計值置信度低,導(dǎo)致噪聲系數(shù)擬合精度不高的問題,提出了基于倒鏡像映射延伸數(shù)據(jù)的動態(tài)總方差法,用總方差替代Allan方差進行截斷窗內(nèi)數(shù)據(jù)的方差估計。將提出的動態(tài)總方差法應(yīng)用于半球諧振陀螺線振動試驗實測數(shù)據(jù)的動態(tài)特征分析,結(jié)果表明:新算法在不降低動態(tài)跟蹤效果的前提下,有效提高了截斷窗內(nèi)方差估計值的置信度,尤其是中、長相關(guān)時間,也就是大平均因子下的方差估計置信度。此外,該算法具有較快的運算速度,可以用于陀螺噪聲系數(shù)的實時在線提取及分析。

[1]方琳,申沖,陳熙源.基于小波多尺度變換的光纖陀螺振動誤差分析與補償[J].傳感技術(shù)學(xué)報,2012,25(7):902-906.

[2]尹杭,張偉,袁琳峰.一種MEMS加速度計誤差分析與校準(zhǔn)方法[J].傳感技術(shù)學(xué)報,2014,27(7):866-869.

[3]Galleani L.The Dynamic Allan Variance III:Confidence and De-tection Surfaces.IEEE Transactions on Ultrasonics,F(xiàn)erroelectrics,and Frequency Control,2011,58(8):1550-1558.

[4]Galleani L,Tavella P.The Dynamic Allan Variance[J].IEEE Transaction on Ultrasonics,F(xiàn)erroelectrics,and Frequency Control,2009,56(3):450-460.

[5]李冀辰,高鳳岐,王廣龍,等.光纖陀螺振動和變溫條件下的DAVAR分析[J].中國激光,2013,40(9):09080041-09080047.

[6]Zhang Chunxi,Wang Lu,Gao Shuang,et al.Dynamic Allan Variance Analysis for Stochastic Errors of Fiber Optic Gyroscope[J]. Infrared and Laser Engineering,2014,43(9):3081-3088.

[7]Li Jintao,F(xiàn)ang Jiancheng.Not Fully Overlapping Allan Variance and Total Variance for Inertial Sensor Stochastic Error Analysis[J].IEEE Transactions on Instrumentation and Measurement,2013,62(10):2659-2672.

[8]Howe D A.The Total Deviation Approach to Long-Term Characterization of Frequency Stability[J].IEEE Transactions on Ultrasonic,F(xiàn)erroelectrics,and Frequency Control,2000,47(5):1102-1110.

[9]程旭維,湯霞清,黃湘遠.基于#1理論方差的光學(xué)陀螺長期隨機誤差分析[J].中國激光,2014,41(10):10005003.1-10005003.8.

[10]Howe D A.ThêoH,A Hybrid,High-Confidence Statistic That Improves on the Allan Deviation[J].Metrologia,2006,43:322-331.

[11]韓軍良,葛升民,沈毅.基于總方差方法的光纖陀螺隨機誤差特性研究[J].哈爾濱工業(yè)大學(xué)學(xué)報,2007,39(5):708-711.

[12]石國祥,陳堅,葉軍,等.總方差方法在光纖陀螺隨機誤差分析中的應(yīng)用[J].光電工程,2012,39(1):63-67.

[13]李燦,汪立新,秦偉偉,等.諧振子偏心對半球諧振陀螺精度的影響[J].傳感技術(shù)學(xué)報,2015,28(6):831-835.

[14]焦月,張升康.頻率穩(wěn)定度測量方法的對比分析[J].宇航計測技術(shù),2013,33(1):35-38.

Dynamic Total Variance and Its Application in Analysis of Dynamic Characteristics of Gyroscope*

ZHU Zhanhui,WANG Lixin*,LI Can
(The Second Artillery Engineering University,Xi'an 710025,China)

The dynamic Allan variance(DAVAR)is a new method for analyzing random error of gyroscope.However,it has defects such as poor confidence on the estimate and great estimation error at long-term τ-values,due to the reduced amount of data captured by the truncation windows.An improvement dynamic Allan variance is proposed to increase the confidence in the case of long-term τ-values.Firstly,the random error of FOG is truncated with the window function.Then the data in the windows are extended by the total variance method to improve the confidence on the estimate,At last the Allan variance of extended data is calculated,and the noise coefficients of gyro can also be identified and extracted at the same time.At last,they are expressed by three or two-dimensional to describe the dynamic characteristics of gyro.The measured data of linear acceleration test for Hemispherical Resonator Gyro is analyzed with proposed algorithm and dynamic Allan variance.The results shows that proposed algorithm can track the non-stationary of gyroscope more effectively and obtain a lower estimation error at long-term τ-values.

gyroscope;dynamic characteristics;dynamic Allan variance;total variance;vibration signal

朱戰(zhàn)輝(1978-),男,河南人,博士研究生,控制科學(xué)與工程專業(yè),主要研究方向為慣性系統(tǒng)及測試,數(shù)字信號處理,zzhhit@sina.com;

汪立新(1966-),男,湖北人,博士學(xué)位,教授,博士生導(dǎo)師,主要研究方向為:慣性技術(shù)、組合導(dǎo)航及信號測試與處理,wlxxian@163.com。

V241.5

A

1004-1699(2015)12-1789-06

項目來源:國家自然科學(xué)基金項目(61503390)

2015-07-17修改日期:2015-09-15

猜你喜歡
振動信號分析
振動的思考
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
隱蔽失效適航要求符合性驗證分析
完形填空二則
振動與頻率
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
基于FPGA的多功能信號發(fā)生器的設(shè)計
電子制作(2018年11期)2018-08-04 03:25:42
中立型Emden-Fowler微分方程的振動性
電力系統(tǒng)及其自動化發(fā)展趨勢分析
基于LabVIEW的力加載信號采集與PID控制
主站蜘蛛池模板: 午夜a级毛片| 亚洲午夜天堂| 亚洲欧美成人在线视频| 波多野结衣在线se| 福利小视频在线播放| 国产麻豆91网在线看| 国产亚洲现在一区二区中文| 尤物午夜福利视频| 亚洲综合色区在线播放2019| 欧美色亚洲| 国产簧片免费在线播放| 日本黄色不卡视频| 美女内射视频WWW网站午夜| 日韩欧美亚洲国产成人综合| 国产精品久久久久久久久久久久| 免费国产无遮挡又黄又爽| 一区二区三区成人| 国产精品亚洲五月天高清| 国产精品久久久久鬼色| 色欲色欲久久综合网| 久久人人妻人人爽人人卡片av| 久久久久亚洲AV成人网站软件| 国产H片无码不卡在线视频| AV天堂资源福利在线观看| 在线观看国产黄色| 国产成人高清亚洲一区久久| 国产区人妖精品人妖精品视频| 在线视频一区二区三区不卡| 欧美国产三级| 国产第三区| 亚洲欧洲美色一区二区三区| 久久国产热| 成人国产精品一级毛片天堂| 久久久久人妻一区精品| 国产白浆一区二区三区视频在线 | 亚洲欧美日韩色图| 青青青伊人色综合久久| 日本在线国产| 日本不卡在线视频| 91精品国产丝袜| 久久久亚洲色| 国产美女精品一区二区| 精品久久久久成人码免费动漫| 玖玖免费视频在线观看| …亚洲 欧洲 另类 春色| 天堂av综合网| 日韩一区二区三免费高清| 天天综合网色中文字幕| 久热这里只有精品6| 一级看片免费视频| 中文字幕乱码二三区免费| 欧美精品H在线播放| 国产白浆视频| 国产99精品久久| 456亚洲人成高清在线| 久久久国产精品无码专区| 五月综合色婷婷| 国产精品综合色区在线观看| 在线观看免费人成视频色快速| 欧美一级色视频| 色综合中文综合网| 五月婷婷精品| 激情乱人伦| 欧美精品三级在线| 九色在线观看视频| 99r在线精品视频在线播放| 视频一区视频二区中文精品| 久久久久九九精品影院| 国产在线精彩视频论坛| 欧美日韩亚洲国产| 欧美一级片在线| 免费人成在线观看视频色| 2021最新国产精品网站| 久久a级片| 五月婷婷激情四射| 日韩精品一区二区三区视频免费看| 免费AV在线播放观看18禁强制| 操国产美女| 狠狠操夜夜爽| 91免费观看视频| av无码一区二区三区在线| 免费 国产 无码久久久|