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

含引力輔助的小推力多任務探測軌道設計

2017-03-06 01:12:35林浩申安喜彬
固體火箭技術 2017年1期
關鍵詞:優化設計

林浩申,何 兵,趙 欣,劉 剛,安喜彬

(第二炮兵工程大學 空間工程系,西安 710025)

含引力輔助的小推力多任務探測軌道設計

林浩申,何 兵,趙 欣,劉 剛,安喜彬

(第二炮兵工程大學 空間工程系,西安 710025)

針對深空探測時間長、成本高的特點,文中以多任務深空探測為背景,建立了包含多模型、多約束、多變量、多階段的探測軌道設計模型。通過融合推力工作模式、分段設計目標函數等策略,克服了模型的內點約束限制,降低了優化模型的復雜度。為了提高優化速度、提升優化精度,結合梯度搜索和粒子群算法的特點,研究提出了一種梯度混合粒子群(GHPSO)算法。將該算法應用到多任務星探測軌道設計模型上,得到了發動機的工作時序,并橫向對比了該文算法與傳統粒子群算法(PSO)和遺傳算法(GA)的優化性能。仿真結果表明,文中提出的算法搜索速度快,第1設計階段GHPSO相對GA提高61.96%,相對PSO提高47.85%,第2設計階段GHPSO相對GA提高61.87%,相對PSO提高43.66%;精度最高,第1設計階段GHPSO相對GA平均精度提高了2.86%,最高精度提高了1.24%,相對PSO平均精度提高了4.19%,最高精度提高了3.97%,第2設計階段GHPSO相對GA平均精度提高了3.33%,最高精度提高了1.63%,相對PSO平均精度提高了4.72%,最高精度提高了3.02%,適用于軌道優化設計類的非線性、多約束全局優化問題。

梯度混合粒子群算法;多任務星際探測;引力輔助

0 引言

自1959年前蘇聯發射第一顆月球探測器以來[1],人類邁開了深空探測的步伐,并隨著空間科學技術的發展不斷向更遠、更深的宇宙進發。1977年,美國宇航局開始醞釀一個旨在研究外太陽系的“大旅行”計劃, “旅行者一號”借助176年一遇的行星連線的契機,途經木星、土星歷經33年飛出太陽系[2]。針對深空探測時間長、成本高的特點,設計更快更高效、更遠更經濟的深空探測軌道是人類探索深空的必然趨勢。

目前的研究主要從三方面提高深空探測飛行任務的效率[3]:一是采用小推力技術提高推進系統的效率[4];二是利用行星引力場的引力輔助技術[5]和利用天體平衡點的不變流形低能轉移技術[6];三是多目標多任務探測[7],是追求低成本、高效率的有效途徑。

小推力轉移軌道優化設計可歸結為非線性常微分方程組的兩點邊值問題。求解該類問題的主要方法包括間接法、直接法、混合法[8]等。間接法實質上通過增加函數乘子和數值乘子作為待求量,把狀態微分方程約束和內點約束引入到廣義指標中,再進一步求解,該方法的優點是待求量少,解的精度高,但其收斂域窄、對初值敏感的特點增加了該算法的操作難度[8];直接法操作簡單,容易得到可行解,但難以保證解的最優性[9];混合法結合了前述兩種方法的優點,但還需要推導一階必要條件和事先假定控制策略[8]。

19世紀,Leverrier和Tisserand在天體軌道攝動方面的研究工作是借力飛行技術的源頭[10],利用該技術探測器可在不消耗燃料的前提下獲得速度增量。借力飛行軌道設計的本質是全局優化問題,通常采用隨機優化方法求解,如遺傳算法(Genetic Algorithm)、進化算法(Evolutionary Algorithm)和模擬退火算法等(Simulated Annealing)。文獻[11]首次將遺傳算法應用在行星探測軌道優化設計問題中。研究結果表明,遺傳算法求解精度較高,但收斂速度較慢。粒子群算法作為智能優化算法的一種,因其參數設置少、收斂速度快等優點,在軌跡優化領域的應用日益受到關注[12],文獻[13]研究了基于改進粒子群算法的航天器燃料最省的多脈沖交匯問題,該算法收斂速度較快,但求解精度較低。

