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

高馬赫數飛行器準平衡飛行段彈道優化方法

2023-10-10 07:25:04叢戎飛丁智堅
兵器裝備工程學報 2023年9期
關鍵詞:優化

溫 杰,周 歡,叢戎飛,丁智堅,張 旭

(1.中國工程物理研究院 總體工程研究所,四川 綿陽 621999;2.中國空氣動力研究與發展中心 空天技術研究所,四川 綿陽 621000)

0 引言

高馬赫數飛行器因具備高機動、快速打擊優勢成為目前各軍事大國競相研究的熱點[1-2],準平衡飛行段彈道優化是突破敵方反導攔截系統的關鍵。高馬赫數飛行器準平衡飛行過程中不僅要考慮動壓、過載、熱流率等約束限制,還需要對多禁飛區(自然環境、政治軍事因素形成的)進行規避[3]。現有針對禁飛區的研究成果相對較淺,主要通過彈道優化或設計橫向制導邏輯規避禁飛區。

在彈道優化方面,高馬赫數飛行器準平衡飛行段彈道優化問題具有較高的非線性和復雜度,經典變分法和Pontryagin極小值原理難以高效求得最優解,因此通常用直接法進行求解[4]。直接打靶法是直接法中典型的代表方法,其通過離散控制變量,將最優控制問題轉化為非線性規劃問題[5],具有通用性好、易實現等優點。但直接打靶法的收斂性能依賴于終端時刻和控制量初末時刻猜想值,在工程實際應用中初值難以預估,需要進行相應的改進[6]。在橫向制導邏輯方面,方法主要包括觸角法[7]、傾側角符號反轉法[8]、人工勢場法[9]等,其通過在線設計制導策略來規避禁飛區,但并未考慮終端時刻性能指標的最優性(如最遠航程等)。

遺傳算法是一種模擬自然選擇和生物進化過程搜索全局最優的方法[10-11],因其具有較強的魯棒性,對初值不敏感,對于大型、復雜非線性系統的全局優化問題具有良好的性能,在彈道優化領域得到了廣泛的應用。如Cai等[12]利用遺傳算法優化了以航程為目標函數的預測模型的飛行路徑;Dancila等[13]提出一種基于遺傳算法的新型飛行軌跡優化方法,選擇優化標準為總成本最小,候選軌跡定義為橫向和縱向兩部分;Ergezer[14]開發了一種新的多目標遺傳算法來解決路徑規劃,優化目標為軌跡長度、完成任務所需的時間以及耗能。但需要強調的是,遺傳算法中適應度函數的選取直接關系到收斂速度以及是否能尋找到最優解[15],文獻[16]在適應度函數中引入自適應調節機制,設計了一種適應度函數可隨種群迭代次數增加而自適應調整的自適應遺傳算法,以此來增加種群多樣性,進一步快速、準確的逼近全局最優解。

基于以上分析,本文中將控制量初末時刻值納入優化設計變量,以終端約束中復合約束項(如速度和高度等)為判斷條件,提出一種基于改進直接打靶法和自適應遺傳算法的混合優化算法,并應用到高馬赫數飛行器準平衡飛行段彈道優化問題中。通過仿真計算控制量及狀態量,評估彈道性能,得出結論。

1 準平衡飛行段彈道優化問題

1.1 高馬赫數飛行器運動學模型

忽略地球自轉以及曲率影響,將地球假設為平面,分別在彈道坐標系和地面坐標系建立動力學方程和運動學方程,得到3自由度高馬赫數飛行器質心運動方程組簡化模型[17]:

(1)

式(1)中:V(t)為飛行器速度,θ(t)和ψv(t)分別為彈道傾角和彈道偏角;y(t)為高度,x(t)和z(t)分別為縱程和橫程;m(t)為飛行器質量,g為重力加速度,α(t)為攻角,γv(t)為速度傾斜角。αB為平衡攻角,βB為平衡側滑角;YB和ZB分別為αB、βB所對應的平衡升力和平衡側向力,X為阻力,計算表達式如下

(2)

式(2)中:Cx、Cy為阻力系數和升力系數,q=0.5*ρ*V2,其中ρ為高馬赫數飛行器所處高度的空氣密度,S為飛行器參考面積。

1.2 約束條件

1) 控制約束

為保證飛行器穩定飛行,控制量攻角α和速度傾斜角γv要限制在一定的范圍內:

(3)

2) 過程約束

