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

基于高斯偽譜法的高超聲速飛行器再入制導(dǎo)研究

2018-01-05 01:10:41仲維昆屈泉酉原勁鵬武文斌
計算機(jī)測量與控制 2017年12期
關(guān)鍵詞:指令系統(tǒng)

仲維昆,屈泉酉,原勁鵬,武文斌

(1.北京空間飛行器總體設(shè)計部,北京 100094; 2.中國兵器工業(yè)第203研究所,西安 710065)

基于高斯偽譜法的高超聲速飛行器再入制導(dǎo)研究

仲維昆1,屈泉酉1,原勁鵬1,武文斌2

(1.北京空間飛行器總體設(shè)計部,北京 100094; 2.中國兵器工業(yè)第203研究所,西安 710065)

針對高超聲速飛行器再入標(biāo)準(zhǔn)軌跡制導(dǎo)方法中存在的制導(dǎo)準(zhǔn)備周期長、彈上需存儲標(biāo)準(zhǔn)軌跡參數(shù)、制導(dǎo)魯棒性較差等缺點,提出一種基于高斯偽譜法與滾動時域控制技術(shù)相結(jié)合的高超聲速飛行器再入預(yù)測-校正制導(dǎo)方案;其中,在線高斯偽譜法采用縱/側(cè)向結(jié)合、全程一體化的制導(dǎo)算法思路,實現(xiàn)了對高超聲速飛行器再入彈道的全程預(yù)測制導(dǎo);同時結(jié)合滾動時域控制技術(shù)從工程上實現(xiàn)了高超聲速飛行器再入制導(dǎo)中對開環(huán)制導(dǎo)信息的閉環(huán)應(yīng)用,完成了飛行器預(yù)測-校正制導(dǎo)方案;通過對高超聲速飛行器再入制導(dǎo)過程進(jìn)行仿真分析,結(jié)果表明應(yīng)用文章設(shè)計的基于高斯偽譜法與滾動時域控制技術(shù)相結(jié)合的高超聲速飛行器再入預(yù)測-校正制導(dǎo)方案,飛行器再入過程中具有良好的制導(dǎo)性能。

高超聲速飛行器; 預(yù)測-校正; 高斯偽譜法; 滾動時域控制

0 引言

高超聲速飛行器再入制導(dǎo)方法根據(jù)是否存儲制導(dǎo)信息可以分為標(biāo)準(zhǔn)軌跡制導(dǎo)方法和預(yù)測-校正制導(dǎo)方法兩大類。預(yù)測-校正制導(dǎo)方法由于具有無需存儲標(biāo)稱軌跡參數(shù)、制導(dǎo)魯棒性強(qiáng)以及制導(dǎo)準(zhǔn)備周期短等優(yōu)點,逐漸成為未來高超聲速飛行器再入制導(dǎo)發(fā)展的重點方向。文章將結(jié)合高斯偽譜法與滾動時域控制方法,研究高超聲速飛行器再入一體化制導(dǎo)方案。這里的一體化是指兩方面,其一是再入飛行器縱/側(cè)向制導(dǎo)的一體化設(shè)計;其二是再入飛行器全程彈道無分段的一體化設(shè)計。飛行器再入過程中基于在線預(yù)測制導(dǎo)模型,應(yīng)用高斯偽譜法快速計算得到的縱/側(cè)向通道制導(dǎo)指令,結(jié)合滾動時域控制技術(shù)不斷實時更新、輸出制導(dǎo)指令,完成有約束情況下的高超聲速飛行器再入制導(dǎo)[1]。

1 高斯偽譜法

