《電子技術(shù)應(yīng)用》
您所在的位置:首頁 > 模擬設(shè)計(jì) > 設(shè)計(jì)應(yīng)用 > 慣性導(dǎo)航系統(tǒng)中浮點(diǎn)計(jì)算加速單元設(shè)計(jì)
慣性導(dǎo)航系統(tǒng)中浮點(diǎn)計(jì)算加速單元設(shè)計(jì)
2019年電子技術(shù)應(yīng)用第8期
田換換1,2,朱曉燕1
1.首都師范大學(xué) 信息工程學(xué)院,北京100048;2.電子系統(tǒng)可靠性技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室,北京100048
摘要: 石英振梁加速度計(jì)采用頻率輸出的形式表示加速度,在慣性導(dǎo)航系統(tǒng)中,需要將頻率值轉(zhuǎn)換為加速度值,再進(jìn)行姿態(tài)解算。采用軟件方法進(jìn)行浮點(diǎn)計(jì)算,需要耗費(fèi)CPU大量的計(jì)算能力。為了優(yōu)化頻率轉(zhuǎn)換的計(jì)算速度,設(shè)計(jì)一種面向頻率轉(zhuǎn)換應(yīng)用的浮點(diǎn)計(jì)算加速單元,并基于FPGA進(jìn)行了實(shí)現(xiàn)與驗(yàn)證。結(jié)果表明,系統(tǒng)從數(shù)據(jù)采樣到頻率轉(zhuǎn)換,然后將頻率值轉(zhuǎn)換成加速度進(jìn)行姿態(tài)解算,陀螺儀測得的角速度進(jìn)行積分,最后完成數(shù)據(jù)融合,使用本文設(shè)計(jì)的浮點(diǎn)加速單元來實(shí)現(xiàn)頻率轉(zhuǎn)換,速度提高了2倍。
中圖分類號: TP33
文獻(xiàn)標(biāo)識碼: A
DOI:10.16157/j.issn.0258-7998.190095
中文引用格式: 田換換,朱曉燕. 慣性導(dǎo)航系統(tǒng)中浮點(diǎn)計(jì)算加速單元設(shè)計(jì)[J].電子技術(shù)應(yīng)用,2019,45(8):79-82,86.
英文引用格式: Tian Huanhuan,Zhu Xiaoyan. Design of floating point calculation acceleration unit in inertial navigation system[J]. Application of Electronic Technique,2019,45(8):79-82,86.
Design of floating point calculation acceleration unit in inertial navigation system
Tian Huanhuan1,2,Zhu Xiaoyan1
1.College of Information Engineering,Capital Normal University,Beijing 100048,China; 2.Beijing Key Laboratory of Electronic System Reliability and Prognostics,Beijing 100048,China
Abstract: Quartz vibrating beam accelerometer expresses acceleration in the form of frequency output. In the inertial navigation system, it is necessary to convert the frequency value into an acceleration value and then perform attitude calculation. Using software method for floating point calculation requires a lot of computing power of CPU. In order to optimize the calculation speed of frequency conversion, the paper designs a floating-point calculation acceleration unit for frequency conversion applications, which is implemented and verified based on FPGA. The results show that the system from data sampling to frequency conversion, then converts the frequency value into acceleration for the attitude calculation. The angular velocity measured by the gyroscope is integrated, and finally the data fusion is completed. Using the floating-point acceleration unit designed in the paper to achieve frequency conversion, the speed is improved by 2 times.
Key words : quartz vibrating beam accelerometer;inertial navigation system;frequency conversion;floating point calculation