本文將梯度搜索的思想融入粒子群算法之中,利用梯度搜索的高效性和粒子群搜索的隨機性解決粒子群算法易“早熟”及收斂精度低等問題,提出了一種改進的混合粒子群算法,并將該算法應用在含引力輔助的小推力軌道設計上。

1 基礎力學模型

1.1 探測器的動力學模型

探測器的推進方式可分為脈沖推進方式和小推力推進。脈沖推進方式可視為探測器在推力作用時刻獲得瞬時脈沖,使得探測器進行Lambert轉移,此種推進方式需要優化的參數少,通過解算Lambert軌道可控制探測器在指定時刻轉移至指定位置;小推力推進方式需要優化的參數多,通過其微分模型可積分得到探測器的實時狀態。文獻[14]給出了脈沖推力模型、小推力模型和Lambert轉移模型的詳細形式,固定時間雙脈沖轉移示意圖如圖1所示。

1.2 引力輔助模型

行星引力輔助是指在探測器飛行過程中,利用行星引力場改變自身速度的過程[5],探測器從無窮遠處以相對速度v∞在tGA時刻進入借力天體P的影響球,引力輔助過程示意圖[15]如圖2所示。

(1)

探測器相對行星的速度大小:

(2)

(3)

式中Rplanet為行星的飛越半徑,取值應不小于借力行星半徑;uplanet為借力行星的引力常數。

進行行星引力輔助時,探測器的位置還應該滿足如下約束條件:

(4)

式中rsc、rplanet分別為航天探測器和借力行星在日心固聯坐標系(HEIRF)中的位置矢量;Rp為行星影響球的半徑。

利用行星引力輔助后,探測器在不消耗任何燃料的情況下獲得的脈沖速度為

(5)

2 多任務深空探測軌道優化模型

2.1 約束條件

深空軌道設計優化的約束條件主要包括發射窗口約束、初始狀態約束、終點狀態約束、微分方程組約束和行星引力輔助約束等。

(1)發射窗口約束

行星在HEIRF內的空間位置是不斷變化的,因此發射窗口的選擇對于軌道優化尤為重要,航天任務規劃通常是在一個相對較長的時間段內,從地球出發時刻t0應滿足如下約束:

(6)

式中 MJDmin、MJDmax分別表示發射任務時間選取的最小值和最大值(改進儒略日)。

(2)探測器初始狀態約束

在從地球出發的時刻( t0),運載火箭能為探測器提供一個最大逃逸速度vΔmax和最大逃逸高程hΔmax,因此t0時刻的探測器位置、速度矢量應該滿足如下約束:

(7)

式中m0為探測器初始質量;mmax為探測器的最大運載質量。

(3)探測器終點狀態約束

根據不同的空間探測任務,終點約束不盡相同,以設計逃離太陽系的空間軌道為例,探測器在終點時刻tend的狀態應該滿足:

(8)

即探測器與太陽的相對距離不小于40個天文單位,相對太陽的軌道偏心率不小于1,以保證探測器能完全脫離太陽的引力束縛,在無動力的情形下飛出太陽系;質量變化量不大于攜帶的燃料質量mfuel。

(4)探測器的運動微分方程約束

探測器在化學推進方式、小推力推進方式及無動力這3種情形下運動方程由文獻[14]已經詳細給出。

(5)行星引力輔助約束

在軌道設計過程中若選擇利用行星引力輔助,引力輔助應該滿足式(1)~式(5)。

2.2 優化變量的設計

