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

RLV末端能量管理段三維魯棒制導方法

2021-11-14 11:29:38穆凌霞張友民
宇航學報 2021年9期
關(guān)鍵詞:指令模型

穆凌霞,余 翔,張友民

(1. 西安理工大學陜西省復雜系統(tǒng)控制與智能信息處理重點實驗室,西安 710048;2. 北京航空航天大學自動化科學與電氣學院,北京 100191;3. 康考迪亞大學機械、工業(yè)與航空工程系,蒙特利爾H3G 1M8)

0 引 言

可重復使用運載器(Reusable launch vehicle, RLV)作為一種新型天地往返空間飛行器,可極大地降低運載成本、提高空間飛行的靈活性[1-2]。RLV返回過程中經(jīng)歷末端能量管理段(Terminal area energy management, TAEM)、高超音速再入段(Hypersonic reentry, HRE)、進場著陸段(Approach and landing, A&L)。其中TAEM段是RLV實現(xiàn)自主水平著陸的關(guān)鍵環(huán)節(jié),在該段RLV無動力滑翔,需實現(xiàn)大范圍內(nèi)的橫側(cè)向機動任務(wù),面臨著在三維空間的高速返回、多約束、多擾動的挑戰(zhàn),如何實現(xiàn)三維高精度魯棒制導成為RLV自主返回的難點之一。

目前RLV末端能量管理段的制導方法主要有兩類。第一類是二維解耦軌跡制導技術(shù)。典型的是美國航天飛機TAEM制導方案[3],將地面軌跡分為S轉(zhuǎn)彎段、捕獲段、航向校準段(HAC)和預著陸段,縱向軌跡則根據(jù)地面航程進行實時調(diào)整。文獻[4]針對亞軌道飛行器Hopper研究了改進的航天飛機TAEM制導方案,通過減少飛行段數(shù)、放寬軌跡幾何約束,提高軌跡的靈活性。文獻[5]設(shè)計了改進的航天飛機TAEM策略,縱向剖面為馬赫數(shù)隨高度變化的三次函數(shù)。文獻[6]通過改變HAC半徑降低了軌跡形狀的約束,迭代三個參數(shù)求解TAEM軌線,增強了軌跡魯棒性。文獻[7]通過分析RLV末端約束條件,提出能量管的概念,若軌跡位于能量管中間位置,則具有較強的魯棒性。文獻[8-11]在航天飛機TAEM制導策略的基礎(chǔ)上,對算法的實時性的提高進行了嘗試。文獻[12]在HAC前增加直線預測捕獲段,并預先規(guī)劃TAEM軌跡,降低實時性的要求。文獻[13]研究了減速板故障下的縱向剖面修正的方法。

第二類是考慮縱橫向耦合的三維空間上的軌跡制導方法。有代表性的為文獻[14]提出的內(nèi)核提取的次優(yōu)節(jié)點應用方法(Sub-optimal nodal application of the kernel extraction, SNAKE),設(shè)計動壓剖面,利用高度域動力學方程迭代計算三維軌跡。文獻[15]對SNAKE迭代求解邏輯和地面軌跡機動方式進行了改進,加快了軌跡推演的速度。文獻[16]將雙回路迭代算法改進為遞推算法,提高了求解速度。文獻[17]提出沿最優(yōu)路徑飛行(Optimum-path-to-go, OPTG)的制導方法,離線構(gòu)造軌跡數(shù)據(jù)庫,并對離線軌跡進行編碼,根據(jù)當前飛行狀態(tài)在線選擇參考軌跡,同時考慮擾動對軌跡進行修正,實現(xiàn)了基于數(shù)據(jù)庫的軌跡重規(guī)劃方法。文獻[18]針對亞軌道飛行器Hopper和X-33,將末端制導問題描述為一定約束條件下的優(yōu)化問題,利用加速梯度投影算法求解該優(yōu)化問題,實現(xiàn)三維參考剖面的更新。文獻[19]利用微分平坦理論對TAEM軌跡優(yōu)化問題進行降維,通過B樣條對問題進行離散化,并采用SNOPT工具箱求解。