高斯偽譜法是求解最優(yōu)控制方法中直接法的一種,又可稱為正交配點法。它首先將原系統(tǒng)狀態(tài)變量以及控制變量在一系列的Gauss點上進(jìn)行離散化。其次應(yīng)用Largrange插值多項式對已離散的系統(tǒng)狀態(tài)變量及控制變量進(jìn)行擬合逼近。以此將原系統(tǒng)微分方程、約束方程等轉(zhuǎn)化為代數(shù)方程約束。而原連續(xù)系統(tǒng)狀態(tài)變量的導(dǎo)數(shù)也通過對離散化后的系統(tǒng)狀態(tài)變量全局插值多項式的求導(dǎo)代替[2]。同時對于原連續(xù)系統(tǒng)最優(yōu)控制問題中積分項,則通過對原系統(tǒng)相應(yīng)狀態(tài)量、控制量的插值多項式進(jìn)行高斯積分求得到。

經(jīng)過上述一系列變換,高斯偽譜法已將原系統(tǒng)最優(yōu)控制問題轉(zhuǎn)化為一種非線性規(guī)劃問題。針對非線性規(guī)劃問題,文章利用成熟的非線性規(guī)劃求解器對原問題進(jìn)行求解。具體表示基于高斯偽譜法求解最優(yōu)控制問題的相關(guān)步驟可得:

1)時域變換

首先對原連續(xù)系統(tǒng)時間進(jìn)行高斯變換,將時間區(qū)間由t∈[t0,tf]轉(zhuǎn)換到τ∈[-1,1],因此引入時域變換:

(1)

則原連續(xù)系統(tǒng)經(jīng)時間歸一化后,性能指標(biāo)函數(shù)即為:

(2)

則原連續(xù)系統(tǒng)經(jīng)時間歸一化后,系統(tǒng)動力學(xué)微分方程即為:

(3)

則原系統(tǒng)經(jīng)時間歸一化后,系統(tǒng)邊界條件和約束即為:

(4)

2)狀態(tài)變量和控制變量的離散化

高斯偽譜法基于K階的Legendre插值多項式Pk(τ)去逼近經(jīng)時間離散化后的系統(tǒng)狀態(tài)變量和控制變量,其中:

(5)

因為在進(jìn)行系統(tǒng)狀態(tài)變量、控制變量離散時,Legendre-Gauss點僅分布在區(qū)間(-1,1)上。所以為了近似原系統(tǒng)狀態(tài)變量、控制變量量的初值,需在離散時間區(qū)間內(nèi)增加一個點τ0=-1。由此構(gòu)成K+1階的Lagrange插值多項式Lagi(τ)(i=0,···,K)將作為基函數(shù)在時間區(qū)間內(nèi)逼近原系統(tǒng)狀態(tài)變量和控制變量。綜上所述,可得原系統(tǒng)狀態(tài)量離散形式為:

(6)

原系統(tǒng)控制變量離散形式為:

(7)

其中:Lagrange插值基函數(shù)為:

(8)

3)微分方程變換

因為原系統(tǒng)狀態(tài)變量已經(jīng)由Lagrange插值多項式在時間區(qū)間內(nèi)離散化,所以原系統(tǒng)狀態(tài)方程亦可以表示為相應(yīng)代數(shù)方程形式,即:

(9)

其中:D∈RK×(K+1)為系統(tǒng)微分矩陣,可以通過離線計算得到。則可得原系統(tǒng)微分方程的代數(shù)表達(dá)形式為:

(k=1,···,K)

(10)

4)性能指標(biāo)變換

原系統(tǒng)狀態(tài)變量、控制變量均已有Lagrange插值多項式在時間區(qū)間內(nèi)離散化。則原連續(xù)系統(tǒng)最優(yōu)控制性能指標(biāo)亦可由相應(yīng)變量離散形式表示。而性能指標(biāo)函數(shù)中存在的積分可由插值多項式的Gauss積分近似。

則原系統(tǒng)性能指標(biāo)函數(shù)可表示為:

(11)

5)終端狀態(tài)變換

(12)

由式(12)可知,系統(tǒng)終端狀態(tài)可由初始狀態(tài)的Gauss積分近似求解。綜上可知系統(tǒng)終端狀態(tài)表達(dá)式為:

(13)