優化變量的設計和優化模型的構造緊密相關,以歐空局Izzo的MGA-1DSM模型[16]為例,該模型描述了探測器在2個行星軌道間的轉移過程,示意圖如圖4所示。圖中,ΔvDSM為探測器進行Lambert轉移時需要的速度增量;探測器借助行星引力輔助時獲得的速度增量為ΔvGA;T0、T1分別為此次轉移的開始時刻和結束時刻;η1為此次轉移過程中Kepler段時間占該次轉移總時間的百分比。

整個轉移過程分為Kepler和Lambert兩段,Kepler段探測器采用小推力沿著速度方向以最大推力持續工作;Lambert段探測器采用脈沖轉移。

深空探測軌道實質上是若干個上述軌道的拼接,因此可設計優化變量為

(9)

式中t0、m0分別表示發射時刻和初始質量;rsc(t0)、αr、βr描述了探測器的初始位置矢量rsc(t0)=r(t0)[cos(βr)cos(αr),cos(βr)sin(αr),sin(βr)];而vsc(t0)、αv、βv則描述了探測器的初始速度矢量vsc(t0)=v(t0)[cos(βv)cos(αv),cos(βv)sin(αv),sin(βv)];探測器第i次行星引力輔助時獲得的速度增量ΔvGAi由δi、γi確定;ΔTi=Ti-Ti-1表示第i次轉移耗費的總時間,ηi表示第i次轉移過程中Kepler段占該次轉移總時間的百分比。

2.3 優化變量的設計

本文研究的是小推力多任務行星探測高速軌道優化設計問題,整體設計軌道時需要優化的變量多,直接得到滿足所有內點約束的全局最優解難度較大,因此合理的設計目標函數至關重要,本文對目標函數的設計主要包括以下幾步:

(1)混合推力模式軌道設計

直接對小推力軌道進行優化設計的難點在于無法瞄準目標大行星,難以滿足引力輔助的內點約束,因此假設探測器能采用化學推進和連續小推力的混合推進方式,以上述MGA-1DSM模型為基礎,探測器通過Lambert轉移能準確抵達目標行星,優化設計過程中無須考慮引力輔助約束。將目標函數設計如下:

(10)

式中 ΔTi=Ti-Ti-1為第i次轉移耗費的總時間;Δt1為探測器最后一次行星借力后發動機工作的時間;Δt2為探測器燃料耗盡后滑行的時間;ΔvDSMi為第i次轉移過程中Lambert轉移需要的速度增量。

式中加入eΔvDSMi項,能引導ΔvDSMi不斷減小,使得各次Lambert脈沖轉移趨于0。通過該目標函數能對發射窗口、探測器初始狀態以及引力輔助策略進行優化,得到一組滿足各項約束的優化結果。

(2)小推力軌道設計

去掉模型中的Lambert轉移軌道的同時引入了新的問題:無法瞄準下一顆借力的大行星。為此重新設計指標函數,新的指標函數設計為

(11)

式中di為Ti時刻探測器與目標行星之間的距離。

通過此指標函數引導探測器在無Lambert轉移軌道的情形下瞄準行星,進而在滿足引力輔助約束的前提下使用引力輔助。

3 梯度混合粒子群算法

3.1 算法思想

傳統粒子群算法通過粒子跟蹤個體最優解和全局最優解更新粒子的位置矢量和速度矢量,這種搜索思想使得粒子群算法的全局搜索能力較弱,同時也決定了該算法容易陷入局部最優出現“早熟”現象。本文結合梯度搜索的思想,對傳統粒子群的搜索方式進行改變,利用相鄰代群的全局最優值得到搜索梯度,通過該梯度指引粒子群更新位置和速度矢量,提高粒子群算法的全局搜索能力。

3.2 算法流程

(1)種群初始化

隨機初始化第1代群體與第2代群體各粒子的位置向量和速度向量,計算各粒子的適應度函數值,同時記錄前兩代群體中的最優個體及全局最優個體。

(2)適應度評價

(3)梯度搜索更新