① 熱流率、動壓、過載約束。一般而言,高馬赫飛行器在準平衡飛行過程中必須滿足熱流率、過載、動壓等硬約束,公式表示如下

(4)

② 禁飛區約束。禁飛區是因自然環境、政治軍事因素形成的不允許飛入的區域,本文中禁飛區選取典型規則形狀圓柱形,中心點為(xi,yi,zi),半徑為ri,具體表達式如下

(5)

設L(x(t),t)是飛行器到禁飛區的最短距離,則禁飛區約束可表示為L(x(t),t)≥εn,其中εn為無窮小正常量[3]。

3) 終端約束。為滿足中末制導交接班需求,需對準平衡飛行段終端時刻的高度、速度、彈道傾角設定約束

(6)

式(6)中:yf、Vf、θf分別為終端時刻的高度、速度和彈道傾角;ydf、Vdf、θdf分別為預定的終端時刻的高度、速度和彈道傾角;εh、εv、εθ分別為3個約束量的誤差上界。

1.3 目標函數

一般而言,目標函數根據具體任務選定,一般可分為最大航程(最大縱程/最大橫程)、飛行時間最短、終端時刻速度最大等。本文中選擇最大航程(最大縱程/最大橫程)為優化目標,目標函數形式如下

minJ=-Ldf

(7)

式(7)中:Ldf為飛行器的航程。

1.4 問題描述

基于以上分析,高馬赫數飛行器準平衡飛行段彈道優化問題可描述為:在滿足質心運動方程組式(1)、控制約束式(3)、過程約束式(4)、式(5)、終端約束式(6)的條件下,尋找控制量攻角α和速度傾斜角γv最優時間序列使得目標函數式(7)最小。更一般性的最優控制問題描述如下:尋找控制變量u(t),使得Bolza型性能指標(即Mayer型和Lagrange型的組合)最小,如式(8)所示:

(8)

且滿足質心運動方程組約束式(9):

(9)

邊界條件約束式(10):

φ(x(t0),t0,x(tf),tf)=0

(10)

以及不等式約束(過程約束)式(11):

C(x(t),u(t),t)≤0

(11)

2 改進直接打靶法——自適應遺傳算法

2.1 改進直接打靶法

直接打靶法是一種僅離散控制變量的方法,描述運動軌跡的狀態變量需要根據參數化后的控制變量對運動方程組數值積分獲得,對于高馬赫數飛行器準平衡飛行段彈道優化問題,一般選擇四階龍格庫塔法(Runge-Kutta)進行求解,以保證數值計算的精度和穩定性。需要注意的是,直接打靶法對初值要求較高,算法的精確性往往依賴于初值猜想的準確程度,可能產生局部極小的問題,為此,對直接打靶進行改進,將控制量初末時刻值一并納入優化設計變量,以終端約束中復合約束項(如速度和高度等)為判斷條件,將彈道結束時刻作為終端時刻,以此達到降低對初值敏感的目的。

1) 時間變量離散

在彈道優化問題中,一般選擇單調變化物理量為參考自變量,直接打靶法中一般選擇對時間進行離散,得到離散時間序列T,如式(12)所示:

t0=t1

(12)

本文中以終端約束中復合約束項(如速度和高度等)為判斷條件,則tN可取一個較大的值,如tN=1 000 s,可令終端約束單調變化量結束時刻為tk=tf。

2) 控制變量離散

以離散節點處的控制變量為設計變量,得到離散控制序列U,如式(13)所示:

[u1,u2,…,uk-1,uk=uf,…,uN-1,uN]

(13)

時間節點之間的控制量通過插值函數ψi(t)獲取,因此,連續空間的控制變量做如下近似

(14)

式(14)中:ui為網格節點處的控制變量。

3) 狀態變量離散

高馬赫數飛行器初始狀態是已知的,即Xi,在通過優化算法得到離散控制序列U的前提下,可通過數值積分的方法,經過多次迭代積分,即可得到與時間序列對應的狀態變量序列X:

[X1,X2,…,Xk-1,Xk=Xf,…,XN-1,XN]

(15)

其中,在高馬赫數飛行器的彈道優化中,數值積分一般選用四階龍格庫塔法(Runge-Kutta)以保證計算精度。式(15)中:Xk=Xf為終端時刻tf對應的終端狀態。

4) 最優控制問題轉化

令設計變量集合如式(16)所示:

