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

一種立方星編隊重構多脈沖機動規(guī)劃方法

2022-12-26 01:49:40劉幸川陳丹鶴廖文和
宇航學報 2022年11期
關鍵詞:規(guī)劃

徐 根,劉幸川,陳丹鶴,廖文和

(南京理工大學機械工程學院,南京 210094)

0 引 言

隨著冷氣推進[1]、電推進[2]等微推進系統(tǒng)向著低功耗、小尺寸、模塊化的方向發(fā)展,使得立方星具備了軌道機動控制的硬件基礎,拓展了低成本立方星的在軌應用領域。國內外多次開展了立方星執(zhí)行軌道機動控制的在軌演示,如加拿大多倫多大學的CAN-4/5雙星編隊任務[3]、Tyvak納衛(wèi)星公司的CPOD交會操作試驗[4]、Aerospace公司的AeroCube-10繞飛觀測試驗[5]等項目,驗證了立方星在復雜空間任務中具備較廣闊的應用潛力,如編隊飛行等任務。

目前,隊形構建、重構和維持過程中的機動規(guī)劃方法仍然是編隊飛行應用的關鍵問題。諸多學者已經對相對運動建模、機動規(guī)劃以及攝動優(yōu)化等方面開展了研究。

針對近圓軌道上的編隊飛行問題,HCW方程提供了最經典的相對運動模型,廣泛應用于相對運動機動規(guī)劃與控制中[6-7]。但是以相對位置速度矢量為狀態(tài)量不能直接反應相對運動軌跡的幾何特性,不利于編隊隊形設計和任務規(guī)劃。H?rting等在靜止軌道衛(wèi)星防碰撞問題中采用了基于相對偏心矢量和相對傾斜矢量的建模方法,D’Amico等[8]在此基礎上提出了相對軌道根數(ROE)模型,應用于低軌衛(wèi)星編隊任務的GNC算法中,并在GRACE[9],TanDEM-X/TerraSAR-X[10],PRISMA[11],AVANTI[12]等任務中獲得了驗證。ROE模型不僅提供了相對軌道運動的幾何解釋,并且基于相對偏心矢量和相對傾斜矢量提供了一種簡便的被動安全編隊隊形的設計方法[13]。Gaias團隊[14-15]分析了J2攝動和大氣阻力攝動對低軌編隊的影響,構建了高精度的相對運動模型,并考慮到星載計算約束,將差分氣動阻力的影響轉換為相對軌道高度的固定漂移,獲得適用于星載計算的ROE狀態(tài)轉移模型[16]。

對于星間距離較大的交會、編隊構建和重構等任務,通常采用基于脈沖推力的機動規(guī)劃方法。孟云鶴等[17]在HCW方程的基礎上,分析了脈沖機動與相對運動之間的關系,研究了燃耗最優(yōu)的組合脈沖機動求解方法。根據相對軌道運動的特性,這些研究中先將面內分量與面外分量的控制解耦,再研究面內各分量之間的控制耦合關系,從而獲得脈沖機動的解析解。而Gaias等[18]、Di Mauro等[19]、Lim等[20]以及Liu等[21]基于ROE模型,從幾何角度闡述了脈沖機動對相對軌道的影響,針對面內各分量的控制耦合問題,給出了多種脈沖組合下的二/三/四脈沖解析解。

上述機動規(guī)劃相關的研究重點,是獲得燃耗最優(yōu)的脈沖機動的解析解,通常不考慮機動規(guī)劃的工程約束。而對于立方星的軌道機動任務,往往會由于姿態(tài)穩(wěn)定控制飽和或推進系統(tǒng)工作時長約束等問題,使得立方星單次機動的速度增量受到約束。對于較大范圍的編隊隊形重構任務,無約束下的脈沖機動解析解往往超出立方星的軌道機動能力,因此立方星軌道控制算法需要具備在機動能力約束下可靠的機動規(guī)劃能力。

針對近地軌道上的立方星編隊飛行面內隊形重構問題,本文提出了一種適用于立方星星載計算的多脈沖機動規(guī)劃算法。首先,梳理了基于相對軌道根數的相對運動理論基礎,包括相對軌道根數定義、考慮近地軌道J2攝動和大氣阻力攝動的線性化狀態(tài)轉移模型、以及脈沖機動對相對軌道根數的控制模型。其次,基于相對運動面內分量和面外分量解耦的特性,提出燃耗最優(yōu)的多脈沖機動規(guī)劃算法,解決了規(guī)劃過程中的單次機動速度增量約束,并提出迭代優(yōu)化策略,提高機動規(guī)劃終點精度。最后,基于“田園一號”立方星平臺開展了編隊面內構型重構任務的仿真驗證和分析。