其中,α為一正常數,α越大,算法全局搜索能力越強,α越小,算法局部搜索能力越強。α的取值策略為

(12)

(4)生成新一代種群

(13)

(5)終止條件判定

判斷優化結果是否達到給定精度或迭代次數是否達到給定值。若滿足給定條件,則終止迭代,并輸出最優解,否則轉(2)。

該算法的優化流程如圖5所示。

4 仿真計算與結果分析

算例問題描述:探測器將于2025年1月1日(儒略日表示為54 000 MJD)至2055年12月31日之間的任意時刻從地球出發,出發時滿足式(12)的初始狀態約束,其中最大逃逸速度vΔmax=3 km/s,最大逃逸高程hΔmax=1 000 km且方向任意。探測器需要完成的任務包括:(1)利用火星的引力輔助效應,同時對火星進行觀察;(2)利用木星的引力輔助效應,同時投放木星探測器;(3)逃離太陽系,進行深空探測。

假設探測器在飛行過程中僅受太陽引力(引力輔助時刻除外)并可在任意時刻在滿足引力輔助約束的前提下利用大行星的引力輔助效應,設計一條能完成上述多任務的探測軌道。

假設設備最大載重為1 800 kg,深空探測設備自重500 kg,攜帶的木星探測設備自重200 kg,攜帶燃料的質量可根據任務自由調整,燃料罐的質量為所攜帶燃料的5 %。

探測器的推進系統采用有限推力的電推進,其比沖為3 000 s,推力最大值為0.5 N,推力的大小和方向均可優化。

(1)參數設置

基于本文提出的混合粒子群算法對多任務深空探測軌道進行優化。問題求解分為混合推力模式軌道設計和小推力軌道設計2個階段。

第1階段,設置粒子群大小為200,最大迭代次數為2 500,慣性權重最大值為0.9,最小值為0.4,學習因子c1的初始值為2.5,終止值為0.5,學習因子c2的初始值為0.5,終止值為2.5,α的初始值為2.5,終止值為0.5。

第2階段以第1階段的結果為初解,假設發動機采取bang-bang控制策略,即推力大小只能為0或最大值0.5 N,方向沿著速度方向或速度反方向。由于具有良好初解,所以設置粒子群大小為20,最大迭代次數為25,其他參數設置不變。

(2)結果分析

給定設計變量的取值范圍,仿真計算50次,得到設計變量優化結果的最優值,如表1所示。

表 1 設計變量優化結果

第2階段混合粒子群算法適應度收斂曲線如圖6所示。

小推力多任務行星探測高速軌道在日心坐標系下的示意圖如圖7所示。圖中圓形軌道由內而外分別為地球軌道、火星軌道和木星軌道,“*”代表探測器從地球離開時地球的空間位置,“?”代表探測器抵達火星借助引力輔助并對火星進行觀測時火星的空間位置,“Δ”代表探測器抵達木星借助引力輔助并拋下木星探測器時木星的空間位置;探測器的軌道中細實線代表小推力沿著速度方向持續工作,粗實線代表探測器發動機不工作。

探測器在木星完成引力輔助后開始逃離太陽系的深空探測任務,軌道示意圖見圖8。圖中圓軌道為木星軌道,探測器的軌道中細實線代表小推力沿著速度方向持續工作,粗實線代表探測器發動機不工作。

經過優化后,得到發動機的工作時序,發動機的工作狀態與探測器所處的狀態如表2所示。

表 2 發動機工作時序

(3)算法對比分析

問題求解分為混合推力模式軌道設計和小推力軌道設計兩個階段,為橫向比較算法的性能,分別采用文獻[18]所述遺傳算法和文獻[19]所述粒子群算法進行小推力軌道優化。

針對混合推力模式軌道設計階段,設置遺傳算法種群大小取200,代數設為2 500;粒子群算法種群大小取200,最大迭代次數取2 500。