0 引言

    隨著微機(jī)電系統(tǒng)(MEMS)技術(shù)的發(fā)展,石英振梁加速度計(jì)作為一種新型MEMS慣性傳感器的全數(shù)字脈沖輸出,越來越廣泛應(yīng)用于小型飛機(jī)、機(jī)器人導(dǎo)航系統(tǒng)等[1-2]。為了使該器件能在慣導(dǎo)系統(tǒng)中發(fā)揮出其優(yōu)良性能,需要對其輸出頻率脈沖信號高精度測量的計(jì)數(shù)值進(jìn)行浮點(diǎn)計(jì)算

    然而隨著超大規(guī)模集成電路設(shè)計(jì)技術(shù)的發(fā)展,在慣性導(dǎo)航系統(tǒng)中,大多數(shù)DSP處理器具有雙精度浮點(diǎn)運(yùn)算能力,但在這其中存在大量的運(yùn)算,會耗用CPU大量的計(jì)算能力,并且采用功能級流水增加了面積開銷和硬件復(fù)雜度[3]。通過軟件設(shè)計(jì)實(shí)現(xiàn),不僅占用大量的軟件資源,而且運(yùn)行速度慢,這些設(shè)計(jì)方法不太適合一些高速數(shù)據(jù)采集和處理應(yīng)用場合[4]

    根據(jù)石英振梁加速度計(jì)的特性,在慣性導(dǎo)航系統(tǒng)中,為了優(yōu)化頻率轉(zhuǎn)換的計(jì)算速度,本文提出一種面向頻率轉(zhuǎn)換應(yīng)用的浮點(diǎn)計(jì)算加速單元設(shè)計(jì),對加速度計(jì)測量的計(jì)數(shù)值進(jìn)行浮點(diǎn)計(jì)算,以加快運(yùn)算的換算速度,減輕了導(dǎo)航計(jì)算機(jī)中處理器的負(fù)擔(dān),從而加快導(dǎo)航姿態(tài)解算的速度。采用FPGA對浮點(diǎn)計(jì)算加速單元進(jìn)行了實(shí)驗(yàn)驗(yàn)證,與傳統(tǒng)的軟件計(jì)算方法進(jìn)行了比較,加速度解算速度提高了2倍,可以很好地提升慣性導(dǎo)航系統(tǒng)的性能。

1 總體設(shè)計(jì)

1.1 石英振梁加速度計(jì)原理

    石英振梁加速度計(jì)是一種基于石英振動(dòng)梁的力-頻特性的新型高精度固態(tài)傳感器。其內(nèi)部交變電場使其推挽方式安裝的兩個(gè)石英梁在伸縮振動(dòng)模式下進(jìn)行振動(dòng)。當(dāng)加速度計(jì)受到外部加速時(shí),敏感質(zhì)量塊產(chǎn)生一個(gè)慣性力以分別作用于這兩個(gè)石英梁上。其中一個(gè)石英梁受到壓縮作用,其諧振頻率會降低;而另一個(gè)受到拉伸作用,其諧振頻率將升高。最終這兩個(gè)石英梁的頻率差與外施加力成比例,即與加速度成比例,進(jìn)而測量出加速度值[5]

    因此每個(gè)加速度計(jì)將輸出兩個(gè)頻率信號,這兩個(gè)頻率信號之差與加速度的關(guān)系如下:

    wdz8-gs1.gif

其中L0≤0.9 g,為零偏值;L1為標(biāo)度因數(shù),典型值為50 Hz/g;L2為二階非線性系數(shù),L2≤20 μg/g2;f1為力敏感石英梁F1的輸出,單位為Hz;f2為力敏感石英梁F2的輸出,單位為Hz;a為即時(shí)加速度計(jì),單位為g。

1.2 系統(tǒng)架構(gòu)設(shè)計(jì)

    傳統(tǒng)的慣性導(dǎo)航系統(tǒng)架構(gòu)設(shè)計(jì)如圖1所示,首先系統(tǒng)對三個(gè)石英振梁加速度計(jì)和三個(gè)陀螺儀進(jìn)行數(shù)據(jù)采集。三個(gè)石英振梁加速度計(jì)的輸出信號和溫度信號均為頻率信號,每個(gè)加速度計(jì)輸出2路頻率信號。根據(jù)式(1)可知,要得到加速度,就要測出2個(gè)石英梁的頻率輸出值,則設(shè)計(jì)的頻率采樣系統(tǒng)就必須能對9路通道的頻率信號進(jìn)行高精度測量與快速計(jì)算。系統(tǒng)對加速度計(jì)進(jìn)行高精度測頻,使用的是等精度測頻法,其公式為:

    wdz8-gs2.gif

其中fs為標(biāo)準(zhǔn)頻率,Ns為標(biāo)準(zhǔn)頻率信號脈沖個(gè)數(shù),fx為被測頻率,Nx為被測頻率信號脈沖個(gè)數(shù)。