Y=[X1,u1,t1,X2,u2,t1,…,Xk-1,uk-1,tk-1,Xf,uf,tf]

(16)

式(8)所示的性能指標可以表示如下:

(17)

式(10)所示的邊界條件可表示為

g1(Y)=0

(18)

式(11)所示的過程約束可表示為

g2(Y)≤0

(19)

綜合式(17)—式(19),可以得到,通過改進直接打靶法彈道優化問題已經轉化為經典的非線性規劃問題,如式(20)所示:

(20)

2.2 自適應遺傳算法

遺傳算法是一類借鑒生物界自然選擇和自然遺傳機制的搜索全局最優的方法,解決了眾多復雜工程問題(如多目標優化等),在科學研究及工程實踐中得到了廣泛的應用[14,18]。遺傳算法經種群初始化編碼,通過選擇、交叉、變異操作,實現優勝劣汰不斷向最優性狀進化。自適應遺傳算法在遺傳算法的基礎上,可在適應度函數、交叉概率、變異概率等幾個方面進行改進。考慮遺傳算法中適應度函數的選取直接關系到收斂速度以及是否能尋找到最優解,故本文中主要針對遺傳算法的適應度函數作出如下改進:① 適應度函數隨著種群迭代次數的增加而自適應調整,以增加種群多樣性,進一步快速、準確的逼近全局最優解。② 適應度函數進行對數變換,縮小性能指標和終端約束條件之間參數量級差距,使其更好的反應個體的優劣。③ 引入多目標分層規劃方法,將多個目標按重要程度分成多個優先層次,以保證終端時刻滿足終端約束的同時使性能指標最優。

具體步驟如下:

1) 種群初始化編碼

自適應遺傳算法首先需要進行初始化編碼,將所求問題可行解表示為遺傳空間的染色體或個體。本文中首先在控制約束范圍內構造控制量數據集,基于改進直接打靶法按照式(12)選取若干時間點,每個節點時刻的控制量攻角和速度傾斜角構成的2個一維數組α=(α1,α2,…,αN-1,αN)、γv=(γv1,γv2,…,γvN-1,γvN)作為彈道個體的染色體,選擇這種編碼有幾下優點:

① 控制約束式(3)可以自動滿足;

② 在3自由度彈道中,當控制量攻角和速度傾斜角歷程確定后可直接得到完整的彈道,若選取其他參數如彈道傾角、高度等為控制量,在迭代時都需要計算出攻角和速度傾斜角,增加了計算的時間和復雜度。

2) 插值函數

在給定的氣動數據中攻角數量是有限的,需要經過插值得到更多的數據量,另外僅有節點時刻攻角值是不夠的,需進行插值計算。插值的方法主要有線性插值、多項式插值、樣條插值等。本文中采用三次樣條插值,原因如下:

① 控制約束式(3)可以自動滿足;

② 算法易實現,與線性插值相比,插值誤差小于線性插值;使用樣條會比使用高階多項式更容易評估,不會受到龍格現象的影響;

③ 通過樣條插值得到的攻角是一條光滑的曲線,不是線性插值的折線,利于飛行器控制系統的設計。

3) 適應度函數

(21)

式(21)中:R1=(cd)λ,fi(x)為違反第i個終端約束的懲罰函數項,pi為第i個終端約束懲罰項的權重系數,c、λ、?均為常數,d為迭代次數。J為彈道優化目標函數項,q為相應的權重系數,K為常數(保證適應度為正即可)。經數值仿真驗證,取c=0.5,λ=?=2是比較可行的方案。另外在綜合適應度函數中,終端約束項之間經過對數處理差距較小,因此主要調節下一層級性能指標權重系數q即可。

4) 選擇操作

從種群中選擇優勝個體,淘汰劣質個體的操作叫做選擇。選擇操作的目的是將優化個體直接遺傳給下一代或者通過配對交叉產生新個體再遺傳給下一代,通常選擇輪盤賭法,設種群規模為Q,其中個體i的適應度為Fi,則第i個個體被選擇的概率如式(22)所示:

(22)

5) 交叉操作

交叉操作是體現自適應遺傳算法全局搜索能力強的核心步驟,其根據交叉概率將兩個染色體交叉組合,以期望父代優秀個體優秀特征遺傳給子代,產生新的優秀個體。本文中選用實數交叉法,即第k個染色體ak和第l個染色體al在j位的交叉操作為

(23)

