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

基于單軸旋轉INS/GPS組合姿態誤差觀測的垂線偏差測量方法

2015-06-05 14:51:33戴東凱王省書戰德軍黃宗升
中國慣性技術學報 2015年2期
關鍵詞:測量方法

戴東凱,王省書,戰德軍,吳 偉,黃宗升

(國防科學技術大學 光電科學與工程學院,長沙 410073)

基于單軸旋轉INS/GPS組合姿態誤差觀測的垂線偏差測量方法

戴東凱,王省書,戰德軍,吳 偉,黃宗升

(國防科學技術大學 光電科學與工程學院,長沙 410073)

單軸旋轉INS/GPS組合導航系統的姿態誤差直接受垂線偏差的影響,因此利用對單軸旋轉INS/GPS組合的姿態誤差觀測也可實現垂線偏差的估計。首先,利用INS中的三個激光陀螺構建了激光陀螺組合體(LGU)并進行自主姿態測量,以LGU作為姿態基準以獲取INS/GPS組合的姿態誤差。然后,建立垂線偏差測量的觀測方程和狀態方程。最終采用Kalman濾波/平滑算法同時實現垂線偏差和其他系統誤差的最優估計。通過對狀態變量精確、合理地建模,并利用全球重力模型補償垂線偏差信號的低頻分量,從而實現垂線偏差與系統誤差的解耦。通過仿真驗證了該方法的可行性, 仿真所用的航跡由實測數據生成。仿真結果表明該方法能夠有效地測量垂線偏差的高頻擾動量。由于該方法的測量精度依賴于所采用的陀螺性能,采用角隨機游走較小的陀螺可以獲得較好的垂線偏差測量結果。船載實驗結果表明,該方法測量得到的垂線偏差數據重復性精度優于0.5″。

垂線偏差;單軸旋轉INS;GPS;激光陀螺組合體

垂線偏差是指地球上某點的真實重力矢量與正常重力矢量在方向上的偏差[1]。垂線偏差含有豐富的重力場高頻信息,相比重力異常,它更能反映重力場的精細結構[2]。在資源勘探、地震和火山監測、衛星精密軌道、重力輔助導航以及高精度慣性導航等應用中,對高精度、高分辨率垂線偏差信息有著極為迫切的需求。

傳統的垂線偏差獲取方法主要有兩種:一是天文大地測量方法,即利用天文經緯儀或數字天頂相機觀測恒星以確定垂線偏差[3]。該方法雖然精度極高,但是操作繁瑣,效率較低,且只能在陸地使用,不適用于海上大動態、大范圍的垂線偏差測量活動。二是重力異常計算方法。該方法是利用衛星重力儀或標量航空重力儀測定區域內的重力異常數據,再通過Vening-Meinesz積分公式[1]計算垂線偏差。其缺點在于所需的高質量重力異常數據較多,數據利用率不高,所獲取的垂線偏差精度和分辨率也較為有限[4]。

基于INS/GPS組合的垂線偏差動態方法是獲取高精度、高分辨率垂線偏差數據的有效手段。該測量方法可以彌補重力異常計算方法精度和分辨率上的不足,克服地面靜態測量方法效率低下的缺陷,是測量精度與測量效率相平衡的一種測量方案。目前,國際上利用INS/GPS組合實現垂線偏差動態測量的主流技術方案是“直接求差法”[5-7]。其基本原理是:利用慣性傳感器或慣性導航系統直接敏感比力信息,同時利用GPS測量干擾加速度,將比力與干擾加速度做差即可得到重力矢量,進而計算得到垂線偏差。該方法面臨兩大難題:一是扣除干擾加速度的影響。20世紀90年代初,載波相位差分GPS技術的發展成功地解決了測量載體運動速度和加速度高精度測量的問題,干擾加速度對慣性測量的影響基本可以消除。二是需要提供高精度的水平基準。高精度水平姿態基準的獲取是制約垂線偏差測量精度的關鍵問題,垂線偏差的存在會引起INS/GPS組合導航系統的姿態誤差,而INS/GPS姿態誤差也會導致相應的垂線偏差測量誤差,因此垂線偏差與慣性姿態基準無法解耦[8]。

