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

基于調頻連續波雷達的物體運動狀態實時檢測算法研究*

2021-11-01 06:10:46屈奎張榮福肖鵬程
物理學報 2021年19期
關鍵詞:測量信號

屈奎 張榮福? 肖鵬程

1) (上海理工大學,光電信息與計算機工程學院,上海 200093)

2) (復旦大學,專用集成電路與系統國家重點實驗室,上海 201203)

微波雷達依靠非接觸、響應速度快、對自然環境的適應性強等特點在物體運動狀態檢測中的應用越來越廣泛.常用的調頻連續波雷達運動檢測算法基于差拍信號頻譜的峰值估計,存在計算量大,抗干擾能力差等缺點.本文通過對運動物體的差拍信號做特定頻率的離散傅里葉變換,將變換后的實部和虛部在互相垂直的兩個方向上進行疊加,其合成軌跡近似為橢圓,求出各軌跡點的相位即可還原物體的運動狀態.該算法無需對每個調頻周期的拍信號做頻譜分析,時間復雜度較低.靜止物體的拍信號被處理成了固定的直流信號,對運動物體的測量不造成影響,具有抗靜止物體干擾的能力.在雷達中心頻率為24 GHz,帶寬為0.15 GHz 的條件下對算法進行了驗證,位移測量精度達到0.27 mm,以500 mm 作為位移的測量范圍,線性度達到0.05%.速度的測量精度為1.11 mm/s.

1 引言

位移和速度是表征物體運動狀態的基本物理量,運動狀態的監測可以應用在精密制造、安防、自動駕駛、健康監控等方面.特別是現代社會向智能化、萬物互聯快速發展的過程中,對各種物體如機器人、無人機、液壓桿、傳送帶等的精確控制提出了越來越高的要求,對物體運動狀態的高精度檢測是精確控制的必要前提.

運動狀態檢測技術可以采用接觸和非接觸兩種方式,接觸式是傳感器與被測物有實體上的連接,該方式受到很多限制,某些情況下不宜采用,比如當被測物與監測儀器距離較遠或者被測物對輕量化要求較高.還有些情況無法采用接觸式,比如自動駕駛中對車道上其他車輛的監測.非接觸測量方式較為靈活,對被測物的運動狀態基本沒有影響,可以在被測物無感知的情況下進行監測,常用的方法有光學[1]、超聲[2]、激光[3]和微波雷達等,光學、超聲、激光方法的共同缺點是測量精度受煙、塵、霧、雨等環境影響較大,有時甚至無法測量,微波雷達受環境的影響很小,是一種實際應用價值很強的運動測量方式.

微波雷達常見的調制方式有單頻[4,5]、脈沖[6,7],調頻[8,9]等,其中單頻雷達發射單一頻率的電磁波,通過計算回波的相位可以獲得很高的位移測量精度,但在獲得物體的絕對距離以及多目標測量中存在缺陷,雖然文獻[10,11]報道了對這些缺陷的彌補,但付出了提高硬件和計算復雜度的代價.脈沖雷達要想獲得較高的測距精度要求有較窄的脈沖寬度和極高的時間分辨能力[12],這些都對硬件提出了很高的要求.調頻連續波(Frequency modulated continuous wave,FMCW)雷達是對發射波的頻率進行調制,最簡單的一種調頻方式是線性調頻,也是較為常用的調頻雷達制式.調頻雷達相比單頻雷達其回波攜帶的信息更加豐富,相比脈沖雷達,其信號收發同時工作,不存在距離盲區,且同樣的探測距離其發射功率峰值較低.