1 基于相對軌道根數的運動模型

1.1 相對軌道根數定義

基于經典軌道根數α=(a,e,i,ω,Ω,M)T,定義在目標星參考系下的無量綱相對軌道根數[8]:

(1)

式中:u=M+ω為平均緯度輻角,下標c和d分別代表目標星和機動星(為簡化描述,后續(xù)公式中不帶下標的絕對軌道根數均為目標星參數)。相對軌道根數δα六個分量可分為軌道面內分量和軌道面外分量,面外分量為相對傾斜量δi,aδi描述了相對運動的法向振幅和初相位,面內分量包括相對高度δa、相對平均緯度幅角δλ以及相對偏心量δe,其中δa和aδe描述了相對運動軌跡在軌道面內投影的形狀,包括平面內運動的振幅、初相位、徑向偏移和跡向距離漂移速度,aδλ描述了兩星之間的跡向距離。

(2)

在二體動力學下的ROE運動模型可表達為:

δα(t)=ΦHCW(t,t0)δα(t0)

(3)

其中狀態(tài)轉移矩陣可表達為:

(4)

1.2 攝動影響

對于近地軌道運行的衛(wèi)星,軌道攝動的主要來源是地球非球形重力場的二階緯向分量J2攝動和大氣阻力攝動。本節(jié)對J2攝動和大氣阻力攝動的影響進行分析和建模。

1.2.1J2攝動影響

通過將J2攝動對絕對軌道根數的影響投影到ROE空間,并通過保留一階項的方式完成線性化,獲得包含J2攝動影響的ROE狀態(tài)轉移矩陣[14]:

ΦHCW+J2(t,t0)=

(5)

式中:

(6)

式中:RE為地球半徑。由式(5)可知,在J2攝動作用下,相對偏心量的方向將會發(fā)生旋轉,旋轉角速度主要取決于目標星軌道的高度和傾角,旋轉方向取決于目標星軌道的傾角。此外,衛(wèi)星的平均運動速度和升交點赤經都會受J2攝動影響產生漂移,從而對編隊衛(wèi)星之間的跡向距離產生了影響,跡向距離漂移的方向取決于目標星軌道傾角和兩星相對高度。

1.2.2大氣阻力攝動影響

低軌稀薄大氣產生的阻力會使衛(wèi)星產生沿負切向的加速度。由于彈道系數和面質比的差異,編隊衛(wèi)星之間存在大氣阻力加速度的差異[22]:

(7)

其中ρ為大氣密度,v=na為目標星的平均運動速度,S為衛(wèi)星迎風面積,m為衛(wèi)星質量,CD為衛(wèi)星彈道系數。

(8)

(9)

其中大氣阻力攝動對ROE的影響表達為:

Φdrag=

(10)

1.3 ROE脈沖控制模型

若機動星執(zhí)行脈沖機動Δv=[ΔvR,ΔvT,ΔvN]T,產生的ROE變化為[18]:

(11)

由式(11)可以看到,ROE的面內分量δa,δλ,δe和面外分量δi的控制是解耦的。對于法向相對運動,可通過在特定相位施加一次法向脈沖實現對δi的控制:

(12)

而徑向/跡向平面內的相對運動顯然存在控制耦合問題:一方面各面內分量之間存在控制耦合,一次脈沖會對多個分量產生影響;另一方面切向脈沖和徑向脈沖的控制效果存在重疊,兩個方向的脈沖控制均會對δe產生影響,但切向脈沖的控制效果是徑向脈沖的控制效果的2倍。

2 相對運動面內重構規(guī)劃算法

由于徑向脈沖無法控制δa,且考慮到燃耗最優(yōu)問題,因此本文所提出的機動規(guī)劃方法僅使用切向脈沖。針對面內各分量之間的耦合關系,由于切向脈沖對δa和δe的控制為瞬時變化,而對δλ的影響需要經過時間的累積。因此在規(guī)劃過程中,按照“先規(guī)劃相對形狀控制機動,后規(guī)劃跡向距離修正機動”的策略進行求解,先完成控制δa和δe的機動規(guī)劃,再進行修正δλ的機動規(guī)劃。