wdz8-t1.gif

    根據(jù)式(2),經(jīng)過頻率轉(zhuǎn)換計(jì)算得到測頻值,然后由式(1)可知,頻率值轉(zhuǎn)換后進(jìn)而得出加速度值。同時(shí)三個(gè)陀螺儀經(jīng)過脈沖測量后得到角速度,最后將得到的加速度值和陀螺儀測得的角速度值輸入到導(dǎo)航計(jì)算機(jī)進(jìn)行后續(xù)的解算。

    系統(tǒng)設(shè)計(jì)實(shí)現(xiàn)流程圖如圖2所示,加速度計(jì)測量的頻率值經(jīng)過轉(zhuǎn)換得到加速度,再通過姿態(tài)解算可得到姿態(tài)角θα,角速度通過積分算出角度θg,使用互補(bǔ)濾波算法進(jìn)行數(shù)據(jù)融合,可得到最終的姿態(tài)角θ。

wdz8-t2.gif

    而在傳統(tǒng)的系統(tǒng)架構(gòu)設(shè)計(jì)中,頻率轉(zhuǎn)換計(jì)算會耗用大量的時(shí)間,為了加快頻率轉(zhuǎn)換的計(jì)算速度,本文設(shè)計(jì)了一種面向頻率轉(zhuǎn)換應(yīng)用的浮點(diǎn)計(jì)算加速單元。如圖3所示,在加速度計(jì)高精度頻率測量模塊后設(shè)計(jì)了FPU模塊,此模塊主要是在FPGA芯片里面完成的。設(shè)計(jì)時(shí),選用了Xilinx公司的XC3S1400AN芯片。在下一節(jié)將對浮點(diǎn)計(jì)算加速單元的實(shí)現(xiàn)進(jìn)行詳細(xì)介紹。

wdz8-t3.gif

2 加速單元的實(shí)現(xiàn)

    浮點(diǎn)計(jì)算是以9路頻率計(jì)數(shù)值的順序執(zhí)行,即按照X1、X2、X3、Y1、Y2、Y3、Z1、Z2和Z3的順序執(zhí)行,在相應(yīng)的計(jì)數(shù)器值鎖定完成后立即開始浮點(diǎn)計(jì)算。浮點(diǎn)計(jì)算過程如圖4所示,主要包括:乘法模塊、定點(diǎn)浮點(diǎn)轉(zhuǎn)換(規(guī)格化)、浮點(diǎn)除法、結(jié)果鎖存4個(gè)過程,最終計(jì)算結(jié)果為滿足IEEE-754標(biāo)準(zhǔn)[6]的規(guī)格化雙精度浮點(diǎn)數(shù)。 

wdz8-t4.gif

2.1 浮點(diǎn)計(jì)算的實(shí)現(xiàn)

2.1.1 乘法模塊

    一個(gè)浮點(diǎn)計(jì)算過程首先是9路通道的計(jì)數(shù)值并行輸入到乘法模塊,然后按照X1,X2,X3,Y1,Y2,Y3,Z1,Z2和Z3的順序?qū)s與Nx,y,z進(jìn)行乘法計(jì)算。其過程在對應(yīng)的計(jì)數(shù)值鎖存完成后,立即開始浮點(diǎn)計(jì)算。9路通道根據(jù)控制器來選擇其中一個(gè)通道的計(jì)數(shù)值進(jìn)入乘法模塊做乘法計(jì)算。本文是根據(jù)工程所需,選用標(biāo)準(zhǔn)頻率fs為80 MHz,通過分步相加實(shí)現(xiàn)的乘法計(jì)算。

2.1.2 規(guī)格化

    在控制器的控制下,對fs與Nx,y,z相乘得到的結(jié)果進(jìn)行規(guī)格化處理,即先判斷計(jì)算的尾數(shù)結(jié)果是否符合規(guī)格化的格式,如果不符合,則將其進(jìn)行移位處理,并對指數(shù)做相應(yīng)的加減操作。9路通道的計(jì)數(shù)值計(jì)算完成后都要進(jìn)行規(guī)格化處理。

2.1.3 級間寄存器

    對fs與Nx,y,z相乘得出的結(jié)果進(jìn)行計(jì)算前導(dǎo)0數(shù)量,即找出最高1的位置。在控制器的控制下,鎖存到寄存器1和寄存器2中;使用最高1位置值加上指數(shù)偏移值,就是規(guī)格化后的指數(shù)。再根據(jù)最高1的位置,把其后的有效數(shù)字鎖存到浮點(diǎn)尾數(shù)從最高位開始的位置,尾數(shù)低位補(bǔ)0。對于擴(kuò)展標(biāo)志,只有0與非0兩種情況,如果計(jì)數(shù)值為全0,擴(kuò)展標(biāo)志為00,否則為01(表示正常數(shù))。