總的來說,上述兩類TAEM制導方法本質(zhì)上為參考軌跡制導,考慮到末端能量管理段任務(wù)的特殊性(三維空間大幅度橫側(cè)向運動),實現(xiàn)過程中,在算法復雜度、制導精度和在線自適應能力上有所折衷。比如第一種航天飛機TAEM制導方案,采用了離線和在線相結(jié)合的思路,地面軌跡形狀為離線設(shè)定,縱向軌跡在線調(diào)整,是一種縱橫向解耦的策略,易于工程實現(xiàn),但是制導精度和自適應能力需進一步提高。第二種SNAKE制導方案,考慮了三維特性,但是軌跡迭代時間較長,實時性和魯棒性有待提高。由于RLV末端返回時需完成三維機動飛行,其制導模型為多輸入多輸出非線性模型,制導問題為強非線性多約束優(yōu)化問題,在線TAEM軌跡制導方法具有一定難度,傳統(tǒng)在線數(shù)值預測-校正制導不適用于多參數(shù)變量和多約束的制導問題,而由最優(yōu)控制問題轉(zhuǎn)化得到的非線性規(guī)劃問題在解的存在性及計算時間方面難以滿足。

本文針對RLV在末端能量管理段所面臨的復雜的制導任務(wù)需求,提出了一種基于滑模控制理論的魯棒制導及指令更新策略,從理論上證明并保證了軌跡生成的可靠性,結(jié)合周期性校正機制,抵抗返回過程的氣動、大氣等干擾。主要貢獻在于:

1)設(shè)計三維滑模面,當?shù)竭_其中兩個滑模面時可滿足TAEM末端狀態(tài)的約束(包括航跡角、航向角、縱向位置和橫向位置),而另一個滑模面則為實際動壓剖面與理想動壓剖面的誤差項,當?shù)竭_該滑模面時,意味著RLV的動壓跟蹤上理想動壓剖面,從而間接滿足TAEM整個過程中的動壓約束;

2)基于高度域制導模型,設(shè)計二階滑模制導律,并證明滑模面在指定高度可達;利用滑模制導律,得到滿足多約束的TAEM軌跡及制導指令,僅與當前系統(tǒng)狀態(tài)相關(guān),無需迭代求解,減少軌跡生成的時間;

3)考慮到RLV飛行中受氣動、大氣和質(zhì)量偏差等干擾,長時間單一制導指令易導致累積誤差無法滿足高精度TAEM末端約束需求,設(shè)計周期性制導指令更新機制,可根據(jù)RLV當前實際狀態(tài)進行軌跡及指令的修正,從而抵抗大氣和氣動等干擾的影響。

1 問題描述

1.1 面向制導的動力學方程

RLV在末端能量管理段的三維滑行動力學方程可表示為:

(1)

式中:ft=[ft1,ft2,ft3,ft4,ft5,ft6]T,具體形式為:

(2)

升力系數(shù)CL和阻力系數(shù)CD為迎角α的多項式函數(shù):

CL=CL0(Ma)+CL1(Ma)α+CL2(Ma)α2

(3)

CDCD0(Ma)+CD1(Ma)α+CD2(Ma)α2+

CD3(Ma)α3+ksb(Ma)δ

(4)

式中:CLi(i=0,1,2)和CDi(i=0,1,2,3)分別為馬赫數(shù)在[0.5,2.5]范圍內(nèi)的多項式系數(shù),隨馬赫數(shù)的變化而變化,ksb為減速板δ對阻力系數(shù)的影響因子。

(5)

式中:fH=[fH 1,fH 2,fH 3,fH 4,fH 5]T,具體形式為:

(6)

1.2 末端能量管理段制導任務(wù)及約束條件分析