使用FMCW 雷達進行精確位移測量,一般核心步驟是使用頻譜的各種估計方法得到待測物體的拍信號對應的精確頻率.大致過程是在每個調頻周期內,對差拍信號的幅度譜通過峰值查找得到目標譜線的粗略位置,再用頻譜細化方法進行高精度修正.齊國清在研究油罐液位的測量中通過對差拍信號做離散傅里葉變換(Discrete Fourier transform,DFT)得到相位,為解決相位隨距離變化的周期性帶來的距離初值的模糊性,將一個調頻帶寬的信號分成兩段,分別求出各段的相位,兩者的差值用來消除距離模糊[13,14];Zwick 研究組[15?17]在對拍信號做DFT 的基礎上,采用插值法,通過離散頻率的振幅比值確定拍頻頻率,該組還嘗試了使用調頻Z 變換(Chirp Z-transform,CZT)進行較為精細的頻率測量;Pohl 研究組[18?20]開發了中心頻率為80 GHz,帶寬達到20 GHz 以上的超寬帶FMCW雷達系統,大的帶寬帶來很高的距離分辨率,但也意味著更復雜的硬件結構.以上研究組針對的測量場景都較為簡單,回避了多徑干擾,把問題歸結為求解單一正弦信號的頻率,在具有多散射物的場景中并不適用.另外,各國以國際通信聯盟制定的ISM(Industrial scientific medical)頻段為基礎,對無線電頻段的使用進行了嚴格的限制,在實際應用中無法通過任意提高帶寬達到提高距離分辨率的目的.

本文所提算法僅對運動物體的差拍信號做特定頻率的DFT,將變換結果中的虛部和實部在互相垂直的兩個方向上疊加,疊加后的軌跡近似為橢圓,求出各軌跡點的相位即可還原物體的運動狀態.在雷達中心頻率為24 GHz,帶寬為0.15 GHz的條件下對理論進行了實驗驗證.

2 FMCW 雷達的差拍信號

FMCW 雷達在一個調頻周期T內發射信號的波形函數為

式中,f0為起始頻率,Bw為調頻帶寬,Φ0為初相.發射波的頻率隨時間的變化關系如圖1 所示.

圖1 發射波的頻率隨時間的變化Fig.1.the time-frequency domain of transmitted signal.

常用FMCW 雷達收發天線的間距遠小于與被測物的間距,因此可將收發天線看成在同一個位置.與天線相距為R的物體,其回波信號相比發射信號時間延遲了td2R/C,C是電磁波在介質中的傳播速度,其回波函數為

回波和發射波的混頻信號為

由于td?T,可以把時間延遲量td的平方項部分省略,得到:

(4)式即是與天線相距為R的物體對應的差拍信號,式中波函數的振幅均做了歸一化處理.

3 算法原理

調頻周期T一般為毫秒甚至微秒級,考慮這樣一種情況,在T的時間內物體的位移很小,以致可以忽略,此情況下,可以認為物體與雷達距離R不變,其差拍信號用(4)式表示.

對一個調頻周期內的差拍信號做特定頻率的N點采樣DFT,表達式為

僅從(5)式很難看出結果的規律性,(6)式是其積分形式,采樣數越密集,(5)式與(6)式的結果越接近.(5)式結果的規律性可以由(6)式的結果來反應:

將(4)式代入(6)式,可得:

R ek和 I mk分別是F(k) 的實部和虛部.令λC/(f0+Bw/2),λ等于發射波的中心波長.(7)式、(8)式可記為

(9)式、(10)式結構相似,均明顯可分為三個因子相乘,即它們的第三項隨R作周期性變化,變化周期為λ/2.當物體的位移不大于第三項的周期,即|?R|≤λ/2 時,如果第一項和第二項乘積的相對變化很小,在該運動過程中就可以近似用一個常數替代它們.

將(11)式分別代入到 R ek和 I mk的前兩項乘積中,可得:

當|?R|≤λ/2 時,所引起的(12)式和(13)式的相對變化量滿足如下關系:

通過分析(14)式和(15)式右邊表達式的變化規律,可以知道,當 ?Rλ/2,k1,δ?0.5 時,(14)式右邊的值最大; ?Rλ/2,k1,δ0.5時,(15)式右邊的值最大.

以在實際應用中常見的24 GHz 雷達系統為例來觀察(14)式、(15)式的具體值.具體參數可以取λ1.25 cm,Bw0.15 GHz,C3×108m/s,估計出(14)式和(15)式的具體值為

(16)式表示在給定波長和帶寬條件下,?R在λ/2 的范圍內變化所引起的,的相對變化量很小,如果把它們看成常數,則(9)式、(10)式可以近似寫為

實際測量環境中,可能會有靜止物體、雜波以及其他運動物體的回波干擾.簡單起見,可以只考慮靜止物體、具有固定偏移量的干擾及頻率穩定的雜波,對于這些干擾的拍信號做特定頻率的DFT后是一個直流信號,相當于在(17)式、(18)式的基礎上增加一個固定的直流分量 R e0和 I m0.