2.1.4 浮點(diǎn)除法 

    在控制器的控制下,將對級間寄存器中的數(shù)據(jù)進(jìn)行浮點(diǎn)除法運(yùn)算。根據(jù)浮點(diǎn)除法運(yùn)算必需的步驟,設(shè)計(jì)浮點(diǎn)除法整體單元如圖5所示[7]。本設(shè)計(jì)由6個(gè)部分構(gòu)成,即預(yù)處理、尾數(shù)除、指數(shù)減、異常處理、規(guī)格化與舍入、溢出判斷與輸出,其中尾數(shù)除是通過SRT基4算法[8]來實(shí)現(xiàn)的。

wdz8-t5.gif

    浮點(diǎn)計(jì)算是在乘法計(jì)算模塊到浮點(diǎn)除法模塊的過程中使用了流水線結(jié)構(gòu),當(dāng)?shù)谝煌ǖ劳瓿沙朔ㄓ?jì)算后,立即送到下一模塊進(jìn)行相應(yīng)的數(shù)據(jù)處理,進(jìn)而將其送到浮點(diǎn)除法模塊做最終計(jì)算。同時(shí)第二通道在處理其他數(shù)據(jù),以此順序執(zhí)行運(yùn)算。

    本文在浮點(diǎn)除法過程中采用非流水線結(jié)構(gòu),此結(jié)構(gòu)是以串行形式的迭代運(yùn)算實(shí)現(xiàn)的,在每次迭代時(shí)都共用同樣的硬件,產(chǎn)生的部分余數(shù)和部分商均在寄存器中更新。此過程是在9路全部完成最終計(jì)算后,再繼續(xù)下一數(shù)據(jù)處理。

2.1.5 結(jié)果鎖存

    在控制器的控制下,對浮點(diǎn)除法模塊計(jì)算的最終結(jié)果進(jìn)行鎖存,然后浮點(diǎn)輸出。在以石英振梁加速度計(jì)為基礎(chǔ)的慣性導(dǎo)航系統(tǒng)中,此結(jié)果經(jīng)過轉(zhuǎn)換為加速度后進(jìn)行姿態(tài)解算。

2.2 控制器的實(shí)現(xiàn)

    本文利用硬件描述語言VHDL實(shí)現(xiàn)對浮點(diǎn)加速單元進(jìn)行邏輯設(shè)計(jì)。在整個(gè)控制過程中,乘法模塊是由乘法變成加法計(jì)算的,使用7個(gè)加法器就可以實(shí)現(xiàn)乘法計(jì)算,3個(gè)周期即可計(jì)算出結(jié)果。

    在乘法計(jì)算模塊到浮點(diǎn)除法模塊的過程中使用了兩級流水線結(jié)構(gòu)。圖6所示為乘法計(jì)算模塊到浮點(diǎn)除法模塊的工作圖,首先是在控制器的控制選擇通道下,第一個(gè)通道輸入計(jì)數(shù)器數(shù)據(jù)有效后進(jìn)行乘法計(jì)算,等待乘法計(jì)算完成后,對其計(jì)算的結(jié)果進(jìn)行規(guī)格化處理也就是定點(diǎn)浮點(diǎn)轉(zhuǎn)換。兩個(gè)周期后,待規(guī)格化數(shù)據(jù)輸出穩(wěn)定后,將其計(jì)算的結(jié)果鎖存到M1。一旦鎖存完成后,立即啟動(dòng)浮點(diǎn)除法計(jì)算單元,進(jìn)行浮點(diǎn)除法運(yùn)算。等待浮點(diǎn)除法完成后,在控制器的控制下,選擇相應(yīng)通道的地址將其計(jì)算的結(jié)果鎖存到D1,同時(shí)系統(tǒng)在控制器的控制下,選擇第二個(gè)通道數(shù)據(jù)輸入到乘法計(jì)算模塊做乘法計(jì)算,并將規(guī)格化后的計(jì)算結(jié)果鎖存到M2,這時(shí)也將第三個(gè)通道的計(jì)數(shù)值輸入到乘法計(jì)算模塊。整個(gè)浮點(diǎn)計(jì)算過程以這種工作方式來完成9路通道數(shù)據(jù)的全部計(jì)算。