RLV在末端能量管理段飛行示意圖如圖1所示。從末端進入點(Terminal entry point, TEP)飛行至進場著陸段的起點ALI處,高度跨越25908 m~3048 m,馬赫數(shù)范圍在Ma2.5~Ma0.5,跨越超音速、跨音速和亞音速飛行。

圖1 RLV末端能量管理段示意圖Fig.1 Schematic diagram of RLV trajectory during TAEM phase

表1 TAEM段飛行約束條件Table 1 Flight constraints during TAEM phase

2 制導策略及流程

設(shè)計TAEM三維魯棒軌跡制導策略如圖2所示,基本思路及實現(xiàn)流程如下:

圖2 TAEM三維魯棒制導策略結(jié)構(gòu)圖Fig.2 Flowchart of the proposed TAEM guidance strategy

1)基于標稱高度域模型設(shè)計滑模制導律,尋找制導指令,使得RLV能夠從當前狀態(tài)(并不需要是起點位置)滑行到末端狀態(tài),并且滿足如表1所示的過程動壓約束和終端約束;將制導律作為輸入代入高度域制導微分方程。利用龍格-庫塔法求解從當前狀態(tài)到末端狀態(tài)的常微分方程初值問題(Initial value problem, IVP),得到離散制導指令序列為u-h。

2)將該制導指令作為時域系統(tǒng)的輸入(時域系統(tǒng)上可增加隨機干擾項),因其隨飛行高度而變化,因此實際的輸入實時地根據(jù)線性插值計算得到。在一個制導周期內(nèi)利用龍格-庫塔法求解由帶干擾的實際時域制導模型及其當前狀態(tài)所決定的常微分方程初值問題。

3)考慮到實際時域系統(tǒng)模型可能受到氣動、大氣和質(zhì)量偏差等干擾,采用理想制導指令不能完全跟蹤期望軌線。設(shè)計了周期性的制導校正機制,一個制導周期后更新RLV系統(tǒng)狀態(tài),并判斷是否到達TAEM終端完成任務(wù)。若到達TAEM終端,則結(jié)束,若未到達,則重新根據(jù)當前狀態(tài)利用滑模制導律,生成滿足約束條件的軌線及制導指令,實現(xiàn)制導更新,抵抗大氣、氣動、質(zhì)量偏差等擾動對TAEM返回過程的影響,完成高精度TAEM制導任務(wù)。

值得注意的是,其中滑模制導律是使得系統(tǒng)(即標稱高度域制導模型)實現(xiàn)期望目標的閉環(huán)狀態(tài)反饋律,因此求解由滑模制導律、高度域模型和初始狀態(tài)構(gòu)建的IVP,得到的狀態(tài)軌線能夠滿足:當RLV到達指定終端高度時,終端約束和飛行過程中動壓滿足要求,這將在第3.1節(jié)得到證明。

此外,為了滿足TAEM飛行過程中的動壓約束條件及終端能量的要求,本文預先設(shè)定了縱向動壓-高度剖面,并將實際剖面與預定剖面之間的誤差作為滑模面之一,利用滑模制導律實現(xiàn)對預定剖面的跟蹤,通過該方式,可間接滿足TAEM過程中的動壓約束及終端能量約束。傳統(tǒng)軌跡生成方法根據(jù)預定的動壓剖面,雙回路或三回路迭代求解其他軌跡參數(shù)[11,14],計算耗時較長,無法適應實時校正的需求。本文通過設(shè)計三維滑模面及二階滑模制導律(詳見第3.1節(jié)),將制導律作為輸入代入高度域制導微分方程,求解IVP問題,從而實現(xiàn)動壓約束和其他終端狀態(tài)約束的同時滿足。無需迭代求解,可避免優(yōu)化無解或時間過長的問題,提高軌跡生成的速度。

3 三維魯棒滑模制導算法實現(xiàn)

3.1 考慮多重約束的滑模制導指令生成