顯然,合成軌跡是一個橢圓,標準形式為

通過橢圓擬合算法,可以得到橢圓參數,也即得到 R e0,I m0,Ak,Bk的值.

為方便描述,以下將這種拍信號在特定頻率下DFT 分量垂直合成軌跡的橢圓化算法(Ellipse algorithm for the vertical synthesis trajectory of the DFT component of the beat signal at a specific frequency)簡記為ETBF 算法.

考慮僅有單個運動物體的情況,特定離散頻率系數k的確定分兩種情況討論.第一種,待測物的初始位置R0已知,則根據k的選取原則,kround(2R0Bw/c),符號round()表示對括號里面的量做四舍五入.初始的k值確定后,可以根據待測物體的位移R′對k值進行更新.更新公式為

第二種情況,待測物初始位置未知.對第一個調頻周期內的拍信號做幅度譜的峰值檢測,在這些峰值對應的離散頻率下將拍信號用ETBF 算法處理,靜止物體的信號被ETBF 算法處理后是一個固定值,運動物體的信號在ETBF 算法處理下會形成橢圓化的合成軌跡,而且當離散頻率最接近被測物的差拍信號頻率時,由于此時(6)式積分中的兩個乘積項具有最強的相關性,得到的橢圓幅度將會最大,所以幅度最大的橢圓對應的離散頻率序數即為初始的k值,再用頻譜細化方法得到較為精確的頻率finitial,根據拍信號頻率和距離的關系可以得到初始位置R0CTfinitial/(2Bw).同樣,采用(22)式對k值進行更新.此情況僅在初始位置的確定中使用了常用的頻譜估計算法.

4 相位解纏(phase unwrapping)算法和運動狀態的計算

得到擬合橢圓的參數后,即可通過(19)式和(20)兩式計算各軌跡點對應相位φ,即,易知:

根據(24)式計算反正切得到的相位φ′被限制在[ 0,2π).

分析(19)式、(20)式可知,當R增大時,Rek和 I mk的合成軌跡點做逆時針運動,反之做順時針運動.相位φ與R成正比,所以,當合成軌跡點做逆時針運動,意味著相位增大,反之減小.考察軌跡點逆時針越過圖2 所示的橫坐標軸的正半軸時的情況,此時實際相位增大,但直接根據(24)式計算出的相位減小;當軌跡點做順時針運動越過橫坐標的正半軸時,實際相位減小,但根據(24)式計算出的相位增大.所以,直接根據(24)式的計算結果無法反應真實的相位變化,需要通過相位解纏算法解決這一問題.而且由于初始時刻的相位值無法獲得,最終得到的是相對初始時刻的相對相位.

圖2 變量R ek和I mk的合成軌跡示意圖.紅色實心點處在橫坐標軸的正半軸上,代表相位φ=2jπ,j∈Z 的位置Fig.2.The synthetic ellipse trajectory diagram of variables Rek and I mk.The red solid point is on the positive half axis of abscissa,which represents the position of φ=2jπ,j∈Z.

將第n個由(24)式計算出的相位采樣值記為φn,解纏后的相對相位記為,則:

jn可以稱為相對相位級數,n=1 時,j10,φ1是n=1 時的φn.

解纏繞的過程就是確定jn的過程.定義 ?φnφn ?φn?1,根據?φn的值可以判斷φn相對φn?1是否越過了橫坐標的正半軸.

1) ?π

后一采樣點相對于前一采樣點在橢圓軌跡上沒有越過橫坐標的正半軸,當前相對相位級數不變

2) π

此情況代表軌跡點做順時針運動且越過了橫坐標的正半軸,相位減小了 2π.所以當前相對相位級數

3) ?2π

此情況代表軌跡點做逆時針運動且越過了橫坐標的正半軸,相位增大了 2π.所以當前相對相位級數

由相對相位可求出其對應的位移為

位移的符號上加了撇號以示和絕對距離相區分.

當位移的測量比較準確時,可以直接對時間求導得到速度

對于離散值,(30)式的數值計算表達式為

實際是用 ?t時間內的平均速度代替瞬時速度,具體計算時,可根據實際情況選取合適的 ?t值.

5 與其他算法的比較