wdz8-t6.gif

    浮點(diǎn)除法模塊是用非流水線結(jié)構(gòu)來實(shí)現(xiàn)的,每次浮點(diǎn)除法計(jì)算需要26個(gè)周期,因此,采用圖4所示的部件級流水結(jié)構(gòu),完成9個(gè)頻率值的浮點(diǎn)計(jì)算共需要260個(gè)周期。

    如果用流水線結(jié)構(gòu)實(shí)現(xiàn)浮點(diǎn)除法運(yùn)算[8],則只需要35個(gè)周期就可以完成9次浮點(diǎn)除法運(yùn)算,但會耗用大量的面積資源。

    本文對這兩種實(shí)現(xiàn)方式都進(jìn)行了實(shí)現(xiàn)與驗(yàn)證,經(jīng)過綜合分析后得出,雖然全流水線結(jié)構(gòu)處理數(shù)據(jù)速度快,但面積開銷大。而部件級流水線結(jié)構(gòu),除法計(jì)算采用非流水結(jié)構(gòu),速度雖然慢,但也能滿足要求,面積開銷以及性能均得到了優(yōu)化折中。

3 結(jié)果與驗(yàn)證

3.1 實(shí)驗(yàn)驗(yàn)證

    本文在慣性導(dǎo)航系統(tǒng)中,對浮點(diǎn)計(jì)算加速單元進(jìn)行程序編寫、綜合及外部端口的輸出狀態(tài)的波形觀測與驗(yàn)證,并在XC3S1400AN芯片上對其進(jìn)行了硬件測試。

    通過驗(yàn)證結(jié)果表明,設(shè)計(jì)的浮點(diǎn)計(jì)算輸出為64位的雙精度浮點(diǎn)數(shù),結(jié)果正確,并且完全符合要求。另外,采用全流水線結(jié)構(gòu)和部件流水結(jié)構(gòu)設(shè)計(jì)的浮點(diǎn)計(jì)算加速單元經(jīng)過綜合、布局布線后,占用資源情況如表1所示。通過兩種方法對比可以看出,全流水線結(jié)構(gòu)實(shí)現(xiàn)浮點(diǎn)計(jì)算,其面積開銷為24%,而部件流水結(jié)構(gòu)的浮點(diǎn)計(jì)算加速單元只需4%的面積開銷,也就是說全流水結(jié)構(gòu)實(shí)現(xiàn)的浮點(diǎn)計(jì)算所消耗的面積資源是部件流水設(shè)計(jì)的6倍。與全流水線結(jié)構(gòu)實(shí)現(xiàn)浮點(diǎn)計(jì)算相比,明顯可以看出本文設(shè)計(jì)的部件流水結(jié)構(gòu)的浮點(diǎn)計(jì)算加速單元極大地減少了面積開銷。

wdz8-b1.gif

3.2 性能驗(yàn)證 

    本文對性能進(jìn)行了測試,系統(tǒng)使用的處理器為TMS320C6713(工作頻率為100 MHz)。首先系統(tǒng)對三個(gè)加速度計(jì)和三個(gè)陀螺儀進(jìn)行數(shù)據(jù)采集,之后對加速度計(jì)做頻率轉(zhuǎn)換計(jì)算,再將頻率值轉(zhuǎn)換成加速度值,經(jīng)過姿態(tài)解算得到姿態(tài)角。同時(shí)三個(gè)陀螺儀再做軟件計(jì)算,經(jīng)過脈沖測量后得到角速度,角速度通過積分得到角度。然后使用互補(bǔ)濾波算法進(jìn)行角度融合,整個(gè)過程運(yùn)行時(shí)間為95.8 μs。而采用本文設(shè)計(jì)的全流水線結(jié)構(gòu)、部件流水結(jié)構(gòu)的浮點(diǎn)計(jì)算加速單元實(shí)現(xiàn)頻率轉(zhuǎn)換計(jì)算,整個(gè)解算過程運(yùn)行時(shí)間均為27.8 μs。與軟件計(jì)算方法相比,運(yùn)行時(shí)間減少了70.9%,速度提高了2倍。