仿真進行20次,對各種算法的平均精度和最高精度對比見表3。

表 3 優化性能對比

仿真結果表明,在求解星際探測軌道時,遺傳算法收斂精度較高,粒子群算法搜索速度較快,梯度混合粒子群算法收斂速度最快,GHPSO相對GA提高61.96%,相對PSO提高47.85% ;梯度混合粒子群算法精度最高,GHPSO相對GA平均精度提高了2.86%,最高精度提高了1.24%,相對PSO平均精度提高了4.19%,最高精度提高了3.97%。

針對小推力軌道設計階段,以GHPSO算法得到混合推力模式軌道設計的結果為初解,設置遺傳算法種群大小取20,代數設為25;粒子群算法種群大小取20,最大迭代次數取25。

仿真進行50次,選擇各算法精度最高的一次作為代表進行對比,遺傳算法、粒子群算法和梯度混合粒子群算法的適應度收斂曲線如圖9所示。

由圖9可知,由于具有良好的初解,各算法收斂較快,但在收斂速度和求解精度上仍然存在差異,表4為50次仿真計算結果的統計值。

仿真結果表明,在求解星際探測軌道時,遺傳算法收斂精度較高,粒子群算法搜索速度較快,梯度混合粒子群算法收斂速度最快,GHPSO相對GA提高61.87%,相對PSO提高43.66%;梯度混合粒子群算法精度最高,GHPSO相對GA平均精度提高了3.33%,最高精度提高了1.63%,相對PSO平均精度提高了4.72%,最高精度提高了3.02%。

本文提出的梯度混合粒子群算法,通過結合梯度搜索的特點使得粒子種群能跳出局部最優解,解決了粒子群算法易“早熟”的現象。仿真結果表明,該算法的全局尋優能力及收斂速度較傳統遺傳算法和粒子群算法都有了顯著提高。

表 4 優化性能對比

5 結論

(1)本文提出了新的“星際高速公路”概念,旨在通過尋找合適的發射窗口,合理利用行星引力場,縮短深空探測需要耗費的時間。

(2)以多任務深空探測為背景,建立了包含多模型、多約束、多變量、多階段的星際高速公路優化設計模型。

(3)通過融合推力工作模式、 分段設計目標函數等策略,克服了模型的內點約束限制,降低了優化模型的復雜度。

(4)為了提高優化速度、提升優化精度,結合梯度搜索和粒子群算法的特點,研究提出了一種梯度混合粒子群算法。

(5)橫向對比了本文算法與傳統粒子群算法和遺傳算法的優化性能,仿真結果表明,本文提出算法搜索速度快,第一設計階段GHPSO相對GA提高61.96%,相對PSO提高47.85%,第二設計階段GHPSO相對GA提高61.87%,相對PSO提高43.66%;精度最高,第一設計階段GHPSO相對GA平均精度提高了2.86%,最高精度提高了1.24%,相對PSO平均精度提高了4.19%,最高精度提高了3.97%,第二設計階段GHPSO相對GA平均精度提高了3.33%,最高精度提高了1.63%,相對PSO平均精度提高了4.72%,最高精度提高了3.02%。

(6)優化得到了發動機的工作時序,發動機從66 844.4 MJD開機,持續工作397.41個儒略日;在67 241.85 MJD關機,持續關機382.89個儒略日;在67 624.74 MJD再次開機直至燃料消耗完畢。

[1] 吳偉仁,劉旺旺,蔣宇平,等.國外月球以遠深空探測的發展及啟示(上)[J].中國航天,2011(7):9-12.

[2] 吳青.飛出太陽系[J].大自然探索,2011(4):20-24.

[3] Cui P Y,Qiao D,Cui H T,et al.Target selection and transfer trajectories design for exploring asteroid mission[J].Science China Technological Sciences,2010,53(4):1150-1158

[4] 李鑒,韓潮.小推力最優軌道轉移問題的UKF估計算法[J].宇航學報,2014,35(2):144-150