INS/GPS組合導航姿態誤差直接受垂線偏差的影響[9],基于此現象文獻[10]提出利用星敏感器和激光陀螺組合體(Laser Gyroscope Unit,LGU)提供姿態基準以測量INS/GPS組合的姿態誤差,進而獲取垂線偏差的方法。該方法能夠有效地解決傳統方法中垂線偏差與組合導航姿態誤差耦合的問題,然而LGU中存在的由初始姿態誤差和陀螺誤差引起的姿態誤差仍無法完全消除。此外,由于該方法需要使用星敏感器,測量只能在晴朗的夜晚進行,且要求載體的運動平穩,嚴重限制了該方法的使用。本文旨在對文獻[10]的方法進行改進,提出一種基于INS/GPS組合姿態誤差觀測的垂線偏差測量新方法。與文獻[10]中的方法相比,新方法有兩方面的改進:一是僅利用LGU構建姿態基準,無需利用星敏感器,因而拓展了測量的適用條件,降低系統的成本和復雜度;二是對LGU姿態誤差進行精確的建模,并利用Kalman濾波實現垂線偏差與LGU姿態誤差的分離。

本文從理論上建立了所提出的垂線偏差測量方法的狀態方程和觀測方程,給出垂線偏差測量的實現方案。通過仿真驗證算法的可行性,研究該測量方法對慣性器件精度的需求。最后通過船載實驗對垂線偏差測量的重復性進行驗證。

1 基于INS/GPS姿態誤差觀測的垂線偏差測量方法

首先定義東北天(E-N-U)坐標系為導航坐標系,即n系。在該坐標系中,原點O位于測量位置,z軸沿著參考橢球的法線指向橢球外部,也就是與O點處正常重力矢量的方向相反(在近地面處忽略正常重力線的彎曲效應);x軸指向東向;y軸指向北向,并與其他兩軸構成正交右手坐標系。

圖1 坐標系與垂線偏差的定義Fig.1 Definitions of the coordinates and DOVs

東向和北向重力擾動分別記為δgE(東向為正方向)和δgN(北向為正方向)。垂線偏差的北-南和東-西角度分量分別記為ξ和η,他們直接與重力擾動矢量的水平分量相關,并定義為:

式中,g是正常重力的大小。在小角度近似條件下, 式(1)和(2) 可近似寫為:

文獻[11]從理論和仿真分析了垂線偏差對INS/ GPS組合導航姿態誤差的影響,分析結果表明垂線偏差與INS/GPS組合導航姿態誤差存在以下關系:

從以上關系可知,通過觀測INS/GPS的姿態誤差即可實現垂線偏差的測量。

1.1 觀測方程的建立

為了獲取INS/GPS的姿態誤差,需要提供一個高精度的姿態基準。文獻[10]提出用INS中的三個激光陀螺構建激光陀螺組合體(Laser Gyroscope Unit,LGU)。在給定初始姿態值后,LGU可以自主地測量本體坐標系(b系)相對于地球坐標系(e系)的姿態[12],利用GPS提供的位置信息可以進一步獲取LGU相對于n系的姿態,以LGU的姿態輸出作為姿態基準獲取INS/GPS組合的姿態誤差。獲取INS/GPS姿態誤差觀測量的流程圖如圖2所示。

圖2 INS/GPS姿態誤差觀測量的獲取流程Fig.2 Procedure for obtaining the observation of INS/GPS integration error

為了便于實施,LGU姿態初始值由INS/GPS組合導航提供。INS/GPS的姿態輸出存在誤差,將引入到LGU的初值中。在LGU姿態解算中,初始誤差和陀螺誤差將引起較大的姿態誤差,因此利用INS/GPS組合與LGU的姿態差中包含了INS/GPS姿態誤差信息和LGU姿態誤差信息。為了獲取垂線偏差信號,關鍵在于如何有效地消除LGU姿態誤差的影響。

