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

月面著陸器與巡視器同波束差分時延相對定位算法

2015-10-28 02:17:43樊敏黃勇李海濤王宏郝萬宏陳少伍
航天器工程 2015年2期
關鍵詞:測量

樊敏黃勇李海濤王宏郝萬宏陳少伍

(1中國科學院上海天文臺,上海 200030)(2北京跟蹤與通信技術研究所,北京 100094)

月面著陸器與巡視器同波束差分時延相對定位算法

樊敏1,2黃勇1李海濤2王宏2郝萬宏2陳少伍2

(1中國科學院上海天文臺,上海 200030)(2北京跟蹤與通信技術研究所,北京 100094)

針對月球著陸巡視探測活動中的月面著陸器與巡視器的相對定位問題,建立了月面雙目標相對運動方程和狀態方程,給出了同波束差分時延測量量關于雙目標相對位置的測量方程,進而實現了基于統計估計方法解算雙目標相對位置的算法。結合嫦娥三號探測器跟蹤測量條件,利用該算法開展仿真分析。結果表明:在測量弧段達到5 min以上、同波束干涉測量(SBI)時延僅有1 ns隨機誤差的情況下,相對定位精度可達20 m;測量數據存在3 ns系統誤差時,相對定位精度為200 m,此時如果增加甚長基線干涉測量(VLBI)時延數據,可將相對定位精度提高到150 m。利用嫦娥三號實測數據處理結果驗證了此算法的正確性和仿真分析的有效性,可為合理制定月面雙目標相對定位策略提供參考。

月面探測;統計估計方法;相對運動方程;同波束干涉測量模型;甚長基線干涉測量

1 引言

2013年12月14日,我國嫦娥三號探測器首次實施了月球軟著陸并開展月面巡視探測。嫦娥三號探測器包含著陸器和巡視器,其中,著陸器在月面軟著陸后固定不動,巡視器與著陸器分離后進行巡視、勘察和采樣分析工作[1]。利用地基測量數據確定著陸器和巡視器在月面的精確位置,向探測器系統和科學應用系統提供高精度的地理位置信息,以保證準確開展巡視勘察活動,是嫦娥三號任務的一項關鍵技術,同時也可為地外天體巡視中經典的視覺導航技術提供參考和輔助驗證。我國喀什、佳木斯深空站和甚長基線干涉測量(VLBI)分系統的4個臺站對嫦娥三號探測器進行跟蹤測量。其中,對于著陸器,深空站可以進行測距、測速,VLBI分系統的臺站可以進行雙差分單向測距(ΔDOR),由此測定著陸器的位置,通過連續跟蹤測量,利用統計方法還可逐步改進著陸器月面位置信息,實現著陸器高精度定位。但是,考慮質量、功耗和科學目標等因素,巡視器上配置的X頻段測控數傳設備不具有利用地面站進行測距和ΔDOR測量的能力,因此無法直接測定巡視器的月面位置。不過,嫦娥三號巡視器僅在著陸器附近區域(相距不超過10 km)進行勘察,由于地月距離遙遠,雙目標相距較近的特點使得它們相對地面站的角距離非常近,可以在地面天線的同一波束內對雙目標進行跟蹤測量。利用兩副地面站天線對雙目標同時進行測量,可以生成高精度的同波束干涉測量(SBI)數據,在美國阿波羅-16、17飛船月球探測任務中,就成功地運用SBI技術確定了月球車相對登月艙的運動軌跡[2]。因此,利用SBI數據是實現著陸器和巡視器相對定位的一種手段,須要進一步研究和實現具體的月面雙目標相對定位。

傳統的單點定位算法利用單歷元時刻的多個測量數據進行定位解算,定位精度主要取決于測量系統的幾何構型和測量數據的系統誤差。由于月面目標距離地球遙遠,利用地基測量系統對其跟蹤測量的幾何構型較差,而且地基測量數據通常有形式復雜的系統誤差,因此,根據現有測量條件和數據精度,利用單點定位算法進行月面目標定位的精度僅能達到數百米量級。為了提高定位精度,可以將長弧段的測量數據歸算到同一定位時刻,再進行統計平差,但歸算過程無法考慮目標運動規律,從而又引入了額外誤差[3]。針對這些問題,本文首先建立了雙目標的相對運動方程和SBI差分時延測量方程;然后利用統計估計方法,建立解算雙目標相對位置的算法并加以實現。在此基礎上,根據嫦娥三號著陸器與巡視器的實際跟蹤測量條件進行仿真,并對嫦娥三號探測器獲取的實測數據應用此算法進行相對定位,驗證了算法的正確性和有效性。