式(23)中:b為[0,1]的隨機數。

6) 變異操作

變異操作的目的有2個,一是增強自適應遺傳算法的局部尋優能力,二是維持種群多樣性,以防止未成熟收斂現象。第i個個體的第j個基因aij變異操作為

(24)

式(24)中:amax、amin分別為基因aij的上界和下界;f(g)=r2(1-g0/Gmax)2,r2是一個隨機數,g0是當前迭代次數,Gmax為最大進化次數,r是[0,1]的隨機數。

2.3 混合優化算法

綜前文所述,本文中提出的高馬赫數飛行器準平衡飛行段彈道優化方法為混合優化算法,即改進直接打靶法-自適應遺傳算法,主要包含以下步驟:

1) 利用改進直接打靶法將原彈道優化問題(原動態最優控制問題)轉化為非線性規劃問題;

2) 在控制約束范圍內構造控制量攻角α和速度傾斜角γv數據集;

3) 將初始狀態X0作為輸入,將控制量初末時刻值一并納入優化設計變量,利用自適應遺傳算法優化改進直接打靶法中離散點處的控制量,并通過三次樣條插值對控制量-時間歷程平滑處理,以終端約束中復合約束項(如速度和高度等)為判斷條件,基于四階龍格庫塔法(Runge-Kutta)進行數值積分,經過若干次迭代計算出最優控制u*,同時可以得到狀態變量曲線最優軌跡X*。

算法思路框圖如圖1、圖2所示。

圖1 改進直接打靶法求解思路

圖2 自適應遺傳算法流程圖

3 仿真算例與分析

3.1 仿真相關參數

飛行器和地球基本參數如表1所示。控制約束如表2所示。過程約束如表3所示。邊界約束(初始狀態和終端狀態)如表4所示。

表1 仿真基本參數

表2 控制約束

表3 過程約束

表4 邊界約束(初始狀態和終端狀態)

自適應遺傳算法主要參數如表5所示。

表5 自適應遺傳算法主要參數

本文中設計的自適應遺傳算法主要針對遺傳算法中適應度函數進行改進,其他參數可在以下范圍內進行設定:種群規模:20~100;迭代次數:100~500;交叉概率:0.4~0.99;變異概率:0.000 1~0.1[19]。為防止早熟,一般采用增大變異率或增加種群數量的方式來解決,故變異率可選范圍內的最大值,種群數量可以適當的增加,迭代次數選擇范圍內最低值以減少計算量,交叉操作是體現自適應遺傳算法全局搜索能力的核心,可取較大值,具體如表5所示。

為了驗證構建的改進直接打靶法-自適應遺傳算法的魯棒性,將實際飛行過程中氣動參數視為不確定性較為合理[20]。不失一般性,在原氣動參數的基礎上,引入正態分布白噪聲,具體公式如下

(25)

3.2 最遠縱程彈道優化

在以最遠縱向航程為目標的彈道優化問題中,縱向航程對應x的變化,則性能指標如式(26)所示:

J0=min(-L0)=min(-x(tf))

(26)

利用原算法(即直接打靶法-遺傳算法,需要不斷試測初值,另外適應度函數需進行對數變換,否則較難收斂)、構建的混合算法:改進直接打靶法-自適應遺傳算法以及在此基礎上引入白噪聲進行仿真,結果如圖3—圖17所示。

圖3 最佳適應度變化曲線

圖5 彈道傾角變化曲線

圖6 彈道偏角變化曲線

圖7 縱程變化曲線

圖8 高度變化曲線

圖9 橫程變化曲線

圖10 攻角變化曲線

圖11 速度傾斜角變化曲線

圖12 高度-縱程變化曲線

圖13 熱流率變化曲線

圖14 動壓變化曲線

圖15 過載變化曲線

圖16 規避禁飛區俯視圖

圖17 規避禁飛區三維立體圖

3.3 最遠橫程彈道優化

在以最遠橫向航程為目標的彈道優化問題中,橫向航程對應z的變化,則性能指標如式(27)所示:

J0=min(-L0)=min(-z(tf))

(27)

利用原算法(即直接打靶法-遺傳算法,需要不斷試測初值,另外適應度函數需進行對數變換,否則較難收斂)、構建的混合算法:改進直接打靶法-自適應遺傳算法以及在此基礎上引入白噪聲進行仿真,結果如圖18—圖32所示。

圖18 最佳適應度變化曲線