INS/GPS與LGU的姿態差可以由下式給出:

式中:ΔΘ為INS/GPS與LGU的姿態差,nδΦ為INS/GPS在n系下的姿態誤差,nδψ為LGU在n系下的姿態誤差。其東向和北向分量可以表示為:

式中:ΔNΘ和ΔEΘ為INS/GPS與LGU姿態差的北向和東向分量;δNφ和δEφ為INS/GPS姿態誤差的北向和東向分量; δNψ和δEψ為LGU姿態誤差 的北向和東向分量。

將式(5)(6)代入式(8)(9)可得:

式中:ΔNΘ和ΔEΘ可以直接利用INS/GPS和LGU的姿態輸出計算得到。式(10)(11)即為垂線偏差估計的觀測方程。

為了利用INS/GPS和LGU姿態差的觀測信息實現垂線偏差的最優估計,下面將對垂線偏差信號和LGU姿態誤差進行狀態空間建模。

1.2 狀態方程的建立

LGU姿態的誤差方程[12]在地球坐標系下可以描述為:

式中:δψe為LGU在e系下的姿態誤差;εb為陀螺在b系下的零偏,為陀螺在b系下的高斯白噪聲。陀螺零偏可以用隨機常值過程建模:

由于GPS可以提供高精度的位置數據,因此LGU相對于e系的姿態可以轉換到n系下。LGU在n系下的姿態誤差方程可以相應的表示為:

由式(14)可以看出,LGU姿態誤差主要是由陀螺誤差和LGU姿態誤差的初始值引起的。由式(14)不難發現,LGU姿態誤差受到地球轉動角速度矢量的調制,而地球的旋轉周期為24 h,因此LGU的姿態誤差將表現為低頻誤差的特性。另一方面,垂線偏差在長波(低頻)區域內也具有較強的功率譜,因此其長波分量很容易與LGU的姿態誤差耦合。

為了避免對垂線偏差進行統一建模,我們將垂線偏差分解為如下兩個部分:

將式(15)和(16)代入式(10)和(11),則式(10)和(11)可以重新整理后得到:

式(17)和(18)為新的觀測方程。

為了使LGU姿態誤差和與垂線偏差擾動不產生耦合,要求對垂線偏差擾動建模的隨機過程其功率譜在低頻區域內有極低的增益,以使低頻系統誤差不會引入到垂線偏差擾動估計結果中。

2階Gauss-Markov過程常用于垂線偏差擾動的建模[15]。然而,該模型在低頻區仍然有一定的增益。本文采用2階Gauss-Markov過程的導數[17]對垂線偏差擾動進行建模。與2階Gauss-Markov模型相比,本文所提出的模型在低頻區域內功率譜密度有更強烈的衰減,因而能實現垂線偏差與系統誤差在頻域上的解耦。δη和δξ的模型可以寫成狀態空間的形式如下:

式中:xE和xN為中間變量,其導數即為垂線偏差擾動;0ω為固有頻率,它與載體的航速有關;?為阻尼系數;qE和qN為相應模型的過程噪聲。

用于垂線偏差擾動估計所選用的狀態空間矢量如式(21)所示,其狀態空間模型可由式(13)(14)(19)以及式(20)給定。

最后,δη和δξ可以由以上給出的狀態空間模型和觀測方程,利用Kalman濾波/平滑算法估計得到。將δη、δξ與、合并即為最終的垂線偏差測量值。

2 仿真與分析

2.1 仿真條件

本節將通過仿真驗證本文方法的可行性。仿真采用導航級的慣性測量單元,且只考慮常值零偏誤差和高斯白噪聲誤差。慣性器件和GPS的誤差參數如表1所示。本文的INS采用的是單軸旋轉調制結構[18],以消除由于慣性器件零偏緩慢變化引起的系統導航誤差。在組合導航中,通過旋轉調制解構可以提高慣性器件誤差的可觀測性,進而抑制測量誤差,提高測量精度。