時間復雜度是對算法的消耗時間進行度量,由于算法消耗時間與工作量成正比,所以也等價于對工作量進行度量.通過大O符號表示法[23]給出計算工作量的漸進函數,可以觀察不同算法在計算工作量上的本質差異.

采用表1 中其他文獻的算法時,雜波會對待測物的測量造成干擾,雜波越強烈,干擾越大,測量誤差越大,如果雜波強度接近或超過待測物的信號強度時,甚至導致待測物無法檢出,所以文獻采用的測量場景都較為簡單,回避了其他散射物的雜波干擾,實際測量環境很難滿足這種嚴格的條件,算法的適應性受到很大的限制.

表1 與其他算法的比較Table 1.Comparison of this work with other methods.

本文所提出的算法回避了頻譜的峰值估計,只需計算一個特定頻率的DFT,時間復雜度顯著低于其他文獻中的算法.由于靜止物體的拍信號被處理成了固定的直流信號,對運動物體的測量不造成影響,所以具有抗靜止物體干擾的能力.由于計算的是相對相位,限制了測量內容為相對運動.

6 實 驗

6.1 硬件系統

得益于集成電路制造技術的飛速發展,雷達的關鍵部件已實現芯片化,大大簡化了小型雷達系統的設計和制造難度.實驗所用雷達系統主要由ADI 公司的ADF4158 芯片和Infineon 公司的BGT24MTR 芯片組成.ADF4158 芯片產生調制信號,控制BGT24MTR 芯片中的壓控振蕩器產生24 GHz 調頻波形,信號經過功率放大器(Power amplifier,PA)放大后,采用微帶天線發送和接收反射回來的信號.接收到的信號經過低噪放大器(Low noise amplifier,LNA)放大濾波后進入混頻器與發射信號的一部分進行混頻,產生的差拍信號由NI 公司的PXI-4461 采集卡進行數據采集,然后輸入到計算機通過MATLAB 進行信號處理.硬件采用了單通道輸入和輸出.雷達前端如圖3 所示,實驗系統的結構如圖4 所示.

圖3 雷達前端 (a)正面為收發天線;(b)反面為電路板Fig.3.Radar front end:(a) The front side is transceiver antenna;(b) the reverse side is circuit board.

圖4 實驗系統結構圖Fig.4.Block diagram of the FMCW radar system.

實驗中調頻信號的參數為f024 GHz,Bw0.15 GHz,T=4 ms.

在校準過的步進電機上固定一塊金屬板作為反射板,步進電機的運動方向垂直雷達的天線面,裝置如圖5 所示.

圖5 測量裝置與載有金屬板的步進電機Fig.5.Measuring set and stepper motor with metal plate.

6.2 R ek 和 I mk 采樣點的合成運動軌跡

由理論公式(19)和(20)可知,λ/2 的運動距離剛好可以使 R ek和 I mk的合成運動軌跡形成一個完整的橢圓.圖6 給出了用計算機仿真的待測物在距離天線800 mm,1200 mm,1600 mm,2000 mm四個起始位置上運動λ/2 時所形成的運動軌跡.仿真采用的調頻信號參數和實驗相同,將待測物拍信號的振幅設置為單位1,在4 個起始距離上設置拍信號振幅為0.5 的干擾物體,并添加功率為–20 dBW 的高斯白噪聲.圖中星號代表采樣點,實線是根據這些采樣點擬合出的橢圓,可以看出,這些采樣點形成的軌跡和擬合橢圓符合得很好.

圖6 計算機模擬在4 種不同起始距離下做 λ/2 位移時的采樣點軌跡變化,起始距離分別為 (a) 800 mm;(b) 1200 mm;(c) 1600 mm;(d) 2000 mmFig.6.Computer simulation of trajectory change by these sampling points at four different starting distances,the displacement is λ/2,the starting distances are:(a) 800 mm;(b) 1200 mm;(c) 1600 mm;(d) 2000 mm.

控制步進電機帶動金屬板分別在距離天線800 mm,1200 mm,1600 mm,2000 mm 的距離上運動λ/2.圖7 給出了實驗測量出的采樣點和根據這些采樣點擬合出的橢圓,采樣點基本落在擬合橢圓上,仿真和實驗結果都與理論分析一致,說明把運動軌跡看成橢圓具有合理性.由于仿真和實驗中拍信號的振幅并不一致,所以圖6 和圖7 的坐標范圍不一致.