2.1 機動約束下的多脈沖機動規(guī)劃算法

對于相對軌道面內分量的控制,由于單次切向脈沖對δa和δe的控制存在耦合,可通過兩次執(zhí)行相位相差180°的切向脈沖實現對δa和δe的分離控制:

(13)

兩次脈沖機動的執(zhí)行相位為

(14)

由于立方星體積和質量的約束,所配備的推進系統(tǒng)和姿態(tài)控制系統(tǒng)往往結構簡單、控制能力有限。一方面,微推力器的推力較小,無法在較短時間內提供較大的速度增量;另一方面,由于推力偏心,會使姿態(tài)穩(wěn)定控制飽和,推進系統(tǒng)單次工作時間受到限制,使得立方星單次軌道機動存在最大速度增量約束。此外,在兩次軌道機動執(zhí)行之間,立方星需要進行動量輪卸載、微推進器調整等工作,使得兩次機動之間存在最小時間間隔的約束。因此,在立方星編隊任務的機動規(guī)劃中必須考慮速度增量與機動時間間隔約束。根據式(13)得到脈沖速度增量往往會超出立方星的單次機動能力。

而考慮到在二體運動下,δa和δe不隨時間發(fā)生改變,根據式(11)可得:

(15)

一次切向脈沖對δa和δe的控制效果等價于由多次在同一相位執(zhí)行的切向脈沖控制效果的線性疊加。因此針對單次機動速度增量約束,可將根據式(13)得到的脈沖機動分解為多個速度增量較小的脈沖機動逐軌執(zhí)行。

(16)

并為分解后的機動逐軌分配執(zhí)行時間:

(17)

(18)

將式(18)的形式轉換為:

aΔδλΔvT=-3ΔvT(tf-tΔvT)

(19)

根據式(19)所建立的切向脈沖機動對跡向距離的控制模型,在施加M1+M2次相對形狀控制機動后,跡向距離產生的變化為:

(20)

(21)

同樣考慮到單次機動速度增量的約束,跡向距離修正機動通過多組脈沖組合實現:

(22)

(23)

(24)

顯然,完成面內分量控制的所有機動需要滿足:

(25)

(26)

(27)

實際任務中期望以較少的機動次數完成任務,因此每一次機動的速度增量都應接近于Δvmax。根據式(27)的目標函數,可將燃耗最優(yōu)問題轉換為使規(guī)劃結果中的脈沖數量最少。注意到:

(28)

可得:

(29)

顯然,機動規(guī)劃結果中∑Qi應為定值。為了使N的值最小,應當使每一組的兩次脈沖執(zhí)行時間的間隔QiT盡可能最大。因此在跡向距離修正脈沖執(zhí)行時間分配時,需要在滿足式(24)的約束下,按時間間隔從大到小依次搜索可行解,圖1給出了機動執(zhí)行時間的分配過程示意。

圖1 多脈沖軌道機動規(guī)劃執(zhí)行時間分配示意Fig.1 Schematic diagram of execution time distribution in multi-impulse orbit maneuver planning

上述面內機動的多脈沖求解算法可總結為:

算法1.機動約束下的多脈沖機動規(guī)劃算法

輸入:目標星軌道根數αc,初始相對軌道δα0,終點目標相對軌道δαf,初始時間t0,終點時間tf,單次機動最大速度增量Δvmax,兩次機動之間最短時間間隔Δtmin。

4) 根據式(25)和(29),確定跡向距離修正需求;

根據上述多脈沖機動規(guī)劃算法,最終共獲得了M1+M2+2N次脈沖機動。由于這些機動需要逐軌執(zhí)行,因此對任務時間提出了要求:首先,任務時間需要滿足相對形狀控制機動的需要,總時間至少要大于M1軌或M2軌;其次,在任務時間內需要搜索到足夠多滿足最小時間間隔約束的Qi,滿足跡向距離修正機動的需要。

2.2 機動規(guī)劃的迭代優(yōu)化策略

對于機動能力受限的立方星來說,執(zhí)行較大范圍的編隊隊形重構任務往往需要較長的時間,而算法1里的機動規(guī)劃是基于二體運動獲得的,對于較長時間的任務,軌道攝動產生的誤差較為明顯,不可忽略。需要對機動規(guī)劃進行優(yōu)化,降低機動規(guī)劃的終點誤差。

(30)