通過上述系統(tǒng)變量離散化、微分方程代數(shù)化等步驟,已把原系統(tǒng)最優(yōu)控制問題轉(zhuǎn)換為非線性規(guī)劃問題。而非線性規(guī)劃問題的求解已經(jīng)具備足夠多的求解方法,例如可以采用序列二次規(guī)劃(SQP)的方法進(jìn)行求解。由于非線性規(guī)劃問題不是文章研究的重點,所以不做詳細(xì)論述。僅應(yīng)用基于Matlab軟件平臺的GPOPS-2軟件包對最優(yōu)控制問題、非線性規(guī)劃問題進(jìn)行求解[3-4]。

2 滾動時域控制

考慮高超聲速飛行器再入制導(dǎo)動力學(xué)微分方程可以通過非線性微分方程加以描述:

(14)

基于高斯偽譜法在線求解高超聲速飛行器再入預(yù)測制導(dǎo)指令后,由于求解得到的制導(dǎo)指令為開環(huán)形式,無法直接作用于飛行器制導(dǎo)控制回路。所以需要對已生成的開環(huán)制導(dǎo)指令進(jìn)行處理,而滾動時域控制恰好解決開環(huán)制導(dǎo)指令如何進(jìn)行閉環(huán)制導(dǎo)的問題。

因為滾動時域控制方法需要在線實時獲得偽譜法計算所需初始狀態(tài)并應(yīng)用其計算所得制導(dǎo)指令,這就對在線高斯偽譜法的軌跡優(yōu)化時間以及滾動時域制導(dǎo)更新周期做出了限制。高超聲速飛行器再入制導(dǎo)過程中的為了加快在線軌跡優(yōu)化算法收斂速率,通常選取第一次軌跡優(yōu)化結(jié)果作為偽譜法軌跡規(guī)劃中狀態(tài)變量、控制變量的參考值,并在再入優(yōu)化過程中不斷迭代更新參考值,以達(dá)到減小在線軌跡規(guī)劃時間的目的。

而對于滾動時域制導(dǎo)更新周期,Ross通過合理的假設(shè)和嚴(yán)謹(jǐn)?shù)臄?shù)學(xué)推倒,建立了滾動時域控制中制導(dǎo)信息更新最大時間間隔和Lipschitz常數(shù)之間的關(guān)系。但是這只是從數(shù)學(xué)理論上推倒得到的滾動時域控制時間間隔條件。由于在實際應(yīng)用中相關(guān)公式所需參數(shù)不易獲得,所以滾動時域控制時間間隔將采用仿真試湊與分析的方法加以確定[6]。

3 預(yù)測-校正再入制導(dǎo)方法

結(jié)合文章前兩小節(jié)內(nèi)容,可知應(yīng)用基于高斯偽譜法與滾動時域控制相結(jié)合的高超聲速飛行器再入預(yù)測-校正制導(dǎo)方法具體步驟如下:

1)初始參數(shù)裝訂:建立高超聲速飛行器再入預(yù)測-校正算法數(shù)學(xué)模型,確定高超聲速飛行器再入點、目標(biāo)點狀態(tài)信息。建立高超聲速飛行器再入過程約束模型及其計算方法,確定高超聲速飛行器再入過程約束、終端狀態(tài)約束參數(shù)值。其中在線軌跡規(guī)劃算法滿足高超聲速飛行器再入過程中的過載、熱流、動壓以及平衡滑翔條件等約束。同時需離線設(shè)計滿足在線計算時間、制導(dǎo)魯棒性要求的滾動時域控制指令更新周期;

2)飛行器在線制導(dǎo)指令計算:高超聲速飛行器再入后需實時測量在線軌跡優(yōu)化算法所需系統(tǒng)狀態(tài)量和控制量,應(yīng)用高斯偽譜進(jìn)行在線軌跡規(guī)劃,生成制導(dǎo)指令。其中,在線軌跡規(guī)劃算法規(guī)劃的彈道為由當(dāng)前高超聲速飛行器所處狀態(tài)至目標(biāo)點的彈道參數(shù),制導(dǎo)指令為縱/側(cè)向一體化計算結(jié)果;