4 結(jié)論

    本文設(shè)計(jì)了一種面向頻率轉(zhuǎn)換應(yīng)用的浮點(diǎn)計(jì)算加速單元,在實(shí)際的導(dǎo)航計(jì)算機(jī)系統(tǒng)中進(jìn)行了實(shí)現(xiàn)驗(yàn)證,功能正確可靠。使用本文設(shè)計(jì)的浮點(diǎn)計(jì)算加速單元,從數(shù)據(jù)采樣到頻率轉(zhuǎn)換,再經(jīng)過頻率值轉(zhuǎn)換成加速度進(jìn)行姿態(tài)解算;陀螺儀測得的角速度再進(jìn)行積分,最后進(jìn)行數(shù)據(jù)融合,速度提高了2倍。相比使用全流水線結(jié)構(gòu),本文設(shè)計(jì)的部件流水結(jié)構(gòu)的浮點(diǎn)計(jì)算加速單元的FPGA面積開銷僅為其1/5。

參考文獻(xiàn)

[1] ZHI W,CHU J,LI J,et al.A novel attitude determination system aided by polarization sensor[J].Sensors,2018(18):158.

[2] 楊挺,楊貴玉,李慶豐.石英振梁加速度計(jì)靜態(tài)輸入輸出特性[J].中國慣性技術(shù)學(xué)報(bào),2014,22(3):386-390.

[3] 苑佳紅.Stratus HLS工具在高性能雙精度浮點(diǎn)乘法設(shè)計(jì)中的應(yīng)用流程[J].電子技術(shù)應(yīng)用,2018,44(8):20-30.

[4] 朱峰,魯征浩,朱青.形式化驗(yàn)證在處理器浮點(diǎn)運(yùn)算單元中的應(yīng)用[J].電子技術(shù)應(yīng)用,2017,43(2):29-32.

[5] 馮麗爽,王文璞,周震,等.石英振梁加速度計(jì)諧振器的結(jié)構(gòu)設(shè)計(jì)[J].中國慣性技術(shù)學(xué)報(bào),2013,21(1):101-105.

[6] ZURAS D,COWLISHAW M,AIKEN A,et al.IEEE standard for floating-point arithmetic[S].IEEE Std 754-2008,2008:1-70.

[7] 陳勛.基于SRT4算法的浮點(diǎn)除法/方根算術(shù)單元設(shè)計(jì)[D].西安:西安電子科技大學(xué),2016.

[8] A VHDL library of parametrisable floating-point and LNS operators for FPGA[DB/OL].(2006-07-31)[2019-01-22].http://www.ens-lyon.fr/LIP/Arenaire/Ware/FPLibrary.



作者信息:

田換換1,2,朱曉燕1

(1.首都師范大學(xué) 信息工程學(xué)院,北京100048;2.電子系統(tǒng)可靠性技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室,北京100048)

此內(nèi)容為AET網(wǎng)站原創(chuàng),未經(jīng)授權(quán)禁止轉(zhuǎn)載。
主站蜘蛛池模板: 狠狠色综合色区| 被公连续侵犯中文字幕| 在线观看成人网站| 中文字字幕码一二区| 日韩欧美亚洲综合久久| 亚洲区中文字幕| 欧美韩国日本在线观看| 女人张腿让男桶免费视频大全| 久久国产乱子伦精品免费强| 欧美一区二区三区久久综合| 亚洲欧美不卡视频| 色婷婷中文字幕| 国产美女被遭强高潮免费网站| jizzjizz成熟丰满舒服| 性欧美视频在线观看| 丰满多毛的陰户视频| 欧美欧美欧美欧美| 人人爽人人爽人人爽人人片av| 麻豆精品传媒一二三区在线视频| 国产精品久久国产精品99| 三级韩国床戏3小时合集| 日本一区二区三区在线视频观看免费| 久久精品国产四虎| 曰批免费视频播放免费| 亚洲中字慕日产2020| 欧美性xxxx极品高清| 亚洲欧美中文日韩在线| 毛片免费vip会员在线看| 亚洲色图第一页| 特级无码毛片免费视频| 免费**毛片在线播放直播| 精品中文字幕在线| 免费观看成年人网站| 精品人妻潮喷久久久又裸又黄| 又粗又大又猛又爽免费视频| 美女视频黄频a免费大全视频| 国产一区二区高清| 色婷婷综合在线| 国产3级在线观看| 羞羞漫画小舞被黄漫免费| 四虎精品在线视频|