此外,注意到式(9)中的狀態(tài)轉移矩陣僅與時間有關,與狀態(tài)變量無關,因此具備以下性質:

Φ(tk,ti)=Φ(tk,tj)Φ(tj,ti)

(31)

因此,若在中間時間tm引入狀態(tài)變化后,得到的狀態(tài)轉移方程可表達為:

X(t)=Φ(t,tm)[Φ(tm,t0)X(t0)+ΔX(tm)]=

Φ(t,t0)X(t0)+Φ(t,tm)ΔX(tm)

(32)

由式(32)可以看到,終點狀態(tài)可被分解為初始狀態(tài)的自然轉移和每次狀態(tài)變化的自然轉移的線性組合。因此依次執(zhí)行根據算法1獲得的脈沖機動序列后,在終點時刻的狀態(tài)量為:

(33)

基于式(33)的線性化終點時刻狀態(tài)計算公式,給出對機動規(guī)劃進行迭代優(yōu)化的步驟:

算法2.考慮攝動影響的機動規(guī)劃迭代優(yōu)化算法

輸入:算法1的輸入變量,終點誤差范圍,最大迭代次數。

2) 根據算法1,獲取脈沖機動序列;

3) 根據式(33)計算機動執(zhí)行的終點位置δα(tf);

5) 重復步驟2),直到滿足終點誤差,或達到最大迭代次數。

3 仿真校驗

為了校驗本文所提出的多脈沖機動方法的有效性,根據“田園一號”立方星的相關參數進行任務仿真,在星箭分離約7日后開始進行跟隨編隊構建任務。任務目標為構建在目標星跡向后方50 km處跟隨飛行的隊形,并降低兩星的相對偏心率,機動規(guī)劃算法需要完成從初始隊形到目標隊形的重構任務。以搭載同一運載火箭入軌的另一顆立方星“金紫荊二號”為目標星,以“田園一號”立方星為機動星,兩顆衛(wèi)星的彈道系數相同,均取CD=2.2。主要仿真參數如表1所示:

表1 主要仿真參數Table 1 Main simulation parameters

“田園一號”立方星所使用的冷氣微推力器的標稱推力為3.6 mN,采用零動量輪系進行姿態(tài)三軸穩(wěn)定控制,動量輪飽和后,需要使用三軸磁力矩器進行動量卸載。單次軌道機動最長時間為4 min,即微推力器單次工作最大速度增量約為0.085 m/s,微推力器兩次工作之間時間間隔至少為20 min。

結合“田園一號”立方星的姿軌控能力,并考慮到低軌衛(wèi)星軌道運行周期,為了便于搜索跡向修正機動的可執(zhí)行時間,設定跡向距離修正機動的執(zhí)行相位取決于相對形狀控制機動的執(zhí)行相位:

(34)

(35)

其中可執(zhí)行機動的時間點需要滿足時間約束:

(36)

在除去執(zhí)行相對形狀控制機動的時間點后,每個相位可對應Ji個可用于執(zhí)行跡向距離修正機動的時間點,將可行時間點一前一后依次組合,得到間隔時間從大到小的N#i個脈沖機動組合:

(37)

考慮到算法1中的任務時間要求,由式(37)獲得的脈沖機動組合需要滿足:

(38)

若式(38)的條件不能得到滿足,則總任務時間過短,無法到達目標跡向位置,需要增加任務時間。若式(38)的條件滿足,則對Q#進行排序,從大到小依次選取Qi,直至滿足式(28)的條件,從而完成對跡向修正機動的規(guī)劃。

3.1 多脈沖機動規(guī)劃算法仿真校驗

根據任務仿真參數設置,軌道機動需要實現的相對軌道高度與相對偏心量變化量為:

計算得到無約束下的相對形狀控制機動:

圖3 任務2:機動約束下在60軌內完成隊形重構Fig.3 Mission 2:Complete formation reconfiguration within 60 orbits with constraints

顯然,“田園一號”立方星無法僅通過兩次脈沖機動完成對相對偏心量和相對軌道高度的控制,需要根據算法1獲得滿足約束的機動序列。

對相對形狀控制機動進行分解:

在執(zhí)行完所有機動序列后,相對運動軌跡如圖2(c)所示。在增加跡向距離修正機動后,機動終點可以達到目標跡向相對位置。針對在40軌內完成隊形重構任務,算法1共規(guī)劃了87次脈沖機動,總速度增量需求約為7.250 m/s。