3)飛行器制導(dǎo)指令的應(yīng)用:基于高斯偽譜法在線計算所得高超聲速飛行器再入縱/側(cè)向制導(dǎo)指令為由當(dāng)前狀態(tài)至終端狀態(tài)的制導(dǎo)參數(shù)。而應(yīng)用滾動時域控制后,在當(dāng)前制導(dǎo)指令更新周期內(nèi),飛行器再入過程中僅應(yīng)用偽譜法計算所得制導(dǎo)指令的前T秒。通過不斷的高超聲速飛行器不斷在線計算、更新,以此達(dá)到閉環(huán)制導(dǎo)效果。

4 仿真驗證與分析

高超聲速飛行器再入大氣層后采用無動力滑翔方式攻擊指定目標(biāo)。而飛行器再入過程中采用傾斜轉(zhuǎn)彎控制方法,因此再入過程中具有較強(qiáng)的側(cè)向機(jī)動能力和較大的側(cè)向位移。這就要求在設(shè)計高超聲速飛行器再入制導(dǎo)方法時,要充分考慮高超聲速飛行器再入過程中地球自轉(zhuǎn)以及飛行器傾側(cè)對再入過程制導(dǎo)效果的影響[7]。因此,建立高超聲速再入軌跡優(yōu)化模型時將考慮地球自轉(zhuǎn)以及飛行器傾側(cè)的影響,所以再入制導(dǎo)三自由度模型如下所示:

(15)

(16)

(17)

其中:ωe為地球旋轉(zhuǎn)角速度;σ為速度偏角,即為飛行器速度方向與正北方向的夾角;υ為飛行器側(cè)傾角,即飛行器升力方向與飛行器所處鉛錘面的夾角;λ為飛行器所處經(jīng)度值;φ為飛行器所處緯度值。

文章在研究高超聲速飛行器再入制導(dǎo)中,對飛行器橫向機(jī)動能力沒有特殊的要求,同時再入過程中未考慮飛行器禁飛區(qū)等因素對再入飛行器射程的影響。因此在高超聲速飛行器其再入過程中,其再入下降段和平衡滑翔段的橫向射程相較其縱向射程而言是一個小量[8]。所以再入飛行器射程可近似按下列公式計算,即:

(18)

其中:R0為地球半徑。

(19)

再入軌跡優(yōu)化過程約束為:

(20)

再入軌跡優(yōu)化性能指標(biāo)為:

(21)

應(yīng)用高斯偽譜法求解最優(yōu)控制問題時,選取飛行器初狀態(tài)參數(shù)為:

初始速度:V0=3 700 m/s;

初始高度:H0=70 km;

初始彈道傾角:Θ0=0°;

初始彈道偏角:ψv 0= 90°。

制導(dǎo)仿真要求末狀態(tài)為:

終端速度:Vf=1 000 m/s;

終端高度:Hf=0 km;

終端速度傾角:Θf=-80°;

終端射程:Lf=3 000 km。

高超聲速飛行器再入過程中為了保證在線軌跡優(yōu)化的快速性、實際制導(dǎo)指令的連續(xù)性以及對參數(shù)不確定的魯棒性,選取的制導(dǎo)指令更新周期應(yīng)本著綜合考慮再入制導(dǎo)指令計算時間、再入制導(dǎo)誤差影響等因素的思路,折中選取具有較高快速性、較小制導(dǎo)誤差、較強(qiáng)制導(dǎo)魯棒性的制導(dǎo)指令更新周期[10]。

因此在綜合考慮以上因素后,文章中選取滾動時域控制中制導(dǎo)指令更新周期為10 s。