圖7 實驗測量在4 種不同起始距離下做 λ/2 位移時的采樣點軌跡變化,起始距離分別為 (a) 800 mm;(b) 1200 mm;(c) 1600 mm;(d) 2000 mmFig.7.Experimental measurement of trajectory changes by these sampling points at four different starting distances,the displacement is λ/2,the starting distances are:(a) 800 mm;(b) 1200 mm;(c) 1600 mm;(d) 2000 mm.

6.3 位移的測量精度

讓金屬板初始時刻距離雷達天線1500 mm,控制步進電機帶動金屬板以10 mm/s,20 mm/s,30 mm/s 三種速度分別勻速運動300 mm,400 mm,500 mm.重復測量20 次,表2 給出了測量的統計結果.平均值的偏差處于(–0.1 mm,0.1 mm)之間,測量的最大偏差為0.27 mm.測量結果的標準差小于0.15 mm.

表2 金屬板運動的測量結果Table 2.Measurement results of metal plate movement.

亞毫米精度的運動測量一般應用在機器人室內定位、機械臂控制、倒車監測等近距探測場景中,范圍一般在5 m 以內,在這樣的距離范圍內,實驗測量未發現有明顯的精度變化,說明算法在近距測量場景中具有較好的精度穩定性.當測量距離進一步增大后,由于待測物的散射立體角減小而導致的信號強度降低越來越明顯,測量精度會逐步下降.

6.4 位移的測量線性度

線性度可以反應系統對待測運動過程還原的準確性.從金屬板以20 mm/s 的速度移動500 mm的20 組實驗數據中任意選取一組,位移隨時間變化的測量結果如圖8(a)所示.勻速運動時,時間與位移應該是線性關系,圖中同時給出了用最小二乘法擬合出的直線作為線性比較的標準.圖8(a)中的子圖是主圖的部分放大圖,從中可看出測量值和擬合直線之間的輕微偏差.圖8(b)為測量偏差變化圖,偏差的平均值為 2.7×10?13mm,標準差為0.22 mm,基本滿足零均值的高斯分布.最大偏差為0.25 mm,以500 mm 作為測量范圍,線性度達到0.05%.

圖8 位移測量結果 (a)位移隨時間的變化;(b)位移偏差隨位移的變化Fig.8.Measurement of displacement of metal plate:(a) Change of displacement with time;(b) change of displacement deviation.

6.5 速度的測量精度

為檢驗速度測量的精確度,設計了如下的運動過程:首先從靜止開始以 1 0 mm/s2的加速度運動3 s,此時速度達到30 mm/s,然后勻速運動8 s,再以 ? 10 mm/s2的加速度運動3 s,此時速度為零,靜止5 s,然后以 ? 5 mm/s2加速度運動6 s.整個過程可分為5 個階段,分別為勻加速、勻速、勻減速、靜止及一個反向勻加速階段.通過對步進電機編程讓其帶動金屬板按照設定運動.

由于采用(31)式計算速度,需要先得到運動過程的位移.圖9(a)同時給出了位移的設定值和實際測量值,兩者的曲線變化幾乎完全一致.圖9(b)是測量值與設定值的差值變化,平均偏差為0.02 mm,最大偏差δR.max為0.29 mm,比線性運動的偏差稍大.偏差較大的部分出現在加速過程中,可能是步進電機的加速不穩定導致.

圖9 (a)測量位移和設定位移的比較;(b)測量位移的偏差變化Fig.9.(a) Change of measured displacement and set displacement with time;(b) deviation of measured displacement.

測量速度和設定速度隨時間的變化如圖10(a)所示.其中計算平均速度所選取的時間間隔 ?t為0.3 s.根據誤差的傳遞公式,測量速度的理論偏差σ應該滿足.圖10(b)是測量值與設定值的差值變化,最大偏差為1.11 mm/s,標準差為0.31 mm/s,符合誤差分析結果.通過適當增加時間間隔可以提高速度的測量精度,但會帶來速度更新率的下降,實際應用中,可根據具體情況靈活設置.

圖10 (a)速度測量值和速度設定值的比較;(b)速度測量值的偏差變化Fig.10.(a) Change of measured speed and set speed with time;(b) deviation of measured speed.