若增加任務時間,在60軌內完成隊形重構任務,仿真結果如圖3所示。可以看到通過調整執(zhí)行時間,僅執(zhí)行相對形狀控制機動即可到達距離目標終點較近的位置,僅需增加一組速度增量較小的跡向距離修正機動即可。與任務1相比,任務2減少了跡向距離修正機動的燃耗,總速度增量需求約為5.443 m/s。

3.2 迭代優(yōu)化策略仿真校驗

基于任務1所獲得的脈沖機動序列,分別根據不同的狀態(tài)轉移模型計算終點相對軌道。其中,從星載機動規(guī)劃的角度,取任務時間內的平均大氣密度為6.569×10-13kg/m3(基于NRLMSISE-00大氣密度模型計算)。

表2 攝動影響下的終點相對軌道誤差Table 2 Terminal ROE errors under perturbation

根據算法2對任務1的編隊重構任務的機動規(guī)劃列進行迭代優(yōu)化,過程如表3所示,經過3次迭代,即可將面內分量的終點誤差都降低至1 m以內。

表3 迭代優(yōu)化過程Table 3 Process of the iterative optimization

表4 迭代優(yōu)化結果對比Table 4 Comparison of the iterative optimization results

圖4 迭代優(yōu)化后的相對偏心量機動軌跡Fig.4 Trajectory of relative eccentricity vector after iterative optimization

表4中對算法1和算法2的機動規(guī)劃結果進行了對比。經迭代優(yōu)化后,機動規(guī)劃的燃耗需求有了較顯著的降低,主要源自相對形狀控制機動的燃耗需求降低。

算法2的優(yōu)化變量實際上是算法1輸入值中的目標相對軌道δαf。如圖4所示,經過迭代優(yōu)化后,算法1需要控制的相對偏心量發(fā)生了改變:

與任務目標相比,經迭代優(yōu)化后,算法1的輸入值中Δδe大小變小,因此降低了機動規(guī)劃的燃耗需求,同時Δδe方向的變化也使機動的執(zhí)行相位發(fā)生了變化。由此可見,在規(guī)劃中考慮J2攝動會主要影響相對偏心量的控制,目標星軌道的傾角會決定相對偏心量的漂移方向和漂移速率,對燃耗的影響需要結合目標星軌道的半長軸、傾角和編隊隊形重構的具體參數進行定量分析。

3.3 機動規(guī)劃算法在復雜攝動下的適應性分析

為了分析算法2的迭代優(yōu)化策略的效果,圖4中的相對偏心量軌跡是通過式(33)的狀態(tài)轉移模型獲得的。本節(jié)基于高精度軌道遞推[23]和瞬/平根數轉換[24]進行機動規(guī)劃算法的驗證,仿真驗證的流程如圖5所示。

圖5 仿真驗證流程Fig.5 Scheme of algorithm verification

圖5中虛線部分表示高精度軌道遞推,包括21×21階的地球引力場模型、大氣阻力模型以及日月第三體引力。仿真初始參數設置如表5所示。

表5 仿真初始參數Table 5 Initial simulation parameters

分別采用算法1和算法2進行機動規(guī)劃,進行仿真。得到如表6所示的終點相對軌道誤差,相對偏心率控制軌跡如圖6所示。

表6 兩種機動規(guī)劃算法的終點相對軌道誤差Table 6 Terminal ROE error of two planning algorithms

圖6 高精度軌道遞推下的相對偏心量機動軌跡Fig.6 Trajectory of relative eccentricity vector under numerical propagation with high-fidelity perturbation

與算法1相比,算法2中的狀態(tài)轉移模型能夠較準確地體現地球非球形攝動的影響,可以預估相對偏心量的方向漂移并加以修正,使得終點位置的相對偏心量與目標相對偏心量方向一致,有利于維持編隊在徑向/法向平面上的投影。

在實際在軌任務中,考慮到測定軌誤差和軌道機動執(zhí)行誤差等因素,不能僅依靠初始時刻的機動規(guī)劃進行完全的開環(huán)控制,需要在過程中修正機動規(guī)劃,因此算法1和算法2在實際在軌任務中均可使用。算法1的計算量較小,在機動過程中可以增加多次重規(guī)劃;而算法2在迭代的過程中需要重復運行算法1,計算量是算法1的4~6倍,但算法2的規(guī)劃誤差更小,可以降低機動重規(guī)劃的頻率。

4 結 論