仿真中測量船的姿態和速度由海上試驗的實測數據給出,利用實測的姿態、速度數據生成用于仿真的陀螺、加表數據。測量船的船姿、船速、船位變化分別如圖3~5所示。

圖3 仿真航跡中的船體姿態Fig.3 Ship attitude on the trajectory in simulation

圖4 仿真航跡中的船體Fig.4 Ship velocity on the trajectory in simulation

圖5 船體的航行軌跡Fig.5 Trajectory of the ship

用于生成仿真數據的重力擾動矢量是從美國公布的其本土的高精度垂線偏差模型DEFLEC12中得到的。該模型的垂線偏差精度和分辨率均高于EGM2008重力模型,能夠更精確地反映重力場的精細結構。如圖4所示,為了能夠利用DEFLEC12模型,仿真時將船位的初始點平移至經度240°,緯度34°的位置。EGM2008相對于DEFLEC12模型的差值如圖6所示,可以看到在某些區域內(山地),EGM2008中包含較大的擾動誤差。

圖6 航跡上EGM2008垂線偏差誤差的分布Fig.6 Distribution of EGM2008 DOV errors along-track

2.2 仿真結果與誤差分析

利用對INS/GPS的姿態誤差的觀測實現垂線偏差測量的可行性是以垂線偏差與INS/GPS組合姿態誤差滿足式(5)和式(6)的關系為基礎的,即δφE-ξ=0,δφN+η=0。圖7給出了仿真得到的δφE-ξ和δφN+η隨時間的變化,可以看到,其數值均在0附近振蕩,且振蕩的幅值較小,因此INS/GPS姿態誤差能夠較好地反映垂線偏差擾動的變化。然而,在垂線偏差變化較為劇烈的區域(第10~14 h),INS/GPS姿態誤差對垂線偏差擾動的跟蹤誤差較大。

圖7 INS/GPS姿態誤差與垂線偏差擾動的關系Fig.7 Relationship between the attitude errors of INS/GPS and DOV disturbances

圖8為LGU姿態基準的誤差變化,由于初始姿態誤差和陀螺誤差的影響,LGU的姿態存在較大的誤差,并表現為低頻趨勢的形式。

圖8 LGU姿態基準誤差Fig.8 Errors of attitude reference provided by LGU

圖9a 垂線偏差擾動的測量結果Fig.9a DOV disturbance measurement result

圖9b 垂線偏差擾動的測量結果Fig.9b DOV disturbance measurement result

利用本文所述的方法對垂線偏差進行估計,最終的垂線偏差測量結果如圖9a、9b所示(前6 h的數據用于INS的精對準,不用于垂線偏差測量)。可以看到,測量得到的垂線偏差擾動為垂線偏差的高頻擾動信息,它能夠更好地反應重力場的細節信息。在本文方法中,垂線偏差的低頻部分由EGM2008提供,在測量時不對它進行估計,然而由于EGM2008在低頻區域內仍然存在一定的誤差,因此這一誤差將直接引入到最終測量結果中。

由圖8可以看到,LGU姿態誤差主要由陀螺零偏誤差和初始姿態誤差引起的,并表現為低頻誤差的特性,因此在頻域上與垂線偏差擾動幾乎沒有混疊。然而,陀螺的角度隨機游走(Angular random walk, ARW)引起的LGU姿態誤差中仍然存在一些高頻成分,這部分高頻擾動誤差有可能會引入到垂線偏差擾動的測量結果中。為了進一步驗證方法的可行性,本文仿真了不同陀螺角度隨機游走條件下的垂線偏差測量結果,仿真結果如圖10所示。由圖10可以看出,隨著陀螺ARW系數的增大,垂線偏差測量的誤差也會相應的增大。因此,本文的垂線偏差測量方法要求所采用的慣導系統具有較低的陀螺角隨機游走誤差。