本節(jié)給出基于高度域制導模型的三維滑模面及二階滑模制導律的設(shè)計過程,旨在保證所生成的狀態(tài)軌線滿足預先設(shè)定的動壓-高度剖面約束,同時滿足TAEM末端位置、航跡角以及航向角的約束,并生成制導指令。

3.1.1滑模面設(shè)計

設(shè)計三維滑模面向量sa[sa1,sa2,sa3]T為:

sa=fa(x(H),H)

(7)

其中:fa=[fa1,fa2,fa3]T設(shè)計為下式的形式:

(8)

滑模面sa對虛擬高度H求一階導和二階導,可得到:

(9)

Guv+F

(10)

式中:矩陣G和F的元素為:

(11)

(12)

G13=G21=G31=F2=F3=0

(13)

式中:uv稱為虛擬輸入,定義為uv[C′D,γ′,′]T。3.1.2滑模制導律綜合設(shè)計

因此,帶約束的TAEM制導問題轉(zhuǎn)變?yōu)槿绾伪WCsa=0和s′a=0。由式(10)可見,對sa求導兩次,顯含控制輸入uv,利用反步法為該二階系統(tǒng)設(shè)計制導律。設(shè)計二階滑模面為:

(14)

式(14)對虛擬高度求一階導后有:

A+F+Guv

(15)

其中:A=(kasa+kas′a(Hf-H))/(Hf-H)2,G和F同式(10)。

定理1.考慮系統(tǒng)(14),在制導律:

uv=-G-1[F+A+ksgn(sb)]

(16)

(17)

的作用下,當H趨于Hr(Hr≤ηHf),sb趨于平衡點0。

(18)

記Hr為到達滑模面sb=0時的高度,此時有Vb(Hr)=0。兩邊積分得到:

(19)

(20)

(21)

式中:Δ為滑模邊界層,取常值。

定理 2.考慮動態(tài)系統(tǒng)

s′a=-kasa/(Hf-H),ka>1

(22)

在區(qū)間[Hr,Hf]內(nèi),當H趨近Hf時,sa和s′a趨近于 0。

(23)

由于Hf總是大于H,且ka>0,有V′a<0。因此當H→Hf有sa→0。

對式(22)的每個分量變換得:

dsai/sai=-kadH/(Hf-H)

(24)

記sa(H=Hr)sar=[sar1,sar2,sar3]T,對上式積分,有:

(25)

因此sai=sari(Hf-H)ka/(Hf-Hr)ka,則s′ai=-kasari(Hf-H)ka-1/(Hf-Hr)ka,寫成矩陣形式有

(26)

因為,ka>1,Hr

根據(jù)虛擬輸入uv,利用方程(5)中的關(guān)系易計算得到實際輸入:

u=fu(xH,uv)

(27)

具體如下:

(28)

(29)

(30)

(31)

fH=[fHi]T|i=1,…,7

(32)

其中:fHi|i=1,…,5由式(6)決定,fH 6=uv1,fH 7=-1/(Vsinγ)。

值得注意的是,上述設(shè)計的高度域滑模制導律將用于生成TAEM軌跡及制導指令,該方法的優(yōu)勢之一是從理論上證明了采用該制導律能夠保證所生成的軌跡滿足約束要求。此外,在生成軌跡時無需求解優(yōu)化問題,避免了優(yōu)化求解時可能因無法找到最優(yōu)解而導致軌跡規(guī)劃失敗、無法獲得有效制導指令的問題,因此本文方法可保證軌跡生成的可靠性及制導指令解的存在。

3.2 實時指令校正機制

理論上采用該制導指令作為時域制導模型的輸入,將能得到與標稱高度域IVP一致的狀態(tài)軌線。但是在實際飛行中,系統(tǒng)受多種擾動影響,如氣動模型的偏差、大氣環(huán)境建模的偏差和RLV質(zhì)量的變化。若采用理想制導指令作為整個返回過程的制導輸入,因擾動引起的累計誤差將越來越大,最終無法到達指定的末端返回點。因此,增加周期性的制導校正環(huán)節(jié),根據(jù)系統(tǒng)實際狀態(tài)進行實時制導修正,抵抗大氣和氣動等干擾的影響。所設(shè)計的RLV末端能量管理段三維魯棒軌跡制導算法步驟如下:

1) 初始化xH(H0)xt(t0)。

2) 標稱軌線生成:考慮標稱高度域模型x′H=fH(xH,u,H),其中u=fu(xH)由式(16)和(27)決定,龍格-庫塔法求解常微分方程初值問題(初始狀態(tài)為RLV當前狀態(tài)),步長ΔH=3m,得到狀態(tài)xH在高度域的軌線,滿足TAEM約束。

3) 制導指令生成:根據(jù)步驟2,計算得到剩余飛行段的離散制導指令相對于高度的序列u=fuh(h)。

4) 制導周期內(nèi)實際軌線生成:針對考慮干擾的實際時域系統(tǒng)模型(1),根據(jù)步驟3中離散制導指令線性插值實時計算任意高度處的制導指令,輸入給時域模型,龍格-庫塔法求解IVP,步長為Δt=0.01 s,持續(xù)一個制導周期(10 s),得到RLV當前的實際狀態(tài)。

5) 狀態(tài)更新:記錄當前狀態(tài),計算下一周期的制導指令,注意自變量H的終端值會發(fā)生變化。時域模型的狀態(tài)更新為xt(t0),高度域模型的狀態(tài)更新為xH(H0),Hf更新為末端初始高度減去實際飛過的高度。

6) 判斷是否到達終端高度,若是結(jié)束,若否返回步驟2。

4 仿真結(jié)果及分析

4.1 模型及其仿真設(shè)置

仿真驗證中采用的RLV模型為X-34驗證機[20],TAEM返回的末端約束如表1所示,初始狀態(tài)設(shè)置為:

(33)

滑模制導律式(14)、(16)、(21)中的參數(shù)設(shè)置為η=0.7,ka=3, Δ=0.1。仿真任務(wù)分為兩部分:1)標稱情形下制導策略有效性及實時性驗證; 2)干擾情形下魯棒性驗證。

4.2 標稱情形下仿真結(jié)果分析

設(shè)定TAEM返回初始狀態(tài)如式(33)所示,根據(jù)標稱高度域制導模型(31)及制導律(式(16)和(27)),在整個TAEM高度范圍內(nèi),計算得到滿足約束條件的狀態(tài)軌線,并計算出相應的理想制導指令。本文利用龍格-庫塔法求解IVP生成軌跡用時約0.55 s,相比于文獻[14]中的基于迭代的三維軌跡規(guī)劃算法平均用時1 min左右、文獻[19]的基于非線性規(guī)劃方法的三維軌跡生成用時約1.39 s,本文方法提高了軌跡生成速度。

所得到的狀態(tài)軌線結(jié)果如圖3和圖4中的實線所示,圖3(a)、(b)分別為TAEM地面軌跡形狀和高度隨時間的變化曲線,其中時間信息是通過式(31)中時間狀態(tài)得到,完成25908 m到3048 m的TAEM返回飛行需要約414.3 s。圖4(a)~(c)分別為RLV的動壓、航跡角和航向角隨高度變化的曲線,可以看到整個TAEM返回過程動壓不超過最大動壓19152.1 Pa的約束,而當高度到達指定的TAEM末端高度3048 m時,航跡角和航向角均達到指定約束條件,即末端航跡角-10°、末端航向角0°。相應的理想制導指令如圖5(a)~(c)所示,迎角、滾轉(zhuǎn)角和減速板,其指令均在合理的范圍內(nèi),并且較為平滑,易于后期內(nèi)環(huán)控制回路的跟蹤。由此可見,所生成的標稱狀態(tài)軌線能夠滿足如表1所示的約束條件,并能夠滿足TAEM返回制導任務(wù)需求。