針對軌道機動能力約束下的編隊隊形面內重構機動規(guī)劃問題,本文基于相對軌道根數提出了一種簡便的多脈沖機動規(guī)劃算法和迭代優(yōu)化策略,結合仿真校驗可得如下結論:

(1) 本文所提出的多脈沖機動算法可有效解決速度增量約束和時間間隔約束,能可靠地獲得燃耗最優(yōu)的機動序列。面內隊形重構任務僅需執(zhí)行切向脈沖,可簡化立方星結構,無需安裝徑向推力器。

(2) 結合考慮J2攝動和大氣阻力攝動的線性化相對狀態(tài)轉移模型,本文所提出的迭代優(yōu)化策略可快速有效地降低攝動影響下的終點誤差。

(3) 攝動因素會影響機動的燃耗需求,本文中對攝動影響的分析可為編隊隊形重構的參數設計提供參考。

鑒于相對軌道面內面外運動解耦的特性,本文所提的算法和任務仿真僅針對面內分量控制,后期將進行考慮面外分量控制的綜合規(guī)劃策略研究,提高該方法在立方星編隊任務中的工程適用性。

猜你喜歡
規(guī)劃
我們的規(guī)劃與設計,正從新出發(fā)!
房地產導刊(2021年6期)2021-07-22 09:12:46
“十四五”規(guī)劃開門紅
“十四五”規(guī)劃建議解讀
發(fā)揮人大在五年規(guī)劃編制中的積極作用
規(guī)劃計劃
規(guī)劃引領把握未來
快遞業(yè)十三五規(guī)劃發(fā)布
商周刊(2017年5期)2017-08-22 03:35:26
基于蟻群算法的3D打印批次規(guī)劃
多管齊下落實規(guī)劃
十三五規(guī)劃
華東科技(2016年10期)2016-11-11 06:17:41
主站蜘蛛池模板: 久久这里只有精品66| 女人18毛片一级毛片在线 | 国产一区在线观看无码| 亚洲制服中文字幕一区二区| 亚洲二区视频| 伊人精品视频免费在线| 国禁国产you女视频网站| 国产成人精品一区二区秒拍1o| 久久久噜噜噜久久中文字幕色伊伊| 天天爽免费视频| 美女无遮挡拍拍拍免费视频| 制服丝袜亚洲| 久久窝窝国产精品午夜看片| 国产成人精品一区二区三区| 伊人成人在线视频| 欧美亚洲激情| 国产成人区在线观看视频| 欧洲亚洲一区| 免费中文字幕在在线不卡| 亚洲国产综合自在线另类| 99爱视频精品免视看| 欧美国产日韩在线| 日韩一区精品视频一区二区| 五月激情婷婷综合| 国产主播一区二区三区| 久久国产V一级毛多内射| 日本草草视频在线观看| 久久久久国产精品免费免费不卡| 天天综合网色| 日本欧美在线观看| 日日拍夜夜操| 又爽又大又光又色的午夜视频| 久久99国产视频| 免费va国产在线观看| 女同久久精品国产99国| 日韩成人在线网站| 久久综合婷婷| 四虎影视无码永久免费观看| 99这里精品| 久久亚洲综合伊人| 午夜日韩久久影院| 91久久偷偷做嫩草影院| 999在线免费视频| 免费一级毛片在线观看| 伊人色天堂| 日本尹人综合香蕉在线观看| 成人精品亚洲| 99久久精品国产麻豆婷婷| 伊人久久大线影院首页| 欧美成人午夜在线全部免费| 免费不卡在线观看av| 久久这里只精品热免费99| 亚洲精品777| 亚洲欧美国产视频| 一级福利视频| 国产在线观看第二页| 欧美综合成人| 国产精品人成在线播放| 国产欧美综合在线观看第七页| 男女男免费视频网站国产| 91外围女在线观看| 日本免费精品| 中文字幕无码制服中字| 被公侵犯人妻少妇一区二区三区| 99久久精品国产自免费| 91啪在线| 天天综合网色| 久操中文在线| 成人无码一区二区三区视频在线观看| 在线观看亚洲国产| 欧美伦理一区| 亚洲永久精品ww47国产| 少妇精品久久久一区二区三区| 2020国产精品视频| www.亚洲色图.com| 国产精品久久久免费视频| 国内精品九九久久久精品| 欧美第二区| 最新亚洲人成网站在线观看| 亚洲一本大道在线| 欧美色99| 色婷婷亚洲十月十月色天|