[5] 譚高威,高揚,楊新.深空探測器多次引力輔助轉移軌道全局搜索[J].航天器工程,2012,21(2):18-27.

[6] David C F,Mark W,Kathleen H,et al.Applications of multi-body dynamical environments:the ARTEMIS transfer trajectory design[J].Acta Astronautica,2012,73(12):237-249.

[7] 孟林智,黃江川,葉培建,等.嫦娥二號衛星多目標多任務設計與經驗[J].中國科學:技術科學,2013,43(6):585-595.

[8] 李俊峰,蔣方華.連續小推力航天器的深空探測軌道優化方法綜述[J].力學與實踐,2011,33(3):1-6.

[9] 尚海濱,崔平遠,徐瑞,等.基于高斯偽光譜的星際小推力轉移軌道快速優化[J].宇航學報,2010,31(4):1005-1011.

[10] 喬棟,崔平遠,徐瑞.星際探測借力飛行軌道的混合設計方法研究[J].宇航學報,2010,31(3):655-660.

[11] Hartmann J W,Coverstone V L,Williams S N.Optimal interplanetary spacecraft trajectories via Pareto genetic algorithm[J].The Journal of the Astronautical Sciences,1998,46(3):267-282.

[12] Mateen-ud-Din Q,He Lin-shu,Tarek E.Rapid trajectory optimization using computational intelligence for guidance and conceptual design of multistage space launch vehicles [C]//AlAA Guidance,Navigation,and Control Conference and Exhibit,San Francisco,California,2005:1-18.

[13] 冉茂鵬,王青.一種基于EPSO的航天器交會軌跡優化方法[J].宇航學報,2013,34(9): 1195-1201.

[14] 唐金國,羅亞中,雍恩米.航天器軌道優化理論、方法及應用[M].北京: 科學出版,2011:169-180.

[15] 趙國強.深空探測飛行任務軌道設計[D].北京:清華大學,2011.

[16] 戴光明,彭雷,羅治情.行星際脈沖轉移軌道設計與優化算法[M].武漢:中國地質大學出版社,2012.

[17] 汪定偉,王俊偉,王洪峰,等.智能優化方法[M].北京:高等教育出版社,2007.

[18] 鮮勇,許立軍.遺傳算法在導彈飛行程序設計中的應用研究[J].系統仿真學報,2009,21(5):1502-1504.

[19] 楊希祥,江振宇,張為華.基于粒子群算法的固體運載火箭上升段彈道優化設計研究[J].宇航學報,2010,31(5):1304-1309.

(編輯:呂耀輝)

Design of low-thrust multitask exploration orbit based on gravity-assist

LIN Hao-shen,HE Bing,ZHAO Xin,LIU Gang,AN Xi-bin