7 總結

本文所提出的算法可以在單通道接收和輸出的硬件條件下實現,單通道的輸出無需考慮正交解調輸出時的信號失衡問題,且硬件結構更簡單.對運動物體的差拍信號做特定頻率的離散傅里葉變換,避免了使用復雜方法進行頻率估算,減少了計算量,易于實時測量.測量中可以隔離其他靜止物體的干擾,直接提取運動物體的信號,降低了對環境的要求.由于算法得到是相對相位,限制了測量內容為相對運動.在一些需要測量相對位移,比如機械的振動頻率,生命信號探測,機械手操控等領域有很大的應用潛力.

在雷達中心頻率為24 GHz,帶寬為0.15 GHz的條件下對算法進行了驗證.位移測量精度達到0.27 mm,以500 mm 作為位移的測量范圍,線性度達到0.05%.速度的測量精度為1.11 mm/s.實驗采用了符合ISM 要求的頻段,更符合實際應用場景,但中心頻率和帶寬低于文獻[16,19?22]中的設置,高的頻率和帶寬可以顯著提高位移的測量精度,由于硬件條件不一致,不能把本實驗的測量精度簡單的與上述文獻比較.研發出更高頻段的雷達硬件對本算法進行驗證及實現絕對距離的測量是下一步要做的工作.

猜你喜歡
測量信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
把握四個“三” 測量變簡單
滑動摩擦力的測量和計算
孩子停止長個的信號
滑動摩擦力的測量與計算
測量的樂趣
測量
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
主站蜘蛛池模板: 国产福利免费观看| 午夜国产大片免费观看| 亚洲女人在线| 成人伊人色一区二区三区| 成人在线欧美| 日本人妻一区二区三区不卡影院 | 久久国产精品波多野结衣| 国产玖玖视频| 亚洲精品麻豆| 免费啪啪网址| 一本久道热中字伊人| 老司国产精品视频| 国产精品欧美在线观看| 波多野吉衣一区二区三区av| 91精品国产自产在线老师啪l| 99国产在线视频| 国产精品毛片一区视频播| 久青草网站| 亚洲乱码精品久久久久..| 成人免费一区二区三区| 国产高清又黄又嫩的免费视频网站| 日本欧美一二三区色视频| 国产夜色视频| 欧美日韩国产在线观看一区二区三区| 免费一级毛片在线观看| 亚洲无码四虎黄色网站| 毛片大全免费观看| 久久狠狠色噜噜狠狠狠狠97视色| 欧美成人午夜在线全部免费| 最新国产高清在线| 大陆国产精品视频| aaa国产一级毛片| 亚洲国产精品无码久久一线| 精品撒尿视频一区二区三区| 亚洲福利视频一区二区| 国产成人a在线观看视频| 中文字幕一区二区人妻电影| 午夜久久影院| 亚洲丝袜第一页| 日本免费高清一区| 日韩经典精品无码一区二区| 亚洲日韩每日更新| 中文字幕久久亚洲一区| 亚洲欧洲日产国码无码av喷潮| 欧美第九页| 亚洲区一区| 一区二区日韩国产精久久| 国产在线拍偷自揄拍精品| 免费高清自慰一区二区三区| 伊人国产无码高清视频| 久久精品这里只有国产中文精品| 国产特一级毛片| 国产精品无码作爱| 国产成人综合久久| 国产一级做美女做受视频| 91成人在线免费观看| 国产va免费精品观看| 国产精品久久国产精麻豆99网站| 亚洲狠狠婷婷综合久久久久| 国产精品免费福利久久播放| 色综合久久88色综合天天提莫| 成人自拍视频在线观看| 亚洲国产精品VA在线看黑人| 91探花国产综合在线精品| 五月婷婷丁香综合| 欧亚日韩Av| 无码乱人伦一区二区亚洲一| 国产成人1024精品下载| 亚洲三级色| 乱人伦视频中文字幕在线| 91九色视频网| 国产精品一区在线麻豆| 亚洲a级在线观看| 日韩第一页在线| 久久精品无码国产一区二区三区| 97久久免费视频| 亚洲色图另类| 91麻豆精品国产高清在线| 国产精品9| 精品视频福利| 波多野衣结在线精品二区| 欧美成人二区|