2 算法原理

利用統計估計方法對月面雙目標進行相對定位,主要借鑒航天器精密定軌的理論[4-5],包括以下3個方面內容。

2.1 相對運動方程和狀態方程

考慮到著陸器在月面著陸后,巡視器行進到指定探測點開展就位探測期間,兩器在月心固連坐標系中均靜止不動,建立的相對運動方程為

式中:Δr為著陸器位置矢量r1和巡視器位置矢量r2之差,即相對位置矢量,可用直角坐標表示,也可用月球地理坐標表示;Δr0為相對位置矢量初始值;為著陸器與巡視器相對速度矢量。

考慮相對位置、速度信息以外的其他待估參數Pg,例如,影響運動狀態的物理參數、地面站坐標的幾何參數和測量數據的系統誤差等,定義狀態矢量X的坐標為由此可得狀態方程為

式中:狀態矢量初值X0的坐標為[Δr00 Pg0],其中,Δr0可取標稱值或根據巡視器遙測數據計算,Pg0可取理論設計值或經驗值。

2.2 SBI測量模型和測量方程

SBI測量原理如圖1所示,當著陸器與巡視器在角度上非常接近時,可在一副地面天線的同一主波束內被觀測,使用兩副天線同時對其進行觀測,即可生成差分干涉測量量[6]。由于SBI的測量量能夠精確確定兩個目標在天平面內的相對位置信息,因此可用于相對定位[7]。

假設著陸器(記為SC1)發射信號的時刻為t1,地面站A、B接收到著陸器信號的時刻分別為t和t2。將t時刻A站接收到著陸器t1時刻發射信號的光行時記為著陸器到A站和B站的幾何時延記為則著陸器t1時刻發射的信號到達B站的時刻那么,可以將信號從著陸器到B站的光行時記為同時,假設在太陽系質心慣性參考坐標系中,t時刻著陸器與巡視器(記為SC2)的位置矢量分別為RSC1和RSC2,A站和B站的位置矢量分別為RA和RB,著陸器t1時刻發射信號到A站和B站的光行時(c表示光速)可以表示為

著陸器t1時刻發射信號到達A站和B站的幾何時延可以表示為

同理可得,巡視器t3時刻發射的信號到達兩站的幾何時延可以表示為

因此,A站和B站對著陸器與巡視器的SBI差分時延測量模型可表示為

式中:ΔR為著陸器與巡視器在太陽系質心慣性參考坐標系中的相對位置矢量。

圖1 著陸器與巡視器同波束干涉測量原理Fig.1 SBI measurement principle of lander and rover

可見,SBI差分時延測量量包含雙目標相對位置信息,通過固定著陸器位置,利用該測量量可以解算相對位置。

由于測量模型是在太陽系質心慣性參考坐標系中建立的,相應的坐標時為太陽系質心動力學時(TDB),在進行光行時解算的時候,要計算著陸器和巡視器在太陽系質心慣性參考坐標系中的相對位置,而且要考慮各大天體引力時延等相對論影響。由于著陸器與巡視器相對狀態方程是在月心固連坐標系中建立的,而且最終的待估參數為雙目標在月心固連坐標系中的相對位置,因此須要考慮月心固連坐標系到太陽系質心慣性參考坐標系的轉換。該轉換過程涉及到的時間系統包括:地面站采用的協調世界時(UTC)、原子時(TAI)、地球時(TT)和TDB[8]。涉及到的坐標系如表1所示。

表1 相對運動方程和測量方程涉及的坐標系Table 1 Coordinate systems used in relative movement equation and measurement equation

此外,月心固連坐標系包括主軸(Principal Axes)和平軸(Mean Rotation Axes)坐標系兩種[9]。其中,月球歷表給出的月球天平動參數和月球重力場采用的是主軸坐標系,而國際天文聯合會(IAU)定義的平軸坐標系主要用于描述月面地形和特征。目前,主軸坐標系和平軸坐標系的差異在月面小于1 km,可以根據不同的JPL DE/LE系列歷表給出的轉換參數進行轉換。對于本文采用的DE 421歷表[10],轉換關系為

式中:P和M分別為主軸坐標系和平軸坐標系中點位置的坐標矢量;Ra(a代表x,y,z)表示繞坐標軸a的旋轉矩陣。