(Department of Space Engineering,The second Artillery Engineering University,Xi'an 710025,China)

Based on long time and high cost features of deep space exploration,an optimization model,containing multiple models,constraints,variables and multi-stage,was established.By integrating thrust mode and step-by-step design of objective function,the complexity of optimization model was reduced.In order to improve the speed and accuracy of optimization algorithms,a Gradient Hybrid Particle Swarm Optimization (GHPSO) algorithm was presented,which combined the characteristic of gradient search and Particle Swarm Optimization algorithm.This algorithm was applied to the trajectory optimization for multitask explorations to obtain engine's timing.Then the performance of GHPSO was compared with traditional particle swarm optimization(PSO)algorithm and genetic algorithm(GA).The results of simulation show that the GHPSO has faster convergence and higher accuracy. In the first stage of simulation,the convergence speed of GHPSO was 61.96% faster than GA and 47.85% than PSO,meanwhile,the average accuracy of GHPSO was 2.86% higher than GA and 4.19% than PSO,and the highest precision of GHPSO was 1.24% higher than GA and 1.63% than PSO.In the second stage,the convergence speed of GHPSO was 61.87% faster than GA and 43.66% than PSO,meanwhile,the average accuracy of GHPSO was 3.33% higher than GA and 4.72% than PSO,and the highest precision of GHPSO was 1.63% higher than GA and 3.02% than PSO.The result shows that the GHPSO can be applied to trajectory optimization design as well other nonlinear constrained global optimization problems.

gradient hybrid particle swarm optimization;multitasking interplanetary exploration;gravity-assist

2015-09-25;

2016-03-07。

基于量子生物地理優化的低空UAV編隊在線協同航跡規劃方法研究(61403399)。

林浩申(1992—),男,碩士,專業方向為空間軌道設計與優化、空間信息融合。E-mail:linhaoshen1@163.com

V412

A

1006-2793(2017)01-0121-07

10.7673/j.issn.1006-2793.2017.01.022

猜你喜歡
優化設計
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
何為設計的守護之道?
現代裝飾(2020年7期)2020-07-27 01:27:42
《豐收的喜悅展示設計》
流行色(2020年1期)2020-04-28 11:16:38
瞞天過海——仿生設計萌到家
藝術啟蒙(2018年7期)2018-08-23 09:14:18
設計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
主站蜘蛛池模板: 国产一级做美女做受视频| 久久99精品国产麻豆宅宅| 五月天在线网站| 国产成人综合久久精品尤物| 黄网站欧美内射| 国产区在线观看视频| 久久国产精品77777| 九九热这里只有国产精品| 亚洲首页在线观看| 中文字幕日韩欧美| 伊人久久久久久久久久| 亚洲色图欧美激情| 精品撒尿视频一区二区三区| 色妞www精品视频一级下载| aⅴ免费在线观看| 亚洲精品国产成人7777| 欧美国产日韩在线播放| 亚洲第一色视频| 亚洲av无码牛牛影视在线二区| 国产乱论视频| 欧美色图久久| 国产成人综合网在线观看| 国产xxxxx免费视频| 国产精品3p视频| 日韩欧美国产另类| 久久国产精品夜色| 亚洲午夜国产精品无卡| 麻豆精品在线视频| 亚卅精品无码久久毛片乌克兰 | 国产精品自在线拍国产电影| AV老司机AV天堂| av在线无码浏览| 国产高清无码麻豆精品| 伊人久久精品亚洲午夜| 欧美精品v日韩精品v国产精品| 久久香蕉国产线看观看式| 极品性荡少妇一区二区色欲| 自拍中文字幕| 国产一级毛片网站| 亚洲熟妇AV日韩熟妇在线| 国产午夜一级淫片| 欧美不卡在线视频| 亚洲无码高清一区| 99伊人精品| 国产乱人伦AV在线A| 精品人妻系列无码专区久久| 国产黄在线免费观看| 2020精品极品国产色在线观看 | 伊在人亞洲香蕉精品區| 巨熟乳波霸若妻中文观看免费| 欧美a级在线| 国产成人调教在线视频| 无码久看视频| 亚洲香蕉在线| 日韩天堂网| 精品成人一区二区三区电影| 动漫精品中文字幕无码| 国产欧美日韩在线一区| 91精品视频网站| 毛片久久网站小视频| 国产福利免费在线观看| 超碰色了色| 精品视频在线一区| 欧美在线三级| 三上悠亚一区二区| 国产97视频在线观看| 毛片a级毛片免费观看免下载| 久久国产高清视频| 久久婷婷综合色一区二区| 精品丝袜美腿国产一区| 欧美精品在线免费| 欧美成人第一页| 国产精品亚洲一区二区三区z| 在线免费亚洲无码视频| 久久人人97超碰人人澡爱香蕉| 日韩视频免费| 8090成人午夜精品| 免费人成黄页在线观看国产| 国产成人亚洲综合A∨在线播放| 美女一级毛片无遮挡内谢| 天天摸夜夜操| 午夜激情福利视频|