圖10 垂線偏差測量誤差與角隨機游走的關系Fig.10 DOV measurement error with respect to ARW

3 實驗結果

圖11 測量船的航跡Fig.11 Trajectory of the survey ship

為驗證本文方法的有效性,課題組于2014年在中國近海進行了船載測量實驗。如圖11所示,在測量的海域內有兩條重復的測線,我們將其命名為測線1 (SL1) 和測線2 (SL2)。在SL1測線上,測量船朝西南方向行駛,在SL2上測量船的行駛方向則朝向東北,每條測線的總長度約為900 km。

圖12 線偏差測量的系統安裝圖Fig.12 System configuration for DOV measurement

圖13a 四組垂線偏差高頻量測量結果的對比(東-西分量)Fig.13a Comparison between four sets of high-frequency DOV measurement (east-west)

圖13b 四組垂線偏差高頻量測量結果的對比(北-南分量)Fig.13b Comparison between four sets of high-frequency DOV measurement (north-south)

如圖12所示,在測量船上安裝有兩套激光陀螺捷聯式慣性導航系統,即SINS1和SINS2,捷聯慣導的數據采樣率為1000 Hz。此外,船上還安裝有GNSS接收機(包括一個GNSS天線)用于獲取速度和位置信息。其定位精度為2 m,速度精度是0.03 m/s。GNSS數據的采樣率為1 Hz。SINS1與GNSS構成的組合導航系統命名為系統1(Sys1),SINS2與GNSS的組合則命名為系統2(Sys2)。

采用本文所提出的垂線偏差測量算法對兩套系統在SL1和SL2測線上采集得到四組實驗數據進行處理,通過四組實驗數據估計得到的垂線偏差擾動結果如圖13a、13b所示(為了更清楚地表現測量細節信息,圖中只顯示了經度從3.5°~7°范圍內垂線偏差擾動的測量結果)。可以看到,所估計得到垂線偏差擾動在0附近波動,其幅值約為2″。四次測量得到的垂線偏差擾動的重復性如表2所示,從表中可以看到,四次測量結果的重復性精度優于0.5″。

表2 四次測量的重復性精度Tab.2 Repeated accuracies of four sets of DOV data

4 結 論

本文提出一種垂線偏差動態測量的新方法。由于單軸旋轉INS/GPS組合導航系統的姿態誤差直接受垂線偏差的影響,本文方法通過對單軸旋轉INS/GPS組合的姿態誤差觀測以實現對垂線偏差的估計。利用真實的船搖數據生成仿真航跡,通過仿真驗證本文方法的可行性。仿真結果表明本文方法能夠有效地測量垂線偏差的高頻擾動量,以獲取更精細的地球重力場信息。該方法的測量精度依賴于所采用的陀螺性能,采用低角隨機游走性能的激光陀螺可以獲得較好的垂線偏差測量結果。通過實驗驗證該方法的有效性,結果表明該方法的垂線偏差測量重復性精度優于0.5″。

(References):

[1] Moritz H. Physical geodesy[M]. Bad V?slau, Austria: SpringerWienNewYork, 2005.

[2] Hirt C, Marti U, Bürki B, et al. Assessment of EGM2008 in Europe using accurate astrogeodetic vertical deflections and omission error estimates from SRTM/DTM2006.0 residual terrain model data[J]. Journal of Geophysical Research, 2010, 115(B10404).

[3] Hirt C, Bürki B, Somieski A, et al. Modern determination of vertical deflections using digital zenith camera[J]. Journal of Surveying Engineering-Asce, 2010, 136(1): 1-12.

[4] 黃謨濤, 翟國君, 管錚, 等. 海洋重力場測定及其應用[M]. 北京: 測繪出版社, 2005.

[5] Kwon J H, Jekeli C. A new approach for airborne vector gravimetry using GPS-INS[J]. Journal of Geodesy, 2001, 74: 690-700.