基于SBI差分時延的測量模型和上述時間坐標系轉換關系,可以建立對著陸器與巡視器進行連續SBI測量的測量方程。將測量值與狀態矢量X之間的函數關系記為G(X,t),考慮測量噪聲,則測量方程為

式中:Yi為ti時刻的測量量;Xi為ti時刻的狀態矢量;εi為ti時刻的測量噪聲。

測量模型中的函數關系是非線性的,須要對其線性化。將測量方程在參考狀態X*(ti)處展開,令

式中:Φ(ti,t0)為狀態轉移矩陣,由于著陸器與巡視器在月心固連坐標系中相對靜止,該矩陣為單位矩陣;矩陣Hi中包含雙目標相對位置矢量從月心固連坐標系到太陽系質心慣性參考坐標系的轉換矩陣及偏導數。

記χ0=X0—X*0,則線性化后的測量方程為

令y=[y1…yk]T,H=[H1…Hk]T,ε=[ε1…εk]T,其中k為測量量的個數,于是可將所有的測量方程總寫為

2.3 統計估計方法

建立測量方程后,需要解決的問題是如何確定上述線性系統的最優估計。通常求解這類問題的直接方法是加權最小二乘批處理算法。根據加權最小二乘估計理論,如果已知待估計參數χ0的先驗估計和先驗估計的加權矩陣,則批處理算法解算的χ0的“最佳”估值為

式中:W為權矩陣。

至此,完成了基于月面目標運動規律的統計估計相對定位算法,實現了對著陸器與巡視器相對位置的解算。對于嫦娥三號著陸器,考慮到其著陸后可以利用較長弧段的深空站測距、測速,以及三向測量和VLBI分系統ΔDOR干涉測量數據,解算出高精度的月心固連坐標系位置信息,因此可增加對著陸器位置的先驗約束來解算著陸器與巡視器的相對位置。同時,還可以考慮利用目前高精度的月面數字高程模型(DEM)提供高程信息,如美國的ULCN2005模型[11]和我國自主研制的DEM模型[12-13]等,增加相對位置中的高程約束來解算雙目標的相對位置。

3 仿真分析與驗證

利用FORTRAN程序實現了上述雙目標定位算法。根據嫦娥三號著陸器與巡視器的跟蹤測量條件,在完成兩器分離后,著陸器定向天線向地面發射數傳信號,而巡視器在開展就位探測期間(每次約20 min),其測控數傳設備也將發射數傳信號,地面站可以利用這些單向數傳信號對雙目標進行SBI干涉測量,以及對單目標進行VLBI干涉測量和單向測速。考慮數傳信號的設計形式和實現方式,目前對雙目標的SBI干涉測量時延精度為1 ns,對單目標的VLBI干涉測量時延精度為1 ns,對單向測速精度為5 cm/s。針對上述實際情況,仿真生成了在一個跟蹤弧段內上海、北京、昆明、烏魯木齊4臺站的SBI和VLBI干涉測量時延數據,以及喀什、佳木斯深空站的測速數據,具體仿真條件如表2所示。在此基礎上,利用本文算法對雙目標相對定位精度進行仿真分析。

表2 仿真分析條件Table 2 Simulation analysis terms

3.1 利用SBI數據相對定位

根據SBI差分時延的測量模型式(6)可知,SBI差分時延數據除包含雙目標相對位置信息外,還包含著陸器本身的位置信息,因此,僅利用SBI差分時延數據進行相對定位時,要增加對著陸器位置的先驗約束來解算著陸器與巡視器的相對位置。表3和表4給出了不同測量弧段和不同解算策略下的仿真計算結果。其中,以著陸器著陸點為坐標原點的北東地坐標系中,X方向為當地正北方向,Y方向為當地正東方向,Z方向與X方向和Y方向構成右手系。

表3 利用SBI數據進行相對定位的仿真分析結果Table 3 Simulation analysis results of relative position determination using SBI data m

表4 解算高程和固定高程的相對定位仿真誤差橢球比較Table 4 Comparison for error ellipsoid of simulation analysis results of relative position determination between estimated and fixed altitudes

可以看出:

(1)利用高精度月面DEM模型增加高程約束,對提高相對定位精度有重要作用。不考慮SBI數據系統誤差的情況下,比較同時解算相對緯度、經度和高程的情況與只解算相對緯度、經度的情況,著陸器與巡視器的相對定位誤差由百米降至10 m,弧段增加到20 min時,還可達到米級。此外,通過對解算時迭代次數的統計發現,解算高程通常要迭代10次以上,而增加高程約束可以有效減少迭代次數至3~4次,從而提高計算效率。