圖19 速度變化曲線

圖20 彈道傾角變化曲線

圖21 彈道偏角變化曲線

圖23 高度變化曲線

圖24 橫程變化曲線

圖3和圖18給出了種群中最佳適應度值隨迭代次數變化的過程,由圖可知,原遺傳算法對應的適應度函數在40代左右趨于收斂,構造的自適應算法對應的最佳適應度值在初期變化較大,隨著迭代次數的增加而逐漸收斂,在10代左右已經趨于平穩,后期在60~90代左右混合算法優化控制量最佳適應度函數值有些許變化,但變化較小,相較于原遺傳算法收斂速度較快。這表明利用自適應遺傳算法對控制變量尋優具有較好的收斂效果。

需要說明的是,原算法(直接打靶法-遺傳算法)是通過大量試測初值得到的,在大部分情況下由于對初值敏感,往往得不到規避禁飛區的軌跡。而構造的混合算法將控制變量初末時刻值納入優化設計變量,以終端約束中單調變化量(如速度)作為軌跡結束時刻的判斷條件,初值不需要不斷試測,效率較高。從圖10—圖11、圖25—圖26中,可以看出,構建的混合算法:改進直接打靶法-自適應遺傳算法在以最遠縱程/最遠橫程為性能指標的仿真中,控制約束在規定的范圍內;圖13—圖15、圖28—圖30中,過程約束中一般約束項熱流率、過載、動壓也在約束范圍內;圖16—圖17、圖31—圖32中,能直觀看出,混合算法能夠實現規避多禁飛區,滿足了實際需求。

圖25 攻角變化曲線

圖26 速度傾斜角變化曲線

圖27 高度-縱程變化曲線

圖29 動壓變化曲線

圖30 過載變化曲線

圖31 規避禁飛區俯視圖

圖32 規避禁飛區三維立體圖

最后,在考慮模型不確定性因素的情況下,即氣動系數引入正態分布的白噪聲后,最佳適應度函數曲線在20~30代左右已經趨于收斂,自適應遺傳算法對控制變量尋優仍具有較好的收斂效果,相應的控制約束、過程約束(熱流率、過載、動壓及多禁飛區)等仍能滿足相應的要求,驗證了構建的混合算法具備一定的魯棒性。

為了更直觀反應構建的混合算法滿足終端約束,建立如表6、表7所示。

表6 終端約束條件對應的終端狀態(性能指標:最遠縱程)

表7 終端約束條件對應的終端狀態(性能指標:最遠橫程)

從圖4、圖5、圖7、圖8和圖19、圖20、圖23、圖24以及表6、表7可以看出,在以最遠縱程/最遠橫程為性能指標的仿真中,混合算法能夠滿足終端約束條件,且相較于原算法性能指標最遠縱程/最遠橫程數值更大,更接近全局最優解,驗證了通過引入多目標分層規劃方法以保證在終端時刻滿足終端約束的同時使性能指標最優的有效性。另外,在考慮氣動參數受到干擾的仿真中,終端約束仍能滿足要求。

綜上所述,構造的基于改進直接打靶法和自適應遺傳算法得到的最優性能指標(最遠縱程/最遠橫程)軌跡滿足控制約束、過程約束(動壓、過載、熱流率及禁飛區)以及終端約束條件,能夠規避多禁飛區,說明了構建的混合算法的有效性。其次,相較于原算法(直接打靶法-遺傳算法),混合算法不需要不斷試測初值,收斂速度較快,性能指標更優,能夠快速、準確的逼近全局最優解。最后,在考慮模型不確定因素,即氣動參數受到正態分布白噪聲干擾時仍能滿足要求,說明了構建的混合算法具備一定的魯棒性。

4 結論

高馬赫數飛行器準平衡飛行段彈道優化方法是高馬赫數飛行器研究中主要的關鍵技術,考慮飛行過程中面臨包含多禁飛區等復雜約束限制,提出了一種基于改進直接打靶法和自適應遺傳算法的混合優化解決方案。

1) 通過對直接打靶法改進,將控制量初末時刻值一并納入優化設計變量,以終端約束中復合約束項(如速度和高度等)為判斷條件,將彈道結束時刻作為終端時刻,以此降低了對初值的敏感程度。

2) 為進一步快速、準確逼近全局最優解,對遺傳算法的適應度函數進行了三點改進,構造了可隨種群迭代次數增加而自適應調整的對數形式的綜合適應度函數,建立了新的自適應遺傳算法。