[6] Li X, Jekeli C. Ground-vehicle INS/GPS vector gravimetry[J]. Geophysics, 2008, 73(2): 11-21.

[7] Cai S, Zhang K, Wu M, et al. Improving airborne strapdown vector gravimetry using stabilized horizontal components[J]. Journal of Applied Geophysics, 2013, 98: 79-89.

[8] Senobari M S. New results in airborne vector gravimetry using strapdown INS/DGPS[J]. Journal of Geodesy, 2010, 84(5): 277-291.

[9] 李勝全, 歐陽永忠, 常國賓, 等. 慣性導航系統重力擾動矢量補償技術[J]. 中國慣性技術學報, 2012. 20(4): 410-413. Li Sheng-quan, Ouyang Yong-zhong, Chang Guo-bin, et al. Compensation technology of gravity disturbance vector in inertial navigation system[J]. Journal of Chinese Inertial Technology, 2012, 20(4): 410-413.

[10] Dai D K, Wang X S, Zhan D J, et al. An improved method for dynamic measurement of deflections of the vertical based on the maintenance of attitude reference[J]. Sensors, 2014, 14(9): 16322-16342.

[11] 戰德軍, 戴東凱, 張忠華, 等. 單軸旋轉INS/GPS組合導航中重力垂線偏差引起的姿態誤差分析[J]. 中國慣性技術學報, 2014, 22(3): 301-305. Zhan De-jun, Dai Dong-kai, Zhang Zhong-hua, et al. Analysis of gravity vertical deflection-induced attitude error in single-axis rotation INS/GPS integrated navigation system[J]. Journal of Chinese Inertial Technology, 2014, 22(3): 301-305.

[12] Wei M, Schwarz K P. A strapdown inertial algorithm using an Earth-fixed Cartesian frame[J]. Navigation, Journal of the Institute of Navigation, 1990, 37(2): 153-167.

[13] Pavlis N, Holmes K A, Kenyon S C, et al. The development and evaluation of the EGM2008[J]. Journal of Geophysical Research, 2012, 117(B04406): 1-38.

[14] Harriman D W, Harrison J C. A statistical analysis of gravity-induced errors in an airborne INS[C]//AIAA Guidance and Control, conference, 1984: 285-295

[15] Kriegsman B A, Mahar K B. Gravity-model errors in mobile inertial-navigation systems[J]. Journal of Guidance, 1986, 9(3): 312-318.

[16] Leonard J M, Nievinski F G, Born G H, et al. Gravity error compensation using second-order gauss Markov processes[J]. Journal of Spacecraft and Rockets, 2013, 50(1): 217-229.

[17] Gibbs B P. Advanced Kalman Filtering, least-squares and modeling[M]. New Jersey: John Wiley & Sons, INC., 2010.

[18] 龍興武, 于旭東, 張鵬飛, 等. 激光陀螺單軸旋轉慣性導航系統[J]. 中國慣性技術學報, 2011, 18(2): 149-153. Long Xing-wu, Yu Xu-dong, Zhang Peng-fei, et al. Single-rotating inertial navigation system with ring laser gyroscope[J]. Journal of Chinese Inertial Technology, 2011, 18(2): 149-153.

Measurement of dynamic vertical deflections by observing attitude errors of single-axis rotation INS/GPS system

DAI Dong-kai, WANG Xing-shu, ZHAN De-jun, WU Wei, HUANG Zong-sheng
(College of Opto-Electronic Science and Engineering, National University of Defense Technology, Changsha 410073, China)