將理想制導指令輸入給時域制導模型,注意到理想制導指令是隨高度變化的離散指令,仿真中采用線性插值實時計算任意高度處所需的指令。得到的狀態(tài)軌線如圖3~圖4中虛線所示,與理想狀態(tài)結(jié)果是一致的。此過程中未考慮時域模型的干擾,因此未加周期性校正環(huán)節(jié)。

圖3 標稱情形下TAEM位置和高度軌線Fig.3 TAEM trajectories in nominal case

圖4 標稱情形下TAEM動壓、航跡角和航向角軌線Fig.4 TAEM states histories with respect to altitude in nominal case

圖5 標稱情形下理想制導指令Fig.5 The predicted guidance commands in nominal case

圖6給出的是滑模面隨虛擬高度H變化的曲線圖。由圖6(a)可見,sb的滑動運動發(fā)生在虛擬高度6400 m處,并保持到終端虛擬高度22860 m,由圖6(b)和6(c)可見,在sb到達零之后,sa和s′a也將趨近于零,正因為如此,系統(tǒng)約束得以滿足。

圖6 標稱情形下滑模面曲線Fig.6 The sliding surface histories with respected to H in nominal case

4.3 干擾情形下仿真結(jié)果分析

考慮干擾因素對所提方法進行驗證,初始條件同式(33),返回過程中的干擾因素見表2,升力系數(shù)、阻力系數(shù)、大氣密度和RLV質(zhì)量出現(xiàn)不超過20%的服從均勻分布的隨機偏差。

表2 TAEM段干擾因素Table 2 Disturbance during TAEM phase

當不加校正環(huán)節(jié)時,在整個TAEM返回過程中采用單一的理想制導指令,仿真結(jié)果如圖7所示,實線為滑模制導計算得到的理想軌線,方點和圓點分別為TAEM起點和期望落點,虛線為實際地軌線,隨著返回過程的持續(xù)進行,干擾影響不斷累加,RLV實際軌線漸漸偏離理想軌跡,實際落點為棱形,偏離期望落點較遠,偏差Δxe=1863.6 m, Δye=-1756.8 m,偏差量級在103m,不能有效完成任務(wù)。

圖7 未加校正的TAEM地軌線Fig.7 The TAEM ground-track trajectory without correction

當加入周期性校正機制,每10 s重新根據(jù)當前狀態(tài)計算理想制導指令,仿真結(jié)果如圖8和9所示。圖8為TAEM地面軌線,其中方點表示RLV在當前制導周期的起點位置,根據(jù)當前位置及其他狀態(tài)信息,利用滑模制導律求解IVP得到的后續(xù)理想軌線為實線所示,通過及時的校正更新可以保證在TAEM末端時滿足約束要求。圖9(a)~(c)分別為加入校正機制之后的動壓、航跡角和航向角軌線,其中圓點為時域制導模型積分10 s時所得到的狀態(tài),代表每飛行10 s后,RLV當前的狀態(tài)。實線表示從當前狀態(tài)(圓點所在處)在剩余高度范圍內(nèi)的理想狀態(tài)軌線。可見,在出現(xiàn)擾動的情況下,RLV利用理想制導指令無法完全跟蹤理想軌線從而導致末端狀態(tài)無法滿足約束,但是引入校正機制后,根據(jù)當前狀態(tài)重新更新飛行狀態(tài)軌線,得到新的制導指令,從而實現(xiàn)魯棒制導,有效抵抗干擾對制導精度的影響。最終的位置偏差Δxe和Δye在3 m以內(nèi),航跡角偏差約0.4°,航向角偏差為0.001°。可見該制導策略精度較高,魯棒性較好。

圖8 干擾情形下校正后的TAEM地軌線Fig.8 The corrected TAEM ground-track trajectory in the presence of disturbances

圖9 干擾情形下校正后的TAEM狀態(tài)軌線Fig.9 The corrected TAEM state histories in the presence of disturbances

5 結(jié) 論