3) 基于建立的自適應遺傳算法優化改進直接打靶法中離散點處的控制量,通過三次樣條插值對控制量-時間歷程平滑處理,并利用四階龍格庫塔法(Runge-Kutta)進行數值積分,經過若干次迭代計算出最優控制時間序列,同時得到了能夠實現規避多禁飛區,滿足約束條件且保證航程最遠的理想軌跡。

試驗表明,建立的混合優化算法不需要不斷試測初值,相較于原算法(直接打靶法-遺傳算法)收斂速度更快,性能指標最遠縱程/最遠橫程數值更大,最接近全局最優解。另外,考慮實際飛行過程中模型存在的不確定性因素,在氣動參數中引入了正態分布白噪聲,結果表明準平衡飛行過程中各類約束條件仍能滿足要求,說明了構造的混合算法具備一定的魯棒性。在應用方面,構建的混合優化算法本質上是一種優化方法,不僅適用于高馬赫數飛行器準平衡飛行段彈道優化,也適用于一般的飛行器彈道優化問題,為下一步研究不確定性量化的魯棒彈道優化問題提供一定的參考價值。

猜你喜歡
優化
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
PEMFC流道的多目標優化
能源工程(2022年1期)2022-03-29 01:06:28
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
圍繞“地、業、人”優化產業扶貧
今日農業(2020年16期)2020-12-14 15:04:59
事業單位中固定資產會計處理的優化
消費導刊(2018年8期)2018-05-25 13:20:08
4K HDR性能大幅度優化 JVC DLA-X8 18 BC
幾種常見的負載均衡算法的優化
電子制作(2017年20期)2017-04-26 06:57:45
主站蜘蛛池模板: 亚洲最黄视频| 国产va在线观看| 青青热久免费精品视频6| 国产欧美日韩va| 波多野结衣爽到高潮漏水大喷| 国产精品护士| 国产精品嫩草影院视频| 九九热视频精品在线| 九九热精品在线视频| 欧美国产菊爆免费观看| 97国产精品视频自在拍| 久久久久国产精品免费免费不卡| 日韩av资源在线| 国产精品福利尤物youwu| 国产成人一区| 一本大道香蕉中文日本不卡高清二区 | 国产SUV精品一区二区| 亚洲精品少妇熟女| 国产精品性| 天天做天天爱夜夜爽毛片毛片| 香蕉久人久人青草青草| 综合社区亚洲熟妇p| 亚洲成aⅴ人在线观看| 国产精品亚洲一区二区三区z| 亚洲精品无码久久毛片波多野吉| 五月丁香伊人啪啪手机免费观看| 少妇露出福利视频| 日韩欧美在线观看| 国产青榴视频| 久热中文字幕在线观看| 国产精品一老牛影视频| 国产精品人成在线播放| 就去色综合| 欧美日韩国产精品va| 99er精品视频| 精品无码日韩国产不卡av| 亚洲欧美自拍中文| 国产成人亚洲无吗淙合青草| 日本亚洲欧美在线| 国产嫩草在线观看| 天堂成人av| 久久久久人妻精品一区三寸蜜桃| 国产精品欧美在线观看| 国产99免费视频| 国产色爱av资源综合区| 国产成人AV综合久久| 全裸无码专区| 在线观看国产精品第一区免费| 亚洲无卡视频| 狠狠色丁香婷婷| 欧美全免费aaaaaa特黄在线| 波多野结衣一区二区三区AV| 国产亚洲视频免费播放| a亚洲视频| a天堂视频| 久久无码免费束人妻| 欧美不卡二区| 欧美在线网| 免费人成视网站在线不卡| 国产玖玖玖精品视频| 五月天在线网站| 久久国产亚洲偷自| 国产簧片免费在线播放| 国产97区一区二区三区无码| 在线欧美a| 日韩国产黄色网站| 99视频在线精品免费观看6| 国模视频一区二区| 91精品啪在线观看国产| 日韩成人免费网站| 黄色在线不卡| 特级aaaaaaaaa毛片免费视频| 男女精品视频| 精品国产Av电影无码久久久| 中字无码精油按摩中出视频| 91麻豆精品视频| 成人在线观看一区| 国产免费羞羞视频| 免费女人18毛片a级毛片视频| 99在线免费播放| 亚洲国产天堂在线观看| 国产精品视频导航|