Since the attitude errors of INS/GPS integrated system are directly related to deflections of the vertical (DOVs), it is expected that DOVs can also be estimated via the observation of INS/GPS attitude errors. In this paper, a laser gyroscope unit (LGU) is constructed by three laser gyroscopes of the INS, and used as an attitude reference for calculating INS/GPS attitude errors. Then, the observation equation and state-space equations for DOV measurement are accurately modeled. Finally, the DOVs and systematic errors are simultaneously estimated by a Kalman filter/smoother. The DOVs and systemic errors are decoupled since the state-space variables are accurately and properly modeled, and the long-wavelength components of DOVs are well compensated by a global gravity model. The feasibility of the proposed method are validated by simulation. Simulation results also show that this method can measure the high-frequency disturbances of DOVs efficiently. As the measurement accuracy of this method depends on the performance of gyroscope, a better accuracy of DOV measurement can be achieved when the gyroscopes with low angular random walk are adopted. The shipborne experiment results show that the repeated accuracy of the DOV measurement method is better than 0.5″.

deflection of the vertical; single-axis rotation INS; GPS; laser gyroscope unit

U666.1

A

1005-6734(2015)02-0172-07

10.13695/j.cnki.12-1222/o3.2015.02.007

2014-11-15;

2015-03-11

國家自然科學基金(61275002)

戴東凱(1986—),男,博士研究生,從事光電儀器與測控技術研究。E-mail:daidongkai@nudt.edu.cn

聯 系 人:王省書(1963—),女,教授,博士生導師,從事光電儀器與測控技術研究。

猜你喜歡
測量方法
把握四個“三” 測量變簡單
學習方法
滑動摩擦力的測量和計算
滑動摩擦力的測量與計算
測量的樂趣
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
測量
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 美女啪啪无遮挡| 精品欧美日韩国产日漫一区不卡| 欧美激情成人网| 国产欧美在线| 欧美日韩北条麻妃一区二区| 亚洲美女高潮久久久久久久| 国产精品午夜电影| 91原创视频在线| 亚洲综合极品香蕉久久网| 青青草国产在线视频| 精品人妻一区无码视频| 免费国产小视频在线观看| 干中文字幕| 免费精品一区二区h| 少妇精品网站| 伊人精品视频免费在线| 天天爽免费视频| 国产激爽爽爽大片在线观看| 狠狠色婷婷丁香综合久久韩国| 国模极品一区二区三区| 国产精品免费露脸视频| 毛片免费高清免费| 内射人妻无套中出无码| 国产精品一老牛影视频| 亚洲第一黄色网址| 国产三级视频网站| 又爽又黄又无遮挡网站| 日本一本在线视频| 女人av社区男人的天堂| 免费观看三级毛片| 91网站国产| 国产成人亚洲无吗淙合青草| 免费无码AV片在线观看国产| 精品少妇人妻无码久久| 国产网友愉拍精品视频| 激情综合婷婷丁香五月尤物 | а∨天堂一区中文字幕| 日韩欧美国产另类| 91在线日韩在线播放| 国产精品部在线观看| 99视频在线看| 亚洲va欧美va国产综合下载| 天堂在线www网亚洲| 亚洲综合在线最大成人| 伊人国产无码高清视频| 久久精品中文字幕少妇| 午夜成人在线视频| 久久午夜夜伦鲁鲁片无码免费| 久久激情影院| 色婷婷天天综合在线| 国产第一页免费浮力影院| 国产v欧美v日韩v综合精品| 熟妇丰满人妻av无码区| 国产综合欧美| 97精品国产高清久久久久蜜芽 | 亚洲第一区在线| 日本一本在线视频| 鲁鲁鲁爽爽爽在线视频观看| 成人噜噜噜视频在线观看| 刘亦菲一区二区在线观看| 国产精品一区二区久久精品无码| 真实国产乱子伦视频| 99国产在线视频| 日韩欧美色综合| 国产区在线看| 成人免费网站久久久| 91热爆在线| 久久一色本道亚洲| 免费观看亚洲人成网站| 欧美激情视频二区| 婷婷午夜影院| 欧美中文字幕在线播放| 精品人妻AV区| 日本欧美视频在线观看| 色视频久久| 狠狠色香婷婷久久亚洲精品| 免费一级毛片在线播放傲雪网| 不卡的在线视频免费观看| 国产精女同一区二区三区久| 看看一级毛片| 91精品日韩人妻无码久久| 欧美黑人欧美精品刺激|