應(yīng)用高超聲速飛行器再入初始狀態(tài)、終端狀態(tài)、過程約束等參數(shù),進(jìn)行高超聲速飛行器再入預(yù)測-校正制導(dǎo)方法仿真驗證。其仿真結(jié)果如圖1~5所示。

圖1 高度隨時間變化、經(jīng)緯度變化示意圖

圖2 速度、彈道傾角隨時間變化示意圖

圖3 彈道偏角、攻角隨時間變化示意圖

圖4 傾側(cè)角、射程隨時間變化示意圖

圖5 攻角導(dǎo)數(shù)、傾側(cè)角導(dǎo)數(shù)隨時間變化示意圖

對仿真結(jié)果進(jìn)行分析可知,應(yīng)用基于高斯偽譜法與滾動時域控制相結(jié)合的高超聲速飛行器再入制導(dǎo)方法,飛行器高度、速度隨時間變化平緩,再入過程中彈道無劇烈抖動情況出現(xiàn)。飛行器高度、速度參數(shù)滿足再入過程中對動壓、熱流約束的要求。而高超聲速飛行器再入過程中攻角、傾側(cè)角導(dǎo)數(shù)指令雖在圖中有較大調(diào)變,但是由于在設(shè)計時已考慮攻角、傾側(cè)角導(dǎo)數(shù)變化率對實際攻角、傾側(cè)角指令的影響,設(shè)計了其變化率最大值約束。所以即使攻角、傾側(cè)角導(dǎo)數(shù)變化劇烈,卻不影響實際再入制導(dǎo)指令攻角、傾側(cè)角的連續(xù)以及可實現(xiàn)性。通過對圖1~圖5彈道仿真結(jié)果的分析可知,高超聲速飛行器再入過程中,其彈道參數(shù)滿足過載約束要求。

綜上可知,應(yīng)用高斯偽譜法與滾動時域控制相結(jié)合的再入飛行器預(yù)測-校正制導(dǎo)方法,其制導(dǎo)指令具有較高的可實現(xiàn)性,飛行器彈道平緩,滿足過載、動壓、熱流以及平衡滑翔約束要求。

高超聲速飛行器再入過程終端狀態(tài)如表1所示。

表1 末端狀態(tài)情況

高超聲速飛行器再入預(yù)測-校正制導(dǎo)方法滿足再入制導(dǎo)過程滿足對過載、動壓以及熱流的約束限制,制導(dǎo)精度滿足高超聲速飛行器再入制導(dǎo)要求,再入過程可是實現(xiàn)對終端速度、速度傾角的控制。其中:速度控制誤差為10/s;速度傾角控制誤差為0.5°;射程控制誤差為2.7 km;可見應(yīng)用基于高斯偽譜法與滾動時域控制相結(jié)合的飛行器再入預(yù)測-校正制導(dǎo)方法,對飛行器再入終端狀態(tài)具有良好的控制效果,其終端參數(shù)控制誤差滿足設(shè)計要求,制導(dǎo)效果良好。

5 結(jié)論

文章以高超聲速飛行器再入制導(dǎo)問題為研究對象,通過建立高超聲速滑翔飛行器再入預(yù)測制導(dǎo)模型,再入飛行器過程約束模型以及飛行器再入終端狀態(tài)約束模型,應(yīng)用高斯偽譜法在線實時規(guī)劃高超聲速飛行器再入剩余彈道,并生成飛行器開環(huán)制導(dǎo)指令。其指令生成算法滿足高超聲速飛行器縱/側(cè)向、全程再入彈道一體化設(shè)計要求,制導(dǎo)指令具有較高連續(xù)性、較強(qiáng)可實現(xiàn)性。同時,文章結(jié)合滾動時域控制技術(shù),應(yīng)用基于偽譜法實時解算出的開環(huán)制導(dǎo)信息更新高超聲速滑翔飛行器再入制導(dǎo)指令,用分段開環(huán)的制導(dǎo)指令的方式達(dá)到高超聲速飛行器再入整體閉環(huán)制導(dǎo)的效果。經(jīng)彈道仿真驗證可知,這種基于高斯偽譜法與滾動時域控制相結(jié)合的高超聲速飛行器再入預(yù)測-制導(dǎo)方法,再入彈道平緩滿足飛行器再入過程約束要求,同時此制導(dǎo)方法對飛行器終端狀態(tài)具有良好的控制效果,其制導(dǎo)性能滿足再入飛行器對落點精度等性能指標(biāo)的要求。