(2)測量數據系統誤差對相對定位精度影響較大。不考慮測量數據系統誤差時,測量弧段增加可有效提高相對定位精度,20 min弧段的相對定位精度優于10 m。但是,當測量數據存在系統誤差時,增加數據弧段對相對定位精度的進一步提高作用不大。測量弧段在5 min以上的情況下,當SBI差分時延數據系統誤差為2 ns時,相對定位誤差小于130 m;當系統誤差為3 ns時,相對定位誤差小于210 m。

綜合考慮上述兩點,在著陸器與巡視器雙目標相對定位時,VLBI分系統應保證獲取5 min以上的SBI數據,并盡可能消除測量系統誤差,在解算相對位置時,考慮DEM模型增加高程約束,以提高著陸器與巡視器的相對定位精度。

3.2 利用VLBI和SBI數據聯合相對定位

考慮到僅利用SBI數據進行相對定位時,要增加對著陸器位置的先驗約束以解算著陸器與巡視器的相對位置,而實際上,地面站在對雙目標進行SBI測量的同時可以獲取單目標的VLBI時延數據。因此,可以綜合利用這兩種數據進行相對定位。表5給出了5 min以上測量弧段、增加高程約束策略下的仿真計算結果。

表5 利用SBI和VLBI數據聯合相對定位的仿真分析結果Table 5 Simulation analysis results of relative position determination using SBI and VLBI data m

比較表5與表3可以看出:①在測量數據沒有系統誤差時,增加VLBI時延測量數據對提高相對定位精度的作用不明顯;②考慮SBI數據的系統誤差,增加VLBI數據可將相對定位誤差由210 m降低到150 m。考慮到SBI數據可能存在系統誤差,因此著陸器與巡視器相對定位時,應綜合利用VLBI數據來提高相對定位精度。

根據地面站對著陸器與巡視器的測量條件,VLBI分系統4個臺站在對雙目標進行SBI測量,以及單目標進行VLBI測量的同時,深空站可以獲取對單目標的單向測速數據。為了分析單向測速數據對相對定位精度的影響,仿真分析了綜合利用單目標測速和VLBI數據,以及雙目標SBI數據進行相對定位的誤差,如表6所示。

表6 利用測速和SBI及VLBI數據聯合相對定位的仿真分析結果Table 6 Simulation analysis results of relative position determination using range rate,SBI and VLBI data m

通過比較表5和表6可見,對于5 min以上測量弧段,采用增加高程約束的相對定位策略,無論是否考慮測量數據的系統誤差,加入測速數據對相對定位精度的影響在厘米級,遠小于相對定位本身能夠實現的精度。因此,在著陸器與巡視器相對定位時,可以不采用深空站的測速數據。

2013年12月14日21:12,嫦娥三號探測器成功實施月面軟著陸,之后著陸器與巡視器分離,測控系統獲取了實測數據。通過對2013年12月15—24日獲取的測量數據進行分析,可得SBI數據隨機誤差約為1 ns,系統誤差約為2 ns。根據本文算法,采用第2.3節增加高程約束的策略,利用實測數據進行相對定位,以精度達到厘米級的視覺導航系統二維相對定位結果為評估基準進行分析,結果如圖2所示。從初步分析的10個弧段的相對定位結果來看,最大偏差約為156 m,最小偏差約為53 m,與表3中仿真分析的結果基本吻合,驗證了本文算法的正確性和仿真分析結果的有效性。

圖2 嫦娥三號著陸器與巡視器實測數據相對定位精度Fig.2 Precision of relative position determination using measurement data of Chang'e-3 lander and rover

4 結束語

本文提出的月面著陸器與巡視器相對定位算法,源于經典的航天器動力學統計定軌方法,根據月球的平動和轉動模型建立雙目標在空間的運動模型,不必考慮復雜的空間飛行動力學建模問題,運動模型的精度僅取決于月球物理天平動參數和月面地理參數的精度,而目前國際通用的月球歷表和地理參數精度已達到厘米級,因此完全滿足相對定位需求。同時,該算法不同于常規的幾何定位歸算算法,可以利用雙目標相對運動方程精確計算任意時刻的測量值,從而提高解算精度。利用此算法對嫦娥三號探測器實測數據進行處理分析,相對定位結果與仿真結果一致,表明算法正確、有效。利用本文算法對著陸器與巡視器實際測量條件開展仿真分析,可以為合理、有效地制定月面雙目標相對定位策略提供參考,后續也可以進一步用于著陸地外天體的兩探測器之間的相對定位。

