• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    行星動力下降多約束最優(yōu)反饋制導(dǎo)方法

    2023-09-22 12:54:46張曉文王天舒王大軼
    宇航學(xué)報 2023年8期
    關(guān)鍵詞:推進(jìn)劑制導(dǎo)校正

    張曉文,王天舒,王大軼,李 驥

    (1. 清華大學(xué)航天航空學(xué)院,北京 100084;2. 北京控制工程研究所,北京 100190;3. 北京空間飛行器總體設(shè)計部,北京 100094)

    0 引 言

    未來行星著陸探測任務(wù)普遍提出定點著陸的工程目標(biāo),一般要求落點精度為100 m。實現(xiàn)此目標(biāo)主要依靠動力制導(dǎo)。以往月球和火星著陸器的動力制導(dǎo)方法研究和工程應(yīng)用多采用多項式制導(dǎo),例如阿波羅制導(dǎo)和E制導(dǎo)。阿波羅制導(dǎo)的優(yōu)點是終端三維位置、速度和加速度全部受控,缺點是固定終端時刻后未對推進(jìn)劑消耗作優(yōu)化。Lu[1-2]為改進(jìn)阿波羅制導(dǎo),先是參考隱式制導(dǎo)在原制導(dǎo)律中增加了一個控制增益作為可調(diào)參數(shù),之后應(yīng)用分?jǐn)?shù)多項式理論又增加了分?jǐn)?shù)指數(shù)參數(shù),并將新方法命名為分?jǐn)?shù)多項式動力下降制導(dǎo)。阿波羅制導(dǎo)中發(fā)動機(jī)推力加速度為時間的二次多項式形式。經(jīng)典的E制導(dǎo)及與其等同的基本型ZEM/ZEV制導(dǎo)[3]的推力加速度的制導(dǎo)指令為時間的一次多項式形式。E制導(dǎo)已被證明是常值引力場假設(shè)下的固定終端時刻的加速度模平方積分最優(yōu)制導(dǎo)[4],即能量最優(yōu)制導(dǎo),雖然不能完全等同于推進(jìn)劑消耗最優(yōu),但是已經(jīng)對推進(jìn)劑消耗做了顯著優(yōu)化。E制導(dǎo)的一項缺點是終端三維加速度不完全受控。E制導(dǎo)和阿波羅制導(dǎo)都未考慮推力范圍限制,因此實際應(yīng)用中指令推力可能觸發(fā)推力飽和。

    E制導(dǎo)或基本型ZEM/ZEV制導(dǎo)盡管是最優(yōu)反饋制導(dǎo),但是面對實際工程應(yīng)用時還是存在一些不足,因此對其改進(jìn)的方法不斷被提出。這些改進(jìn)包括完善最優(yōu)制導(dǎo)時長確定方法,提高制導(dǎo)算法的魯棒性和碰撞規(guī)避能力等。Guo等[5]使用數(shù)值優(yōu)化程序GPOPS選擇中間航路點作為ZEM/ZEV制導(dǎo)的中間目標(biāo),實現(xiàn)了火星動力下降約束下滑角的碰撞規(guī)避。Zhou等[6]針對火星動力下降段制導(dǎo)規(guī)避碰撞地平問題,在ZEM/ZEV制導(dǎo)原性能指標(biāo)中添加了加權(quán)的高度項,然后推導(dǎo)出有地平碰撞規(guī)避能力的改進(jìn)制導(dǎo)律。Wibben等[7]為了提高ZEM/ZEV制導(dǎo)在火星動力下降段的干擾適應(yīng)性,加入最優(yōu)滑動制導(dǎo)并證明了新方法的全局穩(wěn)定性。Zhang等[8]利用平均加速度修正基本型ZEM/ZEV制導(dǎo)加速度的來解決月面著陸地平碰撞規(guī)避問題,還應(yīng)用模型預(yù)測靜態(tài)規(guī)劃(MPSP)方法迭代計算滿足終端狀態(tài)約束的解。Zhang等[9]對Zhou等[6]的方法繼續(xù)改進(jìn),將性能指標(biāo)中高度項的常值加權(quán)系數(shù)改為與高度平方相關(guān)的時變系數(shù),進(jìn)一步提升碰撞規(guī)避能力。Ahn等[10]研究了軌道攔截或交會機(jī)動問題的由最小能量和最小時間組成的性能指標(biāo),利用最優(yōu)控制理論和橫截條件,推導(dǎo)出確定ZEM/ZEV制導(dǎo)最優(yōu)制導(dǎo)時長的代數(shù)方程。Wang等[11]提出一種適用于火星精確著陸的兩段ZEM/ZEV反饋制導(dǎo)方法,通過虛擬終端速度和在線參數(shù)搜索來提高制導(dǎo)的障礙規(guī)避能力。高峰等[12]利用脈寬脈頻(PWPF)調(diào)制器將ZEM/ZEV制導(dǎo)的連續(xù)指令轉(zhuǎn)換為常值推力器的開關(guān)指令。

    碰撞規(guī)避是行星著陸制導(dǎo)必須具備的能力。除了上面提到方法外,還有一種解決思路是控制著陸軌跡的曲率。Cui等[13]先是提出了一種基于多項式的曲率制導(dǎo)策略,通過調(diào)整曲率使著陸器沿凸形軌跡著陸,降低障礙物碰撞的風(fēng)險。之后進(jìn)一步改進(jìn),將曲率約束凸化為二階錐約束,并利用非線性規(guī)劃工具求解凸規(guī)劃問題來優(yōu)化推進(jìn)劑消耗[14]。

    近年來數(shù)值優(yōu)化方法和人工智能技術(shù)被越來越多的應(yīng)用于行星著陸制導(dǎo)問題的理論研究中。林曉輝等[15]應(yīng)用凸優(yōu)化理論將月球動力下降段軌道優(yōu)化問題轉(zhuǎn)化為二階錐問題,并采用內(nèi)點法求解了最優(yōu)標(biāo)稱軌跡。鄧雁鵬等[16]應(yīng)用序列凸優(yōu)化方法解決月面動力下降過程中時變加速度剖面實時估計和補(bǔ)償問題。Furfaro等[17]將深度強(qiáng)化學(xué)習(xí)方法與ZEM/ZEV制導(dǎo)方法相結(jié)合,通過動態(tài)調(diào)節(jié)2個制導(dǎo)增益和制導(dǎo)時長來滿足行星著陸問題的路徑等約束并優(yōu)化推進(jìn)劑消耗。數(shù)值優(yōu)化制導(dǎo)計算量大,依賴高效穩(wěn)定的求解器。人工智能制導(dǎo)目前處于理論研究階段,工程應(yīng)用同樣面臨計算量大和算法穩(wěn)定性方面的挑戰(zhàn)。

    總體而言,傳統(tǒng)最優(yōu)反饋制導(dǎo)方法有著計算量小、穩(wěn)定性好、節(jié)省推進(jìn)劑的優(yōu)勢,但是僅滿足初始和終端位置速度約束。實際工程應(yīng)用中還需要考慮推力范圍受限、連續(xù)變推力、碰撞規(guī)避和終端加速度等多方面約束。單臺變推力發(fā)動機(jī)推力的大小和方向在物理上只能連續(xù)變化,因此要求制導(dǎo)的指令推力連續(xù)變化。著陸終端希望著陸器為豎直姿態(tài),并且合加速度接近為零,這就對制導(dǎo)的終端加速度提出了約束?,F(xiàn)有最優(yōu)反饋制導(dǎo)方法研究欠缺對上述約束的綜合考慮。

    本文考慮推力范圍受限、連續(xù)變推力、碰撞規(guī)避和終端加速度共4方面約束,對動力下降應(yīng)用的傳統(tǒng)最優(yōu)反饋制導(dǎo)方法進(jìn)行改進(jìn)和完善。內(nèi)容安排為,首先建立行星定點著陸動力下降制導(dǎo)問題模型,然后介紹和分析E制導(dǎo)和阿波羅制導(dǎo)方法,接著詳細(xì)闡述提出的多約束最優(yōu)反饋制導(dǎo)方法,最后通過數(shù)學(xué)仿真驗證所提方法有效性。

    1 問題模型

    1.1 制導(dǎo)坐標(biāo)系定義

    制導(dǎo)坐標(biāo)系原點O定義在目標(biāo)著陸點,從原點指向天向定義為Z軸,垂直Z軸從原點指向著陸器初始位置方向為X軸,Y軸按照右手螺旋法則確定。在該坐標(biāo)系下,XOZ平面即為理想著陸軌跡面,或又稱為縱向平面,Y軸即為軌跡面負(fù)法線方向。X軸位移體現(xiàn)著陸器航程變化,因此可以將X向稱之為航向,Z軸位移體現(xiàn)著陸器高度變化,Y軸位移體現(xiàn)著陸器垂直軌跡面位置變化,因此可以將Y向稱之為橫向。

    1.2 行星動力下降動力學(xué)方程

    月球沒有大氣?;鹦谴髿庀”?在著陸器動力下降段氣動力顯著小于制導(dǎo)推力。因此以往分析行星動力下降制導(dǎo)問題時忽略大氣影響,僅考慮中心天體引力和發(fā)動機(jī)推力。制導(dǎo)坐標(biāo)系下的行星動力下降動力學(xué)方程如下:

    (1)

    (2)

    a=F/m

    (3)

    (4)

    1.3 定點著陸制導(dǎo)問題

    定點著陸制導(dǎo)問題即控制著陸器以指定終端速度到達(dá)指定終端位置。研究定點著陸制導(dǎo)問題時一般將終端位置和速度取為零。定點著陸制導(dǎo)問題屬于軌跡優(yōu)化問題,軌跡優(yōu)化問題又屬于最優(yōu)控制問題。設(shè)初始時刻為0,下面依次列出定點著陸制導(dǎo)問題的初始條件、終端條件、性能指標(biāo)和推力約束表述。

    初始條件:

    r(0)=r0

    (5)

    v(0)=v0

    (6)

    m(0)=m0

    (7)

    終端條件:

    r(tf)=rf

    (8)

    v(tf)=vf

    (9)

    性能指標(biāo):

    (10)

    推力約束:

    Fmin≤F≤Fmax

    (11)

    2 多項式制導(dǎo)

    多項式制導(dǎo)即制導(dǎo)指令加速度為時間的多項式形式。常見的多項式制導(dǎo)有E制導(dǎo)和阿波羅制導(dǎo)。下面對E制導(dǎo)和阿波羅制導(dǎo)展開原理介紹和性質(zhì)分析。在做性質(zhì)分析時假設(shè)引力加速度為常值。

    2.1 E制導(dǎo)

    在E制導(dǎo)方法中,制導(dǎo)推力加速度三維矢量取為時間的一次多項式,包含2個待確定系數(shù)矢量,即指令推力加速度為

    a=c0+c1t,t∈[0,tf]

    (12)

    將式(12)代入動力學(xué)方程(2)中,然后連續(xù)積分兩次,并代入位置速度的初值和終值確定制導(dǎo)律系數(shù),結(jié)果如下:

    (13)

    (14)

    將式(13)代入式(12)中后,得到E制導(dǎo)的初始時刻推力加速度計算公式:

    (15)

    定義零控位移偏差rzem和零控速度偏差vzev如下:

    (16)

    vzev?vf-v0-tfg

    (17)

    代入式(15)中,得到零控位移偏差rzem和零控速度偏差vzev形式的制導(dǎo)指令推力加速度計算公式,即基本型ZEM/ZEV制導(dǎo)律公式:

    (18)

    上式說明E制導(dǎo)與基本型ZEM/ZEV制導(dǎo)完全相同。將零控位移偏差rzem和零控速度偏差vzev前的常值控制增益6和-2改為可調(diào)控制增益kr和kv,則得到了一般型的ZEM/ZEV制導(dǎo)律公式[18]:

    (19)

    利用龐特里亞金極小值原理可以證明,E制導(dǎo)是常值引力場假設(shè)下的加速度模平方積分最優(yōu)制導(dǎo),即能量最優(yōu)制導(dǎo),使得如下性能指標(biāo)最小:

    (20)

    E制導(dǎo)性質(zhì)1.E制導(dǎo)的單個方向推力加速度為單調(diào)變化。

    (21)

    E制導(dǎo)性質(zhì)3.已知初始高度大于零,即rz0>0,在t∈[0,tf),滿足rz>0的條件為

    (22)

    2.2 阿波羅制導(dǎo)

    在阿波羅制導(dǎo)方法中,制導(dǎo)加速度三維矢量取為時間的二次多項式,包含3個待確定系數(shù)矢量,即指令推力加速度為

    a=c0+c1t+c2t2,t∈[0,tf]

    (23)

    將式(23)代入動力學(xué)方程(2)中,然后連續(xù)積分兩次,并代入位置、速度和加速度的初值和終值確定制導(dǎo)律系數(shù),結(jié)果如下:

    (24)

    (25)

    (26)

    將式(24)代入式(23)中后,得到E制導(dǎo)的初始時刻推力加速度計算公式:

    (27)

    由阿波羅制導(dǎo)的求解過程可知,求解常值制導(dǎo)參數(shù)時用到了終端位置、速度和加速度條件,因此與E制導(dǎo)相比,阿波羅制導(dǎo)除了可以滿足終端位置和速度約束外,還可以滿足終端加速度約束。

    3 多約束最優(yōu)反饋制導(dǎo)

    3.1 制導(dǎo)時長尋優(yōu)算法

    多約束最優(yōu)反饋制導(dǎo)的基礎(chǔ)是E制導(dǎo),也即基本型ZEM/ZEV制導(dǎo)。E制導(dǎo)給出了終端時刻固定的指令加速度計算公式,沒有給出終端時刻的確定方法。如果制導(dǎo)加速度允許范圍不受約束,則E制導(dǎo)的終端時刻可由E制導(dǎo)性質(zhì)2,即式(21)確定,保證了制導(dǎo)加速度模平方的積分最小。工程應(yīng)用中制導(dǎo)加速度允許范圍一定受限。一般來講,制導(dǎo)時長越短越節(jié)省推進(jìn)劑[6],但是推力變化范圍也會隨之變大,最終指令推力超過了發(fā)動機(jī)實際允許范圍,即指令推力飽和??紤]推力范圍約束后,實際輸出推力矢量滿足如下公式:

    (28)

    E制導(dǎo)有解析的軌跡預(yù)報公式,理論上只要終端指令推力不飽和,則規(guī)劃軌跡滿足終端狀態(tài)約束。因此,將終端指令推力不飽和作為制導(dǎo)算法基本原則。在此基礎(chǔ)上盡量縮短制導(dǎo)時長作為優(yōu)化推進(jìn)劑消耗的方向。最優(yōu)制導(dǎo)時長參數(shù)難以通過解析計算公式直接獲得,因此這里采用數(shù)值方法來尋找。優(yōu)化制導(dǎo)時長首先需要給出其一個合理的初值tf0。該初值應(yīng)盡量接近最優(yōu)值又要便于計算。這里提出將使得初始指令推力正好等于允許最大推力的制導(dǎo)時長作為優(yōu)化初值tf0。設(shè)rf=0和vf=0,由式(15)推導(dǎo)得到E制導(dǎo)的推力加速度模平方的初值為

    (29)

    由最大允許推力Fmax和初始質(zhì)量m0得到最大允許初始加速度為a0max=Fmax/m0。將最大允許初始加速度a0max代入式(29)得到制導(dǎo)時長tf的4次方程:

    (30)

    利用解析公式求解該方程,將其中最大實數(shù)根作為制導(dǎo)時長初值tf0。將制導(dǎo)時長初值tf0代入E制導(dǎo)公式,進(jìn)行首次數(shù)值積分預(yù)報,得到終端時刻標(biāo)稱推力指令。如果終端指令推力飽和,則需要增加制導(dǎo)時長。這里采用預(yù)測校正方法來計算增加量。記數(shù)值積分預(yù)報給出的指令推力在制導(dǎo)末端持續(xù)飽和的時長為tfs,校正系數(shù)記為ka,則給出制導(dǎo)時長第1個校正公式如下:

    tf,i+1=tf,i+katfs

    (31)

    式中:下標(biāo)i為迭代計算拍數(shù)。將校正后制導(dǎo)時長再次代入制導(dǎo)方程,重新進(jìn)行數(shù)值積分預(yù)報。當(dāng)終端指令推力不再飽和時認(rèn)為制導(dǎo)時長優(yōu)化完成。

    如果首次數(shù)值積分預(yù)報的終端指令推力未飽和,則說明制導(dǎo)時長還可以縮短。記制導(dǎo)時長縮短系數(shù)為ks,則制導(dǎo)時長縮短公式如下:

    tf,i+1=kstf,i

    (32)

    應(yīng)保證制導(dǎo)時長縮短足夠充分,使得終端指令推力一定飽和。

    接下來再次進(jìn)行數(shù)值積分預(yù)報。仍然采用預(yù)測校正方法來計算增加量。記校正系數(shù)為kb,則給出制導(dǎo)時長第2個校正公式如下:

    tf,i+1=tf,i+kbtfs

    (33)

    校正系數(shù)ka和kb取值偏小的話,則單次修正量過小,需要的迭代計算次數(shù)多。校正系數(shù)取值偏大的話,則單次修正量過大,計算得到制導(dǎo)時長比其最優(yōu)值偏大較多。

    定義推力俯仰角為推力方向和天向的夾角。為避免著陸器大幅度姿態(tài)機(jī)動帶來多方面不利影響需要約束推力俯仰角。通常不希望推力俯仰角超過90°,對應(yīng)的等效條件是垂向推力加速度非負(fù)。由垂向推力加速度非負(fù)得到的制導(dǎo)時長的不等式約束為:

    (34)

    顯然該式給出了一個制導(dǎo)時長的最小值約束。

    3.2 碰撞規(guī)避算法

    碰撞地面風(fēng)險規(guī)避要求飛行高度不低于地面。由E制導(dǎo)性質(zhì)3得到制導(dǎo)時長需滿足不等式(22),因此定義碰撞規(guī)避時長如下

    (35)

    式中:vzmin為足夠小的速度閾值,避免除零。

    為規(guī)避碰撞地平風(fēng)險,計算垂向制導(dǎo)加速度時,制導(dǎo)時長計算公式如下:

    (36)

    依據(jù)上式得到新的制導(dǎo)加速度,進(jìn)行數(shù)值積分預(yù)報。如果推力沒有飽和,那么用上式方式來制導(dǎo)可以實現(xiàn)臨界碰撞規(guī)避。但是實際推力很可能初始就飽和,且上式取值臨界,此時碰撞規(guī)避就很可能已經(jīng)來不及。這種情況下需要進(jìn)一步優(yōu)化推力加速度的分配。由于碰撞規(guī)避的優(yōu)先級高,因此推力飽和時,優(yōu)先滿足垂向的推力加速度需求。給出碰撞規(guī)避控制的推力優(yōu)化分配律如下:

    (37)

    (38)

    式中:m為著陸器當(dāng)前質(zhì)量;amax為當(dāng)前最大允許推力加速度;ax、ay和az為加速度優(yōu)化分配前指令加速度三軸分量,且az≥0;F為優(yōu)化分配后的指令推力矢量。

    為了進(jìn)一步提高著陸安全性,增加碰撞規(guī)避的裕度,可以再增加下滑角約束,即著陸器與目標(biāo)著陸點連線,與地平的夾角需大于閾值。在滿足地平碰撞規(guī)避的基礎(chǔ)上,下滑角約束通過坐標(biāo)系旋轉(zhuǎn)的方法即可方便的得到滿足[11]。記下滑角閾值為α,則由制導(dǎo)坐標(biāo)系到下滑角坐標(biāo)系的旋轉(zhuǎn)陣為

    (39)

    將初始位置矢量r0和速度矢量v0,終端目標(biāo)位置矢量rf和速度矢量vf,引力加速度矢量g,都左乘坐標(biāo)系旋轉(zhuǎn)陣Cgl,得到下滑角坐標(biāo)系對應(yīng)的分量。將經(jīng)坐標(biāo)系旋轉(zhuǎn)變換后的變量作為新的初始約束、終端約束和引力加速度矢量代入制導(dǎo)計算公式便可得到滿足下滑角約束的解。滿足下滑角約束后既提高了碰撞規(guī)避能力又增加了著陸點可見性。

    3.3 勻速轉(zhuǎn)動制導(dǎo)

    著陸終端一般希望著陸器為豎直姿態(tài)并且合加速度近似為零,即推力加速度應(yīng)為負(fù)的引力加速度。這就對制導(dǎo)的終端加速度提出了約束。但是如前所述,E制導(dǎo)的終端三維加速度不完全受控,即不能滿足該約束。為了滿足終端加速度約束,在E制導(dǎo)結(jié)束之后,引入勻速轉(zhuǎn)動制導(dǎo),迅速調(diào)整推力加速度到目標(biāo)值。

    勻速轉(zhuǎn)動制導(dǎo)的功能是依據(jù)最大角速度和最大推力變化率約束實現(xiàn)最短時間內(nèi)調(diào)整推力加速度矢量的方向和大小。勻速轉(zhuǎn)動制導(dǎo)的指令推力在由初始和終端加速度確定的平面內(nèi)勻速轉(zhuǎn)動。勻速轉(zhuǎn)動制導(dǎo)的輸入為初始制導(dǎo)加速度a0和終端目標(biāo)加速度af,輸出為指令推力F、制導(dǎo)時長tf、位置增量Δr和速度增量Δv。初始制導(dǎo)加速度a0和終端目標(biāo)加速度af的夾角為:

    (40)

    (41)

    定義勻速轉(zhuǎn)動制導(dǎo)的平面坐標(biāo)系為:原點為推力作用點,X軸為初始制導(dǎo)加速度方向,X軸向終端目標(biāo)加速度方向旋轉(zhuǎn)90°得到Y(jié)軸。

    勻速轉(zhuǎn)動制導(dǎo)的原理是首先計算平面坐標(biāo)系下推力加速度矢量,然后將其轉(zhuǎn)換為三維坐標(biāo)系下矢量。最終的勻速轉(zhuǎn)動制導(dǎo)的指令推力計算公式為:

    (42)

    F0=m0||a0||

    (43)

    Ff=mf||af||

    (44)

    θ=θft/tf

    (45)

    式中:M為將平面坐標(biāo)系矢量轉(zhuǎn)為三維坐標(biāo)系矢量的轉(zhuǎn)換陣;mf為著陸器在制導(dǎo)終端的質(zhì)量,可以用初始質(zhì)量m0近似代替。

    對合加速度做一次和二次積分分別得到位置增量Δr和速度增量Δv的解析表達(dá)式:

    (46)

    (47)

    圖1 組合制導(dǎo)參數(shù)預(yù)測校正規(guī)劃方法示意Fig.1 Diagram of predictor-corrector planning method for integrated guidance parameters

    最優(yōu)反饋制導(dǎo)和勻速轉(zhuǎn)動制導(dǎo)交接處的位置rm和速度vm的校正公式如下:

    rm=rf-Δr

    (48)

    vm=vf-Δv

    (49)

    4 仿真校驗

    以火星探測任務(wù)的動力下降段為例,通過數(shù)學(xué)仿真方法來驗證和分析新提出的制導(dǎo)方法。仿真用例選取時參考了文獻(xiàn)[5,7,9,11]中的共同用例。著陸器的初始質(zhì)量取為1 905 kg,發(fā)動機(jī)比沖取為2 205 Ns/kg,制導(dǎo)坐標(biāo)系下的初始位置為[2 000, 0, 1 500]Tm,初始速度的多組工況在后面分別給出;目標(biāo)終端位置為[0, 0, 0]Tm,目標(biāo)終端速度為[0, 0, 0]Tm/s,目標(biāo)終端加速度為[0, 0, 0]Tm/s2。發(fā)動機(jī)推力最大值為13.258 kN,最小值為4.971 kN?;鹦且铀俣热閇0, 0, -3.711 4]Tm/s2。校正系數(shù)ka和kb的取值分別為0.7和0.3。制導(dǎo)時長縮短系數(shù)ks取為0.5。

    4.1 有效性驗證

    設(shè)計制導(dǎo)有效性驗證仿真工況對3.1節(jié)中方法做驗證,設(shè)定了4組初始速度取值,詳見表1。

    表1 有效性驗證仿真工況的初始速度Table 1 Initial velocities of effectiveness verification simulation cases

    有效性驗證仿真工況的軌跡見圖2。從圖中可知,4條軌跡的初始位置相同,由于初始速度不同導(dǎo)致中途軌跡不同,不過最終都匯集到了坐標(biāo)系原點,即目標(biāo)終端位置。比較工況2到3的制導(dǎo)飛行軌跡,當(dāng)航向初始速度為負(fù)值時,制導(dǎo)飛行軌跡偏左上,航向初始速度為正值時,制導(dǎo)飛行軌跡偏右下。比較工況3和4的制導(dǎo)飛行軌跡,當(dāng)垂向和航向初始速度相同,有橫向初始速度時,制導(dǎo)飛行軌跡將偏右下。有效性驗證仿真工況4的橫向位置速度曲線見圖3。從中可知,橫向位置偏差先增大然后逐漸減小到零,橫向速度偏差先由正值變?yōu)樨?fù)值然后逐漸減小到零。有效性驗證仿真工況的推力比見圖4。從中可知,初始航向和橫向速度越大則制導(dǎo)時長越長,初始推力飽和時長也越長。工況1、3和4的推力由飽和段、減小段和增加段組成,即推力大致按照大小大的趨勢變化,近似于最優(yōu)控制理論給出的變推力制導(dǎo)的最大最小最大推力剖面[19]。工況2雖然沒有推力飽和段,但是推力符合大、小、大的變化趨勢。有效性驗證仿真工況的推力俯仰角見圖5。從圖中可發(fā)現(xiàn)推力俯仰角的2個變化規(guī)律。第1個規(guī)律是推力俯仰角先大后小然后再增大,原因是著陸器為了實現(xiàn)水平方向的位移,需要在水平方向上先加速然后再減速。第2個規(guī)律是初始航向和橫向速度越大則初始推力俯仰角越大,原因是制導(dǎo)需要在水平方向分配更多的推力來抵消初始的反向水平速度。

    圖2 有效性驗證仿真工況的軌跡Fig.2 Trajectories of effectiveness verification simulation cases

    圖3 有效性驗證仿真工況4的橫向位置速度Fig.3 Cross position and velocity of effectiveness verification simulation case 4

    圖4 有效性驗證仿真工況的推力比Fig.4 Thrust throttles of effectiveness verification simulation cases

    圖5 有效性驗證仿真工況的推力俯仰角Fig.5 Thrust pitch angles of effectiveness verification simulation cases

    制導(dǎo)有效性驗證仿真工況的制導(dǎo)時長和推進(jìn)劑消耗仿真結(jié)果見表2。從表中可知,初始航向和橫向速度越大則制導(dǎo)時長和推進(jìn)劑消耗越長和越多。工況2由于初始航向速度與目標(biāo)水平移動方向相同,因而制導(dǎo)時長最短并且推進(jìn)劑消耗最少。表中給出了預(yù)測校正法和遍歷法兩種方法的推進(jìn)劑消耗來對比。遍歷法的搜索步長取為1 s??梢钥闯?預(yù)測校正法的推進(jìn)劑消耗十分接近遍歷法,最大僅相差1.1%。表中還給出了兩種預(yù)測校正的實際迭代次數(shù)。工況1實際執(zhí)行式(31)所示的預(yù)測校正0次,執(zhí)行式(33)所示的預(yù)測校正2次。工況2實際執(zhí)行式(31)所示的預(yù)測校正1次,執(zhí)行式(33)所示的預(yù)測校正0次。所有4個工況僅需執(zhí)行1到2次預(yù)測校正便可獲得滿意的結(jié)果,說明所提出的預(yù)測校正方法不僅有效而且高效。

    表2 有效性驗證仿真工況的制導(dǎo)時長和推進(jìn)劑消耗Table 2 Time-to-go and fuel consumption of effectiveness verification simulation cases

    4.2 碰撞規(guī)避驗證

    文獻(xiàn)[5,9,11]中為了測試驗證制導(dǎo)算法碰撞規(guī)避性能,設(shè)計了惡劣的碰撞規(guī)避仿真工況,初始速度取為[100, 0, -75]Tm/s。在此初始速度條件下,分別按照無碰撞規(guī)避、有地平規(guī)避、有8°下滑角約束和有9°下滑角約束進(jìn)行制導(dǎo)方法的測試仿真,對比驗證3.2節(jié)中的碰撞規(guī)避方法。碰撞規(guī)避仿真工況的軌跡見圖6。從中可知,無碰撞規(guī)避的制導(dǎo)飛行軌跡有部分位于了地平線下,該結(jié)果意味著著陸器將與地面碰撞??紤]地平規(guī)避的制導(dǎo)飛行軌跡在著陸后期緊緊貼著地平,實現(xiàn)了臨界的碰撞規(guī)避,該結(jié)果具有理論分析意義,但是不能用于實際飛行,因為其安全裕度不夠高,實際地面必然是具有一定幅度的起伏。有8°和9°下滑角約束的著陸軌跡可以明顯看到著陸后期保持了一定的飛行高度,與地平線形成了固定夾角,顯然該飛行軌跡的碰撞規(guī)避裕度更高,對一定幅度的地形突起也可以實現(xiàn)碰撞規(guī)避,因而更具有工程應(yīng)用價值。工況2到4的高度是單調(diào)變化的,另外著陸后期的軌跡線為直線,原因是著陸后期在下滑角坐標(biāo)系下垂向速度已被控制到零。文獻(xiàn)[5]利用航路點方法作碰撞規(guī)避的后期航跡也近似直線,然而航路點前后制導(dǎo)推力不連續(xù)。文獻(xiàn)[9,11]中碰撞規(guī)避段的高度是先探底然后上升之后再下降,也就是高度非單調(diào)下降,并且文獻(xiàn)[11]中也有制導(dǎo)推力躍變點。碰撞規(guī)避仿真工況的垂向速度見圖7,從中可知,工況1的垂向速度先為負(fù)值后為正值,工況2和3的垂向速度全程為非正值,工況4的垂向速度先為負(fù)值中間為正值后為負(fù)值。工況4的垂向速度出現(xiàn)正值的原因是在碰撞規(guī)避控制段的推力始終飽和。碰撞規(guī)避仿真工況的推力比曲線見圖8,從中可知,全部工況的推力由飽和段、減小段和增加段組成,即推力大致按照大小大的趨勢變化。比較工況2到4可知,碰撞規(guī)避的約束下滑角越大,制導(dǎo)時長越長。另外,考慮碰撞規(guī)避后制導(dǎo)初期的推力飽和段時長增大,并且碰撞規(guī)避裕度越大,制導(dǎo)初期推力飽和段時長增大的越多。碰撞規(guī)避仿真工況的推力俯仰角見圖9,從中可知,碰撞規(guī)避要求越嚴(yán)格則推力的初始俯仰角越小,如果初始推力飽和則初始俯仰角等于負(fù)的下滑角。該現(xiàn)象的原因是,為了實現(xiàn)碰撞規(guī)避就需要減小下降速度,因此制導(dǎo)算法在垂向分配的推力更大,因此推力俯仰角就會相對變小。工況4的推力俯仰角在30 s時發(fā)生3.1°的躍變,原因是碰撞規(guī)避控制結(jié)束時刻已到,但是推力還處于飽和狀態(tài),因此導(dǎo)致制導(dǎo)加速度變化不連續(xù)。工況4說明如果希望制導(dǎo)推力連續(xù),則在碰撞規(guī)避控制結(jié)束時制導(dǎo)推力應(yīng)處于非飽和狀態(tài)。工況3在保證指令推力連續(xù)的前提下滿足了8°下滑角約束,大于文獻(xiàn)[5,11]中相同條件下實現(xiàn)的4°。

    圖6 碰撞規(guī)避仿真工況的軌跡Fig.6 Trajectories of collision avoidance simulation cases

    圖7 碰撞規(guī)避仿真工況的垂向速度Fig.7 Vertical velocities of collision avoidance simulation cases

    圖8 碰撞規(guī)避仿真工況的推力比Fig.8 Thrust throttles of collision avoidance simulation cases

    圖9 碰撞規(guī)避仿真工況的推力俯仰角Fig.9 Thrust pitch angles of collision avoidance simulation cases

    碰撞規(guī)避仿真工況的制導(dǎo)時長和推進(jìn)劑消耗見表3。從中可知,有地平規(guī)避與無碰撞規(guī)避相比制導(dǎo)時長減少3.2 s推進(jìn)劑消耗增加6.5 kg,有8°下滑角約束與無碰撞規(guī)避相比制導(dǎo)時長增加0.7 s推進(jìn)劑消耗增加18.3 kg,有9°下滑角約束與無碰撞規(guī)避相比制導(dǎo)時長增加13.9 s推進(jìn)劑消耗增加75.6 kg。本節(jié)仿真結(jié)果說明新制導(dǎo)方法的碰撞規(guī)避策略有效,并且滿足碰撞規(guī)避要求后付出的推進(jìn)劑代價相對小。

    表3 碰撞規(guī)避仿真工況的制導(dǎo)時長和推進(jìn)劑消耗Table 3 Time-to-go and fuel consumptions of collision avoidance simulation cases

    4.3 有終端加速度約束

    設(shè)計有終端加速度約束仿真工況來驗證3.3節(jié)中方法,初始速度取為[0, 0, -75]Tm/s,目標(biāo)終端加速度為[0, 0, 0]Tm/s2。姿態(tài)機(jī)動最大角速度和最大推力變化率取值見表4。預(yù)測校正迭代計算的終止條件為終端位置偏差小于1 m,速度偏差小于0.1 m/s。

    表4 有終端加速度約束仿真工況的最大角速度和最大推力變化率Table 4 Maximum angular rate and maximum thrust change rate of terminal acceleration constrained simulation cases

    有終端加速度約束仿真工況的軌跡見圖10,從中可知,著陸前中期3條軌跡基本重合,著陸后期3條軌跡略有差異,終端時刻3條軌跡都滿足了終端位置約束。有終端加速度約束仿真工況的推力比見圖11,以工況1為例,前40.2 s為最優(yōu)反饋制導(dǎo),之后為勻速轉(zhuǎn)動制導(dǎo),推力比從0.95降至0.47。有終端加速度約束仿真工況的俯仰角見圖12,以工況1為例,在40.2 s處可以清楚的看到俯仰角發(fā)生轉(zhuǎn)折,這是因為制導(dǎo)律切換為了勻速轉(zhuǎn)動制導(dǎo),俯仰角迅速從54.7°轉(zhuǎn)動到0°。推力比和俯仰角的仿真結(jié)果表明制導(dǎo)終端狀態(tài)滿足了終端加速度約束。

    圖10 有終端加速度約束仿真工況的軌跡Fig.10 Trajectories of terminal acceleration constrained simulation cases

    圖11 有終端加速度約束仿真工況的推力比Fig.11 Thrust throttles of terminal acceleration constrained simulation cases

    圖12 有終端加速度約束仿真工況的俯仰角Fig.12 Thrust pitch angles of terminal acceleration constrained simulation cases

    有終端加速度約束仿真工況的數(shù)值結(jié)果見表5。從表中可知,雖然勻速轉(zhuǎn)動制導(dǎo)時長隨著最大角速度和最大推力變化率增大而減小,但是總制導(dǎo)時長并沒有出現(xiàn)等幅度變化,而是變化很小。勻速轉(zhuǎn)動制導(dǎo)時長越短則推進(jìn)劑消耗越少。工況1的推進(jìn)劑消耗最多,與4.1節(jié)中無終端加速度約束工況相比僅增加了6.3 kg,相對增量為2.7%。工況1到3的預(yù)測校正迭代次數(shù)依次減少,說明勻速轉(zhuǎn)動制導(dǎo)允許的最大角速度和最大推力變化率越大,則預(yù)測校正需要的迭代次數(shù)越少。以上結(jié)果表明最優(yōu)反饋制導(dǎo)和勻速轉(zhuǎn)動制導(dǎo)的組合使用方法有效,使得著陸終端既滿足了位置速度約束,又滿足了加速度約束。

    表5 有終端加速度約束仿真工況的制導(dǎo)時長和推進(jìn)劑消耗Table 5 Time-to-go and fuel consumption of terminal acceleration constrained simulation cases

    5 結(jié) 論

    針對行星定點著陸動力下降段制導(dǎo)問題,從工程實際出發(fā),考慮推力范圍受限、連續(xù)變推力、碰撞規(guī)避和終端加速度共4方面約束,對傳統(tǒng)的最優(yōu)反饋制導(dǎo)方法進(jìn)行改進(jìn)和完善。首先,提出推力范圍受限且連續(xù)變推力約束下的最優(yōu)反饋制導(dǎo)的制導(dǎo)時長高效尋優(yōu)方法,經(jīng)1到2次預(yù)測校正即可獲得推進(jìn)劑消耗量接近遍歷法的結(jié)果,推力剖面為大小大形式。其次,提出基于碰撞規(guī)避時長和推力優(yōu)化分配律的碰撞規(guī)避方法,結(jié)合坐標(biāo)系旋轉(zhuǎn)可滿足下滑角約束,在連續(xù)變推力前提下顯著提高碰撞規(guī)避能力和著陸點可見性,著陸全程高度單調(diào)下降,碰撞規(guī)避控制結(jié)束后的縱向平面軌跡為直線。最后,提出勻速轉(zhuǎn)動制導(dǎo)和基于預(yù)測校正的組合制導(dǎo)參數(shù)規(guī)劃方法,連續(xù)銜接最優(yōu)反饋制導(dǎo)和勻速轉(zhuǎn)動制導(dǎo),充分發(fā)揮推力調(diào)節(jié)能力,使得終端狀態(tài)擴(kuò)充滿足加速度約束。新制導(dǎo)方法改進(jìn)原理清晰,多約束適應(yīng)性強(qiáng),雖然使用了數(shù)值積分預(yù)報,但是總體還屬于解析制導(dǎo)方法,計算量顯著低于凸優(yōu)化等數(shù)值優(yōu)化方法和人工智能方法,適合工程應(yīng)用于行星定點著陸任務(wù)中。

    猜你喜歡
    推進(jìn)劑制導(dǎo)校正
    劉光第《南旋記》校正
    國學(xué)(2020年1期)2020-06-29 15:15:30
    一類具有校正隔離率隨機(jī)SIQS模型的絕滅性與分布
    機(jī)內(nèi)校正
    基于MPSC和CPN制導(dǎo)方法的協(xié)同制導(dǎo)律
    基于在線軌跡迭代的自適應(yīng)再入制導(dǎo)
    帶有攻擊角約束的無抖振滑模制導(dǎo)律設(shè)計
    KNSB推進(jìn)劑最佳配比研究
    復(fù)合制導(dǎo)方式確保精確入軌
    太空探索(2014年1期)2014-07-10 13:41:49
    含LLM-105無煙CMDB推進(jìn)劑的燃燒性能
    無鋁低燃速NEPE推進(jìn)劑的燃燒性能
    天堂8中文在线网| 在线观看免费高清a一片| 91老司机精品| 欧美大码av| 熟女少妇亚洲综合色aaa.| 欧美黄色淫秽网站| 日韩伦理黄色片| 日韩制服骚丝袜av| 久久性视频一级片| 免费高清在线观看视频在线观看| 久久免费观看电影| 国产成人精品无人区| 19禁男女啪啪无遮挡网站| 免费高清在线观看日韩| 久久人妻熟女aⅴ| 亚洲av日韩在线播放| 丝瓜视频免费看黄片| 国产在线一区二区三区精| 精品一品国产午夜福利视频| 精品国产一区二区三区久久久樱花| 青草久久国产| 色综合欧美亚洲国产小说| 精品少妇内射三级| 蜜桃国产av成人99| 蜜桃在线观看..| 精品免费久久久久久久清纯 | 亚洲午夜精品一区,二区,三区| 极品少妇高潮喷水抽搐| 啦啦啦 在线观看视频| 18在线观看网站| 国产男女内射视频| 亚洲欧美一区二区三区国产| 丝瓜视频免费看黄片| 国产欧美亚洲国产| 国产国语露脸激情在线看| 成人影院久久| 丰满饥渴人妻一区二区三| 久久狼人影院| 中国美女看黄片| 熟女少妇亚洲综合色aaa.| av不卡在线播放| 欧美+亚洲+日韩+国产| 亚洲久久久国产精品| 中国美女看黄片| 亚洲欧洲日产国产| 19禁男女啪啪无遮挡网站| 亚洲精品久久成人aⅴ小说| 你懂的网址亚洲精品在线观看| 亚洲国产av影院在线观看| 久久综合国产亚洲精品| 国产精品免费视频内射| 免费看不卡的av| 交换朋友夫妻互换小说| 国产成人欧美在线观看 | 性色av乱码一区二区三区2| 极品人妻少妇av视频| 丁香六月欧美| 韩国精品一区二区三区| 中文字幕另类日韩欧美亚洲嫩草| 菩萨蛮人人尽说江南好唐韦庄| av国产精品久久久久影院| 亚洲成人免费电影在线观看 | 日本av手机在线免费观看| 亚洲av欧美aⅴ国产| 99久久人妻综合| kizo精华| av有码第一页| 欧美成狂野欧美在线观看| 麻豆国产av国片精品| 亚洲精品国产一区二区精华液| 这个男人来自地球电影免费观看| 国产精品 欧美亚洲| 久久精品国产综合久久久| 久久久久久人人人人人| 精品人妻熟女毛片av久久网站| 青草久久国产| 色视频在线一区二区三区| 九色亚洲精品在线播放| 国产一卡二卡三卡精品| 成人黄色视频免费在线看| 久久午夜综合久久蜜桃| 超碰97精品在线观看| 熟女av电影| 国产黄频视频在线观看| 亚洲色图综合在线观看| 亚洲国产欧美日韩在线播放| 久久精品人人爽人人爽视色| 各种免费的搞黄视频| 一区二区三区乱码不卡18| 午夜免费成人在线视频| 色94色欧美一区二区| 久久精品久久久久久噜噜老黄| 脱女人内裤的视频| 国产欧美日韩一区二区三区在线| 色94色欧美一区二区| 精品一区二区三区四区五区乱码 | 国产免费又黄又爽又色| 永久免费av网站大全| 婷婷丁香在线五月| 亚洲精品第二区| 亚洲人成网站在线观看播放| 欧美日韩视频精品一区| 好男人视频免费观看在线| av线在线观看网站| 国产成人精品无人区| 免费人妻精品一区二区三区视频| 亚洲激情五月婷婷啪啪| 香蕉国产在线看| 一二三四社区在线视频社区8| 欧美97在线视频| 两个人免费观看高清视频| 考比视频在线观看| 久久久国产精品麻豆| 伊人亚洲综合成人网| 18禁国产床啪视频网站| 在线天堂中文资源库| 亚洲欧美一区二区三区黑人| 亚洲色图 男人天堂 中文字幕| 色婷婷av一区二区三区视频| 热re99久久国产66热| 亚洲专区中文字幕在线| 你懂的网址亚洲精品在线观看| 成人黄色视频免费在线看| 高潮久久久久久久久久久不卡| 蜜桃国产av成人99| 一级毛片黄色毛片免费观看视频| 18禁裸乳无遮挡动漫免费视频| 免费日韩欧美在线观看| 国产精品国产av在线观看| 国产精品国产三级专区第一集| 成人18禁高潮啪啪吃奶动态图| 久久人人爽人人片av| 亚洲成人国产一区在线观看 | 亚洲国产欧美网| 伊人久久大香线蕉亚洲五| 日韩,欧美,国产一区二区三区| 日本色播在线视频| 王馨瑶露胸无遮挡在线观看| 国产成人系列免费观看| 男的添女的下面高潮视频| svipshipincom国产片| 一级黄色大片毛片| 欧美日本中文国产一区发布| 国产97色在线日韩免费| 日韩一区二区三区影片| 五月开心婷婷网| 两人在一起打扑克的视频| 99久久综合免费| 在线观看人妻少妇| 久久免费观看电影| 久久国产亚洲av麻豆专区| av线在线观看网站| 亚洲精品自拍成人| 国产成人av激情在线播放| 黄色 视频免费看| 老汉色av国产亚洲站长工具| 在线av久久热| 国产精品熟女久久久久浪| 在线观看人妻少妇| 男女免费视频国产| 成人18禁高潮啪啪吃奶动态图| 一区二区三区激情视频| 国产欧美日韩精品亚洲av| 日本一区二区免费在线视频| 午夜激情av网站| 亚洲第一青青草原| 操美女的视频在线观看| 国产黄色视频一区二区在线观看| 久久鲁丝午夜福利片| 精品国产乱码久久久久久小说| 9191精品国产免费久久| 男女无遮挡免费网站观看| 视频在线观看一区二区三区| 大片免费播放器 马上看| 亚洲国产毛片av蜜桃av| 亚洲五月婷婷丁香| 中文字幕av电影在线播放| 女人被躁到高潮嗷嗷叫费观| 成年美女黄网站色视频大全免费| 久久久久久久久免费视频了| av在线app专区| 最近最新中文字幕大全免费视频 | 每晚都被弄得嗷嗷叫到高潮| 欧美xxⅹ黑人| 欧美日韩视频高清一区二区三区二| 飞空精品影院首页| av网站免费在线观看视频| 伦理电影免费视频| 18禁观看日本| 国产在线免费精品| 婷婷丁香在线五月| 91字幕亚洲| 天堂中文最新版在线下载| 伦理电影免费视频| 丝袜在线中文字幕| 国产福利在线免费观看视频| 亚洲精品av麻豆狂野| videos熟女内射| 最近最新中文字幕大全免费视频 | 国产亚洲欧美在线一区二区| av一本久久久久| 亚洲精品中文字幕在线视频| 男女之事视频高清在线观看 | 欧美精品高潮呻吟av久久| 亚洲精品自拍成人| 国产在视频线精品| 欧美黑人欧美精品刺激| 大码成人一级视频| 午夜视频精品福利| 亚洲一区二区三区欧美精品| netflix在线观看网站| 久久亚洲精品不卡| 国产成人一区二区三区免费视频网站 | 亚洲欧美中文字幕日韩二区| 久久女婷五月综合色啪小说| 搡老乐熟女国产| 国产高清视频在线播放一区 | 最新在线观看一区二区三区 | 久久人妻熟女aⅴ| 制服人妻中文乱码| 精品亚洲乱码少妇综合久久| 欧美精品一区二区免费开放| 丝袜脚勾引网站| 日韩,欧美,国产一区二区三区| 国产一区亚洲一区在线观看| 国产精品久久久人人做人人爽| 久久久久网色| 在线 av 中文字幕| 啦啦啦啦在线视频资源| 80岁老熟妇乱子伦牲交| 夫妻性生交免费视频一级片| 免费看十八禁软件| 亚洲精品国产一区二区精华液| 夜夜骑夜夜射夜夜干| 国产精品久久久久久精品古装| 久热这里只有精品99| 日本a在线网址| 少妇人妻久久综合中文| 国产99久久九九免费精品| 日韩一区二区三区影片| 1024香蕉在线观看| 久久天躁狠狠躁夜夜2o2o | 天堂俺去俺来也www色官网| 国产主播在线观看一区二区 | 亚洲国产最新在线播放| 极品少妇高潮喷水抽搐| 最近中文字幕2019免费版| 777米奇影视久久| 国产日韩欧美视频二区| 少妇裸体淫交视频免费看高清 | 制服人妻中文乱码| 国产成人系列免费观看| 亚洲人成电影免费在线| 又紧又爽又黄一区二区| 大陆偷拍与自拍| 啦啦啦在线观看免费高清www| 亚洲熟女毛片儿| 亚洲成色77777| 男女高潮啪啪啪动态图| 久久久国产欧美日韩av| 亚洲天堂av无毛| 亚洲欧美清纯卡通| 捣出白浆h1v1| 多毛熟女@视频| 亚洲图色成人| 后天国语完整版免费观看| 国产精品久久久久成人av| 国产精品麻豆人妻色哟哟久久| 免费高清在线观看视频在线观看| 国产在视频线精品| 五月开心婷婷网| 午夜激情久久久久久久| 成年人午夜在线观看视频| 80岁老熟妇乱子伦牲交| av又黄又爽大尺度在线免费看| 韩国高清视频一区二区三区| 国产成人精品久久二区二区91| 啦啦啦中文免费视频观看日本| 国产亚洲av片在线观看秒播厂| 国产黄色免费在线视频| 国产片特级美女逼逼视频| 亚洲国产精品999| 首页视频小说图片口味搜索 | 制服人妻中文乱码| 国产xxxxx性猛交| 中文字幕人妻熟女乱码| 亚洲熟女毛片儿| 男女之事视频高清在线观看 | 人妻 亚洲 视频| 99香蕉大伊视频| 黄色视频在线播放观看不卡| 国产熟女欧美一区二区| 中文字幕av电影在线播放| 两个人免费观看高清视频| 国产成人av教育| 日韩中文字幕欧美一区二区 | 亚洲精品成人av观看孕妇| 女性生殖器流出的白浆| 色视频在线一区二区三区| 久久精品熟女亚洲av麻豆精品| 国产精品一国产av| 侵犯人妻中文字幕一二三四区| 午夜福利在线免费观看网站| 亚洲欧洲日产国产| 国产在视频线精品| 无遮挡黄片免费观看| 日韩,欧美,国产一区二区三区| 免费看十八禁软件| av福利片在线| 国产成人精品久久二区二区免费| 韩国精品一区二区三区| 国产真人三级小视频在线观看| 九色亚洲精品在线播放| 丰满少妇做爰视频| 精品福利观看| 亚洲专区国产一区二区| 亚洲五月色婷婷综合| 国产精品国产三级国产专区5o| 国产精品99久久99久久久不卡| a级毛片黄视频| 国产亚洲av高清不卡| 精品福利观看| 自线自在国产av| 久久精品熟女亚洲av麻豆精品| 久久天堂一区二区三区四区| 黑丝袜美女国产一区| 18在线观看网站| 两性夫妻黄色片| 午夜激情久久久久久久| 我的亚洲天堂| 欧美性长视频在线观看| 精品国产超薄肉色丝袜足j| 91字幕亚洲| 久久天堂一区二区三区四区| 你懂的网址亚洲精品在线观看| 国产黄频视频在线观看| 日本一区二区免费在线视频| 午夜激情久久久久久久| 中文字幕高清在线视频| 精品少妇内射三级| 免费日韩欧美在线观看| 一二三四在线观看免费中文在| 麻豆乱淫一区二区| 亚洲成人免费av在线播放| 亚洲视频免费观看视频| 肉色欧美久久久久久久蜜桃| 免费在线观看黄色视频的| 男女国产视频网站| 在线观看www视频免费| 岛国毛片在线播放| 两人在一起打扑克的视频| 老司机影院成人| av线在线观看网站| 91麻豆精品激情在线观看国产 | 欧美激情高清一区二区三区| 精品久久久久久电影网| 又大又爽又粗| 国产麻豆69| 男男h啪啪无遮挡| 欧美成人精品欧美一级黄| 久久精品aⅴ一区二区三区四区| 欧美日韩亚洲综合一区二区三区_| 操出白浆在线播放| 久久 成人 亚洲| 国产视频一区二区在线看| 日韩 欧美 亚洲 中文字幕| 国产av精品麻豆| 少妇人妻 视频| 一级毛片女人18水好多 | 人体艺术视频欧美日本| 69精品国产乱码久久久| 久久久久精品国产欧美久久久 | 亚洲精品乱久久久久久| 人人澡人人妻人| 亚洲熟女毛片儿| 国产精品二区激情视频| 一区二区三区精品91| 成人影院久久| 丝袜美足系列| www.自偷自拍.com| 欧美激情极品国产一区二区三区| 亚洲伊人久久精品综合| 国产成人a∨麻豆精品| 亚洲专区中文字幕在线| 99re6热这里在线精品视频| 91国产中文字幕| 搡老乐熟女国产| 一级黄色大片毛片| av欧美777| 精品国产一区二区三区久久久樱花| 国产黄频视频在线观看| 国产高清国产精品国产三级| 男人爽女人下面视频在线观看| 午夜福利视频在线观看免费| 色网站视频免费| 精品国产乱码久久久久久小说| 精品一区二区三区av网在线观看 | 伊人久久大香线蕉亚洲五| 久久久久精品人妻al黑| 亚洲色图 男人天堂 中文字幕| 国产成人av激情在线播放| 一级黄片播放器| 久久久欧美国产精品| a级毛片在线看网站| 国产高清国产精品国产三级| 亚洲欧美精品综合一区二区三区| 国产一区有黄有色的免费视频| 久久中文字幕一级| 18禁裸乳无遮挡动漫免费视频| 精品国产国语对白av| 成人国产一区最新在线观看 | 国产一区二区 视频在线| 欧美大码av| 国产精品一区二区在线不卡| 秋霞在线观看毛片| 成年动漫av网址| 欧美日韩综合久久久久久| 菩萨蛮人人尽说江南好唐韦庄| 国产成人精品久久二区二区免费| av天堂在线播放| 久久女婷五月综合色啪小说| 免费在线观看黄色视频的| 亚洲av片天天在线观看| videos熟女内射| 国产伦人伦偷精品视频| 国产97色在线日韩免费| 午夜福利影视在线免费观看| 国产亚洲av片在线观看秒播厂| 天天添夜夜摸| 久久久精品国产亚洲av高清涩受| 热re99久久国产66热| 亚洲图色成人| 人人妻人人澡人人爽人人夜夜| 一区二区av电影网| 婷婷色综合www| 天天躁狠狠躁夜夜躁狠狠躁| 日韩一本色道免费dvd| 久久影院123| 97精品久久久久久久久久精品| 国产在线视频一区二区| 两人在一起打扑克的视频| 亚洲欧美一区二区三区国产| kizo精华| 日本午夜av视频| 国产老妇伦熟女老妇高清| 熟女av电影| 亚洲伊人久久精品综合| 97精品久久久久久久久久精品| 五月开心婷婷网| 一边亲一边摸免费视频| 蜜桃国产av成人99| 色94色欧美一区二区| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲五月色婷婷综合| 看十八女毛片水多多多| 亚洲精品乱久久久久久| 下体分泌物呈黄色| 老司机影院毛片| 久久ye,这里只有精品| www.自偷自拍.com| 国产精品 国内视频| 丁香六月欧美| 巨乳人妻的诱惑在线观看| 大香蕉久久成人网| 一本久久精品| 午夜免费男女啪啪视频观看| 91九色精品人成在线观看| 成人午夜精彩视频在线观看| 少妇人妻 视频| 欧美成人午夜精品| 99精国产麻豆久久婷婷| 在线观看国产h片| av线在线观看网站| 亚洲激情五月婷婷啪啪| 中文字幕高清在线视频| 国产午夜精品一二区理论片| 又大又黄又爽视频免费| 国产精品99久久99久久久不卡| 视频区欧美日本亚洲| 欧美少妇被猛烈插入视频| 少妇粗大呻吟视频| 亚洲国产精品一区二区三区在线| 黄色怎么调成土黄色| 国产精品麻豆人妻色哟哟久久| 国产精品国产三级专区第一集| 国产在视频线精品| 大香蕉久久网| 男女边吃奶边做爰视频| 男人舔女人的私密视频| 1024视频免费在线观看| 一本综合久久免费| 亚洲成人国产一区在线观看 | 又紧又爽又黄一区二区| 成在线人永久免费视频| 精品第一国产精品| 久久精品熟女亚洲av麻豆精品| 日本wwww免费看| 国产日韩一区二区三区精品不卡| 日本五十路高清| 国产亚洲av片在线观看秒播厂| 免费人妻精品一区二区三区视频| 国产精品熟女久久久久浪| 欧美激情极品国产一区二区三区| 久久久久精品国产欧美久久久 | 国产成人精品久久二区二区免费| 午夜福利一区二区在线看| 狂野欧美激情性xxxx| 久久久久久久国产电影| 欧美在线一区亚洲| 久久 成人 亚洲| 免费看不卡的av| 七月丁香在线播放| 999久久久国产精品视频| 日韩一本色道免费dvd| 热99国产精品久久久久久7| 亚洲成色77777| 我的亚洲天堂| 亚洲视频免费观看视频| 国产激情久久老熟女| tube8黄色片| 婷婷色综合大香蕉| 一边摸一边抽搐一进一出视频| 国产一区二区激情短视频 | 欧美国产精品一级二级三级| 好男人电影高清在线观看| 久久性视频一级片| a 毛片基地| 最近最新中文字幕大全免费视频 | 欧美日本中文国产一区发布| 国产高清国产精品国产三级| 国产成人系列免费观看| 亚洲一码二码三码区别大吗| 搡老乐熟女国产| 亚洲自偷自拍图片 自拍| 嫩草影视91久久| 丝袜喷水一区| 亚洲av电影在线进入| 日韩一卡2卡3卡4卡2021年| 久久国产亚洲av麻豆专区| 新久久久久国产一级毛片| 男人操女人黄网站| 又黄又粗又硬又大视频| 啦啦啦 在线观看视频| 制服人妻中文乱码| 亚洲第一av免费看| 99久久精品国产亚洲精品| 999精品在线视频| 在线观看免费日韩欧美大片| 国产视频一区二区在线看| 麻豆av在线久日| 热99国产精品久久久久久7| 欧美激情高清一区二区三区| 天堂8中文在线网| 美女午夜性视频免费| 制服人妻中文乱码| av电影中文网址| 性高湖久久久久久久久免费观看| 国产一区二区 视频在线| 考比视频在线观看| 2018国产大陆天天弄谢| 欧美av亚洲av综合av国产av| 99热网站在线观看| 韩国精品一区二区三区| 视频区图区小说| 久久人人97超碰香蕉20202| 欧美日韩福利视频一区二区| 蜜桃国产av成人99| 亚洲精品日韩在线中文字幕| 超碰成人久久| 久久人妻熟女aⅴ| 色婷婷久久久亚洲欧美| 亚洲国产日韩一区二区| 男女国产视频网站| 国产一卡二卡三卡精品| 国产视频首页在线观看| 国产91精品成人一区二区三区 | 国产成人系列免费观看| 中文字幕色久视频| 国产亚洲av片在线观看秒播厂| 国产成人91sexporn| 国产亚洲av片在线观看秒播厂| 少妇人妻 视频| 亚洲专区中文字幕在线| 18禁黄网站禁片午夜丰满| 嫩草影视91久久| 老司机午夜十八禁免费视频| 女人高潮潮喷娇喘18禁视频| 国产一区二区在线观看av| 国产在线一区二区三区精| 女警被强在线播放| 欧美日韩亚洲高清精品| 少妇人妻 视频| 丁香六月天网| xxx大片免费视频| 欧美黑人精品巨大| 亚洲欧洲精品一区二区精品久久久| 午夜免费观看性视频| 伦理电影免费视频| 欧美黄色片欧美黄色片| 国产一区二区 视频在线| 国产欧美日韩精品亚洲av| 国产免费视频播放在线视频| 看免费成人av毛片| 国产精品一区二区在线不卡| 欧美黑人欧美精品刺激| 中文字幕av电影在线播放| 国产欧美日韩精品亚洲av| 50天的宝宝边吃奶边哭怎么回事| 深夜精品福利| 国产成人av激情在线播放| 欧美激情极品国产一区二区三区| 亚洲精品av麻豆狂野| 国产男女超爽视频在线观看| 欧美 亚洲 国产 日韩一| 天天操日日干夜夜撸|