[1] 李 強(qiáng). 高超聲速滑翔飛行器再入制導(dǎo)控制技術(shù)研究[D]. 北京: 北京理工大學(xué), 2015.

[2] 李鐵鵬. 基于高斯偽譜法的滑翔彈道優(yōu)化算法研究[J]. 彈箭與制導(dǎo)學(xué)報. 2014, 34(2): 113-116.

[3] 徐明亮, 陳克俊, 劉魯華,等. 高超聲速飛行器準(zhǔn)平衡滑翔自適應(yīng)制導(dǎo)方法[J]. 中國科學(xué):技術(shù)科學(xué). 2014(4).

[4] 張洪倩. 基于高斯偽譜法的彈道優(yōu)化設(shè)計與實現(xiàn)[D]. 南京: 南京理工大學(xué), 2014.

[5] Xie Y, Liu L, Liu J. Rapid Generation of entry Trajectories with Waypoint and No-fly zone Constraints[J]. Acta Astronautica. 2012: 167-181.

[6] Tian B, Zong Q. Optimal Guidance for Reentry Vehicles Based on Indirect Legendre Pseudospectral Method[J]. Acta Astronautica. 2011: 1176-1184.

[7] 陳小慶. 高超聲速滑翔飛行器機(jī)動技術(shù)研究[D]. 長沙: 國防科學(xué)技術(shù)大學(xué), 2011.

[8] 雍恩米. 高超聲速滑翔式再入飛行器軌跡優(yōu)化與制導(dǎo)方法研究[D]. 長沙: 國防科學(xué)技術(shù)大學(xué), 2008.

[9] 雍恩米, 唐國金, 磊 陳. 基于Gauss偽譜方法的高超聲速飛行器再入軌跡快速優(yōu)化[J]. 宇航學(xué)報. 2008, 29(6).

[10] F F, M R I. On discrete-time optimality conditions for pseudo-spectral methods[A]. AIAA/AAS Astrodynamics Specialist Conference and Exhibti[C]. 2006.

[11] M R I, F F. Pseudospectral knotting methods for solving optimal control problems[J]. Journal of Guidance,Control and Dynamics. 2004, 27(3): 397-405.

[12] 趙漢元. 飛行器再入動力學(xué)和制導(dǎo)[M]. 長沙: 國防科技大學(xué)出版社, 1997.

Research on Reentry Integrated Guidance Law of Hypersonic Velocity Aircraft Based on the Gauss Pseudo-spectral Method

Zhong Weikun1,Qu Quanyou1,Yuan Jinpeng1,Wu Wenbin2

(1.Beijing Institute of Spacecraft System Engineering,Beijing 100094, China;2.No.203 Research Institute of China Ordnance Industries, Xi′an 710065, China)

In order to solve the problems of the re-entry process of hypersonic velocity aircraft by using the standard trajectory method. Article gives a Gauss pseudo-spectral method to solve the re-entry problem of hypersonic velocity aircraft. The guidance combines the Gauss Pseudo-spectral Method and the Rolling-time Control Method. For the Gauss Pseudo-spectral Method which based on the integration idea that includes the vertical and the lateral side. And the integration includes All-trajectory integration also. By using the Rolling-time Control method, the guidance could realize the open-loop orders. As the result of the simulation, the data shows that after combine the two methods of Gauss Pseudo-spectral and the Rolling-time Control in hypersonic aircraft guidance control, the aircraft achieve a satisfactory result in re-entry trajectory.