RLV末端能量管理段返回制導時,模型復雜、約束條件眾多,傳統(tǒng)縱橫向解耦制導策略不易滿足高精度要求,而傳統(tǒng)三維制導迭代計算時間不易滿足實時性和魯棒性的要求,針對該問題,本文提出了基于高度域滑模制導律的快速軌跡生成方法及周期性校正策略,提高制導的魯棒性和精度。通過設(shè)計三維滑模面及二階滑模制導指令,生成滿足多重約束要求的狀態(tài)軌線,避免了傳統(tǒng)三維制導方法中迭代求解制導指令的策略,保證了制導指令解的存在以及求解效率,且所生成的制導指令較為平滑易于內(nèi)環(huán)的跟蹤控制。所提制導方法在較大的氣動、大氣和質(zhì)量干擾的情形下,依然能實現(xiàn)高精度制導的需求,具有較好的魯棒性。

猜你喜歡
指令模型
一半模型
聽我指令:大催眠術(shù)
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
ARINC661顯控指令快速驗證方法
LED照明產(chǎn)品歐盟ErP指令要求解讀
電子測試(2018年18期)2018-11-14 02:30:34
殺毒軟件中指令虛擬機的脆弱性分析
電信科學(2016年10期)2016-11-23 05:11:56
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
一種基于滑窗的余度指令判別算法
主站蜘蛛池模板: 国产精品网址你懂的| 日韩小视频网站hq| 无码aⅴ精品一区二区三区| 国产亚洲精品自在线| 欧美一区精品| 欧美日韩高清在线| 亚洲AV无码乱码在线观看裸奔| 超碰aⅴ人人做人人爽欧美| 67194亚洲无码| 亚洲第一黄色网| 国产成人91精品| AV片亚洲国产男人的天堂| 91在线国内在线播放老师| 国产成人av一区二区三区| 国产精品.com| 91久久精品国产| 99热这里只有精品2| 欧美激情综合| 毛片网站免费在线观看| 美女被操91视频| 97精品久久久大香线焦| 国产在线观看91精品亚瑟| 91精品国产一区| 免费又黄又爽又猛大片午夜| 国产三级精品三级在线观看| 色婷婷在线播放| 亚洲va在线∨a天堂va欧美va| 大学生久久香蕉国产线观看| 欧美综合在线观看| 亚洲成a人在线观看| 免费国产无遮挡又黄又爽| 欧美日韩在线观看一区二区三区| 欧美日韩中文国产| 国产91在线|中文| 全裸无码专区| 国产一区二区三区日韩精品| 国产成人8x视频一区二区| 欧美精品伊人久久| 久久精品电影| 青青草一区二区免费精品| 啊嗯不日本网站| 精品国产Ⅴ无码大片在线观看81| 高清免费毛片| 国产欧美在线视频免费| 九色视频一区| 波多野结衣中文字幕一区二区| 中文字幕伦视频| 色偷偷一区| 另类欧美日韩| 日本道综合一本久久久88| 一本大道在线一本久道| 综合色区亚洲熟妇在线| a级毛片在线免费| 波多野结衣第一页| aaa国产一级毛片| 久久综合国产乱子免费| 亚洲人成成无码网WWW| 亚洲欧美一级一级a| 日韩欧美一区在线观看| 国产成人综合网| 秋霞午夜国产精品成人片| 香蕉国产精品视频| 99久久人妻精品免费二区| 综合天天色| 99爱在线| 国产免费福利网站| 欧美国产综合色视频| 亚洲精品你懂的| 免费观看国产小粉嫩喷水| 青青青视频蜜桃一区二区| 91娇喘视频| 中国一级特黄大片在线观看| 国产精品久久久免费视频| 国产小视频a在线观看| 9966国产精品视频| 精品剧情v国产在线观看| 九九九国产| 亚洲国产天堂在线观看| 成人精品午夜福利在线播放| 婷婷六月综合网| 中文成人在线视频| 久久精品人人做人人综合试看|