[1]劉慶會,陳明,熊蔚名,等.基于超高精度多頻點同波束VLBI技術的月球車精密相對定位[J].中國科學:物理學力學天文學,2010,40(2):253-260 Liu Qinghui,Chen Ming,Xiong Weiming,et al.High accurate relative positioning of lunar rover using super high accurate multi point frequencies same beam VLBI[J].Science of China,Scientia Sinica Physical,Mechanical& Astronomical,2010,40(2):253-260(in Chinese)

[2]Salzberg I M.Tracking the Apollo lunar rover with interferometry techniques[J].Proceedings of the IEEE,1973,61(9):1233-1236

[3]黃勇,胡小工,李培佳,等.運動學統計定位方法應用于月球著陸器的精密定位[J].科學通報,2012,57(28/29):2686-2692 Huang Yong,Hu Xiaogong,Li Peijia,et al.Precise positioning of the Chang'e-3 lunar lander using a kinematic statistical method[J].Chinese Science Bull,2012,57(28/29):2686-2692(in Chinese)

[4]李濟生.航天器軌道確定[M].北京:國防工業出版社,2003:255-258 Li Jisheng.Spacecraft orbit determination[M].Beijing: National Defense Industry Press,2003:255-258(in Chinese)

The paper is written in a good language, the logic is clear and the subject and results are discussed graphically and meaningfully.

[5]劉林,胡松杰,王歆.航天動力學引論[M].南京:南京大學出版社,2005:55-58 Liu Lin,Hu Songjie,Wang Xin.An introduction of astrodynamics[M].Nanjing:Nanjing University Press,2005:56-58(in Chinese)

[6]Thornton C L,Border J S.深空導航無線電跟蹤測量技術[M].李海濤,譯.北京:清華大學出版社,2005:3-5 Thornton C L,Border J S.Radiometric tracking techniques for deep space navigation[M].Li Haitao,translated.Beijing:Tsinghua University Press,2005:3-5(in Chinese)

[7]董光亮,郝萬宏,李海濤,等.同波束干涉測量對月面目標相對定位[J].清華大學學報(自然科學版),2010,50(7):1118-1124 Dong Guangliang,Hao Wanhong,Li Haitao,et al.Relative position determination on the lunar surface using same-beam interferometry[J].Journal of Tsinghua University(Science&Technology),2010,50(7): 1118-1124(in Chinese)

[8]Petit G,Luzum B.IERS conventions(2010)[R/OL].[2013-10-01].http://www.iers.org/TN36/

[9]Archinal B A,A'Hearn M F,Bowell E,et al.Report of the IAU Working Group on cartographic coordinates and rotational elements:2009[J].Celestial Mechanics and Dynamical Astronomy,2011,109(2):101-135

[10]Folkner W M,Williams J G,Boggs D H.The planetary and lunar ephemeris DE421[R/OL].[2013-12-01].http://tmo.jpl.nasa.gov/progress_report/42-178/178C.pdf

[11]Archinal B A,Rosiek M R,Kirk R L,et al.Completion of the unified lunar control network 2005 and topographic model[C]//Proceedings of the 37th Annual Lunar and Planetary Science Conference.Houston:Lunar and Planetary Institute,2006

[12]平勁松,黃倩,鄢建國,等.基于嫦娥一號衛星激光測高觀測的月球地形模型CLTM-s01[J].中國科學:物理學力學天文學,2008,38(11):1601-1612 Ping Jinsong,Huang Qian,Yan Jianguo,et al.Lunar topographic model CLTM-s01 from Chang'e-1 laser altimeter[J].Science of China,Scientia Sinica Physical,Mechanical&Astronomical,2008,38(11):1601-1612(in Chinese)

[13]李春來,任鑫,劉建軍,等.嫦娥一號激光測距數據及全月球DEM模型[J].中國科學:地球科學,2010,40(3):281-293 Li Chunlai,Ren Xin,Liu Jianjun,et al.Laser altimetry data of Chang'e-1 and the global lunar DEM model[J].Science of China,Earth Science,2010,40(3): 281-293(in Chinese)

(編輯:夏光)

Algorithm of Relative Positioning for Lander and Rover on Lunar Surface Using Differential Time-delay of SBI

FAN Min1,2HUANG Yong1LI Haitao2WANG Hong2HAO Wanhong2CHEN Shaowu2
(1 Shanghai Astronomical Observatory,Chinese Academy of Sciences,Shanghai 200030,China)
(2 Beijing Institute of Tracking and Telecommunications Technology,Beijing 100094,China)