hypersonic velocity aircraft;prediction-rectification;Gauss pseudo-spectral method;rolling-time control

2017-05-14;

2017-06-08。

仲維昆(1986-),男,山東曹縣人,碩士生,工程師,主要從事航天器項目管理、星務(wù)數(shù)據(jù)管理、姿軌控制等方向的研究。

1671-4598(2017)12-0106-04

10.16526/j.cnki.11-4762/tp.2017.12.028

V411.8

A

猜你喜歡
指令系統(tǒng)
聽我指令:大催眠術(shù)
Smartflower POP 一體式光伏系統(tǒng)
WJ-700無人機(jī)系統(tǒng)
ZC系列無人機(jī)遙感系統(tǒng)
北京測繪(2020年12期)2020-12-29 01:33:58
基于PowerPC+FPGA顯示系統(tǒng)
半沸制皂系統(tǒng)(下)
ARINC661顯控指令快速驗證方法
LED照明產(chǎn)品歐盟ErP指令要求解讀
電子測試(2018年18期)2018-11-14 02:30:34
連通與提升系統(tǒng)的最后一塊拼圖 Audiolab 傲立 M-DAC mini
殺毒軟件中指令虛擬機(jī)的脆弱性分析
主站蜘蛛池模板: 无码中字出轨中文人妻中文中| 国产xx在线观看| 亚洲精品无码日韩国产不卡| 午夜久久影院| 国产亚洲精品自在久久不卡| 欧美日韩另类国产| 99视频在线精品免费观看6| 亚洲动漫h| 99久久人妻精品免费二区| 久久女人网| 婷婷开心中文字幕| 国产色爱av资源综合区| 精品三级在线| 欧美啪啪网| 精品人妻AV区| 自慰高潮喷白浆在线观看| 国产九九精品视频| 欧美一级大片在线观看| 高清免费毛片| 在线亚洲精品福利网址导航| 亚洲床戏一区| 亚洲日本韩在线观看| 国产网站黄| 亚洲成人精品久久| 中文字幕乱码中文乱码51精品| 亚洲欧州色色免费AV| 国产精品香蕉在线观看不卡| 久久久噜噜噜| 国产女同自拍视频| 久久久噜噜噜| 亚洲欧美一区二区三区蜜芽| 成人免费视频一区二区三区 | 亚洲av中文无码乱人伦在线r| 国产午夜福利亚洲第一| 国产波多野结衣中文在线播放| 一级在线毛片| 九色视频线上播放| 一级在线毛片| 少妇高潮惨叫久久久久久| 狠狠躁天天躁夜夜躁婷婷| 国产在线观看人成激情视频| aa级毛片毛片免费观看久| 日韩精品亚洲一区中文字幕| 亚洲免费三区| 秋霞一区二区三区| 无码中文字幕乱码免费2| 欧美日韩动态图| 国产成人无码AV在线播放动漫| 午夜国产理论| 日韩 欧美 小说 综合网 另类| 很黄的网站在线观看| 久草视频精品| 欧美在线伊人| 91久久精品日日躁夜夜躁欧美| 亚洲,国产,日韩,综合一区| 国产网站免费| 再看日本中文字幕在线观看| 97se亚洲综合在线天天| 国产精品密蕾丝视频| 无码区日韩专区免费系列| 亚洲视频欧美不卡| 国产欧美日韩综合在线第一| 91精品情国产情侣高潮对白蜜| 国产精品无码制服丝袜| 久久精品娱乐亚洲领先| 999国产精品| 在线观看精品自拍视频| 亚洲无码A视频在线| 日韩人妻精品一区| 好吊日免费视频| 国产成人精品日本亚洲77美色| 欧洲极品无码一区二区三区| 91成人精品视频| 91精品最新国内在线播放| 无码视频国产精品一区二区| 欧美三级不卡在线观看视频| 综合天天色| 国产精品va| 国内丰满少妇猛烈精品播| 欧美在线观看不卡| 99久久99视频| 激情国产精品一区|