For relative positioning of the lander and the rover in lunar soft landing and surface reconnaissance exploration,an algorithm is presented that can be used in relative positioning for two objects based on the statistical estimation method by establishing the relative kinematic equation,state equation and the measurement model of SBI(same beam interferometry).A simulation of this algorithm is carried out based on emulational measurement data of Chang'e-3 TT&C system.The result shows that the relative positioning error can be 20m when stochastic noise of SBI delay data is 1ns and tracking arc is longer than 5min.Moreover,if there is a bias of 3ns in SBI delay data,the relative positioning error increases to 200m,which can be declined to 150m by adding VLBI delay data.The validity of the algorithm and the effectiveness of the simulation are proved by the analysis of Chang'e-3 measurement data.The algorithm can support effectively determining a strategy of relative positioning of two objects on lunar surface.

exploration on lunar surface;statistical estimation method;equation of relative movement;model of same beam interferometry;very long baseline interferometry

V556

A DOI:10.3969/j.issn.1673-8748.2015.02.003

2014-04-23;

2014-08-13

國家自然科學基金(11473056,11403076)

樊敏,女,工程師,研究方向為航天測控、軌道動力學。Email:fanmin@bittt.cn。

猜你喜歡
測量
測量重量,測量長度……
把握四個“三” 測量變簡單
滑動摩擦力的測量和計算
滑動摩擦力的測量與計算
測量的樂趣
二十四節氣簡易測量
日出日落的觀察與測量
滑動摩擦力的測量與計算
測量
測量水的多少……
主站蜘蛛池模板: 國產尤物AV尤物在線觀看| 国产精品播放| 日韩精品亚洲精品第一页| 99伊人精品| 国产精品自在在线午夜区app| 国产h视频免费观看| 丁香五月亚洲综合在线 | 日韩精品欧美国产在线| 91小视频在线观看| 国产成人精品亚洲77美色| 曰韩免费无码AV一区二区| 极品私人尤物在线精品首页 | 9久久伊人精品综合| 久久99国产精品成人欧美| 欧美全免费aaaaaa特黄在线| 亚洲综合精品香蕉久久网| 久草视频一区| 成人免费午夜视频| 国产黄色片在线看| 国产99视频精品免费观看9e| 亚洲国产欧美国产综合久久 | 亚洲成av人无码综合在线观看| 亚洲国产综合第一精品小说| 日本黄网在线观看| 黄色网页在线播放| 亚洲综合二区| 亚洲欧美日韩中文字幕一区二区三区| 国产精品天干天干在线观看| 国产区福利小视频在线观看尤物| 熟妇人妻无乱码中文字幕真矢织江| 亚洲男女在线| 日韩天堂网| 日本精品一在线观看视频| 亚洲一区二区三区中文字幕5566| 国产人人干| 亚洲精品桃花岛av在线| 亚洲大学生视频在线播放| 亚洲妓女综合网995久久| 九色在线观看视频| 亚洲三级a| 欧美一区二区自偷自拍视频| 免费Aⅴ片在线观看蜜芽Tⅴ| 午夜人性色福利无码视频在线观看| 国产交换配偶在线视频| 日韩无码黄色网站| 成人年鲁鲁在线观看视频| 日韩成人午夜| 天堂成人在线视频| 久爱午夜精品免费视频| 国产高清在线精品一区二区三区| 97精品伊人久久大香线蕉| 97se亚洲综合在线| 久久夜色精品| 国产欧美视频在线| 日本少妇又色又爽又高潮| 无码日韩人妻精品久久蜜桃| 亚洲色图欧美在线| 波多野结衣AV无码久久一区| 欧美国产日产一区二区| 一级看片免费视频| 亚洲成人一区二区三区| 中文国产成人精品久久一| 91外围女在线观看| 亚洲第一福利视频导航| 中文字幕亚洲电影| 一本久道热中字伊人| 国产乱肥老妇精品视频| 国产美女自慰在线观看| 无码高潮喷水在线观看| 日韩欧美91| 国产日韩欧美一区二区三区在线 | 亚洲色图综合在线| 久久这里只有精品国产99| 国产精品久久久精品三级| 亚洲一区第一页| a级毛片免费看| 国产无吗一区二区三区在线欢| 欧美精品啪啪| 无码中文字幕乱码免费2| h视频在线观看网站| 亚洲精品国产首次亮相| 国产黑丝视频在线观看|